CN109446962A - Land cover interannual variance detection method and system - Google Patents
Land cover interannual variance detection method and system Download PDFInfo
- Publication number
- CN109446962A CN109446962A CN201811222022.7A CN201811222022A CN109446962A CN 109446962 A CN109446962 A CN 109446962A CN 201811222022 A CN201811222022 A CN 201811222022A CN 109446962 A CN109446962 A CN 109446962A
- Authority
- CN
- China
- Prior art keywords
- year
- ndvi
- target
- timing
- target area
- 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
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/13—Satellite images
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V20/00—Scenes; Scene-specific elements
- G06V20/10—Terrestrial scenes
- G06V20/194—Terrestrial scenes using hyperspectral data, i.e. more or other wavelengths than RGB
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Astronomy & Astrophysics (AREA)
- Remote Sensing (AREA)
- Multimedia (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The embodiment of the invention provides a kind of Land cover interannual variance detection method and systems, the timing remotely-sensed data by target area in first object year and the 1 target year, determine target area in the NDVI time series data in first object year and the 1 target year respectively；On the basis of default remote sensing temporal model carries out the NDVI time series data reconstruct of two different years, the similarity system design of the timing curve of two different years is carried out using DTW algorithm, obtains the Annual variations testing result of land cover pattern.The use value that timing remotely-sensed data has sufficiently been excavated in the embodiment of the present invention has great application prospect in terms of using remote sensing timing remote sensing imagery change detection.And the displacement problem on time shaft can be overcome by DTW algorithm to a certain extent, it can reduce the influence of the exceptional values such as noise and cloud, obtain better matching effect, improve the detection accuracy of Land cover interannual variance.
Description
Technical field
The present embodiments relate to remote sensing image change detection techniques fields, become more particularly, to land cover pattern year border
Change detection method and system.
Background technique
Land cover pattern variation is the combined reaction of earth's surface lot of essential factors, to Land surface energy budget, climate change, water circulation etc.
It has a direct impact, is numerous subject focus of attention.Remotely-sensed data is that soil covers because of its broad perspectives, periodicity and continuity
The significant data source of lid variation identification.Timing remote sensing image is capable of providing the record of table status long-term sequence over the ground, reflection ground
Dynamic changes of the table in a Long time scale provide reliable data for the space-time identification of land cover pattern variation
Source.The temporal dimension information that remote sensing time series data collection contains sufficiently is excavated, it is automatic to develop the land cover pattern variation based on timing remote sensing
Detection method is one of the hot spot of remote sensing technology area research in recent years.
During sequence remotely-sensed data carries out land cover pattern variation detection when in use, need to eliminate in remotely-sensed data first
The singular value to be formed is influenced by cloud, Yun Ying, mist, haze etc., that is, eliminates the noise of timing curve, this process is exactly remote sensing time sequence
The process of column data filtering and the process of temporal model building.For realize land cover pattern variation automatic detection, with multidate
Remotely-sensed data is information source, have been developed image difference method, Aberrant spectrum method, Principal Component Analysis, pseudo color composing method,
A variety of methods such as wave band Shift Method, classification and predicting method, multiband cross correlation analysis method and Change vector Analysis method, these sides
Method is each advantageous, however land cover pattern variation detection still has difficulty: first is that " the different spectrum of jljl, same object different images " phenomenon makes ground
Table variation cannot be corresponded with image variation characteristic；Second is that land cover pattern variation includes changing between variation and class in class, i.e., together
Land cover pattern variation caused by the land cover pattern variation and vegetation classification that plantation is caused change, variation degree also include
Part variation and all variation, for different goals in research, how to define variation is also a difficult point；Third is that land cover pattern becomes
The unstability of change, so that variation detection is more complicated.
The pass of detection is changed using the time response that land cover pattern embodies in the remote sensing image of long-term sequence
Key is the similarity analysis for comparing two time serieses, i.e. the timing curve of quantitative assessment pixel to be detected and object reference pixel
Similarity between timing curve.The method that time series similarity calculates generally takes Euclidean distance method, Euclidean distance meter
What is calculated is the actual distance of two points of synchronization, but during detection land cover pattern variation, since surface temperature, crop are raw
The reasons such as length, feature changes can occur flexible and be deviated on axis, therefore in the timing curve of different year, atural object in the time
The displacement problem of timing curve on a timeline can not be overcome using Euclidean distance method.
Summary of the invention
In order to overcome the problems referred above or it at least is partially solved the above problem, the embodiment of the invention provides a kind of soils to cover
Lid Annual variations detection method and system.
In a first aspect, the embodiment of the invention provides a kind of Land cover interannual variance detection methods, comprising:
Based on target area in the timing remotely-sensed data in first object year and the 1 target year, the target area is determined respectively
Vegetation index NDVI time series data of the domain in the first object year and the second target year；
Based on default remote sensing temporal model, respectively to the target area in the first object year and second target
Year NDVI time series data be fitted, determine the target area the first object year the first NDVI timing curve with
And the 2nd NDVI timing curve in the second target year；
Based on dynamic time warping DTW algorithm, overture when determining the first NDVI timing curve and two NDVI
DTW distance value between line, and determine the target area in the first object year and described based on the DTW distance value
The Annual variations testing result in 2 target years.
Second aspect, the embodiment of the invention provides a kind of Land cover interannual variance detection systems, comprising: NDVI timing
Data acquisition module, NDVI timing curve obtain module and Annual variations testing result obtains module.Wherein,
NDVI time series data obtains module for distant in the timing in first object year and the 1 target year based on target area
Feel data, determines that normalization difference vegetation of the target area in the first object year and the second target year refers to respectively
Number NDVI time series data；
NDVI timing curve obtains module and is used for based on default remote sensing temporal model, respectively to the target area described
The NDVI time series data in first object year and the second target year is fitted, and determines the target area in first mesh
Mark the first NDVI timing curve in year and the 2nd NDVI timing curve in the second target year；
Annual variations testing result obtains module and is used to be based on dynamic time warping DTW algorithm, determines the first NDVI
DTW distance value between timing curve and the 2nd NDVI timing curve, and the target is determined based on the DTW distance value
Annual variations testing result of the region in the first object year and the second target year.
The third aspect, the embodiment of the invention provides a kind of electronic equipment, comprising:
At least one processor, at least one processor, communication interface and bus；Wherein,
The processor, memory, communication interface complete mutual communication by the bus；
The memory is stored with the program instruction that can be executed by the processor, and the processor calls described program to refer to
It enables, to execute the Land cover interannual variance detection method of first aspect offer.
Fourth aspect, the embodiment of the invention provides a kind of non-transient computer readable storage medium, the non-transient meter
Calculation machine readable storage medium storing program for executing stores computer instruction, and the computer instruction makes the computer execute the soil that first aspect provides
Ground covers Annual variations detection method.
A kind of Land cover interannual variance detection method and system provided in an embodiment of the present invention, by target area
The timing remotely-sensed data in 1 target year and the 1 target year determines target area in first object year and the 1 target year respectively
NDVI time series data；On the basis of default remote sensing temporal model carries out the NDVI time series data reconstruct of two different years, benefit
The similarity system design that the timing curve of two different years is carried out with DTW algorithm obtains the Annual variations detection knot of land cover pattern
Fruit.The use value that timing remotely-sensed data has sufficiently been excavated in the embodiment of the present invention is using the variation inspection of remote sensing timing image
Aspect is surveyed, there is great application prospect.And the displacement on time shaft can be overcome to ask to a certain extent by DTW algorithm
Topic, can reduce the influence of the exceptional values such as noise and cloud, obtains better matching effect, improve the inspection of Land cover interannual variance
Survey precision.
Detailed description of the invention
In order to more clearly explain the embodiment of the invention or the technical proposal in the existing technology, to embodiment or will show below
There is attached drawing needed in technical description to be briefly described, it should be apparent that, the accompanying drawings in the following description is this hair
Bright some embodiments for those of ordinary skill in the art without creative efforts, can be with root
Other attached drawings are obtained according to these attached drawings.
Fig. 1 is the flow diagram of the Land cover interannual variance detection method provided in the embodiment of the present invention；
Fig. 2 is the structural schematic diagram of the Land cover interannual variance detection system provided in the embodiment of the present invention；
Fig. 3 is the structural schematic diagram of a kind of electronic equipment provided in the embodiment of the present invention.
Specific embodiment
In order to make the object, technical scheme and advantages of the embodiment of the invention clearer, below in conjunction with the embodiment of the present invention
In attached drawing, technical scheme in the embodiment of the invention is clearly and completely described, it is clear that described embodiment is
A part of the embodiment of the present invention, instead of all the embodiments.Based on the embodiments of the present invention, those of ordinary skill in the art
Every other embodiment obtained without creative efforts, shall fall within the protection scope of the present invention.
In the description of the embodiment of the present invention, it should be noted that term " center ", "upper", "lower", "left", "right",
The orientation or positional relationship of the instructions such as "vertical", "horizontal", "inner", "outside" is to be based on the orientation or positional relationship shown in the drawings,
It is merely for convenience of the description embodiment of the present invention and simplifies description, rather than the device or element of indication or suggestion meaning must have
There is specific orientation, be constructed and operated in a specific orientation, therefore should not be understood as the limitation to the embodiment of the present invention.In addition,
Term " first ", " second ", " third " are used for descriptive purposes only and cannot be understood as indicating or suggesting relative importance.
In the description of the embodiment of the present invention, it should be noted that unless otherwise clearly defined and limited, term " peace
Dress ", " connected ", " connection " shall be understood in a broad sense, for example, it may be being fixedly connected, may be a detachable connection, or integrally
Connection；It can be mechanical connection, be also possible to be electrically connected；Can be directly connected, can also indirectly connected through an intermediary,
It can be the connection inside two elements.For the ordinary skill in the art, above-mentioned art can be understood with concrete condition
The concrete meaning of language in embodiments of the present invention.
Fig. 1 is the flow diagram of the Land cover interannual variance detection method provided in the embodiment of the present invention, such as Fig. 1 institute
Show, Land cover interannual variance detection method provided in an embodiment of the present invention specifically includes:
S1 determines the mesh based on target area in the timing remotely-sensed data in first object year and the 1 target year respectively
Region is marked in the vegetation index NDVI time series data in the first object year and the second target year；
S2, based on default remote sensing temporal model, respectively to the target area in the first object year and described second
The NDVI time series data in target year is fitted, and determines target area overture in a NDVI in the first object year
Line and the 2nd NDVI timing curve in the second target year；
S3 is based on dynamic time warping DTW algorithm, determines the first NDVI timing curve and the 2nd NDVI timing
DTW distance value between curve, and determine the target area in first object year and described based on the DTW distance value
The Annual variations testing result in the 1 target year.
Specifically, the purpose of the embodiment of the present invention is that according to obtained target area two different years timing remote sensing
Data, determine whether land cover pattern is changed between the two different times.Such as first object year is 2008, the
2 target years were 2009, then obtained target area first in timing remotely-sensed data in 2008 and target area in 2009
Timing remotely-sensed data.Here timing remotely-sensed data refers to the remotely-sensed data being sequentially arranged, in particular to obtains
It the remotely-sensed data of every day and is sequentially arranged in 1 target year (i.e. 2008), obtains the 1 target year (i.e. 2009)
The remotely-sensed data of middle every day is simultaneously equally sequentially arranged.The method that timing remotely-sensed data is obtained in the embodiment of the present invention can
To be determined as needed, the acquisition of Landsat satellite preferably can choose.Remotely-sensed data refers to spectroscopic data,
The reflectivity of earth's surface at each wave band got.
Since NDVI can eliminate partial radiation error, noise, expand differently species it is other between difference, therefore choose
NDVI time series data is studied.According to target area in the timing remotely-sensed data in first object year, target area can be determined
Domain first object year vegetation index (Normalized Difference Vegetation Index,
NDVI) time series data；Similarly, target area can be determined in the timing remotely-sensed data in the 1 target year according to target area
In the NDVI time series data in the 1 target year.Determine that the implementation of NDVI time series data can pass through according to timing remotely-sensed data
The calculation formula of NDVI realizes that NDVI is equal to the difference of the remotely-sensed data of near infrared band and the remotely-sensed data of red spectral band, removes
And, it is embodied as with formula with the remotely-sensed data of the two wave bands:
Wherein, RNIRIndicate the Reflectivity for Growing Season of near infrared band, the i.e. remotely-sensed data of near infrared band, RRedIndicate feux rouges
The Reflectivity for Growing Season of wave band, the i.e. remotely-sensed data of red spectral band.
To different classes of land cover pattern, the value of NDVI is between -1~1: for water body, the value of NDVI is usually
Negative value；For bare area and rock, the value of NDVI is usually 0 or so；The farmland high for vegetation coverage, forest land, NDVI's
Value increases as vegetation coverage increases.
Determining target area in the NDVI timing in first object year and the 1 target year respectively in the embodiment of the present invention
After data, based on default remote sensing temporal model, respectively to target area first object year and the 1 target year NDVI timing
Data are fitted, and determine target area in the first NDVI timing curve in first object year and the second of the 1 target year
NDVI timing curve.Described default remote sensing temporal model in the embodiment of the present invention refers to ordinal number when can be to discrete NDVI
According to the continuous function model being fitted.When by presetting NDVI of the remote sensing temporal model to first object year and the 1 target year
Ordinal number determines that the parameter of the corresponding default remote sensing temporal model of NDVI time series data in first object year takes according to being fitted respectively
The parameter value of the corresponding default remote sensing temporal model of NDVI time series data in value and the 1 target year, when obtaining a NDVI
Overture line and the 2nd NDVI timing curve.
After the first NDVI timing curve and the 2nd NDVI timing curve has been determined, it is based on dynamic time warping (Dynamic
Time Warping, DTW) algorithm, determine the DTW distance value between the first NDVI timing curve and the 2nd NDVI timing curve,
And determine target area in the Annual variations testing result in first object year and the 1 target year based on DTW distance value.DTW algorithm
Thought based on Dynamic Programming (Dynamic Programming, DP) is used to when the time, warping function was met certain condition
The one-to-one relationship of two sequences of input in time is described, guarantees that cumulative distance is minimum when matching the two.It has determined
DTW distance value can determine the similarity between the NDVI timing curve in two target years according to DTW distance value, and then really
Target area is made in the Annual variations testing result in first object year and the 1 target year.Year border in the embodiment of the present invention becomes
Change testing result and generally include two kinds: one is changing, i.e., being sent out in the land cover pattern in first object year and the 1 target year
Changing；One is not changing, i.e., the land cover pattern in first object year and the 1 target year does not change.Similarity
More low, expression changes, and the similarity the high, indicates not change.
The Land cover interannual variance detection method provided in the embodiment of the present invention, by target area in first object year
With the timing remotely-sensed data in the 1 target year, determine target area in the NDVI timing in first object year and the 1 target year respectively
Data；On the basis of default remote sensing temporal model carries out the NDVI time series data reconstruct of two different years, DTW algorithm is utilized
The similarity system design for carrying out the timing curve of two different years obtains the Annual variations testing result of land cover pattern.The present invention
The use value that timing remotely-sensed data has sufficiently been excavated in embodiment, in terms of using remote sensing timing remote sensing imagery change detection, tool
There is great application prospect.And the displacement problem on time shaft can be overcome by DTW algorithm to a certain extent, it can drop
The influence of the exceptional values such as low noise and cloud, obtains better matching effect, improves the detection accuracy of Land cover interannual variance.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention, institute
Stating default remote sensing temporal model is the function model determined by sine term, cosine term, constant term, first order and quadratic term.
Specifically, due in the prior art remote sensing time series data apply during, all multiple timings weights of appearance
Algorithm is built, specifically includes that maximum value synthetic method, at the appointed time in section, by the value of each point in time series, according to big float
Column sequence, takes value of the maximum value as new images corresponding to the position in time series.This method is in certain journey
The interference of cloud and atmosphere can be eliminated on degree, but is handled simply, and the influence of earth's surface bidirectional reflectance is not fully considered, can be damaged
Lose many useful informations.Optimum index slope extraction method, it is believed that coupling relationship be it is stable, when occurring being mutated very in time series
It may be the variation as caused by cloud or sensor visual angle etc., therefore use sliding time window, to judge thresholding, recognition time sequence
Noise in column.This method does not influence gradual change timing, only influences to decline suddenly in timing curve, then gradually rise again
Part, cause higher timing values cannot be removed effectively atmospheric conditions, sliding window threshold value setting also needs to test repeatedly.
Mean iterative filter method is a kind of filtering method of removal time series data medium-high frequency variation.The algorithm smooth effect
Good, algorithm is simple, but has ignored the important change information in timing, and interative computation time-consuming is larger.Savitzky-Golay filter
Wave (S-G filtering), is a kind of convolution algorithm based on least square, the algorithm is not by timing time scale, data space and biography
The limitation of sensor is obtained the siding-to-siding block length of fitting of a polynomial order and filter by experience.The siding-to-siding block length of usual filter is got over
Greatly, the result obtained is more smooth-out, and which results in change informations important in time series to be ignored.Fourier transform
Method is to synthesize a complicated timing curve with a series of mode of sinusoidal and cosine wave superpositions, realizes temporal filtering
Process.Fourier transformation obtains curve with apparent vegetation growth variation characteristic, curve smoothing, but because fitting formula is symmetrical
Property it is stringent, asymmetrical information extract the problem of on be difficult to realize, such as rural activity, settlement place construction either other artificial shadows
Variation caused by ringing can not be applicable in completely in remote sensing time series data is filtered and rebuild.Therefore, it is used in the embodiment of the present invention
Default remote sensing temporal model realizes that the NDVI time series data to target area in first object year and the 1 target year is fitted,
And then determine first NDVI timing curve of the target area in first object year and the 2nd NDVI timing in the 1 target year
Curve.The default remote sensing temporal model used in the embodiment of the present invention is the remote sensing constructed in advance by improved least square method
Temporal model, specifically a kind of function model for containing sine term, cosine term, constant term, first order and quadratic term and determining.It adopts
Apparent vegetation growth variation characteristic can be not only fitted with default remote sensing temporal model, and mankind's activity frequently, soil
Cover type region complicated and changeable also available preferable fitting effect.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention, institute
The specific representation for stating default remote sensing temporal model is as follows:
Wherein, a0For constant term, a1、a2、a3And a4Respectively cosine term coefficient, sinusoidal term coefficient, Monomial coefficient and two
Secondary term coefficient, and be constant；X is the time parameter in the first object year or the second target year as unit of day, T
For the number of days in the first object year or the second target year；For NDVI time series data corresponding with x.
Preferably, the x in the embodiment of the present invention is the number of days of Julian date, i.e., the in 1 year how many days.Cosine term
Coefficient a1, sinusoidal term coefficient a2Express year-end drawdown level, Monomial coefficient a3It expresses and changes between year, two-term coefficient a4Expression
It is disturbed between year.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention,
Determine the target area before the NDVI time series data in the first object year and the second target year, further includes:
The target area is carried out in the timing remotely-sensed data in the first object year and the second target year respectively
Pretreatment；
The pretreatment includes at least: radiation calibration, atmospheric correction and geometric correction.
Specifically, in the embodiment of the present invention, since the target area directlyed adopt is in first object year and the second mesh
The timing remotely-sensed data for marking year, the Annual variations testing results inaccuracy that may result in, it is therefore desirable to obtain when
Sequence remotely-sensed data is pre-processed, and pretreatment operation includes at least radiation calibration, atmospheric correction and geometric correction.Wherein radiation mark
Fixed and atmospheric correction is two processes of radiant correction, and radiation calibration makes the remote sensing image picture element brightness value of remotely-sensed data
(Digital Number, DN) is converted to the radiance value on atmosphere top；Atmospheric correction is to eliminate the factors such as atmosphere and illumination
Influence to clutter reflections.Geometric correction is to eliminate the geometric distortion of image, improves its geo-location precision, can match reality
The position of border atural object.Wherein, radiation calibration by the radiation calibration module in ENVI remote sensing image professional treatment software at
Reason realizes that atmospheric correction disturbs adaptive processing system (Landsat Ecosystem by the Landsat ecosystem
Disturbance Adaptive Processing System, LEDAPS) in atmospheric correction tool carry out processing realization, it is several
What correction carries out processing realization by the geometric correction model in ENVI remote sensing image professional treatment software, makes the mistake of geometric position
Difference control is within a pixel.
It should be noted that the pretreatment operation in the embodiment of the present invention can also include that image is cut, it can be for spy
Determine region to be cut, the remote sensing timing image of the target area after obtaining radiation calibration, atmospheric correction and geometric correction.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention, institute
It states based on dynamic time warping DTW algorithm, determines between the first NDVI timing curve and the 2nd NDVI timing curve
DTW distance value, specifically include:
According to the first NDVI timing curve and the 2nd NDVI timing curve, institute is determined by the DTW algorithm
It states the most short consolidation path of the first NDVI timing curve and the 2nd NDVI timing curve and calculates the most short consolidation path
Length value；
The DTW distance value is calculated based on the length value.
Specifically, the Land cover interannual variance detection method provided in the embodiment of the present invention, it is assumed that calculate similarity
Two time serieses be respectively as follows: the first NDVI timing curve be X={ x1,x2,x3…xm, x1、x2、x3…xmRespectively first
M discrete time coordinate on NDVI timing curve, x1For the starting point of the first NDVI timing curve, xmFor the first NDVI timing
The length of the end point of curve, the first NDVI timing curve is | X |；2nd NDVI timing curve is Y={ y1,y2,y3…yn,
y1、y2、y3…ynN discrete time coordinate on respectively the 2nd NDVI timing curve, y1For rising for the 2nd NDVI timing curve
Initial point, ynLength for the end point of the 2nd NDVI timing curve, the 2nd NDVI timing curve is | Y |.The form in consolidation path
It can indicate are as follows: W={ w1,w2,w3,...,wk,...,wK, wherein Max (| X |, | Y |)≤| W |≤| X |+| Y |.wkForm
For wk(i, j), wherein what i was indicated is i-th of time coordinate in X, and what j was indicated is j-th of time coordinate in Y.Consolidation road
Diameter W starts from (1,1), ends up to (m, n), guarantees that each time coordinate in X and Y occurs in W, w in regular pathk's
Subscript is monotonic increase, i.e. wk=(i, j), wk+1=(i ', j '), i≤i '≤i+1, j≤j '≤j+1.Traverse the first NDVI
One most short regular path finally can be obtained in each time coordinate in timing curve and the 2nd NDVI timing curve.It calculates
The length value in most short consolidation path, i.e., the sum of the distance of all route segments in most short consolidation path, further according to most short consolidation path
Length value DTW distance value is calculated.
In the embodiment of the present invention, calculated by the most short consolidation path of determination, and according to the length value in most short consolidation path
DTW distance value provides a kind of implementable solution for the acquisition of DTW distance value.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention, institute
The length value for stating most short consolidation path is calculated by following formula:
D (i, j)=Dist (i, j)+min { D (i-1, j), D (i, j-1), D (i-1, j-1) }
Wherein, i is i-th of time coordinate in the first NDVI timing curve, and j is the 2nd NDVI timing curve
In j-th of time coordinate, i >=1, j >=1, Dist (i, j)=| xi-yj|, and have D (1,1)=Dist (1,1), xiIt is described
I-th of time coordinate value in first NDVI timing curve,yjFor j-th of time coordinate in the 2nd NDVI timing curve
Value.
Specifically, the length value in most short consolidation path is calculated in the embodiment of the present invention by above-mentioned formula, is returned to be most short
The determination of the length value in whole path provides a kind of feasible scheme.It should be noted that calculating most short return in the embodiment of the present invention
The process of the length value in whole path is actually also the process for determining most short consolidation path, and the two complements each other.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention,
After the length value for determining most short consolidation path, need to calculate DTW distance value according to length value, it is true especially by following formula
It is fixed:
DDTW=D/ (m+n)
Wherein, DDTWFor DTW distance value, D is the length value in the most short consolidation path finally determined.
On the basis of the above embodiments, the Land cover interannual variance detection method provided in the embodiment of the present invention, institute
It states and determines the target area in the Annual variations in the first object year and the second target year based on the DTW distance value
Testing result specifically includes:
The DTW distance value is compared with variation detection threshold value；
If judgement knows that the DTW distance value is greater than the variation detection threshold value, it is determined that the Annual variations detection knot
Fruit is to change.
Specifically, since DTW distance value is bigger, the NDVI timing curve similarity for representing for twice is remoter, i.e. twice
Land cover pattern variation it is bigger.Therefore each Land cover types are comprehensively considered in the embodiment of the present invention in the DTW distance in two times
Value, is set as 1.0 for the size for changing detection threshold value, can extract the changing condition of all Land cover types, finally obtain
The Annual variations testing result of land cover pattern.In the embodiment of the present invention after obtaining DTW distance value, by DTW distance value and variation
Detection threshold value is compared, if judgement knows that DTW distance value is greater than variation detection threshold value, it is determined that Annual variations testing result is
It changes, otherwise determines that Annual variations testing result is not change.
Annual variations testing result is determined by introducing variation detection threshold value in the embodiment of the present invention, so that Annual variations
Testing result is more intuitive.
As shown in Fig. 2, providing a kind of change of border of land cover pattern year on the basis of the above embodiments, in the embodiment of the present invention
Change detection system, comprising: NDVI time series data obtains module 21, NDVI timing curve obtains module 22 and Annual variations detection knot
Fruit obtains module 23.Wherein,
NDVI time series data acquisition module 21 is used for the timing based on target area in first object year and the 1 target year
Remotely-sensed data determines the target area in the normalization difference vegetation in the first object year and the second target year respectively
Index NDVI time series data；
NDVI timing curve obtains module 22 and is used for based on default remote sensing temporal model, respectively to the target area in institute
The NDVI time series data for stating first object year and the second target year is fitted, and determines the target area described first
The first NDVI timing curve in target year and the 2nd NDVI timing curve in the second target year；
Annual variations testing result obtains module 23 and is used to be based on dynamic time warping DTW algorithm, determines described first
DTW distance value between NDVI timing curve and the 2nd NDVI timing curve, and based on described in DTW distance value determination
Annual variations testing result of the target area in the first object year and the second target year.
Specifically, NDVI time series data obtains timing remotely-sensed data of the module 21 according to target area in first object year,
Target area can be determined in vegetation index (the Normalized Difference in first object year
Vegetation Index, NDVI) time series data；It similarly, can according to target area in the timing remotely-sensed data in the 1 target year
To determine target area in the NDVI time series data in the 1 target year.NDVI time series data is determined according to timing remotely-sensed data
Implementation can realize by the calculation formula of NDVI, NDVI be equal near infrared band remotely-sensed data and red spectral band it is distant
Feel data difference, divided by the two wave bands remotely-sensed data and, be embodied as with formula:
Wherein, RNIRIndicate the Reflectivity for Growing Season of near infrared band, the i.e. remotely-sensed data of near infrared band, RRedIndicate feux rouges
The Reflectivity for Growing Season of wave band, the i.e. remotely-sensed data of red spectral band.
NDVI time series data obtains module 21 and is determining target area in first object year and the 1 target year respectively
After NDVI time series data, module 22 is obtained by NDVI timing curve and is based on default remote sensing temporal model, is existed respectively to target area
The NDVI time series data in first object year and the 1 target year is fitted, and determines target area the first of first object year
NDVI timing curve and the 2nd NDVI timing curve in the 1 target year.In the embodiment of the present invention when described default remote sensing
Sequence model refers to the continuous function model that can be fitted to discrete NDVI time series data.By presetting remote sensing timing mould
Type is fitted the NDVI time series data in first object year and the 1 target year, determines the NDVI timing in first object year respectively
The parameter value and the NDVI time series data in the 1 target year of the corresponding default remote sensing temporal model of data be corresponding preset it is distant
The parameter value for feeling temporal model, obtains the first NDVI timing curve and the 2nd NDVI timing curve.
After NDVI timing curve acquisition module 22 has determined the first NDVI timing curve and the 2nd NDVI timing curve, by year
Border, which changes testing result and obtains module 23, is based on dynamic time warping (Dynamic Time Warping, DTW) algorithm, determines the
DTW distance value between one NDVI timing curve and the 2nd NDVI timing curve, and determine that target area exists based on DTW distance value
The Annual variations testing result in first object year and the 1 target year.
The effect of each module and the above method in the Land cover interannual variance detection system provided in the embodiment of the present invention
The operating process of class embodiment is correspondingly that details are not described herein in the embodiment of the present invention.
The Land cover interannual variance detection system provided in the embodiment of the present invention, by target area in first object year
With the timing remotely-sensed data in the 1 target year, determine target area in the NDVI timing in first object year and the 1 target year respectively
Data；On the basis of default remote sensing temporal model carries out the NDVI time series data reconstruct of two different years, DTW algorithm is utilized
The similarity system design for carrying out the timing curve of two different years obtains the Annual variations testing result of land cover pattern.The present invention
The use value that timing remotely-sensed data has sufficiently been excavated in embodiment, in terms of using remote sensing timing remote sensing imagery change detection, tool
There is great application prospect.And the displacement problem on time shaft can be overcome by DTW algorithm to a certain extent, it can drop
The influence of the exceptional values such as low noise and cloud, obtains better matching effect, improves the detection accuracy of Land cover interannual variance.
On the basis of the above embodiments, a kind of Land cover interannual variance detection system provided in the embodiment of the present invention
In further include preprocessing module, the preprocessing module is used for: determining the target area in the first object year and institute
Before the NDVI time series data for stating for the 1 target year, respectively to the target area in the first object year and second mesh
The timing remotely-sensed data in mark year is pre-processed；The pretreatment includes at least: radiation calibration, atmospheric correction and geometric correction.
On the basis of the above embodiments, a kind of Land cover interannual variance detection system provided in the embodiment of the present invention
In, it includes DTW distance value computational submodule, the DTW distance value computational submodule that Annual variations testing result, which obtains module 23,
For: according to the first NDVI timing curve and the 2nd NDVI timing curve, described is determined by the DTW algorithm
The most short consolidation path of one NDVI timing curve and the 2nd NDVI timing curve, and calculate the length in the most short consolidation path
Angle value；The DTW distance value is calculated based on the length value.
On the basis of the above embodiments, a kind of Land cover interannual variance detection system provided in the embodiment of the present invention
In, it further includes that Annual variations testing result determines submodule, the Annual variations inspection that Annual variations testing result, which obtains module 23,
It surveys result and determines that submodule is used for: the DTW distance value is compared with variation detection threshold value；If the DTW is known in judgement
Distance value is greater than the variation detection threshold value, it is determined that the Annual variations testing result is to change.
As shown in figure 3, on the basis of the above embodiments, a kind of electronic equipment is additionally provided in the embodiment of the present invention, wrap
It includes: processor (processor) 301, memory (memory) 302, communication interface (Communications Interface)
303 and bus 304；Wherein,
The processor 301, memory 302, communication interface 303 complete mutual communication by bus 304.It is described to deposit
Reservoir 302 is stored with the program instruction that can be executed by the processor 301, and processor 301 is used to call the journey in memory 302
Sequence instruction, to execute method provided by above-mentioned each method embodiment, for example, S1, based on target area in first object
The timing remotely-sensed data in year and the 1 target year, determines the target area in the first object year and second mesh respectively
Mark the vegetation index NDVI time series data in year；S2, based on default remote sensing temporal model, respectively to the target area
NDVI time series data of the domain in the first object year and the second target year is fitted, and determines the target area in institute
State the first NDVI timing curve in first object year and the 2nd NDVI timing curve in the second target year；S3 is based on
Dynamic time warping DTW algorithm, determine DTW between the first NDVI timing curve and the 2nd NDVI timing curve away from
Determine the target area in the year in the first object year and the second target year from value, and based on the DTW distance value
Border changes testing result.
Logical order in memory 302 can be realized by way of SFU software functional unit and as independent product pin
It sells or in use, can store in a computer readable storage medium.Based on this understanding, technical side of the invention
Substantially the part of the part that contributes to existing technology or the technical solution can be with the shape of software product in other words for case
Formula embodies, which is stored in a storage medium, including some instructions are used so that a calculating
Machine equipment (can be personal computer, server or the network equipment etc.) executes each embodiment the method for the present invention
All or part of the steps.And storage medium above-mentioned include: USB flash disk, mobile hard disk, read-only memory (Read-Only Memory,
ROM), random access memory (Random Access Memory, RAM), magnetic or disk etc. are various can store program
The medium of code.
On the basis of the above embodiments, a kind of non-transient computer readable storage medium is additionally provided in the embodiment of the present invention
Matter, the non-transient computer readable storage medium store computer instruction, and the computer instruction executes the computer
To execute method provided by above-mentioned each method embodiment, for example, S1, based on target area in first object year and second
The timing remotely-sensed data in target year determines the target area returning in the first object year and the second target year respectively
One changes difference vegetation index NDVI time series data；S2, based on default remote sensing temporal model, respectively to the target area described
The NDVI time series data in first object year and the second target year is fitted, and determines the target area in first mesh
Mark the first NDVI timing curve in year and the 2nd NDVI timing curve in the second target year；S3 is based on dynamic time
It is bent DTW algorithm, determines the DTW distance value between the first NDVI timing curve and the 2nd NDVI timing curve, and
Determine that the target area is examined in the Annual variations in the first object year and the second target year based on the DTW distance value
Survey result.
The apparatus embodiments described above are merely exemplary, wherein described, unit can as illustrated by the separation member
It is physically separated with being or may not be, component shown as a unit may or may not be physics list
Member, it can it is in one place, or may be distributed over multiple network units.It can be selected according to the actual needs
In some or all of the modules achieve the purpose of the solution of this embodiment.Those of ordinary skill in the art are not paying creativeness
Labour in the case where, it can understand and implement.
Through the above description of the embodiments, those skilled in the art can be understood that each embodiment can
It realizes by means of software and necessary general hardware platform, naturally it is also possible to pass through hardware.Based on this understanding, on
Stating technical solution, substantially the part that contributes to existing technology can be embodied in the form of software products in other words, should
Computer software product may be stored in a computer readable storage medium, such as ROM/RAM, magnetic disk, CD, including several fingers
It enables and using so that a computer equipment (can be personal computer, server or the network equipment etc.) executes each implementation
Method described in certain parts of example or embodiment.
Finally, it should be noted that the above embodiments are merely illustrative of the technical solutions of the present invention, rather than its limitations；Although
Present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that: it still may be used
To modify the technical solutions described in the foregoing embodiments or equivalent replacement of some of the technical features；
And these are modified or replaceed, technical solution of various embodiments of the present invention that it does not separate the essence of the corresponding technical solution spirit and
Range.
Claims (10)
1. a kind of Land cover interannual variance detection method characterized by comprising
Timing remotely-sensed data based on target area in first object year and the 1 target year determines that the target area exists respectively
The vegetation index NDVI time series data in the first object year and the second target year；
Based on default remote sensing temporal model, respectively to the target area in the first object year and the second target year
NDVI time series data is fitted, determine the target area the first object year the first NDVI timing curve and
The 2nd NDVI timing curve in the second target year；
Based on dynamic time warping DTW algorithm, determine the first NDVI timing curve and the 2nd NDVI timing curve it
Between DTW distance value, and determine the target area in the first object year and second mesh based on the DTW distance value
Mark the Annual variations testing result in year.
2. Land cover interannual variance detection method according to claim 1, which is characterized in that the default remote sensing timing
Model is the function model determined by sine term, cosine term, constant term, first order and quadratic term.
3. Land cover interannual variance detection method according to claim 2, which is characterized in that the default remote sensing timing
The specific representation of model is as follows:
Wherein, a0For constant term, a1、a2、a3And a4Respectively cosine term coefficient, sinusoidal term coefficient, Monomial coefficient and quadratic term
Coefficient, and be constant；X is the time parameter in the first object year or the second target year as unit of day, and T is institute
State the number of days in first object year or the second target year；For NDVI time series data corresponding with x.
4. Land cover interannual variance detection method according to claim 1, which is characterized in that determining the target area
Domain is before the NDVI time series data in the first object year and the second target year, further includes:
The timing remotely-sensed data to the target area in the first object year and the second target year is located in advance respectively
Reason；
The pretreatment includes at least: radiation calibration, atmospheric correction and geometric correction.
5. Land cover interannual variance detection method described in any one of -4 according to claim 1, which is characterized in that the base
In dynamic time warping DTW algorithm, the DTW between the first NDVI timing curve and the 2nd NDVI timing curve is determined
Distance value specifically includes:
According to the first NDVI timing curve and the 2nd NDVI timing curve, described is determined by the DTW algorithm
The most short consolidation path of one NDVI timing curve and the 2nd NDVI timing curve, and calculate the length in the most short consolidation path
Angle value；
The DTW distance value is calculated based on the length value.
6. Land cover interannual variance detection method according to claim 5, which is characterized in that the most short consolidation path
Length value calculated by following formula:
D (i, j)=Dist (i, j)+min { D (i-1, j), D (i, j-1), D (i-1, j-1) }
Wherein, i is i-th of time coordinate in the first NDVI timing curve, and j is in the 2nd NDVI timing curve
J-th of time coordinate, i >=1, j >=1, Dist (i, j)=| xi-yj|, and have D (1,1)=Dist (1,1), xiIt is described first
I-th of time coordinate value in NDVI timing curve, yjFor j-th of time coordinate value in the 2nd NDVI timing curve.
7. Land cover interannual variance detection method described in any one of -4 according to claim 1, which is characterized in that the base
Determine that the target area is detected in the Annual variations in the first object year and the second target year in the DTW distance value
As a result, specifically including:
The DTW distance value is compared with variation detection threshold value；
If judgement knows that the DTW distance value is greater than the variation detection threshold value, it is determined that the Annual variations testing result is
It changes.
8. a kind of Land cover interannual variance detection system characterized by comprising
NDVI time series data obtain module, for based on target area first object year and the 1 target year timing remote sensing number
According to determining the target area in the vegetation index in the first object year and the second target year respectively
NDVI time series data；
NDVI timing curve obtains module, for based on default remote sensing temporal model, respectively to the target area described the
The NDVI time series data in 1 target year and the second target year is fitted, and determines the target area in the first object
Year the first NDVI timing curve and the 2nd NDVI timing curve in the second target year；
Annual variations testing result obtains module, for being based on dynamic time warping DTW algorithm, determines the first NDVI timing
DTW distance value between curve and the 2nd NDVI timing curve, and the target area is determined based on the DTW distance value
In the Annual variations testing result in the first object year and the second target year.
9. a kind of electronic equipment characterized by comprising
At least one processor, at least one processor, communication interface and bus；Wherein,
The processor, memory, communication interface complete mutual communication by the bus；
The memory is stored with the program instruction that can be executed by the processor, and the processor calls described program instruction,
To execute such as Land cover interannual variance detection method of any of claims 1-7.
10. a kind of non-transient computer readable storage medium, which is characterized in that the non-transient computer readable storage medium is deposited
Computer instruction is stored up, the computer instruction executes the computer as soil of any of claims 1-7 covers
Lid Annual variations detection method.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811222022.7A CN109446962B (en) | 2018-10-19 | 2018-10-19 | Land cover annual change detection method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811222022.7A CN109446962B (en) | 2018-10-19 | 2018-10-19 | Land cover annual change detection method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109446962A true CN109446962A (en) | 2019-03-08 |
| CN109446962B CN109446962B (en) | 2020-12-01 |
Family
ID=65547002
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811222022.7A Expired - Fee Related CN109446962B (en) | 2018-10-19 | 2018-10-19 | Land cover annual change detection method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109446962B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111582554A (en) * | 2020-04-17 | 2020-08-25 | 北京农业信息技术研究中心 | Crop growth prediction method and system |
| CN111721711A (en) * | 2020-06-17 | 2020-09-29 | 国家基础地理信息中心 | Method and device for detecting farmland change, electronic equipment and storage medium |
| CN112507858A (en) * | 2020-12-03 | 2021-03-16 | 中国科学院东北地理与农业生态研究所 | Ecological system attribute component composition structure description method based on remote sensing vegetation index |
| CN116912698A (en) * | 2023-09-04 | 2023-10-20 | 中国地质大学（武汉） | Land coverage change monitoring method and device, electronic equipment and storage medium |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102176242A (en) * | 2011-02-01 | 2011-09-07 | 环境保护部卫星环境应用中心 | Method for removing cloud noise effects in normalized difference vegetation index (NDVI) time sequence image |
| CN105512616A (en) * | 2015-11-27 | 2016-04-20 | 吉林大学 | Method for rebuilding remote sensing monitoring data set by eliminating atmospheric noise effect |
| CN106096630A (en) * | 2016-06-03 | 2016-11-09 | 福州大学 | A kind of integrated pixel and the Oryza sativa L. method for detecting automatic variation of object |
| CN106897707A (en) * | 2017-03-02 | 2017-06-27 | 苏州中科天启遥感科技有限公司 | Characteristic image time series synthetic method and device based in multi-source points |
| CN107463775A (en) * | 2017-07-24 | 2017-12-12 | 福州大学 | Vegetation based on more Indices variation tendencies is lost in whereabouts recognition methods |
| US20180293671A1 (en) * | 2015-03-27 | 2018-10-11 | Omniearth, Inc. | System and method for predicting crop yield |
-
2018
- 2018-10-19 CN CN201811222022.7A patent/CN109446962B/en not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102176242A (en) * | 2011-02-01 | 2011-09-07 | 环境保护部卫星环境应用中心 | Method for removing cloud noise effects in normalized difference vegetation index (NDVI) time sequence image |
| US20180293671A1 (en) * | 2015-03-27 | 2018-10-11 | Omniearth, Inc. | System and method for predicting crop yield |
| CN105512616A (en) * | 2015-11-27 | 2016-04-20 | 吉林大学 | Method for rebuilding remote sensing monitoring data set by eliminating atmospheric noise effect |
| CN106096630A (en) * | 2016-06-03 | 2016-11-09 | 福州大学 | A kind of integrated pixel and the Oryza sativa L. method for detecting automatic variation of object |
| CN106897707A (en) * | 2017-03-02 | 2017-06-27 | 苏州中科天启遥感科技有限公司 | Characteristic image time series synthetic method and device based in multi-source points |
| CN107463775A (en) * | 2017-07-24 | 2017-12-12 | 福州大学 | Vegetation based on more Indices variation tendencies is lost in whereabouts recognition methods |
Non-Patent Citations (2)
| Title |
|---|
| 韩晓勇等: "基于DTW距离的时序NDVI数据植被信息提取――以秦巴山区为例 ", 《测绘工程》 * |
| 韩晓勇等: "基于DTW距离的时序NDVI数据植被信息提取――以秦巴山区为例", 《测绘工程》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111582554A (en) * | 2020-04-17 | 2020-08-25 | 北京农业信息技术研究中心 | Crop growth prediction method and system |
| CN111582554B (en) * | 2020-04-17 | 2024-02-23 | 北京市农林科学院信息技术研究中心 | Crop growth prediction method and system |
| CN111721711A (en) * | 2020-06-17 | 2020-09-29 | 国家基础地理信息中心 | Method and device for detecting farmland change, electronic equipment and storage medium |
| CN112507858A (en) * | 2020-12-03 | 2021-03-16 | 中国科学院东北地理与农业生态研究所 | Ecological system attribute component composition structure description method based on remote sensing vegetation index |
| CN112507858B (en) * | 2020-12-03 | 2023-03-24 | 中国科学院东北地理与农业生态研究所 | Ecological system attribute component composition structure description method based on remote sensing vegetation index |
| CN116912698A (en) * | 2023-09-04 | 2023-10-20 | 中国地质大学（武汉） | Land coverage change monitoring method and device, electronic equipment and storage medium |
| CN116912698B (en) * | 2023-09-04 | 2023-12-08 | 中国地质大学（武汉） | Land coverage change monitoring method and device, electronic equipment and storage medium |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109446962B (en) | 2020-12-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Yun et al. | Individual tree crown segmentation from airborne LiDAR data using a novel Gaussian filter and energy function minimization-based approach | |
| de Conto et al. | Performance of stem denoising and stem modelling algorithms on single tree point clouds from terrestrial laser scanning | |
| US11521380B2 (en) | Shadow and cloud masking for remote sensing images in agriculture applications using a multilayer perceptron | |
| CN109446962A (en) | Land cover interannual variance detection method and system | |
| Rizeei et al. | Oil palm counting and age estimation from WorldView‐3 imagery and LiDAR data using an integrated OBIA height model and regression analysis | |
| Chrysafis et al. | Evaluation of single-date and multi-seasonal spatial and spectral information of Sentinel-2 imagery to assess growing stock volume of a Mediterranean forest | |
| Halme et al. | Utility of hyperspectral compared to multispectral remote sensing data in estimating forest biomass and structure variables in Finnish boreal forest | |
| Feng et al. | Robust non-local fuzzy c-means algorithm with edge preservation for SAR image segmentation | |
| CN109146948A (en) | The quantization of crop growing state phenotypic parameter and the correlation with yield analysis method of view-based access control model | |
| Regniers et al. | Supervised classification of very high resolution optical images using wavelet-based textural features | |
| Dao et al. | Improving hyperspectral image segmentation by applying inverse noise weighting and outlier removal for optimal scale selection | |
| US11238307B1 (en) | System for performing change detection within a 3D geospatial model based upon semantic change detection using deep learning and related methods | |
| Wang et al. | Using GF-2 imagery and the conditional random field model for urban forest cover mapping | |
| US11587249B2 (en) | Artificial intelligence (AI) system and methods for generating estimated height maps from electro-optic imagery | |
| Ok et al. | 2-D delineation of individual citrus trees from UAV-based dense photogrammetric surface models | |
| US11747468B2 (en) | System using a priori terrain height data for interferometric synthetic aperture radar (IFSAR) phase disambiguation and related methods | |
| US11636649B2 (en) | Geospatial modeling system providing 3D geospatial model update based upon predictively registered image and related methods | |
| Xing et al. | A novel change detection method using remotely sensed image time series value and shape based dynamic time warping | |
| CN115344815A (en) | Natural runoff change attribution analysis method and system considering vegetation space change | |
| Lei et al. | A novel algorithm of individual tree crowns segmentation considering three-dimensional canopy attributes using UAV oblique photos | |
| Puletti et al. | Enhancing wall-to-wall forest structure mapping through detailed co-registration of airborne and terrestrial laser scanning data in Mediterranean forests | |
| Okubo et al. | Land use/cover classification of a complex agricultural landscape using single-dated very high spatial resolution satellite-sensed imagery | |
| Lasaponara et al. | Pattern recognition and classification using VHR data for archaeological research | |
| US11816793B2 (en) | Geospatial modeling system providing 3D geospatial model update based upon iterative predictive image registration and related methods | |
| Ez-zahouani et al. | Determination of Segmentation Parameters for Object-Based Remote Sensing Image Analysis from Conventional to Recent Approaches: A Review. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201201 Termination date: 20211019 |