A numerical simulation method for scattering wave of viscoacoustic medium based on preconditioned iterative method

Abstract

The application discloses a kind of based on pretreatment iterative method's viscous medium scattering wave numerical simulation method, first in background wave number Small virtual component is introduced into, and then a L-S equation with complex wave number is obtained;Then, a pre-solution factor defined with the aid of the small virtual component is used to pretreat the L-S equation with complex wave number;Finally, the L-S equation after pretreatment is iteratively solved using iterative method, and the Born series improved by complex wave number and pretreatment is obtained;During the whole process, a small virtual component is introduced to improve the divergence problem of Born series in strong disturbance medium, and the calculation efficiency is further accelerated using the pre-solution factor;And in order to avoid solving large iterative matrix, the velocity model is divided into thin plates using the idea of thin plate division, so as to convert global calculation into local calculation, and through the above process, the seismic scattering wave record with high efficiency and good accuracy is finally extracted.

Description

TECHNICAL FIELD

[0001] The present application relates to a geophysical exploration method, in particular to a visco-acoustic medium scattering wave numerical simulation method based on a pre-processing iteration method. BACKGROUND

[0002] When the geological body (visco-acoustic medium) of the exploration target has a complex structure, faults are developed, the stratum inclination is relatively steep, or the lithology is laterally suddenly changed, or different scale non-uniform geological bodies coexist, a very complex seismic wave field with multiple wave groups interfering with each other is formed. In this case, only using reflection wave and diffraction wave information cannot achieve fine imaging in the complex region. Research shows that scattering waves also carry geometric and physical information related to complex structure and complex lithology. Therefore, it is very important to use scattering waves to realize complex structure migration imaging while using conventional reflection waves to realize imaging.

[0003] In order to describe the propagation law of scattered waves, the medium is divided into background medium and perturbed medium by using scattering theory, and the boundary value problem of Helmholtz equation describing the scattering problem is converted into a problem of solving the second Fredholm integral equation by using representation theorem. Since the basic idea of deriving this integral equation is derived from quantum scattering theory, this equation is referred to as Lippmann-Schwinger equation, abbreviated as L-S equation. In this way, the seismic wave scattering problem is converted into the problem of solving the L-S equation. The classical Born series exists slow convergence and is easy to diverge in strong perturbed medium. In order to solve this problem, Osnabruggle et al. published A convergent Born series for solving the inhomogeneous Helmholtz equation in arbitrarily large media in Journal of Computational Physics in 2016. Specifically, the convergence of Born scattering series is improved by using the background Green's function with damping factor (obtained by introducing virtual part in the background wave number) and adopting the pre-solution pre-processing method. In 2020, Huang et al. published the paper On the applicability of a normalized Born series for seismic wavefield modelling in strongly scattering media in Journal of Geophysics and Engineering. The physical significance of the modified Born series proposed by Osnabrugge et al. is reinterpreted from the perspective of renormalization, and the applicability of the method in the seismic wave field forward modeling in strong perturbed scattering medium is verified by numerical examples. However, since the geological body is visco-acoustic medium, it has different physical properties from general acoustic medium, so the method cannot be directly applied to scattered wave migration imaging, and the calculation amount of its global pre-solution factor is relatively large, so the method is not suitable for scattered wave numerical simulation.

[0004] Therefore, how to provide a new scattered wave numerical simulation method, which can quickly and accurately simulate scattered wave in visco-acoustic medium, and further reveal the propagation law of scattered wave in visco-acoustic medium, is the research direction of the present application. SUMMARY

[0005] In view of the problems in the prior art, the present application provides a visco-acoustic medium scattering wave numerical simulation method based on a preconditioned iterative method, which uses a visco-acoustic medium continuation operator based on the preconditioned iterative method to numerically simulate the scattering wave, so that the scattering wave numerical simulation is quickly and accurately realized under the visco-acoustic medium condition, and the scattering wave propagation law in the visco-acoustic medium is revealed.

[0006] In order to achieve the above-mentioned purpose, the technical scheme adopted by the present application is as follows: a visco-acoustic medium scattering wave numerical simulation method based on a preconditioned iterative method, comprising the following steps:

[0007] Step one, a plurality of geophones are arranged in the simulation environment to form an observation system, and the source position and the main frequency of the excited seismic wave are determined, a seismic wave is excited at the source position, and corresponding seismic records are obtained; the seismic records include reflection wave, scattering wave and direct wave data.

[0008] Step two, the velocity model parameters of the seismic scattering wave propagation in the detection area are set, and the absorption and attenuation model of the seismic scattering wave propagation in the detection area is determined.

[0009] Step three, the absorption and attenuation model of step two is divided into a plurality of thin plates along the depth direction by using the thin plate division idea, and the velocity parameters and the quality factor parameters in each thin plate are determined.

[0010] Step four, a Ricker wavelet is loaded at the source position, and Fourier transform is used to convert it to the frequency-wavenumber domain, which is used as an initial value to calculate the downgoing scattering wave field and the first upgoing scattering wave field of the minimum depth thin plate, and the first upgoing scattering wave field is stored; then the downgoing scattering wave field value of the minimum depth thin plate is used as the edge value condition of the adjacent lower thin plate to continue the calculation, and the downgoing scattering wave field and the first upgoing scattering wave field of the adjacent lower thin plate are obtained, and the calculation is repeated until the bottom of the model, and the downgoing scattering wave field and the first upgoing scattering wave field of all thin plates are obtained.

[0011] Step five, the first upgoing scattering wave field of the maximum depth thin plate is used as the edge value condition of the adjacent upper thin plate, the calculation process of step five is repeated, the upgoing scattering wave field of the adjacent upper thin plate is obtained, and the first upgoing scattering wave field of the thin plate calculated in step four is added as the edge value condition of the adjacent upper thin plate, and the calculation process of step five is continued, and the calculation is repeated until the upgoing scattering wave field calculation of the minimum depth thin plate is completed.

[0012] Step six, the frequency domain upgoing scattering wave field of the minimum depth thin plate is converted into the time domain, so that all seismic scattering wave records of each geophone are extracted.

[0013] Further, the step one establishes a coordinate system along the arrangement direction of each detector and the source as X axis and perpendicular to the arrangement direction (i.e. the depth direction) as Z axis, with the coordinate origin at the position (0.0, 0.0); the number of detectors is nx, the distance between adjacent detectors is dr, and the position of the first detector is (0.0, 0.0); the position of the source is (nx / 2, 0.0), the main frequency of the seismic wave is f Hz, the sampling interval is set to l ms, and the sampling time length is nt.

[0014] Further, the velocity model parameter in the step two is , wherein is the coordinate parameter, the grid parameter of X axis is nx, the grid parameter of Z axis is nz, and the grid spacing is dx and dz, respectively; the quality factor parameter is calculated based on the velocity model parameter , and the specific calculation formula is as follows:

[0015] (1)

[0016] The model of absorption attenuation is Kolsky-Futterman model, and the expression of the compensation formula is:

[0017] (2)

[0018] , wherein is the complex velocity, is the reference frequency , and the corresponding velocity approximation value is

[0019] Further, the step four is specifically:

[0020] A. The velocity in each thin plate is decomposed into background velocity and perturbation velocity, i.e. , wherein is the background velocity in the thin plate; is the perturbation velocity in the thin plate; then the quality factor parameter is decomposed into background parameter and perturbation parameter , i.e. .

[0021] B. The imaginary part and the pre-decomposition factor in the thin plate are solved.

[0022] C. is the boundary condition, i.e. , the phase shift compensation processing is performed in the frequency-wavenumber domain to obtain the background wave field , and the calculation formula of the background wave field is:

[0023] (3)

[0024] where is the transverse wavenumber in x direction, is the grid interval, is the vertical wavenumber, .

[0025] D. Transform to frequency space domain Phase shift correction and compensation are performed, and the continuation formula is:

[0026] (4)

[0027] where: .

[0028] E. The down-going scattered wave field of the depth-minimum thin plate is calculated according to the following formula:

[0029] (5)

[0030] where is the background wave field, is the pre-decomposition factor in the depth-minimum thin plate, , , , is the background Green function, , is the background complex velocity.

[0031] F. The first up-going scattered wave field of the depth-minimum thin plate is calculated according to the following formula:

[0032] (6)

[0033] G. The down-going scattered wave field of the depth-minimum thin plate is taken as the boundary value condition of the adjacent lower thin plate, and steps A to F are repeated to calculate the down-going scattered wave field and the first up-going scattered wave field of the adjacent lower thin plate and ; this is repeated until the bottom of the model, and the down-going scattered wave field and the first up-going scattered wave field of all the thin plates are obtained.

[0034] Further, when calculating the up-going scattered wave field of each thin plate from the depth maximum to the depth minimum in step five, all the calculated formulas in step four need to be multiplied by the negative sign before the imaginary part i and then calculated, so as to realize the calculation of the up-going scattered wave field in the direction of the depth minimum.

[0035] Compared with the prior art, the application firstly introduces a small imaginary part into the background wave number, and then obtains an L-S equation with a complex wave number, so that the background Green function corresponding to the L-S equation has a damping factor; then, a pre-decomposition factor defined by means of the small imaginary part is used to preprocess the L-S equation with the complex wave number; finally, the L-S equation after preprocessing is solved by using an iterative method, and a Born series improved by the complex wave number and preprocessing is obtained; in the whole process, in order to avoid the influence of the instability of the background Green function of the strong disturbance medium on the extraction result of the scattered wave, a small imaginary part attenuation component is introduced to improve the divergence problem of the Born series in the strong disturbance medium, and the pre-decomposition factor is used to further speed up the calculation efficiency. In order to avoid solving a large iterative matrix, the velocity model is divided into a plurality of thin plates by using the thin plate division idea, so that the global calculation is converted into local calculation, and then a viscous medium preprocessing iterative method continuation operator based on the thin plate approximation is constructed, and in addition, considering the complexity of the actual medium, the Kolsky-Futterman model is used to equivalent the actual model, which is more in line with the actual situation. Through the above process, the seismic scattered wave record is finally extracted with high efficiency and good precision. BRIEF DESCRIPTION OF DRAWINGS

[0036] Figure 1 is the overall flowchart of the application.

[0037] Figure 2 is the velocity model parameter and the quality factor parameter of the embodiment of the application.

[0038] Figure 3 is the seismic scattered wave record extracted by the embodiment of the application. DETAILED DESCRIPTION

[0039] The application will be further described below.

[0040] As Figure 1 shown, the application comprises the following steps:

[0041] Step one, a plurality of geophones are arranged at equal intervals in a row to form an observation system in a simulation environment, and the position of a seismic source and the main frequency of a seismic wave excited are determined, specifically: a coordinate system is established along the arrangement direction of the geophones and the seismic source as the X axis and perpendicular to the arrangement direction (i.e. the depth direction) as the Z axis, and the coordinate origin is at the (0.0, 0.0) position; the number of geophones is nx, the distance between adjacent geophones is dr, and the position of the first geophone is (0.0, 0.0); the position of the seismic source is (nx / 2, 0.0), the main frequency of the seismic wave is fHz, the sampling interval is set to lms, and the sampling time length is nt; a seismic wave is excited at the position of the seismic source, and the corresponding seismic record is obtained; the seismic record includes reflected wave, scattered wave and direct wave data.

[0042] Step two, set the velocity model parameters of the seismic scattered wave propagation in the detection area as , wherein is the coordinate parameter, the grid parameter of the X axis is nx, the grid parameter of the Z axis is nz, and the grid spacing is dx and dz, respectively; the quality factor parameter is calculated based on the velocity model parameter , and the specific calculation formula is as follows:

[0043] (1)

[0044] The model of absorption attenuation is Kolsky-Futterman model, and the expression of the compensation formula is:

[0045] (2)

[0046] , wherein is the complex velocity, is the reference frequency , and the corresponding velocity approximation value is.

[0047] Step three, the absorption attenuation model of step two is divided into multiple thin plates along the depth direction by adopting the idea of thin plate division, and the velocity parameter and the quality factor parameter in each thin plate are determined as shown in Figure 2 .

[0048] Step four, load the Ricker wavelet at the source position, and use Fourier transform to convert to the frequency wave number domain, which is used as the initial value to calculate the downgoing scattered wave field and the first upgoing scattered wave field of the minimum depth thin plate, and the first upgoing scattered wave field is stored; then the downgoing scattered wave field value of the minimum depth thin plate is used as the edge value condition of the adjacent lower thin plate to continue the calculation, and the downgoing scattered wave field and the first upgoing scattered wave field of the adjacent lower thin plate are obtained, and so on until the bottom of the model, the downgoing scattered wave field and the first upgoing scattered wave field of all thin plates are obtained, and the specific process is as follows:

[0049] A, the velocity in each thin plate is decomposed into background velocity and perturbation velocity, that is, , wherein is the background velocity in the thin plate; is the perturbation velocity in the thin plate; then the quality factor parameter is decomposed into background parameter and perturbation parameter , that is, .

[0050] B, the imaginary part and the pre-decomposition factor in the thin plate are solved.

[0051] C, The boundary condition is The phase shift compensation is performed in the frequency-wavenumber domain to obtain the background wave field The calculation formula of the background wave field is

[0052] (3)

[0053] where is the transverse wavenumber in the x direction, is the grid spacing, is the vertical wavenumber, .

[0054] D. The phase shift correction and compensation are performed in the frequency space domain The continuation formula is

[0055] (4)

[0056] where: .

[0057] E. The calculation formula of the downgoing scattered wave field of the depth-minimum thin plate is as follows:

[0058] (5)

[0059] where is the background wave field, is the pre-decomposition factor in the depth-minimum thin plate, , , , is the background Green function, , is the background complex velocity.

[0060] F. The calculation formula of the first upgoing scattered wave field of the depth-minimum thin plate is as follows:

[0061] (6)

[0062] G. The downgoing scattered wave field of the depth-minimum thin plate is taken as the boundary condition of the adjacent lower thin plate, and steps A to E are repeated to calculate the downgoing scattered wave field and the first upgoing scattered wave field of the adjacent lower thin plate; this is repeated until the bottom of the model, and the downgoing scattered wave field and the first upgoing scattered wave field of all the thin plates are obtained.

[0063] Step five, taking the first up-scattered wave field of the maximum depth plate as the edge value condition of the adjacent upper plate, repeating the calculation process of step five, obtaining the up-scattered wave field of the adjacent upper plate, and adding the first up-scattered wave field of the plate calculated in step four as the edge value condition of the adjacent upper plate, and continuing the calculation process of step five, so repeating (i.e. each time the up-scattered wave field of the upper plate is calculated, the first up-scattered wave field calculated before is added as the edge value condition of the adjacent upper plate for subsequent calculation), until the up-scattered wave field calculation of the minimum depth plate is completed; when calculating the up-scattered wave field of each plate, the imaginary part i of all the calculation formulas in step four needs to be multiplied by a negative sign, so as to realize the up-scattered wave field calculation in the direction of the minimum depth.

[0064] Step six, converting the frequency domain up-scattered wave field of the minimum depth plate into the time domain, as shown in Figure 3 , thereby extracting all the seismic scattered wave records of each geophone.

[0065] The above only describes the preferred embodiments of the present application, and it should be noted that for those skilled in the art, without departing from the principles of the present application, a number of improvements and refinements can be made, and these improvements and refinements should also be considered within the scope of protection of the present application.

Claims

1. A numerical simulation method for scattering waves in viscous acoustic media based on a preprocessing iterative method, characterized in that, Includes the following steps: Step 1: In a simulated environment, multiple geophones are arranged in a row at equal intervals to form an observation system. The location of the earthquake source and the dominant frequency of the excitation seismic wave are determined. A seismic wave is excited at the earthquake source location, and the corresponding earthquake record is obtained. Step 2: Set the velocity model parameters for the propagation of seismic scattered waves in the detection area, and determine the absorption and attenuation model for the seismic scattered waves as they propagate in the detection area; Step 3: Using the idea of ​​thin plate division, the absorption attenuation model in Step 2 is divided into multiple thin plates along the depth direction, and the velocity parameters and quality factor parameters in each thin plate are determined. Step 4: Load the Ricker wavelet at the source location and transform it to the frequency-wavenumber domain using Fourier transform. Use this as the initial value for calculation to obtain the downscattered wave field and the first upscattered wave field of the thinnest plate with the smallest depth, and store the first upscattered wave field. Then, use the downscattered wave field value of the thinnest plate with the smallest depth as the boundary condition of the adjacent thin plates below to continue the calculation, and obtain the downscattered wave field and the first upscattered wave field of the adjacent thin plates below. Repeat this process until the bottom of the model to obtain the downscattered wave field and the first upscattered wave field of all thin plates. Step 5: Use the first upward scattered wave field of the thin plate with the maximum depth as the boundary condition of the adjacent thin plate above. Repeat the calculation process of Step 5 to obtain the upward scattered wave field of the adjacent thin plate above. Add the first upward scattered wave field of the thin plate calculated in Step 4 as the boundary condition of the adjacent thin plate above. Continue the calculation process of Step 5. Repeat this process until the calculation of the upward scattered wave field of the thin plate with the minimum depth is completed. Step 6: Convert the frequency domain upscattered wave field of the thinnest plate with the smallest depth into the time domain, thereby extracting all seismic scattered wave records from each detector.

2. The numerical simulation method for scattering waves in viscous acoustic media based on the preprocessing iterative method according to claim 1, characterized in that, In step one, a coordinate system is established with the X-axis along the direction of each detector and seismic source and the Z-axis perpendicular to the direction of deployment. The origin of the coordinate system is at (0.0, 0.0). The number of detectors is nx, the distance between adjacent detectors is dr, and the position of the first detector is (0.0, 0.0). The location of the seismic source is (nx / 2, 0.0), the dominant frequency of the seismic wave is fHz, the sampling interval is set to 1ms, and the sampling time length is nt.

3. The numerical simulation method for scattering waves in viscous acoustic media based on the preprocessing iterative method according to claim 1, characterized in that, The velocity model parameters in step two are: ,in The coordinate parameters are defined as follows: the grid parameter for the X-axis is nx, the grid parameter for the Z-axis is nz, and the grid spacings are dx and dz, respectively. The quality factor parameters are calculated based on the velocity model parameters. The specific calculation formula is as follows: （1） The absorption attenuation model is the Kolsky-Futterman model, and its compensation formula is expressed as follows: （2） in For complex velocity, Reference frequency The corresponding approximate speed value.

4. The numerical simulation method for scattering waves in viscous acoustic media based on the preprocessing iterative method according to claim 1, characterized in that, Step four specifically involves: A. Within each thin plate, the velocity is decomposed into background velocity and disturbance velocity, i.e. ,in The background velocity within the thin plate; The disturbance velocity within the thin plate is given; then the quality factor parameter is decomposed into background parameters. and disturbance parameters ,Right now ; B. Solving for the imaginary components within the thin plate. and pre-solution factor ; C These are boundary conditions, i.e. Phase shift compensation is performed in the frequency wavenumber domain to obtain the background wavefield. The formula for calculating the background wave field is: （3） in Let x be the transverse wavenumber in the x-direction. For grid spacing, For vertical wavenumber, ; D. Transform to the frequency space domain For phase shift correction and compensation, the continuation formula is: （4） in: ; E. The formula for calculating the downward scattered wave field of the thinnest plate with the minimum depth is as follows: （5） in For the background wave field, The pre-solution factor within the thinnest plate with the minimum depth. , , , For the background Green function, , Background complex velocity; F. The formula for calculating the first-order upward scattered wave field of the thinnest plate with the minimum depth is as follows: （6） G. Using the downward scattered wave field of the thinnest plate with the smallest depth as the boundary condition for the adjacent thin plate below, repeat steps A to F to calculate the downward scattered wave field of the adjacent thin plate below. and a first upscattered wave field Repeat this process until the bottom of the model is reached, obtaining the downscattered wave field and the first upscattered wave field of all the thin plates.

5. The numerical simulation method for scattering waves in viscous acoustic media based on the preprocessing iterative method according to claim 4, characterized in that, When calculating the upward scattering wave field of each thin plate in step five, from the direction of maximum depth to minimum depth, it is necessary to multiply the imaginary part i in all the calculation formulas in step four by a negative sign before performing the calculation.

Images