CN102116869A - High-precision prestack domain least square migration seismic imaging technology - Google Patents

Abstract

The invention provides a high-precision prestack domain least square migration seismic imaging technology. The technology comprises the following steps of: giving an initial speed model; calculating a travel-time table according to a ray tracing method and storing the travel-time table in a memory to be used in a subsequent migration process; calculating a forward modeling wave field by using a high-frequency approximation theory; calculating residual error between the forward modeling wave field and a recording wave field; when the residual error is less than a previously given error standard, stopping iteration and outputting an image result; when deviation is greater than the previously given error standard, shifting the obtained difference data, calculating a correction coefficient and updating the imaging result, and further repeating the operation until the calculation error is less than the given error standard; and outputting the imaging result. The technology can be effectively applied to small-size complex reservoir models, and imaging processing and theoretical research of weak signals on intermediate and deep parts, and has great application prospect in future exploration and development of exploration areas in western China.

Description

High precision prestack territory least square skew seismic imaging technology

Affiliated technical field

The invention belongs to seismic data digital processing field, is to utilize least-squares algorithm to carry out a kind of imaging technique of underground skew.

Background technology

Least square method has obtained extensive utilization as a kind of optimization algorithm in the geophysical research field, and wherein application is comparatively extensive aspect processing filtering, inverting and imaging problem.The least square skew is in early days by (1988) and lambare people such as (1992) propositions such as LeBrad.Nemeth people such as (1999) has further set forth the relative merits of least square migration algorithm from the derivation of method principle subsequently, even and adopt Huygens' principle to explain that least square was offset the mechanism that still can obtain better offset effect when data were imperfect.Then, what Dequet people such as (2000) proposed is handling the relief surface illumination and is having greater advantages, the model boundary zone that especially complicated salt dome bottom and degree of covering are less by the image error that the seismic wave field of irregular thick sampling causes than the Kirchhoff skew based on the least-squares algorithm of ray theory Green equation.Adopting Kirchhoff also is a kind of seismic imaging algorithm that can effectively reduce the false picture of skew as the least square skew of positive propagation operator, in the imaging iterative process, the general conjugated body gradient method that adopts is carried out iteration, part scholar also proposes to adopt the conjugated gradient direction of seeking optimum to come the Approximate Equivalent method of conjugate gradient, further accelerate speed of convergence, improved counting yield.

Kuehl and Sacchi (2001) point out that least square skew can introduce wave equation and just propagate and the anti-spread operator, and have provided least square on this basis and split the step migration algorithm.Jia Xiaofeng etc. have also proposed the calculation process and the algorithm that adopt moving least squares method to find the solution seismic wave equation.Its strong wait (2008) of poplar adopt Fourier finite-difference migration operator to be applied in the poststack least square migration algorithm, and have carried out the examination processing of model.Because the Fourier finite difference operator is that a kind of counting yield is higher, can adapt to violent lateral velocity variation, the adaptability of rate pattern is enhanced.Daiwei etc. (2009) have proposed prestack least square migration algorithm, and use it in the denoising technology of many focus recording geometry seismic imaging, and imaging results shows that the least square skew can be good at suppressing man-made noise and the false picture of skew.

Summary of the invention

The objective of the invention is method principle, and realized prestack territory least square migration algorithm based on the skew of prestack territory high precision least square, and a kind of high precision prestack territory least square skew seismic imaging technology that provides.

The present invention verifies the imaging precision and the practicality of this patent algorithm by the work of three aspects: (1) adopts this algorithm that anomalous body in the three-dimensional experiment model is carried out imaging, and compare with traditional Kirchhoff pre-stack depth migration imaging result, can find by contrast: algorithm of the present invention can better make the diffracted wave energy converge on given anomalous body position in advance, and classic method has higher imaging precision relatively; (2) the complex model imaging tentative calculation of test imaging technique by Mamousi2 and SEG/EAGE-SALT two big standards and algorithm quality as can be seen: from imaging precision and two aspects of imaging resolution, the present invention has imaging results preferably in the complex structure zone, and especially the compacting to man-made noise and the false picture of skew has good effect; (3) carbonatite seam hole, deep preserved the imaging examination processing of body seismogeology equivalent model during the present invention can be used for, and can portray reservoir and the boundary characteristic that carbonatite seam hole preserves body as can be seen preferably from imaging results, and imaging resolution is higher.

The technical solution adopted in the present invention is:

High precision prestack territory least square skew seismic imaging technology, step is:

1) a given initial velocity model is calculated travel timetable and is stored in the internal memory for follow-up migration process use according to ray-tracing scheme;

2) adopt the high-frequency approximation Theoretical Calculation just drilling wave field simultaneously;

3) and then calculate the residual error just drilling wave field and record wave field,, stop iteration and export imaging results when residual error during less than in advance given error criterion;

When 4) deviation is greater than in advance given error criterion, the difference data of gained is offset, and calculates correction factor, the final updated imaging results, and then repeat aforesaid operations up to the error of calculation during, the output imaging results less than given error criterion.

Treatment scheme is:

(1) numerical simulation: the residual error of computational data, utilize formula (a),

r _n＝Lm _n-d

(a)

(2) skew residual error: residual error is offset, and used formula is shown in (b):

g _n＝L ^Tr _n

(b)

(3) calculate step-length: utilize formula (c) to calculate iteration step length,

$k_n=\frac{g_n^T{g}_n}{{\left({\mathrm{Lg}}_n\right)}^T\left({\mathrm{Lg}}_n\right)}$

(c)

(4) new model more: utilize formula (d), calculate the model after upgrading,

m _n+1＝m _n-k _ng _n

(d)

The condition that calculate to stop is divided into two classes: (1) provides certain iterations, and when iterations reached given iterations, iteration stopping was preserved the subsurface model m of this moment _N+1(2) threshold values of assigned error, when error during, iteration stopping less than given threshold values, and the subsurface model m of output after upgrading _N+1

Alternative manner is called the least square migration algorithm, and its content is:

The measured earth model of least square offset method principle supposes that linearity just drilling matrix operator and be

Satisfy following equation (1):

$P=\underset{~}{L}m---\left(1\right)$

Wherein P is the seismologic record vector of simulation, and m is the Reflectivity Model vector of the earth.The seismologic record P that observes ₀Can be expressed as

M wherein ₀Be real albedo of the earth model vector,

For the linearity based on real earth Reflectivity Model is just being drilled matrix operator; In computation process, unless specified otherwise, general supposition

Kirchhoff skew adopts linear transposition of just drilling matrix operator to find the solution underground albedo of the earth model, the result as shown in Equation (2):

$m_k={\underset{~}{L}}^{T}P---\left(2\right)$

M wherein _kBe the Kirchhoff migration result, bring formula (1) into formula (2) and obtain:

$m_k={\underset{~}{L}}^T\underset{~}{L}m---\left(3\right)$

At this moment, Be defined as the Hessian matrix, by formula (3) as can be known: m _kBe equivalent to the result of earth model m through the filtering of Hessian matrix, when

It is a unit matrix vector

The time, by the m of Kirchhoff skew acquisition _kBe exactly real earth model m, yet, in most of the cases,

It or not a unit matrix

It mainly shows as:

Matrix element is unequal on diagonal line; Perhaps element is not 0. in order more clearly to disclose the computing method of Hessian matrix on the off-diagonal, and we have done following equivalence:

Write equation (1) as integrated form, as shown in Equation (4):

P(r，t|s，0)＝∫m(x)W(t)*G(r，t|x，0)*G(x，t|s，0)dx

(4)

In formula (4), G (x, t|s, 0) and G (r, t|x, 0) represent the Green function of focus (s) respectively to scattering point (x) and scattering point (x) to acceptance point (r), W (t) represents source wavelet, people such as Bleistein pointed out in 1984, and Green function can be expressed as formula (5) under the zero-order approximation:

G(x，t|x′，0)＝A _xx′δ(t-τ _xx′)

(5)

Wherein, A _{Xx '}Be separating of transport equation, τ _{Xx '}For between x and the x ' then, consider generalized case, suppose that in recording geometry the sampling function of focus and receiver is respectively h _s(s) and h _r(r),

Bring formula (5) into formula (4), can get:

P(r，t|s，0)＝∫m(x)A _sxA _xrW(t-τ _sx-τ _xr)dx (6)

Consider to excite and receive sampling function h _s(s) and h _r(r), our equation that can obtain Kirchhoff skew can be expressed as:

m _k(x)＝∫h _s(s)ds∫h _r(r)dr×(∫P(r，t|s，0)A _sxA _xrW(t-τ _sx-τ _xr)dt) (7)

Bring formula (6) into formula (7), and the exchange integral order, the Kirchhoff offset equation can be written as:

m _k(x)＝∫m(x′)dx′∫h _s(s)A _sxA _sx′ds∫h _r(r)A _xrA _x′rdr×(∫W(t-τ _sx-τ _xr)W(t-τ _sx′-τ _x′r)dt)

(8)

By formula (8) as can be known, integral kernel K (x, x ') is:

K(x，x′)＝∫h _s(s)ds∫h _r(r)drA _sxA _sx′A _xrA _x′r×R(τ _sx+τ _xr-τ _sx′-τ _x′r) (9)

Wherein, R (τ _Sx+ τ _Xr-τ _{Sx '}-τ _{X ' r}) be source function W in difference simple crosscorrelation constantly, for this reason, (9) are brought into (8), (8) are reduced to:

m _k(x)＝∫K(x，x′)m(x′)dx′ (10)

(10) are write as discrete form, as shown in Equation (11):

$m_k=\underset{~}{K}m---\left(11\right)$

Wherein, by (3) as can be known, matrix

Physical meaning be the Hessian matrix, integral kernel K (x, x ') is equivalent to from m to m _kA mapping, for some scattering point x _f, K (x, x ') represents the system responses of this scattering point.If single track is offset, this response is the oval impulse response of skew; And a large amount of seismic traces is offset, and this response energy will converge to real scattering point position;

When the Hessian matrix was unit matrix, skew model mk equaled spherical model m truly, investigates the result of the diagonal entry of Hessian matrix, and by formula (9) as can be known, the diagonal entry of integral kernel matrix K (x, x ') can be written as following form,

$K(x_0,x_0)=R\left(0\right)\&Integral;h_s\left(s\right)\mathrm{ds}{A}_{{\mathrm{<figure-callout\; id="0"\; label="sx"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">sx</figure-callout>}}_0}^2\&Integral;h_r\left(r\right)\mathrm{dr}{A}_{x_{0}r}^2---\left(12\right)$

By formula (12) as can be known, the element value on the matrix principal diagonal is subjected to

And the noncontinuous sampling function (h of wave field _s(s), h _r(r)) control;

If distributing, focus and receiver caused the image data of non-rule, the Kirchhoff offset method will produce the false picture of skew on migrated section, in order to reduce the false picture of skew to the full extent, obtain comparatively real earth model m, we are given following objective function:

$P\left(m\right)={\left|\right|\underset{~}{L}m-{\mathrm{<figure-callout\; id="0"\; label="P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">P</figure-callout>}}_0\left|\right|}^2+{\mathrm{\&epsiv;}}^2{\left|\right|\underset{~}{C}m-\underset{~}{C}{m}_{\mathrm{apr}}\left|\right|}^2---\left(13\right)$

Wherein, first on formula (13) equation the right is the error function item of data, and second is regularization term, Represent on the model m of a linear matrix operator effect m _AprBe initial reference model, ε ²Be the weight that regularization is handled, equation (13) can be found the solution with the conjugated body gradient method of propositions in 1987 such as Tarantola, and the solving equation of m is:

$m={({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})}^{-1}({\underset{~}{L}}^{\mathrm{<figure-callout\; id="0"\; label="T\; P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">T</figure-callout>}}{P}_{}{}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}})---\left(14\right)$

And equivalent equation, shown in equation (15):

$({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})m={\underset{~}{L}}^{\mathrm{<figure-callout\; id="0"\; label="T\; P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">T</figure-callout>}}{P}_{}{}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}}---\left(15\right)$

Equation (14), the model m in (15) can find the solution with the method for iteration.

Main advantage of the present invention: (1) carries out imaging based on the algorithm of inverting to complex geologic body, can eliminate the false picture of man-made noise and skew preferably; (2) have higher imaging precision and rate pattern adaptability preferably, be useful for especially that carbonatite seam hole, deep preserves body reservoir imaging research in the western exploratory area of China.

Method of the present invention has been considered the normal conditions of Hessian matrix, and unlike the unit matrix that conventional method is supposed, can more accurately model mk be converged to real subsurface model m.In order to compare the difference of the inventive method and conventional offset method imaging effect, we have set up following three-dimensional model to test.

Rate pattern is introduced: the dimension of the three-dimensional model of foundation is respectively nx=ny=nz=75, model three direction sizing grid directions are 20m, adopt 75 gun excitations, 75 roads receive, shot point and acceptance point position as shown in Figure 1,1000 points of single track sampling, source wavelet adopts the ricker wavelet, as shown in Figure 2, dominant frequency is chosen as 15hz.Rate pattern and conventional method and least square migration imaging result are as shown in Figure 3.Contrast as can be known by Fig. 3 (a1) and Fig. 3 (a2), at first investigate horizontal direction, least square skew skew energy is more restrained, energygram is also more concentrated, and not well convergence of diffraction energy in the conventional migration result 3 (a1), there is comparatively significantly not convergent energy ring distribution, will produces the false structure of some artificial skews; In vertical direction, similar imaging results is also arranged, i.e. least square skew has higher imaging precision and resolution with respect to traditional migration algorithm, and the image space of anomalous body also approaches the given anomalous body position of initial model more.

Try result as can be known from the skew of above-mentioned three-dimensional model, the more conventional Kirchhoff prestack depth migration method of the imaging effect of the art of this patent imaging resolution has had tangible improvement, can be effectively applied in the imaging processing and theoretical research of small scale complicated reservoirs model, middle deep feeble signal, in western part of China exploratory area exploratory development from now on, application promise in clinical practice is arranged.

Description of drawings

Fig. 1 is shot point and acceptance point distribution schematic diagram.

Fig. 2 is a Ricker wavelet synoptic diagram.

Fig. 3 is rate pattern and conventional method and least square migration imaging result.Wherein (a) is that level is to initial velocity model, (b) for vertically to initial velocity model; (a1) level is to the Kirchhoff migration result, (b1) vertically to the Kirchhoff migration result; (a2) level is offset migration result to least square, (b2) vertically to least square skew migration result.

Fig. 4 is a least square migration processing process flow diagram.

Fig. 5 is the migration result synoptic diagram of embodiment 1 SEG_EAGE 2D salt dome model velocity field and different least square offset iterations number of times.Wherein (a) is salt dome model synoptic diagram, (b) be 5 results of least square offset iterations, (c) be 10 results of least square offset iterations, (d) be 20 results of least square offset iterations, (e) being 30 results of least square offset iterations, (f) is 50 results of least square offset iterations.

Fig. 6 is an imaging results comparison diagram under the different iterationses of embodiment 2Mamousi2 model.Wherein (a) is 5 results for Mamousi2 model iterations, (b) for Mamousi2 model iterations is 10 results, (c) for Mamousi2 model iterations is 20 results, is 50 results for Mamousi2 model iterations (d).

Fig. 7 is a solution cavity type reservoir-level differential pattern least square migration result at interval.Wherein (a) is solution cavity type reservoir-level differential pattern at interval, (b) is the least square migration result.

Fig. 8 is a solution cavity type reservoir perpendicular separation differential pattern least square migration result.Wherein (a) is solution cavity type reservoir perpendicular separation differential pattern, (b) is the least square migration result.

Fig. 9 is a middle deep carbonate different angle equivalence earthquake geologic model synoptic diagram.Wherein (a) is 30 ° rate pattern for the inclination angle, (b) for the inclination angle being 45 ° rate pattern, (c) is 60 ° rate pattern for the inclination angle, (d) is 70 ° rate pattern for the inclination angle, (e) for the inclination angle being 80 ° rate pattern, (f) is 90 ° rate pattern for the inclination angle.

When Figure 10 is 30 ° for the geology inclination angle, least square migration imaging result schematic diagram under the different iterationses.Wherein (a) (b) is 10 iteration results of least square skew for 5 iteration of least square skew figure as a result, (c) for least square is offset 15 iteration results, is that least square is offset 20 iteration results (d).

Figure 11 is 20 o'clock for iterations, and different geology dip migration results contrast synoptic diagram.Wherein (a) is the migration result at 30 ° of inclination angles, (b) is the migration result at 45 ° of inclination angles, (c) is the migration result at 60 ° of inclination angles, (d) is the migration result at 70 ° of inclination angles, (e) is the migration result at 80 ° of inclination angles, (f) migration result at 90 ° of inclination angles.

Embodiment

High precision prestack territory least square skew seismic imaging technology, step is:

1) a given initial velocity model is calculated travel timetable and is stored in the internal memory for follow-up migration process use according to ray-tracing scheme;

2) adopt the high-frequency approximation Theoretical Calculation just drilling wave field simultaneously;

3) and then calculate the residual error just drilling wave field and record wave field,, stop iteration and export imaging results when residual error during less than in advance given error criterion;

When 4) deviation is greater than in advance given error criterion, the difference data of gained is offset, and calculates correction factor, the final updated imaging results, and then repeat aforesaid operations up to the error of calculation during, the output imaging results less than given error criterion.

Treatment scheme is:

(1) numerical simulation: the residual error of computational data, utilize formula (a),

r _n＝Lm _n-d

(a)

(2) skew residual error: residual error is offset, and used formula is shown in (b):

g _n＝L ^Tr _n

(b)

(3) calculate step-length: utilize formula (c) to calculate iteration step length,

$k_n=\frac{g_n^T{g}_n}{{\left({\mathrm{Lg}}_n\right)}^T\left({\mathrm{Lg}}_n\right)}$

(c)

(4) new model more: utilize formula (d), calculate the model after upgrading,

m _n+1＝m _n-k _ng _n

(d)

The condition that calculate to stop is divided into two classes: (1) provides certain iterations, and when iterations reached given iterations, iteration stopping was preserved the subsurface model m of this moment _N+1(2) threshold values of assigned error, when error during, iteration stopping less than given threshold values, and the subsurface model m of output after upgrading _N+1

Alternative manner is called the least square migration algorithm, and its content is:

The measured earth model of least square offset method principle supposes that linearity just drilling matrix operator and be

Satisfy following equation (1):

$P=\underset{~}{L}m---\left(1\right)$

Wherein P is the seismologic record vector of simulation, and m is the Reflectivity Model vector of the earth.The seismologic record P that observes ₀Can be expressed as

M wherein ₀Be real albedo of the earth model vector,

For the linearity based on real earth Reflectivity Model is just being drilled matrix operator; In computation process, unless specified otherwise, general supposition

Kirchhoff skew adopts linear transposition of just drilling matrix operator to find the solution underground albedo of the earth model, the result as shown in Equation (2):

$m_k={\underset{~}{L}}^{T}P---\left(2\right)$

M wherein _kBe the Kirchhoff migration result, bring formula (1) into formula (2) and obtain:

$m_k={\underset{~}{L}}^T\underset{~}{L}m---\left(3\right)$

At this moment,

Be defined as the Hessian matrix, by formula (3) as can be known: m _kBe equivalent to the result of earth model m through the filtering of Hessian matrix, when

It is a unit matrix vector

The time, by the m of Kirchhoff skew acquisition _kBe exactly real earth model m, yet, in most of the cases, It or not a unit matrix

It mainly shows as:

Matrix element is unequal on diagonal line; Perhaps element is not 0. in order more clearly to disclose the computing method of Hessian matrix on the off-diagonal, and we have done following equivalence:

Write equation (1) as integrated form, as shown in Equation (4):

P(r，t|s，0)＝∫m(x)W(t)*G(r，t|x，0)*G(x，t|s，0)dx

(4)

In formula (4), G (x, t|s, 0) and G (r, t|x, 0) represent the Green function of focus (s) respectively to scattering point (x) and scattering point (x) to acceptance point (r), W (t) represents source wavelet, people such as Bleistein pointed out in 1984, and Green function can be expressed as formula (5) under the zero-order approximation:

G(x，t|x′，0)＝A _xx′δ(t-τ _xx′)

(5)

Wherein, A _{Xx '}Be separating of transport equation, τ _{Xx '}For between x and the x ' then, consider generalized case, suppose that in recording geometry the sampling function of focus and receiver is respectively h _s(s) and h _r(r), bring formula (5) into formula (4), can get:

P(r，t|s，0)＝∫m(x)A _sxA _xrW(t-τ _sx-τ _xr)dx (6)

Consider to excite and receive sampling function h _s(s) and h _r(r), our equation that can obtain Kirchhoff skew can be expressed as:

m _k(x)＝∫h _s(s)ds∫h _r(r)dr×(∫P(r，t|s，0)A _sxA _xrW(t-τ _sx-τ _xr)dt) (7)

Bring formula (6) into formula (7), and the exchange integral order, the Kirchhoff offset equation can be written as:

m _k(x)＝∫m(x′)dx′∫h _s(s)A _sxA _sx′ds∫h _r(r)A _xrA _x′rdr×(∫W(t-τ _sx-τ _xr)W(t-τ _sx′-τ _x′r)dt)

(8)

By formula (8) as can be known, integral kernel K (x, x ') is:

K(x，x′)＝∫h _s(s)ds∫h _r(r)drA _sxA _sx′A _xrA _x′r×R(τ _sx+τ _xr-τ _sx′-τ _x′r) (9)

Wherein, R (τ _Sx+ τ _Xr-τ _{Sx '}-τ _{X ' r}) be source function W in difference simple crosscorrelation constantly, for this reason, (9) are brought into (8), (8) are reduced to:

m _k(x)＝∫K(x，x′)m(x′)dx′ (10)

(10) are write as discrete form, are added shown in the formula (11):

$m_k=\underset{~}{K}m---\left(11\right)$

Wherein, by (3) as can be known, matrix

Physical meaning be the Hessian matrix, integral kernel K (x, x ') is equivalent to from m to m _kA mapping, for some scattering point x _f, K (x, x ') represents the system responses of this scattering point.If single track is offset, this response is the oval impulse response of skew; And a large amount of seismic traces is offset, and this response energy will converge to real scattering point position;

When the Hessian matrix was unit matrix, skew model mk equaled spherical model m truly, investigates the result of the diagonal entry of Hessian matrix, and by formula (9) as can be known, the diagonal entry of integral kernel matrix K (x, x ') can be written as following form,

$K(x_0,x_0)=R\left(0\right)\&Integral;h_s\left(s\right)\mathrm{ds}{A}_{{\mathrm{<figure-callout\; id="0"\; label="sx"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">sx</figure-callout>}}_0}^2\&Integral;h_r\left(r\right)\mathrm{dr}{A}_{x_{0}r}^2---\left(12\right)$

By formula (12) as can be known, the element value on the matrix principal diagonal is subjected to

And the noncontinuous sampling function (h of wave field _s(s), h _r(r)) control;

If distributing, focus and receiver caused the image data of non-rule, the Kirchhoff offset method will produce the false picture of skew on migrated section, in order to reduce the false picture of skew to the full extent, obtain comparatively real earth model m, we are given following objective function:

$P\left(m\right)={\left|\right|\underset{~}{L}m-{\mathrm{<figure-callout\; id="0"\; label="P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">P</figure-callout>}}_0\left|\right|}^2+{\mathrm{\&epsiv;}}^2{\left|\right|\underset{~}{C}m-\underset{~}{C}{m}_{\mathrm{apr}}\left|\right|}^2---\left(13\right)$

Wherein, first on formula (13) equation the right is the error function item of data, and second is regularization term, Represent on the model m of a linear matrix operator effect m _AprBe initial reference model, ε ²Be the weight that regularization is handled, equation (13) can be found the solution with the conjugated body gradient method of propositions in 1987 such as Tarantola, and the solving equation of m is:

$m={({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})}^{-1}({\underset{~}{L}}^{\mathrm{<figure-callout\; id="0"\; label="T\; P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">T</figure-callout>}}{P}_{}{}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}})---\left(14\right)$

And equivalent equation, shown in equation (15):

$({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})m={\underset{~}{L}}^{\mathrm{<figure-callout\; id="0"\; label="T\; P"\; filenames="HSA00000433046300012.png,HSA00000433046300111.png"\; state="\{\{state\}\}">T</figure-callout>}}{P}_{}{}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}}---\left(15\right)$

Equation (14), the model m in (15) can find the solution with the method for iteration.

Based on the said method principle, state in realization on the basis of least square offset method, imaging test by SEG_EAGE salt dome and Mamousi2 model, verified the imaging effect of least square skew, the imaging that simultaneously this method is used for the deep carbonate reservoir is calculated, and preserves effect in the imaging of body seismic reservoir in western carbonatite seam hole to check this method.

Embodiment (one) model 1: the imaging examination result of having tested the 2D salt dome model of SEG_EAGE.As shown in Figure 5, the characteristics of this model are that the speed horizontal change is violent, and the velocity amplitude of salt body is that the twice of sedimentary deposit is many, and the inclination angle on salt body border is big.Testing used salt dome model geological data is the disclosed public in the world normal data of SEG_EAGE, and recording geometry is to blow out in the right, and the left side receives, and has 325 big guns, every big gun has 176 roads, shot interval 48m, track pitch 24m, the smallest offset distance is zero, record length 5s, time-sampling rate 8ms.Least square migration result under the velocity field that Fig. 5 is respectively the salt dome model and the different iterationses, as can be seen along with the increase of iterations, salt dome flank and bottom imaging effect are more and more accurate, salt dome bottom Geological Structural Forms can be sketched the contours of clearly during iteration 20 times, iteration is more than 30 times, and the boundary position imaging effect of salt dome flank and minor fault has also had tangible improvement.

The method of the technology of the present invention in migration process, adopts the method for iteration, constantly revises migration operator, behind the certain number of times of iteration, can realize the playback of reflection wave and the convergence of diffracted wave preferably, has imaging effect preferably.

The imaging examination result of embodiment (two) model 2:Mamousi2 model.Fig. 6 has provided the comparison diagram of least square migration result under another test migration algorithm correctness and adaptive rate pattern (Mamousi2) and the different iterations, the contrast imaging results is as can be known: iterations hour, the least square offset method can provide the profile and the main structural attitude of Mamousi2 rate pattern, increase along with iterations, middle deep energy is further recovered, and imaging resolution improves; Mamousi2 comes out with respect to the initiate low speed body of Mamousi model also imaging comparatively accurately simultaneously, and when iterations reached 50 times, the imaging results of least square skew can reflect the principal character of Mamousi2 rate pattern preferably.In sum: the least square migration algorithm can be used in the migration imaging processing of comparatively complicated rate pattern, and the imaging precision of this method is higher.

We have come out the position with respect to the initiate low-velocity anomal body of Mamousi model among the Mamousi2 with red circle and arrow logo, from imaging examination result as can be seen: the formation method of the art of this patent, even when model is comparatively complicated, also can carry out imaging to small scale low-velocity anomal body reservoir preferably, and the precision of imaging and resolution are higher.The art of this patent will play an important role in the work of exploration and development of western exploratory area complex geological structure environment.

Embodiment (three) model 3: carbonatite seam hole preserves body equivalence earthquake geologic model imaging examination result.Fig. 7 and Fig. 8 two kinds of equivalent earthquake model of geological structure body and least square migration result thereof for preserving based on China's western exploratory area carbonatite solution cavity type that body sets up, iterations is 20 times.Equivalent earthquake geologic model shown in Figure 7, to preserve the body yardstick identical for solution cavity in the group, not solution cavity yardstick on the same group from left to right diameter be respectively 20 meters, 30 meters and 40 meters, the solution cavity spacing is identical in the same depth location place group, minimum spacing is 20 meters, and it is increasing to increase the solution cavity spacing with the degree of depth.Solution cavity equivalent model shown in Figure 8 and Fig. 7 are similar, and vertically depth interval is identical in the same lateral position group of this model that different is, 20 meters of minimum longitudinal separations, and in the group from left to right, longitudinal separation is increasing.From the migration imaging result of two group models as can be known, the least square offset method has all reached 20 meters resolution characteristic in length and breadth to carbonate reservoir to resolution, the diffracted wave energy converges on the solution cavity position preferably, even 20 meters geologic body also can be portrayed the boundary information of anomalous body clearly.

Fig. 9 has provided the equivalent earthquake geologic model of the different fracture dips of middle deep carbonate fracture reservoir, and wherein the width in crack from left to right is respectively 20 meters, 30 meters, 40 metrical scales between the group; Last interval velocity model size is 3000m/s, and lower floor is 4500m/s, and reservoir is 2300m/s.Imaging results when carbonatite seam hole preserved the body least square and is offset different iterations when Figure 10 was 30 ° for the inclination angle, contrasting different iterations results learns, increase iterations and can improve imaging precision, but after greater than 10 iteration imagings, restrain more and more slower, contrast 10 iteration and 20 iterative least square migration result, both are suitable substantially to the degree of convergence that the crack preserves the body diffraction energy.In actual applications, should weigh efficient and effect, determine final iterations.

Figure 11 shows that the least square migration result comparison diagram of 3 slit die type (dominant frequency 50HZ) different angles, by contrasting six migration result as can be known, though the crack angle is very big to the influence of migration imaging resolution, even but least square skew when 90 ° of inclination angles still in the future the diffracted wave of the overwhelming majority at the end, endokinetic fissure top converge picture, reach 20 meters resolution characteristic.When the crack angle is less, not only can restrain the diffracted wave energy of the overwhelming majority of position, the end, top, crack, also can portray the lateral boundaries in every crack clearly.Yet because group internal fissure close together, interference ratio is more serious mutually in the seismic wave propagation process, and because the restriction of migration aperture, the information of the significant wave energy that receives when reservoir is positioned at high angle is less, therefore during high angle, imaging effect is poor during the imaging effect smaller angle of crack lateral boundaries.

We in imaging results (Figure 11) gone out under the different fracture dips portrayal degree on fracture reservoir border with arrow logo.When fracture dip during greater than 60 °, the vertical border portrayal of fracture reservoir is comparatively fuzzy, but the position of the top bottom boundary imaging in crack is comparatively accurate.This shows: the art of this patent, the carbonatite seam hole of burying deeply (greater than 3500 meters), reservoir yardstick little (20 meters) for reservoir preserves body equivalence earthquake geologic model and has higher imaging resolution and imaging precision, can be provided as the picture technical support for the exploratory development of carbonate reservoir under the western China complex structure environment.

Claims (3)

1. high precision prestack territory least square is offset the seismic imaging technology, it is characterized in that step is:

1) a given initial velocity model is calculated travel timetable and is stored in the internal memory for follow-up migration process use according to ray-tracing scheme;

2) adopt the high-frequency approximation Theoretical Calculation just drilling wave field simultaneously;

3) and then calculate the residual error just drilling wave field and record wave field,, stop iteration and export imaging results when residual error during less than in advance given error criterion;

When 4) deviation is greater than in advance given error criterion, the difference data of gained is offset, and calculates correction factor, the final updated imaging results, and then repeat aforesaid operations up to the error of calculation during, the output imaging results less than given error criterion.

2. high precision prestack according to claim 1 territory least square skew seismic imaging technology is characterized in that treatment scheme is:

(1) numerical simulation: the residual error of computational data, utilize formula (a),

r _n＝Lm _n-d

(a)

(2) skew residual error: residual error is offset, and used formula is shown in (b):

g _n＝L ^Tr _n

(b)

(3) calculate step-length: utilize formula (c) to calculate iteration step length,

$k_n=\frac{g_n^T{g}_n}{{\left({\mathrm{Lg}}_n\right)}^T\left({\mathrm{Lg}}_n\right)}$

(c)

(4) new model more: utilize formula (d), calculate the model after upgrading,

m _n+1＝m _n-k _ng _n

(d)

The condition that calculate to stop is divided into two classes: (1) provides certain iterations, and when iterations reached given iterations, iteration stopping was preserved the subsurface model m of this moment _N+1(2) threshold values of assigned error, when error during, iteration stopping less than given threshold values, and the subsurface model m of output after upgrading _N+1

3. high precision prestack according to claim 1 territory least square skew seismic imaging technology is characterized in that alternative manner is called the least square migration algorithm, and its content is:

The measured earth model of least square offset method principle supposes that linearity just drilling matrix operator and be

Satisfy following equation (1):

$P=\underset{~}{L}m---\left(1\right)$

Wherein P is the seismologic record vector of simulation, and m is the Reflectivity Model vector of the earth.The seismologic record P that observes ₀Can be expressed as

M wherein ₀Be real albedo of the earth model vector,

For the linearity based on real earth Reflectivity Model is just being drilled matrix operator; In computation process, unless specified otherwise, general supposition

Kirchhoff skew adopts linear transposition of just drilling matrix operator to find the solution underground albedo of the earth model, the result as shown in Equation (2):

$m_k={\underset{~}{L}}^{T}P---\left(2\right)$

M wherein _kBe the Kirchhoff migration result, bring formula (1) into formula (2) and obtain:

$m_k={\underset{~}{L}}^T\underset{~}{L}m---\left(3\right)$

At this moment,

Be defined as the Hessian matrix, by formula (3) m as can be known _k: be equivalent to the result of earth model m through the filtering of Hessian matrix, when

It is a unit matrix vector

The time, the mk that is obtained by the Kirchhoff skew is exactly real earth model m, yet, in most of the cases,

It or not a unit matrix

It mainly shows as:

Matrix element is unequal on diagonal line; Perhaps element is not 0. in order more clearly to disclose the computing method of Hessian matrix on the off-diagonal, and we have done following equivalence:

Write equation (1) as integrated form, as shown in Equation (4):

P(r，t|s，0)＝∫m(x)W(t)*G(r，t|x，0)*G(x，t|s，0)dx

(4)

In formula (4), G (x, t|s, 0) and G (r, t|x, 0) represent the Green function of focus (s) respectively to scattering point (x) and scattering point (x) to acceptance point (r), W (t) represents source wavelet, people such as Bleistein pointed out in 1984, and Green function can be expressed as formula (5) under the zero-order approximation:

G(x，t|x′，0)＝A _xx′δ(t-τ _xx′)

(5)

Wherein, A _{Xx '}Be separating of transport equation, τ _{Xx '}For between x and the x ' then, consider generalized case, suppose that in recording geometry the sampling function of focus and receiver is respectively h _s(s) and h _r(r),

Bring formula (5) into formula (4), can get:

P(r，t|s，0)＝∫m(x)A _sxA _xrW(t-τ _sx-τ _xr)dx

(6)

Consider to excite and receive sampling function h _s(s) and h _r(r), our equation that can obtain Kirchhoff skew can be expressed as:

m _k(x)＝∫h _s(s)ds∫h _r(r)dr×(∫P(r，t|s，0)A _sxA _xrW(t-τ _sx-τ _xr)dt) (7)

Bring formula (6) into formula (7), and the exchange integral order, the Kirchhoff offset equation can be written as:

m _k(x)＝∫m(x′)dx′∫h _s(s)A _sxA _sx′ds∫h _r(r)A _xrA _x′rdr×(∫W(t-τ _sx-τ _xr)W(t-τ _sx-τ _x′r)dt)

(8)

By formula (8) as can be known, integral kernel K (x, x ') is:

K(x，x′)＝∫h _s(s)ds∫h _r(r)drA _sxA _sx′A _xrA _x′r×R(τ _sx+τ _xr-τ _sx′-τ _x′r)?(9)

Wherein, R (τ _Sx+ τ _Xr-τ _{Sx '}-τ _{X ' r}) be source function W in difference simple crosscorrelation constantly, for this reason, (9) are brought into (8), (8) are reduced to:

m _k(x)＝∫K(x，x′)m(x′)dx′ (10)

(10) are write as discrete form, as shown in Equation (11):

$m_k=\underset{~}{K}m---\left(11\right)$

Wherein, by (3) as can be known, matrix

Physical meaning be the Hessian matrix, integral kernel K (x, x ') is equivalent to from m to m _kA mapping, for some scattering point x _f, K (x, x ') represents the system responses of this scattering point.If single track is offset, this response is the oval impulse response of skew; And a large amount of seismic traces is offset, and this response energy will converge to real scattering point position;

When the Hessian matrix was unit matrix, skew model mk equaled spherical model m truly, investigates the result of the diagonal entry of Hessian matrix, and by formula (9) as can be known, the diagonal entry of integral kernel matrix K (x, x ') can be written as following form,

$K(x_0,x_0)=R\left(0\right)\&Integral;h_s\left(s\right)\mathrm{ds}{A}_{{\mathrm{sx}}_0}^2\&Integral;h_r\left(r\right)\mathrm{dr}{A}_{x_{0}r}^2---\left(12\right)$

By formula (12) as can be known, the element value on the matrix principal diagonal is subjected to

And the noncontinuous sampling function (h of wave field _s(s), h _r(r)) control;

If distributing, focus and receiver caused the image data of non-rule, the Kirchhoff offset method will produce the false picture of skew on migrated section, in order to reduce the false picture of skew to the full extent, obtain comparatively real earth model m, we are given following objective function:

$P\left(m\right)={\left|\right|\underset{~}{L}m-P_0\left|\right|}^2+{\mathrm{\&epsiv;}}^2{\left|\right|\underset{~}{C}m-\underset{~}{C}{m}_{\mathrm{apr}}\left|\right|}^2---\left(13\right)$

Wherein, first on formula (13) equation the right is the error function item of data, and second is regularization term,

Represent on the model m of a linear matrix operator effect m _AprBe initial reference model, ε ²Be the weight that regularization is handled, equation (13) can be found the solution with the conjugated body gradient method of propositions in 1987 such as Tarantola, and the solving equation of m is:

$m={({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})}^{-1}({\underset{~}{L}}^T{P}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}})---\left(14\right)$

And equivalent equation, shown in equation (15):

$({\underset{~}{L}}^T\underset{~}{L}+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C})m={\underset{~}{L}}^T{P}_0+{\mathrm{\&epsiv;}}^2{\underset{~}{C}}^T\underset{~}{C}{m}_{\mathrm{apr}}---\left(15\right)$

Equation (14), the model m in (15) can find the solution with the method for iteration.

Images