CN107103584B - A method of the production high-spatial and temporal resolution NDVI based on space-time weighting - Google Patents
A method of the production high-spatial and temporal resolution NDVI based on space-time weighting Download PDFInfo
- Publication number
- CN107103584B CN107103584B CN201710231237.4A CN201710231237A CN107103584B CN 107103584 B CN107103584 B CN 107103584B CN 201710231237 A CN201710231237 A CN 201710231237A CN 107103584 B CN107103584 B CN 107103584B
- Authority
- CN
- China
- Prior art keywords
- ndvi
- pixel
- increment
- spatial resolution
- moment
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformation in the plane of the image
- G06T3/40—Scaling the whole image or part thereof
- G06T3/4053—Super resolution, i.e. output image resolution higher than sensor resolution
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration by the use of more than one image, e.g. averaging, subtraction
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10028—Range image; Depth image; 3D point clouds
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20212—Image combination
- G06T2207/20221—Image fusion; Image merging
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30181—Earth observation
Abstract
The present invention provides the methods of production high-spatial and temporal resolution NDVI based on space-time weighting a kind of, comprising: to t1The NDVI data of the low spatial resolution at moment and the NDVI data and t of high spatial resolutionpThe low spatial resolution NDVI data at moment obtain the two kinds of NDVI increments of high spatial resolution data in time by Decomposition of Mixed Pixels and space interpolation；Weight combination is carried out to described two NDVI increments using Bayesian model averaging method to calculate, and obtains the NDVI synthesis increment of each pixel；NDVI synthesis increment is finally added to t1It can be obtained required high-spatial and temporal resolution NDVI on the NDVI of the high spatial resolution at moment.The present invention carries out Bayes's synthesis to two kinds of NDVI increments, takes full advantage of the space-time relationship between geodata, provides a kind of effective method to produce the vegetation index of high-spatial and temporal resolution.
Description
Technical field
The present invention relates to a kind of method of multi-sources RS data fusion, it is mainly used to produce the normalization of high-spatial and temporal resolution
Vegetation index (NDVI) belongs to the scope of space-time data fusion.
Background technique
Normalized differential vegetation index (Normalized Difference Vegetation Index, be abbreviated as NDVI) is one
The widely used vegetation index of kind.NDVI numerical value is able to reflect the upgrowth situation of vegetation and green with the physiological parameter of vegetation such as leaf
Cellulose content, water content and leaf area index etc. have direct connection.
CN 102831310B discloses a kind of method of time series data for constructing high spatial resolution, background skill
The relevant knowledge of NDVI data is described in detail in art part, wherein referring to two kinds of existing NDVI data: Landsat
NDVI data and MODIS NDVI data.Wherein, the spatial resolution of the NDVI numerical value in Landsat NDVI data is 30m*
30m, temporal resolution is 16 days, and the maximum space resolution ratio of the NDVI numerical value in MODIS NDVI data is 250m*250m,
Temporal resolution is 1 day.That is, Landsat NDVI data spatial resolution with higher and lower time resolution
Rate, and MODIS NDVI data have lower spatial resolution and higher temporal resolution.
In order to obtain the NDVI data of high spatial resolution and high time resolution simultaneously, the prior art is developed
A variety of pairs of Landsat NDVI data and MODIS NDVI data merge, to obtain the NDVI data of high spatial resolution
Method mainly includes three types: (1) based on the method for weight function；(2) based on the method for Decomposition of Mixed Pixels；(3) based on mixed
The method of molding type.
Wherein, (1) with STARFM, (room and time adaptively reflects Fusion Model, Spatial based on the method for weight function
And Temporal Adaptive Reflectance Fusion Model) it is Typical Representative, referring to Gao, F., Masek,
J., Schwaller, M. , &Hall, F. (2006), " merge to Landsat with the Reflectivity for Growing Season of MODIS: prediction is daily
Landsat Reflectivity for Growing Season ", the geoscience of institute of Electrical and Electronic Engineers and remote sensing journal, 44,2207-2218.
(Gao, F., Masek, J., Schwaller, M. , &Hall, F. (2006) .On the blending of the Landsat
and MODIS surface reflectance:Predicting daily Landsat surface
Reflectance.IEEE Transactions on Geoscience and Remote Sensing, 44,2207-2218).
Such method is conceived to the search of similar pixel and its determination of weight, by current pixel and surrounding pixel in time, space
And the similitude in the features such as spectrum determines similar pixel and calculates corresponding weight, using the average increment of similar pixel as
The increment of current pixel.Such methods often just have relatively good effect in the relatively high region of degree of purity.
(2) based on the method for Decomposition of Mixed Pixels with NDVI-LMGM (NDVI- linear hybrid growth model, NDVI
Linear Mixing Growth Model) it is representative, referring to Rao, Y., Zhu, X., Chen, J. , &Wang, J. (2015),
" a kind of ground that upgraded uses multidate MODIS NDVI data and Landsat TM/ETM+ image production high-spatial and temporal resolution
The method of 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.Remote
Sensing, 7,7865-7891).Such method calculates the increment of each pixel in time by Decomposition of Mixed Pixels, from
And the NDVI data of predicted time are obtained, the particular content in relation to this method also can be found in aforementioned prior art CN
102831310B, the two are substantially the same.The such methods region relatively high for heterogeneous degree, which can also have, compares
Good effect, but it assumes that land cover pattern situation does not change, and this point not necessarily can satisfy.
(3) based on the method for mixed model with FSDAF (flexible space-time data fusion method, Flexible
Spatiotemporal Data Fusion method) it is representative, referring to Zhu, X., Helmer, E.H., Gao, F., Liu,
D., Chen, J. , &Lefsky, M.A. (2016), it is " a kind of that the satellite image of different resolution is neatly subjected to temporal-spatial fusion
Method ", environmental remote sensing, 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).Such method will be more
Kind model, such as weight function model, Pixel Unmixing Models and spatial interpolation model combine, and high-resolution is calculated
The increment of each pixel on rate image.Such methods are often more complicated, but precision is relatively high, and can cover to soil
Lid variation carries out a degree of capture.
The above three classes method all achieves relatively good effect in the field of each self application.But they are not all accounted for
To two phase low spatial resolution data interpolatings are used, the increment between two phase interpolation results can to the mode of high spatial resolution
Reflect the increment of grid cell size to a certain extent, while also containing the information of land cover pattern variation.The present invention is based on pattra leaves
This model averaging method integrates the increment of the increment of Decomposition of Mixed Pixels and space interpolation, can effectively improve data and melt
The precision of conjunction.New approach is provided to produce the NDVI of high-spatial and temporal resolution.
Summary of the invention
The technical problem to be solved in the present invention is to provide a kind of production high-spatial and temporal resolution NDVI's based on space-time weighting
Method, the problem of being formerly mentioned is reduced or avoided.
In order to solve the above technical problems, the invention proposes a kind of production high-spatial and temporal resolution NDVI based on space-time weighting
Method, for by the Landsat NDVI data of high spatial resolution and low temporal resolution and low spatial resolution and Gao Shi
Between resolution ratio MODIS NDVI data fusion, from reference instant t1Obtain prediction time tpHigh-spatial and temporal resolution NDVI number
According to described method includes following steps:
Step A: with same reference instant t1Divided for starting point using MODIS NDVI data and Landsat NDVI data
It Huo get not prediction time tpFirst increment △ NDVI of each pixel on NDVIdWith the second increment △ NDVIs；Wherein,
The first increment △ NDVIdIt refers to from reference instant t1To prediction time tp, it is assumed that MODIS NDVI data
Each pixel NDVI increment, equal to the other pixel of differently species of Landsat NDVI data under same area
The average value of NDVI increment then utilizes t1The MODIS NDVI data and Landsat NDVI data and t at momentpMoment
MODIS NDVI data calculate and obtain t1Moment is to tpThe NDVI increment of each pixel of the Landsat NDVI data at moment, should
Increment is the first increment △ NDVId；
The second increment △ NDVIsIt refers to using the thin plate spline method work integrated in existing ENVI/IDL software
Tool, directly by t1Moment and tpThe low spatial resolution MODIS pixel at moment is interpolated into the Landsat of high spatial resolution respectively
Pixel gets on, and then uses tpThe interpolation result and t at moment1The interpolation result at moment subtracts each other, and result is the institute of each pixel
State the second increment △ NDVIs；
Step B is existed the pixel of the step A each high spatial resolution obtained using Bayesian model averaging method
First increment on NDVI carries out weight with second increment and combines calculating, obtains the comprehensive increasing of NDVI of each pixel
Amount；
Step C, the t being calculated by step B1Moment is to tpMoment, the NDVI synthesis increment of each pixel add
To t1Numerically, t can be obtained in moment corresponding high spatial resolution NDVIpThe high spatial resolution NDVI numerical value at moment.
Preferably, in the step B, each high spatial resolution pixel (xi, yj) the NDVI synthesis increment formula
It is expressed as,
ΔNDVIcom(xi,yj)=ws*ΔNDVIs(xi,yj)+wd*ΔNDVId(xi,yj) (1)
Wherein, in formula (1), wdAnd wsRespectively indicate the weight of first increment and second increment, △ NDVIcom
(xi, yj) indicate the NDVI synthesis increment that first increment obtains in conjunction with second increment.
Preferably, the weight w of each pixel in the formula (1)dAnd wsIt is asked by Bayesian model averaging method
Solution obtains.
Preferably, the step B further comprises the step of carrying out error correction to the NDVI synthesis increment.
Preferably, the step of error correction includes:
The error for calculating low spatial resolution pixel (x, y) entirety first is,
In formula (2), m indicates the number for the high spatial resolution pixel that low spatial resolution pixel (x, y) is covered；R
The error of (x, y) expression low spatial resolution pixel (x, y) entirety；△NDVIC(x, y) indicates tpMoment and t1Moment low spatial
The difference of the NDVI of resolution ratio pixel (x, y)；
Then the homogenieity for calculating high spatial resolution pixel is as follows,
HI(xi,yj)=nsmae/ z, (3)
In formula (3), HI (xi, yj) it is with high spatial resolution pixel (xi, yj) centered on, the window of the 8*8 of building
Homogenieity parameter；Z indicates the number of resolution pixel between this window high and medium, nsameIt indicates in this window, with center pel
(xi, yj) with same type high spatial resolution pixel number；
It is later, as follows to be allocated to error according to the homogenieity of pixel,
In formula (4), r (xi, yj) it is high spatial resolution pixel (xi, yj) error that should be assigned to；
Finally, the synthesis increment that this error is added in formula (1) is got on, corrected high spatial resolution picture is obtained
The NDVI synthesis increment of member is as follows,
In formula (5), △ NDVIcom new(xi, yj) it is the high spatial resolution pixel (x being finally calculatedi, yj) from
t1Moment is to tpThe corrected NDVI synthesis increment at moment.
Preferably, in the step C, t is calculated using the corrected NDVI synthesis incrementpThe high-space resolution at moment
The formula of rate NDVI are as follows:
In formula (6), NDVIp(xi, yj) it is tpEach pixel (x at momenti, yj) high spatial resolution NDVI, NDVI1
(xi, yj) it is t1Each pixel (x at momenti, yj) high spatial resolution NDVI.
The present invention by data to low spatial resolution and high spatial resolution data carry out Decomposition of Mixed Pixels and
Space interpolation obtains the two kinds of increments of high spatial resolution data in time, then uses Bayesian model averaging method by institute
It states two kinds of increments and carries out weight combination calculating, the synthesis increment optimized improves the accuracy of incremental computations, sufficiently benefit
With the space-time connection between geodata, fusion accuracy is improved.
Detailed description of the invention
The following drawings are only intended to schematically illustrate and explain the present invention, not delimit the scope of the invention.Wherein,
The schematic illustration of the first increment of a specific embodiment according to the present invention is shown in Fig. 1；
Land cover pattern region of variation and corresponding Landsat NDVI data and MODIS is shown in Fig. 2 a-2f respectively
NDVI data；Certain wherein indicated catchment basin area the case where, be respectively as follows: Fig. 2 a indicate on November 26th, 2004 Landsat
The Landsat NDVI in 12 days December in 2004 that NDVI, Fig. 2 b are indicated；It is simulated respectively from above-mentioned image and obtains Fig. 2 d expression
The MODIS NDVI in 12 days December in 2004 that on November 26th, 2004 and Fig. 2 e are indicated；What Fig. 2 c was indicated is on November 26th, 2004
Landsat image, Fig. 2 f indicates corresponding unsupervised classification result.
Several different data fusion methods are shown in the fusion knot of land cover pattern region of variation in Fig. 3 a-3e respectively
Fruit；What wherein Fig. 3 a-3d was respectively indicated is Landsat NDVI on December 12 in 2004 by NDVI-LMGM fusion, STARFM
Fusion, FSDAF fusion, method of the invention fusion as a result, what Fig. 3 e indicated is true value.
Specific embodiment
For a clearer understanding of the technical characteristics, objects and effects of the present invention, this hair of Detailed description of the invention is now compareed
Bright specific embodiment.But it will be appreciated by those skilled in the art that following embodiment is not to technical solution of the present invention
The unique restriction made, all any equivalents or change done under technical solution of the present invention Spirit Essence are regarded as
It belongs to the scope of protection of the present invention.
The following detailed description of the method for the production high-spatial and temporal resolution NDVI based on space-time weighting provided according to the present invention
Principle Landsat NDVI data and MODIS NDVI data are carried out with three types described in background parts certainly
As the method for fusion, the present invention is also the method merged to the remotely-sensed data of above-mentioned different data sources, also that is, this hair
The method of the bright production high-spatial and temporal resolution NDVI based on space-time weighting is used for high spatial resolution and low temporal resolution
Landsat NDVI data and low spatial resolution and high time resolution MODIS NDVI data fusion, from reference instant
t1Obtain prediction time tpHigh-spatial and temporal resolution NDVI data.
Satellite due to obtaining Landsat NDVI data and MODIS NDVI data has similar orbit parameter and not
30 minutes satellite transit time intervals (interval can be ignored for daily 24 hours) of foot, therefore, in theory
On assume that the same place of both data and same reference instant t1NDVI data there is only the differences in spatial resolution
Not, remaining is negligible.Wherein, the spatial resolution of Landsat NDVI data is 30m*30m, temporal resolution 16
It；The spatial resolution of MODIS NDVI data is 250m*250m, and temporal resolution is 1 day；It then produces and obtains after data fusion
The NDVI data of the high-spatial and temporal resolution obtained, maximum space resolution ratio can be equal to the spatial discrimination of Landsat NDVI data
Rate (30m*30m), maximum time resolution ratio can be equal to the temporal resolution (1 day) of MODIS NDVI data.That is,
The MODIS NDVI data and spatial resolution that by known spatial resolution ratio be 250m*250m, temporal resolution is 1 day
The Landsat NDVI data for being 16 days for 30m*30m, temporal resolution, with same reference instant t1It is of the invention for starting point
Method can be predicted to obtain 15 days thereafter any 1 day prediction time tpHigh-spatial and temporal resolution NDVI data.
The following detailed description of the specific step of the method for the production high-spatial and temporal resolution NDVI of the invention based on space-time weighting
Suddenly.
Step A: with same reference instant t1Divided for starting point using MODIS NDVI data and Landsat NDVI data
It Huo get not prediction time tpFirst increment △ NDVI of each pixel on NDVIdWith the second increment △ NDVIs。
Wherein, the first increment △ NDVIdIt refers to from reference instant t1To prediction time tp, it is assumed that MODIS NDVI
The NDVI increment of each pixel of data, equal to the other pixel of differently species of the Landsat NDVI data under same area
NDVI increment average value, then utilize t1The MODIS NDVI data and Landsat NDVI data and t at momentpMoment
MODIS NDVI data calculate and obtain t1Moment is to tpThe NDVI increment of each pixel of the Landsat NDVI data at moment, should
Increment is the first increment △ NDVId。
Specifically, reference instant t1The low spatial resolution MODIS of (time interval that there are in fact 30 minutes)
There is only the differences in spatial resolution for same place for NDVI data and the Landsat NDVI data of high spatial resolution
Not, the position alignment of two kinds of data is facilitated into calculating when in order to merge, it can be by the MODIS NDVI data of low spatial resolution
From 250m*250m resampling to 240*240m, so that 8*8=64 Landsat picture can be completely covered in a MODIS pixel
Member.
As shown in Figure 1, the schematic illustration of its first increment that a specific embodiment according to the present invention is shown,
Wherein what the top of time shaft indicated is the MODIS NDVI data of low spatial resolution, and what is indicated below time shaft is high spatial
The Landsat NDVI data of resolution ratio.With same reference instant t1For starting point, t is utilized1The MODIS NDVI data at moment and
Landsat NDVI data and tpThe MODIS NDVI data at moment can obtain the original of time shaft lower right dotted line expression
Originally the t being not presentpThe Landsat NDVI data of the high spatial resolution at moment.
At this time, it is assumed that from reference instant t1To prediction time tp, the NDVI increment of each pixel of MODIS NDVI data,
Equal to the average value of the NDVI increment of the other pixel of differently species of the Landsat NDVI data under same area.That is, Fig. 1
The heart-shaped 1 MODIS pixel indicated is from t above middle time shaft1To tpThe NDVI increment at moment, can be by by tpMoment it is low
The MODIS NDVI data and t of spatial resolution1The MODIS NDVI data of the low spatial resolution at moment directly subtract each other to obtain.
Corresponding to the area of 1 MODIS pixel above time shaft, equal to 64 Landsat pixels below time shaft
Area, on an average according to its area by the NDVI increment of this 64 Landsat pixels, can regard as be equal to time shaft on
The NDVI increment of 1 MODIS pixel of the same area of side.Its calculating process is, according to the Landsat of high spatial resolution
The terrain classification of NDVI data calculates every kind of type of ground objects in the pixel of the MODIS NDVI data of each low spatial resolution
Abundance.Classification data can derive from existing sort product, such as GlobeLand30, can also be by high spatial point
The satellite image data of the Landsat NDVI data of resolution are classified to obtain.According to the Landsat of high spatial resolution
The terrain classifications of NDVI data is as a result, calculating in the pixels of the MODIS NDVI data of each low spatial resolution includes not
With the area ratio of type of ground objects, i.e. abundance.Therefore, by the product addition of the NDVI increment of every kind of type of ground objects and its abundance,
It can be obtained by the average value of the NDVI increment of whole Landsat pixels (64 pixels) under the area, this average value is just
Well equal to the MODIS pixel (1 pixel) of aforementioned same area from t1To tpThe NDVI increment at moment.
For example, 64 Landsat pixels below Fig. 1 time shaft include three kinds of types of ground objects, then the 1st class type of ground objects
24 pixels have a NDVI increment, and 20 pixels of the 2nd class type of ground objects have another NDVI increment, the 3rd class atural object
20 pixels of type have another NDVI increment, these three NDVI increments, are it for each pixel from t1
To tpThe first increment △ NDVI at momentd.In order to calculate these three increments, needs to construct at least three equations, assume at this time most
The NDVI of 8 pixels (that is, around 8 pixels of heart-shaped label pixel) of heart-shaped label pixel above neighbouring time shaft increases
Amount is identical as the heart label NDVI increment of pixel, then can construct 3*3's centered on each low spatial resolution pixel
Local window.Assuming that each type of ground objects is in t in this local window1Moment is to tpNDVI increment in moment is identical.
According to the NDVI increment of 9 low spatial resolution pixels in this local window, following system of linear equations, formula can be constructed
It is expressed as follows,
In formula (1-1), △ NDVIC(x, y) is indicated from t1To tpThe NDVI of moment low spatial resolution pixel (x, y) increases
Amount, fl(x, y) indicates abundance of the l seed type in low spatial resolution pixel (x, y), △ NDVIlThe office of (x, y) expression 3*3
The NDVI increment of l seed type in portion's window.
It is carried out solving available △ NDVI by least square according to formula (1-1)l(x, y).For illustrating local window
For low spatial resolution pixel (x, y) in mouthful, the high-resolution pixel (x under coverage areai, yj) increasing in time
Amount is △ NDVId(xi, yj), this increment is known as the first increment, and subscript (i, j) indicates under the covering of low spatial resolution pixel
(i, j) a high spatial resolution pixel.It is hypothesized in this local window, the NDVI increment of every kind of type of ground objects is solid
Fixed.Therefore, for high spatial resolution pixel (xi, yj) for, increment △ NDVI in timed(xj, yj) be equal to
High spatial resolution pixel (xj, yj) belonging to that type of ground objects NDVI increment.If high spatial resolution pixel (xi,
yj) belong to l class, then △ NDVId(xi, yj) it is equal to △ NDVIl(x, y).
For the first increment △ NDVIdFor, due to tpThe type of ground objects and t at moment1The type of ground objects at moment may have compared with
It is different, such as soil erosion, deforestation, fire, flood etc. can cause the type of ground objects of different time to change.Cause
This, the abundance f of every kind of type of ground objects in formula (1-1)l(x, y) is possible to be different under two time, will lead to ground
The NDVI increment △ NDVI of species typel(x, y) calculates error.And by the NDVI increment △ NDVI of every kind of type of ground objectsl(x,
Y) t is distributed to by typepIt can also malfunction when the high spatial resolution pixel at moment.
In order to avoid the influence of this unfavorable factor, therefore, present invention further proposes use the second increment △ NDVIs
Come the design being modified to the first increment.
Wherein, the second increment △ NDVIsIt refers to using the thin plate spline integrated in existing ENVI/IDL software
(Thin Plate Spline, TPS) methods and instruments, directly by t1Moment and tpThe low spatial resolution MODIS pixel at moment point
Be not interpolated into the Landsat pixel of high spatial resolution up (the interpolation method tool has been integrated in ENVI/IDL software,
For well known method, function is called to can be completed).Then, using tpThe interpolation result and t at moment1The interpolation result phase at moment
Subtract, result is the second increment △ NDVI of each pixels.Due to the second increment △ NDVIsHave ignored type of ground objects
Difference, can not only provide information spatially, and not by the interference of type of ground objects, therefore by following step by
One and second increment merged, the first increment △ NDVI can be eliminateddDefect and adverse effect.
That is, 1 MODIS pixel of heart-shaped label shown in FIG. 1,64 Landsat pixels under corresponding to, each pixel
First increment △ NDVIdIt is related to type of ground objects, it is indicated in figure there are three types of type of ground objects, then 24 under the 1st class type of ground objects
First increment △ NDVI of pixeldIt is all identical, the first increment △ NDVI of 20 pixels under the 2nd class type of ground objectsdIt is also
It is identical, the first increment △ NDVI of 20 pixels under the 3rd class type of ground objectsdAlso be all it is identical, certainly, these three first
Increment △ NDVIdNumerical value calculate separately acquisition, be independent from each other.Above-mentioned first increment △ NDVIdIt is it is assumed that same
What the atural object classification of pixel obtained in the case that no change has taken place, this will obviously have certain problems.
And the second increment △ NDVIsAcquisition be the methods and instruments according to thin plate spline directly by t1Moment and tpMoment
The MODIS pixel interpolation of low spatial resolution obtains tpLandsat pixel (the side of thin plate spline of the high spatial resolution at moment
Method principle is to find out multiple match points in two low spatial resolution MODIS images, can be by this using thin plate spline function
Multiple points are deflected into corresponding position, while giving the interpolation of whole image), the t that interpolation is obtainedpThe high-space resolution at moment
The NDVI numerical value of each pixel (such as corresponding to 64 Landsat pixels in Fig. 1) of the Landsat image of rate with it is known
t1The NDVI numerical value of each pixel of the Landsat image at moment subtracts each other, and can be obtained each pixel from t1To tpThe second of moment
Increment △ NDVIs, i.e. the second increment △ NDVIsUnrelated with type of ground objects, each pixel all has an independent second increment △
NDVIs.Certainly, since the Landsat of high spatial resolution can be obtained from the direct interpolation of MODIS pixel of low spatial resolution
Pixel then individually calculates and obtains the second increment △ NDVIsBe it is nugatory, therefore, propose the second increment △ NDVI hereins's
Purpose is exactly the fusion for next step, individually uses the first increment △ NDVI to eliminatedWhen obtaining high-spatial and temporal resolution NDVI
Brought problem.
Step B, using Bayesian model averaging method, by first increasing of the pixel of each high spatial resolution on NDVI
Amount carries out weight with the second increment and combines calculating, obtains the NDVI synthesis increment of each pixel.
Come used here as two weights to each pixel (xi, yj) corresponding to the first increment △ NDVId(xi, yj) and second
Increment △ NDVIs(xi, yj) carry out weight combination.Obtain each high spatial resolution pixel (xi, yj) NDVI synthesis increment be,
ΔNDVIcom(xi,yj)=ws*ΔNDVIs(xi,yj)+wd*ΔNDVId(xi,yj) (1)
In formula (1), wdAnd wsThe weight of the first increment and the second increment is respectively indicated, the two weights are not with pixel
Variation, only one w of entire imagedWith a ws。△NDVIcom(xi, yj) indicate the first increment in conjunction with the second increment
The NDVI synthesis increment arrived.
(Bayesian Model Averaging, BMA) method will be needed to be averaged using Bayesian model below to solve
The two weights of each pixel.BMA method is a kind of classical multi-model integrated method of comparison, it by likelihood function into
Row optimizes to calculate the weight of final each model.Detailed content may refer to bibliography Raftery, A.E., Gneiting,
T., Balabdaoui, F. , &Polakowski, M. (2005), " using Bayesian model averaging method to the DATA PROCESSING IN ENSEMBLE PREDICTION SYSTEM of weather
Correction ", monthly weather is commented on, 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).First increment in the present invention and second are increased
The maximum likelihood function of amount is expressed as follows,
In formula (1-2), L indicates maximum likelihood value, wsAnd wdRespectively indicate the weight of two kinds of increments, △ NDVId,iIt indicates
First increment of i-th of pixel, △ NDVIs,iIndicate second increment of i-th of pixel, △ NDVIiIndicate the true of i-th of pixel
Real increment.gsAnd gdRespectively indicate the probability density function of two increment sets, it is believed that it is normal distribution, and N indicates training sample
Number.In the present invention, due to the true NDVI incremental data of not no high spatial resolution, t is used1Moment and tpMoment
The true increment of low resolution pixel all carries out 8* as training sample, by the first increment of high spatial resolution and the second increment
8 synthesis rises on scale to low spatial resolution.Can be obtained in this way low spatial resolution scale the true NDVI increment of pixel,
First increment of low spatial resolution and the second increment of low spatial resolution.Gone out according to the first incremental computations that low spatial is differentiated
Its Mean Square Error (Mean Square Error, MSE) between the true NDVI increment of low spatial resolution, this is flat
Mean square error is exactly probability density function gdVariance, can similarly obtain probability density function gsVariance.gs(△NDVIi
|△NDVIs,i) indicate probability density function gsUsing its Mean Square Error as variance, with △ NDVIs,iFor mean value, △ can be generated
NDVIiProbability.Equally, gd(△NDVIi|△NDVId,i) indicate probability density function gdUsing its Mean Square Error as variance,
With △ NDVId,iFor mean value, △ NDVI can be generatediProbability.By maximizing to likelihood function (3), available weight
wsAnd wd。
The explanation of value, although NDVI synthesis increment △ NDVIcom(xi, yj) have been able to it is relatively good close to high spatial
Resolution ratio pixel is from t1Moment is to tpThe true NDVI increment at moment, but it still may include error, therefore in order to obtain more
Accurately as a result, in a preferred embodiment, further error correction can also be carried out to the NDVI synthesis increment of acquisition,
Its amendment step are as follows:
The error for calculating low spatial resolution pixel (x, y) entirety first is,
In formula (2), m indicates the number for the high spatial resolution pixel that low spatial resolution pixel (x, y) is covered,
What is shown in Fig. 1 is 64；The error of R (x, y) expression low spatial resolution pixel (x, y) entirety；△NDVIC(x, y) indicates tp
Moment and t1The difference of the NDVI of moment low spatial resolution pixel (x, y).
Herein, it needs the error R (x, y) that low spatial resolution grid cell size is included being assigned to its high-altitude covered
Between resolution ratio pixel (xi, yj) on.In general, (such as a sheet of water body or one big of the relatively high region of homogenieity on image
Piece meadow) often pre- time determination error is smaller, and region that homogenieity is low (such as trees and exposed soil are interspersed sparse wood)
Often prediction error is larger.Therefore, the present invention is according to the homogenieity (x of high spatial resolution pixeli, yj) come to error R (x, y)
It is allocated.
That is, the homogenieity for then calculating high spatial resolution pixel is as follows,
HI(xi,yj)=nsmae/ z, (3)
In formula (3), HI (xi, yj) it is with high spatial resolution pixel (xi, yj) centered on, the window of the 8*8 of building
Homogenieity parameter.Z indicates the number that pixel is differentiated between this window high and medium, i.e. 64 here, nsameIt indicates in this window,
With center pel (xi, yj) with same type high spatial resolution pixel number.
It is later, as follows to be allocated to error according to the homogenieity of pixel,
In formula (4), r (xi, yj) it is high spatial resolution pixel (xi, yj) error that should be assigned to.
Finally, the synthesis increment that this error is added in formula (1) is got on, corrected high spatial resolution picture is obtained
The NDVI synthesis increment of member is as follows,
In formula (5), △ NDVIcom new(xi, yj) it is the high spatial resolution pixel (x being finally calculatedi, yj) from
t1Moment is to tpThe corrected NDVI synthesis increment at moment.
Step C calculates tpThe high spatial resolution NDVI at moment；
The t calculated in the step of front1Moment is to tpMoment, each high spatial resolution pixel (xi, yj)
NDVI synthesis increment, is added to t for synthesis increment below1On the high spatial resolution NDVI at moment, t can be obtainedpThe high-altitude at moment
Between resolution, N DVI data.Calculating is as follows,
In formula (6), NDVIp(xi, yj) it is tpThe high spatial resolution NDVI, NDVI at moment1(xi, yj) it is t1Moment
High spatial resolution NDVI.
As described above, after the processing of step A, B and C, prediction time tpAll high spatial resolution pixels can
More accurate fusion results are obtained by the prediction process that space-time weights.For existing data fusion method, this
Invention has more fully used information included in low spatial resolution data by BMA (Bayesian model is average) method,
The problem of Land cover types variation causes can preferably be coped with.
Technical effect in order to better illustrate the present invention, three that the method in the present invention is referred to background technology part
The data fusion method of seed type is compared, that is, is directed to STARFM (room and time adaptively reflects Fusion Model), NDVI-
LMGM (NDVI- linear hybrid growth model) and FSDAF (flexible space-time data fusion method) carries out the survey of time series data respectively
The test of examination and single scape data.With NDVI to merge object in test, low spatial resolution data are MODIS data, high spatial
Resolution data is Landsat data.Since the transit time of both data acquisitions is identical, and orbit parameter is close, because
This is relatively more suitable for for carrying out fusion treatment.
Land cover pattern region of variation and corresponding Landsat NDVI data and MODIS is shown in Fig. 2 a-2f respectively
NDVI data；Certain wherein indicated catchment basin area the case where, the test data of specific choice is located at the new Nan Weier of Australia
The downstream the Gwydir basin research area (Lower Gwydir Catchment site) in scholar the north, longitude and latitude is respectively 29 ° of 07 ' S
With 149 ° of 04 ' E.The region is a typical crop planting area, and primary bigger flood has occurred in December, 2004
Water event has a large amount of farmlands to be submerged, and obvious land cover pattern variation has occurred.Soil is covered in order to verify the present invention
The detection effect for covering variation, from Australia Commonwealth Scientific and Industrial Research Organisation (Commonwealth Scientific and
Industrial Research Organization, CSIR O) official website: https: //data.csiro.au/dap/
Landingpage? execution=e2s2&_eventId=viewDescription downloads to flood and occurs first 2004
The Landsat TM image on December 12nd, 2004, cuts it, calculates correspondence after November 26 and flood occur
NDVI (Fig. 2 a and 2b).It is carried out according to the image (Fig. 2 c) that on November 26th, 2004 cut based on the non-supervisory of isodata
Classification, obtains classification results (Fig. 2 f).In order to avoid between different sensors systematic error and geometric correction and atmosphere school
Bring error during just waiting, the present invention synthesize to obtain MODIS data (figure using the liter scale that Landsat data carry out 8*8
2d and 2e).
Several different data fusion methods are shown in the fusion knot of land cover pattern region of variation in Fig. 3 a-3e respectively
Fruit melts the data on December 12nd, 2004 using the Landsat NDVI data on November 26th, 2004 as reference data
It closes.By the fusion results (Fig. 3 a, Fig. 3 b, Fig. 3 c and Fig. 3 d) of four kinds of methods and true NDVI (Fig. 3 e) on December 12nd, 2004
It is compared.It can see from the entirety of partial enlarged view, the fusion results of NDVI-LMGM method have apparent patch problem,
And the flatness of the result of method fusion of the invention is relatively good.Lesser box can see from Fig. 3 a-3e；STARFM melts
The detectability that conjunction changes land cover pattern is inadequate, and the black region in lesser box is not identified, and this hair
Bright method can well identify the variation of land cover pattern；It can see from biggish box, the knot of FSDAF fusion
Fruit by classify influenced it is relatively larger, in fusion results for flood occur before water boundary have apparent reservation, make
The transition that region must be covered with water is unsmooth, and method of the invention is influenced relatively small knowledge to land cover pattern by classifying
Other ability is stronger, and fusion results and true value are closer.
In conclusion the present invention is based on Bayesian model averaging methods by the increment of Decomposition of Mixed Pixels and space interpolation
Increment is integrated, and the calculating process of increment is optimized.Space-time data is integrated, data is effectively increased and melts
The precision of conjunction.The present invention takes full advantage of the space-time relationship between geodata, refers to produce the vegetation of high-spatial and temporal resolution
Number provides a kind of effective method.
It will be appreciated by those skilled in the art that although the present invention is described in the way of multiple embodiments,
It is that not each embodiment only contains an independent technical solution.So narration is used for the purpose of for the sake of understanding in specification,
The skilled in the art should refer to the specification as a whole is understood, and by technical solution involved in each embodiment
Regard as and can be combined with each other into the modes of different embodiments to understand protection scope of the present invention.
The foregoing is merely the schematical specific embodiment of the present invention, the range being not intended to limit the invention.It is any
Those skilled in the art, made equivalent variations, modification and combination under the premise of not departing from design and the principle of the present invention,
It should belong to the scope of protection of the invention.
Claims (6)
1. a kind of method of the production high-spatial and temporal resolution NDVI based on space-time weighting, is used for high spatial resolution and low time
The Landsat NDVI data of resolution ratio and the MODIS NDVI data fusion of low spatial resolution and high time resolution, from base
It is punctual to carve t1Obtain prediction time tpHigh-spatial and temporal resolution NDVI data, which is characterized in that the method includes walking as follows
It is rapid:
Step A: with same reference instant t1It is obtained respectively for starting point using MODIS NDVI data and Landsat NDVI data
Prediction time tpFirst increment △ NDVI of each pixel on NDVIdWith the second increment △ NDVIs；Wherein,
The first increment △ NDVIdIt refers to from reference instant t1To prediction time tp, it is assumed that MODIS NDVI data it is every
The NDVI increment of a pixel, the NDVI equal to the other pixel of differently species of the Landsat NDVI data under same area increase
The average value of amount then utilizes t1The MODIS NDVI data and Landsat NDVI data and t at momentpThe MODIS at moment
NDVI data calculate and obtain t1Moment is to tpThe NDVI increment of each pixel of the Landsat NDVI data at moment, the increment are
For the first increment △ NDVId；
The second increment △ NDVIsIt refers to using the thin plate spline methods and instruments integrated in existing ENVI/IDL software, directly
It connects t1Moment and tpThe low spatial resolution MODIS pixel at moment is interpolated into respectively on the Landsat pixel of high spatial resolution
It goes, then uses tpThe interpolation result and t at moment1The interpolation result at moment subtracts each other, and result is described the second of each pixel
Increment △ NDVIs；
Step B, using Bayesian model averaging method, by the pixel of the step A each high spatial resolution obtained on NDVI
First increment and second increment carry out weight and combine calculating, obtain the NDVI synthesis increment of each pixel；
Step C, the t being calculated by step B1Moment is to tpMoment, the NDVI synthesis increment of each pixel are added to t1
Numerically, t can be obtained in moment corresponding high spatial resolution NDVIpThe high spatial resolution NDVI numerical value at moment.
2. the method as described in claim 1, which is characterized in that in the step B, each high spatial resolution pixel (xi, yj)
The NDVI synthesis increment be formulated as,
ΔNDVIcom(xi,yj)=ws*ΔNDVIs(xi,yj)+wd*ΔNDVId(xi,yj) (1)
Wherein, in formula (1), wdAnd wsRespectively indicate the weight of first increment and second increment, △ NDVIcom(xi,
yj) indicate the NDVI synthesis increment that first increment obtains in conjunction with second increment.
3. method according to claim 2, which is characterized in that the weight w of each pixel in the formula (1)dAnd ws
It is solved and is obtained by Bayesian model averaging method.
4. method as claimed in claim 2 or claim 3, which is characterized in that the step B further comprises increases comprehensive to the NDVI
Amount carries out the step of error correction.
5. method as claimed in claim 4, which is characterized in that the step of error correction includes:
The error for calculating low spatial resolution pixel (x, y) entirety first is,
In formula (2), m indicates the number for the high spatial resolution pixel that low spatial resolution pixel (x, y) is covered；R (x, y)
Indicate the error of low spatial resolution pixel (x, y) entirety；△NDVIC(x, y) indicates tpMoment and t1Moment low spatial resolution
The difference of the NDVI of pixel (x, y)；
Then the homogenieity for calculating high spatial resolution pixel is as follows,
HI(xi,yj)=nsmae/ z, (3)
In formula (3), HI (xi, yj) it is with high spatial resolution pixel (xi, yj) centered on, the homogeneous of the window of the 8*8 of building
Property parameter；Z indicates the number of resolution pixel between this window high and medium, nsameIt indicates in this window, with center pel (xi,
yj) with same type high spatial resolution pixel number；
It is later, as follows to be allocated to error according to the homogenieity of pixel,
In formula (4), r (xi, yj) it is high spatial resolution pixel (xi, yj) error that should be assigned to；
Finally, the synthesis increment that this error is added in formula (1) is got on, corrected high spatial resolution pixel is obtained
NDVI synthesis increment is as follows,
In formula (5), △ NDVIcom new(xi, yj) it is the high spatial resolution pixel (x being finally calculatedi, yj) from t1When
It is carved into tpThe corrected NDVI synthesis increment at moment.
6. method as claimed in claim 5, which is characterized in that in the step C, increase using the corrected NDVI is comprehensive
Amount calculates tpThe formula of the high spatial resolution NDVI at moment are as follows:
In formula (6), NDVIp(xi, yj) it is tpEach pixel (x at momenti, yj) high spatial resolution NDVI, NDVI1(xi,
yj) it is t1Each pixel (x at momenti, yj) high spatial resolution NDVI.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231237.4A CN107103584B (en) | 2017-04-11 | 2017-04-11 | A method of the production high-spatial and temporal resolution NDVI based on space-time weighting |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231237.4A CN107103584B (en) | 2017-04-11 | 2017-04-11 | A method of the production high-spatial and temporal resolution NDVI based on space-time weighting |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107103584A CN107103584A (en) | 2017-08-29 |
| CN107103584B true CN107103584B (en) | 2019-06-14 |
Family
ID=59675152
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710231237.4A Expired - Fee Related CN107103584B (en) | 2017-04-11 | 2017-04-11 | A method of the production high-spatial and temporal resolution NDVI based on space-time weighting |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107103584B (en) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109033543B (en) * | 2018-06-29 | 2020-06-30 | 北京师范大学 | Method, device and equipment for estimating vegetation coverage of heterogeneous land areas |
| CN109685108B (en) * | 2018-11-23 | 2022-05-03 | 电子科技大学 | Method for generating high-space-time resolution NDVI long-time sequence |
| CN109612588A (en) * | 2019-01-02 | 2019-04-12 | 中国科学院新疆生态与地理研究所 | LST image data prediction technique, device and electronic equipment |
| CN110147832B (en) * | 2019-05-08 | 2021-03-05 | 北京师范大学 | Method model for obtaining high space-time sequence EVI (error vector index) based on Bayesian theory |
| CN110363246B (en) * | 2019-07-18 | 2023-05-09 | 滨州学院 | Fusion method of vegetation index NDVI with high space-time resolution |
| CN111125277B (en) * | 2019-11-14 | 2023-04-11 | 国家气候中心 | Landsat remote sensing vegetation index restoration method based on cube technology |
| CN110909821B (en) * | 2019-12-03 | 2020-07-28 | 中国农业科学院农业资源与农业区划研究所 | Method for carrying out high-space-time resolution vegetation index data fusion based on crop reference curve |
| CN111369483B (en) * | 2020-03-05 | 2020-11-13 | 北京师范大学 | Method for generating high-spatial-temporal-resolution remote sensing data by fusing multi-source remote sensing data |
| CN111523451B (en) * | 2020-04-22 | 2023-05-30 | 重庆邮电大学 | Method for constructing high space-time resolution NDVI data |
| CN112052720B (en) * | 2020-07-16 | 2021-04-20 | 贵州师范学院 | High-space-time normalization vegetation index NDVI fusion model based on histogram clustering |
| CN111898494B (en) * | 2020-07-16 | 2022-09-27 | 大同煤矿集团有限责任公司 | Mining disturbed land boundary identification method |
| CN114301905B (en) * | 2020-09-23 | 2023-04-04 | 华为技术有限公司 | Resolution conversion method and device |
| CN112819697B (en) * | 2021-02-04 | 2023-04-14 | 北京师范大学 | Remote sensing image space-time fusion method and system |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102831310A (en) * | 2012-08-17 | 2012-12-19 | 北京师范大学 | Method for building high-spatial resolution NDVI (normalized difference vegetation index) time series data |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9408342B2 (en) * | 2010-10-25 | 2016-08-09 | Trimble Navigation Limited | Crop treatment compatibility |
-
2017
- 2017-04-11 CN CN201710231237.4A patent/CN107103584B/en not_active Expired - Fee Related
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102831310A (en) * | 2012-08-17 | 2012-12-19 | 北京师范大学 | Method for building high-spatial resolution NDVI (normalized difference vegetation index) time series data |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107103584A (en) | 2017-08-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107103584B (en) | A method of the production high-spatial and temporal resolution NDVI based on space-time weighting | |
| Danner et al. | Efficient RTM-based training of machine learning regression algorithms to quantify biophysical & biochemical traits of agricultural crops | |
| Zhu et al. | An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions | |
| Verhegghen et al. | Mapping Congo Basin vegetation types from 300 m and 1 km multi-sensor time series for carbon stocks and forest areas estimation | |
| dela Torre et al. | Remote sensing-based estimation of rice yields using various models: A critical review | |
| Bagan et al. | Land cover classification from MODIS EVI times-series data using SOM neural network | |
| Wang et al. | On deep learning-based bias correction and downscaling of multiple climate models simulations | |
| CN106918816A (en) | crop growth monitoring system and its data processing method and module | |
| CN111652443A (en) | Method for predicting ocean fishery resource abundance by comprehensive multi-source satellite remote sensing and application thereof | |
| Gandharum et al. | Remote sensing versus the area sampling frame method in paddy rice acreage estimation in Indramayu regency, West Java province, Indonesia | |
| Ozdarici Ok et al. | A segment-based approach to classify agricultural lands by using multi-temporal optical and microwave data | |
| CN112836725A (en) | Weak supervision LSTM recurrent neural network rice field identification method based on time sequence remote sensing data | |
| CN114926743A (en) | Crop classification method and system based on dynamic time window | |
| CN114529097A (en) | Multi-scale crop phenological period remote sensing dimensionality reduction prediction method | |
| Kenduiywo et al. | Crop type mapping from a sequence of TerraSAR-X images with dynamic conditional random fields | |
| CN103034988B (en) | A kind of space-time quantitative remote sensing fusion method of any sensor quantity | |
| Zhang et al. | Merging multisatellite precipitation products using stacking method and the censored-shifted gamma ensemble model output statistics in china's Beimiaoji basin | |
| Shiu et al. | Mapping paddy rice agriculture in a highly fragmented area using a geographic information system object-based post classification process | |
| Kumari et al. | Object-based machine learning approach for soybean mapping using temporal sentinel-1/sentinel-2 data | |
| Guo et al. | Panicle Ratio Network: streamlining rice panicle measurement by deep learning with ultra-high-definition aerial images in the field | |
| Chawla | Possibilistic c-means-spatial contextual information based sub-pixel classification approach for multi-spectral data | |
| CN114782825B (en) | Crop identification method and device based on incomplete remote sensing data and electronic equipment | |
| CN112052720B (en) | High-space-time normalization vegetation index NDVI fusion model based on histogram clustering | |
| Misra et al. | Mapping specific crop-A multi sensor temporal approach | |
| Han et al. | Combining Sentinel-1 and-3 imagery for retrievals of regional multitemporal biophysical parameters under a deep learning framework |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190614 Termination date: 20200411 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |