fracspy.location.migration.diffstack#
- fracspy.location.migration.diffstack(data, n_xyz, x, y, z, tt, dt, nforhc=10, output_type=None, stack_type=None, swsize=0, polcor_type=None, recs=None)[source]#
Diffraction stacking for microseismic source location.
This routine performs imaging of microseismic data by diffraction stacking with various stacking approaches and optional polarity correction.
- Parameters:
- data
numpy.ndarray Data of shape \(n_r \times n_t\)
- n_xyz
tuple Number of grid points in X-, Y-, and Z-axes for the imaging area
- x
numpy.ndarray Imaging area grid vector in X-axis
- y
numpy.ndarray Imaging area grid vector in Y-axis
- z
numpy.ndarray Imaging area grid vector in Z-axis
- tt
numpy.ndarray Traveltime table of size \(n_r \times n_x \times n_y \times n_z\)
- dt
float Time step
- nforhc
int, optional, default: 10 Number of points for hypocenter
- output_type
str, optional, default: None Defines the output type. In the default case of None it is set to
max. Types:max(output 3D volume as a maximum of the image function through time),mean(output 3D volume as an average of the image function through time),sumsq(output 3D volume as a stack of the squared image function values through time),full(output full 4D image function: x,y,z,t). In this case hc is set to None.- stack_type
str, optional, default: None Diffraction stacking type (imaging condition), default None is the same as
squared. Types:absolute(absolute-value),squared(squared-value),semblance(semblance-based).- swsize
int, optional, default: 0 Sliding time window size for semblance-based type, amount of time steps
- polcor_type
str, optional, default: None Polarity correction type to be used for amplitudes. None is default: no polarity correction,
mtiis for polarity correction based on moment tensor inversion.- recs
numpy.ndarray, optional, default: None Array of shape (3, nrec) containing receiver coordinates. Must be provided if polcor_type is not None
- data
- Returns:
- ds_im_vol
numpy.ndarray Diffraction stack volume, returned if output_type is not “full”
- hc
numpy.ndarray Estimated hypocentral location, if output_type is not “full” is set to None
- ds_full
numpy.ndarray Diffraction stack full array of shape (nt, n_xyz), returned if output_type is “full”
- ds_im_vol
- Raises:
- ValueError
if data and grid sizes do not fit the traveltimes array
- ValueError
if output_type value is unknown
- ValueError
if stack_type value is unknown
- ValueError
if polcor_type value is unknown
- ValueError
if recs is None and polcor_type is not None
Notes
The subsurface volume is discretised and each grid node is considered to be a potential source position or a so-called image point. In other words, each image point represents a possible diffraction point from which seismic energy radiates. The term “diffraction stacking” dates back to the works of Claerbout (1971, 1985) [3], [4] and Timoshin (1972) [7]. The concept was initially related to seismic migration and imaging of reflectors as a set of diffraction points (in the context of the exploding reflector principle). Here, the diffraction stacking is related to imaging of real excitation sources.
Diffraction stacking is applied to data with corrected event moveout (EMO), and can be done in several ways:
Absolute-value based
\[F(\mathbf{r},t) = \left| \sum_{R=1}^{N_R} A_R^{EMO}(t,\mathbf{r}) \right|,\]where \(\mathbf{r}\) is a vector that defines a spatial position \((x, y, z)\) of the image point, \(A_R^{EMO}(t,\mathbf{r})\) represents the EMO-corrected data at the receiver \(R\), and \(N_R\) is a number of receivers [1].
Squared-value based
Similarly, but the stack is taken to the power of 2:
\[E(\mathbf{r},t) = \left[ \sum_{R=1}^{N_R} A_R^{EMO}(t,\mathbf{r}) \right]^2,\]Semblance-based
Simple stacking can be further improved using a semblance-based approach.
The semblance is a coherency or similarity measure and can be understood as the ratio of the total energy (the square of sum of amplitudes) to the energy of individual traces (the sum of squares) [5]:
\[S(\mathbf{r},t) = \frac{\left[ \sum_{R=1}^{N_R} A_R^{EMO}(t,\mathbf{r}) \right]^2} {N_R \sum_{R=1}^{N_R} \left[ A_R^{EMO}(t,\mathbf{r}) \right]^2} \]In order to suppress the effect of noise even better, it is possible to extend the semblance-based approach by introducing a sliding time window \(W\) over which the energy measures are summed:
\[S_W(\mathbf{r},t,W) = \frac{\sum_{k=it-W}^{it+W}\left[ \sum_{R=1}^{N_R} A_R^{EMO}(t,\mathbf{r}) \right]^2} {N_R \sum_{k=it-W}^{it+W}\sum_{R=1}^{N_R} \left[ A_R^{EMO}(t,\mathbf{r}) \right]^2}\]where \(k\) an index of the time-discretised signal within a sliding time interval consisting of the \(2W + 1\) samples, and \(it\) is the index of time \(t\) [8].
Depending on the desired output, the function can output either a 3D volume and a location, or a full 4D array (time dimension + 3 spatial dimensions).
By default, a 3D image volume obtained as a maximum of the 4D image function, e.g.:
\[F_s(\mathbf{r}) = \max_t F(\mathbf{r},t).\]Optionally, one can output a 3D image volume obtained as a mean of the 4D image function over time, e.g.:
\[F_s(\mathbf{r}) = \underset{t}{\text{mean}}~F(\mathbf{r}, t).\]Another option is to output a 3D image volume obtained as a sum of squared values of the 4D image function over time, e.g.:
\[F_s(\mathbf{r}) = \sum_t F^2(\mathbf{r},t).\]The full 4D output is useful for subsequent detection of multiple events.
Polarity correction, if requested, is done using moment tensor inversion [2]:
\[F_{pc}(\mathbf{r},t) = \left|\sum_{R=1}^{N_R} \phi_R(\mathbf{r}) A_R^{EMO+PC} \right|,\]where
\[A_R^{EMO+PC}(t,\mathbf{r}) = \phi_R(\mathbf{r}) A_R \left(t + T_R(\mathbf{r})\right),\]and \(\phi_R(\mathbf{r})\) are the defined through the \(\mathrm{sign}\) function:
\[\phi_R(\mathbf{r})=\mathrm{sign}(\mathbf{M}(\mathbf{r})\mathbf{G}_R^{\phantom{T}}(\mathbf{r})),\quad R=1\dots N_R,\]where \(\mathbf{M}(\mathbf{r})\) is the vectorized moment tensor derived [6] as
\[\mathbf{M}(\mathbf{r}) = \left(\sum_{R=1}^{N_R} \mathbf{G}_R^T(\mathbf{r})\mathbf{G}_R^{\phantom{T}}(\mathbf{r})\right)^{-1} \mathbf{d}(\mathbf{r}),\]where \(\mathbf{G}_R^{\phantom{T}}(\mathbf{r})\) is the vectorized derivative of Green’s function for receiver \(R\) and \(\mathbf{d}(\mathbf{r})\) is a vector obtained from amplitudes of the corresponding waveforms:
\[\mathbf{d}(\mathbf{r}) = \sum_{R=1}^{N_R} A_R(t+T_R(\mathbf{r})) \mathbf{G}_R^{\phantom{T}}(\mathbf{r}).\]More details on the described inversion including stability issues can be found in Zhebel & Eisner (2014) [9].
References
[1]Anikiev, D. (2015). Joint detection, location and source mechanism determination of microseismic events (Doctoral dissertation). St. Petersburg State University. St. Petersburg. https://disser.spbu.ru/files/phdspsu2015/Anikiev_PhD_web_final.pdf
[2]Anikiev, D., Valenta, J., Staněk, F. & Eisner, L. (2014). Joint location and source mechanism inversion of microseismic events: Benchmarking on seismicity induced by hydraulic fracturing. Geophysical Journal International, 198(1), 249–258. https://doi.org/10.1093/gji/ggu126
[3]Claerbout, J. F. (1971). Toward a unified theory of reflector mapping. Geophysics, 36(3), 467-481. https://doi.org/10.1190/1.1440185
[4]Claerbout, J. (1985). Imaging the earth’s interior. Oxford, England: Blackwell Scientific Publications. https://sepwww.stanford.edu/sep/prof/iei2/
[5]Neidell, N. S., & Taner, M. T. (1971). SEMBLANCE AND OTHER COHERENCY MEASURES FOR MULTICHANNEL DATA. Geophysics, 36(3), 482–497. https://doi.org/10.1190/1.1440186
[6]Sipkin, S. A. (1982). Estimation of earthquake source parameters by the inversion of waveform data: synthetic waveforms. Physics of the Earth and Planetary Interiors, 30(2–3), 242–259. https://doi.org/10.1016/0031-9201(82)90111-x
[7]Timoshin Yu. V. (1972). Fundamentals of diffraction conversion of seismic recordings. Moscow: Nedra. [In Russian].
[8]Trojanowski, J., & Eisner, L. (2016). Comparison of migration‐based location and detection methods for microseismic events. Geophysical Prospecting, 65(1), 47–63. https://doi.org/10.1111/1365-2478.12366
[9]Zhebel, O., & Eisner, L. (2014). Simultaneous microseismic event localization and source mechanism determination. Geophysics, 80(1), KS1–KS9. https://doi.org/10.1190/geo2014-0055.1