A seismic imaging method, system, terminal, and storage medium based on an approximate Hessian array.

Abstract

This invention discloses a seismic imaging method, system, terminal, and storage medium based on an approximate Hessian array. The method includes: acquiring seismic data and migration parameters, and performing pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and a dip field; determining the target imaging region, dividing it into multiple one-dimensional subspaces, acquiring the target parameters and cutoff frequency in each one-dimensional subspace, and generating frequency integral results based on the target parameters, cutoff frequency, and dip field; determining a preset approximate Hessian array analytical expression, and inputting the frequency integral results, cutoff frequency, and target parameters into the preset approximate Hessian array analytical expression to obtain the target approximate Hessian array; constructing a linear inversion equation system based on the conventional migration results and the target approximate Hessian array, and solving it to obtain the seismic imaging results. This invention can significantly improve the efficiency and accuracy of least-squares migration seismic imaging in the imaging domain.

Description

Technical Field

[0001] This invention relates to the field of seismic imaging processing technology, and in particular to a seismic imaging method, system, terminal, and computer-readable storage medium based on an approximate Hessian array. Background Technology

[0002] In seismic imaging processing, the conventional migration operator is used as an adjoint operator of the forward modeling operator rather than an inverse operator. Under non-ideal conditions such as incomplete observation system, uneven sampling, or limited coverage, it is prone to causing a series of imaging defects such as amplitude distortion, migration arcing, and uneven illumination.

[0003] To address the aforementioned issues, the imaging domain least squares migration method was proposed. Its core lies in explicitly constructing and storing the Hessian array, and then compensating for imaging errors through the inversion process, recovering true amplitude information, and effectively suppressing migration noise.

[0004] However, the imaging results obtained by the existing least squares migration method are computationally inefficient and have poor imaging quality, which cannot meet the needs of seismic imaging.

[0005] Therefore, existing technologies still need to be improved and developed. Summary of the Invention

[0006] The main objective of this invention is to provide a seismic imaging method, system, terminal, and computer-readable storage medium based on an approximate Hessian array, aiming to solve the problems of low computational efficiency and poor imaging results obtained by the current least squares migration method, which cannot meet the requirements of seismic imaging.

[0007] To achieve the above objectives, the present invention provides a seismic imaging method based on an approximate Hessian array, the seismic imaging method based on an approximate Hessian array comprising the following steps: Seismic data and migration parameters are acquired, and the seismic data are processed by pre-stack time migration according to the migration parameters to obtain conventional migration results and dip field. The target imaging region is determined, and the target imaging region is divided into multiple one-dimensional subspaces. The target parameters and cutoff frequencies in each one-dimensional subspace are obtained, and frequency integral results are generated based on the target parameters, the cutoff frequencies, and the tilt field. Determine a preset approximate Hessian matrix analytical expression, and input the frequency integral result, the cutoff frequency and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace; Based on the conventional migration results and all the target approximate Hessian arrays, a set of linear inversion equations is constructed, and the set of linear inversion equations is solved to obtain the seismic imaging results.

[0008] Optionally, in the seismic imaging method based on the approximate Hessian array, the dip field includes a first dip field and a second dip field; The process of acquiring seismic data and migration parameters, and performing pre-stack time migration processing on the seismic data based on the migration parameters to obtain conventional migration results and dip fields, specifically includes: The raw seismic data is acquired and preprocessed to obtain seismic data. The preprocessing includes static correction, deconvolution, and amplitude compensation. The migration parameters are obtained, including the time sampling interval and number of time sampling points of the original seismic data, the spatial sampling interval, number of sampling points and spatial starting position of the target imaging area, the root mean square velocity field of the migration corresponding to the imaging space, and the sampling interval, number of samples, starting migration distance and migration aperture parameters of the pre-stack migration distance domain. Based on the principle of diffraction stacking, the seismic data is processed by pre-stack time migration according to the migration parameters to obtain conventional migration results and dip angle gathers. The reflection Fresnel zone range corresponding to each point in the imaging space is extracted based on the dip gather, and the actual stratum dip angle corresponding to each point in the imaging space is spatially mapped based on the reflection Fresnel zone range to obtain a first dip field along the lateral line and a second dip field along the direction perpendicular to the lateral line.

[0009] Optionally, in the seismic imaging method based on the approximate Hessian array, the target parameters include a first target parameter and a second target parameter; The process of determining the target imaging region involves dividing the target imaging region into multiple one-dimensional subspaces, obtaining the target parameters and cutoff frequency within each one-dimensional subspace, and generating a frequency integral result based on the target parameters, the cutoff frequency, and the tilt field. Specifically, this includes: Determine the target imaging region and divide the target imaging region into multiple one-dimensional subspaces extending along the depth direction according to a preset horizontal coordinate. Obtain the first target parameters in each of the one-dimensional subspaces, wherein the first target parameters include the gun point coordinates, the detector coordinates, the imaging point coordinates, the diffraction point coordinates, and the root mean square velocity field; A preset travel time formula is determined, and the first target parameter, the first dip angle field, and the second dip angle field are input into the preset travel time formula to obtain the second target parameter; The second target parameters include the travel time from the shot point to the imaging point, the travel time from the shot point to the diffraction point, the travel time from the detector point to the imaging point, the travel time from the detector point to the diffraction point, and the shortest travel time from the shot point through the reflection plane to the detector. The seismic spectrum is acquired, and the cutoff frequency corresponding to the seismic spectrum is determined. Numerical integration is performed based on the cutoff frequency and the second target parameter to obtain the frequency integration result.

[0010] Optionally, in the seismic imaging method based on the approximate Hessian array, the calculation expression for the second target parameter is: ; ; ; ; ; in, The travel time from the shot point to the imaging point. The travel time from the firing point to the diversion point. The travel time from the receiver point to the imaging point, The travel time from the receiver point to the diffraction point, The travel time of the reflected light from the shot point, through the reflecting plane, and towards the detector. The root mean square velocity field, Let the coordinate vector of the imaging point be... Let be the coordinate vector of the diffraction point. For the first A diffraction point, The x-coordinate of the shot point The x-coordinate of the receiver point The x-coordinate of the imaging point. Let x be the x-coordinate of the diffraction point. The vertical coordinate of the shot point. The vertical coordinate of the receiver point. The vertical coordinate of the imaging point. Let be the vertical coordinate of the diffraction point. and Let be the ordinate of the coordinate vector of the shot-receiver pair. The vertical coordinate of the imaging point. The ordinate of the diffraction point is . , as well as These are intermediate variables used in the derivation.

[0011] Optionally, in the seismic imaging method based on the approximate Hessian array, the expression for calculating the frequency integral result is: ; in, For the frequency integral result, The cutoff frequency, Angular frequency, It is the imaginary unit.

[0012] Optionally, in the seismic imaging method based on the approximate Hessian array, the calculation expression for the target approximate Hessian array is: ; in, For the target approximation Hessian array, The aperture is [not specified].

[0013] Optionally, the seismic imaging method based on the approximate Hessian array, wherein the step of constructing a linear inversion equation system based on the conventional migration results and all the target approximate Hessian arrays, and solving the linear inversion equation system to obtain the seismic imaging results, specifically includes: A set of linear inversion equations is constructed based on the target approximation Hessian matrix corresponding to each one-dimensional subspace and the conventional migration results, and a sparse regularization strategy is used to construct an objective function based on the set of linear inversion equations. The objective function is solved using an iterative shrinkage threshold algorithm to obtain the least squares offset result for each one-dimensional subspace. The least squares migration results corresponding to all the one-dimensional subspaces are combined according to the preset horizontal coordinates to obtain the least squares migration results of the whole space, and the least squares migration results are used as the seismic imaging results.

[0014] Furthermore, to achieve the above objectives, the present invention also provides a seismic imaging system based on an approximate Hessian array, wherein the seismic imaging system based on an approximate Hessian array comprises: The time migration processing module is used to acquire seismic data and migration parameters, and to perform pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and dip field. The frequency integration result generation module is used to determine the target imaging region, divide the target imaging region into multiple one-dimensional subspaces, obtain the target parameters and cutoff frequency in each one-dimensional subspace, and generate frequency integration results based on the target parameters, the cutoff frequency and the tilt field. An approximate Hessian matrix generation module is used to determine a preset approximate Hessian matrix analytical expression, and input the frequency integration result, the cutoff frequency and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace; The seismic imaging results output module is used to construct a set of linear inversion equations based on the conventional migration results and all the target approximate Hessian arrays, and to solve the set of linear inversion equations to obtain the seismic imaging results.

[0015] Furthermore, to achieve the above objectives, the present invention also provides a terminal, wherein the terminal includes: a memory, a processor, and a seismic imaging program based on an approximate Hessian array stored in the memory and executable on the processor, wherein when the seismic imaging program based on an approximate Hessian array is executed by the processor, it implements the steps of the seismic imaging method based on an approximate Hessian array as described above.

[0016] Furthermore, to achieve the above objectives, the present invention also provides a computer-readable storage medium, wherein the computer-readable storage medium stores a seismic imaging program based on an approximate Hessian array, which, when executed by a processor, implements the steps of the seismic imaging method based on an approximate Hessian array as described above.

[0017] In this invention, seismic data and migration parameters are acquired, and pre-stack time migration processing is performed on the seismic data according to the migration parameters to obtain conventional migration results and dip fields. A target imaging region is determined, and the target imaging region is divided into multiple one-dimensional subspaces. Target parameters and cutoff frequencies are acquired within each one-dimensional subspace, and frequency integration results are generated based on the target parameters, cutoff frequencies, and dip fields. A preset approximate Hessian array analytical expression is determined, and the frequency integration results, cutoff frequencies, and target parameters are input into the preset approximate Hessian array analytical expression to obtain the target approximate Hessian array corresponding to each one-dimensional subspace. A linear inversion equation system is constructed based on the conventional migration results and all the target approximate Hessian arrays, and the linear inversion equation system is solved to obtain the seismic imaging results. This invention constructs an approximate Hessian array by superimposing Hessian array elements within the same reflecting surface, and then constructs and solves the linear inversion equation system to obtain the seismic imaging results, significantly improving the efficiency and accuracy of seismic imaging. Attached Figure Description

[0018] Figure 1 This is a flowchart of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention; Figure 2 This is a schematic diagram illustrating the specific implementation process of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention; Figure 3 This is a schematic diagram of an approximate Hessian array, representing a preferred embodiment of the seismic imaging method based on an approximate Hessian array according to the present invention. Figure 4 This is a schematic diagram comparing the theoretical model migration profile of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention. Figure 5 This is a schematic diagram of the theoretical model of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention, showing the amplitude comparison along the phase axis. Figure 6 This is a schematic diagram comparing the amplitude spectra of the theoretical model of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention; Figure 7 This is a schematic diagram of the actual model dip gather and dip field of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention; Figure 8 This is a schematic diagram comparing the actual model migration profiles of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention. Figure 9 This is a schematic diagram comparing actual model CRP gathers of a preferred embodiment of the seismic imaging method based on the approximate Hessian array of the present invention; Figure 10 This is a structural diagram of a preferred embodiment of the seismic imaging system based on an approximate Hessian array according to the present invention; Figure 11 This is a structural diagram of a preferred embodiment of the terminal of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer and more explicit, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

[0020] In seismic imaging processing, the conventional migration operator is used as an adjoint operator of the forward modeling operator rather than an inverse operator. Under non-ideal conditions such as incomplete observation systems, uneven sampling, or limited coverage, it is prone to a series of imaging defects such as amplitude distortion, migration arcing, and uneven illumination. To overcome these problems, the imaging domain least-squares migration method was proposed. Its core lies in explicitly constructing and storing the Hessian array, then compensating for imaging errors through an inversion process, recovering true amplitude information, and effectively suppressing migration noise. Current calculations of the Hessian array mainly rely on two types of methods: one is an indirect calculation method based on the point spread function, which requires constructing the point spread function point by point through migration and inverse migration; the other is a direct calculation method based on high-frequency asymptotic approximation. This method has certain advantages in computational efficiency, but it does not approximate the size of the Hessian array, resulting in relatively low efficiency.

[0021] Nevertheless, due to the quadratic relationship between the dimension of the Hessian array and the size of the imaging domain, its storage requirements in practical 3D applications are extremely large, becoming a key bottleneck restricting the industrial-scale promotion of this method. Although the storage capacity of computing hardware has been continuously improving in recent years, traditional imaging domain least squares migration still cannot be effectively applied within acceptable storage and computing costs when dealing with large-scale 3D seismic data from actual work areas. Therefore, existing imaging domain least squares migration techniques have not yet been able to fully realize their theoretical advantages in actual production, making it difficult to stably obtain high-fidelity, high-resolution true amplitude imaging results.

[0022] In summary, there is still a significant gap between the theoretical feasibility and the engineering feasibility of the current technology. Systematic improvements are urgently needed in areas such as algorithm design, storage optimization, and computing architecture to promote the transition of imaging domain least squares offset from theoretical methods to practical industrial applications.

[0023] To address the issue of excessive storage requirements for Hessian arrays in traditional least-squares migration in the imaging domain, this invention proposes an improved method based on the concept of localization. This method first constructs an approximate Hessian array by superimposing Hessian array elements within the same reflective surface. Then, based on the principle of phase stability, its analytical expression is derived and integrated into the least-squares migration process in the imaging domain. Since the size of the approximate Hessian array is much smaller than the original Hessian array, and the analytical expression avoids time-consuming explicit superposition operations, this invention significantly reduces storage requirements and computational complexity. Experimental results show that this method effectively improves amplitude distortion and migration arcing caused by defects in the observation system while significantly enhancing the computational efficiency of the imaging process.

[0024] This invention relates to the field of seismic exploration. In view of the above-mentioned shortcomings of traditional imaging domain least squares migration, it provides a seismic imaging method based on an approximate Hessian array, aiming to solve the problem that existing imaging domain least squares migration imaging cannot efficiently obtain imaging results.

[0025] The seismic imaging method based on the approximate Hessian array described in the preferred embodiment of the present invention, such as... Figure 1 and Figure 2 As shown, the seismic imaging method based on the approximate Hessian array includes the following steps: Step S10: Acquire seismic data and migration parameters, and perform pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and dip fields. The dip fields include a first dip field and a second dip field.

[0026] This invention comprises a migration processing module, an approximate Hessian array generation module, and an inversion module. The migration processing module performs conventional pre-stack time migration on the original seismic data according to migration processing parameters, generating migration results and dip gathers. The migration processing parameters include the temporal sampling interval and duration of the original seismic data, the spatial sampling interval and size of the imaging target region, the root mean square velocity field used for migration, and related pre-stack domain parameters. The pre-stack domain parameters mainly include the migration distance sampling interval and the migration distance range. Since the influence of the Hessian array is not considered, the migration results may contain some migration arcing and amplitude distortion. The dip gathers are used to estimate the range of the Fresnel reflection zone corresponding to each point in the imaging space, and then derive the mapping relationship of the true strata dip angle at each point in the time domain.

[0027] Specifically, raw seismic data is acquired and preprocessed to obtain seismic data. The preprocessing includes static correction, deconvolution, and amplitude compensation. Migration parameters are acquired, including the time sampling interval and number of time sampling points of the raw seismic data, the spatial sampling interval, number of sampling points, and spatial starting position of the target imaging area, the root mean square velocity field of the migration corresponding to the imaging space, and the sampling interval, number of samples, starting migration distance, and migration aperture parameters of the pre-stack migration distance domain. Based on the diffraction stacking principle, pre-stack time migration processing is performed on the seismic data according to the migration parameters to obtain conventional migration results and dip gathers. The reflection Fresnel zone range corresponding to each point in the imaging space is extracted from the dip gathers, and the spatial domain mapping processing is performed on the actual stratum dip angle corresponding to each point in the imaging space based on the reflection Fresnel zone range to obtain a first dip field along the lateral direction and a second dip field perpendicular to the lateral direction.

[0028] like Figure 2 As shown, based on seismic data and migration parameters obtained from seismic exploration, this invention generates conventional migration results and dip gathers through pre-stack time migration based on the principle of diffraction stacking. Seismic data refers to digital signals generated through artificial seismic exploration, recorded by ground geophones, and preprocessed. The preprocessing includes static correction, deconvolution, and amplitude compensation. The migration parameters set in this invention include: the time sampling interval dt and the number of time sampling points nt of the original seismic data; the spatial sampling intervals dx and dy, the number of sampling points nx and ny, and the spatial starting positions xbeg and ybeg of the target imaging area; the root mean square velocity field Vrms(nx, ny, nt) corresponding to the imaging space; and the sampling interval doff, the number of samples noff, the starting migration distance offbeg, and the migration aperture parameters of the pre-stack migration distance domain.

[0029] This invention utilizes the aforementioned data and parameters for pre-stack time migration to obtain conventional migration results mig(nx, ny, nt) and dip gathers. Since Hessian array correction is not introduced, conventional migration results may contain migration arcs and amplitude distortion. Dip gathers can be obtained by rearranging the migration results according to specific parameters. They can be used to extract the Fresnel reflection zone range corresponding to each point in the imaging space through manual picking, algorithmic picking, or artificial intelligence methods, and then estimate the mapping of the true stratigraphic dip angle at each point in space to the spatial domain, obtaining the dip field dipx(nx, ny, nt) along the lateral direction and the dip field dipy(nx, ny, nt) perpendicular to the lateral direction.

[0030] Step S20: Determine the target imaging region, divide the target imaging region into multiple one-dimensional subspaces, obtain the target parameters and cutoff frequency in each one-dimensional subspace, and generate a frequency integral result based on the target parameters, the cutoff frequency, and the tilt field. The target parameters include a first target parameter and a second target parameter.

[0031] Specifically, the target imaging region is determined and divided into multiple one-dimensional subspaces extending along the depth direction according to a preset horizontal coordinate; a first target parameter is obtained in each of the one-dimensional subspaces, wherein the first target parameter includes the shot point coordinates, the detector coordinates, the imaging point coordinates, the diffraction point coordinates, and the root mean square velocity field; a preset travel time formula is determined, and the first target parameter, the first tilt field, and the second tilt field are input into the preset travel time formula to obtain the second target parameter; wherein the second target parameter includes the travel time from the shot point to the imaging point, the travel time from the shot point to the diffraction point, the travel time from the detector point to the imaging point, the travel time from the detector point to the diffraction point, and the shortest travel time from the shot point through the reflection plane to the detector.

[0032] The expression for calculating the second target parameter is: ; ; ; ; ; in, The travel time from the shot point to the imaging point. The travel time from the firing point to the diversion point. The travel time from the receiver point to the imaging point, The travel time from the receiver point to the diffraction point, The travel time of the reflected light from the shot point, through the reflecting plane, and towards the detector. The root mean square velocity field, Let the coordinate vector of the imaging point be... Let be the coordinate vector of the diffraction point. For the first A diffraction point, The x-coordinate of the shot point The x-coordinate of the receiver point The x-coordinate of the imaging point. Let x be the x-coordinate of the diffraction point. The vertical coordinate of the shot point. The vertical coordinate of the receiver point. The vertical coordinate of the imaging point. Let be the vertical coordinate of the diffraction point. and Let be the ordinate of the coordinate vector of the shot-receiver pair. The vertical coordinate of the imaging point. The ordinate of the diffraction point is . , as well as These are intermediate variables used in the derivation.

[0033] To implement a localization strategy, this invention first divides the imaging space (i.e., the target imaging area in this invention) into several one-dimensional subspaces extending along the depth direction according to different horizontal coordinates (nx, ny) (i.e., the preset horizontal coordinates in this invention). Within each one-dimensional subspace, based on the shot coordinates (Xs, Ys, 0), detector coordinates (Xg, Yg, 0), imaging point coordinates (Xp, Yp, Zp), diffraction point coordinates (Xq, Yq, Zq), root mean square velocity field Vrms (nx, ny, nt), tilt field dipx (nx, ny, nt) along the lateral direction, and tilt field dipy (nx, ny, nt) perpendicular to the lateral direction, the travel times tsp and tsq from the shot to the imaging and diffraction points, the travel times tgp and tgq from the detector to the imaging and diffraction points, and the shortest travel time tsg from the shot to the detector via the reflection plane are calculated using a preset travel time formula. The calculation expression is as follows: ; ; ; ; .

[0034] The seismic spectrum is acquired, and the cutoff frequency corresponding to the seismic spectrum is determined. Numerical integration is performed based on the cutoff frequency and the second target parameter to obtain the frequency integration result.

[0035] The expression for calculating the frequency integral result is as follows: ; in, For the frequency integral result, The cutoff frequency, Angular frequency, It is the imaginary unit.

[0036] It is understood that after obtaining the aforementioned travel time (i.e., the second target parameter in this invention), the frequency integral result is obtained by numerical integration based on the cutoff frequency fc determined from the seismic spectrum and the aforementioned travel time. The expression for calculating the frequency integral result is as follows: .

[0037] Step S30: Determine the preset approximate Hessian matrix analytical expression, and input the frequency integral result, the cutoff frequency, and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace.

[0038] This invention also includes an approximate Hessian array generation module, which calculates the approximate Hessian array in each subspace based on the shot-receiver pair coordinates of each seismic data, the cutoff frequency of the seismic data, the aperture range corresponding to each imaging point, and the dip field obtained by the migration processing module, for subsequent inversion use. The shot-receiver pair coordinates are mainly used to calculate the travel times from the shot point and receiver to the imaging point and diffraction point (i.e., the travel times from the shot point to the imaging point, from the shot point to the diffraction point, from the receiver to the imaging point, from the receiver to the diffraction point, and the shortest travel time from the shot point through the reflection plane to the receiver) and the first half of the amplitude term of the approximate Hessian array; the cutoff frequency of the seismic data... The aperture range used to determine the frequency integration range must be consistent with the aperture range used in the offset process.

[0039] Specifically, the calculation expression for the target approximate Hessian matrix is ​​as follows: ; in, For the target approximation Hessian array, The aperture is [not specified].

[0040] Substituting the aforementioned intermediate parameters into the above approximate Hessian matrix analytical expression (i.e., inputting the frequency integral result, the cutoff frequency, and the target parameter into the preset approximate Hessian matrix analytical expression), the value of the approximate Hessian matrix is ​​calculated. The specific calculation expression is as follows: .

[0041] It is understood that this invention provides an approximate analytical expression for the Hessian matrix based on the principle of stable phase, avoiding the numerical superposition of the Hessian matrix and significantly reducing the computational cost of the aforementioned localization method. The specific process is as follows: Starting from the three-dimensional Hessian array based on deconvolution imaging conditions within the pre-stack time migration framework, namely: ; in, and For the coordinate vector of the shot-receiver pair, sum the signs. This indicates the aperture corresponding to each imaging point. Internal artillery inspection and summation, Let the travel times from the shot point and detector to the imaging point and diffraction point be respectively. The expressions for the above travel times are: ; ; ; ; ; in, The travel time from the shot point to the imaging point. The travel time from the firing point to the diversion point. The travel time from the receiver point to the imaging point, The travel time from the receiver point to the diffraction point, The travel time of the reflected light from the shot point, through the reflecting plane, and towards the detector. The root mean square velocity field, Let the coordinate vector of the imaging point be... Let be the coordinate vector of the diffraction point. For the first A diffraction point, The x-coordinate of the shot point The x-coordinate of the receiver point The x-coordinate of the imaging point. Let x be the x-coordinate of the diffraction point. The vertical coordinate of the shot point. The vertical coordinate of the receiver point. The vertical coordinate of the imaging point. Let be the vertical coordinate of the diffraction point. and Let be the ordinate of the coordinate vector of the shot-receiver pair. The vertical coordinate of the imaging point. The ordinate of the diffraction point is . , as well as These are intermediate variables used in the derivation.

[0042] Substituting the above parameters into the approximate Hessian array expression, considering the spectral characteristics of seismic data, this expression is actually an oscillatory integral, which can be approximated based on the steady-phase principle, thus obtaining the corrected approximate Hessian array analytical expression, namely: ; in, , , and These are intermediate variables used in the derivation. The travel time of the reflected light from the shot point, through the reflecting plane, and towards the detector. and The tilt angle of the reflecting plane along the lateral line and perpendicular to the lateral line (e.g.) Figure 3 As shown), the expression for the above parameters is as follows: ; ; ; ; .

[0043] Step S40: Construct a linear inversion equation set based on the conventional migration results and all the target approximate Hessian arrays, and solve the linear inversion equation set to obtain the seismic imaging results.

[0044] This invention also includes an inversion module, used to construct a system of linear equations within each subspace based on the regular migration result vector generated by the migration processing module and the approximate Hessian matrix provided by the approximate Hessian matrix generation module, and to obtain the least-squares migration result by solving this system of equations. Considering that in practice, the Hessian matrix and its approximate form often have large condition numbers, which can easily lead to unstable or unsolvable solutions to the linear equations, the inversion module introduces a hybrid sparse regularization method based on TV and L1 norms in the inversion process to enhance the numerical stability of the inversion process.

[0045] Specifically, a set of linear inversion equations is constructed based on the target approximation Hessian matrix corresponding to each one-dimensional subspace and the conventional migration results. A sparse regularization strategy is then used to construct an objective function based on the set of linear inversion equations. An iterative shrinkage threshold algorithm is used to solve the objective function to obtain the least squares migration results corresponding to each one-dimensional subspace. The least squares migration results corresponding to all one-dimensional subspaces are combined according to the preset horizontal coordinates to obtain the least squares migration results for the entire space. The least squares migration results are then used as the seismic imaging results.

[0046] It is understood that this invention constructs a linear equation system h0*r0=mig0 by combining the conventional offset result mig0(nt) in the current subspace with the approximate Hessian matrix h0(nt, nt), where r0(nt) represents the least-squares offset result of the subspace, specifically: ; Furthermore, this invention introduces a sparse regularization strategy to construct an objective function, and uses an iterative shrinking threshold algorithm to solve the equations to obtain the seismic imaging result r0(nt). Finally, the least squares migration results in all subspaces are combined according to different horizontal coordinates (nx, ny) to obtain the full-space least squares migration imaging result r(nx, ny, nt).

[0047] It is understood that this invention provides a Hessian matrix localization method based on diagonal superposition, thereby obtaining a modified approximate Hessian matrix expression. The specific process is as follows: 1. Based on the least squares migration theory of the imaging domain, the migration result... It can be viewed as an underground reflection coefficient matrix After Heisenberg array The filtered result is: ; in, Let the coordinate vector of the imaging point be... Let be the coordinate vector of the diffraction point. For the entire imaging space.

[0048] 2. Considering that the size of the entire Hessian array is the square of the total imaging space, it generally cannot be explicitly calculated and used for inversion. To reduce the size of the imaging space, this invention only considers several diffraction points on the reflection plane and makes the assumption that the reflection coefficient is the same along the reflection plane. Then, the migration result can be approximated as: ; in, For the diffraction point The reflecting plane, To consider the number of reflecting planes, we can treat the integral term on the right-hand side of the above equation as an approximate Hessian matrix, then we have: ; in, for The diffraction point on the reflecting plane relative to the imaging point The influence of this, namely the approximate Hessian array proposed in this invention.

[0049] 3. For a scale of The three-dimensional problem (where, (These represent the number of sampling points in the three directions of the imaging space, respectively). To implement the localization strategy, a system of linear equations is established by taking the imaging points and diffraction points with the same horizontal coordinates and solving the equations (total number of...). The influence of diffraction points with different horizontal coordinates on the imaging point is transformed into a corrected approximate Hessian array through superposition. This successfully reduces the size of the Hessian array from (a few points) to (a few points). Down to Furthermore, by repeating this process for different horizontal coordinates, the least squares migration result for the entire space can be obtained, and the least squares migration result can be used as the seismic imaging result.

[0050] The imaging effect of this invention on the theoretical model is as follows: Figure 4 As shown, compared with the conventional offset results (i.e. Figure 4 Compared to (a) in the previous example, the least squares offset result (i.e. Figure 4 (b) effectively eliminates the offset arc and amplitude distortion phenomena. Furthermore, the normalized amplitude along the phase axis (such as...) Figure 5 The normalized amplitude spectrum (as shown) is closer to the true value than the actual value. Figure 6 As shown in the figure, significant compensation was also obtained, verifying the correctness and rationality of the proposed method. The application effect of this method on a practical model is as follows: Figure 7 , Figure 8 and Figure 9 As shown. Least squares offset results (i.e. Figure 8 and Figure 9 (b) in the middle is different from the conventional offset result (i.e. Figure 8 and Figure 9 (a) is closer to the migration result when the observation system is complete (i.e. Figure 8 and Figure 9 (c) indicates that the present invention can effectively compensate for the problem of insufficient lighting caused by the incomplete observation system.

[0051] Therefore, this invention calculates the approximate Hessian arrays corresponding to the theoretical and actual models respectively, and constructs a system of linear equations based on these matrices and conventional migration results to solve the inverse problem, ultimately obtaining the least-squares migration results in the imaging domain. Because the constructed approximate Hessian array fully considers the influence of the seismic wave acquisition system and propagation process, it effectively eliminates amplitude distortion and migration arcing, improving the amplitude fidelity of the migration data and thus enhancing the accuracy of lithology inversion based on this data.

[0052] Furthermore, this invention provides a GPU parallelization process and comprehensive acceleration scheme based on the above-mentioned system. This process covers the entire process from data preprocessing to result aggregation, specifically including: first, data preprocessing, transferring relevant data such as offset results, shot-receiver coordinates, and root mean square velocity field from the host memory to the GPU; then, parallel execution of the approximate Hessian matrix calculation and inversion process on the GPU, generating least-squares offset results for each sub-block; finally, transferring the results of each sub-block back from the GPU to the host, where they are combined into a complete three-dimensional imaging domain least-squares offset result. Regarding acceleration optimization, this invention achieves efficient parallel processing of multiple sub-blocks through batch computation optimization, maximizing the utilization of GPU computing resources; with operator fusion optimization, inversion is performed directly after the approximate Hessian matrix calculation is completed on the GPU, avoiding redundant data transmission between the host and the device; simultaneously, memory access optimization is implemented, redesigning the mapping relationship between threads and matrix elements to optimize memory access patterns, significantly improving data read / write efficiency, thereby systematically improving the overall computational performance.

[0053] Beneficial effects: This invention provides a Hessian array localization method based on superimposed reflective surfaces, and derives a modified approximate analytical expression for the Hessian array using the phase-stable principle, achieving effective dimensionality reduction and efficient computation for large-scale three-dimensional least-squares migration problems. By decomposing the global Hessian array into a series of approximate Hessian arrays, this invention significantly reduces the dimensionality of the Hessian array, greatly reducing the storage space requirements and making three-dimensional least-squares migration based on real-world large-scale data feasible on conventional computing platforms. Furthermore, the approximate analytical expression proposed in this invention is constructed based on the phase-stable principle, avoiding the explicit superposition operations required in the aforementioned localization strategies, thereby significantly reducing computational complexity. Compared with existing imaging domain least-squares migration techniques, this invention has significant advantages in both storage efficiency (space complexity) and computational speed (time complexity), while maintaining good imaging accuracy control.

[0054] To fully realize the engineering application potential of this method, this invention also includes a targeted GPU parallel computing architecture and multi-level optimization strategies. Through batch computation optimization, operator fusion optimization, and memory access optimization, the system improves the overall computational efficiency. This invention provides stable, reliable, and highly scalable technical support for realizing least-squares migration imaging of large-scale 3D seismic data in the imaging domain.

[0055] Furthermore, such as Figure 10 As shown, based on the above-described seismic imaging method based on the approximate Hessian array, the present invention also provides a seismic imaging system based on the approximate Hessian array, wherein the seismic imaging system based on the approximate Hessian array includes: The time migration processing module 51 is used to acquire seismic data and migration parameters, and to perform pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and dip field. The frequency integration result generation module 52 is used to determine the target imaging region, divide the target imaging region to obtain multiple one-dimensional subspaces, obtain the target parameters and cutoff frequency in each one-dimensional subspace, and generate frequency integration results based on the target parameters, the cutoff frequency and the tilt field. The approximate Hessian matrix generation module 53 is used to determine a preset approximate Hessian matrix analytical expression, and input the frequency integral result, the cutoff frequency and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace; The seismic imaging result output module 54 is used to construct a linear inversion equation set based on the conventional migration results and all the target approximate Hessian arrays, and to solve the linear inversion equation set to obtain the seismic imaging results.

[0056] Furthermore, such as Figure 11 As shown, based on the above-mentioned seismic imaging method and system based on the approximate Hessian array, the present invention also provides a terminal, which includes a processor 10, a memory 20 and a display 30. Figure 11 Only some of the terminal components are shown; however, it should be understood that it is not required to implement all of the components shown, and more or fewer components may be implemented instead.

[0057] In some embodiments, the memory 20 may be an internal storage unit of the terminal, such as a hard disk or memory. In other embodiments, the memory 20 may be an external storage device of the terminal, such as a plug-in hard disk, smart media card (SMC), secure digital card (SD), flash card, etc. Further, the memory 20 may include both internal and external storage devices. The memory 20 is used to store application software and various types of data installed on the terminal, such as the program code installed on the terminal. The memory 20 can also be used to temporarily store data that has been output or will be output. In one embodiment, the memory 20 stores a seismic imaging program 40 based on an approximate Hessian array, which can be executed by the processor 10 to implement the seismic imaging method based on an approximate Hessian array as described in this application.

[0058] In some embodiments, the processor 10 may be a central processing unit (CPU), a microprocessor, or other data processing chip, used to run program code stored in the memory 20 or process data, such as executing the seismic imaging method based on the approximate Hessian array.

[0059] In some embodiments, the display 30 may be an LED display, a liquid crystal display, a touch-sensitive liquid crystal display, or an OLED (Organic Light-Emitting Diode) touchscreen. The display 30 is used to display information on the terminal and to display a visual user interface.

[0060] In one embodiment, when the processor 10 executes the seismic imaging program 40 based on the approximate Hessian array in the memory 20, it implements the steps of the seismic imaging method based on the approximate Hessian array as described above.

[0061] The present invention also provides a computer-readable storage medium, wherein the computer-readable storage medium stores a seismic imaging program based on an approximate Hessian array, which, when executed by a processor, implements the steps of the seismic imaging method based on an approximate Hessian array as described above.

[0062] In summary, this invention provides a seismic imaging method, system, terminal, and storage medium based on an approximate Hessian array. The method includes: acquiring seismic data and migration parameters, and performing pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and a dip field; determining a target imaging region, dividing the target imaging region into multiple one-dimensional subspaces, acquiring target parameters and cutoff frequencies in each one-dimensional subspace, and generating frequency integration results based on the target parameters, the cutoff frequencies, and the dip field; determining a preset approximate Hessian array analytical expression, and inputting the frequency integration results, the cutoff frequencies, and the target parameters into the preset approximate Hessian array analytical expression to obtain the target approximate Hessian array corresponding to each one-dimensional subspace; constructing a linear inversion equation system based on the conventional migration results and all the target approximate Hessian arrays, and solving the linear inversion equation system to obtain the seismic imaging result. This invention constructs an approximate Hessian array by superimposing Hessian array elements within the same reflecting surface, and then constructs and solves the linear inversion equation system to obtain the seismic imaging result, significantly improving the efficiency and accuracy of seismic imaging.

[0063] It should be noted that, in this document, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or terminal that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or terminal. Unless otherwise specified, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or terminal that includes that element.

[0064] Of course, those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware (such as a processor, controller, etc.). The program can be stored in a computer-readable storage medium, and when executed, it can include the processes described in the above method embodiments. The computer-readable storage medium can be a memory, magnetic disk, optical disk, etc.

[0065] It should be understood that the application of the present invention is not limited to the examples above. Those skilled in the art can make improvements or modifications based on the above description, and all such improvements and modifications should fall within the protection scope of the appended claims.

Claims

1. A seismic imaging method based on an approximate Hessian array, characterized in that, The seismic imaging method based on the approximate Hessian array includes: Seismic data and migration parameters are acquired, and the seismic data are processed by pre-stack time migration according to the migration parameters to obtain conventional migration results and dip field. The target imaging region is determined, and the target imaging region is divided into multiple one-dimensional subspaces. The target parameters and cutoff frequencies in each one-dimensional subspace are obtained, and frequency integral results are generated based on the target parameters, the cutoff frequencies, and the tilt field. Determine a preset approximate Hessian matrix analytical expression, and input the frequency integral result, the cutoff frequency and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace; Based on the conventional migration results and all the target approximate Hessian arrays, a set of linear inversion equations is constructed, and the set of linear inversion equations is solved to obtain the seismic imaging results.

2. The seismic imaging method based on an approximate Hessian array according to claim 1, characterized in that, The tilt field includes a first tilt field and a second tilt field; The process of acquiring seismic data and migration parameters, and performing pre-stack time migration processing on the seismic data based on the migration parameters to obtain conventional migration results and dip fields, specifically includes: The raw seismic data is acquired and preprocessed to obtain seismic data. The preprocessing includes static correction, deconvolution, and amplitude compensation. The migration parameters are obtained, including the time sampling interval and number of time sampling points of the original seismic data, the spatial sampling interval, number of sampling points and spatial starting position of the target imaging area, the root mean square velocity field of the migration corresponding to the imaging space, and the sampling interval, number of samples, starting migration distance and migration aperture parameters of the pre-stack migration distance domain. Based on the principle of diffraction stacking, the seismic data is processed by pre-stack time migration according to the migration parameters to obtain conventional migration results and dip angle gathers. The reflection Fresnel zone range corresponding to each point in the imaging space is extracted based on the dip gather, and the actual stratum dip angle corresponding to each point in the imaging space is spatially mapped based on the reflection Fresnel zone range to obtain a first dip field along the lateral line and a second dip field along the direction perpendicular to the lateral line.

3. The seismic imaging method based on an approximate Hessian array according to claim 2, characterized in that, The target parameters include a first target parameter and a second target parameter; The process of determining the target imaging region involves dividing the target imaging region into multiple one-dimensional subspaces, obtaining the target parameters and cutoff frequency within each one-dimensional subspace, and generating a frequency integral result based on the target parameters, the cutoff frequency, and the tilt field. Specifically, this includes: Determine the target imaging region and divide the target imaging region into multiple one-dimensional subspaces extending along the depth direction according to a preset horizontal coordinate. Obtain the first target parameters in each of the one-dimensional subspaces, wherein the first target parameters include the gun point coordinates, the detector coordinates, the imaging point coordinates, the diffraction point coordinates, and the root mean square velocity field; A preset travel time formula is determined, and the first target parameter, the first dip angle field, and the second dip angle field are input into the preset travel time formula to obtain the second target parameter; The second target parameters include the travel time from the shot point to the imaging point, the travel time from the shot point to the diffraction point, the travel time from the detector point to the imaging point, the travel time from the detector point to the diffraction point, and the shortest travel time from the shot point through the reflection plane to the detector. The seismic spectrum is acquired, and the cutoff frequency corresponding to the seismic spectrum is determined. Numerical integration is performed based on the cutoff frequency and the second target parameter to obtain the frequency integration result.

4. The seismic imaging method based on an approximate Hessian array according to claim 3, characterized in that, The expression for calculating the second target parameter is: ； ； ； ； ； in, The travel time from the shot point to the imaging point. The travel time from the firing point to the diversion point. The travel time from the receiver point to the imaging point, The travel time from the receiver point to the diffraction point, The travel time of the reflected light from the shot point, through the reflecting plane, and towards the detector. The root mean square velocity field, Let the coordinate vector of the imaging point be... Let be the coordinate vector of the diffraction point. For the first A diffraction point, The x-coordinate of the shot point The x-coordinate of the receiver point The x-coordinate of the imaging point. Let x be the x-coordinate of the diffraction point. The vertical coordinate of the shot point. The vertical coordinate of the receiver point. The vertical coordinate of the imaging point. Let be the vertical coordinate of the diffraction point. and Let be the ordinate of the coordinate vector of the shot-receiver pair. The vertical coordinate of the imaging point. The ordinate of the diffraction point is . , as well as These are intermediate variables used in the derivation.

5. The seismic imaging method based on an approximate Hessian array according to claim 4, characterized in that, The expression for calculating the frequency integral result is as follows: ； in, For the frequency integral result, The cutoff frequency, Angular frequency, It is the imaginary unit.

6. The seismic imaging method based on an approximate Hessian array according to claim 5, characterized in that, The calculation expression for the target approximate Hessian matrix is ​​as follows: ； in, For the target approximation Hessian array, The aperture is [not specified].

7. The seismic imaging method based on an approximate Hessian array according to claim 3, characterized in that, The process of constructing a linear inversion equation set based on the conventional migration results and all the target approximate Hessian arrays, and solving the linear inversion equation set to obtain the seismic imaging results, specifically includes: A set of linear inversion equations is constructed based on the target approximation Hessian matrix corresponding to each one-dimensional subspace and the conventional migration results, and a sparse regularization strategy is used to construct an objective function based on the set of linear inversion equations. The objective function is solved using an iterative shrinkage threshold algorithm to obtain the least squares offset result for each one-dimensional subspace. The least squares migration results corresponding to all the one-dimensional subspaces are combined according to the preset horizontal coordinates to obtain the least squares migration results of the entire space, and the least squares migration results are used as the seismic imaging results.

8. A seismic imaging system based on an approximate Hessian array, characterized in that, The seismic imaging system based on the approximate Hessian array includes: The time migration processing module is used to acquire seismic data and migration parameters, and to perform pre-stack time migration processing on the seismic data according to the migration parameters to obtain conventional migration results and dip field. The frequency integration result generation module is used to determine the target imaging region, divide the target imaging region into multiple one-dimensional subspaces, obtain the target parameters and cutoff frequency in each one-dimensional subspace, and generate frequency integration results based on the target parameters, the cutoff frequency and the tilt field. An approximate Hessian matrix generation module is used to determine a preset approximate Hessian matrix analytical expression, and input the frequency integral result, the cutoff frequency and the target parameter into the preset approximate Hessian matrix analytical expression to obtain the target approximate Hessian matrix corresponding to each one-dimensional subspace; The seismic imaging results output module is used to construct a set of linear inversion equations based on the conventional migration results and all the target approximate Hessian arrays, and to solve the set of linear inversion equations to obtain the seismic imaging results.

9. A terminal, characterized in that, The terminal includes: a memory, a processor, and a seismic imaging program based on an approximate Hessian array stored in the memory and executable on the processor. When executed by the processor, the seismic imaging program based on an approximate Hessian array implements the steps of the seismic imaging method based on an approximate Hessian array as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a seismic imaging program based on an approximate Hessian array, which, when executed by a processor, implements the steps of the seismic imaging method based on an approximate Hessian array as described in any one of claims 1-7.

Images