A two-way wave migration imaging method based on a Muller model and a VTI medium

Abstract

The application discloses a two-way wave migration imaging method based on a Muller model and a VTI medium, adopts a thin plate division idea to divide a model into thin plates, obtains a positive transmission down-going migration wave field and a positive transmission up-going migration wave field of each thin plate based on a seismic source calculation, then obtains a reverse transmission down-going migration wave field and a reverse transmission up-going migration wave field based on a seismic record, finally obtains all imaging results by using a cross-correlation imaging condition, and in addition, the Muller model is used as an absorption attenuation model, which is different from the existing Kolsky-Futterman model, can further avoid the problem that the existing model cannot accurately satisfy the minimum time delay criterion, can utilize down-going migration wave field information and up-going migration wave field information through the above process, and can correct the influence of anisotropic parameters on migration imaging results and compensate the amplitude attenuation problem of seismic waves through two one-way migration operators, so that the migration imaging precision of steep dip structures is effectively improved.

Description

Technical Field

[0001] This invention relates to a geophysical exploration method, specifically a multi-forward and single-backward two-way wave pre-stack migration imaging method based on the Müller model and VTI medium. Background Technology

[0002] As the objects of reflection seismic exploration become increasingly complex, the requirements for medium models also increase. To better characterize the absorption attenuation and anisotropic properties of real rocks, among numerous medium models, viscoacoustic anisotropic medium models (Müller model and VTI medium) offer the following advantages: 1) They can better approximate the main characteristics of viscous porous elastic media; 2) They can provide relatively good migration velocity models for reflection seismic exploration of deep-water sedimentary strata; 3) They effectively reduce the problems of numerous parameters, complex algorithms, and high computational costs associated with viscoelastic anisotropic media. Compared to differential equation-based micro-classification methods, integral classification methods have advantages such as the ability to filter operators and establish step-by-step algorithms, thus enabling their application in seismic migration imaging studies of viscoacoustic anisotropic media.

[0003] To improve the accuracy of migration imaging in viscoacoustic anisotropic media, some scholars have proposed a one-way wave migration method for viscoacoustic VTI media based on the DeWolf approximation. This method utilizes the one-way wave DeWolf approximation to achieve pre-stack migration imaging of VTI media. Although the one-way wave operator can compensate for the amplitude loss of seismic waves and eliminate the influence of anisotropic parameters on the migration imaging results to some extent, the operator has dip angle limitations, resulting in relatively low imaging accuracy for steep dip structures. At the same time, the computational efficiency of existing migration imaging methods is low.

[0004] To overcome this problem, the research direction of this invention is to provide a new offset imaging method that can effectively improve the imaging accuracy of steeply tilted structures while ensuring computational efficiency. Summary of the Invention

[0005] To address the problems existing in the prior art, this invention provides a two-way wave migration imaging method based on the Müller model and VTI medium. It combines down-travel migration wavefield information and one-way up-travel migration wavefield information, and constructs a two-way wave migration operator by using two unidirectional migration operators and cross-correlating them, thereby effectively improving the imaging accuracy of steep-tilt structures.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a two-way wave migration imaging method based on the Müller model and VTI medium, comprising the following steps:

[0007] Step 1: Arrange multiple geophones in a row at equal intervals in the area to be detected to form an observation system, and arrange multiple seismic sources in a row at equal intervals and number them.

[0008] Step 2: Excite each seismic source in sequence according to its number, and the observation system acquires the seismic record at the time of excitation of each seismic source.

[0009] Step 3: Set up the VIT medium model and Müller absorption attenuation model for the propagation of seismic scattered waves in the detection area.

[0010] Step 4: Using the idea of ​​thin plate division, the VIT medium model and Müller absorption attenuation model from Step 3 are divided into multiple thin plates along the depth direction, and the velocity parameters, quality factor parameters, and anisotropy parameters within each thin plate are determined.

[0011] Step 5: Load the Rack wavelet at the first source location and convert it to the frequency-wavenumber domain using Fourier transform. Use this as the initial value to calculate and obtain the forward down-propagation offset wavefield and forward up-propagation offset wavefield for all thin plates.

[0012] Step 6: Perform Fourier transform on the seismic record corresponding to the excitation of the first source to convert it to the frequency-wavenumber domain, and use it as the initial value to perform calculations using the calculation process in Step 5 to obtain the reverse propagation downlink migration wavefield and reverse propagation uplink migration wavefield of all thin plates.

[0013] Step 7: Using cross-correlation imaging conditions, perform migration imaging on the forward and reverse propagation downlink and uplink migration wavefields obtained in Steps 5 and 6 to obtain the imaging results of the first source.

[0014] Step 8: Repeat steps 5 to 7 for each seismic source in sequence according to the source number to obtain the imaging results of each source.

[0015] Step 9: Superimpose the imaging results of all seismic sources to obtain the final offset imaging result of the detection area.

[0016] Furthermore, in step one, a coordinate system is established with the X-axis along the deployment direction of each detector and seismic source and the Z-axis perpendicular to the deployment direction (i.e., the depth direction of the detection area), with the origin 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 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 dominant frequency of the wavelet is f Hz; the sampling interval is set to 1ms; and the sampling time length is nt.

[0017] Furthermore, step three is as follows:

[0018] I. Determine the speed parameters as follows ,in The coordinate parameters are nx for the X-axis and nz for the Z-axis, with grid spacings of dx and dz, respectively.

[0019] II. The expression for the Müller absorption attenuation model is as follows:

[0020] (1)

[0021] in , , It is the angular frequency. , and These are the reference frequencies.

[0022] III. Anisotropic parameters of the VTI medium model and The expression is:

[0023] (2)

[0024] in These are the maximum speed and the minimum speed, respectively.

[0025] Furthermore, step five specifically includes:

[0026] A. The velocity parameters within the j-th thin plate Quality factor parameters and anisotropy parameters and Decomposed into background parameters and disturbance parameters, where the background parameters include background velocity. Background quality factors and background anisotropy parameters and The disturbance parameters include velocity disturbances. Quality factor perturbation and anisotropic parameter perturbation and .

[0027] B. Loading a Ricker wavelet at the location of the first hypocenter. And it is transformed to the frequency wavenumber domain using Fourier transform. .

[0028] C. Perform phase-shift calculations in the frequency wavenumber domain to obtain... item.

[0029] (3)

[0030] in , , , .

[0031] D. Perform phase-shifting calculations in the frequency-wavenumber domain, including both frequency-wavenumber domain calculations and phase-shifting calculations. , , , .

[0032] (4)

[0033] (5)

[0034] (6)

[0035] (7)

[0036] in:

[0037] (8)

[0038] in for The corresponding transverse wavenumber, , , .

[0039] E. , , and Transformed into the frequency space domain as follows , , and After correction, the forward downshifting offset field of all thin plates is obtained. (down indicates a downward wave) (circular frequency) and forward propagation first upshift offset wave field The specific formula is as follows:

[0040] (9)

[0041] in:

[0042] (10)

[0043] in Represents Fourier transform, Represents the inverse Fourier transform. For position Vertical wavenumber at that location Represent , , Thin plate spacers, For velocity perturbation, and For anisotropic parameter perturbation, and These are the background anisotropy parameters; , For background speed, This is a virtual part unit.

[0044] First upshifted wave field The expression for (up represents an upward wave) is:

[0045] (11)

[0046] in .

[0047] F. After calculating to the bottom of the model, the first upward migration wave field of the second thin plate, sorted from deepest to shallowest. As its upward offset wave field The initial value is calculated according to formula (9) (when calculating the upward wave, the imaginary part i needs to be multiplied by a sign), and the previous upward offset wave field of the thin plate is added. The calculation continues using the initial values ​​of the adjacent thin plates above them; until the top of the model, the forward propagation upward migration wavefield of all thin plates is obtained. .

[0048] Furthermore, the calculation formula for cross-correlation imaging in step seven is as follows:

[0049] (8)

[0050] Where * represents the conjugate, and W represents the cutoff frequency.

[0051] Compared to the traditional viscoacoustic anisotropy De Wolf approximation one-way wave method, this invention first employs the idea of ​​thin-plate partitioning, dividing the model parameters into individual thin plates, which effectively improves computational efficiency. Then, based on the seismic source, the forward down-propagation migration wavefield and the first forward up-propagation migration wavefield are calculated sequentially for each thin plate from the shallowest to the deepest point. Next, the forward up-propagation migration wavefield is calculated for each thin plate at different locations from the deepest to the shallowest point. The first forward up-propagation migration wavefield of each thin plate is added to the previous one as the initial value for the next thin plate, and this process is repeated to obtain the forward up-propagation migration wavefield for all thin plates. Then, based on the seismic record, the reverse down-propagation migration wavefield and the reverse up-propagation migration wavefield are obtained according to the above calculation process. Finally, using… The cross-correlation imaging conditions are used to obtain all imaging results. In addition, the present invention uses the Müller model of VIT medium as the absorption attenuation model. Unlike the existing Kolsky-Futterman model, the Müller model of the present invention can further avoid the problem that the Kolsky-Futterman model cannot accurately satisfy the minimum time delay criterion. Through the above process, not only down-going migration wavefield information but also up-going migration wavefield information can be used. Two one-way migration operators are used to form a two-way wave migration operator, which can correct the influence of anisotropic parameters on migration imaging results while compensating for the amplitude attenuation problem of seismic waves, and finally effectively improve the migration imaging accuracy of steep-dipping structures. Attached Figure Description

[0052] Figure 1 These are the velocity model parameters, quality factor, and anisotropy parameters in embodiments of the present invention.

[0053] Where (a) are the velocity model parameters; (b) are the quality factor Q model; and (c) are the anisotropy parameters. Parameters; (d) Anisotropy parameter.

[0054] Figure 2 This is the offset imaging result of an embodiment of the present invention. Detailed Implementation

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

[0056] This invention includes the following steps:

[0057] Step 1: Arrange multiple geophones in a row at equal intervals to form an observation system in the required detection area, and arrange multiple seismic sources in a row at equal intervals and number them. Establish a coordinate system with the X-axis along the direction of each geophone and seismic source and the Z-axis perpendicular to the direction of arrangement (i.e., the depth direction of the detection area), with the origin at (0.0, 0.0). 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 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 dominant frequency of the wavelet is fHz, the sampling interval is 1ms, and the sampling time length is nt.

[0058] Step 2: Each seismic source is excited sequentially according to its number. The observation system acquires the seismic records at the time of excitation of each source. .

[0059] Step 3: Establish the VIT medium model and Müller absorption attenuation model for the propagation of seismic scattered waves in the detection area, specifically as follows:

[0060] I. Determine the speed parameters as follows ,in The coordinate parameters are nx for the X-axis and nz for the Z-axis, with grid spacings of dx and dz, respectively.

[0061] II. The expression for the Müller absorption attenuation model is as follows:

[0062] (1)

[0063] in , , It is the angular frequency. , and These are the reference frequencies.

[0064] III. Anisotropic parameters of the VTI medium model and The expression is:

[0065] (2)

[0066] in These are the maximum speed and the minimum speed, respectively.

[0067] Step 4: Using the idea of ​​thin plate partitioning, both the VIT medium model and the Müller absorption attenuation model from Step 3 are divided into multiple thin plates along the depth direction, totaling j plates. The thickness of each thin plate is dp, which is usually the thickness of one mesh, dz. and These represent the velocity parameter and quality factor parameter within a thin plate, respectively. and Anisotropic parameters within a thin plate, such as Figure 1 As shown, the anisotropy parameters and The calculation formula is as follows:

[0068] (2)

[0069] in and These are the maximum speed and the minimum speed, respectively.

[0070] The formula for calculating the quality factor parameter is as follows:

[0071] (3)

[0072] in The unit is km / s.

[0073] Step 5: Load the Ricker wavelet at the location of the first seismic source. (t is a time parameter), and then transformed to the frequency and wavenumber domain using Fourier transform. ( Using the angular frequency as an initial value, the forward downshifting offset wavefield and forward upshifting offset wavefield of all thin plates are obtained through calculation, specifically:

[0074] A. The velocity parameters within the j-th thin plate Quality factor parameters and anisotropy parameters and Decomposed into background parameters and disturbance parameters, where the background parameters include background velocity. Background quality factors and background anisotropy parameters and The disturbance parameters include velocity disturbances. Quality factor perturbation and anisotropic parameter perturbation and .

[0075] B. Load the Ricker wavelet at the location of the first hypocenter. And it is transformed to the frequency wavenumber domain using Fourier transform. .

[0076] C. Perform phase-shift calculations in the frequency wavenumber domain to obtain... item.

[0077] (3)

[0078] in , , , .

[0079] D. Perform phase-shifting calculations in the frequency-wavenumber domain, including both frequency-wavenumber domain calculations and phase-shifting calculations. , , , .

[0080] (4)

[0081] (5)

[0082] (6)

[0083] (7)

[0084] in:

[0085] (8)

[0086] in for The corresponding transverse wavenumber, , , .

[0087] E. , , and Transformed into the frequency space domain as follows , , and After correction, the forward downshifting offset field of all thin plates is obtained. (down indicates a downward wave) (circular frequency) and forward propagation first upshift offset wave field The specific formula is as follows:

[0088] (9)

[0089] in:

[0090] (10)

[0091] in Represents Fourier transform, Represents the inverse Fourier transform. For position Vertical wavenumber at that location Represent , , Thin plate spacers, For velocity perturbation, and For anisotropic parameter perturbation, and These are the background anisotropy parameters; , For background speed, This is a virtual part unit.

[0092] First upshifted wave field The expression for (up represents an upward wave) is:

[0093] (11)

[0094] in .

[0095] F. After calculating to the bottom of the model, the first upward migration wave field of the second thin plate, sorted from deepest to shallowest. As its upward offset wave field The initial value is calculated according to formula (9) (when calculating the upward wave, the imaginary part i needs to be multiplied by a sign), and the previous upward offset wave field of the thin plate is added. The calculation continues using the initial values ​​of the adjacent thin plates above them; until the top of the model, the forward propagation upward migration wavefield of all thin plates is obtained. .

[0096] Step 6: Perform a Fourier transform on the seismic record corresponding to the excitation of the first hypocenter to convert it to the frequency-wavenumber domain. Using this as the initial value, the calculation process in step five is performed to obtain the corresponding reverse propagation downlink offset wavefield. and rising wave field .

[0097] Step 7: Using cross-correlation imaging conditions, perform migration imaging on the forward and reverse propagation downlink and uplink migration wavefields obtained in Steps 5 and 6 to obtain the imaging results of the first source. The specific calculation formula is as follows:

[0098] (12)

[0099] Where * represents the complex conjugate, and W represents the cutoff frequency.

[0100] Step 8: Repeat steps 5 to 7 for each seismic source, arranged sequentially by source number, to obtain the imaging results for each source, as follows: , … .

[0101] Step 9: Superimpose the imaging results of all seismic sources, as shown below. Figure 2 As shown, the final offset imaging result of the detection area is obtained. The specific formula is as follows:

[0102] (13).

[0103] 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 two-way wave migration imaging method based on the Müller model and VTI medium, characterized in that, Includes the following steps: Step 1: Arrange multiple geophones in a row at equal intervals in the area to be detected to form an observation system, and arrange multiple seismic sources in a row at equal intervals and number them. Step 2: Excite each seismic source in sequence according to its number, and the observation system acquires the seismic record at the time of excitation of each seismic source; Step 3: Set up the VTI medium model and Müller absorption attenuation model for the propagation of seismic scattered waves in the detection area; Step 4: Using the idea of ​​thin plate division, the VTI medium model and Müller absorption attenuation model from Step 3 are divided into multiple thin plates along the depth direction, and the velocity parameters, quality factor parameters and anisotropy parameters within each thin plate are determined. Step 5: Load the Rack wavelet at the first source location and convert it to the frequency wavenumber domain using Fourier transform. Use this as the initial value to calculate and obtain the forward down-propagation offset wavefield and forward up-propagation offset wavefield of all thin plates. Step 6: Perform Fourier transform on the seismic record corresponding to the excitation of the first source to convert it to the frequency-wavenumber domain, and use it as the initial value to perform calculations using the calculation process in Step 5 to obtain the reverse propagation downlink migration wavefield and reverse propagation uplink migration wavefield of all thin plates. Step 7: Using cross-correlation imaging conditions, perform migration imaging on the forward and reverse propagation downlink and uplink migration wavefields obtained in Steps 5 and 6 to obtain the imaging results of the first source. Step 8: Repeat steps 5 to 7 for each seismic source in sequence according to the source number to obtain the imaging results of each source. Step 9: Superimpose the imaging results of all seismic sources to obtain the final offset imaging result of the detection area.

2. The two-way wave migration imaging method based on the Müller model and VTI medium 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, with the origin at (0.0, 0.0). There are nx detectors, the distance between adjacent detectors is dr, and the position of the first detector is (0.0, 0.0). There are ns seismic sources, the distance between adjacent seismic sources is ds, and the position of the first seismic source is (0.0, 0.0). The dominant frequency of the wavelet is f Hz, the sampling interval is set to l ms, and the sampling time length is nt.

3. The two-way wave migration imaging method based on the Müller model and VTI medium according to claim 1, characterized in that, Step three is as follows: I. Determine the speed parameters as follows ,in The coordinate parameters are nx for the X-axis and nz for the Z-axis, with grid spacings of dx and dz, respectively. II. The expression for the Müller absorption attenuation model is as follows: （1） in , , It is the angular frequency. , and These are the reference frequencies; III. Anisotropic parameters of the VTI medium model and The expression is: （2） in These are the maximum speed and the minimum speed, respectively.

4. The two-way wave migration imaging method based on the Müller model and VTI medium according to claim 1, characterized in that, Step five specifically involves: A. The velocity parameters within the j-th thin plate Quality factor parameters and anisotropy parameters and Decomposed into background parameters and disturbance parameters, where the background parameters include background velocity. Background quality factors and background anisotropy parameters and The disturbance parameters include velocity disturbances. Quality factor perturbation and anisotropic parameter perturbation and ; B. Load the Ricker wavelet at the location of the first hypocenter. And it is transformed to the frequency wavenumber domain using Fourier transform. ; C. Perform phase-shift calculations in the frequency wavenumber domain to obtain... item: （3） in , , , ; D. Perform phase-shifting calculations in the frequency-wavenumber domain, including both frequency-wavenumber domain calculations and phase-shifting calculations. , , , ; （4） （5） （6） （7） in: （8） in for The corresponding transverse wavenumber, , , E. , , and Transformed into the frequency space domain as follows , , and After correction, the forward downshifting offset field of all thin plates is obtained. And the forward transmission of the first upward offset wave field The specific formula is as follows: （9） in: （10） in Represents Fourier transform, Represents the inverse Fourier transform. For position Vertical wavenumber at that location Represent , , Thin plate spacers, For velocity perturbation, and For anisotropic parameter perturbation, and These are the background anisotropy parameters; , For background speed, For imaginary units; First upshifted wave field The expression is: （11） in ; F. After calculating to the bottom of the model, the first upward migration wave field of the second thin plate, sorted from deepest to shallowest. As its upward offset wave field The initial value is calculated according to formula (9), and is added to the previous upward offset wave field of the thin plate. The calculation continues using the initial values ​​of the adjacent thin plates above them; until the top of the model, the forward propagation upward migration wavefield of all thin plates is obtained. .

5. The two-way wave migration imaging method based on the Müller model and VTI medium according to claim 1, characterized in that, The calculation formula for cross-correlation imaging in step seven is as follows: （12） Where * represents the conjugate, and W represents the cutoff frequency.

Images