CN102426387A - Seismic scattering wave imaging method - Google Patents

Abstract

The invention discloses a seismic scattering wave imaging method, overcoming the defects that the traditional reflection wave imaging method has narrow application range and can not process multiple practical geologic structures and terrain imaging and capable of imaging geologic bodies such as a complex fault block, a dike, a lenticulated ore body and the like and enhancing imaging effect. The obtained scattering trace gather is smoother, energy distribution is reasonable, signal to noise ratio is high, migration velocity is easy to obtain, and stepped noise of a converted wave is eliminated. Further, the invention also discloses a seismic scattering wave imaging system on the basis of the seismic scattering wave imaging method, used for seismic scattering wave imaging; and the seismic scattering wave imaging system has the advantages of wide application range, good imaging effect, smooth scattering trace gather and high signal to noise ratio, and the stepped noise of the converted wave can be eliminated aimed at a common scattered point scattering trace gather of the converted wave.

Description

A kind of seismic scattering wave imaging method

Technical field

The invention discloses a kind of seismic event scattering formation method, belong to seismic event imaging processing technical field.

Background technology

Even to this day; The seismic imaging technology mainly is to utilize the reflected energy of seismic event; According to the Fermat principle in the seismic motion: along the shortest perpendicular to the path time of wavefront surface, promptly ripple is along ray propagates, and ripple is the strongest along the energy of ray propagates. and ray propagates meets the Si Naier theorem; Reflective wave method only just is suitable under the mild situation of underground medium interface occurrence. and actual geological problem is often very complicated; Like metallic ore exploration, complex area petroleum prospecting, lithology complex area engineering investigation etc., exploration targets tectonic structure complicacy, mature fault, stratigraphic dip laterally suddenly change the anisotropically plastid symbiosis of different scale than steep or lithology; Can form seismic wave field very complicated, that multiple phase of wave is interfered mutually; In these cases, the reflected signal that receives on the face of land is weak, signal to noise ratio (S/N ratio) is lower, even does not receive reflection wave; On the other hand; Consider the complicacy of the geologic prospecting geologic condition of practical application; Usually in the seismic event composition reflection wave ratio be seldom, and the seismic wave propagation rule of reflection wave and other compositions is also inequality, if the seismic event of all the components is all handled according to reflection wave; The imaging results of gained and actually bigger discrepancy will occur, no Practical significance.

Although development along with technology; People on basis, proposed based on reflection wave imaging improved during based on reflection wave apart from hyp CMP gather treatment technology or common-conversion point gather treatment technology; Yet CMP, CCP are a kind of imaging processing technology that develops out to horizontal layer uniform dielectric model; This model is only better to the geologic condition treatment effect that is the horizontal layer uniform dielectric substantially strictly speaking; Yet this kind geologic condition often is difficult to run in reality, perhaps only individually has just now at the utmost point, and there is the technology barrier that is difficult to overcome in therefore traditional seismic event imaging technique.

Therefore there is following Technology Need in this area: deposit into inhomogeneous area, tectonic structure complex area (rotten body, solution cavity, metallic ore etc.) and belt of weathering etc. and seek a kind of new imaging seismic event formation method more accurately to the complicated region of underground structure (like massif zone of fracture, broken area, edge, basin, the movable area of tectonic plates), geologic medium.

Summary of the invention

The invention discloses a kind of seismic scattering wave imaging method; Overcome traditional reflection wave imaging method narrow application range, can't handle the actual multiple tectonic structure that faces and the defective of terrain imaging; Can increase imaging effect to the imaging of geologic bodies such as complex fault block, vein, lenticular ore body; And the scattering road collection of gained is more level and smooth, energy distribution is reasonable, signal to noise ratio (S/N ratio) is high, obtain migration velocity easily, and the formation method that is provided has been eliminated the stepped noise of common scattering point road collection of transformed wave.

For realizing the foregoing invention purpose, the present invention realizes through following technical proposals:

A kind of seismic scattering wave imaging method comprises the steps:

1) reads in geological data;

2) calculate each sampled point of each input channel and be mapped to the concentrated equivalent offset that belongs in scattering point road, extract scattering point road collection apart from the position; Described scattering point road collection altogether forms and adopts following method to handle:

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

H wherein: half offset distance; X: scattering point road collection CSP is to the distance of common midpoint (CMP) altogether; He: equivalent offset distance; V: scattering point speed; T: ripple from the focus to the scattering point again to the whilst on tour of acceptance point; I=1,2,3 ..., the earthquake Taoist monastic name is concentrated in the scattering point road altogether;

According to above-mentioned formula with the input channel data map to scattering point road collection altogether.

3) to step 2) the result altogether carrying out the stack velocity analysis of spectrum on the scattering point road collection;

4) step 3) velocity analysis data are carried out horizontal and vertical interpolation, speed parameter is provided for speed change is total to the collection formation of scattering point road;

5) to original input data, the velocity information according to step 4) obtains extracts scattering point road collection altogether again, carries out step 3) then, up to obtaining rational migration velocity information;

6) according to Kirchoff pre-stack time migration formula, to common scattering point road collection and the velocity analysis result that step 5) obtains, carry out the energy playback, the time section after obtaining squinting promptly gets and forms images;

Wherein for scattering wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range;

For compressional wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

In above-mentioned steps; It can be to read in prestack big gun collection or common midpoint collection scale-of-two SEGY form geological data that said step 1) is read in data, also can be the geological data of other types, for example GRISYS, SEG2 etc.; Consider versatility, preferred geological data form is the SEGY form.As is known to the person skilled in the art, the geological data file of SEGY form generally comprises three parts, and first is EBCDIC file header (3200 byte), is made up of 40 cards, is used for preserving the information that some are described seismic data volume; Second portion is that binary file head (400 byte) is used for storing some key messages of describing the SEGY file; Comprise the information such as data layout, sampling number, SI, measuring unit of SEG-Y file, these information generally are stored on the fixed position of binary file head; Third part is actual seismic trace, and every seismic trace all comprises the road header and the seismic channel data of 240 bytes.The general information such as the corresponding wire size of this seismic trace, Taoist monastic name, sampling number, terrestrial coordinate of preserving in the trace header data; Seismic channel data is that the waveform of seismic signal is taken a sample by the certain hour interval of delta t, notes this a series of discrete amplitude with certain mode again.The present invention does not limit the concrete type that data storage is become the SEGY form, and IBM or IEEE all can.

In above-mentioned steps; Step 2) ultimate principle is single square root equation of variable for the two square roots with focus-scattering point-acceptance point whilst on tour are converted into the equivalent offset distance, concentrates scattared energy to press the hyperbolic curve rule by the common scattering point road of single square root equation formation and distributes.Single square root equation can be converted into

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

(each symbol implication as stated)

According to this formula with the input channel data map to scattering point road collection altogether.

In above-mentioned steps, step 2) form formula 1 according to scattering point road collection altogether) calculate each sampled point of each input channel be mapped to concentrate the place in the scattering point road equivalent offset apart from the position, extract scattering point road collection.How being mapped to the scattering road for the seismic amplitude without time shift concentrates; Amplitude processing when wherein scattering point road collection forms altogether: in the compressional wave data, adopt amplitude by distance weighted, scattering degree of covering normalization, the processing of phase place is on average carried out phase correction according to front and back sampled point amplitude by the sampling time.During altogether scattering point road collection forms, do not run cautious the exchange at transformed wave, thus keep shot point to geophone station direction consistance compacting eliminate stepped noise.

To utmost point low signal-to-noise ratio data processing, method of the present invention is in step 2 for further) the basis on also comprise the steps:

2.1) in steps 2) the common scattering point road collection that obtains is according to the relation stack of scattering point position and equivalent offset distance, obtain the poststack section, and the skew playback handled.

2.2) judge the position of scatterer: to step 2.1) position of preliminary identification scatterer in the section that obtains, the grown form characteristic is for velocity analysis provides distinguishing rule.

In above-mentioned steps, step 3) is utilized the conventional speeds analytical approach to step 2) the result carry out the stack velocity analysis of spectrum.In the enterprising line speed analysis of scattering point road collection altogether; Its pickup velocity principle is different from conventional road collection; It is according to being: the effective scattering wave of 1. correct identification; Apart from being to be exactly effectively on 0 the road, scattering wave hyperbolic curve summit is not to be exactly invalid on 0 the road in the equivalent offset distance in equivalent offset on scattering wave hyperbolic curve summit; 2. velocity pick criterion, with the summit in equivalent offset apart from being that the formed energy group of scattering wave hyperbolic curve on 0 the road is as the criterion, and even up and the stacking image effect improves the velocity pick precision according to effective scattering wave lineups; When the scattering wave energy is very weak, can not form energy group, velocity pick can not only be as the criterion with energy group.

Above-mentioned conventional speeds analytical approach can adopt present seismic imaging field existing available method arbitrarily; As is known to the person skilled in the art; This type of velocity analysis method is widely used; Actual product is all arranged in multiple business system; The for example velocity analysis method that adopted of the velocity analysis module VELDEF among the business software FOCUS, the velocity analysis method that CGG Geoclusteur Systems Analysis Module VESPA is adopted; The principle of speed method is also widely known, and " Review of seismic velocity analysis methods " literary composition that for example is published in " exploration geophysics progress " 2006.29 (5) has promptly carried out detailed introduction to ultimate principle, development and the application conditions etc. of earthquake research field stack velocity analysis, migration velocity analysis and the typical velocity analysis method of speed tomographic inversion three major types commonly used, and is like that; Those skilled in the art can repeat no more from the conventional velocity analysis method of available technology adopting here.

In above-mentioned steps, step 4) speed interpolation is conventional speeds to be analyzed data carry out horizontal and vertical interpolation, for speed change is total to scattering point road collection formation module speed parameter is provided.Its objective is to obtaining smooth two-dimension speed model; Therefore step 4) is carried out time domain and spatial domain linear interpolation and level and smooth; Main process is following: 4.1) velocity analysis is carried out one-dimensional linear speed interpolation in the scattering point position altogether; Carry out translation at first time domain velocity analysis point constantly to t=0; In the end a time domain velocity analysis point writes down cut-off time to t=, carries out translation, then the velocity field after the time domain linear interpolation (velocity variations and time are varied to direct ratio) is carried out smoothly; 4.2): carry out the space two-dimensional linear interpolation; Obtain profile direction two-dimension speed field; Carry out translation at first spatial domain velocity analysis point to the section reference position; In the end a spatial domain velocity analysis point carries out translation to the section final position, then the velocity field after the spatial domain linear interpolation (velocity variations and space length are varied to direct ratio) is carried out smoothly.

Although above-mentionedly disclose the used speed difference approach of the present invention, however said method be merely of the present invention one preferred, and unrestrictedly only carry out the speed difference with the method.Those skilled in the art adopt other concrete grammar according to invention connotation of the present invention, are the present invention and cover as long as played the effect of speed difference.

In above-mentioned steps, step 5) is to original input data, and the velocity information according to step 4) obtains extracts scattering point road collection altogether again, repeats step 3) then, up to obtaining rational migration velocity information.

In above-mentioned steps, step 6) is used for migration imaging, and common scattering point road collection and velocity analysis result to a last step obtains carry out the energy playback, the time section after obtaining squinting.

The CSP road collection migration imaging formula of step 6) can be write as:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range.

When said ripple is transformed wave (PS), only need root-mean-square velocity v with compressional wave _RmsBe changed to the equivalent offset speed v _SemBe that transformed wave CSP road collection migration imaging formula is:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

The principle of above-mentioned transformed wave is:

Under the nonhomogeneous media condition, transformed wave is a kind of form of expression of scattering wave.The P ripple propagates into scattering point from focal point and converts the S ripple into, and the total whilst on tour that propagates into acceptance point with the S ripple again is T=T _p+ T _s(with reference to figure 7).Seismic wave is with velocity equivalent v hypothetically _eFrom focus arrival scattering point, arrive acceptance point then, used whilst on tour equals T, introduces equivalent offset simultaneously apart from h _e, then have

$T=\frac{{({z_0}^2+{h_s}^2)}^{\frac{1}{2}}}{V_p}+\frac{{({z_0}^2+{h_r}^2)}^{\frac{1}{2}}}{V_s}=\frac{{({z_0}^2+{h_e}^2)}^{\frac{1}{2}}}{V_e}+\frac{{({z_0}^2+{h_e}^2)}^{\frac{1}{2}}}{V_e}$

Wherein: x is that the mid point (CMP) of focal point-acceptance point arrives the distance of scattering point at floor projections point (O), and h is the half the of geophone offset, h _eBe equivalent offset distance, h _sFor focal point S to the distance of scattering point at floor projections point, h _rFor acceptance point R to the distance of scattering point at floor projections point.V _pBe p wave interval velocity, V _sBe S wave velocity h _s=x+h, h _r=x-h.

By above-mentioned visible, the relation between TWT and the equivalent offset distance just is reduced to single square root equation, is double curve.It is thus clear that in the method for the invention, the abstracting method of transformed wave CSP road collection can be identical with compressional wave.

Seismic scattering wave imaging method of the present invention to recording geometry variation, the unequal problem of degree of covering, is taked the collection degree of covering normalization of scattering point road; The strong reflection amplitude according to the weighting of scattering wave propagation distance, is played the scattering curve smooth interaction; For utmost point low signal-to-noise ratio data, formation is pressed the zero-offset section processing with all common scattering point road collection after being total to scattering point road collection; Utilize and interfere superposition principle; Can improve signal to noise ratio (S/N ratio), and not need velocity information, tentatively judge the scatterer position in initial the interference when superposeing; To the sorting of input channel set reason, do not carry out the cautious exchange of big gun, give up the stepped energy of dispersing noise of formation, improve the signal to noise ratio (S/N ratio) of transformed wave CSP road collection.

Further; On the said method basis; The invention also discloses a kind of seismic scattering ripple imaging system; Comprise that common scattering point Dao Ji road collection forms module, interferes laminating module, velocity analysis module, speed interpolating module, migration imaging module, said scattering point Dao Ji road collection altogether form module comprise the normal speed of P ripple altogether scattering point Dao Ji road collection form submodule, P ripple speed change altogether scattering point Dao Ji road collection form submodule, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule.

Said velocity analysis module is carried out on existing business software; Belong to prior art; Repeat no more here; Those skilled in the art can select any existing proper speed analysis module as required, well-known FOCUS velocity analysis module VELDEF, the CGG Geoclusteur Systems Analysis Module VESPA etc. of comprising of this generic module.

In said system; The said module of scattering point Dao Ji road collection formation altogether is used for after reading in geological data, calculating each sampled point of each input channel and is mapped to the concentrated equivalent offset that belongs in scattering point road apart from the position; Extract scattering point road collection, and according to the following equation with the input channel data map to scattering point road collection altogether:

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

H wherein: half offset distance; X: scattering point road collection CSP is to the distance of common midpoint (CMP) altogether; He: equivalent offset distance; V: scattering point speed; T: ripple from the focus to the scattering point again to the whilst on tour of acceptance point; I=1,2,3 ..., the earthquake Taoist monastic name is concentrated in the scattering point road altogether.

As stated; Altogether scattering point Dao Ji road of the present invention collection form module comprise the normal speed of P ripple altogether scattering point Dao Ji road collection form submodule, P ripple speed change altogether scattering point Dao Ji road collection form submodule, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule; All can adopt scattering point road collection altogether of the present invention to form module to different scattering wave types handles the geological data that reads in

In seismic scattering ripple imaging system of the present invention, be total to scattering point road collection, transformed wave speed change scattering point road collection altogether to the normal speed of transformed wave, can it be adopted with p phase of wave abstracting method together and handle, it is according to being:

Under the nonhomogeneous media condition, transformed wave is a kind of form of expression of scattering wave.The P ripple propagates into scattering point from focal point and converts the S ripple into, and the total whilst on tour that propagates into acceptance point with the S ripple again is T=T _p+ T _s(with reference to figure 2).Seismic wave is with velocity equivalent v hypothetically _eFrom focus arrival scattering point, arrive acceptance point then, used whilst on tour equals T, introduces equivalent offset simultaneously apart from h _e, then have

$T=\frac{{({z_0}^2+{h_s}^2)}^{\frac{1}{2}}}{V_p}+\frac{{({z_0}^2+{h_r}^2)}^{\frac{1}{2}}}{V_s}=\frac{{({z_0}^2+{h_e}^2)}^{\frac{1}{2}}}{V_e}+\frac{{({z_0}^2+{h_e}^2)}^{\frac{1}{2}}}{V_e}$

Wherein: x is that the mid point (CMP) of focal point-acceptance point arrives the distance of scattering point at floor projections point (O), and h is the half the of geophone offset, h _eBe equivalent offset distance, h _sFor focal point S to the distance of scattering point at floor projections point, h _rFor acceptance point R to the distance of scattering point at floor projections point.V _pBe p wave interval velocity, V _sBe S wave velocity h _s=x+h, h _r=x-h.

Relation between TWT and the equivalent offset distance just is reduced to single square root equation, is double curve.By above-mentioned visible, in seismic scattering ripple imaging system of the present invention, the abstracting method of transformed wave CSP road collection can be identical with compressional wave.

Know based on above-mentioned, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule and adopt with corresponding p ripple and be total to the same abstracting method of scattering point Dao Ji road collection formation submodule.

Above-mentioned interference laminating module is used to handle utmost point low signal-to-noise ratio data processing; Utilizing altogether, scattering point road collection has self excitation and self receiving zero-offset section characteristic and relevant principle; Do not need velocity information; The common scattering point road collection that obtains of step on all is obtained the poststack section according to the relation stack of scattering point position and equivalent offset distance, and skew playback processing.

Above-mentioned speed interpolating module is used for the data of conventional speeds analysis module are carried out horizontal and vertical interpolation; For speed change is total to scattering point road collection formation module speed parameter is provided; And with resulting velocity information; Again extract scattering point road collection altogether, then data are reached the velocity analysis module, repeat up to obtaining rational migration velocity information.

Above-mentioned migration imaging module is used for carrying out the energy playback to being total to scattering point road collection and the resultant result of velocity analysis module, the time section after obtaining squinting, thereby imaging, and it forms images according to following method: CSP road collection migration imaging formula can be write as:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range.

When said ripple is transformed wave (PS), only need root-mean-square velocity v with compressional wave _RmsBe changed to the equivalent offset speed v _SemBe that transformed wave CSP road collection migration imaging formula is:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

Seismic scattering wave imaging method of the present invention and be a kind of based on the imaging solution of scattering point road collection altogether based on this seismic scattering ripple imaging system; Overcome the defective that reflection wave can't widespread use; Solved the scattering point road collection problem altogether that extracts in the scatterer imaging; Can be applied to geologic body imagings such as complex fault block, vein, lenticular ore body, increase imaging effect, and scattering road collection is more level and smooth, energy distribution is reasonable; And signal to noise ratio (S/N ratio) is high, has eliminated the transformed wave stepped noise of scattering point road collection altogether.

Description of drawings

Fig. 1 reads in the input channel data map to the synoptic diagram that is total on the scattering point road collection for method of the present invention;

Fig. 2 is a transformed wave whilst on tour synoptic diagram;

Fig. 3 a-3d is for extracting face of land difference 200m (csp40), 300m (csp60), 500m (csp100), the common scattering point road collection of 600m (csp120);

Fig. 4 representes the process synoptic diagram with a plurality of scattering point road collection altogether stack of Fig. 3;

Fig. 5 is Fig. 4 gained stacked section;

Fig. 6 is transformed wave CSP road collection big gun, the cautious back principle schematic of the generation whilst on tour time difference of exchanging;

The scattering wave hyperbolic curve of Fig. 7 for shining upon on the CSP of the transformed wave as a result road collection based on Fig. 6;

Fig. 8 a, 8b are for adopting the equivalent offset based on CSP road collection of the present invention apart from migrated section and traditional conventional post-stack migration section based on CMP road collection respectively;

Fig. 9 is the common scattering point road collection to inclination tabular molding middle part;

Figure 10 a, 10b are for being CMP gather velocity spectrum, common scattering point road collection velocity spectrum;

Figure 11 a, 11b are that the common scattering point road of the transformed wave collection that contains stepped noise, the transformed wave of eliminating stepped noise are total to scattering point road collection;

Figure 12 is the workflow diagram of seismic scattering ripple imaging system of the present invention.

Embodiment

Can better understand connotation of the present invention in conjunction with accompanying drawing; Yet the accompanying drawing of giving only is used to the result to a certain concrete geology exemplary application the inventive method gained is described, is not to be used to limit the present invention or to limit the present invention only be used to handle appended embodiment.

The invention discloses a kind of seismic scattering wave imaging method, comprise the steps

1. read in geological data;

2. calculate each sampled point of each input channel and be mapped to the concentrated equivalent offset that belongs in scattering point road, extract scattering point road collection apart from the position; Described scattering point road collection altogether forms and adopts following method to handle:

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

H wherein: half offset distance; X: scattering point road collection CSP is to the distance of common midpoint (CMP) altogether; He: equivalent offset distance; V: scattering point speed; T: ripple from the focus to the scattering point again to the whilst on tour of acceptance point; I=1,2,3 ...., the earthquake Taoist monastic name;

According to above-mentioned formula with the input channel data map to scattering point road collection altogether.

3. the result to step 2 is carrying out the stack velocity analysis of spectrum on the scattering point road collection altogether; The principle of its foundation is the effective scattering wave of 1. correct identification: apart from being to be exactly effectively on 0 the road, scattering wave hyperbolic curve summit is not to be exactly invalid on 0 the road in the equivalent offset distance in equivalent offset on scattering wave hyperbolic curve summit; 2. velocity pick criterion: with the summit in equivalent offset apart from being that the formed energy group of scattering wave hyperbolic curve on 0 the road is as the criterion, and even up with the stacking image effect decision according to effective scattering wave lineups and pick up precision.

4. step 3 velocity analysis data are carried out horizontal and vertical interpolation, speed parameter is provided for speed change is total to the collection formation of scattering point road;

5. to original input data, the velocity information according to step 4 obtains extracts scattering point road collection altogether again, carry out step 3 then, up to obtaining rational migration velocity information;

6. the common scattering point road collection and the velocity analysis result that step 5 are obtained carry out the energy playback, and the time section after obtaining squinting promptly gets and forms images;

Wherein for scattering wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range.

For compressional wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

To utmost point low signal-to-noise ratio data processing, method of the present invention also comprises the steps: on the basis of step 2 for further

2.1 in steps the 2 common scattering point road collection that obtain obtain the poststack section according to the relation stack of scattering point position and equivalent offset distance, and the skew playback is handled.

2.2 judge the position of scatterer: the position of preliminary identification scatterer in the section that step 2.1 is obtained, the grown form characteristic is for velocity analysis provides distinguishing rule.

In conjunction with accompanying drawing 3-5, be appreciated that formation method of the present invention carries out the stack of CSP road collection based on step (2), (2.1), (2.2), (3), form one have high letter make than the process of stacked section.

As shown in Figure 3, underground have scattering point A (200m, 300m) and B (500m, 100m); Medium velocity position 2000m/s extracts face of land difference 200m (csp40) respectively, 300m (csp60), 500m (csp100); The common scattering point road collection (Fig. 3 a-d) of 600m (csp120), from A, the scattared energy that B is 2 is according to separately hyperbola distribution, and the hyperbolic curve summit is positioned at 100ms and 300ms place respectively; Equal 300 meters of hyperbolic curve apex horizontal spacing, each common scattering point road is concentrated, and the zero-offset road is the 161st road; Single-side shift is 160 roads apart from number (EON), and equivalent offset is 5 meters apart from interval (EOS), and maximum equivalent offset distance (MEO) is 800 meters.According to CSP road collection diverse location; With all CSP roads collection stacks (Fig. 4), form one have high letter make than stacked section (STACK), as shown in Figure 5; Whole additive process does not need velocity information, therefore visible seismic scattering wave imaging method additive process of the present invention do not rely on velocity information.

With reference to figure 6,7, be appreciated that traditional shot point and acceptance point exchange and cause the whilst on tour time difference, thereby causing that the scattering wave hyperbolic curve is stepped on the transformed wave CSP road collection disperses.As far as the PS ripple, incident wave is different with outgoing wave speed, has produced the whilst on tour time difference after causing big gun, the cautious exchange, and has increased along with the increase of geophone offset.As shown in Figure 6, on the CMP road of transformed wave collection, after shot point and acceptance point exchanged, whilst on tour became t2 by t1.T1 is shot point coordinate 100m, the seismic event whilst on tour when the geophone station coordinate is 300m, and t2 is the seismic event whilst on tour after the cautious exchange of big gun.If directly form the CSP road collection that formula obtains transformed wave according to the scattering point road collection of simplifying, then twice mapping back produces the time difference (T2-T1), T2; (shown in Figure 6 for the scattering point position be (0m, 200m), velocity equivalent=2000m/s, vp=2000m/s, vp/vs=1.8, t1 be shot point (100m, 0m), acceptance point (300m, the seismic event whilst on tour in the time of 0m), t2 to be that big gun is cautious exchange after seismic event whilst on tour).To depart from the scattering wave hyperbolic curve that the CSP road is concentrated, this time difference increases along with the increase of geophone offset, and stepped the dispersing of scattering wave hyperbolic curve of on transformed wave CSP road collection, shining upon is as shown in Figure 7.It is thus clear that, produced the whilst on tour time difference after big gun, the cautious exchange, be to cause the stepped basic reason of dispersing of scattering wave hyperbolic curve on the transformed wave CSP road collection.

With reference to figure 8a, 8b, Fig. 8 a is for adopting the equivalent offset based on CSP road collection of the present invention apart from migrated section, and Fig. 8 b is for adopting traditional conventional post-stack migration section based on CMP road collection.What it was directed against is that metallic ore intrusive mass model designs, and owing to the metal ore body is lain concealed under intrusive vein, and by the intrusive mass shielding, it is big that it surveys difficulty, in metallic ore is surveyed, has general representativeness.Therefore generally adopt this model to be used to weigh imaging effect in the seismic imaging field; Can find out by figure; Seismic scattering wave imaging method based on common scattering point of the present invention is compared with traditional formation method based on common midpoint; To complicated geological, imaging is concentrated, more clear more, and imaging effect is better.

Fig. 9 is used to tilt the common scattering point road collection at tabular molding middle part for adopting formation method of the present invention, and this model is abstract in mathematics physics model designs to geological phenomenons such as subsurface fault shatter belt, inclined thin rock mass, and is widely known in the seismic imaging field.Can find out that by figure formation method scattering of the present invention road collection is more level and smooth, energy distribution is reasonable.

Figure 10 a is the CMP gather velocity spectrum, and Figure 10 b is for being total to scattering point road collection velocity spectrum, to same geologic condition.By finding out among the figure, scattering point Dao Ji road collection signal to noise ratio (S/N ratio) is obviously heightened altogether, obtains migration velocity easily with the conventional speeds analysis.From figure, can see that compare with Figure 10 a, Figure 10 b degree of covering has increased by 1 times than Figure 10 a, signal to noise ratio (S/N ratio) obviously improves; The lineups continuity strengthens, and smoother; From velocity spectrum, can find out obviously that the following multiple reflection of 3s is by significantly compacting.

With reference to figure 11a, 11b; To subsurface anomaly body yardstick when close with the seismic event wavelength; Can the anomalous body boundary scattering be approximately the point scattering model; The point scattering model is the basic model of research scattering problem; Institute's diagrammatic sketch is total to scattering point road collection (11b) for adopting method of the present invention to the transformed wave that the transformed wave that contains stepped noise is total to the stepped noise of elimination that obtains after scattering point road collection (11a) is handled, and visible by the contrast of 11a and 11b, method of the present invention has been eliminated the transformed wave stepped noise of scattering point road collection altogether.

With reference to Figure 12, be appreciated that the principle of work of seismic scattering ripple imaging system of the present invention, Figure 12 has showed that the course of work of corresponding each module of imaging system of the present invention and streams data each other connect.Seismic scattering ripple imaging system disclosed by the invention; Comprise that common scattering point Dao Ji road collection forms module, interferes laminating module, velocity analysis module, speed interpolating module, migration imaging module; Said altogether scattering point Dao Ji road collection form module comprise the normal speed of P ripple altogether scattering point Dao Ji road collection form submodule, P ripple speed change altogether scattering point Dao Ji road collection form submodule, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule; Said velocity analysis module is carried out on existing business software; Belong to prior art; Repeat no more, those skilled in the art can select any existing proper speed analysis module as required here.

In said system; The said module of scattering point Dao Ji road collection formation altogether is used for after reading in geological data, calculating each sampled point of each input channel and is mapped to the concentrated equivalent offset that belongs in scattering point road apart from the position; Extract scattering point road collection, and with the input channel data map to scattering point road collection altogether:

Wherein, Altogether scattering point Dao Ji road of the present invention collection form module comprise the normal speed of P ripple altogether scattering point Dao Ji road collection form submodule, P ripple speed change altogether scattering point Dao Ji road collection form submodule, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule; All can adopt scattering point road collection altogether of the present invention to form module to different scattering wave types handles the geological data that reads in

Interfere laminating module to be used to handle utmost point low signal-to-noise ratio data processing; Utilizing altogether, scattering point road collection has self excitation and self receiving zero-offset section characteristic and relevant principle; Do not need velocity information; The common scattering point road collection that obtains of step on all is obtained the poststack section according to the relation stack of scattering point position and equivalent offset distance, and skew playback processing.

The speed interpolating module is used for the data of conventional speeds analysis module are carried out horizontal and vertical interpolation; For speed change is total to scattering point road collection formation module speed parameter is provided; And with resulting velocity information; Again extract scattering point road collection altogether, then data are reached the velocity analysis module, repeat up to obtaining rational migration velocity information.

The migration imaging module is used for carrying out the energy playback to being total to scattering point road collection and the resultant result of velocity analysis module, the time section after obtaining squinting, thereby imaging.

Claims (8)

1. a seismic scattering wave imaging method comprises the steps:

1) reads in geological data;

2) calculate each sampled point of each input channel and be mapped to the concentrated equivalent offset that belongs in scattering point road, extract scattering point road collection apart from the position; Described scattering point road collection altogether forms and adopts following method to handle:

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

H wherein: half offset distance; X: scattering point road collection CSP is to the distance of common midpoint CMP altogether; He: equivalent offset distance; V: scattering point speed; T: ripple from the focus to the scattering point again to the whilst on tour of acceptance point; I=1,2,3 ..., the earthquake Taoist monastic name is concentrated in the scattering point road altogether;

According to above-mentioned formula with the input channel data map to scattering point road collection altogether;

3) result to step (2) is carrying out the stack velocity analysis of spectrum on the scattering point road collection altogether;

4) step (3) velocity analysis data are carried out horizontal and vertical interpolation, speed parameter is provided for speed change is total to the collection formation of scattering point road;

5) to original input data, the velocity information according to step (4) obtains extracts scattering point road collection altogether again, carries out step (3) then, up to obtaining rational migration velocity information;

6) the common scattering point road collection and the velocity analysis result that step (5) are obtained carry out the energy playback, and the time section after obtaining squinting promptly gets and forms images;

Wherein for scattering wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range; For compressional wave, its formation method is according to following:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

2. method according to claim 1 is characterized in that said step 1) geological data is a scale-of-two SEGY form geological data.

3. method according to claim 1; It is characterized in that said step 2) amplitude processing adopted following method when scattering point road collection formed altogether: in the compressional wave data, adopt amplitude by distance weighted, scattering degree of covering normalization, the processing of phase place is on average carried out phase correction according to front and back sampled point amplitude by the sampling time.

4. method according to claim 1 is characterized in that said step 2) during altogether scattering point road collection forms, do not run cautious the exchange at transformed wave, keep shot point to the compacting of geophone station direction consistance, eliminate stepped noise.

5. method according to claim 1 is characterized in that said step 2) be used for utmost point low signal-to-noise ratio data processing and also comprise the steps:

2.1) in steps 2) the common scattering point road collection that obtains is according to the relation stack of scattering point position and equivalent offset distance, obtain the poststack section, and the skew playback handled.

2.2) judge the position of scatterer: to step 2.1) position of preliminary identification scatterer in the section that obtains, the grown form characteristic is for velocity analysis provides distinguishing rule.

6. method according to claim 1; It is characterized in that said step 3) adopts following principle to carry out velocity analysis: the effective scattering wave of 1. correct identification; Apart from being to be exactly effectively on 0 the road, scattering wave hyperbolic curve summit is not to be exactly invalid on 0 the road in the equivalent offset distance in equivalent offset on scattering wave hyperbolic curve summit; 2. velocity pick criterion, with the summit in equivalent offset apart from being that the formed energy group of scattering wave hyperbolic curve on 0 the road is as the criterion, and even up and the stacking image effect improves the velocity pick precision according to effective scattering wave lineups.

7. method according to claim 1; It is characterized in that said step 4) is carried out time domain and the spatial domain linear interpolation is also level and smooth; Comprise the steps: 4.1) velocity analysis is carried out one-dimensional linear speed interpolation in the scattering point position altogether, carry out translation at first time domain velocity analysis point constantly to t=0, in the end a time domain velocity analysis point writes down cut-off time to t=; Carry out translation, then the velocity field after the time domain linear interpolation is carried out smoothly; 4.2): carry out the space two-dimensional linear interpolation; Obtain profile direction two-dimension speed field; Carry out translation at first spatial domain velocity analysis point to the section reference position; In the end a spatial domain velocity analysis point carries out translation to the section final position, then the velocity field after the spatial domain linear interpolation is carried out smoothly.

8. seismic scattering ripple imaging system; Comprise that common scattering point Dao Ji road collection forms module, interferes laminating module, velocity analysis module, speed interpolating module, migration imaging module; Said altogether scattering point Dao Ji road collection form module comprise the normal speed of P ripple altogether scattering point Dao Ji road collection form submodule, P ripple speed change altogether scattering point Dao Ji road collection form submodule, the normal speed of transformed wave altogether scattering point Dao Ji road collection form submodule, transformed wave speed change altogether scattering point Dao Ji road collection form submodule; It is characterized in that the said module of scattering point Dao Ji road collection formation altogether is used for after reading in geological data, calculating each sampled point of each input channel and is mapped to the concentrated equivalent offset that belongs in scattering point road apart from the position; Extract scattering point road collection, and according to the following equation with the input channel data map to scattering point road collection altogether:

$t_i=\frac{2\mathrm{xh}}{v_{\mathrm{mig}}{(x^2+h^2-{h_{\mathrm{ei}}}^2)}^{\frac{1}{2}}}$

H wherein: half offset distance; X: scattering point road collection CSP is to the distance of common midpoint (CMP) altogether; He: equivalent offset distance; V: scattering point speed; T: ripple from the focus to the scattering point again to the whilst on tour of acceptance point; I=1,2,3 ..., the earthquake Taoist monastic name is concentrated in the scattering point road altogether;

Described interference laminating module is used to handle utmost point low signal-to-noise ratio data processing; Utilizing altogether, scattering point road collection has self excitation and self receiving zero-offset section characteristic and relevant principle; Do not need velocity information; The common scattering point road collection that obtains of step on all is obtained the poststack section according to the relation stack of scattering point position and equivalent offset distance, and skew playback processing;

Described speed interpolating module is used for the data of conventional speeds analysis module are carried out horizontal and vertical interpolation; For speed change is total to scattering point road collection formation module speed parameter is provided; And with resulting velocity information; Again extract scattering point road collection altogether, then data are reached the velocity analysis module, repeat up to obtaining rational migration velocity information;

Described migration imaging module is used for carrying out the energy playback to being total to scattering point road collection and the resultant result of velocity analysis module, the time section after obtaining squinting, thereby imaging, and it forms images according to following method:

CSP road collection migration imaging formula can be write as:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{rms}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{rms}}})]$

Wherein: P _Out(x, τ t=0) are imaging results,

Be the input wave field of CSP road collection, ρ (t) expression is to the time-derivative of input wave field, and symbol * is a convolution operator, and x is the position of scattering point road collection, h _eFor concentrating and put equivalent focus and equivalent received is put the distance of scattering point at surface projection, Δ h in the scattering point road _eFor equivalent offset apart from interval, ∑ is illustrated in summation in the optimized migration pore diameter range;

When said ripple is transformed wave PS, only need root-mean-square velocity v with compressional wave _RmsBe changed to the equivalent offset speed v _SemBe that transformed wave CSP road collection migration imaging formula is:

$P_{\mathrm{out}}(x,\mathrm{\τ},t=0)=\frac{{\mathrm{\Δh}}_e}{2\mathrm{\π}}\underset{h_e}{\mathrm{\Σ}}[\frac{\cos \mathrm{\θ}}{\sqrt{v_{\mathrm{sem}}r}}\mathrm{\ρ}\left(t\right)*P_{\mathrm{in}}(x,h_e,t-\frac{r}{v_{\mathrm{sem}}})].$

Images