METHOD AND SYSTEM FOR SEISMIC DATA COMPRESSION AND NOISE REDUCTION
20220381933 · 2022-12-01
Assignee
Inventors
Cpc classification
International classification
Abstract
In the field of seismic data collection and analysis, the problem of data compression is considered. It is desired to reduce the storage and transmission requirements of seismic data, such as associated with microseismic monitoring and processing, VSP (vertical seismic profile) surveys, and the like, for instance using a distributed acoustic sensor (DAS), which can generate in excess of 50 GB of data per hour on a survey that lasts days. A method and system for data compression that separates the data collected into additive signal and noise components, and compresses the estimated signal component for transmission, storage, and analysis, is described. The idea is that the signal component, which exhibits clear structure across the traces, may be accurately described with relatively few parameters, and therefore may be significantly compressed without loss of important information.
Claims
1. A method for processing a set of seismic traces from an array of sensors so as to reduce its storage or transmission requirements, while reducing the presence of seismic background or noise in a reconstructed set of seismic signals, comprising: identifying one or more component seismic signal sources appearing additively in the set, each component source defined by one or more component source time waveforms and related amplitudes; and storing or transmitting data representing the component source time waveforms and related sensor amplitudes.
2. A method for processing a set of seismic traces from an array of sensors so as to reduce its storage or transmission requirements comprising: identifying a component seismic signal of interest, the component seismic signal of interest comprising a time waveform and sensor amplitudes; and storing or transmitting data representing the component seismic signal of interesting time waveform and sensor amplitudes.
3. The method of claim 2, wherein the component seismic signal component of interest also comprises one or more sensor time offsets.
4. The method of claim 2, further comprising: forming basis vectors that describe a subspace of a matrix or tensor of data derived from the set of seismic traces.
5. A method for reducing the storage or transmission requirements of a set of seismic traces from an array of sensors comprising: forming basis vectors that describe a subspace of a matrix or tensor of data derived from the set of seismic traces.
Description
BRIEF DESCRIPTION OF THE DRAWINGS
[0016] 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:
[0017]
[0018]
[0019]
[0020]
[0021]
[0022]
[0023]
[0024]
[0025]
[0026]
[0027]
DETAILED DESCRIPTION
[0028] 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.
[0029] According to certain general aspects, an approach to data compression according to embodiments flows from the observation that the desired arriving wavefronts have structure which can be expressed with relatively few parameters, while the additive seismic background and sensor noise are difficult to compress, requiring a significant amount of data to describe. Here, we process the recorded seismic traces to produce an encoded signal component that has been separated from an additive noise component. In this way, the signal component retains its clear structure across the traces, and is compactly described without loss of important information.
[0030] One example system according to embodiments is shown in
[0031] It should be pointed out that while the following discussion is relative to transient seismic events, the teachings described herein apply just as well to other sources of seismic energy arriving at the sensor array in a way that is spatially structured, that is structured in some way across sensors.
[0032] Consider first a seismic record from an array of M sensors containing a single seismic event. Each sensor receives signals g.sub.m(t), measured along one or more sensor components as a function of time t. The measured signal at sensor m is assumed to be an event waveform s.sub.m(t−τ.sub.m) arriving at time τ.sub.m with polarization β.sub.m in additive noise n.sub.m(t) consisting of seismic background and sensor noise,
g.sub.m(t)=s.sub.m(t−τ.sub.m)β.sub.m.sup.τ+m.sub.m(t), (1)
where.Math..sup.τ denotes matrix transpose. The goal is to estimate the arriving event waveforms s.sub.m(t) at each sensor m from the measured seismic traces g.sub.m(t), m=1, . . . , M. It should be noted that for P-wave arrivals, s.sub.m(t) and β.sub.m will have one column, while for S-wave arrivals, s.sub.m(t) and β.sub.m can have two columns. For clarity in the following, s.sub.m(t) and β.sub.m are treated as having one 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 principles described herein to the case of s.sub.m(t) and β.sub.m having multiple (e.g. two or more) columns.
[0033] A set of measured traces 202 is shown in
[0034] In an embodiment, trace data (band-pass filtered, e.g.) is time aligned according to an estimate of the arrival moveout so as to have the event wavefront arrive at all sensors simultaneously. Additionally, if the sensors contain multiple components, these components are combined in such a way as to point each sensor in the direction at which it achieves its maximum SNR, thereby generating a single trace 402 for each sensor. (Other methods, such as pointing the sensor according to the source direction, may be used to combine the sensor traces into a single component.) This set of aligned, oriented traces is then time windowed to select the duration over which the wavefront appears. The traces of
[0035] 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 pointed out that the direction at which the maximum SNR is achieved is not necessarily determined by only 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 perpendicular to a dominant noise source than along the source direction.
[0036] Denote by h.sub.m(t) the time-aligned, oriented, and windowed sensor signals, evaluated at sample time t,
h.sub.m(t)=w(t).Math.g.sub.m(t+τ.sub.m)β.sub.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 h.sub.m the column of T samples of h.sub.m(t),
h.sub.m=[h.sub.m(1)h.sub.m(2) . . . h.sub.m(T)].sup.T, (3)
and by H the T×M matrix of time-aligned, oriented, and windowed sensor signals,
H=[h.sub.1 h.sub.2 . . . h.sub.M]. (4)
[0037] The columns of the arrival waveform matrix H consist of sampled event waveforms s.sub.m(t) in a background of noise n.sub.m(t+τ.sub.m).
[0038] 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, one can decompose the arrival waveform matrix into signal and noise subspaces (say, 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.
[0039] Denoting by S the matrix of arrival waveforms,
S=[s.sub.1 s.sub.2 . . . s.sub.M], (5)
where s.sub.m is the vector of sensor m event waveform samples s.sub.m(t), t=1, 2, . . . , T, one can model the arrival waveforms as
Ŝ=U.sub.LW.sub.L, (6)
where U.sub.L is an L-column matrix of T-long orthonormal waveform basis shapes, and W.sub.L is an M-column matrix of L-element weights. The basis waveform shapes and weights can be found via singular value decomposition of H, for instance using the MATLAB function call svd(H,′econ′), to write
UDV.sup.1=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. To estimate the sensor waveforms during the arrival, the order L is selected, and (6) is used to construct the estimated waveforms, with U.sub.L being the first L columns of U, and W.sub.L being the first L rows of the product DV.sup.1. (Note that the singular values D could be apportioned in any manner between U and V.sup.1 in computing the waveform estimates.) In some sense, Ŝ reconstructs the event arrival waveforms without using the noise-like components.
[0040] As an example,
[0041] While estimates of the arriving waveforms are useful in determining arrival parameters, it may also useful to generate “cleaned” versions of the original trace data that has been similarly processed to preserve event arrival waveforms while reducing background noise. Note that the singular value decomposition above produces a time basis U and a sensor basis V. To “clean” an entire seismic record containing the wavefront arrival, the entire record may be projected onto the sensor basis. This may be done as follows: First generate time-aligned, oriented traces (essentially h.sub.m(t), but not windowed),
k.sub.m(t)=g.sub.m(t+τ.sub.m)β.sub.m, m=1,2, . . . ,M, (8)
and form the matrix K using k.sub.m(t) as its t, mth row-column element. Then project the traces onto the first L columns of V, denoted V.sub.L,
{circumflex over (R)}−KV.sub.LV.sub.L.sup.1. (9)
[0042] In this way, the entire duration of the traces K is processed, not just the time window used to find V. Note that the estimated traces {circumflex over (R)} will match the estimated arrival waveforms Ŝ in the arrival time window. Finally, the estimated waveforms {circumflex over (R)} are time shifted back to their original moveout, and each sensor signal is rotated according to its polarization to generate multi-component “cleaned” traces. As an example, processed versions of the raw trace data appearing in
[0043] Trace estimates 902 made using L=5 are shown in
[0044] One question concerns the nature of the waveform estimation errors in Ŝ or {circumflex over (R)}. It is expected that the estimation errors to 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
[0045] The estimated arrival waveforms can be used to generate improved event arrival-time and polarization estimates for each sensor. Event arrival polarizations may be estimated using the projection of the estimated arrival waveforms Ŝ onto the corresponding sensor's multi-component trace data:
β.sub.m=[g.sub.m(1+τ.sub.m).sup.1 . . . g.sub.m(T+τ.sub.m).sup.1]Ŝ.sub.m, (10)
β.sub.m=ρ.sub.m/∥ρ.sub.m∥, (11)
where ŝ.sub.m is the estimated arrival waveform at the mth sensor. Similarly, 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.
[0046] The case of planar motion (two-dimensional s.sub.m(t) and β.sub.m) may be accommodated by first solving for the first column of β.sub.m1, then subtracting that direction ŝβ.sub.1.sup.T from g.sub.m, and finally repeating the 1-D process with that first dimension of motion removed.
[0047] Improved event arrival waveform estimates may be found iteratively by first using the process above to estimate the arriving waveforms, and then using those waveforms to form new arrival-time and polarization estimates. The new arrival times and polarizations are then used to generate a new set of arrival waveform estimates, and so on until the process converges, or for a given number of iterations.
[0048] Another embodiment 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 {tilde over (h)}.sub.m(t) the time-aligned and windowed component signals of sensor m, evaluated at sample time t,
{tilde over (h)}.sub.m(t)=w(t).Math.g.sub.m(t+τ.sub.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 {tilde over (H)} the T×NM matrix of time-aligned and windowed multi-component sensor signals,
where the multi-component arrival waveform matrix {tilde over (H)} is N-times as wide as the oriented arrival waveform matrix H.
[0049] 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,
UDV.sup.1={tilde over (H)}, (14)
with Ũ, {tilde over (D)} and {tilde over (V)} and 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,
{tilde over (S)}=Ũ.sub.L{tilde over (W)}.sub.L, (15)
where Ũ.sub.L is the matrix made from the first L columns of Ũ, and {tilde over (W)}.sub.l is the first L rows of {tilde over (D)}{tilde over (V)}.sup.τ. Similarly, any duration of trace data may be projected onto the columns of {tilde over (V)}.sub.l (as above) to estimate the arrival signal from an entire record.
[0050] Note that the N columns of {tilde over (W)}.sub.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 Ũ.sub.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 {tilde over (W)}.sub.L associated with each sensor by either a rank one or rank two approximation, depending on the arrival, for instance found by singular value decomposition of each sensor's L×N portion of {tilde 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.
[0051] Data compression may be achieved in a number of ways. In one method, arrivals of seismic energy contained in the trace data are detected, their parameters such as polarizations and moveouts estimated, and their waveforms modeled using a small, but sufficient number of basis shapes as described above. The resulting basis shapes and other parameters, which represent a potentially significant reduction in data needed to approximate the original sensor data or desired components therein, are then compressed using a lossless or nearly lossless data compression scheme.
[0052] The trace data may be compressed by retaining only L signal subspace components of the original trace data, that is by storing V.sub.L and KV.sub.L along with the associated arrival times and polarizations. Doing so requires only about L/(qN) the amount of storage needed for the original data, q being the number of axes measured for each trace, e.g. q−1 for a DAS. For the example here, one can achieve a data reduction by a factor of 10 for the 20 traces of data by retaining 2 components. Furthermore, standard lossless or nearly lossless data compression methods may be used to achieve further data reduction by encoding the weights and basis waveforms V.sub.L and KV.sub.L using a lossless or nearly lossless data compression scheme.
[0053] Note that in certain scenarios, it may be desired to transmit or store a “low-resolution” version of the trace data using little bandwidth or space, and then transmit additional data to achieve a more complete rendering of the seismic data. This can be achieved by first transmitting an initial number of basis vectors and weights (the most energetic dimensions of the signal subspace, e.g.), and then transmitting more of the basis vectors and weights.
[0054] Alternatively, data compression may be achieved using a convolutional model to separate multiple simultaneously arriving waveforms from a noise background. In these cases, there will be a number of estimated sources, each source formed according to a set of basis waveforms in time that are weighted according to the sensor number, and delayed according to an arrival time and moveout.
[0055] Data compression is also possible using SVD, NNMF, PLCA or other matrix or tensor factorization methods to extract a low-dimensional signal space from a set of traces directly without explicitly estimating features such as moveouts of separate seismic sources.
[0056] Another method of achieving data compression takes the matrix of traces over time, and forms, for each trace, a time-frequency analysis such as a spectrogram or wavelet transform as a pre-processing step to the above SVD, PLCA, or related method. The time-frequency analysis reveals frequencies at which energy appears across neighboring traces, and can lead to compact representations of the underlying trace data.
[0057] Finally, only certain frequency bands may be of interest, or different frequency bands may be useful for different purposes. In compressing the trace data for transmission or storage, the different frequency bands may be separately processed using the techniques described herein via band-splitting or similar filters, such as those used in perceptual audio coding or subband image coding. Note that each frequency band processed may be downsampled according to its bandwidth, further reducing the data storage space and transmission bandwidth requirements.
[0058] Additionally, different subbands may have different associated rates of compression, e.g., using different subspace dimensions L. An example of such a system is shown in
[0059] An example flowchart of the above-described methodology is shown in
[0060] The herein described subject matter sometimes illustrates different components contained within, or connected with, different other components. It is to be understood that such depicted architectures are illustrative, and that in fact many other architectures can be implemented which achieve the same functionality. In a conceptual sense, any arrangement of components to achieve the same functionality is effectively “associated” such that the desired functionality is achieved. Hence, any two components herein combined to achieve a particular functionality can be seen as “associated with” each other such that the desired functionality is achieved, irrespective of architectures or intermedial components. Likewise, any two components so associated can also be viewed as being “operably connected,” or “operably coupled,” to each other to achieve the desired functionality, and any two components capable of being so associated can also be viewed as being “operably coupleable,” to each other to achieve the desired functionality. Specific examples of operably coupleable include but are not limited to physically mateable and/or physically interacting components and/or wirelessly interactable and/or wirelessly interacting components and/or logically interacting and/or logically interactable components.
[0061] With respect to the use of plural and/or singular terms herein, those having skill in the art can translate from the plural to the singular and/or from the singular to the plural as is appropriate to the context and/or application. The various singular/plural permutations may be expressly set forth herein for sake of clarity.
[0062] It will be understood by those within the art that, in general, terms used herein, and especially in the appended claims (e.g., bodies of the appended claims) are generally intended as “open” terms (e.g., the term “including” should be interpreted as “including but not limited to,” the term “having” should be interpreted as “having at least,” the term “includes” should be interpreted as “includes but is not limited to,” etc.).
[0063] Although the figures and description may illustrate a specific order of method steps, the order of such steps may differ from what is depicted and described, unless specified differently above. Also, two or more steps may be performed concurrently or with partial concurrence, unless specified differently above. Such variation may depend, for example, on the software and hardware systems chosen and on designer choice. All such variations are within the scope of the disclosure. Likewise, software implementations of the described methods could be accomplished with standard programming techniques with rule-based logic and other logic to accomplish the various connection steps, processing steps, comparison steps, and decision steps.
[0064] It will be further understood by those within the art that if a specific number of an introduced claim recitation is intended, such an intent will be explicitly recited in the claim, and in the absence of such recitation, no such intent is present. For example, as an aid to understanding, the following appended claims may contain usage of the introductory phrases “at least one” and “one or more” to introduce claim recitations. However, the use of such phrases should not be construed to imply that the introduction of a claim recitation by the indefinite articles “a” or “an” limits any particular claim containing such introduced claim recitation to inventions containing only one such recitation, even when the same claim includes the introductory phrases “one or more” or “at least one” and indefinite articles such as “a” or “an” (e.g., “a” and/or “an” should typically be interpreted to mean “at least one” or “one or more”); the same holds true for the use of definite articles used to introduce claim recitations. In addition, even if a specific number of an introduced claim recitation is explicitly recited, those skilled in the art will recognize that such recitation should typically be interpreted to mean at least the recited number (e.g., the bare recitation of “two recitations,” without other modifiers, typically means at least two recitations, or two or more recitations).
[0065] Furthermore, in those instances where a convention analogous to “at least one of A, B, and C, etc.” is used, in general such a construction is intended in the sense one having skill in the art would understand the convention (e.g., “a system having at least one of A, B, and C” would include but not be limited to systems that have A alone, B alone, C alone, A and B together, A and C together, B and C together, and/or A, B, and C together, etc.). In those instances where a convention analogous to “at least one of A, B, or C, etc.” is used, in general, such a construction is intended in the sense one having skill in the art would understand the convention (e.g., “a system having at least one of A, B, or C” would include but not be limited to systems that have A alone, B alone, C alone, A and B together, A and C together, B and C together, and/or A, B, and C together, etc.). It will be further understood by those within the art that virtually any disjunctive word and/or phrase presenting two or more alternative terms, whether in the description, claims, or drawings, should be understood to contemplate the possibilities of including one of the terms, either of the terms, or both terms. For example, the phrase “A or B” will be understood to include the possibilities of “A” or “B” or “A and B.”
[0066] Further, unless otherwise noted, the use of the words “approximate,” “about,” “around,” “substantially,” etc., mean plus or minus ten percent.
[0067] Although the present embodiments have been particularly described with reference to preferred examples 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.