Gradient preconditioned image domain least squares reverse time migration method

Abstract

The application discloses a gradient precondition-based image domain least square reverse time migration method and belongs to the field of deep-sea oil and gas field exploration and development. The method aims to solve the problems of imaging amplitude distortion, resolution decline and slow convergence caused by the traditional gradient optimization method due to the neglect of the second-order curvature information of the Hessian matrix. The technical scheme points are as follows: an image domain least square reverse time migration objective function is constructed; the Hessian matrix is approximated based on a point spread function; a training pair of the gradient and the migration gradient is constructed by using the Born scattering theory; an improved U-Net network with a channel-space double attention mechanism is built to learn the mapping relationship of the Hessian inverse matrix in an unsupervised manner; and the trained network is embedded into the optimization framework of a fast iterative shrinkage threshold algorithm as a gradient precondition operator to correct the model updating direction. The application effectively improves the imaging quality of uneven illumination areas and improves the amplitude fidelity and inversion convergence speed.

Description

Technical Field

[0001] This invention belongs to the field of deep-sea oil and gas field exploration and development technology, and particularly relates to an image domain least squares reverse time migration method based on gradient preconditions. Background Technology

[0002] In the field of seismic exploration for deep-sea oil and gas fields, reverse time migration (RTM) has become a key technology for high-precision imaging of complex structures. To further improve imaging resolution, image-domain least-squares RTM has been extensively studied. Existing typical techniques first obtain the full-space point spread function, which reflects spatial variation characteristics, using numerical methods to approximate the Hessian matrix characterizing the blurring effect of the imaging system. This allows for the construction of a forward and inverse framework describing the imaging process. For the resulting ill-conditioned inverse problem, gradient-based optimization algorithms are typically used for efficient solution, aiming to improve imaging quality at a computational cost close to that of conventional RTM.

[0003] However, this existing technical solution still has significant drawbacks in practical applications. Traditional gradient-based optimization methods only utilize the first-order gradient information of the objective function during model updates, completely ignoring the second-order curvature information contained in the Hessian matrix. In complex deep-sea geological structures, this neglect makes it difficult for the method to effectively compensate for the complex wavefield propagation paths and severe illumination inhomogeneity caused by salt domes, steeply dipping faults, and high-speed shielding layers. The direct consequences are distortion of imaging amplitude, discontinuity of phase axes, and decreased spatial resolution in deep and structural boundary regions. Furthermore, optimization paths relying solely on first-order information are often not optimal, resulting in slow convergence speed and a large number of iterations required for the entire inversion process, thus affecting imaging efficiency.

[0004] While theoretically, explicitly constructing and storing the Hessian inverse matrix and applying it to the gradient can introduce second-order curvature information to correct the aforementioned defects, it faces insurmountable obstacles in practical large-scale 3D seismic imaging tasks. The explicit computational complexity of the full-dimensional Hessian matrix and its inverse is extremely high, requiring computational resources that grow superlinearly, while its storage requirements reach massive levels, far exceeding the capabilities of current commercial computing hardware. This high computational and storage cost makes preconditioning techniques based on explicit Hessian inverse matrices completely impractical for industrial production. Therefore, finding a technical solution in seismic imaging that can efficiently introduce second-order optimization information while keeping computational costs under control has long been a major challenge for researchers. Summary of the Invention

[0005] To address the aforementioned technical problems, this invention proposes an image domain least squares reverse time migration method based on gradient preconditions, thereby resolving the issues present in the prior art.

[0006] Firstly, to achieve the above objectives, the present invention provides an image domain least-squares reverse time migration method based on gradient preconditioning, comprising the following steps:

[0007] Construct an image domain least-squares reverse time migration objective function based on the Hessian matrix;

[0008] The Hessian matrix is ​​approximated based on the point spread function;

[0009] Based on Born scattering theory and an approximate Hessian matrix, gradient data pairs are constructed for network training.

[0010] A deep neural network is built and trained to learn the mapping relationship from the offset gradient to the original gradient, thereby obtaining a surrogate model of the Hessian inverse matrix.

[0011] The objective function is solved iteratively using a fast iterative shrinkage threshold algorithm;

[0012] In each iteration of the model update step, the surrogate model is used as a preconditioning operator on the gradient of the objective function to correct the model update direction.

[0013] Optionally, the process of constructing the image domain least-squares reverse time migration objective function includes:

[0014] Establish a linear forward modeling relationship between observational data and reflectance models;

[0015] The observed data is migrated and imaged using the adjoint operator of the forward modeling operator to obtain the migrated image;

[0016] An imaging model is established that correlates the offset image with the true reflectance model via a Hessian matrix;

[0017] Based on the imaging model, an objective function is constructed that aims to minimize the L2 norm of the difference between the offset image and the Hessian matrix acting on the reflectivity model.

[0018] Optionally, the process of approximating the Hessian matrix based on the point spread function includes:

[0019] Sparse control points are set up in the model space, and point scatterers are set at the control points;

[0020] Calculate the point spread function at each control point, the point spread function characterizing the offset response of the imaging system to the scatterer at that point;

[0021] Based on the point spread function at all control points, the point spread function distribution across the entire space is reconstructed through interpolation;

[0022] The complete Hessian matrix is ​​approximated by using the reconstructed full-space point spread function distribution.

[0023] Optionally, the process of constructing gradient data pairs for network training includes:

[0024] During the iteration process, obtain the gradient of the objective function corresponding to the current reflectivity model;

[0025] Using the point spread function obtained from Born inverse offset-offset calculation, a non-stationary convolution operation is performed on the gradient of the objective function to obtain a simulated offset gradient image;

[0026] The simulated offset gradient image is used as the network input, and the gradient of the objective function is used as the label of the network's expected output to form a training data pair.

[0027] Optionally, the deep neural network is an improved U-Net network with an encoder-decoder structure;

[0028] The improved U-Net network includes dual attention modules with both channel attention and spatial attention mechanisms in both the encoder and decoder.

[0029] Optionally, the process of training the deep neural network adopts an unsupervised learning method;

[0030] The loss function used for training is the L2 error between the gradient predicted by the network and the gradient of the objective function used as the label.

[0031] Optionally, the iterative solution process using the fast iterative shrinking threshold algorithm includes:

[0032] Introduce the L1 norm regularization term into the objective function;

[0033] In each iteration, the gradient of the objective function with respect to the current model is calculated;

[0034] A shrinkage threshold operation is applied to the intermediate variables after gradient descent and preconditioning to promote the sparsity of the solution;

[0035] The iteration step size and reflectivity model are updated according to the algorithm rules.

[0036] Optionally, the process of using the surrogate model as a preconditioning operator in the model update step is as follows:

[0037] The output of the trained deep neural network is used as an approximation of the gradient of the objective function after filtering with the Hessian inverse matrix.

[0038] The approximate result is multiplied by the iteration step size and then used to update the current reflectivity model.

[0039] In a second aspect, the present invention also provides a computer terminal device, comprising:

[0040] One or more processors;

[0041] A memory, coupled to the processor, for storing one or more programs;

[0042] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the gradient-preconditioned image domain least-squares inverse time migration method in the first aspect described above.

[0043] Thirdly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the computer program is executed by a processor, it implements the steps of the gradient-preconditioned image domain least squares reverse time migration method described in the first aspect above.

[0044] Compared with the prior art, the present invention has the following advantages and technical effects:

[0045] This invention provides an image-domain least-squares reverse time migration method based on gradient preconditioning. By introducing a deep learning mechanism to implicitly approximate the Hessian inverse matrix, it effectively overcomes the problem of insufficient adaptability of traditional methods to complex geological structures. It significantly improves the imaging amplitude distortion and resolution reduction defects in deep and boundary regions caused by uneven wavefield illumination, effectively maintaining the integrity and continuity of the imaging phase axis. By embedding the trained network as a gradient preconditioning operator into a fast iterative optimization process, second-order curvature information is introduced to correct the model update path, thereby significantly accelerating the inversion convergence speed and enhancing the stability of the iterative process. Furthermore, the physics-driven unsupervised learning framework avoids the dependence on large amounts of labeled data in traditional supervised learning. While ensuring improved imaging quality, it significantly reduces the computational complexity and cost in practical applications, providing an efficient and feasible technical path for large-scale seismic imaging tasks. Attached Figure Description

[0046] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0047] Figure 1 This is a flowchart illustrating an image domain least squares reverse time migration method based on gradient preconditioning according to an embodiment of the present invention.

[0048] Figure 2 This is a schematic diagram of the improved U-Net network structure according to an embodiment of the present invention;

[0049] Figure 3This is a horizontal velocity model for an embodiment of the present invention. The model size is 256×192, the grid spacing is 10 meters, and the velocity range is 2000-3000m / s.

[0050] Figure 4 The PSF generated based on the acoustic wave equation and the corresponding Born scattering equation in this embodiment of the invention. A PSF B Schematic diagram, where (a) is PSF A Schematic diagram, (b) shows PSF B Schematic diagram;

[0051] Figure 5 This is a schematic diagram of the gradient, offset gradient, network prediction gradient, and network correction gradient in the DL-IDLSRTM process of this invention, where (a) is the gradient, (b) is the offset gradient, (c) is the network prediction gradient, and (d) is the network correction gradient.

[0052] Figure 6 The loss function of the network during the training process in this embodiment of the invention. and IDLSRTM objective function The convergence curve is shown in the figure, where (a) is the loss function of the network. (b) is the objective function The convergence curve;

[0053] Figure 7 The following are schematic diagrams of imaging results obtained by different methods in embodiments of the present invention, wherein (a) is the reflectance model, (b) is the result of reverse time migration (RTM), (c) is the result of image domain least squares reverse time migration (IDLSRTM), and (d) is the result of gradient preconditioned image domain least squares reverse time migration (DL-IDLSRTM).

[0054] Figure 8 The following are schematic diagrams of vertical offset profiles at positions of 1km and 1.9km for different imaging results of embodiments of the present invention, where (a) is a vertical offset profile at position 1km and (b) is a vertical offset profile at position 1.9km. Detailed Implementation

[0055] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other. The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0056] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0057] This embodiment provides an image domain least squares reverse time migration method based on gradient preconditions, including: constructing an image domain least squares reverse time migration objective function based on the Hessian matrix;

[0058] The Hessian matrix is ​​approximated based on the point spread function;

[0059] Based on Born scattering theory and an approximate Hessian matrix, gradient data pairs are constructed for network training.

[0060] A deep neural network is built and trained to learn the mapping relationship from the offset gradient to the original gradient, thereby obtaining a surrogate model of the Hessian inverse matrix.

[0061] The objective function is solved iteratively using a fast iterative shrinkage threshold algorithm;

[0062] In each iteration of the model update step, the surrogate model is used as a preconditioning operator on the gradient of the objective function to correct the model update direction.

[0063] Furthermore, the process of constructing the image domain least-squares reverse time migration objective function includes:

[0064] Establish a linear forward modeling relationship between observational data and reflectance models;

[0065] The observed data is migrated and imaged using the adjoint operator of the forward modeling operator to obtain the migrated image;

[0066] An imaging model is established that correlates the offset image with the true reflectance model via a Hessian matrix;

[0067] Based on the imaging model, an objective function is constructed that aims to minimize the L2 norm of the difference between the offset image and the Hessian matrix acting on the reflectivity model.

[0068] Specifically, the implementation process of this embodiment includes:

[0069] Construct the IDLSRTM objective function:

[0070] In the frequency domain, suppose the excitation source is located at... Its source signal is Then at the location of the detector point The observed reflection seismic records at this location are expressed as follows:

[0071] (1);

[0072] in, It is the imaginary unit. Angular frequency, For the imaging point in space, For reflectivity model, For the source wavelet, and These represent the locations of the earthquake source and receiver, respectively. Describes the epicenter To the imaging point The incident wave field, Then describe from the imaging point to the receiving point The reflected wave field, Indicates that it is caused by the epicenter Excitation, at the receiving point The frequency domain seismic wavefield data recorded at [location]. The model is represented in compact matrix form:

[0073] (2);

[0074] Where d is the observed seismic data, L is the linear forward modeling operator, and m is the reflectivity model.

[0075] However, in actual imaging processes, due to limitations in computational complexity and practical observation systems, conventional migration imaging typically employs the adjoint operator of the forward modeling operator rather than the inverse operator. The imaging results are expressed as follows:

[0076] (3);

[0077] in It is a time-shifted image. It is the adjoint operator of the linear forward modeling operator. Substituting equation (2) into equation (3), we obtain the reverse time offset. With underground reflectivity model The relationship between them:

[0078] (4);

[0079] in, The Hessian matrix characterizes the blurring effect on reverse-time migrated images caused by factors such as limited bandwidth seismic sources, acquisition system configuration, and background velocity. Therefore, conventional migration imaging results can be considered as the true reflectance after being processed by the Hessian matrix. To eliminate the blurring effect of the Hessian matrix and obtain higher resolution imaging results, this embodiment establishes the following objective function:

[0080] (5).

[0081] Furthermore, the process of approximating the Hessian matrix based on the point spread function includes:

[0082] Sparse control points are set up in the model space, and point scatterers are set at the control points;

[0083] Calculate the point spread function at each control point, the point spread function characterizing the offset response of the imaging system to the scatterer at that point;

[0084] Based on the point spread function at all control points, the point spread function distribution across the entire space is reconstructed through interpolation;

[0085] The complete Hessian matrix is ​​approximated by using the reconstructed full-space point spread function distribution.

[0086] Specifically, the implementation process of this embodiment includes:

[0087] Approximate Hessian matrix:

[0088] The exact Hessian matrix cannot be directly obtained and applied due to its huge dimension and computational cost.

[0089] Therefore, this embodiment proposes a Hessian approximation strategy based on PSF. PSF characterizes the offset response of the imaging system to a spatial unit point scatterer, and its mathematical essence corresponds to the column vector (i.e., local kernel) describing that spatial location in the Hessian matrix. It fully characterizes the imaging blurring effect jointly introduced by the finite bandwidth source, the non-uniform observation system, and the velocity non-uniformity of the background medium, reflecting the spatial resolution and local illumination characteristics of the imaging system, thus providing a physically clear solution path for efficiently approximating the Hessian matrix. In the frequency domain, the Hessian matrix integral kernel in the model space is expressed as:

[0090] (6);

[0091] in, This indicates taking the real part of a complex number. and These are the locations of the seismic source and the receiver, respectively. Angular frequency; To enhance weights at higher frequencies, The spectrum of the source wavelet. Indicates its energy. and These are two locations in the imaging space; Green's function is the function used in an acoustic medium to describe the propagation response of a wave between two points. This is the conjugate of the Green's function. This equation has a dual physical meaning: firstly, it describes the scattering point... For imaging points The contribution weight; secondly, when fixed When the wave field is scattered at a point, the equation degenerates into the PSF, which describes the system's scattering response at that point. The four Green's function terms correspond to the forward and reverse propagation paths of the wave field, respectively, collectively reflecting the wave field's ability to focus on the subsurface scattering point. Based on this theory, this invention achieves efficient Hessian matrix approximation through the following steps:

[0092] (1) Sparse control point layout: Point scatterers are uniformly distributed in the offset model space to construct the scatterer model;

[0093] (2) PSF calculation: The local PSF at each control point is obtained by calculating the offset response of the scatterer model;

[0094] (3) Interpolation and reconstruction: For any non-control point location, interpolation is performed using the PSF at the nearest control point to quickly obtain the PSF distribution across the entire space.

[0095] Furthermore, the process of constructing gradient data pairs for network training includes:

[0096] During the iteration process, obtain the gradient of the objective function corresponding to the current reflectivity model:

[0097] Using the point spread function obtained from Born inverse offset-offset calculation, a non-stationary convolution operation is performed on the gradient of the target function to obtain a simulated offset gradient image;

[0098] The simulated offset gradient image is used as the network input, and the gradient of the objective function is used as the label of the network's expected output to form a training data pair.

[0099] Furthermore, the process of training the deep neural network adopts an unsupervised learning method;

[0100] The loss function used for training is the L2 error between the gradient predicted by the network and the gradient of the objective function used as the label.

[0101] Specifically, the implementation process of this embodiment includes:

[0102] Based on Born scattering theory, a training pair of gradient and offset gradient is constructed:

[0103] The key to constructing an efficient surrogate model that can accurately approximate the Hessian inverse matrix lies in obtaining training samples that can fully characterize the spatial fuzzy properties of the Hessian matrix. Based on the wave field linearization propagation mechanism revealed by Born scattering theory, this invention proposes to construct "input-output" training data pairs with clear physical meaning by sequentially performing inverse and offset operations on the current gradient.

[0104] Offset gradient image expression for:

[0105] (7);

[0106] Among them, the linear forward modeling operator Defined as:

[0107] (8);

[0108] In the formula, These are spatial coordinates. It is time. It is the location of the epicenter. It's the background speed. It is a disturbance term. It is the scattered wave field caused by the perturbation term. It is the background wave field. It is the Laplacian operator. Since each iteration requires an inverse offset-offset process to calculate the effect of the Hessian matrix on the gradient, the computational cost is high. Therefore, this embodiment proposes an approximation method to model the offset gradient image as:

[0109] (9);

[0110] in, The point spread function is calculated based on Born's inverse offset-offset. This is a non-stationary two-dimensional convolution operation. Once the model is established... arrive Given the end-to-end mapping relationship, deep neural networks (DNNs) can serve as efficient surrogate models for the Hessian inverse matrix. The training process is formalized as an optimization problem, with its loss function defined as:

[0111] (10).

[0112] in, Indicates the learnable parameters Parameterized proxy model.

[0113] Furthermore, the deep neural network is an improved U-Net network with an encoder-decoder structure;

[0114] The improved U-Net network includes dual attention modules with both channel attention and spatial attention mechanisms in both the encoder and decoder.

[0115] Specifically, the implementation process of this embodiment includes:

[0116] An improved U-Net network with a dual channel-space attention mechanism was constructed;

[0117] This invention constructs the task of learning the Hessian inverse matrix as a nonlinear mapping problem from an offset gradient image to a gradient image, which is implemented by a convolutional neural network (CNN) based on the U-Net architecture capable of simulating such complex transformations. Figure 2 As shown, the network adopts an encoder-decoder structure. The encoding stage is used to capture multi-scale features, while the decoding stage is dedicated to reconstructing a high-fidelity target image.

[0118] The encoder section comprises five cascaded convolutional blocks, each consisting of two consecutive 3×3 convolutional layers (Conv), batch normalization (BatchNorm), dropout regularization, and a Leaky-ReLU activation function. The negative slope of Leaky-ReLU is set to 0.1 to preserve crucial negative phase information in the input features. A dropout method with a probability of 0.1 is used during training to mitigate overfitting. To further enhance feature representation capabilities, each module's output is refined using a dual attention mechanism: the Channel Attention Module (CAM) calculates channel weights through adaptive pooling and a shared multilayer perceptron (MLP); the Spatial Attention Module (SAM) focuses attention by concatenating the mean and maximum features at the channel level, followed by 3×3 convolution and element-wise multiplication. After attention refinement, each encoder block undergoes downsampling via 2×2 max pooling.

[0119] The decoder consists of five convolutional blocks, mirroring the encoder's structure. Each block first upsamples the feature map using nearest-neighbor interpolation, then concatenates it with high-resolution features from the corresponding encoder layer to preserve rich spatial details. The concatenated features are further processed by the convolutional blocks, progressively compressing the channel dimensions and restoring spatial resolution. Finally, a separate convolutional layer generates an output prediction that matches the target dimension; this output is the network's approximate representation of the Hessian inverse matrix.

[0120] Furthermore, the iterative solution process using the fast iterative shrinkage threshold algorithm includes:

[0121] Introduce the L1 norm regularization term into the objective function;

[0122] In each iteration, the gradient of the objective function with respect to the current model is calculated;

[0123] A shrinkage threshold operation is applied to the intermediate variables after gradient descent and preconditioning to promote the sparsity of the solution;

[0124] The iteration step size and reflectivity model are updated according to the algorithm rules.

[0125] Specifically, the implementation process of this embodiment includes:

[0126] Solving the inverse problem based on the FISTA algorithm;

[0127] Based on formula (5), to stabilize the solution to the inverse problem, L1 norm regularization can be introduced to constrain the objective function. This leverages the sparse representation characteristics of reflectance images to effectively improve the stability and interpretability of the solution. The constructed regularized objective function expression is as follows:

[0128] (11);

[0129] in, For regularization parameters, It is the L1 norm of m. Within this framework, the update formula before shrinkage thresholding is:

[0130] (8);

[0131] in, It is the first Model update amount in each iteration To update the step size, k is the iteration index. This represents the gradient of the objective function at the k-th iteration.

[0132] In each iteration, the sparse solution Generated by the shrinkage threshold operator, its expression is:

[0133] (12);

[0134] in, This is a non-negative threshold used to control the sparsity of the iteration results; It is a sparse solution; This represents the shrinkage threshold operator; For symbolic functions, This means that each component is compared with... Compare, if less than Set to zero if necessary, otherwise subtract... The step size is calculated using the following formula:

[0135] (13);

[0136] The reflectivity model is updated as follows:

[0137] (14);

[0138] The initial condition for the iteration is set to , Threshold Defined as ,in for The reciprocal of the largest eigenvalue, The parameter is used to control sparsity and adjust the sparsity level of the solution.

[0139] Furthermore, the process of using the surrogate model as a preconditioning operator in the model update step is as follows:

[0140] The output of the trained deep neural network is used as an approximation of the gradient of the objective function after filtering with the Hessian inverse matrix.

[0141] The approximate result is multiplied by the iteration step size and then used to update the current reflectivity model.

[0142] Specifically, the implementation process of this embodiment includes:

[0143] A gradient preconditioner is constructed by approximating the Hessian inverse matrix using a network.

[0144] During the iterative inversion process, in order to find the objective function in the current model The update direction at that point is derived based on its truncated second-order Taylor expansion:

[0145] (15);

[0146] in, For the objective function in gradient at, Represented as the corresponding Hessian matrix, This represents the model update amount. For formula (15)... Finding the first variation, we get:

[0147] (16);

[0148] when When taking the minimum value:

[0149] (17);

[0150] This relationship indicates that the gradient is essentially the result of the model update amount after filtering with the Hessian inverse matrix. Therefore, substituting equation (14) into the iterative equation, the solution after Hessian inverse matrix preprocessing is expressed as:

[0151] (18).

[0152] Since the Hessian matrix characterizes the local curvature properties of the objective function, its inverse matrix can compensate for amplitude distortion caused by geometric diffusion and wave propagation transmission loss, significantly improving the spatial resolution and amplitude fidelity of the reconstructed image. However, explicitly constructing and storing the complete Hessian inverse matrix incurs extremely high computational and storage costs, making it difficult to implement in practical applications. Therefore, this invention employs a data-adaptive surrogate model based on a deep neural network (DNN) to achieve an efficient approximation of the Hessian inverse matrix through unsupervised training. Based on this surrogate model, the pre-threshold solution can be calculated using the following formula:

[0153] (19);

[0154] in, Indicates the learnable parameters A parameterized surrogate model. This network does not rely on manually labeled training data, but is trained using gradient data implicitly generated during the DL-IDLSRTM iteration process. The input data is constructed by first inversely shifting the gradients and then performing the shift operation, thereby effectively capturing the response characteristics of the imaging system.

[0155] In this embodiment, a computer terminal device is provided, including:

[0156] One or more processors;

[0157] A memory, coupled to the processor, for storing one or more programs;

[0158] When the one or more programs are executed by the one or more processors, the one or more processors implement the steps of the above-described gradient-preconditioned image domain least-squares inverse time migration method.

[0159] In this embodiment, a computer-readable storage medium is also provided, on which a computer program is stored. When the computer program is executed by a processor, it implements the steps of the above-described gradient-preconditioned image domain least-squares inverse time-shifting method.

[0160] An application example of this invention is as follows:

[0161] The workflow and improved network structure of this embodiment are respectively... Figure 1 and Figure 2 Displayed in China. Figure 3The experiment used a velocity model as an example, with a grid size of 256×192 and velocities ranging from 2000 to 3000 m / s. Synthetic seismic data was generated from 10 seismic sources spaced 200 m apart, with detectors evenly distributed along the surface at 10 m intervals. The seismic sources used a Ricker wavelet with a dominant frequency of 25 Hz, and the total recording time was 3 seconds. The experiment employed a finite difference scheme with second-order temporal precision and twelfth-order spatial precision to achieve seismic wavefield extension. Perfectly matched layer (PML) absorbing boundary conditions were applied to suppress spurious boundary reflections. All synthetic seismic data were computed on a uniform grid with a spatial sampling interval of 10 m and a temporal sampling interval of 1 ms. The training of the convolutional neural network (CNN) was completed within the PyTorch framework, with accelerated computation utilizing an NVIDIA RTX 2080 Ti GPU.

[0162] The experiment first transformed the calculation of the Hessian matrix into a process of calculating the PSF at the control points and performing fast interpolation. Due to differences in the theoretical formulas, two different PSFs were obtained. Figure 4 The PSF generated based on the acoustic wave equation and the corresponding Born scattering equation in this embodiment of the invention. A PSF B Schematic diagram, where (a) is PSF A Schematic diagram, (b) shows PSF B Schematic diagram. The initial gradient image is shown below. Figure 5 As shown in (a). Figure 5 (b) shows the offset gradient, which is achieved through... Figure 4 PSF in (b) B The results generated by non-stationary convolution with the initial gradient show the blurring effect in the deep region caused by the inverse-offset process.

[0163] During network training, gradient images were used as the ground truth labels, and offset gradient images were used as the input data. The network parameters were optimized using the Adam optimizer, with an initial learning rate of 0.001 and a loss function defined by formula (10). Training lasted for 200 epochs. The gradients predicted by the network were as follows: Figure 5 As shown in (c), the high consistency between the two indicates that the network successfully captured the key features of the Hessian inverse matrix. The network's training convergence curve is shown in Figure 1. Figure 6 As shown in (a), the loss function rapidly decreases from 0.8 to 0.1 within the first 75 cycles, demonstrating the model's efficient optimization capability in the initial stage. After approximately 150 cycles, the loss enters a stable phase with negligible oscillations, indicating that the model has essentially converged to a stable solution. Due to the adoption of a parameter inheritance strategy, the initial state of the network is optimized through continuous iterations; therefore, as the number of DL-IDLSRTM iterations increases, the initial loss value gradually decreases.

[0164] After network training is complete, the gradient is input into the trained model for forward propagation to obtain the predicted gradient, such as... Figure 5 As shown in (d), the effect of equalization amplitude compensation from shallow to deep regions is demonstrated, which corrects the attenuation problem with depth in the original gradient. Figure 6 (b) compares the convergence curves of IDLSRTM and DL-IDLSRTM. Both show a decreasing trend, but DL-IDLSRTM, which uses gradient preprocessing, achieves higher convergence accuracy. The reflectivity model is as follows... Figure 7 As shown in (a). Compared to RTM ( Figure 7 Compared to (b), IDLSRTM ( Figure 7 While (c) reduces offset artifacts and wavelet sidelobe effects, it results in some loss of effective reflection information due to uneven gradient energy distribution and shrinkage threshold processing. DL-IDLSRTM ( Figure 7 (d) not only retains the advantages of traditional methods, but also significantly improves the imaging quality and amplitude fidelity in weakly illuminated areas. These conclusions are in... Figure 8 This was further verified in the vertical offset profile, where Figure 8 (a) is a vertical offset profile at a location of 1 km. Figure 8 (b) shows the vertically offset profile at a location of 1.9 km. Overall, the improved DL-IDLSRTM provides higher quality imaging results.

[0165] This invention provides an image-domain least-squares reverse time migration method based on gradient preconditioning. By introducing a deep learning mechanism to implicitly approximate the Hessian inverse matrix, it effectively overcomes the problem of insufficient adaptability of traditional methods to complex geological structures. It significantly improves the imaging amplitude distortion and resolution reduction defects in deep and boundary regions caused by uneven wavefield illumination, effectively maintaining the integrity and continuity of the imaging phase axis. By embedding the trained network as a gradient preconditioning operator into a fast iterative optimization process, second-order curvature information is introduced to correct the model update path, thereby significantly accelerating the inversion convergence speed and enhancing the stability of the iterative process. Furthermore, the physics-driven unsupervised learning framework avoids the dependence on large amounts of labeled data in traditional supervised learning. While ensuring improved imaging quality, it significantly reduces the computational complexity and cost in practical applications, providing an efficient and feasible technical path for large-scale seismic imaging tasks.

[0166] The above are merely preferred embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A gradient-preconditioned image domain least-squares reverse-time migration method, characterized in that, include: Construct an image domain least-squares reverse time migration objective function based on the Hessian matrix; The Hessian matrix is ​​approximated based on the point spread function; Based on Born scattering theory and an approximate Hessian matrix, gradient data pairs are constructed for network training. Build and train a deep neural network to learn the mapping relationship from the offset gradient to the original gradient, thereby obtaining a surrogate model of the Hessian inverse matrix; The objective function is solved iteratively using a fast iterative shrinkage threshold algorithm; In each iteration of the model update step, the surrogate model is used as a preconditioning operator to apply to the gradient of the objective function in order to correct the model update direction; The process of constructing the image domain least-squares reverse time migration objective function includes: Establish a linear forward modeling relationship between observational data and reflectivity models; The observation data represents data from the earthquake source. Excitation, at the receiving point Frequency domain seismic wavefield data recorded at the location; (1); where d is the observed data, L is the linear forward modeling operator, and m is the reflectance model; The adjoint operator of the linear forward modeling operator is used to perform migration imaging on the observed data to obtain a reverse time-migrated image; (2); among which It is a time-shifted image. It is the adjoint operator of linear forward operators; An imaging model is established that links the reverse-time migrated image and the reflectance model through the Hessian matrix; substituting equation (1) into equation (2), the reverse-time migrated image is obtained. With reflectivity model The relationship between them: (3); among them, It is a Hessian matrix; Based on the imaging model, an objective function is constructed that aims to minimize the L2 norm of the difference between the reverse time-shifted image and the Hessian matrix acting on the reflectivity model.

2. The method according to claim 1, characterized in that, The process of approximating the Hessian matrix based on the point spread function includes: Sparse control points are set up in the model space, and point scatterers are set at the control points; Calculate the point spread function at each control point, the point spread function characterizing the offset response of the imaging system to the scatterer at that point; Based on the point spread function at all control points, the distribution of the point spread function in the entire space is reconstructed by interpolation; The complete Hessian matrix is ​​approximated by using the reconstructed full-space point spread function distribution.

3. The method according to claim 1, characterized in that, The process of constructing gradient data pairs for network training includes: During the iteration process, obtain the gradient of the objective function corresponding to the current reflectivity model; Using the point spread function obtained from Born inverse offset-offset calculation, a non-stationary convolution operation is performed on the gradient of the target function to obtain a simulated offset gradient image; The simulated offset gradient image is used as the network input, and the gradient of the objective function is used as the label of the network's expected output to form a training data pair.

4. The method according to claim 1, characterized in that, The deep neural network is an improved U-Net network with an encoder-decoder structure; The improved U-Net network includes dual attention modules with both channel attention and spatial attention mechanisms in both the encoder and decoder.

5. The method according to claim 1, characterized in that, The process of training the deep neural network adopts an unsupervised learning method; The loss function used for training is the L2 error between the gradient predicted by the network and the gradient of the objective function used as the label.

6. The method according to claim 1, characterized in that, The iterative solution process using the fast iterative shrinking threshold algorithm includes: Introduce the L1 norm regularization term into the objective function; In each iteration, the gradient of the objective function with respect to the current model is calculated; A shrinkage threshold operation is applied to the intermediate variables after gradient descent and preconditioning to promote the sparsity of the solution; The iteration step size and reflectivity model are updated according to the algorithm rules.

7. The method according to claim 1, characterized in that, The process of using the surrogate model as a preconditioning operator in the model update step is as follows: The output of the trained deep neural network is used as an approximation of the gradient of the objective function after filtering with the Hessian inverse matrix. The approximate result is multiplied by the iteration step size and then used to update the current reflectivity model.

8. A computer terminal device, characterized in that, include: One or more processors; A memory, coupled to the processor, for storing one or more programs; When the one or more programs are executed by the one or more processors, the one or more processors perform the steps of the method as described in any one of claims 1-7.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the steps of the method as described in any one of claims 1-7.

Images