Gaussian beam form optimization method based on minimum mean square width constraint

Abstract

The invention belongs to the technical field of seismic migration imaging, and particularly relates to a Gaussian beam form optimization method based on minimum mean square width constraint. According to the method, symmetric matrix parameters representing the optimal initial form of a Gaussian beam in a plane perpendicular to a central ray are calculated, a smooth objective function represented by the initial parameters is minimized, and a novel Gaussian beam is constructed; a green function is analyzed through integral superposition fitting of a novel Gaussian beam to solve an integral weight coefficient, then the green function and an uplink and downlink continuation wave field are constructed, and an imaging value at an underground imaging point is extracted by applying a deconvolution imaging condition. Compared with a traditional Gaussian beam migration method, the Gaussian beam form optimization method based on the minimum mean square width constraint solves the contradiction between the imaging precision of a near-surface layer and the imaging precision of a middle-deep layer, enables a migration result profile to be more accurate and reliable, and facilitates subsequent interpretation processing.

Description

Technical Field

[0001] The invention belongs to the technical field of seismic migration imaging, and in particular relates to a Gaussian beam morphology optimization method constrained by minimum mean square width. Background Art

[0002] As the difficulty of exploration and development in my country gradually increases, exploration targets have begun to shift to complex structural exploration areas such as fractured oil and gas reservoirs and lithologic oil and gas reservoirs. Complex undulating surfaces have gradually become normalized, and the lateral sudden changes in near-surface and mid-deep velocities have made the processing of seismic data more difficult. Research on high-precision and highly adaptable prestack depth migration algorithms is of great significance in subsequent seismic interpretation, attribute analysis and reservoir prediction.

[0003] Several commonly used depth domain imaging methods in the industry are reverse time migration, Fourier finite difference migration, Kirchhoff migration and Gaussian beam migration. The above four methods have their own advantages and disadvantages: reverse time migration based on the two-way exact solution of the wave equation is the migration method with the highest imaging accuracy so far, but this method has low computational efficiency and consumes a lot of memory; Fourier finite difference migration can act on both the time domain and the frequency domain, but its imaging accuracy cannot be compared with reverse time migration; Kirchhoff migration is widely used in commercial-scale production due to its unique computational efficiency advantage and matching ability for irregular observation data, but it cannot handle complex structural imaging problems in geophysical exploration and cannot meet the requirements of fine exploration of target blocks. The Gaussian beam migration method can not only solve the theoretical defects of Kirchhov migration (such as the singular amplitude in the caustic zone, the lack of wave field energy in the shadow zone, and the problem of multiple wave arrivals), but also inherit its flexibility and high efficiency, and has imaging adaptability comparable to that of wave equation migration for complex high-steep structures and block areas; in addition, its dependence on the initial depth domain velocity field is weaker than that of the wave equation migration method, and the velocity model obtained by the conventional travel time tomography method can meet the imaging requirements.

[0004] The Chinese invention patent application with publication number CN111210495A discloses a three-dimensional model driving method, device, terminal and computer-readable storage medium, and the method includes: determining the target skeleton of the three-dimensional model, the target skeleton is obtained from at least two candidate skeletons; binding the target skeleton to the three-dimensional model, wherein the target skeleton is used to drive the movement of the three-dimensional model. The Chinese invention patent application with publication number CN110138426A discloses a panoramic beam construction method, system, device and medium based on satellite communication, which is applied to a ground station, and the method includes: simulating historical downlink data modulated by satellite; demodulating the downlink data modulated by simulated satellite to obtain demodulation parameters corresponding to a set number of spacecraft; obtaining real-time downlink data modulated by satellite, demodulating the real-time modulated downlink data based on the demodulation parameters, and obtaining signal characteristic data of a set number of spacecraft carried in the real-time modulated downlink data. It obtains the demodulation parameters required by the ground station equipment by simulating actual scenarios, and improves the corresponding ground station equipment based on the obtained demodulation parameters. Based on the improved ground station equipment, it can realize the deployment of static panoramic beams, so that the communication satellite communicating with the ground station equipment can cover all spacecraft within its coverage area, and realize real-time communication between the user center communicating with the ground station and each spacecraft.

[0005] In the prior art, the traditional Gaussian beam migration method is usually implemented in the frequency domain, among which the most typical one is the shot domain Gaussian beam migration algorithm proposed by Hill and Gray, etc.; since then, the algorithm has been gradually extended to true amplitude, undulating surface, anisotropic media, and many new seismic beam migration methods, such as control beam, fast beam, focused beam, Fresnel beam, complex valued beam, adaptive focused beam, etc. Although Gaussian beam migration overcomes the shortcomings of Kirchhoff migration that cannot handle multiple arrivals and one-way wave wave equation migration that cannot accurately image steep structures, under the dual complex conditions of complex surface and complex structure in land exploration, its imaging accuracy depends on the selected initial beam width, that is, when the initial beam width is small, the near-surface imaging accuracy is high, but the imaging quality of the middle and deep layers is poor; on the contrary, when the initial beam width is large, the imaging quality of the middle and deep layers is improved, but the near-surface imaging accuracy is reduced. Therefore, in order to obtain a more ideal migration result profile, that is, to ensure the imaging quality of the mid-depth layer while taking into account the imaging accuracy of the near-surface, a new Gaussian beam morphology optimization method is needed. Summary of the invention

[0006] The purpose of the present invention is to provide a Gaussian beam morphology optimization method with minimum mean square width constraint to solve the contradiction in imaging accuracy between near-surface and mid-depth layers of the Gaussian beam migration method, so as to make the migration result profile more accurate and reliable, which is conducive to subsequent interpretation and processing.

[0007] In order to achieve the above purpose, the technical solution adopted by the present invention is:

[0008] A Gaussian beam morphology optimization method with minimum mean square width constraint comprises the following steps:

[0009] (1) Calculate and store the ray propagation matrix, slowness vector Cartesian components, and transformation matrix based on the ray path, departure angle, and initial point position;

[0010] (2) Calculating the symmetric matrix parameters that characterize the optimal initial shape of the Gaussian beam in the plane perpendicular to the central ray, and storing the initial parameters of the Gaussian beam on regular grid points;

[0011] (3) Selecting a smooth objective function that minimizes the initial parameter representation to obtain a smooth distribution of the initial parameters and then construct a new Gaussian beam;

[0012] (4) The integral weight coefficient is solved by analyzing the Green's function through the integral superposition fitting of the new Gaussian beam, and then the Green's function and the uplink and downlink extended wave fields are constructed, and the imaging value at the underground imaging point is extracted using the deconvolution imaging condition.

[0013] Compared with the traditional Gaussian beam migration method, the Gaussian beam morphology optimization method with minimum mean square width constraint provided by the present invention solves the contradiction between the imaging accuracy of the near-surface and the middle-depth layers, making the migration result profile more accurate and reliable, which is conducive to subsequent interpretation and processing.

[0014] Preferably, in step (2), the symmetric matrix parameters are obtained using the target functional under the minimum mean square width constraint, and the target functional under the minimum mean square width constraint is:

[0015]

[0016] Among them, s 1 ,s 2 Represent any two points of the fixed part of the central ray, tr represents the trace of the matrix, Re represents the real part, G is the weighting matrix, and Ψ(s) is the intermediate matrix variable, that is,

[0017]

[0018] Among them, the weight matrix G is in the form of R(s) and Y(s) represent the real and imaginary parts of the matrix Ψ(s), respectively.

[0019] Further preferably, in step (3), the minimization is solved iteratively by formula (3):

[0020]

[0021] Where O represents the minimum value of equation (3), T is the target functional in equation (1), and R 0 and Y 0are the initial matrices of R(s) and Y(s) at the point where the central ray emerges, respectively. The operator |||| 2 It means taking the square of the matrix's modulus.

[0022] Preferably, in step (3), the construction of the new Gaussian beam is to obtain the new Gaussian beam function U(s,n,ω) by calculation using the following formula:

[0023]

[0024] Among them, τ(s) is the travel time obtained by kinematic ray tracing, v(s) is the value of the input model velocity field; p(s) and q(s) are the dynamic parameters of the seismic beam obtained during ray tracing, (s,n) is the ray center coordinate system, n represents the distance of the paraxial ray in the direction perpendicular to the central ray; ω represents the beam frequency.

[0025] Further preferably, in step (4), the Green's function G(x, x′, ω) is calculated using the following formula:

[0026]

[0027]

[0028] Among them, Φ(φ) is its integral weight coefficient; U φ (x, x′, ω) is a beam emitted at an angle φ; x′, x represent the center ray emission position of the beam and the underground target imaging point respectively; A(x, x′) and T(x, x′) are the amplitude and travel time of the seismic beam respectively; π is pi; p 1 (s),p 2 (s) are the two components of the complex variable p(s), q 1 (s),q 2 (s) are the two components of the complex variable q(s), subscript 1 represents the imaginary part, and subscript 2 represents the real part; ε(x) represents the complex-valued beam parameter of the underground target imaging point x; v(x′), v(x) are the values ​​of the input model velocity field at the central ray exit position x′ and the underground target imaging point x, respectively; τ(x, x′) is the travel time of the central ray from the central ray exit position x′ to the underground target imaging point x.

[0029] More preferably, in step (4), the shot point x s , detection point x r Perform ray tracing downwards to construct the up-and-down extension wave field of the gun inspection point:

[0030]

[0031] Among them, v s 、v rThey are respectively the gun point x s , detection point x r Speed ​​at R(x r ,x s ,ω) is the shot point x s Excitation detection point x r Received earthquake records; and

[0032] φ s =α s -β s ,φ r =α r -β r (8)

[0033] Among them, α s and β s They are respectively the gun point x s The ray departure angle and the surface inclination at α r and β r They are detection points x r The ray departure angle and the surface inclination at P D represents the downlink forward propagation wave field of the shot point, P U It indicates the upward back propagation extension wave field of the detection point.

[0034] Preferably, in step (4), the imaging value at the underground imaging point is obtained using the following formula:

[0035]

[0036] Where I(x) represents the imaging value at the underground imaging point x, and the superscript * represents the operator conjugate. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a flowchart of the Gaussian beam morphology optimization method with minimum mean square width constraint of the present invention;

[0038] Figure 2 The beam in the present invention moves from the initial position P(s) along the central ray Ω 0 ,0) to point Q(s,0);

[0039] Figure 3 The figure is a beam half-width distribution diagram of the smoothed Marmousi velocity model in Example 1 of the present invention; the upper left is a smoothed Marmousi velocity model, the upper right is a conventional Gaussian beam half-width distribution diagram, and the lower row is a new Gaussian beam half-width distribution diagram, the left is a first iteration, and the right is a second iteration;

[0040] Figure 4It is the Green's function characterized by the new Gaussian beam at different angular intervals in the homogeneous medium in Example 1 of the present invention; 2° on the upper left, 5° on the upper right, 10° on the lower left, and 20° on the lower right;

[0041] Figure 5 The depth domain migration velocity model of the work area 1 in Example 1 of the present invention;

[0042] Figure 6 The depth domain migration result of the work area 1 in Example 1 of the present invention; the upper one is Gaussian beam migration, and the lower one is new Gaussian beam migration;

[0043] Figure 7 The depth domain migration velocity model of the work area 2 in Example 1 of the present invention;

[0044] Figure 8 1 and 1 are the depth domain migration results of work area 2 in Example 1 of the present invention; the upper one is Gaussian beam migration, and the lower one is new Gaussian beam migration. DETAILED DESCRIPTION

[0045] The implementation process of the present invention is described in detail below in conjunction with specific embodiments.

[0046] Example 1

[0047] The Gaussian beam morphology optimization method with minimum mean square width constraint in this embodiment has a flow chart as shown in FIG. Figure 1 As shown, the following steps are included:

[0048] (1) Based on the ray path, the departure angle, and the initial point position, the ray propagation matrix, the slowness vector Cartesian components, and the transformation matrix are calculated and stored.

[0049] This step is the data preparation process, which prepares for step (2) to determine the optimal shape of the Gaussian beam.

[0050] Specifically, a smooth model velocity field is input, and ray tracing is performed on the velocity field to obtain the travel time, amplitude, and dynamic parameters p(s) and q(s) of the seismic beam. As the ray arc length s changes continuously, the various parameters such as the ray propagation matrix in step (1) are continuously updated during the ray tracing process. Different from the conventional Gaussian beam method, two additional parameters, the Cartesian component of the slowness vector and the transformation matrix, need to be stored here.

[0051] Figure 2 The beam is shown along the central ray from the initial position (s 0 ,0) to point (s,0).

[0052] (2) Calculate the symmetric matrix parameters that characterize the optimal initial shape of the Gaussian beam in the plane perpendicular to the central ray, and store the initial parameters of the Gaussian beam on regular grid points.

[0053] In this step, when the beam width takes the minimum mean square width, it can be considered as the best initial form. The "symmetric matrix parameters" are included in the initial parameters of the Gaussian beam, and the Gaussian beam is subsequently constructed according to formula (4); regular grid points refer to the equally spaced grids that divide the underground imaging area.

[0054] This step is to obtain the symmetric matrix parameters using the target functional under the minimum mean square width constraint. The target functional under the minimum mean square width constraint is:

[0055]

[0056] Among them, s 1 ,s 2 Represent any two points of the fixed part of the central ray, tr represents the trace of the matrix, Re represents the real part, G is the weighting matrix, and Ψ(s) is the intermediate matrix variable, that is,

[0057]

[0058] Among them, the weight matrix G is in the form of

[0059] R(s) and Y(s) represent the real and imaginary parts of the matrix Ψ(s), respectively. The superscript [] -1 Represents the inverse operation.

[0060] (3) Select a smooth objective function that minimizes the initial parameters to obtain a smooth distribution of the initial parameters, and then construct a new Gaussian beam.

[0061] In this step, the smoothing objective function represented by the initial parameters is minimized by iteratively solving formula (3):

[0062]

[0063] Where O represents the minimum value of equation (3), T is the target functional in equation (1), and R 0 and Y 0 are the initial matrices of R(s) and Y(s) at the point where the central ray emerges, respectively. The operator || || 2 It means taking the square of the matrix's modulus.

[0064] The Gaussian beam function U(s,n,ω) is constructed using the following formula:

[0065]

[0066] Among them, τ(s) is the travel time obtained by kinematic ray tracing, v(s) is the value of the input model velocity field, p(s) and q(s) are the dynamic parameters of the seismic beam obtained during ray tracing, (s,n) is the ray center coordinate system, and n represents the distance of the paraxial ray in the direction perpendicular to the central ray ( Figure 2 ), exp represents the exponential function of the natural constant e, ω represents the beam frequency,

[0067] In order to verify the effectiveness of the new beam building strategy, numerical simulations from simple models to complex models were carried out.

[0068] First consider the smoothed Marmousi velocity model, such as Figure 3 As shown in the upper left figure, the color RGB values ​​in the other three figures describe the half-width of the Gaussian beam, not the mean square half-width. Yellow corresponds to the half-width of the Gaussian beam at all frequencies of 0 meters. The half-width increases from yellow to green and from blue to red, corresponding to the standard half-width of 1010 meters and above (frequency is 35 Hz) and 1890 meters and above (frequency is 10 Hz). By comparison, it is found that the beam half-width of the conventional Gaussian beam increases rapidly ( Figure 3 This will result in low calculation accuracy when constructing the Green function, which will in turn affect the construction accuracy and imaging quality of the seismic wave field. The new seismic beam is generated by one iteration ( Figure 3 Bottom left) and two iterations ( Figure 3 The optimal initial parameters characterizing the Gaussian beam morphology are smoothed to a constant value, and the half-width distribution results are in line with expectations. The smoothing constraint can not only minimize the width of the Gaussian beam, but also minimize the secondary changes along any surface, structural interface, or wavefront section. The iteration can also converge quickly. Smoothing the initial beam parameters three or more times will no longer change the half-width distribution, optimize the Gaussian beam morphology, and effectively reduce the calculation error of constructing the Green's function. Therefore, it is beneficial to improve the accuracy of seismic wave field construction and the final migration quality, especially in complex surface exploration areas with drastic changes in elevation and near-surface velocity.

[0069] (4) The integral weight coefficient is solved by analyzing the Green's function through the integral superposition fitting of the new Gaussian beam, and then the Green's function and the uplink and downlink extended wave fields are constructed, and the imaging value at the underground imaging point is extracted using the deconvolution imaging condition.

[0070] In this step, the Green's function is characterized by the superposition form of the new Gaussian beam in step (3), and the pre-stack depth domain migration imaging formula of the double complex condition is derived by applying the deconvolution imaging condition.

[0071] Specifically, the Green's function G(x, x′, ω) is calculated using the following formula:

[0072]

[0073]

[0074] Among them, Φ(φ) is its integral weight coefficient; U φ (x, x′, ω) is a beam emitted at an angle φ; x′, x represent the center ray emission position of the beam and the underground target imaging point respectively; A(x, x′) and T(x, x′) are the amplitude and travel time of the seismic beam respectively; π is pi; p 1 (s),p 2 (s) are the two components of the complex variable p(s), q 1 (s),q 2 (s) are the two components of the complex variable q(s), subscript 1 represents the imaginary part, and subscript 2 represents the real part; ε(x) represents the complex-valued beam parameter of the underground target imaging point x; v(x′), v(x) are the values ​​of the input model velocity field at the central ray exit position x′ and the underground target imaging point x, respectively; τ(x, x′) is the travel time of the central ray from the central ray exit position x′ to the underground target imaging point x.

[0075] For ease of explanation, the system of dynamic ray tracing equations is given here:

[0076] Where v(s) is the propagation speed mentioned in formula (4). 1 (s),p 1 (s)], [q 2 (s),p 2 (s)] is the complex linear combination solution of the dynamic ray tracing equations. If the dynamic ray tracing equations choose [q 1 (s),p 1 (s)], we can get the plane wave solution of the wave equation; if the dynamic ray tracing equation is chosen [q 2 (s),p 2 (s)], the point source solution of the wave equation can be obtained; in order to make the seismic wave field regular everywhere, the complex-valued initial beam parameter ε = -iε is introduced i (ε r =0,ε i > 0) The point source solution and the plane wave solution are combined into a complex linear combination, that is, p(s) = ε(s)p 1 (s)+p 2 (s),q(s)=ε(s)q 1 (s)+q 2 (s), further solving the dynamic ray tracing equation can cleverly avoid the amplitude singularity problem.

[0077] By gun point x s , detection point x rPerform ray tracing downwards to construct the up-and-down extension wave field of the gun inspection point:

[0078]

[0079] Among them, v s 、v r They are respectively the gun point x s , detection point x r Speed ​​at R(x r ,x s ,ω) is the shot point x s Excitation detection point x r Received earthquake records; and

[0080] φ s =α s -β s ,φ r =α r -β r (8)

[0081] Among them, α s and β s They are respectively the gun point x s The ray departure angle and the surface inclination at α r and β r They are detection points x r The ray departure angle and the surface inclination at P D represents the downlink forward propagation wave field of the shot point, P U It indicates the upward back propagation extension wave field of the detection point.

[0082] The imaging value at the underground imaging point is obtained using the following formula:

[0083]

[0084] In formula (9), I(x) represents the imaging value at the underground imaging point x, the superscript * represents the operator conjugate, and cos is the cosine function.

[0085] Figure 4 The Green's function constructed by the new Gaussian beam under different angle intervals in a uniform medium at 2000m / s. As the ray angle interval increases, the approximate accuracy of the Green's function represented by the new Gaussian beam gradually decreases, especially at far sources. This is because the new Gaussian beam has the same beam width on the entire ray in a uniform medium. Therefore, when the ray angle interval is large, the beam cannot achieve sufficient illumination at far sources. Therefore, in practical applications, the new Gaussian beam should pay attention to the selection of ray angle intervals. Too large angle intervals should not be used. Generally, 1 to 5 degrees can achieve sufficient illumination of the underground medium.

[0086] In order to test the robustness of the new Gaussian beam migration method, two actual work areas are tested.

[0087] The conventional Gaussian beam method does not optimize the initial morphological parameters. During data preparation, the conventional Gaussian beam method does not need to store the two additional parameters of the slowness vector Cartesian component and the transformation matrix, and there are no steps (2) and (3). It is directly constructed by formula (4).

[0088] Figure 5 is the migration velocity model for work area 1, Figure 6 The upper sub-figure shows the migration result of the conventional Gaussian beam method in work area 1, and the lower sub-figure shows the migration result of the new Gaussian beam method in work area 1.

[0089] Figure 7 is the migration velocity model for work area 2, Figure 8 The upper sub-figure shows the migration result of the conventional Gaussian beam method in work area 2, and the lower sub-figure shows the migration result of the new Gaussian beam method in work area 2.

[0090] Since the Gaussian beam uses constant initial beam parameters and has a larger beam width in the middle and deep layers, its imaging results have larger offset noise. However, the new Gaussian beam method of the present invention constrains the main energy of the seismic beam to the neighborhood of the central ray, so it has a smaller wave field energy mapping range in the imaging space, which helps to effectively focus the imaging energy, reduce offset noise, and improve the imaging quality of near-surface and middle and deep layers.

[0091] The embodiments of the present invention have been described above, and the above description is exemplary, not exhaustive, and is not limited to the disclosed embodiments. Many modifications and changes will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments.

Claims

1. A Gaussian beam morphology optimization method with minimum mean square width constraint, characterized in that: The following steps are involved: (1) Calculate and store the ray propagation matrix, slowness vector Cartesian components, and transformation matrix based on the ray path, departure angle, and initial point position; (2) Calculating the symmetric matrix parameters that characterize the optimal initial shape of the Gaussian beam in the plane perpendicular to the central ray, and storing the initial parameters of the Gaussian beam on regular grid points; (3) Selecting a smooth objective function that minimizes the initial parameter representation to obtain a smooth distribution of the initial parameters and then construct a new Gaussian beam; (4) The integral weight coefficient is solved by analyzing the Green's function through the integral superposition fitting of the new Gaussian beam, and then the Green's function and the uplink and downlink extended wave fields are constructed, and the imaging value at the underground imaging point is extracted using the deconvolution imaging condition.

2. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 1, characterized in that: In step (2), the symmetric matrix parameters are obtained using the target functional under the minimum mean square width constraint, and the target functional under the minimum mean square width constraint is: Among them, s1 and s2 represent any two points of the fixed part of the central ray, tr represents the trace of the matrix, Re represents the real part, G is the weighting matrix, and Ψ(s) is the intermediate matrix variable, that is, Among them, the weight matrix G is in the form of R(s) and Y(s) represent the real and imaginary parts of the matrix Ψ(s), respectively.

3. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 2, characterized in that: In step (3), the minimization is performed by iteratively solving formula (3): Where O represents the minimum value of equation (3), T is the target functional in equation (1), R0 and Y0 are the initial matrices of R(s) and Y(s) at the exit point of the central ray, respectively, and the operator || || 2 It means taking the square of the matrix's modulus.

4. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 1, characterized in that: In step (3), the construction of the new Gaussian beam is to obtain the new Gaussian beam function U(s,n,ω) by calculation using the following formula: Among them, τ(s) is the travel time obtained by kinematic ray tracing, v(s) is the value of the input model velocity field; p(s) and q(s) are the dynamic parameters of the seismic beam obtained during ray tracing, (s,n) is the ray center coordinate system, n represents the distance of the paraxial ray in the direction perpendicular to the central ray; ω represents the beam frequency.

5. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 4, characterized in that: In step (4), the Green's function G(x, x′, ω) is calculated using the following formula: Among them, Φ(φ) is its integral weight coefficient; U φ (x, x′, ω) is a beam emitted at an angle φ; x′, x represent the center ray emission position of the beam and the underground target imaging point, respectively; A(x, x′) and T(x, x′) are the amplitude and travel time of the seismic beam, respectively; π is pi; p1(s) ​​and p2(s) are two components of the complex variable p(s), q1(s) and q2(s) are two components of the complex variable q(s), subscript 1 represents the imaginary part, and subscript 2 represents the real part; ε(x) represents the complex-valued beam parameter of the underground target imaging point x; v(x′) and v(x) are the values ​​of the input model velocity field at the center ray emission position x′ and the underground target imaging point x, respectively; τ(x, x′) is the travel time of the center ray from the center ray emission position x′ to the underground target imaging point x.

6. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 5, characterized in that: In step (4), the gun point x s , detection point x r Perform ray tracing downwards to construct the up-and-down extension wave field of the gun inspection point: Among them, v s 、v r They are respectively the gun point x s , detection point x r Speed ​​at R(x r ,x s ,ω) is the shot point x s Excitation detection point x r Received earthquake records; and f s =a s -b s ,f r =a r -b r (8) Among them, α s and β s They are respectively the gun point x s The ray departure angle and the surface inclination at α r and β r They are detection points x r The ray departure angle and the surface inclination at P D represents the downlink forward propagation wave field of the shot point, P U It indicates the upward back propagation extension wave field of the detection point.

7. The Gaussian beam morphology optimization method with minimum mean square width constraint as claimed in claim 6, characterized in that: In step (4), the imaging value at the underground imaging point is obtained using the following formula: Where I(x) represents the imaging value at the underground imaging point x, and the superscript * represents the operator conjugate.

Images