Viscous-acoustic medium preprocessing iteration method one-way wave migration method based on thin plate approximation

Abstract

The invention discloses a viscous sound medium preprocessing iteration method one-way wave migration method based on thin plate approximation. The method comprises the following steps: firstly, obtaining an actual seismic record; calculating a forward scattering migration wave field by using a viscous-acoustic medium preprocessing iteration method based on thin plate approximation, and calculating a back scattering migration wave field of the seismic record by taking observation data as a boundary value condition; and finally, scattering wave imaging is carried out by using a cross-correlation imaging condition, and a final migration imaging field is obtained. In the whole process, in order to avoid the influence of the instability of a strong disturbance medium background Green function on an imaging result, a small imaginary part attenuation component is introduced to improve the divergence problem of the Born series in a strong disturbance medium, and the calculation efficiency is further accelerated by utilizing a pre-solution factor. In order to avoid the problem of solving a large iterative matrix, a speed model is divided into thin plates by utilizing the thought of thin plate division, so that global calculation is converted into local calculation, and migration imaging with relatively high efficiency and relatively high precision is finally realized through the process.

Description

Technical Field

[0001] The invention relates to a geophysical exploration method, in particular to a one-way wave migration method based on a thin plate approximation viscoacoustic medium preprocessing iteration method. Background Art

[0002] In contrast to differential operators, integral operators typically utilize perturbation theory to decompose the medium into a superposition of background and perturbation parameters, transforming the scattering problem into a boundary value problem for solving the second-kind Fredholm integral equation. Solving this problem using Green's function yields the Lippmann-Schwinger (LS) integral equation. Due to the nonlinearity of the LS equation, an iterative method can be used to obtain the Born scattering series. However, the traditional Born series converges slowly in strongly perturbed media, and the first-order Born approximation suffers from relatively large computational errors in strongly perturbed media. Therefore, the Born approximation is only applicable to weakly transversely inhomogeneous media (velocity perturbations less than 10%). To address this issue, industry researchers are currently introducing the renormalization concept from quantum scattering theory. By renormalizing the Born series, they can improve the convergence of the scattering series in strongly perturbed media, thereby obtaining a more accurate numerical solution with fewer iterations. In the field of seismic scattering, the renormalization method arranges the scattering series into a series of subseries by splitting the operator and summing up some of the subseries to eliminate the divergent terms of the series.

[0003] The purpose of the preconditioned iteration method is to improve the convergence of the Born series, which is similar to the preconditioned method for solving a system of linear equations. In order to obtain a converged Born series, a small imaginary component is first introduced into the background wavenumber, thereby obtaining an LS equation with a complex wavenumber. The background Green's function corresponding to the LS equation has a damping factor. Then, the LS equation with the complex wavenumber is preconditioned using a pre-solution factor defined by a small imaginary component. Finally, the preconditioned LS equation is iteratively solved using an iterative method to obtain the Born series improved by the complex wavenumber and preconditioning. Through experimental verification, the Born series improved in this way also converges for strong scattering. However, due to different preconditioning factors, the convergence speed of the Born series is different.

[0004] Therefore, how to provide a new migration imaging method, by selecting appropriate preprocessing factors and applying the above-mentioned preprocessing iterative method to viscoacoustic medium one-way wave migration imaging, and ultimately achieve migration imaging with high efficiency and good accuracy, is the research direction required by the present invention. Summary of the Invention

[0005] In response to the problems existing in the above-mentioned prior art, the present invention provides a one-way wave migration method of viscoacoustic medium preprocessing iteration method based on thin plate approximation, which can effectively solve the above-mentioned technical problems.

[0006] To achieve the above object, the present invention adopts a technical solution: a viscoacoustic medium preprocessing iterative method one-way wave migration method based on thin plate approximation, comprising the following steps: Step 1: Arrange multiple detectors in a row at equal intervals in the required detection area to form an observation system, and arrange multiple seismic sources in a row at equal intervals and number them.

[0007] Step 2: Stimulate each earthquake source in sequence according to the number, and the observation system obtains the earthquake record when each earthquake source is excited.

[0008] Step 3: Set the velocity model parameters of the seismic scattered wave propagating in the detection area, and determine the absorption and attenuation model of the seismic scattered wave when propagating in the detection area.

[0009] Step 4: Using the idea of ​​thin plate division, the absorption attenuation model of step 3 is divided into multiple thin plates along the depth direction, and the velocity parameters and quality factor parameters in each thin plate are determined.

[0010] Step 5: Load the Ricker wavelet at the first source position and convert it into the frequency wavenumber domain using Fourier transform. Use it as the initial value for calculation to obtain the forward scattering migration wave field of all thin plates.

[0011] Step 6: Perform Fourier transform on the seismic record corresponding to the first source excitation, convert it into the frequency-wavenumber domain, and use it as the initial value to perform calculations using the calculation process of step 5 to obtain the corresponding backscattering migration wave field.

[0012] Step 7: Use cross-correlation imaging conditions to perform migration imaging on the forward scattering migration wavefield and the backscattering migration wavefield obtained in steps 5 and 6 to obtain the imaging result of the first earthquake source.

[0013] Step 8: Repeat steps 5 to 7 for each earthquake source in order of earthquake source number, thereby obtaining imaging results for each earthquake source.

[0014] Step 9: Superimpose the imaging results of all earthquake sources to finally obtain the offset imaging results of the detection area.

[0015] Furthermore, in step 1, a coordinate system is established along the arrangement direction of each detector and source as the X axis and perpendicular to the arrangement direction (i.e., the depth direction of the detection area) as the 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 number of sources is ns, the distance between adjacent sources is ds, and the position of the first source is (0.0, 0.0), the main frequency of the wavelet is fHz, the sampling interval is set to 1ms, and the sampling time length is nt.

[0016] Furthermore, the velocity model parameters in step 3 are ,in 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 parameters , the specific calculation formula is as follows: (1) in The unit is km / s; the absorption attenuation model is the Kolsky-Futterman model, and its compensation formula is expressed as follows: (2) in is the complex velocity, is the reference frequency The corresponding approximate speed.

[0017] Furthermore, the step five is specifically as follows: A. Decompose the velocity in each thin plate into background velocity and disturbance velocity, i.e. ,in is the background velocity in the thin plate; is the perturbation velocity in the thin plate; then the quality factor parameter is decomposed into the background parameter and the disturbance parameter ,Right now .

[0018] B. Solve for the imaginary component in the thin plate and pre-solution factors .

[0019] C. As the boundary condition, 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: (3) in is the transverse wave number in the x direction, is the grid spacing, is the vertical wave number, .

[0020] D. Convert to frequency space domain Perform phase shift correction and compensation, and the extension formula is: (4) in: .

[0021] E. The forward scattering migration wave field calculation formula of the j+1th thin plate is as follows: (5) in is the background wave field, is the pre-solution factor in the j+1th thin plate, , , , is the background Green function.

[0022] Substitute the result obtained in step D into formula (5) to obtain the forward scattering migration wave field after correction of the j+1th thin plate: (6) F. Use the forward scattering offset wave field calculated by formula (6) as the boundary condition of the next thin plate, repeat steps A to E, and calculate the forward scattering offset wave field of the next thin plate. ; Repeat this process until the forward scattering offset wave fields of all thin plates are obtained.

[0023] Furthermore, the calculation formula for the cross-correlation imaging in step seven is: (7) Where * is the conjugate and W is the cutoff frequency.

[0024] Compared with the existing technology, the present invention first obtains the actual seismic record; secondly, it uses the thin-plate approximation viscoacoustic medium preprocessing iterative method extension operator to calculate the forward scattering migration wave field; then, the observation data is used as the boundary condition to realize the backscattering migration wave field calculation of the seismic record; finally, the scattered wave imaging is performed using the cross-correlation imaging condition to obtain the final migration imaging field. Throughout this process, to avoid the influence of the instability of the background Green function of the strongly disturbed medium on the imaging results, a small imaginary attenuation component is introduced to improve the divergence of the Born series in the strongly disturbed medium, and the pre-solution factor is used to further accelerate the computational efficiency. To avoid the problem of solving large iterative matrices, the velocity model is divided into thin plates using the thin-plate partitioning concept, thereby converting the global calculation into a local calculation, and then constructing the thin-plate approximation viscoacoustic medium preprocessing iterative method extension operator. In addition, considering the complexity of the actual medium, the actual model is equivalent to the Kolsky-Futterman model, which is more in line with the actual situation. Through the above process, efficient and accurate migration imaging is ultimately achieved. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 It is the overall flow chart of the present invention.

[0026] Figure 2 are velocity model parameters and quality factor parameters of the embodiment of the present invention.

[0027] Figure 3 This is the migration imaging result of the embodiment of the present invention. DETAILED DESCRIPTION

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

[0029] like Figure 1 As shown, the present invention includes the following steps: Step 1: Arrange multiple detectors in a row at equal intervals in the required detection area to form an observation system, and arrange multiple seismic sources in a row at equal intervals and number them; establish a coordinate system along the arrangement direction of each detector and seismic source as the X-axis and perpendicular to the arrangement direction as the 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 number of seismic sources is ns, the distance between adjacent seismic sources is ds, and the position of the first seismic source is (0.0, 0.0); the main frequency of the wavelet is fHz, the sampling interval is set to 1ms, and the sampling time length is nt.

[0030] Step 2: Stimulate each source in sequence according to the number. The observation system obtains the earthquake records when each source is excited. .

[0031] Step 3: Set the velocity model parameters of the seismic scattered wave propagating in the detection area, and determine the absorption and attenuation model of the seismic scattered wave when propagating in the detection area, such as Figure 2 As shown, specifically: the velocity model parameters are ,in 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 parameters , the specific calculation formula is as follows: (1) in The unit is km / s; the absorption attenuation model is the Kolsky-Futterman model, and its compensation formula is expressed as follows: (2) in is the complex velocity, is the reference frequency The corresponding approximate speed.

[0032] Step 4: Using the idea of ​​thin plate division, the absorption attenuation model in step 3 is divided into multiple thin plates along the depth direction, a total of j thin plates, with a thickness of dp, which is usually the thickness of a grid dz; and represent the velocity parameter and quality factor parameter in a thin plate respectively.

[0033] Step 5: Load the Ricker wavelet at the first source location (t is the time parameter), and is converted to the frequency wavenumber domain using Fourier transform ( is the circular frequency), and the forward scattering migration wave field of all thin plates is obtained after calculation using it as the initial value, specifically: A. Decompose the velocity in each thin plate into background velocity and disturbance velocity, i.e. ,in is the background velocity in the thin plate; is the perturbation velocity in the thin plate; then the quality factor parameter is decomposed into the background parameter and the disturbance parameter ,Right now .

[0034] B. Solve for the imaginary component in the thin plate and pre-solution factors .

[0035] C. As the boundary condition, 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: (3) in is the transverse wave number in the x direction, is the grid spacing, is the vertical wave number, .

[0036] D. Convert to frequency space domain Perform phase shift correction and compensation, and the extension formula is: (4) in: .

[0037] E. The forward scattering migration wave field calculation formula of the j+1th thin plate is as follows: (5) in is the background wave field, is the pre-solution factor in the j+1th thin plate, , , , is the background Green function.

[0038] Substitute the result obtained in step D into formula (5) to obtain the forward scattering migration wave field after correction of the j+1th thin plate: (6).

[0039] F. Use the forward scattering offset wave field calculated by formula (6) as the boundary condition of the next thin plate, repeat steps A to E, and calculate the forward scattering offset wave field of the next thin plate. ; Repeat this process until the maximum depth of the model is reached and the forward scattering migration wave fields of all thin plates are obtained.

[0040] Step 6: Perform Fourier transform on the seismic record corresponding to the first source excitation and convert it into the frequency-wavenumber domain , and use it as the initial value to calculate using the calculation process of step 5 to obtain the corresponding backscattering migration wave field .

[0041] Step 7: Use the cross-correlation imaging condition to perform migration imaging on the forward scattering migration wavefield and the reverse scattering migration wavefield obtained in steps 5 and 6 to obtain the imaging result of the first earthquake source. ; The calculation formula of cross-correlation imaging is: (7) Where * is the conjugate and W is the cutoff frequency.

[0042] Step 8: Repeat steps 5 to 7 for each earthquake source in sequence according to the earthquake source number, so as to obtain the imaging results of each earthquake source, which are 、 … .

[0043] Step 9: Superimpose the imaging results of all earthquake sources, such as Figure 3 As shown, the final offset imaging result of the detection area is obtained , the specific formula is: (8).

[0044] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements and modifications can be made without departing from the principles of the present invention. These improvements and modifications should also be regarded as the scope of protection of the present invention.

Claims

1. A one-way wave migration method based on thin plate approximation and viscoacoustic medium preprocessing iterative method, characterized in that: The following steps are involved: Step 1: Arrange multiple geophones in a row at equal intervals in the desired detection area to form an observation system, and arrange multiple seismic sources in a row at equal intervals and number them; Step 2: Stimulate each seismic source in sequence according to the number, and the observation system obtains the earthquake record when each seismic source is excited; Step 3: Set the velocity model parameters of the seismic scattered wave propagating in the detection area, and determine the absorption and attenuation model of the seismic scattered wave when propagating in the detection area; Step 4: Using the thin plate partitioning concept, the absorption attenuation model of step 3 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 5: Load the Ricker wavelet at the first source location and convert it into the frequency wavenumber domain using Fourier transform. Use it as the initial value to calculate and obtain the forward scattering migration wavefield of all thin plates. Step 6: Perform Fourier transform on the seismic record corresponding to the first source excitation to convert it into the frequency-wavenumber domain, and use it as the initial value to perform calculations using the calculation process of step 5 to obtain the corresponding backscattering migration wave field; Step 7: Using cross-correlation imaging conditions, perform migration imaging on the forward scattering migration wavefield and the reverse scattering migration wavefield obtained in steps 5 and 6 to obtain an imaging result of the first earthquake source; Step 8: Repeat steps 5 to 7 for each earthquake source in order of earthquake source number, thereby obtaining imaging results for each earthquake source; Step 9: Superimpose the imaging results of all earthquake sources to finally obtain the offset imaging results of the detection area.

2. The one-way wave migration method of viscoacoustic medium preprocessing iteration method based on thin plate approximation according to claim 1 is characterized in that: In the step 1, a coordinate system is established along the X-axis along the layout direction of each detector and source and perpendicular to the layout direction as the 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 number of sources is ns, the distance between adjacent sources is ds, the position of the first source is (0.0, 0.0), the main frequency of the sub-wave is fHz, the sampling interval is set to 1ms, and the sampling time length is nt.

3. The one-way wave migration method of viscoacoustic medium preprocessing iteration method based on thin plate approximation according to claim 1 is characterized in that: The velocity model parameters in step 3 are: ,in 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 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 is the complex velocity, is the reference frequency The corresponding approximate speed.

4. The one-way wave migration method of viscoacoustic medium preprocessing iteration method based on thin plate approximation according to claim 1 is characterized in that: The step five is specifically as follows: A. Decompose the velocity in each thin plate into background velocity and disturbance velocity, i.e. ,in is the background velocity in the thin plate; is the perturbation velocity in the thin plate; then the quality factor parameter is decomposed into the background parameter and the disturbance parameter ,Right now ; B. Solve for the imaginary component in the thin plate and pre-solution factors ; C. As the boundary condition, 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: （3） in is the transverse wave number in the x direction, is the grid spacing, is the vertical wave number, ; D. Convert to frequency space domain Perform phase shift correction and compensation, and the extension formula is: （4） in: ; E. The forward scattering migration wave field calculation formula of the j+1th thin plate is as follows: （5） in is the background wave field, is the pre-solution factor in the j+1th thin plate, , , , is the background Green function; Substitute the result obtained in step D into formula (5) to obtain the forward scattering migration wave field after correction of the j+1th thin plate: （6） F. Use the forward scattering offset wave field calculated by formula (6) as the boundary condition of the next thin plate, repeat steps A to E, and calculate the forward scattering offset wave field of the next thin plate. ; Repeat this process until the forward scattering offset wave fields of all thin plates are obtained.

5. The one-way wave migration method of viscoacoustic medium preprocessing iteration method based on thin plate approximation according to claim 1 is characterized in that: The calculation formula for the cross-correlation imaging in step seven is: （7） Where * is the conjugate and W is the cutoff frequency.

Images