Least square reverse time migration method and computer readable storage medium

Abstract

The invention provides a least square reverse time migration method and a computer readable storage medium, and belongs to the technical field of seismic wave imaging. The method comprises the steps of firstly obtaining a background wave field and a back propagation wave field, performing vector decomposition to obtain upper, lower, left and right traveling waves of the background wave field and the back propagation wave field, then performing gradient solution and suppressing low-frequency noise according to a vector decomposition imaging condition, and obtaining a final imaging result by using a target functional of least square reverse time migration based on an L-2 norm. According to the invention, low-frequency noise can be suppressed, the convergence speed of iterative calculation is improved, and the imaging quality is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic wave imaging, and in particular relates to a least squares reverse time migration method and a computer-readable storage medium. Background Art

[0002] The least squares reverse time migration method is an inversion imaging method under the least squares framework. In essence, it is to deblur the conventional reverse time migration imaging results by solving the Hessian matrix, thereby achieving high-resolution imaging. However, due to the huge amount of calculation of the Hessian matrix, it cannot be calculated directly.

[0003] Therefore, most of the least squares reverse time migration is currently implemented in the data domain, that is, by constructing the target functional for iterative solution, gradually approximating the Hessian matrix without explicitly solving the Hessian matrix. However, in the actual imaging process, the conventional cross-correlation imaging conditions used will cause severe low-frequency noise in the least squares gradient solution process, which on the one hand reduces the convergence speed of the calculation, and on the other hand, it is likely to cause severe low-frequency noise in the final imaging result, thereby affecting the final imaging quality. Therefore, a least squares reverse time migration method that can suppress severe low-frequency noise is needed. Summary of the invention

[0004] The purpose of the present invention is to solve the above-mentioned difficulties in the prior art, and to provide a least squares reverse time migration method and a computer-readable storage medium. Under the Born approximation assumption, the shot point wave field and the detection point wave field are decomposed into four directions of up, down, left and right by using the traveling wave decomposition method, and then the imaging conditions are modified again to suppress low-frequency noise, increase the convergence speed of iterative calculation, and improve the imaging quality.

[0005] The present invention is achieved through the following technical solutions:

[0006] A first aspect of the present invention provides a least squares reverse time migration method, which first obtains a background wave field and a back propagation wave field, and performs vector decomposition to obtain the up, down, left, and right traveling waves of the background wave field and the back propagation wave field, then performs gradient solution according to the vector decomposition imaging condition, suppresses low-frequency noise, and obtains the final imaging result using the target functional of the least squares reverse time migration based on the L-2 norm.

[0007] A further improvement of the present invention is:

[0008] The method comprises:

[0009] (1) Input velocity field, original shot record, and seismic wavelet;

[0010] (2) Obtain the opposite number of the original shot record and use it as residual data;

[0011] (3) Calculate the background wave field and the back-propagation wave field, and separate the background wave field and the back-propagation wave field into upward, downward, left, and right traveling waves;

[0012] (4) Using the four directional traveling waves separated from the background wave field and the back-propagation wave field, a new imaging condition of vector decomposition is constructed;

[0013] (5) Calculate the iteration step size;

[0014] (6) Update imaging results;

[0015] (7) Obtain residual data and obtain the value of the target functional;

[0016] (8) Determine whether the termination condition is met. If yes, proceed to step (9); if no, return to step (3);

[0017] (9) Output the final imaging result.

[0018] A further improvement of the present invention is:

[0019] The operations of calculating the background wave field and the back propagation wave field in step (3) include:

[0020] The background wave field is obtained by calculation using the following formula:

[0021]

[0022] Among them, u 0 (x,z;t) represents the background wave field, v 0 (x,z) represents the velocity field, f represents the seismic wavelet, t represents time, and x and z represent the directions of the X and Z axes of space, respectively;

[0023] The back propagation wave field is calculated using the following formula:

[0024]

[0025] where u s (x,z;t) represents the back propagation wave field, v 0 (x,z) represents the velocity field, d res Represents residual data, t represents time, and x and z represent the X and Z axis directions of space respectively.

[0026] A further improvement of the present invention is:

[0027] In step (3), the Poynting vector is used to separate the upward, downward, left and right traveling waves of the background wave field and the back propagation wave field. The specific operations include:

[0028] The following formula is used to perform vector decomposition of the background wave field:

[0029]

[0030]

[0031]

[0032]

[0033]

[0034]

[0035] in, Indicates the up, down, left, and right traveling waves of the background wave field, u 0 (x,z;t) is the background wave field, Respectively represent the horizontal and vertical components of the Poynting vector;

[0036] The vector decomposition of the back propagation wave field is performed using the following formula:

[0037]

[0038]

[0039]

[0040]

[0041]

[0042]

[0043] in, Indicates the up, down, left, and right traveling waves of the background wave field, u s (x,z;t) is the back propagation wave field, Represent the horizontal and vertical components of the Poynting vector respectively.

[0044] A further improvement of the present invention is:

[0045] The imaging conditions of the new vector decomposition in step (4) are as follows:

[0046]

[0047] Where t represents time, i represents the shot number, N represents the total number of shots in the seismic channel, and g represents the total number of shots in the seismic channel. k represents the gradient of the kth iteration, v 0(x,z) represents the velocity field.

[0048] A further improvement of the present invention is:

[0049] The operation of step (5) includes:

[0050] Using the conjugate gradient method, calculate the iteration step size:

[0051]

[0052] z k+1 =g k+1 +βz k

[0053]

[0054] Among them, g k represents the gradient of the kth iteration, L is the forward operator, α is the iteration step size, β and z are intermediate variables.

[0055] A further improvement of the present invention is:

[0056] The operation of step (6) includes:

[0057] Update the imaging results using the following formula:

[0058] m k+1 =m k -αg k

[0059] Among them, m is the imaging result, α is the iteration step, g k represents the gradient of the kth iteration, where k is the number of iterations.

[0060] A further improvement of the present invention is:

[0061] In step (7), the following formula is used to obtain the residual data:

[0062] d res =d cal -d obs

[0063] Among them, d cal represents simulated data, d cal =Lm, L represents the forward operator, d obs represents the original gun record, d res represents residual data;

[0064] The operation of obtaining the value of the target functional in step (7) includes:

[0065] Construct the target functional of the least squares reverse time migration based on the L-2 norm:

[0066] g min =||d cal -d obs || 2 (11)

[0067] Calculate the value of the target functional g min .

[0068] A further improvement of the present invention is:

[0069] The termination condition in step (8) is as follows:

[0070] g min The value of is less than the set threshold, or k is equal to the set number of iterations.

[0071] According to a second aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps in the least squares reverse time migration method described above.

[0072] Compared with the prior art, the beneficial effects of the present invention are: utilizing the present invention can suppress low-frequency noise, increase the convergence speed of iterative calculation, and improve imaging quality. BRIEF DESCRIPTION OF THE DRAWINGS

[0073] Figure 1 The velocity model of the embodiment;

[0074] Figure 2 The original gun records entered;

[0075] Figure 3 Theoretical scattering intensity;

[0076] Figure 4 Reverse time migration results;

[0077] Figure 5 Conventional least squares reverse time migration results after 40 iterations;

[0078] Figure 6 Vector decomposition least squares reverse time migration results after 40 iterations;

[0079] Figure 7 Offset distance 1500m Figure 3 , Figure 5 and Figure 6 Amplitude curve comparison chart of ;

[0080] Figure 8 It is the normalized curve graph of conventional least squares reverse time migration and vector decomposition least squares reverse time migration;

[0081] Fig. 9 A flowchart of the steps of the method of the present invention. DETAILED DESCRIPTION

[0082] The present invention is further described in detail below in conjunction with the accompanying drawings:

[0083] The least squares reverse time migration technique directly solves the two-way wave equation without approximate processing, and is a high-precision, high-resolution migration imaging method. However, the conventional least squares migration method generally uses cross-correlation imaging conditions, which leads to a large amount of low-frequency noise in the gradient term, affecting the imaging quality and iterative convergence speed, and hindering the application of the least squares reverse time migration technique on actual data.

[0084] The present invention provides a least squares reverse time migration method for vector decomposition imaging conditions, which is still based on the least squares inversion framework. The method of the present invention first obtains the background wave field and the back propagation wave field, and performs vector decomposition to obtain the up, down, left and right traveling waves of the background wave field and the back propagation wave field, and then performs gradient solution according to the vector decomposition imaging conditions to suppress low-frequency noise, and obtains the final imaging result using the target functional of the least squares reverse time migration based on the L-2 norm.

[0085] The embodiments of the method of the present invention are as follows:

[0086] Embodiment 1:

[0087] like Fig. 9 As shown, the method includes:

[0088] (1) Input velocity field, original shot record (i.e., observation data), and seismic wavelet;

[0089] (2) Obtain the opposite number of the original shot record and use it as residual data;

[0090] (3) Calculate the background wave field and the back-propagation wave field, and separate the background wave field and the back-propagation wave field into upward, downward, left, and right traveling waves;

[0091] (4) Using the four directional traveling waves separated from the background wave field and the back propagation wave field, the low-frequency noise term is removed and a new imaging condition of vector decomposition is constructed;

[0092] (5) Calculate the iteration step size;

[0093] (6) Update imaging results;

[0094] (7) Obtain residual data and obtain the value of the target functional;

[0095] (8) Determine whether the termination condition is met. If yes, proceed to step (9); if no, return to step (3);

[0096] (9) Output the final imaging result.

[0097] Embodiment 2:

[0098] The operation of obtaining the opposite number of the original shot record in step (2) includes: directly taking the negative sign of the original shot record, or multiplying it by -1, to obtain the opposite number of the original shot record

[0099] Embodiment three:

[0100] The operations of calculating the background wave field and the back propagation wave field in step (3) include:

[0101] The background wave field is obtained by calculation using the following formula:

[0102]

[0103] Among them, u 0 (x,z;t) represents the background wave field, v 0 (x,z) represents the velocity field (the background velocity field directly input), f represents the seismic wavelet, t represents time, and x and z represent the directions of the spatial X and Z axes, respectively. Is a mathematical symbol representing the second-order derivative of space.

[0104] The back propagation wave field is calculated using the following formula:

[0105]

[0106] where u s (x,z;t) represents the back propagation wave field, v 0 (x,z) represents the velocity field, d res Represents residual data, t represents time, and x and z represent the X and Z axis directions of space respectively.

[0107] Embodiment 4:

[0108] In step (3), the Poynting vector is used to separate the upward, downward, left and right traveling waves of the background wave field and the back propagation wave field. The specific operations include:

[0109] The following formula is used to perform vector decomposition of the background wave field:

[0110]

[0111]

[0112]

[0113]

[0114]

[0115]

[0116] in, Indicates the up, down, left, and right traveling waves of the background wave field, u 0 (x,z;t) is the background wave field, Represent the horizontal and vertical components of the Poynting vector respectively.

[0117] The vector decomposition of the back propagation wave field is performed using the following formula:

[0118]

[0119]

[0120]

[0121]

[0122]

[0123]

[0124] in, Indicates the up, down, left, and right traveling waves of the background wave field, u s (x,z;t) is the back propagation wave field, Represent the horizontal and vertical components of the Poynting vector respectively.

[0125] Embodiment five:

[0126] The imaging condition in the conventional algorithm is as follows:

[0127]

[0128] Among them, u s 、u 0 They refer to the back-propagation wave field and the background wave field, respectively. 0 (x,z) represents the velocity field. After the wave field is decomposed in step (3), the above formula is completely equivalent to:

[0129]

[0130] Because when the background wave field and the back wave field have the same traveling direction, low-frequency noise is displayed, and when the background wave field and the back wave field have opposite directions, effective imaging results are displayed. Therefore, the present invention removes the low-frequency noise items with the same traveling direction and only retains the items with the opposite traveling direction. The imaging conditions for reconstructing the new vector decomposition in step (4) are as follows:

[0131]

[0132] Where t represents time, i represents the shot number, N represents the total number of shots in the seismic channel, and g represents the total number of shots in the seismic channel. k represents the gradient of the kth iteration, v 0 (x,z) represents the velocity field.

[0133] u in the above formula s 、u 0 They refer to the back propagation wave field and the background wave field respectively. Up, down, light and right represent the results calculated after decomposition in the up, down, left and right directions respectively. The above formula adds i to the formula of two vector decompositions, which represents the back propagation wave field and background wave field of the ith shot.

[0134] There are 8 items in the conventional algorithm, but only 4 items are retained in the present invention, and 4 low-frequency noise items are removed, thus suppressing the low-frequency noise in the gradient.

[0135] Embodiment six:

[0136] The operation of step (5) includes:

[0137] Using the conjugate gradient method, calculate the iteration step size:

[0138]

[0139] z k+1 =g k+1 +βz k

[0140]

[0141] Among them, g k represents the gradient of the kth iteration, L is the forward operator, α is the iteration step, β and z are the intermediate variables of the calculation, and T represents the transpose.

[0142]

[0143]

[0144] d lg =u 1 (x r ,z r ; t)

[0145] where d lg It represents the calculation result of Lg, that is, to obtain the specific position (x r ,z r ) and record it as the value of Lg. g represents the gradient, g k Represents the gradient of the kth iteration. During the iteration, g k Substitute g into the above equation; u 0(x, z; t) represents the background wave field, u 1 (x,z;t) represents the scattered wave field, v 0 (x,z) represents the velocity field, f represents the seismic wavelet, t represents time, x and z represent the X and Z axis directions of space respectively, r 、z r They respectively represent the detection point positions, and the detection point position information is read from the header of the input original shot record.

[0146]

[0147]

[0148] d lz =u 1 (x r ,z r ; t)

[0149] Among them, d lz It represents the calculation result of Lz, that is, obtaining the specific position (x r ,z r ) and record it as the value of Lz. z represents the intermediate result, z k Indicates the intermediate result at the kth iteration. During the iteration, z k Substitute z into the above equation. 0 (x, z; t) represents the background wave field, u 1 (x,z;t) represents the scattered wave field, v 0 (x,z) represents the velocity field, f represents the seismic wavelet, t represents time, x and z represent the X and Z axis directions of space respectively, r 、z r They respectively represent the detection point positions, and the detection point position information is read from the header of the input original shot record.

[0150] Embodiment seven:

[0151] The operation of step (6) includes:

[0152] Update the imaging results using the following formula:

[0153] m k+1 =m k -αg k (6)

[0154] Among them, m is the imaging result, α is the iteration step, g k represents the gradient of the kth iteration, where k is the number of iterations.

[0155] Embodiment eight:

[0156] The operation of obtaining the residual data in step (7) includes:

[0157] The residual data is the difference between the simulated data and the original shot record, and is obtained using the following formula:

[0158] d res =d cal -d obs (13)

[0159] Among them, d cal represents simulated data, d obs represents the original gun record, d res represents residual data;

[0160] d cal =Lm

[0161] Among them, L represents the forward operator, which represents a forward process. The formula is as follows, m represents the imaging result, d obs Represents the original gun record:

[0162]

[0163] Among them, d cal It represents the calculation result of Lm, that is, obtaining the specific position (x r ,z r ) and record it as the value of Lm. 0 (x, z; t) represents the background wave field, u 1 (x,z;t) represents the scattered wave field, v 0 (x,z) represents the velocity field, f represents the seismic wavelet, t represents time, x and z represent the X and Z axis directions of space respectively, d cal represents the calculated simulation data, m(x,z) represents the imaging result, x r 、z r They respectively represent the detection point positions, and the detection point position information is read from the header of the input original shot record.

[0164] The operation of obtaining the value of the target functional in step (7) includes:

[0165] Construct the target functional of the least squares reverse time migration based on the L-2 norm:

[0166] g min =||d cal -d obs || 2 (11)

[0167] Calculate the value of the target functional g min .

[0168] Embodiment nine:

[0169] The termination condition in step (8) is as follows:

[0170] g min The value of is less than the set threshold, or k is equal to the set number of iterations.

[0171] Embodiment ten:

[0172] The operation of step (9) includes:

[0173] Output the final least squares reverse time migration result, that is, the imaging result m after the last iteration k+1 .

[0174] Embodiment eleven:

[0175] The present invention will be further described below by taking the depression model as an example in combination with the accompanying drawings and specific implementation methods.

[0176] (1) Reading speed model v as Figure 1 As shown, and the original gun records, such as Figure 2 As shown;

[0177] (2) Finite difference is used to calculate the background wave field and the back propagation wave field, and the Poynting vector is used to separate the upward, downward, left and right traveling waves of the background wave field and the back propagation wave field respectively;

[0178] (3) Using the four directional traveling waves separated from the background wave field and the back propagation wave field, a new imaging condition for vector decomposition is constructed to solve the gradient;

[0179] (4) Using the conjugate gradient method, the iteration step length is calculated and the imaging results are updated;

[0180] (5) Determine whether the termination condition is met. The termination condition given here is g min If it is less than 0.0001 or the number of iterations is 40, the condition is met and the process goes directly to step (6). If the condition is not met, the difference between the simulated record and the original record is calculated and input into step (2) for iteration.

[0181] (6) Output the final least squares reverse time migration result, such as Figure 6 As shown;

[0182] Figure 4 The reverse time migration result is equivalent to the result of only one iteration of the conventional least squares migration. It can be seen that there is still a lot of noise in the imaging result, and the phase axis is not clear. Figure 5 and Figure 6 The results of conventional least squares reverse time migration and least squares reverse time migration based on vector decomposition imaging conditions after 40 point spread functions are shown. Figure 4 The reverse time migration results of both methods have a good noise suppression effect, but Figure 6 The suppression effect is more obvious, and the energy of the phase axis is more balanced. Figure 3 The theoretical value and Figure 5 and Figure 6 At the offset of 1500 meters, the amplitude curve is extracted, such as Figure 7 As shown ( Figure 7 The red line in the figure represents the least squares reverse time migration result of the vector decomposition imaging condition of the present invention, and the blue line represents the conventional least squares reverse time migration result). It can be seen that Figure 6 The results are more consistent with the theoretical results. Figure 8 It is a normalized curve diagram of the two methods. The convergence speed of the vector decomposition least squares reverse time migration is faster, indicating that the method proposed in the present invention can obtain better imaging effect and improve the convergence speed.

[0183] Least squares reverse time migration has great advantages in high-precision imaging in complex areas. The implementation method in the data domain is generally to construct a target functional and solve it iteratively. Then, due to the conventional cross-correlation imaging conditions, a large amount of low-frequency noise will be present in the calculated gradient, which will affect the convergence speed of the iteration and the final imaging effect, and affect its practical process. In response to this problem, the present invention constructs new imaging conditions by vector decomposing the back-propagation wave field and the background wave field, suppresses the low-frequency noise in the gradient, improves the imaging effect, and accelerates the convergence speed, which is conducive to the production application of the technology in actual work areas.

[0184] The present invention suppresses low-frequency noise in the gradient, improves imaging quality, and accelerates the convergence speed of iteration by changing the imaging conditions and modifying the gradient term in the least squares migration iteration process. Meanwhile, the increased amount of calculation is almost negligible, which is helpful for the subsequent application of the least squares reverse time migration technology and promotes the further development of the technology.

[0185] In the description of the present invention, it should be noted that, unless otherwise clearly specified and limited, the terms "connected" and "connection" should be understood in a broad sense, for example, it can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0186] In the description of the present invention, unless otherwise specified, the terms "upper", "lower", "left", "right", "inside", "outside", etc. indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings. They are only for the convenience of describing the present invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific direction, be constructed and operated in a specific direction. Therefore, they cannot be understood as limitations on the present invention.

[0187] The above technical solution is only one implementation mode of the present invention. For those skilled in the art, it is easy to make various types of improvements or modifications based on the principles disclosed in the present invention, and it is not limited to the technical solution described in the above specific embodiments of the present invention. Therefore, the above description is only preferred and does not have a restrictive meaning.

Claims

1. A least squares reverse time migration method, Features: The method first obtains the background wave field and the back propagation wave field, and performs vector decomposition to obtain the up, down, left and right traveling waves of the background wave field and the back propagation wave field, then performs gradient solution according to the vector decomposition imaging condition, suppresses low-frequency noise, and obtains the final imaging result by using the target functional of the least squares reverse time migration based on the L-2 norm.

2. The least squares reverse time migration method according to claim 1, Features: The method comprises: (1) Input velocity field, original shot record, and seismic wavelet; (2) Obtain the opposite number of the original shot record and use it as residual data; (3) Calculate the background wave field and the back-propagation wave field, and separate the background wave field and the back-propagation wave field into upward, downward, left, and right traveling waves; (4) Using the four directional traveling waves separated from the background wave field and the back-propagation wave field, a new imaging condition of vector decomposition is constructed; (5) Calculate the iteration step size; (6) Update imaging results; (7) Obtain residual data and obtain the value of the target functional; (8) Determine whether the termination condition is met. If yes, proceed to step (9); if no, return to step (3); (9) Output the final imaging result.

3. The least squares reverse time migration method according to claim 2, Features: The operations of calculating the background wave field and the back propagation wave field in step (3) include: The background wave field is obtained by calculation using the following formula: Among them, u 0 (x,z;t) represents the background wave field, v 0 (x,z) represents the velocity field, f represents the seismic wavelet, t represents time, and x and z represent the directions of the X and Z axes of space, respectively; The back propagation wave field is calculated using the following formula: where u s (x,z;t) represents the back propagation field, v 0 (x,z) represents the velocity field, d res Represents residual data, t represents time, and x and z represent the X and Z axis directions of space respectively.

4. The least squares reverse time migration method according to claim 3, Features: In step (3), the Poynting vector is used to separate the upward, downward, left and right traveling waves of the background wave field and the back propagation wave field. The specific operations include: The following formula is used to perform vector decomposition of the background wave field: in, Indicates the up, down, left, and right traveling waves of the background wave field, u 0 (x,z;t) is the background wave field, Respectively represent the horizontal and vertical components of the Poynting vector; The vector decomposition of the back propagation wave field is performed using the following formula: in, Indicates the up, down, left, and right traveling waves of the background wave field, u s (x,z;t) is the back propagation wave field, Represent the horizontal and vertical components of the Poynting vector respectively.

5. The least squares reverse time migration method according to claim 4, Features: The imaging conditions of the new vector decomposition in step (4) are as follows: Where t represents time, i represents the shot number, N represents the total number of shots in the seismic channel, and g represents the total number of shots in the seismic channel. k represents the gradient of the kth iteration, v 0 (x,z) represents the velocity field.

6. The least squares reverse time migration method according to claim 5, Features: The operation of step (5) includes: Using the conjugate gradient method, calculate the iteration step size: z k+1 =g k+1 +βz k Among them, g k represents the gradient of the kth iteration, L is the forward operator, α is the iteration step size, β and z are intermediate variables.

7. The least squares reverse time migration method according to claim 6, Features: The operation of step (6) includes: Update the imaging results using the following formula: m k+1 =m k -αg k Among them, m is the imaging result, α is the iteration step, g k represents the gradient of the kth iteration, where k is the number of iterations.

8. The least squares reverse time migration method according to claim 7, Features: In step (7), the following formula is used to obtain the residual data: d res =d cal -d obs Among them, d cal represents simulated data, d cal =Lm, L represents the forward operator, d obs represents the original gun record, d res represents residual data; The operation of obtaining the value of the target functional in step (7) includes: Construct the target functional of the least squares reverse time migration based on the L-2 norm: g min =||d cal -d obs || 2 Calculate the value of the target functional g min .

9. The least squares reverse time migration method according to claim 8, Features: The termination condition in step (8) is as follows: g min The value of is less than the set threshold, or k is equal to the set number of iterations.

10. A computer-readable storage medium, It is characterized in that The computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps in the least squares reverse time migration method according to any one of claims 1 to 9.

Images