CN106959467A - Seimic wave velocity inversion method and device - Google Patents
Seimic wave velocity inversion method and device Download PDFInfo
- Publication number
- CN106959467A CN106959467A CN201710168052.3A CN201710168052A CN106959467A CN 106959467 A CN106959467 A CN 106959467A CN 201710168052 A CN201710168052 A CN 201710168052A CN 106959467 A CN106959467 A CN 106959467A
- Authority
- CN
- China
- Prior art keywords
- data
- wave
- seismic
- seismic data
- common imaging
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 44
- 238000003384 imaging method Methods 0.000 claims abstract description 84
- 238000013508 migration Methods 0.000 claims abstract description 43
- 230000005012 migration Effects 0.000 claims abstract description 43
- 238000012545 processing Methods 0.000 claims abstract description 40
- 238000012937 correction Methods 0.000 claims abstract description 17
- 238000000926 separation method Methods 0.000 claims description 9
- 230000008859 change Effects 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000003068 static effect Effects 0.000 claims description 6
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000013075 data extraction Methods 0.000 claims description 3
- 238000004321 preservation Methods 0.000 claims description 3
- 230000007423 decrease Effects 0.000 claims description 2
- 238000009826 distribution Methods 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 description 5
- 230000006870 function Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 230000002411 adverse Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000004587 chromatography analysis Methods 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 239000010410 layer Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention provides a kind of seimic wave velocity inversion method and device；Wherein, this method includes：Obtain the seismic data of earth's surface in default geologic province；Wherein, the seismic data of the earth's surface includes shear wave data and P wave data；Using initial velocity model, pre-stack depth migration processing is carried out to seismic data, the common imaging gather of seismic wave is generated；The remaining depth data of lineups is extracted from common imaging gather；Wherein, the remaining depth data includes remaining depth correction amount and the reflecting interface position of seismic wave；According to remaining depth data, the seismic data of underground medium in inverting geologic province；Wherein, the seismic data of the underground medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters.The seismic data accuracy that the mode for the seimic wave velocity inverting that the present invention is provided is obtained is higher, and then improves the resolution ratio of seismic imaging.
Description
Technical field
The present invention relates to seismic exploration technical field, more particularly, to a kind of seimic wave velocity inversion method and device.
Background technology
With the development and the raising of computer software and hardware level of oilfield prospecting developing technology, pre-stack depth migration
(PSDM) it is increasingly becoming the key technology of the regional high-resolution seismic exploration processing of complicated structure.In PSDM, Wave velocity calculation
Inaccurate follow-up imaging and explanation can be caused more significantly to adversely affect.Therefore, speed point of the research suitable for PSDM
Analysis and inversion method have important theory significance and practical value.
In the prior art, the PSDM guiding chromatographies based on ray theory develop to migration velocity analysis.These methods include
Carry out renewal speed by the way that interpretation horizon, as constraint, is solved into interlayer velocity gradient；And utilize the data-driven based on grid
Algorithm realizes velocity analysis.However, the optimal speed model for the seismic wave that aforesaid way is obtained has multi-solution, and then cause speed
Degrees of data is inaccurate, have impact on seismic tomography effect.
For above-mentioned seimic wave velocity inverting mode accuracy it is poor the problem of, not yet propose effective solution.
The content of the invention
In view of this, it is an object of the invention to provide a kind of seimic wave velocity inversion method and device, to improve earthquake
The seismic data accuracy that the mode of wave velocity inverting is obtained, and then improve the resolution ratio of seismic imaging.
In a first aspect, the embodiments of the invention provide a kind of seimic wave velocity inversion method, including：Obtain default geological province
The seismic data of earth's surface in domain；Wherein, the seismic data of the earth's surface includes shear wave data and P wave data；Utilize initial speed
Model is spent, pre-stack depth migration processing is carried out to seismic data, the common imaging gather of seismic wave is generated；Wherein, this is initial
Rate pattern is set up according to historical earthquake wave number according to the corresponding relation of historical earthquake wave velocity；Extracted from common imaging gather
The remaining depth data of lineups；Wherein, the remaining depth data includes the remaining depth correction amount and reflecting interface of seismic wave
Position；According to remaining depth data, the seismic data of underground medium in inverting geologic province；Wherein, the ground of the underground medium
Seismic wave data include shear wave velocity, velocity of longitudinal wave and anisotropic parameters.
With reference in a first aspect, the embodiments of the invention provide the possible embodiment of the first of first aspect, wherein, on
Stating the seismic data of earth's surface in the default geologic province of acquisition includes：Subsurface geology letter is carried in the default geologic province of collection
The earthquake shot gather data of breath；Wherein, the subsurface geology information includes geological structure information and geologic lithology change information；By dissipating
Degree and curl operator, wave field separation processing, generation shear wave data and P wave data are carried out to earthquake shot gather data.
With reference in a first aspect, the embodiments of the invention provide the possible embodiment of second of first aspect, wherein, on
State and utilize initial velocity model, pre-stack depth migration processing is carried out to seismic data, generate the common imaging gather of seismic wave
Including：Seismic data is pre-processed, seismic data of the generation available for skew；Wherein, the pretreatment is included over the ground
Seismic wave data carry out noise remove, static correction and to seismic data and the historical earthquake wave number pre-saved according to progress one
One correspondence；Using initial velocity model, pre-stack depth migration processing is carried out to the seismic data available for skew, earthquake is generated
The common imaging gather of ripple；Wherein, the initial velocity model includes constant velocity model；The abscissa of common imaging gather is skew
Away from；The ordinate of common imaging gather is depth；The numerical value of common imaging gather is seismic amplitude；Generate seismic wave altogether into
After picture point road collection, the above method also includes：Summed along the X direction of common imaging gather, generate the folded of seismic wave
Preceding depth migration data.
With reference to second of possible embodiment of first aspect, the embodiments of the invention provide the third of first aspect
Possible embodiment, wherein, the above-mentioned remaining depth data that lineups are extracted from common imaging gather includes：From common imaging
Corresponding first depth value of each offset distance is extracted in point road collection；Wherein, the offset distance be geophone and earthquake shot point it
Between distance；First depth value is handled through pre-stack depth migration by the corresponding geological data of offset distance and obtained；From common imaging
Corresponding second depth value of zero-offset is extracted in point road collection；The difference of the first depth value and the second depth value is calculated, is total to
The remaining depth correction amount of image gather；From the reflecting interface of pre-stack depth migration extracting data seismic wave, reflected
Interface location；Wherein, the reflecting interface position includes the locus coordinate of reflecting interface.
With reference to the third possible embodiment of first aspect, the embodiments of the invention provide the 4th of first aspect kind
Possible embodiment, wherein, above-mentioned according to remaining depth data, the seismic data bag of underground medium in inverting geologic province
Include：Grid discretization processing, formation speed mesh point are carried out to initial velocity model and reflecting interface position；Pass through 2 rays
The mode of tracking, calculates shot point and geophone station position to the ray of pip respectively；Wherein, the pip is located at reflecting interface position
Put；The information of ray at least includes ray traveltime and propagation path；By the point distribution on ray to putting closest speed
Spend on mesh point；The corresponding parameter of speed mesh point is set to the seismic data of underground medium in geologic province；Pass through ground
Seismic velocity inverse model, solves the seismic data of underground medium in geologic province；Wherein, the seimic wave velocity inverse model
Set up according to remaining depth data.
With reference to the 4th kind of possible embodiment of first aspect, the embodiments of the invention provide the 5th of first aspect kind
Possible embodiment, wherein, above-mentioned seimic wave velocity inverse model includes：
Wherein, J (Δ λ) is object function；A is M1×M2The matrix of row N row；The element of matrix A is excursions depths to model
The derivative of parameterz0For remaining depth；xjFor the abscissa positions of inverted parameters；hkFor offset distance
Half；λiFor inverted parameters；M1For the quantity of the offset distance of common imaging gather；M2For the common imaging for renewal speed model
The quantity of point road collection；P represents compressional wave；S represents shear wave；Vectorial b includes M1×M2Individual element,Z is
Excursions depths；For the average value of different offset distance excursions depths；Vectorial Δ λ is model modification amount；Δλ0For initial model to
Amount；β is regularization parameter；| | | | norm is sought in representative.
With reference to the 4th kind of possible embodiment of first aspect, the embodiments of the invention provide the 6th of first aspect kind
Possible embodiment, wherein, above by seimic wave velocity inverse model, solve the seismic wave of underground medium in geologic province
Data include：Initial temperature T is set0, initial solution Δ λ0, outer circulation iterations outn, interior circulation Markov chain length inn, temperature
Spend attenuation coefficient dec, the excursion a and minimum temperature t of model parametermin, i=0；Carry out following outer circulation operations：Adopt
Δ λ is calculated with Memoryless Quasi-Neton Methodi, preserve Δ λi；J (Δ λ are calculated by seimic wave velocity inverse modelι), preserve J (Δs
λι)；K=0, t are setk=T0With scale factor b=0；Carry out following interior circulation operations：Judge tk<tminWhether set up；If
It is to perform following update to operate；If not, setting Δ λk=Δ λi；Calculate Δ λk+1=Δ λk+ξka；Wherein, ξk=tk*sign
(r-0.5)*[(1-1/tk)|2r-1|- 1], r is to obey the random number uniformly divided between [0,1]；The difference of calculating target function value
Δ E=J (Δ λk+1)-J(Δλk)；Judge whether Δ E≤0 sets up；If Δ E≤0, Δ λ is setk+1=Δ λk；If Δ E>0,
JudgeWhether set up；IfΔ λ is setk+1=Δ λk；IfJudging k=inn is
It is no to set up；If not, continuing executing with above-mentioned interior circulation operation；If it is, setting tk+1=tk*(dec)b, b=b+1, k=k+1；
Continue executing with above-mentioned interior circulation operation；Following update is carried out to operate：I=i+1 is set；Judge i<Whether outn sets up；If it is,
Continue executing with above-mentioned outer circulation operation；If not, from J (the Δ λ of preservationι) in value, filter out J (the Δ λ of minimumι) value；By minimum
J (Δ λι) the corresponding velocity amplitude of value is defined as seismic data.
Second aspect, the embodiments of the invention provide a kind of seimic wave velocity inverting device, including：Seismic data is obtained
Module, the seismic data for obtaining earth's surface in default geologic province；Wherein, the seismic data of the earth's surface includes shear wave number
According to and P wave data；Pre-stack depth migration processing module, for utilizing initial velocity model, carries out prestack deep to seismic data
Migration processing is spent, the common imaging gather of seismic wave is generated；Wherein, the initial velocity model according to historical earthquake wave number according to going through
The corresponding relation of history seimic wave velocity is set up；Remaining depth data extraction module, for extracting same phase from common imaging gather
The remaining depth data of axle；Wherein, the remaining depth data includes remaining depth correction amount and the reflecting interface position of seismic wave；
Seismic data inverting module, for according to remaining depth data, the seismic data of underground medium in inverting geologic province；Its
In, the seismic data of the underground medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters.
With reference to second aspect, the embodiments of the invention provide the possible embodiment of the first of second aspect, wherein, on
Stating seismic data acquisition module includes：Collecting unit, subsurface geology information is carried for gathering in default geologic province
Earthquake shot gather data；Wherein, the subsurface geology information includes geological structure information and geologic lithology change information；Wave field separation
Processing unit, for by divergence and curl operator, carrying out wave field separation processing to earthquake shot gather data, generation shear wave data and
P wave data.
With reference to second aspect, the embodiments of the invention provide the possible embodiment of second of second aspect, wherein, on
Stating pre-stack depth migration processing module includes：Pretreatment unit, for being pre-processed to seismic data, generation can be used for inclined
The seismic data of shifting；Wherein, the pretreatment includes carrying out seismic data noise remove, static correction and to seismic wave
Data and the historical earthquake wave number that pre-saves are according to being corresponded；Pre-stack depth migration processing unit, for using initially
Rate pattern, carries out pre-stack depth migration processing to the seismic data available for skew, generates the common imaging point road of seismic wave
Collection；Wherein, the initial velocity model includes constant velocity model；The abscissa of common imaging gather is offset distance；Common imaging point road
The ordinate integrated is depth；The numerical value of common imaging gather is seismic amplitude；Device also includes：Sum unit, for along
The X direction of common imaging gather is summed, and generates the pre-stack depth migration data of seismic wave.
The embodiment of the present invention brings following beneficial effect：
Seimic wave velocity inversion method and device provided in an embodiment of the present invention, by initial velocity model, to earth's surface
Seismic data carries out pre-stack depth migration processing, can obtain the common imaging gather of seismic wave；By from common imaging point road
The remaining depth data for extracting lineups is concentrated, and according to the remaining depth data, underground in geologic province can be finally inversed by and be situated between
The seismic data of matter；The seismic data accuracy that the mode of above-mentioned seimic wave velocity inverting is obtained is higher, and then improves
The resolution ratio of seismic imaging.
Other features and advantages of the present invention will be illustrated in the following description, also, partly be become from specification
Obtain it is clear that or being understood by implementing the present invention.The purpose of the present invention and other advantages are in specification, claims
And specifically noted structure is realized and obtained in accompanying drawing.
To enable the above objects, features and advantages of the present invention to become apparent, preferred embodiment cited below particularly, and coordinate
Appended accompanying drawing, is described in detail below.
Brief description of the drawings
, below will be to specific in order to illustrate more clearly of the specific embodiment of the invention or technical scheme of the prior art
The accompanying drawing to be used needed for embodiment or description of the prior art is briefly described, it should be apparent that, in describing below
Accompanying drawing is some embodiments of the present invention, for those of ordinary skill in the art, before creative work is not paid
Put, other accompanying drawings can also be obtained according to these accompanying drawings.
Fig. 1 is a kind of flow chart of seimic wave velocity inversion method provided in an embodiment of the present invention；
Fig. 2 be a kind of seimic wave velocity inversion method provided in an embodiment of the present invention in, using initial velocity model over the ground
Seismic wave data carry out pre-stack depth migration processing, generate the particular flow sheet of the common imaging gather of seismic wave；
Fig. 3 be a kind of seimic wave velocity inversion method provided in an embodiment of the present invention in, pass through seimic wave velocity inverting mould
Type, solves the particular flow sheet of the seismic data of underground medium in geologic province；
Fig. 4 is a kind of structural representation of seimic wave velocity inverting device provided in an embodiment of the present invention.
Embodiment
To make the purpose, technical scheme and advantage of the embodiment of the present invention clearer, below in conjunction with accompanying drawing to the present invention
Technical scheme be clearly and completely described, it is clear that described embodiment is a part of embodiment of the invention, rather than
Whole embodiments.Based on the embodiment in the present invention, those of ordinary skill in the art are not making creative work premise
Lower obtained every other embodiment, belongs to the scope of protection of the invention.
In view of existing seimic wave velocity inverting mode accuracy it is poor the problem of, the embodiments of the invention provide one
Seimic wave velocity inversion method and device are planted, the technology can apply in high-resolution seismic imaging, can also be applied to
In seismic exploration；The technology can be realized using related software and hardware, be described below by embodiment.
Embodiment one：
A kind of flow chart of seimic wave velocity inversion method shown in Figure 1, this method comprises the following steps：
Step S102, obtains the seismic data of earth's surface in default geologic province；Wherein, the seismic data bag of the earth's surface
Include shear wave data and P wave data；
Generally, seismic wave is, from earthquake source to the vibration blazed about, to refer to the elasticity for producing and being radiated to surrounding from focus
Ripple.Seismic wave can be divided into compressional wave (also referred to as P ripples), shear wave (also referred to as S ripples) and face ripple (L ripples) three types by circulation way；
Wherein, compressional wave and shear wave belong to bulk wave.When earthquake occurs, rupture and motion rapidly occur for the medium of focal area, this to disturb
Dynamic to constitute a wave source, due to the continuity of ball medium, just to earth interior and top layer, spread out is gone everywhere for fluctuation, is formd
Elastic wave in continuous media.
Step S104, using initial velocity model, pre-stack depth migration processing is carried out to seismic data, seismic wave is generated
Common imaging gather；Wherein, initial velocity model according to historical earthquake wave number according to and historical earthquake wave velocity corresponding relation
Set up；
Step S106, extracts the remaining depth data of lineups from common imaging gather；Wherein, the remaining depth data
Remaining depth correction amount and reflecting interface position including seismic wave；
Step S108, according to above-mentioned remaining depth data, the seismic data of underground medium in inverting geologic province；Its
In, the seismic data of the underground medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters., can when actually realizing
With the earthquake for the equations set up according to above-mentioned remaining depth data in the Mathematical Modeling related to seimic wave velocity, the model
Wave velocity and the direct relation of remaining depth data；By the Mathematical Modeling, it can obtain optimal by the way of iterative cycles
Wave speed of the earthquake angle value.
Seimic wave velocity inversion method provided in an embodiment of the present invention, by initial velocity model, to the seismic wave of earth's surface
Data carry out pre-stack depth migration processing, can obtain the common imaging gather of seismic wave；By being carried from common imaging gather
The remaining depth data of lineups is taken, and according to the remaining depth data, the ground of underground medium in geologic province can be finally inversed by
Seismic wave data；The seismic data accuracy that the mode of above-mentioned seimic wave velocity inverting is obtained is higher, so improve earthquake into
The resolution ratio of picture.
Poor in view of the primary data handlability that is gathered by sensing device, ground in geologic province is preset in above-mentioned acquisition
The seismic data of table, comprises the following steps：(1) the earthquake big gun for carrying subsurface geology information in default geologic province is gathered
Collect data；Wherein, the subsurface geology information includes geological structure information and geologic lithology change information；Specifically, the earthquake big gun
Collect data by being obtained after being pre-processed to 3-D seismics shot gather data, the earthquake shot gather data is available for direct imaging
Earthquake shot gather data；The pretreatment include to earthquake shot gather data carry out noise removal process, and to earthquake shot gather data with
The historical earthquake data pre-saved are corresponded；Further, the pretreatment can also include observation system loading；
(2) by divergence and curl operator, wave field separation processing, generation shear wave data and P wave data are carried out to earthquake shot gather data.
The above method carries out preliminary treatment by the earthquake shot gather data that will be got, and improves the follow-up handlability of data.
In a kind of seimic wave velocity inversion method shown in Figure 2, seismic data is entered using initial velocity model
The processing of row pre-stack depth migration, generates the particular flow sheet of the common imaging gather of seismic wave, and this method comprises the following steps：
Step S202, is pre-processed to seismic data, seismic data of the generation available for skew；Wherein, this is pre-
Processing includes carrying out seismic data noise remove, static correction and to seismic data and the historical earthquake pre-saved
Wave number is according to being corresponded；
Step S204, using initial velocity model, is carried out at pre-stack depth migration to the seismic data available for skew
Reason, generates the common imaging gather of seismic wave；Wherein, initial velocity model includes constant velocity model；The horizontal stroke of common imaging gather
Coordinate is offset distance；The ordinate of common imaging gather is depth；The numerical value of common imaging gather is seismic amplitude.
Above-mentioned steps S202 and step S204, can also be accomplished in the following manner：According to compressional wave and shear wave pretreatment ground
Data and the initial compressional wave and shear wave velocity model that provide are provided, pre-stack depth migration is carried out, compressional wave and shear wave imaging point altogether is obtained
Road collection.
After step s 204, above-mentioned seimic wave velocity inversion method also includes：Along the transverse axis side of common imaging gather
To being summed, the pre-stack depth migration data of seismic wave are generated.Wherein, asked along the X direction of common imaging gather
The mode of sum is referred to as being overlapped along the X direction of common imaging gather.
Further, the above-mentioned remaining depth data that lineups are extracted from common imaging gather, comprises the following steps：
(1) corresponding first depth value of each offset distance is extracted from common imaging gather；Wherein, the offset distance be geophone with
The distance between earthquake shot point；First depth value is handled by the corresponding geological data of offset distance through pre-stack depth migration
Arrive；(2) corresponding second depth value of zero-offset is extracted from common imaging gather；(3) the first depth value and second is calculated deeply
The difference of angle value, obtains the remaining depth correction amount of common imaging gather；(4) from pre-stack depth migration extracting data seismic wave
Reflecting interface, obtain reflecting interface position；Wherein, the reflecting interface position includes the locus coordinate of reflecting interface.
When implementing, the step in above-mentioned seimic wave velocity inversion method can respectively be carried out to compressional wave and shear wave
Processing；For example, in above-mentioned steps (1)-(4), the P wave data in common imaging gather can be handled first, the surplus of compressional wave is generated
Remaining depth correction amount and reflecting interface position；The shear wave data in common imaging gather are reprocessed, the remaining depth of shear wave is generated
Correcting value and reflecting interface position.
Further, it is above-mentioned according to remaining depth data, the seismic data of underground medium in inverting geologic province, including
Following steps：(1) grid discretization processing, formation speed mesh point are carried out to initial velocity model and reflecting interface position；(2)
By way of 2 ray traces, shot point and geophone station position are calculated respectively to the ray of pip；Wherein, pip position
In on reflecting interface position；The information of the ray at least includes ray traveltime and propagation path；(3) by the point on ray distribute to
With putting on closest speed mesh point；(4) the corresponding parameter of speed mesh point is set to underground medium in geologic province
Seismic data；(5) by seimic wave velocity inverse model, the seismic data of underground medium in geologic province is solved；Its
In, the seimic wave velocity inverse model is set up according to remaining depth data.When implementing, the earthquake of compressional wave can be first generated
Wave number evidence, regenerates the seismic data of shear wave.Then the remaining depth correction of the common imaging gather of compressional wave and shear wave is combined
Amount, the velocity amplitude on ray path grid and anisotropic parameters value set up seimic wave velocity inverse model and (are also referred to as elastic wave
Regularization velocity inversion model)；By non-linearization iterative seimic wave velocity inverse model, the speed of each mesh point is obtained
Degree and anisotropic parameters value.
Specifically, the seimic wave velocity inverse model can be expressed as：
Wherein, J (Δ λ) is object function；A is M1×M2The matrix of row N row；The element of matrix A is excursions depths to model
The derivative of parameterz0For remaining depth；xjFor the abscissa positions of inverted parameters；hkFor offset distance
Half；λiFor inverted parameters；M1For the quantity of the offset distance of common imaging gather；M2For the common imaging for renewal speed model
The quantity of point road collection；P represents compressional wave；S represents shear wave；Vectorial b includes M1×M2Individual element,Z is
Excursions depths；For the average value of different offset distance excursions depths；Vectorial Δ λ is model modification amount；Δλ0For initial model to
Amount；β is regularization parameter；| | | | norm is sought in representative.
In a kind of seimic wave velocity inversion method shown in Figure 3, by seimic wave velocity inverse model, geology is solved
The particular flow sheet of the seismic data of underground medium in region, this method comprises the following steps：
Step S302, sets initial temperature T0, initial solution Δ λ0, outer circulation iterations outn, interior circulation Markov chain it is long
Spend inn, temperature decline coefficient dec, the excursion a and minimum temperature t of model parametermin, i=0；
Step S304, carries out following outer circulation operations：Δ λ is calculated using Memoryless Quasi-Neton Methodi, preserve Δ λi；
Step S306, J (Δ λ are calculated by seimic wave velocity inverse modelι), preserve J (Δ λι)；
Step S308, sets k=0, tk=T0With scale factor b=0；
Step S310, carries out following interior circulation operations：Judge tk<tminWhether set up；If it is, performing step S326；Such as
It is really no, perform step S312；
Step S312, sets Δ λk=Δ λi；Calculate Δ λk+1=Δ λk+ξka；
Wherein, ξk=tk*sign(r-0.5)*[(1-1/tk)|2r-1|- 1], r is that uniformly divide random is obeyed between [0,1]
Number；Difference DELTA E=J (the Δ λ of calculating target function valuek+1)-J(Δλk)；
Step S318, judges whether Δ E≤0 sets up；If Δ E≤0, step S322 is performed；If Δ E>0, perform step
Rapid S320；
Step S320, judgesWhether set up；IfPerform step S322；IfPerform step S324；
Step S322, is set
Step S324, judges whether k=inn sets up；If not, performing step S310；If it is, performing step S325；
Step S325, sets tk+1=tk*(dec)b, b=b+1, k=k+1, perform step S310；
Step S326, carries out following update and operates：I=i+1 is set；
Step S328, judges i<Whether outn sets up；If it is, performing step S304；If not, performing step S330；
Step S330, from J (the Δ λ of preservationι) in value, filter out J (the Δ λ of minimumι) value；By minimum J (Δ λι) value pair
The velocity amplitude answered is defined as seismic data.
Optimal seismic wave velocity of longitudinal wave and shear wave velocity, earthquake wave number can be obtained above by the mode of loop iteration
It is higher according to accuracy, and then improve the resolution ratio of seismic imaging.
Embodiment two：
Corresponding to above method embodiment, a kind of structural representation of seimic wave velocity inverting device shown in Figure 4,
The device includes following part：
Seismic data acquisition module 402, the seismic data for obtaining earth's surface in default geologic province；Wherein, should
The seismic data of earth's surface includes shear wave data and P wave data；
Pre-stack depth migration processing module 404, for utilizing initial velocity model, prestack depth is carried out to seismic data
Migration processing, generates the common imaging gather of seismic wave；Wherein, the initial velocity model is according to historical earthquake wave number evidence and history
The corresponding relation of seimic wave velocity is set up；
Remaining depth data extraction module 406, the remaining depth data for extracting lineups from common imaging gather；
Wherein, the remaining depth data includes remaining depth correction amount and the reflecting interface position of seismic wave；
Seismic data inverting module 408, for according to remaining depth data, the ground of underground medium in inverting geologic province
Seismic wave data；Wherein, the seismic data of the underground medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters.
Seimic wave velocity inverting device provided in an embodiment of the present invention, by initial velocity model, to the seismic wave of earth's surface
Data carry out pre-stack depth migration processing, can obtain the common imaging gather of seismic wave；By being carried from common imaging gather
The remaining depth data of lineups is taken, and according to the remaining depth data, the ground of underground medium in geologic province can be finally inversed by
Seismic wave data；The seismic data accuracy that the mode of above-mentioned seimic wave velocity inverting is obtained is higher, so improve earthquake into
The resolution ratio of picture.
Poor, the above-mentioned seismic data acquisition module bag in view of the primary data handlability that is gathered by sensing device
Include following part：(1) collecting unit, for gathering the earthquake big gun collection number for carrying subsurface geology information in default geologic province
According to；Wherein, the subsurface geology information includes geological structure information and geologic lithology change information；(2) wave field separation processing unit,
For by divergence and curl operator, wave field separation processing, generation shear wave data and P wave data to be carried out to earthquake shot gather data.
The above method carries out preliminary treatment by the earthquake shot gather data that will be got, and improves the follow-up handlability of data.
Further, above-mentioned pre-stack depth migration processing module includes：(1) pretreatment unit, for seismic data
Pre-processed, seismic data of the generation available for skew；Wherein, the pretreatment includes removing seismic data progress noise
Remove, static correction and to seismic data and the historical earthquake wave number that pre-saves according to corresponding；(2) prestack depth
Migration processing unit, for utilizing initial velocity model, is carried out at pre-stack depth migration to the seismic data available for skew
Reason, generates the common imaging gather of seismic wave；Wherein, the initial velocity model includes constant velocity model；The common imaging gather
Abscissa be offset distance；The ordinate of the common imaging gather is depth；The numerical value of the common imaging gather shakes for seismic wave
Width；Above-mentioned seimic wave velocity inverting device also includes：Sum unit, for being asked along the X direction of common imaging gather
With the pre-stack depth migration data of generation seismic wave.
Seimic wave velocity inversion method and the computer program product of device that the embodiment of the present invention is provided, including storage
The computer-readable recording medium of program code, the instruction that described program code includes can be used for performing previous methods embodiments
Described in method, implement can be found in embodiment of the method, will not be repeated here.
If the function is realized using in the form of SFU software functional unit and is used as independent production marketing or in use, can be with
It is stored in a computer read/write memory medium.Understood based on such, technical scheme is substantially in other words
The part contributed to prior art or the part of the technical scheme can be embodied in the form of software product, the meter
Calculation machine software product is stored in a storage medium, including some instructions are make it that a computer equipment (can be individual
People's computer, server, or network equipment etc.) perform all or part of step of each of the invention embodiment methods described.
And foregoing storage medium includes：USB flash disk, mobile hard disk, read-only storage (ROM, Read-Only Memory), arbitrary access are deposited
Reservoir (RAM, Random Access Memory), magnetic disc or CD etc. are various can be with the medium of store program codes.
Finally it should be noted that：Embodiment described above, is only the embodiment of the present invention, to illustrate the present invention
Technical scheme, rather than its limitations, protection scope of the present invention is not limited thereto, although with reference to the foregoing embodiments to this hair
It is bright to be described in detail, it will be understood by those within the art that：Any one skilled in the art
The invention discloses technical scope in, it can still modify to the technical scheme described in previous embodiment or can be light
Change is readily conceivable that, or equivalent is carried out to which part technical characteristic；And these modifications, change or replacement, do not make
The essence of appropriate technical solution departs from the spirit and scope of embodiment of the present invention technical scheme, should all cover the protection in the present invention
Within the scope of.Therefore, protection scope of the present invention described should be defined by scope of the claims.
Claims (10)
1. a kind of seimic wave velocity inversion method, it is characterised in that including：
Obtain the seismic data of earth's surface in default geologic province；Wherein, the seismic data of the earth's surface includes shear wave data
And P wave data；
Using initial velocity model, pre-stack depth migration processing is carried out to the seismic data, the common imaging of seismic wave is generated
Point road collection；Wherein, the initial velocity model is set up according to historical earthquake wave number according to the corresponding relation of historical earthquake wave velocity；
The remaining depth data of lineups is extracted from the common imaging gather；Wherein, the remaining depth data includes institute
State remaining depth correction amount and the reflecting interface position of seismic wave；
According to the remaining depth data, described in inverting in geologic province underground medium seismic data；Wherein, the underground
The seismic data of medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters.
2. according to the method described in claim 1, it is characterised in that the earthquake wave number for obtaining earth's surface in default geologic province
According to including：
The earthquake shot gather data for carrying subsurface geology information in the default geologic province of collection；Wherein, the subsurface geology letter
Breath includes geological structure information and geologic lithology change information；
By divergence and curl operator, wave field separation processing, generation shear wave data and compressional wave are carried out to the earthquake shot gather data
Data.
3. according to the method described in claim 1, it is characterised in that the utilization initial velocity model, to the earthquake wave number
According to pre-stack depth migration processing is carried out, the common imaging gather of generation seismic wave includes：
The seismic data is pre-processed, seismic data of the generation available for skew；Wherein, the pretreatment includes
Noise remove, static correction are carried out to the seismic data and to the seismic data and the historical earthquake pre-saved
Wave number is according to being corresponded；
Using initial velocity model, pre-stack depth migration processing, generation ground are carried out to the seismic data that can be used for skew
The common imaging gather of seismic wave；Wherein, the initial velocity model includes constant velocity model；The horizontal seat of the common imaging gather
It is designated as offset distance；The ordinate of the common imaging gather is depth；The numerical value of the common imaging gather is seismic amplitude；
After the common imaging gather of the generation seismic wave, methods described also includes：
Summed along the X direction of the common imaging gather, generate the pre-stack depth migration data of the seismic wave.
4. method according to claim 3, it is characterised in that described that lineups are extracted from the common imaging gather
Remaining depth data includes：
Corresponding first depth value of each offset distance is extracted from the common imaging gather；Wherein, the offset distance is earthquake
The distance between wave detector and earthquake shot point；First depth value is deep through prestack by the corresponding geological data of the offset distance
Degree migration processing is obtained；
Corresponding second depth value of zero-offset is extracted from the common imaging gather；
The difference of first depth value and second depth value is calculated, the remaining depth school of the common imaging gather is obtained
Positive quantity；
From the reflecting interface of seismic wave described in the pre-stack depth migration extracting data, reflecting interface position is obtained；Wherein, institute
Stating reflecting interface position includes the locus coordinate of the reflecting interface.
5. method according to claim 4, it is characterised in that described according to the remaining depth data, described in inverting
The seismic data of underground medium includes in matter region：
Grid discretization processing, formation speed mesh point are carried out to the initial velocity model and the reflecting interface position；
By way of 2 ray traces, shot point and geophone station position are calculated respectively to the ray of pip；Wherein, it is described anti-
Exit point is located on the reflecting interface position；The information of the ray at least includes ray traveltime and propagation path；
By on point distribution to the speed mesh point closest with the point on the ray；
The corresponding parameter of the speed mesh point is set to the seismic data of underground medium in the geologic province；
By seimic wave velocity inverse model, the seismic data of underground medium in the geologic province is solved；Wherein, describedly
Seismic velocity inverse model is set up according to the remaining depth data.
6. method according to claim 5, it is characterised in that the seimic wave velocity inverse model includes：
Wherein, J (Δ λ) is object function；A is M1×M2The matrix of row N row；The element of matrix A is excursions depths to model parameter
Derivativez0For the remaining depth；xjFor the abscissa positions of inverted parameters；hkFor offset distance
Half；λiFor inverted parameters；M1For the quantity of the offset distance of the common imaging gather；M2For being total to for renewal speed model
The quantity of image gather；P represents compressional wave；S represents shear wave；Vectorial b includes M1×M2Individual element,Z is excursions depths；For the average value of different offset distance excursions depths；Vectorial Δ λ be model more
New amount；Δλ0For initial model vector；β is regularization parameter；| | | | norm is sought in representative.
7. method according to claim 5, it is characterised in that described by seimic wave velocity inverse model, is solved described
The seismic data of underground medium includes in geologic province：
Initial temperature T is set0, initial solution Δ λ0, outer circulation iterations outn, interior circulation Markov chain length inn, temperature decline
Subtract coefficient d ec, the excursion a and minimum temperature t of model parametermin, i=0；
Carry out following outer circulation operations：
Δ λ is calculated using Memoryless Quasi-Neton Methodi, preserve the Δ λi；
J (Δ λ are calculated by the seimic wave velocity inverse modelι), preserve J (the Δ λι)；
K=0, t are setk=T0With scale factor b=0；
Carry out following interior circulation operations：
Judge tk<tminWhether set up；Operated if it is, performing following update；If not, setting Δ λk=Δ λi；
Calculate Δ λk+1=Δ λk+ξka；Wherein, ξk=tk*sign(r-0.5)*[(1-1/tk)|2r-1|- 1], r is clothes between [0,1]
From the random number uniformly divided；
Difference DELTA E=J (the Δ λ of calculating target function valuek+1)-J(Δλk)；
Judge whether Δ E≤0 sets up；
If Δ E≤0, Δ λ is setk+1=Δ λk；If Δ E>0, judgeWhether set up；
IfΔ λ is setk+1=Δ λk；IfJudge whether k=inn sets up；If not, continuing
Perform above-mentioned interior circulation operation；If it is, setting tk+1=tk*(dec)b, b=b+1, k=k+1；Continue executing with above-mentioned interior circulation
Operation；
Following update is carried out to operate：
I=i+1 is set；
Judge i<Whether outn sets up；If it is, continuing executing with above-mentioned outer circulation operation；If not, from J (the Δ λ of preservationι) value
In, filter out J (the Δ λ of minimumι) value；
By minimum J (the Δ λι) the corresponding velocity amplitude of value is defined as seismic data.
8. a kind of seimic wave velocity inverting device, it is characterised in that including：
Seismic data acquisition module, the seismic data for obtaining earth's surface in default geologic province；Wherein, the earth's surface
Seismic data includes shear wave data and P wave data；
Pre-stack depth migration processing module, for utilizing initial velocity model, carries out prestack depth inclined to the seismic data
Shifting is handled, and generates the common imaging gather of seismic wave；Wherein, the initial velocity model is according to historical earthquake wave number evidence and history
The corresponding relation of seimic wave velocity is set up；
Remaining depth data extraction module, the remaining depth data for extracting lineups from the common imaging gather；Its
In, the remaining depth data includes remaining depth correction amount and the reflecting interface position of the seismic wave；
Seismic data inverting module, for according to the remaining depth data, underground medium in geologic province described in inverting
Seismic data；Wherein, the seismic data of the underground medium includes shear wave velocity, velocity of longitudinal wave and anisotropic parameters.
9. device according to claim 8, it is characterised in that the seismic data acquisition module includes：
Collecting unit, for gathering the earthquake shot gather data for carrying subsurface geology information in default geologic province；Wherein, institute
Stating subsurface geology information includes geological structure information and geologic lithology change information；
Wave field separation processing unit, for by divergence and curl operator, being carried out to the earthquake shot gather data at wave field separation
Reason, generation shear wave data and P wave data.
10. device according to claim 8, it is characterised in that the pre-stack depth migration processing module includes：
Pretreatment unit, for being pre-processed to the seismic data, seismic data of the generation available for skew；Its
In, the pretreatment includes carrying out the seismic data noise remove, static correction and to the seismic data and in advance
The historical earthquake wave number first preserved is according to being corresponded；
Pre-stack depth migration processing unit, for utilizing initial velocity model, enters to the seismic data that can be used for skew
The processing of row pre-stack depth migration, generates the common imaging gather of seismic wave；Wherein, the initial velocity model includes constant velocity mould
Type；The abscissa of the common imaging gather is offset distance；The ordinate of the common imaging gather is depth；The common imaging
The numerical value of point road collection is seismic amplitude；
Described device also includes：Sum unit, is summed for the X direction along the common imaging gather, generates institute
State the pre-stack depth migration data of seismic wave.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710168052.3A CN106959467B (en) | 2017-03-20 | 2017-03-20 | Seimic wave velocity inversion method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710168052.3A CN106959467B (en) | 2017-03-20 | 2017-03-20 | Seimic wave velocity inversion method and device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106959467A true CN106959467A (en) | 2017-07-18 |
| CN106959467B CN106959467B (en) | 2019-02-19 |
Family
ID=59470528
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710168052.3A Expired - Fee Related CN106959467B (en) | 2017-03-20 | 2017-03-20 | Seimic wave velocity inversion method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106959467B (en) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108710148A (en) * | 2018-05-29 | 2018-10-26 | 中国科学院地质与地球物理研究所 | The steady phase prestack depth migration method in three-dimensional dip domain and device |
| CN111596346A (en) * | 2019-02-20 | 2020-08-28 | 中国石油天然气集团有限公司 | Elastic wave velocity inversion method and device |
| CN111971586A (en) * | 2018-03-30 | 2020-11-20 | Bp北美公司 | Seismic velocity derived hydrocarbon indication |
| CN113031073A (en) * | 2019-12-24 | 2021-06-25 | 中国石油天然气股份有限公司 | Method, device and equipment for determining prestack depth migration speed |
| CN113126149A (en) * | 2018-12-28 | 2021-07-16 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102841375A (en) * | 2012-09-06 | 2012-12-26 | 中国石油大学(华东) | Method for tomography velocity inversion based on angle domain common imaging gathers under complicated condition |
| CN104268404A (en) * | 2014-09-25 | 2015-01-07 | 北京航空航天大学 | Automatic pick-up method of residual depth |
| CN104749631A (en) * | 2015-03-11 | 2015-07-01 | 中国科学院地质与地球物理研究所 | Sparse inversion based migration velocity analysis method and device |
| CN105572728A (en) * | 2014-10-16 | 2016-05-11 | 中国石油化工股份有限公司 | Reverse illumination speed inversion method based on least-square objective functional |
| US20160341839A1 (en) * | 2015-05-20 | 2016-11-24 | Conocophillips Company | Surface wave tomography using sparse data acquisition |
-
2017
- 2017-03-20 CN CN201710168052.3A patent/CN106959467B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102841375A (en) * | 2012-09-06 | 2012-12-26 | 中国石油大学(华东) | Method for tomography velocity inversion based on angle domain common imaging gathers under complicated condition |
| CN104268404A (en) * | 2014-09-25 | 2015-01-07 | 北京航空航天大学 | Automatic pick-up method of residual depth |
| CN105572728A (en) * | 2014-10-16 | 2016-05-11 | 中国石油化工股份有限公司 | Reverse illumination speed inversion method based on least-square objective functional |
| CN104749631A (en) * | 2015-03-11 | 2015-07-01 | 中国科学院地质与地球物理研究所 | Sparse inversion based migration velocity analysis method and device |
| US20160341839A1 (en) * | 2015-05-20 | 2016-11-24 | Conocophillips Company | Surface wave tomography using sparse data acquisition |
Non-Patent Citations (2)
| Title |
|---|
| CAIXIA YU 等: ""Seismic detection method for small-scale discontinuities based on dictionary learning and sparse representation"", 《JOURNAL OF APPLIED GEOPHYSICS》 * |
| 于彩霞 等: ""VTI介质快速偏移速度分析"", 《石油地球物理勘探》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111971586A (en) * | 2018-03-30 | 2020-11-20 | Bp北美公司 | Seismic velocity derived hydrocarbon indication |
| CN108710148A (en) * | 2018-05-29 | 2018-10-26 | 中国科学院地质与地球物理研究所 | The steady phase prestack depth migration method in three-dimensional dip domain and device |
| CN108710148B (en) * | 2018-05-29 | 2019-05-24 | 中国科学院地质与地球物理研究所 | The steady phase prestack depth migration method in three-dimensional dip domain and device |
| CN113126149A (en) * | 2018-12-28 | 2021-07-16 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
| CN113126149B (en) * | 2018-12-28 | 2024-04-09 | 中国石油化工股份有限公司 | Method and system for seismic image processing to enhance geological structure fidelity |
| CN111596346A (en) * | 2019-02-20 | 2020-08-28 | 中国石油天然气集团有限公司 | Elastic wave velocity inversion method and device |
| CN113031073A (en) * | 2019-12-24 | 2021-06-25 | 中国石油天然气股份有限公司 | Method, device and equipment for determining prestack depth migration speed |
| CN113031073B (en) * | 2019-12-24 | 2024-03-01 | 中国石油天然气股份有限公司 | Pre-stack depth migration speed determination method, device and equipment |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106959467B (en) | 2019-02-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106772583B (en) | A kind of earthquake diffracted wave separation method and device | |
| CN106959467A (en) | Seimic wave velocity inversion method and device | |
| CN105277978B (en) | A kind of method and device for determining near-surface velocity model | |
| CN108369289A (en) | The method for analyzing design geophysical reconnaissance using full wave field inversion point spread function | |
| CN104237940A (en) | Diffracted wave imaging method and diffracted wave imaging device based on dynamical features | |
| CN106526674A (en) | Three-dimensional full waveform inversion energy weighted gradient preprocessing method | |
| US6018499A (en) | Three-dimensional seismic imaging of complex velocity structures | |
| CN102944896A (en) | Model method static correction method for surface survey data | |
| CN101984366A (en) | Fluctuating surface pre-stack time migration method and device based on common aperture surface | |
| CN105319589A (en) | Full-automatic three-dimensional tomography inversion method using a local event slope | |
| Tran et al. | 3-D time-domain Gauss–Newton full waveform inversion for near-surface site characterization | |
| CN110471113A (en) | Bearing calibration, device and storage medium are moved in inverting based on unstable state seismic data | |
| CN103926619A (en) | Reverse time migration method of three-dimensional VSP data | |
| CN107462924A (en) | A kind of absolute wave impedance inversion method independent of well-log information | |
| CN103119472B (en) | Utilize simultaneously and order source method carries out the mixed method of full waveform inversion | |
| CN109633749A (en) | Non-linear Fresnel zone seismic traveltime tomography method based on scattering integral method | |
| CN106772593A (en) | The imaging method and device of diffracted wave | |
| CN102053269A (en) | Analysis method of speed in seismic data | |
| Guo et al. | Topography-dependent eikonal tomography based on the fast-sweeping scheme and the adjoint-state technique | |
| US5764514A (en) | Method for modelling kinematic seismic data processed with at least one motion operator | |
| Sambolian et al. | Consistent seismic event location and subsurface parameters inversion through slope tomography: a variable-projection approach | |
| CN111399037B (en) | Method and apparatus for high speed top interface extraction | |
| CN108957539B (en) | Ray tracing method and device in chromatography migration velocity analysis | |
| CN104950330B (en) | The velocity modeling method and system of gas cap reservoir depth domain imaging | |
| CN107918147A (en) | Diffraction wave imaging method and device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20181116 Address after: 100120 Xicheng District six paw Kang in Beijing City Applicant after: CHINA NATIONAL PETROLEUM Corp. Applicant after: INSTITUTE OF GEOLOGY AND GEOPHYSICS, CHINESE ACADEMY OF SCIENCES Applicant after: BGP INC., CHINA NATIONAL PETROLEUM Corp. Address before: No. 19, North Tu Cheng West Road, Chaoyang District, Beijing Applicant before: Institute of Geology and Geophysics, Chinese Academy of Sciences |
|
| 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: 20190219 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |