A comprehensive method to suppress the strong shielding of igneous rocks and recover the effective weak signal energy of the underlying strata

Abstract

The present invention provides a comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata, comprising: step 1, establishing a corresponding forward model for the actual geological conditions of the igneous rock development area, quantitatively analyzing the seismic response characteristics of the igneous rocks and the factors affecting the underlying strata, and establishing an identification plate; step 2, according to the results of step 1, comparing and analyzing the changes in the main geophysical characteristics of the actual seismic data of the igneous rocks and their influence on the underlying strata; step 3, according to the results of step 2, conducting an analysis of the dominant frequency bands of the igneous rocks and their underlying strata, and optimizing the two-dimensional mixed phase graded wavelet in the dominant frequency band of the underlying target layer for time-frequency domain deconvolution, and performing energy recovery processing on the convolution results to obtain optimized seismic data. This method weakens the interference of igneous rocks, enhances the reflection characteristics of sand bodies, improves the identification accuracy of faults, and has high application value and good promotion prospects.

Description

Technical Field

[0001] The invention relates to the field of geophysical exploration and comprehensive research of clastic rock reservoirs in igneous rock development areas, and in particular to a comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of underlying strata. Background Art

[0002] Igneous rocks are widely developed in the exploration area of ​​Jiyang Depression in Shengli Oilfield: they have the characteristics of wide distribution layers and range, and diverse lithology. As a special lithology, igneous rocks are significantly different from sedimentary rocks. The reflection characteristics that appear on conventional seismic profiles are relatively special, and the amplitude and frequency seismic reflection characteristics vary greatly with different lithofacies. In the process of seismic excitation, igneous rocks are not only an energy shielding layer, but also cause distortion interference and a large amount of scattering to the seismic signal, resulting in the weakening of the energy of the down-transmitted waves and the complexity of the seismic wavelength, which causes the seismic reflection energy to weaken, the seismic signal-to-noise ratio to decrease, and it is difficult to effectively identify faults and clastic reservoirs, which seriously restricts the oil and gas exploration in the igneous rock development area.

[0003] At present, there are two main exploration challenges facing igneous rock development areas:

[0004] (1) Due to the abnormally high velocity and strong shielding effect of the overlying igneous rocks, the effective seismic signal of the underlying strata becomes weaker and the accuracy of fault identification is reduced. At the same time, due to the abnormally high velocity of the igneous rocks, the travel time of the seismic waves is reduced, and structural illusions appear in the time domain seismic data.

[0005] (2) In areas where igneous rocks are developed, igneous rock formations often have strong reflections because igneous rocks have strong shielding and absorption effects on seismic waves. This will weaken the seismic reflections of the underlying formations due to the strong shielding effects of igneous rocks, making them difficult to identify. How to extract and remove strong reflections on seismic data, eliminate the shielding effects of strong reflections of igneous rocks, and highlight the reflection information of effective reservoirs is a difficult point in reservoir prediction in areas where igneous rocks are developed.

[0006] The existing literature on igneous rocks mainly focuses on the research results of reservoirs, accumulation and lithofacies, but the impact on seismic data is less studied and most of them focus on single aspects such as velocity and signal-to-noise ratio. It is unable to provide effective guidance for the effective identification and interpretation of faults and clastic reservoirs in igneous rock development areas.

[0007] In the Chinese patent application with application number: CN201910271274.7, a method for amplitude compensation of underlying strata based on quantitative analysis of volcanic rock shielding is involved. First, the wells that have been drilled are selected; second, the time-depth calibration of the wells is selected; third, synthetic seismic records of seismic data of different frequency bands are produced; fourth, the reflection amplitude ratio of volcanic rock and underlying strata in seismic data of different frequency bands and synthetic records is calculated; fifth, the frequency band with a ratio close to 1 is selected as the reference frequency band for amplitude compensation of subsequent seismic data; sixth, the volcanic rock in the actual data is calculated. stratum and underlying strata; seven, decompose the seismic data into different frequency subsets in the time-frequency domain; eight, perform weighted amplitude compensation on seismic data in different frequency bands to make the amplitude ratio of the overlying volcanic strata and the underlying strata in each frequency band consistent; nine, complete the time-frequency domain frequency division amplitude compensation; ten, compare the compensation energy changes of seismic data in different frequency bands and the full frequency band and conduct quality control analysis, and repeat four to nine when the compensation effect is not ideal; simple and efficient amplitude preservation eliminates the problem of weakening the reflected energy of the underlying strata caused by energy shielding of the overlying high-speed volcanic rocks.

[0008] In the Chinese patent application with application number: CN201610798160.4, a method for eliminating the influence of coal-bearing strata and restoring the seismic reflection amplitude of the underlying strata is involved. This method studies the amplitude shielding effect of coal-bearing strata of different thicknesses on the sand bodies in the underlying strata of the same set of strata, and combines drilling data with seismic profile data to obtain the shielding amount of the unit coal-bearing strata in the study area on the underlying target strata - the shielding influence rate of the coal-bearing strata, so as to obtain the shielding amount of different coal-bearing strata on the amplitude value of the same target stratum, and the shielding amount plus the apparent amplitude value of the existing target stratum, that is, the true amplitude value of the target stratum is obtained.

[0009] In the Chinese patent application with application number: CN201811155579.3, a method for depth domain velocity modeling of igneous rock is involved, including: obtaining the initial velocity field and logging stratification data of the igneous rock; inverting the background velocity field based on the initial velocity field and logging stratification data; based on the background velocity field, improving the resolution and performing network tomography inversion iteration to obtain the initial igneous rock depth domain velocity model; adding well seismic errors and structural errors to the objective function of the initial igneous rock depth domain velocity model to obtain the final igneous rock depth domain velocity model.

[0010] The above existing technologies are greatly different from the present invention and fail to solve the technical problem we want to solve. For this reason, we have invented a new comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata. Summary of the invention

[0011] The purpose of the present invention is to provide a new comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of underlying strata by reconstructing multi-information fusion after seismic signal decomposition based on matching pursuit (MP) algorithm to improve the recovery of effective signal energy in igneous rock development areas.

[0012] The purpose of the present invention can be achieved by the following technical measures: a comprehensive method for suppressing the strong shielding of igneous rocks to restore the effective weak signal energy of the underlying strata, the comprehensive method for suppressing the strong shielding of igneous rocks to restore the effective weak signal energy of the underlying strata comprises:

[0013] Step 1: Establish a corresponding forward model based on the actual geological conditions of the igneous rock development area, quantitatively analyze the seismic response characteristics of igneous rocks and the factors affecting the underlying strata, and establish an identification plate;

[0014] Step 2: Based on the results of step 1, compare and analyze the changes in the main geophysical characteristics of the actual seismic data of the igneous rock and its underlying strata;

[0015] Step 3, based on the results of step 2, analyze the dominant frequency bands of igneous rocks and their underlying strata, select a two-dimensional mixed phase-graded wavelet for time-frequency domain deconvolution within the dominant frequency band of the underlying target layer, and perform energy recovery processing on the convolution result to obtain optimized seismic data.

[0016] The purpose of the present invention can also be achieved by the following technical measures:

[0017] In step 1, forward modeling quantitatively analyzes the seismic response characteristics of igneous rocks including frequency, velocity, energy and attenuation.

[0018] In step 1, the impacts on the underlying strata include: the impact of different velocities between igneous rocks and surrounding rocks on the underlying surrounding rocks and faults; the impact of thickness changes of a single set of igneous rocks with the same velocity on the underlying surrounding rocks and faults; the impact of the occurrence of igneous rocks, including the impact of different dips on the underlying surrounding rocks and faults; and the impact of the distance between igneous rocks and underlying strata on the underlying surrounding rocks and faults.

[0019] In step 2, the seismic geophysical characteristics of igneous rocks and their underlying strata include changes in kinematic, dynamic and statistical seismic characteristic quantities such as frequency, energy and waveform, which provide a basis for the strong shielding stripping of igneous rocks and the recovery of weak signal energy.

[0020] In step 3, based on the results of step 2, the dominant frequency bands of igneous rocks and their underlying strata are analyzed, and a two-dimensional mixed phase-graded wavelet is selected to perform time-varying frequency deconvolution within the dominant frequency band of the underlying target layer. The convolution results are subjected to energy recovery processing such as dip-guided filtering, amplitude equalization and AGC to obtain optimized seismic data.

[0021] In step 3, based on the results of step 2, different dominant frequency bands with better energy, amplitude and waveform morphology of igneous rock and underlying target layer are preferably distinguished; frequency band optimization is implemented based on matching pursuit MP algorithm; the dominant frequency band of the underlying target layer is selected for time-frequency domain deconvolution, and the result is subjected to amplitude equalization and AGC and other energy recovery processing, and finally an optimized seismic data body is obtained; the core idea of ​​the matching pursuit algorithm is to represent the signal as a linear combination of a series of time-frequency atoms that best match the local structural characteristics of the signal; the essence of finding the best matching time-frequency atom is to calculate the local instantaneous properties of the seismic signal, and the superposition of the time-frequency distribution of each time-frequency atom is the matching pursuit time-frequency.

[0022] Step 3 includes:

[0023] Step 3a, performing sparse expression of seismic signals;

[0024] Step 3b, obtaining the time-frequency spectra of the decomposed different atomic signals;

[0025] Step 3c, extracting multi-well mixed phase wavelets;

[0026] Step 3d, using the optimal wavelet to perform deconvolution, the convolution result is subjected to amplitude equalization, dip-guided filtering and AGC amplitude gain and other energy recovery and enhancement processing, that is, the final seismic data after removing the strong shielding effect of igneous rocks is obtained.

[0027] In step 3a, the sparse expression of the seismic signal is: Where s(t) is the band-limited seismic signal, Rs (m) (t) is the residual signal after matching, m is the number of atoms in iteration, and the termination condition of iteration is judged by the set threshold, thereby determining the number of atoms in iteration; atoms are determined by four parameters t j Delay, f j Main frequency, a j Correlation coefficient, φ j Phase controlled;

[0028] Given an atomic library D = {w k ,k=1,2,...K}, and satisfy ||w k ||=1, these vectors can be expanded into N-dimensional Hilbert space; for any signal s(t), find the atom w that best matches it from D kr ; The following matching conditions are met: Among them, <…> represents the inner product, and Rs (1) (t) represents the difference between the projection value and the signal, that is, the residual signal of the first iteration, so the signal can be expressed as: s(t) = <s,w k0 >w k0 +Rs (1) (t);

[0029] In the initial state, Rs (0) (t) = s(t); obviously w k0 and Rs. (1) (t) is orthogonal, we can get: ||Rs (0) (t)|| 2 =|〈s,w k0 >| 2 +||Rs (1) (t)|| 2 , successive iterations approach the optimization, which must minimize the energy of the residual Rs(t) and make the projection |〈s,w k0 >|Maximization; after n iterations, Rs (n) (t) = a n w n +Rs (n+1) (t), The signal is finally broken down as shown below:

[0030] In step 3b, the analytical expression of the time-frequency spectrum is obtained according to Liu Jianlei:

[0031] j W R (f,f j )env[w R (tt j ,f j )]

[0032] Where: a j is the wavelet factor, is the frequency distribution, is the instantaneous envelope property in the frequency domain, is the phase change under Fourier transform; tt j is the time difference, f j is the frequency, For phase.

[0033] ① Perform Hilbert transformation on s(t) to construct a complex analysis signal: S(t) = s(t) + jHT(s(t)), and determine the instantaneous envelope and instantaneous phase formula of the signal as follows:

[0034] Determine the delay t according to the corresponding time of the instantaneous envelope extreme value and the peak and valley position j , and find the instantaneous phase θ(t j); where S(t) is the complex analysis signal of the residual signal, s(t) is the initial residual signal. jHT(s(t)) is the Hilbert change of the initial signal, A(t) is the instantaneous envelope property, and θ(t) is the instantaneous phase property.

[0035] ② Obtain the instantaneous frequency of the signal, construct a matching sub-wave near the delay, and select the sub-wave W (t j ,f j ,θ j ) Its corresponding main frequency is the optimal main frequency of the matching sub-wave; W(t j ,f j ,θ j ) is the optimal sub-wave under different time, frequency and phase.

[0036] ③According to the correlation, the factors of the wavelet are obtained under the condition of normalization of the wavelet energy. The matching atom with the largest correlation with the signal is subtracted from the seismic trace to obtain the residual signal, which is then iterated recursively within a certain threshold range; thus, the time-frequency spectra of all matching sub-waves are superimposed to obtain the time-frequency spectra of the initial signal.

[0037] In step 3c, ignoring the noise term, the relationship between the seismic record and the seismic wavelet and reflection coefficient sequence can be expressed by the convolution model: in the frequency domain, it can be expressed as: x(ω)=w(ω)ξ(ω); where x(ω), w(ω), ξ(ω) are the seismic record, seismic wavelet and reflection coefficient in the frequency domain, respectively. After taking the logarithm, it can be transformed into a linear system lnx(ω)=lnw(ω)+lnξ(ω), and then the anti-rich transformation can be obtained: They are called the complex spectrum sequences of x(t), w(t), ξ(t) respectively; in the complex spectrum domain, the minimum and maximum phase components of the wavelet can be determined, and then a set of wavelets with different phase spectra can be determined by scanning;

[0038] Based on this, multiple wells are selected to scan the wellside seismic traces in a small window in the reciprocal spectrum domain, and the optimal wavelet is selected and calculated by comparing the synthetic records with the actual time-frequency volume amplitude spectrum.

[0039] The comprehensive method of suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata in the present invention is to establish a seismic identification panel of the clastic reservoir underlying the igneous rocks on the basis of quantitatively analyzing the seismic response characteristics of the igneous rocks and the factors affecting the surrounding rocks by forward simulation; combine the original seismic data to perform effective signal spectrum analysis, use the wellside seismic trace to optimize the two-dimensional mixed phase graded wavelet in the time-frequency domain, and perform time-frequency domain deconvolution reconstruction in combination with the local dip angle information of the target layer, suppress the strong interference of the igneous rocks, and recover the energy of the weak signals of the underlying strata, thereby improving the effective reflection signals of the underlying strata in the igneous rock development area, forming a complete set of comprehensive methods for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata. This method weakens the interference effect of igneous rocks, enhances the reflection characteristics of sand bodies, effectively improves the quality of seismic data in the igneous rock development area, improves the recognition accuracy of faults, and provides a reliable basis for the description of effective reservoirs and faults in the igneous rock development area, which has high application value and good promotion prospects.

[0040] The present invention is proposed to solve the problem of chaotic seismic reflection signals and low identification accuracy of faults and clastic reservoirs in igneous rock development areas. It provides a comprehensive method for suppressing the strong shielding of igneous rocks and restoring the effective weak signal energy of the underlying strata, effectively improving the quality of seismic data, the identification accuracy of faults and effective reservoirs, and guiding exploration and development.

[0041] Aiming at the characteristics of chaotic seismic reflection characteristics and strong shielding of seismic reflection signals of underlying strata in igneous rock development areas, this invention proposes a method of suppressing igneous rock shielding in the time-frequency domain and enhancing the seismic reflection characteristics of underlying strata by mixed phase deconvolution. Through the fusion of multiple information from seismic and drilling, the quality of seismic data in igneous rock development areas is effectively improved. Compared with the existing technology, it fully utilizes various information of seismic data and reasonable geophysical means to improve the identification accuracy of faults and reservoirs under the influence of igneous rock shielding, and has good application effect and promotion prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 A flowchart of a specific embodiment of the comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying stratum of the present invention;

[0043] Figure 2 It is a schematic diagram of the geological forward model and seismic response of the underlying strata affected by igneous rocks in the specific embodiment 1 of the present invention;

[0044] Figure 3 A comparison diagram of geophysical characteristics of different facies of igneous rocks in specific embodiment 1 of the present invention;

[0045] Figure 4 It is a schematic diagram of the identification plate of the seismic response characteristics of the underlying strata affected by igneous rocks in the specific embodiment 1 of the present invention;

[0046] Figure 5 A comparison diagram of the geophysical characteristics of igneous rock and surrounding rock in the specific embodiment 1 of the present invention;

[0047] Figure 6 It is a schematic diagram of decomposition of actual seismic data signals in specific embodiments 1 and 3 of the present invention;

[0048] Figure 7 It is a flow chart of time-frequency spectrum calculation by the matching pursuit algorithm (MP) in specific embodiments 1 and 3 of the present invention and a time-frequency spectrum diagram of an actual seismic signal;

[0049] Figure 8 It is a diagram of the multi-well mixed phase wavelet calculation process in specific embodiments 1 and 3 of the present invention;

[0050] Fig. 9 It is a comparison diagram of the original seismic profile, the intermediate profile of the signal recovery process and the final profile in the specific embodiment 1 of the present invention;

[0051] Fig.10 It is the reservoir plane prediction map (seismic data volume after signal energy recovery processing) in the specific embodiment 1 of the present invention;

[0052] Fig.11 It is a comparison diagram of the coherence and dip superposition plane of the original seismic and multi-body fusion attribute volumes after processing in the specific embodiment 2 of the present invention;

[0053] Fig.12 It is a comparison diagram between the original seismic profile and the processed and reconstructed profile in specific embodiment 3 of the present invention. DETAILED DESCRIPTION

[0054] It should be noted that the following detailed descriptions are exemplary and are intended to provide further explanation of the present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present invention belongs.

[0055] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit exemplary embodiments according to the present invention. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprising" and / or "including" are used in this specification, it indicates the presence of features, steps, operations and / or combinations thereof.

[0056] The comprehensive method of suppressing igneous rock strong shielding and recovering the effective weak signal energy of the underlying stratum of the present invention comprises the following steps:

[0057] Step 1: Establish a corresponding forward model based on the actual geological conditions of the igneous rock development area, quantitatively analyze the seismic response characteristics of igneous rocks and the characteristics of the factors affecting the underlying strata, and establish an identification plate;

[0058] The forward simulation quantitative analysis of the seismic response characteristics of igneous rocks mainly includes: frequency, velocity, energy and attenuation. The impact on the underlying strata mainly includes: the different velocities of igneous rocks and surrounding rocks affect the underlying surrounding rocks and faults; the thickness change of a single set of igneous rocks with the same velocity affects the underlying surrounding rocks and faults; the occurrence of igneous rocks mainly affects the underlying surrounding rocks and faults due to different dip angles; the distance between igneous rocks and underlying strata affects the underlying surrounding rocks and faults.

[0059] Step 2: Based on the results of step 1, compare and analyze the changes in the main geophysical characteristics of the actual seismic data of the igneous rock and its underlying strata;

[0060] The seismic geophysical characteristics of igneous rocks and their underlying strata mainly include changes in kinematic, dynamic and statistical seismic characteristic quantities such as frequency, energy and waveform, which provide a basis for the strong shielding stripping of igneous rocks and the recovery of weak signal energy.

[0061] Step 3: According to the result of step 2, the dominant frequency band of igneous rock and its underlying strata is analyzed, and a two-dimensional mixed phase graded wavelet is selected to perform time-varying frequency deconvolution in the dominant frequency band of the underlying target layer. The convolution result is subjected to a series of energy recovery processing such as dip-guided filtering, amplitude equalization and AGC to obtain optimized seismic data.

[0062] According to the results of step 2, it is preferred to distinguish different dominant frequency bands with better energy, amplitude and waveform morphology of igneous rocks and underlying target layers. Frequency band optimization is implemented based on the matching pursuit (MP) algorithm. The dominant frequency band of the underlying target layer is selected for time-frequency domain deconvolution, and the result is subjected to a series of energy recovery processing such as amplitude equalization and AGC, and finally the optimized seismic data body is obtained. The core idea of ​​the matching pursuit algorithm: the signal is represented as a linear combination of a series of time-frequency atoms that best match the local structural characteristics of the signal. The essence of finding the best matching time-frequency atom: calculating the local instantaneous properties of the seismic signal, the superposition of the time-frequency distribution of each time-frequency atom is the matching pursuit time-frequency.

[0063] In general, seismic wave s(t) is regarded as the convolution of sub-wave ω(t) and reflection coefficient sequence γ(t), plus noise n(t), that is, s(t) = ω(t)*γ(t)+n(t). Ignoring the influence of noise, as long as the best atom of the target layer is found in the time-frequency domain, the strong shielding effect of igneous rock can be suppressed to the maximum extent and the signal energy of the target layer can be restored. The whole process is realized by time-frequency analysis based on MP algorithm.

[0064] Given an atomic library D = {g r,r=1,2,…R}, and satisfy ||g r ||=1, these vectors can be expanded into N-dimensional Hilbert space. For any vector signal f in the space, its best matching atom can be found from D The signal f is first decomposed into: Here Rf is The residual after decomposition in the direction is obviously Orthogonal to Rf, so In order to minimize ||Rf||, we must find make The energy is the largest. Under this condition, we only need to find the best atomic vector in the atomic library. make Here α is the optimization factor and satisfies 0<α≤1.

[0065] Let R 0 f = f, assuming that the residual R of the nth time is known n f, n ≥ 0. There is a good match for the remainder R n Elements of f Then the residual R n f can be decomposed into This formula defines the residual of the n+1th order. n+1 f Orthogonal: Continue this decomposition until the mth time, decomposing f into a cascade summation form From formula 1, we can get The matching pursuit (MP) algorithm is convergent. As the decomposition proceeds, the error energy gradually decreases. If the atom library is complete enough and the number of iterations is unlimited, the linear combination of the decomposed atoms can be close enough to the original signal. The frequency band optimization is completed by obtaining the time-frequency spectrum comparison analysis for each time-frequency atom after decomposition.

[0066] The specific implementation process is as follows:

[0067] (1) The sparse expression of seismic signal is: Where s(t) is the band-limited seismic signal, Rs (m) (t) is the residual signal after matching, m is the number of atoms in iteration, and the termination condition of iteration is judged by setting the threshold, thereby determining the number of atoms in iteration. Atoms are determined by four parameters t j (delay), f j (main frequency), a j (correlation coefficient), φ j (phase) to control. Given the atomic library D = {w k ,k=1,2,...K}, and satisfy ||wk ||=1, these vectors can be expanded into N-dimensional Hilbert space. For any signal s(t), find the atom w that best matches it from D. kr . The following matching conditions are met: Among them, <…> represents the inner product, and Rs (1) (t) represents the difference between the projection value and the signal, that is, the residual signal of the first iteration, so the signal can be expressed as: s(t) = <s,w k0 >w k0 +Rs (1) (t). In the initial state, Rs (0) (t) = s(t). Obviously w k0 and Rs. (1) (t) is orthogonal, we can get: ||Rs (0) (t)|| 2 =| <s,w k0 >| 2 +||Rs (1) (t)|| 2 , successive iterations approach the optimization, which must minimize the energy of the residual Rs(t) and make the projection | <s,w k0 >|Maximize. After n iterations, Rs (n) (t) = a n w n +Rs (n+1) (t), The signal is finally broken down as shown below:

[0068] (2) Obtain the time-frequency spectrum of different decomposed atomic signals: According to the analytical expression for obtaining the time-frequency spectrum proposed by Liu Jianlei:

[0069] W R (f,f j )env[w R (tt j ,f j )]

[0070] (Where: a j is the wavelet factor, is the frequency distribution, is the instantaneous envelope property in the frequency domain, (Phase change under Fourier transform)① Perform Hilbert transformation on s(t) to construct a complex analysis signal: S(t) = s(t) + jHT(s(t)), and determine the instantaneous envelope and instantaneous phase formula of the signal as follows: Determine the delay t according to the corresponding time of the instantaneous envelope extreme value and the peak and valley position j , and find the instantaneous phase θ(tj ). ② Obtain the instantaneous frequency of the signal, construct a matching sub-wave near the delay, and select the sub-wave W (t j ,f j ,θ j ) Its corresponding main frequency is the optimal main frequency of the matching sub-wave. ③According to the correlation, the factor of the sub-wave is obtained under the condition of normalization of the sub-wave energy. The matching atom with the largest correlation with the signal is subtracted from the seismic trace to obtain the residual signal, which is iterated recursively within a certain threshold range. In this way, the time-frequency spectrum of all matching sub-waves is superimposed to obtain the time-frequency spectrum of the initial signal.

[0071] (3) Extraction of multi-well mixed phase wavelets: The above two steps can obtain a time-frequency volume that highlights different signal energies. By comparing and analyzing the geophysical characteristics of the actual seismic data of igneous rocks and their underlying strata statistically analyzed in step 2, a time-frequency volume that eliminates the shielding effect of igneous rocks on the underlying strata signal energy can be obtained. However, there is still a problem here: the wavelets in the atomic library during signal decomposition and reconstruction are generally standard wavelets with zero phase, while changes in the thickness and occurrence of strata in actual earthquakes will cause changes in the wavelet phase, energy and frequency. In order to obtain a time-frequency volume that can more accurately reflect the actual geological conditions, it is necessary to further optimize the wavelets. Ulrych (1971) introduced deconvolution technology into seismic exploration. Ignoring the noise term, the relationship between the seismic record and the seismic wavelet and reflection coefficient sequence can be expressed by the convolution model: In the frequency domain, it is expressed as: x(ω) = w(ω)ξ(ω). Where x(ω), w(ω), ξ(ω) are the seismic record, seismic wavelet and reflection coefficient in the frequency domain, respectively. After taking the logarithm, it can be transformed into a linear system lnx(ω)=lnw(ω)+lnξ(ω), and then an anti-rich transformation can be obtained: They are called the reciprocal spectrum sequences x(t), w(t), ξ(t) respectively. In the reciprocal spectrum domain, we can determine the minimum and maximum phase components of the wavelet, and then determine a set of wavelets with different phase spectra by scanning. Based on this, we can select multiple wells to scan the wellside seismic traces in the reciprocal spectrum domain with a small window, and select and calculate the optimal wavelet by comparing the synthetic record with the actual time-frequency volume amplitude spectrum.

[0072] (4) Deconvolution is performed using the optimal wavelet, and a series of energy recovery and enhancement processes such as amplitude equalization, dip-guided filtering, and AGC amplitude gain are performed on the convolution results to obtain the final seismic data after removing the strong shielding effect of igneous rocks.

[0073] Embodiment 1:

[0074] In the specific embodiment 1 of the present invention, Figure 1As shown: This embodiment provides a comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata, including the following steps:

[0075] Step 1: Establish a corresponding forward model based on the actual geological conditions of the igneous rock development area, quantitatively analyze the seismic response characteristics of igneous rocks and the characteristics of the factors affecting the underlying strata, and establish an identification plate;

[0076] In areas where igneous rocks are developed, strong reflections often occur in igneous rock strata because igneous rocks have extremely strong shielding and absorption effects on seismic waves. Since the energy of these strong reflection axes is too strong, the weaker reflections near the strong reflection axes will be shielded, making them difficult to identify. According to the actual geological conditions, four types of igneous rock geological models are designed. Based on the models, the seismic response characteristics of igneous rocks and the factors affecting the underlying strata are analyzed, and a quantitative identification template for the influencing factors is established.

[0077] ① A single-layer volcanic rock geological model was designed, with the igneous rock velocity unchanged and the thickness increasing from 5m to 60m, with a thickness interval of 5m;

[0078] ② An igneous rock model with the same thickness was designed, with the speed increased from 4000m / s to 6000m / s, with a speed interval of 200m / s;

[0079] ③ An igneous rock model with the same thickness and velocity was designed, and its inclination angle with the underlying strata increased from 5 degrees to 50 degrees; the inclination interval was 5 degrees.

[0080] ④ An igneous rock model with the same thickness and velocity was designed, and the distance between it and the underlying strata increased from 10m to 100m, with a distance interval of 10m; ( Figure 2 )

[0081] Based on the established geological model, the seismic response characteristics of igneous rocks and the quantitative identification of the seismic response characteristics of the underlying strata are summarized; Figure 3 , Figure 4 )

[0082] Step 2: Based on the results of the forward model, compare and analyze the changes in the main geophysical characteristics of the actual seismic data of the igneous rock and its underlying strata: mainly including the changes in kinematic and statistical seismic characteristic quantities such as velocity, frequency, energy and waveform, to provide a basis for the strong shielding stripping of igneous rocks and the recovery of weak signal energy. Figure 5 )

[0083] Step 3, based on the analysis results of step 2, carry out analysis of the dominant frequency bands of igneous rocks and their underlying strata, select a multi-dimensional mixed phase-graded wavelet for time-varying frequency deconvolution within the dominant frequency band of the underlying target layer, and perform a series of energy recovery processing such as dip-guided filtering, amplitude equalization and AGC on the convolution results to obtain optimized seismic data.

[0084] The specific implementation process in step 3 is as follows:

[0085] (1) Actual seismic data signal decomposition

[0086] The sparse expression of the seismic signal is: Where s(t) is the band-limited seismic signal, Rs (m) (t) is the residual signal after matching, m is the number of atoms in iteration, and the termination condition of iteration is judged by setting the threshold, thereby determining the number of atoms in iteration. Atoms are determined by four parameters t j (delay), f j (main frequency), a j (correlation coefficient), φ j (phase) to control. Given the atomic library D = {w k ,k=1,2,...K}, and satisfy ||w k ||=1, these vectors can be expanded into N-dimensional Hilbert space. For any signal s(t), find the atom w that best matches it from D. kr . The following matching conditions are met: Among them, <…> represents the inner product, and Rs (1) (t) represents the difference between the projection value and the signal, that is, the residual signal of the first iteration, so the signal can be expressed as: s(t) = <s,w k0 >w k0 +Rs (1) (t). In the initial state, Rs (0) (t) = s(t). Obviously w k0 and Rs. (1) (t) is orthogonal, we can get: ||Rs (0) (t)|| 2 =| <s,w k0 >| 2 +||Rs (1) (t)|| 2 , successive iterations approach the optimization, which must minimize the energy of the residual Rs(t) and make the projection | <s,w k0 >|Maximize. After n iterations, Rs (n) (t) = a n w n +Rs (n+1) (t), The signal is finally broken down as shown below: ( Figure 6 )

[0087] (2) Obtain the time-frequency spectrum of a series of atomic signals obtained by decomposition

[0088] Obtain the time-frequency spectrum of different decomposed atomic signals: According to the analytical expression for obtaining the time-frequency spectrum proposed by Liu Jianlei: (Where: a j is the wavelet factor, W R (f,f j ) is the frequency distribution, env[w R (tt j ,f j )] is the instantaneous envelope property in the frequency domain, (Phase change under Fourier transform)① Perform Hilbert transformation on s(t) to construct a complex analysis signal: S(t) = s(t) + jHT(s(t)), and determine the instantaneous envelope and instantaneous phase formula of the signal as follows: Determine the delay t according to the corresponding time of the instantaneous envelope extreme value and the peak and valley position j , and find the instantaneous phase θ(t j ). ② Obtain the instantaneous frequency of the signal, construct a matching sub-wave near the delay, and select the sub-wave W (t j ,f j ,θ j ) Its corresponding main frequency is the optimal main frequency of the matching sub-wave. ③According to the correlation, the factor of the sub-wave is obtained under the condition of normalization of the sub-wave energy. The matching atom with the largest correlation with the signal is subtracted from the seismic trace to obtain the residual signal, and the residual signal is iterated recursively within a certain threshold range. In this way, the time-frequency spectrum of all matching sub-waves is superimposed to obtain the time-frequency spectrum of the initial signal. ( Figure 7 )

[0089] (3) Calculation of multi-well mixed phase wavelets: The above two steps can obtain a time-frequency volume that highlights different signal energies. By comparing and analyzing the geophysical characteristics of the actual seismic data of igneous rocks and their underlying strata statistically analyzed in step 2, a time-frequency volume that eliminates the shielding effect of igneous rocks and reflects the signal energy of the underlying strata can be obtained. However, there is still a problem here: the wavelets in the atomic library during signal decomposition and reconstruction are generally standard wavelets with zero phase, while changes in the thickness and occurrence of strata in actual earthquakes will cause changes in the phase, energy and frequency of the wavelets. In order to obtain a time-frequency volume that can more accurately reflect the actual geological conditions, it is necessary to further optimize the wavelets. Ulrych (1971) introduced deconvolution technology into seismic exploration. Ignoring the noise term, the relationship between the seismic record and the seismic wavelet and reflection coefficient sequence can be expressed by the convolution model: In the frequency domain, it is expressed as: x(ω) = w(ω)ξ(ω). Where x(ω), w(ω), ξ(ω) are the seismic record, seismic wavelet and reflection coefficient in the frequency domain, respectively. After taking the logarithm, it can be transformed into a linear system lnx(ω)=lnw(ω)+lnξ(ω), and then an anti-rich transformation can be obtained: They are called the complex spectrum sequence of x(t), w(t), ξ(t) respectively. In the complex spectrum domain, we can determine the minimum and maximum phase components of the wavelet, and then determine a set of wavelets with different phase spectra by scanning. Based on this: we can select multiple wells to scan the wellside seismic traces in the complex spectrum domain with a small window, and select and calculate the optimal wavelet by comparing the synthetic record with the actual time-frequency volume amplitude spectrum. ( Figure 8 )

[0090] (4) Deconvolution of seismic data in the time-frequency domain: Deconvolution is performed using the optimal wavelet, and a series of energy recovery and enhancement processes such as amplitude equalization, dip-guided filtering, and AGC amplitude gain are performed on the convolution results to obtain the final seismic data after removing the strong shielding effect of igneous rocks. Fig. 9 )

[0091] After removing the strong shielding effect of igneous rocks, the amplitude attributes of the seismic data volume after the signal energy of the underlying target layer is restored can well reflect the distribution range of the reservoir in the igneous rock development area, and has a good application effect. Fig.10 )

[0092] Embodiment 2:

[0093] In this embodiment, the implementation process of step 1 and step 2 is the same as that of embodiment 1, and the only difference is step 3, which is specifically described as follows:

[0094] Step 3, perform wavelet time-frequency decomposition on the original seismic data, compare and analyze the dominant frequency band that can more clearly reflect the underlying stratum fault under the influence of strong shielding of igneous rocks based on the results of step 2, select the wavelet time-frequency decomposition body for fault boundary enhancement processing such as variance and ant tracking, perform multi-body fusion, and finally obtain a seismic data body that can better reflect fault changes.

[0095] The specific implementation process is as follows:

[0096] (1) Perform multi-wavelet seismic trace decomposition on the original data. The specific principles and steps are the same as those in step 3 (1) and (2) of embodiment 1. Fig.11 )

[0097] (2) By comparing and analyzing the statistical results of step 2, a new seismic gather that is more conducive to directly reflecting the changes in faults under the strong shielding of igneous rocks is screened and synthesized to extract the seismic components.

[0098] (3) Perform dip-guided filtering on the seismic components, extract variance volumes, ant tracking volumes, and other methods based on the seismic component data to better reflect the distribution of discontinuities in the data volume in order to enhance the response to the fault boundary.

[0099] (4) Volume fusion: The attributes that can better reflect the changes in the fracture boundary of the earthquake component volume are subjected to PCA-RGB volume fusion. The PCA principal component analysis technology is used to reduce the dimension of multiple earthquake attribute volumes of the same type to obtain the corresponding principal components. Then, the RGB color fusion technology is used to fuse the principal components of different types of earthquake attribute volumes in order to obtain an attribute volume that can more accurately and intuitively reflect the fracture boundary.

[0100] ( Fig.12 )

[0101] Example 3

[0102] In this embodiment, the purpose is to study the problem that the pinch-out point of the underlying sandstone reservoir is inaccurate and difficult to track due to the shielding effect of the overlying igneous rock in the igneous rock development area. In order to eliminate the influence of the igneous rock, the pinch-out range of the sand body can be accurately determined. The implementation process of step 1 and step 2 is the same as that of embodiment 1, and the only difference is step 3, which is specifically introduced as follows:

[0103] Step 3, perform synthetic record calibration and time-frequency spectrum analysis on the actual wells, and combine the model statistical results of step 2 to select the dominant frequency band range of the average thickness of the sandstone reservoir in the target well area as the basis, and then perform spectrum decomposition processing on the original seismic data body to extract the dominant time-frequency body of the sandstone reservoir in the well area, and calculate the mixed phase wavelet that can better reflect the average thickness of the sandstone in the target area according to the wavelet of each single well in the synthetic record calibration process of the actual wells. The specific process is the same as part (3) of step 3 in embodiment 1. Deconvolution is performed, and the convolution result is subjected to energy balance and AGC amplitude gain processing to obtain the final seismic data after removing the strong shielding effect of igneous rock.

[0104] The specific implementation process is as follows:

[0105] (1) Perform synthetic record calibration and time-frequency spectrum analysis on the actual well, and select the dominant frequency band range and sub-waves of the target sand body based on the statistical results of the model in step 2.

[0106] (2) Perform time-frequency decomposition on the original seismic data volume and select the time-frequency volume that can highlight the target sand body. This process uses the MP matching pursuit method, which represents the seismic signal as a linear combination of a series of time-frequency atoms that best match the local structural characteristics of the signal. Figure 7 )

[0107] (3) Calculate the mixed phase wavelet based on the multiple wavelets selected from the single well synthetic record. The specific process is the same as part (3) of step 3 in Example 1. Figure 8 )

[0108] (4) Finally, the time-frequency volume calculated in (2) and the mixed phase wavelet calculated in (3) are deconvolved, and a series of interpretative processing is performed on the convolution results to finally obtain the seismic data volume with the igneous rock shielding effect removed.

[0109] The comprehensive method of suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata of the present invention establishes a seismic identification panel of the clastic rock reservoir underlying the igneous rocks on the basis of quantitatively analyzing the seismic response characteristics of the igneous rocks and the factors affecting the surrounding rocks through forward simulation; combines the original seismic data to perform effective signal spectrum analysis, uses the wellside seismic trace to optimize the two-dimensional mixed phase graded wavelet in the time-frequency domain, combines the local dip angle information of the target layer to perform time-frequency domain deconvolution reconstruction, suppresses the strong interference of the igneous rocks, and recovers the energy of the weak signals of the underlying strata, thereby improving the effective reflection signals of the underlying strata in the igneous rock development area, forming a complete set of comprehensive methods for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata. The method weakens the interference effect of the igneous rocks, enhances the reflection characteristics of the sand body, improves the recognition accuracy of the fault, and provides a reliable basis for the description of the effective reservoirs and faults in the igneous rock development area, and has high application value and good promotion prospects.

[0110] Finally, it should be noted that the above is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions recorded in the aforementioned embodiments or replace some of the technical features therein with equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

[0111] Except for the technical features described in the specification, all other technical features are known technologies to those skilled in the art.

Claims

1. A comprehensive method for suppressing the strong shielding of igneous rocks and recovering the effective weak signal energy of the underlying strata, characterized in that: The comprehensive method of suppressing the strong shielding of igneous rock to recover the effective weak signal energy of the underlying strata includes: Step 1: Establish a corresponding forward model based on the actual geological conditions of the igneous rock development area, quantitatively analyze the seismic response characteristics of igneous rocks and the factors affecting the underlying strata, and establish an identification plate; Step 2: Based on the results of step 1, compare and analyze the changes in the main geophysical characteristics of the actual seismic data of the igneous rock and its underlying strata; Step 3: According to the result of step 2, the dominant frequency band of the igneous rock and its underlying strata is analyzed, a two-dimensional mixed phase graded wavelet is selected in the dominant frequency band of the underlying target layer to perform time-frequency domain deconvolution, and energy recovery processing is performed on the deconvolution result to obtain optimized seismic data; In step 3, according to the result of step 2, different dominant frequency bands are selected to distinguish the energy, amplitude and waveform morphology of igneous rock and underlying target layer; the frequency band selection is implemented based on the matching pursuit MP algorithm; the dominant frequency band of the underlying target layer is selected for time-frequency domain deconvolution, and the result is subjected to amplitude equalization and AGC and other energy recovery processing, and finally the optimized seismic data body is obtained; the core idea of ​​the matching pursuit MP algorithm is to represent the signal as a linear combination of a series of time-frequency atoms that best match the local structural characteristics of the signal; the essence of finding the best matching time-frequency atom is to calculate the local instantaneous properties of the seismic signal, and the superposition of the time-frequency distribution of each time-frequency atom is the matching pursuit time-frequency.

2. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 1 is characterized in that: In step 1, forward modeling quantitatively analyzes the seismic response characteristics of igneous rocks including frequency, velocity, energy and attenuation.

3. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 2 is characterized in that: In step 1, the impacts on the underlying strata include: the impact of different velocities between igneous rocks and surrounding rocks on the underlying surrounding rocks and faults; the impact of thickness changes of a single set of igneous rocks with the same velocity on the underlying surrounding rocks and faults; the impact of the occurrence of igneous rocks, including the impact of different dips on the underlying surrounding rocks and faults; and the impact of the distance between igneous rocks and underlying strata on the underlying surrounding rocks and faults.

4. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 1 is characterized in that: In step 2, the seismic geophysical characteristics of igneous rocks and their underlying strata include changes in frequency, energy and waveform kinematics, dynamics and statistics of seismic characteristics, which provide a basis for strong shielding stripping of igneous rocks and weak signal energy recovery.

5. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 1 is characterized in that: In step 3, based on the results of step 2, the dominant frequency bands of igneous rocks and their underlying strata are analyzed, and a two-dimensional mixed phase-graded wavelet is selected in the dominant frequency band of the underlying target layer for time-varying frequency deconvolution. The convolution results are subjected to energy recovery processing such as dip-guided filtering, amplitude equalization and AGC to obtain optimized seismic data.

6. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 1 is characterized in that: Step 3 includes: Step 3a, performing sparse expression of seismic signals; Step 3b, obtaining the time-frequency spectra of the decomposed different atomic signals; Step 3c, extracting multi-well mixed phase wavelets; Step 3d, using the selected wavelet for deconvolution, the convolution result is subjected to amplitude equalization, dip-guided filtering and AGC amplitude gain, which are energy recovery and enhancement processes, to obtain the final seismic data after removing the strong shielding effect of igneous rocks.

7. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 6 is characterized in that: In step 3a, the sparse expression of the seismic signal is: Where s(t) is the band-limited seismic signal, Rs (m) (t) is the residual signal after matching, m is the number of atoms in iteration, and the termination condition of iteration is judged by the set threshold, thereby determining the number of atoms in iteration; atoms are determined by four parameters t j Delay, f j Main frequency, a j Correlation coefficient, φ j Phase controlled; Given an atomic library D = {w k ,k=1,2,...K}, and satisfy ||w k ||=1, these vectors can be expanded into N-dimensional Hilbert space; for any signal s(t), find the atom w that best matches it from D kr ; The following matching conditions are met: Among them, <"> represents the inner product, and Rs (1) (t) represents the difference between the projection value and the signal, that is, the residual signal of the first iteration, so the signal can be expressed as: s(t) = <s,w k0 >w k0 +Rs (1) (t); s(t) is the initial signal, <s,w k0 > is the projection value of the first iteration, w k0 The best matching atom in the first iteration; in the initial state, Rs (0) (t) = s(t); obviously w k0 and Rs. (1) (t) is orthogonal, we can get: ||Rs (0) (t)|| 2 =| <s,w k0 >| 2 +||Rs (1) (t)|| 2 , successive iterations approach the optimization, which must minimize the energy of the residual Rs(t) and make the projection | <s,w k0 >|Maximization; after n iterations, Rs (n) (t) = a n w n +Rs (n+1) (t), a n Representing the matching factor for n iterations, the signal is finally decomposed as shown below:

8. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 7 is characterized in that: In step 3b, the analytical expression of the time-frequency spectrum is obtained: Where: a j is the wavelet factor, j W R (f,f j ) is the frequency distribution, env[w R (tt j ,f j )] is the instantaneous envelope property in the frequency domain at different time periods, is the phase change under Fourier transform, tt j is the time difference, f j is the frequency, is the phase; ① Perform Hilbert transformation on s(t) to construct a complex analysis signal: S(t) = s(t) + jHT(s(t)), and determine the instantaneous envelope and instantaneous phase formula of the signal as follows: Determine the delay t according to the corresponding time of the instantaneous envelope extreme value and the peak and valley position j , and find the instantaneous phase θ(t j ); where S(t) is the complex analysis signal of the residual signal, s(t) is the initial residual signal; jHT(s(t)) is the Hilbert change of the initial signal, A(t) is the instantaneous envelope attribute, and θ(t) is the instantaneous phase attribute; ② Obtain the instantaneous frequency of the signal, construct a matching sub-wave near the delay, and select the sub-wave W (t j, f j ,θ j ) Its corresponding main frequency is the optimal main frequency of the matching sub-wave; W(t j, f j ,θ j ) is the optimal sub-wave under different time, frequency and phase; ③According to the correlation, the wavelet factor is obtained under the condition of normalization of the wavelet energy. The matching atom with the largest correlation with the signal is subtracted from the seismic trace to obtain the residual signal, which is then iterated recursively within a certain threshold range; thus, the time-frequency spectra of all matching sub-waves are superimposed to obtain the time-frequency spectra of the initial signal.

9. The comprehensive method for suppressing igneous rock strong shielding to recover the effective weak signal energy of the underlying strata according to claim 8 is characterized in that: In step 3c, ignoring the noise term, the relationship between the seismic record and the seismic wavelet and reflection coefficient sequence can be expressed in the frequency domain by the convolution model: x(ω)=w(ω)ξ(ω); where x(ω), w(ω), ξ(ω) are the seismic record, seismic wavelet and reflection coefficient in the frequency domain respectively; after taking the logarithm, it can be transformed into a linear system lnx(ω)=lnw(ω)+lnξ(ω), and then the anti-rich transformation is performed to obtain: They are respectively called the complex spectrum sequences of x(t), w(t), ξ(t); the minimum and maximum phase components of the wavelet can be determined in the complex spectrum domain, and then a set of wavelets with different phase spectra can be determined by scanning; Based on this, multiple wells are selected to scan the wellside seismic traces in a small window in the reciprocal spectrum domain, and the optimal wavelet is selected and calculated by comparing the synthetic records with the actual time-frequency volume amplitude spectrum.

Images