CN101923175B - Method for directly generating angle gathers by using wave-equation migration - Google Patents

Abstract

The invention relates to a method for directly generating angle gather by using wave-equation migration, which is applied to seismic data processing in seismic exploration and is a method for forming the wave-equation pre-stack depth migration of the angle gather of seismic migration. The method directly utilizes shot-domain migration results for generating the angle gather imaged by an underground reflecting interface by calculating the incident angle of a seismic source outgoing wave at an underground imaging point, and compensates the change of the nonuniform coverage of an observing system for amplitude. The angle gather directly corresponds to the parameter changes of physical properties of target layers and provides accurate angular-domain seismic reflecting strength for pre-stack inversion. The method has the core that because the incident angle of a seismic source wave field is determined by directly using a single-route wave algorithm, the method is better matched with a wave-equation migration method while avoiding the instability of a ray method. Compared with migration calculation, because only the simple-harmonic wave field of the single-route wave is utilized for the calculation of the incident angle, the metering amount in the method can be ignored. The method can be applied to two-dimensional and three-dimensional wave-equation pre-stack depth migration treatment and has an important application value in the explorations of oil gas resources and mineral resources.

Description

Wave equation migration directly produces road, angle diversity method

Technical field

The invention belongs to reflected seismic information processing technology field in the seismic prospecting, relate to migration before stack imaging and prestack inversion technology category in the seismic data processing procedure, is a kind of method that generates road, the angle collection of migration before stack at wave equation pre stack depth migration.

Background technology

In the seismic data treatment scheme, the migration before stack imaging is the link of most critical, and its main task is with the actual position of the seismic reflection signals of ground table record " focusing " to underground generation reflection.To one of migration result more general application is with the migration result stack of different big guns or different offset distances, obtains stacked section, thereby obtains the image of underground structure; And being migration result by different big guns or different offset distances, more deep application obtains the reflection strength of subsurface reflection point, in view of the above the oily situation of Direct Recognition underground medium to the different angles incident wave.

Existing pre-stack time migration and integral method pre-stack depth migration can obtain the migration before stack result of branch offset distance easily, thereby simply obtain the reflection strength of subsurface reflection point, promptly obtain common reflection point (CRP) the road collection that changes with offset distance different offset distances incident field.The simple hour offset distance of velocity variations has better simply transformational relation with incident angle, but the corresponding relation of following offset distance of complicated speed structure and incident angle becomes complicated, transforms to incident angle by offset distance and has become very difficult.

When the overlying strata of underground structure velocity variations is violent, based on the pre-stack time migration of pyramid face (pyramid) whilst on tour hypothesis and based on progressive (ray) theoretical integral method pre-stack depth migration just can not correct response the communication process of seismic wave practically, therefore can not exactly seismic reflection signals " focusings " be arrived the actual position that underground generation is reflected.Like this, even obtain common reflection point (CRP) the road collection that changes with offset distance, because the identical reflection spot that does not come is actually nonsensical.

Wave equation pre stack depth migration is to describe the propagation of seismic wave process with wave theory more accurately, the ripple propagation phenomenon that the repeatedly arrival that the correct Simulation of Complex overlying strata of energy causes etc. are complicated and the diffusing phenomenon of seismic wave propagation, therefore can be expected to obtain the imaging amplitude of correct response underground structure reflection strength again to the correct imaging of complex structure.But mostly the wave equation pre stack depth migration of serving the complex structure imaging is what the big gun territory skew by minute gun realized; It is simple obtaining stacked section by big gun territory migration result, but the big gun territory migration result of minute gun can not obtain common reflection point (CRP) the road collection with the offset distance variation simply as integral method.Since what really react the variation of underground structures parameter is the reflection strength (Jiao Daoji) that changes with incident angle, people just wish and can directly obtain Jiao Daoji by big gun territory migration result.

Developed the method for two classes by big gun territory Wave equation depth migration method generation road, angle collection.The one, before using image-forming condition, the anti-pass wave field of the incident field of shot point and geophone station is carried out the part plan wavelength-division and separate, to the imaging respectively of every pair of plane wave, the 2nd, in the image-forming condition of big gun territory skew, introduce " local offset distance " parameter, form the multiplanar imaging body, obtain Jiao Daoji by beam steering equal angles transform method.First method needs that imaging point is one by one carried out the part plan wavelength-division and separates, and converts the imaging of single-point to an imaging array, and calculated amount was difficult to bear when three-D migration calculated; Although the second method calculated amount reduces to some extent than first kind, introduce poly-the increasing of memory requirements that " local offset distance " parameter still causes imaging to be calculated.Two kinds of methods can not consider that all shot point covers the inhomogeneous imaging amplitude error that causes, and the effective frequency upper limit of road, the angle collection that these two class methods are generated has been controlled in the spatial sampling of calculations of offset.

The idea of a nature is directly to calculate the incident angle of the outgoing wave of shot point to imaging point, directly obtains Jiao Daoji by partial stack by big gun territory migration result then.But the method for existing calculating incident angle is based on ray-tracing algorithm, and this high-frequency approximation method is responsive to rate pattern, be difficult to obtain correct incident angle under the complicated velocity variations situation, so this method runs into many difficulties in actual applications.

The Jiao Daoji that is generated migration before stack by wave equation pre stack depth migration is still the problem that needs are furtherd investigate, has only accurate acquisition angle Dao Ji, could bring into play the advantage of wave equation pre stack depth migration more all sidedly, thereby for Direct Recognition reservoir under the complicated overlying strata contain fluid and oily lays the foundation, reduce the exploration risk of oil and gas reservoir under the complex structure situation.

Summary of the invention

The objective of the invention is: the method that a kind of road, angle collection by the wave equation pre stack depth migration output migration before stack of big gun territory is provided, the Jiao Daoji that reservoir physical parameter changes can be under the situation of complicated overlying strata, directly reacted, and the change of the non-homogeneous covering of recording geometry can be compensated amplitude; Road, the angle collection of output provides the seismic reflection intensity that changes with incident angle accurately for the prestack inversion technology.Compare with the big gun territory wave equation pre stack depth migration method of routine, this method does not increase too many calculated amount when forming road, angle collection; It has overcome existing all kinds of wave equation migrations and has generated calculated amount and the huge dual problem of memory requirements that road, angle collection technology exists.This method generates the wave equation migration of 3-D seismics data road, angle collection and is advanced to the industrial applications stage based on the inverting of road, angle collection.

The technical solution used in the present invention is: Wave equation depth migration directly produces road, angle diversity method, and concrete steps comprise:

(1) with monolateral or bilateral acquisition mode, obtains the seismic reflection signals that man-made explosion excites, record on the tape by wave detector;

(2) read seismic reflection signals from tape, be stored in the computing machine by the ordering of prestack big gun collection;

(3) utilize the one way ripple algorithm of wave equation migration respectively the source wavefield of shot point and the reception wave field of geophone station to be carried out the continuation of the wave field degree of depth along degree of depth axle, shot point after the continuation and geophone station wave field application deconvolution image-forming condition are carried out imaging, specifically: in the territory, frequency space, remove the geophone station wave field with the shot point wave field after the continuation, and with the real part weighted accumulation of the multiple wave field value of different frequency, division calculation adopts can guarantee the algorithm that result of calculation is stable;

(4) utilize the one way wave method to calculate the incident angle of the outgoing wave of shot point to underground each imaging point, specifically: the dominant frequency component to shot point carries out degree of depth continuation, with two kinds of algorithms the continuation result is calculated two spatial field, calculate the merchant of two fields and carry out medium filtering, level and smooth and normalization, try to achieve the incident angle of shot point outgoing wave underground each imaging point by arc cosine;

(5) according to the incident angle of each the shot point outgoing wave that obtains, will add up and compensate the influence of non-homogeneous covering by angle by each single big gun migration result that the deconvolution image-forming condition obtains, form angle Dao Ji amplitude;

(6) diagonal angle Dao Ji does the moving school of residue, and the corresponding numerical value that obvious stretched portion occurs is changed to zero, suppresses random noise with singular value decomposition method, promptly obtains the Jiao Daoji of the migration before stack of reaction medium surface physics parameter variation.

(7) read the variation of imaging amplitude in different spatial place with angle by the road, angle is concentrated, can be obtained the petrophysical parameter on stratum by inversion method, the petrophysical parameter image can be indicated the oily situation of underground structure.

Describedly shot point after the continuation and geophone station wave field are used the deconvolution image-forming condition carry out imaging and be achieved in that and make T ₀Be the duration of seismologic record, frequency sampling Δ ω=2 π/T then ₀(radian per second); If Δ z is the spatial sampling of depth direction, n is the Ser.No. of degree of depth up-sampling point, and x and y are the horizontal level coordinates, and m is the Ser.No. of frequency sampling point; If shot point and geophone station wave field are through being respectively P after the degree of depth continuation on the n Δ z degree of depth _D(x, y, n Δ z, m Δ ω) and P _U(x, y, n Δ z, m Δ ω) calculates

A(x，y，nΔz，mΔω)＝P _D(x，y，nΔz，mΔω){P _D(x，y，nΔz，mΔω)} ^*

Subscript * is total to volume, shot point and geophone station wave field P again in the formula _DAnd P _UUnit can be speed or pressure according to the wave detector type, the unit of frequency is a radian per second, the unit of the degree of depth and horizontal coordinate be rice.Make A ₀(n Δ z) is the mean value of A (x, y, n Δ z, m Δ ω) on degree of depth n Δ z, as A (x, y, n Δ z, m Δ ω) 〉=1.2A ₀(n Δ z) calculates

$B(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})=\frac{P_U(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω}){\left\{P_D(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})\right\}}^{*}}{A(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})}$

As A (x, y, n Δ z, m Δ ω)＜1.2A ₀(n Δ z) calculates

$B(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})=\frac{P_U(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω}){\left\{P_D(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})\right\}}^{*}}{A_0\left(\mathrm{n\Δz}\right)}(2-a)(1+{(1-a)}^2),$ Wherein $a=\frac{A(x,y,\mathrm{n\Δz},\mathrm{m\Δ\ω})}{A_0\left(\mathrm{n\Δz}\right)}$

The upper and lower limit that makes the seismic data effective band is respectively l ₁Δ ω and l ₂Δ ω, l here ₁And l ₂Be positive integer, then the imaging amplitude φ of position, space deconvolution image-forming condition (x, y, n Δ z) is

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\underset{m=l_1}{\overset{{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="DEST\_PATH\_HSB00000018701500031.png,DEST\_PATH\_HSB00000018701500051.png"\; state="\{\{state\}\}">l</figure-callout>}}_2}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{\mathrm{\&beta;B}(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})\right\}$

If j is a unit imaginary number, β=jm Δ ω when doing three-D migration calculating in the following formula is when two-dimensional migration calculates $\mathrm{\&beta;}=\sqrt{\mathrm{jm\&Delta;\&omega;}};$ The real part of multiple wave field value is asked in symbol Re representative.

Describedly utilize outgoing wave that the one way wave method calculates shot point that the incident angle of underground each imaging point is achieved in that to make that x and y are horizontal coordinates, z is a depth coordinate, and the dominant frequency of focus is ω _p(radian per second), the locus of focus are (x _s, y _s, 0), the speed of medium be c (x, y, z); Utilize the one way ripple algorithm of wave equation migration, in frequencies omega _pTo space impulse δ (x-x _s, y-y _s) carry out degree of depth continuation, obtaining with Δ z is multiple wave field value P (x, y, n Δ z, ω on each degree of depth discrete point of spacing _p, x _s, y _s); Calculate respectively

I ₀(x，y，nΔz)＝|Re{P(x，y，nΔz，ω _p，x _s，y _s)}|

$I_1(x,y,\mathrm{n\&Delta;z})=\mathrm{Re}\left\{j\frac{c(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&omega;}}_p}P(x,y,\mathrm{n\&Delta;z},{\mathrm{\&omega;}}_p,x_s,y_s)\right\}$

$I_2(x,y,\mathrm{n\&Delta;z})=\left|\frac{2}{3\mathrm{\&Delta;z}}\right[I_1(x,y,\mathrm{n\&Delta;z}+\mathrm{\&Delta;z})-I_1(x,y,\mathrm{n\&Delta;z}-\mathrm{\&Delta;z})]$

$-\frac{1}{12\mathrm{\&Delta;z}}[I_1(x,y,\mathrm{n\&Delta;z}+2\mathrm{\&Delta;z})-I_1(x,y,\mathrm{n\&Delta;z}-2\mathrm{\&Delta;z})]|$

X in the following formula _s, y _sAnd ω _pBe parameter, so in the function of left end, it is ignored; Symbol || it is the Ser.No. of degree of depth up-sampling point that absolute value, n are asked in representative; (unit z) is a meter per second to speed c for x, y, and the unit of frequency is a radian per second, and the unit of the degree of depth and horizontal coordinate is a rice, and the unit of space impulse is a pressure.If I ₂(x, y, n Δ z)＞1.05I ₀(x, y, n Δ z) then makes I ₂(x, y, n Δ z)=0.Make φ ₀(n Δ z) is I ₀(x, y, n Δ z) mean value on degree of depth n Δ z is worked as I ₀(x, y, n Δ z) 〉=1.2 φ ₀(n Δ z) calculates

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\frac{I_2(x,y,\mathrm{n\&Delta;z})}{I_0(x,y,\mathrm{n\&Delta;z})}$

Work as I ₀(x, y, n Δ z)＜1.2 φ ₀(n Δ z) calculates

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\frac{I_2(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&phi;}}_0\left(\mathrm{n\&Delta;z}\right)}(2-a)(1+{(1-a)}^2),$ Wherein $a=\frac{I_0(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&phi;}}_0\left(\mathrm{n\&Delta;z}\right)}$ φ (x, y, n Δ z) value to the space each point is carried out following processing:

(1) carries out medium filtering along depth direction;

(2) it is level and smooth to carry out space three-dimensional;

(3) to each degree of depth n Δ z, choose maximal value,, realize the maximal value normalization on this degree of depth with the numerical value of this value except that the every bit of this degree of depth.

Making the value after the above-mentioned processing is φ (x, y, n Δ z), then can get focus to the incident angle θ (unit is degree) of underground each imaging point to be:

θ＝acos(φ(x，y，nΔz))

(x, y, n Δ z) promptly is the volume coordinate of imaging point in the formula.。

The incident angle of the outgoing wave of each shot point that described foundation obtains, to add up by angle by each single big gun migration result that the deconvolution image-forming condition obtains and compensate the influence of non-homogeneous covering amplitude, form angle Dao Ji and be achieved in that definition angle sampling interval Δ θ (degree), to one-dimension array of each imaging point definition, array length is less than 90/ θ Δ; To whole big gun circulations, single big gun migration result to every big gun, read the incident angle θ of the outgoing wave of this shot point that obtains by step (4) to each imaging point, calculate l=1+int{ θ/Δ θ }, the amplitude of this imaging point is added on l the element of corresponding array, and statistics is added to the big gun number on this element simultaneously; Finish the circulation of whole big guns, the imaging amplitude array of the different depth point of same horizontal level is formed two bit array by the degree of depth, angle, promptly obtain the Jiao Daoji that compensates without non-homogeneous covering; The big gun number that adds up that obtains with statistics is removed the numerical value of the different angles of each imaging point, just obtains finishing the Jiao Daoji of non-homogeneous covering compensation.

Wave equation depth migration of the present invention directly produces road, angle diversity method, can obtain the Jiao Daoji of the migration before stack of subsurface reflection point through the part stack directly by the big gun territory migration result of minute gun.

Wave equation depth migration of the present invention directly produces road, angle diversity method, can reject the influence of complicated overlying strata, is directly reacted the migration before stack angle Dao Ji that the target zone physical parameter changes.

Wave equation depth migration of the present invention directly produces road, angle diversity method, can compensate the change of the non-homogeneous covering of recording geometry to amplitude.

Wave equation depth migration of the present invention directly produces road, angle diversity method, compares with the big gun territory wave equation pre stack depth migration method of routine, does not increase too many calculated amount.

Wave equation depth migration of the present invention directly produces road, angle diversity method, can be applicable to the two and three dimensions wave equation pre stack depth migration and handles.

Wave equation depth migration of the present invention directly produces road, angle diversity method, can be used for the petrophysical parameter on inverting stratum, directly obtains the oily situation of underground structure.

Specific implementation principle of the present invention is as follows:

Core concept of the present invention is directly to calculate the incident angle of the outgoing wave of shot point at underground imaging point with the one way wave method, has avoided rays method that the instability of rate pattern is also more mated with the wave equation migration method.Incident angle calculates the simple harmonic quantity wave field that only utilizes the one way ripple, compares calculated amount with migration algorithm and can ignore.

The principle of calculating incident angle is, by one way wave equation as can be known, and descending frequency field wave field

Satisfy

$\frac{\&PartialD;{\tilde{P}}^{+}(x,y,z;\mathrm{\&omega;})}{\&PartialD;z}=-j{k}_z{\tilde{P}}^{+}(x,y,z;\mathrm{\&omega;})---\left(1\right)$

Wherein j is a unit imaginary number, k _zBe vertical wave number, ω is an angular frequency, and x and y are the horizontal level coordinates, and z is a depth coordinate, is got by dispersion relation

$k_z=\frac{\mathrm{\&omega;}}{c(x,y,z)}\cos \mathrm{\&theta;}---\left(2\right)$

(x, y z) are the medium velocity of wave of position, space to c in the formula, and angle θ promptly is that (x, y z) locate the angle of wavefront surface normal and degree of depth axle z in the locus.With (2) formula substitution (1) formula and do discrete time domain Fu Shi inverse transformation, and make time t=0, can get

$\frac{\&PartialD;\left[\underset{\mathrm{\&omega;}>0}{\mathrm{\&Sigma;}}\mathrm{Re}\left(j\frac{c(x,y,z)}{\mathrm{\&omega;}}{\tilde{P}}^{+}(x,y,z,\mathrm{\&omega;})\right)\right]}{\&PartialD;z}=\cos \mathrm{\&theta;}\left[\underset{\mathrm{\&omega;}>0}{\mathrm{\&Sigma;}}\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,\mathrm{\&omega;})\right)\right]---\left(3\right)$ The real part of multiple wave field value is asked in symbol Re representative in the formula.If order ${\tilde{P}}^{+}(x,y,z=0,\mathrm{\&omega;})=\mathrm{\&delta;}(x-x_s,y-y_s)f\left(\mathrm{\&omega;}\right),$ Wherein f (ω) is the spectrum of seismic wavelet, then in (3) formula

Place (x corresponding to focal point _s, y _s, 0) and impulse response when locating, and the wavefront surface normal is exactly (x _s, y _s, 0) and locate the outgoing wave of shot point in (x, y, incident wave direction z).

Impulse response

In the space, only there is a wavefront surface,, can not determines shot point incident wave line of propagation with wavefront surface to there not being the locus of wavefront surface.If only the dominant frequency component is calculated impulse response, promptly only calculates Re

Then wavefront surface will be distributed in whole zone, be easy to try to achieve the incident wave direction of all imaging points, and calculate the also significantly minimizing of calculated amount of impulse response this moment.Only consider the dominant frequency component, can try to achieve by (3) formula

$\cos \mathrm{\&theta;}=\frac{\&PartialD;\left[\mathrm{Re}\left(j\frac{c(x,y,z)}{{\mathrm{\&omega;}}_p}{\tilde{P}}^{+}(x,y,z,{\mathrm{\&omega;}}_p)\right)\right]/\&PartialD;z}{\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,{\mathrm{\&omega;}}_p)\right)}---\left(4\right)$

Will there be a plurality of zero points in denominator in the formula (4), locate this zero point to obtain correct incident wave angle; Since the angle of incident wave does not have remarkable conversion in a wavelength coverage, the incident wave angle interpolation at the wave crest of available wavefront surface and trough place obtains the incident wave angle of close position.We the design following smoothing algorithm, realize this thought.At first calculate

$I_0(x,y,z)=\left|\mathrm{Re}\left({\tilde{P}}^{+}(x,y,z,{\mathrm{\&omega;}}_p)\right)\right|---\left(5\right)$

$I_1(x,y,z)=\mathrm{Re}\left(j\frac{c(x,y,z)}{{\mathrm{\&omega;}}_p}{\tilde{P}}^{+}(x,y,z,{\mathrm{\&omega;}}_p)\right)---\left(6\right)$

With the first order derivative in quadravalence difference algorithm calculating (4) formula, the spatial sampling that the depth direction of impulse response is calculated in order is Δ z, can get

$I_2(x,y,z)=\left|\frac{2}{3\mathrm{\&Delta;z}}\right[I_1(x,y,z+\mathrm{\&Delta;z})-I_1(x,y,z-\mathrm{\&Delta;z})]---\left(7\right)$

$-\frac{1}{12\mathrm{\&Delta;z}}[I_1(x,y,z+2\mathrm{\&Delta;z})-I_1(x,y,z-2\mathrm{\&Delta;z})]$

Adopting absolute value in formula (5) and (7), is because for descending ripple cos θ always just.By (4) formula as can be known, I ₂(x, y z) should be greater than I ₀(x, y, z).But when the skew noise is strong (in the more weak shadow region of incident wave), do not satisfy this condition sometimes; This shows that there is not incident wave accurately in this place, therefore can not obtain correct incident angle.In calculating, if I ₂(x, y, z)＞1.05I ₀(x, y z), can make I ₂(x, y z)=0, have avoided the difficulty of shadow region angle calculation like this.

For avoiding the instability of division calculation in the formula (4), we have developed the high-precision stable algorithm that is divided by, and the core of this algorithm is to utilize the series expansion approximate expression:

$\frac{1}{x}=\frac{1}{1-(1-x)}\&ap;1+(1-x)+{(1-x)}^2+{(1-x)}^3=(2-x)(1+{(1-x)}^2)---\left(8\right)$

When x goes to zero, formula (8) still can obtain stable result, and works as | and during x|≤1.2, formula (8) approximate has excellent precision.Utilize formula (8), we have developed the stable algorithm that is divided by,

$\mathrm{\&phi;}(x,y,z)=\frac{I_2(x,y,z)}{I_0(x,y,z)}=\left\{\begin{array}{c}\begin{array}{cc}\frac{I_2(x,y,z)}{I_0(x,y,z)},& I_0(x,y,z)\&GreaterEqual;1.2{\mathrm{\&phi;}}_0\left(z\right)\end{array}\\ \frac{I_2(x,y,z)}{{\mathrm{\&phi;}}_0\left(z\right)}(2-\frac{I_0(x,y,z)}{{\mathrm{\&phi;}}_0\left(z\right)})(1+{(1-\frac{I_0(x,y,z)}{{\mathrm{\&phi;}}_0\left(z\right)})}^2)\end{array}\right.---\left(9\right)$

φ in the formula ₀(z) be I ₀(x, y, z) mean value on depth z.

To the result of (9) formula, do medium filtering earlier along depth direction, it is level and smooth to carry out space three-dimensional then, can be obtained the incident angle of close position by the incident wave angle interpolation at wave crest and trough place; Be to avoid (7) formula difference etc. to calculate the error that produces, will (x, y z) carry out normalized, guarantee the existence of vertical incidence to the φ after the smoothing processing on each degree of depth.So just can pass through the degree of depth continuation of the simple harmonic quantity wave field of single frequency, try to achieve the incident angle of shot point outgoing wave underground each point

θ＝acos(φ(x，y，z))(10)

(x, y z) promptly are the volume coordinate of underground each point in the formula.

The present invention adopts the deconvolution image-forming condition to the skew of big gun territory, rather than adopts the dependent imaging condition; The deconvolution image-forming condition can compensate the diffusional effect of seismic event, obtains correct amplitude.

We utilize formula (8) to solve the stable problem of division in the deconvolution image-forming condition once more.The deconvolution image-forming condition of big gun territory skew can be expressed as

$\mathrm{\&phi;}(x,y,z)=\underset{m=l_1}{\overset{{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="DEST\_PATH\_HSB00000018701500031.png,DEST\_PATH\_HSB00000018701500051.png"\; state="\{\{state\}\}">l</figure-callout>}}_2}{\mathrm{\&Sigma;}}}\mathrm{Re}\left(\mathrm{\&beta;}\frac{P_U(x,y,z,\mathrm{m\&Delta;\&omega;})}{P_D(x,y,z,\mathrm{m\&Delta;\&omega;})}\right)---\left(11\right)$

The real part of multiple wave field value is asked in symbol Re representative in the formula, and Δ ω is a frequency sampling, and m is the Ser.No. of frequency sampling point, P _U(x, y, z, m Δ ω) and P _D(x, y, z, m Δ ω) is respectively the shot point and the geophone station wave field in depth z upper frequency territory after the degree of depth continuation; l ₁And l ₂Be positive integer, l ₁Δ ω and l ₂Δ ω is the upper and lower limit of corresponding seismic data effective band respectively; It is consistent with the geophone station wave field for the waveform that is held in picture introducing parameter beta, if j is a unit imaginary number, β=jm Δ ω when doing three-D migration calculating is when two-dimensional migration calculates $\mathrm{\&beta;}=\sqrt{\mathrm{jm\&Delta;\&omega;}}.$

Division in the formula (11) can further be expressed as

$\frac{P_U(x,y,z,\mathrm{m\&Delta;\&omega;})}{P_D(x,y,z,\mathrm{m\&Delta;\&omega;})}=\frac{P_U(x,y,z,\mathrm{m\&Delta;\&omega;}){\left\{P_D(x,y,z,\mathrm{m\&Delta;\&omega;})\right\}}^{*}}{P_D(x,y,z,\mathrm{m\&Delta;\&omega;}){\left\{P_D(x,y,z,\mathrm{m\&Delta;\&omega;})\right\}}^{*}}---\left(12\right)$

Subscript * representative is total to volume again in the formula.Making the denominator of real number in the formula (12) is A (x, y, z, m Δ ω), by (8) Shi Kede

$\frac{1}{A(x,y,z,\mathrm{m\&omega;})}=\left\{\begin{array}{c}\begin{array}{cc}\frac{1}{A(x,y,z,\mathrm{m\&omega;})},& A(x,y,z,\mathrm{m\&omega;})\&GreaterEqual;1.2A_0\left(z\right)\end{array}\\ \frac{1}{A_0\left(z\right)}(2-\frac{A(x,y,z,\mathrm{m\&omega;})}{A_0\left(z\right)})(1+{(1-\frac{A(x,y,z,\mathrm{m\&omega;})}{A_0\left(z\right)})}^2)\end{array}\right.---\left(13\right)$

A in the formula ₀(z) be the mean value of A (x, y, z, m ω) on depth z.With formula (13) substitution (12) substitution again (11), deconvolution image-forming condition that can be stable.

The skew amplitude that the deconvolution image-forming condition obtains has been reacted the pairing correct reflection strength of incident angle of imaging point place shot point wave field.When forming road, angle collection, although the sampling of angle is identical, the big gun number that adds up on the different angles is different, and angle is big more, and the quantity that adds up is many more; And the shot point covering is close more, and the quantity that adds up is also just many more.For obtaining correctly to be reacted into the imaging amplitude of picture point physical parameter contrast, we add up the big gun number that adds up in each angular interval in cumulative process, remove accumulation result with this quantity.So promptly considered imaging results angle domain distribute inhomogeneous, also compensated shot point and covered uneven influence.

Beneficial effect of the present invention: this method can be applicable to the seismic migration imaging of oil and gas reservoir under the complicated overlying strata, can generate the directly Jiao Daoji of the migration before stack of reaction zone of interest physical parameter variation.Oil gas and the direct detection techniques of fluid (having avoided various compensation based on well-log information) such as prestack inversion are served in this road, angle energy collecting better, have reduced the risk of no well or few well area oil-gas exploration.This method is advanced to lithology identification with wave equation pre stack depth migration by structure imaging, has expanded and promoted the range of application and the effect of wave equation pre stack depth migration.This method has significant application value to oil gas, the mineral resources exploration in a large amount of zones of having complicated superstratum.

Description of drawings

Fig. 1 is the typical single shot record of grand celebration LMD regional development earthquake, and collection is 16 lines, bilateral observation, track pitch 40m; This adopts wave equation migration directly to generate road, angle diversity method and has used 3600 big guns altogether, completely covers 11.2 square kilometres.

Fig. 2 a is the two dimension slicing along survey line 322 of the three-dimensional velocity model of pre-stack depth migration use, and Fig. 2 b was the two-dimension speed section of the vertical line direction in edge of CDP400.

Fig. 3 is the common image gather that wire size 322 and CDP number 400 place are made of all single big gun migration result, and deconvolution image-forming condition of the present invention has been adopted in single big gun skew.

Fig. 4 a and Fig. 4 b are the incident angle distribution plan of shot point at the shot point wave field at wire size 318 and CDP number 415 places, and incident angle is calculated and adopted new method of the present invention; Fig. 4 a is the two dimension angular section along line direction, and Fig. 4 b is the two dimension angular section along vertical line direction.

Fig. 5 is the Jiao Daoji that does not carry out non-homogeneous covering compensation at wire size 318 and CDP number 415 places, and this collects together is according to the incident angle of each shot point at wire size 318 and CDP number 415 places, does local stack by each single big gun migration result by angle and obtains.

Fig. 6 uses result behind the nonuniform compensation to road, Fig. 5 angle collection, comparison diagram 5 and Fig. 6 as can be known, low-angle amplitude is enhanced.

Fig. 7 at first uses the moving school of residue to Fig. 6 angle Dao Ji, and the corresponding numerical value that will occur obvious stretched portion again is changed to zero, then the final angle Dao Ji that obtains with singular value decomposition method compacting random noise.

Fig. 8 a is the partial enlarged drawing of road, Fig. 7 angle collection, and Fig. 8 b is that the amplitude that changes with angle at place, two typical layers positions contrasts.The partial enlarged drawing correspondence of Fig. 8 a the ground layer segment of well-log information arranged.Solid line in the comparison diagram of Fig. 8 b is the variation with angle of the degree of depth 970m that read by Fig. 7 and 1030m place amplitude, dotted line is to locate P ripple, S wave velocity and the density that sound wave and density logging data read according to this position (wire size 318 and CDP number 415), the reflection strength that changes with angle that is obtained by Theoretical Calculation.As seen both matches are fine in contrast from Fig. 8 b, and the Jiao Daoji that this surperficial the present invention tries to achieve does not need " correction " incident angle and reflection magnitude relation (AVA) feature of match reflecting interface preferably.

Fig. 9 a is the Jiao Daoji at the well location place of wire size 323 and CDP numbers 360, and Fig. 9 b is that amplify the part of road, angle collection, and Fig. 9 c is that two typical layers positions (degree of depth 840m and 975m) are located amplitude and changed contrast with theoretical AVA feature with angle.The meaning of each figure is identical with Fig. 8 b with Fig. 7, Fig. 8 a.

Figure 10 a is the Jiao Daoji at the well location place of wire size 363 and CDP numbers 448, and Figure 10 b is that amplify the part of road, angle collection, and Figure 10 c is that two typical layers positions (degree of depth 995m and 1095m) are located amplitude and changed contrast with theoretical AVA feature with angle.The meaning of each figure is identical with Fig. 8 b with Fig. 7, Fig. 8 a.

Figure 11 a and Figure 11 b are the grand celebration LMD block skew stacked sections that this migration processing obtains.Figure 11 a is the two dimension slicing along survey line 322, and Figure 11 b was the two dimension slicing of the vertical line direction in edge of CDP400.Figure 11 a and Figure 11 b are obtained along the angle direction stack by road, the angle collection of each imaging point.

Embodiment

Embodiment 1: Wave equation depth migration directly produces road, angle diversity method, is example at grand celebration LMD area, is specially following steps:

(1) use bilateral acquisition mode, obtain the seismic reflection signals that man-made explosion excites by wave detector, record on the tape, every big gun is laid 16 and is gathered line, and the track pitch of gathering on the line is 40m, and every big gun contains 16 * 168=2688 road altogether; Excite and write down 3600 big guns altogether, Fig. 1 is typical single shot record.

(2) read seismic reflection signals from tape, be stored in the PC-Cluster computing machine by the ordering of prestack big gun collection, calculations of offset will circulate to big gun and carry out by reading each single shot record.

(3) according to the migration velocity model of Fig. 2 a and 2b, utilize the one way ripple algorithm of wave equation migration respectively the source wavefield (being the room and time impulse function) and the big gun record (seismic signal that receives) of shot point to be carried out wave field extrapolation along degree of depth axle, shot point after the continuation and geophone station wave field application deconvolution image-forming condition are carried out imaging.

Make T ₀Be the duration of seismologic record, frequency sampling Δ ω=2 π/T then ₀(radian per second); If Δ z is the spatial sampling of depth direction, n is the Ser.No. of degree of depth up-sampling point, and x and y are the horizontal level coordinates, and m is the Ser.No. of frequency sampling point; If shot point and geophone station wave field are through being respectively P after the degree of depth continuation on the n Δ z degree of depth _D(x, y, n Δ z, m Δ ω) and P _U(x, y, n Δ z, m Δ ω) calculates

A (x, y, n Δ z, m Δ ω)=P _D(x, y, n Δ z, m Δ ω) { P _D(x, y, n Δ z, m Δ ω) } ^*Subscript * is total to volume, shot point and geophone station wave field P again in the formula _DAnd P _UUnit can be speed or pressure according to the wave detector type, the unit of frequency is a radian per second, the unit of the degree of depth and horizontal coordinate be rice.Make A ₀(n Δ z) is the mean value of A (x, y, n Δ z, m Δ ω) on degree of depth n Δ z, as A (x, y, n Δ z, m Δ ω) 〉=1.2A ₀(n Δ z) calculates

$B(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})=\frac{P_U(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;}){\left\{P_D(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})\right\}}^{*}}{A(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})}$

As A (x, y, n Δ z, m Δ ω)＜1.2A ₀(n Δ z) calculates

$B(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})=\frac{P_U(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;}){\left\{P_D(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})\right\}}^{*}}{A_0\left(\mathrm{n\&Delta;z}\right)}(2-a)(1+{(1-a)}^2),$ Wherein

$a=\frac{A(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})}{A_0\left(\mathrm{n\&Delta;z}\right)}$

The upper and lower limit that makes the seismic data effective band is respectively l ₁Δ ω and l ₂Δ ω, l here ₁And l ₂Be positive integer, then the imaging amplitude φ of position, space deconvolution image-forming condition (x, y, n Δ z) is

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\underset{m=l_1}{\overset{{\mathrm{<figure-callout\; id="2"\; label="l"\; filenames="DEST\_PATH\_HSB00000018701500031.png,DEST\_PATH\_HSB00000018701500051.png"\; state="\{\{state\}\}">l</figure-callout>}}_2}{\mathrm{\&Sigma;}}}\mathrm{Re}\left\{\mathrm{\&beta;B}(x,y,\mathrm{n\&Delta;z},\mathrm{m\&Delta;\&omega;})\right\}$

If j is a unit imaginary number, β=jm Δ ω when doing three-D migration calculating in the following formula is when two-dimensional migration calculates $\mathrm{\&beta;}=\sqrt{\mathrm{jm\&Delta;\&omega;}};$ The real part of multiple wave field value is asked in symbol Re representative.

Imaging results φ (x, y, n Δ z) is stored in the hard disk of computing machine.Fig. 3 is the common image gather that wire size 322 and CDP number 400 place are made of all single big gun migration result; The parameter that all single big gun skews are adopted all is: Δ z=5m, the total number of sample points n=1000 on the degree of depth, Δ ω=1.40 (radian per second), l ₁=13, l ₂=337.

(4) to each shot point (focus), according to the migration velocity model of Fig. 2 a and 2b, utilize the one way wave method that frequency is carried out degree of depth continuation for the space impulse of 25Hz, get Δ z=5m, try to achieve the incident angle of the outgoing wave of shot point to underground each imaging point.

Computation process is: make that x and y are horizontal coordinates, z is a depth coordinate, and the dominant frequency of focus is ω _p(radian per second), the locus of focus are (x _s, y _s, 0), the speed of medium be c (x, y, z); Utilize the one way ripple algorithm of wave equation migration, in frequencies omega _pTo space impulse δ (x-x _s, y-y _s) carry out degree of depth continuation, obtaining with Δ z is multiple wave field value P (x, y, n Δ z, ω on each degree of depth discrete point of spacing _p, x _s, y _s); Calculate respectively

I ₀(x，y，nΔz)＝|Re{P(x，y，nΔz，ω _p，x _s，y _s)}|

$I_1(x,y,\mathrm{n\&Delta;z})=\mathrm{Re}\left\{j\frac{c(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&omega;}}_p}P(x,y,\mathrm{n\&Delta;z},{\mathrm{\&omega;}}_p,x_s,y_s)\right\}$

$I_2(x,y,\mathrm{n\&Delta;z})=\left|\frac{2}{3\mathrm{\&Delta;z}}\right[I_1(x,y,\mathrm{n\&Delta;z}+\mathrm{\&Delta;z})-I_1(x,y,\mathrm{n\&Delta;z}-\mathrm{\&Delta;z})]$

$-\frac{1}{12\mathrm{\&Delta;z}}[I_1(x,y,\mathrm{n\&Delta;z}+2\mathrm{\&Delta;z})-I_1(x,y,\mathrm{n\&Delta;z}-2\mathrm{\&Delta;z})|$

X in the following formula _s, y _sAnd ω _pBe parameter, so in the function of left end, it is ignored; Symbol || it is the Ser.No. of degree of depth up-sampling point that absolute value, n are asked in representative; (unit z) is a meter per second to speed c for x, y, and the unit of frequency is a radian per second, and the unit of the degree of depth and horizontal coordinate is a rice, and the unit of space impulse is a pressure.If I ₂(x, y, n Δ z)＞1.05I ₀(x, y, n Δ z) then makes I ₂(x, y, n Δ z)=0.Make φ ₀(n Δ z) is I ₀(x, y, n Δ z) mean value on degree of depth n Δ z is worked as I ₀(x, y, n Δ z) 〉=1.2 φ ₀(n Δ z) calculates

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\frac{I_2(x,y,\mathrm{n\&Delta;z})}{I_0(x,y,\mathrm{n\&Delta;z})}$

Work as I ₀(x, y, n Δ z)＜1.2 φ ₀(n Δ z) calculates

$\mathrm{\&phi;}(x,y,\mathrm{n\&Delta;z})=\frac{I_2(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&phi;}}_0\left(\mathrm{n\&Delta;z}\right)}(2-a)(1+{(1-a)}^2),$ Wherein $a=\frac{I_0(x,y,\mathrm{n\&Delta;z})}{{\mathrm{\&phi;}}_0\left(\mathrm{n\&Delta;z}\right)}$

φ (x, y, n Δ z) value to the space each point is carried out following processing:

1) carries out medium filtering along depth direction;

2) it is level and smooth to carry out space three-dimensional;

3) to each degree of depth n Δ z, choose maximal value,, realize the maximal value normalization on this degree of depth with the numerical value of this value except that the every bit of this degree of depth.

Making the value after the above-mentioned processing is φ (x, y, n Δ z), then can get focus to the incident angle θ (unit is degree) of underground each imaging point to be:

θ＝acos(φ(x，y，nΔz))

(x, y, n Δ z) promptly is the volume coordinate of imaging point in the formula.

Fig. 4 a and Fig. 4 b are the space incident angle distribution plans of certain shot point outgoing wave of calculating.

(5) read single big gun migration result of each big gun, according to the incident angle value of this shot point outgoing wave that calculates, migration result is added up by angle and compensate the influence of non-homogeneous covering to amplitude, form angle Dao Ji, the angle sampling interval is taken as Δ θ=2 °.

Computation process is: definition angle sampling interval Δ θ (degree), and to one-dimension array of each imaging point definition, array length is less than 90/ Δ θ; To whole big gun circulations, single big gun migration result to every big gun, read the incident angle θ of the outgoing wave of this shot point that obtains by step (4) to each imaging point, calculate l=1+int{ θ/Δ θ }, the amplitude of this imaging point is added on l the element of corresponding array, and statistics is added to the big gun number on this element simultaneously; Finish the circulation of whole big guns, the imaging amplitude array of the different depth point of same horizontal level is formed two bit array by the degree of depth, angle, promptly obtain the Jiao Daoji that compensates without non-homogeneous covering; The big gun number that adds up that obtains with statistics is removed the numerical value of the different angles of each imaging point, just obtains finishing the Jiao Daoji of non-homogeneous covering compensation.

Fig. 5 is the Jiao Daoji that certain horizontal position is not carried out non-homogeneous covering compensation.Fig. 6 is the Jiao Daoji that certain horizontal position is finished non-homogeneous covering compensation.

(6) to the circulation of road, whole angle collection, do the moving school of residue, the corresponding numerical value that obvious stretched portion occurs is changed to zero, suppress random noise with singular value decomposition method, Fig. 7 is the Jiao Daoji that finishes subsequent treatment.

(7) those skilled in the art can be concentrated by the road, angle and read the variation of different spatial place amplitude with angle, and obtain the petrophysical parameter on stratum by inversion method, and the petrophysical parameter image can be indicated the oily situation of underground structure.Consult Fig. 8 a and 8b, Fig. 9 a, 9b and 9c, Figure 10 a, 10b and 10c and Figure 11 a and 11b.

Claims (1)

1. a Wave equation depth migration directly produces road, angle diversity method, it is characterized in that adopting following steps: A) with monolateral or bilateral acquisition mode, obtain the seismic reflection signals that man-made explosion excites by wave detector, record on the tape; B) read seismic reflection signals from tape, be stored in the computing machine by the ordering of prestack big gun collection; C) utilize the one way ripple algorithm of wave equation migration respectively the source wavefield of shot point and the reception wave field of geophone station to be carried out the continuation of the wave field degree of depth along degree of depth axle, shot point after the continuation and geophone station wave field application deconvolution image-forming condition are carried out imaging; D) utilize the one way wave method to calculate the incident angle of the outgoing wave of shot point to underground each imaging point; E) incident angle of each the shot point outgoing wave that obtains according to step D will be added up and compensate the influence of non-homogeneous covering to amplitude by angle by each single big gun migration result that step C obtains, and form angle Dao Ji; F) diagonal angle Dao Ji does the moving school of residue, and the corresponding numerical value that obvious stretched portion occurs is changed to zero, suppresses random noise with singular value decomposition method, promptly obtains the Jiao Daoji of the migration before stack of reaction medium surface physics parameter variation; G) concentrated by the road, angle and read the variation of imaging amplitude in different spatial place with angle, obtain the petrophysical parameter on stratum by inversion method, the petrophysical parameter image can be indicated the oily situation of underground structure;

Describedly shot point after the continuation and geophone station wave field are used the deconvolution image-forming condition carry out imaging and be achieved in that and make T ₀Be the duration of seismologic record, frequency sampling Δ ω=2 π/T then ₀Radian per second; If Δ z is the spatial sampling of depth direction, n is the Ser.No. of degree of depth up-sampling point, and x and y are the horizontal level coordinates, and m is the Ser.No. of frequency sampling point; If shot point and geophone station wave field are through being respectively P after the degree of depth continuation on the n Δ z degree of depth _D(x, y, n Δ z, m Δ ω) and P _U(x, y, n Δ z, m Δ ω) calculates

A(x，y，nΔz，mΔω)＝P _D(x，y，nΔz，mΔω){P _D(x，y，nΔz，mΔω)} ^*

Subscript * is a complex conjugate in the formula, shot point and geophone station wave field P _DAnd P _UUnit be speed or pressure according to the wave detector type, the unit of frequency is a radian per second, the unit of the degree of depth and horizontal coordinate be rice; Make A ₀(n Δ z) is the mean value of A (x, y, n Δ z, m Δ ω) on degree of depth n Δ z, as A (x, y, n Δ z, m Δ ω) 〉=1.2A ₀(n Δ z) calculates

As A (x, y, n Δ z, m Δ ω)＜1.2A ₀(n Δ z) calculates

Wherein

The upper and lower limit that makes the seismic data effective band is respectively l ₁Δ ω and l ₂Δ ω, l here ₁And l ₂Be positive integer,

Then the imaging amplitude φ of position, space deconvolution image-forming condition (x, y, n Δ z) is

If j is a unit imaginary number, β=jm Δ ω when doing three-D migration calculating in the following formula is when two-dimensional migration calculates

The real part of multiple wave field value is asked in symbol Re representative;

Describedly utilize outgoing wave that the one way wave method calculates shot point that the incident angle of underground each imaging point is achieved in that to make that x and y are horizontal coordinates, z is a depth coordinate, and the dominant frequency of focus is ω _pRadian per second, the locus of focus are (x _s, y _s, 0), the speed of medium be c (x, y, z); Utilize the one way ripple algorithm of wave equation migration, in frequencies omega _pTo space impulse δ (x-x _s, y-y _s) carry out degree of depth continuation, obtaining with Δ z is multiple wave field value P (x, y, n Δ z, ω on each degree of depth discrete point of spacing _p, x _s, y _s); Calculate respectively

I ₀(x，y，nΔz)＝|Re{P(x，y，nΔz，ω _p，x _s，y _s)}|

X in the following formula _s, y _sAnd ω _pBe parameter, so in the function of left end, it is ignored; Symbol || it is the Ser.No. of degree of depth up-sampling point that absolute value, n are asked in representative; (unit z) is a meter per second to speed c for x, y, and the unit of frequency is a radian per second, and the unit of the degree of depth and horizontal coordinate is a rice, and the unit of space impulse is a pressure, if I ₂(x, y, n Δ z)＞1.05I ₀(x, y, n Δ z) then makes I ₂(x, y, n Δ z)=0 makes φ ₀(n Δ z) is I ₀(x, y, n Δ z) mean value on degree of depth n Δ z is worked as I ₀(x, y, n Δ z) 〉=1.2 φ ₀(n Δ z) calculates

Work as I ₀(x, y, n Δ z)＜1.2 φ ₀(n Δ z) calculates

Wherein

φ (x, y, n Δ z) value to the space each point is carried out following processing:

(1) carries out medium filtering along depth direction;

(2) it is level and smooth to carry out space three-dimensional;

(3) to each degree of depth n Δ z, choose maximal value,, realize the maximal value normalization on this degree of depth with the numerical value of this value except that the every bit of this degree of depth;

Make the value after the above-mentioned processing be

Then can get focus to the incident angle θ degree of underground each imaging point is:

(x, y, n Δ z) promptly is the volume coordinate of imaging point in the formula;

The incident angle of each the shot point outgoing wave that obtains according to step D, to add up by angle by each single big gun migration result that step C obtains and compensate the influence of non-homogeneous covering amplitude, form angle Dao Ji and be achieved in that definition angle sampling interval Δ θ degree, to one-dimension array of each imaging point definition, array length is less than 90/ Δ θ; To whole big gun circulations, single big gun migration result to every big gun, read the incident angle θ of the outgoing wave of this shot point that obtains by step D to each imaging point, calculate l=1+int{ θ/Δ θ }, the amplitude of this imaging point is added on l the element of corresponding array, and statistics is added to the big gun number on this element simultaneously; Finish the circulation of whole big guns, the imaging amplitude array of the different depth point of same horizontal level is formed two bit array by the degree of depth, angle, promptly obtain the Jiao Daoji that compensates without non-homogeneous covering; The big gun number that adds up that obtains with statistics is removed the numerical value of the different angles of each imaging point, just obtains finishing the Jiao Daoji of non-homogeneous covering compensation.

Images