A method for compensating low-frequency information of seismic data and its application

Abstract

This study discloses a low-frequency compensation method for seismic data in desert areas, and establishes a new low-frequency compensation and expansion method for seismic data and an optimization method for processing parameters. It is used to enhance the energy of low-frequency components at the low-frequency end of seismic data, expand the effective frequency band range at the low-frequency end, and optimize processing parameters according to the noise level of seismic data and the original frequency band characteristics, so as to enhance the low-frequency effective information and resolution of seismic data as much as possible. This method has good applicability and promotion based on existing data, and the research ideas and technical methods are highly feasible, providing reliable seismic success data for the evaluation and deployment of fault zones and carbonate fault-controlled reservoirs in desert areas.

Description

Technical Field

[0001] The invention belongs to the field of ultra-deep carbonate fracture reservoir imaging in the desert area of ​​the Tarim Basin, and specifically relates to a seismic data low-frequency information compensation method and application thereof. Background Art

[0002] In the Shunbei area of ​​the Tarim Basin, the Ordovician carbonate reservoir is generally buried at a depth of more than 7000m. Due to the influence of surface sand dunes and stratum absorption, the energy of the high-frequency components of seismic waves decays quickly, and the signal-to-noise ratio of seismic data is generally low. Due to the small vertical fault throw of the strike-slip fault, the amplitude energy of the diffraction wave is relatively weak. However, since the deep and large strike-slip fault has a certain width and extension length, it has the amplitude-frequency response characteristics of low frequency, weak amplitude, and low wave number on the seismic profile. The low-frequency component of the seismic signal plays a vital role in the imaging of deep reservoirs (Zhang et al., 2015; Wang et al., 2017; Yang et al., 2018). During the propagation of seismic waves, the high-frequency components of the seismic signal are absorbed due to stratum absorption and energy attenuation factors. In deep or ultra-deep layers, the main component of the seismic signal is low-frequency information.

[0003] In the study of low-frequency compensation of seismic data, Lines (1977) first proposed a standard autoregressive method for spectrum extrapolation; Guan Luping (1990) proved that the spectrum of the signal can be expressed as an autoregressive process, and the low-frequency band can be predicted by using a reliable intermediate frequency band (such as 10-50Hz) to apply autoregression to the low-frequency band, and the spectrum can be extrapolated according to the principle of minimizing the forward prediction error and the reverse error. Han Liguo et al. (2012) proposed low-frequency compensation of seismic data based on compressed sensing and sparse inversion, only intercepting the low-frequency band as the low-frequency component of the expanded spectrum, and retaining the high-frequency band of the original signal. The envelope inversion proposed by Wu Rushan et al. (2014) shows that the amplitude envelope of the seismic signal contains low-frequency information. Sajid (2014) proposed an improved algorithm to expand the bandwidth of seismic data without introducing high-frequency noise. The algorithm uses a set of three cascaded difference operators to increase the high frequency and a simple smoothing operator to increase the low frequency. Woodburn et al. (2012) proposed a frequency extension method based on the deconvolution operator, and Wang Peng (2016) used heterodyne technology to achieve the recovery and compensation of low-frequency components. However, the above low-frequency extension methods all have certain limitations. When there is strong low-speed band interference, the low-frequency extension effect processed by the autoregressive method is not obvious; the frequency stabilization method based on spatial filtering only performs frequency division and mixing processing on the low-frequency band information, and the meaning of frequency compensation is unclear; and the low-frequency extension method based on the deconvolution operator only processes the sub-waves, without considering the transformation of the stratum to the sub-waves. The low-frequency compensation method based on extrapolation technology changes the spectral characteristics of the original seismic record, and the physical meaning is unclear. Zhang Junhua et al. (2015) studied the frequency extension method based on compressed sensing and sparse constraints, and processed the actual CIP gathers, achieving good application results, but the method has limited ability to recover and compensate low frequencies. Summary of the invention

[0004] Based on the above-mentioned deficiencies in the prior art in the art, the object of the present invention is to provide a method for compensating low-frequency information of desert seismic data, so as to enhance the low-frequency effective information and resolution of desert seismic data.

[0005] The technical solution of the present invention is as follows:

[0006] A method for compensating low-frequency information of seismic data, characterized by comprising:

[0007] A. Extract seismic wavelets from seismic data;

[0008] B. Obtain the reflection coefficient of seismic data through sparsity regularized inversion;

[0009] C. reshaping the seismic wavelet extracted in step A to obtain a low-frequency seismic wavelet;

[0010] D. Convolve the low-frequency seismic wavelet obtained in step C with the reflection coefficient obtained by inversion in step B.

[0011] Extracting the seismic wavelet in step A means: obtaining the frequency domain analytical expression of the seismic wavelet by deducing the Ricker wavelet, and substituting the determined order u and reference frequency ω0 into the frequency domain analytical expression of the seismic wavelet;

[0012] Preferably, the analytical expression of the Ricker wavelet is subjected to Fourier transformation, multiplied by a frequency factor, and then normalized to obtain the frequency domain analytical expression of the seismic wavelet.

[0013] The analytical expression of the Ricker wavelet is:

[0014]

[0015] Among them, τ is the time variable with the dimension of s, ω0 is the reference frequency with the dimension of πrad / s;

[0016] Preferably, the analytical expression of the Ricker wavelet is transformed by Fourier transform to obtain the following equation (2):

[0017]

[0018] Preferably, the analytical expression of the Ricker wavelet is subjected to Fourier transformation and multiplied by the frequency factor to obtain the following equation (3):

[0019]

[0020] In the formula, (iω) u is the frequency factor;

[0021] Preferably, the normalization refers to formula (3) multiplied by a normalization factor;

[0022] Preferably, the normalization factor is

[0023] Preferably, the frequency domain analytical expression of the seismic wavelet is:

[0024]

[0025] Where u is the order and dimensionless.

[0026] The order u is determined by the following formula (7):

[0027]

[0028] Preferably, the reference frequency ω0 is determined by the following formula (8):

[0029]

[0030] In formula (7) and formula (8), ω m is the mean frequency, the dimension is Hertz, Hz; ω σ is the frequency deviation, the dimension is Hertz, Hz;

[0031] Preferably, the mean frequency ω m and frequency deviation ω σ It is determined by fitting the power spectrum of the seismic record with a Gaussian function;

[0032] Preferably, the Gaussian function is

[0033] The reflection coefficient in step B is obtained by the following formula (15):

[0034]

[0035] Wherein, d is the seismic data in matrix form, dimensionless, W is the seismic wavelet in matrix form extracted in step A, dimensionless, and λ is the regularization parameter, dimensionless;

[0036] Preferably, the sparsity regularized inversion refers to sparsity regularized inversion based on the least squares method.

[0037] The reshaping process in step C refers to: performing frequency extension processing on the extracted seismic wavelet.

[0038] Step C further obtains a sufficiently low-frequency seismic wavelet by reshaping;

[0039] Preferably, step D convolves the sufficiently low-frequency seismic wavelet obtained in step C with the reflection coefficient in step B, so as to obtain seismic data after low-frequency information compensation;

[0040] Preferably, the reshaping process refers to performing frequency extension processing on the seismic wavelet extracted by equation (4).

[0041] The seismic data are desert seismic data.

[0042] The method for compensating low-frequency information of seismic data is used in analyzing seismic signals, oil exploration, and reservoir imaging.

[0043] The present invention provides a technical method for enhancing the low-frequency effective information and resolution of seismic data in desert areas, the method comprising:

[0044] First: Analytical seismic wavelet extraction technology

[0045] Since actual seismic data often cannot meet some basic assumptions of seismic wavelet extraction, and some methods have their own shortcomings, the seismic wavelets actually extracted are not smooth and have low accuracy, which makes it difficult to meet the requirements of seismic exploration technology. Therefore, it is necessary to combine the existing new methods and effective tools in the field of signal processing to further study the seismic wavelet extraction algorithm and improve the accuracy and stability of seismic wavelet estimation. In the conventional seismic wavelet extraction method, each sampling point of the seismic wavelet is essentially regarded as an unknown number, so the conventional seismic wavelet extraction method actually solves a lot of unknown numbers at the same time. Compared with solving a large number of unknown numbers, we can turn to solve the analytical seismic wavelet controlled by several key parameters. The analytical seismic wavelet has a clear expression formula. Therefore, by using the characteristics of the analytical seismic wavelet, if the key control parameters can be estimated from the actual seismic data, the corresponding seismic wavelet can be determined by the expression. Thereby, the efficiency and stability of seismic wavelet extraction can be improved.

[0046] Second: Instantaneous frequency-converting seismic wavelet reshaping technology:

[0047] Since the seismic record in the time domain can be regarded as the result of the convolution of the reflection coefficient and the seismic wavelet plus noise (i.e., the content expressed by formula (9), which is well known in the art), in the frequency domain, it is the result of the unequal weighted multiplication of the reflection coefficient spectrum by the seismic wavelet spectrum (similar to bandpass filtering). Therefore, the optimized generalized S transform or the continuous wavelet transform matching the seismic wavelet is used to map the time-frequency distribution of the seismic data, fit the seismic wavelet spectrum trend from the seismic record spectrum, highlight the reflection coefficient spectrum, and then use the instantaneous variable frequency seismic wavelet that can directly control the characteristic parameters such as the signal energy, bandwidth, and attenuation at the low frequency end to perform reshaping processing, enhance the energy of the low frequency component in the low frequency direction, and expand the effective frequency band range of the low frequency end, so as to obtain a seismic record with enhanced low frequency components and expanded low frequency effective frequency band.

[0048] Third: Sparse regularization inversion low-frequency compensation technology:

[0049] When the seismic wavelet is known, the least square method can be used to obtain the reflection coefficient sequence of the seismic data. Due to the sparsity of the seismic reflection coefficient, the L1 norm is used for constraint, and the least square sparse regularized inversion of the seismic data is performed to obtain the reflection coefficient sequence. Extracting the optimal seismic wavelet from the seismic data is also the key to the technology. We use the ASSW method proposed by Gao Jinghuai et al. (2017) to extract the optimal amplitude spectrum of the wavelet from the seismic amplitude spectrum and obtain the optimal seismic wavelet. While performing sparse inversion, the optimal parameters are optimized to obtain the optimal reflection coefficient sequence. Then, the instantaneous variable frequency seismic wavelet that can directly control the characteristic parameters such as low-frequency signal energy, bandwidth, and attenuation is used for reshaping processing, and the seismic record with enhanced low-frequency components and expanded low-frequency effective frequency band can be obtained.

[0050] The technical effects of the present invention are as follows:

[0051] (1) The low-frequency cutoff frequency and bandwidth of low-frequency compensation can be determined according to the characteristics of actual seismic data, and the high-frequency components can be expanded;

[0052] (2) It is used to enhance the energy of low-frequency components at the low-frequency end of seismic data, expand the effective frequency band range at the low-frequency end, and optimize processing parameters according to the noise level of seismic data and the original frequency band characteristics, thereby enhancing the low-frequency effective information and resolution of seismic data, significantly improving the imaging accuracy and resolution of geological bodies, and reducing the generation of reflection illusions.

[0053] The present invention can determine the low-frequency cutoff frequency and bandwidth of the low-frequency compensation according to the characteristics of the actual seismic data, and can expand the high-frequency components; it is used to enhance the energy of the low-frequency components at the low-frequency end of the seismic data, expand the effective frequency band range of the low-frequency end, and can optimize the processing parameters according to the noise level of the seismic data and the original frequency band characteristics, enhance the low-frequency effective information and resolution of the seismic data, significantly improve the imaging accuracy and resolution of the geological body, and reduce the generation of reflection illusions.

[0054] The present invention establishes a new method for low-frequency compensation and expansion of seismic data and a method for optimizing its processing parameters, which is used to enhance the energy of low-frequency components at the low-frequency end of the seismic data and expand the effective frequency band range of the low-frequency end. The processing parameters can be optimized according to the noise level of the seismic data and the original frequency band characteristics, thereby enhancing the low-frequency effective information and resolution of the seismic data as much as possible.

[0055] Aiming at the deep burial of Ordovician carbonate reservoirs in the desert area of ​​Tarim Basin, serious influence of complex surface layers and absorption attenuation on seismic signals, low signal-to-noise ratio and weak amplitude of seismic data, the present invention establishes a new low-frequency compensation and expansion method for seismic data and its processing parameters, which are used to enhance the energy of low-frequency components at the low-frequency end of seismic data, expand the effective frequency band range of the low-frequency end, improve the imaging accuracy and resolution of geological bodies, and significantly improve the imaging quality of fault-controlled reservoirs in medium and deep layers.

[0056] In Shunbei area, the Ordovician carbonate reservoir is deeply buried, and the seismic signal is seriously affected by the complex surface layer and absorption attenuation, so the signal-to-noise ratio of the seismic data is low and the amplitude is weak. In the process of seismic data processing, it is found that the low-frequency-weak amplitude response characteristics of deep and large faults are more sensitive to the low-frequency information in the seismic data. By making full use of the effective information at the low-frequency end of the seismic data, the imaging quality of the fault system can be significantly improved. The present invention is mainly used for the technology of compensating and expanding the low-frequency information of actual seismic data. According to the geological characteristics of Shunbei area, a numerical model of the seismic wave field under different low-frequency characteristics is established, and the influence of the low-frequency component of the seismic signal on the resolution of seismic data and the seismic wave field of different geological targets in deep strata is analyzed. Then, by extracting the original seismic data wavelet, a variable-frequency seismic wavelet analysis method is established, which can directly control the characteristic parameters such as the energy, bandwidth, and attenuation of the low-frequency end seismic signal. By selecting the optimal widened wavelet low cutoff frequency, the sparse regularization inversion method is adopted to realize the expansion of the low-frequency signal of the seismic data, thereby enhancing the low-frequency effective information of the seismic data and improving the resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] Figure 1 The generalized seismic wavelet extracted from seismic data in the experimental example of the present invention. (a) Original seismic data (the red arrow indicates the seismic trace used for seismic wavelet extraction); (b) Power spectrum of seismic data and its fitted Gaussian function, where the red curve is the amplitude spectrum of the frequency domain seismic wavelet shown in equation (4) of the present invention; the black curve is the amplitude spectrum of the seismic record instead of the power spectrum; (c) Extracted generalized seismic wavelet; (d) Comparison of the seismic trace (red) synthesized using the extracted seismic wavelet and the surrounding seismic traces.

[0058] Figure 2 The experimental example of the present invention is a time domain seismic wavelet and its spectrum before and after the seismic wavelet shaping process; wherein (a) is the time domain variable frequency seismic wavelet under different control parameters, and (b) is the seismic wavelet spectrum before and after the shaping (both the low frequency and the high frequency are expanded at the same time).

[0059] Figure 3 FIG. 1 is a processing flow chart of a method for compensating low-frequency information of seismic data provided by an embodiment of the present invention; the "seismic wavelet spectrum" obtained by S2 in the figure refers to the frequency domain seismic wavelet in equation (4), but the "amplitude spectrum of the seismic record" in S1 is as follows: Figure 1 The black curve in (b) is not a smooth curve. In this field, "spectrum" and "amplitude spectrum" are not the same concept. "Spectrum" is a complex number. After taking the modulus of "spectrum", we get "amplitude spectrum". For example, the modulus of complex number a+bi is

[0060] Figure 4This is an application example of the instantaneous frequency-converting seismic wavelet reshaping method. Left: before processing (noise suppression preprocessing has been performed), right: after processing; the stripes in the figure are called "even phase axes" in this field. The stripe density in the left figure is relatively sparse, and the stripe density in the right figure is denser. This difference in density shows that the seismic data has a higher resolution, which is the effect brought about by the broadening of the seismic wavelet spectrum.

[0061] Figure 5 The frequency spectrum of seismic records before and after processing by the method for compensating for low-frequency information of seismic data provided by an embodiment of the present invention is compared; wherein, the black curve is the amplitude spectrum of the seismic record before low-frequency compensation (i.e., the amplitude spectrum of the original seismic record), and the red curve is the amplitude spectrum of the seismic record after low-frequency compensation. The area of ​​the red curve after low-frequency compensation is larger than the area of ​​the black curve before low-frequency compensation, which is the result of the widening of the frequency band of the low-frequency information; the left figure and the right figure are actually different representations of the same content, the ordinate of the left figure is amplitude, and the ordinate of the right figure is dB, and the conversion method is: log10 (Amplitude) * 20. The black curve in the figure is the amplitude spectrum of the seismic record before low-frequency compensation, and the red curve is the amplitude spectrum of the seismic record after low-frequency compensation; it should be noted that, Figure 5 It is shown that the compensation method of the present invention not only expands the low-frequency information of the seismic data (the left side of the peak of the black curve belongs to the low-frequency area), but also expands the high-frequency information of the seismic data (the right side of the peak of the black curve belongs to the high-frequency area). DETAILED DESCRIPTION

[0062] The content of the present invention is further described in detail below in conjunction with specific embodiments, but the protection scope of the present invention is not limited thereto.

[0063] In this article, "seismic data" has the conventional technical meaning commonly understood by those skilled in the art, specifically referring to the vibration waveform curves recorded by the detectors placed on the surface of the earth as a wave signal generated by the vibration of the underground crust. The amplitude spectrum of seismic records, the power spectrum of seismic data, the spectrum and amplitude spectrum of seismic wavelets are all derived data generated by processing seismic data.

[0064] The “low-frequency information” of the seismic data refers to the region with a frequency lower than 10 Hz in the amplitude spectrum of the seismic record.

[0065] The first group of embodiments, the method for compensating low-frequency information of seismic data of the present invention

[0066] This group of embodiments provides a method for compensating low-frequency information of seismic data. All embodiments of this group have the following common features: the method for compensating low-frequency information of seismic data includes:

[0067] A. Extract seismic wavelets from seismic data;

[0068] B. Obtain the reflection coefficient of seismic data through sparsity regularized inversion;

[0069] C. reshaping the seismic wavelet extracted in step A to obtain a low-frequency seismic wavelet;

[0070] D. Convolve the low-frequency seismic wavelet obtained in step C with the reflection coefficient obtained by inversion in step B.

[0071] In a specific embodiment, extracting the seismic wavelet in step A means: obtaining a frequency domain analytical expression of the seismic wavelet by Ricker wavelet deduction, and substituting the determined order u and reference frequency ω0 into the frequency domain analytical expression of the seismic wavelet;

[0072] Preferably, the analytical expression of the Ricker wavelet is subjected to Fourier transformation, multiplied by a frequency factor, and then normalized to obtain the frequency domain analytical expression of the seismic wavelet.

[0073] The analytical expression of the Ricker wavelet is:

[0074]

[0075] Among them, τ is the time variable with the dimension of s, ω0 is the reference frequency with the dimension of πrad / s;

[0076] Preferably, the analytical expression of the Ricker wavelet is transformed by Fourier transform to obtain the following equation (2):

[0077]

[0078] Preferably, the analytical expression of the Ricker wavelet is subjected to Fourier transformation and multiplied by the frequency factor to obtain the following equation (3):

[0079]

[0080] In the formula, (iω) u is the frequency factor;

[0081] Preferably, the normalization refers to formula (3) multiplied by a normalization factor;

[0082] Preferably, the normalization factor is

[0083] Preferably, the frequency domain analytical expression of the seismic wavelet is:

[0084]

[0085] Where u is the order and dimensionless.

[0086] The order u is determined by the following formula (7):

[0087]

[0088] Preferably, the reference frequency ω0 is determined by the following formula (8):

[0089]

[0090] In formula (7) and formula (8), ω m is the mean frequency, the dimension is Hertz, Hz; ω σ is the frequency deviation, the dimension is Hertz, Hz;

[0091] Preferably, the mean frequency ω m and frequency deviation ω σ It is determined by fitting the power spectrum of the seismic record with a Gaussian function;

[0092] Preferably, the Gaussian function is

[0093] The reflection coefficient in step B is obtained by the following formula (15):

[0094]

[0095] Wherein, d is the seismic data in matrix form, dimensionless, W is the seismic wavelet in matrix form extracted in step A, dimensionless, and λ is the regularization parameter, dimensionless;

[0096] Preferably, the sparsity regularized inversion refers to sparsity regularized inversion based on the least squares method.

[0097] The reshaping process in step C refers to: performing frequency extension processing on the extracted seismic wavelet.

[0098] Step C further obtains a sufficiently low-frequency seismic wavelet by reshaping;

[0099] Preferably, step D convolves the sufficiently low-frequency seismic wavelet obtained in step C with the reflection coefficient in step B, so as to obtain seismic data after low-frequency information compensation;

[0100] Preferably, the reshaping process refers to performing frequency extension processing on the seismic wavelet extracted by equation (4).

[0101] The seismic data are desert seismic data.

[0102] The second group of embodiments, the use of the low-frequency information compensation method of seismic data of the present invention

[0103] This group of embodiments provides the use of a method for compensating low-frequency information of seismic data provided by any one of the first group of embodiments in analyzing seismic signals, oil exploration, and reservoir imaging.

[0104] Experimental Example 1: Specific operation steps of the method of the present invention

[0105] 1. Analytical seismic wavelet extraction technology

[0106] Analytical seismic wavelets have clear expression formulas. Therefore, by using the characteristics of analytical seismic wavelets, if the key control parameters can be estimated from the actual seismic data, the corresponding seismic wavelets can be determined by the expression. This can improve the efficiency and stability of seismic wavelet extraction. At present, several commonly used analytical wavelets are: Ricker wavelet, bandpass wavelet, damped cosine wavelet, controllable source wavelet, broadband Ricker wavelet (Yu wavelet), damped stretched sine wavelet, explosion seismic wavelet, etc. The most commonly used Ricker wavelet is theoretically a solution to the Stokes differential equation. Mathematically, the Ricker wavelet is a second-order derivative of a Gaussian function, and therefore, is symmetrical in the time domain. However, the seismic signals recorded at different depths in actual vertical boreholes are often asymmetric, and are closer to the first-order derivative or 1.5-order derivative of a Gaussian function. In order to better represent the actually observed seismic signals in non-Ricker wavelet form, Wang (2015) proposed the concept of generalized seismic wavelet, that is, modifying the order of the derivative to an arbitrary positive real number to obtain a generalized seismic wavelet. The Ricker wavelet is a special case when the order is 2.

[0107] Ricker (1943; 1944) pointed out that wavelet functions can be expressed as polynomials of various derivatives of potential functions. Therefore, referring to the analytical expression of Ricker wavelet, a Gaussian potential function can be set:

[0108]

[0109] Where τ is the time variable (unit: s), and ω0 is the reference frequency (unit: πrad / s). The reference angular frequency ω0 is inversely proportional to the deviation from the Gaussian distribution, which physically reflects the viscoelastic characteristics of the underground medium.

[0110] Performing Fourier transform on equation (1) yields:

[0111]

[0112] For any positive real-order derivative of a Gaussian function, its spectrum can be expressed as G(ω) multiplied by a frequency factor (iω) u :

[0113]

[0114] Where u is the order. Equation (3) is the unnormalized spectrum. To normalize the spectrum, you need to multiply it by a normalization factor Right now:

[0115]

[0116] Equation (4) is the frequency domain analytical formula of the generalized seismic wavelet. It can be seen that the generalized seismic wavelet is determined by three parameters u, ω0 and τ0. (u) By performing an inverse Fourier transform on (ω), we can transform the seismic wavelet in the frequency domain into the seismic wavelet Φ in the time domain. (u) (τ):

[0117]

[0118] In fact, the reference frequency ω0 and the order u have the most important influence on the seismic wavelet. The time symmetric position τ0 is usually set to 0 in actual use.

[0119] Wang (2015) pointed out that the statistical characteristics of the discrete Fourier spectrum of an actual seismic signal can be described by the mean frequency and frequency deviation. m and frequency deviation ω σ The analytical expression of can be written as:

[0120]

[0121] From the above formula, we can see that the mean frequency ω m and frequency deviation ω σ is a function of the reference frequency and order. If the mean frequency ω is known m and frequency deviation ω σ , then the order u and the reference frequency ω0 can be found. In addition, the mean frequency ω m and frequency deviation ω σ You can also define a Gaussian function:

[0122]

[0123] Therefore, based on this idea, the Gaussian function can be used to fit the power spectrum of the seismic record. The use of the Gaussian function to fit the power spectrum of the seismic record is a conventional operation in the art, and specific reference can be made to the fitting method described in the article "High-frequency recovery of seismic data constrained by well data".

[0124] The mean frequency ω is determined based on the result obtained by fitting the Gaussian function shown in the above equation (17): m and frequency deviation ω σ , and then determine the order u according to the following formula:

[0125]

[0126] In the above steps, the power spectrum of the seismic record is obtained by the following process: first perform Fourier transform on the seismic record, take the modulus of the Fourier transform result, and then square it. Mean frequency ω m and frequency deviation ω σ The process of determining these two parameters belongs to curve fitting, which is a known method in the art. Assuming that the power spectrum of the seismic record is x, the Gaussian function fitting is as follows: Conduct a series of m and ω σ A series of Gaussian curves can be constructed by using the value, and the Gaussian curve closest to the power spectrum x of the earthquake record can be found. The corresponding ω m and ω σ The value is the best result to be found. Using the power spectrum of the seismic record, a Gaussian function is constructed with different parameters for curve fitting. The set of parameters corresponding to the Gaussian function curve that is closest to the shape of the power spectrum of the seismic record is the mean frequency ω. m and frequency deviation ω σ Two parameters.

[0127] Once the order u is determined, the reference frequency ω0 can be determined as follows:

[0128]

[0129] After determining the order u and the reference frequency ω0, substitute them into equation (4) to extract the seismic wavelet.

[0130] Figure 1 The application example of extracting generalized seismic wavelets from actual seismic data is demonstrated. Figure 1 The seismic wavelet is extracted from the seismic trace pointed by the red arrow in a (time range: 2306-2505ms). The comparison between the amplitude spectrum of the seismic record and the amplitude spectrum of the extracted seismic wavelet is shown in Figure 1. Figure 1 b. Using the extracted seismic wavelet ( Figure 1 c) Synthetic seismic records such as Figure 1 As shown in the red seismic record in d, its correlation coefficient with the actual well bypass is 0.836.

[0131] 2. Instantaneous frequency conversion seismic wavelet reshaping technology

[0132] Since the seismic record in the time domain can be regarded as the result of the convolution of the reflection coefficient and the seismic wavelet plus noise, in the frequency domain it is the result of the unequally weighted multiplication of the reflection coefficient spectrum by the seismic wavelet spectrum (similar to bandpass filtering). Therefore, by fitting the seismic wavelet spectrum trend from the seismic record spectrum and highlighting the reflection coefficient spectrum, and then using the instantaneous variable frequency seismic wavelet that can directly control the characteristic parameters such as the low-frequency signal energy, bandwidth, and attenuation to perform reshaping processing, we can obtain seismic records that enhance the low-frequency components and expand the low-frequency effective frequency band.

[0133] In the above steps, the instantaneous frequency-converted seismic wavelet is the result of frequency extension processing of the seismic wavelet extracted by formula (4). "Reshaping processing" specifically refers to frequency extension processing of the seismic wavelet extracted by formula (4). Reshaping processing is a conventional method in this field, for example, the method described in Chinese patent 201910153993.9 "A seismic resolution enhancement technology based on wide-spectrum seismic wavelet shaping", or the method described in the document "Seismic Resolution Enhancement by Frequency Dependent Wavelet Scaling", (Shuangquan Chen, Yanghua Wang, 2018).

[0134] The amplitude spectrum of the seismic record after low-frequency compensation is obtained by multiplying the extracted reflection coefficient spectrum by the low-frequency seismic wavelet spectrum (i.e., the result obtained after the instantaneous frequency-converted seismic wavelet is subjected to the above-mentioned reshaping process);

[0135] By performing an inverse Fourier transform on the amplitude spectrum of the seismic record after low-frequency compensation, the seismic data after low-frequency compensation can be obtained. Figure 2 The results shown, Figure 2 (b) shows that the spectrum has been broadened. At this step, the low-frequency compensation purpose of the present invention has been basically achieved. However, the full effect of this low-frequency compensation still needs to be combined with step 3. The result of this step is reflected in the seismic data after the low-frequency compensation.

[0136] Figure 2 The comparison of wavelets after expanding different frequency band components is shown in Figure . Figure 2 The original seismic wavelet (black line) of a has a wide time duration and strong side lobes. Its spectrum ( Figure 2 b) has a narrow bandwidth and lacks low-frequency energy, but after the frequency band is expanded at the low and high frequency ends ( Figure 2 b) The corresponding time-domain seismic wavelet (blue line) has significantly weakened side lobes and a narrower time duration.

[0137] 3. Sparse regularization inversion low-frequency compensation technology

[0138] The basic process of inverting seismic reflection coefficients using the sparse regularization method is as follows:

[0139] Assume that the earthquake data is represented as:

[0140] d=w*r+n (9)

[0141]

[0142] Among them, d, w, r represent seismic records, seismic wavelets and reflection coefficients respectively, and n represents seismic noise.

[0143] The above-mentioned seismic record d is a component of the "seismic data" in this article. The "seismic data" is composed of a series of seismic records d. Generally speaking, seismic data is three-dimensional, and seismic records are one-dimensional; seismic record d is dimensionless.

[0144] The seismic wavelet w here refers to the seismic wavelet in the time domain, that is, the seismic wavelet Φ in the time domain expressed by equation (5) (u) (τ), rather than Φ in (4) (u) (ω) represents the seismic wavelet in the frequency domain.

[0145] The seismic wavelet w is used to construct the n×m dimensional Toeplitz structure matrix, which is as follows:

[0146]

[0147] The above seismic data can be expressed in matrix form as:

[0148] d=W·r+n (12)

[0149] The reflection coefficient is inverted as follows:

[0150]

[0151] Then we have:

[0152] r=(W Τ W) -1 W Τ d (14)

[0153] Considering the sparseness of the reflection coefficient, the above formula is sparsely constrained and the L1 norm regularization term is added to solve it:

[0154]

[0155] Among them, ||r||1 represents the L1 norm of r; d is the earthquake record, which is known; λ is the regularization parameter, which has no fixed calculation method. Generally, different parameters are tried and the parameter with the best effect is selected as the regularization parameter. The method for selecting and determining the regularization parameter is a conventional method in the field, for example, the regularization method described in the article "Regularization parameter selection for downward extension of potential field based on L-curve method" can be referred to.

[0156] Based on sparse regularized inversion, the reflection coefficient of seismic data is obtained, and then a sufficiently low-frequency seismic wavelet is further constructed and convolved with the obtained reflection coefficient to obtain low-frequency compensated seismic data.

[0157] To obtain the reflection coefficient, we must first fit the seismic wavelet spectrum trend from the seismic record spectrum, then use the method in step 1 to obtain the frequency domain seismic wavelet of formula (4), and after inverse Fourier transform, obtain the time domain seismic wavelet of formula (5), and finally use formula (15) to calculate the reflection coefficient; the frequency domain seismic wavelet obtained in this process, as shown in formula (4), is further used to construct a sufficiently low-frequency seismic wavelet. The above process needs to be performed for different seismic records.

[0158] "Sufficiently low frequency" in earthquakes means as low as 2-3 Hz. To achieve sufficiently low frequency is related to the selection of characteristic parameters such as low-frequency signal energy, bandwidth, and attenuation of the instantaneous variable frequency seismic wavelet.

[0159] Experimental Example 2

[0160] First, seismic wavelets are fitted from seismic data, and reflection coefficients are extracted using sparse regularization inversion. Then, parameter-controllable wide low-frequency analytical seismic wavelets that can directly control characteristic parameters such as low-frequency signal energy, bandwidth, and attenuation are used for reshaping. In this way, seismic records with enhanced low-frequency components and expanded low-frequency effective frequency bands can be obtained. The implementation process is as follows: Figure 3 shown.

[0161] Experimental Example 3

[0162] This experimental example uses the low-frequency information compensation method of seismic data provided by an embodiment of the present invention to perform low-frequency compensation on seismic data in a certain area. The results are as follows: Figure 5 shown. Figure 5 This is a comparison of data from a certain area before and after processing. It can be seen that the resolution of the seismic profile after processing is significantly improved. Compared with the spectrum of the seismic data before processing ( Figure 4 The black line in the figure shows the frequency spectrum of the processed seismic data ( Figure 5 (as shown by the red line in the figure) is significantly expanded at both the low and high frequency ends, and the frequency band is significantly widened.

Claims

1. A method for compensating low-frequency information of seismic data, characterized in that: include: A. Extract seismic wavelets from seismic data; B. Obtain the reflection coefficient of seismic data through sparsity regularized inversion; C. reshaping the seismic wavelet extracted in step A to obtain a low-frequency seismic wavelet; D. convolving the low-frequency seismic wavelet obtained in step C with the reflection coefficient obtained by inversion in step B; Extracting the seismic wavelet in step A means: obtaining the frequency domain analytical expression of the seismic wavelet by deducing the Ricker wavelet, and substituting the determined order u and reference frequency ω0 into the frequency domain analytical expression of the seismic wavelet; The analytical expression of the Ricker wavelet is: Among them, τ is the time variable with the dimension of s, ω0 is the reference frequency with the dimension of πrad / s; And / or, the analytical expression of Ricker wavelet is transformed by Fourier transform to obtain the following equation (2):

2. A method for compensating low-frequency information of seismic data according to claim 1, characterized in that: The frequency domain analytical expression of the seismic wavelet can be obtained by Fourier transforming the analytical expression of the Ricker wavelet, multiplying it by the frequency factor, and then normalizing it.

3. A method for compensating low-frequency information of seismic data according to claim 2, characterized in that: The analytical expression of Ricker wavelet is transformed by Fourier and multiplied by the frequency factor to obtain the following equation (3): In the formula, (iω) u is the frequency factor; And / or, the normalization refers to multiplying formula (3) by a normalization factor; And / or, the normalization factor is And / or, the frequency domain analytical expression of the seismic wavelet is: Where u is the order and dimensionless.

4. A method for compensating low-frequency information of seismic data according to claim 2 or 3, characterized in that: The order u is determined by the following formula (7): And / or, the reference frequency ω0 is determined by the following equation (8): In formula (7) and formula (8), ω m is the mean frequency, the dimension is Hertz, Hz; ω σ is the frequency deviation, the dimension is Hertz, Hz; and / or, the mean frequency ω m and frequency deviation ω σ It is determined by fitting the power spectrum of the seismic record with a Gaussian function; And / or, the Gaussian function is 5. A method for compensating low-frequency information of seismic data according to claim 1, characterized in that: The reflection coefficient in step B is obtained by the following formula (15): Wherein, d is the seismic data in matrix form, dimensionless, W is the seismic wavelet in matrix form extracted in step A, dimensionless, and λ is the regularization parameter, dimensionless; And / or, the sparsity regularized inversion refers to the sparsity regularized inversion based on the least squares method.

6. A method for compensating low-frequency information of seismic data according to claim 1, characterized in that: The reshaping process in step C refers to: performing frequency extension processing on the extracted seismic wavelet.

7. A method for compensating low-frequency information of seismic data according to claim 6, characterized in that: Step C obtains a sufficiently low-frequency seismic wavelet by reshaping; And / or, step D convolves the sufficiently low-frequency seismic wavelet obtained in step C with the reflection coefficient in step B to obtain seismic data after low-frequency information compensation.

8. A method for compensating low-frequency information of seismic data according to claim 7, characterized in that: The reshaping process refers to performing frequency extension processing on the seismic wavelet extracted by equation (4).

9. A method for compensating low-frequency information of seismic data according to any one of claims 1-3 and 5-8, characterized in that: The seismic data are desert seismic data.

10. Use of the method for compensating low-frequency information of seismic data as claimed in any one of claims 1 to 9 in analyzing seismic signals, oil exploration and reservoir imaging.

Images