CN103617629B - High-resolution remote sensing image vegetation index time series bearing calibration based on MODIS remote sensing image - Google Patents

Abstract

The invention discloses a kind of high-resolution remote sensing image vegetation index time series bearing calibration based on MODIS remote sensing image, does the method comprise HJ? is between resolution remote sense image NDVI seasonal effect in time series structure, HJ and low-to-medium altitude, the time series of the corresponding pixel of resolution remote sense image assisted whole inspection, HJ between NDVI and low-to-medium altitude? the abnormal section inspection of NDVI time series, HJ? the abnormal section of NDVI time series is proofreaied and correct four-stage. Can the method be revised HJ? NDVI time series, because climate reasons causes the abnormal of sequence curve form, improves its correctness and availability, provides valid data for study biological geochemical model under high-spatial and temporal resolution. This law is also applicable to the time series of certain vegetation index that comes from other Moderate-High Spatial Resolution Remote Sensing Image simultaneously.

Description

High-resolution remote sensing image vegetation index time sequence correction method based on MODIS remote sensing image

Technical Field

The invention relates to a high-resolution remote sensing image vegetation index time sequence correction method based on MODIS remote sensing images, and belongs to the technical field of geographic information systems.

Background

The remote sensing image time sequence is an image sequence which utilizes a series of remote sensing image data original images or index images derived from the original images to be arranged according to the time sequence, so that the remote sensing image time sequence can express the time variation characteristic of a remote sensing object. A data source is provided for analyzing from the time dimension using the remote sensing images.

The method is most widely applied at present and is a vegetation index time sequence of the remote sensing image with medium and low spatial resolution. The NDVI (normalized vegetation index, English is normalized Difference VegetationIndex) time series product provided by medium and low spatial resolution remote sensing satellite ModerateResolutionImagingSpectrophotometer (MODIS) is a vegetation index product widely used in the world. The product is derived from daily remote sensing image data which is subjected to atmospheric correction and surface bidirectional reflection, a targeted method is used for the data, the stability and reliability of the data quality are guaranteed, and a continuous high-time-resolution NDVI time sequence product covering the whole world can be provided. However, the maximum spatial resolution of modistvi is 250m, and the lower spatial resolution brings difficulty in studying finer terrain conditions. High spatial resolution and high temporal resolution are a pair of spears, and satellites with higher spatial resolution (such as Landsat and SPOT series satellites) cannot provide a high enough re-turn period and cannot construct a high temporal resolution time sequence. In contrast, many scholars have developed data fusion research, and use different sensor images to perform spatial and temporal fusion to make up for the contradiction, so that the fused images retain the dual characteristics of high temporal resolution and high spatial resolution as much as possible. Relevant studies have been reviewed by PohlandVanGenderen (1998), AmolinscandZhang (2007), ThomaandRanchen (2008). Although new research efforts continue to emerge, it is still not possible to fully preserve the fused image with all the details of the parent image.

With the emission of an environment disaster reduction small satellite constellation (HJ-1A/1 Bsmallsatellitectionstellation, HJ-1A/1B) emitted by China in 9 months of 2008, another solution is provided for constructing a remote sensing image time sequence with high spatial resolution and high temporal resolution. HJ-1A/1B can provide multispectral remote sensing images with spatial resolution of 30m and revisit period of 2 days, so that NDVI time sequence (HJNDDVI time sequence) with higher spatial-temporal resolution can be generated by using the HJ-1A/1B remote sensing images without fusion processing of high temporal resolution and high spatial resolution images. Theoretically, such a temporal resolution should be able to meet the requirements of NDVI time series for vegetation research, which some scholars have also applied to. However, the HJNDIVI time sequence still has the following problems: (1) the number of images acquired by the HJ satellite in unit time is far smaller than that of MODIS, and the influence of noise on the reduction of NDVI (normalized difference of absolute difference) values cannot be eliminated by using a maximum value synthesis (MVC) method; (2) the HJ-1A/1B remote sensing image does not provide a product for evaluating the pixel quality, and the generated NDVI time sequence cannot be adjusted according to the platform of the HJ-1A/1B remote sensing image, so that the influence of noise-polluted time phase is weakened. (3) The downloaded video is not radiation corrected and atmospheric corrected, and the HJNDDVI time sequence is not stable. The above problems cause the phenomena of low NDVI value persistence, abnormal fluctuation and the like in the application of the HJNDVI time sequence. In order to be able to accurately use the HJNDVI time series, it must be corrected. For the correction research of vegetation index time series, various methods for different data sources exist, such as a best exponential slope extraction algorithm (BISE), a Fourier fitting method, a Savitzky-Golay filtering and an asymmetric Gaussian function fitting method. But none of these solves the above-mentioned problem of the HJNDVI time series well.

NDVI, although not an intrinsic physical quantity, is related to the physical properties of vegetation (e.g., leaf area index, vegetation coverage, vegetation growth conditions, biomass, etc.). Therefore, regardless of which sensor calculates the NDVI time series, the NDVI time series will often exhibit seasonal variations in a seasonal or yearly cycle, and the NDVI time series of different sensors should have the same variation trend over the same time period in the same region. Because the quality of the NDVI time sequence of the medium-low spatial resolution remote sensing image is relatively reliable, the abnormal section of the HJNDVI time sequence can be corrected by referring to the NDVI time sequence of the medium-low spatial resolution remote sensing image, the time sequence abnormality caused by the climate can be corrected, the application requirement can be met, and effective data can be provided for researching a biological geochemical model under higher spatial resolution.

Disclosure of Invention

The technical problem solved by the invention is as follows: the method can correct abnormal sections of the HJNDVI time sequence, correct the time sequence abnormality caused by climate reasons, and provide effective data for researching a biological geochemical model under higher spatial resolution.

In order to solve the technical problems, the technical scheme provided by the invention is as follows: a vegetation index time sequence correction method for a high-resolution remote sensing image based on an MODIS remote sensing image comprises the following steps:

1) acquiring HJ-1A/1B remote sensing image data of a preset region and a preset time period and MODIS remote sensing data of the region in the same time period;

2) reading MODIS remote sensing data, and obtaining an NDVI value, an actual imaging day and a quality evaluation value of each pixel in an MODIS remote sensing image, wherein the quality evaluation values are divided into the following categories: "1" indicates no evaluation result, "0" indicates that the data quality is very good, "1" indicates that the data is available but the quality is general, "2" indicates that there is snow ice covering, "3" indicates that there is cloud interference;

3) respectively combining NDVI values of pixels with the same coordinate value in all MODIS remote sensing images into a first MODISNDVI time sequence according to the time sequence of the actual imaging day, and replacing the NDVI values of the pixels with quality evaluation values not equal to 0 and 1 in the first MODISNDVI time sequence by interpolation of adjacent NDVI values in the sequence to obtain a second MODISNDVI time sequence;

4) carrying out normalization processing on the vegetation index of each HJ-1A/1B remote sensing image pixel to obtain a HJNDDVI image;

5) respectively combining NDVI values of pixels with the same coordinate value in all HJNDDVI images into a first HJNDDVI time sequence corresponding to the coordinate value according to the time sequence of the actual imaging day, and interpolating the first HJNDVI time sequence according to the most adjacent NDVI value according to the actual imaging day of a second MODISNDVI time sequence to obtain a second HJNDVI time sequence consistent with the actual imaging day of the second MODISNDVI time sequence;

6) for each pixel of the MODIS remote sensing image in the preset area, finding out a pixel in the HJ-1A/1B remote sensing image which is spatially covered by the pixel, forming a time sequence pair by a second MODISNDVI time sequence and a second HJNDDVI time sequence which have a covering relation, and carrying out a co-integration test on the time sequence pair by taking the second MODISNDVI time sequence as a standard sequence;

7) subtracting a corresponding second HJNDDVI time sequence from a second MODISNDVI time sequence of the time sequence pair passing the coordination check to obtain a fluctuation sequence; carrying out abnormal section detection on the fluctuation sequence to obtain a normal fluctuation section and an abnormal fluctuation section, and dividing the second HJNDDVI time sequence into a normal HJNDDVI time section and an abnormal HJNDDVI time section according to the actual imaging days of the normal fluctuation section and the abnormal fluctuation section;

8) and correcting the abnormal HJNDVI time section in the second HJNDVI time sequence to obtain a high-resolution remote sensing image time sequence with higher time resolution in the preset region and the preset time section.

The further improvement of the technical scheme is as follows: the second MODISNDVI time series and the second HJNDVI time series are first Savitzky-Golay filtered before performing step 6).

The further improvement of the technical scheme is as follows: and 8) correcting the abnormal HJNDDVI time section by using the normal HJNDDVI time section and the second MODISNDVI time sequence of the corresponding time section.

The further improvement of the technical scheme is as follows: and 4) before normalizing each HJ-1A/1B remote sensing image in the step 4), performing radiation correction and geometric fine correction on each HJ-1A/1B remote sensing image.

Preferably, the method for performing the co-integration test on the time series pairs in the step 6) is as follows:

1A) judging whether a second MODISNDVI time sequence and a second HJNDDVI time sequence in the corresponding time sequence pair are of the same order and are single, and turning to the step 1B) if the second MODISNDVI time sequence and the second HJNDVI time sequence are of the same order and are single; if the order is not the same, the time sequence pair does not have the coordination relation;

1B) carrying out co-integration regression on the time sequence pairs of the same order single integration, using an ADF method to check the stationarity of regression residual errors, and if the residual errors are stable, the time sequence pairs are co-integrated and have a dependent relationship; otherwise, the two have no coordination relation.

Preferably, the method for performing the abnormal segment inspection on the fluctuation sequence in the step 7) is as follows:

2A) assuming that the second MODISNDVI time sequence is X and the second hjnndvi time sequence is Y, the central bound CL and the upper bound UCL of the fluctuation sequence are calculated as follows:

$\left\{\begin{array}{c}\mathrm{CL}=\stackrel{\‾}{U}\\ \mathrm{UCL}=\stackrel{\‾}{U}+\mathrm{A\σ}\end{array}\right.$

wherein, $\stackrel{\‾}{U}=\mathrm{\Σ}{u}_t/n,$

u_t=X_t-Y_t，

$\mathrm{\σ}=\sqrt{\frac{\mathrm{\Σ}(u_t/\stackrel{\‾}{U})}{n-1}},$

a is a control parameter with the value range of [1,3],

X_tis the NDVI value of the t-th element in the second MODISNDVI time series,

Y_tfor the NDVI value of the tth element in the second HJNDVI temporal sequence,

n is the number of elements of a second MODISNDVI time sequence or a second HJNDDVI time sequence;

2B) traversing the fluctuation sequence to obtain abnormal points of which the point values exceed the upper control limit UCL;

2C) if the number of the abnormal points continuously appearing between the central limit CL and the upper limit UCL is larger than L, the continuously appearing abnormal points form an abnormal section, and the value range of L is [5,30 ].

Preferably, the specific method for correcting the abnormal HJNDVI time section in step 8) is as follows:

3A) setting a second MODISNDVI time sequence corresponding to the normal fluctuation section as Xc and a second HJNDVI time sequence corresponding to the normal fluctuation section as Yc, performing a co-integration regression on the Xc and the Yc, and obtaining estimated values of regression parameters alpha and beta and a residual error sequence ec by a least square method;

wherein the co-integration regression formula is:

$\hat{Y}{c}_t=\mathrm{\α}+\mathrm{\βX}{c}_t$

$e{c}_t=Y{c}_t-\hat{Y}{c}_t$

in the formula, Xc_tThe NDVI value of the t-th element in the second MODISNDVI time series corresponding to the normal fluctuation section,

Yc_tthe NDVI value of the t-th element in the second HJNDVI temporal sequence corresponding to the normal fluctuation zone,

ec_tis the value of the t-th element in the residual error sequence;

3B) establishing an error correction model, and estimating β in the error correction model by a least square method₁、β₂And lambda; the formula of the error correction model is as follows:

$\mathrm{\Δ}\hat{Y}{c}_t={\mathrm{\β}}_1\mathrm{\ΔX}{c}_t+{\mathrm{\β}}_2\mathrm{\ΔY}{c}_{t-1}+\mathrm{\λ}(Y{c}_{t-1}-a-\mathrm{\βX}{c}_{t-1})$

in the formula,. DELTA.Xc_tThe difference between NDVI values of t +1 th element and t th element of the second MODISNDVI time sequence corresponding to the normal fluctuation section, i.e., Δ Xc_t=Xc_t+1-Xc_t；

ΔYc_t-1The difference between the NDVI values of the tth element and the tth-1 element of the second HJNDDVI time series corresponding to the normal fluctuation segment, Δ Yc_t-1=Yc_t-Yc_t-1；

3C) Let the second MODISNDVI time sequence be X, the second HJNDDVI time sequence be Y, X_tIs the NDVI value of the t-th element in the second MODISNDVI time series, Y_tThe NDVI value of the t-th element in the second HJNDV time sequence, the number of abnormal points of a certain abnormal HJNDV time section is N, and the NDVI value of the first element of the abnormal HJNDV time section is Y_KThe NDVI value of the last element being Y_K+N-1Then Y is_K-1、Y_K-2Belongs to the normal HJNDDVI time segment, and Y is obtained according to the following formula_KCorrected value Y'_K：

Y'_K=Y_K-1+β₁(X_K-X_K-1)+β₂(Y_K-1-Y_K-2)+λ·(Y_K-1-βX_K-1-α)

3D) From Y'_KReplacement of Y_KAnd calculating Y 'from the above formula'_K+1Recursion is done until Y 'is obtained by calculation'_K+N-1And is of Y'_K+N-1Replacement of Y_K+N-1Completing the NDVI value correction of the elements in the abnormal HJNDDVI time section;

3E) and after all abnormal HJNDDVI time sections are corrected, obtaining a corrected HJNDDVI time sequence.

Preferably, the MODIS remote sensing data are MODISMOD13Q1 data and MODISMYD13Q1 data.

The invention has the following beneficial effects:

1) in the present invention, NDVI is not an intrinsic physical quantity, but is related to the physical properties of vegetation (e.g., leaf area index, vegetation coverage, vegetation growth conditions, biomass, etc.). Therefore, regardless of which sensor calculates the NDVI time series, the NDVI time series will often exhibit seasonal variations in a seasonal or yearly cycle, and the NDVI time series of different sensors should have the same variation trend over the same time period in the same region. Because the quality of the NDVI time sequence of the medium-low spatial resolution remote sensing image is relatively reliable, the abnormal section of the HJNDVI time sequence is corrected by referring to the NDVI time sequence of the medium-low spatial resolution remote sensing image, the time sequence abnormality caused by the climate is corrected, the time sequence abnormality accords with the application requirement, and effective data is provided for researching a biological geochemical model under higher spatial resolution.

2) The previous research application of the HJNDIVI time sequence ignores the defect of an abnormal section in the sequence or finds that no compensation scheme is provided. The method provides a method for judging and correcting the abnormal segment of the HJNDDVI, and fills up the technical blank.

3) In the correction process of the HJNDDVI time sequence, the MODISNDVI time sequence with more reliable quality is introduced, so that a more accurate result is brought to the correction of the HJNDDVI time sequence. The invention provides a prerequisite and a correction model for correcting the HJNDDVI time sequence by using the MODISNDVI time sequence.

4) The HJNDDVI time sequence data corrected by the method can correct the abnormality of sequence curve form caused by climate reasons in the time sequence data with medium and high resolution, provides effective data for researching a biological geochemical model under high temporal and spatial resolution, and promotes the vegetation time sequence research to develop in a refined way.

5) The invention is not only suitable for HJ-1A/1B satellite data, but also has reference significance for other vegetation index time sequences formed by other medium-high resolution remote sensing images.

In summary, the invention provides a HJNDDVI time sequence correction method taking MODISNDVI time sequences as reference aiming at the data characteristics of HJ-1A/1B satellites, designs a set of step methods for finding and correcting abnormal sections of the HJNDDVI time sequences, can effectively improve the correctness and the usability of the HJ-1A/1BNDVI time sequences, and brings a data source with higher spatial resolution for the high-time-resolution NDVI time sequence correlation research.

Drawings

The vegetation index time series correction method of the high-resolution remote sensing image based on the MODIS remote sensing image is further explained with reference to the attached drawings.

FIG. 1 is a schematic flow chart of an embodiment of the present invention.

Fig. 2 is a remote sensing image coverage map of the embodiment.

FIG. 3 shows the result of the calibration verification of this embodiment.

Fig. 4 is a result of the HJNDVI time-series abnormal section inspection of the sample area A, B, C of fig. 2.

Fig. 5 is a result of the HJNDVI time-series correction of the sample area A, B, C of fig. 2.

Detailed Description

Examples

The coverage area of the remote sensing image adopted in the embodiment is shown in fig. 2, the vegetation coverage of the forest park in the old mountain country of Nanjing China is a massive mixed forest with deciduous forests, evergreen broad-leaved forests and coniferous forests, and the land types such as bare land, reservoir and residential land exist in other areas. This embodiment is implemented by MATLABR2012 a.

As shown in fig. 1, the method for correcting vegetation index time series of high-resolution remote sensing images based on MODIS remote sensing images in this embodiment includes the following steps:

1) and acquiring HJ-1A/1B remote sensing image data of a preset region and a preset time period and MODIS remote sensing data of the region and the time period.

The MODIS remote sensing data of the embodiment comprises MODISMOD13Q1 data and MODISMYD13Q1 data, firstly MOD13Q1 data and MYD13Q1 data (band number h28v05, both of which are 16 days Level3 vegetation index data) of Terra and Aqua in 1 month to 12 months in 2012 are obtained from NASA (website address: http:// ladssweb. nascom. NASA. gov/data/search. html), and the size of a single image grid is 231m after the single image grid is re-projected by MRT MOD (ISrojectionTool).

The HJ-1A/1B remote sensing image data of the embodiment is obtained from a chinese resource satellite application center (http:// www.cresda.com), and specifically is 216 scene multispectral remote sensing image data (cloud amount is less than 20%) covering the study area from 1 month to 12 months in 2009, and the size of a single image grid after sampling is 30 m.

2) Reading MODIS remote sensing data, and obtaining an NDVI value, an actual imaging day and a quality evaluation value of each pixel in an MODIS remote sensing image, wherein the quality evaluation values are divided into the following categories: "1" indicates no evaluation result, "0" indicates that the data quality is very good, "1" indicates that the data is available but the quality is general, "2" indicates that there is snow ice covering, and "3" indicates that there is cloud disturbance.

Recorded on the actual imaging day of this example is the value used for the NDVI "image on day N, 365/366 of the year". MODISMOD13Q1 and MYD13Q1 obtain a scene of data each day 16, and through the actual imaging day, the two types of data can be combined according to the time sequence of the actual imaging day to synthesize an NDVI time sequence with the time resolution being almost 8 days.

3) And respectively combining the NDVI values of the pixels with the same coordinate value in all the MODIS remote sensing images into a first MODISNDVI time sequence according to the time sequence of the actual imaging day, and replacing the NDVI values of the pixels with the quality evaluation values not equal to 0 and 1 in the first MODISNDVI time sequence by interpolation of adjacent NDVI values in the sequence to obtain a second MODISNDVI time sequence.

In this embodiment, after removing the NDVI values of MOD/MYD13Q1 in the modistvi time series, which are not equal to 0 and 1, the values obtained by performing linear interpolation on the two nearest NDVI are used instead. The linear interpolation method adopted in the embodiment is as follows:

$X_i=\frac{X_{i^{+}}-X_{i^{-}}}{i^{+}-i^{-}}*(i-i^{-})$

wherein X_tFor the NDVI value of the t-th element in the modistvi time series,

i is a time when the quality evaluation value is not "0" or "1",

i⁺the most adjacent time instant of i equals "0" or "1",

i^-i is the most adjacent previous time instant equal to "0" or "1".

4) And carrying out normalization processing on the vegetation index of the pixel of each HJ-1A/1B remote sensing image to obtain the HJNDDVI image.

The conversion of the HJNDVI image obtained by normalizing the HJ-1A/1B remote sensing image pixels in this embodiment is obtained by the following formula (each pixel in the image is calculated by the following formula):

$\mathrm{NDVI}=\frac{\mathrm{band}4-\mathrm{band}3}{\mathrm{band}4+\mathrm{band}3}$

wherein the band4 is the apparent reflectivity (rho) of the fourth wave band in the HJ-1A/1B remote sensing image₄），

band3 is the value of the apparent reflectivity (rho) of the third band in the HJ-1A/1B remote sensing image₃）。

5) And respectively combining NDVI values of pixels with the same coordinate value in all the HJNDDVI images into a first HJNDDVI time sequence corresponding to the coordinate value according to the time sequence of the actual imaging day, and interpolating the first HJNDVI time sequence according to the most adjacent NDVI value according to the actual imaging day of a second MODISNDVI time sequence to obtain a second HJNDVI time sequence consistent with the actual imaging day of the second MODISNDVI time sequence.

In this embodiment, a linear interpolation method is adopted when the first hjnndvi time sequence is interpolated according to the nearest NDVI value on the actual imaging day of the second modistvi time sequence, and the basic linear interpolation method is the same as the linear interpolation method in step 3).

Thus, the obtained second HJNDVI time-series and second MODISNDVI time-series have the same imaging date distribution.

6) And for each pixel of the MODIS remote sensing image in the preset area, finding the pixel in the HJ-1A/1B remote sensing image covered by the pixel, forming a time sequence pair by a second MODISNDVI time sequence and a second HJNDVI time sequence which have a covering relation, and carrying out a co-integration test on the time sequence pair by taking the second MODISNDVI time sequence as a standard sequence.

In the embodiment, HJ-1A/1B remote sensing image pixels which can be covered by each MODIS remote sensing image pixel are searched, pixels having an inclusion relationship are called as "corresponding pixels", the corresponding relationship is recorded, and the "corresponding pixels" in the research area of the embodiment are shown in fig. 2. The second MODISNDVI time sequence and the second HJNDVI time sequence of the corresponding picture element are referred to as "corresponding sequences". In this embodiment, the method for performing the co-integration test on all the "corresponding sequences" includes the following steps:

1A) judging whether a second MODISNDVI time sequence and a second HJNDDVI time sequence in the corresponding time sequence pair are of the same order and are single, and turning to the step 1B) if the second MODISNDVI time sequence and the second HJNDVI time sequence are of the same order and are single; and if the order is not the same, the time sequence pair has no coordination relation.

The time sequence with the statistical characteristics not changing along with the time translation is called a stable time sequence, if a certain time sequence is stable after N-order difference, the time sequence is called an N-order single integer sequence, and N is a single integer order. And acquiring a second HJNDDVI time sequence of one HJ pixel and a second MODISNDVI time sequence covering the HJDVI time sequence, and respectively checking the single integer order of the two time sequences.

The method for judging the single integer order adopted by the embodiment comprises the following steps: AugmentDickey-Fullertest (ADF) test was used to determine whether the sequence was smooth, and the ADF test was implemented in statistical software Eviews 5.0. If a group of corresponding sequences (a second modistvi time sequence and a second hjnndvi time sequence) are stable after N-order difference, the sequence is called an N-order single integer sequence, and the single integer numbers of the two time sequences are the same, the two time sequences are the same order single integer.

1B) Carrying out co-integration regression on the time sequence pairs of the same order single integration, using an ADF method to check the stationarity of regression residual errors, and if the residual errors are stable, the time sequence pairs are co-integrated and have a dependent relationship; otherwise, the two have no coordination relation.

In the embodiment, firstly, the co-integration regression is carried out on the corresponding sequences of the same order single integration, then the stability of the regression residual is checked by using the ADF, if the residual is stable, the two time sequences are co-integrated, and the two time sequences have a dependent relationship; if the residual error is not stable, the two do not have a co-integration relationship. The formula of the co-integration regression is:

${\hat{Y}}_t=\mathrm{\α}+\mathrm{\β}{X}_i$

$e_t=Y_t-{\hat{Y}}_t$

wherein, X_tIs the NDVI value of the t-th element in the second MODISNDVI time series,

Y_tfor the NDVI value of the tth element of the second HJNDVI time series,

e_tis the value of the t-th element in the estimated residual error sequence;

regression parameters α and β are estimated using the ordinary least squares method (OLS).

In this embodiment, the following formula of non-intercept term and non-trend term is selected, and H is assumed to be_t:<Residual error e at 0 using t-test_tWhether or not zero hypothesis H is satisfied₀:=0：

$\mathrm{\&Delta;}{e}_t=\mathrm{\&delta;}{e}_{t-1}+\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{\mathrm{\&theta;}}_i\mathrm{\&Delta;}{e}_{t-1}+{\mathrm{\&epsiv;}}_t$

Wherein, Δ e_tRepresenting the magnitude of the change of the residual between time t +1 and time t,

、θ_ifor the coefficients to be estimated of the polynomial,

p is the sequence correlation order set by the user,

_tis the adjustment term of the polynomial.

To e_tThe test of (2) was completed with a t-test, and the test thresholds for different sample volumes in the bivariate case are shown in table 1. Rejecting null hypothesis H if, at a certain hysteresis order, the value of the t statistic is less than the critical value at the corresponding sample capacity₀=0, consider e_tThere is no unit root and it is stable, at this time the two time series are coordinated.

TABLE 1 bivariate synergistic ADF test thresholds

This example performs a co-integration regression on the residual e_tJudgment of stationarity and Pair E_tAll are tested inThe statistical software Eviews5.0 is implemented, and the result of the collaborative inspection is shown in FIG. 3.

7) Subtracting a corresponding second HJNDDVI time sequence from a second MODISNDVI time sequence of the time sequence pair passing the coordination check to obtain a fluctuation sequence; and performing abnormal section detection on the fluctuation sequence to obtain a normal fluctuation section and an abnormal fluctuation section, and dividing the second HJNDDVI time sequence into a normal HJNDDVI time section and an abnormal HJNDDVI time section according to the actual imaging days of the normal fluctuation section and the abnormal fluctuation section.

In this embodiment, the fluctuation sequence is denoted as U, and a sigma quality control method is adopted as a method for performing abnormal section inspection on the fluctuation sequence U. The method comprises the following specific steps:

2A) assuming that the second MODISNDVI time sequence is X and the second hjnndvi time sequence is Y, the central bound CL and the upper bound UCL of the fluctuation sequence are calculated as follows:

$\left\{\begin{array}{c}\mathrm{CL}=\stackrel{\&OverBar;}{U}\\ \mathrm{UCL}=\stackrel{\&OverBar;}{U}+\mathrm{A\&sigma;}\end{array}\right.$

wherein, $\stackrel{\&OverBar;}{U}=\mathrm{\&Sigma;}{u}_t/n,$

u_t=X_t-Y_t，

$\mathrm{\&sigma;}=\sqrt{\frac{\mathrm{\&Sigma;}(u_t/\stackrel{\&OverBar;}{U})}{n-1}},$

a is a control parameter with the value range of [1,3],

X_tis the NDVI value of the t-th element in the second MODISNDVI time series,

Y_tfor the NDVI value of the tth element in the second HJNDVI temporal sequence,

n is the number of elements of a second MODISNDVI time sequence or a second HJNDDVI time sequence;

2B) traversing the fluctuation sequence to obtain abnormal points of which the point values exceed the upper control limit UCL;

2C) if the number of the abnormal points continuously appearing between the central limit CL and the upper limit UCL is larger than L, the continuously appearing abnormal points form an abnormal section, and the value range of L is [5,30 ].

When the value of a certain point on the second HJNDVI time sequence exceeds the upper control limit UCL, it indicates that the second HJNDVI time sequence may be interfered by stronger pollution at this point. The points at which the point values appear consecutively on the side of the boundary are called "chains", and the number of points contained therein is called the chain length. If the chain length of a segment between CL and UCL is greater than L, then it is assumed that the second HJNDDVI time-series may have a lower persistence of NDVI values during this time due to climatic factors. According to the above-mentioned criterion several abnormal sections in the time sequence can be obtained, and the sections which do not meet the abnormal section judgement criterion are normal sections.

The result of the abnormal HJNDVI temporal section check of the sample area A, B, C in fig. 2 in this embodiment is shown in fig. 4.

8) And correcting the abnormal HJNDVI time section in the second HJNDVI time sequence to obtain a high-resolution remote sensing image time sequence with higher time resolution in the preset region and the preset time section.

In this embodiment, the abnormal HJNDVI time segment is corrected by the normal HJNDVI time segment and the second modindvi time sequence of the corresponding time segment, and the specific method is as follows:

3A) setting a second MODISNDVI time sequence corresponding to the normal fluctuation section as Xc and a second HJNDVI time sequence corresponding to the normal fluctuation section as Yc, performing a co-integration regression on the Xc and the Yc, and obtaining estimated values of regression parameters alpha and beta and a residual error sequence ec by a least square method;

wherein the co-integration regression formula is:

$\hat{Y}{c}_t=\mathrm{\&alpha;}+\mathrm{\&beta;X}{c}_t$

$e{c}_t=Y{c}_t-\hat{Y}{c}_t$

in the formula, Xc_tThe NDVI value of the t-th element in the second MODISNDVI time series corresponding to the normal fluctuation section,

Yc_tthe NDVI value of the t-th element in the second HJNDVI temporal sequence corresponding to the normal fluctuation zone,

ec_tis the value of the t-th element in the residual error sequence;

3B) establishing an error correction model, and estimating β in the error correction model by a least square method₁、β₂And lambda; the formula of the error correction model is as follows:

$\mathrm{\&Delta;}\hat{Y}{c}_t={\mathrm{\&beta;}}_1\mathrm{\&Delta;X}{c}_t+{\mathrm{\&beta;}}_2\mathrm{\&Delta;Y}{c}_{t-1}+\mathrm{\&lambda;}(Y{c}_{t-1}-a-\mathrm{\&beta;X}{c}_{t-1})$

in the formula,. DELTA.Xc_tThe difference between NDVI values of t +1 th element and t th element of the second MODISNDVI time sequence corresponding to the normal fluctuation section, i.e., Δ Xc_t=Xc_t+1-Xc_t；

ΔYc_t-1The difference between the NDVI values of the tth element and the tth-1 element of the second HJNDDVI time series corresponding to the normal fluctuation segment, Δ Yc_t-1=Yc_t-Yc_t-1；

3C) Let the second MODISNDVI time sequence be X, the second HJNDDVI time sequence be Y, X_tIs the NDVI value of the t-th element in the second MODISNDVI time series, Y_tThe NDVI value of the t-th element in the second HJNDV time sequence, the number of abnormal points of a certain abnormal HJNDV time section is N, and the NDVI value of the first element of the abnormal HJNDV time section is Y_KThe NDVI value of the last element being Y_K+N-1Then Y is_K-1、Y_K-2Belongs to the normal HJNDDVI time segment, and Y is obtained according to the following formula_KCorrected value Y'_K：

Y'_K=Y_K-1+β₁(X_K-X_K-1)+β₂(Y_K-1-Y_K-2)+λ·(Y_K-1-βX_K-1-α)

3D) From Y'_KReplacement of Y_KAnd calculating Y 'from the above formula'_K+1Recursion is done until Y 'is obtained by calculation'_K+N-1And is of Y'_K+N-1Replacement of Y_K+N-1Completing the NDVI value correction of the elements in the abnormal HJNDDVI time section;

3E) and after all abnormal HJNDDVI time sections are corrected, obtaining a corrected HJNDDVI time sequence.

The HJNDVI time-series correction results for the sample area A, B, C in fig. 2 in this embodiment are shown in fig. 5.

The embodiment can be further modified as follows:

a) before normalizing each HJ-1A/1B remote sensing image in the step 4), firstly carrying out radiation correction and geometric fine correction on each high-rate remote sensing image.

In the embodiment, the radiation correction is divided into two steps of converting the pixel value of the original image into the radiation value at the sensor and calculating the apparent reflectivity. The formula for converting the pixel value of the HJ-1A/1B image into the radiation value at the sensor is as follows:

$L_{\mathrm{\&lambda;}}=\frac{{\mathrm{DN}}_{\mathrm{\&lambda;}}}{A_{\mathrm{\&lambda;}}}+\mathrm{bia}{s}_{\mathrm{\&lambda;}}$

in the formula, L_λFor the spectral radiance value of the band lambda at the sensor,

A_λin order to scale the gain of the coefficient absolutely,

DN_λis the pixel value of the image in the lambda band,

bias_λis the absolute scaling coefficient offset.

The spectral radiance values must be converted to the apparent reflectance at the sensor for comparison, the conversion equation is (each pixel in the image is calculated as follows):

${\mathrm{\&rho;}}_{\mathrm{\&lambda;}}=\frac{\mathrm{\&pi;}\&CenterDot;L_{\mathrm{\&lambda;}}\&CenterDot;D^2}{{\mathrm{ESUN}}_{\mathrm{\&lambda;}}\&CenterDot;\cos {\mathrm{\&theta;}}_s}$

in the formula,_ρλthe picture element is the apparent reflectivity of the top layer of the atmosphere, D is the distance between the sun and the ground, ESUN_λIs a top layer of atmosphereThe average solar irradiance, theta, is the zenith angle of the sun.

The parameter data in the above formula can be obtained from the image header file and the China resource satellite application center. The geometric fine correction of each HJ-1A/1B remote sensing image can be completed by the function of image geometric correction of version ENVI5.0 of remote sensing image processing software.

b) Performing Savitzk on the second MODISNDVI time sequence and the second HJNDDVI time sequence before performing step 6)_y-Gola_yAnd (6) filtering.

Savitzk selected for use in this example_y-Gola_yThe filter formula is as follows:

in the formula, Y_j ^*Is the NDVI value of the jth element of the filtered second modistvi time series or second HJNDVI time series,

Y_j+1the NDVI value of the j +1 th element of the second modistvi time series or the second HJNDVI time series,

m is a positive integer, 2m +1 is the size of the sliding window,

C_iis the coefficient of the ith element of the second MODISNDVI time sequence or the second HJNDDVI time sequence.

The averaging method can retain the alternating characteristics of seasonal peaks and troughs of the NDVI sequence.

The invention is not limited to the specific technical solutions described in the above embodiments, and all technical solutions formed by equivalent substitutions are within the scope of the claims of the invention.

Claims (7)

1. A vegetation index time sequence correction method for a high-resolution remote sensing image based on an MODIS remote sensing image comprises the following steps:

1) acquiring HJ-1A/1B remote sensing image data of a preset region and a preset time period and MODIS remote sensing data of the region in the same time period;

2) reading MODIS remote sensing data, and obtaining an NDVI value, an actual imaging day and a quality evaluation value of each pixel in an MODIS remote sensing image, wherein the quality evaluation values are divided into the following categories: "1" indicates no evaluation result, "0" indicates that the data quality is very good, "1" indicates that the data is available but the quality is general, "2" indicates that there is snow ice covering, "3" indicates that there is cloud interference;

3) respectively combining NDVI values of pixels with the same coordinate value in all MODIS remote sensing images into a first MODISNDVI time sequence according to the time sequence of the actual imaging day, and replacing the NDVI values of the pixels with quality evaluation values not equal to 0 and 1 in the first MODISNDVI time sequence by interpolation of adjacent NDVI values in the sequence to obtain a second MODISNDVI time sequence;

4) carrying out normalization processing on the vegetation index of each HJ-1A/1B remote sensing image pixel to obtain a HJNDDVI image;

5) respectively combining NDVI values of pixels with the same coordinate value in all HJNDDVI images into a first HJNDDVI time sequence corresponding to the coordinate value according to the time sequence of the actual imaging day, and interpolating the first HJNDVI time sequence according to the most adjacent NDVI value according to the actual imaging day of a second MODISNDVI time sequence to obtain a second HJNDVI time sequence consistent with the actual imaging day of the second MODISNDVI time sequence;

6) for each pixel of the MODIS remote sensing image in the preset area, finding out a pixel in the HJ-1A/1B remote sensing image which is spatially covered by the pixel, forming a time sequence pair by a second MODISNDVI time sequence and a second HJNDDVI time sequence which have a covering relation, and carrying out a co-integration test on the time sequence pair by taking the second MODISNDVI time sequence as a standard sequence;

7) subtracting a corresponding second HJNDDVI time sequence from a second MODISNDVI time sequence of the time sequence pair passing the coordination check to obtain a fluctuation sequence; carrying out abnormal section detection on the fluctuation sequence to obtain a normal fluctuation section and an abnormal fluctuation section, and dividing the second HJNDDVI time sequence into a normal HJNDDVI time section and an abnormal HJNDDVI time section according to the actual imaging days of the normal fluctuation section and the abnormal fluctuation section;

8) and correcting the abnormal HJNDVI time section in the second HJNDVI time sequence to obtain a high-resolution remote sensing image time sequence with higher time resolution in the preset region and the preset time section.

2. The MODIS remote sensing image-based vegetation index time series correction method for the high-resolution remote sensing image according to claim 1, characterized in that: the second MODISNDVI time sequence and the second HJNDVI time sequence are first Savitzky-Golay filtered before performing step 6).

3. The MODIS remote sensing image-based vegetation index time series correction method for the high-resolution remote sensing image according to claim 1, characterized in that: and 8) correcting the abnormal HJNDDVI time section by using the normal HJNDDVI time section and the second MODISNDVI time sequence of the corresponding time section.

4. The MODIS remote sensing image-based vegetation index time series correction method for the high-resolution remote sensing image according to claim 1, characterized in that: and 4) before normalizing each HJ-1A/1B remote sensing image in the step 4), performing radiation correction and geometric fine correction on each HJ-1A/1B remote sensing image.

5. The MODIS remote sensing image-based vegetation index time series correction method for the high-resolution remote sensing image according to claim 1, wherein the method for performing the collaborative inspection on the time series pair in the step 6) is as follows:

1A) judging whether a second MODISNDVI time sequence and a second HJNDDVI time sequence in the corresponding time sequence pair are of the same order and are single, and turning to the step 1B) if the second MODISNDVI time sequence and the second HJNDVI time sequence are of the same order and are single; if the order is not the same, the time sequence pair does not have the coordination relation;

1B) carrying out co-integration regression on the time sequence pairs of the same order single integration, using an ADF method to check the stationarity of regression residual errors, and if the residual errors are stable, the time sequence pairs are co-integrated and have a dependent relationship; otherwise, the two have no coordination relation.

6. The MODIS remote sensing image-based vegetation index time series correction method for the high resolution remote sensing image according to claim 3, wherein the specific method for correcting the abnormal HJNDI time section in the step 8) is as follows:

3A) setting a second MODISNDVI time sequence corresponding to the normal fluctuation section as Xc and a second HJNDVI time sequence corresponding to the normal fluctuation section as Yc, performing a co-integration regression on the Xc and the Yc, and obtaining estimated values of regression parameters alpha and beta and a residual error sequence ec by a least square method;

wherein the co-integration regression formula is:

$\hat{Y}{c}_t=\mathrm{\&alpha;}+{\mathrm{\&beta;Xc}}_t$

${\mathrm{ec}}_t={\mathrm{Yc}}_t-\hat{Y}{c}_t$

in the formula, Xc_tThe NDVI value of the t-th element in the second MODISNDVI time series corresponding to the normal fluctuation section,

Yc_tthe NDVI value of the t-th element in the second HJNDVI temporal sequence corresponding to the normal fluctuation zone,

ec_tis the value of the t-th element in the residual error sequence;

3B) establishing an error correction model, and estimating β in the error correction model by a least square method₁、β₂And lambda; the formula of the error correction model is as follows:

$\mathrm{\&Delta;}\hat{Y}{c}_t={\mathrm{\&beta;}}_1{\mathrm{\&Delta;Xc}}_t+{\mathrm{\&beta;}}_2{\mathrm{\&Delta;Yc}}_{t-1}+\mathrm{\&lambda;}({\mathrm{Yc}}_{t-1}-a-{\mathrm{\&beta;Xc}}_{t-1})$

in the formula, △ Xc_tThe difference between NDVI values of t +1 th element and t th element of the second MODISNDVI time sequence corresponding to the normal fluctuation segment, i.e., △ Xc_t＝Xc_t+1-Xc_t；

△Yc_t-1The difference between the NDVI values of the t-th element and the t-1-th element of the second HJNDDVI time series corresponding to the normal fluctuation segment, △ Yc_t-1＝Yc_t-Yc_t-1；

3C) Let the second MODISNDVI time sequence be X, the second HJNDDVI time sequence be Y, X_tIs the NDVI value of the t-th element in the second MODISNDVI time series, Y_tThe NDVI value of the t-th element in the second HJNDV time sequence, the number of abnormal points of a certain abnormal HJNDV time section is N, and the NDVI value of the first element of the abnormal HJNDV time section is Y_KThe NDVI value of the last element being Y_K+N-1Then Y is_K-1、Y_K-2Belongs to the normal HJNDDVI time segment, and Y is obtained according to the following formula_KCorrected value Y'_K：

Y'_K＝Y_K-1+β₁(X_K-X_K-1)+β₂(Y_K-1-Y_K-2)+λ·(Y_K-1-βX_K-1-α)

3D) From Y'_KReplacement of Y_KAnd using the above formula to calculateY'_K+1Recursion is done until Y 'is obtained by calculation'_K+N-1And is of Y'_K+N-1Replacement of Y_K+N-1Completing the NDVI value correction of the elements in the abnormal HJNDDVI time section;

3E) and after all abnormal HJNDDVI time sections are corrected, obtaining a corrected HJNDDVI time sequence.

7. The MODIS remote sensing image-based vegetation index time series correction method for the high-resolution remote sensing image according to claim 1, characterized in that: the MODIS remote sensing data are MODISMOD13Q1 and MODISMYD13Q1 data.