CN111208568B - Time domain multi-scale full waveform inversion method and system - Google Patents
Time domain multi-scale full waveform inversion method and system Download PDFInfo
- Publication number
- CN111208568B CN111208568B CN202010045893.7A CN202010045893A CN111208568B CN 111208568 B CN111208568 B CN 111208568B CN 202010045893 A CN202010045893 A CN 202010045893A CN 111208568 B CN111208568 B CN 111208568B
- Authority
- CN
- China
- Prior art keywords
- velocity model
- observation data
- iteration
- seismic observation
- model
- 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 78
- 238000004364 calculation method Methods 0.000 claims abstract description 20
- 238000003384 imaging method Methods 0.000 claims abstract description 19
- 238000001914 filtration Methods 0.000 claims abstract description 12
- 238000003325 tomography Methods 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 claims description 99
- 238000004088 simulation Methods 0.000 claims description 20
- 238000010276 construction Methods 0.000 claims description 4
- 238000009795 derivation Methods 0.000 claims description 3
- 230000036039 immunity Effects 0.000 abstract description 6
- 230000002708 enhancing effect Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 17
- 230000008569 process Effects 0.000 description 10
- 239000011159 matrix material Substances 0.000 description 9
- 108010001267 Protein Subunits Proteins 0.000 description 6
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 208000010392 Bone Fractures Diseases 0.000 description 2
- 206010017076 Fracture Diseases 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 241001482237 Pica Species 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000012795 verification 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
-
- 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/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
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 discloses a time domain multi-scale full waveform inversion method and a system, wherein the method comprises the following steps: step S1: acquiring earthquake observation data by using a geophone; step S2: determining an initial velocity model by using a tomography method; step S3: filtering the seismic observation data to different frequency bands by using a wiener low-pass filter; step S4: determining an optimal velocity model according to the initial velocity model and the seismic observation data of different frequency bands; step S5: and performing seismic imaging based on the optimal velocity model, improving the convergence velocity and the calculation efficiency of full-waveform inversion, and enhancing the noise immunity of the inversion.
Description
Technical Field
The invention relates to the technical field of seismic exploration, in particular to a time domain multi-scale full-waveform inversion method and a time domain multi-scale full-waveform inversion system.
Background
With the rapid development of social industry, the demand for resources is increasing, so the demand for deep resource detection in thick covering zones, severe undulating terrains (steep mountainous regions) and complex underground geological structures (underground fracture development, structural fracture, serious stratum deformation and complex lithology) is more and more strong. Seismic exploration is a system project, and whether a seismic method can be successfully used for resource exploration depends on seismic data processing technology to a great extent. The seismic imaging is a core link in seismic data processing and is the centralized embodiment of final results and benefits. With the development of computer hardware, seismic imaging techniques, represented by reverse time migration imaging (RTM), have been widely developed and applied. And a Full waveform inversion (Full waveform inversion) technology taking Full wavefield information fitting medium models as means can provide a high-resolution and high-precision velocity model for constructing the most complex area imaging.
Full waveform inversion starts from a wave equation, information such as amplitude, phase and the like in a pre-stack seismic wave field is comprehensively utilized, and an initial assumed velocity model is corrected by fitting an actual wave field and a predicted wave field to obtain high-precision underground velocity information. The full waveform inversion provides reliable basis for deep large-scale structural evolution analysis, accurate seismic imaging, velocity modeling and the like, so that the full waveform inversion is regarded as a high-precision velocity model construction method and is a research hotspot of the international geophysical field at present. Full waveform inversion is a highly nonlinear problem, and is sensitive to the quality of actual observed data (offset, effective low frequency and noise), technical means adopted in the inversion process (inversion algorithm, step size search method, objective function and inversion strategy) and the accuracy of an initial input model. At the same time, these factors all affect the application of full waveform inversion to practical seismic data to varying degrees.
The main objective of the full-waveform inversion is to minimize an objective function formed by observation data and synthetic seismic data, calculate gradient through an accompanying operator, calculate step length by using a certain search method, construct a search direction by combining a proper inversion algorithm, continuously update model parameters, and finally obtain a high-resolution parameter model. The objective function and the step search method are two factors that are crucial in full waveform inversion.
The objective function determines the degree of non-linearity and noise immunity of the full waveform inversion. The objective function of the conventional full waveform inversion is based on the L2 norm of the calculated data and observed data residuals, but the L2 norm has the problem of being sensitive to non-gaussian noise, and the full waveform inversion easily falls into local minima in the face of various strong noises present in the actual data.
The optimized step size search method can accelerate the convergence of full waveform inversion to a minimum value, few extra forward modeling is involved, and the calculation efficiency of full waveform inversion is improved, but the prior Pica et al (1990) proposes a method for directly solving the optimal step size from an objective function, but the method is only suitable for the L2 norm with poor noise immunity. Vigh et al (2009) propose to use a parabolic method, solve a parabolic coefficient by using three step values satisfying the relationship and objective function values corresponding to the three step values, and obtain an optimal step value. The non-precise linear search method ensures that the number of times of function estimation is small, but the search process usually needs multiple times of additional forward calculation, and the calculation cost is high.
Disclosure of Invention
Based on this, the invention aims to provide a time domain multi-scale full waveform inversion method and a time domain multi-scale full waveform inversion system, so as to improve the noise immunity, the convergence speed and the calculation efficiency of full waveform inversion.
To achieve the above object, the present invention provides a time domain multi-scale full waveform inversion method, including:
step S1: acquiring earthquake observation data by using a geophone;
step S2: determining an initial velocity model by using a tomography method;
step S3: filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, wherein the low-frequency band seismic observation data, the medium-frequency band seismic observation data and the high-frequency band seismic observation data are respectively low-frequency band seismic observation data, medium-frequency band seismic observation data and high-frequency band seismic observation data;
step S4: determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data;
step S5: and performing seismic imaging based on the optimal velocity model.
Optionally, the determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data, and the high-frequency band seismic observation data includes:
step S41: carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model;
step S42: carrying out inversion according to the middle-frequency-band seismic observation data and the first velocity model to obtain a second velocity model;
step S43: and carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model.
Optionally, the performing inversion according to the low-frequency band seismic observation data and the initial velocity model to obtain a first velocity model includes:
step S411: calculating simulated wave field data corresponding to the initial velocity model;
step S412: constructing an objective function based on the simulated wavefield data and the seismic observation data;
step S413: calculating the gradient of the objective function by adopting a adjoint state method;
step S414: determining the searching direction of the target function according to the gradient;
step S415: determining an optimal step size of the objective function;
step S416: determining a speed model of the next iteration according to the search direction and the optimal step length;
step S417: judging whether the iteration times are larger than or equal to a first iteration time threshold value or not; if the iteration number is smaller than the first iteration number threshold, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the iteration number is greater than or equal to the first iteration number threshold, taking the speed model of the next iteration as the first speed model, and executing step S42;
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 first set coefficient, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the relative difference of the successive objective functions is smaller than the first set coefficient, the velocity model of the next iteration is taken as the first velocity model, and "step S42" is executed.
Optionally, the calculating the simulated wave field data corresponding to the initial velocity model includes:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
The invention also provides a time domain multi-scale full waveform inversion system, which comprises:
the earthquake observation data acquisition module is used for acquiring earthquake observation data by using the geophone;
an initial velocity model determination module for determining an initial velocity model using a tomography method;
the filtering module is used for filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, and the seismic observation data are low-frequency-band seismic observation data, middle-frequency-band seismic observation data and high-frequency-band seismic observation data respectively;
the optimal velocity model determining module is used for determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data;
and the seismic imaging module is used for performing seismic imaging based on the optimal velocity model.
Optionally, the optimal speed model determining module includes:
the first velocity model determining unit is used for carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model;
the second velocity model determining unit is used for carrying out inversion according to the intermediate frequency band seismic observation data and the first velocity model to obtain a second velocity model;
and the optimal velocity model determining unit is used for carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model.
Optionally, the first speed model determining unit includes:
the first simulation wave field data determining subunit is used for calculating simulation wave field data corresponding to the initial velocity model;
an objective function construction subunit for constructing an objective function based on the simulated wavefield data and the seismic observation data;
a gradient determining subunit, configured to calculate a gradient of the objective function by using a adjoint state method;
a search direction determining subunit, configured to determine a search direction of the objective function according to the gradient;
an optimal step length determining subunit, configured to determine an optimal step length of the objective function;
the next iteration speed model determining subunit is used for determining the next iteration speed model according to the search direction and the optimal step length;
the first judgment subunit is used for judging whether the iteration times are greater than or equal to a first iteration time threshold value; if the iteration times are smaller than a first iteration time threshold value, taking a velocity model of the next iteration as the initial velocity model, and returning to the 'first simulation wave field data determination subunit'; if the iteration times are larger than or equal to the first iteration time threshold, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit;
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 first set coefficient, taking the velocity model of the next iteration as the initial velocity model, and returning to 'first simulation wave field data determination subunit'; and if the relative difference value of the continuous objective functions is smaller than a first set coefficient, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit.
Optionally, the determining the subunit by the first simulated wave field data specifically includes:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
According to the specific embodiment provided by the invention, the invention discloses the following technical effects:
the invention discloses a time domain multi-scale full waveform inversion method and a system, wherein the method comprises the following steps: step S1: acquiring earthquake observation data by using a geophone; step S2: determining an initial velocity model by using a tomography method; step S3: filtering the seismic observation data to different frequency bands by using a wiener low-pass filter; step S4: determining an optimal velocity model according to the initial velocity model and the seismic observation data of different frequency bands; step S5: and performing seismic imaging based on the optimal velocity model, improving the convergence velocity and the calculation efficiency of full-waveform inversion, and enhancing the noise immunity of the inversion.
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 time domain multi-scale full waveform inversion method according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of an initial velocity model according to an embodiment of the present invention;
FIG. 3 is an optimal velocity model obtained by inversion according to an embodiment of the present invention;
FIG. 4 is a schematic diagram of the optimal velocity model depth profiles obtained by inversion according to the embodiment of the present invention in horizontal positions, respectively;
FIG. 5 is a graph of the convergence of the objective function of the inversion process according to an embodiment of the present invention;
FIG. 6 is a schematic diagram of seismic survey data distribution according to an embodiment of the invention;
FIG. 7 is an optimal velocity model obtained by inversion after noise is added in an embodiment of the present invention;
FIG. 8 is a schematic diagram of the optimal velocity model depth profiles obtained by inversion after noise is added, respectively in the horizontal position in the embodiment of the present invention;
FIG. 9 is a graph of the convergence of the objective function of the inversion process after noise is added in accordance with an embodiment of the present invention;
fig. 10 is a structural diagram of a time domain multi-scale full waveform inversion system according to an embodiment of the 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 time domain multi-scale full waveform inversion method and a time domain multi-scale full waveform inversion system, so as to improve the noise immunity, the convergence speed and the calculation efficiency 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 time domain multi-scale full waveform inversion method according to an embodiment of the present invention, and as shown in fig. 1, the present invention provides a time domain multi-scale full waveform inversion method, where the method includes:
step S1: and acquiring seismic observation data by using a detector.
Step S2: an initial velocity model is determined using a tomographic method.
Step S3: and filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, wherein the low-frequency band seismic observation data, the medium-frequency band seismic observation data and the high-frequency band seismic observation data are respectively low-frequency band seismic observation data, medium-frequency band seismic observation data and high-frequency band seismic observation data.
Step S4: and determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data.
Step S5: and performing seismic imaging based on the optimal velocity model.
The individual steps are discussed in detail below:
step S1: and acquiring seismic observation data by using a detector.
And exciting a seismic source in the field, and acquiring seismic observation data by using a detector, wherein the seismic observation data can be seismic observation data added with noise or seismic observation data without noise.
Step S2: an initial velocity model is determined using a tomographic method.
The strong nonlinear relation between the data residual and the model parameters in the full-waveform inversion causes that the inversion process is easy to fall into a local minimum value, and the initial velocity model directly determines whether the inversion is successful or not, so that the chromatographic imaging method is used for determining the initial velocity model, and the phenomenon that the inversion process falls into the local minimum value is prevented.
Step S3: and filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, wherein the low-frequency band seismic observation data, the medium-frequency band seismic observation data and the high-frequency band seismic observation data are respectively low-frequency band seismic observation data, medium-frequency band seismic observation data and high-frequency band seismic observation data.
The range of each frequency band is determined according to the frequency spectrum distribution of actual data, and low frequency bands and high frequency bands of different regions or different data can be different without a uniform standard.
According to the invention, the low-frequency-band seismic observation data can construct a macrostructural velocity model, then the low-frequency-band seismic observation data is input into the middle-frequency-band inversion process and then input into the high-frequency-band inversion process, and with the continuous improvement of the inversion frequency band, the model can be finely depicted, so that a high-precision velocity model is obtained. The method comprises the following specific steps:
step S4: the determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data comprises the following steps:
step S41: and carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model.
Step S42: and carrying out inversion according to the intermediate frequency band seismic observation data and the first velocity model to obtain a second velocity model.
Step S43: and carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model.
Wherein, step S41: performing inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model, wherein the inversion comprises the following steps:
step S411: and calculating the simulated wave field data corresponding to the initial velocity model.
The method solves the equation (1) by using a finite difference method, is a widely used numerical simulation method, and is simple in solving, easy to realize programming, small in required memory and high in calculation efficiency by using a local operator. The seismic source forward transmission and residual reverse transmission related by the method are carried out by adopting a finite difference format of 2-order time and 10-order space, and the strong reflection generated by artificial boundaries is absorbed by adopting a complete matching layer absorption boundary condition (PML) in numerical simulation to finish the seismic source forward transmission, and the method comprises the following specific steps of:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
The acoustic wave equation is:
wherein t is the wave field propagation time, x and z are the horizontal coordinate and the vertical coordinate of the wave field respectively,
for the simulated wavefield data of the kth iteration, vkFor the velocity model for the kth iteration, s (x, z, t) is the source term.
Step S412: constructing an objective function based on the simulated wavefield data and the seismic observation data; the objective function is an L1 norm type objective function; the objective function is:
for convenience of calculation, let in equation (2):
wherein s is the number of shots, l is the number of detectors,
for the simulated wavefield data of the kth iteration, pk obs(x, z, t) is seismic observation data for the kth iteration.
Step S413: calculating the gradient of the objective function by adopting a adjoint state method; the gradient formula is:
wherein, gkAnd the gradient of the kth iteration, T is the maximum calculation time, and lambda is the residual back-propagation wave field at the detector.
The quasi-Newton algorithm is characterized in that the inverse of the Hessian matrix does not need to be directly calculated (the calculation cost is high), but a certain means is adopted to approximately construct the inverse of the Hessian matrix. The finite memory Newton method (L-BFGS) avoids storing an approximate hessian matrix inverse matrix, and constructs a search direction by storing finite gradient variation and speed model updating quantity (3-20), so that the following steps are disclosed:
step S414: determining the searching direction of the target function according to the gradient; specifically, a double-loop iteration method is adopted, and the search direction of the target function is determined according to the gradient; the search direction formula is as follows:
dk＝-Hkgk (5)
wherein, gkGradient for the k-th iteration, dkFor the search direction of the kth iteration, HkIs the inverse of the Hessian matrix in the kth iteration, HkThe expression of (a) is as follows:
wherein the content of the first and second substances,
i is an identity matrix and is a matrix of the identity,
ykgradient change, y, for the k-th iterationk＝gk+1-gk，skFor the speed of the k-th iterationThe amount of the degree update is calculated,
the initial matrix is approximate to the hessian inverse matrix, and m is the number of gradient change amounts stored, and is usually 3-20.
Step S415: and determining the optimal step size of the objective function.
Taylor expansion is performed on the simulated wave field data to obtain:
wherein p iscal(vk) To simulate the observed wave field, αkFor the optimal step size of the kth iteration, vkThe velocity model representing the kth iteration, i.e. the optimal velocity model,
is the laplacian operator.
The objective function (2) can thus be written as:
wherein p isobs(vk) For the actual observed wavefield;
to find the optimal objective function, the objective function is made to derive the variable step size and the derivative is made equal to 0, i.e.:
through the calculation and derivation of the above series of formulas, the optimal step calculation formula for the L1 norm type objective function can be obtained as follows:
wherein alpha istFor measuring step length, c1Is a value of the weight, and is,
step S416: and determining a speed model of the next iteration according to the search direction and the optimal step length, wherein the specific formula is as follows:
vk+1＝vk+αkdk (11)
wherein v isk+1Velocity model for the k +1 th iteration, vkVelocity model for the kth iteration, αkFor the optimal step size of the kth iteration, dkThe search direction for the kth iteration.
Step S417: judging whether the iteration times are larger than or equal to a first iteration time threshold value or not; if the iteration number is smaller than the first iteration number threshold, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the iteration number is greater than or equal to the first iteration number threshold, taking the speed model of the next iteration as the first speed model, and executing step S42;
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 first set coefficient, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the relative difference of the successive objective functions is smaller than the first set coefficient, the velocity model of the next iteration is taken as the first velocity model, and "step S42" is executed.
Step S42: and carrying out inversion according to the intermediate frequency band seismic observation data and the first velocity model to obtain a second velocity model.
Step S421: calculating the simulated wave field data corresponding to the first velocity model, specifically comprising:
and calculating the simulated wave field data corresponding to the first speed model by using an acoustic wave equation by adopting a finite difference method.
Steps S422 to S426 are completely the same as steps S412 to S416, and are not described in detail herein.
Step S427: judging whether the iteration times are larger than or equal to a second iteration time threshold value or not; if the iteration number is smaller than the second iteration number threshold, taking the speed model of the next iteration as the first speed model, and returning to the step S421; if the number of iterations is greater than or equal to the second iteration number threshold, the speed model of the next iteration is taken as the second speed model, and "step S43" is performed.
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 second set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a second set coefficient, taking the speed model of the next iteration as the first speed model, and returning to the step S421; if the relative difference of the successive objective functions is smaller than the second set coefficient, the speed model of the next iteration is taken as the second speed model, and step S43 is executed.
Step S43: performing inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model, wherein the inversion comprises the following steps:
step S431: calculating the simulated wave field data corresponding to the second velocity model, specifically comprising:
and calculating the simulated wave field data corresponding to the second velocity model by using an acoustic wave equation by using a finite difference method.
Steps S432 to S436 are completely the same as steps S412 to S416, and are not described again.
Step S437: judging whether the iteration times are larger than or equal to a third iteration time threshold value or not; if the iteration times are smaller than the third iteration time threshold, taking the speed model of the next iteration as the second speed model, and returning to the step S431; if the number of iterations is greater than or equal to the third iteration number threshold, the velocity model of the next iteration is taken as the optimal velocity model, and "step S5" is performed.
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 third set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a third set coefficient, taking the speed model of the next iteration as the second speed model, and returning to the step S431; if the relative difference of the successive objective functions is smaller than the third setting coefficient, the velocity model of the next iteration is taken as the optimal velocity model, and step S5 is executed.
Fig. 10 is a structural diagram of a time domain multi-scale full waveform inversion system according to an embodiment of the present invention, and as shown in fig. 10, the present invention further provides a time domain multi-scale full waveform inversion system, where the system includes:
the earthquake observation data acquisition module 1 is used for acquiring earthquake observation data by using a wave detector;
an initial velocity model determining module 2, for determining an initial velocity model by using a tomography method;
the filtering module 3 is used for filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, wherein the frequency bands are low-frequency band seismic observation data, middle-frequency band seismic observation data and high-frequency band seismic observation data;
the optimal velocity model determining module 4 is used for determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data;
and the seismic imaging module 5 is used for performing seismic imaging based on the optimal velocity model.
As an embodiment, the optimal speed model determining module 4 of the present invention includes:
and the first velocity model determining unit is used for carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model.
And the second velocity model determining unit is used for carrying out inversion according to the intermediate frequency band seismic observation data and the first velocity model to obtain a second velocity model.
And the optimal velocity model determining unit is used for carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model.
As an embodiment, the first velocity model determining unit of the present invention includes:
the first simulation wave field data determining subunit is used for calculating simulation wave field data corresponding to the initial velocity model;
an objective function construction subunit for constructing an objective function based on the simulated wavefield data and the seismic observation data;
a gradient determining subunit, configured to calculate a gradient of the objective function by using a adjoint state method;
a search direction determining subunit, configured to determine a search direction of the objective function according to the gradient;
an optimal step length determining subunit, configured to determine an optimal step length of the objective function;
the next iteration speed model determining subunit is used for determining the next iteration speed model according to the search direction and the optimal step length;
the first judgment subunit is used for judging whether the iteration times are greater than or equal to a first iteration time threshold value; if the iteration times are smaller than a first iteration time threshold value, taking a velocity model of the next iteration as the initial velocity model, and returning to the 'first simulation wave field data determination subunit'; if the iteration times are larger than or equal to the first iteration time threshold, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit;
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 first set coefficient, taking the velocity model of the next iteration as the initial velocity model, and returning to 'first simulation wave field data determination subunit'; and if the relative difference value of the continuous objective functions is smaller than a first set coefficient, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit.
As an embodiment, the determining subunit of the first simulated wave field data specifically includes:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
As an embodiment, the second velocity model determining unit of the present invention includes:
the second simulation wave field data determining subunit is used for calculating the simulation wave field data corresponding to the first speed model;
the second judgment subunit is used for judging whether the iteration times are greater than or equal to a second iteration time threshold value; if the iteration times are smaller than a second iteration time threshold value, taking a speed model of the next iteration as the first speed model, and returning to a 'second simulation wave field data determining subunit'; if the iteration times are larger than or equal to the second iteration time threshold, taking the speed model of the next iteration as a second speed model, and executing a third speed model determining unit;
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 second set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a second set coefficient, taking the velocity model of the next iteration as the first velocity model, and returning to a 'second simulation wave field data determining subunit'; and if the relative difference value of the continuous objective functions is smaller than a second set coefficient, taking the speed model of the next iteration as a second speed model, and executing a third speed model determining unit.
The remaining sub-units in the second velocity model determining unit are completely the same as the target function constructing sub-unit, the gradient determining sub-unit, the search direction determining sub-unit, the optimal step length determining sub-unit and the velocity model determining sub-unit of the next iteration in the first velocity model determining unit, and are not described in detail herein.
As an embodiment, the second simulated wave field data determining subunit specifically includes:
and calculating the simulated wave field data corresponding to the first speed model by using an acoustic wave equation by adopting a finite difference method.
As an embodiment, the optimal velocity model determining unit according to the present invention includes:
the third simulated wave field data determining subunit is used for calculating the simulated wave field data corresponding to the second velocity model;
a third judging subunit, configured to judge whether the iteration number is greater than or equal to a third iteration number threshold; if the iteration times are smaller than a third iteration time threshold value, taking the speed model of the next iteration as the second speed model, and returning to a third simulation wave field data determination subunit; if the iteration times are larger than or equal to a third iteration time threshold value, taking the velocity model of the next iteration as an optimal velocity model, and executing a seismic imaging 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 third set coefficient; if the relative difference value of the continuous objective functions is larger than or equal to a third set coefficient, taking the speed model of the next iteration as the second speed model, and returning to 'third simulation wave field data determination subunit'; and if the relative difference value of the continuous objective functions is smaller than a third set coefficient, taking the velocity model of the next iteration as an optimal velocity model, and executing a seismic imaging module.
The remaining sub-units in the optimal speed model determining unit and the target function constructing sub-unit, the gradient determining sub-unit, the search direction determining sub-unit, the optimal step length determining sub-unit and the speed model determining sub-unit of the next iteration in the first speed model determining unit are complete, and are not described in detail herein.
As an embodiment, the determining subunit of the third simulated wave field data specifically includes:
and calculating the simulated wave field data corresponding to the second velocity model by using an acoustic wave equation by using a finite difference method.
The time domain multi-scale full-waveform inversion method and the time domain multi-scale full-waveform inversion system disclosed by the invention are adopted for experimental verification, and FIG. 2 is a schematic diagram of an initial velocity model in the embodiment of the invention; FIGS. 2(a) and 2(b) are an Overthrast true model and a smoothed initial model, respectively; FIG. 3 is an optimal velocity model obtained by inversion according to an embodiment of the present invention, and it can be found that the velocity model obtained by inversion according to the present invention has a clear structure, and both the propagation volume and the layered model are clear; FIG. 4 is a schematic diagram of horizontal positions of three depth profiles arbitrarily extracted from an optimal velocity model obtained by inversion according to an embodiment of the present invention; (a) the method comprises the following steps of (a) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 3.0km for inversion, (b) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 4.8km for inversion, and (c) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 6.0km for inversion, wherein the depth profile obtained by the method is well matched with a real depth profile; FIG. 5 is a graph of the convergence of the objective function of the inversion process according to an embodiment of the present invention; FIG. 6 is a schematic diagram of distribution of seismic observation data according to an embodiment of the present invention, which shows that the convergence of the objective function corresponding to the present invention is fast, and the inversion error is small and the accuracy is high under the condition of fewer iteration times; FIG. 6 (a) and (b) are schematic diagrams showing the distribution of seismic observation data of shot 30 of the Overhaust model and seismic observation data added with strong noise, respectively; FIG. 7 is an optimal velocity model obtained by inversion after noise is added in the embodiment of the present invention, and it can be found that in the case of strong noise mixed in data, although the inversion accuracy is reduced compared with the case of no noise, the present invention can still obtain a velocity model with high accuracy; FIG. 8 is a schematic diagram of horizontal positions of three depth profiles arbitrarily extracted from an optimal velocity model obtained by inversion after noise is added in the embodiment of the invention; (a) the method comprises the following steps of (a) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 3.0km after noise is added, (b) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 4.8km after noise is added, and (c) obtaining a schematic diagram of an optimal velocity model depth profile at a horizontal position of 6.0km after noise is added, wherein the method can find that the integral matching degree with a real velocity model is good, and the velocities of a deep layer and a shallow layer are well restored and reconstructed; FIG. 9 is a graph of the convergence curve of the objective function in the inversion process after noise is added, and it can be found that, under the condition of adding strong noise, the convergence speed of the objective function curve is faster on the premise of ensuring the accuracy; fig. 2-9 demonstrate the robustness of the method of the invention.
The invention provides a time domain multi-scale full waveform inversion method and a time domain multi-scale full waveform inversion system, wherein an L1 norm type target function is used in full waveform inversion, the anti-noise performance of inversion is improved, and a satisfactory inversion result can still be obtained under the condition that seismic data contain strong noise. Aiming at the L1 norm type target function, an applicable optimal step size searching method is provided, so that the convergence speed and the calculation efficiency of full waveform inversion are improved, and the anti-noise property of the inversion is enhanced.
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 (4)
1. A method of time domain multi-scale full waveform inversion, the method comprising:
step S1: acquiring earthquake observation data by using a geophone;
step S2: determining an initial velocity model by using a tomography method;
step S3: filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, wherein the low-frequency band seismic observation data, the medium-frequency band seismic observation data and the high-frequency band seismic observation data are respectively low-frequency band seismic observation data, medium-frequency band seismic observation data and high-frequency band seismic observation data;
step S4: determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data;
step S5: performing seismic imaging based on the optimal velocity model;
the determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data comprises the following steps:
step S41: carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model;
step S42: carrying out inversion according to the middle-frequency-band seismic observation data and the first velocity model to obtain a second velocity model;
step S43: carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model;
and performing inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model, wherein the inversion comprises the following steps:
step S411: calculating simulated wave field data corresponding to the initial velocity model;
step S412: constructing an objective function based on the simulated wavefield data and the seismic observation data;
step S413: calculating the gradient of the objective function by adopting a adjoint state method;
step S414: determining the searching direction of the target function according to the gradient;
step S415: determining an optimal step size of the objective function;
taylor expansion is performed on the simulated wave field data to obtain:
wherein p iscal(vk) To simulate the observed wave field, αkFor the optimal step size of the kth iteration, vkThe velocity model representing the kth iteration, i.e. the optimal velocity model,
is Laplace operator, dkThe search direction for the kth iteration;
thus the objective function
Can be written as:
wherein p isobs(vk) For the actual observed wavefield, s is the number of shots, l is the number of detectors,
for the simulated wavefield data of the kth iteration, pk obs(x, z, T) is seismic observation data of the kth iteration, and T is maximum calculation time;
to find the optimal objective function, the objective function is made to derive the variable step size and the derivative is made equal to 0, i.e.:
through the calculation and derivation of the above series of formulas, the optimal step calculation formula for the L1 norm type objective function can be obtained as follows:
wherein alpha istTo testStep size, c1Is a value of the weight, and is,
step S416: determining a speed model of the next iteration according to the search direction and the optimal step length;
step S417: judging whether the iteration times are larger than or equal to a first iteration time threshold value or not; if the iteration number is smaller than the first iteration number threshold, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the iteration number is greater than or equal to the first iteration number threshold, taking the speed model of the next iteration as the first speed model, and executing step S42;
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 first set coefficient, taking the speed model of the next iteration as the initial speed model, and returning to the step S411; if the relative difference of the successive objective functions is smaller than the first set coefficient, the velocity model of the next iteration is taken as the first velocity model, and "step S42" is executed.
2. The time domain multi-scale full waveform inversion method according to claim 1, wherein the calculating the simulated wavefield data corresponding to the initial velocity model comprises:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
3. A time domain multi-scale full waveform inversion system, the system comprising:
the earthquake observation data acquisition module is used for acquiring earthquake observation data by using the geophone;
an initial velocity model determination module for determining an initial velocity model using a tomography method;
the filtering module is used for filtering the seismic observation data to different frequency bands by using a wiener low-pass filter, and the seismic observation data are low-frequency-band seismic observation data, middle-frequency-band seismic observation data and high-frequency-band seismic observation data respectively;
the optimal velocity model determining module is used for determining an optimal velocity model according to the initial velocity model, the low-frequency band seismic observation data, the middle-frequency band seismic observation data and the high-frequency band seismic observation data;
the seismic imaging module is used for carrying out seismic imaging based on the optimal velocity model;
the optimal speed model determination module comprises:
the first velocity model determining unit is used for carrying out inversion according to the low-frequency-band seismic observation data and the initial velocity model to obtain a first velocity model;
the second velocity model determining unit is used for carrying out inversion according to the intermediate frequency band seismic observation data and the first velocity model to obtain a second velocity model;
the optimal velocity model determining unit is used for carrying out inversion according to the high-frequency-band seismic observation data and the second velocity model to obtain an optimal velocity model;
the first velocity model determination unit includes:
the first simulation wave field data determining subunit is used for calculating simulation wave field data corresponding to the initial velocity model;
an objective function construction subunit for constructing an objective function based on the simulated wavefield data and the seismic observation data;
a gradient determining subunit, configured to calculate a gradient of the objective function by using a adjoint state method;
a search direction determining subunit, configured to determine a search direction of the objective function according to the gradient;
an optimal step length determining subunit, configured to determine an optimal step length of the objective function;
taylor expansion is performed on the simulated wave field data to obtain:
wherein p iscal(vk) To simulate the observed wave field, αkFor the optimal step size of the kth iteration, vkThe velocity model representing the kth iteration, i.e. the optimal velocity model,
is Laplace operator, dkThe search direction for the kth iteration;
thus the objective function
Can be written as:
wherein p isobs(vk) For the actual observed wavefield, s is the number of shots, l is the number of detectors,
for the simulated wavefield data of the kth iteration, pk obs(x, z, T) is seismic observation data of the kth iteration, and T is maximum calculation time;
to find the optimal objective function, the objective function is made to derive the variable step size and the derivative is made equal to 0, i.e.:
through the calculation and derivation of the above series of formulas, the optimal step calculation formula for the L1 norm type objective function can be obtained as follows:
wherein alpha istFor measuring step length, c1Is a value of the weight, and is,
the next iteration speed model determining subunit is used for determining the next iteration speed model according to the search direction and the optimal step length;
the first judgment subunit is used for judging whether the iteration times are greater than or equal to a first iteration time threshold value; if the iteration times are smaller than a first iteration time threshold value, taking a velocity model of the next iteration as the initial velocity model, and returning to the 'first simulation wave field data determination subunit'; if the iteration times are larger than or equal to the first iteration time threshold, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit;
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 first set coefficient, taking the velocity model of the next iteration as the initial velocity model, and returning to 'first simulation wave field data determination subunit'; and if the relative difference value of the continuous objective functions is smaller than a first set coefficient, taking the speed model of the next iteration as a first speed model, and executing a second speed model determining unit.
4. The time domain multi-scale full waveform inversion system of claim 3, wherein the first simulated wavefield data determining subunit comprises:
and calculating the simulated wave field data corresponding to the initial velocity model by using an acoustic wave equation by adopting a finite difference method.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010045893.7A CN111208568B (en) | 2020-01-16 | 2020-01-16 | Time domain multi-scale full waveform inversion method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010045893.7A CN111208568B (en) | 2020-01-16 | 2020-01-16 | Time domain multi-scale full waveform inversion method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111208568A CN111208568A (en) | 2020-05-29 |
| CN111208568B true CN111208568B (en) | 2021-04-09 |
Family
ID=70787266
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010045893.7A Expired - Fee Related CN111208568B (en) | 2020-01-16 | 2020-01-16 | Time domain multi-scale full waveform inversion method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111208568B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112379415A (en) * | 2020-11-06 | 2021-02-19 | 中国科学技术大学 | Multi-scale full-waveform inversion method and device for reconstructed low-frequency data based on downsampling |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105319581B (en) * | 2014-07-31 | 2018-01-16 | 中国石油化工股份有限公司 | A kind of efficient time-domain full waveform inversion method |
| WO2016165064A1 (en) * | 2015-04-14 | 2016-10-20 | 中国科学院自动化研究所 | Robust foreground detection method based on multi-view learning |
| CN105005076B (en) * | 2015-06-02 | 2017-05-03 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105549079A (en) * | 2016-01-12 | 2016-05-04 | 中国矿业大学(北京) | Method and device for establishing full-waveform inversion model for geophysics parameters |
| US10436926B2 (en) * | 2016-08-17 | 2019-10-08 | Pgs Geophysical As | Marine vibrator source acceleration and pressure |
| CN106501852B (en) * | 2016-10-21 | 2018-06-08 | 中国科学院地质与地球物理研究所 | A kind of multiple dimensioned full waveform inversion method of three-dimensional acoustic wave equation arbitrarily-shaped domain and device |
| CN107765302B (en) * | 2017-10-20 | 2018-06-26 | 吉林大学 | Inversion method when the time-domain single-frequency waveform of source wavelet is walked is not depended on |
| CN107894618B (en) * | 2017-11-10 | 2018-08-21 | 中国海洋大学 | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm |
-
2020
- 2020-01-16 CN CN202010045893.7A patent/CN111208568B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN111208568A (en) | 2020-05-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106526674B (en) | Three-dimensional full waveform inversion energy weighting gradient preprocessing method | |
| CA2726462C (en) | Estimation of soil properties using waveforms of seismic surface waves | |
| KR101549388B1 (en) | Prestack elastic generalized-screen migration method for seismic multicomponent data | |
| CN107843925B (en) | A kind of reflection wave inversion method based on orrection phase place | |
| Vigh et al. | Developing earth models with full waveform inversion | |
| Raknes et al. | Time-lapse full-waveform inversion of limited-offset seismic data using a local migration regularization | |
| Li et al. | Elastic reflection waveform inversion with variable density | |
| CA2795340A1 (en) | Artifact reduction in iterative inversion of geophysical data | |
| CN113740901B (en) | Land seismic data full-waveform inversion method and device based on complex undulating surface | |
| US20130258810A1 (en) | Method and System for Tomographic Inversion | |
| EA032186B1 (en) | Seismic adaptive focusing | |
| KR20170118185A (en) | Multi-stage full wave inversion process to generate multiple sets of nonreflective data sets | |
| CN105093278A (en) | Extraction method for full waveform inversion gradient operator based on excitation main energy optimization algorism | |
| CN111580163B (en) | Full waveform inversion method and system based on non-monotonic search technology | |
| CN109655890B (en) | Depth domain shallow-medium-deep layer combined chromatography inversion speed modeling method and system | |
| CN109946742A (en) | The pure rolling land qP shakes digital simulation method in a kind of TTI medium | |
| CN104749625B (en) | A kind of geological data inclination angle method of estimation based on Regularization Technique and device | |
| Hu et al. | Ray-illumination compensation for adjoint-state first-arrival traveltime tomography | |
| Gao et al. | An efficient vector elastic reverse time migration method in the hybrid time and frequency domain for anisotropic media | |
| CN111208568B (en) | Time domain multi-scale full waveform inversion method and system | |
| Li et al. | Waveform inversion of seismic first arrivals acquired on irregular surface | |
| Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
| CN111435174B (en) | Method and device for compensating amplitude of seismic data in strong reflection area | |
| CN114814944B (en) | Elastic longitudinal and transverse wave field separation and reverse time migration imaging method based on divergence and rotation | |
| Wang et al. | Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20210409 Termination date: 20220116 |