A high-precision three-dimensional elastic wave full waveform inversion method

Abstract

This invention discloses a high-precision three-dimensional elastic wave full-waveform inversion method. First, several discrete frequency points are selected from the observation data and processed using a conventional band-limited wavelet source to form a new source. Then, the forward and reverse propagating wave fields are obtained using the new source, and a gradient is constructed using Lamé parameters as model parameters. Finally, the L-BFGS method, which has second-order convergence, is used to calculate the inverse of the Hessian matrix and continuously update the initial velocity field model, obtaining the final velocity field model parameters, which represent the high-precision three-dimensional imaging of the desired region. This method enhances the energy of low-frequency signals during inversion, effectively increasing the contribution of low-frequency signals to the low wavenumber information of the model, thereby effectively reducing the dependence of the inversion on the initial velocity field model and improving the inversion accuracy, achieving the goal of high-precision inversion.

Description

Technical Field

[0001] This invention relates to an imaging method for underground three-dimensional elastic media, and more particularly to a three-dimensional elastic wave full waveform inversion method for resource exploration such as oil and coal mines, belonging to the field of seismic exploration velocity modeling technology. Background Technology

[0002] The exploitation of underground resources relies heavily on the accurate identification of target areas. This includes precise reservoir identification and appropriate well placement in oil and gas exploration, and accurate identification of key strata and the delineation of hidden hazardous factors in coal mining. To ensure the safe exploitation of underground resources, it is essential and urgent to ascertain the geological structure distribution and lithological variations of the target area. Seismic exploration, a branch of geophysical exploration, possesses the capability to detect fine structures and plays a dominant role in resource exploration. As a high-precision seismic exploration technique, full-waveform inversion analyzes the propagation of seismic waves underground to infer the structure and properties of the subsurface medium.

[0003] Because geological bodies are objectively distributed in three dimensions, two-dimensional full-waveform inversion cannot accurately depict their spatial distribution. Three-dimensional full-waveform inversion technology has become one of the research hotspots in the field of geophysical exploration. Compared to two-dimensional full-waveform inversion, three-dimensional full-waveform inversion requires solving more complex equations and their discrete equations, demanding tens or even hundreds of times more computational resources and efficiency than two-dimensional problems. Simultaneously, the complexity and ambiguity of the inversion increase.

[0004] Current 3D inversion methods are slower than 2D inversion. Due to the dominance of the dominant frequency in time-domain inversion, frequency information is often underutilized. This significantly impacts the accuracy of 3D full-waveform inversion, especially since it is highly dependent on the initial model, leading to reduced computational accuracy. The reason for this is that low-frequency information in the data cannot be effectively used to update the high wavenumber information in the model. For iterative algorithms like full-waveform inversion, inaccurate low-frequency inversion results in significant deviations in high-wavenumber updates. The reasons for this ineffective utilization of low-frequency information are twofold: 1. The data itself contains very limited and weak low-frequency information; 2. Current full-waveform inversion methods cannot effectively utilize this weak low-frequency signal. Although commonly used multi-scale inversion strategies based on convolution alleviate these problems to some extent, they still cannot solve them completely. Furthermore, poor reference trace quality during inversion leads to unsatisfactory results, all of which reduce inversion accuracy.

[0005] Therefore, the question is how to implement a high-precision three-dimensional full-waveform inversion method that can effectively use low-frequency information from observation data to update high-wavenumber information in the model, and ultimately ensure the accuracy of three-dimensional inversion through multi-scale inversion. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention provides a high-precision three-dimensional elastic wave full waveform inversion method, which can effectively use low-frequency information in the observation data to update the high wavenumber information in the model, and ultimately ensure the accuracy of the three-dimensional inversion through multi-scale inversion.

[0007] To achieve the above objectives, the technical solution adopted by this invention is: a high-precision three-dimensional elastic wave full waveform inversion method, the specific steps of which are as follows:

[0008] Step 1: Obtain the observation data collected by the observation system. First, give an initial velocity field model, select several discrete frequency points from the observation data, and combine them with multiple frequency points in the conventional band-limited wavelet source to perform amplitude normalization processing to form a new source.

[0009] Step 2: Using the new seismic source obtained in Step 1, calculate the forward propagation wave field using the three-dimensional elastic wave method, and combine it with the wave field data recorded by the observation system to obtain synthetic data;

[0010] Step 3: Establish a least squares objective function based on the synthetic data obtained in Step 2 and perform calculations. The calculated data residuals are used as the back-transmission source.

[0011] Step 4: Calculate the anti-propagating wave field under the anti-propagating source condition based on the anti-propagating source information obtained in Step 3.

[0012] Step 5: Construct the gradient of the model parameters using the forward propagation wave field obtained in Step 2 and the reverse propagation wave field obtained in Step 4;

[0013] Step 6: Based on the gradient obtained in Step 5, the inverse matrix of the Hessian matrix is ​​calculated using the L-BFGS method with second-order convergence, and the initial velocity field model in Step 1 is continuously updated until the data residual reaches the set threshold or the number of iterations reaches the preset value. Then, the iteration update is stopped to obtain the final velocity field model parameters. The velocity field model parameters at this time are the three-dimensional high-precision imaging of the required area.

[0014] Furthermore, the specific formula for normalization to form a new seismic source in step one is as follows:

[0015]

[0016] Where F represents the Fourier transform, Indicates frequency ω i The inverse Fourier transform of the spectrum, x s Indicates the location of the earthquake source. Let ω represent the normalization of signal A, N represent the number of frequencies used, and ω = (ω1, ω2, ..., ω...). N ), s(xs ,t) is a conventional band-limited wavelet source, where t represents time and i represents a cyclic subscript.

[0017] Furthermore, the specific calculation formula for calculating the propagating wave field using the three-dimensional elastic wave method in step two is as follows:

[0018]

[0019] Where ρ represents the density of the medium, and the vibration velocity of a particle in the forward propagation wave field along the three directions is v = [v x v y v z The normal stress and shear stress wave field of a particle in the propagating wave field in three directions are τ = [τ]. xx τ yy τ zz τ xy τ xz τ yz ]; Let t denote the time partial derivative operator, and t denote time. The link matrix is... The Hooke matrix related to the elasticity parameter is:

[0020] The specific formula for generating synthetic data is:

[0021] ψ syn (x r ,t;m)=ψ(x,t;m)δ(x r ),

[0022] Where x represents the computational region, δ(x) r ) represents the impulse function, and the synthesized data ψ(x,t;m)=[v τ] T At detector point x r Wave field data recorded at the location.

[0023] Furthermore, the objective function of least squares in step three is:

[0024]

[0025] Where the model parameters m = [α β] T α and β are the longitudinal and transverse wave velocities, respectively, and ψ syn Represents synthetic data, ψ obs Representing the observed data, the corresponding data residuals after calculating the above objective function are:

[0026]

[0027] Furthermore, the calculation formula for the backpropagating wavefield in step four is as follows:

[0028]

[0029] Among them, the reverse propagation wave field The vibrational velocities of a particle in the anti-propagating wave field along the three directions are υ = [υ x υ y υ z The normal stress and shear stress wavefield of a particle in the three directions in the reverse propagation wavefield are σ = [σ xx σ yy σ zz σ xy σ yz σ xz ].

[0030] Furthermore, in step five, the gradient of the model parameters is constructed using the forward and reverse propagation wavefields. The model parameters used are Lamé parameters, and the gradient of the Lamé parameters is expressed as follows:

[0031]

[0032]

[0033] Where λ and μ are Lamé parameters, and nt represents the number of time steps recorded. and Represents the spatial partial derivative operators with respect to the x, y, and z directions, respectively; the relationship between the Lamé parameters and the P-wave and S-wave velocities is λ+2μ=ρα. 2 and μ=ρβ 2 Let α and β be the P-wave and S-wave velocities, respectively. Calculate the gradient of the objective function with respect to the P-wave and S-wave velocities:

[0034]

[0035]

[0036] Furthermore, the specific process of calculating the inverse of the Hessian matrix using the L-BFGS method with second-order convergence in step six and continuously updating the initial velocity field model in step one is as follows:

[0037] The inverse matrix D of the Hessian matrix is ​​obtained using the L-BFGS method. k+1 The specific formula is as follows:

[0038]

[0039] Where I is a unit diagonal matrix, and y is the gradient change. k =g k+1 -g k The change s in the velocity field model k =m k+1 -mk , The gradient in this formula is obtained from step five; the inverse matrix D k+1 The initial D in the formula k Let k be the identity matrix, where the index k represents the iteration number. The model update is as follows:

[0040] m k+1 =m k -D k+1 g k+1 .

[0041] Using the above formula, the initial velocity field model is continuously updated until the data residual reaches a set threshold or the number of iterations reaches a preset value, at which point the iteration update stops and the final velocity field model parameters are obtained.

[0042] Compared with the prior art, the present invention has the following advantages:

[0043] 1. This invention first selects several discrete frequency points from the observation data, and then normalizes them by combining them with conventional band-limited wavelet sources and their locations to form a new source. This new wavelet source is then used for inversion, thereby overcoming the defect of insufficient low wavenumber updates caused by using conventional band-limited wavelet source data in the inversion method.

[0044] 2. This invention obtains the forward propagation wavefield and the reverse propagation wavefield by combining the new seismic source, and then constructs the gradient of the model parameters. By using several frequency groups from low to high in this way, multi-scale inversion can be effectively achieved, avoiding the problem of poor reference trace quality leading to unsatisfactory inversion results in convolutional multi-scale inversion.

[0045] 3. This invention uses the L-BFGS method combined with the acquired gradient to continuously update the initial velocity field model, and then uses the updated model for inversion. This method can further enhance the energy of low-frequency signals in the inversion, and enhance the contribution of low-frequency signals to the low wavenumber information of the model, thereby effectively reducing the dependence of the inversion on the initial velocity field model, thus improving the inversion accuracy and achieving the goal of high-precision inversion.

[0046] 4. The three-dimensional high-precision full waveform inversion scheme of the present invention, compared with the two-dimensional inversion imaging method, can obtain accurate spatial distribution of underground three-dimensional media and provide data guidance for the exploration and development of underground resources. Attached Figure Description

[0047] Figure 1 This is a spectral distribution diagram of the observation data in an embodiment of the present invention;

[0048] Figure 2 This is a spectrum diagram obtained through conventional inversion;

[0049] Figure 3 This is a spectrum diagram of an embodiment of the present invention;

[0050] Figure 4 This is the spectrum of a conventional band-limited wavelet source;

[0051] Figure 5 This is the spectrum diagram of the new seismic source formed by the present invention;

[0052] Figure 6 This is an image of the P-wave velocity results obtained through conventional inversion;

[0053] Figure 7 This is an image of the longitudinal wave velocity results obtained by the inversion of this invention;

[0054] Figure 8 This is an image of the shear wave velocity results obtained through conventional inversion;

[0055] Figure 9 This is an image of the transverse wave velocity results obtained by the inversion of this invention. Detailed Implementation

[0056] The present invention will be further described below.

[0057] The specific steps of this invention are as follows:

[0058] Step 1: Acquire observation data collected by the observation system. First, given an initial velocity field model, select several discrete frequency points from the observation data. Combine these with multiple frequency points from a conventional band-limited wavelet source and perform amplitude normalization to form a new source. The specific process is as follows: The specific formula for the new source is:

[0059]

[0060] Where F represents the Fourier transform, Indicates frequency ω i The inverse Fourier transform of the spectrum, x s Indicates the location of the earthquake source. Let ω represent the normalization of signal A, N represent the number of frequencies used, and ω = (ω1, ω2, ..., ω...). N ), s(x s ,t) is a conventional band-limited wavelet source, where t represents time and i represents a cyclic subscript.

[0061] Step 2: Using the new seismic source obtained in Step 1, calculate the propagating wavefield using the three-dimensional elastic wave method, and combine it with the wavefield data recorded by the observation system to obtain synthetic data; the specific calculation formula for the propagating wavefield is as follows:

[0062]

[0063] Where ρ represents the density of the medium, and the vibration velocity of a particle in the forward propagation wave field along the three directions is v = [v x v y v z The normal stress and shear stress wave field of a particle in the propagating wave field in three directions are τ = [τ]. xx τ yy τ zz τ xy τ xz τ yz ]; Let t denote the time partial derivative operator, and t denote time. The link matrix is... The Hooke matrix related to the elasticity parameter is:

[0064] The specific formula for generating synthetic data is:

[0065] ψ syn (x r ,t;m)=ψ(x,t;m)δ(x r ),

[0066] Where x represents the computational region, δ(x) r ) represents the impulse function, and the synthesized data ψ(x,t;m)=[v τ] T At detector point x r Wave field data recorded at the location.

[0067] Step 3: Based on the synthetic data obtained in Step 2, establish and calculate the least squares objective function. The calculated data residuals are used as the back-transmission source. The least squares objective function is:

[0068]

[0069] Where the model parameters m = [α β] T α and β are the longitudinal and transverse wave velocities, respectively, and ψ syn Represents synthetic data, ψ obs Representing the observed data, the corresponding data residuals after calculating the above objective function are:

[0070]

[0071] The result of the above calculation is used as the back-transmission source.

[0072] Step 4: Based on the anti-reverse source obtained in Step 3, calculate the anti-reverse wavefield under the anti-reverse source condition; the formula for calculating the anti-reverse wavefield is:

[0073]

[0074] Among them, the reverse propagation wave field The vibrational velocities of a particle in the anti-propagating wave field along the three directions are υ = [υ x υ y υ z The normal stress and shear stress wavefield of a particle in the three directions in the reverse propagation wavefield are σ = [σ xx σ yy σ zz σ xy σ yz σ xz ].

[0075] Step 5: Construct the gradient of the model parameters using the forward propagation wavefield obtained in Step 2 and the reverse propagation wavefield obtained in Step 4. Specifically, the model parameters used are Lamé parameters, and the gradient of the Lamé parameters is expressed as:

[0076]

[0077]

[0078] Where λ and μ are Lamé parameters, and nt represents the number of time steps recorded; and Represents the spatial partial derivative operators with respect to the x, y, and z directions, respectively; the relationship between the Lamé parameters and the P-wave and S-wave velocities is λ+2μ=ρα. 2 and μ=ρβ 2 Let α and β be the P-wave and S-wave velocities, respectively. Calculate the gradient of the objective function with respect to the P-wave and S-wave velocities:

[0079]

[0080]

[0081] Step Six: Based on the gradient obtained in Step Five, calculate the inverse of the Hessian matrix using the L-BFGS method, which has second-order convergence, and continuously update the initial velocity field model from Step One. The specific process is as follows:

[0082] The inverse matrix D of the Hessian matrix is ​​obtained using the L-BFGS method. k+1 The specific formula is as follows:

[0083]

[0084] Where I is a unit diagonal matrix, and y is the gradient change. k =g k+1 -g k , velocity model change s k =m k+1 -m k gradient of the objective function The gradient in this formula is obtained from step five; the inverse matrix Dk+1 The initial D in the formula k Let k be the identity matrix, where the index k represents the iteration number. The model update is as follows:

[0085] m k+1 =m k -D k+1 g k+1 .

[0086] Using the above formula, the initial velocity field model is continuously updated until the data residual reaches a set threshold or the number of iterations reaches a preset value. Then, the iteration updates are stopped to obtain the final velocity field model parameters. At this point, the velocity field model parameters are the three-dimensional high-precision imaging of the required area.

[0087] Experiments have shown that:

[0088] For a certain mining area, both conventional three-dimensional full waveform inversion and the three-dimensional full waveform inversion of this invention use the same observation data, such as... Figure 1 As shown, the collected observation data exhibits a dominance of dominant frequency energy, with weaker low-frequency signals. However, low-frequency data is crucial for overcoming the low-frequency dependence of full-waveform inversion and improving its accuracy. Therefore, effectively utilizing this low-frequency data to update model parameters in 3D inversion is critical. Figure 2 The diagram shown illustrates the spectral matching process of conventional 3D full-waveform inversion on observed data. This diagram demonstrates that even with multi-scale inversion strategies, existing methods often suffer from insufficient utilization of weak amplitude signals within each frequency band due to the dominant frequency characteristic of wavelets within each band. In contrast, this invention first generates a new seismic source and then performs the inversion. Figure 3 As shown, when the present invention combines new seismic sources for inversion, the selected frequency amplitudes in each inversion frequency band are used with equal weight. In particular, in low-frequency inversion, weak amplitude signals are effectively used, thereby ensuring low wavenumber updates of model parameters.

[0089] like Figure 4 As shown, this is the spectrum of a conventional band-limited wavelet source. Figure 5 The present invention selects several discrete frequency points from the collected observation data, and combines them with conventional band-limited wavelet sources and their locations for normalization processing to form a new source spectrum. Through these two figures, it can be seen that the spectrum of conventional band-limited wavelet sources has a wide frequency band and large amplitude differences at different frequencies, which is not conducive to the use of low-frequency information with weak amplitude. At the same time, wideband data is not conducive to multi-scale inversion.

[0090] like Figures 6 to 9As shown, when both conventional and the methods of this invention select poor initial models, conventional inversion methods cannot obtain high-precision inversion results, especially for shear wave inversion, which is more dependent on the initial model, where the inversion accuracy is significantly reduced. However, the method of this invention, by incorporating a new seismic source and calculating the forward and reverse propagation wave fields, obtains the gradient of the model parameters. Finally, it uses the L-BFGS method to continuously update the initial velocity field model using the obtained gradients. After the update, high-precision P- and S-wave inversion results are obtained. The inversion results are more clearly constructed, have stronger three-dimensionality, and show significant improvement, especially in deep inversion, resulting in higher inversion accuracy.

[0091] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A high-precision three-dimensional elastic wave full waveform inversion method, characterized in that, The specific steps are as follows: Step 1: Obtain the observation data collected by the observation system. First, give an initial velocity field model, select several discrete frequency points from the observation data, and combine them with multiple frequency points in the conventional band-limited wavelet source to perform amplitude normalization processing to form a new source. Step 2: Using the new seismic source obtained in Step 1, calculate the forward propagation wave field using the three-dimensional elastic wave method, and combine it with the wave field data recorded by the observation system to obtain synthetic data; Step 3: Establish a least squares objective function based on the synthetic data obtained in Step 2 and perform calculations. The calculated data residuals are used as the back-transmission source. Step 4: Calculate the anti-propagating wave field under the anti-propagating source condition based on the anti-propagating source information obtained in Step 3. Step 5: Construct the gradient of the model parameters using the forward propagation wave field obtained in Step 2 and the reverse propagation wave field obtained in Step 4; Step 6: Based on the gradient obtained in Step 5, the inverse matrix of the Hessian matrix is ​​calculated using the L-BFGS method with second-order convergence, and the initial velocity field model in Step 1 is continuously updated until the data residual reaches the set threshold or the number of iterations reaches the preset value. Then, the iteration update is stopped to obtain the final velocity field model parameters. The velocity field model parameters at this time are the three-dimensional high-precision imaging of the required area.

2. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 1, characterized in that, The specific formula for normalization to form a new seismic source in step one is as follows: Where F represents the Fourier transform, Indicates frequency ω i The inverse Fourier transform of the spectrum, x s Indicates the location of the earthquake source. Let ω represent the normalization of signal A, N represent the number of frequencies used, and ω = (ω1, ω2, ..., ω...). N ), s(x s ,t) is a conventional band-limited wavelet source, where t represents time and i represents a cyclic subscript.

3. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 2, characterized in that, The specific calculation formula for calculating the propagating wave field using the three-dimensional elastic wave method in step two is as follows: Where ρ represents the density of the medium, and the vibration velocity of a particle in the forward propagation wave field along the three directions is v = [v x v y v z The normal stress and shear stress wave field of a particle in the propagating wave field in three directions are τ = [τ]. xx τ yy τ zz τ xy τ xz τ yz ]; Let t denote the time partial derivative operator, and t denote time. The link matrix is... The Hooke matrix related to the elasticity parameter is: The specific formula for generating synthetic data is: ψ syn (x r ,t;m)=ψ(x,t;m)δ(x r ), Where x represents the computational region, δ(x) r ) represents the impulse function, and the synthesized data ψ(x,t;m)=[vτ] T At detector point x r Wave field data recorded at the location.

4. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 3, characterized in that, The objective function of least squares in step three is: Where, the model parameter m = [αβ] T α and β are the longitudinal and transverse wave velocities, respectively, and ψ syn Represents synthetic data, ψ obs Representing the observed data, the corresponding data residuals after calculating the above objective function are:

5. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 4, characterized in that, The formula for calculating the backpropagating wave field in step four is as follows: Among them, the reverse propagation wave field The vibrational velocities of a particle in the anti-propagating wave field along the three directions are υ = [υ x υ y υ z The normal stress and shear stress wavefield of a particle in the three directions in the reverse propagation wavefield are σ = [σ xx σ yy σ zz σ xy σ yz σ xz ].

6. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 5, characterized in that, In step five, the gradient of the model parameters is constructed using the forward propagation wavefield and the reverse propagation wavefield. The model parameters used are Lamé parameters, and the gradient of the Lamé parameters is expressed as follows: Where λ and μ are Lamé parameters, and nt represents the number of time steps recorded. and Represents the spatial partial derivative operators with respect to the x, y, and z directions, respectively; the relationship between the Lamé parameters and the P-wave and S-wave velocities is λ+2μ=ρα. 2 and μ=ρβ 2 Let α and β be the P-wave and S-wave velocities, respectively. Calculate the gradient of the objective function with respect to the P-wave and S-wave velocities:

7. The high-precision three-dimensional elastic wave full waveform inversion method according to claim 6, characterized in that, The specific process of calculating the inverse of the Hessian matrix and continuously updating the initial velocity field model in step one using the L-BFGS method with second-order convergence in step six is ​​as follows: The inverse matrix D of the Hessian matrix is ​​obtained using the L-BFGS method. k+1 The specific formula is as follows: Where I is a unit diagonal matrix, and y is the gradient change. k =g k+1 -g k , velocity model change s k =m k+1 -m k gradient of the objective function The gradient in this formula is obtained from step five; the inverse matrix D k+1 The initial D in the formula k Let k be the identity matrix, where the index k represents the iteration number. The model update is as follows: m k+1 =m k -D k+1 g k+1 Using the above formula, the initial velocity field model is continuously updated until the data residual reaches a set threshold or the number of iterations reaches a preset value, at which point the iteration update stops and the final velocity field model parameters are obtained.

Images