阅读说明：本技术 基于点扩散函数的Hessian矩阵的计算与存储方法 (Calculation and storage method of Hessian matrix based on point spread function ) 是由 倪文军 刘少勇 韦俊坤 于 2021-06-17 设计创作，主要内容包括：本发明提供基于点扩散函数的Hessian矩阵的计算与存储方法,包括基于高斯束进行射线追踪获得射线参数,通过射线参数计算旅行时,振幅,用点扩散函数近似替代Hessian矩阵,利用WKBJ近似导出Hessian矩阵解析公式,计算Hessian矩阵,利用Hessian矩阵对偏移结果进行校正。本发明引入短波逼近法(WKBJ)下的近似格林函数以及在局部窗中对Hessian矩阵进行构建,提高了计算效率,降低了运算的规模。(The invention provides a calculation and storage method of a Hessian matrix based on a point spread function, which comprises the steps of carrying out ray tracing based on a Gaussian beam to obtain ray parameters, calculating amplitude during travel through the ray parameters, approximately replacing the Hessian matrix with the point spread function, deriving a Hessian matrix analytic formula by utilizing WKBJ approximation, calculating the Hessian matrix, and correcting an offset result by utilizing the Hessian matrix. The invention introduces approximate Green's function under short wave approximation (WKBJ) and constructs Hessian matrix in local window, thus improving calculation efficiency and reducing operation scale.)

1. The calculation and storage method of the Hessian matrix based on the point spread function is characterized by comprising the following steps of:

s1, carrying out ray tracing and calculating the travel time of the underground space based on the Gaussian beam and according to the background speed;

s2, obtaining PSFs by adopting an asymptotic Green function under the short wave approximation method WKBJ approximation and approximate estimation in a local window;

s3, the offset image is corrected by using the obtained PSFs, and a corrected image is obtained, where the formula is as follows:

m＝H^-1I_mig

wherein, the offset imaging result is marked as I_migThe ideal image is denoted as m.

2. The method for calculating and storing the Hessian matrix based on the point spread function as claimed in claim 1, wherein the S1 is specifically as follows:

s11, applying for the construction structure array Cells (Nz, Nx, Npx) to store travel time in different directions, and applying for the travel time array to store single shot time (Nz, Nx, Nshot); defining Cells (Nz, Nx, Npx), time (Nz, Nx, Nshot) and related cycle-related variables including Nz, Nx, Npx, Nshot, Ix, Iz, Ipx;

s12, shot point circulation, namely Nshot circulation is carried out on shot points;

s13, shot direction circulation, namely, Npx times of circulation is carried out on shot direction shot, forward modeling is carried out on a Gaussian beam in one direction Ipx to obtain ray parameters, and Cells (Nz, Nx, Npx) in the direction are obtained;

s14, circulating model space for 1-Nx NZ times; accumulating traveling time in different directions for the grid points of Iz and Ix in the model space to obtain single-cannon traveling time (Iz, Ix and Ishot);

where z and x refer to the two-dimensional coordinates of a point in space, where I ═ 1, N.

3. The method for calculating and storing the Hessian matrix based on the point spread function as claimed in claim 1, wherein the S2 is specifically as follows:

s21, applying for a Hessian array H (Ni, Nj), wherein Ni and Nj are the sizes of a horizontal model and a vertical model of the Hessian matrix respectively;

s22, shot point circulation, namely performing Nshot circulation on the shot point Ishot;

s23, wave detection point circulation, namely processing any wave detector Ireceiver, and performing 1-Ntrace times of circulation on the wave detection points, wherein Nshot is the number of shot points;

s24, model space circulation, circulation is performed for 1-Ni Nj times, and model space integral output is as follows: h (i, j)

Wherein, Ni Nj is the size of model space of the Hessian matrix.

4. The method for calculating and storing the Hessian matrix based on the point spread function as claimed in claim 3, wherein the model space integration outputs: the specific formula of H (i, j) is as follows:

wherein R is_ssFor the autocorrelation function, T represents the travel time, xi represents the index of the shot-check pair, f (omega; s)_x) Representing the source wavelet, A_rayRepresentative amplitude can be obtained by travel time, ω represents frequency, x_rRepresenting the coordinates of the demodulator probe, x_sRepresenting the coordinates of the shot point, x_j、x_iRepresenting the subsurface space point coordinates.

Technical Field

The invention relates to the technical field of oil and gas exploration, in particular to a calculation and storage method of a Hessian matrix based on a point spread function.

Background

The Hessian matrix is a function describing the imaging resolution at a point in the model, which is usually a function of space and contains all the information that affects the imaging resolution. The Hessian matrix is a large-scale sparse matrix, and can be effectively replaced by a set of Point Spread Functions (PSFs). In 2005 Xie et al, in the "Wave-equalization-based diagnosis analysis" published by Geophysics, studied the relationship between point spread function and Wave number domain illumination. "An effective method for broadband semiconductor irradiation and resolution analysis" published in SEG by Chen and Xie2015 studies the influence of the point spread function of the acquisition system in a given macroscopic model on the imaging resolution by using a method of calculating point-spread source imaging.

A paper 'Fast least-square migration with a deblocking filter' published by Aoki, Schuster2009 in Geophysics obtains a traditional Hessian matrix by performing a round of migration on a group of pulse scatterers, but the Hessian matrix in the mode is large in scale and high in calculation cost. A method for obtaining a Hessian matrix by utilizing PSFs (particle swarm optimization) approximation is proposed in a paper 'Inversion after depth imaging' published by Fletcher in 2012, wherein the method greatly reduces the scale of the Hessian matrix. The paper "Image domain least-squares simulation with a fluidic matrix estimated by non-stationary simulation filters" published in Journal of Geophysics and Engineering by Guo and Wang2019 also resulted in PSFs by performing a round of migration on the conventional migration results. However, the traditional PSFs estimation method has higher calculation cost, and the invention provides an effective PSFs estimation method on the basis of the prior art, so that the calculation amount can be effectively reduced.

Disclosure of Invention

In view of the above, the present invention provides a method for calculating and storing a Hessian matrix based on a point spread function, including the following steps:

s1, carrying out ray tracing and calculating the travel time of the underground space based on the Gaussian beam and according to the background speed;

s2, obtaining PSFs by adopting an asymptotic Green function under the short wave approximation method WKBJ approximation and approximate estimation in a local window;

s3, the offset image is corrected by using the obtained PSFs, and a corrected image is obtained, where the formula is as follows:

m＝H^-1I_mig

wherein, the offset imaging result is marked as I_migThe ideal image is recorded as m, H^-1Is the inverse of the Hessian matrix.

The technical scheme provided by the invention has the beneficial effects that: the invention introduces approximate Green's function under short wave approximation (WKBJ) and calculates Hessian matrix in local window, thus improving calculation efficiency and reducing calculation scale.

Drawings

FIG. 1 is a flow chart of a calculation and storage method of a Hessian matrix based on a point spread function according to the present invention;

FIG. 2 is a PSF estimated under an asymptotic Hessian/PSF calculation strategy;

FIG. 3 is a schematic diagram of Hessian calculation;

FIG. 4 is a PSF calculation diagram;

FIG. 5 asymptotic Hessian/PSF calculation strategy and estimated PSF;

FIG. 6 is a flow chart of the present invention for calculating a travel in a subterranean space;

FIG. 7 is a flow chart of the present invention for obtaining PSFs by approximate estimation;

FIG. 8 is a diagram of a real reflectivity image of the sigma 2a model;

FIG. 9 is a sigma 2a model offset image;

fig. 10 is an image of the sigma 2a model corrected by the PSF.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, embodiments of the present invention will be further described with reference to the accompanying drawings.

Referring to fig. 1, the invention provides a calculation and storage method of a Hessian matrix based on a point spread function, which is applied to the aspects of illumination analysis, least square offset and the like in seismic exploration data processing;

in the imaging problem of seismic exploration, the relationship of seismic data to a model of reflection coefficients is non-linear.

d＝L(m) (1)

Wherein d is seismic data, m is a model of reflection coefficients, and L is a propagation operator.

This non-linear process can be approximated as linear under a first-order Born approximation, resulting in

d＝Lm (2)

The model for estimating the reflection coefficient of the given seismic observation data is a typical seismic inverse problem and can be solved by the following objective function;

in the formula (d)_obsFor seismic observation data, the first-order partial derivative of the objective function to the model is the gradient, the second-order partial derivative is the Hessian operator, let the gradient be zero, we can get the following equation:

L^Hd＝L^HLm (4)

wherein L is^Hd is the conventional offset imaging result and is denoted as I_mig，L^HL is Hessian operator and is denoted as H, and the above equation can be written as follows:

I_mig＝Hm (5)

the Hessian matrix can be regarded as a fuzzy operator, and then the ideal image m and the offset result I are combined_migIn connection with the above formula, in order to find the ideal image m, the above formula can be:

m＝H^-1I_mig (6)

therefore, the calculation of the Hessian matrix is very important, and when the linearization imaging problem is considered, the Hessian matrix can be expressed as a double integral of the frequency ω and the shot detection ξ, wherein the kernel of the integral is a scattering green function, see formula (2).

H(i,j)＝∫_ω∫ξ[-ω²G(ω,x_r；x_i)G(ω,x_i；x_s)f(ω；s_x)]^*×[-ω²G(ω,x_r；x_j)G(ω,x_j；x_s)f(ω；s_x)]dξdω

＝∫_ω∫ξω⁴|f(ω；s_x)|²G^* _scatter(ω,x_r；x_i；x_s)G_scatter(ω,x_r；x_j；x_s)dξdω (7)

Wherein G represents the Green function, G_scatterRepresents the scattered wave green function, and T represents the travel time; xi represents the index of the shot-to-shot pair, f (omega; s)_x) Representing source wavelet, ω frequency, x_rRepresenting the coordinates of the demodulator probe, x_sRepresenting the coordinates of the shot point, x_j，x_iRepresenting the subsurface space point coordinates.

The ideal image m can be recovered according to the formula, but the inversion including the Hessian matrix is non-trivial and faces a challenge with huge calculation amount. As can be seen from equation (2), the calculation cost is derived from two aspects:

the Green function is calculated in a large quantity, the frequency domain Green function is a high-dimensional function with frequency, shot point and underground space as independent variables, and the Green function is shown in a formula:

wherein the content of the first and second substances,representing the initial wave vector direction, u, of the local plane wave at the shot point_GB(x,x_s,p^sω) represents a gaussian beam in a cartesian coordinate system, where x represents the imaging point coordinates and x represents the imaging point coordinates_sIndicating shot coordinates, ω being the frequency of the seismic waves, and i being a complex number.

Large-scale calculation of the Hessian matrix, because the Hessian matrix is a large-scale sparse matrix, and the size of the Hessian matrix is the model dimension, for example, the dimension of the underground space model is 500 × 500, refer to fig. 2;

the present invention addresses the above difficulties in two main respects:

1) the Hessian is calculated in the local window, and can be obtained from the formula (2), wherein one element of the Hessian is an inner product of products of all shot-inspection pairs of Green functions of a certain point in the model and products of all shot-inspection pairs of Green functions of the certain point, and reference is made to fig. 3. In the model, when calculating the Hessian of a certain point, the point farther from the point has a small value of the dot inner product with the point, and can be ignored, so that n surrounding points closer to the point can be selected to be the dot inner product with the point, the larger element in the Hessian of the row can be calculated, and the size of the Hessian matrix is obviously reduced by ignoring the element close to zero, referring to fig. 4.

2) The integral equation (2) can be simplified by using the asymptotic green's function under WKBJ approximation, see equation (9), and thus the analytic form of Hessian is expressed in equation (4) as the integral of the ray amplitude and source wavelet autocorrelation by the approximation. Due to the autocorrelation function R_ssThe attenuation characteristic of (2), the intensity of which is a function of x_iAnd x_jThe increase in distance is attenuated, which also coincides with the local support of the PSF. The Hessian matrix calculation formula under the WKBJ approximation is converted from double integration to single integration, an integration kernel is changed into an autocorrelation function for controlling the ray amplitude and the propagation time, and the propagation time and the amplitude of the ray are low-dimensional functions relative to a Green function, so that the operation scale is greatly reduced. The arrows in diagram (a) of FIG. 5 represent the ray paths (x) in the desired equation (10)_r←x_i←x_s) And (x)_r←x_j←x_s)。

Therefore, the method comprises the following specific steps:

s1, carrying out ray tracing and calculating the travel time of the underground space based on the Gaussian beam and according to the background speed; the method comprises the following specific steps:

s11, applying for the construction structure array Cells (Nz, Nx, Npx) to store travel time in different directions, and applying for the travel time array to store single shot time (Nz, Nx, Nshot); defining Cells (Nz, Nx, Npx), time (Nz, Nx, Nshot) and related cycle-related variables (such as Nz, Nx, Npx, Nshot, Ix, Iz, Ipx, etc.);

s12, shot point circulation, namely Nshot circulation is carried out on shot points;

s13, shot direction circulation, namely, Npx times of circulation is carried out on shot direction shot, forward modeling is carried out on a Gaussian beam in one direction ipx to obtain ray parameters, and Cells (Nz, Nx, Npx) in the direction are obtained;

s14, circulating model space for 1-Nx NZ times; accumulating traveling time in different directions for IZ and Ix grid points in the model space to obtain traveling time (IZ, Ix and Ishot) of a single cannon, and referring to FIG. 6;

wherein z and x refer to two-dimensional coordinates of a space point, Nshot refers to the cycle number of a shot point, and shot refers to a shot point; wherein I ═ 1, N.

S2, obtaining PSFs by adopting an asymptotic Green' S function under the short wave approximation method WKBJ approximation and approximation estimation in a local window, wherein the method comprises the following steps:

s21, application Hessian array H (Ni, Nj)

S22, shot point circulation, namely performing Nshot circulation on the shot point Ishot;

s23, circulating the detection points, processing any detector Ireceiver, and circulating the detection points for 1-Ntrace times;

s24, model space circulation, wherein circulation is performed for 1-Ni Nj times, Ni Nj is the size of the Hessian matrix model space, and the model space integral is output:

wherein R is_ssAs an autocorrelation function, A_rayRepresentative amplitude can be obtained by travel time, T represents travel time, xi represents index of shot-geophone pair, f (omega; s)_x) Representing source wavelet, ω frequency, x_rRepresenting the coordinates of the demodulator probe, x_sRepresenting the coordinates of the shot point, x_j，x_iRepresenting the subsurface space point coordinates, see FIG. 7;

s3, the offset image is corrected by using the obtained PSFs, and a corrected image is obtained, where the formula is as follows:

m＝H^-1I_mig

wherein, the offset imaging result is marked as I_migThe ideal image is denoted as m.

In order to test the effectiveness of the invention, a sigbee 2a model is selected for testing; fig. 8 shows the true reflection coefficient of the sigma 2a model, fig. 9 shows an offset image obtained by the sigma 2a model with a background velocity, and fig. 10 shows an image obtained by correcting the offset image by the sigma 2a model with PSFs. It can be seen that the image position obtained by correcting the offset image by the Hessian matrix calculated by the method is clear, and the energy is relatively balanced.

The previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.