CN115201914B - Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction - Google Patents
Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction Download PDFInfo
- Publication number
- CN115201914B CN115201914B CN202210784961.0A CN202210784961A CN115201914B CN 115201914 B CN115201914 B CN 115201914B CN 202210784961 A CN202210784961 A CN 202210784961A CN 115201914 B CN115201914 B CN 115201914B
- Authority
- CN
- China
- Prior art keywords
- data
- matching
- time
- waveform inversion
- convolution
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000000605 extraction Methods 0.000 title claims abstract description 21
- 238000005070 sampling Methods 0.000 claims abstract description 64
- 238000004088 simulation Methods 0.000 claims abstract description 14
- 238000004364 calculation method Methods 0.000 claims description 19
- 239000013598 vector Substances 0.000 claims description 14
- 238000005457 optimization Methods 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 10
- 238000012937 correction Methods 0.000 claims description 9
- 230000003068 static effect Effects 0.000 claims description 6
- 230000001427 coherent effect Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 19
- 239000011159 matrix material Substances 0.000 description 11
- 230000008569 process Effects 0.000 description 8
- 230000001965 increasing effect Effects 0.000 description 6
- 230000036039 immunity Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 230000002829 reductive effect Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000002238 attenuated effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- 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
-
- 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/34—Displaying seismic recordings or visualisation of seismic data or attributes
- G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/53—Statics correction, e.g. weathering layer or transformation to a datum
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
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)
Abstract
The invention provides a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction, which is characterized in that convolution operation is carried out on single-channel observation data and simulation data through one-dimensional Gaussian convolution cores with different lengths to extract comprehensive attribute features of wave fields of each sampling point in different time ranges; calculating the Euclidean distance of the characteristics between each sampling point of the single-channel analog data and all sampling points of the corresponding channel of the analog data, obtaining the unique optimal matching point with each analog data sampling point in the observed data, breaking the limit of point-to-point matching in time of conventional full waveform inversion, and realizing dynamic matching; introducing an attenuation factor, and carrying out amplitude exponential attenuation on the analog data subjected to dynamic matching according to the absolute value of the optimal matching time difference; a global cross-correlation objective function is employed. The invention breaks through the limit of the point-to-point subtraction in time of the traditional full waveform inversion, realizes the self-adaptive cross-time matching according to the characteristics of the data, and improves the convergence of the full waveform inversion.
Description
Technical Field
The invention belongs to the technical field of seismic exploration, in particular to a seismic exploration imaging method, and particularly relates to a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction.
Background
Seismic exploration is the only geophysical method with large exploration depth and high resolution, and is recognized as the preferred exploration method technology in the field of deep resource exploration. At present, the demand of China for resources such as petroleum and mineral products is increased year by year, and under the condition that the position and reserve of shallow resources are basically ascertained, the center of gravity of resource exploration is gradually changed to deep complex geological structures. In this context, countries place higher demands on depth and accuracy of seismic exploration. The main function of the seismic exploration is to realize the exploration of the underground geological structure and lithology information after a series of data processing according to the seismic signals acquired by the detectors arranged on the earth surface. Therefore, the distribution of physical parameters of the underground medium obtained through the seismic data is particularly important.
Among the many physical parameters, the propagation velocity information of seismic waves in different media is a key factor in determining the offset imaging quality. At present, underground medium speed modeling methods widely used at home and abroad mainly comprise three types of speed analysis, tomography and full waveform inversion. Wherein velocity analysis and tomography use only travel time information of seismic wave propagation and ignore many other important information, and therefore, the velocity modeling result cannot meet the requirement of high-precision seismic exploration. Full waveform inversion (Full Waveform Inversion, FWI) techniques utilize all information of the pre-stack seismic wavefield, including amplitude, phase, travel time, etc., kinematics and dynamics information, to enable high-precision velocity modeling of subsurface media. In addition to modeling speed parameters, full waveform inversion can also model parameters such as density, quality factors, and anisotropic coefficients with high accuracy.
Lailly (1983) and Tarantola (1984) convert the strong nonlinearity problem of subsurface medium velocity inversion into an optimization problem in the least squares sense. By minimizing the observed data and the simulated data under the constraint of the L2 norm, the initial velocity model is continuously iteratively updated to gradually approximate the simulated data to the observed data, thereby finding the optimal solution for the inversion of the velocity parameters of the underground medium. Since directly calculating the derivative of friechet results in a large calculation amount, full waveform inversion is difficult in case of insufficient calculation power. The gradient of model update can be obtained under the condition that the Frectchet derivative is not directly calculated by the proposal of the accompanying state method, so that the calculation efficiency of full waveform inversion is greatly improved, and the enthusiasm of domestic and foreign scientific researchers on full waveform inversion research is greatly promoted. At present, full waveform inversion is a great topic in the field of exploration geophysics by virtue of the advantage of high-precision modeling of various physical parameters in the underground, and countless scientific researchers at home and abroad perform intensive research on all aspects of full waveform inversion, so that the development of the method is greatly promoted.
Meanwhile, as the full waveform inversion is a strong nonlinear problem, the full waveform inversion is approximately a local optimization problem under the constraint of a traditional objective function, so that the full waveform inversion is very dependent on an accurate initial model and low-frequency components in observed data. However, for seismic data of complex geologic formation areas, it is very difficult to provide a good initial model for full waveform inversion; in addition, since it is difficult to generate a low frequency signal when the source is excited and it is difficult for the detectors to respond to the low frequency information in the seismic signal, the low frequency component is difficult to exist or the energy is weak in the seismic data obtained by actual work. Thus, cycle skip issues follow, which is also a major problem faced by full waveform inversion. In addition, noise immunity and reduced wavelet error effects on inversion are also problems that are currently under investigation and addressed. The traditional full waveform inversion based on the L2 norm objective function is obtained by point-to-point subtraction of observed data and simulated data in the time direction, has high dependence on the amplitude precision of the simulated data, and has poor noise resistance. The correlation-based objective function reduces the dependence of FWI on the accuracy of amplitude information, highlights the matching of phase information with better linearity, solves the problem of interference of amplitude information errors or uneven energy of signals with different offsets on inversion results caused by the seismic acquisition process, and has certain noise immunity.
In order to alleviate the cycle skip phenomenon, a plurality of methods are developed currently, for example, bunks (1995) propose to decompose the seismic data into different frequency scales by using a low-pass filter according to frequency bands, and an inversion strategy from a low frequency band to a high frequency band in a time domain is adopted to enable the FWI objective function to be better converged, so that the occurrence of local extremum is reduced, and the cycle skip problem is relieved. Shin and Cha (2008) change the seismic data into Laplace domain, and data containing abundant low-frequency information can be obtained through e exponential decay. Laplace domain waveform inversion can construct long wavelength components of subsurface media. Zhang et al (2017) propose constructing observation data containing low frequency information by using a sparse blind deconvolution algorithm in combination with artificially designed low frequency wavelets based on a convolution model of seismic data. Wang et al (2016) propose to utilize the dynamic normalization method in the signal recognition field, find out the best matching point between two ways according to the similarity between identical track numbers of data, this method has jumped out the traditional defect that point-to-point subtracts on sampling time based on L2 norm FWI, have succeeded in relieving the cycle skip problem. The frequency expansion of the seismic data can be realized by using compressed sensing and sparse inversion in digital signal processing. However, these methods either rely on assistance from other techniques or have poor computational stability, and have some drawbacks in terms of computational efficiency or inversion results.
Disclosure of Invention
The invention aims to provide a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction, which breaks the limit of point-to-point matching of conventional full waveform inversion on sampling time, utilizes one-dimensional convolution kernel to extract the features of single-channel seismic data to realize dynamic matching of a simulation record and an observation record, solves the problem of cycle skip, improves the convergence of full waveform inversion and the precision of inversion results, and does not obviously increase extra calculation amount.
The invention aims at realizing the following technical scheme:
a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction comprises the following steps:
a. Performing preprocessing such as de-compiling, high Cheng Jing correction, first arrival static correction, earth surface consistency static correction, random noise and coherent noise removal, wavelet estimation, missing seismic channel compensation and the like on an actual seismic observation record;
b. Firstly, a linear incremental initial model is established in an estimated speed range, and time domain full waveform inversion related parameters including a seismic wavelet main frequency f, a model size nz x nx, grid distances dx and dz, a sampling total time T, a time sampling interval dt, a maximum iteration number itermax, an iteration step q of an optimization algorithm, an objective function required precision tol, and a maximum value vmax and a minimum value vmin of model speed estimation are set according to requirements.
C. performing finite difference forward modeling calculation on the initial model by using the estimated wavelets according to an actual observation system to obtain a simulated seismic record;
d. designing a one-dimensional Gaussian convolution kernel:
k g denotes a gaussian convolution kernel, n denotes an element in the convolution kernel, l denotes a length of the convolution kernel, exp ()' denotes an exponential function based on a natural constant e, and n m denotes a center point of the convolution kernel;
e. Extracting wave field characteristics of single-channel seismic data of observation data and simulation data in different convolution kernel length time ranges by using Gaussian convolution kernels of different lengths:
* A convolution operation is represented and is performed, And/>Respectively representing simulation data and observation data of the jth channel of the ith gun, ns and nr respectively representing the total gun number and the total channel number of each gun, T and T respectively representing time variable and maximum sampling time, and/(respectively) calculating the time variable and the maximum sampling timeAnd/>Features of single-pass analog data and observed data extracted by a w-th length gaussian convolution kernel are respectively represented, nw represents the total number of used different-length gaussian convolution kernels, the convolution kernel length l monotonically increases from a first value l 1 to a last value l nw, || represents an absolute value taken on the data, norm (-) represents an operator for normalizing the data, M [ ] represents taking the middle part of the convolution result so that the characteristic time length of the clipped data is the same as that of the original seismic record, namely, feature data in a time period from n m to (t+1-n m) are extracted from the convolution result;
f. constructing a feature vector of each sampling point in the single-channel seismic data:
And/> Respectively representing the characteristic vector of each sampling point;
g. Calculating Euclidean distance between each sampling point of single-channel analog data and all sampling point characteristics corresponding to the observed data, and finding out the observed data sampling point which is optimally matched with each analog data sampling point according to the minimum Euclidean distance, namely realizing dynamic matching:
t s represents a sampling point in the analog data, and t o represents a sampling point in the observation data that best matches the analog data;
h. calculating the best matching time difference:
Δt＝ts-to,
Δt represents the best match time difference;
i. According to the size of the absolute value of the optimal matching time difference, amplitude attenuation is carried out on the analog data sampling points so as to realize constraint on the reliability of dynamic matching, and the greater the absolute value of the optimal matching time difference, the higher the attenuation degree of the analog data sampling points is:
Analog data after dynamic matching and amplitude attenuation are represented, gamma represents an attenuation factor, and the attenuation degree of the analog data amplitude along with the increase of the absolute value of the optimal matching time difference is controlled;
j. Establishing a global cross-correlation full-waveform inversion objective function:
J represents an objective function, v represents the speed information of the underground medium, and the gradient expression obtained by taking the derivative of the speed at two ends of the objective function is as follows:
P f represents the time domain forward wavefield, P b represents the inverse wavefield;
k. And carrying out iterative updating on the velocity model by using an L-BFGS optimization algorithm, firstly, reversing the large-scale structure of the model by using a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction, then carrying out full-frequency-band conventional full waveform inversion by taking the inversion result as an initial model, inverting the detailed structure of the model, and finally obtaining the high-precision underground model.
Compared with the prior art, the invention has the beneficial effects that:
The invention provides a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction, which is simple in principle and easy to program, and the comprehensive attribute features (amplitude, phase, travel time and the like) of wave fields of each sampling point of data in different time ranges are extracted by carrying out convolution operation on single-channel observation data and analog data through one-dimensional Gaussian convolution check with different lengths; the method breaks through the limit of point-to-point subtraction in time of the traditional full-waveform inversion, realizes self-adaptive cross-time matching according to the characteristics of data, effectively reduces the occurrence of cycle skipping under the condition of keeping the calculation efficiency unchanged, is very suitable for actual seismic exploration, and can obviously improve the accuracy of the inversion result of the seismic exploration; the method can be combined with a mixed seismic source coding strategy to further improve the calculation efficiency, and compared with the conventional full waveform inversion, the calculation amount is not obviously increased; the invention utilizes the one-dimensional convolution kernel to extract the single-channel seismic data characteristics to realize the dynamic matching of the simulation data and the observation data, and has the following advantages:
1. Extracting comprehensive characteristics of wave fields of single-channel observation data and simulation data in different time ranges by utilizing a plurality of one-dimensional convolution kernels, finding out an optimal matching point pair according to the minimum Euclidean distance of all characteristics of each sampling point of a corresponding channel of the simulation data and the observation data, breaking through the limit of point-to-point matching of conventional full waveform inversion in the same time, realizing self-adaptive dynamic matching according to the characteristic change of the wave field data, effectively eliminating cycle jump phenomenon caused by waveform matching errors, remarkably enhancing the convergence of full waveform inversion and improving the accuracy of inversion results;
2. the dependence of inversion on amplitude information is reduced by using a global cross-correlation objective function. The amplitudes of partial sampling points of the analog records are attenuated, so that the phase information of each sampling point is more similar to the observation records, but the corresponding amplitude information is not corrected, so that the influence of amplitude errors on inversion is reduced by adopting a global cross-correlation objective function, the weight of a remote offset record is increased by cross-correlation calculation, and the remote offset record often carries information of a large-scale structure, so that the inversion precision of the large-scale structure of the model can be further enhanced, and a more accurate initial model is provided for the subsequent full-band inversion;
3. The convolution kernel used for feature extraction is a Gaussian convolution kernel with strong noise immunity, the feature of the data can still be accurately extracted under the condition that the observed data contains strong noise, and in addition, the cross-correlation-based objective function also has strong noise immunity, so that the time domain dynamic wave field matching full waveform inversion method based on the convolution feature extraction has strong noise immunity, and an accurate inversion result can be obtained under the condition that the observed data contains strong noise.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings that are needed in the embodiments will be briefly described below, it being understood that the following drawings only illustrate some embodiments of the present invention and therefore should not be considered as limiting the scope, and other related drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flowchart of a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction;
FIG. 2a removes the 20Hz dominant frequency Rake wavelet for information below 11 Hz;
FIG. 2b is the wavelet spectrum of FIG. 2 a;
FIG. 3a is a real model diagram;
FIG. 3b is a linearly increasing initial model diagram;
FIG. 4a is a graph of inversion results based on conventional time domain full waveform inversion;
FIG. 4b is a graph of the results of a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction.
Detailed Description
The invention is described in further detail below with reference to the drawings and examples. It is to be understood that the specific embodiments described herein are merely illustrative of the invention and are not limiting thereof. It should be further noted that, for convenience of description, only some, but not all of the structures related to the present invention are shown in the drawings.
It should be noted that: like reference numerals and letters denote like items in the following figures, and thus once an item is defined in one figure, no further definition or explanation thereof is necessary in the following figures. Meanwhile, in the description of the present invention, the terms "first", "second", and the like are used only to distinguish the description, and are not to be construed as indicating or implying relative importance.
The time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction is realized on a MATLAB2016b platform, a program is written under the MATLAB2016b software framework, an installation environment of a MATLAB parallel working library is built according to corresponding parallel computing requirements, and a MATLAB parallel computing tool box (Parallel Computing Toolbox) is installed.
The core is to provide a wave field dynamic matching method based on convolution feature extraction. Because low-frequency information is often lost in the actually acquired seismic signals, and an accurate initial model is difficult to provide for full waveform inversion, the problem of cycle skip caused by direct matching of medium-high frequency signals is the most important reason for influencing the accuracy of full waveform inversion results. The feature of each sampling point is defined as the convolution value when the middle element of the convolution kernel is right against a certain sampling point in the convolution calculation, namely the feature value of the sampling point, so that the length selection of the convolution kernel is odd. And extracting the characteristics of each sampling point in different time window ranges by using a plurality of convolution kernels with different lengths so as to improve the accuracy of dynamic matching. The Euclidean distance of all the characteristics between each sampling point of the single-channel analog data and all the sampling points of the corresponding single-channel observation data is calculated, the sampling point of the observation data which is optimally matched with each sampling point of the analog data can be found out according to the minimum principle of the Euclidean distance of the total sum of the characteristics, and the optimal matching time difference is obtained, and the difference value can be positive or negative or zero during the optimal matching. In order to evaluate the reliability of dynamic matching and restrict unreliable dynamic matching to avoid more cycle skip phenomenon, the invention proposes to exponentially attenuate the amplitude of the analog data after dynamic matching according to the absolute value of the optimal matching time difference, i.e. the larger the absolute value of the optimal matching time difference is, the lower the possibility that the matching is correct is, and the greater the attenuation of the amplitude is required. The invention effectively solves the cycle skip problem on the premise of not obviously increasing the calculated amount.
The invention discloses a time domain dynamic wave field matching full-waveform inversion method based on convolution feature extraction, which comprises the following steps:
a. And (3) performing preprocessing such as de-compiling, high Cheng Jing correction, first arrival static correction, earth surface consistency static correction, random noise and coherent noise removal, wavelet estimation, missing seismic channel compensation and the like on the actual seismic observation record.
B. Firstly, a linear incremental initial model is established in an estimated speed range, and time domain full waveform inversion related parameters including a seismic wavelet main frequency f, a model size nz x nx, grid distances dx and dz, a sampling total time T, a time sampling interval dt, a maximum iteration number itermax, an iteration step q of an optimization algorithm, an objective function required precision tol, and a maximum value vmax and a minimum value vmin of model speed estimation are set according to requirements.
C. and carrying out finite difference forward modeling calculation on the initial model by using the estimated wavelets according to an actual observation system to obtain a simulated seismic record.
D. designing a one-dimensional Gaussian convolution kernel:
k g represents a gaussian convolution kernel. n represents an element in the convolution kernel. l denotes the length of the convolution kernel. exp (-) represents an exponential function based on a natural constant e. n m denotes the center point of the convolution kernel.
E. Extracting wave field characteristics of single-channel seismic data of observation data and simulation data in different convolution kernel length time ranges by using Gaussian convolution kernels of different lengths:
* Representing a convolution operation. And/>Respectively representing simulation data and observation data of the jth channel of the ith gun. ns and nr represent the total number of shots and the total number of shots per shot, respectively. T and T represent the time variable and the maximum sampling time, respectively. /(I)And/>Features of single-pass simulation data and observation data extracted by a w-th length gaussian convolution kernel are represented, respectively. nw represents the total number of gaussian convolution kernels of different lengths used. The convolution kernel length l increases monotonically from the first value l 1 to the last value l nw. |. | represents taking absolute value on the data. Firstly, taking absolute values of single-channel seismic data to avoid that the difference of features extracted by different convolution kernels is not obvious due to the fact that the subtracting interference occurs between positive and negative amplitudes during the convolution feature extraction. norm (-) represents an operator that normalizes data. And normalizing the single-channel observation data and the analog data after taking the absolute value to prevent the disturbance of the abnormal amplitude value on convolution characteristic extraction. M </i > ] represents taking the middle part of the convolution result so that the clipped data feature duration is the same as the original seismic record, i.e., feature data in the time period from n m to (T+1-n m) are extracted from the convolution result. And finally, performing normalization processing for the second time so that the features extracted by the convolution kernels with different lengths have the same weight in dynamic matching.
F. constructing a feature vector of each sampling point in the single-channel seismic data:
And/> Representing the feature vector for each sample point separately.
G. Calculating Euclidean distance between each sampling point of single-channel analog data and all sampling point characteristics corresponding to the observed data, and finding out the observed data sampling point which is optimally matched with each analog data sampling point according to the minimum Euclidean distance, namely realizing dynamic matching:
t s represents a sampling point in the analog data, and t o represents a sampling point in the observation data that best matches the analog data.
H. calculating the best matching time difference:
Δt＝ts-to,
Δt represents the best match time difference.
I. and according to the absolute value of the optimal matching time difference, carrying out amplitude attenuation on the analog data sampling points to realize the constraint on the reliability of dynamic matching. The greater the absolute value of the best match time difference, the higher the degree of attenuation of the analog data sampling point is:
Representing the analog data after dynamic matching and amplitude decay. Gamma denotes the attenuation factor, controlling the degree to which the analog data amplitude decays as the absolute value of the best match moveout increases.
J. Establishing a global cross-correlation full-waveform inversion objective function:
j represents an objective function. v represents velocity information of the subsurface medium. The derivative of speed at both ends of the objective function can be given by the gradient expression:
The simplification can be expressed as:
wherein λ is the companion source, and its expression is:
The time domain full waveform inversion gradient formula based on the accompanying state method can be obtained through a chain law:
Where P f is the time domain forward wavefield, L -1 lambda denotes the inverse wavefield, where L -1 is the inverse operator.
K. the inversion is based on the L-BFGS optimization algorithm to carry out iterative updating, and a proper step length is searched through a Wolfe convergence criterion; the iterative formula is expressed as:
mk+1＝mk-αkHkgk
Wherein H k is an inverse matrix of an approximate Hessian matrix, m k is a model update parameter, g k is a model update gradient, alpha k is a step size obtained by Wolfe linear search, and k represents the number of iterations.
The L-BFGS optimization algorithm has few matrixes to be saved in the iterative calculation process, has small requirement on the memory of a computer, and simultaneously has the following iterative formula for updating the Hessian matrix:
Hk+1＝Vk THkVk+ρksksk T
Vk＝I-ρkyksk T
sk＝mk+1-mk,yk＝gk+1-gk
Wherein V k is an intermediate variable in the calculation process, s k is an intermediate variable in the calculation process, ρ k is an intermediate variable, I represents an identity matrix, H k+1 is calculated according to vector pairs { s k,yk } and H k, and only n vector pairs are stored to implicitly express the inverse matrix of the Hessian matrix. The product of H kgk can be obtained by summing the inner product of the series of vectors between the gradient g k and the vector pair { s k,yk }, and the vector. After each update is completed, the last vector pair will be replaced by the current new vector pair { s k+1,yk+1 }. Thus, the set of vector pairs contains the curvature signal of the last n iterations. In practice, more satisfactory results are obtained when n.gtoreq.5. Initial approximate Hessian matrix of L-BFGS optimization algorithm Each iteration is different. The inverse matrix H k of the approximate Hessian matrix needs to satisfy the following update formula:
the updating direction of the model is realized by the following method:
(1) Let q=g k be the number, Q=q- α iyi, where i=k-1, k-2, … k-n; alpha i is an intermediate variable in the calculation process;
(2) Order the R=r-s i(αi - β), where i=k-n, k-n+1, …, k-1; /(I)Approximating the Hessian matrix for the initial;
(3) The update amount H kgk =r is obtained by the above steps.
The updating quantity H kgk of the model is obtained through the method, and then the proper step length alpha k is obtained through Wolfe linear search for updating iteration. The L-BFGS optimization algorithm effectively overcomes the difficulty of explicitly solving the Hessian matrix and the inverse thereof, has the advantages of high super-linear convergence speed, high calculation efficiency, small occupied memory, high precision and the like, and is suitable for solving the large-scale nonlinear optimization problem.
And judging whether the inversion result meets the set termination condition, namely, if so, stopping iteration, and if not, continuing to calculate until the set termination condition is met, and ending calculation.
In order to prevent positive and negative amplitude offset in the convolution feature extraction process, firstly taking an absolute value of original data, and respectively carrying out normalization processing on features extracted by convolution kernels with different lengths after convolution operation extraction so that the weights of the features extracted by the convolution kernels with different lengths on dynamic matching contributions are the same; calculating the Euclidean distance of the characteristics between each sampling point of the single-channel analog data and all sampling points of the corresponding channel of the analog data, obtaining the unique optimal matching point with each analog data sampling point in the observed data according to the minimum Euclidean distance of the characteristics, breaking the limit of point-to-point matching in time of conventional full waveform inversion, and realizing dynamic matching; in order to evaluate and restrict the accuracy of dynamic matching, an attenuation factor is introduced, and amplitude exponential attenuation is carried out on the analog data after dynamic matching according to the size of the absolute value of the optimal matching time difference, namely the greater the absolute value of the dynamic matching time difference is, the higher the attenuation degree of the corresponding sampling point of the analog data after dynamic matching is; the amplitude of the analog data is attenuated in the process, so that the original carried underground medium information is destroyed; therefore, in order to reduce the dependence of inversion on amplitude information, and highlight the effect of phase information with higher linearity, the invention adopts a global cross-correlation objective function.
Example 1
According to the exploration requirements, parallel Computing Toolbox and MATLAB Distributed Computing Server (R2016 b) are installed under a Windows 10 professional system, and a MATLAB parallel platform is built.
And carrying out full waveform inversion testing on the missing low-frequency information by using Marmousi. A real model (fig. 3 a) and a linearly increasing initial model (fig. 3 b).
Table 1 missing low frequency information full waveform inversion test parameters
| Grid size | Grid distance | Length of measuring line | Longitudinal depth | Frequency bandwidth |
| 69×192 | 12.5m | 2400m | 862.5m | Full frequency band |
The size of the grid of the model is 69 multiplied by 192, the grid distance is 12.5, the transverse distance is 2400m, the longitudinal depth is 862.5m, the seismic wave speed in the model ranges from 1500m/s to 4000m/s, the geophones are arranged on the surface of the model, each grid point is one geophone, the interval between the geophones is 12.5m, and 12 seismic sources are excited in total. The seismic source adopts a Rake wavelet (figure 2 a) which lacks information below 11Hz and has 20Hz main frequency, the sampling interval is 0.001s, the actual sampling total length is 2s, and the frequency range is full-band data.
The total number nw of convolution kernels is 40, the maximum length l nw of convolution kernels is 401, and the attenuation factor gamma is 5.
The inversion results are shown in fig. 4. Fig. 4a shows that the conventional full waveform inversion result is greatly affected by the loss of low frequency information, and obvious cycle skipping occurs at a plurality of positions. For a linear incremental initial model, the initial model is greatly different from a real model, in the waveform matching process, the phase difference between a plurality of observation records and a simulation record is more than half a period due to the lack of low-frequency information, and the cycle skipping occurs, so that the gradient update quantity calculation error of the corresponding position is caused. Along with the increase of iteration times in the inversion process, the error speed update quantity is continuously accumulated, and finally the inversion result is greatly different from the real model. Fig. 4b is a full waveform inversion result based on the method of the present invention. From the figure it can be seen that the inversion results can invert substantially all the configurations of the real model without velocity anomaly regions. After the dynamic data extracted based on convolution characteristics are matched, the matching degree of the simulation record and the observation record is improved, so that the model speed value is updated to the correct direction, and finally, a relatively accurate inversion result is obtained without a speed abnormal region.
Note that the above is only a preferred embodiment of the present invention and the technical principle applied. It will be understood by those skilled in the art that the present invention is not limited to the particular embodiments described herein, but is capable of various obvious changes, rearrangements and substitutions as will now become apparent to those skilled in the art without departing from the scope of the invention. Therefore, while the invention has been described in connection with the above embodiments, the invention is not limited to the embodiments, but may be embodied in many other equivalent forms without departing from the spirit or scope of the invention, which is set forth in the following claims.
Claims (1)
1. The time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction is characterized by comprising the following steps of:
a. performing preprocessing of de-compiling, high Cheng Jing correction, first arrival static correction, earth surface consistency static correction, random noise and coherent noise removal, wavelet estimation and missing seismic channel compensation on an actual seismic observation record;
b. firstly, establishing a linear incremental initial model in an estimated speed range, and setting time domain full waveform inversion related parameters according to requirements, wherein the parameters comprise a main frequency f of a seismic wavelet, a model size nz x nx, grid distances dx and dz, a total sampling time T, a time sampling interval dt, a maximum iteration number itermax, an iteration step length q of an optimization algorithm, a target function required precision tol, and a maximum value vmax and a minimum value vmin of model speed estimation;
c. performing finite difference forward modeling calculation on the initial model by using the estimated wavelets according to an actual observation system to obtain a simulated seismic record;
d. designing a one-dimensional Gaussian convolution kernel:
k g denotes a gaussian convolution kernel, n denotes an element in the convolution kernel, l denotes a length of the convolution kernel, exp ()' denotes an exponential function based on a natural constant e, and n m denotes a center point of the convolution kernel;
e. Extracting wave field characteristics of single-channel seismic data of observation data and simulation data in different convolution kernel length time ranges by using Gaussian convolution kernels of different lengths:
* A convolution operation is represented and is performed, And/>Respectively representing simulation data and observation data of the jth channel of the ith gun, ns and nr respectively representing the total gun number and the total channel number of each gun, T and T respectively representing time variable and maximum sampling time, and/(respectively) calculating the time variable and the maximum sampling timeAnd/>Features of single-pass analog data and observed data extracted by a w-th length gaussian convolution kernel are respectively represented, nw represents the total number of used different-length gaussian convolution kernels, the convolution kernel length l monotonically increases from a first value l 1 to a last value l nw, || represents an absolute value taken on the data, norm (-) represents an operator for normalizing the data, M [ ] represents taking the middle part of the convolution result so that the characteristic time length of the clipped data is the same as that of the original seismic record, namely, feature data in a time period from n m to (t+1-n m) are extracted from the convolution result;
f. constructing a feature vector of each sampling point in the single-channel seismic data:
And/> Respectively representing the characteristic vector of each sampling point;
g. Calculating Euclidean distance between each sampling point of single-channel analog data and all sampling point characteristics corresponding to the observed data, and finding out the observed data sampling point which is optimally matched with each analog data sampling point according to the minimum Euclidean distance, namely realizing dynamic matching:
t s represents a sampling point in the analog data, and t o represents a sampling point in the observation data that best matches the analog data;
h. calculating the best matching time difference:
Δt＝ts-to,
Δt represents the best match time difference;
i. According to the size of the absolute value of the optimal matching time difference, amplitude attenuation is carried out on the analog data sampling points so as to realize constraint on the reliability of dynamic matching, and the greater the absolute value of the optimal matching time difference, the higher the attenuation degree of the analog data sampling points is:
Analog data after dynamic matching and amplitude attenuation are represented, gamma represents an attenuation factor, and the attenuation degree of the analog data amplitude along with the increase of the absolute value of the optimal matching time difference is controlled;
j. Establishing a global cross-correlation full-waveform inversion objective function:
J represents an objective function, v represents the speed information of the underground medium, and the gradient expression obtained by taking the derivative of the speed at two ends of the objective function is as follows:
P f represents the time domain forward wavefield, P b represents the inverse wavefield;
k. And carrying out iterative updating on the velocity model by using an L-BFGS optimization algorithm, firstly, reversing the large-scale structure of the model by using a time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction, then carrying out full-frequency-band conventional full waveform inversion by taking the inversion result as an initial model, inverting the detailed structure of the model, and finally obtaining the high-precision underground model.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210784961.0A CN115201914B (en) | 2022-06-30 | 2022-06-30 | Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210784961.0A CN115201914B (en) | 2022-06-30 | 2022-06-30 | Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115201914A CN115201914A (en) | 2022-10-18 |
| CN115201914B true CN115201914B (en) | 2024-05-28 |
Family
ID=83578598
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210784961.0A Active CN115201914B (en) | 2022-06-30 | 2022-06-30 | Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115201914B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20160095690A (en) * | 2015-02-03 | 2016-08-12 | 한국해양대학교 산학협력단 | Seismic waveform inversion method and apparatus for underground velocity information derived from the single-channel seismic data |
| CN107356968A (en) * | 2017-08-25 | 2017-11-17 | 四川易诚智讯科技有限公司 | Three-dimensional level set tomography curved surface extraction method based on crease |
| CN109459789A (en) * | 2018-12-18 | 2019-03-12 | 吉林大学 | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation |
| CN113093272A (en) * | 2021-03-29 | 2021-07-09 | 吉林大学 | Time domain full waveform inversion method based on convolutional coding |
| US11079509B1 (en) * | 2021-01-19 | 2021-08-03 | Emerson Paradigm Holding Llc. | Building accurate training images for automatic seismic interpretation |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP3168653B1 (en) * | 2015-11-05 | 2019-07-17 | CGG Services SAS | Device and method for full waveform inversion |
-
2022
- 2022-06-30 CN CN202210784961.0A patent/CN115201914B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20160095690A (en) * | 2015-02-03 | 2016-08-12 | 한국해양대학교 산학협력단 | Seismic waveform inversion method and apparatus for underground velocity information derived from the single-channel seismic data |
| CN107356968A (en) * | 2017-08-25 | 2017-11-17 | 四川易诚智讯科技有限公司 | Three-dimensional level set tomography curved surface extraction method based on crease |
| CN109459789A (en) * | 2018-12-18 | 2019-03-12 | 吉林大学 | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation |
| US11079509B1 (en) * | 2021-01-19 | 2021-08-03 | Emerson Paradigm Holding Llc. | Building accurate training images for automatic seismic interpretation |
| CN113093272A (en) * | 2021-03-29 | 2021-07-09 | 吉林大学 | Time domain full waveform inversion method based on convolutional coding |
Non-Patent Citations (3)
| Title |
|---|
| Full waveform inversion based on dynamic data matching of convolutional wavefields;董士琦等;《Frontiers in Earth Science》;20230419;全文 * |
| 基于波场局部尺度匹配的时间域全波形反演研究;董士琦;《中国优秀博士论文全文数据库（基础科学辑）》;20220115;全文 * |
| 梁煌 ; 韩立国 ; 许卓 ; 胡勇 ; 邹佳儒 ; .互相关与最小二乘加权目标函数全波形反演.世界地质.2017,(02),全文. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115201914A (en) | 2022-10-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109407151B (en) | Time-domain full waveform inversion method based on wave field local correlation time shift | |
| CN111239802B (en) | Deep learning speed modeling method based on seismic reflection waveform and velocity spectrum | |
| US8923093B2 (en) | Determining the quality of a seismic inversion | |
| CN105467444B (en) | A kind of elastic wave full waveform inversion method and device | |
| CN113093272A (en) | Time domain full waveform inversion method based on convolutional coding | |
| CN106405651B (en) | Full waveform inversion initial velocity model construction method based on logging matching | |
| CN109459789B (en) | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation | |
| CN111239806B (en) | Time domain full waveform inversion method based on amplitude increment coding | |
| CN102099705A (en) | System and method for processing seismic data for interpretation | |
| CN113687433A (en) | Bi-LSTM-based magnetotelluric signal denoising method and system | |
| CN110531410A (en) | A kind of least square reverse-time migration gradient Preconditioning method based on through wave field | |
| CN111025388B (en) | Multi-wave combined prestack waveform inversion method | |
| CN109738950A (en) | The noisy-type data primary wave inversion method of domain inverting is focused based on sparse 3 D | |
| CN107942374A (en) | Diffracted wave field extracting method and device | |
| Romero et al. | Seeing through the CO2 plume: Joint inversion-segmentation of the Sleipner 4D seismic data set | |
| CN115201914B (en) | Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction | |
| CN111239833B (en) | K-value robust YPD (pre-stack simultaneous inversion) method based on Poisson ratio decomposition | |
| Sheng et al. | Arrival-time picking of microseismic events based on MSNet | |
| CN111273346B (en) | Method, device, computer equipment and readable storage medium for removing deposition background | |
| CN108680957B (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting | |
| Zhou et al. | Optimizing orthogonal-octahedron finite-difference scheme for 3D acoustic wave modeling by combination of Taylor-series expansion and Remez exchange method | |
| CN112630839B (en) | Logging curve standardization method and system | |
| CN111208568B (en) | Time domain multi-scale full waveform inversion method and system | |
| CN110888158B (en) | Full waveform inversion method based on RTM constraint | |
| CN113138412A (en) | Deep shale porosity earthquake prediction method and device |
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 |