CN107704807A - A kind of dynamic monitoring method based on multi-source remote sensing sequential images - Google Patents
A kind of dynamic monitoring method based on multi-source remote sensing sequential images Download PDFInfo
- Publication number
- CN107704807A CN107704807A CN201710788089.6A CN201710788089A CN107704807A CN 107704807 A CN107704807 A CN 107704807A CN 201710788089 A CN201710788089 A CN 201710788089A CN 107704807 A CN107704807 A CN 107704807A
- Authority
- CN
- China
- Prior art keywords
- time series
- remote sensing
- trend
- sensing images
- time
- 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
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/25—Fusion techniques
- G06F18/251—Fusion techniques of input or preprocessed data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
- G06V10/00—Arrangements for image or video recognition or understanding
- G06V10/20—Image preprocessing
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Multimedia (AREA)
- Astronomy & Astrophysics (AREA)
- Remote Sensing (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
The present invention relates to the dynamic monitoring method suitable for multi-source time series remote sensing images, this method effectively can improve fusion method using the spatial information of data, reduce the feature difference between multi-source Remote Sensing Images, improve the effect of dynamic monitoring.First, multi-source Remote Sensing Images sequence is built, sequence remote sensing images are pre-processed using appropriate image pre-processing method, make remote sensing images sequence has the features such as coherent image radiation characteristic, higher image registration accuracy；Next, the enhancing multi-source Remote-sensing Image Fusion eMulTiFuse based on spatial and temporal distributions merges to remote sensing images sequence, make time series that there are increasingly similar space characteristics；Finally, based on pixel feature, dynamic monitoring is carried out to time series remote sensing images using Mann Kendall trend-monitorings methods, obtains the dynamic monitor result.The present invention has the higher degree of accuracy.
Description
First, technical field
The present invention relates to the dynamic monitoring method suitable for multi-source time series remote sensing images, belongs to Spatial Information Technology neck
Domain.
2nd, background technology
At present, dynamic monitoring based on time series remote sensing images is in, low resolution (especially using MODIS time serieses as
It is most) in the majority, and wherein especially based on the dynamic monitoring of single source time sequence remote sensing images.Single source time sequence remote sensing images or more
Or it is few exist without contaminated regions such as value, bad value, more clouds and mists, especially in the area of warm and moist (such as Chongqing, Sichuan, expensive
The ground such as state, Yunnan) interference of cloud and mist etc. is frequently subjected to, the simple remote sensing images obtained by a satellite are often difficult to satisfaction pair
The dynamic monitoring in this place.There are some researches show be longer than 2 years, the accuracy of dynamic monitoring is bright when without cloud and mist view data time interval
It is aobvious to reduce.Multi-source data can supplement missing and contaminated data in single source images in time, greatly improve the knot of dynamic monitoring
Fruit.Remotely-sensed data image quality and processing procedure difference from different sensors are larger, it is impossible to it is mixed to be used directly to generation
The research that time series carries out dynamic monitoring is closed, therefore it is particularly heavy to study the dynamic monitoring method based on multi-source Remote Sensing Image Fusion
Will.
With the variation in remotely-sensed data source and going deep into for remote sensing technology research, the dynamic monitoring based on remote sensing is special
It is not because gradually being used by people with many merits, therefore more and more based on the dynamic monitoring of Multitemporal Remote Sensing Images
People achieves great successes to having carried out extensive research based on this method.Existing Multitemporal Remote Sensing Images dynamic is supervised
Survey method can not only study on monitoring head of district's period variation tendency, and be able to detect that catastrophe and the crop of generation
Deng real-time growing way dynamic etc., multi-source time series remote sensing images, which are introduced into the technique study of dynamic monitoring, has bright prospects
And significance, but it may also make because solar illumination angle, atmospheric conditions, surface humidity, sensor accuracy etc. change
Monitoring result misalignment.Based on this, the present invention have studied a kind of multi-source time series remote sensing image fusion side of combination space-time characteristic
Method, this method effectively can improve fusion method using the spatial information of data, and the feature reduced between multi-source Remote Sensing Images is poor
It is different, and then improve the effect of dynamic monitoring.
3rd, the content of the invention
1st, purpose：Disturb etc. and ask to caused by time series analysis herein for cloud layer in single source time sequence remote sensing images
Topic, proposes a kind of multi-source time series remote sensing image fusion method of combination space-time characteristic, and this method can effectively utilize data
Spatial information improves fusion method, reduces the feature difference between multi-source Remote Sensing Images, and then improve the effect of dynamic monitoring.
2nd, technical scheme：
Suitable for the dynamic monitoring method of multi-source time series remote sensing images, it is characterised in that comprise the following steps (as schemed
1)：Step 1：The structure of multi-source Remote Sensing Images sequence
The feature that the non-pending remote sensing images of remotely-sensed data feature referred to herein have in itself, and refer to possess the progress time
Some features that the remote sensing images of sequence satellite image analysis should have：
First, time-series image radiation characteristic is consistent.Multitemporal Remote Sensing Images are often due to external factor (such as light
According to, atmospheric radiation, sensor platform etc.) influence and cause atural object that " puppet change " occurs, therefore image radiation characteristic is unanimously
Critically important.
Secondly, image registration accuracy will height.Analysis to time series remote sensing images at present is based on middle low-resolution image
More, the analysis method of use is mostly the analysis method based on pixel, this precision to geometric registration of imagery propose it is higher will
Ask.
Finally, image picture elements value will be the reflection of atural object real features.By image bad pixel, image bad track and complicated gas
Remote sensing images pixel caused by pixel distortion etc. caused by waiting phenomenon (such as rain, cloud, mist or snow etc.) can not reflect truly earth's surface
State, the presence of this feature make the accuracy of sequential remote Sensing Image Analysis.
Therefore, on the basis of suitable remote sensing images are chosen, remote sensing images are entered using appropriate image pre-processing method
Row pretreatment, makes remote sensing images sequence meet that features described above is vital.
Further, since remote sensing images are influenceed in imaging by factors such as sensor, flying platform, atmospheric radiations, make
Distortion in terms of geometry and radiation to some extent be present in image, influence picture quality, the use to image brings certain
It is difficult.Using associated picture preprocess method, geometry, radiation error are reduced and eliminated as far as possible, is strengthened remote sensing images,
Important leverage is provided for remote Sensing Image Analysis, interpretation, understanding, identification and information extraction.
Step 2：Multi-source Remote Sensing Image Fusion
(1) multi-source Remote-sensing Image Fusion-MulTiFuse based on pixel
When being merged using MulTiFuse fusion methods to two time series remote sensing images from different sensors,
First, the weight temporal sequence relation of two time serieses is derived, in order to maximize the statistical nature of this relation and consideration
Abnormal conditions, weight optimization is associated to it before regression analysis is weighted to two time serieses；Then basis
Optimum regression model merges to multi-source Remote Sensing Images sequence.
The acquisition time of discrete-time series remote sensing images X time series pixel value is set as TX, time series X length
Spend for kx；The acquisition time of discrete-time series remote sensing images Y time series pixel value is similarly set as TY, time series Y's
Length is ky.Assuming that above-mentioned two discrete-time series data length is inconsistent, image data acquiring time interval is different, and
Two time series acquisition time scopes have lap.Be primarily based on linear interpolation method to X time series Y acquisition time
Row interpolation is clicked through, obtains interpolation time sequence Xint, row interpolation is similarly entered to time series Y and obtains interpolation time sequence Yint, obtain
To the time series as shown in following formula (1) and (2)：
Wherein, XlinAnd YlinX the and Y time serieses formed by interpolation are represented respectively.
Singular value may be obtained by above-mentioned simple linear interpolation., particularly in acquisition time interval longer time
Sequence, influenceing can be very big.In order to solve this problem, the observation difference (i.e. amplitude) of two time is introduced, its absolute value is got over
Greatly, weight is returned with regard to smaller；On the contrary, the absolute value of amplitude is smaller, it is bigger to return weight.
Next, definition is interpolated amplitude and normalizes the amplitude.First, it is defined on the timeI=1 ..., kxIt is exhausted
To amplitude mY, such as following formula (3)：
Wherein,WithRepresent that time series Y originates in time series X time respectivelyFront and rear.WithRespectively away from the X timesNearest Y observation.
To amplitude mYIt is normalized, sees below formula (4) and (5)：
Similar, define mX、mX lWithBecause overlapping period X and Y length of time series is different, it is therefore desirable to by the time
Sequence length is taken into account, rightWithIt is normalized again, such as following formula (6) and (7)：
Wherein, lxRefer to time tX iQuantity, lYRefer to time tY iQuantity.Return weightIt is defined as follows formula
(8)：
Exponential weighting function ewf is introduced herein to adjust weight, such as following formula (9)：
Wherein, ewf=0,1,2, no weight, a weight, square weight are indicated respectively.Interpolation time sequence XintAnd Yint
Fit correlation by Seber's and Lee《Linear regression analysis(2nded.)》Weighting in one book is minimum
Two, which multiply the fitting Return Law, is calculated, and linear fit model is obtained by formula (10)：
Yint=a+bXint (10)
Wherein, a and b is constant.
(2) enhancing multi-source Remote-sensing Image Fusion-eMulTiFuse based on spatial and temporal distributions
Above-mentioned algorithm provides good thinking to multisource data fusion from the angle of temporal characteristics, and we are on this basis
Carry out unified, selection and amendment to unusual pixel on space characteristics to similar pixel again, strengthen multi-source Remote Sensing Images number
According to syncretizing effect.
This research carries out unsupervised classification to time series remote sensing images, remote sensing image data is clustered, as right
The foundation that unusual pixel is modified.Carrying out the conventional algorithm of cluster analysis has ISODATA (iteration self-organizing data analysis), K-
Means, chain method etc..ISODATA cluster centre is determined that this point is calculated with K-Means by the interative computation of this average
Method is similar, but the algorithm can draw the experience obtained by intermediate result, has self-organization, it can be said that compared with K-Means
In principle further, it is the comparatively ideal algorithm of classifying quality, this research is using this method progress cluster analysis.
According to the result of multispectral image unsupervised classification, the interpolation image that time prediction obtains is corresponded to the image and is carried out
Classification, then judge whether interpolation image pixel belongs to such with mahalanobis distance (Mahalanobis distance), geneva
Distance algorithm such as following formula (11)：
Wherein, XiPixel to be discriminated is represented,Refer to classification r average pixel value,Reflect pixel i and generic
The diversity of remaining pixel.
Because mahalanobis distance assumes that data are multi-source normalities in itself, therefore the approximate card side that obeys of above-mentioned mahalanobis distance value divides
Cloth.We select the quantile of chi square distribution 90% criterion whether abnormal as pixel value is judged herein, more than the quantile value
Pixel value exception must be thought, average is obtained by the use of its generic pixel value and replace the pixel value and be used as correction result.
Step 3：Dynamic monitoring based on time series remote sensing images
The present invention carries out the dynamic monitoring of sequence remote sensing image, Mann- using Mann-Kendall trend tests method
Which type of regularity of distribution is Kendall trend test methods known pixel value need not have in advance without input parameter value,
The interference of noise can be reduced, convenience of calculation, can be used for the hypothesis testing of variable trends.Mann-Kendall trend tests are general
Based on pixel or single grid, by assuming that having inspected the variation tendency monitoring of atural object index.
In Mann-Kendall trend tests, make the following assumptions：Assume H in source0For：Time series data (X to be detected1,
X2,...,Xn), meet n independent, stochastic variables with being distributed, i.e. random distribution, without notable trend；Alternative hypothesis H1For：It is right
In all k, j≤n, and k ≠ j, XkAnd XjDistribution differ, sequence have enhancing or attenuation trend.Test statistics S is such as
Following formula (14) calculates：
Work as n>When 10, S is approximate to obey standardized normal distribution, and its average is 0, and variance can be using approximate representation as Var(S)=n
(n-1)(2n+5)/18.The normally distributed variable Z of standard is test statistics, for trend test, is calculated by following formula (16)：
In bilateral trend test, in given α confidence levels, if-Z1-α/2≤Z≤Z1-α/2, then it is assumed that null hypothesis
H0 is reliable, and receiving should be it is assumed that illustrates time series without significant change trend.If | Z | ＞ Z1-α/2, then it is assumed that assuming that H1 is correct,
Illustrate that sequence has obvious enhancing or attenuation trend.For test statistics Z, Z>0 shows that time series has enhancing trend, Z<
0 shows that time series has attenuation trend.Confidence level be 90% significance test, Z0.95Equal to 1.28；Confidence level is 95%
Significance test when, Z0.975Equal to 1.64；Confidence level be 99% significance test, Z0.9995Equal to 2.32.
Work as n<When 10, bilateral trend test is carried out based on test statistics S.When determining confidence level α, significance test
In, if | S |≤Sα/2, then null hypothesis is received, the time series is without significant change trend；On the contrary, refusal null hypothesis, receives standby
Select it is assumed that the time series is in significantly to rise or fall trend.For test statistics S, if S>0, show that time series has
It is on the rise, if S=0, no visible trend, if S<0, show that time series has downward trend.
The intermediate value of the ratio of difference and time difference of the Theil Sen slopes based on different year gray value or exponential quantity and
, noise disturbance can be preferably avoided, this is worth calculating no significance test.The computational methods of the value such as following formula (17)：
Wherein, 1<j<i<N (N represents the length of time series).Work as Slope>When 0, gray scale or index value sequence present liter
Trend, work as Slope<When 0, the sequence is on a declining curve.The value of Theil Sen slopes is bigger, then illustrates that the degree of its change is got over
Greatly, it is on the contrary then smaller.
The trend analysis of vegetation coverage is carried out using Mann-Kendall trend tests method, with 95% confidence during inspection
Degree, that is, take level of signifiance α=0.05, work as Slope>0 and | Z | when≤1.96, time series in rise but unconspicuous trend, can
It is in less ascendant trend to think vegetation coverage；Work as Slope>0 and | Z | during ＞ 1.96, time series becomes in obvious rise
Gesture, it is believed that vegetation coverage is in obvious ascendant trend；Work as Slope<0 and | Z | when≤1.96, time series in rise but unobvious
Trend, it is believed that vegetation coverage is in less downward trend；Work as Slope<0 and | Z | during ＞ 1.96, time series is in obvious
Downward trend, it is believed that vegetation coverage declines obvious.
3rd, advantage and effect：The present invention is for cloud layer in single source time sequence remote sensing images to caused by time series analysis
The problems such as interference, propose a kind of multi-source time series remote sensing image fusion method of combination space-time characteristic, first, multi-source remote sensing figure
As the structure of sequence, sequence remote sensing images are pre-processed using appropriate image pre-processing method, make remote sensing images sequence
There is the feature such as coherent image radiation characteristic, higher image registration accuracy；Next, using the enhancing based on spatial and temporal distributions
Multi-source Remote-sensing Image Fusion eMulTiFuse merges to remote sensing images sequence, has time series increasingly similar
Space characteristics；Finally, based on pixel feature, action is entered to time series remote sensing images using Mann-Kendall trend-monitorings method
State monitors, and obtains the dynamic monitor result.The present invention can improve the precision of dynamic monitoring.
4th, illustrate
Regional dynamics monitoring route maps of the Fig. 1 based on multi-source remote sensing sequential images
5th, embodiment
The present invention relates to the dynamic monitoring method suitable for multi-source time series remote sensing images, it is characterised in that including such as
Lower step (such as Fig. 1)：
Step 1：The structure of multi-source Remote Sensing Images sequence
The feature that the non-pending remote sensing images of remotely-sensed data feature referred to herein have in itself, and refer to possess the progress time
Some features that the remote sensing images of sequence satellite image analysis should have：
First, time-series image radiation characteristic is consistent.Multitemporal Remote Sensing Images are often due to external factor (such as light
According to, atmospheric radiation, sensor platform etc.) influence and cause atural object that " puppet change " occurs, therefore image radiation characteristic is unanimously
Critically important.
Secondly, image registration accuracy will height.Analysis to time series remote sensing images at present is based on middle low-resolution image
More, the analysis method of use is mostly the analysis method based on pixel, this precision to geometric registration of imagery propose it is higher will
Ask.
Finally, image picture elements value will be the reflection of atural object real features.By image bad pixel, image bad track and complicated gas
Remote sensing images pixel caused by pixel distortion etc. caused by waiting phenomenon (such as rain, cloud, mist or snow etc.) can not reflect truly earth's surface
State, the presence of this feature make the accuracy of sequential remote Sensing Image Analysis.
Therefore, on the basis of suitable remote sensing images are chosen, remote sensing images are entered using appropriate image pre-processing method
Row pretreatment, makes remote sensing images sequence meet that features described above is vital.
Further, since remote sensing images are influenceed in imaging by factors such as sensor, flying platform, atmospheric radiations, make
Distortion in terms of geometry and radiation to some extent be present in image, influence picture quality, the use to image brings certain
It is difficult.Using associated picture preprocess method, geometry, radiation error are reduced and eliminated as far as possible, is strengthened remote sensing images,
Important leverage is provided for remote Sensing Image Analysis, interpretation, understanding, identification and information extraction.
Step 2：Multi-source Remote Sensing Image Fusion
(1) multi-source Remote-sensing Image Fusion-MulTiFuse based on pixel
When being merged using MulTiFuse fusion methods to two time series remote sensing images from different sensors,
First, the weight temporal sequence relation of two time serieses is derived, in order to maximize the statistical nature of this relation and consideration
Abnormal conditions, weight optimization is associated to it before regression analysis is weighted to two time serieses；Then basis
Optimum regression model merges to multi-source Remote Sensing Images sequence.
The acquisition time of discrete-time series remote sensing images X time series pixel value is set as TX, time series X length
Spend for kx；The acquisition time of discrete-time series remote sensing images Y time series pixel value is similarly set as TY, time series Y's
Length is ky.Assuming that above-mentioned two discrete-time series data length is inconsistent, image data acquiring time interval is different, and
Two time series acquisition time scopes have lap.Be primarily based on linear interpolation method to X time series Y acquisition time
Row interpolation is clicked through, obtains interpolation time sequence Xint, row interpolation is similarly entered to time series Y and obtains interpolation time sequence Yint, obtain
To the time series as shown in following formula (1) and (2)：
Wherein, XlinAnd YlinX the and Y time serieses formed by interpolation are represented respectively.
Singular value may be obtained by above-mentioned simple linear interpolation., particularly in acquisition time interval longer time
Sequence, influenceing can be very big.In order to solve this problem, the observation difference (i.e. amplitude) of two time is introduced, its absolute value is got over
Greatly, weight is returned with regard to smaller；On the contrary, the absolute value of amplitude is smaller, it is bigger to return weight.
Next, definition is interpolated amplitude and normalizes the amplitude.First, it is defined on the timeI=1 ..., kxIt is exhausted
To amplitude mY, such as following formula (3)：
Wherein,WithRepresent that time series Y originates in time series X time respectivelyFront and rear.WithRespectively away from the X timesNearest Y observation.
To amplitude mYIt is normalized, sees below formula (4) and (5)：
Similar, define mX、mX lWithBecause overlapping period X and Y length of time series is different, it is therefore desirable to by the time
Sequence length is taken into account, rightWithIt is normalized again, such as following formula (6) and (7)：
Wherein, lxRefer to time tX iQuantity, lYRefer to time tY iQuantity.Return weightIt is defined as follows formula
(8)：
Exponential weighting function ewf is introduced herein to adjust weight, such as following formula (9)：
Wherein, ewf=0,1,2, no weight, a weight, square weight are indicated respectively.Interpolation time sequence XintAnd Yint
Fit correlation by Seber's and Lee《Linear regression analysis(2nded.)》Weighting in one book is minimum
Two, which multiply the fitting Return Law, is calculated, and linear fit model is obtained by formula (10)：
Yint=a+bXint (10)
Wherein, a and b is constant.
(2) enhancing multi-source Remote-sensing Image Fusion-eMulTiFuse based on spatial and temporal distributions
Above-mentioned algorithm provides good thinking to multisource data fusion from the angle of temporal characteristics, and we are on this basis
Carry out unified, selection and amendment to unusual pixel on space characteristics to similar pixel again, strengthen multi-source Remote Sensing Images number
According to syncretizing effect.
This research carries out unsupervised classification to time series remote sensing images, remote sensing image data is clustered, as right
The foundation that unusual pixel is modified.Carrying out the conventional algorithm of cluster analysis has ISODATA (iteration self-organizing data analysis), K-
Means, chain method etc..ISODATA cluster centre is determined that this point is calculated with K-Means by the interative computation of this average
Method is similar, but the algorithm can draw the experience obtained by intermediate result, has self-organization, it can be said that compared with K-Means
In principle further, it is the comparatively ideal algorithm of classifying quality, this research is using this method progress cluster analysis.
According to the result of multispectral image unsupervised classification, the interpolation image that time prediction obtains is corresponded to the image and is carried out
Classification, then judge whether interpolation image pixel belongs to such with mahalanobis distance (Mahalanobis distance), geneva
Distance algorithm such as following formula (11)：
Wherein, XiPixel to be discriminated is represented,Refer to classification r average pixel value,Reflect pixel i and generic
The diversity of remaining pixel.
Because mahalanobis distance assumes that data are multi-source normalities in itself, therefore the approximate card side that obeys of above-mentioned mahalanobis distance value divides
Cloth.We select the quantile of chi square distribution 90% criterion whether abnormal as pixel value is judged herein, more than the quantile value
Pixel value exception must be thought, average is obtained by the use of its generic pixel value and replace the pixel value and be used as correction result.
Step 3：Dynamic monitoring based on time series remote sensing images
The present invention carries out the dynamic monitoring of sequence remote sensing image, Mann- using Mann-Kendall trend tests method
Which type of regularity of distribution is Kendall trend test methods known pixel value need not have in advance without input parameter value,
The interference of noise can be reduced, convenience of calculation, can be used for the hypothesis testing of variable trends.Mann-Kendall trend tests are general
Based on pixel or single grid, by assuming that having inspected the variation tendency monitoring of atural object index.
In Mann-Kendall trend tests, make the following assumptions：Assume H in source0For：Time series data (X to be detected1,
X2,...,Xn), meet n independent, stochastic variables with being distributed, i.e. random distribution, without notable trend；Alternative hypothesis H1For：It is right
In all k, j≤n, and k ≠ j, XkAnd XjDistribution differ, sequence have enhancing or attenuation trend.Test statistics S is such as
Following formula (14) calculates：
Work as n>When 10, S is approximate to obey standardized normal distribution, and its average is 0, and variance can be using approximate representation as Var(S)=n
(n-1)(2n+5)/18.The normally distributed variable Z of standard is test statistics, for trend test, is calculated by following formula (16)：
In bilateral trend test, in given α confidence levels, if-Z1-α/2≤Z≤Z1-α/2, then it is assumed that null hypothesis
H0Reliably, receiving should be it is assumed that illustrates time series without significant change trend.If | Z | ＞ Z1-α/2, then it is assumed that assuming that H1 is correct,
Illustrate that sequence has obvious enhancing or attenuation trend.For test statistics Z, Z>0 shows that time series has enhancing trend, Z<
0 shows that time series has attenuation trend.Confidence level be 90% significance test, Z0.95Equal to 1.28；Confidence level is 95%
Significance test when, Z0.975Equal to 1.64；Confidence level be 99% significance test, Z0.9995Equal to 2.32.
Work as n<When 10, bilateral trend test is carried out based on test statistics S.When determining confidence level α, significance test
In, if | S |≤Sα/2, then null hypothesis is received, the time series is without significant change trend；On the contrary, refusal null hypothesis, receives standby
Select it is assumed that the time series is in significantly to rise or fall trend.For test statistics S, if S>0, show that time series has
It is on the rise, if S=0, no visible trend, if S<0, show that time series has downward trend.
The intermediate value of the ratio of difference and time difference of the Theil Sen slopes based on different year gray value or exponential quantity and
, noise disturbance can be preferably avoided, this is worth calculating no significance test.The computational methods of the value such as following formula (17)：
Wherein, 1<j<i<N (N represents the length of time series).Work as Slope>When 0, gray scale or index value sequence present liter
Trend, work as Slope<When 0, the sequence is on a declining curve.The value of Theil Sen slopes is bigger, then illustrates that the degree of its change is got over
Greatly, it is on the contrary then smaller.
The trend analysis of vegetation coverage is carried out using Mann-Kendall trend tests method, with 95% confidence during inspection
Degree, that is, take level of signifiance α=0.05, work as Slope>0 and | Z | when≤1.96, time series in rise but unconspicuous trend, can
It is in less ascendant trend to think vegetation coverage；Work as Slope>0 and | Z | during ＞ 1.96, time series becomes in obvious rise
Gesture, it is believed that vegetation coverage is in obvious ascendant trend；Work as Slope<0 and | Z | when≤1.96, time series in rise but unobvious
Trend, it is believed that vegetation coverage is in less downward trend；Work as Slope<0 and | Z | during ＞ 1.96, time series is in obvious
Downward trend, it is believed that vegetation coverage declines obvious.
Embodiment 1：
Remote sensing images are carried out with radiant correction, relative radiometric normalization, geometrical registration and correction, image cloud removing, figure
Slice, which is taken, to be removed and image mosaic pre-processes with the series such as cutting, and obtains spectrum consistent, high registration accuracy, the true atural object of reflection
Spectral value.To differing greatly from different sensors spectral value, can not handle to obtain desired light by relative radiometric normalization
The remote sensing images of spectrum, using the enhancing multi-source Remote-sensing Image Fusion-eMulTiFuse based on spatial and temporal distributions of proposition, enter
Row multi-source Remote Sensing Image Fusion, obtain the remote sensing images with larger similarity spectral value.Time series remote sensing images are calculated
Normalized differential vegetation index, normalized differential vegetation index image sequence is obtained, be then based on Mann-Kendall trend-monitorings method and putting
Reliability is to make bilateral trend-monitoring in 95% fiducial range, obtains the degree of vegetation dynamic changes in study period, is based on
Theil Sen slopes, ask intermediate value to estimate variation tendency according to spectral value difference and time difference, obtain coupling relationship trend with becoming
The dynamic change result that change degree combines.Finally ground acquisition, high spatial resolution remote sense are used for the result of dynamic monitoring
The sampled results such as data acquisition carry out the calculating and analysis of precision.
Using Chongqing City Fuling District as survey region, the Dynamic Changes Analysis of sequence remote sensing image is carried out with this method, it is fixed
Property and the quantitative checking present invention performance.By 1999-2008 Landsat time serieses remote sensing images and 2009-2013
Year synthesizes a long period sequence remote sensing images through the environment star time series remote sensing images obtained by image co-registration, is based on
Mann-Kendall trend-monitorings method carries out trend analysis to Forest Resources Condition of the research area over 15 years and carried out corresponding
Checking, overall accuracy reach 87.18%, it is believed that can carry out dynamic monitoring to research area using the method.
The front and rear statistics value list of certain region of table 1 fusion
Obtained via Landsat remote sensing images sequence with environment No.1 remote sensing images sequence according to fusion method eMulTiFuse
To fused images, table 1 is the statistical result such as the average of a certain area image pixel, standard deviation before and after fusion, with reference to visual effect
Fruit and statistical result, it is known that remote sensing images have the space characteristics increasingly similar with Landsat after fusion.
Claims (1)
1. suitable for the dynamic monitoring method of multi-source time series remote sensing images, it is characterised in that comprise the following steps (as schemed
1)：
Step 1：The structure of multi-source Remote Sensing Images sequence
The feature that the non-pending remote sensing images of remotely-sensed data feature referred to herein have in itself, and refer to possess carry out time series
Some features that the remote sensing images of satellite image analysis should have：
First, time-series image radiation characteristic is consistent.Multitemporal Remote Sensing Images are (such as illumination, big often due to external factor
Gas radiation, sensor platform etc.) influence and cause atural object that " puppet change " occurs, therefore image radiation characteristic is unanimously critically important
's.
Secondly, image registration accuracy will height.Analysis to time series remote sensing images at present is more based on middle low-resolution image,
The analysis method of use is mostly the analysis method based on pixel, and this precision to geometric registration of imagery proposes higher requirement.
Finally, image picture elements value will be the reflection of atural object real features.Showed by image bad pixel, image bad track and complicated weather
Remote sensing images pixel caused by pixel distortion etc. as caused by (such as rain, cloud, mist or snow etc.) can not reflect truly table status,
The presence of this feature makes the accuracy of sequential remote Sensing Image Analysis.
Therefore, on the basis of suitable remote sensing images are chosen, remote sensing images are carried out using appropriate image pre-processing method pre-
Processing, makes remote sensing images sequence meet that features described above is vital.
Further, since remote sensing images are influenceed in imaging by factors such as sensor, flying platform, atmospheric radiations so that figure
As the distortion in terms of geometry and radiation being to some extent present, picture quality is influenceed, the use to image brings certain difficulty.
Using associated picture preprocess method, geometry, radiation error are reduced and eliminated as far as possible, is strengthened remote sensing images, also to be distant
Feel graphical analysis, interpretation, understanding, identification and information extraction and important leverage is provided.
Step 2：Multi-source Remote Sensing Image Fusion
(1) multi-source Remote-sensing Image Fusion-MulTiFuse based on pixel
When being merged using MulTiFuse fusion methods to two time series remote sensing images from different sensors, first,
The weight temporal sequence relation of two time serieses is derived, in order to maximize the statistical nature of this relation and consider abnormal feelings
Condition, weight optimization is associated to it before regression analysis is weighted to two time serieses；Then according to optimal time
Model is returned to merge multi-source Remote Sensing Images sequence.
The acquisition time of discrete-time series remote sensing images X time series pixel value is set as TX, time series X length is
kx；The acquisition time of discrete-time series remote sensing images Y time series pixel value is similarly set as TY, time series Y length
For ky.Assuming that above-mentioned two discrete-time series data length is inconsistent, image data acquiring time interval is different, and two
Time series acquisition time scope has lap.Be primarily based on linear interpolation method to X time series Y acquisition time
Enter row interpolation, obtain interpolation time sequence Xint, row interpolation is similarly entered to time series Y and obtains interpolation time sequence Yint, obtain
Time series as shown in following formula (1) and (2)：
Wherein, XlinAnd YlinX the and Y time serieses formed by interpolation are represented respectively.
Singular value may be obtained by above-mentioned simple linear interpolation., particularly in acquisition time interval longer time sequence,
Influenceing can be very big.In order to solve this problem, the observation difference (i.e. amplitude) of two time is introduced, its absolute value is bigger, returns
Return weight with regard to smaller；On the contrary, the absolute value of amplitude is smaller, it is bigger to return weight.
Next, definition is interpolated amplitude and normalizes the amplitude.First, it is defined on the timeI=1 ..., kxAbsolute width
Value mY, such as following formula (3)：
Wherein,WithRepresent that time series Y originates in time series X time respectivelyFront and rear.WithRespectively away from the X timesNearest Y observation.
To amplitude mYIt is normalized, sees below formula (4) and (5)：
Similar, define mX、WithBecause overlapping period X and Y length of time series is different, it is therefore desirable to by time series
Length is taken into account, rightWithIt is normalized again, such as following formula (6) and (7)：
Wherein, lxRefer to the timeQuantity, lYRefer to the timeQuantity.Return weightIt is defined as follows formula (8)：
Exponential weighting function ewf is introduced herein to adjust weight, such as following formula (9)：
Wherein, ewf=0,1,2, no weight, a weight, square weight are indicated respectively.Interpolation time sequence XintAnd YintPlan
Conjunction relation is by Seber's and Lee《Linear regression analysis(2nded.)》Weighted least-squares in one book
The fitting Return Law is calculated, and linear fit model is obtained by formula (10)：
Yint=a+bXint (10)
Wherein, a and b is constant.
(2) enhancing multi-source Remote-sensing Image Fusion-eMulTiFuse based on spatial and temporal distributions
Above-mentioned algorithm provides good thinking to multisource data fusion from the angle of temporal characteristics, and we are on this basis again
Unified, selection and amendment to unusual pixel are carried out on space characteristics to similar pixel, enhancing multi-source Remote Sensing Images data
Syncretizing effect.
This research carries out unsupervised classification to time series remote sensing images, remote sensing image data is clustered, as to unusual
The foundation that pixel is modified.Carrying out the conventional algorithm of cluster analysis has ISODATA (iteration self-organizing data analysis), K-
Means, chain method etc..ISODATA cluster centre is determined that this point is calculated with K-Means by the interative computation of this average
Method is similar, but the algorithm can draw the experience obtained by intermediate result, has self-organization, it can be said that compared with K-Means
In principle further, it is the comparatively ideal algorithm of classifying quality, this research is using this method progress cluster analysis.
According to the result of multispectral image unsupervised classification, the interpolation image that time prediction obtains is corresponded to the image and is divided
Class, then with mahalanobis distance (Mahalanobis distance) judge interpolation image pixel whether belong to such, geneva away from
From algorithm such as following formula (11)：
Wherein, XiPixel to be discriminated is represented,Refer to classification r average pixel value,Reflect pixel i and generic remaining
The diversity of pixel.
Because mahalanobis distance assumes that data are multi-source normalities in itself, therefore above-mentioned mahalanobis distance value approximation obeys chi square distribution.
We select the quantile of chi square distribution 90% criterion whether abnormal as pixel value is judged herein, are worth recognizing more than the quantile
It is abnormal for the pixel value, obtain average by the use of its generic pixel value and replace the pixel value and be used as correction result.
Step 3：Dynamic monitoring based on time series remote sensing images
The present invention carries out the dynamic monitoring of sequence remote sensing image using Mann-Kendall trend tests method, and Mann-Kendall becomes
Which type of regularity of distribution gesture method of inspection known pixel value need not have in advance without input parameter value, can reduce and make an uproar
The interference of sound, convenience of calculation, it can be used for the hypothesis testing of variable trends.Mann-Kendall trend tests typically with pixel or
Based on the single grid of person, by assuming that having inspected the variation tendency monitoring of atural object index.
In Mann-Kendall trend tests, make the following assumptions：Assume H in source0For：Time series data (X to be detected1,
X2,...,Xn), meet n independent, stochastic variables with being distributed, i.e. random distribution, without notable trend；Alternative hypothesis H1For：It is right
In all k, j≤n, and k ≠ j, XkAnd XjDistribution differ, sequence have enhancing or attenuation trend.Test statistics S is such as
Following formula (14) calculates：
Work as n>When 10, S is approximate to obey standardized normal distribution, and its average is 0, and variance can be using approximate representation as Var(S)=n (n-1)
(2n+5)/18.The normally distributed variable Z of standard is test statistics, for trend test, is calculated by following formula (16)：
In bilateral trend test, in given α confidence levels, if-Z1-α/2≤Z≤Z1-α/2, then it is assumed that null hypothesis H0Can
Lean on, receiving should be it is assumed that illustrates time series without significant change trend.If | Z | ＞ Z1-α/2, then it is assumed that assuming that H1Correctly, explanation
Sequence has obvious enhancing or attenuation trend.For test statistics Z, Z>0 shows that time series has enhancing trend, Z<0 table
Bright time series has attenuation trend.Confidence level be 90% significance test, Z0.95Equal to 1.28；Confidence level be 95% it is aobvious
When work property is examined, Z0.975Equal to 1.64；Confidence level be 99% significance test, Z0.9995Equal to 2.32.
Work as n<When 10, bilateral trend test is carried out based on test statistics S.When determining confidence level α, in significance test, such as
Fruit | S |≤Sα/2, then null hypothesis is received, the time series is without significant change trend；On the contrary, refusal null hypothesis, receives alternative vacation
If the time series is in significantly to rise or fall trend.For test statistics S, if S>0, it is upper to show that time series has
The trend of liter, if S=0, no visible trend, if S<0, show that time series has downward trend.
The intermediate value of the ratio of difference and time difference of the Theil Sen slopes based on different year gray value or exponential quantity and obtain, energy
Enough preferably to avoid noise disturbance, this is worth calculating no significance test.The computational methods of the value such as following formula (17)：
Wherein, 1<j<i<N (N represents the length of time series).Work as Slope>When 0, gray scale or index value sequence are in rising trend,
Work as Slope<When 0, the sequence is on a declining curve.The value of Theil Sen slopes is bigger, then illustrates that the degree of its change is bigger, instead
It is then smaller.
The trend analysis of vegetation coverage is carried out using Mann-Kendall trend tests method, with 95% confidence level during inspection,
Level of signifiance α=0.05 is taken, works as Slope>0 and | Z | when≤1.96, time series in rise but unconspicuous trend, can recognize
It is in less ascendant trend for vegetation coverage；Work as Slope>0 and | Z | during ＞ 1.96, time series is in obvious ascendant trend,
It is in obvious ascendant trend to think vegetation coverage；Work as Slope<0 and | Z | when≤1.96, time series in rise but it is unconspicuous
Trend, it is believed that vegetation coverage is in less downward trend；Work as Slope<0 and | Z | during ＞ 1.96, time series is under obvious
Drop trend, it is believed that vegetation coverage declines obvious.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710788089.6A CN107704807B (en) | 2017-09-05 | 2017-09-05 | Dynamic monitoring method based on multi-source remote sensing sequence image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710788089.6A CN107704807B (en) | 2017-09-05 | 2017-09-05 | Dynamic monitoring method based on multi-source remote sensing sequence image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107704807A true CN107704807A (en) | 2018-02-16 |
| CN107704807B CN107704807B (en) | 2020-12-29 |
Family
ID=61171988
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710788089.6A Expired - Fee Related CN107704807B (en) | 2017-09-05 | 2017-09-05 | Dynamic monitoring method based on multi-source remote sensing sequence image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107704807B (en) |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108564002A (en) * | 2018-03-22 | 2018-09-21 | 中国科学院遥感与数字地球研究所 | A kind of remote sensing image time series variation detection method and system |
| CN109035150A (en) * | 2018-04-24 | 2018-12-18 | 湖南省有色地质勘查研究院 | A kind of remote sensing image data defogging method based on middle-low alloy steels |
| CN109063663A (en) * | 2018-08-10 | 2018-12-21 | 武汉大学 | A kind of spissatus detection of timing remote sensing image and minimizing technology by slightly to essence |
| CN110765369A (en) * | 2019-09-11 | 2020-02-07 | 安徽先兆科技有限公司 | Real-time monitoring data processing method and system based on time-space attributes |
| CN110926430A (en) * | 2019-11-22 | 2020-03-27 | 海南省林业科学研究所 | Air-ground integrated mangrove forest monitoring system and control method |
| CN111742329A (en) * | 2020-05-15 | 2020-10-02 | 安徽中科智能感知产业技术研究院有限责任公司 | Mining typical ground object dynamic monitoring method and platform based on multi-source remote sensing data fusion and deep neural network |
| CN112150555A (en) * | 2020-08-27 | 2020-12-29 | 北京空间机电研究所 | In-orbit relative radiation calibration method for geosynchronous orbit area-array camera |
| CN112418554A (en) * | 2020-12-07 | 2021-02-26 | 江苏高速公路工程养护技术有限公司 | Method and device for determining influence factors of pavement performance decay |
| CN113553697A (en) * | 2021-06-24 | 2021-10-26 | 中国矿业大学(北京) | Long-time-sequence multi-source data-based vegetation disturbance analysis method for coal mining |
| CN114035021A (en) * | 2021-10-08 | 2022-02-11 | 北京航空航天大学 | Circuit fault prediction method based on EEMD-Prophet |
| CN114202705A (en) * | 2022-02-16 | 2022-03-18 | 清华大学 | Spectral feature time sequence construction method and system |
| CN114494909A (en) * | 2022-02-16 | 2022-05-13 | 中国科学院空天信息创新研究院 | Method and system for generating spatial distribution diagram of soybean growing season |
| CN116916166A (en) * | 2023-09-12 | 2023-10-20 | 湖南湘银河传感科技有限公司 | Telemetry terminal based on AI image analysis |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105184048B (en) * | 2015-07-23 | 2018-04-13 | 河南大学 | A kind of preparation method of the data of land utilization space year by year based on multi-source data |
-
2017
- 2017-09-05 CN CN201710788089.6A patent/CN107704807B/en not_active Expired - Fee Related
Cited By (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108564002A (en) * | 2018-03-22 | 2018-09-21 | 中国科学院遥感与数字地球研究所 | A kind of remote sensing image time series variation detection method and system |
| CN108564002B (en) * | 2018-03-22 | 2020-08-21 | 中国科学院遥感与数字地球研究所 | Method and system for detecting time sequence change of remote sensing image |
| CN109035150A (en) * | 2018-04-24 | 2018-12-18 | 湖南省有色地质勘查研究院 | A kind of remote sensing image data defogging method based on middle-low alloy steels |
| CN109035150B (en) * | 2018-04-24 | 2019-05-10 | 湖南省有色地质勘查研究院 | A kind of remote sensing image data defogging method based on middle-low alloy steels |
| CN109063663B (en) * | 2018-08-10 | 2021-08-03 | 武汉大学 | Thick cloud detection and removal method for time sequence remote sensing image from coarse to fine |
| CN109063663A (en) * | 2018-08-10 | 2018-12-21 | 武汉大学 | A kind of spissatus detection of timing remote sensing image and minimizing technology by slightly to essence |
| CN110765369A (en) * | 2019-09-11 | 2020-02-07 | 安徽先兆科技有限公司 | Real-time monitoring data processing method and system based on time-space attributes |
| CN110765369B (en) * | 2019-09-11 | 2022-03-29 | 安徽先兆科技有限公司 | Real-time monitoring data processing method and system based on time-space attributes |
| CN110926430A (en) * | 2019-11-22 | 2020-03-27 | 海南省林业科学研究所 | Air-ground integrated mangrove forest monitoring system and control method |
| CN111742329A (en) * | 2020-05-15 | 2020-10-02 | 安徽中科智能感知产业技术研究院有限责任公司 | Mining typical ground object dynamic monitoring method and platform based on multi-source remote sensing data fusion and deep neural network |
| CN111742329B (en) * | 2020-05-15 | 2023-09-12 | 安徽中科智能感知科技股份有限公司 | Mining typical feature dynamic monitoring method and platform based on multi-source remote sensing data fusion and deep neural network |
| CN112150555A (en) * | 2020-08-27 | 2020-12-29 | 北京空间机电研究所 | In-orbit relative radiation calibration method for geosynchronous orbit area-array camera |
| CN112150555B (en) * | 2020-08-27 | 2024-02-09 | 北京空间机电研究所 | On-orbit relative radiation calibration method for geosynchronous orbit area array camera |
| CN112418554A (en) * | 2020-12-07 | 2021-02-26 | 江苏高速公路工程养护技术有限公司 | Method and device for determining influence factors of pavement performance decay |
| CN113553697A (en) * | 2021-06-24 | 2021-10-26 | 中国矿业大学(北京) | Long-time-sequence multi-source data-based vegetation disturbance analysis method for coal mining |
| CN114035021A (en) * | 2021-10-08 | 2022-02-11 | 北京航空航天大学 | Circuit fault prediction method based on EEMD-Prophet |
| CN114202705A (en) * | 2022-02-16 | 2022-03-18 | 清华大学 | Spectral feature time sequence construction method and system |
| CN114494909A (en) * | 2022-02-16 | 2022-05-13 | 中国科学院空天信息创新研究院 | Method and system for generating spatial distribution diagram of soybean growing season |
| CN114202705B (en) * | 2022-02-16 | 2022-05-20 | 清华大学 | Spectral feature time sequence construction method and system |
| CN116916166A (en) * | 2023-09-12 | 2023-10-20 | 湖南湘银河传感科技有限公司 | Telemetry terminal based on AI image analysis |
| CN116916166B (en) * | 2023-09-12 | 2023-11-17 | 湖南湘银河传感科技有限公司 | Telemetry terminal based on AI image analysis |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107704807B (en) | 2020-12-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107704807A (en) | A kind of dynamic monitoring method based on multi-source remote sensing sequential images | |
| CN110378196B (en) | Road visual detection method combining laser point cloud data | |
| CN104751478B (en) | Object-oriented building change detection method based on multi-feature fusion | |
| Braaten et al. | Automated cloud and cloud shadow identification in Landsat MSS imagery for temperate ecosystems | |
| CN102799903B (en) | High-spectrum automatic cloud detection method based on space and spectral information | |
| Zhang et al. | Mapping individual tree species in an urban forest using airborne lidar data and hyperspectral imagery | |
| Morel et al. | Evaluating the potential to monitor aboveground biomass in forest and oil palm in Sabah, Malaysia, for 2000–2008 with Landsat ETM+ and ALOS-PALSAR | |
| Pignatti et al. | Evaluating Hyperion capability for land cover mapping in a fragmented ecosystem: Pollino National Park, Italy | |
| CN102540165B (en) | Method and system for preprocessing MODIS (Moderate-Resolution Imaging Spectroradiometer) surface albedo data | |
| CN106443701B (en) | Method for early warning before flood and waterlog based on sequential water range remote sensing image | |
| Behrangi et al. | Daytime precipitation estimation using bispectral cloud classification system | |
| CN109961460A (en) | A kind of multiple target method for inspecting based on improvement YOLOv3 model | |
| CN106886760B (en) | A kind of EO-1 hyperion Ship Detection combined based on empty spectrum information | |
| Su et al. | Improving MODIS sea ice detectability using gray level co-occurrence matrix texture analysis method: A case study in the Bohai Sea | |
| CN103049751A (en) | Improved weighting region matching high-altitude video pedestrian recognizing method | |
| CN115170979B (en) | Mining area fine land classification method based on multi-source data fusion | |
| Szantoi et al. | A tool for rapid post-hurricane urban tree debris estimates using high resolution aerial imagery | |
| CN103164695A (en) | Fruit identification method based on multi-source image information fusion | |
| CN102540271B (en) | Semi-supervised hyperspectral sub-pixel target detection method based on enhanced constraint sparse regression method | |
| CN106778629A (en) | Greenhouse recognition methods and device | |
| US20210286970A1 (en) | Cloud detection method based on landsat 8 snow-containing image | |
| CN107292920B (en) | A kind of multidate full-polarization SAR remote sensing imagery change detection method of joint classification | |
| Fu et al. | Evaluation of LAI estimation of mangrove communities using DLR and ELR algorithms with UAV, hyperspectral, and SAR images | |
| US7058511B2 (en) | Sub-visible cloud cover assessment: VNIR-SWIR | |
| Kasimati et al. | Investigation of the similarities between NDVI maps from different proximal and remote sensing platforms in explaining vineyard variability |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20201229 Termination date: 20210905 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |