CN102262245A - Self-adaptive weighted SIRT inversion method and system in seismic tomography processing - Google Patents

Abstract

The invention provides a self-adaptive weighted SIRT inversion method and a system thereof in seismic tomography processing. The inversion method and the inversion system add the grid slowness correction value iteration calculation formula

The weight correction part provides ray length-based correction for the slowness correction value of the grid, and solves the problem that the travel time errors caused by different ray lengths are different. By applying a weighting factor p_iThe weighted power exponent h in the calculation formula adopts a weighted power exponent h calculation formula based on ray density, so that the slowness correction value c of the grid_jDifferent degrees of weighted correction can be obtained according to different grid ray densities. Interpolation calculation modules and smoothing processing modules are added in the slowness value calculation units to perform interpolation calculation and smoothing processing on each grid slowness value result, so that the seismic tomography result is more reasonable.

Description

Adaptive weighted SIRT inversion method and system thereof during a kind of seismic tomography is handled

Technical field

The present invention relates to the geophysical exploration technology field, adaptive weighted SIRT during particularly a kind of seismic tomography is handled (while iterative reconstruction technology, simultaneous iterative reconstruction techniques) inversion method and system thereof.

Background technology

Ray tracing technique based on Fermat (Fermat) principle and Huygens-Fresnel (Huygens Fresnel) principle has a wide range of applications in the seismic tomography treatment technology.The positive algorithm of whilst on tour linear interpolation ray tracing that is based on ray theory that present technique adopts, by analysis, and then be finally inversed by a kind of geophysical method of the important informations such as structure, velocity distribution of the underground medium that covers by a large amount of rays to the seimic travel time data that observe.

Fig. 1 is seismic tomography process flow figure.As shown in the figure, this seismic tomography disposal route can be divided into following four steps: first arrival is picked up, ray tracing, equation and inverting were found the solution when foundation was walked.When at first, walking as actual shock wave by the first break time of picking up earthquake big gun collection record.Secondly, the ray tracing step utilizes ray-tracing scheme to find the solution ray tracing and whilst on tour that seismic event passes imaging region by under the condition of known media velocity distribution model.Here, described medium velocity distributed model generally is to treat the exploratory area to set up the grid cell model, supposes that the slowness in each grid is constant.Like this, on this net boundary arbitrarily the whilst on tour of any all can calculate by the whilst on tour linear interpolation of two discrete points that are adjacent on this net boundary.Then, the travel-time difference distance that data and forward simulation calculate when walking according to actual shock wave judges whether it satisfies certain condition, then sets up inverting equation when walking as not satisfying.The medium velocity distributed model of hypothesis before the equation correction when walking at last, according to the inverting of setting up.So circulate, finally obtain the structure and the velocity profile information of believable underground medium.

Wherein, described inverting solution procedure is when utilizing walking that known media velocity distribution model carries out that Theoretical Calculation obtains and the difference of real data when picking up walking of acquisition, and the medium velocity distributed model is revised.At present, prior art normally adopts while iterative reconstruction technology (being called for short SIRT) that the medium velocity distributed model is revised in this step.The treatment step of this method is as follows:

Suppose that the slowness in j the grid of medium velocity distributed model is s _j, j=1,2 ..., M (M is the model meshes sum);

Utilize known media velocity distribution model to carry out Theoretical Calculation, the T when theory that obtains every ray is walked _Ic:

$T_{\mathrm{ic}}=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}{s}_j,$ I=1,2 ..., N (N is the ray total number) (1)

G in the formula _IjFor the i bar ray that obtains by ray tracing j grid inner rays length.T when obtaining every actual shock wave of ray and walk according to this formula _ImT when just drilling theory and walk _IcDifference Δ t _i:

Δt _i＝T _im-T _ic，i＝1，2…，N (2)

If the slowness modified value c of i bar ray in j grid _Ij, then:

$\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}{c}_{\mathrm{ij}}=\mathrm{\Δ}{t}_i,$ i＝1，2…，N (3)

If modified value c _IjBe proportional to the path g that j grid inner rays passed through _IjWith the long R of this ray _iThe ratio, that is:

$c_{\mathrm{ij}}={\mathrm{\α}}_i\frac{g_{\mathrm{ij}}}{R_i},$ i＝1，2…，N (4)

α in the formula _iBe the proportionality constant of i bar ray, R _iBe the total length of i bar ray:

$R_i=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{g}_{\mathrm{ij}},$ i＝1，2，…，N (5)

With (4), (5) formula substitution (3) formula, and readjusting and simplifying, can get:

$c_{\mathrm{ij}}=\frac{\mathrm{\Δ}{t}_i{g}_{\mathrm{ij}}}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\left(g_{\mathrm{ij}}\right)}^2},$ i＝1，2，…，N (6)

If the ray total number through j grid is Y _j, then the average modified value of each grid is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{1}{Y_j}{c}_{\mathrm{ij}}---\left(7\right)$

With this average modified value

Slowness value s to j grid _jRevise, that is:

$s_j^{\left(q\right)}=s_j^{(q-1)}+c_j^q---\left(8\right)$

Yet data can be subjected to noise when walking owing to the shock wave of actual acquisition, and it is irrational to model parameter correction that above-mentioned existing SIRT algorithm adopts average weights to each bar ray simultaneously.Because can be subjected to more noise during than the walking of long shot thread path, so the modified value that equation produced of this ray correspondence should give less weights; Otherwise data can be reliable during the walking of short raypath correspondence, and should give bigger correction weights.Simultaneously, this revises weights also should be relevant with the density of ray, and the big more correction weights of radiographic density are more little, and the positive weights of light maintenance are big more more for radiographic density.So-called radiographic density is meant the number of lines of penetrating that passes in the model meshes.In addition, owing to the cautious quantity of big gun set in the seismic tomography process is very limited, having many grids does not have the ray process.These do not have the grid of ray process, and the grid slowness can not get revising in the process of inverting, can produce bigger adverse effect to the seismic tomography effect.

Based on above reason, we are necessary to design the adaptive weighted SIRT inversion method in a kind of seismic tomography processing, so that the problems referred to above are carried out the self-adaptation adjustment.

Summary of the invention

Fundamental purpose of the present invention is to solve problems of the prior art, provide a kind of seismic tomography adaptive weighted SIRT inversion method and the system thereof in handling, with to carrying out the self-adaptation adjustment because of caused each round-off error of noise, radiographic density and vacant grid.

The objective of the invention is to be achieved by following technical proposals:

Adaptive weighted SIRT inversion method during a kind of seismic tomography is handled is characterized in that, comprising:

Data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _Ij

Data T when the actual shock wave that is obtained being walked based on grid slowness modified value iterative computation formula _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}^2}\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; Q is the number of times of iterative computation;

Be the q time iteration slowness modified value of j grid; g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor; p _iBe the adaptive weighted factor, this adaptive weighted factor

Wherein, h is the weighting power exponent, and this weighting power exponent h is the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly $L_p=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{p}_i;$

Slowness modified value c according to the grid that calculates _jCalculate the slowness value s of this grid _j:

Described data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjStep comprises following concrete steps:

In actual earthquake big gun collection record, pick up first break time T when walking as actual shock wave _Im

By the ray tracing module known media velocity distribution model is carried out forward simulation, calculate the path g of i bar ray in j grid _Ij, just drilling theoretical shock wave data when walking

The equation when inverting of setting up ray is walked.

Described weighting power exponent h value is 1.

At described weighting factor p _iIn also introduced an equilibrium factor δ _i, i.e. this weighting factor p _iCalculating formula be $p_i={(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{g}_{\mathrm{ij}}+{\mathrm{\δ}}_i)}^{-h};$ This equilibrium factor δ _iIt is a constant.

Described equilibrium factor δ _iValue is 0.01.

Described weighting power exponent h calculates and gets based on radiographic density; This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\ρ}}_j-{\mathrm{\ρ}}_{\min }}{{\mathrm{\ρ}}_{\max }-{\mathrm{\ρ}}_{\min }}+1$

Wherein, ρ _jBe the radiographic density of j grid, promptly pass the number of lines of penetrating of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

Calculating the slowness value s of grid _jAfterwards, also include the step of carrying out interpolation calculation at unknown grid; Should carry out based on following calculating formula at the interpolation calculation of unknown grid:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_{\mathrm{<figure-callout\; id="2"\; label="j\; d\; kj"\; filenames="HDA0000068191830000031.png"\; state="\{\{state\}\}">j</figure-callout>}}}{d_{\mathrm{kj}}^{}}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]$

Wherein, s _kIt is the slowness value of k unknown grid; s _jBe slowness value through the revised j of slowness modified value known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

The known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

Calculating the slowness value s of grid _jAfterwards, also include the step of smoothing processing:

Each treats that it is the mean value of all the grid slowness values in certain smoothing processing scope at center that the slowness value of smoothing processing grid equals with this grid.

Described smoothing processing scope is for being the rectangular area of 3 * 3 grids at center with the grid for the treatment of smoothing processing.

Adaptive weighted SIRT inverting system during a kind of seismic tomography is handled is characterized in that, comprising: first arrival pickup unit, ray tracing unit, adaptive weighted inverting unit and slowness value computing unit;

Described first arrival pickup unit is when picking up that actual shock wave is walked in the earthquake big gun collection record;

Described ray tracing unit is in order to pass through the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjJust drilling theoretical shock wave data T when walking _Ic

Described adaptive weighted inverting unit, data T when walking in order to receive the actual shock wave that sends by the first arrival pickup unit _Im, and the path g of i bar ray in j grid that sends by the ray tracing unit _IjJust drilling theoretical shock wave data T when walking _Ic, and data T when walking according to the actual shock wave that is obtained _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; The number of times of q iterative computation;

Be the q time iteration slowness modified value of j grid; g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor; p _iBe the adaptive weighted factor, this adaptive weighted factor

Wherein, h is the weighting power exponent, and this weighting power exponent h is the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly $L_p=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{p}_i;$

Described slowness value computing unit is in order to receive each grid slowness modified value c that sends adaptive weighted inverting unit _j, and according to slowness modified value c _jCalculate the slowness value s of this grid _j:

Described weighting power exponent h value is 1.

At described weighting factor p _iIn also introduced an equilibrium factor δ _i, i.e. this weighting factor p _iCalculating formula be $p_i={(\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}+{\mathrm{\&delta;}}_i)}^{-h};$ This equilibrium factor δ _iIt is a constant.

Described equilibrium factor δ _iValue is 0.01.

Described weighting power exponent h calculates and gets based on radiographic density; This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\&rho;}}_j-{\mathrm{\&rho;}}_{\min }}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1$

Wherein, ρ _jIt is the radiographic density of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

Also be provided with the interpolation calculation module in the described slowness value computing unit;

Described interpolation calculation module is used to receive the revised slowness value of the slowness modified value s that known grid obtains through Inversion Calculation _j, and according to the revised slowness value of this known grid s _jSlowness value to unknown grid is carried out the interpolation calculation correction; This interpolation calculation formula is:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_{\mathrm{<figure-callout\; id="2"\; label="j\; d\; kj"\; filenames="HDA0000068191830000031.png"\; state="\{\{state\}\}">j</figure-callout>}}}{d_{\mathrm{kj}}^{}}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]$

Wherein, s _kIt is the slowness value of k unknown grid; s _jBe slowness value through the revised j of slowness modified value known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

The known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

Also be provided with the smoothing processing module in the described slowness value computing unit;

Described smoothing processing module is used to receive the slowness value s of grid _j, and according to treat that the smoothing processing grid is the slowness value after all grid slowness values in certain smoothing processing scope at center are calculated this and treated smoothing processing grid smoothing processing.

Described smoothing processing scope is for being the rectangular area of 3 * 3 grids at center with the grid for the treatment of smoothing processing.

The invention has the beneficial effects as follows:

(1) the present invention is by having added in grid slowness modified value iterative computation formula

The weight retouch, provide correction to the slowness modified value of grid based on ray length, solved because ray length different introduce walk the also different error problem of time error.

(2) pass through weighting factor p _iWeighting power exponent h in the calculating formula adopts the weighting power exponent h calculating formula based on radiographic density, makes the slowness modified value c of grid _jCan obtain weighting correction in various degree with the difference of grid radiographic density.

(3) by in slowness value computing unit, increasing interpolation calculation module and level and smooth processing module, each grid slowness value result is carried out interpolation calculation and smoothing processing processing, make that this seismic tomography result is more reasonable.

Description of drawings

Accompanying drawing described herein is used to provide further understanding of the present invention, constitutes the application's a part, does not constitute limitation of the invention.In the accompanying drawings:

Fig. 1 is seismic tomography process flow figure;

Fig. 2 is the process flow diagram of the adaptive weighted SIRT inversion method in the seismic tomography processing;

Fig. 3 is the system construction drawing of the adaptive weighted SIRT inverting system of seismic tomography in handling.

Embodiment

For making the purpose, technical solutions and advantages of the present invention clearer,, the present invention is described in further details below in conjunction with embodiment and accompanying drawing.At this, exemplary embodiment of the present invention and explanation thereof are used to explain the present invention, but not as a limitation of the invention.

As previously mentioned, data can be subjected to noise when walking owing to the shock wave of actual acquisition, and the SIRT inversion method employing during therefore existing seismic tomography is handled is irrational to the method that each bar ray is averaged weights to model parameter correction.Because can be subjected to more noise during than the walking of long shot thread path, so the modified value that equation produced of this ray correspondence should give less weights; Otherwise data can be reliable during the walking of short raypath correspondence, and should give bigger correction weights.Therefore, the SIRT inversion technique during the present invention handles existing seismic tomography is improved, and has added the adaptive weighted inversion algorithm relevant with ray length.

Fig. 2 is the process flow diagram of the adaptive weighted SIRT inversion method in the designed seismic tomography processing of the present invention.As shown in the figure, this adaptive weighted SIRT inversion method comprises:

Data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _Ij

Data T when the actual shock wave that is obtained being walked based on grid slowness modified value iterative computation formula _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)---\left(10\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; Q is the number of times of iterative computation;

Be the q time iteration slowness modified value of j grid, wherein g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor,

p _iBe the adaptive weighted factor, this adaptive weighted factor

Wherein, h is the weighting power exponent, and its value can be determined according to the weight size, obtains excessive weight for preventing extremely short ray, and generally this weighting power exponent h gets the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly

Slowness modified value c according to the grid that calculates _jCalculate the slowness value s of this grid _j:

Therefore, the adaptive weighted SIRT inversion method that the present invention is designed, its key is slowness modified value c in each grid _jCalculating on.Wherein, the computing formula of existing slowness modified value algorithm is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{1}{Y_j}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)$

On this basis, the present invention is based on the aforementioned reason of analyzing, in this formula, added

The weight retouch, iterative computation is partly revised.Because

Therefore this weight partly is equivalent to

Wherein, the adaptive weighted factor

Because, g _IjBe the length of i bar ray in j grid, therefore wherein

Be the total length of i bar ray.Therefore, this weight retouch

Correction based on ray length is provided, solved aforementioned point out because ray length different introduce walk the also different error problem of time error.

Wherein, we can pass through the numerical value of the weighting power exponent h of this adaptive weighted factor of adjusting, reach the purpose of regulating this weight retouch weighting size, so that its adjusting factor is tending towards reasonable.Generally, this weighting power exponent h should value be 1.

Described data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjStep comprises following concrete steps:

In actual earthquake big gun collection record, pick up first break time T when walking as actual shock wave _Im

By the ray tracing module known media velocity distribution model is carried out forward simulation, calculate the path g of i bar ray in j grid _Ij, just drilling theoretical shock wave data when walking

The equation when inverting of setting up ray is walked.

In addition, because the present invention's set weight retouch in grid slowness modified value iterative computation formula

In the adaptive weighted factor be

Based on this calculating formula, if the length of ray is too short, then its weighted value in formula is excessive, can influence final result of calculation.For avoiding the appearance of such situation, we are also at weighting factor p _iCalculating formula in introduced an equilibrium factor δ _i, be about to this weighting factor p _iCalculating formula change into $p_i={(\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}+{\mathrm{\&delta;}}_i)}^{-h}.$ This equilibrium factor δ _iBe a very little constant, given experience value is 0.01 in the present embodiment.By at weighting factor p _iCalculating formula in introduce such equilibrium factor δ _i, make the weight retouch that the present invention introduced

Can not produce the bigger error of calculation because of the too short ray of individual lengths.Certainly, the user can be according to actual needs to this equilibrium factor δ _iValue regulate, to correspond to actual needs.But no matter which kind of value the user adopts in actual applications, and it all should be considered as within protection scope of the present invention.

In addition, described weighting factor

In weighting power exponent h be that fixed value also is inappropriate.Because, the slowness modified value c of grid _jNot only relevant with the length of ray, also should be relevant with the density of ray.Here the radiographic density of indication is the bar number that passes ray in the grid.Radiographic density is big more, and the chance that this grid model parameter obtains revising is many more, and slowness correction weights should suitably reduce, otherwise radiographic density is more little, and slowness correction weights are big more.Therefore, we are to weighting factor p _iWeighting power exponent h in the calculating formula has designed a calculating formula based on radiographic density.This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\&rho;}}_j-{\mathrm{\&rho;}}_{\min }}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1---\left(11\right)$

Wherein, ρ _jBe the radiographic density of j grid, promptly pass the number of lines of penetrating of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

By adopting the calculating formula of above-mentioned weighting power exponent h based on radiographic density, weighting factor p _iWeighting power exponent h in the calculating formula no longer is a fixed value, but the variable quantity based on radiographic density, makes the slowness modified value c of grid _jCan obtain weighting correction in various degree with the difference of grid radiographic density.

Have, in order to improve the effect of chromatography static correction, the present invention is also calculating the slowness value s of grid again _jAfterwards, increased interpolation algorithm to velocity field.This is that the grid slowness can not get revising in aforesaid refutation process because the cautious quantity of the big gun of seismic tomography is very limited, and having many grids does not have the raypath process, and these do not have the grid of ray process.Test shows, if the slowness value of these grids is not handled, can produce bigger adverse effect to the effect of chromatography static correction.Therefore, we have designed the method for speed interpolation, by around have the speed interpolation of modified value grid not revised the slowness value of grid.

In the present embodiment, our the concrete interpolation algorithm that adopts is anti-square distance method, promptly the distance of unknown grid and known grid as weight factor, the distance of unknown grid and known grid is near more, its weight is big more, otherwise then more little, weight is provided by the inverse ratio of square distance, and concrete interpolation algorithm expression formula is:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_{\mathrm{<figure-callout\; id="2"\; label="j\; d\; kj"\; filenames="HDA0000068191830000031.png"\; state="\{\{state\}\}">j</figure-callout>}}}{d_{\mathrm{kj}}^{}}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]---\left(12\right)$

Wherein, s _kBe the slowness value of k unknown grid, promptly need to carry out the grid slowness of interpolation calculation; s _jSlowness value for the revised j of slowness modified value that obtains through Inversion Calculation known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

By above-mentioned interpolation algorithm, we can pass through the slowness value of the revised known grid of inversion algorithm according to the front, the correction that the slowness value of unknown grid is rationalized, thus obtain better chromatography static correction effect.

Here, the known grid distribution range that participates in interpolation calculation should be too not big, otherwise can increase search time, reduces counting yield, and limited to the influence of unknown grid from unknown grid known grid slowness value far away.For this reason, can set the known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

And, interpolation processing algorithm given here is not that any situation all is suitable for, if known grid is more than unknown number of grid, and unknown grid is to be distributed in discretely in the middle of the known grid, and then these unknown grids are higher through the result reliability that interpolation calculation obtains.If unknown grid is more than known grid, perhaps unknown grid be continuous distribution in a bulk of zone, then the result's that obtains through interpolation calculation of these unknown grids reliability is lower.

In addition,, more approach actual central velocity field characteristic distributions for the velocity field that makes renewal changes too fierceness, and the stability that guarantees inverting.The present invention is also calculating the slowness value s of grid _jAfterwards, increased the step of velocity field being made smoothing processing.Specifically the smoothing processing method that we taked is as follows:

Calculating the slowness value s of grid _jAfterwards, also comprise the smoothing processing step:

Each treats that it is the mean value of all the grid slowness values in certain smoothing processing scope at center that the slowness value of smoothing processing grid equals with this grid.

Medium velocity distributed model according to seismic tomography is generally divided with rectangular node, and therefore the smoothing processing scope of above-mentioned delineation also should be the rectangular area.And the zone of the smoothing processing scope of this delineation should not be too big, and too conference reduces precision, and speed of convergence slows down.The rectangular area is too little, level and smooth poor effect.In the present embodiment, actual what adopt is to be that the rectangular area of 3 * 3 grids at center is this smoothing processing scope with the grid for the treatment of smoothing processing.

Above-mentioned smoothing processing method not only can be arranged on the slowness value s that calculates grid by inversion algorithm _jAfterwards, equally also can be arranged on through after the above-mentioned interpolation algorithm processing.At this, such treatment step is not just being given unnecessary details.

Fig. 3 is the system construction drawing of the adaptive weighted SIRT inverting system of the designed seismic tomography of the present invention in handling.As shown in the figure, this adaptive weighted inverting system comprises: first arrival pickup unit 1, ray tracing unit 2, adaptive weighted inverting unit 3 and slowness value computing unit 4.

Described first arrival pickup unit 1, T when picking up that actual shock wave is walked in the earthquake big gun collection record _Im

Described ray tracing unit 2 is in order to pass through the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjJust drilling theoretical shock wave data T when walking _Ic

Described adaptive weighted inverting unit 3, data T when walking in order to receive the actual shock wave that sends by first arrival pickup unit 1 _Im, and the path g of i bar ray in j grid that sends by ray tracing unit 2 _IjJust drilling theoretical shock wave data T when walking _Ic, and data T when walking according to the actual shock wave that is obtained _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is: $c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)---\left(10\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; Q is the number of times of iterative computation;

Be the q time iteration slowness modified value of j grid, wherein

g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor, p _iBe the adaptive weighted factor, this adaptive weighted factor

Wherein, h is the weighting power exponent, and its value can be determined according to the weight size, obtains excessive weight for preventing extremely short ray, and generally this weighting power exponent h gets the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly

Described slowness value computing unit 4 is in order to receive each grid slowness modified value c that sends adaptive weighted inverting unit 3 _j, and according to slowness modified value c _jCalculate the slowness value s of this grid _j:

In the adaptive weighted SIRT inverting system of above-mentioned design, its key is by adaptive weighted inverting unit 3 slowness modified value c in each grid _jCalculating on.Wherein, the computing formula of existing slowness modified value algorithm is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{1}{Y_j}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)$

On this basis, the present invention is based on the aforementioned reason of analyzing, in this formula, added

The weight retouch, iterative computation is partly revised.Because

Therefore this weight partly is equivalent to Wherein, the adaptive weighted factor

Because, g _IjBe the length of i bar ray in j grid, therefore wherein

Be the total length of i bar ray.Therefore, this weight retouch

Correction based on ray length is provided, solved aforementioned point out because ray length different introduce walk the also different error problem of time error.

Wherein, we can pass through the numerical value of the weighting power exponent h of this adaptive weighted factor of adjusting, reach the purpose of regulating this weight retouch weighting size, so that its adjusting factor is tending towards reasonable.Generally, this weighting power exponent h should value be 1.

In addition, because the present invention's set weight retouch in grid slowness modified value iterative computation formula

In the adaptive weighted factor be

Based on this calculating formula, if the length of ray is too short, then its weighted value in formula is excessive, can influence final result of calculation.For avoiding the appearance of such situation, we are also at weighting factor p _iCalculating formula in introduced an equilibrium factor δ _i, be about to this weighting factor p _iCalculating formula change into $p_i={(\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}+{\mathrm{\&delta;}}_i)}^{-h}.$ This equilibrium factor δ _iBe a very little constant, given experience value is 0.01 in the present embodiment.By at weighting factor p _iCalculating formula in introduce such equilibrium factor δ _i, make the weight retouch that the present invention introduced

Can not produce the bigger error of calculation because of the too short ray of individual lengths.Certainly, the user can be according to actual needs to this equilibrium factor δ _iValue regulate, to correspond to actual needs.But no matter which kind of value the user adopts in actual applications, and it all should be considered as within protection scope of the present invention.

In addition, described weighting factor

In weighting power exponent h be that fixed value also is inappropriate.Because, the slowness modified value c of grid _jNot only relevant with the length of ray, also should be relevant with the density of ray.Here the radiographic density of indication is the bar number that passes ray in the grid.Radiographic density is big more, and the chance that this grid model parameter obtains revising is many more, and slowness correction weights should suitably reduce, otherwise radiographic density is more little, and slowness correction weights are big more.Therefore, we are to weighting factor p _iWeighting power exponent h in the calculating formula has designed a calculating formula based on radiographic density.This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\&rho;}}_j-{\mathrm{\&rho;}}_{\min }}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1---\left(11\right)$

Wherein, ρ _jBe the radiographic density of j grid, promptly pass the number of lines of penetrating of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

By adopting the calculating formula of above-mentioned weighting power exponent h based on radiographic density, weighting factor p _iWeighting power exponent h in the calculating formula no longer is a fixed value, but the variable quantity based on radiographic density, makes the slowness modified value c of grid _jCan obtain weighting correction in various degree with the difference of grid radiographic density.

Have, because the cautious quantity of the big gun of seismic tomography is very limited, having many grids does not have the raypath process again, and these do not have the grid of ray process, and the grid slowness can not get revising in aforesaid refutation process.Test shows, if the slowness value of these grids is not handled, can produce bigger adverse effect to the effect of chromatography static correction.Therefore, in order to improve the effect of chromatography static correction, as shown in Figure 3, the present invention also is provided with interpolation calculation module 41 in the slowness value computing unit 4 of above-mentioned adaptive weighted SIRT inverting system.

Described interpolation calculation module 41 is used to receive the revised slowness value of the slowness modified value s that known grid obtains through Inversion Calculation _j, and according to the revised slowness value of this known grid s _jSlowness value to unknown grid is carried out the interpolation calculation correction; This interpolation calculation formula is:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_{\mathrm{<figure-callout\; id="2"\; label="j\; d\; kj"\; filenames="HDA0000068191830000031.png"\; state="\{\{state\}\}">j</figure-callout>}}}{d_{\mathrm{kj}}^{}}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]---\left(12\right)$

Wherein, s _kBe the slowness value of k unknown grid, promptly need to carry out the grid slowness of interpolation calculation; s _jSlowness value for the revised j of slowness modified value that obtains through Inversion Calculation known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

By above-mentioned interpolation algorithm, we can pass through the slowness value of the revised known grid of inversion algorithm according to the front, the correction that the slowness value of unknown grid is rationalized, thus obtain better chromatography static correction effect.

Here, the known grid distribution range that participates in interpolation calculation should be too not big, otherwise can increase search time, reduces counting yield, and limited to the influence of unknown grid from unknown grid known grid slowness value far away.For this reason, can set the known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

In addition,, more approach actual central velocity field characteristic distributions for the velocity field that makes renewal changes too fierceness, and the stability that guarantees inverting.As shown in Figure 3, the present invention also is provided with smoothing processing module 42 in the slowness value computing unit 4 of above-mentioned adaptive weighted SIRT inverting system.

Described smoothing processing module 42 is used to receive the slowness value s of grid _j, and according to treat that the smoothing processing grid is the slowness value after all grid slowness values in certain smoothing processing scope at center are calculated this and treated smoothing processing grid smoothing processing.

Medium velocity distributed model according to seismic tomography is generally divided with rectangular node, and therefore the smoothing processing scope of above-mentioned delineation also should be the rectangular area.And the zone of the smoothing processing scope of this delineation should not be too big, and too conference reduces precision, and speed of convergence slows down.The rectangular area is too little, level and smooth poor effect.In the present embodiment, actual what adopt is to be that the rectangular area of 3 * 3 grids at center is this smoothing processing scope with the grid for the treatment of smoothing processing.

In sum, the present invention is by having added in grid slowness modified value iterative computation formula

The weight retouch, provide correction to the slowness modified value of grid based on ray length, solved because ray length different introduce walk the also different error problem of time error.Simultaneously, by to adaptive weighted factor p wherein _iA series of adjustment, make the correction result of this weight retouch better accurately.In addition, by increase interpolation calculation module and level and smooth processing module in slowness value computing unit, the result adjusts processing to each grid slowness value, makes that this seismic tomography result is more reasonable.Persons skilled in the art any not creative transformation of doing under this design philosophy all should be considered as within protection scope of the present invention.

Claims (19)

1. the adaptive weighted SIRT inversion method during a seismic tomography is handled is characterized in that, comprising:

Data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _Ij

Data T when the actual shock wave that is obtained being walked based on grid slowness modified value iterative computation formula _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; Q is the number of times of iterative computation;

Be the q time iteration slowness modified value of j grid; g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor; p _iBe the adaptive weighted factor, this adaptive weighted factor

Wherein, h is the weighting power exponent, and this weighting power exponent h is the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly $L_p=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{p}_i;$

Slowness modified value c according to the grid that calculates _jCalculate the slowness value s of this grid _j:

2. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 1 is handled is characterized in that: described data T when obtaining the cautious actual shock wave of each big gun that seismic tomography obtains in handling and walking _ImJust drilling theoretical shock wave data T when walking _IcAnd by the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjStep comprises following concrete steps:

In actual earthquake big gun collection record, pick up first break time T when walking as actual shock wave _Im

By the ray tracing module known media velocity distribution model is carried out forward simulation, calculate the path g of i bar ray in j grid _Ij, just drilling theoretical shock wave data when walking

The equation when inverting of setting up ray is walked.

3. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 1 is handled, it is characterized in that: described weighting power exponent h value is 1.

4. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 1 is handled is characterized in that: at described weighting factor p _iIn also introduced an equilibrium factor δ _i, i.e. this weighting factor p _iCalculating formula be $p_i={(\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}+{\mathrm{\&delta;}}_i)}^{-h};$ This equilibrium factor δ _iIt is a constant.

5. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 4 is handled is characterized in that: described equilibrium factor δ _iValue is 0.01.

6. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 1 is handled is characterized in that:

Described weighting power exponent h calculates and gets based on radiographic density; This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\&rho;}}_j-{\mathrm{\&rho;}}_{\min }}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1$

Wherein, ρ _jBe the radiographic density of j grid, promptly pass the number of lines of penetrating of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

7. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 1 is handled is characterized in that: calculating the slowness value s of grid _jAfterwards, also include the step of carrying out interpolation calculation at unknown grid; Should carry out based on following calculating formula at the interpolation calculation of unknown grid:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_j}{d_{\mathrm{kj}}^2}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]$

Wherein, s _kIt is the slowness value of k unknown grid; s _jBe slowness value through the revised j of slowness modified value known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

8. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 7 is handled, it is characterized in that: the known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

9. as the adaptive weighted SIRT inversion method in claim 1 or the 7 described seismic tomographies processing, it is characterized in that: calculating the slowness value s of grid _jAfterwards, also include the step of smoothing processing:

Each treats that it is the mean value of all the grid slowness values in certain smoothing processing scope at center that the slowness value of smoothing processing grid equals with this grid.

10. the adaptive weighted SIRT inversion method during seismic tomography as claimed in claim 9 is handled is characterized in that: described smoothing processing scope is for being the rectangular area of 3 * 3 grids at center with the grid for the treatment of smoothing processing.

11. the adaptive weighted SIRT inverting system during a seismic tomography is handled is characterized in that, comprising: first arrival pickup unit, ray tracing unit, adaptive weighted inverting unit and slowness value computing unit;

Described first arrival pickup unit is when picking up that actual shock wave is walked in the earthquake big gun collection record;

Described ray tracing unit is in order to pass through the path g of known media velocity distribution Model Calculation i bar ray in j grid _IjJust drilling theoretical shock wave data T when walking _Ic

Described adaptive weighted inverting unit, data T when walking in order to receive the actual shock wave that sends by the first arrival pickup unit _Im, and the path g of i bar ray in j grid that sends by the ray tracing unit _IjJust drilling theoretical shock wave data T when walking _Ic, and data T when walking according to the actual shock wave that is obtained _ImJust drilling theoretical shock wave data T when walking _IcCarry out iterative computation, obtain the slowness modified value c in each grid _jDescribed grid slowness modified value iterative computation formula is:

$c_j^{\left(q\right)}=c_j^{(q-1)}+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\left(\frac{p_i}{L_p}\frac{g_{\mathrm{ij}}{r_i}^{(q-1)}}{\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}^2}\right)$

Wherein, j is the model meshes sequence number; M is the model meshes sum; I is the sequence number of ray; N is the total number of ray; The number of times of q iterative computation;

Be the q time iteration slowness modified value of j grid; g _IjBe the length of i bar ray in j grid; r _iData T when walking for actual shock wave _ImJust drilling theoretical shock wave data T when walking _IcPoor; p _iBe the adaptive weighted factor, this adaptive weighted factor Wherein, h is the weighting power exponent, and this weighting power exponent h is the real number more than or equal to 1; L _pBe all adaptive weighted factor p _iAnd, promptly $L_p=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{p}_i;$

Described slowness value computing unit is in order to receive each grid slowness modified value c that sends adaptive weighted inverting unit _j, and according to slowness modified value c _jCalculate the slowness value s of this grid _j:

The adaptive weighted SIRT inverting system during 12. seismic tomography as claimed in claim 11 is handled, it is characterized in that: described weighting power exponent h value is 1.

13. the adaptive weighted SIRT inverting system during seismic tomography as claimed in claim 11 is handled is characterized in that: at described weighting factor p _iIn also introduced an equilibrium factor δ _i, i.e. this weighting factor p _iCalculating formula be $p_i={(\underset{j=1}{\overset{M}{\mathrm{\&Sigma;}}}{g}_{\mathrm{ij}}+{\mathrm{\&delta;}}_i)}^{-h};$ This equilibrium factor δ _iIt is a constant.

14. the adaptive weighted SIRT inverting system during seismic tomography as claimed in claim 13 is handled is characterized in that: described equilibrium factor δ _iValue is 0.01.

The adaptive weighted SIRT inverting system during 15. seismic tomography as claimed in claim 11 is handled, it is characterized in that: described weighting power exponent h calculates and gets based on radiographic density; This calculating formula based on the weighting power exponent h of radiographic density is:

$h=\frac{{\mathrm{\&rho;}}_j-{\mathrm{\&rho;}}_{\min }}{{\mathrm{\&rho;}}_{\max }-{\mathrm{\&rho;}}_{\min }}+1$

Wherein, ρ _jIt is the radiographic density of j grid; ρ _MaxAnd ρ _MinBe respectively the maximal value and the minimum value of radiographic density.

16. the adaptive weighted SIRT inverting system during seismic tomography as claimed in claim 11 is handled is characterized in that: also be provided with the interpolation calculation module in the described slowness value computing unit;

Described interpolation calculation module is used to receive the revised slowness value of the slowness modified value s that known grid obtains through Inversion Calculation _j, and according to the revised slowness value of this known grid s _jSlowness value to unknown grid is carried out the interpolation calculation correction; This interpolation calculation formula is:

$s_k=[\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{s_j}{d_{\mathrm{kj}}^2}/\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{d_{\mathrm{kj}}^2}]$

Wherein, s _kIt is the slowness value of k unknown grid; s _jBe slowness value through the revised j of slowness modified value known grid; d _KjBe the distance of k unknown grid to j known grid; M is the number of the known grid of k unknown grid interpolation calculation of participation.

The adaptive weighted SIRT inverting system during 17. seismic tomography as claimed in claim 16 is handled, it is characterized in that: the known grid that participates in unknown grid interpolation calculation is four the known grids nearest apart from this unknown grid.

18. the adaptive weighted SIRT inverting system in handling as claim 11 or 16 described seismic tomographies is characterized in that: also be provided with the smoothing processing module in the described slowness value computing unit;

Described smoothing processing module is used to receive the slowness value s of grid _j, and according to treat that the smoothing processing grid is the slowness value after all grid slowness values in certain smoothing processing scope at center are calculated this and treated smoothing processing grid smoothing processing.

19. the adaptive weighted SIRT inverting system during seismic tomography as claimed in claim 18 is handled is characterized in that: described smoothing processing scope is for being the rectangular area of 3 * 3 grids at center with the grid for the treatment of smoothing processing.

Images