CN107103584A

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

Info

Publication number: CN107103584A
Application number: CN201710231237.4A
Authority: CN
Inventors: 曹鑫; 刘梦; 陈晋; 陈学泓; 杨林清
Original assignee: Beijing Normal University
Current assignee: Beijing Normal University
Priority date: 2017-04-11
Filing date: 2017-04-11
Publication date: 2017-08-29
Anticipated expiration: 2037-04-11
Also published as: CN107103584B

Abstract

The invention provides a kind of production high-spatial and temporal resolution NDVI weighted based on space-time method, including：To t₁The NDVI data of the low spatial resolution of time and the NDVI data and t of high spatial resolution_pThe low spatial resolution NDVI data of time obtain the two kinds of NDVI increments of high spatial resolution data in time by Decomposition of Mixed Pixels and space interpolation；Described two NDVI increments are carried out with weight combination using Bayesian model averaging method to calculate, the NDVI synthesis increments of each pixel are obtained；NDVI synthesis increments are finally added to t₁Required high-spatial and temporal resolution NDVI can be obtained on the NDVI of the high spatial resolution of time.The present invention carries out Bayes's synthesis to two kinds of NDVI increments, takes full advantage of the space-time relationship between geodata, and the vegetation index for production high-spatial and temporal resolution provides a kind of effective method.

Description

Method for producing high-space-time resolution NDVI (normalized difference vector magnitude) based on space-time weighting

Technical Field

The invention relates to a method for fusing multi-source remote sensing data, which is mainly used for producing a normalized vegetation index (NDVI) with high space-time resolution and belongs to the field of space-time data fusion.

Background

The Normalized Difference Vegetation Index (abbreviated NDVI) is a widely used Vegetation Index. The NDVI value can reflect the growth condition of the vegetation and is directly related to the physiological parameters of the vegetation, such as chlorophyll content, water content, leaf area index and the like.

CN 102831310B discloses a method for constructing time series data with high spatial resolution, and its background section details the relevant knowledge of NDVI data, wherein two existing NDVI data are mentioned: LandsatNDVI data and MODIS NDVI data. The spatial resolution of NDVI values in Landsat NDVI data was 30m × 30m and the temporal resolution was 16 days, while the maximum spatial resolution of NDVI values in MODIS NDVI data was 250m × 250m and the temporal resolution was 1 day. That is, the Landsat NDVI data has a higher spatial resolution and a lower temporal resolution, while the MODIS NDVI data has a lower spatial resolution and a higher temporal resolution.

In order to obtain NDVI data with high spatial resolution and high temporal resolution simultaneously, the prior art develops various methods for fusing Landsat NDVI data and MODIS NDVI data to obtain NDVI data with high spatial resolution, which mainly include three types: (1) a weight function based approach; (2) a method based on mixed pixel decomposition; (3) a hybrid model based approach.

Among them, (1) the weight function-based method is typically represented by STARFM (spatial and Temporal Adaptive reflection Fusion Model), see Gao, F., mask, j., Schwaller, m., & Hall, F. (2006), "Fusion of land and MODIS surface reflectivities: prediction of daily Landsat surface reflectance ", the institute of Electrical and electronics Engineers, science and remote sensing, 44, 2207-. (Gao, F., Masek, J., Schwaller, M., & Hall, F. (2006.) the blending of the bonded modified surface reflection. IEEE Transactions On Geoscience and remove Sensing, 44, 2207- & 2218). The method focuses on the search of similar pixels and the determination of weights thereof, determines the similar pixels and calculates corresponding weights according to the similarity of the current pixel and the surrounding pixels on the characteristics of time, space, spectrum and the like, and takes the average increment of the similar pixels as the increment of the current pixel. Such methods tend to work well in areas of relatively high purity.

(2) The Method based on mixed pixel decomposition is represented by NDVI-LMGM (NDVI-linear mixed Growth Model), see Rao, Y., Zhu, X., Chen, J., & Wang, J. (2015), "a Method for Producing High Spatial and Temporal Resolution NDVI Time series Data by using Multi-Temporal MODIS NDVI Data and Landsat TM/ETM + images" remote sensing, 7, 7865 + 7891(Rao, Y., Zhu, X., & Chen, J., & Wang, J. (2015). The method calculates the increment of each image element in time through mixed image element decomposition so as to obtain NDVI data of the predicted time, and specific contents of the method can also be referred to in the aforementioned prior art CN 102831310B, which are substantially the same. Such methods can also be effective in areas with a relatively high degree of heterogeneity, but this is not necessarily satisfactory assuming that the land cover conditions do not change.

(3) The method based on the hybrid model is represented by FSDAF (flexible spatiotemporal Data Fusion method), see Zhu, x., helper, e.h., Gao, f., Liu, d., Chen, j., & Lefsky, M.A. (2016), "a method for flexibly spatiotemporal Fusion of satellite images of different resolutions", environmental remote Sensing, 172, 165-. The method combines a plurality of models, such as a weight function model, a mixed pixel decomposition model, a spatial interpolation model and the like, and calculates to obtain the increment of each pixel on the high-resolution image. Such methods tend to be relatively complex, but relatively accurate, and can capture land cover variations to some extent.

The three methods have good effect in respective applicable fields. However, the method of interpolation from two-stage low-spatial-resolution data to high-spatial-resolution data is not considered, and the increment between two-stage interpolation results can reflect the increment of pixel dimension to a certain extent and also contains information of land coverage change. The invention integrates the increment of mixed pixel decomposition and the increment of spatial interpolation based on the Bayesian model average method, and can effectively improve the precision of data fusion. Provides a new way for producing NDVI with high space-time resolution.

Disclosure of Invention

The technical problem underlying the present invention is to provide a method for producing high spatial and temporal resolution NDVI based on spatial and temporal weighting to reduce or avoid the aforementioned problems.

In order to solve the technical problem, the invention provides a method for producing high-space-time resolution NDVI based on space-time weighting, which 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_pIs at a first 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 MODIS NDVI 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.

Preferably, in step B, each high spatial resolution image element (x)_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.

Preferably, said weight w of each pixel of said formula (2)_dAnd w_sAnd solving by a Bayesian model averaging method.

Preferably, the step B further comprises the step of error correcting the NDVI composite delta.

Preferably, the error correction step includes:

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

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,

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,

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.

Preferably, in said step CCalculating t using said modified NDVI synthetic delta_pThe formula for the temporal high spatial resolution NDVI is:

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.

According to the method, the two time increments of the high-spatial-resolution data are obtained by performing mixed pixel decomposition and spatial interpolation on the low-spatial-resolution data and the high-spatial-resolution data, and then the two increments are subjected to weight combination calculation by adopting a Bayesian model averaging method to obtain the optimized comprehensive increment, so that the accuracy of increment calculation is improved, the spatio-temporal relation among the geographic data is fully utilized, and the fusion precision is improved.

Drawings

The drawings are only for purposes of illustrating and explaining the present invention and are not to be construed as limiting the scope of the present invention. Wherein,

FIG. 1 is a schematic diagram illustrating a first increment according to an embodiment of the present invention;

FIGS. 2a-2f show a land cover change area and corresponding Landsat NDVI data and MODISNDVI data, respectively; the conditions of a certain catchment basin area shown therein are respectively: FIG. 2a shows Landsat NDVI at 26/11/2004, and FIG. 2b shows Landsat NDVI at 12/2004; respectively simulating MODIS NDVI represented by 26 days at 11 months in 2004 and 12 days at 12 months in 2004 in FIG. 2d and 2e from the images; fig. 2c shows Landsat image at 26/11/2004, and fig. 2f shows the corresponding unsupervised classification result.

FIGS. 3a-3e show the fusion results of several different data fusion methods in a land cover change area, respectively; in which FIGS. 3a-3d show the results of NDVI-LMGM fusion, STARFM fusion, FSDAF fusion, and the method of the present invention, respectively, for Landsat NDVI at 12 months and 12 days 2004, and FIG. 3e shows the true values.

Detailed Description

In order to more clearly understand the technical features, objects, and effects of the present invention, embodiments of the present invention will now be described with reference to the accompanying drawings. However, it should be understood by those skilled in the art that the following examples are not intended to limit the scope of the present invention, and any equivalent changes or modifications made within the spirit of the present invention should be considered as falling within the scope of the present invention.

The principle of the method for producing high spatial-temporal resolution NDVI based on spatial-temporal weighting according to the present invention is described in detail below, and it is understood that, like the three types of method 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 of the above-mentioned different data sources, that is, the method for producing high spatial-temporal resolution NDVI based on spatial-temporal weighting according to the present invention is used for fusing Landsat NDVI data with high spatial resolution and low temporal resolution with MODIS NDVI data with low spatial resolution and high temporal resolution from the reference time t₁Obtaining a predicted time t_pHigh spatial and temporal resolution NDVI data.

Since the satellites from which the Landsat NDVI and MODIS NDVI data are obtained have similar orbital parameters and a satellite transit time interval of less than 30 minutes (which is negligible for 24 hours per day), it can be theoretically assumed that the two data are co-located and at the same reference time t₁The NDVI data of (A) only differ in spatial resolutionThe rest can be ignored. Wherein the spatial resolution of the Landsat NDVI data is 30m by 30m, and the time resolution is 16 days; the spatial resolution of the MODIS NDVI data is 250m by 250m, and the temporal resolution is 1 day; the data fusion then produces high spatial and temporal resolution NDVI data with a maximum spatial resolution that may be equal to the spatial resolution of Landsat NDVI data (30m x 30m) and a maximum temporal resolution that may be equal to the temporal resolution of MODIS NDVI data (1 day). That is, the MODIS NDVI data with spatial resolution of 250m and temporal resolution of 1 day and the Landsat NDVI data with spatial resolution of 30m and temporal resolution of 16 days are known, and the same reference time t is used₁As a starting point, the method of the present invention may predict the predicted time t to be obtained for any 1 day 15 days after the prediction_pHigh spatial and temporal resolution NDVI data.

The specific steps of the method of the present invention for producing high space-time resolution NDVI based on space-time weighting are described in detail below.

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_pIs at a first increment △ NDVI on NDVI_dAnd a second increment △ NDVI_s。

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 MODIS NDVI 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。

Specifically, the reference time t₁Low spatial resolution modistvi data and high spatial resolution Landsat NDVI data (actually, 30-minute time intervals)Only the difference in spatial resolution exists for the same site, and in order to align the positions of the two data for convenient calculation during fusion, the MODIS NDVI data with low spatial resolution can be resampled from 250m by 250m to 240m by 240m, so that one MODIS pixel can completely cover 8 x 8 to 64 Landsat pixels.

Referring to fig. 1, a schematic diagram of a first increment according to an embodiment of the present invention is shown, wherein the upper part of the time axis represents the low spatial resolution MODIS NDVI data, and the lower part of the time axis represents the high spatial resolution Landsat NDVI data. At the same reference time t₁As starting point, using t₁MODIS NDVI and Landsat NDVI of time, and t_pThe MODIS NDVI data of time can obtain t which does not exist originally and is indicated by a dotted line at the lower right of a time axis_pTemporally high spatial resolution Landsat NDVI data.

At this time, assume that the reference time t is counted from₁To the predicted time t_pThe NDVI increment of each pel of MODIS NDVI data is equal to the average of the NDVI increments of pels of different terrain categories of landmass NDVI data under the same area. That is, 1 MODIS pel, represented by a heart above the time axis in FIG. 1, starts from t₁To t_pNDVI of time can be increased by adding t_pMODIS NDVI data and t with low spatial resolution of time₁The MODIS NDVI data with low spatial resolution of time are directly obtained by subtraction.

The area of 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 the 64 Landsat pixels are averaged according to the areas, so that the NDVI increment of 1 MODIS pixel above the time axis and having the same area can be regarded as the NDVI increment of the 1 MODIS pixel above the time axis. The calculation process comprises the step of calculating the abundance of each ground object type in the pixel of each MODIS NDVI data with low spatial resolution according to the ground object classification of the LandsatNDVI data with high spatial resolution. The classification data can be derived from existing classification products, such as GlobeLand30, or can be obtained by counting satellite images of Landsat NDVI data with high spatial resolutionObtained by classification. And according to the ground object classification result of the LandsatNDVI data with high spatial resolution, calculating the area proportion, namely the abundance, of different ground object types contained in the pixel of each MODIS NDVI data with low spatial resolution. Therefore, the average value of the NDVI increment of all Landsat image elements (64 image elements) in the area can be obtained by adding the product of the NDVI increment of each ground type and the abundance of the NDVI increment, and the average value is exactly equal to the average value of the NDVI increment of the MODIS image elements (1 image element) in the same area from t₁To t_pNDVI increment of time.

For example, the 64 Landsat pels below the time axis of FIG. 1 include three feature types, then 24 pels of type 1 feature type have one NDVI increment, 20 pels of type 2 feature type have another NDVI increment, and 20 pels of type 3 feature type have yet another NDVI increment, with respect to each pel it is from t₁To t_p△ NDVI as a first increment of time_d. To calculate these three increments, at least three equations need to be constructed, and assuming that the NDVI increment of the 8 pixels nearest the heart marked pixel above the time axis (i.e., the 8 pixels surrounding the heart marked pixel) is the same as the NDVI increment of the heart marked pixel, a 3 × 3 local window can be constructed centering on each low spatial resolution pixel. Suppose that each type of terrain in the local window is at t₁Time to t_pThe NDVI increment over time is the same. From the NDVI increments for the 9 low spatial resolution pixels in this local window, a system of linear equations can be constructed, formulated as follows,

in equation (1), △ NDVI^C(x, y) denotes the number from t₁To t_pNDVI increment, f, of temporal low spatial resolution pixel (x, y)_l(x, y) denotes the abundance of the l-th type in the low spatial resolution picture element (x, y), △ NDVI_l(x, y) denotes a 3 x 3 partial windowNDVI increment of type i.

△ NDVI can be obtained by solving with least squares according to equation (1)_l(x, y). For a low spatial resolution pel (x, y) in the illustrated local window, it covers a high resolution pel (x) under its coverage_i，y_j) The increment in time is △ NDVI_d(x_i，y_j) This increment is called the first increment, and the index (i, j) indicates the (i, j) th high spatial resolution pel under the coverage of the low spatial resolution pel. It has been assumed that the NDVI increment for each terrain type is fixed in this local window. Thus, for high spatial resolution picture elements (x)_i，y_j) In other words, its increment in time is △ NDVI_d(x_j，y_j) Is equal to the high spatial resolution picture element (x)_j，y_j) The NDVI increment of the type of terrain to which it belongs. If the high spatial resolution pel (x)_i，y_j) Belongs to class i, then △ NDVI_d(x_i，y_j) Is equal to △ NDVI_l(x，y)。

△ NDVI for the first increment_dIn terms of t_pType of ground object and t of time₁The types of the ground features of different time periods can be greatly different, for example, soil erosion, forest felling, fire, flood and the like can cause the types of the ground features of different time periods to change. Therefore, the abundance f of each terrain type in formula (1)_l(x, y) may be different at two times, resulting in an increase △ NDVI in the NDVI of the terrain type_l(x, y) is calculated in error, and the NDVI for each terrain type is incremented △ NDVI_l(x, y) is assigned to t by type_pErrors also occur in temporally high spatial resolution pels.

To avoid the effect of such adverse factors, the present invention further proposes to use a second increment △ NDVI_sThe idea of modifying the first increment.

Wherein the second increment △ NDVI_sThat is, the existing ENVI/IDL software centralization is adoptedFormed Thin Plate Spline (TPS) method tool, directly cutting t₁Time and t_pAnd respectively interpolating the time low-spatial-resolution MODIS pixels to the high-spatial-resolution Landsat pixels (the interpolation method tool is integrated in ENVI/IDL software and can be completed by calling a function for a known method). Then, use 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△ NDVI due to the second increment_sIgnoring the difference in ground object type, not only can provide spatial information, but also is not interfered by the ground object type, so the first increment △ NDVI can be eliminated by fusing the first and second increments through the following steps_dDefects and adverse effects.

That is, 1 MODIS pel of the cardioid mark shown in FIG. 1, corresponding to 64 Landsat pels, the first increment △ NDVI of each pel_dWith respect to the feature type, three feature types are shown in the figure, and the first increment △ NDVI of 24 image elements under the feature type 1_dAll the same, first increment △ NDVI of 20 picture elements in type 2 terrain type_dAlso the same, first increment △ NDVI of 20 picture elements under type 3 terrain type_dAnd are all the same, of course, these three first increments △ NDVI_dAre calculated separately and independently of each other, the first increment △ NDVI described above_dIs obtained under the assumption that the ground object class of the same image element is not changed, which obviously has certain problems.

And a second increment △ NDVI_sIs obtained by directly using t according to a thin-plate spline method tool₁Time and t_pT is obtained by the interpolation of MODIS pixel elements with low spatial resolution of time_pThe method for obtaining the time high-spatial-resolution Landsat pixel (thin plate spline is based on the principle that a plurality of matching points are found out in two low-spatial-resolution MODIS images, the points can be deformed to corresponding positions by applying a thin plate spline function, and meanwhile, the interpolation of the whole image is given), and t obtained by interpolation is used_pThe NDVI value of each pel of the temporally high spatial resolution Landsat image (e.g., corresponding to 64 Landsat pels in fig. 1) is compared to the known t₁The NDVI values of each pixel of the time Landsat image are subtracted, and the difference between t and t of each pixel is obtained₁To t_p△ NDVI A second increment of time_sI.e. the second increment △ NDVI_sIndependent of the type of terrain, each pixel has an independent second increment △ NDVI_sOf course, since the Landsat pixel with high spatial resolution can be directly interpolated from the MODIS pixel with low spatial resolution, the second increment △ NDVI is obtained by independent calculation_sIs of no value, therefore, a second delta △ NDVI is proposed here_sThe purpose of (2) is to eliminate the need to use the first increment △ NDVI alone for further fusion_dProblems that arise when high spatial and temporal resolution NDVI is obtained.

And step 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.

Here two weights are used to apply to each pel (x)_i，y_j) Corresponding first increment △ NDVI_d(x_i，y_j) And a second increment △ NDVI_s(x_i，y_j) And performing weight combination. Obtaining each high spatial resolution pixel (x)_i，y_j) The integrated increment of NDVI of (a) is,

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

in the formula (2), w_dAnd w_sWeights representing the first increment and the second increment respectively, the weights being invariant with respect to picture element, the whole image having only one w_dAnd a w_s。△NDVI_com(x_i，y_j) Representing the NDVI composite gain obtained by combining the first increment and the second incrementAmount of the compound (A).

These two weights for each pel will need to be solved below using a Bayesian Model Averaging (BMA) method. The BMA method is a relatively classical multi-model synthesis method, which calculates the final weights of each model by optimizing a likelihood function. Details can be found in references rafter, a.e., gneting, t., Balabdaoui, f, & Polakowski, M. (2005), "correction of collective forecast of Weather Using Bayesian Model Averaging", monthly Weather Review, 133, 1155-. The maximum likelihood function for the first and second increments of the present invention is expressed as follows,

in the formula (3), L represents a maximum likelihood value, w_sAnd w_dRepresenting the weight of two increments, △ NDVI, respectively_d,iRepresenting a first increment of the ith pixel, △ NDVI_s,iA second increment, △ NDVI, representing the ith pixel_iRepresenting the true increment of the ith pixel element. g_sAnd g_dThe probability density functions of the two increment sets are respectively expressed and considered to be normal distribution, and N represents the number of training samples. In the present invention, t is used because there is no real NDVI incremental data of high spatial resolution₁Time and t_pAnd taking the real increment of the time low-resolution pixel as a training sample, and performing 8-by-8 synthesis upscaling on the first increment and the second increment of the high spatial resolution to the low spatial resolution. Thus, the real NDVI increment of the pixel with low spatial resolution scale, the first increment of the low spatial resolution scale and the second increment of the low spatial resolution scale can be obtained. The Mean Square Error (MSE) between the first increment of low spatial resolution and the true NDVI increment of low spatial resolution is calculated from the first increment of low spatial resolutionThe mean square error is the probability density function g_dThe same way to obtain the probability density function g_sThe variance of (c). g_s(△NDVI_i|△NDVI_s,i) Representing a probability density function g_sUsing its mean squared error as variance, and △ NDVI_s,iAs an average, △ NDVI can be generated_iThe probability of (c). Likewise, g_d(△NDVI_i|△NDVI_d,i) Representing a probability density function g_dUsing its mean squared error as variance, and △ NDVI_d,iAs an average, △ NDVI can be generated_iThe probability of (c). By maximizing the likelihood function (3), the weight w can be obtained_sAnd w_d。

The values are illustrated although the NDVI integrated delta △ NDVI_com(x_i，y_j) The sub-pixel from t which is close to high spatial resolution and can be better₁Time to t_pThe actual NDVI increment of time, but it may still contain errors, so in order to obtain more accurate results, in a preferred embodiment, further error correction may be performed on the obtained NDVI integrated increment, and the correction steps are:

in equation (4), m represents the number of high spatial resolution pixels covered by the low spatial resolution pixel (x, y), which is shown as 64 in FIG. 1, 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).

Here, the error R (x, y) contained in the low spatial resolution pel scale needs to be assigned to the high spatial resolution pel (x, y) it covers_i，y_j) The above. Generally, homogeneity is imagewiseRegions of relatively high quality (e.g., a large body of water or a large piece of grass) tend to have relatively small prediction errors, while regions of low homogeneity (e.g., sparse forest areas with trees and bare earth staggered) tend to have relatively large prediction errors. The invention is thus based on the homogeneity (x) of the high spatial resolution picture elements_i，y_j) To assign the error R (x, y).

I.e. 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 window. m denotes the number of high spatially resolved picture elements in this window, here 64, n_sameIndicating that in this window, the pixel is centered on (x)_i，y_j) The number of high spatial resolution picture elements having the same type.

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.

in equation (7), △ NDVI_com ^new(x_i，y_j) I.e. the final calculationHigh spatial resolution picture element (x)_i，y_j) From t₁Time to t_pA corrected NDVI integrated increment of time.

Step C, calculating t_pA high spatial resolution NDVI of time;

t already calculated in the previous step₁Time to t_pTemporal, each high spatial resolution pel (x)_i，y_j) NDVI synthetic increment of (1), hereinafter the synthetic increment is added to t₁T can be obtained on high spatial resolution NDVI of time_pHigh spatial resolution NDVI data in time. The calculation is as follows,

in formula (8), NDVI_p(x_i，y_j) Is t_pHigh spatial resolution NDVI in time₁(x_i，y_j) Is t₁High spatial resolution NDVI in time.

As described above, after the processing of steps A, B and C, the time t is predicted_pAll the high spatial resolution pixels can obtain a relatively accurate fusion result through a prediction process of space-time weighting. Compared with the existing data fusion method, the method has the advantages that the information contained in the low-spatial-resolution data is more fully used through a BMA (Bayesian model average) method, and the problems caused by the change of the land cover type can be better solved.

To better illustrate the technical effect of the present invention, the method of the present invention is compared with the three types of data fusion methods mentioned in the background section, namely, the time series data test and the single scene data test are respectively performed for STARFM (spatial and temporal adaptive reflection fusion model), NDVI-LMGM (NDVI-linear mixed growth model) and FSDAF (flexible spatio-temporal data fusion method). In the test, NDVI is taken as a fusion object, low spatial resolution data is MODIS data, and high spatial resolution data is Landsat data. Since the transit times of the two data acquisitions are the same and the orbit parameters are close, it is better to use them for the fusion process.

FIGS. 2a-2f show a land cover change area and corresponding Landsat NDVI data and MODISNDVI data, respectively; in the case of a certain Catchment basin area shown therein, the specifically selected test data is located in the Lower Gwydir catch site in northern New south Wales, Australia, with a latitude and longitude of 29 DEG 07 'S and 149 DEG 04' E, respectively. The area is a typical crop planting area, a relatively large flood accident occurs in 2004 for 12 months, a large amount of farmlands are submerged, and relatively obvious land coverage changes occur. To verify the detection of land cover changes by the present invention, the results were obtained from the official website of the australian federal Scientific and industrial Research Organization (CSIRO): https:// data. csiro. au/dap/landingpage? The execution is e2s2& _ evenld ═ viewDescription, and Landsat TM images on 26 days 11 and 26 days 2004 and 12 months 12 and 12 after the flood were downloaded and cropped to calculate the corresponding NDVI (fig. 2a and 2 b). Unsupervised classification by isodata was performed on images cut on day 26/11/2004 (fig. 2c) to obtain classification results (fig. 2 f). In order to avoid systematic errors among different sensors and errors caused in the processes of geometric correction, atmospheric correction and the like, the invention uses Landsat data to perform 8 × 8 upscaling synthesis to obtain MODIS data (fig. 2d and 2 e).

Fig. 3a to 3e respectively show the fusion results of several different data fusion methods in the land cover change area, and the Landsat NDVI data of 26/11/2004 is taken as reference data to fuse the data of 12/2004. The fusion results of the four methods (fig. 3a, 3b, 3c and 3d) were compared to the real NDVI at 12 months and 12 days 2004 (fig. 3 e). From the whole of the partial enlarged view, the fusion result of the NDVI-LMGM method has obvious plaque problems, while the fusion result of the method of the invention has better smoothness. As can be seen from the smaller boxes in fig. 3a-3 e; the detection capability of the STARFM fusion on the land cover change is not enough, the black area in the smaller square box is not identified, and the method of the invention can well identify the land cover change; from a larger square frame, the FSDAF fusion result is relatively greatly influenced by classification, the fusion result obviously reserves the water body boundary before flood occurs, so that the transition of the water covered area is not smooth, the method disclosed by the invention is relatively slightly influenced by classification, the recognition capability of the land coverage is stronger, and the fusion result is closer to the true value.

In conclusion, the increment of the mixed pixel decomposition and the increment of the spatial interpolation are integrated based on the Bayesian model averaging method, and the incremental calculation process is optimized. The spatio-temporal data are integrated, and the precision of data fusion is effectively improved. The invention fully utilizes the space-time correlation among the geographic data and provides an effective method for producing the vegetation index with high space-time resolution.

It should be appreciated by those of skill in the art that while the present invention has been described in terms of several embodiments, not every embodiment includes only a single embodiment. The description is given for clearness of understanding only, and it is to be understood that all matters in the embodiments are to be interpreted as including technical equivalents which are related to the embodiments and which are combined with each other to illustrate the scope of the present invention.

The above description is only an exemplary embodiment of the present invention, and is not intended to limit the scope of the present invention. Any equivalent alterations, modifications and combinations can be made by those skilled in the art without departing from the spirit and principles of the invention.

Claims

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；

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)

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:

<mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mi>&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>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <msub> <mi>&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>

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

<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>&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>&rsqb;</mo> <mo>/</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>8</mn> </munderover> <mo>&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>&rsqb;</mo> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msubsup> <mi>&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>&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>

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>&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>