阅读说明：本技术 一种混合域地震偏移海森矩阵估计方法 (Mixed domain seismic migration hessian matrix estimation method ) 是由 杨继东 张东林 黄建平 李振春 徐洁 秦善源 于 2021-08-23 设计创作，主要内容包括：本发明公开了一种混合域地震偏移海森矩阵估计方法,涉及勘探地球物理技术领域,该方法首先利用常规偏移方法获得初始成像结果,然后将该偏移结果作为反射率模型,进行线性正演和偏移获得二次成像结果,使用这两个偏移结果在空间-波数域估计非稳相滤波器,用于近似汉森矩阵,并利用混合域反褶积计算汉森矩阵的逆,再使用估计的汉森矩阵的逆在最小二乘偏移中作为泛函的梯度,以提高反演的收敛速率,使得最小二乘反演在4-5次迭代内收敛至准确解。本发明的有益效果是,该混合域地震偏移汉森矩阵估计方法既具有较高的计算效率,又可加快最小二乘偏移收敛速度,降低最小二乘反射率反演的计算成本,具有较好的应用前景。(The invention discloses a mixed domain seismic migration hessian matrix estimation method, which relates to the technical field of exploration geophysics, and is characterized in that an initial imaging result is obtained by utilizing a conventional migration method, then the migration result is used as a reflectivity model, linear forward modeling and migration are carried out to obtain a secondary imaging result, the two migration results are used for estimating a non-stationary phase filter in a space-wave number domain to approximate a hessian matrix, the inverse of the hessian matrix is calculated by utilizing mixed domain deconvolution, and the estimated inverse of the hessian matrix is used as a functional gradient in least square migration to improve the convergence rate of inversion, so that the least square inversion converges to an accurate solution in 4-5 iterations. The mixed domain seismic migration Hansen matrix estimation method has the advantages of being high in calculation efficiency, capable of accelerating the least square migration convergence speed, reducing the calculation cost of least square reflectivity inversion and good in application prospect.)

1. A mixed domain seismic migration Hessian matrix estimation method is characterized by comprising the following steps:

step one, acquiring input data;

calculating a seismic source wave field by solving a wave equation, calculating a wave field of a wave detection point, and calculating a preliminary imaging result by adopting a zero-delay cross-correlation imaging condition;

step three, taking the preliminary imaging result as a reflectivity model, and calculating and synthesizing seismic data by using linear forward modeling;

step four, using the seismic data and the migration velocity model synthesized in the step three, adopting the imaging processing process of the step two to carry out migration again, and obtaining a secondary imaging result;

fifthly, applying space-wave number domain transformation to the primary imaging result and the secondary imaging result to obtain a mixed domain imaging result;

sixthly, calculating a non-stationary phase mixed domain filter approximate to a Hansen matrix;

step seven, calculating the inverse of the Hansen matrix by using mixed domain deconvolution;

step eight, the inverse of the estimated Hansen matrix is used as a preconditioner and applied to the functional gradient of least square offset to improve the convergence rate of the inverse;

and step nine, applying the preconditioner in the step eight in each least square iteration to ensure that the least square offset is rapidly converged to an accurate solution.

2. The method of mixed-domain seismic migrated hessian matrix estimation as claimed in claim 1, wherein in step one, the input data comprises: offset velocity model v (x), source function wavelet f (t), field observation data d_obs(x_S,x_RT), mixed domain filter window length l, x_RIs the position of the detection point.

3. The method for mixed-domain seismic migrated hessian matrix estimation as claimed in claim 2, wherein step two comprises:

(2.1) calculating the source wavefield U by solving the wave equation using the source function wavelets f (t) and the offset velocity model v (x)_S(x, t), the wave equation is as follows:

wherein x is the underground imaging spatial position, x_SRepresenting the position of a seismic source, t representing a time variable, and delta being a dirac delta function;

(2.2) Using field observations d_obs(x_S,x_rT) and a migration velocity model v (x), calculating a wave field U of the wave detection point_R(x, t), the specific expression is as follows:

(2.3) calculating a preliminary imaging result m by adopting a zero-delay cross-correlation imaging condition₀(x) And cross-correlation imagingThe conditions are as follows:

wherein T is the length of the seismic record, U_S(x, t) is the source wavefield, U_R(x, t) the demodulator probe wavefield.

4. The mixed-domain seismic migration hessian matrix estimation method of claim 3, wherein in step three, the preliminary imaging result m is processed₀(x) As a reflectivity model, synthetic seismic data d is calculated_syn(x_S,x_rT) the linear forward modeling expression employed is as follows:

wherein, U₀(x, t) is the background wavefield, U₁(x, t) is the perturbing wavefield, δ (x-x)_R) Is a dirac delta function.

5. The mixed-domain seismic migrated hessian matrix estimation method of claim 4, where in step five, the primary imaging result m is subjected to₀(x) And secondary imaging result m₁(x) Applying a space-wavenumber domain transform, the space-wavenumber domain transform expression being:

wherein x is [ x, y, z ═ x, y, z]^TFor space vector, k ═ k_x,k_y,k_z]^TIs a wavenumber vector, x '═ x', y ', z']^TFor space vectors after space-wavenumber domain transformation, M₀(x',k)、M₁(x', k) are each m₀(x)、m₁(x) Result of the space-wavenumber domain transformation of σ_x、σ_yAnd σ_zSet σ for the window lengths in the x, y and z directions, respectively_x、σ_yAnd σ_zThe value of (d) is the input mixed domain filter window length l.

6. The method of mixed-domain seismic migrated hessian matrix estimation as claimed in claim 5, where in step six, M is used₀(x', k) and M₁(x ', k) computing a mixture of non-stationary domain filter F (x', k) approximating the Hansen matrix, expressed as follows:

wherein epsilon is 0.001 max | M₁(x',k)|。

7. The method of claim 6, wherein in step seven, the inverse H of the Hansen matrix_inv(x', k) is:

8. the mixed-domain seismic migration hessian matrix estimation method of claim 7, wherein in step eight, the preconditioner is as follows:

wherein g (x', k) is the space-wavenumber domain spectrum of the functional gradient, g_pre(x) Is the gradient after applying the preconditions.

Technical Field

The invention relates to the technical field of exploration geophysics, in particular to a mixed domain seismic migration hessian matrix estimation method.

Background

At present, high-precision seismic migration imaging is one of important key technologies in oil and gas exploration, and a conventional seismic migration method comprises the following steps: kirchhoff migration, gaussian beam migration, one-way wave migration, reverse time migration, etc. can be viewed mathematically as a adjoint operator of a linear forward modeling, not an inverse operator in the true sense. Therefore, conventional migration methods based on adjoint operators have difficulty obtaining high quality imaging results when seismic acquisition is irregular, the source and detector bands are narrow, and subsurface illumination is not uniform. The least square migration based on the inversion framework can gradually approximate the inversion of a forward operator through data fitting, so that the amplitude balance and the spatial resolution of a migration result are greatly improved, and the method is a hotspot of research in the field of seismic imaging at present.

The core of least squares offset imaging is how to efficiently compute and invert the hansen matrix. The conventional data domain least square migration firstly constructs an error functional of observation data and prediction data, then uses linear forward modeling to calculate the prediction data, uses adjoint migration to calculate functional gradient, and approaches the inverse of a Hansen matrix step by step in an iterative mode. The method usually needs dozens of times to hundreds of forward modeling and accompanying deviation, has low calculation efficiency and is difficult to be applied to large-scale actual production processing. The convergence rate of least square deviation can be accelerated by designing a reasonable preconditioner, but the conventional seismic source illumination compensation can only improve the relative strength relation of amplitude and cannot improve the spatial resolution.

In order to make the least square offset move to the productive application, it is necessary to develop precise and efficient preconditioner so that it can converge to the vicinity of the true solution with a few iterations.

Disclosure of Invention

In order to solve the technical problems of low convergence rate and low calculation efficiency in the least square migration imaging, the invention discloses a mixed domain seismic migration hessian matrix estimation method.

In order to achieve the purpose, the invention adopts the following technical scheme:

a mixed domain seismic migration Hessian matrix estimation method specifically comprises the following steps:

step one, acquiring input data;

calculating a seismic source wave field by solving a wave equation, calculating a wave field of a wave detection point, and calculating a preliminary imaging result by adopting a zero-delay cross-correlation imaging condition;

step three, taking the preliminary imaging result as a reflectivity model, and calculating and synthesizing seismic data by using linear forward modeling;

step four, using the seismic data and the migration velocity model synthesized in the step three, adopting the imaging processing process of the step two to carry out migration again, and obtaining a secondary imaging result;

fifthly, applying space-wave number domain transformation to the primary imaging result and the secondary imaging result to obtain a mixed domain imaging result;

sixthly, calculating a non-stationary phase mixed domain filter approximate to a Hansen matrix;

step seven, calculating the inverse of the Hansen matrix by using mixed domain deconvolution;

step eight, applying the inverse of the estimated Hansen matrix as a preconditioner to a functional gradient of least square offset to improve the convergence rate of the inversion;

and step nine, applying the preconditioner in the step eight in each least square iteration to ensure that the least square offset is rapidly converged to an accurate solution.

Further, in the first step, inputting data includes: offset velocity model v (x), source function wavelet f (t), field observation data d_obs(x_S,x_RT), mixed domain filter window length l, x_RIs the position of the detection point.

Further, the second step comprises:

(2.1) calculating the source wavefield U by solving the wave equation using the source function wavelets f (t) and the offset velocity model v (x)_S(x, t), wave equations such asThe following:

wherein x is the underground imaging spatial position, x_SRepresenting the source position, t representing a time variable, δ (x-x)_s) Is a dirac delta function.

(2.2) Using field observations d_obs(x_S,x_rT) and a migration velocity model v (x), calculating a wave field U of the wave detection point_R(x, t), the specific expression is as follows:

(2.3) calculating a preliminary imaging result m by adopting a zero-delay cross-correlation imaging condition₀(x) And the cross-correlation imaging conditions are:

wherein T is the length of the seismic record, U_S(x, t) is the source wavefield, U_R(x, t) the demodulator probe wavefield.

Further, in step three, the preliminary imaging result m is obtained₀(x) As a reflectivity model, synthetic seismic data d is calculated_syn(x_S,x_rT) the linear forward modeling expression employed is as follows:

wherein, U₀(x, t) is the background wavefield, U₁(x, t) is the perturbing wavefield.

Further, in step five, the primary imaging result m is subjected to₀(x) And secondary imaging result m₁(x) Applying a space-wavenumber domain transform, the space-wavenumber domain transform expression being:

wherein x is [ x, y, z ═ x, y, z]^TFor space vector, k ═ k_x,k_y,k_z]^TIs a wavenumber vector, x '═ x', y ', z']^TFor space vectors after space-wavenumber domain transformation, M₀(x',k)、M₁(x', k) are each m₀(x)、m₁(x) Result of the space-wavenumber domain transformation of σ_x、σ_yAnd σ_zSet σ for the window lengths in the x, y and z directions, respectively_x、σ_yAnd σ_zThe value of (d) is the input mixed domain filter window length l.

Further, in step six, M is used₀(x', k) and M₁(x ', k) computing a mixture of non-stationary domain filter F (x', k) approximating the Hansen matrix, expressed as follows:

wherein, to avoid the small value of zero, let ε be 0.001 max | M₁(x',k)|。

Further, in step seven, the inverse H of the Hansen matrix_inv(x', k) is:

further, in step eight, the preconditioner is as follows:

wherein g (x', k) is the space-wavenumber domain spectrum of the functional gradient, g_pre(x) Is the gradient after applying the preconditions.

The method has the advantages that firstly, the Hansen matrix analyzed in the uniform medium is calculated, and the analysis finds that the point spread function formed when the column vector of the matrix is projected to the imaging space has the localization characteristic. And then the inverse of the Hansen matrix is used as a preconditioner of least square offset, so that the convergence rate of inversion can be greatly improved, and the least square inversion can be converged to an accurate solution in 4-5 iterations. The mixed domain seismic migration Hansen matrix estimation method has high calculation efficiency, can accelerate the least square migration convergence rate, reduces the calculation cost of least square reflectivity inversion, and has good application prospect.

Drawings

FIG. 1 is a diagram of a partial Sigsbee velocity model;

FIG. 2 is an imaging result m obtained by a conventional offset method₀(x)；

FIG. 3 is a graph consisting of m₀(x) As a reflectance model, an imaging result m obtained by performing linear forward modeling and re-shifting₁(x)；

FIG. 4 is a pair m₀(x) Applying an inverse of the hansen matrix to obtain a corrected offset result;

FIG. 5 is an inversion imaging result obtained after 5 times of iteration updating of conventional least square offset

FIG. 6 shows inversion imaging results obtained after the optimized least squares offset is iteratively updated 5 times using the Hansen matrix of the present invention as a preconditioner.

Detailed Description

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.

Based on a linear inversion framework, least squares migration can eliminate migration artifacts due to acquisition irregularities and can improve deep amplitude equalization and overall spatial resolution. However, in the conventional data domain least square offset, the inverse of the hansen matrix is gradually approximated by iterative data fitting, in practical application, dozens of forward simulations and accompanying offsets are often required, the calculation amount is large, and the practical production application is still difficult to be realized under the current calculation condition.

The invention discloses a mixed domain seismic migration Hessian matrix estimation method aiming at the problems of low convergence rate and low calculation efficiency in least square migration imaging, which specifically comprises the following steps:

(1) input data is acquired.

The input data includes: offset velocity model v (x), source function wavelet f (t), field observation data d_obs(x_S,x_RT), mixed domain filter window length l, x_RIs the position of the detection point.

(2) And calculating a seismic source wave field by solving a wave equation, calculating a wave field of a wave detection point, and calculating a preliminary imaging result by adopting a zero-delay cross-correlation imaging condition.

Specifically, the method comprises the following steps:

(2.1) calculating the source wavefield U by solving the wave equation using the source function wavelets f (t) and the offset velocity model v (x)_S(x, t), the wave equation is as follows:

wherein x is the underground imaging spatial position, x_SRepresenting the source position, t representing a time variable, δ (x-x)_s) Is a dirac delta function.

(2.2) Using field observations d_obs(x_S,x_rT) and a migration velocity model v (x), calculating a wave field U of the wave detection point_R(x, t), the specific expression is as follows:

(2.3) calculating a preliminary imaging result m by adopting a zero-delay cross-correlation imaging condition₀(x) And the cross-correlation imaging conditions are:

wherein T is the length of the seismic record, U_S(x, t) is the source wavefield, U_R(x, t) the demodulator probe wavefield.

(3) And taking the preliminary imaging result as a reflectivity model, and calculating and synthesizing the seismic data by using linear forward modeling.

The preliminary imaging result m₀(x) As a reflectivity model, synthetic seismic data d is calculated_syn(x_S,x_rT) the linear forward modeling expression employed is as follows:

wherein, U₀(x, t) is the background wavefield, U₁(x, t) is the perturbing wavefield.

(4) Using the seismic data d synthesized in step (3)_syn(x_S,x_rT) and a migration velocity model, and adopting the imaging processing process of the step (2) to carry out secondary migration to obtain a secondary imaging result m₁(x)；

(5) Applying space-wave number domain transformation to the primary imaging result and the secondary imaging result to obtain a mixed domain imaging result M₀(x', k) and M₁(x',k)。

For the primary imaging result m₀(x) And secondary imaging result m₁(x) Applying a space-wavenumber domain transform, the space-wavenumber domain transform expression being:

wherein x is [ x, y, z ═ x, y, z]^TFor space vector, k ═ k_x,k_y,k_z]^TIs a wavenumber vector, x '═ x', y ', z']^TFor space vectors after space-wavenumber domain transformation, M₀(x',k)、M₁(x', k) are each m₀(x)、m₁(x) Result of the space-wavenumber domain transformation of σ_x、σ_yAnd σ_zThe window lengths in the x, y and z directions, respectively, are set in the method_x、σ_yAnd σ_zThe value of (d) is the input mixed domain filter window length l.

(6) A non-stationary phase mixed domain filter approximating a hansen matrix is calculated.

Using M₀(x', k) and M₁(x ', k) computing a mixture of non-stationary domain filter F (x', k) approximating the Hansen matrix, expressed as follows:

wherein, to avoid the small value of zero, let ε be 0.001 max | M₁(x',k)|。

(7) The inverse of the hansen matrix is calculated using mixed domain deconvolution.

Inverse H of Hansen matrix_inv(x', k) is:

(8) and applying the inverse of the estimated Hansen matrix as a preconditioner to the functional gradient of least square offset to improve the convergence rate of the inversion.

The preconditioner is as follows:

wherein g (x', k) is the space-wavenumber domain spectrum of the functional gradient, g_pre(x) Is composed ofThe pre-conditioned gradient is applied.

(9) The preconditioner in step eight is applied in each least squares iteration so that the least squares offset converges quickly to an accurate solution.

The method is applied to the seismic migration imaging of the international standard Sigsbee model, and a better imaging result is obtained.

FIG. 1 is a partial Sigsbee velocity model diagram, the seismic source function of the observed data is a 15Hz Rake wavelet, the observed data set is 50 shot sets, the shot sets are uniformly distributed on the ground surface, each shot set comprises 500 wave detection points, and the interval is 10 m. FIG. 2 is a graph of offset imaging results m obtained using observation data₀(x) It can be seen that the conventional offset is good for shallow reflector patterning, but for deep sub-salt deposits, it is difficult to achieve high quality imaging results due to the salt dome energy shielding effect. FIG. 3 is a graph showing the result m of imaging using the image shown in FIG. 2₀(x) As a reflectivity model, a quadratic shift result m obtained by performing linear forward modeling and shifting₁(x) Under the influence of the hansen matrix, the amplitude of the deep part of the secondary offset result is weaker and the resolution is poorer compared with the initial offset result. The hansen matrix is estimated in the space-wave number domain by using the offset results shown in fig. 2 and 3, and is applied to the initial result to obtain the imaging result of fig. 4, so that it can be seen that the deep energy is obviously compensated; in addition, the overall spatial resolution is also greatly improved. Figures 5 and 6 compare the imaging results obtained after 5 iterative updates of conventional least squares and least squares offset using the hansen matrix as a precondition proposed by the present patent. It can be seen that the conventional least squares, due to the slower convergence, does not improve the imaging results much after 5 iterations compared to the imaging results in fig. 2. However, after the hansen matrix is used as the preconditioner, the imaging quality of the salt dome and the salt deposit layer is greatly improved after 5 times of iteration.

The invention provides a mixed domain seismic migration Hansen matrix estimation method, which mainly solves the problems of low convergence rate and low calculation efficiency in conventional least square migration, and can make the least square migration converge to an accurate solution after 4-5 iterations by using an estimated Hansen matrix as a preconditioner of a functional gradient in the least square migration.

The method comprises the steps of firstly obtaining an initial imaging result by using a conventional migration method, then using the migration result as a reflectivity model, carrying out linear forward modeling and migration to obtain a secondary imaging result, using the two migration results to estimate a non-stationary phase filter in a space-wavenumber domain for approximating a Hansen matrix, calculating the inverse of the Hansen matrix by using a mixed domain deconvolution, and using the estimated inverse of the Hansen matrix as the gradient of a functional in least square migration, so that the convergence rate can be greatly improved, the convergence rate can be converged to an accurate solution in 4-5 iterations, and finally the practical process of the least square migration is promoted.

It is to be understood that the above description is not intended to limit the present invention, and the present invention is not limited to the above examples, and those skilled in the art may make modifications, alterations, additions or substitutions within the spirit and scope of the present invention.