Angle domain illumination compensation elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition
Through the wave mode decoupling and angle decomposition methods, the problems of vertical and transverse wave crosstalk and insufficient deep structure accuracy in elastic wave full waveform inversion are solved, and the decoupling and illumination compensation of vertical and transverse wave velocity gradients are realized, and the inversion accuracy and resolution of deep structures are improved.
Patent Information
- Application Number
- CN202510870247.7
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-06-26
- Publication Date
- 2025-08-15
AI Technical Summary
In elastic wave full waveform inversion, there are problems of vertical and transverse wave crosstalk and insufficient deep structural inversion accuracy. The existing methods have failed to effectively solve the lateral complexity of multi-parameter crosstalk and ignoring the wave field.
The full waveform inversion method of angle-domain illumination compensation elastic wave based on wave mode decoupling and angle decomposition is adopted. The wave field mode separation and angle decomposition are realized through local decomposition transformation, the local illumination matrix and local partial resolution function are calculated, and the vertical and horizontal wave velocity gradient is decoupled and illuminated compensation is performed, and the vertical and horizontal wave velocity is iteratively optimized inversion.
Effectively solve the crosstalk problem of vertical and transverse waves, improve the inversion accuracy and resolution of deep structures, avoid multi-parameter crosstalk, and improve the inversion effect of deep structures.
Smart Images
Figure CN120491171A_ABST
Abstract
Description
Technical Field
[0001] The present invention belongs to the technical field of geophysical exploration, and in particular relates to an angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition. Background Art
[0002] With the deepening of oil and gas exploration and development, the target has shifted from structural reservoirs to lithologic reservoirs, and the difficulty of exploration has gradually increased. To meet this demand, full waveform inversion has rapidly developed and become a research hotspot in the geophysical community. Full waveform inversion is an effective method for estimating formation parameters. It not only provides high-resolution and high-precision velocity models for seismic imaging, but also has the ability to simultaneously invert a variety of different parameters.
[0003] Currently, full waveform inversion has been widely developed in the field of acoustic waves. Acoustic full waveform inversion is relatively simple to implement, but the acoustic wave approximation only considers the longitudinal wave component of the wave field and ignores the transverse wave component of the wave field. In comparison, the use of elastic wave equations can more accurately characterize the propagation characteristics of the wave field, which is conducive to a more precise understanding of the formation medium information. Therefore, the study of elastic wave full waveform inversion is of great significance. However, in the case of elastic media, there is mutual coupling of multiple wave fields, which leads to the mutual crosstalk problem of multi-parameter inversion, bringing difficulties to elastic wave full waveform inversion. Scholars have studied different methods to solve the multi-parameter crosstalk problem in elastic wave inversion, such as alternating inversion, hierarchical inversion, and decoupling using the Hessian matrix.
[0004] While effectively addressing the crosstalk between P- and S-wave velocities, improving the inversion accuracy of deep structures is a crucial issue. Numerous methods have been studied and developed in this area. To compensate for the attenuation of wave energy with depth, one strategy is to simply assume the attenuation is a function of depth and use this to scale the calculated gradient. However, simple one-dimensional functions often cannot accurately model the attenuation of energy with depth. Subsequently, layered inversion methods were proposed, which first focus on shallow structure recovery using shallow reflections and then gradually recover deeper regions based on the retrieved shallow structure. These methods only consider the attenuation of wave energy with depth and ignore the lateral complexity of the wavefield. Recently, illumination-compensated waveform inversion methods based on the wave equation have also been proposed for deep structure enhancement. However, these methods do not consider angular information. More recently, an angle-domain illumination-compensated acoustic full-waveform inversion method has been proposed to improve the inversion accuracy of deep structures. This method enhances the inversion of deep structures by leveraging local angular information around a local target. However, these methods have not yet considered the case of elastic waves, which more accurately describe wave characteristics than the acoustic wave equation. Summary of the Invention
[0005] The purpose of the present invention is to provide an angle-domain illumination-compensated elastic wave full-waveform inversion method based on wave mode decoupling and angular decomposition. This method, by adopting a local decomposition transformation in the slowness domain, can not only achieve the separation of different wave field modes, but also achieve the angular decomposition of the wave field to obtain a local illumination matrix and a local resolution function, thereby simultaneously solving the problem of longitudinal and transverse wave crosstalk and improving the inversion accuracy of deep structures.
[0006] The technical solution adopted by the present invention is an angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition, specifically: Step 1: Obtain multi-component observation data ; Step 2: Calculate frequency domain seismic data ; Step 3: Construct a rectangular grid geological model; Step 4: Given the objective function ; Step 5: Given the frequency components for inversion ; Step 6: Calculate the frequency component to be Frequency domain seismic data , then decompose; Step 7: Obtain longitudinal wave field and shear wave fields ; Step 8: P-wave velocity gradient after wave mode decoupling and shear wave velocity gradient ; Step 9: Calculate the local illumination matrix and the local resolution function; Step 10, obtaining the longitudinal wave velocity gradient and the shear wave velocity gradient after illumination compensation; Step 11: For other frequency components in step 5 , repeat steps 6-10 until the longitudinal wave velocity is obtained and shear wave velocity .
[0007] The present invention is also characterized in that: Step 1 is as follows: set the shot points and receiver points, collect the original seismic data and obtain the time variable The detection point position is The gun point position is Multi-component observation data , To include Quantity and The multi-component wave field vector of the component: ; Step 2 is as follows: Calculate the frequency variable according to Fourier transform: The detection point position is The gun point position is Frequency domain seismic data , Still includes Quantity and The multi-component wave field vector of the component: ;Frequency domain seismic data It is calculated according to the following Fourier transform formula: (1) In formula (1), is the maximum observation time of the wave field; The size of the rectangular grid geological model constructed in step 3 is: horizontal grid points, vertically grid points, and set the spatial sampling interval and time sampling interval of the forward simulation , maximum sampling time ; Among them, the spatial sampling interval includes the horizontal sampling interval and longitudinal sampling interval .
[0008] Step 4 is as follows: Given the longitudinal wave velocity model and shear wave velocity model Initial P-wave velocity model and the initial shear wave velocity model , given the objective function ; In step 4, set the objective function The bi-norm of the residual between the frequency domain observation wavefield and the frequency domain calculation wavefield has the following form: (2) In formula (2), and Respectively represent the frequency domain observation wave field and the frequency domain calculation wave field; represents transpose, Indicates taking conjugate; Among them, the frequency domain observation wave field By multi-component observation data Perform the Fourier transform described in step 2 to obtain; Frequency domain wavefield calculation It is obtained by the following method: According to the time domain elastic wave equation, based on the initial longitudinal wave velocity model and the initial shear wave velocity model Calculate the time domain wave field , then Perform the Fourier transform described in step 2 to obtain ; Among them, the time domain calculates the wave field The following time domain elastic wave equation is satisfied: (3) In formula (3), represents the time domain wave field, To include Quantity and The multi-component wave field vector of the component: ; and Represent the time variables as ,lie in The wave field Quantity and component source functions; Indicates the formation density, which is set as a constant; and represents the Lame coefficient, which is related to the initial longitudinal wave velocity model A specific location P-wave velocity at and the initial shear wave velocity model A specific location Shear wave velocity at The following relationship exists: (4).
[0009] In step 5, the frequency components used for inversion are given , N Take 15, The corresponding settings are 6Hz, 7Hz, 8Hz, 9Hz, 10Hz, 11Hz, 12Hz, 13Hz, 14Hz, 15Hz, 16Hz, 17Hz, 18Hz, 19Hz, and 20Hz respectively.
[0010] In step 6, the decomposed frequency domain seismic data Obtain the longitudinal wave field in the frequency angle domain and shear wave fields ,in represents the angle component; Among them, the longitudinal wave field in the frequency angle domain Including source side wave field and the receiving point side wave field , shear wave field in frequency-angle domain Also includes the source side wave field and the receiving point side wave field ,in and Represents spatial position The wave field propagation direction on the source side and the wave field propagation direction on the receiving point side; The specific implementation of step 6 is as follows: First, use the method in step 2 to calculate the frequency component Frequency domain seismic data Decompose to obtain the longitudinal wave field in the frequency angle domain and shear wave fields , the specific formula is: (5) In formula (5), represents the slowness vector, where is the unit vector pointing to the direction of wave field propagation, is the slowness vector The absolute value of Indicates the coordinate point The sampling window centered at "and" " represent cross product operation and dot product operation respectively; Then Substitution In Substitution The angle domain longitudinal wave field is obtained and shear wave fields , the specific formula is: (6) In formula (6), and represents the absolute value of the longitudinal wave slowness and the absolute value of the shear wave slowness, that is, , ,in and Indicates the average longitudinal wave velocity and the average shear wave velocity within the sampling window; Combining equations (5) and (6), we can get the angle domain longitudinal wave field: and shear wave fields ; Finally, according to equations (5) and (6), the source side wave field is and the receiving point side wave field Decomposition is performed, where the two superscript arrows are used to distinguish the source side wave field and the receiver side wave field, and the source side angle domain longitudinal wave field is obtained. and shear wave fields , and the longitudinal wave field in the angle domain at the receiving point and shear wave fields .
[0011] Step 7 is as follows: the frequency-angle domain longitudinal wave field obtained in step 6 is and shear wave fields Perform angle superposition and obtain the frequency component without angle information: Local location P-wave field at and shear wave fields , the expression is as follows: (7) (8) In formula (7) and formula (8), and They represent the frequency components without angle information on the source side. Local location The longitudinal wave field and the shear wave field at and They represent the frequency components without angle information at the receiving point side. Local location The longitudinal wave field and the shear wave field at and They represent the sum of the wave field propagation angles on the source side and the sum of the wave field propagation angles on the receiving point side. Abbreviated as , accordingly, the shear wave field in the angle domain on the source side Abbreviated as , longitudinal wave field in the angle domain at the receiving point side Abbreviated as , Shear wave field in the angle domain at the receiving point side Abbreviated as .
[0012] Step 8 is implemented as follows: For the P-wave velocity model and shear wave velocity model , using the longitudinal wave field without angle information on the source side obtained in step 7 and shear wave fields And the longitudinal wave field without angle information at the receiving point and shear wave fields Obtain the longitudinal wave velocity gradient after wave mode decoupling and shear wave velocity gradient , specifically: First, the velocity gradient formulas for different wave field modes without angle information are obtained by the chain rule as follows: (9) (10) (11) (12) In equations (9) to (12), “Re” represents the real part. After using equations (9) to (12) to obtain the velocity gradients of different wave field modes without angle information, the P-wave velocity gradient after wave mode decoupling is obtained: and shear wave velocity gradient , the specific expression is: (13).
[0013] Step 9 is as follows: Use the frequency angle domain longitudinal wave field and shear wave fields Obtain PP local illumination matrix and SS local lighting matrix ; Then use the PP local lighting matrix and SS local lighting matrix Obtaining the PP local resolution function and SS local resolution function ,in Indicates spatial location The interface dip wave number at ; In step 9, the PP local lighting matrix and SS local lighting matrix Calculated as follows: (14) In formula (14), represents the sum of the source and the receiving points, Indicates spatial location From the epicenter location The angle of propagation is The frequency is The incident wave field, Indicates spatial location From the receiving point The angle of propagation is The frequency is The scattered wave field.
[0014] In step 9, the PP local resolution function and SS local resolution function Calculated as follows: (15) In formula (15), represents the interface inclination angle, represents the background wave number, Indicates when , PP local illumination matrix when Indicates when , SS local lighting matrix when .
[0015] Step 10 is as follows: Using the PP local resolution function P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated longitudinal wave velocity gradient ; Using SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated shear wave velocity gradient ; For the P-wave velocity model and shear wave velocity model , respectively using the longitudinal wave velocity gradient after illumination compensation based on wave mode decoupling and shear wave velocity gradient The objective function is optimized by iterative optimization algorithm. Optimize and obtain the inverted longitudinal wave velocity based on wave mode decoupling and illumination compensation and shear wave velocity ; In step 10, the PP local resolution function is used P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the longitudinal wave velocity gradient after illumination compensation and using the SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the shear wave velocity gradient after illumination compensation The specific implementation is as follows: (16) In formula (16), For The wavenumber domain gradient value obtained by local Fourier transform is For The wavenumber domain gradient value obtained by local Fourier transform is calculated as follows: (17) In formula (17), represents the spatial integration range, represents the imaginary unit; Equation (16) represents the illumination compensation of velocity gradient in the wavenumber domain, Through Perform inverse Fourier transform to obtain, Can be achieved through Perform inverse Fourier transform and the calculation formula is as follows: (18) In formula (18), Indicates the wavenumber integration range; Calculate the longitudinal wave velocity gradient based on wave mode decoupling and illumination compensation and shear wave velocity gradient Then, the longitudinal wave velocity is calculated using the following criterion: and shear wave velocity Perform iterative updates: (19) In formula (19), Respectively Step and The longitudinal wave velocity value of the step iteration, Respectively Step and The shear wave velocity value of the step iteration; is the iteration step size, which is a positive constant.
[0016] The beneficial effects of the present invention are: The method of the present invention: First, starting from a smooth initial model, the lowest frequency component in the seismic data is decomposed and transformed to obtain the P-wave and S-wave wavefields at different angles. Second, based on the decomposed wavefield, the P-wave and S-wave wavefields without angle information are obtained, thereby obtaining the P-wave velocity gradient and S-wave velocity gradient after wave mode decoupling. The P-wave and S-wave wavefields with angle information are used to calculate the local illumination matrix and local resolution function. Third, based on the P-wave velocity gradient and S-wave velocity gradient after wave mode decoupling, angle domain illumination compensation is used to improve the inversion resolution of deep structures. Finally, all frequency components are traversed until the final inverted fine P-wave velocity structure and S-wave velocity structure are obtained. Compared with commonly used inversion methods, this method can effectively avoid the multi-parameter crosstalk problem in elastic wave full waveform inversion and the problem of poor resolution in deep inversion. The reasons why this method has the above advantages are: first, because the wave field is separated and the longitudinal wave velocity and shear wave velocity are decoupled according to their scattering characteristics, crosstalk between the two can be avoided; second, because the local resolution function is taken into account in the inversion process, the relationship between the target position inclination and the incident wave and the scattered wave can be fully considered, and the inversion effect of deep structures can be effectively improved through angle domain illumination compensation. BRIEF DESCRIPTION OF THE DRAWINGS
[0017] Figure 1 is a flow chart of the method of the present invention; Figure 2is the longitudinal wave velocity model of the Marmousi elastic velocity model; Figure 3 is the shear wave velocity model of the Marmousi elastic velocity model; Figure 4 is the smoothed initial P-wave velocity model; Figure 5 is the smoothed initial shear wave velocity model; Figure 6 It is the P-wave velocity model obtained by inversion using the traditional method; Figure 7 It is the shear wave velocity model obtained by inversion using the traditional method; Figure 8 is the final refined P-wave velocity model obtained by inversion using the method of the present invention; Figure 9 It is the final refined shear wave velocity model obtained by inversion using the method of the present invention. DETAILED DESCRIPTION
[0018] The present invention will be described in detail below with reference to the accompanying drawings and specific embodiments.
[0019] The present invention provides an angle domain illumination compensation elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition, such as Figure 1 The specific steps are as follows: Step 1: Set the shot points and detection points, collect the original seismic data and obtain the time variable as The detection point position is The gun point position is Multi-component observation data , To include Quantity and The multi-component wave field vector of the component: ; Step 2: Calculate the frequency variable according to Fourier transform: The detection point position is The gun point position is Frequency domain seismic data , Still includes Quantity and The multi-component wave field vector of the component: ;Frequency domain seismic data It is calculated according to the following Fourier transform formula: (1) In formula (1), is the maximum observation time of the wave field.
[0020] Step 3: Construct a rectangular grid geological model and set the model size to: horizontal grid points, vertically grid points, and set the spatial sampling interval and time sampling interval of the forward simulation , maximum sampling time ; Among them, the spatial sampling interval includes the horizontal sampling interval and longitudinal sampling interval ; In practice, relevant parameters can be determined based on the effective frequency band of seismic data, the length of seismic recording time, and the seismic data observation system. The standard for selecting grid parameters is to ensure that when performing finite difference forward modeling based on the grid, it can not only meet the stability conditions but also effectively suppress numerical dispersion. Step 4: Given the longitudinal wave velocity model and shear wave velocity model Initial P-wave velocity model and the initial shear wave velocity model , given the objective function ; In step 4, set the objective function The bi-norm of the residual between the frequency domain observation wavefield and the frequency domain calculation wavefield has the following form: (2) In formula (2), and Respectively represent the frequency domain observation wave field and the frequency domain calculation wave field; represents transpose, Indicates taking conjugate; Among them, the frequency domain observation wave field By multi-component observation data Perform the Fourier transform described in step 2 to obtain; Frequency domain wavefield calculation It is obtained by the following method: According to the time domain elastic wave equation, based on the initial longitudinal wave velocity model and the initial shear wave velocity model Calculate the time domain wave field , then Perform the Fourier transform described in step 2 to obtain ; Among them, the time domain calculates the wave field The following time domain elastic wave equation is satisfied: (3) In formula (3), represents the time domain wave field, To include Quantity and The multi-component wave field vector of the component: ; and Represent the time variables as ,lie in The wave field Quantity and component source functions; Indicates the formation density, which is set as a constant; and represents the Lame coefficient, and the longitudinal wave velocity model A specific location P-wave velocity at and shear wave velocity model A specific location Shear wave velocity at The following relationship exists: (4); Step 5: Given the frequency components for inversion , N Take 15, The corresponding settings are 6Hz, 7Hz, 8Hz, 9Hz, 10Hz, 11Hz, 12Hz, 13Hz, 14Hz, 15Hz, 16Hz, 17Hz, 18Hz, 19Hz, 20Hz respectively; Step 6: Use the method in step 2 to calculate the frequency component: Frequency domain seismic data , and then decompose it to get the longitudinal wave field in the frequency angle domain and shear wave fields ,in represents the angle component; Among them, the longitudinal wave field in the frequency angle domain Including source side wave field and the receiving point side wave field , shear wave field in frequency-angle domain Also includes the source side wave field and the receiving point side wave field ,in and Represents spatial position The wave field propagation direction on the source side and the wave field propagation direction on the receiving point side; The specific implementation of step 6 is as follows: First, use the method in step 2 to calculate the frequency component Frequency domain seismic data Decompose to obtain the longitudinal wave field in the frequency angle domain and shear wave fields , the specific formula is: (5) In formula (5), represents the slowness vector, where is the unit vector pointing to the direction of wave field propagation, is the slowness vector The absolute value of Indicates the coordinate point The sampling window centered at "and" " represent cross product operation and dot product operation respectively; Then Substitution In Substitution The angle domain longitudinal wave field is obtained and shear wave fields , the specific formula is: (6) In formula (6), and represents the absolute value of the longitudinal wave slowness and the absolute value of the shear wave slowness, that is, , ,in and Indicates the average longitudinal wave velocity and the average shear wave velocity within the sampling window; Combining equations (5) and (6), we can get the angle domain longitudinal wave field: and shear wave fields ; Finally, according to equations (5) and (6), the source side wave field is and the receiving point side wave field Decomposition is performed, where the two superscript arrows are used to distinguish the source side wave field and the receiver side wave field, and the source side angle domain longitudinal wave field is obtained. and shear wave fields , and the longitudinal wave field in the angle domain at the receiving point and shear wave fields .
[0021] Step 7: The frequency-angle domain longitudinal wave field obtained in step 6 and shear wave fields Perform angle superposition and obtain the frequency component without angle information: Local location P-wave field at and shear wave fields , the expression is as follows: (7) (8) In formula (7) and formula (8), and They represent the frequency components without angle information on the source side. Local location The longitudinal wave field and the shear wave field at and They represent the frequency components without angle information at the receiving point side. Local location The longitudinal wave field and the shear wave field at and They represent the sum of the wavefield propagation angles on the source side and the sum of the wavefield propagation angles on the receiving point side respectively; To simplify the description, the longitudinal wave field in the source side angle domain is Abbreviated as , accordingly, the shear wave field in the angle domain on the source side Abbreviated as , longitudinal wave field in the angle domain at the receiving point side Abbreviated as , Shear wave field in the angle domain at the receiving point side Abbreviated as .
[0022] Step 8: For the longitudinal wave velocity model and shear wave velocity model , using the longitudinal wave field without angle information on the source side obtained in step 7 and shear wave fields And the longitudinal wave field without angle information at the receiving point and shear wave fields Obtain the longitudinal wave velocity gradient after wave mode decoupling and shear wave velocity gradient , the specific implementation is as follows: During elastic wave propagation, longitudinal and transverse waves will generate conversion waves after scattering. There are multiple modes of waves in the entire wave field, namely PP waves (longitudinal waves-longitudinal waves), PS waves (longitudinal waves-transverse waves), SP waves (transverse waves-longitudinal waves), and SS waves (transverse waves-transverse waves). Different modes of waves have different sensitivities to different medium parameters, that is, different modes of waves will be generated after scattering with different parameters. Mainly produces PP waves, while the shear wave speed Generate SS wave, SP wave, PS wave and PP wave, among which SS wave is the main one. and The mutual crosstalk problem, inversion Only PP waves are used when inverting When only SS, SP and PS waves are used, the velocity gradient formulas of different wave field modes without angle information are first obtained by the chain rule under the condition of wave mode decoupling: (9) (10) (11) (12) In equations (9) to (12), “Re” represents the real part. After using equations (9) to (12) to obtain the velocity gradients of different wave field modes without angle information, the P-wave velocity gradient after wave mode decoupling is obtained: and shear wave velocity gradient , the specific expression is: (13); Step 9: Use the frequency-angle domain longitudinal wave field and shear wave fields Obtain PP local illumination matrix and SS local lighting matrix ; Then use the PP local lighting matrix and SS local lighting matrix Obtaining the PP local resolution function and SS local resolution function ,in Indicates spatial location The interface dip wave number at ; In step 9, the PP local lighting matrix and SS local lighting matrix Calculated as follows: (14) In formula (14), represents the sum of the source and the receiving points, Indicates spatial location From the epicenter location The angle of propagation is The frequency is The incident wave field, Indicates spatial location From the receiving point The angle of propagation is The frequency is The scattered wave field.
[0023] In step 9, the PP local resolution function and SS local resolution function Calculated as follows: (15) In formula (15), represents the interface inclination angle, represents the background wave number, Indicates when , PP local illumination matrix when Indicates when , SS local lighting matrix when .
[0024] Step 10: Use the PP local resolution function P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated longitudinal wave velocity gradient ; Using SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated shear wave velocity gradient ; For the P-wave velocity model and shear wave velocity model , respectively using the longitudinal wave velocity gradient after illumination compensation based on wave mode decoupling and shear wave velocity gradient The objective function is optimized by iterative optimization algorithm. Optimize and obtain the inverted longitudinal wave velocity based on wave mode decoupling and illumination compensation and shear wave velocity ; In step 10, the PP local resolution function is used P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the longitudinal wave velocity gradient after illumination compensation and using the SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the shear wave velocity gradient after illumination compensation The specific implementation is as follows: (16) In formula (16), For The wavenumber domain gradient value obtained by local Fourier transform is For The wavenumber domain gradient value obtained by local Fourier transform is calculated as follows: (17) In formula (17), represents the spatial integration range, represents the imaginary unit; Equation (16) represents the illumination compensation of velocity gradient in the wavenumber domain, Through Perform inverse Fourier transform to obtain, Can be achieved through Perform inverse Fourier transform and the calculation formula is as follows: (18) In formula (18), Indicates the wavenumber integration range; Calculate the longitudinal wave velocity gradient based on wave mode decoupling and illumination compensation and shear wave velocity gradient Then, the longitudinal wave velocity is calculated using the following criterion: and shear wave velocity Perform iterative updates: (19) In formula (19), Respectively Step and The longitudinal wave velocity value of the step iteration, Respectively Step and The shear wave velocity value of the step iteration; is the iteration step size, which is a positive constant.
[0025] The above method is used to calculate the longitudinal wave velocity gradient after illumination compensation based on wave mode decoupling. and shear wave velocity gradient And iterate, after the iteration converges, the longitudinal wave velocity obtained by angle domain illumination compensation elastic wave full waveform inversion based on wave mode decoupling and angle decomposition can be obtained and shear wave velocity .
[0026] Step 11: For other frequency components in step 5 Repeat steps 6 to 10 until the longitudinal wave velocity is obtained. and shear wave velocity , which is the final inversion result.
[0027] Example 1 The specific implementation process of the present invention is applied to the Marmousi elastic velocity model. The model has a total of 4600 meters in the horizontal direction and 1500 meters in the vertical direction. The longitudinal wave velocity model and the shear wave velocity model of the Marmousi elastic velocity model are respectively as follows: Figure 2 and Figure 3As shown. Set the number of shot points to 45 and the number of receiver points to 230. The shot points and receiver points are evenly spaced at the top of the Marmousi elastic velocity model. The Marmousi elastic velocity model is gridded at equal intervals in the horizontal and vertical directions, with 460 grids in the horizontal direction and 150 grids in the vertical direction. The horizontal sampling interval and the vertical sampling interval are both 10m. The time sampling interval for recording wave field data is 0.001s. According to the maximum velocity of the Marmousi elastic velocity model of 4700m / s, the relationship between the time step and the grid interval, , meeting the stability condition; according to the relationship between the minimum velocity of the Marmousi elastic velocity model of 1260 m / s, the main frequency of the earthquake source of 10 Hz and the grid spacing, , satisfying the dispersion condition.
[0028] The initial models of the P-wave velocity model and S-wave velocity model are set as smooth initial models, respectively. Figure 4 and Figure 5 shown.
[0029] The P-wave velocity model and S-wave velocity model obtained by the traditional inversion method are as follows: Figure 6 and Figure 7 As shown in Figure 2, it can be seen that the velocity model obtained by the traditional inversion method is affected by the multi-parameter crosstalk problem and the poor resolution of deep inversion.
[0030] The final refined P-wave velocity model and S-wave velocity model obtained by inversion using the method proposed in the present invention are as follows: Figure 8 and Figure 9 As shown. In step 5, N Take 15, The corresponding settings are 6Hz, 7Hz, 8Hz, 9Hz, 10Hz, 11Hz, 12Hz, 13Hz, 14Hz, 15Hz, 16Hz, 17Hz, 18Hz, 19Hz, and 20Hz respectively. Figure 6 and Figure 7 It can be seen that during the inversion process, wave mode decoupling avoids the impact of multi-parameter crosstalk on the inversion results. By considering the local resolution function during the inversion process, the use of angle-domain illumination compensation effectively improves the inversion resolution of deep structures. It can be seen that the obtained P-wave and S-wave velocities are close to the true model, demonstrating the effectiveness of the proposed method.
[0031] Example 2 Step 1: Obtain multi-component observation data ; Step 2: Calculate frequency domain seismic data ; Step 3: Construct a rectangular grid geological model; Step 4: Given the objective function ; Step 5: Given the frequency components for inversion ; Step 6: Calculate the frequency component to be Frequency domain seismic data , then decompose; Step 7: Obtain longitudinal wave field and shear wave fields ; Step 8: P-wave velocity gradient after wave mode decoupling and shear wave velocity gradient ; Step 9: Calculate the local illumination matrix and the local resolution function; Step 10, obtaining the longitudinal wave velocity gradient and the shear wave velocity gradient after illumination compensation; Step 11: For other frequency components in step 5 , repeat steps 6-10 until the longitudinal wave velocity is obtained and shear wave velocity .
[0032] Example 3 The difference from Example 2 is that: Step 1 specifically includes: setting shot points and detection points, collecting original seismic data to obtain the time variable The detection point position is The gun point position is Multi-component observation data , To include Quantity and The multi-component wave field vector of the component: .
[0033] Example 4 The difference from Example 3 is that step 2 is specifically: the frequency variable is calculated according to Fourier transform as The detection point position is The gun point position is Frequency domain seismic data , Still includes Quantity and The multi-component wave field vector of the component: ;Frequency domain seismic data It is calculated according to the following Fourier transform formula: (1) In formula (1), is the maximum observation time of the wave field.
[0034] Example 5 The difference from Example 4 is that the size of the rectangular grid geological model constructed in step 3 is: grid points, vertically grid points, and set the spatial sampling interval and time sampling interval of the forward simulation , maximum sampling time ; Among them, the spatial sampling interval includes the horizontal sampling interval and longitudinal sampling interval .
[0035] Example 6 The difference from Example 5 is that step 4 is specifically: given a longitudinal wave velocity model and shear wave velocity model Initial P-wave velocity model and the initial shear wave velocity model , given the objective function .
Claims
1. An angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition, characterized by: Specifically: Step 1: Obtain multi-component observation data ; Step 2: Calculate frequency domain seismic data ; Step 3: Construct a rectangular grid geological model; Step 4: Given the objective function ; Step 5: Given the frequency components for inversion ; Step 6: Calculate frequency domain seismic data , then decompose; Step 7: Obtain longitudinal wave field and shear wave fields ; Step 8: P-wave velocity gradient after wave mode decoupling and shear wave velocity gradient ; Step 9: Calculate the local illumination matrix and the local resolution function; Step 10, obtaining the longitudinal wave velocity gradient and the shear wave velocity gradient after illumination compensation; Step 11: For other frequency components , repeat steps 6-10 until the longitudinal wave velocity is obtained and shear wave velocity .
2. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 1 is characterized in that: Step 1 is as follows: set the shot points and receiver points, collect the original seismic data and obtain the time variable The detection point position is The gun point position is Multi-component observation data , To include Quantity and The multi-component wave field vector of the component: ; Step 2 is as follows: Calculate the frequency variable according to Fourier transform: The detection point position is The gun point position is Frequency domain seismic data , Still includes Quantity and The multi-component wave field vector of the component: ;Frequency domain seismic data It is calculated according to the following Fourier transform formula: (1) In formula (1), is the maximum observation time of the wave field; The size of the rectangular grid geological model constructed in step 3 is: horizontal grid points, vertically grid points, and set the spatial sampling interval and time sampling interval of the forward simulation , maximum sampling time ; Among them, the spatial sampling interval includes the horizontal sampling interval and longitudinal sampling interval .
3. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 2 is characterized in that: Step 4 is as follows: Given the P-wave velocity model and shear wave velocity model Initial P-wave velocity model and the initial shear wave velocity model , given the objective function ; In step 4, set the objective function The bi-norm of the residual between the frequency domain observation wavefield and the frequency domain calculation wavefield has the following form: (2) In formula (2), and Respectively represent the frequency domain observation wave field and the frequency domain calculation wave field; represents transpose, Indicates taking conjugate; Among them, the frequency domain observation wave field By multi-component observation data Perform the Fourier transform described in step 2 to obtain; Frequency domain wavefield calculation It is obtained by the following method: According to the time domain elastic wave equation, based on the initial longitudinal wave velocity model and the initial shear wave velocity model Calculate the time domain wave field , then Perform the Fourier transform described in step 2 to obtain ; Among them, the time domain calculates the wave field The following time domain elastic wave equation is satisfied: (3) In formula (3), represents the time domain wave field, To include Quantity and The multi-component wave field vector of the component: ; and Represent the time variables as ,lie in The wave field Quantity and component source functions; Indicates the formation density, which is set as a constant; and represents the Lame coefficient, which is related to the initial longitudinal wave velocity model A specific location P-wave velocity at and the initial shear wave velocity model A specific location Shear wave velocity at The following relationship exists: (4)。 4. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 3 is characterized in that: In step 5, the frequency components used for inversion are given , N Take 15, The corresponding settings are 6Hz, 7Hz, 8Hz, 9Hz, 10Hz, 11Hz, 12Hz, 13Hz, 14Hz, 15Hz, 16Hz, 17Hz, 18Hz, 19Hz, and 20Hz respectively.
5. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 4 is characterized in that: In step 6, the decomposed frequency domain seismic data Obtain the longitudinal wave field in the frequency angle domain and shear wave fields ,in represents the angle component; Among them, the longitudinal wave field in the frequency angle domain Including source side wave field and the receiving point side wave field , frequency-angle domain shear wave field Also includes the source side wave field and the receiving point side wave field ,in and Represents spatial position The wave field propagation direction on the source side and the wave field propagation direction on the receiving point side; The specific implementation of step 6 is as follows: First, use the method in step 2 to calculate the frequency component Frequency domain seismic data Decompose to obtain the longitudinal wave field in the frequency angle domain and shear wave fields , the specific formula is: (5) In formula (5), represents the slowness vector, where is the unit vector pointing to the direction of wave field propagation, is the slowness vector The absolute value of Indicates the coordinate point The sampling window centered at " "and" " represent cross product operation and dot product operation respectively; Then Substitution In Substitution The angle domain longitudinal wave field is obtained and shear wave fields , the specific formula is: (6) In formula (6), and represents the absolute value of the longitudinal wave slowness and the absolute value of the shear wave slowness, that is, , ,in and Indicates the average longitudinal wave velocity and the average shear wave velocity within the sampling window; Combining equations (5) and (6), we can get the angle domain longitudinal wave field: and shear wave fields ; Finally, according to equations (5) and (6), the source side wave field is and the receiving point side wave field Decomposition is performed, where the two superscript arrows are used to distinguish the wave field on the source side and the wave field on the receiving point side, and the longitudinal wave field in the angle domain on the source side is obtained. and shear wave fields , and the longitudinal wave field in the angle domain at the receiving point and shear wave fields .
6. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 5 is characterized in that: Step 7 is as follows: the frequency-angle domain longitudinal wave field obtained in step 6 is and shear wave fields Perform angle superposition and obtain the frequency component without angle information: Local location P-wave field at and shear wave fields , the expression is as follows: (7) (8) In formula (7) and formula (8), and They represent the frequency components without angle information on the source side. Local location The longitudinal wave field and the shear wave field at and They represent the frequency components without angle information at the receiving point side. Local location The longitudinal wave field and the shear wave field at and They represent the sum of the wave field propagation angles on the source side and the sum of the wave field propagation angles on the receiving point side respectively; in order to simplify the expression, the longitudinal wave field in the angle domain on the source side will be Abbreviated as , accordingly, the shear wave field in the angle domain on the side of the earthquake source is Abbreviated as , longitudinal wave field in the angle domain at the receiving point side Abbreviated as , Shear wave field in the angle domain at the receiving point side Abbreviated as .
7. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 6 is characterized in that: Step 8 is implemented as follows: For the P-wave velocity model and shear wave velocity model , using the longitudinal wave field without angle information on the source side obtained in step 7 and shear wave fields And the longitudinal wave field without angle information at the receiving point and shear wave fields Obtain the longitudinal wave velocity gradient after wave mode decoupling and shear wave velocity gradient , specifically: First, the velocity gradient formulas for different wave field modes without angle information are obtained by the chain rule as follows: (9) (10) (11) (12) In equations (9) to (12), "Re" represents the real part. After using equations (9) to (12) to obtain the velocity gradients of different wave field modes without angle information, the longitudinal wave velocity gradient after wave mode decoupling is obtained: and shear wave velocity gradient , the specific expression is: (13)。 8. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 7 is characterized in that: Step 9 is as follows: Use the frequency angle domain longitudinal wave field and shear wave fields Obtain PP local illumination matrix and SS local lighting matrix ; Then use the PP local lighting matrix and SS local lighting matrix Obtaining the PP local resolution function and SS local resolution function ,in Indicates spatial location The interface dip wave number at ; In step 9, the PP local lighting matrix and SS local lighting matrix Calculated as follows: (14) In formula (14), represents the sum of the source and the receiving points, Indicates spatial location From the epicenter location The angle of propagation is The frequency is The incident wave field, Indicates spatial location From the receiving point The angle of propagation is The frequency is The scattered wave field; In step 9, the PP local resolution function and SS local resolution function Calculated as follows: (15) In formula (15), represents the interface inclination angle, represents the background wave number, Indicates when , PP local illumination matrix when Indicates when , SS local lighting matrix when .
9. The angle-domain illumination-compensated elastic wave full waveform inversion method based on wave mode decoupling and angle decomposition according to claim 8, characterized in that: Step 10 is as follows: Using the PP local resolution function P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated longitudinal wave velocity gradient ; Using SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the compensated shear wave velocity gradient ; For the P-wave velocity model and shear wave velocity model , respectively using the longitudinal wave velocity gradient after illumination compensation based on wave mode decoupling and shear wave velocity gradient The objective function is optimized by iterative optimization algorithm. Optimize and obtain the inverted longitudinal wave velocity based on wave mode decoupling and illumination compensation and shear wave velocity ; In step 10, the PP local resolution function is used P-wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the longitudinal wave velocity gradient after illumination compensation and using the SS local resolution function Shear wave velocity gradient after wave mode decoupling Perform illumination compensation to obtain the shear wave velocity gradient after illumination compensation The specific implementation is as follows: (16) In formula (16), For The wavenumber domain gradient value obtained by local Fourier transform is For The wavenumber domain gradient value obtained by local Fourier transform is calculated as follows: (17) In formula (17), represents the spatial integration range, represents the imaginary unit; Equation (16) represents the illumination compensation of velocity gradient in the wavenumber domain, Through Perform inverse Fourier transform to obtain, Can be achieved through Perform inverse Fourier transform and the calculation formula is as follows: (18) In formula (18), Indicates the wavenumber integration range; Calculate the longitudinal wave velocity gradient based on wave mode decoupling and illumination compensation and shear wave velocity gradient Then, the longitudinal wave velocity is calculated using the following criterion: and shear wave velocity Perform iterative updates: (19) In formula (19), Respectively Step and The longitudinal wave velocity value of the step iteration, Respectively Step and The shear wave velocity value of the step iteration; is the iteration step size, which is a positive constant.