CN106950596A - A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate - Google Patents
A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate Download PDFInfo
- Publication number
- CN106950596A CN106950596A CN201710231308.0A CN201710231308A CN106950596A CN 106950596 A CN106950596 A CN 106950596A CN 201710231308 A CN201710231308 A CN 201710231308A CN 106950596 A CN106950596 A CN 106950596A
- Authority
- CN
- China
- Prior art keywords
- contrast
- source
- wavelet
- finite difference
- initial
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 83
- 238000004364 calculation method Methods 0.000 claims abstract description 16
- 239000011159 matrix material Substances 0.000 claims description 11
- 239000004576 sand Substances 0.000 claims description 6
- 239000013598 vector Substances 0.000 claims description 5
- 239000002131 composite material Substances 0.000 claims 1
- 238000000354 decomposition reaction Methods 0.000 description 8
- 238000004088 simulation Methods 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 3
- 230000036039 immunity Effects 0.000 description 3
- 238000007796 conventional method Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000002360 explosive Substances 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
Source full waveform inversion method is contrasted the invention discloses a kind of finite difference based on wavelet iterative estimate, be the described method comprises the following steps：S1, inputs initial velocity model, builds finite difference operator, and calculate background wave field；S2, calculates approximate focus ratio, thus calculates and obtains scattered field, and obtains contrast source initial value, contrast initial value and total wave field initial value；S3, each big gun of calculating initially contrasts the Grad corresponding to source, builds conjugated gradient direction；S4, asks for the step-length that every big gun contrast source updates, and be iterated renewal to contrast source；S5, judges that the size of target function value now and allowable error determines whether iteration stops；S6, calculates according to contrast value now and obtains spread speed of the ripple in medium.The method of the present invention is effectively overcome when used wavelet differs larger with true wavelet, the difficult point of inverting acquired results substantial deviation true model, and amount of calculation is significantly reduced.
Description
Technical Field
The invention relates to a finite difference contrast source full waveform inversion method based on wavelet iterative estimation, and belongs to the technical field of frequency domain acoustic full waveform inversion.
Background
Full Waveform Inversion (FWI) extracts velocity parameters from the dynamics of seismic data, and is expected to obtain a more accurate subsurface velocity structure than methods based on data kinematics, such as travel time tomography. At present, there are two main categories from the implementation point of view, time domain and frequency domain. The inversion result with higher precision can be obtained only by a plurality of frequencies in the frequency domain, and the multi-scale calculation is easy to realize and is rapidly popularized. The most important calculation amount of the frequency domain full waveform inversion is a large amount of forward modeling, a large linear equation set needs to be solved, and the direct method (such as LU decomposition) has the advantage of calculating multiple shots, so the direct method is widely used, and the main calculation amount is the decomposition operation of a matrix, and is large. Meanwhile, due to the lack of accurate source wavelet information, when the difference between the simulated wavelet and the real wavelet is too large, the inversion result of the simulated wavelet is completely deviated from the real model.
Disclosure of Invention
The invention aims to provide a finite difference contrast source full waveform inversion method based on wavelet iterative estimation, which can solve the problems in the prior art, remarkably reduce the calculated amount and can still invert to obtain a correct velocity model result when a seismic source wavelet cannot be accurately obtained.
In order to solve the technical problems, the invention adopts the following technical scheme: a finite difference contrast source full waveform inversion method based on wavelet iterative estimation comprises the following steps:
s1, constructing a finite difference operator according to the initial velocity model, and calculating a background wave field;
s2, calculating an initial approximate seismic source ratio by using the simulated wave field and the real wave field at the receiving point corresponding to any one shot, thereby obtaining a scattered field, calculating a contrast source initial value according to the scattered field, and calculating a contrast initial value and a total wave field initial value by using the contrast source initial value;
s3, calculating a gradient value corresponding to each shot initial contrast source according to the inversion target function, and constructing a conjugate gradient direction;
s4, calculating the updating step length of each shot of comparison source, and performing iterative updating on the comparison source;
s5, updating the total wave field, the contrast and the approximate seismic source ratio according to the comparison source value after iterative updating, and judging the relation between the objective function value and the allowable error, if the objective function value is larger than the allowable error, the step S3 is carried out by taking the updated comparison source, the contrast, the total wave field and the approximate seismic source ratio as initial values, otherwise, the step S6 is carried out;
and S6, calculating the propagation speed of the wave in the medium according to the contrast value at the moment.
In the foregoing method for inverting full waveforms of finite difference contrast sources based on wavelet iterative estimation, the approximate source ratio λ in step S2 is
r ∈ S, whereinAnd urefAnd (r) column vectors respectively representing simulated and observed wavefields at the receiving point corresponding to the one shot.
In the foregoing method for full waveform inversion of finite difference contrast source based on wavelet iterative estimation, the formula of the fringe field in step S2 is as follows:r∈S。
in the foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimation, in step S2, a back propagation method is used to obtain a contrast source initial value, and the calculation method includes:
in the foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimation, the contrast value calculation method in step S2 includes:
wherein Re represents the real part of the complex number,
uj,nis the reconstructed total wave field
In the foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimation, in step S3, the inversion target function Fn is formed by a data domain and an inversion domain:
wherein the linear operator LbActing on the contrast source to generate a scattered wave field; projection operator MSAnd MDWherein M isSIndicating that the data in the total area T is projected onto the data field S, and MDProjecting the data on the T to an inversion domain D; djRepresenting the received scattered field, ηSAndare respectively shown as followsηSAndthe method plays a role in balancing the data domain and the inversion domain weight.
In the foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimation, the step size of the nth iteration in step S4Let the conjugate gradient direction of the nth iteration be pj,nComparison source wjThe update formula of (2) is:
in the foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimation, the calculation formula of the propagation velocity of the wave calculated in step S6 is as follows:r ∈ D, wherein v (r) and vb(r) is the propagation velocity of the wave in the real model and the background model, respectively.
Compared with the prior art, the invention has the following advantages:
1. when the difference between the used wavelet and the real wavelet is large, the inversion result in the prior art is seriously deviated from a real model, and the inversion effect basically consistent with the use of the real wavelet can be obtained by adopting the method, namely the dependence of the method on the wavelet is low.
2. Compared with the conventional full waveform inversion method, the method provided by the invention has the advantages that the updated object is the contrast source of each shot, and the background model is kept unchanged in the single-frequency inversion process, so that only one time of decomposition of the impedance matrix is needed for one frequency, and the forward simulation calculation amount in the inversion process is obviously reduced.
3. The method adopted by the invention has certain noise resistance; when the signal-to-noise ratio is 30dB, the structure has a discontinuous phenomenon; when the signal-to-noise ratio is 25dB, the speed model is seriously interfered, but the speed model does not deviate from the actual model seriously, and the method also has better stability on data containing random noise.
Drawings
FIG. 1 is a flow chart of one embodiment of the present invention;
FIG. 2 is an initial velocity model in an embodiment of the present invention, where a is an initial velocity model obtained by co-shot recording and b is an initial velocity model obtained by true model smoothing;
FIG. 3 is a wavelet shape and amplitude spectrum of an embodiment of the present invention, wherein (a) a real wavelet; (b) simulating a wavelet; (c) an amplitude spectrum of the real wavelet; (d) simulating an amplitude spectrum of the wavelet;
FIG. 4 is a schematic representation of a 24Hz single-shot real wavefield and the real part of the reconstructed wavefield in an embodiment of the present invention, wherein (a) the real wavefield; (b) a reconstructed wavefield for the KSS method; (c) a reconstructed wavefield of the IES method;
FIG. 5 is a real part of a 24Hz contrast source and contrast for an embodiment of the present invention, wherein (a) the single shot contrast source of the KSS method; (b) a single shot contrast source for the IES method; (c) contrast for the KSS method; (d) contrast of the IES method;
FIG. 6 is a final inversion result in an embodiment of the present invention, wherein (a) the KSS method; (b) an IES method; (c) analog wavelets are used and no processing is done;
fig. 7 shows a comparison of a single speed profile of the present invention with (a) x being 1.5 km; (b) x is 3.0 km;
FIG. 8 is a graph of relative error curves for wavelets at various frequencies in an embodiment of the present invention;
FIG. 9 is a chart of shot records containing 35dB of random noise in an embodiment of the present invention;
FIG. 10 shows the inversion results for different SNR for the example of the present invention, (a)40 dB; (b)35 dB; (c)30 dB; (d)25 dB;
FIG. 11 shows an example of a gradient initial model (a) and an inversion result (b).
FIG. 12 is a flow chart of an embodiment of the present invention.
The invention is further described with reference to the following figures and detailed description.
Detailed Description
In practical situations, a common shot gather record (time domain, observation data) acquired in the field is used for obtaining an observation wave field of each frequency at the position of the detector by utilizing Fourier transform, and the wavelet excited by the explosive source is the real wavelet. In the inversion, a simulation wavelet is also needed to generate a simulation wave field, but the real wavelet is not known generally, so that inversion failure is easily caused when the difference between the artificially set simulation wavelet and the real wavelet is large, and the problem is solved by using a wavelet iteration method.
One embodiment of the invention: a finite difference contrast source full waveform inversion method based on wavelet iterative estimation is disclosed, as shown in FIG. 1, and comprises the following steps:
s1, constructing a finite difference operator according to the initial velocity model, and calculating a background wave field; in this embodiment, the grid size is 101 × 251, and the horizontal and vertical spatial sampling intervals are all 15 m. The seismic source and the detector are located on the earth surface, the first shot is located at 60m, the shot distance is 75m, the number of shots is 49, the detector is set to be fully received on the earth surface, so that a common shot set record is obtained, an initial velocity model is obtained through velocity analysis, and then Fourier transform is performed on the common shot set record to obtain an observation wave field of each frequency at the position of the detector as shown in figure 2 a. The observed wavefield is the true wavefield received at the location of the geophone. In seismic exploration, only the wavefield at the receiver location is obtained, and the true wavefield is unknown in most parts of the subsurface.
In an isotropic medium, the propagation of seismic waves satisfies the frequency-space domain Helmholtz equation
Wherein u isj(r) representation by the true seismic source f1The total wave field that is generated is,representing the spatial location of the source, and k (r) is the wavenumber at spatial location r. Impedance matrix H for constructing background model by adopting mixed 9-point differential formatbAnd further calculate the linear operatorkbRepresenting the wave number of the background model, from which k is calculatedbFurther obtain the background wave fieldWherein the initial velocity model is used as background model in the first inversion, and the velocity obtained by the inversion of the previous frequency in the later inversionThe model serves as a background model.
S2, calculating initial approximate source ratio by using the simulated wave field and the observed wave field at the receiving point corresponding to any one of the cannons, thereby obtaining a scattered field, calculating a contrast source initial value according to the scattered field, calculating a contrast initial value and a total wave field initial value by using the contrast source initial value, wherein the simulated wave field can be obtained by the input simulated wavelet.
In actual seismic exploration, the real wavelet is usually difficult to obtain accurately, and the real wavelet f1And analog wavelet f2Are not equal, by f2Generated background wavefieldThe scatter field cannot be obtained directly from the total wave field, for which purpose an approximate source ratio λ is defined as
Wherein in the formulaA column vector representing the simulated wavefield at the receiver corresponding to one of the shots, obtained from the forward modeling, urefAnd (r) is a column vector formed by the observed wave fields, an asterisk "+" represents a conjugate transpose, r represents a space position, S is an area where the detector is located, and lambda obtained by calculation before iteration starts is an initial approximate seismic source ratio.
When the background model is close to the real model, λ ≈ f2/f1. Therefore, the scattered field d received at this timejIs composed of
Where j denotes the shot location, uj(r) represents the observed wavefield and,representing the background wavefield.
At the beginning of the inversion at each frequency, we also need to provide an initial contrast source for each shot, i.e., the initial value of the contrast source. The initial value of the contrast source is obtained by using a back propagation method:
wherein, the contrast source wj(r) from the contrast χ (r) and the total wave field uj(r) constitution
wj(r)＝χ(r)uj(r) (5)
The contrast is expressed as follows
Wherein v (r) and vb(r) is the propagation velocity of the wave in the real case and background models, respectively.
Then the initial contrast value is obtained by using the following contrast updating formula (namely when n is equal to 0),
wherein Re represents the real part of the complex number, uj,nThe reconstructed total wave field is obtained, when n is 0, the initial value of the total wave field is obtained,
the initial contrast source calculated by equation (4) and the contrast source calculated by equation (8)The total wavefield is equivalent to the sum of the simulated seismic sources f2And (4) generating. The initial approximate source ratio only provides a substantially correct direction for the inversion, and after model iteration, λ needs to be recalculated to continue model updating, and the method adopts the method after model iterationLambda is recalculated.
S3, under the framework of the contrast source inversion method, the target region of the exploration is called D, the region where the detector is located is called a data field S, and the two regions form a total region T.
The inverted target function is composed of a data domain and an inversion domain
In the above formula, the linear operator LbActing on the contrast source to generate a scattered wave field; projection operator MSAnd MDWherein M isSIndicating that the data in the total area T is projected onto the data field S, and MDProjecting the data on the T to an inversion domain D; djRepresenting the received scattered field, ηSAndhas the function of balancing the weight of the data domain and the inversion domain to ensure that the inversion is stably and effectively carried out, and are respectively expressed as follows
It can be seen from the objective function that the contrast source inversion method is very different from the conventional full waveform inversion, and after a frequency is determined, the parameters to be updated are the contrast source and the contrast, and the update of the velocity parameter is not involved. That is, only one LU decomposition is required to be performed on the impedance matrix after a background model is selected, and the forward calculation in all iterations at the current frequency can use the result of the first decomposition, thereby significantly reducing the calculation amount and improving the inversion efficiency.
The objective function contains two kinds of variables (x, w)j) First, let χ be constant, for wjAnd (6) updating. To simplify the representation, the data field error for the nth iteration is defined as ρj,nAnd the inversion domain error of the nth iteration is rj,nThen the objective function is paired with wjDerivation is carried out to obtain each cannon wjUpdate the required gradient:
in the above formula, the upper dash line represents taking complex conjugate,represents LbAssociated operator of MS*Representing the projection of the matrices in the data field into the total area T, MD*Representing the projection of the matrix in the inversion domain to the total area T; it is noted that the result of the decomposition of the impedance matrix can also be used to construct an operatorAnd does not require re-decomposition.
In this embodiment, after the gradient direction is obtained by the conjugate gradient method, the Polak-Ribiere direction is constructed as the conjugate gradient direction Pj(i.e., CG direction P)j)。
S4, calculating the updating step length of each shot of comparison source, and performing iterative updating on the comparison source;
is the step size of the nth iteration, orderThen
Let the conjugate gradient direction of the nth iteration be pj,nThen w isjIs as follows
S5, updating the total wave field by a formula (8) according to the comparison source value after iterative updating; the selection criterion of contrast is to minimize the overall objective function, but onlyIncluding contrast, byMolecular item of (1)Minimizing to 0 to obtain an updated formula (7) of the contrast, updating the approximate seismic source ratio according to the formula (2), judging the relation between the objective function value and the allowable error at the moment after updating the contrast and the approximate seismic source ratio, if the objective function value is greater than the allowable error, taking the updated contrast source, the contrast, the total wave field and the approximate seismic source ratio as initial values to step S3, otherwise, entering step S6, wherein the allowable error as an empirical value can be set according to different conditions.
S6, calculating the propagation velocity of the wave in the medium according to the contrast value at the moment, wherein the calculation formula is as follows:
wherein v (r) and vb(r) is the propagation velocity of the wave in the medium and in the background model, respectively.
When the difference between the used wavelet and the real wavelet is large, the inversion result in the prior art is seriously deviated from the real model. The method of the invention provides an updated direction for inversion by using the initial approximate seismic source ratio, so that the inversion is carried out towards the basically correct direction, and the generated simulated wave field is closer to the real situation after the model is updated, thereby the approximate seismic source ratio is also closer to the real seismic source ratio. By approximating the source ratio more accurately, the model obtained by inversion can be more accurate, thereby realizing a benign cycle. Even under the condition that wavelets are different greatly, a correct inversion result can be obtained, and therefore dependency on the wavelets is reduced.
Example 2
As another embodiment of the present invention, the inventors further verified the method of the present invention by applying the scheme of the present invention to the overhrust model (i.e. the inferred volume model), as shown in FIG. 12, where an initial smooth velocity model (as in FIG. 2b) is smoothed from the overhrust model while the first derivative of the Gaussian function is used as the true wavelet (FIG. 3a) for generating the observed wavefield, and FIG. 3c is its corresponding amplitude spectrum. The simulated wavelet generating the background wave field in the inversion is a Ricker wavelet with the dominant frequency of 10Hz (figure 3b), the frequency band of the real wavelet is obviously wider than that of the simulated wavelet from the amplitude spectrum (figure 3d), the contrast source inversion method with larger amplitude difference of the real wavelet and the simulated wavelet is different from the conventional full waveform inversion, the updated object is the contrast source of each cannon, then the scattered field is updated according to the contrast source, and the higher the coincidence degree of the reconstructed wave field synthesized by the scattered field and the background field and the real wave field is, the better the inversion result is. Constructing a finite difference operator by inputting an initial velocity model, and calculating a background wave field; then, calculating an approximate seismic source ratio, thereby obtaining a scattering field through calculation, and obtaining a contrast source initial value, a contrast initial value and a total wave field initial value; calculating a gradient value corresponding to each shot initial contrast source, and constructing a conjugate gradient direction; further solving the updating step length of each shot of the comparison source, and carrying out iterative updating on the comparison source; determining whether the iteration is stopped or not according to the objective function value and the allowable error; and finally, calculating the propagation speed of the wave in the medium according to the contrast value when iteration stops. And the obtained inversion result is compared with the inversion method commonly used at present, and the detailed comparison scheme is as follows:
three different inversion methods were applied to the overhrust model: 1) using real wavelets for inversion, denoted as KSS (knock source signature); 2) selecting analog wavelets, and using a finite difference contrast source full waveform inversion strategy based on wavelet iterative estimation, and recording the strategy as IES (iterative estimation of source signature), namely the scheme of the invention; 3) analog wavelets were selected and no processing was performed. The frequencies used are all from 3Hz to 30Hz, inversion is carried out on 10 frequencies with 3Hz as an interval, and the inversion result of the low frequency is used as the input of the adjacent high frequency.
FIG. 4a is the real part of the true wavefield for the 25 th shot at 24Hz, and FIGS. 4b and 4c are the real parts of the reconstructed wavefields for the KSS and IES methods, respectively, noting that the reconstructed wavefield for the IES method is obtained by dividing the updated total wavefield by the approximate source ratio, unlike when known wavefields are used. It is not difficult to find that the reconstructed wavefields of both methods are very close to the true wavefield, which indicates that the inversion is proceeding smoothly. FIG. 5(a-b) is a comparison of the real source parts of the same cannon for the KSS and IES methods at 24Hz, which are irregular in shape, but still show traces of the wavefield of FIG. 4. When the contrast sources for all the shots and the reconstructed wavefield are obtained, the contrast shown in FIG. 5(c-d) can be obtained (FIG. 5 c: KKS; FIG. 5 d: IES).
Fig. 6(a-c) are inversion results of three cases, and the inversion result is the best when the wavelet is correct (fig. 6a), which can accurately reflect the construction information in the real model, and this also shows the effectiveness of the contrast source inversion method. FIG. 6b is the inversion result of the IES method, which is substantially the same as the inversion result when the wavelets are correct, but the convergence is slightly worse at the deep position. Because the difference between the real wavelet and the simulated wavelet is too large, when an inversion strategy of wavelet iterative estimation is not used, the inversion result (figure 6c) completely deviates from a real model, and the importance of the wavelet in full waveform inversion can be seen.
The speed profiles at 1.5km and 3.0km were extracted from each of FIGS. 6(a-b), and a comparison of the single-pass speed profiles as shown in FIGS. 7(a-b) was plotted. The velocity profiles of the two methods substantially coincide and are well matched to the true velocity, especially at levels where the structural changes are not very drastic. We define the relative error function of the wavelets
Wherein λ istrueIs the ratio of the simulated wavelet to the real wavelet. FIG. 8 is a plot of the relative error for the initial and final wavelets for each frequency. From the initial relative error, it can be seen that the approximate source ratio defined by equation (7) is relatively close to the true source ratio, which is also a prerequisite for proper inversion by the IES method. As the iteration is carried out, the model is closer to the real situation, so that the wavelet is updated, and the relative error of the wavelet is smaller and smaller.
In addition, the inventor also performed other verification works:
1. the noise immunity of the method is verified, random white noise is added into shot records by using an AWGN function in Matlab, synthetic seismic data with signal-to-noise ratios of 40dB, 35dB, 30dB and 25dB are respectively obtained, and FIG. 9 shows that a single shot record with the signal-to-noise ratio of 35dB is obtained, wherein a same-phase axis with a small deep amplitude is basically annihilated by the noise.
Inversion results of the IES method are tested by utilizing the noisy data, wherein the inversion results (shown in FIGS. 10a-b) at 40dB and 35dB basically can restore the real situation except that the convergence of the high-speed fault position in the middle is slightly poor, and the IES method has good noise immunity. When the signal-to-noise ratio is 30dB, due to the influence of noise, a discontinuous phenomenon appears in a plurality of structures (FIG. 10c), but the inversion result still does not deviate from the actual model in general. When the signal-to-noise ratio is 25dB, the phenomenon of constructing discontinuity is serious, but compared with the inversion result in FIG. 6c, FIG. 10d does not fall into a local minimum value seriously, and the inversion strategy embodying wavelet iterative estimation is not sensitive to random noise.
2. And (5) verifying the adaptability of the method to the initial model. The initial models selected in the above test are all smooth models in fig. 2b, the full waveform inversion has strong dependence on the initial models, and particularly when an inversion strategy of wavelet iterative estimation is used, the initial models not only affect the nonlinearity degree of the inversion problem, but also determine the accuracy of the initial approximate source ratio. Next we used a more nonlinear gradient model as the initial model (FIG. 11a), with a top velocity of 2700m/s and a bottom velocity of 6000m/s, with uniform variation in the middle. The inversion result of the IES method is shown in fig. 11b, and the inversion effect is substantially the same as that in fig. 6b, which not only reflects that the dependence of the contrast source inversion method on the initial model is relatively low, but also reflects that the inversion strategy of wavelet iterative estimation still has good stability when the initial model is poor.
On the basis of finite difference contrast source full waveform inversion, the invention provides a frequency domain full waveform inversion method with high iteration efficiency and small dependence on the wavelet by introducing an inversion strategy of wavelet iteration estimation, and the method has the following advantages: (1) compared with the conventional method, the method has the advantages that the updated object is the contrast source of each shot, and the background model is kept unchanged in the single-frequency inversion process, so that the decomposition of the impedance matrix is only needed for one frequency, and the forward simulation calculation amount in the inversion process is obviously reduced; (2) when the difference between the simulated wavelet and the real wavelet used in the conventional technology is large, the inversion result obtained by the conventional method deviates from the real model seriously, and the inversion strategy of the invention can obtain the inversion effect basically consistent with that of the real wavelet. (3) By comparing the inversion results of the synthetic data with different signal-to-noise ratios, the method finds that: the influence on the inversion result is small when the signal-to-noise ratio is 40dB and 35dB, and the method is verified to have certain noise immunity; at 30dB, the structure is discontinuous; when the signal-to-noise ratio is 25dB, the speed model is seriously interfered, but the speed model does not deviate from the actual model seriously, which also shows that the method has better stability on data containing random noise; (4) the inversion results of the gradient initial model and the smooth initial model are basically consistent, namely the method has low dependency on the initial model.
Claims (8)
1. A finite difference contrast source full waveform inversion method based on wavelet iterative estimation is characterized by comprising the following steps:
s1, constructing a finite difference operator according to the initial velocity model, and calculating a background wave field;
s2, calculating an initial approximate seismic source ratio by using the simulated wave field and the observed wave field at the receiving point corresponding to any one shot, thereby obtaining a scattered field by calculation, calculating a contrast source initial value according to the scattered field, and obtaining a contrast initial value and a total wave field initial value by calculation of the contrast source initial value;
s3, calculating a gradient value corresponding to each shot initial contrast source according to the inversion target function, and constructing a conjugate gradient direction;
s4, calculating the updating step length of each shot of comparison source, and performing iterative updating on the comparison source;
s5, updating the total wave field, the contrast and the approximate seismic source ratio according to the comparison source value after iterative updating, and judging the relation between the objective function value and the allowable error, if the objective function value is larger than the allowable error, the step S3 is carried out by taking the updated comparison source, the contrast, the total wave field and the approximate seismic source ratio as initial values, otherwise, the step S6 is carried out;
and S6, calculating the propagation speed of the wave in the medium according to the contrast value at the moment.
2. The wavelet iterative estimation-based finite difference contrast source full waveform inversion method according to claim 1, wherein the approximate source ratio λ in the step S2 is
r ∈ S, whereinColumn vectors, u, representing the components of the simulated wavefield at the receiver point corresponding to a particular one of the shotsrefAnd (r) is a column vector formed by the observation wave field of a certain cannon, an asterisk ". indicates a conjugate transpose, r indicates a space position, and S is the area where the detector is located.
3. The wavelet iterative estimation-based finite difference contrast source full waveform inversion method of claim 2, wherein the scatterfield d is calculated in step S2 by using the following methodj(r)：r ∈ S, where j denotes the shot location, uj(r) representing the observed wavefield，Representing the background wavefield.
4. The wavelet iterative estimation-based finite difference contrast source full waveform inversion method of claim 3, wherein in step S2, a contrast source w is obtained by using a back propagation methodjThe initial values are specifically calculated as follows:
wherein,representing a linear operator LbAssociated operator of MSRepresenting the projection of the matrix in the total area T onto the data field S, MS*Representing the projection of the matrix in the data field S to the total area T; l isbFrom the background model impedance matrix HbInverse of (2) and background wave number kbThe components are combined to form the composite material,
5. the wavelet iterative estimation-based finite difference contrast source full waveform inversion method of claim 3, wherein in step S2, when n is 0, the initial value of the total wavefield is calculatedAnd contrast initial valueWhere Re represents the real part of the complex number, MDRepresenting the projection of the matrix in the total area T onto the inversion area D, the complex conjugate is shown by the upper line, and n refers to the number of iterations.
6. The wavelet iterative estimation-based finite difference contrast source full waveform inversion method of claim 5Method, characterized in that in said step S3, the inverse objective function Fn is formed from a data domain functional FSSum-inversion domain functionalJointly constitute:wherein the linear operator LbActing on the contrast source to generate a scattered wave field; djRepresenting the received scattered field, ηSAndrespectively, as follows:
7. the wavelet iterative estimation-based finite difference contrast source full waveform inversion method of claim 6, wherein the step size of the nth iteration in the step S4Wherein the conjugate gradient direction of the nth iteration is pj,n，And (3) updating a contrast source by adopting the following method:
8. the finite difference contrast source full waveform inversion method based on wavelet iterative estimation as claimed in claim 7, wherein said step ofIn step S6, the method for calculating the propagation velocity of the wave in the medium is:r ∈ D, wherein v (r) and vb(r) is the propagation velocity of the wave in the medium and in the background model, respectively.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231308.0A CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231308.0A CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106950596A true CN106950596A (en) | 2017-07-14 |
| CN106950596B CN106950596B (en) | 2019-02-26 |
Family
ID=59474578
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710231308.0A Expired - Fee Related CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106950596B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108508482A (en) * | 2018-05-25 | 2018-09-07 | 中国海洋石油集团有限公司 | A kind of subterranean fracture seismic scattering response characteristic analogy method |
| CN109239771A (en) * | 2018-08-10 | 2019-01-18 | 杭州电子科技大学 | A kind of elastic wave imaging method based on non-homogeneous background media |
| CN109541681A (en) * | 2017-09-22 | 2019-03-29 | 中国海洋大学 | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate |
| CN112444870A (en) * | 2019-08-30 | 2021-03-05 | 中国石油化工股份有限公司 | Multi-scale finite difference contrast source inversion method based on backscattering theory |
Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103091711A (en) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | Method and device for full-wave-shape inversion |
| CN103135132A (en) * | 2013-01-15 | 2013-06-05 | 中国科学院地质与地球物理研究所 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
| CN103713315A (en) * | 2012-09-28 | 2014-04-09 | 中国石油化工股份有限公司 | Seismic anisotropy parameter full waveform inversion method and device |
| US20150012221A1 (en) * | 2013-07-08 | 2015-01-08 | Reeshidev Bansal | Full-Wavefield Inversion of Primaries and Multiples in Marine Environment |
| CN104570082A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on green function characterization |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | Efficient time domain full waveform inversion method |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod and system based on gradient processing |
| CN105467444A (en) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | An elastic wave full-waveform inversion method and apparatus |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
| CN105676277A (en) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | Full waveform joint inversion method for improving high steep structure velocity inversion efficiency |
| US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
| CN106501852A (en) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | A kind of multiple dimensioned full waveform inversion method of three-dimensional acoustic wave equation arbitrarily-shaped domain and device |
| CN106526674A (en) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | Three-dimensional full waveform inversion energy weighted gradient preprocessing method |
-
2017
- 2017-04-11 CN CN201710231308.0A patent/CN106950596B/en not_active Expired - Fee Related
Patent Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103713315A (en) * | 2012-09-28 | 2014-04-09 | 中国石油化工股份有限公司 | Seismic anisotropy parameter full waveform inversion method and device |
| CN103135132A (en) * | 2013-01-15 | 2013-06-05 | 中国科学院地质与地球物理研究所 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
| CN103091711A (en) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | Method and device for full-wave-shape inversion |
| US20150012221A1 (en) * | 2013-07-08 | 2015-01-08 | Reeshidev Bansal | Full-Wavefield Inversion of Primaries and Multiples in Marine Environment |
| CN104570082A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on green function characterization |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | Efficient time domain full waveform inversion method |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod and system based on gradient processing |
| US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105467444A (en) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | An elastic wave full-waveform inversion method and apparatus |
| CN105676277A (en) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | Full waveform joint inversion method for improving high steep structure velocity inversion efficiency |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
| CN106501852A (en) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | A kind of multiple dimensioned full waveform inversion method of three-dimensional acoustic wave equation arbitrarily-shaped domain and device |
| CN106526674A (en) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | Three-dimensional full waveform inversion energy weighted gradient preprocessing method |
Non-Patent Citations (4)
| Title |
|---|
| A ABUBAKAR 等: ""A finite-difference contrast source inversion method"", 《INVERSION PROBLEMS》 * |
| 李庆洋 等: ""不依赖子波、基于包络的WFI初始模型建立方法研究"", 《地球物理学报》 * |
| 段晓亮 等: ""基于逆散射理论的地震波速度正则化反演"", 《物理学报》 * |
| 邵婕 等: ""地震波散射理论及应用研究进展"", 《地球物理学进展》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109541681A (en) * | 2017-09-22 | 2019-03-29 | 中国海洋大学 | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate |
| CN108508482A (en) * | 2018-05-25 | 2018-09-07 | 中国海洋石油集团有限公司 | A kind of subterranean fracture seismic scattering response characteristic analogy method |
| CN109239771A (en) * | 2018-08-10 | 2019-01-18 | 杭州电子科技大学 | A kind of elastic wave imaging method based on non-homogeneous background media |
| CN109239771B (en) * | 2018-08-10 | 2020-01-31 | 杭州电子科技大学 | elastic wave imaging method based on non-uniform background medium |
| CN112444870A (en) * | 2019-08-30 | 2021-03-05 | 中国石油化工股份有限公司 | Multi-scale finite difference contrast source inversion method based on backscattering theory |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106950596B (en) | 2019-02-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104597490B (en) | Multi-wave AVO reservoir elastic parameter inversion method based on accurate Zoeppritz equations | |
| US9910174B2 (en) | Seismic imaging apparatus and method for performing iterative application of direct waveform inversion | |
| CN108873066B (en) | Elastic medium wave equation reflected wave travel time inversion method | |
| US10520619B2 (en) | FWI model domain angle stacks with amplitude preservation | |
| CN106950596A (en) | A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate | |
| US20110131020A1 (en) | Dip guided full waveform inversion | |
| KR20090075843A (en) | Iterative inversion of data from simultaneous geophysical sources | |
| US20150293246A1 (en) | Frequency-domain augmented time-domain full wavefield inversion | |
| CN110187384B (en) | Bayes time-lapse seismic difference inversion method and device | |
| CN106526674A (en) | Three-dimensional full waveform inversion energy weighted gradient preprocessing method | |
| CN106569262A (en) | Background speed model reconstructing method in absence of low frequency earthquake data | |
| KR101820850B1 (en) | Seismic imaging apparatus and method using iterative direct waveform inversion | |
| CN111077567B (en) | Method for double-pass wave prestack depth migration based on matrix multiplication | |
| US11340366B2 (en) | Accurate velocity model estimation and imaging in the presence of localized attenuation (Q) anomalies | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| Jun et al. | Weighted pseudo-Hessian for frequency-domain elastic full waveform inversion | |
| CN108680968B (en) | Evaluation method and device for seismic exploration data acquisition observation system in complex structural area | |
| CN107479091A (en) | A kind of method for extracting reverse-time migration angle gathers | |
| CN105607123A (en) | Method and device for calculating seismic wave characteristic information of random pore medium model | |
| CN111077566B (en) | Method for double-pass wave prestack depth migration based on matrix decomposition | |
| CN108680957A (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting | |
| Re | Building the initial velocity model for refraction tomography: A data-driven approach | |
| Zhang et al. | Got noisy impedance? Try this | |
| CN111665549A (en) | Inversion method of stratum acoustic wave attenuation factor | |
| Zhao et al. | Double plane wave reverse time migration in the time domain |
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: 20200411 |