CN115201914A - 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
- CN115201914A CN115201914A CN202210784961.0A CN202210784961A CN115201914A CN 115201914 A CN115201914 A CN 115201914A CN 202210784961 A CN202210784961 A CN 202210784961A CN 115201914 A CN115201914 A CN 115201914A
- Authority
- CN
- China
- Prior art keywords
- data
- matching
- representing
- waveform inversion
- time
- 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 43
- 238000000605 extraction Methods 0.000 title claims abstract description 21
- 238000005070 sampling Methods 0.000 claims abstract description 58
- 238000004088 simulation Methods 0.000 claims abstract description 21
- 238000004364 calculation method Methods 0.000 claims description 17
- 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 9
- 238000010606 normalization Methods 0.000 claims description 4
- 230000001427 coherent effect Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 230000006870 function Effects 0.000 description 20
- 239000011159 matrix material Substances 0.000 description 11
- 230000008569 process Effects 0.000 description 8
- 230000036039 immunity Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 230000002238 attenuated effect Effects 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 230000000670 limiting effect Effects 0.000 description 2
- 230000001902 propagating 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
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects 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
- 230000002452 interceptive effect Effects 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 230000008707 rearrangement Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
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 characteristic extraction, which comprises the steps of performing convolution operation on single-channel observation data and simulation data through one-dimensional Gaussian convolution checks with different lengths to extract the comprehensive attribute characteristics of a wave field of each sampling point of the data 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, then obtaining the unique optimal matching point of each analog data sampling point in the observation data, breaking the limitation of point-to-point matching in time of the conventional full waveform inversion, and realizing dynamic matching; introducing attenuation factors, and carrying out amplitude exponential attenuation on the dynamically matched analog data according to the magnitude of the absolute value of the optimal matching time difference; a global cross-correlation objective function is employed. The method breaks through the limitation 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, and improves the convergence of the full waveform inversion.
Description
Technical Field
The invention belongs to the technical field of seismic exploration, particularly relates 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 detection depth and high resolution, and is acknowledged 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 the reserve of shallow resources are basically ascertained, the gravity center of resource exploration gradually turns to a deep complex geological structure. Against this background, countries put higher demands on the depth and accuracy of seismic exploration. The main function of seismic exploration is to realize the exploration of underground geological structure and lithology information after a series of data processing according to seismic signals acquired by detectors arranged on the earth surface. Therefore, it is important to obtain the distribution of the physical parameters of the subsurface medium through seismic data.
Among the physical parameters, the information of the propagation velocity of seismic waves in different media is a key factor determining the quality of offset imaging. At present, the underground medium velocity modeling method widely used at home and abroad mainly comprises three methods of velocity analysis, tomography and full waveform inversion. The velocity analysis and tomography only use travel time information of seismic wave propagation and neglect many other important information, so that the velocity modeling result cannot meet the requirement of high-precision seismic exploration. Full Waveform Inversion (FWI) techniques utilize all the information of the prestack seismic wavefield (including the kinematic and kinetic information such as amplitude, phase, travel time, etc.), and achieve high-precision velocity modeling of the subsurface medium. Besides modeling the velocity parameters, full waveform inversion can also model parameters such as density, quality factors and anisotropy coefficients with high precision.
Laily (1983) and Tarantola (1984) transform the strong nonlinear problem of subsurface medium velocity inversion into an optimization problem in the least squares sense. By minimizing the observation data and the simulation data under the constraint of the L2 norm, the initial velocity model is continuously updated in an iterative manner to enable the simulation data to gradually approach the observation data, so that the optimal solution for the velocity parameter inversion of the underground medium is found. The direct calculation of the Frechet derivative results in a large amount of calculation, which makes it difficult to perform full waveform inversion under the condition of insufficient calculation capability. With the adoption of the state method, the gradient of model updating can be obtained under the condition that the Frechet derivative is not directly calculated, so that the calculation efficiency of full waveform inversion is greatly improved, and the enthusiasm of domestic and foreign researchers on full waveform inversion research is greatly promoted. At present, full waveform inversion becomes a hot topic in the field of exploration geophysical by virtue of the advantage of high-precision modeling of various underground physical parameters, and numerous scientific researchers at home and abroad deeply and meticulously research various aspects of full waveform inversion, thereby greatly promoting the development of the method.
Meanwhile, full waveform inversion is a strong nonlinear problem and is approximate to 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 geological formation regions, 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 seismic source is excited and the detector responds to low frequency information in the seismic signal, the low frequency component is difficult to exist or weak in energy in the seismic data obtained by actual operation. Therefore, the cycle skip problem arises, which is also the main problem facing full waveform inversion. In addition, the noise interference resistance and the influence of wavelet error reduction on inversion are problems to be researched and solved at present. The traditional full waveform inversion based on the L2 norm objective function is that observed data and simulated data are subjected to point-to-point subtraction in the time direction, the amplitude accuracy dependency on the simulated data is high, and the noise immunity is poor. The target function based on the cross correlation reduces the dependence of the FWI on the amplitude information precision, highlights the matching of the phase information with better linearity, solves the problem of interference of amplitude information errors or uneven energy of signals with different offset distances to an inversion result caused by the seismic acquisition process, and has certain noise immunity.
In order to alleviate the cycle skip phenomenon, many methods are currently developed, such as Bunks (1995) which proposes to decompose seismic data into different frequency scales by dividing frequency bands by a low-pass filter, and an inversion strategy from a low frequency band to a high frequency band in a time domain can be adopted to enable a FWI objective function to be better converged, reduce the occurrence of local extrema, and alleviate the cycle skip problem. Shin and Cha (2008) change seismic data into a Laplace domain, and data containing rich low-frequency information can be obtained through e exponential attenuation. Laplace domain waveform inversion can construct long-wavelength components of the subsurface medium. Zhang et al (2017) propose to construct 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 a dynamic warping method using the signal identification field, and find out the best matching point between two channels according to the similarity between the same channel number of data, which overcomes the defect of point-to-point subtraction in sampling time based on the L2 norm FWI in the prior art, and successfully alleviates the problem of cycle skip. Frequency extension of seismic data can be achieved by utilizing compressed sensing and sparse inversion in digital signal processing. However, these methods rely on other technical means for assistance or have poor computational stability, and have some disadvantages 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 through the limitation of point-to-point matching on sampling time of conventional full waveform inversion, utilizes a one-dimensional convolution kernel to extract the features of single-channel seismic data to realize dynamic matching of simulation records and observation records, solves the problem of cycle jump, improves the convergence of full waveform inversion and the accuracy of inversion results, and does not obviously increase extra calculated amount.
The purpose of the invention is realized by 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. preprocessing actual earthquake observation records, such as de-coding, elevation static correction, first arrival static correction, surface consistency static correction, random noise and coherent noise removal, wavelet estimation, missing seismic channel compensation and the like;
b. firstly, establishing a linear increasing initial model in a pre-estimated speed range, and setting time domain full waveform inversion related parameters including seismic wavelet main frequency f, model size nz multiplied by nx, grid distances dx and dz, total sampling time T, time sampling interval dt, maximum iteration times itermax, iteration step length q of an optimization algorithm, target function required precision tol, and maximum value vmax and minimum value vmin of model speed estimation according to requirements.
c. Performing finite difference forward modeling calculation on the initial model by using the estimated wavelet according to an actual observation system to obtain a simulated seismic record;
d. designing a one-dimensional Gaussian convolution kernel:
k g representing a gaussian convolution kernel, n representing an element in the convolution kernel, l representing the length of the convolution kernel, exp (.) representing an exponential function based on a natural constant e, n m Representing the center point of the convolution kernel;
e. extracting wave field characteristics of single-channel seismic data of observation data and simulation data in time ranges with different convolution kernel lengths by using Gaussian convolution kernels with different lengths:
* Which represents a convolution operation, is a function of,
and
respectively representing the simulation data and observation data of the jth track of the ith cannon, ns and nr respectively representing the total cannon number and the total track number of each cannon, T and T respectively representing time variable and maximum sampling time,
and
respectively representing the features of the single-channel analog data and observed data extracted by the w-th length of Gaussian convolution kernel, nw representing the total number of different-length Gaussian convolution kernels used, the length l of the convolution kernel being from a first value l 1 To the last value l nw Monotonically increasing, | denotes taking the absolute value of the data, norm (.) denotes the operator normalizing the data, M [.]Meaning that the convolution result is centered so that the cropped data feature duration is the same as the original seismic record, i.e. n is extracted from the convolution result m To (T + 1-n) m ) Characteristic data over a period of time;
f. constructing a characteristic vector of each sampling point in single-channel seismic data:
g. calculating Euclidean distances between each sampling point of the single-channel analog data and the characteristics of all sampling points corresponding to the observed data, and finding out the observation 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 representing sample points, t, in analogue data o Representing a sampling point in the observed data which is best matched with the simulated data;
h. calculating the optimal matching time difference:
Δt＝t s -t o ,
Δ t represents the best match time difference;
i. according to the magnitude of the absolute value of the optimal matching time difference, carrying out amplitude attenuation on the analog data sampling points to realize the restraint on the reliability of dynamic matching, wherein the larger the absolute value of the optimal matching time difference is, the higher the attenuation degree of the analog data sampling points is:
j. establishing a global cross-correlation full-waveform inversion target function:
j represents an objective function, v represents velocity information of the underground medium, and the gradient expression obtained by solving the derivative of the velocity at two ends of the objective function is as follows:
P f representing the time domain forward wave field, P b Representing a back-propagation wave field;
k. and carrying out iterative updating on the velocity model by using an L-BFGS optimization algorithm, inverting the large-scale structure of the model by using a time domain dynamic wave field matching full waveform inversion method based on convolution characteristic extraction, carrying out full-band conventional full waveform inversion by using 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, the principle is simple and easy to realize by programming, single-channel observation data and simulation data are subjected to convolution operation through one-dimensional Gaussian convolution kernels with different lengths to extract the comprehensive attribute features (amplitude, phase, travel time and the like) of the wave field of each sampling point of the data in different time ranges; the method breaks through the limitation of point-to-point subtraction on time in 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 seismic exploration inversion results; 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 characteristics of the single-channel seismic data to realize the dynamic matching of the simulation data and the observation data, and has the following advantages:
1. the comprehensive characteristics of wave fields of single-channel observation data and simulation data in different time ranges are extracted by utilizing a plurality of one-dimensional convolution kernels, the optimal matching point pair can be found according to the minimum Euclidean distance of all characteristics of each sampling point corresponding to the simulation data and the observation data, the limitation of point-to-point matching of the conventional full waveform inversion on the same time is broken through, the self-adaptive dynamic matching according to the characteristic change of the wave field data is realized, the cycle jump phenomenon caused by waveform matching errors is effectively eliminated, the convergence of the full waveform inversion is obviously enhanced, and the accuracy of the inversion result is improved;
2. the global cross-correlation objective function is adopted to reduce the dependence of inversion on amplitude information. Because the amplitude of part of sampling points of the analog record is attenuated, the phase information of each sampling point is closer to the observation record, but the corresponding amplitude information is not corrected, the influence of amplitude error on inversion is reduced by adopting a global cross-correlation target function, meanwhile, the weight of the far offset record is increased by cross-correlation calculation, and the far offset record always 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 features of data can still be accurately extracted under the condition that the observed data contains strong noise, in addition, the target function based on cross-correlation also has strong noise immunity, so that the time domain dynamic wave field matching full waveform inversion method based on 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 needed to be used in the embodiments will be briefly described below, it should be understood that the following drawings only illustrate some embodiments of the present invention and therefore should not be considered as limiting the scope, and for those skilled in the art, other related drawings can be obtained according to the drawings without inventive efforts.
FIG. 1 is a flow chart of a time domain dynamic wavefield matching full waveform inversion method based on convolution feature extraction;
FIG. 2a is a graph of a 20Hz dominant frequency Rake wavelet with information below 11Hz removed;
FIG. 2b the wavelet spectrum of FIG. 2 a;
FIG. 3a is a true model diagram;
FIG. 3b is a diagram of a linear incremental initial model;
FIG. 4a is a graph of inversion results based on a conventional time domain full waveform inversion;
FIG. 4b is a time domain dynamic wavefield matching full waveform inversion method result diagram based on convolution feature extraction.
Detailed Description
The present invention will be described in further detail with reference to the accompanying drawings and examples. It is to be understood that the specific embodiments described herein are merely illustrative of the invention and are not to be construed as limiting the invention. It should be further noted that, for the convenience of description, only some structures related to the present invention are shown in the drawings, not all of them.
It should be noted that: like reference numbers and letters refer to like items in the following figures, and thus, once an item is defined in one figure, it need not be further defined or explained in subsequent figures. Meanwhile, in the description of the present invention, the terms "first", "second", and the like are used only for distinguishing 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 an MATLAB2016b platform, a program is written and completed under an MATLAB2016b software framework, an installation environment of an MATLAB Parallel working library is built according to corresponding Parallel Computing requirements, and an MATLAB Parallel Computing Toolbox (Parallel Computing Toolbox) is installed.
The core of the method 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 it is difficult to provide an accurate initial model for full waveform inversion, the cycle jump problem caused by direct matching of medium-high frequency signals is the most important reason influencing the accuracy of the full waveform inversion result. And carrying out convolution operation on the single-channel observation data and the analog data through a one-dimensional Gaussian convolution kernel, and defining the characteristic of each sampling point as a convolution value when a convolution kernel directly faces a certain sampling point during convolution calculation, namely the characteristic value of the sampling point, so that the length of the convolution kernel is selected to be 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 precision of dynamic matching. And calculating Euclidean distances of all characteristics between each sampling point of the single-channel simulation data and all sampling points of the corresponding single-channel observation data, finding out an observation data sampling point which is optimally matched with each sampling point of the simulation data according to the Euclidean distance minimum principle of the sum of the characteristics, and obtaining the optimal matching time difference, wherein the optimal matching time difference can be positive or negative, and can also be zero. In order to evaluate the reliability of dynamic matching and constrain unreliable dynamic matching to avoid artificially generating more cycle jumps, the invention provides amplitude exponential attenuation on the dynamically matched analog data according to the magnitude of the optimal matching time difference absolute value, namely, the larger the optimal matching time difference absolute value is, the lower the possibility that the matching is correct is, the more the amplitude of the analog data needs to be attenuated to a greater extent. The invention effectively solves the problem of cycle jump without obviously increasing the calculated amount.
The invention relates to a time domain dynamic wave field matching full waveform inversion method based on convolution characteristic extraction, which comprises the following steps of:
a. and (3) preprocessing actual seismic observation records such as de-coding, elevation static correction, first arrival static correction, surface consistency static correction, random noise and coherent noise removal, wavelet estimation, missing seismic channel compensation and the like.
b. Firstly, establishing a linear increasing initial model in a pre-estimated speed range, and setting time domain full waveform inversion related parameters including seismic wavelet main frequency f, model size nz multiplied by nx, grid distances dx and dz, total sampling time T, time sampling interval dt, maximum iteration times itermax, iteration step length q of an optimization algorithm, target function required precision tol, and maximum value vmax and minimum value vmin of model speed estimation according to requirements.
c. And performing finite difference forward calculation on the initial model according to an actual observation system by using the estimated wavelet to obtain a simulated seismic record.
d. Designing a one-dimensional Gaussian convolution kernel:
k g representing 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 with a natural constant e as the base. n is a radical of an alkyl radical m Representing the center point of the convolution kernel.
e. Extracting wave field characteristics of single-channel seismic data of observation data and simulation data in time ranges of different convolution kernel lengths by using Gaussian convolution kernels of different lengths:
* Representing a convolution operation.
And
respectively representing the simulation data and the observation data of the jth path of the ith gun. ns and nr represent the total number of shots and total number of tracks per shot, respectively. T and T represent the time variable and the maximum sampling time, respectively.
And
respectively representing the characteristics of the single-channel simulation data and the observation data extracted by the w-th length of the Gaussian convolution kernel. nw represents the total number of different length gaussian convolution kernels used. The length l of the convolution kernel is derived from the first value l 1 To the last value l nw Monotonically increasing. |. | represents taking the absolute value of the data. The method comprises the steps of firstly, taking an absolute value of single-channel seismic data to avoid that the feature difference of different convolution kernels is not obvious due to the fact that positive and negative amplitudes are subjected to subtraction interference during convolution feature extraction. norm (.) represents the operator that normalizes the data. The single-channel observation data and the analog number after the absolute value is takenAnd the abnormal amplitude value is prevented from interfering the convolution characteristic extraction according to the normalization. M.]Meaning that the convolution result is centered so that the cropped data feature duration is the same as the original seismic record, i.e. n is extracted from the convolution result m To (T + 1-n) m ) Characteristic data over a period of time. And finally, carrying out second normalization processing to ensure 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:
g. Calculating Euclidean distances between each sampling point of the single-channel analog data and the characteristics of all sampling points corresponding to the observed data, and finding out the observation 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 representing sample points, t, in analogue data o Representing the sample points in the observed data that best match the simulated data.
h. Calculating the optimal matching time difference:
Δt＝t s -t o ,
Δ 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 point to realize the constraint on the reliability of dynamic matching. The greater the absolute value of the best match time difference, the higher the attenuation of the analog data sampling point to it:
j. Establishing a global cross-correlation full-waveform inversion target function:
j denotes an objective function. v represents velocity information of the subsurface medium. The gradient expression obtained by taking the derivative of the velocity at both ends of the objective function is as follows:
after simplification, can be expressed as:
wherein λ is a satellite source, and the expression is:
the time domain full waveform inversion gradient formula based on the adjoint state method can be obtained through a chain rule:
wherein P is f For time domain forward propagating wave field, L -1 λ represents the backward propagating wave field, where L -1 Is a passback operator.
k. Performing iterative updating on the inversion based on an L-BFGS optimization algorithm, and searching a proper step length through a Wolfe convergence criterion; the iterative formula is expressed as:
m k+1 ＝m k -α k H k g k
wherein H k Is an inverse of the approximate Hessian matrix, m k Updating parameters for the model, g k Update gradient, α, for model k For the step size obtained by Wolfe linear search, k represents the number of iterations.
The L-BFGS optimization algorithm has few matrixes to be stored in the iterative calculation process, has low requirement on the memory of a computer, and simultaneously has the following iterative formula for updating the Hessian matrix:
H k+1 ＝V k T H k V k +ρ k s k s k T
s k ＝m k+1 -m k ，y k ＝g k+1 -g k
wherein, V k To calculate intermediate variables of the process, s k To calculate intermediate variables of the process, p k For intermediate variables, I denotes the identity matrix, H k+1 Is based on vector pairs s k ,y k H and H k And calculating to obtain the vector, and only storing n vector pairs to implicitly express an inverse matrix of the Hessian matrix. H k g k May pass through a gradient g k And vector pair s k ,y k The inner product of a series of vectors between and the sum of the vectors. After each update, the last vector pair is updated by the current new vector pair s k+1 ,y k+1 And (4) substituting. Thus, the curvature signal of the last n iterations is contained in the set of vector pairs. In practice, when n is more than or equal to 5, a more satisfactory result can be obtained. Initial approximation Hessian matrix of L-BFGS optimization algorithm
Different in each iteration. Inverse matrix H of approximate Hessian matrix k The following update formula needs to be satisfied:
the updating direction of the model is realized by the following method:
(1) Let q = g k ，
Q = q- α i y i Wherein i = k-1,k-2, … k-n; alpha is alpha i Intermediate variables in the calculation process;
(2) Order to
Then r = r-s i (α i - β), wherein i = k-n, k-n + 1, …, k-1;
an initial approximate Hessian matrix;
(3) Obtaining the updated quantity H through the steps k g k ＝r。
Obtaining the update quantity H of the model by the method k g k Then obtaining a suitable step size alpha through Wolfe linear search k An update iteration is performed. The L-BFGS optimization algorithm effectively overcomes the difficulty of explicitly solving the Hessian matrix and the inverse of the Hessian matrix, has the advantages of super-linear convergence speed, high calculation efficiency, small occupied memory, high precision and the like, and is more suitable for solving the large-scale nonlinear optimization problem.
And judging whether the inversion result meets the set termination condition, namely the difference between the inversion result and the real model is very small, if so, stopping iteration, if not, continuing to calculate until the termination condition is met, and finishing the 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 carrying out convolution operation extraction, so that the weight of the feature extracted by each length convolution kernel on the dynamic matching contribution is the same; calculating characteristic Euclidean distances between each sampling point of the single-channel analog data and all sampling points of the corresponding channel of the analog data, obtaining a unique optimal matching point in the observed data and each analog data sampling point according to the minimum characteristic Euclidean distance, breaking the limitation 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 dynamically matched analog data according to the magnitude of the optimal matching time difference absolute value, namely the larger the dynamic matching time difference absolute value is, the higher the attenuation degree of a corresponding sampling point of the dynamically matched analog data is; because the amplitude of the analog data is attenuated in the process, the originally carried underground medium information is damaged; 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 requirement, a Parallel Computing Toolbox and an MATLAB Distributed Computing Server (R2016 b) are installed under a Windows 10 professional version system, and an MATLAB Parallel platform is built.
And carrying out full waveform inversion test 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 full waveform inversion test parameters for missing low frequency information
| Mesh size | Grid distance | Measuring length of wire | Longitudinal depth | Frequency bandwidth |
| 69×192 | 12.5m | 2400m | 862.5m | Full frequency band |
The size of the model grid 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 range in the model is 1500m/s to 4000m/s, the geophones are arranged on the surface of the model, each grid point is a geophone, the interval between the geophones is 12.5m, and 12 seismic sources are excited in total. The seismic source selects a Rake wavelet (figure 2 a) lacking information below 11Hz with the dominant frequency of 20Hz, 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, and the maximum length l of the convolution kernels nw Take 401 and the attenuation factor gamma takes 5.
The inversion results are shown in fig. 4. Fig. 4a shows that the conventional full waveform inversion result is greatly affected by the missing of low frequency information, and a plurality of obvious cycle skipping occurs. For a linear incremental initial model, the initial model has a large phase difference with a real model, and in the waveform matching process, the phase difference between multiple observation records and simulation records is more than half a cycle due to the absence of low-frequency information, cycle skipping occurs, and the gradient update quantity calculation error of the corresponding position is caused. With the increase of iteration times in the inversion process, the error speed updating amount is continuously accumulated, and finally the difference between the inversion result and the real model is larger. FIG. 4b shows the result of a full waveform inversion based on the method of the present invention. It can be seen from the figure that the inversion results can substantially invert all constructions of the real model without velocity anomaly regions. After the dynamic data extracted based on the convolution characteristics are matched, the matching degree of the simulation record and the observation record is improved, the model velocity value is updated to the correct direction, and a relatively accurate inversion result is finally obtained without a velocity abnormal area.
It is to be noted that the foregoing description is only exemplary of the invention and that the principles of the technology may be employed. 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, although the present invention has been described in greater detail by the above embodiments, the present invention is not limited to the above embodiments, and may include other equivalent embodiments without departing from the spirit of the present invention, and the scope of the present invention is determined by the scope of the appended claims.
Claims (1)
1. A time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction is characterized by comprising the following steps:
a. preprocessing actual earthquake observation records such as de-coding, elevation static correction, first arrival static correction, surface consistency static correction, random noise and coherent noise removal, wavelet estimation, missing seismic channel compensation and the like;
b. firstly, establishing a linear increasing initial model in a pre-estimated speed range, and setting time domain full waveform inversion related parameters including seismic wavelet main frequency f, model size nz multiplied by nx, grid distances dx and dz, total sampling time T, time sampling interval dt, maximum iteration times itermax, iteration step length q of an optimization algorithm, target function required precision tol, and maximum value vmax and minimum value vmin of model speed estimation according to requirements.
c. Performing finite difference forward modeling calculation on the initial model by using the estimated wavelet according to an actual observation system to obtain a simulated seismic record;
d. designing a one-dimensional Gaussian convolution kernel:
k g representing a gaussian convolution kernel, n representing an element in the convolution kernel, l representing the length of the convolution kernel, exp (.) representing an exponential function based on a natural constant e, n m Representing the center point of the convolution kernel;
e. extracting wave field characteristics of single-channel seismic data of observation data and simulation data in time ranges of different convolution kernel lengths by using Gaussian convolution kernels of different lengths:
* Which represents a convolution operation, is a function of,
and
respectively representing the simulation data and observation data of the jth path of the ith shot, ns and nr respectively representing the total number of shots and the total number of paths of each shot, T and T respectively representing a time variable and maximum sampling time,
and
respectively representing the features of the single-channel analog data and observed data extracted by the w-th length of Gaussian convolution kernel, nw representing the total number of different-length Gaussian convolution kernels used, the length l of the convolution kernel being from a first value l 1 To the last value l nw Monotonically increasing, | denotes taking the absolute value of the data, norm (|) denotes pairsOperator of data normalization, M.]Meaning that the convolution result is centered so that the cropped data feature duration is the same as the original seismic record, i.e. n is extracted from the convolution result m To (T + 1-n) m ) Characteristic data over a period of time;
f. constructing a characteristic vector of each sampling point in single-channel seismic data:
g. calculating Euclidean distances between each sampling point of the single-channel analog data and the characteristics of all sampling points corresponding to the observed data, and finding out the observation 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 representing sample points, t, in analogue data o Representing a sampling point in the observed data which is best matched with the simulated data;
h. calculating the optimal matching time difference:
Δt＝t s -t o ,
Δ t represents the best match time difference;
i. according to the magnitude of the absolute value of the optimal matching time difference, carrying out amplitude attenuation on the analog data sampling points to realize the restraint on the reliability of dynamic matching, wherein the larger the absolute value of the optimal matching time difference is, the higher the attenuation degree of the analog data sampling points is:
j. establishing a global cross-correlation full-waveform inversion target function:
j represents an objective function, v represents velocity information of the underground medium, and the gradient expression obtained by solving the derivative of the velocity at two ends of the objective function is as follows:
P f representing the time domain forward wave field, P b Representing a back-propagation wave field;
k. and carrying out iterative updating on the velocity model by using an L-BFGS optimization algorithm, inverting the large-scale structure of the model by using a time domain dynamic wave field matching full waveform inversion method based on convolution characteristic extraction, carrying out full-band conventional full waveform inversion by using 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 true CN115201914A (en) | 2022-10-18 |
| CN115201914B 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 (6)
| 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 |
| US20170131418A1 (en) * | 2015-11-05 | 2017-05-11 | Cgg Services Sa | Device and method for full waveform inversion |
| 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 |
-
2022
- 2022-06-30 CN CN202210784961.0A patent/CN115201914B/en active Active
Patent Citations (6)
| 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 |
| US20170131418A1 (en) * | 2015-11-05 | 2017-05-11 | Cgg Services Sa | Device and method for full waveform inversion |
| 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 (2)
| Title |
|---|
| 董士琦: "基于波场局部尺度匹配的时间域全波形反演研究", 《中国优秀博士论文全文数据库（基础科学辑）》, 15 January 2022 (2022-01-15) * |
| 董士琦等: "Full waveform inversion based on dynamic data matching of convolutional wavefields", 《FRONTIERS IN EARTH SCIENCE》, 19 April 2023 (2023-04-19) * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115201914B (en) | 2024-05-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109407151B (en) | Time-domain full waveform inversion method based on wave field local correlation time shift | |
| CN113093272A (en) | Time domain full waveform inversion method based on convolutional coding | |
| US8923093B2 (en) | Determining the quality of a seismic inversion | |
| US8880384B2 (en) | Artifact reduction in iterative inversion of geophysical data | |
| US8509027B2 (en) | Continuous adaptive surface wave analysis for three-dimensional seismic data | |
| CN109459789B (en) | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation | |
| CN102099705A (en) | System and method for processing seismic data for interpretation | |
| CN103954992B (en) | Deconvolution method and device | |
| CN111239806B (en) | Time domain full waveform inversion method based on amplitude increment coding | |
| US8139440B2 (en) | Spectral conditioning for surface seismic data | |
| CN111045076A (en) | Multi-mode Rayleigh wave frequency dispersion curve parallel joint inversion method | |
| CN111665556B (en) | Stratum acoustic wave propagation velocity model construction method | |
| CA2412995C (en) | Seismic survey system | |
| Romero et al. | Seeing through the CO2 plume: Joint inversion-segmentation of the Sleipner 4D seismic data set | |
| Zhou et al. | Optimizing orthogonal-octahedron finite-difference scheme for 3D acoustic wave modeling by combination of Taylor-series expansion and Remez exchange method | |
| CN115201914B (en) | Time domain dynamic wave field matching full waveform inversion method based on convolution feature extraction | |
| CN108680957B (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting | |
| MacBeth et al. | Upper crustal shear velocity models from higher mode Rayleigh wave dispersion in Scotland | |
| CN111208568B (en) | Time domain multi-scale full waveform inversion method and system | |
| CN110888158B (en) | Full waveform inversion method based on RTM constraint | |
| CN112630839B (en) | Logging curve standardization method and system | |
| CN112363222A (en) | Post-stack adaptive broadband constraint wave impedance inversion method and device | |
| CN113640880B (en) | Sandstone type uranium ore seismic exploration method and device | |
| CN112558147B (en) | Polarization analysis method and system for microseism data in well | |
| CN115993632A (en) | Method and system for optimizing key acquisition parameters |
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 |