CN110133718A - A kind of attenuation anisotropy elasticity of fluid impedance inversion approach - Google Patents
A kind of attenuation anisotropy elasticity of fluid impedance inversion approach Download PDFInfo
- Publication number
- CN110133718A CN110133718A CN201910407041.5A CN201910407041A CN110133718A CN 110133718 A CN110133718 A CN 110133718A CN 201910407041 A CN201910407041 A CN 201910407041A CN 110133718 A CN110133718 A CN 110133718A
- Authority
- CN
- China
- Prior art keywords
- formula
- parameter
- crack
- attenuation
- decaying
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000012530 fluid Substances 0.000 title claims abstract description 54
- 238000013459 approach Methods 0.000 title claims abstract description 14
- 239000011435 rock Substances 0.000 claims abstract description 54
- 238000000034 method Methods 0.000 claims abstract description 23
- 229920006395 saturated elastomer Polymers 0.000 claims abstract description 13
- 238000012512 characterization method Methods 0.000 claims description 29
- 239000011159 matrix material Substances 0.000 claims description 20
- 239000013598 vector Substances 0.000 claims description 15
- 230000006698 induction Effects 0.000 claims description 14
- 230000008859 change Effects 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 102000008297 Nuclear Matrix-Associated Proteins Human genes 0.000 claims description 3
- 108010035916 Nuclear Matrix-Associated Proteins Proteins 0.000 claims description 3
- 239000000203 mixture Substances 0.000 claims description 3
- 210000000299 nuclear matrix Anatomy 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 230000010287 polarization Effects 0.000 claims description 3
- 230000008569 process Effects 0.000 claims description 3
- 239000007787 solid Substances 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 230000008030 elimination Effects 0.000 abstract description 4
- 238000003379 elimination reaction Methods 0.000 abstract description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 abstract description 3
- 208000010392 Bone Fractures Diseases 0.000 description 5
- 206010017076 Fracture Diseases 0.000 description 5
- 230000015572 biosynthetic process Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 241000208340 Araliaceae Species 0.000 description 2
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 2
- 235000003140 Panax quinquefolius Nutrition 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 235000008434 ginseng Nutrition 0.000 description 2
- 238000001615 p wave Methods 0.000 description 2
- 239000011148 porous material Substances 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000007423 decrease Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
-
- 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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
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 kind of attenuation anisotropy elasticity of fluid impedance inversion approach, comprising the following steps: step 1: obtaining the linear slide HTI model in view of the intrinsic decaying of background rock and crack induced attenuation；Step 2: the relation equation between the dry elasticity modulus of fissure rock and saturated flow bulk modulus is established；Step 3: saturation fluid apertures gap fissuted medium Complex modes real part is obtained using the formula in step 1 and step 2；Step 4: the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is obtained；Step 5: longitudinal wave reflection coefficient equation is obtained；Step 6: decaying elastic impedance equation in orientation is obtained by the longitudinal wave reflection coefficient equation of step 5；Step 7: inverting is carried out to crack " weakness " parameter and crack induced attenuation parameter by the decaying elastic impedance difference of different direction.Inversion method of the present invention realizes the pre-stack seismic inversion of crack elimination background fluid parameter, fracture water flow parameter and attenuation parameter.
Description
Technical field
The present invention relates to petroleum resources geophysics exploratory techniques fields, and in particular to a kind of attenuation anisotropy fluid bullet
Property impedance inversion approach.
Background technique
As to the continuous growth of Demand of Oil & Gas, the oil-gas exploration direction of Geophysicist is also gradually in global range
Some special oil and gas reservoirs are turned to by conventional gas and oil reservoir, wherein the exploration and development of slit formation oil and gas reservoir is ground as oil-gas exploration
The focus studied carefully.Due to the combined influence by factors such as strata pressures, most cracks are all high angles or close in actual formation
Vertical fracture, when fracture development is to certain scale, the propagation of seismic wave and its seismic response will appear cross in crack elimination
Wavefront splitting, azimuthal velocity anisotropic character etc..The list along advantage horizontal direction is developed in homogeneous isotropism background media
Group vertical fracture can be considered the equivalent HTI medium of long wavelength.Therefore, underground medium crack seismic response spy is described using seismic data
Sign, there is important directive significance and real value to the exploration and development of crack type reservoir.
Stronger velocity dispersion and decaying can occur when practical fractured reservoir is propagated for seismic wave, and elastic response has frequency
Rate dependence.Previous research or the only intrinsic decaying of consideration background media or only consideration matrix pores and fluid flowing
The influence of effect；And the application will consider that matrix pores and fluid flowing effect and the intrinsic decaying of background and crack induction decline simultaneously
Subtract the influence to anisotropy poroelasticity approximation theory, establishes a kind of anisotropic fluid elastic impedance inversion method.
Summary of the invention
The purpose of the present invention is to overcome above-mentioned the deficiencies in the prior art, provide a kind of attenuation anisotropy elasticity of fluid resistance
Anti- inversion method.
To achieve the above object, the present invention adopts the following technical solutions: a kind of attenuation anisotropy elasticity of fluid impedance is anti-
Drill method, comprising the following steps:
Step 1: the linear slide HTI model in view of the intrinsic decaying of background rock and crack induced attenuation is obtained；
Step 2: the relation equation between the dry elasticity modulus of fissure rock and saturated flow bulk modulus is established；
Step 3: it is real that saturation fluid apertures gap fissuted medium Complex modes are obtained using the formula in step 1 and step 2
Portion；
Step 4: the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is obtained；
Step 5: longitudinal wave reflection coefficient equation is obtained；
Step 6: decaying elastic impedance equation in orientation is obtained by the longitudinal wave reflection coefficient equation of step 5；
Step 7: by different direction decaying elastic impedance difference to crack " weakness " parameter and crack induced attenuation parameter into
Row inverting.
Preferably, in the step 1, isotropism background modulus parameter and crack ginseng in the linear slide HTI model
Number is represented as plural form, and formula is as follows:
In formula (1), WithRespectively indicate the multiple longitudinal wave mould of Isotropic Decay background rock
Amount, the first and second Lame constants, and
WithIt indicates
The normal direction of plural form and tangential crack " weakness " parameter can be used along the crack of symmetrical axis direction induction P-wave And S decaying
WithIt is represented simply as
Under the glutinous assumed condition for playing background rock of isotropism, the model parameter of plural formWithIt can indicate
For
In formula (3),WithBackground P-wave And S is respectively indicated against quality factor；
Publicity (2), (3) are brought into formula (1) simultaneously and are obtained in view of the intrinsic decaying of background rock and crack
The linear slide HTI model of induced attenuation.
Preferably, the pass in the step 2, between the dry elasticity modulus and saturated flow bulk modulus of the fissure rock
It is that equation is as follows:
In formula (4),It is the saturation fluid fissure rock effective rigidity elastic matrix of plural number characterization,It is plural number
The desiccation crack rock effective rigidity elastic matrix of characterization；It is the anisotropic class Biot coefficient of plural number characterization, can indicates
For
Wherein, KgIndicate the effective bulk modulus of composition rock solid particle；It is the anisotropic class of plural number characterization
Gassmann interstitial space modulus, is represented by
In formula,It is the dry rock volume modulus of anisotropic class of plural number characterization；φ indicates hole
Porosity；κfIndicate pore-fluid effective bulk modulus.
Preferably, in the step 3, the method for obtaining saturation fluid apertures gap fissuted medium Complex modes real part is as follows:
Wherein, the real part of the desiccation crack rock weak anisotropy stiffness matrix of plural number characterization is represented by
In formula (7),
WithThe intrinsic P-wave And S decaying of characterization background and crack induction are vertical respectively
The overall attenuation factor of wave attenuation product；
Formula (5), (6), (7) are substituted into formula (4), and close based on weak anisotropy approximation and the intrinsic decaying of weak background
Seemingly, then within the scope of seismic frequency, the dry Modulus of Rocks of background, desiccation crack flexibility, fluid modulus, the intrinsic P-wave And S of background are obtained
Decay and the approximate expression of the saturation fluid apertures gap fissuted medium Complex modes real part of crack induction attenuation of P-wave characterization is
In formula (8),
β0=1-Kdry/KgIt is Biot coefficient；
KdryIndicate the bulk modulus of dry rock.
Preferably, in the step 4, the method for obtaining PP wave reflection coefficient equation is as follows:
Based on scattering theory, the longitudinal wave reflection coefficient of HTI medium be may be expressed as:
Wherein, θ indicates incidence angle, and ρ indicates the density item of homogeneous isotropism background media,Indicate saturated rock
Disturbance stiffness parameters, ξmnIt is related with slowness vector sum polarization vector；
Therefore, convolution (8) and (9) can derive the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium, table
It is up to formula
In formula,
In formula (10),Indicate azimuth, Δ indicates the difference of each attribute of two layers of medium, f=fluid/hole
Gap is a fluid factor；Crack " weakness " parameter can be estimated according to log data and petrophysical model；
Preferably, in the step 5, according to horizontal interface longitudinal wave reflection coefficient and azimuthal anisotropy elastic impedance it
Between relationship, longitudinal wave reflection coefficient equation may be expressed as:
In formula (11),Indicate that orientation decays elastic impedance, Δ EI (θ, φ) indicates upper one layer and next
The difference of layer elastic impedance.
Preferably, in the step 6, orientation decaying elastic impedance side is obtained by the longitudinal wave reflection coefficient equation of step 5
The method of journey is as follows:
Weak elastic parameter it is poor (| Δ f/f | < < 1, | Δ μ/μb| < < 1 and | Δ ρ/ρb| < < 1), rimala " weakness " (And δT< < 1) and underdamp (And) under assumed condition,
The relative difference of background elasticity modulus can be made to replace with lower aprons in formula (11):
Δ f/f ≈ Δ (lnf),
Δμ/μb≈Δ(lnμb),
Δρ/ρb≈Δ(lnρb), and
In the case where elastic parameter, parameters of fissure, attenuation parameter and elastic impedance parametric continuity change assumed condition, Δ (ln
F) ≈ d (ln f), Δ (ln μb)≈d(lnμb), Δ (ln ρb)≈d(lnρb),ΔδT≈dδT,And
Above approximately replacement is brought into formula (11), following formula is obtained:
Integral is taken to formula (12), and passes through corresponding operation, can get the orientation decaying elastic impedance equation of log-domain,
It is expressed as
Logarithm operation is taken to formula (13), then decaying elastic impedance equation in final orientation is represented by
Preferably, in the step 7, by the decaying elastic impedance difference of different direction to crack " weakness " parameter and crack
It is as follows that induced attenuation parameter carries out the step of inverting:
Step 7-1: difference linearization process is carried out to log-domain different direction decaying elastic impedance, obtains difference linearisation
Expression formula is shown below；
In formula (15), M is the quantity of incidence angle, and N is the quantity of azimuth difference, and L is the quantity of reflecting interface；It is the elastic impedance difference of log-domain different direction, mcQIt is the crack " weakness " and crack induced attenuation parameter to inverting,
Δ X is forward operator matrix relevant to different direction reflection coefficient weight difference operator matrix Δ A, can be indicated respectively
For
Step 7-2: decorrelation and normalized, just asking after obtaining decorrelation are carried out to difference linearized expression
Topic；
The decorrelation and normalized between model parameter are considered in the elastic impedance difference inverting of log-domain orientation, are gone
Nuclear matrix Δ X after correlation becomesDynamic domain model parameter vector mcQBecomeDecorrelation
Direct problem afterwards is represented by
Step 7-3: objective function is converted by the direct problem after decorrelation；
Use Cauchy probability distribution as priori probability density function, uses Gaussian Profile as likelihood function, then posteriority
Probability density function is solved using the joint probability density function of priori probability density function and likelihood function, i.e.,
In formula (17),It is noise variance,It is the variance of quasi-static model parameter vectors.To formula (17)
Maximization posterior probability density function is taken, and combines initial model low-frequency information regularization constraint item, by taking logarithmic transformation, then
Objective function can be ultimately expressed as
In formula (18),It is the regularization coefficient of quasi-static model parameter；
WhereinIndicate the first of model parameter
Value.
Step 7-4: objective function is solved；
Target function type (18) are solved, can be obtained
In formula (18),
QCauchyIndicate Cauchy sparse matrix,
Indicate reflection coefficient；
Step 7-5: solution is iterated to formula (19) using iteration weight weighted least-squares optimization algorithm, obtains dynamic
Domain model parameter vector, formula are as follows:
The invention has the following advantages:
The present invention comprehensively considers the intrinsic decaying of fissure rock background and crack induced attenuation, and combines the crack of plural number characterization
Relation equation between the dry elasticity modulus and saturated flow bulk modulus of rock, and based on weak anisotropy approximation and weak background sheet
Sign decaying is approximate, and within the scope of seismic frequency, it is intrinsic to obtain the dry Modulus of Rocks of background, desiccation crack flexibility, fluid modulus, background
The approximate expression of the saturation fluid apertures gap fissuted medium Complex modes real part of P-wave And S decaying and crack induction attenuation of P-wave characterization
Formula；Based on scattering theory, the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is derived, and by longitudinal wave reflection coefficient
Equation weak elastic parameter is poor, under rimala " weakness " and underdamp assumed condition and elastic parameter, parameters of fissure, attenuation parameter
Change under assumed condition with elastic impedance parametric continuity, is deduced orientation decaying elastic impedance equation；Propose it is sparse about
The orientation decaying elastic impedance Inverse iteration method of beam regularization and low-frequency information constraint regularization, realizes crack elimination back
The pre-stack seismic inversion of scape fluid parameter, fracture water flow parameter and attenuation parameter.
Detailed description of the invention
The accompanying drawings constituting a part of this application is used to provide further understanding of the present application, and the application's shows
Meaning property embodiment and its explanation are not constituted an undue limitation on the present application for explaining the application.
Fig. 1 is the flow chart of attenuation anisotropy elasticity of fluid impedance inversion approach of the present invention；
Specific embodiment
It is noted that following detailed description is all illustrative, it is intended to provide further instruction to the application.Unless another
It indicates, all technical and scientific terms used herein has usual with the application person of an ordinary skill in the technical field
The identical meanings of understanding.
It should be noted that term used herein above is merely to describe specific embodiment, and be not intended to restricted root
According to the illustrative embodiments of the application.As used herein, unless the context clearly indicates otherwise, otherwise singular
Also it is intended to include plural form, additionally, it should be understood that, when in the present specification using term "comprising" and/or " packet
Include " when, indicate existing characteristics, step, operation, device, component and/or their combination.
Present invention will be further explained below with reference to the attached drawings and examples.
A kind of attenuation anisotropy elasticity of fluid impedance inversion approach, as shown in Figure 1, comprising the following steps:
Step 1: the linear slide HTI model in view of the intrinsic decaying of background rock and crack induced attenuation is obtained；
Preferably, in the step 1, isotropism background modulus parameter and crack ginseng in the linear slide HTI model
Number is represented as plural form, and formula is as follows:
In formula (1), WithRespectively indicate the multiple longitudinal wave of Isotropic Decay background rock
Modulus, the first and second Lame constants, and
WithIt indicates
The normal direction of plural form and tangential crack " weakness " parameter can be used along the crack of symmetrical axis direction induction P-wave And S decaying
WithIt is represented simply as
Publicity (2) is brought into the linear slide HTI obtained in formula (1) in view of crack induction crack induced attenuation
Model；
Under the glutinous assumed condition for playing background rock of isotropism, the model parameter of plural formWithIt can indicate
For
In formula (3),WithBackground P-wave And S is respectively indicated against quality factor；
Publicity (3) is brought into and obtains the linear slide HTI mould for considering the intrinsic decaying of background rock in formula (1)
Type；
Publicity (2), (3) are brought into formula (1) simultaneously and are obtained in view of the intrinsic decaying of background rock and crack
The linear slide HTI model of induced attenuation.
Step 2: the relation equation between the dry elasticity modulus of fissure rock and saturated flow bulk modulus is established；
Preferably, the pass in the step 2, between the dry elasticity modulus and saturated flow bulk modulus of the fissure rock
It is that equation is as follows:
In formula (4),It is the saturation fluid fissure rock effective rigidity elastic matrix of plural number characterization,It is plural number
The desiccation crack rock effective rigidity elastic matrix of characterization；It is the anisotropic class Biot coefficient of plural number characterization, can indicates
For
Wherein, KgIndicate the effective bulk modulus of composition rock solid particle；It is the anisotropic class of plural number characterization
Gassmann interstitial space modulus, is represented by
In formula,It is the dry rock volume modulus of anisotropic class of plural number characterization；φ indicates hole
Porosity；κfIndicate pore-fluid effective bulk modulus.
Step 3: it is real that saturation fluid apertures gap fissuted medium Complex modes are obtained using the formula in step 1 and step 2
Portion；
Preferably, in the step 3, the method for obtaining saturation fluid apertures gap fissuted medium Complex modes real part is as follows:
Existing research shows that the reflection coefficient imaginary part of plural number characterization will be far smaller than real part, therefore, in the present invention only
The derivation of reflection coefficient real part is focused on, and ignores the imaginary part of reflection coefficient；
Wherein, the real part of the desiccation crack rock weak anisotropy stiffness matrix of plural number characterization is represented by
In formula (7),
WithThe intrinsic P-wave And S decaying of characterization background and crack induction are vertical respectively
The overall attenuation factor of wave attenuation product；
Indicate that background is done under rocky conditionWithRatio；
Wherein the parameter in formula of the present invention with subscript dry indicates in the case of dry rock (not having fluid in hole)
Parameter, such as:Indicate the multiple longitudinal wave modulus of Isotropic Decay background rock, thenThen indicate that Isotropic Decay is dry
The multiple longitudinal wave modulus of background rock；
Formula (5), (6), (7) are substituted into formula (4), and close based on weak anisotropy approximation and the intrinsic decaying of weak background
Seemingly, then within the scope of seismic frequency, the dry Modulus of Rocks of background, desiccation crack flexibility, fluid modulus, the intrinsic P-wave And S of background are obtained
Decay and the approximate expression of the saturation fluid apertures gap fissuted medium Complex modes real part of crack induction attenuation of P-wave characterization is
In formula (8),
β0=1-Kdry/KgIt is Biot coefficient；
KdryIndicate the bulk modulus of dry rock.
Step 4: the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is obtained；
Preferably, in the step 4, the method for obtaining PP wave reflection coefficient equation is as follows:
Based on scattering theory, the longitudinal wave reflection coefficient of HTI medium be may be expressed as:
Wherein, θ indicates incidence angle, and ρ indicates the density item of homogeneous isotropism background media,Indicate saturated rock
Disturbance stiffness parameters, ξmnIt is related with slowness vector sum polarization vector；
Therefore, convolution (8) and (9) can derive the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium, specifically
Ground, in formula (9)Carry out derivation and then and ξmnEvent was ignored in multiplication can be obtained by formula (10), expression formula
For
In formula,
In formula (10),Indicate azimuth, Δ indicates the difference of each attribute of two layers of medium, f=fluid/hole
Gap is a fluid factor；Crack " weakness " parameter can be estimated according to log data and petrophysical model；
Indicate that background is saturated in the case of fluidWithRatio；
Wherein the parameter in the present invention with subscript sat indicates the parameter in the case of saturation fluid.
Step 5: longitudinal wave reflection coefficient equation is obtained；
Preferably, in the step 5, according to horizontal interface longitudinal wave reflection coefficient and azimuthal anisotropy elastic impedance it
Between relationship, longitudinal wave reflection coefficient equation may be expressed as:
In formula (11),Indicate that orientation decays elastic impedance, Δ EI (θ, φ) indicates upper one layer and next
The difference of layer elastic impedance；
Wherein longitudinal wave is exactly p wave, and it is all p wave that pp, which indicates incident and reflection,.
Step 6: decaying elastic impedance equation in orientation is obtained by the longitudinal wave reflection coefficient equation of step 5；
Preferably, in the step 6, orientation decaying elastic impedance side is obtained by the longitudinal wave reflection coefficient equation of step 5
The method of journey is as follows:
Weak elastic parameter it is poor (| Δ f/f | < < 1, | Δ μ/μb| < < 1 and | Δ ρ/ρb| < < 1), rimala " weakness " (And δT< < 1) and underdamp (And) under assumed condition,
The relative difference of background elasticity modulus can be made to replace with lower aprons in formula (11):
Δ f/f ≈ Δ (lnf),
Δμ/μb≈Δ(lnμb),
Δρ/ρb≈Δ(lnρb), and
In the case where elastic parameter, parameters of fissure, attenuation parameter and elastic impedance parametric continuity change assumed condition, Δ (ln
F) ≈ d (ln f), Δ (ln μb)≈d(lnμb), Δ (ln ρb)≈d(lnρb),ΔδT≈dδT,And
Above approximately replacement is brought into formula (11), following formula is obtained:
Integral is taken to formula (12), and passes through corresponding operation, can get the orientation decaying elastic impedance equation of log-domain,
It is expressed as
Logarithm operation is taken to formula (13), then decaying elastic impedance equation in final orientation is represented by
Step 7: by different direction decaying elastic impedance difference to crack " weakness " parameter and crack induced attenuation parameter into
Row inverting.
Preferably, in the step 7, by the decaying elastic impedance difference of different direction to crack " weakness " parameter and crack
It is as follows that induced attenuation parameter carries out the step of inverting:
Step 7-1: difference linearization process is carried out to log-domain different direction decaying elastic impedance, obtains difference linearisation
Expression formula is shown below；
In formula (15), M is the quantity of incidence angle, and N is the quantity of azimuth difference, and L is the quantity of reflecting interface；It is the elastic impedance difference of log-domain different direction, mcQIt is the crack " weakness " and crack induced attenuation parameter to inverting,
Δ X is forward operator matrix relevant to different direction reflection coefficient weight difference operator matrix Δ A, can be indicated respectively
For
Step 7-2: decorrelation and normalized, just asking after obtaining decorrelation are carried out to difference linearized expression
Topic；
The decorrelation and normalized between model parameter are considered in the elastic impedance difference inverting of log-domain orientation, are gone
Nuclear matrix Δ X after correlation becomesDynamic domain model parameter vector mcQBecomeDecorrelation
Direct problem afterwards is represented by
Step 7-3: objective function is converted by the direct problem after decorrelation；
Use Cauchy probability distribution as priori probability density function, uses Gaussian Profile as likelihood function, then posteriority
Probability density function is solved using the joint probability density function of priori probability density function and likelihood function, i.e.,
Herein apply Bayes principle: Posterior distrbutionp=prior probability × likelihood function, to the direct problem after decorrelation into
Row processing；
In formula (17),It is noise variance,It is the variance of quasi-static model parameter vectors.To formula (17)
Maximization posterior probability density function is taken, and combines initial model low-frequency information regularization constraint item, by taking logarithmic transformation, then
Objective function can be ultimately expressed as
In formula (18),It is the regularization coefficient of quasi-static model parameter；
WhereinIndicate the first of model parameter
Value.
Step 7-4: objective function is solved；
Target function type (18) are solved, can be obtained
In formula (18),
QCauchyIndicate Cauchy sparse matrix,
Indicate reflection coefficient；
Step 7-5: solution is iterated to formula (19) using iteration weight weighted least-squares optimization algorithm, obtains dynamic
Domain model parameter vector, formula are as follows:
Orientation decaying elastic impedance inverting side proposed by the invention is verified using the orientation prestack seismic gather of synthesis
The reliability of method: when synthesizing trace gather without noise or containing appropriate noise, using the model parameter of inversion method inverting of the present invention
Inversion result and real data coincide it is preferable.The orientation angles trace gather and original trace gather synthesized by model parameter inversion result
Comparison, comparison further demonstrate the feasibility of orientation decaying elastic impedance inversion method it is found that the two difference is smaller.
The present invention comprehensively considers the intrinsic decaying of fissure rock background and crack induced attenuation, and combines the crack of plural number characterization
Relation equation between the dry elasticity modulus and saturated flow bulk modulus of rock, and based on weak anisotropy approximation and weak background sheet
Sign decaying is approximate, and within the scope of seismic frequency, it is intrinsic to obtain the dry Modulus of Rocks of background, desiccation crack flexibility, fluid modulus, background
The approximate expression of the saturation fluid apertures gap fissuted medium Complex modes real part of P-wave And S decaying and crack induction attenuation of P-wave characterization
Formula；Based on scattering theory, the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is derived, and by longitudinal wave reflection coefficient
Equation weak elastic parameter is poor, under rimala " weakness " and underdamp assumed condition and elastic parameter, parameters of fissure, attenuation parameter
Change under assumed condition with elastic impedance parametric continuity, is deduced orientation decaying elastic impedance equation；Propose it is sparse about
The orientation decaying elastic impedance Inverse iteration method of beam regularization and low-frequency information constraint regularization, realizes crack elimination back
The pre-stack seismic inversion of scape fluid parameter, fracture water flow parameter and attenuation parameter.
Above-mentioned, although the foregoing specific embodiments of the present invention is described with reference to the accompanying drawings, not to limit of the invention
System, those skilled in the art should understand that, based on the technical solutions of the present invention, those skilled in the art do not need to pay
The various modifications or changes that creative work can be made out are still within protection scope of the present invention.
Claims (8)
1. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach, characterized in that the following steps are included:
Step 1: the linear slide HTI model in view of the intrinsic decaying of background rock and crack induced attenuation is obtained；
Step 2: the relation equation between the dry elasticity modulus of fissure rock and saturated flow bulk modulus is established；
Step 3: saturation fluid apertures gap fissuted medium Complex modes real part is obtained using the formula in step 1 and step 2；
Step 4: the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium is obtained；
Step 5: longitudinal wave reflection coefficient equation is obtained；
Step 6: decaying elastic impedance equation in orientation is obtained by the longitudinal wave reflection coefficient equation of step 5；
Step 7: crack " weakness " parameter and crack induced attenuation parameter are carried out by the decaying elastic impedance difference of different direction anti-
It drills.
2. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as shown in claim 1, characterized in that the step
In one, isotropism background modulus parameter and parameters of fissure are represented as plural form, formula in the linear slide HTI model
It is as follows:
In formula (1), WithRespectively indicate Isotropic Decay background rock multiple longitudinal wave modulus,
First and second Lame constants, and
WithIndicate plural shape
The normal direction of formula and tangential crack " weakness " parameter can be used along the crack of symmetrical axis direction induction P-wave And S decayingWithLetter
Singly it is expressed as
Under the glutinous assumed condition for playing background rock of isotropism, the model parameter of plural formWithIt is represented by
In formula (3),WithBackground P-wave And S is respectively indicated against quality factor；
Publicity (2), (3) are brought into formula (1) simultaneously and are obtained in view of the intrinsic decaying of background rock and crack induction
The linear slide HTI model of decaying.
3. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as shown in claim 2, characterized in that the step
In two, the relation equation between the dry elasticity modulus and saturated flow bulk modulus of the fissure rock is as follows:
In formula (4),It is the saturation fluid fissure rock effective rigidity elastic matrix of plural number characterization,It is plural characterization
Desiccation crack rock effective rigidity elastic matrix；It is the anisotropic class Biot coefficient of plural number characterization, is represented by
Wherein, KgIndicate the effective bulk modulus of composition rock solid particle；It is the anisotropic class of plural number characterization
Gassmann interstitial space modulus, is represented by
In formula,It is the dry rock volume modulus of anisotropic class of plural number characterization；φ indicates hole
Degree；κfIndicate pore-fluid effective bulk modulus.
4. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as stated in claim 3, characterized in that the step
In three, the method for obtaining saturation fluid apertures gap fissuted medium Complex modes real part is as follows:
Wherein, the real part of the desiccation crack rock weak anisotropy stiffness matrix of plural number characterization is represented by
In formula (7),
WithThe intrinsic P-wave And S decaying of characterization background declines with crack induction longitudinal wave respectively
Subtract the overall attenuation factor of product；
Formula (5), (6), (7) are substituted into formula (4), and approximate based on weak anisotropy approximation and the intrinsic decaying of weak background, then
Within the scope of seismic frequency, obtain the decaying of the dry Modulus of Rocks of background, desiccation crack flexibility, fluid modulus, background intrinsic P-wave And S and
The approximate expression of saturation fluid apertures gap fissuted medium Complex modes real part of crack induction attenuation of P-wave characterization is
In formula (8),
β0=1-Kdry/KgIt is Biot coefficient；
KdryIndicate the bulk modulus of dry rock.
5. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as shown in claim 4, characterized in that the step
In four, the method for obtaining PP wave reflection coefficient equation is as follows:
Based on scattering theory, the longitudinal wave reflection coefficient of HTI medium be may be expressed as:
Wherein, θ indicates incidence angle, and ρ indicates the density item of homogeneous isotropism background media,Indicate the disturbance of saturated rock
Stiffness parameters, ξmnIt is related with slowness vector sum polarization vector；
Therefore, convolution (8) and (9) can derive the PP wave reflection coefficient equation of dynamic equivalent HTI attenuation medium, expression formula
For
In formula,
In formula (10),Indicate azimuth, Δ indicates the difference of each attribute of two layers of medium, and f=fluid/hole, is one
A fluid factor；Crack " weakness " parameter can be estimated according to log data and petrophysical model.
6. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as stated in claim 5, characterized in that the step
In five, according to the relationship between horizontal interface longitudinal wave reflection coefficient and azimuthal anisotropy elastic impedance, longitudinal wave reflection coefficient side
Journey may be expressed as:
In formula (11),Indicate orientation decaying elastic impedance, Δ EI (θ, φ) indicates upper one layer and next layer of bullet
The difference of property impedance.
7. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as shown in claim 6, characterized in that the step
In six, the method for obtaining orientation decaying elastic impedance equation by the longitudinal wave reflection coefficient equation of step 5 is as follows:
Weak elastic parameter it is poor (| Δ f/f | < < 1, | Δ μ/μb| < < 1 and | Δ ρ/ρb| < < 1), rimala " weakness " (And δT< < 1) and underdamp (And) under assumed condition,
The relative difference of background elasticity modulus can be made to replace with lower aprons in formula (11):
Δ f/f ≈ Δ (ln f),
Δμ/μb≈Δ(lnμb),
Δρ/ρb≈Δ(lnρb), and
In the case where elastic parameter, parameters of fissure, attenuation parameter and elastic impedance parametric continuity change assumed condition, Δ (ln f) ≈
D (ln f), Δ (ln μb)≈d(lnμb), Δ (ln ρb)≈d(lnρb),ΔδT≈dδT,And
Above approximately replacement is brought into formula (11), following formula is obtained:
Integral is taken to formula (12), and passes through corresponding operation, the orientation decaying elastic impedance equation of log-domain is can get, indicates
For
Logarithm operation is taken to formula (13), then decaying elastic impedance equation in final orientation is represented by
8. a kind of attenuation anisotropy elasticity of fluid impedance inversion approach as shown in claim 7, characterized in that the step
In seven, the step of inverting is carried out to crack " weakness " parameter and crack induced attenuation parameter by the decaying elastic impedance difference of different direction
It is rapid as follows:
Step 7-1: difference linearization process is carried out to log-domain different direction decaying elastic impedance, obtains difference linearisation expression
Formula is shown below；
In formula (15), M is the quantity of incidence angle, and N is the quantity of azimuth difference, and L is the quantity of reflecting interface；It is
The elastic impedance difference of log-domain different direction, mcQThe crack " weakness " and crack induced attenuation parameter to inverting, Δ X be with not
With orientation reflection coefficient weight difference operator matrix Δ A, relevant forward operator matrix, can be expressed as respectively
[δT]L×1=[δT(t1),...,δT(tL)]T,
Step 7-2: decorrelation and normalized are carried out to difference linearized expression, the direct problem after obtaining decorrelation；
The decorrelation and normalized between model parameter, decorrelation are considered in the elastic impedance difference inverting of log-domain orientation
Nuclear matrix Δ X afterwards becomesDynamic domain model parameter vector mcQBecomeAfter decorrelation
Direct problem is represented by
Step 7-3: objective function is converted by the direct problem after decorrelation；
Use Cauchy probability distribution as priori probability density function, uses Gaussian Profile as likelihood function, then posterior probability
Density function is solved using the joint probability density function of priori probability density function and likelihood function, i.e.,
In formula (17),It is noise variance,It is the variance of quasi-static model parameter vectors.Formula (17) is taken most
Bigization posterior probability density function, and combine initial model low-frequency information regularization constraint item, by taking logarithmic transformation, then target
Function can be ultimately expressed as
In formula (18),It is the regularization coefficient of quasi-static model parameter；
WhereinIndicate the initial value of model parameter；
Step 7-4: objective function is solved；
Target function type (18) are solved, can be obtained
In formula (18),
QCauchyIndicate Cauchy sparse matrix,
Indicate reflection coefficient；
Step 7-5: solution is iterated to formula (19) using iteration weight weighted least-squares optimization algorithm, obtains dynamic domain mould
Shape parameter vector, formula are as follows:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910407041.5A CN110133718B (en) | 2019-05-16 | 2019-05-16 | Attenuation anisotropic fluid elastic impedance inversion method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910407041.5A CN110133718B (en) | 2019-05-16 | 2019-05-16 | Attenuation anisotropic fluid elastic impedance inversion method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110133718A true CN110133718A (en) | 2019-08-16 |
| CN110133718B CN110133718B (en) | 2020-12-01 |
Family
ID=67574530
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910407041.5A Expired - Fee Related CN110133718B (en) | 2019-05-16 | 2019-05-16 | Attenuation anisotropic fluid elastic impedance inversion method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110133718B (en) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110515126A (en) * | 2019-09-12 | 2019-11-29 | 中国石油大学(华东) | A kind of velocity of sound calculation method of the transverse isotropy of crack containing random distribution rock |
| CN110687601A (en) * | 2019-11-01 | 2020-01-14 | 中南大学 | Inversion method for fluid factor and fracture parameter of orthotropic medium |
| CN111077568A (en) * | 2019-12-20 | 2020-04-28 | 中国石油大学(北京) | Method and equipment for detecting oil and gas reservoir by fluid factor of tight oil and gas reservoir |
| CN111077573A (en) * | 2019-12-30 | 2020-04-28 | 中国石油大学(北京) | Method, device and system for determining stratum elastic parameters |
| CN111551990A (en) * | 2020-04-07 | 2020-08-18 | 西安科技大学 | Method and system for calculating seismic wave reflection coefficient of HTI (HTI) coal seam |
| CN112230276A (en) * | 2020-09-30 | 2021-01-15 | 甘肃省地震局（中国地震局兰州地震研究所） | Fracture type tight reservoir fluid identification method, system, identification instrument, medium and application |
| CN112630829A (en) * | 2019-10-08 | 2021-04-09 | 中国石油化工股份有限公司 | Method and system for analyzing tight sandstone elastic wave attenuation property |
| CN112649871A (en) * | 2020-12-18 | 2021-04-13 | 中国矿业大学(北京) | Longitudinal wave reflection coefficient determining method and device, electronic equipment and storage medium |
| CN113219536A (en) * | 2021-06-25 | 2021-08-06 | 成都理工大学 | Pre-stack seismic inversion method of longitudinal and transverse wave attenuation parameters depending on frequency |
| CN114063163A (en) * | 2021-11-11 | 2022-02-18 | 中南大学 | Fracture type reservoir monoclinic equivalent medium seismic characterization and inversion method and system |
| CN116522810A (en) * | 2023-04-07 | 2023-08-01 | 中国地质调查局油气资源调查中心 | Reservoir frequency dispersion attenuation analysis method based on Norris-KG equivalent medium modeling |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070280048A1 (en) * | 2006-06-06 | 2007-12-06 | Baker Hughes Incorporated | P-wave anisotropy evaluation by measuring acoustic impedance of the rock by beam-steering from within the borehole at different angles |
| CN102147478A (en) * | 2010-12-29 | 2011-08-10 | 中国海洋大学 | Pre-stack low frequency signal recognition method of complex oil pool |
| CN102879800A (en) * | 2011-07-15 | 2013-01-16 | 中国石油天然气集团公司 | Method for detecting shear-wave splitting fracture |
| CN104820239A (en) * | 2015-05-13 | 2015-08-05 | 中国石油大学(华东) | Azimuth pre-stack seismic attribution decoupling extraction method |
| CN105005074A (en) * | 2015-06-23 | 2015-10-28 | 成都理工大学 | Method for identifying gas reservoir by using frequency-variable seismic reflection coefficient |
-
2019
- 2019-05-16 CN CN201910407041.5A patent/CN110133718B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20070280048A1 (en) * | 2006-06-06 | 2007-12-06 | Baker Hughes Incorporated | P-wave anisotropy evaluation by measuring acoustic impedance of the rock by beam-steering from within the borehole at different angles |
| CN102147478A (en) * | 2010-12-29 | 2011-08-10 | 中国海洋大学 | Pre-stack low frequency signal recognition method of complex oil pool |
| CN102879800A (en) * | 2011-07-15 | 2013-01-16 | 中国石油天然气集团公司 | Method for detecting shear-wave splitting fracture |
| CN104820239A (en) * | 2015-05-13 | 2015-08-05 | 中国石油大学(华东) | Azimuth pre-stack seismic attribution decoupling extraction method |
| CN105005074A (en) * | 2015-06-23 | 2015-10-28 | 成都理工大学 | Method for identifying gas reservoir by using frequency-variable seismic reflection coefficient |
Non-Patent Citations (1)
| Title |
|---|
| 潘新朋，等: "裂缝-孔隙型含气储层流体与裂缝参数贝叶斯地震反演方法", 《中国科学：地球科学》 * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110515126A (en) * | 2019-09-12 | 2019-11-29 | 中国石油大学(华东) | A kind of velocity of sound calculation method of the transverse isotropy of crack containing random distribution rock |
| CN112630829A (en) * | 2019-10-08 | 2021-04-09 | 中国石油化工股份有限公司 | Method and system for analyzing tight sandstone elastic wave attenuation property |
| CN110687601A (en) * | 2019-11-01 | 2020-01-14 | 中南大学 | Inversion method for fluid factor and fracture parameter of orthotropic medium |
| CN111077568A (en) * | 2019-12-20 | 2020-04-28 | 中国石油大学(北京) | Method and equipment for detecting oil and gas reservoir by fluid factor of tight oil and gas reservoir |
| CN111077568B (en) * | 2019-12-20 | 2021-04-23 | 中国石油大学(北京) | Method and equipment for detecting oil and gas reservoir by fluid factor of tight oil and gas reservoir |
| CN111077573A (en) * | 2019-12-30 | 2020-04-28 | 中国石油大学(北京) | Method, device and system for determining stratum elastic parameters |
| CN111551990A (en) * | 2020-04-07 | 2020-08-18 | 西安科技大学 | Method and system for calculating seismic wave reflection coefficient of HTI (HTI) coal seam |
| CN111551990B (en) * | 2020-04-07 | 2023-05-23 | 西安科技大学 | Acquisition system for HTI type coal seam seismic wave reflection coefficient |
| CN112230276A (en) * | 2020-09-30 | 2021-01-15 | 甘肃省地震局（中国地震局兰州地震研究所） | Fracture type tight reservoir fluid identification method, system, identification instrument, medium and application |
| CN112230276B (en) * | 2020-09-30 | 2023-07-07 | 兰州石化职业技术学院 | Crack type tight reservoir fluid identification method, system, identification instrument, medium and application |
| CN112649871A (en) * | 2020-12-18 | 2021-04-13 | 中国矿业大学(北京) | Longitudinal wave reflection coefficient determining method and device, electronic equipment and storage medium |
| CN112649871B (en) * | 2020-12-18 | 2021-08-24 | 中国矿业大学(北京) | Longitudinal wave reflection coefficient determining method and device, electronic equipment and storage medium |
| CN113219536A (en) * | 2021-06-25 | 2021-08-06 | 成都理工大学 | Pre-stack seismic inversion method of longitudinal and transverse wave attenuation parameters depending on frequency |
| CN113219536B (en) * | 2021-06-25 | 2022-03-01 | 成都理工大学 | Pre-stack seismic inversion method of longitudinal and transverse wave attenuation parameters depending on frequency |
| CN114063163A (en) * | 2021-11-11 | 2022-02-18 | 中南大学 | Fracture type reservoir monoclinic equivalent medium seismic characterization and inversion method and system |
| CN114063163B (en) * | 2021-11-11 | 2022-11-11 | 中南大学 | Fracture type reservoir monoclinic equivalent medium seismic characterization and inversion method and system |
| CN116522810A (en) * | 2023-04-07 | 2023-08-01 | 中国地质调查局油气资源调查中心 | Reservoir frequency dispersion attenuation analysis method based on Norris-KG equivalent medium modeling |
| CN116522810B (en) * | 2023-04-07 | 2024-01-30 | 中国地质调查局油气资源调查中心 | Reservoir frequency dispersion attenuation analysis method based on Norris-KG equivalent medium modeling |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110133718B (en) | 2020-12-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110133718A (en) | A kind of attenuation anisotropy elasticity of fluid impedance inversion approach | |
| US10620331B2 (en) | Reverse time migration in anisotropic media with stable attenuation compensation | |
| Zhang et al. | Amplitude-preserving reverse time migration: From reflectivity to velocity and impedance inversion | |
| Regone et al. | Geologic model building in SEAM Phase II—Land seismic challenges | |
| Wang et al. | Accurate porosity prediction for tight sandstone reservoir: A case study from North China | |
| Beckwith et al. | Estimating frequency-dependent attenuation quality factor values from prestack surface seismic data | |
| Behura et al. | Estimation of interval anisotropic attenuation from reflection data | |
| Thiel et al. | Comparison of acoustic and elastic full‐waveform inversion of 2D towed‐streamer data in the presence of salt | |
| Rabben et al. | AVA inversion of the top Utsira Sand reflection at the Sleipner field | |
| Zhang | Simultaneous inversion for S-wave velocity and density from the SV-SV wave | |
| Wang et al. | Analysis and estimation of an inclusion-based effective fluid modulus for tight gas-bearing sandstone reservoirs | |
| Jiang et al. | Rock-physics and seismic-inversion based reservoir characterization of the Haynesville Shale | |
| Zhu et al. | Recent applications of turning-ray tomography | |
| US11604299B2 (en) | Mixed-phase source wavelet estimation from recorded seismic data | |
| CN105929453B (en) | State estimation method for infinite distribution time lag neural network system with channel fading | |
| Zhang et al. | Wave equation tomographic velocity inversion method based on the Born/Rytov approximation | |
| Schwenk | Constrained parameterization of the multichannel analysis of surface waves approach with application at Yuma Proving Ground, Arizona | |
| Wang et al. | Challenges in Brazil pre-salt carbonate reservoir characterization: a realistic 3D synthetic study | |
| Zhang et al. | One-way wave propagation in the ray-centred coordinate system for vertical transversely isotropic media | |
| Shekar et al. | Kirchhoff modeling for attenuative anisotropic media using Gaussian beams | |
| Chen et al. | Observation of a mesoscale warm eddy impacts acoustic propagation in the slope of the South China Sea | |
| Pei et al. | Research and application of 5D seismic prediction technology | |
| Yang et al. | Micro-seismic monitoring using sparse planar array and a weak signal enhancement method | |
| Glushchenko et al. | A broadband full azimuth land seismic case study Onshore Abu Dhabi | |
| Wang et al. | DIRECT PREDICTION METHOD OF FRACTURING ABILITY IN SHALE FORMATIONS BASED ON PRE-STACK SEISMIC 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 | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20191101 Address after: 410083 Hunan province Changsha Lushan Road No. 932 Applicant after: Central South University Applicant after: China University of Petroleum (East China) Address before: 266580 Qingdao economic and Technological Development Zone, Changjiang Road, No. 66, Shandong Applicant before: China University of Petroleum (East China) |
|
| TA01 | Transfer of patent application right | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201201 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |