A waveform correction method and device based on time-frequency domain adaptive filtering

Abstract

The present invention discloses a waveform correction method and device based on time-frequency domain adaptive filtering. The method comprises: obtaining the P-wave velocity ratio of an underground medium, selecting a reference P-wave velocity ratio, calculating the wavelet and the reference wavelet of the matched converted wave data at each time sampling point, obtaining the amplitude spectrum of the wavelet and the reference wavelet at each time sampling point, transforming the matched converted wave data into the time-frequency domain, constructing a time-frequency domain adaptive filter using the amplitude spectrum of the wavelet and the amplitude spectrum of the reference wavelet, multiplying the adaptive filter with the data in the time-frequency domain to obtain a filtered time-frequency spectrum, transforming the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace, thereby correcting the non-steady-state converted wave data matched to the P-wave time domain into steady-state data, thereby making the subsequent joint inversion of the P-wave and converted waves more accurate.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration, and in particular to a waveform correction method and device based on time-frequency domain adaptive filtering. Background Art

[0002] Seismic exploration refers to a geophysical exploration method that uses the differences in velocity and density of underground media caused by artificially excited elastic waves to observe and analyze the propagation patterns of seismic waves generated by artificial earthquakes underground to infer the properties and morphology of underground rock formations.

[0003] Seismic inversion is the process of imaging (solving) the spatial structure and physical properties of underground rock formations using surface seismic data, constrained by known geological laws and drilling and logging data. In a broad sense, seismic inversion encompasses the entire content of seismic processing and interpretation.

[0004] Multiwave and multicomponent exploration, also known as vector exploration, refers to an exploration technique that utilizes a combination of P- and S-wave sources and multicomponent geophones to observe various wave fields, revealing more information about subsurface structure, lithology, and oil and gas. With the advancement of multicomponent seismic acquisition technology, converted-wave exploration is increasingly being used in seismic exploration. Compared to P-waves, converted waves have greater penetration into strata. Therefore, combined inversion of P- and converted waves yields more accurate reservoir elastic parameter information than P-wave inversion alone.

[0005] However, due to the different P- and S-wave velocities in the subsurface, the P-wave and converted-wave reflected from the same interface have different travel times. Conventional joint inversion based on the Zoeppritz equation and its approximations requires matching the P-wave and converted-wave data to the same time domain, typically matching the converted-wave data to the P-wave time domain. This matching causes stretching or squeezing of the converted-wave data, rendering it no longer a steady-state seismic data, leading to inaccurate inversion results. Summary of the Invention

[0006] The present invention provides a waveform correction method based on time-frequency domain adaptive filtering to solve the problem in the prior art that when converted wave data is matched to the longitudinal wave time domain, the converted wave data is no longer steady-state seismic data, which leads to inaccurate inversion results. The method includes:

[0007] Extract converted wave seismic wavelet;

[0008] Matching converted wave data to the P-wave time domain;

[0009] Obtaining the ratio of the longitudinal and shear wave velocities of the underground medium, and based on Fourier scaling theory, using the ratio of the longitudinal and shear wave velocities of the underground medium and the converted wave seismic wavelet, calculating the wavelet of the matched converted wave data at each time sampling point, and obtaining the amplitude spectrum of the wavelet at each time sampling point;

[0010] Selecting a reference P-wave and S-wave velocity ratio, and based on Fourier scaling transform theory, using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio, calculating a reference wavelet of the matched converted wave data, and obtaining an amplitude spectrum of the reference wavelet;

[0011] Transforming the matched converted wave data into the time-frequency domain to obtain data in the time-frequency domain, and constructing a time-frequency domain adaptive filter using the amplitude spectrum of the wavelet and the amplitude spectrum of the reference wavelet at each time sampling point obtained;

[0012] The adaptive filter is multiplied by the data in the time-frequency domain to obtain a filtered time-frequency spectrum, and the filtered time-frequency spectrum is converted back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

[0013] In one embodiment, matching the converted wave data to the longitudinal wave time domain includes:

[0014] Match the converted wave data to the longitudinal wave time domain to obtain the signal s(t);

[0015] The signal s(t) satisfies the following relationship:

[0016] FT f S(t / β)=|β|S0(βf)

[0017] Among them, FT f S(t / β) is the spectrum of the compressed or stretched signal s(t), β is the compression or stretching coefficient, and S0(f) is the Fourier transform of the original signal s0(t).

[0018] In one embodiment, the amplitude spectrum S(f) of the wavelet of the compressed or stretched signal s(t) is:

[0019] S(f)=|β|S0(βf)

[0020] Where β is the compression or stretch coefficient, and S(f) is the amplitude spectrum of the wavelet.

[0021] In one embodiment, the compression or tension coefficient β is obtained by the following formula:

[0022]

[0023] Where γ is the ratio of the longitudinal and transverse wave velocities of the underground medium.

[0024] In one embodiment, the longitudinal and transverse wave velocity ratios include a reference longitudinal and transverse wave velocity ratio γ ref , the compression or tension coefficient includes a reference compression or tension coefficient β ref , the reference longitudinal and transverse wave velocity ratio γ ref and the reference compression or tension coefficient β ref Satisfy the formula

[0025] In one embodiment, selecting a reference longitudinal and transverse wave velocity ratio includes:

[0026] The P-wave velocity ratio at the target layer or the average P-wave velocity ratio in the exploration area is selected as the reference P-wave velocity ratio.

[0027] In one embodiment, a wavelet shaping filter is obtained based on the obtained amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet:

[0028]

[0029] Among them, S ref is the amplitude spectrum of the reference wavelet, S pred is the amplitude spectrum of the wavelet, and ε is a constant.

[0030] In one embodiment, after the matched converted wave data is transformed into the time-frequency domain, the time-frequency domain adaptive filter is:

[0031]

[0032] Where α is the scale factor, 0≤α≤1, N is the filter order, is the Hadamard product of matrices, is a box function.

[0033] In one embodiment, the filtered converted wave seismic trace is:

[0034]

[0035] Among them, ifft is the inverse Fourier transform, is the filtered time-frequency domain data.

[0036] The embodiment of the present invention further provides a waveform correction device based on time-frequency domain adaptive filtering, comprising:

[0037] Extraction module, extracting converted wave seismic wavelet;

[0038] Conversion module, used to match converted wave data to the longitudinal wave time domain;

[0039] an acquisition calculation module for obtaining the P-wave and S-wave velocity ratio of the underground medium, and calculating the wavelet of the matched converted wave data at each time sampling point using the P-wave and S-wave velocity ratio of the underground medium and the converted wave seismic wavelet based on the Fourier scaling theory, and obtaining the amplitude spectrum of the wavelet at each time sampling point; selecting a reference P-wave and S-wave velocity ratio, and calculating the reference wavelet of the matched converted wave data based on the Fourier scaling theory using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio, and obtaining the amplitude spectrum of the reference wavelet;

[0040] A construction module is used to transform the matched converted wave data into the time-frequency domain to obtain time-frequency domain data, and construct a time-frequency domain adaptive filter using the obtained amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet;

[0041] The processing module is used to multiply the adaptive filter by the time-frequency domain data to obtain a filtered time-frequency spectrum, and transform the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

[0042] In an embodiment of the present invention, a waveform correction method and device based on time-frequency domain adaptive filtering are provided. Compared with the technical solutions in the prior art, the method and device obtain the P-wave velocity ratio of the underground medium, select a reference P-wave velocity ratio, calculate the wavelet and reference wavelet of the matched converted wave data at each time sampling point, obtain the amplitude spectrum of the wavelet and the reference wavelet at each time sampling point, transform the matched converted wave data into the time-frequency domain, construct a time-frequency domain adaptive filter using the amplitude spectrum of the wavelet and the amplitude spectrum of the reference wavelet, multiply the adaptive filter with the time-frequency domain data to obtain a filtered time-frequency spectrum, and transform the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace. In this way, the non-steady-state converted wave data matched to the P-wave time domain is corrected into steady-state data, thereby making the subsequent joint inversion of P-waves and converted waves more accurate. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative work. In the drawings:

[0044] Figure 1 Flowchart of a waveform correction method based on time-frequency domain adaptive filtering in one embodiment of the present invention;

[0045] Figure 2Flowchart of a waveform correction method based on time-frequency domain adaptive filtering in one embodiment of the present invention;

[0046] Figure 3 The input actual logging data and the estimated P-wave and S-wave velocity ratios in one embodiment of the present invention;

[0047] Figures 4A-4E According to the embodiment of the present invention Figure 2 The well logging data shown are angle gathers synthesized using the exact Zoeppritz equation;

[0048] Figures 5A-5C PS seismic wavelet used in one embodiment of the invention and reference wavelet when the reference P-wave and S-wave velocity ratios are 2.5 and 4, respectively;

[0049] Figures 6A-6C The filtering results obtained when the reference longitudinal and transverse wave velocity ratio is 2.5 are selected for one embodiment of the present invention;

[0050] Figures 7A-7C The filtering result obtained when the reference longitudinal and transverse wave velocity ratio is 4 is selected for an embodiment of the present invention;

[0051] Figure 8 The inversion results of P-wave and S-wave velocities and densities obtained by performing a joint inversion based on waveform-corrected data according to an embodiment of the present invention;

[0052] Figure 9 Schematic diagram of a waveform correction device based on time-frequency domain adaptive filtering in one embodiment of the present invention. DETAILED DESCRIPTION

[0053] To make the purpose, technical solutions and advantages of the embodiments of the present invention more clear, the embodiments of the present invention are further described in detail below with reference to the accompanying drawings. Here, the exemplary embodiments of the present invention and their descriptions are used to explain the present invention, but are not intended to limit the present invention.

[0054] Multiwave and multicomponent exploration, also known as vector exploration, refers to an exploration technique that utilizes a combination of P- and S-wave sources and multicomponent geophones to observe various wave fields, revealing more information about subsurface structure, lithology, and oil and gas. With the advancement of multicomponent seismic acquisition technology, converted-wave exploration is increasingly being used in seismic exploration. Compared to P-waves, converted waves have greater penetration into strata. Therefore, combined inversion of P- and converted waves yields more accurate reservoir elastic parameter information than P-wave inversion alone.

[0055] However, due to the different P- and S-wave velocities of the subsurface medium, the P-wave and converted-wave reflected from the same interface underground have different travel times. Conventional joint inversion based on the Zoeppritz equation and its approximations requires matching the P-wave and converted-wave data to the same time domain, typically matching the converted-wave data to the P-wave time domain. This matching causes stretching or squeezing of the converted-wave data, rendering it no longer a steady-state seismic data. Extracting a statistical wavelet from the entire seismic data and then using it for inversion will affect the inversion results.

[0056] One approach to addressing the non-stationary nature of the matched converted wave data is to perform inversion using time windows. However, when the P-wave and S-wave velocity ratios of the subsurface medium vary dramatically, it is difficult to find a satisfactory time window within which the seismic data meet the steady-state assumption. Another approach is to use Fourier scaling theory to process the data. First, a steady-state seismic wavelet is extracted from the original converted wave data. Then, based on the estimated P-wave and S-wave velocity ratios, the seismic wavelet at each time after compression is calculated using Fourier scaling theory. A reference P-wave and S-wave velocity ratio is also selected, and a reference seismic wavelet is also calculated. A zero-phase frequency-domain filter is then constructed using the reference wavelet and the amplitude spectrum of the estimated seismic wavelet at each time. Finally, this filter is applied to the compressed converted wave data to obtain a filtered data volume. However, this method is susceptible to random noise.

[0057] However, the current approach to joint inversion of P-wave and converted-wave data often makes it difficult to select an appropriate time window when the P-wave and S-wave velocity ratios in the subsurface vary dramatically. Furthermore, constructing a spectral balance filter based on Fourier scaling theory can amplify random noise in the raw data.

[0058] The present invention provides a waveform correction method based on time-frequency domain adaptive filtering, such as Figure 1 As shown, the following steps may be included:

[0059] Step 101: extracting converted wave seismic wavelet;

[0060] It should be noted that the method for extracting the converted wave seismic wavelet can be a deterministic wavelet extraction method or a statistical wavelet extraction method, and this application does not specifically limit the seismic wavelet extraction method.

[0061] Step 102: Match the converted wave data to the longitudinal wave time domain;

[0062] It should be noted that because the speed of converted waves is lower than that of P-waves, the arrival times of the same reflector on converted and P-wave profiles are inconsistent. Therefore, time registration is required during seismic data interpretation. Time registration involves compressing the time on the converted wave profile to the P-wave time based on horizon comparison, using the arrival time of the P-wave profile as a standard and referencing several key target reflection horizons. This allows the P-wave and converted waves to be on the same timeline, facilitating comparative interpretation of wave groups.

[0063] According to Fourier scaling theory, when the time domain signal s0(t) is compressed or stretched, its Fourier transform corresponds to stretching or squeezing, and the amount of squeezing or stretching corresponds to the following relationship:

[0064] FT f S(t / β)=|β|S0(βf) (1)

[0065] Among them, FT f S(t / β) is the spectrum of the compressed or stretched signal s(t), β is the compression or stretching coefficient, and S0(f) is the Fourier transform of the signal s0(t).

[0066] To calculate the wavelet and its amplitude spectrum at each time sampling point of the matched converted wave data, it is necessary to extract the seismic wavelet of the converted wave from the seismic data. A seismic wavelet is a signal with a definite start time, limited energy, and a certain duration; it is the basic unit of a seismic record. Seismic wavelets can be extracted using either deterministic or statistical wavelet extraction methods; this application does not specifically limit the extraction method.

[0067] Step 103: Obtain the P-wave and S-wave velocity ratio of the underground medium, and based on Fourier scaling theory, use the P-wave and S-wave velocity ratio of the underground medium and the converted wave seismic wavelet to calculate the wavelet of the matched converted wave data at each time sampling point, and obtain the amplitude spectrum of the wavelet at each time sampling point;

[0068] It should be noted that the present application does not specifically limit the method for obtaining the ratio of longitudinal and transverse wave velocities and the method for obtaining the velocity.

[0069] After obtaining the P-wave and S-wave velocity ratio of the underground medium, based on the Fourier scaling theory, the P-wave and S-wave velocity ratio of the underground medium and the converted wave seismic wavelet are used to calculate the wavelet of the matched converted wave data at each time sampling point, thereby obtaining the amplitude spectrum of the wavelet at each time sampling point. And from formula (1), it can be concluded that the Fourier transform S(f) of the compressed or stretched signal s(t), that is, the amplitude spectrum of the wavelet, is:

[0070] S(f)=|β|S0(βf) (2)

[0071] Where β is the compression or stretch coefficient, and S0(f) is the Fourier transform of the signal s0(t).

[0072] According to this relationship, the spectrum of the compressed or stretched signal can be obtained in the frequency domain, and then the compressed or stretched signal can be obtained by using inverse Fourier transform.

[0073] It should be noted that the compression or stretch coefficient β is obtained by the following formula:

[0074]

[0075] Where γ is the ratio of the longitudinal and transverse wave velocities of the underground medium.

[0076] Step 104: Select a reference P-wave and S-wave velocity ratio, and based on Fourier scaling theory, use the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio to calculate a reference wavelet of the matched converted wave data, and obtain an amplitude spectrum of the reference wavelet.

[0077] It should be noted that the selection of the reference P-wave and S-wave velocity ratio includes: selecting the P-wave and S-wave velocity ratio at the target layer or the average P-wave and S-wave velocity ratio in the exploration area as the reference P-wave and S-wave velocity ratio.

[0078] After selecting the reference P-wave and S-wave velocity ratio, based on the Fourier scaling theory, the reference wavelet of the converted wave data after matching is calculated using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio to obtain the amplitude spectrum of the reference wavelet. And from formula (1), it can be concluded that the reference wavelet s ref The Fourier transform S of (t) ref (f) That is, the spectrum of the wavelet is:

[0079] S ref (f)=|β ref |S0(β ref f) (4)

[0080] Among them, β ref is the reference compression or stretch coefficient, and S0(f) is the Fourier transform of the signal s0(t).

[0081] According to this relationship, the spectrum of the compressed or stretched signal can be obtained in the frequency domain, and then the compressed or stretched signal can be obtained by using inverse Fourier transform.

[0082] It should be noted that the reference compression or tension coefficient β ref Obtained by the following formula:

[0083]

[0084] Among them, γ ref is the reference P-wave velocity ratio.

[0085] Step 105: transforming the matched converted wave data into the time-frequency domain to obtain data in the time-frequency domain, and constructing a time-frequency domain adaptive filter using the amplitude spectrum of the wavelet and the amplitude spectrum of the reference wavelet at each time sampling point.

[0086] It should be noted that after obtaining the amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet, the following wavelet shaping filter can be constructed:

[0087]

[0088] Among them, S ref is the amplitude spectrum of the reference wavelet, S pred is the amplitude spectrum of the wavelet at each time sampling point calculated by formula (2), ε is a constant, which is used to prevent S from being pred is 0, which leads to numerical instability of the sub-shaping filter F. Usually, this value is taken as S pred The maximum value is 0.000001.

[0089] It should be noted that, as can be seen from Equation (6), the filtering process is equivalent to a spectrum balancing process. When the P-wave velocity ratio γ of the underground medium is greater than the reference P-wave velocity ratio γ re f region, the filter is similar to a low-pass filter, and the high-frequency components in the original data will be reduced accordingly. When the longitudinal and transverse wave velocity ratio γ of the underground medium is less than the reference longitudinal and transverse wave velocity ratio γ re When f, the filter is similar to a high-pass filter, and the high-frequency components in the original data will increase.

[0090] Step 106: multiply the adaptive filter by the time-frequency domain data to obtain a filtered time-frequency spectrum, and transform the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

[0091] Considering the need to achieve a higher resolution of the filtered data, reduce artifacts, and improve the quality of the filtered data, as shown in equation (6), when we select a very high reference P-wave and S-wave velocity ratio, the increase in high-frequency components will amplify the noise energy in the original signal, especially random noise. Even when the noise energy in the original data is very weak, an unreasonably high reference P-wave and S-wave velocity ratio will produce some artifacts, thereby affecting the quality of the filtered data. Therefore, we usually select the P-wave and S-wave velocity ratio at the target layer as the reference P-wave and S-wave velocity ratio.

[0092] It should be noted that even if we select a reasonable reference P-wave velocity ratio, the high-frequency noise in the original data may be amplified. When the reference P-wave velocity ratio is greater than the P-wave velocity ratio at the target layer, the filter defined by Equation (6) will amplify the random noise in the original data. Therefore, it is very necessary to propose a filter that can suppress random noise and perform wavelet shaping.

[0093] To solve this problem, we first transform the matched PS wave data into the time-frequency domain and then filter the PS wave data in the time-frequency domain. In this embodiment, the time-frequency analysis technique we selected is short-time Fourier transform. Of course, other higher-resolution time-frequency analysis techniques can also be selected to improve the filtering effect. The present invention does not specifically limit the time-frequency analysis technique selected.

[0094] According to the short-time Fourier transform, the short-time Fourier transform of the time domain signal x(t) is:

[0095]

[0096] Among them, w(t) is the window function, which satisfies

[0097] After transforming the signal into the time-frequency domain, we can set an adaptive filter in the time-frequency domain. Assuming that the PS wave data matched to the PP time domain is s(t), its short-time Fourier transform is S(ω,τ), then we can set the following filter in the time-frequency domain:

[0098]

[0099] Among them, α is the scale factor, and 0≤α≤1, N is the filter order, is the Hadamard product of matrices, is a box function, and its expression is:

[0100]

[0101] Among them, WL is the time length of the box function. It can be seen from formula (8) that compared with the traditional wavelet shaping filter shown in formula (6), this filter adds an adaptive term, which is similar to a Butterworth filter, that is, the amplitude in the time spectrum of the original data that is less than α times the maximum value of the amplitude spectrum is considered to be random noise and is then suppressed. In this way, by selecting a suitable scale factor α and filter order, it is possible to suppress the amplification of random noise while shaping the wavelet. In addition, when the filter order N is 0, the filter degenerates into a traditional wavelet shaping filter, except for a constant 2. After the above filter is applied to the PS wave data in the time-frequency domain, the final filtered and corrected data fs(t) can be obtained by the following formula:

[0102]

[0103] Among them, ifft is the inverse Fourier transform, is the filtered time-frequency domain data, which can be expressed as:

[0104]

[0105] like Figure 2 As shown, the present invention provides a waveform correction method based on time-frequency domain adaptive filtering, which may include the following steps:

[0106] Step 201: extracting seismic wavelets from PS data in the PS time domain, and converting the PS data into the PP time domain using the acquired subsurface P-wave and S-wave velocity ratios;

[0107] Step 202: using the extracted PS seismic wavelet and the underground P-wave velocity ratio, calculate the seismic wavelet of the matched PS data at each time sampling point according to Fourier scaling theory;

[0108] Step 203: Calculate a reference wavelet using the selected reference P-wave and S-wave velocity ratios, and construct a wavelet shaping filter at each sampling point using the calculated seismic wavelet at each sampling point;

[0109] Step 204: Convert the matched PS data into the time-frequency domain using a time-frequency transformation technique, and construct an adaptive filter in the time-frequency domain;

[0110] Step 205: Multiply the constructed adaptive filter and the wavelet shaping filter, and apply the result to the PS data in the time-frequency domain. Then, perform inverse time-frequency transformation to obtain the filtered result.

[0111] like Figure 3As shown, it is the actual logging data used in the embodiment of the present invention, from left to right are P-wave velocity (Vp), S-wave velocity (Vs), density (Density), P-wave velocity ratio (Vp / Vs) and estimated P-wave velocity ratio (Estimated Vp / Vs). The vertical axis Time represents time. In the embodiment of the present invention, the PS and PP data matching is completed using the dynamic time warping algorithm, and the P-wave velocity ratio can be estimated, but this does not belong to the research content of the embodiment of the present invention, and the present invention does not specifically limit the PS and PP data matching algorithm. From the comparison of the estimated P-wave velocity ratio and the true P-wave velocity ratio, the trend of the estimated P-wave velocity ratio is the same as the true P-wave velocity ratio, so the P-wave velocity ratio can be used for subsequent processing.

[0112] Figures 4A-4E The embodiment of the present invention is based on Figure 3 The well log data shown are synthesized seismic data using the exact Zoeppritz equation, where Figure 4A is the PP wave data, Figure 4B PS wave data in the PS time domain, Figure 4C To match the PS wave data in the PP time domain, Figure 4D is the standard PS wave data when the reference P-wave velocity ratio is 2.5, Figure 4E is the standard PS wave data when the reference longitudinal-to-short wave velocity ratio is 4, and the standard PS wave data is obtained by convolving the reference wavelet with the PS reflection coefficient in the PP time domain. Figures 4A-4E In the figure, the vertical axis Time represents time, and the horizontal axis Angle represents angle.

[0113] Figures 5A-5C is the seismic wavelet used in this embodiment, where Figure 5A is the original PS wavelet, Figure 5B Select the seismic wavelet when the reference P-wave velocity ratio is 2.5, Figure 5C Seismic wavelets are selected for a P-wave to S-wave velocity ratio of 4. As can be seen from the figure, the main events of the PS wave data after matching to the PP time domain essentially correspond to those of the PP wave data. In this example, the PS wave reflection coefficient and the PP wave reflection coefficient have opposite polarities. After matching to the PP time domain, the wavelets of the PS wave data vary significantly from shallow to deep layers. Directly using this data for joint inversion and extracting a statistical wavelet would inevitably affect the inversion results, necessitating wavelet correction. Figures 5A-5C In the figure, the vertical axis Amplitude represents the amplitude, and the horizontal axis Time represents the time.

[0114] Figures 6A-6Cis the filtering result obtained when the reference P-wave velocity ratio is 2.5. Figures 6A-6C Only the superposition data volume is shown in FIG, and the PS wave superposition data matched to the PP time domain is added with random noise with a signal-to-noise ratio of 10 to verify the effect of the invention solution. Figure 3 Judging from the estimated actual P-wave and S-wave velocity ratio of the strata, this value is almost smaller than the P-wave and S-wave velocity ratio of all layers in the study area. Therefore, theoretically, sampling this value for wavelet shaping filtering will not amplify the random noise in the original data.

[0115] Figure 6A The PS wave data matched to the PP time domain is compared with the standard trace. As can be seen from the figure, due to the compression of the original wavelet, the matched seismic trace is quite different from the standard seismic trace. Figure 6B and Figure 6C The filtering results obtained by the traditional method and the embodiment of the present invention are respectively. It can be seen from the figure that when the reference P-wave and S-wave velocity ratio is selected as 2.5, both methods obtain a good filtering result. Compared with the unfiltered data (Perfect), the PS wave data after filtering (Filtered) is more consistent with the standard trace.

[0116] Figures 7A-7C The filtering result is obtained when the reference P-wave velocity ratio is 4. Because this reference P-wave velocity ratio is greater than the P-wave velocity ratios in the studied area, random noise is clearly amplified in the traditional filtering results. If this data is used for subsequent joint inversion, the inversion result will inevitably contain excessive numerical noise. The results obtained by the present invention correct the wavelet waveform without excessively amplifying the random noise in the original data.

[0117] Figure 8 The P-wave and S-wave velocity and density inversion results obtained by performing a joint P-wave and S-wave inversion based on the filtered data in an embodiment of the present invention, where the black solid line represents the true P-wave and S-wave velocity and density values, and the black dotted line represents the inversion result. The inversion result is very consistent with the true result and does not introduce obvious numerical noise.

[0118] The present invention provides a waveform correction method based on time-frequency domain adaptive filtering. Compared with the traditional waveform correction method, which uses the amplitude spectrum of the corrected wavelet and the amplitude spectrum of the reference wavelet to construct a wavelet shaping filter and then performs filtering correction on the matched PS wave data, this method will amplify the random noise contained in the original data when the selected reference P-wave velocity ratio is greater than the actual P-wave velocity ratio of the formation. In addition, even if there is no random noise in the original data, this method will introduce a certain amount of numerical noise due to the amplification of the high-frequency components. The embodiment of the present invention adds an adaptive filtering item to the traditional wavelet shaping filter. This adaptive filtering is performed in the time-frequency domain and is adaptively filtered according to the amplitude of the original signal in the time-frequency domain, thereby suppressing the random noise in the original data. In this way, the obtained filtering result corrects the waveform of the wavelet without amplifying the high-frequency noise in the original data, so that a high-quality filtering result can be obtained for subsequent joint inversion of P-waves and S-waves.

[0119] Based on the above invention concept, Figure 9 As shown, the present invention also proposes a waveform correction device based on time-frequency domain adaptive filtering, comprising:

[0120] Extraction module 901, extracting converted wave seismic wavelet;

[0121] Conversion module 902, for matching converted wave data to the longitudinal wave time domain;

[0122] The acquisition calculation module 903 is used to obtain the P-wave and S-wave velocity ratio of the underground medium. Based on the Fourier scaling transform theory, the P-wave and S-wave velocity ratio of the underground medium and the converted wave seismic wavelet are used to calculate the wavelet at each time sampling point of the converted wave data after matching, and obtain the amplitude spectrum of the wavelet at each time sampling point.

[0123] Selecting a reference P-wave and S-wave velocity ratio, and based on Fourier scaling transform theory, using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio, calculating a reference wavelet of the matched converted wave data, and obtaining an amplitude spectrum of the reference wavelet;

[0124] A construction module 904 is configured to transform the matched converted wave data into the time-frequency domain, and construct a time-frequency domain adaptive filter using the obtained amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet;

[0125] The processing module 905 is configured to multiply the adaptive filter by the time-frequency domain data to obtain a filtered time-frequency spectrum, and convert the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

[0126] It should be noted that the above-mentioned device may also include other implementations according to the description of the method embodiment. Specific implementations can refer to the description of the relevant method embodiment and will not be described in detail here.

[0127] The methods or devices described in the above embodiments provided in this specification can implement business logic through computer programs and record them on a storage medium, which can be read and executed by a computer to achieve the effects of the solutions described in the embodiments of this specification. Therefore, this specification also provides an electronic device for use in a server, which may include a processor and a memory storing processor-executable instructions. When the instructions are executed by the processor, the steps of the method described in any of the above embodiments are implemented.

[0128] The storage medium may include a physical device for storing information, typically digitizing the information and then storing it in a medium utilizing electrical, magnetic, or optical means. Examples of such storage media include: devices that use electrical energy to store information, such as various types of memory, such as RAM and ROM; devices that use magnetic energy to store information, such as hard disks, floppy disks, magnetic tapes, magnetic core memories, bubble memories, and USB flash drives; and devices that use optical means to store information, such as CDs or DVDs. Of course, there are also other types of readable storage media, such as quantum memories and graphene memories.

[0129] The embodiments of this specification are not limited to situations that must conform to standard data models / templates or the embodiments of this specification. Certain industry standards or slightly modified implementation plans based on the implementation described in the custom methods or embodiments can also achieve the same, equivalent or similar implementation effects as the above embodiments, or the expected implementation effects after deformation. The embodiments obtained by applying these modified or deformed data acquisition, storage, judgment, processing methods, etc. can still fall within the scope of the optional implementation plans of this specification.

[0130] The various embodiments in this specification are described in a progressive manner. Similar or identical parts between the various embodiments can be referenced to each other, and each embodiment focuses on the differences from other embodiments. In particular, since the device embodiments are generally similar to the method embodiments, the description is relatively simple, and relevant parts can be referenced to the partial description of the method embodiments. In the description of this specification, reference to the terms "one embodiment," "some embodiments," "examples," "specific examples," or "some examples" means that the specific features, structures, materials, or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of this specification. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in any appropriate manner in any one or more embodiments or examples. In addition, those skilled in the art may combine and combine the different embodiments or examples described in this specification, as well as the features of different embodiments or examples, without conflict.

[0131] The foregoing is merely an example of the present invention and is not intended to limit the present invention. Various modifications and variations are possible for those skilled in the art. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention are intended to be included within the scope of the claims of the present invention.

Claims

1. A waveform correction method based on time-frequency domain adaptive filtering, characterized in that: include: Extract converted wave seismic wavelet; Matching converted wave data to the P-wave time domain; Obtaining the ratio of the longitudinal and shear wave velocities of the underground medium, and based on Fourier scaling theory, using the ratio of the longitudinal and shear wave velocities of the underground medium and the converted wave seismic wavelet, calculating the wavelet of the matched converted wave data at each time sampling point, and obtaining the amplitude spectrum of the wavelet at each time sampling point; Selecting a reference P-wave and S-wave velocity ratio, and based on Fourier scaling transform theory, using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio, calculating a reference wavelet of the matched converted wave data, and obtaining an amplitude spectrum of the reference wavelet; The matched converted wave data is transformed into the time-frequency domain to obtain the time-frequency domain data. The amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet are used to construct a time-frequency domain adaptive filter. The time-frequency domain adaptive filter is: α is the scale factor, 0≤α≤1, N is the filter order, o is the Hadamard product of the matrix, is a box function; WL is the time length of the box function; The adaptive filter is multiplied by the data in the time-frequency domain to obtain a filtered time-frequency spectrum, and the filtered time-frequency spectrum is converted back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

2. The waveform correction method based on time-frequency domain adaptive filtering according to claim 1, characterized in that: Matching the converted wave data to the longitudinal wave time domain includes: Match the converted wave data to the longitudinal wave time domain to obtain the signal s(t); The signal s(t) satisfies the following relationship: FT f S(t / β)=|β|S0(βf) Among them, FT f S(t / β) is the spectrum of the compressed or stretched signal s(t), β is the compression or stretching coefficient, and S0(βf) is the signal after Fourier transform of the original signal s0(t).

3. The waveform correction method based on time-frequency domain adaptive filtering according to claim 2, characterized in that: The compression or tension coefficient β is obtained by the following formula: Where γ is the ratio of the longitudinal and transverse wave velocities of the underground medium.

4. The waveform correction method based on time-frequency domain adaptive filtering according to claim 3, characterized in that: The longitudinal and transverse wave velocity ratios include a reference longitudinal and transverse wave velocity ratio γ ref , the compression or tension coefficient includes a reference compression or tension coefficient β ref , the reference longitudinal and transverse wave velocity ratio γ ref and the reference compression or tension coefficient β ref Satisfy the formula 5. The waveform correction method based on time-frequency domain adaptive filtering according to claim 4, characterized in that: The selecting of a reference longitudinal and transverse wave velocity ratio comprises: The P-wave velocity ratio at the target layer or the average P-wave velocity ratio in the exploration area is selected as the reference P-wave velocity ratio.

6. The waveform correction method based on time-frequency domain adaptive filtering according to claim 3, characterized in that: Based on the obtained amplitude spectrum of the wavelet at each time sampling point and the amplitude spectrum of the reference wavelet, a wavelet shaping filter is obtained: Among them, S ref is the amplitude spectrum of the reference wavelet, S pred is the amplitude spectrum of the wavelet, and ε is a constant.

7. The waveform correction method based on time-frequency domain adaptive filtering according to claim 1, characterized in that: The converted wave seismic trace after filtering is: Among them, ifft is the inverse Fourier transform, is the filtered time-frequency domain data.

8. A waveform correction device based on time-frequency domain adaptive filtering, characterized in that: include: Extraction module, extracting converted wave seismic wavelet; Conversion module, used to match converted wave data to the longitudinal wave time domain; an acquisition calculation module for obtaining the P-wave and S-wave velocity ratio of the underground medium, and calculating the wavelet of the matched converted wave data at each time sampling point using the P-wave and S-wave velocity ratio of the underground medium and the converted wave seismic wavelet based on the Fourier scaling theory, and obtaining the amplitude spectrum of the wavelet at each time sampling point; selecting a reference P-wave and S-wave velocity ratio, and calculating the reference wavelet of the matched converted wave data based on the Fourier scaling theory using the converted wave seismic wavelet and the reference P-wave and S-wave velocity ratio, and obtaining the amplitude spectrum of the reference wavelet; A construction module is used to transform the matched converted wave data into the time-frequency domain to obtain time-frequency domain data, and construct a time-frequency domain adaptive filter using the amplitude spectrum of the wavelet and the amplitude spectrum of the reference wavelet at each time sampling point. The time-frequency domain adaptive filter is: α is the scale factor, 0≤α≤1, N is the filter order, o is the Hadamard product of the matrix, is a box function; WL is the time length of the box function; The processing module is used to multiply the adaptive filter by the time-frequency domain data to obtain a filtered time-frequency spectrum, and transform the filtered time-frequency spectrum back to the time domain through an inverse time-frequency transform to obtain a filtered converted wave seismic trace.

Images