CN111580163B - Full waveform inversion method and system based on non-monotonic search technology - Google Patents
Full waveform inversion method and system based on non-monotonic search technology Download PDFInfo
- Publication number
- CN111580163B CN111580163B CN202010466886.4A CN202010466886A CN111580163B CN 111580163 B CN111580163 B CN 111580163B CN 202010466886 A CN202010466886 A CN 202010466886A CN 111580163 B CN111580163 B CN 111580163B
- Authority
- CN
- China
- Prior art keywords
- iteration
- equal
- objective function
- monotonic
- seismic
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 92
- 230000003121 nonmonotonic effect Effects 0.000 title claims abstract description 51
- 238000005516 engineering process Methods 0.000 title claims abstract description 32
- 238000004364 calculation method Methods 0.000 claims abstract description 15
- 238000003384 imaging method Methods 0.000 claims abstract description 12
- 230000006870 function Effects 0.000 claims description 117
- 238000004088 simulation Methods 0.000 claims description 12
- 238000010276 construction Methods 0.000 claims description 6
- 238000001514 detection method Methods 0.000 claims description 6
- 238000010586 diagram Methods 0.000 description 9
- 238000005457 optimization Methods 0.000 description 5
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000011159 matrix material Substances 0.000 description 1
- 230000000750 progressive effect Effects 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/16—Receiving elements for seismic signals; Arrangements or adaptations of receiving elements
- G01V1/18—Receiving elements, e.g. seismometer, geophone or torque detectors, for localised single point measurements
- G01V1/181—Geophones
-
- 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/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
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)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention aims to provide a full waveform inversion method and a full waveform inversion system based on a non-monotonic search technology, wherein the method comprises the following steps: firstly, dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; secondly, determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length; then judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; and finally, performing seismic imaging based on the optimal velocity model. According to the method, the non-monotonic search technology is combined with the initial step length to determine the step length, and the method can obtain a high-precision inversion result on the premise of ensuring the inversion calculation efficiency; and when the initial velocity model is poor, the global convergence of inversion is ensured to a certain extent, and a better inversion result is obtained.
Description
Technical Field
The invention relates to the technical field of seismic exploration, in particular to a full waveform inversion method and system based on a non-monotonic search technology.
Background
Full Waveform Inversion (FWI) comprehensively utilizes Full wave information in a seismic wave field, and quantitatively corrects an initial model by fitting an actual observed wave field and a predicted wave field to reconstruct a parameter model of the underground medium, such as velocity, density and the like, so that high-precision elastic parameter information of longitudinal and transverse wave velocity, density and the like can be provided for seismic exploration and exploration of the internal structure of the earth.
Full waveform inversion is generally a nonlinear data fitting process that is repeated in an iterative fashion until the objective function formed by the observed data and the synthetic seismic data is sufficiently small. Firstly, an initial estimation value of underground stratum parameters is required to be given, an objective function and the gradient thereof are obtained, a search direction is constructed by combining an inversion optimization algorithm, the step length is calculated, the objective function is gradually reduced by continuously updating model parameters, and a satisfactory result is finally obtained. For the nonlinear seismic waveform inversion problem, the iteration step size must be obtained by adopting a linear search method along the search direction. The linear search method mainly comprises two main methods of non-precise linear search and precise linear search. The precise linear search method is a calculation method for solving from an objective function to obtain an optimal precise step length, and different objective functions correspond to different precise step lengths. The parabolic method and the direct solution method are widely used accurate linear search methods, but both are obtained by solving based on a two-norm type objective function. And a precise step length calculation method for adapting to the mixed norm is also provided for the Huber norm and the Hybrid norm. Although the accurate linear search method has certain advantages in inversion accuracy and convergence speed, different accurate step lengths need to be calculated for different target functions, the universality is poor, target functions with good development prospects such as cross-correlation target functions and envelope type target functions at present, and the accurate step length calculation method has not been researched.
Compared with the accurate linear search method, the non-accurate linear search method has wide applicability to the target function. The non-precise linear search method mainly utilizes the judgment condition and the initial step length to ensure that the target function has enough descending amount in the search direction in each iteration process, can generally make the target function descend quickly, does not require the target function to reach the precise minimum in each step, and only needs proper descending amount. At present, the traditional non-precise linear search method in seismic full waveform inversion generally uses a monotone search technology to require that an objective function is gradually decreased, such as a classic Wolfe judgment condition. Although the method can also make inversion converge, the objective function is forced to be gradually reduced, the global convergence property of the inversion optimization method is damaged, and when the inversion falls into the vicinity of a local minimum value, small step length and zigzag phenomena occur, the calculation efficiency and accuracy of the inversion are reduced, and the gradual reduction of the objective function is not beneficial to protecting the global convergence of the inversion.
Disclosure of Invention
Based on this, the invention aims to provide a full waveform inversion method and a full waveform inversion system based on a non-monotonic search technology, which combine the non-monotonic search technology with an initial step size to determine the step size so as to improve the calculation efficiency and accuracy of full waveform inversion.
In order to achieve the above object, the present invention provides a full waveform inversion method based on a non-monotonic search technique, the method comprising:
step S1: collecting seismic observation data by using a geophone;
step S2: dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1;
step S3: determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length;
step S4: judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; if j is less than or equal to N, making j equal to j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3;
step S5: and performing seismic imaging based on the optimal velocity model.
Optionally, the determining, by using a non-monotonic search technique in combination with the initial step length, a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model specifically includes:
step S31: solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method, and calculating the simulation wave field data of the kth iteration of the jth frequency group;
step S32: constructing an objective function value according to the seismic observation data and the simulated wave field data;
step S33: calculating a gradient from the simulated wavefield data using an adjoint method;
step S34: constructing a search direction according to the gradient by using a finite memory quasi-Newton method;
step S35: determining the step size of seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology in combination with the initial step size;
step S36: determining a speed model of the next iteration corresponding to the jth frequency group according to the search direction and the step length;
step S37: judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, let k be k +1, and return to "step S31"; if the iteration number k is greater than or equal to the iteration number threshold, executing step S4;
and/or calculating a relative difference value of the continuous objective function, and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient; if the relative difference of the continuous objective functions is greater than or equal to a set coefficient, making k equal to k +1, and returning to step S31; if the continuous objective function relative difference is smaller than the set coefficient, "step S4" is performed.
Optionally, the step length of the seismic full waveform inversion is determined according to the objective function value, the search direction, and the gradient by using a non-monotonic search technique in combination with an initial step length, and a specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkIs a non-linear degree, takes a value of 0 or 1,Qkfor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
representing the initial step size of the k-th iteration.
Optionally, the gradient is calculated according to the simulated wave field data by using an adjoint state method, and a specific formula is as follows:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
Optionally, the objective function value is constructed according to the seismic observation data and the simulated wave field data, and the specific formula is as follows:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
The invention also provides a full waveform inversion system based on a non-monotonic search technique, the system comprising:
the acquisition module is used for acquiring earthquake observation data by using a wave detector;
the grouping module is used for dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1;
the velocity model determining module is used for determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length;
the judging module is used for judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; if j is less than or equal to N, making j equal to j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3;
and the seismic imaging module is used for performing seismic imaging based on the optimal velocity model.
Optionally, the speed model determining module specifically includes:
the simulation wave field data determining unit is used for solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method and calculating simulation wave field data of the kth iteration of the jth frequency group;
the objective function value construction unit is used for constructing an objective function value according to the seismic observation data and the simulated wave field data;
a gradient calculation unit for calculating a gradient from the simulated wave field data using an adjoint method;
a search direction construction unit for constructing a search direction according to the gradient by using a finite memory quasi-Newton method;
the step length determining unit is used for determining the step length of the seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology and combining an initial step length;
the speed model determining unit is used for determining a speed model of the next iteration corresponding to the jth frequency group according to the searching direction and the step length;
the judging unit is used for judging whether the iteration times are larger than or equal to the iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, making k equal to k +1, and returning to the simulated wave field data determining unit; if the iteration number k is greater than or equal to the iteration number threshold, executing a judgment module;
and/or calculating a relative difference value of the continuous objective function, and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a set coefficient, enabling k to be k +1, and returning to a simulation wave field data determining unit; and if the relative difference value of the continuous objective functions is smaller than a set coefficient, executing a judgment module.
Optionally, the step length of the seismic full waveform inversion is determined according to the objective function value, the search direction, and the gradient by using a non-monotonic search technique in combination with an initial step length, and a specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkFor non-linearity, values of 0 or 1, QkFor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
representing the initial step size of the k-th iteration.
Optionally, the gradient is calculated according to the simulated wave field data by using an adjoint state method, and a specific formula is as follows:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
Optionally, the objective function value is constructed according to the seismic observation data and the simulated wave field data, and the specific formula is as follows:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention aims to provide a full waveform inversion method and a full waveform inversion system based on a non-monotonic search technology, wherein the method comprises the following steps: firstly, dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; secondly, determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length; then judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; and finally, performing seismic imaging based on the optimal velocity model. According to the method, the non-monotonic search technology is combined with the initial step length to determine the step length, and the method can obtain a high-precision inversion result on the premise of ensuring the inversion calculation efficiency; and when the initial velocity model is poor, the global convergence of inversion is ensured to a certain extent, and a better inversion result is obtained.
Drawings
In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings needed to be used in the embodiments will be briefly described below, and it is obvious that the drawings in the following description are only some embodiments of the present invention, and it is obvious for those skilled in the art to obtain other drawings without inventive exercise.
FIG. 1 is a flow chart of a full waveform inversion method based on a non-monotonic search technique according to an embodiment of the present invention;
FIG. 2 is a block diagram of a full waveform inversion system based on a non-monotonic search technique according to an embodiment of the present invention;
FIG. 3 is a diagram of a Marmousi model according to an embodiment of the present invention;
FIG. 4 is a schematic illustration of a smooth model according to an embodiment of the present invention;
FIG. 5 is a diagram illustrating a first inversion result according to an embodiment of the present invention;
FIG. 6 is a schematic diagram of an Overthrast model according to an embodiment of the present invention;
FIG. 7 is a schematic view of a homogeneous model according to an embodiment of the present invention;
FIG. 8 is a diagram illustrating a second inversion result according to an embodiment of the present invention.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
The invention aims to provide a full waveform inversion method and a full waveform inversion system based on a non-monotonic search technology, which combine the non-monotonic search technology with an initial step size to determine the step size so as to improve the calculation efficiency and accuracy of full waveform inversion.
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in further detail below.
Fig. 1 is a flowchart of a full waveform inversion method based on a non-monotonic search technique according to an embodiment of the present invention, and as shown in fig. 1, the present invention discloses a full waveform inversion method based on a non-monotonic search technique, where the method includes:
step S1: and collecting seismic observation data by using a detector.
Step S2: dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1.
Step S3: and determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology and combining the initial step length.
Step S4: judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; and if j is less than or equal to N, taking j as j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3.
Step S5: and performing seismic imaging based on the optimal velocity model.
Before step S2, the method further includes: and denoising and amplitude energy compensation processing are carried out on the seismic observation data.
The individual steps are discussed in detail below:
step S3: determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology and combining an initial step length, wherein the method specifically comprises the following steps:
step S31: solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method, and calculating the simulation wave field data of the kth iteration of the jth frequency group; the simulated wavefield data is a forward wavefield.
Step S32: and constructing an objective function value according to the seismic observation data and the simulated wave field data, wherein a specific formula is as follows:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
Step S33: calculating the gradient according to the simulated wave field data by using an adjoint state method, wherein the concrete formula is as follows:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
Step S34: and constructing a search direction according to the gradient by using a finite memory quasi-Newton method, wherein a specific formula is as follows:
dk＝-Hk·gk (3)；
wherein d iskDenotes the search direction of the kth iteration, gkGradient, H, representing the k-th iterationkThe inverse of the hessian matrix for the kth iteration is represented.
Step S35: determining the step length of seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology and combining with an initial step length, wherein the specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkFor non-linearity, values of 0 or 1, QkFor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
representing the initial step size of the k-th iteration.
When etakThe judgment condition of equation (4) degenerates to the conventional monotone judgment condition, ηkWhen 1 is CkRepresents near QkAverage of the objective function in the case of sub-iterations. Of course, instead of selecting a certain number of target function averages under continuous iteration to perform non-monotonic condition determination, the maximum value of these target functions may be tried as Ck. The non-monotonic search technology does not require the monotonous decrease of the target function, allows the local increase of the target function in the inversion process, and not only fully utilizes the products in the iteration processThe generated related target function information also protects the global convergence property of inversion to a certain extent.
Step S36: determining a speed model of the next iteration corresponding to the jth frequency group according to the search direction and the step length, wherein the specific formula is as follows:
mk+1＝mk+αkdk (5)；
wherein m isk+1Velocity model for the k +1 th iteration, mkVelocity model for the kth iteration, αkStep size of the kth iteration, dkThe search direction for the kth iteration.
Step S37: judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, let k be k +1, and return to "step S31"; if the iteration number k is greater than or equal to the iteration number threshold, "step S4" is performed.
Or, step S37: calculating a relative difference value of a continuous target function, and judging whether the relative difference value of the continuous target function is smaller than a first set coefficient; if the relative difference of the continuous objective functions is greater than or equal to a set coefficient, making k equal to k +1, and returning to step S31; if the continuous objective function relative difference is smaller than the set coefficient, "step S4" is performed.
Or, step S37: judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, let k be k +1, and return to "step S31"; if the iteration number k is greater than or equal to the iteration number threshold, "step S38" is performed.
Step S38: calculating a relative difference value of a continuous target function, and judging whether the relative difference value of the continuous target function is smaller than a first set coefficient; if the relative difference of the continuous objective functions is greater than or equal to a set coefficient, making k equal to k +1, and returning to step S31; if the continuous objective function relative difference is smaller than the set coefficient, "step S4" is performed.
Step S4: judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration as the optimal speed model; if j is less than or equal to N, let j be j +1, take the speed model of the next iteration as the initial speed model, and return to step S3.
Step S5: and performing seismic imaging based on the optimal velocity model.
As shown in fig. 2, the present invention further provides a full waveform inversion system based on a non-monotonic search technique, the system comprising:
and the acquisition module 1 is used for acquiring seismic observation data by using a geophone.
The grouping module 2 is used for dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1.
And the velocity model determining module 3 is used for determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology and combining the initial step length.
The judging module 4 is used for judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; and if j is less than or equal to N, taking j as j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3.
And the seismic imaging module 5 is used for performing seismic imaging based on the optimal velocity model.
As an embodiment, the speed model determining module 3 of the present invention specifically includes:
and the simulated wave field data determining unit is used for solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method and calculating the simulated wave field data of the kth iteration of the jth frequency group.
And the objective function value construction unit is used for constructing an objective function value according to the seismic observation data and the simulated wave field data.
And the gradient calculation unit is used for calculating the gradient according to the simulated wave field data by using an adjoint state method.
And the search direction construction unit is used for constructing a search direction according to the gradient by using a finite memory quasi-Newton method.
And the step length determining unit is used for determining the step length of the seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology and combining with the initial step length.
And the speed model determining unit is used for determining a speed model of the next iteration corresponding to the jth frequency group according to the searching direction and the step length.
The first judging unit is used for judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, making k equal to k +1, and returning to the simulated wave field data determining unit; if the iteration number k is greater than or equal to the iteration number threshold, a judgment module is executed.
Or, the second judging unit is used for calculating the relative difference value of the continuous objective function and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a set coefficient, enabling k to be k +1, and returning to a simulation wave field data determining unit; and if the relative difference value of the continuous objective functions is smaller than a set coefficient, executing a judgment module.
Or, a third judging unit, configured to judge whether the iteration count is greater than or equal to an iteration count threshold; if the iteration number k is smaller than the iteration number threshold, making k equal to k +1, and returning to the simulated wave field data determining unit; if the iteration number k is greater than or equal to the iteration number threshold, executing a 'fourth judgment unit';
the fourth judging unit is used for calculating the relative difference value of the continuous objective function and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient or not; if the relative difference value of the continuous objective functions is larger than or equal to a set coefficient, enabling k to be k +1, and returning to a simulation wave field data determining unit; and if the relative difference value of the continuous objective functions is smaller than a set coefficient, executing a judgment module.
As an embodiment, the non-monotonic search technique is combined with an initial step size, and a step size of seismic full waveform inversion is determined according to the objective function value, the search direction, and the gradient, and a specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkFor non-linearity, values of 0 or 1, QkFor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
representing the initial step size of the k-th iteration.
In one embodiment, the method of calculating the gradient according to the simulated wave field data by using the adjoint state method includes:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
As an embodiment, the method for constructing an objective function value according to the seismic observation data and the simulated wave field data includes:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
The invention solves the optimization problem by utilizing the non-monotonic search technology, does not require the monotonous decrease of the objective function, allows the local increase of the objective function in the inversion process, and fully utilizes the related objective function information generated in the iteration process. The optimization problem is solved by matching the non-monotonic search technology with an inversion optimization algorithm, the calculation times of the objective function and the gradient in the inversion iteration process can be reduced, and the calculation efficiency is improved. Meanwhile, the global convergence property of inversion can be protected to a certain extent.
Fig. 3 is a schematic diagram of a Marmousi model according to an embodiment of the present invention, and fig. 4 is a schematic diagram of a smooth model according to an embodiment of the present invention, in which the Marmousi model is used as a true parameter model, the smooth model is used as an initial velocity model to perform inversion, and a specific inversion result is shown in fig. 5, and it can be seen from fig. 5 that both the overall structure and the fine structure of the Marmousi model are well inverted. The test can obtain high-precision inversion results by using 95 iterations and 105 additional objective function calculations. The example proves that the method can obtain the inversion result with high precision under the condition of ensuring the inversion calculation efficiency.
Fig. 6 is a schematic diagram of an overhhrust model according to an embodiment of the present invention, and fig. 7 is a schematic diagram of a uniform model according to an embodiment of the present invention, where the invention uses the overhhrust model as a real parameter model, and uses the uniform model as an initial velocity model for inversion, and a specific inversion result is shown in fig. 8, and it can be seen from fig. 8 that a non-monotonic search technique can still recover a main structure for reconstructing a middle-shallow part of the overhhrust model. This example illustrates to some extent that non-monotonic search techniques can preserve the global convergence of full waveform inversion when the initial model is poor.
The embodiments in the present description are described in a progressive manner, each embodiment focuses on differences from other embodiments, and the same and similar parts among the embodiments are referred to each other.
The principles and embodiments of the present invention have been described herein using specific examples, which are provided only to help understand the method and the core concept of the present invention; meanwhile, for a person skilled in the art, according to the idea of the present invention, the specific embodiments and the application range may be changed. In view of the above, the present disclosure should not be construed as limiting the invention.
Claims (6)
1. A method for full waveform inversion based on a non-monotonic search technique, the method comprising:
step S1: collecting seismic observation data by using a geophone;
step S2: dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1;
step S3: determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length;
step S4: judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; if j is less than or equal to N, making j equal to j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3;
step S5: performing seismic imaging based on the optimal velocity model;
determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology and combining an initial step length, wherein the method specifically comprises the following steps:
step S31: solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method, and calculating the simulation wave field data of the kth iteration of the jth frequency group;
step S32: constructing an objective function value according to the seismic observation data and the simulated wave field data;
step S33: calculating a gradient from the simulated wavefield data using an adjoint method;
step S34: constructing a search direction according to the gradient by using a finite memory quasi-Newton method;
step S35: determining the step length of seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology and combining with an initial step length, wherein the specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkFor non-linearity, values of 0 or 1, QkFor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
represents the initial step size of the kth iteration;
step S36: determining a speed model of the next iteration corresponding to the jth frequency group according to the search direction and the step length;
step S37: judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, let k be k +1, and return to "step S31"; if the iteration number k is greater than or equal to the iteration number threshold, executing step S4;
or, step S37: calculating a relative difference value of a continuous target function, and judging whether the relative difference value of the continuous target function is smaller than a first set coefficient; if the relative difference of the continuous objective functions is greater than or equal to a set coefficient, making k equal to k +1, and returning to step S31; if the relative difference of the continuous objective functions is smaller than a set coefficient, executing step S4;
or, step S37: judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, let k be k +1, and return to "step S31"; if the iteration number k is greater than or equal to the iteration number threshold, executing step S38;
step S38: calculating a relative difference value of a continuous target function, and judging whether the relative difference value of the continuous target function is smaller than a first set coefficient; if the relative difference of the continuous objective functions is greater than or equal to a set coefficient, making k equal to k +1, and returning to step S31; if the continuous objective function relative difference is smaller than the set coefficient, "step S4" is performed.
2. The full waveform inversion method based on the non-monotonic search technique as claimed in claim 1, wherein the step of calculating the gradient from the simulated wavefield data using the adjoint method comprises the following steps:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
3. The full waveform inversion method based on the non-monotonic search technique according to claim 1, wherein an objective function value is constructed according to the seismic observation data and the simulated wavefield data, and the objective function value is defined by the following formula:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
4. A full waveform inversion system based on a non-monotonic search technique, the system comprising:
the acquisition module is used for acquiring earthquake observation data by using a wave detector;
the grouping module is used for dividing the multi-scale inversion frequency into N frequency groups according to a set frequency group range; wherein N is a positive integer greater than or equal to 1;
the velocity model determining module is used for determining a velocity model of the next iteration corresponding to the jth frequency group according to the seismic observation data and the initial velocity model by adopting a non-monotonic search technology in combination with the initial step length;
the judging module is used for judging whether j is larger than N; if j is larger than N, outputting the speed model of the next iteration corresponding to the jth frequency group as the optimal speed model; if j is less than or equal to N, making j equal to j +1, taking the speed model of the next iteration corresponding to the j-th frequency group as the initial speed model, and returning to the step S3;
the seismic imaging module is used for carrying out seismic imaging based on the optimal velocity model;
the speed model determining module specifically comprises:
the simulation wave field data determining unit is used for solving an acoustic wave equation according to the initial velocity model by adopting a finite difference method and calculating simulation wave field data of the kth iteration of the jth frequency group;
the objective function value construction unit is used for constructing an objective function value according to the seismic observation data and the simulated wave field data;
a gradient calculation unit for calculating a gradient from the simulated wave field data using an adjoint method;
a search direction construction unit for constructing a search direction according to the gradient by using a finite memory quasi-Newton method;
the step length determining unit is used for determining the step length of the seismic full waveform inversion according to the objective function value, the search direction and the gradient by adopting a non-monotonic search technology and combining an initial step length;
the speed model determining unit is used for determining a speed model of the next iteration corresponding to the jth frequency group according to the searching direction and the step length;
the first judging unit is used for judging whether the iteration times are larger than or equal to an iteration time threshold value or not; if the iteration number k is smaller than the iteration number threshold, making k equal to k +1, and returning to the simulated wave field data determining unit; if the iteration number k is greater than or equal to the iteration number threshold, executing a judgment module;
or, the second judging unit is used for calculating the relative difference value of the continuous objective function and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a set coefficient, enabling k to be k +1, and returning to a simulation wave field data determining unit; if the relative difference value of the continuous objective functions is smaller than a set coefficient, executing a judgment module;
or, a third judging unit, configured to judge whether the iteration count is greater than or equal to an iteration count threshold; if the iteration number k is smaller than the iteration number threshold, making k equal to k +1, and returning to the simulated wave field data determining unit; if the iteration number k is greater than or equal to the iteration number threshold, executing a 'fourth judgment unit';
the fourth judging unit is used for calculating the relative difference value of the continuous objective function and judging whether the relative difference value of the continuous objective function is smaller than a first set coefficient or not; if the relative difference value of the continuous objective functions is larger than or equal to a set coefficient, enabling k to be k +1, and returning to a simulation wave field data determining unit; if the relative difference value of the continuous objective functions is smaller than a set coefficient, executing a judgment module;
the step length of the seismic full waveform inversion is determined according to the objective function value, the search direction and the gradient by adopting the non-monotonic search technology and combining the initial step length, and the specific formula is as follows:
wherein m isk+1Velocity model for the k +1 th iteration, αkStep size of the kth iteration, dkSearch direction for the k-th iteration, c1Is a constant number, Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηkFor non-linearity, values of 0 or 1, QkFor the number of successive objective functions selected, CkIs QkAverage of the sub-objective function, f (m)k+1) Value of the objective function for k +1 iterations, gkThe gradient of the k-th iteration is represented,
representing the initial step size of the k-th iteration.
5. The non-monotonic search technique based full waveform inversion system of claim 4 wherein the computing the gradient from the simulated wavefield data using an adjoint method is formulated as:
wherein, gkGradient, m, representing the k-th iterationkVelocity model for the kth iteration, T denotes the length of observation time, λkRepresenting the back-propagation residual wavefield, p, of the kth iterationk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wavefield data.
6. The non-monotonic search technique based full waveform inversion system of claim 4, wherein the objective function value is constructed from the seismic observation data and the simulated wavefield data by the following formula:
wherein, f (m)k) Representing the value of the objective function of the k-th iteration, mkVelocity model for the kth iteration, xsRepresenting the s-th source position, xrDenotes the location of the r-th detection point, t denotes time, pk(xr,t|xs) Indicating that time t is defined by being at xsAt the seismic source xrThe calculated simulated wave field data is processed,
indicating that time t is defined by being at xsAt the seismic source xrObserved wavefield data of (a).
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010466886.4A CN111580163B (en) | 2020-05-28 | 2020-05-28 | Full waveform inversion method and system based on non-monotonic search technology |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010466886.4A CN111580163B (en) | 2020-05-28 | 2020-05-28 | Full waveform inversion method and system based on non-monotonic search technology |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111580163A CN111580163A (en) | 2020-08-25 |
| CN111580163B true CN111580163B (en) | 2021-08-10 |
Family
ID=72125676
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010466886.4A Expired - Fee Related CN111580163B (en) | 2020-05-28 | 2020-05-28 | Full waveform inversion method and system based on non-monotonic search technology |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111580163B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112084655A (en) * | 2020-09-08 | 2020-12-15 | 南京众诚土地规划设计咨询有限公司 | Ground penetrating radar parameter inversion method based on non-monotonic line search |
| CN113552620A (en) * | 2021-09-07 | 2021-10-26 | 中国地震局地球物理研究所 | Step length calculation method and system suitable for waveform inversion |
| CN115184986B (en) * | 2022-06-28 | 2023-09-29 | 吉林大学 | Global envelope cross-correlation full waveform inversion method independent of seismic source |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107450102A (en) * | 2017-07-28 | 2017-12-08 | 西安交通大学 | Multiple dimensioned full waveform inversion method based on the controllable envelope generating operator of resolution ratio |
| CN108549100A (en) * | 2018-01-11 | 2018-09-18 | 吉林大学 | The multiple dimensioned full waveform inversion method of time-domain of frequency is opened up based on non-linear high order |
| CN110058302A (en) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | A kind of full waveform inversion method based on pre-conditional conjugate gradient accelerating algorithm |
| CN110058307A (en) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | A kind of full waveform inversion method based on quick quasi-Newton method |
| WO2019186286A1 (en) * | 2018-03-27 | 2019-10-03 | King Abdullah University Of Science And Technology | Robust full waveform inversion of seismic data method and device |
| WO2020009752A1 (en) * | 2018-07-02 | 2020-01-09 | Exxonmobil Upstream Research Company | Full wavefield inversion with an image-gather-flatness constraint |
-
2020
- 2020-05-28 CN CN202010466886.4A patent/CN111580163B/en not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107450102A (en) * | 2017-07-28 | 2017-12-08 | 西安交通大学 | Multiple dimensioned full waveform inversion method based on the controllable envelope generating operator of resolution ratio |
| CN108549100A (en) * | 2018-01-11 | 2018-09-18 | 吉林大学 | The multiple dimensioned full waveform inversion method of time-domain of frequency is opened up based on non-linear high order |
| WO2019186286A1 (en) * | 2018-03-27 | 2019-10-03 | King Abdullah University Of Science And Technology | Robust full waveform inversion of seismic data method and device |
| WO2020009752A1 (en) * | 2018-07-02 | 2020-01-09 | Exxonmobil Upstream Research Company | Full wavefield inversion with an image-gather-flatness constraint |
| CN110058302A (en) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | A kind of full waveform inversion method based on pre-conditional conjugate gradient accelerating algorithm |
| CN110058307A (en) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | A kind of full waveform inversion method based on quick quasi-Newton method |
Non-Patent Citations (1)
| Title |
|---|
| 基于初始模型优化的稀疏高效全波形反演方法;段超然;《中国博士学位论文全文数据库 基础科学辑》;20181215(第12期);第2,4章 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111580163A (en) | 2020-08-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111580163B (en) | Full waveform inversion method and system based on non-monotonic search technology | |
| CN106526674B (en) | Three-dimensional full waveform inversion energy weighting gradient preprocessing method | |
| CN107843925B (en) | A kind of reflection wave inversion method based on orrection phase place | |
| WO2023087451A1 (en) | Observation data self-encoding-based multi-scale unsupervised seismic wave velocity inversion method | |
| CN105319589B (en) | A kind of fully automatic stereo chromatography conversion method using local lineups slope | |
| SG184803A1 (en) | Artifact reduction in method of iterative inversion of geophysical data | |
| CN109188519B (en) | System and method for inverting longitudinal and transverse wave speeds of elastic waves under polar coordinates | |
| CN106483559B (en) | A kind of construction method of subsurface velocity model | |
| CN110618453A (en) | Wave impedance inversion method based on improved damping least square method | |
| CN103454677B (en) | Based on the earthquake data inversion method that population is combined with linear adder device | |
| CN104237937B (en) | Pre-stack seismic inversion method and system thereof | |
| CN105319581A (en) | Efficient time domain full waveform inversion method | |
| CN104965222B (en) | Three-dimensional longitudinal wave impedance full-waveform inversion method and apparatus | |
| CN110879412A (en) | Underground transverse wave velocity inversion method, device, computing equipment and storage medium | |
| CN113486591B (en) | Gravity multi-parameter data density weighted inversion method for convolutional neural network result | |
| CN109541691B (en) | Seismic velocity inversion method | |
| CN109239776B (en) | Seismic wave propagation forward modeling method and device | |
| CN104730572A (en) | Diffracted wave imaging method and device based on L0 semi-norm | |
| CN117331119A (en) | Rapid earthquake wave travel time calculation method for tunnel detection | |
| CN111273346B (en) | Method, device, computer equipment and readable storage medium for removing deposition background | |
| CN111208568B (en) | Time domain multi-scale full waveform inversion method and system | |
| Wang et al. | Memory optimization in RNN-based full waveform inversion using boundary saving wavefield reconstruction | |
| CN110161565A (en) | A kind of Reconstruction of seismic data method | |
| CN113866827A (en) | Method, system, medium and device for explanatory velocity modeling seismic imaging | |
| CN110764146B (en) | Space cross-correlation elastic wave reflection waveform inversion method based on acoustic wave operator |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210810 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |