GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression

Abstract

The present invention discloses a GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression, comprising the following steps: S1: Select data sources to obtain a first dataset of remote sensing images; the first dataset includes GIMMS NDVI 3g dataset and MOD13A2 NDVI dataset; S2: Perform filtering preprocessing on the first dataset to obtain a second dataset after filtering preprocessing; S3: For the data in the overlapping time period of the two remote sensing images in the second dataset, perform maximum value composition processing using the maximum value composition method to obtain a monthly resolution third dataset of the two types of data; S4: Establish a fusion model, and according to the third dataset, use point spread function simulation and time series regression to determine the unknown parameters of the fusion model, so as to obtain the final fusion model; S5: Use the final fusion model to process the MOD13A2 NDVI dataset after 2015, so as to obtain the GIMMS NDVI 3g dataset with an extended time period. The present invention can extend the time period of the GIMMS NDVI 3g dataset, and the extension accuracy is high and the effect is good.

Description

Technical Field

[0001] The invention relates to the technical field of ecological remote sensing, and in particular to a GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression. Background Art

[0002] Vegetation, as an important component of terrestrial ecosystems, plays a key role in the connection between the pedosphere, atmosphere and hydrosphere. It directly or indirectly affects the carbon cycle, water cycle and energy exchange of different ecosystems. With the development of satellite sensor technology, remote sensing technology has been widely used in vegetation growth monitoring due to its large coverage area and long time series. NDVI has been proven to have a good relationship with vegetation coverage, and it has been widely used in fields such as vegetation dynamic monitoring.

[0003] At present, the main NDVI data products include SPOT NDVI, MODIS NDVI and GIMMS NDVI. 3g Among them, SPOTNDVI is obtained by the VEGETATION sensor carried by the SPOT-4 satellite, with a spatial resolution of 1km, a temporal resolution of 10 days, and a time period from 1998 to the present; MODIS NDVI is a MODIS data series regularly released by NASA in the United States, with spatial resolutions of 250m, 500m, 1000m, etc., a temporal resolution of 16 days, etc., and a time period from 2000 to the present; GIMMS NDVI 3g It is the third generation NDVI dataset synthesized by the NASA Global Monitoring and Simulation Research Group using NOAA satellites, with a spatial resolution of 8km, a temporal resolution of half a month, and a time period of 1982-2015. Although these three NDVI data products are widely used around the world, their application is always limited by the time span. Summary of the invention

[0004] In view of the above problems, the present invention aims to provide a GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression.

[0005] The technical solution of the present invention is as follows:

[0006] A GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression includes the following steps:

[0007] S1: Select a data source to obtain a remote sensing image dataset 1; the dataset 1 includes GIMMS NDVI 3g dataset and MOD13A2 NDVI dataset;

[0008] S2: performing filtering preprocessing on the data set 1 to obtain a data set 2 after filtering preprocessing;

[0009] S3: For the data of the overlapping time periods of the two remote sensing images in the second data set, the maximum value synthesis method is used to perform maximum value synthesis processing to obtain the monthly resolution data set 3 of the two data;

[0010] S4: establishing a fusion model, according to the data set 3, using point spread function simulation and time series regression to determine unknown parameters of the fusion model, so as to obtain a final fusion model;

[0011] S5: Use the final fusion model to process the MOD13A2 NDVI dataset after 2015 to obtain the GIMMS NDVI after the time period is extended 3g Dataset.

[0012] Preferably, in step S2, Savitzky-Golay filtering is used for filtering preprocessing, and the Savitzky-Golay filtering is:

[0013]

[0014] Where: Y is the synthetic sequence data; f is half of the smoothing window size; C q is the filter coefficient; Y h+1 is the original sequence data; D is the size of the smoothing window.

[0015] Preferably, in step S3, the maximum value synthesis method is:

[0016]

[0017] Where: MNDVI is the maximum NDVI value of a certain pixel each month; e is the serial number of the first image of each month; n is the serial number of the last image of each month; NDVI g is the NDVI value of the g-th image.

[0018] Preferably, in step S4, the fusion model is:

[0019]

[0020] Where: GIMMS NDVI predicted by the fusion model 3g image; and are the values ​​of linear correction parameters a and b obtained through optimization; P n is the normalized two-dimensional Gaussian kernel; M is the original 1km MOD13A2 NDVI image;

[0021] The unknown parameters include With P n ,in and Calculate using the following formula:

[0022]

[0023] Where: and are the values ​​of linear correction parameters a and b obtained by optimization at the (x, y) coordinates; N represents the number of samples participating in the time series regression; represents the upscaled 8kmMOD13A2 NDVI image, (x, y) is its horizontal and vertical coordinates; G i (x,y) represents the monthly maximum GIMMS NDVI 3g data;

[0024] P n Calculate using the following formula:

[0025]

[0026] In the formula: r is the multiple of the scale increase; K is the excess Part of the space covered The number of fine pixels in Part of the space covered The index of the fine pixel in M; (i, j) is the horizontal and vertical coordinates of M pixels, (i o ,j o )for The corresponding position at fine spatial resolution; (i k ,j k ) is the kth The position of MOD13A2 pixel in the image; For The portion of space covered; To exceed The part of space covered; M p is an image block of image M.

[0027] As a preference, the formula used for time series regression is:

[0028]

[0029] Preferably, when simulating the point spread function, the upscaling model used is:

[0030]

[0031] K=(2S+r) 2 -r 2 (8)

[0032] Where: M t For the tth MOD13A2 pixel value in; M k For the kth MOD13A2 pixel value in; S is beyond The width of the portion of the space covered.

[0033] Preferably, the width S is determined by the following formula:

[0034]

[0035]

[0036] Where: is the value of S obtained by optimization; SSIM is the structural similarity index; G represents GIMMS NDVI 3g ;μ G and Represents G and The mean of; C1 and C2 are both constants; σ G , as well as G and The variance and covariance of .

[0037] The beneficial effects of the present invention are:

[0038] Based on point spread function simulation and time series regression, the present invention can solve the problem of MOD13A2 NDVI to GIMMS NDVI. 3g The spatial structure difference (i.e., scale difference) and radiation difference in the conversion process make the present invention extend the GIMMSNDVI 3g The time period of the data set is of high accuracy and good effect.

[0039] In order to make the upscaled MOD13A2 NDVI closer to the GIMMS NDVI in spatial structure 3g The present invention uses a normalized two-dimensional Gaussian kernel to simulate the point spread function, which not only considers the contribution of the fine pixels spatially covered by each coarse pixel, but also fully considers the contribution of the fine pixels around its spatial coverage range. Therefore, it can obtain a better quality upscaled MOD13A2 NDVI than the traditional method.

[0040] Eliminating MOD13A2 NDVI and GIMMS NDVI 3gWhen there is a radiation difference between the pixels, the present invention uses pixel-level time series regression to coordinate the data. This regression can eliminate the mutual interference between pixels and perform special regression learning for specific pixels, thereby effectively eliminating the radiation difference. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0042] Figure 1 It is a flow chart of the GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression of the present invention;

[0043] Figure 2 For a specific example, the target area is the GIMMS NDVI in January 2001. 3g Maximum spatial distribution map;

[0044] Figure 3 The spatial distribution map of the maximum value of MOD13A2 NDVI in the target area in January 2001 is a specific embodiment;

[0045] Figure 4 For a specific example, the GIMMS NDVI simulated for the target area on January 17, 2016 3g Spatial distribution map. DETAILED DESCRIPTION

[0046] The present invention is further described below in conjunction with the accompanying drawings and embodiments. It should be noted that, in the absence of conflict, the embodiments in this application and the technical features in the embodiments can be combined with each other. It should be noted that, unless otherwise specified, all technical and scientific terms used in this application have the same meanings as those generally understood by those of ordinary skill in the art to which this application belongs. The words "including" or "comprising" and the like used in the disclosure of the present invention mean that the elements or objects appearing before the word cover the elements or objects listed after the word and their equivalents, without excluding other elements or objects.

[0047] like Figure 1 As shown, the present invention provides a GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression, comprising the following steps:

[0048] S1: Select a data source to obtain a remote sensing image dataset 1; the dataset 1 includes GIMMS NDVI3g dataset and MOD13A2 NDVI dataset.

[0049] S2: Perform filtering preprocessing on the data set 1 to obtain a data set 2 after filtering preprocessing.

[0050] In a specific embodiment, Savitzky-Golay filtering is used for filtering preprocessing, and the Savitzky-Golay filtering is:

[0051]

[0052] Where: Y is the synthetic sequence data; f is half of the smoothing window size; C q is the filter coefficient; Y h+1 is the original sequence data; D is the size of the smoothing window, D=2f+1, which is also the number of convolutions.

[0053] It should be noted that the purpose of filtering preprocessing is to remove outliers in a data set. In addition to the filtering method used in the above embodiment, other filtering methods in the prior art may also be applicable to the present invention.

[0054] S3: The data of the overlapping time periods of the two remote sensing images in the second data set are subjected to maximum synthesis processing by using the maximum synthesis method to obtain a monthly resolution data set three of the two data.

[0055] In a specific embodiment, the maximum value synthesis method is:

[0056]

[0057] Where: MNDVI is the maximum NDVI value of a certain pixel each month; e is the serial number of the first image of each month; n is the serial number of the last image of each month; NDVI g is the NDVI value of the g-th image.

[0058] S4: Establish a fusion model, and determine the unknown parameters of the fusion model by using point spread function simulation and time series regression according to the data set 3, so as to obtain the final fusion model.

[0059] In a specific embodiment, the fusion model is:

[0060]

[0061] Where: GIMMS NDVI predicted by the fusion model 3g image; and are the values ​​of linear correction parameters a and b obtained through optimization; P n is the normalized two-dimensional Gaussian kernel; M is the original 1km MOD13A2 NDVI image;

[0062] The unknown parameters include With P n The present invention uses point spread function simulation and time series regression to determine the unknown parameters. MOD13A2 NDVI and GIMMS NDVI 3g The differences include scale differences, and because the two come from different sensors, there will be radiation differences caused by a series of complex factors such as imaging time differences, imaging angle differences, band setting differences, etc. In the present invention, point spread function simulation is used to eliminate scale differences, and time series regression is used to eliminate radiation differences.

[0063] The point spread function is a description of the complex spatial response relationship in the process of remote sensing image scale change. Its purpose is to truly describe the spatial correspondence between remote sensing images of different scales. Fully considering the influence of PSF in the upscaling task is conducive to improving the spatial structure accuracy of the remote sensing image after upscaling. In the upscaling task, the pixels (coarse pixels) of the coarse-scale image will be affected by the pixels (fine pixels) of the fine-scale image in its spatial coverage area, and will also be affected by the fine pixels around the range. And this influence is often large and cannot be ignored. The existing method uses mean filtering to complete upscaling, which essentially oversimplifies the real spatial response relationship, that is, oversimplifies the point spread function. In the present invention, a normalized two-dimensional Gaussian kernel is used to simulate the point spread function for upscaling of MOD13A2 NDVI, specifically:

[0064] First, a normalized two-dimensional Gaussian kernel is used to simulate the point spread function for upscaling the MOD13A2 NDVI. The process is shown in the following formula:

[0065]

[0066] Where: represents the upscaled 8kmMOD13A2 NDVI image, (x, y) is its horizontal and vertical coordinates; P n is the normalized two-dimensional Gaussian kernel; * is the convolution operation; M represents the original 1km MOD13A2 NDVI image, M p is an image block of image M, which is n Consistent size.

[0067] The 1kmMOD13A2 NDVI pixels in this image block will affect the upscaling process. The fine pixels beyond the block are considered to have no impact on it. p It consists of two parts, one of which is The portion of the space covered is defined here as The second is the part beyond its spatial coverage, defined as Its width is defined as S. Correspondingly, P is defined as:

[0068]

[0069] Where: r is the multiple of upscaling; (i, j) is the horizontal and vertical coordinates of M pixels, (i o ,j o )for The corresponding positions at fine spatial resolution;

[0070] Then the normalized form of P can be written as:

[0071]

[0072] Where: K is the excess Part of the space covered The number of fine pixels in Part of the space covered The index of the fine pixel in k ,j k ) is the kth The position of the MOD13A2 pixel in the image.

[0073] Combining formula (5) with formula (11), we can obtain the upscaling model shown below:

[0074]

[0075] K=(2S+r) 2 -r 2 (8)

[0076] Where: M t For the tth MOD13A2 pixel value in; M k For the kth The MOD13A2 pixel value in .

[0077] From formulas (7)-(8), it can be seen that there is only one unknown variable S in the upscaling model, which controls By determining the range of this parameter, the upscaling can be completed.

[0078] The purpose of upscaling is to unify MOD13A2 NDVI and GIMMS NDVI 3gThe spatial resolution between the data is increased to achieve the effect of similar spatial structure. Compared with the real GIMMS NDVI 3g The spatial structural similarity of (denoted as G) is taken as the goal and optimized to adaptively determine S. The adaptive determination of the parameters in the simulation by optimization enables the results to be flexibly applied to various research areas.

[0079] In a specific embodiment, the structural similarity index (Structure Similarity, SSIM) is used to measure the spatial structural similarity between and , and its formula is as follows:

[0080]

[0081] Where: μ G and Represents G and The mean of; C1 and C2 are both constants, which can avoid the SSIM value being 0; σ G , as well as G and The variance and covariance of .

[0082] The optimization process can be expressed as:

[0083]

[0084] Where: is the value of S obtained through optimization.

[0085] By optimizing the above formula for data set 3, we can finally determine Thus, an upscaling model with all parameters known is obtained, thus completing the upscaling.

[0086] In a specific embodiment, when performing time series regression, the formula used is:

[0087]

[0088] Using the above formula to perform time series regression on data set 3, the radiation difference can be eliminated through linear correction, so that Closer to G. When performing time series regression, the linear correction parameters a and b are determined by the following formula:

[0089]

[0090] Where: and are the values ​​of linear correction parameters a and b obtained by optimization at the (x, y) coordinates; N represents the number of samples participating in the time series regression; represents the upscaled 8kmMOD13A2 NDVI image, (x, y) is its horizontal and vertical coordinates; G i (x,y) represents the monthly maximum GIMMS NDVI 3g data;

[0091] S5: Use the final fusion model to process the MOD13A2 NDVI dataset after 2015 to obtain the GIMMS NDVI after the time period is extended 3g Dataset.

[0092] In a specific embodiment, taking the Qinghai-Tibet Plateau as an example, the GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression of the present invention is used to fusion the GIMMS NDVI of the study area. 3g The data time period is extended.

[0093] The Qinghai-Tibet Plateau (26°00′12″-39°46′50″N, 73°18′52″-104°46′59″E) is located in southwestern China, with an area of ​​about 2.57×10 6 km 2 The Qinghai-Tibet Plateau stretches from the Pamir Plateau and the Karakoram Mountains in the west to the Yulong Snow Mountain, the Daxue Mountain, the Jiajin Mountain, the Qionglai Mountain and the Min Mountain in the east, the southern edge of the Himalayas in the south, and the Kunlun Mountain, the Altun Mountain and the Qilian Mountain in the north. The average altitude of the Qinghai-Tibet Plateau is about 4,000 meters. There are many high mountains in the region, so it is also called the "Roof of the World", the Asian Water Tower and the Third Pole of the Earth. In terms of administrative divisions, the Qinghai-Tibet Plateau involves six provinces and regions, namely the Tibet Autonomous Region, Qinghai Province, Sichuan Province, Gansu Province, Yunnan Province and the Xinjiang Uygur Autonomous Region, with a total of 201 counties (cities).

[0094] In this embodiment, the data set 1 includes the GIMMS NDVI with a spatial resolution of 8 km, a temporal resolution of half a month, and a time period of 1982-2015. 3g The data set has a spatial resolution of 1 km, a temporal resolution of 16 days, and a time period of 2001 to the present MOD13A2 NDVI data set. Thus, when the maximum value synthesis method is used to perform maximum value synthesis processing on the data of the overlapping time period of the two remote sensing images in step S3, the maximum value synthesis method is used to perform maximum value synthesis processing on the data from 2001 to 2015. In this embodiment, GIMMS NDVI 3g The maximum composite processing results of NDVI of MOD13A2 and MOD13A2 are as follows: Figure 2 and Figure 3The fusion model obtained in this embodiment is used to process the MOD13A2 NDVI dataset after 2015, and the GIMMS NDVI 3g The time period is extended, and the GIMMS NDVI simulated on January 17, 2016 3g The results are as follows Figure 4 shown.

[0095] In addition, in this embodiment, the GIMMS NDVI with a temporal resolution of half a month is also converted according to the overlapping time period (2001-2015). 3g The maximum value of NDVI for each month is calculated by the maximum value synthesis method with the 16-day simulation data of the time resolution, and then the two are compared and verified to obtain the simulation accuracy of the present invention. The experimental results show that the present invention can extend the GIMMS NDVI with high accuracy. 3g time period.

[0096] In summary, the present invention can accurately measure the GIMMS NDVI 3g The data set is extended over a time period. Compared with the prior art, the present invention has significant improvements.

[0097] The above description is only a preferred embodiment of the present invention and does not constitute any form of limitation to the present invention. Although the present invention has been disclosed as a preferred embodiment as above, it is not intended to limit the present invention. Any technician familiar with the profession can make some changes or modifications to equivalent embodiments of equivalent changes using the technical contents disclosed above without departing from the scope of the technical solution of the present invention. However, any simple modification, equivalent change and modification made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solution of the present invention still falls within the scope of the technical solution of the present invention.

Claims

1. A GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression, characterized in that: The following steps are involved: S1: Select a data source to obtain a remote sensing image dataset 1; the dataset 1 includes GIMMS NDVI 3g dataset and MOD13A2 NDVI dataset; S2: performing filtering preprocessing on the data set 1 to obtain a data set 2 after filtering preprocessing; S3: For the data of the overlapping time periods of the two remote sensing images in the second data set, the maximum value synthesis method is used to perform maximum value synthesis processing to obtain the monthly resolution data set 3 of the two data; S4: Establish a fusion model. According to the data set 3, use point spread function simulation and time series regression to determine the unknown parameters of the fusion model, so as to obtain the final fusion model; the fusion model is: Where: GIMMS NDVI predicted by the fusion model 3g image; and are the values ​​of linear correction parameters a and b obtained through optimization; P n is the normalized two-dimensional Gaussian kernel; M is the original 1km MOD13A2 NDVI image; The unknown parameters include With P n ,in and Calculate using the following formula: Where: and are the values ​​of linear correction parameters a and b obtained by optimization at the (x, y) coordinates; N represents the number of samples participating in the time series regression; represents the upscaled 8km MOD13A2 NDVI image, (x, y) is its horizontal and vertical coordinates; G i (x,y) represents the monthly maximum GIMMS NDVI 3g data; P n Calculate using the following formula: In the formula: r is the multiple of the scale increase; K is the excess Part of the space covered The number of fine pixels in Part of the space covered The index of the fine pixel in ; (i,j) is the horizontal and vertical coordinates of M pixels, (i o ,j o )for The corresponding position at fine spatial resolution; (i k ,j k ) is the kth The position of MOD13A2 pixel in the image; For The portion of space covered; To exceed The part of space covered; M p is an image block of image M; When simulating the point spread function, the upscaling model used is: K=(2S+r) 2 -r 2 (8) Where: M t For the tth MOD13A2 pixel value in; M k For the kth MOD13A2 pixel value in; S is beyond the width of the portion of the space covered; S5: Use the final fusion model to process the MOD13A2 NDVI dataset after 2015 to obtain the GIMMS NDVI after the time period is extended 3g Dataset.

2. The GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression according to claim 1, characterized in that: In step S2, Savitzky-Golay filtering is used for filtering preprocessing, and the Savitzky-Golay filtering is: Where: Y is the synthetic sequence data; f is half of the smoothing window size; C q is the filter coefficient; Y h+1 is the original sequence data; D is the size of the smoothing window.

3. The GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression according to claim 1, characterized in that, In step S3, the maximum value synthesis method is: Where: MNDVI is the maximum NDVI value of a certain pixel each month; e is the serial number of the first image of each month; n is the serial number of the last image of each month; NDVI g is the NDVI value of the g-th image.

4. The GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression according to claim 1, characterized in that: When performing time series regression, the formula used is:

5. The GIMMS-MODIS NDVI data fusion method based on point spread function simulation and time series regression according to claim 1, characterized in that, The width S is determined by the following formula: Where: is the value of S obtained through optimization; SSIM is the structural similarity index; G stands for GIMMS NDVI 3g ;μ G and Represents G and The mean of; C1 and C2 are both constants; σ G , as well as G and The variance and covariance of .

Images