CN107103584A - A kind of production high-spatial and temporal resolution NDVI weighted based on space-time method - Google Patents

Abstract

The present invention provides a kind of method based on spatio-temporal weighting of producing high spatio-temporal resolution NDVI, comprising: the NDVI data of the low spatial resolution of _t1 time and the NDVI data of high spatial resolution and the low spatial resolution of t _p time Two kinds of NDVI increments of high spatial resolution data in time are obtained through mixed pixel decomposition and spatial interpolation of NDVI data; the Bayesian model averaging method is used to calculate the weight combination of the two NDVI increments to obtain each image The NDVI synthesis increment of the element; finally, add the NDVI synthesis increment to the NDVI with high spatial resolution at time t ₁ to obtain the required high spatio-temporal resolution NDVI. The invention performs Bayesian synthesis on two kinds of NDVI increments, fully utilizes the time-space correlation between geographic data, and provides an effective method for producing vegetation indices with high time-space resolution.

Description

A method of producing high spatiotemporal resolution NDVI based on spatiotemporal weighting

technical field

The invention relates to a multi-source remote sensing data fusion method, which is mainly used to produce a normalized difference vegetation index (NDVI) with high temporal and spatial resolution, and belongs to the category of temporal and spatial data fusion.

Background technique

The Normalized Difference Vegetation Index (NDVI for short) is a widely used vegetation index. NDVI values can reflect the growth status of vegetation and are directly related to physiological parameters of vegetation such as chlorophyll content, water content and leaf area index.

CN 102831310 B discloses a method for constructing high-spatial-resolution time-series data. Its background technology section describes the relevant knowledge of NDVI data in detail, and mentions two existing NDVI data: Landsat NDVI data and MODIS NDVI data data. Among them, the spatial resolution of NDVI values in Landsat NDVI data is 30m*30m, and the time resolution is 16 days, while the maximum spatial resolution of NDVI values in MODIS NDVI data is 250m*250m, and the time resolution is 1 day. That is to say, Landsat NDVI data has higher spatial resolution and lower temporal resolution, while MODIS NDVI data has lower spatial resolution and higher temporal resolution.

In order to obtain NDVI data with high spatial resolution and high temporal resolution at the same time, the existing technology has developed a variety of methods for fusing Landsat NDVI data and MODIS NDVI data to obtain NDVI data with high spatial resolution, mainly including Three types: (1) method based on weight function; (2) method based on mixed pixel decomposition; (3) method based on mixed model.

Among them, (1) the weight function-based method is typically represented by STARFM (Spatial and Temporal Adaptive Reflectance Fusion Model), see Gao, F., Masek, J., Schwaller, M., & Hall , F. (2006), "Fusing Landsat and MODIS Surface Reflectance: Predicting Daily Landsat Surface Reflectance", IEEE Transactions on Earth Science and Remote Sensing, 44, 2207-2218. (Gao, F., Masek, J., Schwaller, M., & Hall, F. (2006). On the blending of the Landsat and MODIS surface reflection: Predicting daily Landsat surface reflection. IEEE Transactions on Geoscience and Remote Sensing, 44, 2207 -2218). This type of method focuses on the search for similar pixels and the determination of their weights. The similarity between the current pixel and surrounding pixels in terms of time, space, and spectrum is used to determine similar pixels and calculate the corresponding weights. The average increment of is used as the increment of the current cell. This type of method often has better results in areas with higher purity.

(2) The method based on mixed pixel decomposition is represented by NDVI-LMGM (NDVI-Linear Mixing Growth Model, NDVILinear Mixing Growth Model), see Rao, Y., Zhu, X., Chen, J., & Wang, J. (2015), "An upgraded method of using multi-temporal MODIS NDVI data and Landsat TM/ETM+ images to produce high temporal and spatial resolution NDVI time series data", Remote Sensing, 7, 7865-7891 (Rao, Y., Zhu, X., Chen, J., & Wang, J.(2015). An Improved Method for Producing High Spatial-Resolution NDVI Time Series Datasets with Multi-Temporal MODIS NDVI Data and Landsat TM/ETM+Images. RemoteSensing, 7, 7865-7891) . This type of method calculates the time increment of each pixel by decomposing the mixed pixels, so as to obtain the NDVI data of the predicted time. For the specific content of this method, please refer to the aforementioned prior art CN 102831310B. The two are essentially Are the same. This type of method can also have better results for areas with a relatively high degree of heterogeneity, but it assumes that the land cover status does not change, which may not be satisfied.

(3) The method based on the mixed model is represented by FSDAF (Flexible Spatiotemporal Data Fusion method), see Zhu, X., Helmer, E.H., Gao, F., Liu, D., Chen, J., & Lefsky, M.A. (2016), "A method for flexible spatiotemporal fusion of satellite imagery with different resolutions", Remote Sensing of the Environment, 172, 165-177 (Zhu, X., Helmer, E.H., Gao, F., Liu , D., Chen, J., & Lefsky, M.A. (2016). A flexible spatiotemporal method for fusing satellite images with different resolutions. Remote Sensing of Environment, 172, 165-177). This type of method combines a variety of models, such as weight function model, mixed pixel decomposition model, and spatial interpolation model, to calculate the increment of each pixel on the high-resolution image. Such methods are often complex, but have relatively high precision and can capture land cover changes to a certain extent.

The above three types of methods have achieved relatively good results in their respective applicable fields. However, none of them considered the method of using two periods of low spatial resolution data to interpolate to high spatial resolution. The increment between the two periods of interpolation results can reflect the increment of the pixel scale to a certain extent, and it also includes land Override information on changes. The invention synthesizes the increment of mixed pixel decomposition and the increment of space interpolation based on the Bayesian model averaging method, and can effectively improve the accuracy of data fusion. It provides a new way to produce NDVI with high spatio-temporal resolution.

Contents of the invention

The technical problem to be solved by the present invention is to provide a method for producing NDVI with high spatiotemporal resolution based on spatiotemporal weighting, so as to reduce or avoid the aforementioned problems.

In order to solve the above technical problems, the present invention proposes a method for producing high spatio-temporal resolution NDVI based on spatiotemporal weighting, which is used to combine Landsat NDVI data with high spatial resolution and low temporal resolution with low spatial resolution and high temporal resolution The MODIS NDVI data fusion of rate obtains the NDVI data of the high _spatio -temporal resolution of prediction time t from reference time t ₁ , described method comprises the following steps:

Step A: Starting from the same reference time _t1 , using MODIS NDVI data and Landsat NDVI data, respectively obtain the first increment △NDVI _d and the second increment △NDVI of each pixel on NDVI at the prediction time t _{p s} ; among them,

The first increment ΔNDVI _d refers to, from the reference time t ₁ to the forecast time t _p , assuming that the NDVI increment of each pixel of the MODIS NDVI data is equal to the different features of the Landsat NDVI data under the same area The average value of the NDVI increment of the pixels of the category is calculated by using the MODIS NDVI data and Landsat NDVI data at time t ₁ and the MODIS NDVI data at time t _p to obtain the Landsat NDVI data from time t _p to time t ₁ The NDVI increment of pixels, which is the first increment ΔNDVI _d ;

The second increment ΔNDVI _s refers to directly interpolating MODIS pixels with low spatial resolution at time t ₁ and time t _p to high The Landsat pixel of spatial resolution goes up, and then the interpolation result of time t _p is subtracted from the interpolation result of time t ₁ , and the result is the second increment ΔNDVI _s of each pixel;

Step B, using the Bayesian model averaging method, performing weighted combination calculations on the first increment and the second increment of each pixel with high spatial resolution obtained in step A on NDVI, to obtain each The NDVI comprehensive increment of the pixel;

Step C, from the _t1 time to the _tp time calculated by the step B, the NDVI comprehensive increment of each pixel is added to the high spatial resolution NDVI value corresponding to the _t1 time, and the _tp time can be obtained High spatial resolution NDVI values for .

Preferably, in the step B, the NDVI comprehensive increment of each high spatial resolution pixel ( _xi , y _j ) is expressed as,

ΔNDVI _com ( _xi ,y _j )=w _s *ΔNDVI _s ( _xi ,y _j )+w _d *ΔNDVI _d ( _xi ,y _j ) (2)

Wherein, in formula (2), w _d and w _s represent the weights of the first increment and the second increment respectively, and △NDVI _com ( _xi , y _j ) represents the weight of the first increment and the weight of the second increment The NDVI composite increment obtained by combining the second increment.

Preferably, the weights w _d and w _s of each pixel in the formula (2) are obtained by solving the Bayesian model averaging method.

Preferably, said step B further includes the step of performing error correction on said NDVI synthesis increment.

Preferably, the step of error correction includes:

Firstly, the overall error of the low spatial resolution pixel (x, y) is calculated as,

In formula (4), m represents the number of high spatial resolution pixels covered by low spatial resolution pixels (x, y); R(x, y) represents the overall low spatial resolution pixels (x, y) △ NDVI ^C (x, y) represents the difference between the NDVI of the low spatial resolution pixel (x, y) at time t _p and time t ₁ ;

The homogeneity of the high spatial resolution pixel is then calculated as follows,

HI(x _i ,y _j )=n _smae /m, (5)

In formula (5), HI( _xi , y _j ) is the homogeneity parameter of an 8*8 window constructed with the high spatial resolution pixel ( _xi , y _j ) as the center; m represents the height of the window The number of spatially resolved pixels, n _same represents the number of pixels with the same type of high spatial resolution as the central pixel ( _xi , _yj ) in this window;

After that, the error is distributed according to the homogeneity of the pixel as follows,

In formula (6), r(x _i , y _j ) is the error that the high spatial resolution pixel (x _i , y _j ) should be allocated to;

Finally, this error is added to the comprehensive increment in formula (2), and the corrected NDVI increment of the high spatial resolution pixel is as follows,

In formula (7), △NDVI _com ^new ( _xi , y _j ) is the corrected value of the finally calculated high spatial resolution pixel ( _xi , y _j ) from time t ₁ to time t _p NDVI composite increment.

Preferably, in the step C, the formula for calculating the high spatial resolution NDVI at time _tp using the revised NDVI comprehensive increment is:

In formula (8), NDVI _p ( _xi , y _j ) is the high spatial resolution NDVI of each pixel ( _xi , y _j ) at time t _p , and NDVI ₁ ( _xi , y _j ) is High spatial resolution NDVI for each pixel ( _xi ,y _j ) at time t ₁ .

The present invention obtains two increments of high spatial resolution data in time by performing mixed pixel decomposition and spatial interpolation on low spatial resolution data and high spatial resolution data, and then adopts the Bayesian model averaging method to The weight combination calculation of the two kinds of increments is carried out to obtain the optimized comprehensive increment, which improves the accuracy of the increment calculation, makes full use of the spatio-temporal relationship between geographic data, and improves the fusion accuracy.

Description of drawings

The following drawings are only intended to illustrate and explain the present invention schematically, and do not limit the scope of the present invention. in,

Figure 1 shows a schematic diagram of the principle of the first increment according to a specific embodiment of the present invention;

Figures 2a-2f show the land cover change area and the corresponding Landsat NDVI data and MODISNDVI data respectively; the situation of a catchment area represented therein is: Figure 2a shows the Landsat NDVI on November 26, 2004, and Figure 2b Landsat NDVI on December 12, 2004; Figure 2d shows November 26, 2004 and Figure 2e shows MODIS NDVI on December 12, 2004; Figure 2c shows November 2004 For the Landsat image on the 26th, Figure 2f shows the corresponding unsupervised classification results.

Figures 3a-3e respectively show the fusion results of several different data fusion methods in land cover change areas; Figures 3a-3d respectively represent Landsat NDVI on December 12, 2004 after NDVI-LMGM fusion, STARFM fusion, The result of FSDAF fusion and the method fusion of the present invention, Fig. 3e shows the real value.

detailed description

In order to have a clearer understanding of the technical features, purposes and effects of the present invention, the specific implementation manners of the present invention will now be described with reference to the accompanying drawings. But those skilled in the art should know that the following examples are not the sole limitation to the technical solution of the present invention, and any equivalent transformation or modification made under the spirit of the technical solution of the present invention should be considered as belonging to the protection of the present invention scope.

The principle of the method for producing high spatio-temporal resolution NDVI based on spatio-temporal weighting provided by the present invention will be described in detail below. Of course, it is the same as the three types of methods for fusing Landsat NDVI data and MODIS NDVI data described in the background section. The present invention is also a method for fusing remote sensing data from the above-mentioned different data sources, that is, the method for producing NDVI with high spatial-temporal resolution based on time-space weighting of the present invention is used to combine Landsat NDVI data with high spatial resolution and low temporal resolution Fused with MODIS NDVI data at low spatial resolution and high temporal resolution, NDVI data with high temporal and spatial resolution at prediction time _tp are obtained from reference time _t1 .

Since the satellites acquiring Landsat NDVI data and MODIS NDVI data have similar orbital parameters and satellite transit time intervals of less than 30 minutes (this interval is negligible for a 24-hour day), it can theoretically be assumed that both data The NDVI data at the same location and the same reference time t ₁ only have differences in spatial resolution, and the rest can be ignored. Among them, the spatial resolution of Landsat NDVI data is 30m*30m, and the temporal resolution is 16 days; the spatial resolution of MODIS NDVI data is 250m*250m, and the temporal resolution is 1 day; The maximum spatial resolution of the NDVI data at a high rate can be equal to the spatial resolution of Landsat NDVI data (30m*30m), and its maximum temporal resolution can be equal to the temporal resolution of MODIS NDVI data (1 day). That is to say, through the known MODIS NDVI data with a spatial resolution of 250m*250m and a time resolution of 1 day, and the Landsat NDVI data with a spatial resolution of 30m*30m and a time resolution of 16 days, the same reference time t ₁ is the starting point, and the method of the present invention can predict and obtain NDVI data with high temporal and spatial resolution at the prediction time t _p of any one day in the next 15 days.

The specific steps of the method for producing NDVI with high spatiotemporal resolution based on spatiotemporal weighting of the present invention will be described in detail below.

Step A: Starting from the same reference time _t1 , using MODIS NDVI data and Landsat NDVI data, respectively obtain the first increment △NDVI _d and the second increment △NDVI of each pixel on NDVI at the prediction time t _{p s} .

Wherein, the first increment ΔNDVI _d refers to, from the reference time t ₁ to the predicted time t _p , assuming that the NDVI increment of each pixel of MODIS NDVI data is equal to the difference of Landsat NDVI data under the same area The average value of the NDVI increment of the pixel of the ground object category is calculated by using the MODIS NDVI data and Landsat NDVI data at time _t1 , and the MODIS NDVI data at time _tp to obtain the Landsat NDVI data from time _tp to time _t1 The NDVI increment of each pixel of , which is the first increment ΔNDVI _d .

Specifically, the low-spatial-resolution MODISNDVI data and the high-spatial-resolution Landsat NDVI data at the reference time t ₁ (actually there is a time interval of 30 minutes) only have a difference in spatial resolution for the same location. Aligning the positions of the two data is convenient for calculation, and the MODIS NDVI data with low spatial resolution can be resampled from 250m*250m to 240*240m, so that one MODIS pixel can completely cover 8*8=64 Landsat pixels .

As shown in Figure 1, what it shows is the schematic diagram of the principle of the first increment according to a specific embodiment of the present invention, and what the top of wherein the time axis represents is the MODIS NDVI data of low spatial resolution, and what represent below the time axis is Landsat NDVI data with high spatial resolution. Starting from the same reference time _t1 , using the MODIS NDVI data and Landsat NDVI data at time _t1 and the MODIS NDVI data at time _tp , the high space at time _tp that did not exist originally indicated by the dotted line at the lower right of the time axis can be obtained Resolution Landsat NDVI data.

At this time, it is assumed that from the reference time t ₁ to the prediction time t _p , the NDVI increment of each pixel of the MODIS NDVI data is equal to the average of the NDVI increments of the pixels of different object categories in the Landsat NDVI data under the same area value. That is, the _NDVI increment of _{one MODIS} pixel represented by the heart shape above the time axis in _Fig . The MODIS NDVI data with spatial resolution were directly subtracted.

The area corresponding to 1 MODIS pixel above the time axis is equal to the area of 64 Landsat pixels below the time axis, and the NDVI increments of these 64 Landsat pixels are averaged according to their areas, which can be regarded as equal to time NDVI increment of 1 MODIS cell of the same area above the axis. The calculation process is to calculate the abundance of each surface object type in each pixel of MODIS NDVI data with low spatial resolution according to the object classification of LandsatNDVI data with high spatial resolution. Classification data can come from existing classification products, such as GlobeLand30, or it can be obtained by classifying satellite image data of Landsat NDVI data with high spatial resolution. According to the object classification results of LandsatNDVI data with high spatial resolution, the area proportion of different object types contained in each pixel of MODIS NDVI data with low spatial resolution is calculated, that is, the abundance. Therefore, by adding the product of the NDVI increment of each surface object type and its abundance, the average value of the NDVI increment of all Landsat pixels (64 pixels) under the area can be obtained, which is exactly equal to NDVI increment of MODIS pixel (1 pixel) with the same area from t ₁ to t _p time.

For example, the 64 Landsat pixels below the time axis in Figure 1 include three types of objects, then the 24 pixels of the first type of objects have an NDVI increment, and the 20 pixels of the second type of objects have Another NDVI increment, 20 pixels of the third type of ground object type has another NDVI increment, these three NDVI increments, relative to each pixel, are its time from t ₁ to t _p The first increment ΔNDVI _d . In order to calculate these three increments, at least three equations need to be constructed, assuming that the 8 cells closest to the heart-shaped marker cell above the time axis (that is, the 8 cells surrounding the heart-shaped marker cell) The NDVI increment of is the same as the NDVI increment of the heart-shaped marker pixel, then a 3*3 local window can be constructed centered on each low spatial resolution pixel. Assume that the NDVI increment of each object type in this local window is the same from time _t1 to time _tp . According to the NDVI increments of the 9 low spatial resolution pixels in this local window, the following linear equations can be constructed, and the formulas are as follows,

In formula (1), △NDVI ^C (x, y) represents the NDVI increment of the low-spatial resolution pixel (x, y) from t ₁ to time t _p , and f _l (x, y) represents the type l Abundance in low spatial resolution pixel (x, y), △ NDVI _l (x, y) represents the NDVI increment of type l in a 3*3 local window.

△NDVI _l (x, y) can be obtained by solving the least squares according to formula (1). For the low spatial resolution pixel (x, y) in the local window shown in the figure, the time increment of the high resolution pixel ( _xi , y _j ) under its coverage is △NDVI _d (x _i , y _j ), this increment is called the first increment, and the subscript (i, j) indicates the (i, j)th high spatial resolution pixel covered by the low spatial resolution pixel. It has been assumed that in this local window, the NDVI increment of each object type is fixed. Therefore, for a high spatial resolution pixel ( _xi , y _j ), its time increment △NDVI _d (x _j , y _j ) is equal to the high spatial resolution pixel (x _j , y _j ) is the NDVI increment of the feature type to which it belongs. If the high spatial resolution pixel ( _xi , y _j ) belongs to class l, then △NDVI _d ( _xi , y _j ) is equal to △NDVI _l (x, y).

For the first increment △NDVI _d , since the type of ground objects at time _tp may be quite different from that at time _t1 , for example, soil erosion, deforestation, fire, flood, etc. The feature type changes. Therefore, the abundance f _l (x, y) of each object type in formula (1) may be different at two times, which will lead to the NDVI increment of the object type △ NDVI _l (x, y ) calculation error. And it will also make mistakes when assigning the NDVI increment △NDVI _l (x, y) of each surface object type to the high spatial resolution pixel at time t _p according to the type.

In order to avoid the influence of such unfavorable factors, the present invention further proposes the idea of using the second increment ΔNDVI _s to correct the first increment.

Wherein, the second increment △NDVI _s refers to using the thin plate spline (Thin Plate Spline, TPS) method tool integrated in the existing ENVI/IDL software to directly convert the low space between the time t ₁ and the time t _p The resolution MODIS pixels are interpolated to the high spatial resolution Landsat pixels respectively (this interpolation method tool has been integrated in the ENVI/IDL software, which is a known method and can be completed by calling a function). Then, the interpolation result at time t _p is subtracted from the interpolation result at time t ₁ , and the result is the second increment ΔNDVI _s of each pixel. Since the second increment △NDVI _s ignores the difference of object types, it can not only provide spatial information, but also not be disturbed by the types of objects, so the first and second increments are fused through the following steps, The defects and adverse effects of the first increment △NDVI _d can be eliminated.

That is, one MODIS pixel of the heart-shaped mark shown in Figure 1 corresponds to the 64 Landsat pixels below, and the first increment △NDVI _d of each pixel is related to the type of ground object. object type, the first increment △NDVI _d of the 24 pixels under the first type of object type is the same, and the first increment △NDVI _d of the 20 pixels under the second type of object type is also Similarly, the first increment △NDVI _d of the 20 pixels under the third type of ground object is also the same. Of course, the values of these three first increments △NDVI _d are calculated separately, and are mutually independant. The above-mentioned first increment △NDVI _d is obtained under the assumption that the object category of the same pixel has not changed, which obviously has certain problems.

The acquisition of the second increment △ NDVI _s is based on the method of thin plate splines, directly interpolating MODIS pixels with low spatial resolution at time _t1 and time _tp to obtain Landsat pixels with high spatial resolution at time _tp (The principle of the thin-plate spline method is to find multiple matching points in two low-spatial-resolution MODIS images, and apply the thin-plate spline function to deform these multiple points to the corresponding positions, and at the same time give the interpolation of the entire image ), the NDVI value of each pixel (for example, corresponding to 64 Landsat pixels in Figure 1) of the Landsat image with high spatial resolution at time t _p obtained by interpolation and the known Landsat image at time t ₁ Subtract the NDVI value of each pixel to obtain the second increment △NDVI _s of each pixel from t ₁ to t _p , that is, the second increment △NDVI _s has nothing to do with the type of ground objects, and each image Each element has an independent second increment ΔNDVI _s . Of course, since Landsat pixels with high spatial resolution can be directly interpolated from MODIS pixels with low spatial resolution, it is worthless to obtain the second increment △ NDVI _s by separate calculation. Therefore, the second increment △ NDVI s is proposed here The purpose of △NDVI _s is for the next step of fusion, so as to eliminate the problems caused by using the first increment △NDVI _d alone to obtain NDVI with high spatio-temporal resolution.

In step B, the Bayesian model averaging method is used to combine the weights of the first increment and the second increment of each pixel with high spatial resolution on NDVI to obtain the comprehensive NDVI increment of each pixel.

Here two weights are used to adjust the first increment △NDVI _d ( _xi , y _j ) and the second increment △NDVI _s ( _xi , y _j ) corresponding to each pixel ( _xi , y _j ) Combine weights. The NDVI comprehensive increment of each high spatial resolution pixel (x _i , y _j ) is obtained as,

ΔNDVI _com ( _xi ,y _j )=w _s *ΔNDVI _s ( _xi ,y _j )+w _d *ΔNDVI _d ( _xi ,y _j ) (2)

In formula (2), w _d and w _s represent the weights of the first increment and the second increment respectively, these two weights do not vary with the pixel, and the whole image has only one w _d and one w _s . ΔNDVI _com ( _xi , y _j ) represents the NDVI comprehensive increment obtained by combining the first increment and the second increment.

Next, the Bayesian Model Averaging (BMA) method will be used to solve the two weights for each pixel. The BMA method is a more classic multi-model synthesis method, which calculates the final weight of each model by optimizing the likelihood function. For details, see references Raftery, A.E., Gneiting, T., Balabdaoui, F., & Polakowski, M. (2005), "Ensemble Forecast Correction for Weather Using Bayesian Model Averaging", Monthly Weather Review, 133 , 1155-1174 (Raftery, A.E., Gneiting, T., Balabdaoui, F., & Polakowski, M. (2005). Using Bayesian Model Averaging to Calibrate Forecast Ensembles. Monthly Weather Review, 133, 1155-1174). The maximum likelihood function for the first increment and the second increment in the present invention is expressed as follows,

In formula (3), L represents the maximum likelihood value, w _s and w _d represent the weights of the two increments respectively, △NDVI _d,i represents the first increment of the i-th pixel, △NDVI _s,i Indicates the second increment of the i-th pixel, and △NDVI _i represents the real increment of the i-th pixel. g _s and g _d respectively represent the probability density functions of the two incremental sets, which are considered to be normal distribution, and N represents the number of training samples. In the present invention, since there is no real NDVI increment data of high spatial resolution, the real increment of the low-resolution pixel at time _t1 and time _tp is used as the training sample, and the first increment of high spatial resolution and The second increments are all 8*8 synthetically upscaled to low spatial resolution. In this way, the real NDVI increment of the pixel at the low spatial resolution scale, the first increment of the low spatial resolution and the second increment of the low spatial resolution can be obtained. Calculate the mean square error (Mean Square Error, MSE) between the first increment of low spatial resolution and the real NDVI increment of low spatial resolution. This mean square error is the variance of the probability density function g _d , similarly The variance of the probability density function g _s can be obtained. g _s (△NDVI _i |△NDVI _s,i ) means the probability that the probability density function g _s can produce △NDVI _i with its mean squared error as the variance and △NDVI _s,i as the mean. Similarly, g _d (△NDVI _i |△NDVI _d,i ) represents the probability that the probability density function g _d can produce △NDVI _i with its average squared error as the variance and △NDVI _d,i as the mean. By maximizing the likelihood function (3), the weights w _s and w _d can be obtained.

The value shows that although the NDVI comprehensive increment △ NDVI _com ( _xi , y _j ) can be relatively close to the real NDVI increment of the high spatial resolution pixel from time t ₁ to time t _p , it is still May contain errors, so in order to obtain more accurate results, in a preferred embodiment, further error correction can also be carried out to the obtained NDVI comprehensive increment, and its correction steps are:

Firstly, the overall error of the low spatial resolution pixel (x, y) is calculated as,

In formula (4), m represents the number of high spatial resolution pixels covered by low spatial resolution pixels (x, y), which is 64 in Fig. 1; R(x, y) represents low spatial resolution The overall error of the resolution pixel (x, y); △ NDVI ^C (x, y) represents the difference between the NDVI of the low spatial resolution pixel (x, y) at time t _p and time t ₁ .

Here, it is necessary to distribute the error R(x, y) contained in the pixel scale of low spatial resolution to the pixel (x _i , y _j ) of high spatial resolution covered by it. Generally speaking, areas with relatively high homogeneity on the image (such as a large body of water or a large grassland) tend to have less error in prediction, while areas with low homogeneity (such as sparse forest land interlaced with trees and bare soil) tend to The prediction error is large. Therefore, the present invention distributes the error R(x,y) according to the homogeneity ( _xi , _yj ) of the high spatial resolution pixel.

That is, the homogeneity of the high spatial resolution pixel is then calculated as follows,

HI(x _i ,y _j )=n _smae /m, (5)

In formula (5), HI( _xi , y _j ) is a homogeneity parameter of an 8*8 window constructed with the high spatial resolution pixel ( _xi , y _j ) as the center. m represents the number of high spatial resolution pixels in this window, 64 here, and n _same represents the number of high spatial resolution pixels of the same type as the central pixel ( _xi , y _j ) in this window.

After that, the error is distributed according to the homogeneity of the pixel as follows,

In formula (6), r( _xi , y _j ) is the error to which the high spatial resolution pixel ( _xi , y _j ) should be assigned.

Finally, this error is added to the comprehensive increment in formula (2), and the corrected NDVI increment of the high spatial resolution pixel is as follows,

In formula (7), △ NDVI _com ^new ( _xi , _{y j )} is the _{corrected NDVI} comprehensive increment.

Step C, calculating the high spatial resolution NDVI at time _tp ;

From time t ₁ to time t _p calculated in the previous steps, the NDVI comprehensive increment of each high spatial resolution pixel ( _xi , y _j ) will be added to the high space of time t ₁ below The high spatial resolution NDVI data of time t _p can be obtained on the resolution NDVI. Calculated as follows,

In formula (8), NDVI _p ( _xi , y _j ) is the high spatial resolution NDVI at time t _p , and NDVI ₁ ( _xi , y _j ) is the high spatial resolution NDVI at time t ₁ .

As mentioned above, after the processing of steps A, B and C, all the pixels with high spatial resolution at the prediction time _tp can obtain more accurate fusion results through the prediction process of spatiotemporal weighting. Compared with the existing data fusion method, the present invention uses the information contained in the low spatial resolution data more fully through the BMA (Bayesian model average) method, and can better deal with the changes caused by land cover type changes. The problem.

In order to better illustrate the technical effects of the present invention, the method in the present invention is compared with the three types of data fusion methods mentioned in the background technology section, namely for STARFM (Space and Time Adaptive Reflective Fusion Model), NDVI- LMGM (NDVI-Linear Mixed Growth Model) and FSDAF (Flexible Spatio-Temporal Data Fusion Method) are used to test time-series data and single-scene data respectively. In the test, NDVI is used as the fusion object, the low spatial resolution data is MODIS data, and the high spatial resolution data is Landsat data. Since the transit time of the two data acquisitions is the same, and the orbital parameters are close, they are more suitable for fusion processing.

Figures 2a-2f show the land cover change area and the corresponding Landsat NDVI data and MODISNDVI data; they represent the situation of a catchment area, and the specific test data is located in the Gwydir lower watershed research area in northern New South Wales, Australia (Lower Gwydir Catchment site), the latitude and longitude are 29°07′S and 149°04′E respectively. This area is a typical crop planting area, and a relatively large flood event occurred in December 2004, a large number of farmlands were submerged, and relatively obvious land cover changes occurred. In order to verify the detection effect of the present invention for land cover change, from the official website of the Commonwealth Scientific and Industrial Research Organization (CSIRO): https://data.csiro.au/dap/landingpage? execution=e2s2&_eventId=viewDescription, downloaded the Landsat TM images of November 26, 2004 before the flood and December 12, 2004 after the flood, clipped them, and calculated the corresponding NDVI (Fig. 2a and 2b). According to the clipped image on November 26, 2004 (Fig. 2c), the non-supervised classification based on isodata was carried out, and the classification result was obtained (Fig. 2f). In order to avoid systematic errors between different sensors, as well as errors caused by geometric correction and atmospheric correction, the present invention uses Landsat data to perform 8*8 upscaling synthesis to obtain MODIS data (Figure 2d and 2e).

Figures 3a-3e respectively show the fusion results of several different data fusion methods in the land cover change area. Taking the Landsat NDVI data on November 26, 2004 as the reference data, the data on December 12, 2004 were fused . The fusion results of the four methods (Fig. 3a, Fig. 3b, Fig. 3c, and Fig. 3d) were compared with the real NDVI on December 12, 2004 (Fig. 3e). It can be seen from the whole of the partially enlarged picture that the fusion result of the NDVI-LMGM method has obvious plaque problems, while the smoothness of the fusion result of the method of the present invention is relatively good. It can be seen from the smaller boxes in Figure 3a-3e; STARFM fusion is not enough for the detection of land cover change, the black area in the smaller box has not been identified, and the method of the present invention can be well It can be seen from the larger box that the result of FSDAF fusion is relatively greatly affected by classification, and the fusion result clearly retains the boundary of the water body before the flood, making the area covered by water The transition is not smooth, but the method of the present invention is relatively less affected by classification, has a stronger ability to identify land cover, and the fusion result is closer to the real value.

To sum up, the present invention synthesizes the increment of mixed pixel decomposition and the increment of spatial interpolation based on the Bayesian model averaging method, and optimizes the calculation process of the increment. The spatio-temporal data are synthesized, which effectively improves the accuracy of data fusion. The invention makes full use of the time-space correlation between geographic data, and provides an effective method for producing vegetation indices with high time-space resolution.

Those skilled in the art should understand that although the present invention is described in terms of multiple embodiments, not each embodiment only includes an independent technical solution. The description in the description is only for the sake of clarity, and those skilled in the art should understand the description as a whole, and understand the present invention by considering the technical solutions involved in each embodiment as being able to be combined with each other to form different embodiments scope of protection.

The above descriptions are only illustrative specific implementations of the present invention, and are not intended to limit the scope of the present invention. Any equivalent changes, modifications and combinations made by those skilled in the art without departing from the concept and principle of the present invention shall fall within the protection scope of the present invention.

Claims (6)

1. A method for producing high space-time resolution NDVI based on space-time weighting is used for fusing Landsat NDVI data with high space resolution and low time resolution with MODIS NDVI data with low space resolution and high time resolution from a reference time t₁Obtaining a predicted time t_pThe high spatial and temporal resolution NDVI data of (1), the method comprising the steps of:

step A: at the same reference time t₁As a starting point, obtaining the predicted time t by using MODIS NDVI data and Landsat NDVI data respectively_pEach pixel ofFirst increment △ NDVI on NDVI_dAnd a second increment △ NDVI_s(ii) a Wherein,

the first increment △ NDVI_dIs to mean from the reference time t₁To the predicted time t_pAssuming that the NDVI increment of each pixel of the MODIS NDVI data is equal to the average value of the NDVI increments of the pixels of different ground object types of the Landsat NDVI data under the same area, t is used₁MODIS NDVI and Landsat NDVI of time, and t_pObtaining t by calculating MODISNDVI data of time_pTime to t₁The NDVI increment of each pixel of the Landsat NDVI data of the time is the first increment △ NDVI_d；

The second increment △ NDVI_sThe method is characterized in that t is directly processed by adopting a thin plate spline method tool integrated in the existing ENVI/IDL software₁Time and t_pRespectively interpolating the time low-spatial resolution MODIS pixels to the high-spatial resolution Landsat pixels, and then using t_pInterpolation result of time and t₁The interpolation results of the time are subtracted, and the result is the second increment △ NDVI of each pixel_s；

B, performing weight combination calculation on the first increment and the second increment of each pixel with high spatial resolution on the NDVI by using a Bayes model average method to obtain an NDVI comprehensive increment of each pixel;

step C, t obtained by calculating in step B₁Time to t_pTime, the NDVI integrated increment for each pel, added to t₁On the high spatial resolution NDVI value corresponding to time, t can be obtained_pHigh spatial resolution NDVI values in time.

2. The method of claim 1, wherein in step B, each high spatial resolution pel (x) is_i，y_j) Is formulated as,

ΔNDVI_com(x_i,y_j)＝w_s*ΔNDVI_s(x_i,y_j)+w_d*ΔNDVI_d(x_i,y_j) (2)

wherein, in the formula (2), w_dAnd w_sWeights representing said first increment and said second increment, respectively, △ NDVI_com(x_i，y_j) Representing the NDVI composite delta resulting from the combination of the first delta and the second delta.

3. A method according to claim 2, characterized in that the weight w for each picture element of said formula (2)_dAnd w_sAnd solving by a Bayesian model averaging method.

4. The method of any one of claims 1-3, wherein said step B further comprises the step of error correcting said NDVI composite delta.

5. The method of claim 4, wherein the step of error correction comprises:

firstly, the error of the whole low-spatial resolution pixel (x, y) is calculated as,

<mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&amp;Delta;NDVI</mi> <mi>C</mi> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mn>1</mn> <mi>m</mi> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <msub> <mi>&amp;Delta;NDVI</mi> <mrow> <mi>c</mi> <mi>o</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

in equation (4), m represents the number of high spatial resolution pixels covered by the low spatial resolution pixel (x, y), R (x, y) represents the error of the low spatial resolution pixel (x, y) as a whole, △ NDVI^C(x, y) denotes t_pTime and t₁Difference in NDVI for temporal low spatial resolution pixel (x, y);

the homogeneity of the high spatial resolution picture elements is then calculated as follows,

HI(x_i,y_j)＝n_smae/m， (5)

in the formula (5), HI (x)_i，y_j) Is a high spatial resolution picture element (x)_i，y_j) As a center, homogeneity parameter of the constructed 8 x 8 windows; m represents the number of high spatial resolution picture elements in this window, n_sameIndicating that in this window, the pixel is centered on (x)_i，y_j) A number of high spatial resolution picture elements having the same type;

the error is then assigned according to the homogeneity of the picture elements as follows,

<mrow> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mi>m</mi> <mo>*</mo> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>*</mo> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mi>I</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>/</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <mo>&amp;lsqb;</mo> <mn>1</mn> <mo>-</mo> <mi>H</mi> <mi>I</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

in the formula (6), r (x)_i，y_j) I.e. high spatial resolution picture element (x)_i，y_j) The error that should be assigned;

finally, this error is added to the integrated increment in equation (2) to obtain the NDVI increment for the corrected high spatial resolution pixel as follows,

<mrow> <msubsup> <mi>&amp;Delta;NDVI</mi> <mrow> <mi>c</mi> <mi>o</mi> <mi>m</mi> </mrow> <mrow> <mi>n</mi> <mi>e</mi> <mi>w</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&amp;Delta;NDVI</mi> <mrow> <mi>c</mi> <mi>o</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <mi>r</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

in equation (7), △ NDVI_com ^new(x_i，y_j) Namely the high spatial resolution pixel (x) obtained by final calculation_i，y_j) From t₁Time to t_pA corrected integrated increment of the NDVI in time.

6. The method of claim 5, wherein in step C, t is calculated using the modified NDVI synthetic delta_pThe formula for the temporal high spatial resolution NDVI is:

<mrow> <msub> <mi>NDVI</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>NDVI</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msubsup> <mi>&amp;Delta;NDVI</mi> <mrow> <mi>c</mi> <mi>o</mi> <mi>m</mi> </mrow> <mrow> <mi>n</mi> <mi>e</mi> <mi>w</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

in formula (8), NDVI_p(x_i，y_j) Is t_pEach picture element of time (x)_i，y_j) High spatial resolution NDVI, NDVI₁(x_i，y_j) Is t₁Each picture element of time (x)_i，y_j) High spatial resolution NDVI.