Remote sensing vegetation index time series data reconstruction method based on improved space-time SG
By improving the spatiotemporal SG method, using the autocorrelation of the VI time series and the spatiotemporal data of the current year, the problems of overcorrection and missing value filling in the reconstruction of the remote sensing vegetation index timing data are solved, and more accurate vegetation phenological characteristics simulation and stability are achieved performance.
Patent Information
- Application Number
- CN202411209955.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-08-30
- Publication Date
- 2025-01-03
AI Technical Summary
The prior art is prone to over-correcting the true local low VI value when reconstructing remote sensing vegetation index timing data, and performance in missing value filling depends on multi-year data consistency, especially when crop rotation or insufficient data is insufficient for reconstruction accuracy.
Using the improved spatiotemporal SG method, the autocorrelation of the VI time series and the spatiotemporal data of the current year are used to reconstruct the vegetation index time series through steps such as mass label reclassification, spatiotemporal fusion, and iterative SG filtering to avoid over-correction and be applicable to time. Filling of consecutive missing values.
It improves the simulation accuracy, stability and computational efficiency of vegetation phenological characteristics, and is suitable for VI time series with time-consecutive missing values. It only requires data from the current year, and its performance is better than traditional methods.
Smart Images
Figure CN119251329A_ABST
Abstract
Description
Technical Field
[0001] The invention relates to the technical field of remote sensing image processing, and in particular to a remote sensing vegetation index time series data reconstruction method for improving spatiotemporal SG. Background Art
[0002] Remote sensing vegetation index (VI) time series data is one of the most widely used remote sensing datasets. Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) time series data have been widely used in the study of various surface parameters, including leaf area index, vegetation phenology, and land cover. Although high-quality remote sensing NDVI / EVI time series data are crucial for these applications, these data are easily contaminated by adverse atmospheric conditions (such as cloud cover).
[0003] A variety of methods have been developed to remove residual noise from remote sensing VI time series products. Early methods were time-based, which used various temporal filters to smooth the original VI time series or adopted different mathematical functions to simulate the annual trajectory of the VI time series. In recent years, more complex methods integrating spatiotemporal information have been proposed to reconstruct remote sensing VI time series. These methods consider both the spatial correlation between pixels in adjacent windows and the temporal correlation of multiple years of VI data, thus achieving significantly better performance in reconstructing VI time series with continuous data missing. Typical spatiotemporal-based methods include the spatiotemporal Savitzky-Golay (STSG) method, the spatiotemporal tensor completion method, and the spatiotemporal pre-filling method with time series harmonic analysis.
[0004] In order to improve the quality of remote sensing vegetation index time series data, most previous methods have focused on correcting the negative bias VI noise by approximating the upper envelope of the VI time series. However, this method may cause over-correction of some real local low VI values in the time series, resulting in inaccurate simulation of vegetation phenological characteristics. In addition, another challenge in reconstructing remote sensing VI time series is to fill in the continuous missing values in time. The early spatiotemporal Savitzky-Golay (STSG) method solved this problem by utilizing multi-year VI data, but its performance is highly dependent on the consistency of data from different years. If some crop rotations are obvious or multi-year data cannot be provided, the reconstruction accuracy will not be guaranteed. Therefore, it is necessary to further improve the remote sensing vegetation index time series image reconstruction method. Summary of the invention
[0005] The present invention proposes a method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG. The method of the present invention takes into account the autocorrelation in the VI time series, and uses the spatiotemporal VI data of the current year to fill in the missing values in the VI time series to reconstruct the vegetation index time series image.
[0006] The technical solution adopted by the present invention is as follows:
[0007] A method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG, characterized in that the method uses a single-year VI (vegetation index) time series as raw data, and reconstructs a VI time series image by using temporal correlation and spatial correlation. Specifically, the method comprises the following steps:
[0008] S1. Reclassify the original data according to the quality label band, determine the quality index (Reliability Index, RI) and perform preprocessing to obtain the preprocessed VI time series.
[0009] S2. For the VI time series obtained in step S1, the spatial correlation of the VI remote sensing image and the temporal correlation of the VI time series are used to open spatial windows and temporal windows for spatiotemporal fusion, respectively, to generate an initial VI estimation value.
[0010] S3. Correct the preprocessed VI time series according to the quality index and the initial VI estimation value to obtain a corrected VI time series.
[0011] S4. Perform iterative SG filtering on the corrected VI time series and keep the phenological characteristics to obtain the reconstructed VI time series.
[0012] Furthermore, in step S1, the original data is reclassified according to the quality label band as follows: each pixel in the remote sensing image is marked, and the VI value with good quality is marked with RI (quality index) as good; the VI value with uncertain quality is marked with RI as uncertain; and the VI value with vacancy, cloud, ice and snow is marked with RI as contaminated.
[0013] The preprocessing method is as follows: retain the VI values marked as good and uncertain, set the VI values of the pixels marked as contaminated by RI as missing, and exclude the values with VI values less than a given threshold (such as NDVI < 0.15, which means that the pixel position is not vegetation and does not need to be reconstructed).
[0014] Furthermore, in step S2, define is the contaminated or uncertain VI value of the target pixel p0 on date Ti, and its initial value is recorded as The initial value is estimated based on the VI value in the time window centered on date Ti, using the following formula:
[0015]
[0016] in, represents the VI value at date Tj, and are the slope and intercept in linear regression; tw represents the length of the time half-window;
[0017] Use the VI values marked as good and uncertain within the time window to estimate and
[0018]
[0019] in, and Respectively represent adjacent pixels p in a space n The good or uncertain VI values at dates Ti and Tj in the VI time series. N represents the number of similar pixels (adjacent pixels with correlation coefficients greater than 0.85) within the spatial window centered on the target pixel p0, and the spatial half-window length is denoted as sw.
[0020] Furthermore, in step S3, the preprocessed VI time series is corrected according to the quality index and the initial VI estimate value as follows: the VI values marked as contaminated in the preprocessed VI time series are directly replaced by the corresponding initial VI estimate values; the VI values marked as uncertain in the preprocessed VI time series are compared with the corresponding initial VI estimate values, and the larger value is selected to be retained, because the contaminated VI (NDVI and EVI) values are usually low; the corrected VI time series is recorded as VI (corrected).
[0021] Furthermore, in step S4, the method of iterative SG filtering of the corrected VI time series and maintaining the phenological characteristics is as follows: applying the SG (Savitzky-Golay) filter to smooth VI (corrected), setting the window length and order of the polynomial in the SG filter to 7 and 4 respectively, and obtaining the VI time series SG[VI (corrected)]1 after the first filtering.
[0022] Compare each value of the corrected VI time series VI(corrected) with the filtered VI time series SG[VI(corrected)]1, retain the larger value between VI(corrected) and SG[VI(corrected)]1 to obtain the fused filtered fused VI time series; repeat the filtering and fusion process to obtain
[0023]
[0024] in, represents the corrected VI value of the target pixel p0 at time Ti, Represents the VI value of the target pixel p0 after the k-th filtering and fusion at time Ti. The recommended value of k is 3.
[0025] Define normalized to the interval [0,1] The expression is:
[0026]
[0027]
[0028] in, Represents the vegetation phenology feature weight of the target pixel p0 in the local low-value segment at time Ti. This weight is related to the vegetation phenology feature of the local low-value segment. The more obvious the feature, the higher the weight value. It represents the corrected VI value of the target pixel p0 in the time window centered at time Ti, and meanV represents the mean of the VI values in the time window excluding the center time Ti.
[0029] right and Perform weighted summation to generate the reconstructed
[0030]
[0031] in, Represents the reconstructed VI value of the target pixel p0 at date Ti.
[0032] Compared with the prior art, the present invention has the following beneficial effects: the method of the present invention better simulates the vegetation phenological characteristics of certain key stages in the reconstructed VI time series by automatically quantifying possible overcorrections in the VI time series and making corresponding adjustments, thereby avoiding overcorrection of some real local low VI values caused by excessive pursuit of the upper envelope. The present method fills the missing values in the VI time series by considering the autocorrelation of the VI time series, which makes the new method applicable to VI time series with temporally continuous missing values. In addition, the method of the present invention only requires the VI time series data of the current year, and compared with the STSG method that requires many years of historical VI data, it has more stable performance and can achieve a good balance between reconstruction accuracy and computational efficiency. BRIEF DESCRIPTION OF THE DRAWINGS
[0033] Figure 1 The flowchart of the remote sensing vegetation index time series data reconstruction method for improving the spatiotemporal SG of the present invention;
[0034] Figure 2 is the polluted or uncertain NDVI value of the target pixel p0 on date Ti in the embodiment of the present invention Its initial value the determination process;
[0035] Figure 3 Comparison of the remote sensing vegetation index time series data reconstruction method of the improved spatiotemporal SG proposed in the embodiment of the present invention with the Chen's SG and STSG methods in reconstructing the real NDVI time series results in 2020 in four study areas;
[0036] Figure 4 The average absolute difference between the NDVI values reconstructed by the embodiment of the present invention on the post-harvest date (161st day) of winter wheat pixels in the North China Plain and the reference NDVI values is compared with the Chen's SG and STSG methods;
[0037] Figure 5 The embodiment of the present invention reconstructs NDVI time series data in Northeast Asia to detect the vegetation greening period (GOD) and the vegetation withering period (DOD), and compares the difference histogram of directly detecting GOD and DOD with the reference curve and the histogram of Chen's SG and STSG methods. DETAILED DESCRIPTION
[0038] In order to make the purpose and technical solution of the embodiments of the present application clearer, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.
[0039] Taking the reconstruction of the 2020 NDVI time series in four 1000km*1000km study areas in Northeast Asia, North China Plain, Southeast China, and Southeast Asia as an example, in this embodiment, the remote sensing vegetation index time series data reconstruction method of the improved spatiotemporal SG is as follows: Figure 1 As shown, the following steps are included:
[0040] S1. Obtain MODIS NDVI data in the study area, reclassify it according to the quality label band, determine the quality index (Reliability Index, RI) and preprocess it to obtain the preprocessed NDVI time series.
[0041] S2. For the NDVI time series obtained in step S1, the spatial correlation of the NDVI remote sensing image and the temporal correlation of the NDVI time series are used to open spatial windows and temporal windows for spatiotemporal fusion to generate an initial NDVI estimate.
[0042] S3. Correct the preprocessed NDVI time series according to the quality index and the initial NDVI estimation value to obtain a corrected NDVI time series.
[0043] S4. Perform iterative SG filtering on the corrected NDVI time series and keep the phenological characteristics to obtain the reconstructed NDVI time series.
[0044] Specifically, in step S1, the MODIS NDVI time series data for 2020 in the study area is obtained. In the MODIS quality label band, NDVI values with good quality (i.e., no clouds and no snow) have a pixel label of 0, while uncertain NDVI values are labeled 1, and contaminated NDVI values are labeled -1 (no data), 2 (snow and ice), or 3 (clouds). The NDVI values with a quality label value of 0 are marked as good by RI; the NDVI values with a quality label of 1 are marked as uncertain by RI; and the NDVI values with quality labels of -1, 2, and 3 are marked as contaminated by RI. The NDVI values marked as good and uncertain are retained, the contaminated NDVI values are set to missing, and the NDVI values less than 0.15 are excluded to obtain the preprocessed NDVI time series.
[0045] In step S2, the time half-window length is 4, and the space half-window length is 5 km. In the case of 16-day time resolution, the time window is 64 days before and after the time point; the space window size is 11 km*11 km. Adjacent pixels in the space window are considered similar if the correlation coefficient with the target pixel is greater than 0.85. Figure 2 is the polluted or uncertain NDVI value of the target pixel p0 on date Ti Its initial value If there are not enough pixels available in the time window to generate an initial estimate The time window and space window can be appropriately expanded, or filled using linear interpolation.
[0046] is the polluted or uncertain NDVI value of the target pixel p0 on date Ti, and its initial value is recorded as The initial NDVI value is estimated based on the NDVI value in the time window centered on date Ti, and the formula is as follows:
[0047]
[0048] in, represents the NDVI value at date Tj, and are the slope and intercept in the linear regression; tw represents the length of the time half-window; the NDVI values marked as good and uncertain within the time window are used to estimate the equation and The expression is as follows:
[0049]
[0050] in, and Respectively represent adjacent pixels p in a space n good or uncertain NDVI values at dates Ti and Tj in the NDVI time series, and N represents the number of similar pixels within the spatial window centered on the target pixel p0.
[0051] In step S3, the preprocessed NDVI time series is corrected according to the quality index and the initial NDVI estimate. The NDVI values marked as contaminated in the preprocessed NDVI time series are directly replaced by the corresponding initial NDVI estimate. The NDVI values marked as uncertain in the preprocessed NDVI time series are compared with the corresponding initial NDVI estimate, and the larger value is selected to be retained because the contaminated NDVI values are usually low. The corrected NDVI time series is recorded as VI (corrected).
[0052] In step S4, the SG (Savitzky-Golay) filter is applied to smooth the corrected NDVI time series VI (corrected). The window length and order of the polynomial in the SG filter are set to 7 and 4, respectively, and the NDVI time series SG[VI (corrected)]1 after the first filtering is obtained.
[0053] Compare each value in the corrected NDVI time series VI(corrected) with the filtered NDVI time series SG[VI(corrected)]1, retain the larger value between VI(corrected) and SG[VI(corrected)]1 to obtain the fused filtered fused NDVI time series, repeat the filtering and fusion process, and obtain
[0054]
[0055] in, represents the corrected NDVI value of the target pixel p0 at time Ti, It represents the NDVI value of the target pixel p0 after the k-th filtering and fusion at time Ti. In this embodiment, k is 3.
[0056] Define normalized to the interval [0,1] The expression is:
[0057]
[0058]
[0059] in, Represents the vegetation phenology feature weight of the target pixel p0 in the local low-value segment at time Ti. This weight is related to the vegetation phenology feature of the local low-value segment. The more obvious the feature, the higher the weight value. It represents the corrected NDVI value of the target pixel p0 in the time window centered at time Ti, and meanV represents the mean of the NDVI values in the time window excluding the center time Ti.
[0060] After three iterations, and Perform weighted summation to generate the reconstructed
[0061]
[0062] in, Represents the reconstructed NDVI value of the target pixel p0 on date Ti.
[0063] like Figure 3 As shown in the figure, the remote sensing vegetation index time series data reconstruction method of the improved spatiotemporal SG of the present invention was tested in four study areas, and compared with Chen's SG and STSG methods, which are also based on SG filtering. The four test areas are representative of different climates and vegetation cover types. The first site is located in northeastern China and eastern Russia (referred to as Northeast Asia), mainly covered by deciduous forests and grasslands. The second site is the North China Plain, which is mainly dominated by double-season crops (winter wheat and corn) and grasslands. Winter wheat is harvested from May to June, followed by corn planting, resulting in two vegetation index peaks in a year. The third site is located in southeastern China, with a mixture of evergreen forests and deciduous forests. The fourth site is located in Southeast Asia, covering most of Thailand, where rice is the main crop. About 75% of rice is grown in the rainy season (May to October), while the remaining 25% is grown in both the dry and rainy seasons, which is double-season rice.
[0064] The following three experiments are provided to evaluate the effectiveness of the improved spatiotemporal SG remote sensing vegetation index time series data reconstruction method of the present invention.
[0065] Since the real VI time series data itself has noise and no true value, for quantitative evaluation, simulated VI time series data for 2020 and reference VI time series data as the true value were generated. The mean absolute difference (MAE) between the simulated VI time series and the reference VI time series data after smoothing by the method of the present invention is used for quantitative evaluation to quantify the overall fidelity. The method for generating simulated VI time series data and reference VI time series data is as follows (taking NDVI as an example): The reference NDVI time series for each pixel is generated by averaging all valid NDVI values on the same day from 2001 to 2020. The reference NDVI time series is regarded as the true NDVI time series of the pixel and is regarded as the true value. Then, according to the NDVI data quality label band and data conditions in 2020, the reference NDVI time series data is simulated for missing values and NDVI negative bias noise, and the reference NDVI time series data after noise addition is regarded as NDVI time series data with noise, which is called simulated NDVI time series data.
[0066] The design and implementation results of each experiment are as follows:
[0067] Experiment 1: Three methods including the method of the present invention were quantitatively evaluated based on simulated NDVI time series data and EVI time series (simulation method is consistent with NDVI) in four test areas. The test methods include the remote sensing vegetation index time series data reconstruction method of the present invention with improved spatiotemporal SG, Chen's SG and STSG methods.
[0068] Tables 1 and 2 show the mean absolute difference (MAE) values between the VI time series reconstructed by the three SG-based methods at the four test sites and the reference VI time series. The proposed method performed best on the five test sites with the smallest MAE value, followed by the STSG and Chen's SG methods. Specifically, the average MAE values of NDVI at the four test sites were 0.012 for the proposed method, 0.018 for the STSG method, and 0.020 for the Chen's SG method; the average MAE values of EVI were 0.011 for the proposed method, 0.021 for the STSG method, and 0.016 for the Chen's SG method. In addition, at the Southeast Asian site, due to the significant interannual land cover changes caused by frequent human activities in the region, the STSG method performed the worst using multi-year NDVI time series data, which explains its poor performance at the site. This experiment confirms that the proposed method is superior to the previous SG basic method in terms of the improvement in the fidelity of reconstructed VI values.
[0069] Table 1. Average MAE between the reconstructed NDVI time series and the reference NDVI time series
[0070]
[0071] Table 2. Average MAE between the reconstructed EVI time series and the reference EVI time series
[0072]
[0073] Experiment 2: Evaluate the maintenance of low NDVI values caused by the critical harvest period of winter wheat in the North China Plain.
[0074] Pixels with continuous winter wheat cultivation from 2008 to 2020 in the North China Plain of the test area were identified as evaluation pixels. The three SG-based methods were applied to the simulated NDVI time series data, and the MAE between the reconstructed NDVI value and the reference NDVI value at the approximate wheat post-harvest date (day 161) for each winter wheat pixel was calculated, as shown in Figure 4 As shown in the figure, the MAE value (0.020) of the proposed method on this date is smaller than that of the other two methods (STSG is 0.045, Chen's SG is 0.132), which indicates that the proposed method performs better in simulating the phenological characteristics of winter wheat harvest and can better simulate the vegetation characteristics of the key harvest stage of the reconstructed NDVI time series.
[0075] Experiment 3: The difference between the real phenological date and the estimated phenological date based on the simulated NDVI time series data in the northeastern East Asia in 2020.
[0076] Three SG-based methods including the method of the present invention were used to smooth the simulated NDVI time series data, and the vegetation greening period (GOD) and vegetation withering period (DOD) were detected, and the differences between the vegetation greening period (GOD) and vegetation withering period (DOD) were directly detected with the reference curve. Figure 5 As shown, the GOD values generated by Chen's SG method are significantly advanced (i.e., the GOD-d values are positive), and the DOD values are significantly delayed (i.e., the DOD-d values are negative), highlighting its limitations in vegetation phenology applications. In contrast, the method of the present invention effectively solves this systematic error and achieves the smallest GOD-d (1.1 day) and DOD-d (0.3 day) values among the three methods. In previous methods, the principle of approaching the upper envelope of the NDVI time series may lead to overcorrection of some true local low NDVI values. The method of the present invention solves this problem by automatically quantifying the possible overcorrection in the NDVI time series and making corresponding adjustments. Compared with previous methods, the method of the present invention can more accurately extract the growth start date (GOD) and dormancy start date (DOD) in the NDVI time series.
[0077] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.
Claims
1. A method for reconstructing remote sensing vegetation index time series data based on improved spatiotemporal SG, characterized in that: The method uses a single-year VI time series as raw data and reconstructs a VI time series image using temporal correlation and spatial correlation. Specifically, the method includes the following steps: S1, reclassify the original data according to the quality label band, determine the quality index and preprocess it to obtain the preprocessed VI time series; S2, for the VI time series obtained in step S1, using the spatial correlation of the VI remote sensing image and the temporal correlation of the VI time series, respectively open a spatial window and a temporal window for spatiotemporal fusion to generate an initial VI estimation value; S3, correcting the preprocessed VI time series according to the quality index and the initial VI estimation value to obtain a corrected VI time series; S4. Perform iterative SG filtering on the corrected VI time series and keep the phenological characteristics to obtain the reconstructed VI time series.
2. The method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG according to claim 1, characterized in that: In step S1, the original data is reclassified according to the quality label bands as follows: Each pixel in the remote sensing image is marked, and the VI value with good quality is marked as good; The VI value with uncertain mass has its RI marked as uncertain; The VI values with quality of vacancy, cloud presence, and ice and snow presence have their RI marked as contaminated; The preprocessing method is as follows: retain the VI values marked as good and uncertain by RI, set the VI values marked as contaminated pixels by RI as missing, and exclude the values with VI values less than a given threshold.
3. The method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG according to claim 2, characterized in that: In step S2, define is the contaminated or uncertain VI value of the target pixel p0 on date Ti, and its initial value is recorded as The initial value is estimated based on the VI value in the time window centered on date Ti, using the following formula: in, represents the VI value at date Tj, a Tj and b Tj are the slope and intercept in linear regression; tw represents the length of the time half-window; Use the VI values marked as good and uncertain within the time window to estimate a Tj and b Tj : in, and Respectively represent adjacent pixels p in a space n The good or uncertain VI values at dates Ti and Tj in the VI time series; N represents the number of similar pixels in the spatial window centered on the target pixel p0, and the spatial half-window length is denoted as sw.
4. The method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG according to claim 3, characterized in that: In step S3, the preprocessed VI time series is corrected according to the quality index and the initial VI estimate value as follows: The VI values marked as contaminated by RI in the preprocessed VI time series are directly replaced by the corresponding initial VI estimates; the VI values marked as uncertain by RI in the preprocessed VI time series are compared with the corresponding initial VI estimates, and the larger value is selected to be retained; The corrected VI time series is recorded as VI(corrected).
5. The method for reconstructing remote sensing vegetation index time series data by improving spatiotemporal SG according to claim 4, characterized in that: In step S4, the method of iterative SG filtering of the corrected VI time series and maintaining the phenological characteristics is as follows: applying the SG filter to smooth VI(corrected), setting the window length and order of the polynomial in the SG filter to 7 and 4 respectively, and obtaining the VI time series SG[VI(corrected)]1 after the first filtering; Compare each value of the corrected VI time series VI(corrected) with the filtered VI time series SG[VI(corrected)]1, and retain the larger value between VI(corrected) and SG[VI(corrected)]1 to obtain the fused filtered fused VI time series; Repeat the filtering and fusion process to obtain in, represents the corrected VI value of the target pixel p0 at time Ti, Represents the VI value of the target pixel p0 after the kth filtering and fusion at time Ti. Define normalized to the interval [0,1] The expression is: in, Represents the vegetation phenology feature weight of the target pixel p0 in the local low-value segment at time Ti. This weight is related to the vegetation phenology feature of the local low-value segment. The more obvious the feature, the higher the weight value. represents the corrected VI value of the target pixel p0 in the time window centered at time Ti, and meanV represents the mean of the VI values in the time window excluding the center time Ti; right and Perform weighted summation to generate the reconstructed in, Represents the reconstructed VI value of the target pixel p0 at date Ti.