CN107255832B - A kind of inversion method of subsurface structure - Google Patents
A kind of inversion method of subsurface structure Download PDFInfo
- Publication number
- CN107255832B CN107255832B CN201710554305.0A CN201710554305A CN107255832B CN 107255832 B CN107255832 B CN 107255832B CN 201710554305 A CN201710554305 A CN 201710554305A CN 107255832 B CN107255832 B CN 107255832B
- Authority
- CN
- China
- Prior art keywords
- wave field
- residual error
- data
- reflected wave
- data residual
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
Abstract
The embodiment of the present application discloses a kind of inversion method of subsurface structure, which comprises calculates the first data residual error according to observation data, obtains diffracted wave field and reflected wave field in the first data residual error；According to preset reflected wave field weighting function and the reflected wave field, diffracted wave field, the second data residual error is determined；According to the second data residual error, it is iterated inversion imaging.The present processes utilize weighting function by separation diffracted wave field and reflected wave field, are weighted to reflected wave field, are iterated inverting to generate the second data residual error, realize the high-precision imaging of heterogeneous geology objective body；The weight of back wave, the resolution of the inhomogeneous structures such as very high fracture hole of lower opening over the ground can also be constantly reduced in an iterative process by changing weighting function.
Description
Technical field
The present invention relates to a kind of inversion methods of subsurface structure, belong to oil gas physical prospecting engineering field.
Background technique
In carbonate reservoir, hole, seam, hole construction be rich accumulation of oil and gas important place, the Important Problems of exploration be as
What realizes the identification and imaging of the heterogeneous constructions such as hole, seam, hole.The seismic response of heterogeneous construction is mainly shown as diffracted wave
?.
But diffracted wave is weaker relative to reflected energy, in imaging, the corresponding heterogeneous construction of diffracted wave would generally
It is submerged in Background Construction, it is difficult to identify.Further, since diffracted wave field energy is too weak, it is imaged individually extracting diffracted wave
When, the result of conventional image forming process is more vulnerable to noise jamming, exists simultaneously that diffracted wave is not restrained, imaging resolution is low etc. and asks
Topic.
Based on this, it would be desirable to more accurate subsurface structure inversion method.
Summary of the invention
The embodiment of the present application provides a kind of inversion method of subsurface structure, for solving at least one of following problems:
It is heterogeneous to construct the problem of being difficult to when inverting subsurface structure, or, heterogeneous structure imaging is inaccurate when inverting subsurface structure
True problem.
Based on this, the embodiment of the present application provides a kind of inversion method of subsurface structure, which comprises
The first data residual error is calculated according to observation data, obtains diffracted wave field and reflected wave field in the first data residual error；
According to preset reflected wave field weighting function and the reflected wave field, diffracted wave field, the second data residual error is determined；
According to the second data residual error, it is iterated inversion imaging.
Compared to existing technologies, the embodiment of the present application is by separating diffracted wave field and reflection in the first data residual error
Wave field, and weighting function is utilized, reflected wave field is weighted, to generate the second data residual error (namely weighting reflected wave field
Data residual error afterwards) it is iterated inversion imaging.The heterogeneous geologic structure such as lower opening fracture hole over the ground may be implemented in the method for the present invention
High-precision imaging；The weight of back wave can also be constantly reduced in an iterative process by changing weighting function, improve over the ground
The resolution of the inhomogeneous structures such as lower opening fracture hole.
Detailed description of the invention
Fig. 1 is the flow diagram of one embodiment provided by the embodiments of the present application；
Fig. 2 is the migration velocity field of one embodiment of the present of invention；
Fig. 3 is conventional least square reverse-time migration imaging results；
Fig. 4 is conventional diffracted wave reverse-time migration imaging results；
Fig. 5 is the diffracted wave least-squares migration imaging results of one embodiment of the present of invention；
Fig. 6 is the all-wave field least square reverse-time migration imaging results of one embodiment of the present of invention；
Fig. 7 is the embodiment of the present invention and the wave-number spectrum comparison of other two methods results；
Fig. 8 is the normalization data residual error convergence of the embodiment of the present invention and conventional least square reverse-time migration methods and results
Curve；
The present invention is further illustrated with reference to the accompanying drawings and detailed description.
Specific embodiment
To keep the purposes, technical schemes and advantages of the application clearer, below in conjunction with the application specific embodiment and
Technical scheme is clearly and completely described in corresponding attached drawing.Obviously, described embodiment is only the application one
Section Example, instead of all the embodiments.Based on the embodiment in the application, those of ordinary skill in the art are not doing
Every other embodiment obtained under the premise of creative work out, shall fall in the protection scope of this application.
The embodiment of the present application provides a kind of inversion method of subsurface structure, as shown in Figure 1, comprising the following steps:
S101 calculates the first data residual error according to observation data, obtains the diffracted wave field in the first data residual error and reflection
Wave field.Specifically, include the following steps:
Step a obtains model data and observes the difference of data, wherein the model data is according to observation data
Determined by imaging results and the number of iterations；
It is 0 that initial imaging results, which are arranged, and obtaining first imaging results according to observation data, (i.e. imaging results need not
Stop iteration), inverse time inverse migration is carried out to current imaging results, obtains model data.
Specifically, being realized by inverse time inverse migration computation model data by solving following wave equation:
Wherein, v0It (x) is migration velocity field, p (x, t) is the model data that inverse time inverse migration is calculated,It is general to draw
Laplacian operater, m (x) are imaging results, p0(x, t) is background wave field.Solve the wave equation of background wave field are as follows:
Wherein, f (t) is source wavelet.
It obtains model data and observes the difference of data, calculating process is indicated in the form of matrix operation are as follows:
r(k)=Lm(k-1)-d (3)
Wherein, r(k)For the difference corresponding in kth iteration, that is, the first data residual error.
L is inverse time inverse migration operator, m(k-1)For the imaging results obtained when -1 iteration of kth, that is, to current imaging
As a result inverse time inverse migration is carried out, is obtained model data p (x, t), d is observation data, and inverse time inverse migration passes through shown in solution formula (1)
Wave equation carry out.
Due to being carried out by the way of iteration in refutation process, each iteration can all change current imaging results, from
And the model data and and the first data residual error be also required to correspondence and be iterated.
Step b isolates diffracted wave field included in the difference.
After obtaining the difference, the difference is transformed to common offset trace gather by common-shot-gather, to its into
The conventional diffracted wave separation based on plane wave destructing filter of row, obtains diffracted wave field:
D(k)=Sr(k) (4)
Wherein, D(k)For diffracted wave field included in kth time iteration time difference value, S is the plane wave for carrying out diffracted wave separation
Deconstruct filter.
Since in common offset trace gather, back wave corresponds to the pseudo-linear lineups of more continuous and derivable, diffracted wave pair
It should be in discontinuous, rough hyperbolic lineups.It is estimated that the part of lineups is inclined by deconstructing filter with plane wave
Angle removes the reflected wave field of continuous and derivable, obtains diffracted wave field.
Step c determines reflected wave field according to the difference and diffracted wave field.
It is determined by following formula:
R(k)=r(k)-D(k) (5)
Wherein, R(k)Reflected wave field when iteration secondary for kth in the difference.
S102 determines that the second data are residual according to preset reflected wave field weighting function and the reflected wave field, diffracted wave field
Difference.Specifically, the second data residual error can be determined in the following manner:
Wherein,For the second data residual error, D(k)For diffracted wave field,For reflected wave field weighting function, R(k)For reflection
Wave field, k are the number of iterations；Wherein, the reflected wave field weighting function is equal to 0, or, the reflected wave field weighting function is greater than 0
And less than 1.
If imageable target is the heterogeneous bodies such as hole, seam, hole, the tomography of underground, it is arrangedIt is 0, then it is only right at this time
Diffracted wave field imaging.To improve the image quality of diffracted wave field subsequently through iteration more new algorithm, underground is provided anisotropically
The exact image of plastid.
If imageable target is all subsurface structures, it is arrangedThe reflected wave field weighting functionGreater than 0 and it is less than
1, andIt is gradually reduced with the increase of the number of iterations.There are many structural forms for weighting function at this time, can be according to imaging mesh
Mark is artificial given.For example, settingFor reflected wave field weighting function, b is control reflected wave field weighting coefficient size
Parameter, b > 0, so as to the specific gravity that diffracted wave field prominent in successive iterations refutation process accounts in objective function, and then improve around
The imaging resolution of the corresponding heterogeneous body in ejected wave field.
S103 is iterated inversion imaging according to the second data residual error.
At this time according to the second data residual error, objective function and more new direction are determined, objective function is iterated instead
It drills.
It can be by constructing following objective function:It is iterated inverting.Wherein J
It (m) is objective function, m is imaging results, and D is diffracted wave field, R and wRRespectively reflected wave field and its reflected wave field weight letter
Number, TλFor carry out L1 modular constraint threshold function table, to suppress in imaging results be less than threshold value signal.
Determine the formula of more new direction are as follows:
Wherein, J (m) is objective function, D*(x, t) and R*(x, t) is diffracted wave field D (x in data residual errorg, t) and reflection
Wave field R (xg, t) adjoint wave field, p0(x, t) is background wave field, and x is underground any point, xgFor receiving point.Diffracted wave field
It can be calculated by following wave equation with wave field:
Similarly, the adjoint wave field of reflected wave field can be calculated by following wave equation:
Then judge whether more new direction meets error condition, error condition are as follows: 2 norms of more new direction are less than ε, and ε is
One, close to 0 real number, can rule of thumb be set.Continue iteration if being unsatisfactory for error condition, meets error condition
When iteration ends,
When 2 norms of more new direction are greater than ε, i.e., when error condition is unsatisfactory for, imaging results are carried out according to more new direction
Iteration updates.Specifically, imaging results are iterated with the formula of update are as follows:
Wherein, m(k)For the imaging results that kth iterates to calculate, m(k-1)Obtained imaging knot is iterated to calculate for kth -1
Fruit,Update step-length when iteration secondary for kth, g(k)More new direction when iteration secondary for kth, update side
It is solved to by formula (7).The threshold function table of L1 modular constraint is carried out when kth time iteration, it can be according to the signal-to-noise ratio of observation data
Artificial given, the signal-to-noise ratio for observing data is higher, and threshold value is smaller.L1 modular constraint is realized by threshold value shrinkage method, in the number of iterations
When smaller, biggish threshold value is selected rapidly to suppress the offset noise in imaging results, with the increase of the number of iterations, threshold value by
It is decrescence small, to update the corresponding geological structure of echo signal.
When 2 norms of more new direction be less than ε when, determine current imaging results be ultimately imaged as a result, export finally at
As result.
Using method of the invention, back wave is weighted by the weighting function for being arranged different, to determine the second number
According to residual error, when diffracted wave is imaged according to the second data residual error (when namely only the inhomogeneous structures such as slit being imaged), this
Inventive method can improve the image quality of diffracted wave field by iteration more new algorithm, provide the accurate of underground heterogeneous body
Image；When according to the second data residual error to full Seismic imaging (when namely whole underground structures being imaged), the method for the present invention
During reflected wave field weighting function is gradually reduced, the corresponding large scale homophaneous structure of reflected wave field is updated first, is improved
Specific gravity of the reflected wave field in objective function is gradually reduced while its image quality, prominent diffracted wave field improves diffracted wave
The resolution of the corresponding inhomogeneous structure in field, may finally obtain the high precision image of underground possessive construction.
To make technical solution of the present invention is more plain to understand, the example for providing a concrete application below is said
It is bright:
Method of the invention is applied for Sigsbee 2A model, input observation data, including migration velocity field is (such as
Shown in Fig. 2), observation seismic data, observation system, reflected wave field weighting function, L1 modular constraint threshold function table, and being arranged initial
Imaging results be 0, observation system distribution are as follows: be evenly arranged 160 big guns on sea with 100 meters of intervals, it is equal in seabed by 1600
The geophone station of even arrangement receives, and is divided into 10 meters between geophone station, offset parameter is as follows: migration velocity field transverse direction sampled point be 1600 and
Longitudinal sampled point is 600, and spatial sampling interval is 10 meters, and time sampling interval is 0.6 millisecond, and the sampling time is 5.4 seconds, dominant frequency
It is 15 hertz, L1 modular constraint threshold function table is gradually decreased as 1% by the 5% of imaging results maximum value；Based on current imaging knot
Fruit carries out inverse time inverse migration computation model data, subtracts each other with observation data and obtains the first data residual error；First data are residual
Difference is transformed to common offset trace gather by common-shot-gather, and isolates reflected wave field and diffracted wave field.
Reflected wave field weighting function is set and is equal to 0, inverting is iterated according to the second data residual error at this time, output is final
Imaging results, as shown in Figure 5.From imaging results, using diffraction wave imaging method of the invention, relative to conventional diffracted wave
Reverse-time migration imaging results (as shown in Figure 3) have suppressed imaging noise, and the convergence of diffraction wave energy, amplitude is more balanced, improves
The imaging resolution of the heterogeneous construction in underground.
Reflected wave field weighting function is setIt is iterated inverting according to the second data residual error at this time, is obtained
Imaging results it is as shown in Figure 6.Using imaging method of the invention, relative to conventional least square reverse-time migration imaging results
(as shown in Figure 4) be more clear the structure imagings such as tomography, hole under salt dome, improve to the possessive construction of underground at
As precision.
Fig. 7 is precedent when weighting function is 0 and weighting functionWhen, method of the invention with other two kinds
The wave-number spectrum of imaging method result compares.Wave of the wave-number spectrum of imaging results of the invention relative to conventional diffracted wave imaging results
Number spectrum, enhances high wavenumber information, compensates for the imaging energy of heterogeneous construction.It is conventional minimum since diffraction wave energy is weaker
Two multiply in reverse-time migration methods and results only middle lower wave number information rich in, lack high wavenumber information, i.e., in its imaging results
The imaging energy of heterogeneous construction is weaker.The wave-number spectrum of all-wave field imaging results of the invention is relative to the conventional least square inverse time
The wave-number spectrum of migration imaging result, frequency band is wider, and high wave number energy is stronger, to underground possessive construction all accurate imagings.
Fig. 8 is in precedent when weighting function is 0 and weighting functionWhen, with conventional least square reverse-time migration
The normalization data residual error convergence curve of methods and results compares.The data residual error of conventional least square reverse-time migration method is with repeatedly
The increase convergence rate of generation number is slack-off, and the second data residual error of method of the invention remains to continue quickly to receive with iteration increase
It holds back.Therefore, the method for the present invention is opposite and art methods, can improve the heterogeneous construction such as hole, seam, hole and tomography of underground
Imaging precision, while can to underground possessive construction high-precision be imaged, be of great significance to the exploration of carbonate reservoir.
Claims (3)
1. a kind of inversion method of subsurface structure, which is characterized in that the described method includes:
The first data residual error is calculated according to model data and observation data, obtains the diffracted wave field in the first data residual error and reflection
Wave field；
According to preset reflected wave field weighting function and the reflected wave field, diffracted wave field, the second data residual error is determined, it is specific to wrap
It includes:Wherein,For the second data residual error, D(k)For diffracted wave field,Letter is weighted for reflected wave field
Number, R(k)For reflected wave field, k is the number of iterations, wherein the reflected wave field weighting function is equal to 0, or, the reflected wave field adds
Weight function is greater than 0 and less than 1, and, increase with the number of iterations and successively decreases；
According to the second data residual error, it is iterated inversion imaging.
2. the inversion method of subsurface structure as described in claim 1 calculates the first data according to model data and observation data
Residual error obtains diffracted wave field and reflected wave field in the first data residual error, comprising:
Obtain model data and observe data difference, wherein the model data be according to observation data imaging results and
Determined by the number of iterations；
Isolate diffracted wave field included in the difference；
Reflected wave field is determined according to the difference and diffracted wave field.
3. the inversion method of subsurface structure as described in claim 1, if the second data residual error is determined, according to second data
Residual error is iterated inversion imaging, comprising:
According to following objective function determined by the second data residual error:Carry out inverting at
Picture；Wherein, J (m) is objective function, and m is imaging results, TλFor threshold function table.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710554305.0A CN107255832B (en) | 2017-07-10 | 2017-07-10 | A kind of inversion method of subsurface structure |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710554305.0A CN107255832B (en) | 2017-07-10 | 2017-07-10 | A kind of inversion method of subsurface structure |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107255832A CN107255832A (en) | 2017-10-17 |
| CN107255832B true CN107255832B (en) | 2019-02-26 |
Family
ID=60026350
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710554305.0A Expired - Fee Related CN107255832B (en) | 2017-07-10 | 2017-07-10 | A kind of inversion method of subsurface structure |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107255832B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109490951B (en) * | 2018-12-03 | 2019-11-15 | 中国矿业大学(北京) | Diffraction wave imaging method, device and electronic equipment |
| CN113608261B (en) * | 2021-07-30 | 2021-12-28 | 中国矿业大学(北京) | Diffracted wave imaging method and device and electronic equipment |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102455439A (en) * | 2010-11-02 | 2012-05-16 | 中国石油大学(北京) | Diffracted wave field separation method based on Kirchhoff integral method |
| CN102520444A (en) * | 2011-12-13 | 2012-06-27 | 中国科学院地质与地球物理研究所 | Diffraction wave information extraction method in post-stack seismic wave |
-
2017
- 2017-07-10 CN CN201710554305.0A patent/CN107255832B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102455439A (en) * | 2010-11-02 | 2012-05-16 | 中国石油大学(北京) | Diffracted wave field separation method based on Kirchhoff integral method |
| CN102520444A (en) * | 2011-12-13 | 2012-06-27 | 中国科学院地质与地球物理研究所 | Diffraction wave information extraction method in post-stack seismic wave |
Non-Patent Citations (1)
| Title |
|---|
| 基于平面波记录的绕射目标成像方法研究;孔雪 等;《石油地球物理勘探》;20120831;第47卷(第4期);第674-681页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107255832A (en) | 2017-10-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109425896B (en) | Dolomite oil and gas reservoir distribution prediction method and device | |
| CN106154323B (en) | The thin method for predicting reservoir of phased stochastic inverse of frequency processing is opened up based on earthquake | |
| CN104597490B (en) | Multi-wave AVO reservoir elastic parameter inversion method based on accurate Zoeppritz equations | |
| CN104570125B (en) | A kind of method utilizing well data to improve image taking speed model accuracy | |
| CA2940406C (en) | Characterizing a physical structure using a multidimensional noise model to attenuate noise data | |
| CN108369289A (en) | The method for analyzing design geophysical reconnaissance using full wave field inversion point spread function | |
| CN109444956B (en) | Three-dimensional undulation observation surface seismic slope tomography method | |
| CN107065013B (en) | A kind of interval velocity under earthquake scale determines method and device | |
| Wu et al. | A deep learning estimation of the earth resistivity model for the airborne transient electromagnetic observation | |
| CN110007340A (en) | Salt dome speed density estimation method based on the direct envelope inverting of angle domain | |
| CN109633749A (en) | Non-linear Fresnel zone seismic traveltime tomography method based on scattering integral method | |
| CN113945982A (en) | Method and system for removing low frequency and low wave number noise to generate enhanced image | |
| CN107255832B (en) | A kind of inversion method of subsurface structure | |
| CN104977615B (en) | A kind of multiple ripple drawing method of deep water OBC data based on modeling statistics pickup | |
| CN104181586A (en) | Inversion method of waterland detector data seabed reflection coefficient | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| CN106443776B (en) | A kind of submarine seismograph method for relocating based on isochronous surface method | |
| CN109143336B (en) | A method of overcome the period in full waveform inversion to jump | |
| CN104316961A (en) | Method for obtaining geological parameters of weathered layer | |
| EP3749987B1 (en) | Systems and methods to enhance 3-d prestack seismic data based on non-linear beamforming in the cross-spread domain | |
| CN104536041B (en) | Optimization method of seismological observation system parameters | |
| CN108375794A (en) | Based on the VSP fracture hole Diffraction Imaging technical methods symmetrically observed | |
| US11662493B2 (en) | Methodology for enhancing properties of geophysical data with deep learning networks | |
| CN109387872A (en) | Surface-related multiple prediction technique | |
| Nivlet | Low-frequency constrain in a priori model building for stratigraphic inversion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190226 Termination date: 20190710 |