Systems and methods for non-parametric autopicking of seismic wave features from seismic data

Abstract

Systems and methods for automatically detecting and identifying seismic wave features in seismic data are provided. In general, the systems and methods utilize a nonparametric time series classification method to detect seismic wave features that may otherwise be difficult to automatically identify in seismic data. Instead of building a model by estimating parameters from the seismic data, the data itself is used to define a trained model. The systems and methods described here provide the ability to detect and identify seismic wave features with reasonably fast and extremely accurate results without needing to compute parameters.

Claims

1. A computer system for detecting and identifying a seismic wave feature in seismic data using a non-parametric time series classification model, the computer system comprising: a memory having instructions stored therein; a processor configured to access the instructions stored in the memory and execute the instructions to cause the computer system to: (a) store a set of training seismic data comprising a positive reference containing an S wave and a negative reference in the memory; (b) generate a detection window having a detection window size corresponding to a temporal characteristic of the S-wave and store the detection window in the memory; (c) provide seismic data based on a seismic event that includes seismic traces; (d) scale the seismic data based on information contained in the seismic data; (e) perform a sliding window analysis by moving the detection window through a plurality of different time points in the seismic data and computing a similarity score for each of the plurality of different time points, wherein computing the similarity score includes: computing a positive reference similarity score between the seismic data and the positive reference; computing a negative reference similarity score between the seismic data and the negative reference; and computing the similarity score based on a ratio between the positive reference similarity score and the negative reference similarity score; (f) select at least one of the similarity scores from the plurality of different points in the seismic data as corresponding to the S-wave; and (g) generate a report that includes a location or identification of the seismic event based on the S-wave.

2. The computer system as recited in claim 1, wherein step (d) includes computing an absolute value of an amplitude of the seismic data and scaling the seismic data based on the computed absolute value.

3. The computer system as recited in claim 2, wherein step (d) includes computing an average amplitude over a portion of the seismic data and scaling the seismic data by dividing the computed absolute value by the computed average.

4. The computer system as recited in claim 1, wherein the negative reference is defined as a portion of the training seismic data that does not include an S-wave.

5. The computer system as recited in claim 1, wherein step (e) includes computing a distance metric between the training seismic data and the seismic data.

6. The computer system as recited in claim 5, wherein the distance metric is a squared Euclidean distance.

7. The computer system as recited in claim 1, wherein selecting the portion of the provided seismic data as containing the S-wave in step (f) includes determining the time point at which the detection window is associated with a maximum similarity score.

8. The computer system as recited in claim 1, wherein computing the similarity score includes computing a vote for each detection window, wherein the vote is based in part on a distance metric between the training seismic data and the seismic data.

9. The computer system as recited in claim 8, wherein the vote includes weighting the distance metric by a weighting value.

10. The computer system as recited in claim 1, further comprising filtering the provided seismic data using a band-pass filter to remove selected frequency components from the seismic data.

11. The computer system as recited in claim 10, further comprising filtering the provided training seismic data using the band-pass filter to also remove the selected frequency components from the training seismic data.

12. The computer system as recited in claim 1, wherein the provided training seismic data includes a numerical model that represents the S-wave.

13. The computer system as recited in claim 1, wherein the positive reference in the provided set of training seismic data further includes at least one identified P-wave, P wave to S wave conversion, or surface wave, and wherein step (d) further includes processing the seismic data and the training seismic data to identify a portion of the seismic data as containing a P-wave, P wave to S wave conversion, or surface wave based on a computed similarity score between the seismic data and the training seismic data.

14. The computer system as recited in claim 1, further comprising determining a quality of the S-wave based on a Wadati diagram.

15. The computer system as recited in claim 7, wherein the S-wave is selected as corresponding to a center position of the detection window having the maximum similarity score.

16. The computer system as recited in claim 1, wherein the report that includes the location or identification of the seismic event based on the S-wave is generated in real-time.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flowchart setting forth the steps of an example method for detecting a seismic wave feature of interest in seismic data using appropriately selected training data;

(2) FIG. 2 is an example of a Wadati diagram and L.sub.1 fit for an example training data set, which can be used to improve the accuracy of the training data by identifying outliers to be removed from the training data;

(3) FIG. 3 is an example of positive and negative reference data sets, which include four horizontal component records containing S-wave analyst picks and four horizontal component records containing noise, and in which each window is scaled using the absolute amplitude of the data in the window and divided by the mean of the second half of the window;

(4) FIG. 4 is an example of two horizontal traces (solid lines) and their associated R.sub.pair score curve (dotted line);

(5) FIG. 5 is an example illustrating a re-scaling scheme used in some embodiments of the present invention, whereby an un-scaled example of an S-wave arrival at 0.33 seconds (left) is rescaled by computing the absolute value of the data's amplitude (middle) and then dividing that absolute amplitude by the mean of the second half of the S-wave portion of the signal (0.34 seconds to 0.67 seconds) (right); and

(6) FIG. 6 is a block diagram of an example computer system that can implement some embodiments of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

(7) Described here are systems and methods for automatically detecting and identifying seismic wave features in seismic data. Seismic wave features that can be detected and identified include P-waves, S-waves, P-wave to S-wave conversions, surface waves, and so on.

(8) In general, the systems and methods utilize a nonparametric time series classification method to detect seismic wave features that may otherwise be difficult to automatically identify in seismic data. Instead of building a model by estimating parameters from the seismic data, the data are used to define a model. This approach runs contrary to previous methods for detecting seismic wave features, which relied on computing parameters from the seismic data and then basing the feature detection on those parameters. As an example, some machine learning methods, such as neural networks and neural trees, have been used for phase identification, but the majority of these previous methods are trained on parameters, such as AR-model coefficients or short-term average/long-term average (STA/LTA) values, instead of on the raw data. These previous methods have limited reliability because detection is based on parameters, whose computation could be affected by error sources, rather than the raw data. The systems and methods of the present invention, however, provide the ability to detect and identify seismic wave features with reasonably fast and extremely accurate results without needing to compute parameters.

(9) In some embodiments, to address the limitations of previous seismic wave feature detection techniques, the systems and methods of the present invention are based on computing the similarity of input signals to a reference, or training, data set containing the seismic wave features of interest. Potential seismic wave features of interest, which may include S-wave arrivals, are identified based on their similarity to a set of training data that includes positive and negative reference sets. As an example, the positive set contains examples of analyst-identified seismic wave features and the negative set contains noise, which may include P-wave or S-wave coda.

(10) Referring now to FIG. 1, a flowchart is illustrated as setting forth the steps of an example method for automatically detecting and identifying a seismic wave feature in seismic data. A seismic wave feature may include a P-wave, and S-wave, a P-wave to S-wave conversion, a surface wave, and so on.

(11) The method includes providing training data, as indicated at step 102. In some embodiments, the training data may include seismic data containing at least one seismic wave feature that is similar to the seismic wave feature of interest to be detected and identified in the acquired data. For instance, if the detection task is to identify S-waves in seismic data, then the training data will include at least one S-wave.

(12) In general, the training data includes a positive reference and a negative reference. The positive reference can be a portion of seismic data that has been identified as corresponding to the seismic wave feature of interest. For instance, the positive reference can include a portion of seismic data manually identified as containing an S-wave. For example, the positive reference can include horizontal traces identified as containing the seismic wave feature of interest. The negative reference can be a portion of seismic data that does not contain the seismic wave feature of interest. As one example, the negative reference can be selected as a noisy portion of the seismic data occurring after the positive reference. For instance, the negative reference can be selected as S-wave coda windows extracted from the positive reference. As another example, the negative reference can be selected as a noisy portion of the seismic data occurring before the positive reference.

(13) The training data does not need to be geographically or geologically specific. That is, training data based on seismic data acquired from one geographic or geologic region can be used to detect seismic wave features in seismic data acquired from a completely different geographic or geologic region.

(14) In some embodiments, the training data can include a simple model of the desired seismic wave feature. As one example, the training data can comprise a step function that models the initiation of an S-wave.

(15) As one specific, and non-limiting, example, the training data can be established as follows. For this example, two different datasets containing S-wave picks from two different analysts were used as the basis for two different training data sets. The first data was acquired from the Parkfield region in central California, and the second was acquired from the Alpine Fault region on the South Island of New Zealand. The Alpine Fault dataset includes 209 analyst picks from over 30 different broadband and short period stations along approximately 50 km of the Alpine Fault from the WIZARD array operated in 2012-2014. The Parkfield dataset includes 721 analyst picks from over 50 stations of the PASO array along the San Andreas Fault operated in 2001-2002.

(16) These training data sets were then processed, as will be described in more detail below. In this example, the horizontal traces for each station-event pair were band-pass filtered between 3 Hz and 30 Hz and cut 0.5 seconds after the known P-wave arrival. In addition, all signals were down-sampled to the minimum sampling rate of all signals. Next, an expected S-wave arrival was calculated assuming a P-wave speed to S-wave speed (.sub.P/.sub.S) ratio of 1.7, and by defining an S-wave arrival region as 2 seconds around the expected S-wave arrival. Signals for which the analyst S-wave pick fell outside this 4 second region were removed from the analysis, considering them to be inconsistent and therefore unreliable.

(17) In addition, erroneous analyst picks can be identified and removed using a Wadati diagram, excluding outliers from a least absolute deviation (L.sub.1) fit. An example of such a Wadati diagram is illustrated in FIG. 2, which shows a Wadati diagram and L.sub.1 fit for Parkfield event 2002009235200. Analyst pick outliers (solid black dots) are identified as erroneous picks and removed from the analysis if points fall outside of lines with slope +/0.2 that of the L.sub.1 fit. Using these quality control methods, analyst picks can be reasonably removed from the provided training data sets to improve their accuracy.

(18) Referring again to FIG. 1, after the training data have been provided, seismic data are acquired or otherwise provided, as indicated at step 104. In some embodiments, the seismic data can be acquired and processed in real-time.

(19) A detection window size is then selected, as indicated at step 106. This detection window size defines a detection window that will be used for automatically detecting and identifying the seismic wave feature of interest in provided seismic data. By way of example, the detection window size can be selected based on known temporal characteristic of the seismic wave feature of interest. For instance, if the seismic wave feature of interest is an S-wave arrival, then the window size can be selected based on the known temporal characteristics of S-wave arrivals. As one example, the window length can be selected as 0.67 seconds. It will be appreciated that other window sizes can be suitably selected depending at least in part on the seismic wave feature of interest.

(20) The seismic data are then filtered, as indicated at step 108. As an example, the seismic data are filtered using a bandpass filter that is designed to remove parts of the seismic data that are not important for the detection task. When the training data includes seismic data, the seismic data are preferably also filtered, as indicated at step 110, and using the same filter used on the acquired seismic data, as described above. By using the same filter, the frequency content of the seismic data and the training data will be equalized, which can improve the detection efficiency of the seismic wave feature picker.

(21) As another pre-processing step, the seismic data are scaled, as indicated at step 112. In general, the seismic data can be scaled using a scaling factor, which in some embodiments is determined from the seismic data. As one example, the scaling factor can be determined as the scaling factor that normalizes a selected portion of the seismic data. For instance, a portion of the seismic data that is thought to contain the seismic wave feature of interest can be selected and normalized, and the scaling factor that normalizes that portion of the seismic data can be used to scale the seismic data as a whole.

(22) Using the training data and the selected detection window, the scaled seismic data are processed to detect and identify seismic wave features, as indicated at step 114. In general, the seismic wave feature of interest can be detected and identified using a sliding window analysis, whereby seismic data within the selected detection window is compared with the training data.

(23) In some embodiments, a k-nearest neighbors based approach is used to detect and identify seismic wave features with the appropriately prepared training data. As one example, similarity between the seismic data and the training data can be computed. In some embodiments, similarity can be defined as the Euclidean distance between the scaled absolute values of the amplitudes of the observed signal and a given reference example. It will be appreciated by those skilled in the art, however, that other distance or similarity metrics can also be used.

(24) By way of example, a seismic wave feature of interest can be detected and identified as follows. A single pass is done through the band-passed, scaled seismic data to produce a seismic wave feature score function that is defined, for each window, as the ratio of the sum of similarity to positive reference examples over the sum of similarity to negative reference examples. The seismic wave feature of interest is chosen as the center position of the window that maximizes this score function.

(25) Using an S-wave arrival as the seismic wave feature of interest, S-wave onsets can be identified in seismic data by taking successive windows of a given horizontal trace and computing a similarity metric to the set of positive and negative reference set examples. For both horizontal traces in a given station-event combination, an S-wave score function is computed for the two traces at each time step, and these scores are multiplied together to produce a final score. The S-wave arrival is chosen at the center of the window corresponding to the maximum score.

(26) In this example, each horizontal trace is processed independently and the following analysis is applied to each un-rotated horizontal trace. A moving-window procedure is performed, using the previously selected window size. For each time window, a vote, V, is calculated by comparing the observed data, s, to the reference trace, r, in the training data:
V(r,s)=e.sup.d(r,s) (1);

(27) where is a weighting parameter and d(r,s) is a distance metric between the reference signal and the observed data. As one example, the distance metric can be the square of the Euclidean distance,

(28) $\begin{array}{cc}d\left(r,s\right)=\underset{i=1}{.Math.}{\left(r_i-s_i\right)}^2.& \left(2\right)\end{array}$

(29) Each vote represents the degree of similarity between the observed data and a given reference signal in the training data. As described above, the training data, R, contains equal numbers of positive and negative reference signals. An example of a positive reference set and a negative reference set is illustrated in FIG. 3. The votes are aggregated for the positive and negative comparisons, and a score, R(s), is calculated for each window by taking the ratio of the positive to negative vote sums:

(30) $\begin{array}{cc}R\left(s\right)=\frac{\underset{r{R}_{+}}{.Math.}V\left(r,s\right)}{\underset{r{R}_{-}}{.Math.}V\left(r,s\right)}.& \left(3\right)\end{array}$

(31) FIG. 4 shows an example of two horizontal traces and their associated R-curves. For a given station-event pair, two scores R.sub.H1 and R.sub.H2 are computed for the two horizontal traces, and a single score, R.sub.pair, is computed by combining the two horizontal R-score values, which can be accomplished using multiplication or addition as follows:
R.sub.pair(s)=R.sub.H1(s).Math.R.sub.H2(s) (4);

(32) or,
R.sub.pair(s)=R.sub.H1(s)+R.sub.H2(s) (5).

(33) The S-wave window is then determined as the window, s.sub.max, that maximizes R.sub.pair,

(34) $\begin{array}{cc}{s}_{\max }=\underset{s}{\max }\left(R_{\mathrm{pair}}\left(s\right)\right).& \left(6\right)\end{array}$

(35) Finally, the S-wave pick is selected as the time of the center of the s.sub.max window.

(36) As one example implementation of the methods described above, the weighting parameter, , can be selected as 0.0001, the moving-window length can be selected as 0.67 seconds, and the data can be band-pass filtered between 3 Hz and 30 Hz. Positive reference waveforms were created using the portion of the training data 0.33 seconds before and after the analyst S-wave pick, whereas negative reference waveforms were created by cutting a 0.67 second window several seconds after the S-wave analyst pick. Both the positive and negative reference waveforms were scaled by taking the absolute value of the amplitude and dividing by the mean amplitude of the second half of the 0.67 second window of the positive examples. An example of seismic data scaled in this manner is illustrated in FIG. 5.

(37) A quality rating that includes two classes can be used to rate the quality of a seismic wave feature pick. For example, a class 1 can be defined based on a signal-to-noise ratio (SNR) greater than 1, and a class 0 can be defined based on an SNR greater than 10. SNR can be calculated by selecting noise as the RMS value from the first half of the window before the S-wave onset and selecting the signal as the RMS value in the second half of the window after the S-wave onset.

(38) In some embodiments, the detected and identified seismic wave features can also be processed to remove erroneous picks, similar to those methods that can be applied to improve the accuracy of the training data. For instance, a Wadati diagram and L.sub.1 fit, such as those illustrated in FIG. 2, can be used to this end. As one example, after automatic S-wave picks are identified, a Wadati diagram can be constructed for each event using automatic S-wave picks and analyst P-wave picks. Outliers in this Wadati diagram can then be marked as incorrect picks and are therefore not used.

(39) Thus, a simple, fast, and accurate method for automatically identifying seismic wave features, such as S-wave onsets, in seismic data has been provided. In some embodiments, a k-nearest neighbor based technique is applied to band-passed horizontal seismograms. The method identifies potential S-wave onsets by comparing the similarity of a given window to a set of known S-wave onsets and to a set of examples containing noise. An S-wave onset score function is calculated using the ratio of similarity to the S-wave examples divided by the similarity to the noise examples. The method only needs to use the horizontal component of seismograms, although other components can also be included in the analysis, and needs only one pass through the band-passed data with a fixed window size.

(40) The method has few parameters to tune, and the method has been found to be invariant to the training set. It is contemplated that the method requires only two reference picks to achieve an acceptable classification rate. The success of the method using training sets from different geographic locations suggests that, at least in some cases, the method can work as an out-of-the-box auto-picker using a provided training set.

(41) Referring now to FIG. 6, a block diagram is illustrated of an example computer system 600 that can be used for detecting and identifying seismic wave features in seismic data in accordance with some embodiments of the present invention. The system 600 generally may include an input 602, at least one processor 604, a memory 606, and an output 608. The system 600 may be, for example, a workstation, a notebook computer, a personal digital assistant (PDA), a multimedia device, a network server, a mainframe, or any other general-purpose or application-specific computing device. The computer system 600 may operate autonomously or semi-autonomously, or may read executable software instructions from a computer-readable medium (such as a hard drive, a CD-ROM, flash memory, and the like), or may receive instructions from a user, or any another source logically connected to a computer or device, such as another networked computer or server, via the input 602.

(42) The input 602 may take any shape or form, as desired, for operation of the computer system 600, including the ability for selecting, entering, or otherwise specifying parameters consistent with operating the computer system 600. In some instances, the input 602 may be designed to receive training data that represents a seismic wave feature of interest and seismic data.

(43) Among the processing tasks for operating the computer system 600, the at least one processor 604 may be configured to perform the method described above with respect to FIG. 1.

(44) The memory 606 may contain software 608, and may be configured for storage and retrieval of processed information and data to be processed by the processor 604. In some aspects, the software 608 may contain instructions directed to performing the method described above with respect to FIG. 1.

(45) The present invention has been described in terms of one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention.