MODIS image filtering method based on spatiotemporal multi-dimensional fusion technology of vegetation index

Abstract

The present invention provides a MODIS image filtering method based on the spatiotemporal multi-dimensional fusion technology of vegetation index, which solves the defects of irregular changes in NDVI time series, unreliable information extraction during long-term downward trend, continuous mutation points, and inability to consider the spatial domain compared with the prior art. The present invention includes the following steps: acquisition of MODIS pixel data; preprocessing of MODIS pixel data; time dimension filtering of the preprocessed MODIS NDVI band I; space dimension filtering of the preprocessed MODIS NDVI band I; and acquisition of MODIS image filtering results. The present invention improves the smoothness of time series data through the time dimension, eliminates the uncertainty of single-dimensional filtering data, eliminates the influence of continuous mutation values, increases the credibility of data, eliminates the influence of irregularly changing NDVI time series, and provides more realistic filtering data.

Description

Technical Field

[0001] The invention relates to remote sensing image processing technology, in particular to a MODIS image filtering method based on vegetation index spatiotemporal multi-dimensional fusion technology. Background Art

[0002] The time series data of Normalized Difference Vegetation Index (NDVI) reflects the dynamic change process of vegetation. However, due to the influence of factors such as atmosphere, cloud cover, and snow cover, there is a lot of noise in the remote sensing vegetation index image, which limits the in-depth analysis and application of the data. At present, there are more than 20 methods for time series data reconstruction (filtering, denoising), mainly including Best Index Slope Extraction (BISE), Median Filtering (MF), SG method (Savitzky Golay, SG), Fourier Transform (FT), Wavelet Transform (WT), Interpolation for Data Reconstruction (IDR), Mean Value Iteration filter (MVI), Maximum Value Composite (MVC), etc.

[0003] The above time series data reconstruction methods can generally be divided into three types, each of which has advantages and disadvantages in reconstructing NDVI time series data. The first is a local filtering method in the time domain, such as BISE. This method is based on the predictability of vegetation changes. The filtered results are often close to the original NDVI time series, but it requires a sliding period to determine the NDVI regeneration percentage, and it is also necessary to determine a smooth area to indicate whether the increase in the area is within an acceptable range. Therefore, when the NDVI time series shows a long-term downward trend, this method makes the extracted information unreliable. The second is a frequency domain denoising method, such as FT. This method fits the first three harmonics with a least squares solver. The NDVI curve obtained by this method is relatively smooth, and the vegetation growth characteristics are obvious. However, it is not suitable for irregular or asymmetric NDVI data. The third is a function-based method, which is based on the least squares of the piecewise function of the NDVI time series, such as the dual logistic function fitting method. However, this method is difficult to fit the NDVI time series with irregular seasonal time distribution, and its application is limited when the area is large. It should be pointed out that the vegetation index time series is a kind of time series data with strong periodicity, but due to cloud and snow occlusion, a certain number of mutation points will be generated. These mutation points are high-frequency components in the frequency domain. Therefore, proper smoothing of such high-frequency components is the key to denoising the vegetation index time series data. Wavelet transform filtering is an effective frequency domain filtering algorithm. It inherits and develops the localization idea of ​​short-time Fourier transform and overcomes certain shortcomings. However, the filtering effect is still not ideal for mutation points or continuous mutation points with obvious changes in vegetation index, and further optimization is needed. In addition, the above method can only focus on the change characteristics of NDVI in the time-frequency dimension, and cannot consider the spatial correlation of NDVI.

[0004] Therefore, how to design a new MODIS image filtering method to solve the problems of irregular changes in NDVI time series, unreliable information extraction during long-term downward trends, continuous mutation points, and inability to consider the spatial domain. Summary of the invention

[0005] The purpose of the present invention is to solve the defects of the prior art, such as irregular changes in NDVI time series, unreliable information extraction during long-term downward trends, continuous mutation points, and inability to consider the spatial domain, and to provide a MODIS image filtering method based on spatiotemporal multidimensional fusion technology of vegetation index to solve the above problems.

[0006] In order to achieve the above object, the technical solution of the present invention is as follows:

[0007] A MODIS image filtering method based on vegetation index spatiotemporal multi-dimensional fusion technology includes the following steps:

[0008] Acquisition of MODIS pixel data: Extract MODIS NDVI band I and pixel quality band II data from NASA;

[0009] Preprocessing of MODIS pixel data: Grouping and preprocessing of MODIS NDVI band I and pixel quality band II data;

[0010] The pre-processed MODIS NDVI band I is subjected to time dimension filtering: weighted interpolation and wavelet transform filtering are performed on the grouped data of band I in the time dimension;

[0011] The pre-processed MODIS NDVI band I is subjected to spatial filtering: in the spatial dimension, the low-quality pixels of the grouped data of band I are spatially filtered using the pixel quality information of band II and the neighboring pixels;

[0012] Obtaining MODIS image filtering results: Weighted fusion of time dimension and space dimension filtering data is performed to obtain MODIS image filtering results.

[0013] The acquisition of the MODIS pixel data comprises the following steps:

[0014] Band I and pixel quality band II of NASA's MOD13Q1 data product were extracted. The extraction methods included image stitching, projection, conversion, and cropping to obtain the monthly NDVI images of the study period. The pixels of each period of the extracted pixel quality band II corresponded to each period of NDVI band I. The quality of each pixel was divided into five levels: -1 level means no value, 0 or 1 level means valid value, and 2 or 3 level means invalid value.

[0015] All image pixel values ​​are multiplied by a coefficient of 0.0001 to control the NDVI interval to [-1,1], where the [-1,0] interval represents bare surface and water areas with no vegetation; the NDVI in the [0,1] interval is retained.

[0016] The preprocessing of the MODIS pixel data includes the following steps:

[0017] Preprocessing of MODIS NDVI band I is used to group NDVI images during the study period:

[0018] NDVI group one is a two-dimensional table of one-dimensional time series of the same coordinate pixel of each period of research image. The number of research image periods is equal to the length of the one-dimensional time series, and the number of one-dimensional time series is equal to the number of research image pixels.

[0019] NDVI group 2 is a three-dimensional table containing all research images. Each layer of the table is a research image of each period. The height of the table is equal to the number of research image periods.

[0020] Among them, the two-dimensional table of one-dimensional time series is the result of flattening all the research images. The flattening method is to arrange each pixel in each period of image from top to bottom and from left to right to form a column of data. Each column of data is a period of research image. After all the research images are flattened, a two-dimensional table is formed. Each row of the two-dimensional table is the one-dimensional time series mentioned.

[0021] Pixel quality band II is divided into two groups. Pixel quality group one is a two-dimensional table of one-dimensional time series containing all pixel quality band images, and pixel quality group two is a three-dimensional table containing all pixel quality band images.

[0022] The time dimension filtering process for the pre-processed MODIS NDVI band I comprises the following steps:

[0023] Traverse the pixel quality group one corresponding to the NDVI group one;

[0024] If the quality level of a pixel in pixel quality group 1 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 1 is retained;

[0025] If the quality level of a pixel in pixel quality group 1 is 2 or 3, the interpolated pixel value is obtained by weighting the pixel values ​​before and after the current pixel in the corresponding NDVI group 1. The calculation formula is as follows:

[0026]

[0027] m=|NDVI j，i -NDVI j+1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0028] m=|NDVI j，i -NDVI j-1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0029] In the formula, is the interpolated NDVI of the i-th pixel in the j-th period, NDVI j，i is the original NDVI of the i-th pixel in the j-th period, NDVI j+1，i is the original NDVI of the i-th pixel in the j+1th period, NDVI j-1，iis the original NDVI of the i-th pixel in the j-1th period, PixRe is the pixel quality level, m and n are the weights during interpolation, j is an image of a certain period, and i is a pixel in the image;

[0030] For the interpolated one-dimensional time series NDVI group one, wavelet transform filtering is performed. The calculation formula of wavelet transform filtering is:

[0031]

[0032] Where WT(α, τ) is the result of wavelet transform filtering under a certain scale and translation, α is the scale, τ is the translation, is the wavelet function, t is the time domain, and f(t) is the signal;

[0033] The NDVI group that has completed time dimension filtering is marked with TNDVI.

[0034] The spatial dimension filtering process of the pre-processed MODIS NDVI band I comprises the following steps:

[0035] Traverse the pixel quality group 2 corresponding to NDVI group 2:

[0036] If the quality level of a pixel in pixel quality group 2 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 2 is retained;

[0037] If the quality level of a pixel in pixel quality group 2 is 2 or 3, the 8 neighboring pixels corresponding to the current pixel position in NDVI group 2 are sorted from small to large, and the median value of the 8 pixel values ​​is used to replace the current pixel value;

[0038] The traversal method for NDVI group 2 is set as follows: generate a 3x3 filter, start from the upper left corner of each image, slide from left to right and from top to bottom, and retrieve the value of each pixel in the center of the 3x3 filter;

[0039] The NDVI group 2 that has completed spatial filtering is marked with SNDVI.

[0040] The acquisition of the MODIS image filtering result comprises the following steps:

[0041] Calculate the corresponding weights p and 1-p of TNDVI and SNDVI,

[0042] That is, if the quality level corresponding to a pixel at a certain position in a certain period is 2 or 3, then the number of periods num with quality levels 2 or 3 corresponding to the pixels at the same position in all periods snum is counted, and the p value is num / snum. The corresponding 1-p is 1-num / snum, and the weighted sum calculation formula is:

[0043] NDVI st =(1-p)·SNVDVI+p·TNDVI

[0044] Where NDVI st is the NDVI time series obtained by filtering;

[0045] Output NDVI st This is the filtered NDVI time series image.

[0046] Beneficial Effects

[0047] Compared with the prior art, the MODIS image filtering method based on the spatiotemporal multidimensional fusion technology of vegetation index of the present invention improves the smoothness of time series data through the time dimension, eliminates the problem of being unable to consider the spatial domain through the space dimension, improves the authenticity of time series data by combining time and space, eliminates the uncertainty of single-dimensional filtering data, eliminates the influence of continuous mutation values, increases the credibility of data, eliminates the influence of irregularly changing NDVI time series, and provides more realistic filtering data.

[0048] The present invention can restore the original value of the data to the greatest extent while ensuring the periodicity and trend of the vegetation index data, reduce the fluctuation of the time series data, weaken the fault phenomenon of the image, and make the data smoother and more continuous. Compared with the complex parameter requirements of other reconstruction methods, the present invention is more operational. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 It is a schematic diagram of the process of the present invention;

[0050] Figure 2a NDVI time series graph of checkpoints for reconstructing NDVI images using other reconstruction methods;

[0051] Figure 2b A time series diagram of the NDVI checkpoints for reconstructing the NDVI image of the present invention;

[0052] Figure 3a This is the original image of the NDVI image;

[0053] Figure 3b It is the noise map of NDVI image;

[0054] Figure 3c The NDVI image map is reconstructed for the present invention. DETAILED DESCRIPTION

[0055] In order to have a further understanding and recognition of the structural features and the effects achieved by the present invention, a preferred embodiment and accompanying drawings are used for detailed description as follows:

[0056] like Figure 1As shown, the MODIS image filtering method based on spatiotemporal multi-dimensional fusion technology of vegetation index described in the present invention has superior effect, lower requirements on input data than existing methods, simple data preprocessing, simple parameter setting, no complex expression in mathematical calculation, fast calculation and high execution efficiency. It includes the following steps:

[0057] The first step is to obtain MODIS pixel data: extract MODIS NDVI band I and pixel quality band II data from NASA.

[0058] (1) Band I and pixel quality band II of the MOD13Q1 data product from NASA were extracted by stitching, projecting, converting, and cropping the images to obtain the monthly NDVI images of the study period. The pixels of each period of the extracted pixel quality band II corresponded to the NDVI band I of each period. The quality of each pixel was divided into five levels: -1 for no value, 0 or 1 for valid values, and 2 or 3 for invalid values.

[0059] (2) Multiply all image pixel values ​​by a coefficient of 0.0001 to keep the NDVI interval within [-1,1], where the interval [-1,0] represents bare surface and water with no vegetation; retain the NDVI in the interval [0,1].

[0060] The second step is the preprocessing of MODIS pixel data: grouping and preprocessing of MODIS NDVI band I and pixel quality band II data.

[0061] (1) Preprocessing of MODIS NDVI band I to group NDVI images during the study period:

[0062] NDVI group one is a two-dimensional table of one-dimensional time series of the same coordinate pixel of each period of research image. The number of research image periods is equal to the length of the one-dimensional time series, and the number of one-dimensional time series is equal to the number of research image pixels.

[0063] NDVI group 2 is a three-dimensional table containing all research images. Each layer of the table is a research image of each period. The height of the table is equal to the number of research image periods.

[0064] Among them, the two-dimensional table of the one-dimensional time series is the result of flattening all the research images. The flattening method is to arrange each pixel in each period of image in order from top to bottom and from left to right to form a column of data. Each column of data is a period of research image. After all the research images are flattened, a two-dimensional table is formed. Each row of the two-dimensional table is the one-dimensional time series mentioned.

[0065] (2) Pixel quality band II is divided into two groups. Pixel quality group 1 is a two-dimensional table of one-dimensional time series containing all pixel quality band images, and pixel quality group 2 is a three-dimensional table containing all pixel quality band images.

[0066] The third step is to perform time dimension filtering on the pre-processed MODIS NDVI band I: weighted interpolation and wavelet transform filtering are performed on the grouped data of band I in the time dimension. The specific steps are as follows:

[0067] (1) Traverse the pixel quality group 1 corresponding to the NDVI group 1;

[0068] If the quality level of a pixel in pixel quality group 1 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 1 is retained;

[0069] If the quality level of a pixel in pixel quality group 1 is 2 or 3, the interpolated pixel value is obtained by weighting the pixel values ​​before and after the current pixel in the corresponding NDVI group 1. The calculation formula is as follows:

[0070]

[0071] m=|NDVI j，i -NDVI j+1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0072] m=|NDVI j，i -NDVI j-1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0073] In the formula, is the interpolated NDVI of the i-th pixel in the j-th period, NDVI j，i is the original NDVI of the i-th pixel in the j-th period, NDVI j+1，i is the original NDVI of the i-th pixel in the j+1th period, NDVI j-1，i is the original NDVI of the i-th pixel in the j-1th period, PixRe is the pixel quality level, m and n are the weights during interpolation, j is an image of a certain period, and i is a pixel in the image.

[0074] (2) For the interpolated one-dimensional time series NDVI group one, wavelet transform filtering is performed. The calculation formula of wavelet transform filtering is:

[0075]

[0076] Where WT(α, τ) is the result of wavelet transform filtering under a certain scale and translation, α is the scale, τ is the translation, is the wavelet function, t is the time domain, and f(t) is the signal.

[0077] (3) The NDVI grouping after time dimension filtering is marked with TNDVI.

[0078] The fourth step is to perform spatial filtering on the pre-processed MODIS NDVI band I: In the spatial dimension, the low-quality pixels of the band I grouped data are spatially filtered using the band II pixel quality information and neighboring pixels. The specific steps are as follows:

[0079] (1) Traverse the pixel quality group 2 corresponding to NDVI group 2:

[0080] If the quality level of a pixel in pixel quality group 2 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 2 is retained;

[0081] If the quality level of a pixel in pixel quality group 2 is 2 or 3, the 8 neighboring pixels corresponding to the current pixel position in NDVI group 2 are sorted from small to large, and the median value of the 8 pixel values ​​is used to replace the current pixel value.

[0082] (2) The traversal method for NDVI group 2 is set as follows: generate a 3x3 filter, start from the upper left corner of each image, slide from left to right and from top to bottom, and retrieve the value of each pixel at the center of the 3x3 filter.

[0083] (3) The NDVI grouping after spatial filtering is marked with SNDVI.

[0084] The fifth step is to obtain the MODIS image filtering results: weighted fusion of the time dimension and space dimension filtering data to obtain the MODIS image filtering results.

[0085] (1) Calculate the corresponding weights p and 1-p of TNDVI and SNDVI,

[0086] That is, if the quality level corresponding to a pixel at a certain position in a certain period is 2 or 3, then the number of periods num with quality levels 2 or 3 corresponding to the pixels at the same position in all periods snum is counted, and the p value is num / snum. The corresponding 1-p is 1-num / snum, and the weighted sum calculation formula is:

[0087] NDVI st =(1-p)·SNVDVI+p·TNDVI

[0088] Where NDVI st is the NDVI time series obtained by filtering.

[0089] (2) Output NDVI st This is the filtered NDVI time series image.

[0090] Here, the following embodiments are provided:

[0091] 1. Extracted data: The extracted MODIS NDVI band I and pixel quality band II are from the MOD13Q1 data product of NASA (Nationalaeronautics and space administration) (http: / / search.earthdata.nasa.gov / search). The product has a temporal resolution of 16 days and a spatial resolution of 250 m. It includes 12 bands, the first band being the NDVI band and the twelfth band being the pixel credibility band (pixel quality band). The research area of ​​the present invention is the southern Anhui Province, and the research period is 2010-2011.

[0092] The extraction method is to use MRT (Modis reprojection tools) and ArcGIS tools to stitch, project, convert and crop the images, and use the MVC (Maximum value composites) maximization synthesis method to obtain the monthly NDVI images of the study period, so there are 46 images in 2010-2011. According to the MOD13Q1 product description, all image pixel values ​​are multiplied by a coefficient of 0.0001 to control the NDVI interval to [-1, 1], where the [-1, 0] interval is the exposed surface and water area, without vegetation distribution, so only the NDVI in the [0, 1] interval is retained.

[0093] In order to better demonstrate the model effect, the NDVI extracted here is used as the true value, and random noise points are added on the basis of this true value as the data for adding noise. The pixels of each period of pixel quality band II correspond to each period of NDVI band I. The quality of each pixel is divided into 5 levels: -1 level is no value, generally for bare surface; 0 or 1 level is a valid value and can be used directly; 2 or 3 level is an invalid value, generally for the surface covered by clouds and snow. When adding noise points, the quality level of the same position of the noise point on the extracted pixel quality image is randomly changed to 2 or 3.

[0094] 2. Preprocessing of MODIS pixel data: Group and preprocess MODIS NDVI band I and pixel quality band II data.

[0095] The band I images within the study period are grouped, and group one is a two-dimensional table of one-dimensional time series of the same coordinate pixels of the study images of each period. The number of study image periods is equal to the length of the one-dimensional time series, and the number of one-dimensional time series is equal to the number of study image pixels.

[0096] Group 2 is a three-dimensional table containing band I images. Each layer of the table is each period of research images. The height of the table is equal to the number of research image periods. This embodiment is implemented through MATLAB's reshape(A, sz) method.

[0097] The band II images during the study period were grouped in the same way as band I.

[0098] 3. Perform time dimension filtering on the preprocessed MODIS NDVI band I: perform weighted interpolation and wavelet transform filtering on the grouped data of band I in the time dimension.

[0099] Traverse the pixel quality group 1 corresponding to NDVI group 1. If the quality level of a pixel in pixel quality group 1 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 1 is retained. If the quality level of a pixel in pixel quality group 1 is 2 or 3, the pixel value before and after the current pixel in the corresponding NDVI group 1 is weighted and calculated. The calculation formula is:

[0100]

[0101] m=|NDVI j，i -NDVI j+1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0102] m=|NDVI j，i -NDVI j-1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |)

[0103] In the formula is the interpolated NDVI of the i-th pixel in the j-th period, NDVI j，i is the original NDVI of the i-th pixel in the j-th period, PixRe is the pixel quality level, and m and n are the weights during interpolation.

[0104] For the interpolated one-dimensional time series NDVI group one, wavelet transform filtering is performed. The calculation formula of wavelet transform filtering is:

[0105]

[0106] Where α is the scale, τ is the translation, is a wavelet function, f(t) is a signal, and this embodiment uses the wdenoise(X, LEVEL, 'DenoisingMethod', DMETHOD) method in MATLAB to complete the wavelet transform filtering of group 1. The NDVI group 1 that completes the time dimension filtering is marked with TNDVI.

[0107] 4. Perform spatial filtering on the preprocessed MODIS NDVI band I: In the spatial dimension, use the band II pixel quality information and neighboring pixels to perform spatial filtering on the low-quality pixels of the band I grouped data.

[0108] Traverse the pixel quality group 2 corresponding to NDVI group 2. If the quality level of a pixel in pixel quality group 2 is -1, 0 or 1, retain the original value of the pixel in the corresponding NDVI group 2. If the quality level of a pixel in pixel quality group 2 is 2 or 3, sort the 8 neighboring pixels corresponding to the current pixel position in NDVI group 2 from small to large, and take the median value of the 8 pixel values ​​to replace the current pixel value.

[0109] For the traversal of NDVI group 2, the specific traversal method is: generate a 3x3 filter, start from the upper left corner of each image, slide from left to right and from top to bottom, and retrieve each pixel value at the center of the 3x3 filter.

[0110] The filter described above can also be set to 4x4, 5x5 and other widths, and the number of neighboring pixels will change accordingly. It should be pointed out that the filter width directly affects the image filtering effect and filtering speed. The larger the width, the worse the filtering effect and the faster the filtering speed. The smaller the width, the better the filtering effect and the slower the filtering speed. In practical applications, the optimal filter size should be selected according to the situation.

[0111] The NDVI group 2 that has completed spatial filtering is marked with SNDVI.

[0112] 5. Obtaining MODIS image filtering results: weighted fusion of time dimension and space dimension filtering data.

[0113] Calculate the corresponding weights p and 1-p of TNDVI and SNDVI. Specifically, if the quality level corresponding to the pixel at a certain position in a certain period is 2 or 3, then count the number of periods num with quality levels 2 or 3 corresponding to the pixels at the same position in all periods snum. The p value is num / snum, and the corresponding 1-p can also be obtained. The weighted sum calculation formula is:

[0114] NDVI st=(1-p)·SNVDVI+p·TNDVI

[0115] Where NDVI st is the NDVI time series obtained by filtering.

[0116] Output NDVI st This is the filtered NDVI time series image.

[0117] Below Figure 2a , Figure 2b , Figure 3a to Figure 3c Table 1 is used as an example to illustrate the superiority of this method:

[0118] Select a checkpoint to compare the NDVI time series after filtering by the five reconstruction methods. Figure 2a As shown in the figure, NDVI is the true value, NDVINOISE is the noise value, FT is the Fourier transform method, MVC is the maximum value synthesis method, MF is the median filter method, SG is the Savitzky Golay method, and WT is the wavelet transform method. It can be seen that the effects of the five reconstruction methods are not very rational, and there is a large gap with the real NDVI curve. The NDVI time series reconstructed by using the present invention after selecting the same checkpoints as above is shown in the figure. Figure 2b As shown, it can be seen that the NDVI time series curve after filtering by the present invention (STFW) is closer to the real NDVI data.

[0119] The original image of a certain period in the study area is as follows Figure 3a As shown, the noise image is Figure 3b As shown, the image reconstructed using the present invention is as follows Figure 3c As shown, in order to reflect the NDVI image reconstruction ability of the present invention, it is compared with other 5 commonly used reconstruction methods, including Fourier transform method FT, maximum value synthesis method MVC, median filter method MF, SG method, and wavelet transform method WT. Two commonly used evaluation indicators are used, namely peak signal-to-noise ratio PNSR and structural similarity SSIM. The larger the PNSR and SSIM values, the better the reconstruction effect.

[0120] The reconstruction effects of various methods are shown in Table 1. STFW is the method proposed by the present invention. From Table 1, it can be seen that the STFW method is superior to other reconstruction methods in various indicators.

[0121] Table 1 is a comparison table of the reconstruction effects of the present invention and other methods

[0122] Reconstruction methods PNSR SSIM FT 67.4290 0.7009 MVC 67.1105 0.7561 MF 72.7094 0.6559 SG 83.8940 0.8966 WT 60.9354 0.3273 STFW 87.4731 0.9686

[0123] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions only describe the principles of the present invention. The present invention may be subject to various changes and improvements without departing from the spirit and scope of the present invention. These changes and improvements fall within the scope of the present invention. The scope of protection claimed by the present invention is defined by the attached claims and their equivalents.

Claims

1. A MODIS image filtering method based on spatiotemporal multidimensional fusion technology of vegetation index, characterized in that: The following steps are involved: 11) Acquisition of MODIS pixel data: Extract MODIS NDVI band I and pixel quality band II data from NASA; 12) Preprocessing of MODIS pixel data: Grouping and preprocessing of MODIS NDVI band I and pixel quality band II data; The preprocessing of the MODIS pixel data includes the following steps: 121) Preprocessing of MODIS NDVI band I to group NDVI images during the study period: NDVI group one is a two-dimensional table of one-dimensional time series of the same coordinate pixel of each period of research image. The number of research image periods is equal to the length of the one-dimensional time series, and the number of one-dimensional time series is equal to the number of research image pixels. NDVI group 2 is a three-dimensional table containing all research images. Each layer of the table is a research image of each period. The height of the table is equal to the number of research image periods. Among them, the two-dimensional table of one-dimensional time series is the result of flattening all the research images. The flattening method is to arrange each pixel in each period of image from top to bottom and from left to right to form a column of data. Each column of data is a period of research image. After all the research images are flattened, a two-dimensional table is formed. Each row of the two-dimensional table is the one-dimensional time series mentioned. 122) Divide the pixel quality band II into two groups, the pixel quality group 1 is a two-dimensional table of one-dimensional time series containing all pixel quality band images, and the pixel quality group 2 is a three-dimensional table containing all pixel quality band images; 13) Perform time dimension filtering on the pre-processed MODIS NDVI band I: perform weighted interpolation and wavelet transform filtering on the grouped data of band I in the time dimension; The time dimension filtering process for the pre-processed MODIS NDVI band I comprises the following steps: 131) traverse the pixel quality group 1 corresponding to the NDVI group 1; If the quality level of a pixel in pixel quality group 1 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 1 is retained; If the quality level of a pixel in pixel quality group 1 is 2 or 3, the interpolated pixel value is obtained by weighting the pixel values ​​before and after the current pixel in the corresponding NDVI group 1. The calculation formula is as follows: m=|NDVI j，i -NDVI j+1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |) m=|NDVI j，i -NDVI j-1，i | / (|NDVI j，i -NDVI j-1，i |+|NDVI j，i -NDVI j+1，i |) In the formula, is the interpolated NDVI of the i-th pixel in the j-th period, NDVI j，i is the original NDVI of the i-th pixel in the j-th period, NDVI j+1，i is the original NDVI of the i-th pixel in the j+1th period, NDVI j-1，i is the original NDVI of the i-th pixel in the j-1th period, PixRe is the pixel quality level, m and n are the weights during interpolation, j is an image of a certain period, and i is a pixel in the image; 132) For the interpolated one-dimensional time series NDVI group one, wavelet transform filtering is performed. The calculation formula of wavelet transform filtering is: Where WT(α,τ) is the result of wavelet transform filtering under a certain scale and translation, α is the scale, τ is the translation, is the wavelet function, t is the time domain, and f(t) is the signal; 133) completing the NDVI grouping after time dimension filtering and marking it with TNDVI; 14) Perform spatial filtering on the pre-processed MODIS NDVI band I: Use the band II pixel quality information and neighboring pixels in the spatial dimension to perform spatial filtering on the low-quality pixels of the band I grouped data; 15) Obtaining MODIS image filtering results: Weighted fusion of time dimension and space dimension filtering data to obtain MODIS image filtering results.

2. The MODIS image filtering method based on the spatiotemporal multidimensional fusion technology of vegetation index according to claim 1 is characterized in that: The acquisition of the MODIS pixel data comprises the following steps: 21) Extract band I and pixel quality band II from NASA's MOD13Q1 data product. The extraction method includes image stitching, projection, conversion, and cropping to obtain the monthly NDVI images of the study period. The pixels of each period of the extracted pixel quality band II correspond to each period of NDVI band I. The quality of each pixel is divided into 5 levels: -1 level means no value, 0 or 1 level means valid value, and 2 or 3 level means invalid value. 22) Multiply all image pixel values ​​by a coefficient of 0.0001 to control the NDVI interval to [-1,1], where the [-1,0] interval represents bare surface and water areas with no vegetation; retain the NDVI in the [0,1] interval.

3. The MODIS image filtering method based on the spatiotemporal multidimensional fusion technology of vegetation index according to claim 1 is characterized in that: The spatial dimension filtering process of the pre-processed MODIS NDVI band I comprises the following steps: 31) Traverse the pixel quality group 2 corresponding to NDVI group 2: If the quality level of a pixel in pixel quality group 2 is -1, 0 or 1, the original value of the pixel in the corresponding NDVI group 2 is retained; If the quality level of a pixel in pixel quality group 2 is 2 or 3, the 8 neighboring pixels corresponding to the current pixel position in NDVI group 2 are sorted from small to large, and the median value of the 8 pixel values ​​is used to replace the current pixel value; 32) The traversal method for NDVI group 2 is set as follows: generate a 3x3 filter, start from the upper left corner of each image, slide from left to right and from top to bottom, and retrieve the value of each pixel in the center of the 3x3 filter; 33) The NDVI grouping after spatial filtering is marked with SNDVI.

4. The MODIS image filtering method based on the spatiotemporal multidimensional fusion technology of vegetation index according to claim 1 is characterized in that: The acquisition of the MODIS image filtering result comprises the following steps: 41) Calculate the corresponding weights p and 1-p of TNDVI and SNDVI, That is, if the quality level corresponding to a pixel at a certain position in a certain period is 2 or 3, then the number of periods num with quality levels 2 or 3 corresponding to the pixels at the same position in all periods snum is counted, and the p value is num / sum. The corresponding 1-p is 1-num / sum. The weighted sum calculation formula is: NDVI st =(1-p)·SNVDVI+p·TNDVI Where NDVI st is the NDVI time series obtained by filtering; 42) Output NDVI st This is the filtered NDVI time series image.

Images