The present disclosure relates generally to analyzing microseismic events, and more particularly to techniques for estimating features of a wavefront associated with a microseismic event.
Passive localization of microseismic events is commonly used to monitor resource extraction processes such as hydraulic fracture stimulation, or “fracking.” In a typical scenario, an array of geophones is positioned in a monitor well nearby the well undergoing treatment. The array records seismic energy released impulsively from induced failure events as fractures form. Arrival times and polarizations of P-waves and S-waves impinging the array, among other features, are used to estimate the location of each detected event. Additional information about the event, treatment, and seismic propagation environment is contained in the arriving waveforms. For example, the measure of seismic propagation loss, Q, may be estimated by analyzing the spectral content of arrivals as the wavefront propagates across the sensors of the array. Additionally, the event moment magnitude and its moment tensor may be estimated from arriving waveform features.
Co-pending U.S. patent application Ser. Nos. 13/598,580, 14/340,356 and 14/536,985, the contents of which are all incorporated by reference herein in their entirety, have dramatically advanced the state of the art by providing various techniques for analyzing and characterizing microseismic events utilizing sensor data such as data from an array of geophones as described above. However, certain challenges remain.
Embodiments of the present invention provide a system and method for extracting event arrival waveforms and estimating sensor wavefront features such as arrival times and polarizations, which can operate in the presence of significant seismic background energy and which do not result in any systematic distortion of the arrival waveforms. According to certain aspects, an approach taken in the present invention is to recast the problem of suppressing the seismic background and sensor noise as one of estimating the arriving waveforms in a “noise” background. An aspect of this approach is to simultaneously estimate at each sensor pod wavefront features such as arrival times and polarizations, along with the arriving waveform itself.
In an embodiment of the present invention, recorded traces are time-aligned according to an estimate of the waveform moveout, and projected onto a low-dimensional space to estimate the arriving waveforms. In another embodiment, the estimated waveforms are correlated with the trace data to re-estimate wavefront arrival times, which in turn, are used in an iterative process to update the waveform estimates.
In other embodiments, the component traces for each sensor pod are combined according to an estimated event direction, forming a single oriented trace for each pod. In yet another embodiment, the component traces for each sensor pod are combined according to a direction at which the SNR is approximately maximized. In still another embodiment, the estimated waveforms are expanded to multi-component data according to the orientations used to generate the single-trace data used in estimating the waveforms and the event direction.
In further embodiments, the low-dimensional projection is an average of the time-aligned waveforms, scaled according to the inner product of the average and the portion of the trace data thought to contain the arrival. In another embodiment the average is weighted according to the SNR or energy of the time-aligned waveforms.
Yet another embodiment estimates wavefront arrival times and polarizations by correlating for each sensor pod a single estimated waveform with the pod trace data.
These and other aspects and features of the present embodiments will become apparent to those ordinarily skilled in the art upon review of the following description of specific embodiments in conjunction with the accompanying figures, wherein:
The present embodiments will now be described in detail with reference to the drawings, which are provided as illustrative examples of the embodiments so as to enable those skilled in the art to practice the embodiments and alternatives apparent to those skilled in the art. Notably, the figures and examples below are not meant to limit the scope of the present embodiments to a single embodiment, but other embodiments are possible by way of interchange of some or all of the described or illustrated elements. Moreover, where certain elements of the present embodiments can be partially or fully implemented using known components, only those portions of such known components that are necessary for an understanding of the present embodiments will be described, and detailed descriptions of other portions of such known components will be omitted so as not to obscure the present embodiments. Embodiments described as being implemented in software should not be limited thereto, but can include embodiments implemented in hardware, or combinations of software and hardware, and vice-versa, as will be apparent to those skilled in the art, unless otherwise specified herein. In the present specification, an embodiment showing a singular component should not be considered limiting; rather, the present disclosure is intended to encompass other embodiments including a plurality of the same component, and vice-versa, unless explicitly stated otherwise herein. Moreover, applicants do not intend for any term in the specification or claims to be ascribed an uncommon or special meaning unless explicitly set forth as such. Further, the present embodiments encompass present and future known equivalents to the known components referred to herein by way of illustration.
According to certain general aspects, embodiments of the invention provide techniques for extracting event arrival waveforms and estimating sensor wavefront features such as arrival times and polarizations, which techniques can operate in the presence of significant seismic background energy and which do not result in any systematic distortion of the arrival waveforms.
Referring to
According to certain aspects, the present inventors recognize that the wavefronts arriving at the sensors 104 from source 102 are often corrupted by significant seismic background energy and sensor noise, which obscures important event wavefront features. Processing of the sensor traces such as bandpass filtering or Wiener filtering is commonly performed to suppress frequency bands in which the event wavefront signal-to-noise ratio is poor (note that here “noise” refers to both the seismic background and sensor noise); however, in many cases such processing does not remove a sufficient amount of the seismic background to effectively understand details of the arriving wavefront. Other non-linear processing techniques, such as adaptive filtering methods, suppress noise effectively, but can alter the arriving waveforms in doing so.
According to certain aspects, therefore, embodiments provide a seismic trace processing method which can remove a significant amount of seismic background during an event wavefront arrival, while preserving important features of the arriving waveforms.
As in
gm(t)=sm(t−τm)βmT+nm(t) (1)
where •T denotes a matrix transpose. The goal is to estimate the arriving event waveforms sm(t) at each sensor m from the measured seismic traces gm(t); m=1 to M. It should be noted that for P-wave arrivals, sm(t) and βm will have one dimension or column, while for S-wave arrivals, sm(t) and βm can have two dimensions or columns. For clarity in the following, sm(t) and βm are treated as having one dimension or column, representing P-waves or S-waves polarized along a single direction. It will be apparent to those skilled in the art how to extend the invention described herein to the case of sm(t) and βm having two dimensions or columns.
An example set of raw measured traces is shown in
A functional block diagram of an example method and system for implementing aspects of embodiments of the invention is shown in
Returning to
Additionally in block 402, if the sensors contain multiple components, the respective traces from these components are combined into a single trace for each sensor. In one example, using an estimate for each respective polarization βm, the data for each sensor is “pointed” in the source direction, so as to combine the multiple traces for each sensor into a single component and single trace. The initial estimate for βm (i.e. estimate of source location relative to the position of the sensor in the array) can be determined in various ways, such as those described in the co-pending applications.
In another possible example, a single trace for each sensor is generated by pointing each sensor in the direction at which it achieves its maximum SNR. It should be noted that SNR may be estimated by examining the signal level before and during the event arrival, and that the direction of maximum SNR may be found by a search over arrival directions. It should be further noted that the direction at which the maximum SNR is achieved is not necessarily determined by the source direction. This is because while the source energy is maximized when the sensor is oriented according to the source direction, the maximum is rather broad. By contrast, energy from a noise source is minimized when the sensor is oriented perpendicular to the noise source polarization, but the minimum is relatively narrow. As a result, the SNR maximizing orientation will often be more aligned relative to a dominant noise source than according to the source direction.
Finally, in block 402, the set of time-aligned, oriented single traces is time-windowed to select the duration of trace data over which the wavefront appears at each sensor. An example of how block 402 processes the traces of
Denote by hm(t) the time-aligned, oriented, and windowed sensor signals as produced by block 402, and evaluated at sample time t,
hm(t)=w(t)·gm(t+τm)βm; m=1,2, . . . M (2)
where w(t) represents the time window, which is zero outside the T-sample long duration of the aligned arrivals. Denote by Hm the column of T samples of hm(t),
Hm=[hm(1)hm(2) . . . hm(T)]T (3)
and by H the T×M matrix of time-aligned, oriented, and windowed sensor signals,
H=[H1H2 . . . HM] (4)
The columns of the arrival waveform matrix H consist of sampled event waveforms sm(t) in a background of noise nm(t+τm).
The event waveforms are expected to be similar sensor to sensor, and should be accurately represented by the weighted sum of a relatively small number of basis shapes. The noise, on the other hand, is expected to vary in an uncorrelated manner from sensor to sensor. In this way, embodiments decompose the arrival waveform matrix into signal and noise subspaces (for example, via singular value decomposition), and estimate the event waveforms at each sensor as that portion of the arrival waveform matrix H lying in the subspace associated with its few larger singular values.
Denoting by S the matrix of arrival waveforms,
S=[s1s2 . . . sM] (5)
where sm is the vector of sensor m event waveform samples sm(t); t=1, 2, . . . T, embodiments model the arrival waveforms as
Ŝ=ULWL (6)
where UL is an L-column matrix of T-long orthonormal waveform basis shapes, and WL is an M-column matrix of L-element weights. The basis waveform shapes and weights (i.e. the matrices U and V) can be found via singular value decomposition of H as performed in block 404. In one non-limiting example, singular value decomposition is performed using the MATLAB function call svd(H, ‘econ’). The outputs of this function call using the matrix H formed using the outputs of block 402 can be expressed as:
UDVT=H (7)
where U is the T×M matrix of orthonormal basis shapes, D is an M×M diagonal matrix of singular values dl; l=1, 2, . . . M, arranged largest to smallest, and V is an M×M orthonormal matrix of normalized weights.
In forming the model waveform basis shapes U and sensor weights W, it is desired to take into consideration the relative sensor SNRs, in particular, it is desired to give more weight to those sensors having larger SNRs. A weighted decomposition may be formed as follows. Denote by Σ the sensor noise covariance matrix, e.g. a diagonal matrix of sensor noise variances. The matrix HΣ−1 then weighs the columns of H by the respective sensor inverse variances. Performing a singular value decomposition on HΣ−1 gives (after multiplying through by Σ),
UΣDΣVΣTΣ=H,
with the orthonormal UΣ being the waveform basis shapes, and DΣVΣTΣ being the sensor weights.
Using these matrices, in block 406, to estimate the sensor waveforms during the arrival, the order L is selected, and Equation (6) is used to construct the estimated waveforms, with UL being the first L columns of U, and WL being the first L rows of the product DV. (Note that the singular values D could be apportioned in any manner between U and VT in computing the waveform estimates.) In some sense, Ŝ reconstructs the event arrival waveforms without using the noise-like components.
As an example,
Returning to the example embodiment shown in
First, in block 408, the first L columns of matrix V are obtained from the output of block 404.
Next, in block 408, embodiments generate time-aligned, oriented traces (essentially hm(t), but not windowed),
km(t)=gm(t+τm)βm; m=1,2, . . . M (8)
and form the matrix K using km(t) as its t, mth row-column element. In this way, the entire duration of the traces is processed, not just the time window used to find V.
Then, in block 408, embodiments project the traces onto the first L columns of V, denoted VL to form the estimated traces Ř, using the following equation:
Ř=KVLVLT (9)
Note that the estimated traces Ř will match the estimated arrival waveforms Ŝ in the arrival time window.
Finally in block 410, the estimated waveforms Ř are time shifted back to their original moveout using τm, and each sensor signal is rotated according to its polarization using βm to generate multicomponent “cleaned” traces. As an example, processed versions of the raw trace data appearing in
Trace estimates made using L=5 are shown in
One question concerns the nature of the waveform estimation errors in Ŝ or Ř. An expectation is that the estimation errors will be additive and noise-like as the estimates are constructed additively, and their errors are expected to result from including unwanted portions of additive noise in the signal basis or including noise subspace components by choosing an order L which is too large. Therefore, this process is not expected introduce systematic changes in the phase or magnitude of the arriving wavefronts. This seems to be the case, as evidenced by the similarity between the “cleaned” and bandpass filtered arrivals as the wavefront traverses the array in.
According to certain additional aspects, the estimated arrival waveforms can be used to generate improved event arrival-time and polarization estimates for each sensor in embodiments of the invention.
First, in block 1102, the estimated arrival waveforms Ŝ are generated as described above (e.g. in connection with blocks 402, 404 and 406).
Next, in block 1104, event arrival polarizations may be estimated using the projection of the estimated arrival waveforms S onto the corresponding sensor's multi-component trace data:
where ŝm is the estimated arrival waveform at the mth sensor.
Similarly, in block 1104, arrival times can be estimated by cross-correlating the estimated arrival waveforms with the corresponding sensor's trace data, preferably oriented according to the arrival polarization or the direction of maximum SNR, and finding the correlation lag at which the cross correlation is maximized. Note that these cross-correlations could also be performed using multi-component input trace data and the multi-component signal estimates made using the method described below directly, without reducing the data to a 1-dimensional signal per sensor.
The case of planar motion (two-dimensional sm(t) and m) may be accommodated by first solving for the first column of βm1, then subtracting that direction ŝβm1T from gm, and finally repeating the 1-D process with that first dimension of motion removed.
Improved event arrival waveform estimates may be found iteratively by first using the process above to estimate the arriving waveforms in block 1102, and then using those waveforms to form new arrival-time and polarization estimates in block 1104. As shown in
According to still further aspects, another embodiment of the invention forms the arrival waveform matrix without first combining multi-component sensor signals to produce a reduced-dimension signal for each sensor such as a single, directionally oriented trace for each sensor. Denote by ĥm(t) the time-aligned and windowed component signals of sensor m, evaluated at sample time t,
ĥm(t)=w(t)·gm(t+τm); m=1,2, . . . M (12)
where w(t) represents the time window, which is zero outside the T-sample long duration of the aligned arrivals. For a collection of sensors each having N components, denote by Ĥ the T×NM matrix of time-aligned and windowed multi-component sensor signals,
where the multi-component arrival waveform matrix Ĥ is N-times as wide as the oriented arrival waveform matrix H.
To estimate the event arrival waveforms, the multi-axis waveform matrix is split into signal and subspace matrices. As above, a singular value decomposition can be used,
{circumflex over ({dot over (U)})}Ď{tilde over (V)}T={hacek over (H)} (14)
with {circumflex over ({dot over (U)})}, Ď and {tilde over (V)} being the orthonormal basis shapes, diagonal matrix of singular values (arranged largest to smallest), and orthonormal matrix of normalized weights. As above, the multi-component arriving wavefront signal may be estimated using a given number of singular values L,
Ŝ={circumflex over ({dot over (U)})}L{acute over (W)}L (15)
where {circumflex over ({dot over (U)})}L is the matrix made from the first L columns of {circumflex over ({dot over (U)})}, and {acute over (W)}L is the first L rows of Ď{circumflex over (V)}T. Similarly, any duration of trace data may be projected onto the columns of {circumflex over (V)}L (as above) to estimate the arrival signal from an entire record.
Note that the N columns of {acute over (W)}L associated with a given sensor are likely to have a rank which is the lesser of N and L. Accordingly, when L>1, the different waveform shapes specified by {circumflex over ({dot over (U)})}L will be added in different proportions to estimate the arrival waveform along each component of a sensor. However, a P-wave arrival has motion along only the arrival direction; an S-wave arrival can involve motion only in the plane perpendicular to the arrival direction. As a result, it is preferred to replace the columns of {acute over (W)}L associated with each sensor by either a rank-one or rank-two approximation, depending on the arrival type, for instance found by singular value decomposition of each sensor's L×N portion of {acute over (W)}. Note that the row space of this rank-one or rank-two approximation can be used as an estimate of the arrival polarization. Finally, note that this procedure may be applied in the case that different sensors within the array have varying number of components.
Although the present embodiments have been particularly described with reference to preferred ones thereof, it should be readily apparent to those of ordinary skill in the art that changes and modifications in the form and details may be made without departing from the spirit and scope of the present disclosure. It is intended that the appended claims encompass such changes and modifications.
The present application claims priority to U.S. Prov. Appln. No. 62/211,893, filed Aug. 31, 2015, the contents of which are incorporated herein by reference in their entirety.
Number | Name | Date | Kind |
---|---|---|---|
3453614 | Bose | Jul 1969 | A |
3460013 | Gaylor | Aug 1969 | A |
4683556 | Willis | Jul 1987 | A |
5583825 | Carrazzone et al. | Dec 1996 | A |
5996726 | Sorrells | Dec 1999 | A |
7391675 | Drew | Jun 2008 | B2 |
7643374 | Plona et al. | Jan 2010 | B2 |
8107317 | Underhill et al. | Jan 2012 | B2 |
8649980 | Gulati | Feb 2014 | B2 |
8705316 | Thornton et al. | Apr 2014 | B2 |
8800652 | Bartko et al. | Aug 2014 | B2 |
20050060099 | Sorrells et al. | Mar 2005 | A1 |
20080159075 | Underhill et al. | Jul 2008 | A1 |
20090010104 | Leaney | Jan 2009 | A1 |
20090125240 | Den Boer et al. | May 2009 | A1 |
20100262373 | Khadhraoui et al. | Oct 2010 | A1 |
20100312480 | Hansteen et al. | Dec 2010 | A1 |
20110091078 | Kherroubi et al. | Apr 2011 | A1 |
20110120702 | Craig | May 2011 | A1 |
20110317518 | Guigne et al. | Dec 2011 | A1 |
20130144532 | Williams et al. | Jun 2013 | A1 |
20130215717 | Hofland et al. | Aug 2013 | A1 |
20140019057 | Diller | Jan 2014 | A1 |
20140112099 | Hofland et al. | Apr 2014 | A1 |
20140278120 | Kahn et al. | Sep 2014 | A1 |
20140288840 | Vermilye | Sep 2014 | A1 |
20150226868 | Cieplicki | Aug 2015 | A1 |
Number | Date | Country |
---|---|---|
1 485 733 | Dec 2011 | EP |
2 784 552 | Oct 2014 | EP |
2 450 707 | Jan 2009 | GB |
2 525 996 | Nov 2015 | GB |
WO-2004070424 | Aug 2004 | WO |
WO-2005106533 | Nov 2005 | WO |
WO-2008124759 | Oct 2008 | WO |
WO-2009026979 | Mar 2009 | WO |
WO-2010116236 | Oct 2010 | WO |
WO-2011077223 | Jun 2011 | WO |
WO-2011077227 | Jun 2011 | WO |
WO-2012085848 | Jun 2012 | WO |
Entry |
---|
Coffin, et al., “On the Accuracy of Microseismic Event Location Estimates”, undated white paper, pp. 1-18. |
Bahorich et al., “3-D Seismic Discontinuity for Faults and Stratigraphic Features: The Coherence Cube”, The Leading Edge, Oct. 1995, pp. 1053-1058. |
Chopra et al., “Curvature Attributes Delineate Subtleties”, E & P, Sep. 2006, pp. 72-73. |
Culver et al., “Interpretation of Multipath-Analysis Localizations of Microseismic Data: An Alberta Montney Shale Example”, undated white paper, 5 pages. |
Eisner, et al., “Comparison of Surface and Borehole Locations of Induced Seismicity”, Geophysical Prospecting, 2010, 58, 809-820. |
Eisner, et al., “Determination of S-wave Slowness from a Linear Array of Borehole Receivers”, Geophys. J. Int. (2009) 176, 31-39. |
Eisner, et al., “Uncertainties in Passive Seismic Monitoring”, The Leading Edge, Jun. 2009, pp. 648-655. |
Fuller et al., “Seismic Wave Phenomena and Implications for Accuracy of Microseismic Results”, CSEG Recorder, Nov. 2010, pp. 36-37. |
Grenchka, et al., “Data-Acquisition Design for Microseismic Monitoring”, The Leading Edge, Mar. 2010, pp. 1-5. |
Hur, et al., “An Analytic Model for Microseismic Event Location Estimate Accuracy”, Oct. 2011, 25 pages. |
Hur et al., “An Analytic Model for Microseismic Event Location Estimate Accuracy”, First Break, undated, 7 pages. |
International Search Report and Written Opinion dated Sep. 25, 2013 for International Patent Application No. PCT/US2013/040191, 11 pages. |
Kidney, et al., “Impact of Distance Dependent Location Dispersion Error on Interpretation of Microseismic Event Distributions”, The Leading Edge, Mar. 2010, pp. 1-5. |
Maxwell, et al., “Key Criteria for a Successful Microseismic Project”, SPE Paper 134695, Sep. 2010, pp. 1-16. |
Maxwell, et al., “Microseismic Location Uncertainty”, CSEG Recorder, Apr. 2009, pp. 41-46. |
International Search Report dated May 27, 2013 for International Patent Application No. PCT/US2013/022950, 1 page. |
Rutledge, et al., “Hydraulic Stimulation of Natural Fractures as Revealed by Induced Microearthquakes”, Carthage Cotton Valley Gas Field, East Texas, Geophysics, vol. 68, No. 2, 2003, pp. 1-37. |
Roth et al., “Fracture Interpretation in the Barnett Shale, Using Macro and Microseismic Data”, Frontiers and Innovation—CSPG CSEG SWLS Convention, 2009, pp. 1-4. |
Zimmer et al., “Localization of Microseismic Events Using Headwaves and Direct Waves”, SEG Expanded Abstracts 2010, pp. 2196-2200. |
Number | Date | Country | |
---|---|---|---|
62211893 | Aug 2015 | US |