High-resolution normalized difference vegetation index time-series reconstruction method based on adaptive uncertainty

Abstract

This invention belongs to the field of remote sensing science and technology, and relates to a high-resolution normalized difference vegetation index (NDVI) time-series reconstruction method based on adaptive uncertainty. The method includes collecting data to obtain a priori NDVI image time series and a coarse-resolution NDVI time series; establishing a dynamic state-space model based on the priori NDVI image time series; estimating the high-resolution observation uncertainty, the dynamic state-space model uncertainty, and the coarse-resolution observation uncertainty; constructing an ensemble Kalman smoother for multi-scale data assimilation to obtain the reconstructed high-resolution NDVI time series. This invention solves the technical problems in existing technologies, such as the inability to capture rapid changes due to data scarcity, the failure of overly simplified assumptions in mixed cropping areas, and the lack of flexibility for specific sensors.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing science and technology, specifically relating to a high-resolution normalized differential vegetation index time-series reconstruction method based on adaptive uncertainty. Background Technology

[0002] Monitoring terrestrial ecosystems requires consistent and fine-scale observations of vegetation dynamics. Time series of vegetation indices provide crucial information for ecosystem processes such as productivity estimation, phenological tracking, and land cover change detection. However, a persistent trade-off between spatial and temporal resolution, particularly in historical satellite archives, limits the ability to monitor vegetation at long-term, fine-resolution scales. Atmospheric disturbances, cloud pollution, and sensor limitations further exacerbate this challenge, collectively hindering applications that require frequent and fine-resolution observations.

[0003] To bridge data gaps, multi-sensor integration has become a practical solution. In recent years, satellites like Sentinel-2, equipped with dual-satellite constellations, can theoretically achieve 5-day revisit observations at a resolution of 10-20 meters. Combined with Landsat, the Landsat-Sentinel-2 (HLS) dataset now provides global land observations every 2-3 days at a resolution of 30 meters. However, the practical usability of Sentinel-2 data remains limited by persistent cloud cover, limited swath width, and a short historical coverage area, highlighting the need for alternative methods to generate dense, high-quality VI (Vegetation Index) time series. Spatiotemporal fusion methods address this issue by integrating coarse-resolution but frequent observations with sparse high-resolution imagery to generate dense, high-resolution time series. Over the past two decades, various spatiotemporal fusion techniques have been developed, including unmixing-based, weighted function-based, Bayesian-based, learning-based, and hybrid methods. These methods typically rely on coarse-resolution and high-resolution image pairs (base pairs) to predict high-resolution data for a target date, enabling detailed vegetation monitoring. Among these methods, unmixing-based approaches were the earliest, while weight function-based methods, such as the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM) and its enhanced version (ESTARFM), are widely used. Weight function-based methods utilize temporal interpolation and spatial similarity constraints to predict high-resolution observations for target dates; recent improvements, such as STARFM based on similarity weight learning (SWL-STARFM), optimize the similarity metric. Hybrid models, such as Flexible Spatiotemporal Data Fusion (FSDAF) and its extended version SFSDAF, combine unmixing and weighting strategies to better recover sub-pixel information. Furthermore, dynamic assimilation methods based on Kalman filtering have been introduced, utilizing linear mixture models or prior high-resolution information to integrate coarse-resolution and high-resolution data.

[0004] Spatiotemporal fusion can be performed at the reflectance level or directly on vegetation indices. For vegetation monitoring, two strategies are commonly used: (1) fusing reflectance first and then calculating indices (fusion before indexing), or (2) directly fusing indices (indexing before fusion). The latter is generally preferred in phenological studies because it reduces error propagation, lowers computational costs, and is less sensitive to radiation inconsistencies and algorithmic assumptions.

[0005] Despite these methodological advances, existing spatiotemporal fusion methods still face fundamental limitations in heterogeneous agricultural regions, particularly under data-scarce conditions. In practice, data scarcity manifests not only in the limited number of high-resolution images but also in the inability to capture rapid surface changes. Due to the combined effects of cloud contamination, revisit cycles, and sensor swath width limitations, available high-resolution observations are often temporally distant from the prediction date. This problem is particularly acute during periods of rapid vegetation growth or aging, when critical phenological transitions occur within short time windows that are difficult to observe with sparse imagery. Most spatiotemporal fusion methods were initially designed for landscapes dominated by single-season vegetation, and their foundation rests on strong assumptions about the dynamic spatiotemporal behavior of the surface, such as the temporal stationarity of spatial patterns (i.e., time-stable coarse-to-fine resolution relationships) or the spatial invariance of temporal models (i.e., uniform temporal variation across spatial scales). These assumptions break down under sparse high-resolution observations, and some methods, including ESTARFM (Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model) and LMGM (Local Multi-scale Gradient Maximum), further assume linear temporal variations between observation dates.

[0006] While these assumptions may be acceptable in relatively homogeneous landscapes with consistent vegetation composition evolution, they often fail in heterogeneous agricultural areas where single- and double-cropping systems coexist (e.g., winter wheat-summer maize rotation). In such mixed-cropping systems, vegetation dynamics become phenologically asynchronous at the sub-pixel scale, with different crop components potentially exhibiting different or even opposite temporal variations during key transition periods. Consequently, coarse-resolution vegetation indices like MODIS often mask multi-seasonal signals, and the assumed temporal similarity between the base image pair and the predicted date becomes highly unstable. Previous studies have shown that spatiotemporal fusion accuracy deteriorates rapidly with increasing temporal distance between base image pairs precisely because this assumed temporal similarity no longer holds. In heterogeneous agricultural landscapes, this similarity fluctuates nonlinearly with seasonal transitions, making the reconstruction of multi-seasonal phenological curves particularly difficult—especially when high-resolution observations are sparse. Furthermore, many existing spatiotemporal fusion algorithms are tightly coupled to specific sensor pairs, limiting their flexibility in incorporating additional high-resolution observations as the availability of high-resolution observation data varies with sensor and time. Although many approaches were developed around the MODIS-Landsat combination, they are often limited to a single high-resolution data source when extended to newer sensors such as Sentinel-2, which limits their portability and scalability in multi-sensor environments.

[0007] Meanwhile, the performance of spatiotemporal fusion methods in reconstructing phenology in mixed seasonal landscapes, especially in system evaluation under data-scarce conditions, remains underdeveloped. Therefore, how to solve the technical problems in existing technologies, such as the inability to capture rapid changes due to data scarcity, the failure of overly simplified assumptions in mixed-cropping areas, and the lack of flexibility in sensor binding, has become a key bottleneck that urgently needs to be overcome. Summary of the Invention

[0008] The purpose of this invention is to provide a high-resolution normalized differential vegetation index time-series reconstruction method based on adaptive uncertainty. This method solves the problem that existing spatiotemporal fusion methods are generally based on strong assumptions such as single-season vegetation, linear changes, and specific sensor coupling, which makes it difficult to effectively reconstruct the complex, asynchronous, nonlinear phenological dynamics in single- and double-season mixed crop areas under sparse high-resolution observation conditions, thus limiting their accuracy and applicability in actual heterogeneous agricultural landscape monitoring.

[0009] To solve the above-mentioned technical problems, the technical solution adopted by this invention is a high-resolution normalized differential vegetation index time-series reconstruction method based on adaptive uncertainty, comprising the following steps: S1: Collect historical high-resolution normalized difference vegetation index (NDVI) images and coarse-resolution NDVI time series, preprocess them respectively to obtain a priori NDVI time series with fixed time intervals and a coarse-resolution NDVI time series. S2: Establish a dynamic state-space model based on the time series of prior normalized vegetation index with fixed time intervals; S3: Estimate the uncertainty of high-resolution observations, the uncertainty of dynamic state-space models, and the uncertainty of coarse-resolution observations, respectively; S4: Construct an ensemble Kalman smoother to perform multi-scale data assimilation and obtain the reconstructed high-resolution normalized differential vegetation index time series.

[0010] The beneficial effects of this invention are: 1) The FSATS (Flexible and Self-Adaptive Trend Surface) algorithm proposed in this invention integrates adaptive seasonal decomposition derived from long-term historical phenological trends with a unified ensemble Kalman Smoother (EnKS), thereby achieving uncertainty-aware assimilation of multi-source NDVI (Normalized Difference Vegetation Index) time series. Extensive experiments across multiple fusion scenarios and benchmark algorithms in this invention demonstrate the robustness and adaptability of FSATS. When jointly assimilating Landsat and MODIS observations, FSATS achieves the highest reconstruction accuracy (RMSE = 0.092, R = 0.838), outperforming three representative STF methods, namely ESTARFM (Enhanced Spatialand Temporal Adaptive Reflectance Fusion Model), LMGM (Linear Mixing Growth Model), and PK-STEM (Prior Knowledge-based Spatio-Temporal Enhancement Method), in terms of spatial detail preservation and temporal continuity.

[0011] 2) In terms of methodology, the FSATS algorithm proposed in this invention is specifically designed to address two persistent challenges in spatiotemporal fusion: sparse high-resolution observations and phenological complexity in heterogeneous agricultural systems. By adaptively characterizing single-season and double-season vegetation using long-term phenological signals and embedding this information into an uncertainty-aware fusion framework, the FSATS algorithm remains robust even when the temporal similarity between observations fails under mixed cropping conditions. This enables reliable fine-scale phenological reconstruction in scenarios where the performance of traditional STF methods degrades due to data scarcity and asynchronous vegetation dynamics. Attached Figure Description

[0012] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0013] Figure 1 This is a flowchart of the FSATS (Flexible and Self-Adaptive Trend Surface) algorithm of this invention; Figure 2 This is a comparison of HANTS (Harmonic Analysis of Time Series) and DL (Dual Logic Method) filtering on single and double-season MODIS NDVI time series. Among them, (a) is a comparison of the smoothing results of MODIS NDVI time series of HANTS and DL filtering on double-season vegetation 1, (b) is a comparison of the smoothing results of MODIS NDVI time series of HANTS and DL filtering on double-season vegetation 2, (c) is a comparison of the smoothing results of MODIS NDVI time series of HANTS and DL filtering on single-season vegetation 1, and (d) is a comparison of the smoothing results of MODIS NDVI (Moderate Resolution Imaging Spectroradiometer Normalized Difference Vegetation Index) time series of HANTS and DL filtering on single-season vegetation 2. Figure 3 It is a map of the study area; Figure 4 The time distribution of Landsat NDVI imagery from 1995 to 2007 based on day sequence is shown. Figure 5 The spatial distribution of the uncertainty estimates of the dynamic model and MODIS NDVI observations for DOY (Day of Year) 001, 065, 129, and 289 is shown in (a) for DOY=001; (b) for DOY=065; (c) for DOY=129; and (d) for DOY=289. Figure 6This study compares the spatial and quantitative predictions of FSATS across dates under three operating modes within the research area. (a) shows the actual Landsat NDVI; (b) shows the prediction results for the DOY 289 sub-region under Mode I; (c) shows the prediction results for the DOY 289 sub-region under Mode II; (d) shows the prediction results for the DOY 289 sub-region under Mode III; (e) shows the correlation coefficient histograms for the three modes at DOY 65, 129, and 289 and on the three dates; (f) shows the root mean square error histograms for the three modes at DOY 65, 129, and 289 and on the three dates; (g) shows the mean absolute error histograms for the three modes at DOY 65, 129, and 289 and on the three dates; (h) shows the scatter plot of predicted and actual values ​​on the validation date DOY 065; (i) shows the scatter plot of predicted and actual values ​​on the validation date DOY 129; and (j) shows the prediction results for the three validation dates (DOY 065, ... Scatter plot of predicted and actual values ​​(065, 129, 289); Figure 7 This is a comparison of the NDVI time series predicted by FSATS for different sample points under modes I–III, where (a) is cultivated land, (b) is forest land, and (c) is grassland. Figure 8 The comparison shows the NDVI prediction results of DOY 289, where (a) is input image A, (b) is input image B, (c) is input image C, (d) is input image D, (e) is input image E, (f) is the true Landsat NDVI, (g) is the prediction result obtained by ESTARFM, (h) is the prediction result obtained by LMGM, (i) is the prediction result obtained by PK-STEM, and (j) is the prediction result obtained by FSATS. Figure 9These are scatter plots of actual and predicted NDVI for the ESTARFM, LMGM, PK-STEM, and FSATS methods in DOY 065, 129, and 289. (a) shows the actual and predicted NDVI for the ESTARFM method in DOY 065, 129, and 289; (b) shows the actual and predicted NDVI for the LMGM method in DOY 065, 129, and 289; (c) shows the actual and predicted NDVI for the PK-STEM method in DOY 065, 129, and 289; (d) shows the actual and predicted NDVI for the FSATS method in DOY 065, 129, and 289; (e) shows the actual and predicted NDVI for the ESTARFM (Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model) method in DOY 065; and (f) shows the actual and predicted NDVI for the LMGM (LinearMixing Growth Model) method in DOY 065. Scatter plots of actual and predicted NDVI in DOY 065, (g) scatter plot of actual and predicted NDVI in DOY 065 using the PK-STEM (Prior Knowledge-based Spatio-Temporal Enhancement Method), (h) scatter plot of actual and predicted NDVI in DOY 065 using the FSATS method, (i) scatter plot of actual and predicted NDVI in DOY 129 using the ESTARFM method, (j) scatter plot of actual and predicted NDVI in DOY 129 using the LMGM method, (k) scatter plot of actual and predicted NDVI in DOY 129 using the PK-STEM method, (l) scatter plot of actual and predicted NDVI in DOY 129 using the FSATS method, (m) scatter plot of actual and predicted NDVI in DOY 289 using the ESTARFM method, (n) scatter plot of actual and predicted NDVI in DOY 289 using the LMGM method, and (o) scatter plot of actual and predicted NDVI in DOY 065 using the PK-STEM method. Scatter plot of actual and predicted NDVI for DOY 289, (p) is the scatter plot of actual and predicted NDVI for DOY 289 using the FSATS method; Figure 10 This is a comparison of the fused NDVI prediction time series (mean ± standard deviation) generated by four fusion methods; Figure 11 The image shows a comparison of pixel-level NDVI curves generated by four fusion methods, where (a) represents cultivated land, (b) represents woodland, and (c) represents grassland. Detailed Implementation

[0014] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0015] like Figure 1 As shown, the high-resolution normalized difference vegetation index time-series reconstruction method based on adaptive uncertainty proposed in this invention is carried out according to the following steps: S1: Collect historical high-resolution normalized difference vegetation index (NDVI) images and coarse-resolution NDVI time series, preprocess them respectively to obtain high-resolution prior normalized difference vegetation index time series and coarse-resolution prior normalized difference vegetation index time series at fixed time intervals. In a specific embodiment of the present invention, step S1 is as follows: S1.1: Collect historical high-resolution normalized vegetation index images available for many years, and compile the same DOY images for each year according to day sequence (DOY) to obtain a dense but irregular historical time series; In a specific embodiment of the present invention, this DOY-based synthesis method provides empirical priors for phenological dynamics. However, due to irregular collection intervals and interannual variability, the generated time series exhibits temporal discontinuity and abrupt fluctuations. To alleviate these problems, the present invention applies temporal filtering to smooth the synthesized NDVI trajectory curve.

[0016] S1.2: Harmonic analysis was performed on the historical time series using the HANTS algorithm to obtain the prior normalized differential vegetation index time series. As shown in equation (1): (1) in: For the prior normalized differential vegetation index time series. For the first The time points of each observation sample Here are the indices of the observed time points in the time series, j=1,2,...,N. The coefficient at zero frequency, The number of periodic terms in the series. For harmonic component indices, n=1, 2, ..., N f N is the number of observation samples in the time observation series. Harmonic frequency The corresponding coefficients of the first triangular component, For harmonic frequencies, Harmonic frequency The corresponding coefficients of the second triangular component, For frequency, parameter , , It belongs to real numbers and has no fixed range of values. From the perspective of the normalized difference vegetation index constraint, the parameters... Scope ,and and Each amplitude does not have its own fixed range, but the total amplitude must satisfy the following formula: (2) A significant feature of the FSATS algorithm proposed in this invention is its ability to simulate changing vegetation seasonality using a fixed-order HANTS filter (e.g., three harmonics). The advantage lies in the fact that this order is sufficient to capture single-season and double-season phenological patterns in different landscapes. This harmonic order is chosen after a trade-off between model flexibility and robustness, and previous studies have shown that it is sufficient to represent annual and semi-annual vegetation cycles without overfitting. This approach avoids the need for predefined land cover maps or crop calendars while maintaining robustness in heterogeneous agricultural areas.

[0017] like Figure 2 (a)-(d) show the comparison of smoothing results of MODIS NDVI (Medium Resolution Imaging Spectroradiometer Normalized Differential Vegetation Index) time series on double-season vegetation 1, double-season vegetation 2, single-season vegetation 1, and single-season vegetation 2, respectively. Compared with classical filters such as the dual logic (DL) method (A. Fischer, "A model for the seasonal variations of vegetation indices in coarse resolution data and its inversion to extract crop parameters," Remote Sensing of Environment, vol. 48, no. 2, pp. 220-230, 1994.), the HANTS algorithm used in this invention shows better adaptability to both single-season and double-season NDVI profiles, recovers peak vegetation dynamics, and maintains a smooth seasonal envelope. The filtered historical high-resolution NDVI dataset (spatial resolution of 30 meters) is called the prior NDVI time series.

[0018] S1.3: Temporally resample the prior normalized difference vegetation index time series at fixed time intervals (e.g., 8 days) to generate a prior normalized difference vegetation index time series with fixed time intervals (8 days). The purpose is to align with the time step of MODIS (Moderate Resolution Imaging Spectroradiometer) data (8 days) to facilitate the calculation of subsequent downscaling coarse resolution observation errors.

[0019] S1.4: Preprocess the coarse-resolution normalized differential vegetation index time series to obtain the coarse-resolution prior normalized differential vegetation index time series.

[0020] S2: Establish a dynamic state-space model based on the time series of prior normalized vegetation index images; The core of the dynamic state-space model is to use the data of two adjacent moments in the time series of the normalized differential vegetation index with a fixed time interval obtained in S1.3, and calculate the rate of change (i.e. the transition matrix) between adjacent moments through equations (3) and (4). Based on this rate of change, the normalized vegetation index at any moment can be obtained by forward reasoning from the initial state, or the state at past moments can be inferred backward from the new observation data. The FSATS algorithm captures the temporal evolution of the normalized vegetation index through the state-space model of physical information.

[0021] S2.1: Define the state equations for the dynamic state-space model, expressed as: (3) in: For the high-resolution normalized vegetation index field at time t, For the high-resolution normalized vegetation index field at time t-1, Here is the transition matrix; S2.2: Prior Normalized Differential Vegetation Index Value at Time t in the Prior Normalized Differential Vegetation Index Time Series The prior normalized difference vegetation index value at time t-1 Calculate the transition matrix As shown in equation (4), the dynamic state space model is obtained; (4) in: Let be the prior NDVI value at time t. Let be the prior NDVI value at time t-1. It is the first constant. In this invention, 0.001 is used to prevent division by zero. Let be the transition matrix.

[0022] This invention combines the prior normalized differential vegetation index value at time t. The prior normalized differential vegetation index value at time t-1 Calculate the transition matrix and with the transition matrix The values ​​of each row and column encode the relative growth rate and decay rate as relative changes. Compared with traditional additive models or linear interpolation, the method proposed in this invention is closer to the vegetation biophysical process. The dynamic model is particularly suitable for phenological modeling because the relative changes in vegetation conditions are more stable than absolute reflectance values ​​on sensor and spatial scales. This allows the FSATS algorithm proposed in this invention to realistically simulate the normalized differential vegetation index changes of asymptotic and abrupt changes during vegetation transition.

[0023] S3: Estimate the uncertainty of high-resolution observations, the uncertainty of dynamic models, and the uncertainty of coarse-resolution observations, respectively; A key innovation of the FSATS proposed in this invention lies in its dynamic and data-driven multi-scale uncertainty quantification, which enables adaptive and context-aware observation weighting during the assimilation process. Unlike traditional methods that assume fixed or empirical errors, FSATS explicitly simulates uncertainty from three different sources: (i) high-resolution observations, (ii) dynamic models, and (iii) coarse-resolution observations. Each component is based on the normalized vegetation index trend surface formula and estimated using forward simulation and back projection strategies. The dynamic model uncertainty characterizes the prediction error of the state-space model, while the coarse-resolution observation uncertainty characterizes the error introduced by the coarse-resolution observation data after spatial downscaling. In this invention, high-resolution observations are considered the most reliable observation source and assigned a constant uncertainty. =0.05.

[0024] The specific steps for S3 are as follows: S3.1: Addressing the uncertainty of high-resolution observations Let it be a constant =0.05, used to account for sensor radiated noise, atmospheric correction error and geometric misalignment. This parameter also serves as the minimum threshold for other uncertainty components to prevent overconfidence assimilation. High-resolution observations refer to remote sensing observation data with high spatial resolution that can provide more detailed surface spatial information.

[0025] S3.2: Estimating the uncertainty of the dynamic model This includes the following sub-steps: S3.2.1: Collect observations for each high-resolution normalized differential vegetation index ; To assess the inherent uncertainty of the normalized differential vegetation index (FSATS) dynamic model, the FSATS algorithm introduces a novel back projection strategy, the specific implementation of which is shown in S3.2.2: S3.2.2: Using the transition matrix Each high-resolution normalized differential vegetation index observation By projecting backwards to the initial time, we obtain the implicit initial state. , is represented as: (5) in: For high-resolution normalized differential vegetation index observation The implicit initial state For high-resolution normalized differential vegetation index observation, Let be the inverse matrix of the normalized differential vegetation index at time t.

[0026] For each high-resolution normalized difference vegetation index observation, operation S3.2.2 was performed to obtain a set of implicit initial state values. ,in, Let be the set of implicit initial state values ​​corresponding to all high-resolution observation times t. This refers to the data collection time for high-resolution normalized differential vegetation index observations.

[0027] S3.2.3: Calculate each high-resolution normalized differential vegetation index observation The corresponding implicit initial state The standard deviation of is used to obtain the uncertainty of the dynamic state-space model, which is expressed as: (6) in: For the uncertainty of dynamic models, For high-resolution normalized differential vegetation index observation Quantity, For high-resolution normalized differential vegetation index observation The collection time, The mean of the normalized differential vegetation index values ​​at the initial time (t=1) is reconstructed for all high-resolution normalized differential vegetation index observations. This refers to a specific observation time within the data collection period for high-resolution normalized difference vegetation index (NDVI) observations.

[0028] The mean of the initial time-normalized differential vegetation index values ​​reconstructed from all high-resolution normalized differential vegetation index observations is shown in Equation (7): (7) In S3.2.2, the projection reconstruction estimates the initial normalized difference vegetation index (NDVI) value at time t=1 if the modeling transformation perfectly captures the temporal dynamics of the high-resolution NDVI. Due to imperfect dynamics, the NDVI value is estimated from multiple high-resolution NDVI observations. Reconstructed They will deviate from each other, and this variability is used to calculate the model uncertainty of the initial state.

[0029] After completing the initial high-resolution normalized differential vegetation index prior field construction and uncertainty estimation, the uncertainty of the dynamic model... Forward propagation is performed using a dynamic process model; specifically, it utilizes the transition matrix. During the state prediction process, the mean of each set member is zero and the variance is calculated by superimposing a sum of zero values. The random noise is used to characterize the uncertainty brought about by the model structure error and the unmodeled process. Equation (6) produces the spatiotemporal adaptive model uncertainty that reflects the transformation noise and model shortcomings.

[0030] S3.3: Estimating the uncertainty of coarse-resolution observations This includes the following sub-steps: Coarse-resolution observations refer to remote sensing observation data with low spatial resolution but high temporal observation frequency. For example, the coarse-resolution normalized difference vegetation index provides high temporal resolution, but its spatial mismatch and scale difference introduce significant uncertainty when resampled to 30 meters. To address this, the FSATS algorithm proposed in this invention estimates this error by projecting high-resolution observations onto coarse-resolution time series.

[0031] S3.3.1: Construct a coarse-resolution time transition matrix using the time series of the coarse-resolution normalized differential vegetation index. , is represented as: (8) in: This is a coarse-resolution time transition matrix. The value of the coarse-resolution normalized vegetation index at time t is [value]. The value of the coarse-resolution normalized vegetation index at time t-1 is... It is the second constant. In this invention, 0.001 is used to avoid the denominator being 0.

[0032] S3.3.2: Utilizing a coarse-resolution transition matrix Each high-resolution normalized differential vegetation index observation Back-projecting to the initial time, we obtain the normalized differential vegetation index observations derived from the coarse-resolution time series. Implicit initial state , represented as: (9) in: For high-resolution normalized vegetation index observation The initial state, Let be the inverse of the coarse-resolution time transition matrix at time t. It is the inverse of the coarse-resolution time transition matrix at time t-1.

[0033] Equation (9) approximates the initial normalized vegetation index state implied by the normalized vegetation index dynamics derived from coarse-resolution time series, and all high-resolution observations Corresponding initial state The variability is used to define the uncertainty of coarse-resolution observations.

[0034] S3.3.3: Calculate the rate-normalized differential vegetation index derived from each coarse-resolution time series observation. Implicit initial state The standard deviation is used to obtain the coarse-resolution observation uncertainty. , represented as: (10) in: The mean of the coarse-resolution normalized vegetation index values ​​at the initial time (t=1) for dynamic reconstruction of all high-resolution normalized vegetation index observations; All high-resolution normalized vegetation index observations The mean value of the coarse-resolution normalized vegetation index for the initial date of dynamic reconstruction is shown in Equation (11): (11) Coarse resolution observation uncertainty σ c Using the coarse-resolution time transition matrix G t Forward propagation generates a temporally dense but spatially adaptive error field. Equation (11) takes into account potential sources of coarse resolution error, such as mixed pixels, cloud contamination, and angle effects.

[0035] S3.4: Addressing the uncertainties of the dynamic model respectively and coarse resolution observation uncertainty Apply a lower bound constraint to ensure that both are not less than the uncertainty of high-resolution observations. .

[0036] To ensure numerical stability and avoid overconfident assimilation in the low variance region, this invention addresses the uncertainty of the dynamic model obtained through dynamic estimation. Coarse-resolution observation uncertainty Apply a lower bound constraint condition as follows .

[0037] S4: Construct an ensemble Kalman smoother to perform multi-scale data assimilation and obtain the reconstructed high-resolution normalized differential vegetation index time series.

[0038] In a specific embodiment of the present invention, a high-resolution normalized vegetation index set is used as the system state variable, and MODIS low-resolution products are introduced as observation constraints. The state set is updated by forward filtering and backward smoothing within the time window. By fusing model prediction and multi-source observation information, the joint correction of the high-resolution remote sensing state field is achieved, and finally the posterior high-resolution NDVI time series under full-time observation conditions is obtained.

[0039] Since NDVI is the only state variable in this model, the system state is determined by the set matrix. Let N be the number of members in the set matrix, and let N be the mean of the set matrix. Set perturbation matrix Defined as .

[0040] S4.1: Use the dynamic state-space model to derive the normalized differential vegetation index value corresponding to the initial time step, and add an unbiased Gaussian distribution perturbation to this initial value of the normalized differential vegetation index value. Generate a set and the set As the initial time analysis set ; In a specific embodiment of the present invention, at the initial time step, the forward and backward rates of change of adjacent times are first calculated using prior normalized vegetation index time series data through equations (3) and (4), and then the initial normalized vegetation index value is obtained by combining equation (5). Then, the standard deviation is... Unbiased Gaussian distribution perturbation Add to the initial normalized vegetation index value to generate a set. .

[0041] S4.2: For each subsequent time step, utilize the dynamic transition matrix. The analysis set from the previous moment Forward propagation generates the prediction set for the current time step. ; In a specific embodiment of the present invention, the prediction set The former refers to the result obtained by extrapolating from time t-1 to time t using only the transition matrix, while the latter is the set obtained by correcting the prediction set using the Kalman gain matrix when there is observation data at time t. The prediction set provides the forward model prediction of the normalized vegetation index before absorbing new observations, and it is not the same as the forward analysis set.

[0042] S4.3: Determine whether there are any new observations at the current time. If it is a high-resolution observation or a coarse-resolution observation, calculate the Kalman gain using the corresponding observation uncertainty and update the prediction set to obtain the current analysis set. If there are no new observations, recursively obtain the normalized differential vegetation index value at the current time based on the predicted value at the previous time and the dynamic transition matrix. Integrate the results of each time to obtain the high spatial resolution normalized differential vegetation index time series. When performing high-resolution observations, the specific steps to obtain a high spatial resolution normalized differential vegetation index time series are as follows: S4.3.1: Selecting High-Resolution Observation Uncertainties And calculate its square. As observation variance; S4.3.2: Adding variance to high-resolution observations Gaussian perturbation, generating a set of perturbation observations. ; S4.3.3: Using measurement operators The prediction set at the current moment Mapping to a high-resolution observation space yields the predicted observation set. ; S4.3.4: Calculate the prediction error covariance matrix And construct the diagonal observation error covariance matrix. ; S4.3.5: Calculate the Kalman gain matrix; S4.3.6: Update the prediction set at the current time step according to the Kalman filter formula. Updated to forward analytic set ; Based on prediction set Calculate the prediction error covariance matrix As shown in equation (12), combined with the observation error covariance matrix Calculate the Kalman gain matrix with the measurement operator As shown in Equation (13), uncertainty perturbation is added to the observation data to generate a perturbation observation set.

[0043] (12) in: To predict the error covariance, The total number of members in the set. Let be the perturbation deviation matrix of the predicted set at time t. This is the transpose of the perturbation deviation matrix of the prediction set at time t.

[0044] (13) in: Here is the Kalman gain matrix. For the transpose of the prediction error covariance, A measurement operator that maps model states to the observation space. To map the model state to the transpose of the measurement operator, The covariance matrix of diagonal observation errors is used to represent the coarse-resolution observation uncertainty. and high-resolution observation uncertainty Construction, diagonal elements correspond to high-resolution observation uncertainties The value of the observation error covariance matrix at that location is Diagonal elements correspond to coarse-resolution observation uncertainty The value of the observation error covariance matrix at that location is .

[0045] In a specific embodiment of the present invention, the Kalman gain matrix is ​​utilized. By correcting the prediction set with the perturbation observation set, we obtain the forward analysis set at time t. As shown in equation (14).

[0046] (14) in: Let be the forward analysis set at time t. For the perturbation observation set, , For the real number domain, Gaussian noise is added to the observed NDVI values ​​according to their respective variances.

[0047] S4.3.7: Using the backward propagation algorithm, perform backward smoothing updates on all historical times less than t to generate backward analysis sets for each historical time. ; S4.3.8: Calculate the mean of the forward analysis set and the mean of the backward analysis set for each historical time point; S4.3.9: The forward analysis set mean and the backward smoothing set mean are weighted and fused to obtain a high-resolution normalized differential vegetation index time series.

[0048] In a specific embodiment of the present invention, high-resolution data is used... Coarse resolution data High-resolution observation error By constructing a set of disturbance observations and the diagonal observation error covariance matrix Proceed to the analysis step and further influence the Kalman gain. The calculation is used to adjust the magnitude of the model state correction to high-resolution observation, thereby achieving weighted assimilation of high-resolution NDVI information. The low-resolution observation error covariance matrix is ​​calculated by the deviation between the 30-meter MODIS obtained by downscaling and the high-resolution observation, and is transmitted over time t through the rate of change. The high-resolution observation error covariance is assigned a value of 0.05 when high-resolution data is available.

[0049] To enhance temporal consistency and high reliability of utilizing high-resolution data, the proposed algorithm FSATS includes a backward smoothing step. When a high-resolution observation is available at time t, it is used not only to update the current state but also to refine the previous state by performing a backward pass using equation (15). The state at any given moment, This refers to any historical moment before the current high-resolution observation time t.

[0050] When using coarse-resolution observations, the specific method for obtaining high spatial resolution normalized differential vegetation index time series is as follows: When coarse-resolution observations are acquired, the current prediction set is updated based on the uncertainty of the coarse-resolution observations to obtain the analysis set. The mean of the analysis set is taken as the current high spatial resolution normalized differential vegetation index value. No backward smoothing is performed. The analysis set is then integrated step by step to form the high spatial resolution normalized differential vegetation index time series.

[0051] In a specific embodiment of the present invention, the backward analysis set This supplements the forward analysis set obtained during the forward filtering process. (The NDVI state set is updated hourly using ensemble Kalman filtering to obtain an analysis set that integrates observation information. This analysis set is then used as the initial condition for model advancement at the next time step, thus propagating observation constraints over time and providing a consistent state trajectory for subsequent smoothing stages.) The smoothed ensemble mean is then calculated using a time-weighted average of the forward and backward means. As shown in equation (15): (15) in: The smoothed set mean , As a weighting factor, for Forward analysis of the set mean at each time step for Backward analysis set mean at time step.

[0052] In Equation 15, the set refers to the set of state members (Ensemble) under the Kalman Filter (EnKF) framework. The value ranges from 0 to 1, depending on... With available high-resolution observation time and ( The temporal proximity of ) is shown in equation (16): (16) in: For any historical moment before the current high-resolution observation time t, For the first available high-resolution observation time, This is the second available high-resolution observation time.

[0053] When using coarse-resolution observations, the specific method for obtaining high spatial resolution normalized differential vegetation index time series is as follows: When coarse-resolution observations are acquired, the current prediction set is updated based on the uncertainty of the coarse-resolution observations to obtain the analysis set. The mean of the analysis set is taken as the current high spatial resolution normalized differential vegetation index value. No backward smoothing is performed. The analysis set is then integrated step by step to form the high spatial resolution normalized differential vegetation index time series.

[0054] To accurately monitor complex phenological signals in heterogeneous agricultural landscapes and overcome the limitations of temporal sparsity in high-resolution remote sensing data, this invention... Figure 3The study shown is located in Ruzhou, Henan Province, China (33°56′–-34°20′N, 112°31′–-113°07′E), situated in the core region of the Huang-Huai-Hai Plain. This is one of the most agriculturally active and phenologically diverse regions in China, exhibiting a heterogeneous landscape composed of single-cropping (e.g., winter wheat or summer maize) and double-cropping systems (e.g., winter wheat-summer maize rotation), occasionally interspersed with cash crops such as cotton. These planting patterns are closely coupled with seasonal climate variability and irrigation practices, resulting in complex and spatially dispersed phenological signals. Importantly, this region is also one of the few in China with publicly available high-quality field phenological observations from 2008, enabling rigorous validation of fine-scale phenological reconstructions. This combination of phenological complexity and reliable ground observations makes Ruzhou an ideal benchmark location for evaluating the spatiotemporal consistency and phenological accuracy of spatiotemporal fusion algorithms. A key challenge in monitoring such heterogeneous agricultural landscapes is the limited temporal availability of high-resolution remote sensing data. Although the nominal revisit interval for the Landsat series is 16 days, cloud contamination and atmospheric interference significantly reduce the effective number of available images for the region. For example, in 2008, only 5 out of 17 Landsat images had cloud cover below 10%, while 12 had cloud cover above 18%, and 10 had cloud cover above 40%. Similarly, in 2007, only 2 out of 13 available images met the <20% cloud cover threshold. These statistics highlight the severe temporal sparsity of high-resolution optical imagery over Ruzhou, underscoring the need for robust spatiotemporal fusion techniques to fill observational gaps and support time-sensitive phenological analysis. To overcome these limitations, this invention employs a multi-source data fusion strategy for 2008, combining historical Landsat observations with medium-resolution MODIS data. 2008 was chosen primarily because contemporaneous field phenological records were available, which is crucial for quantitatively assessing the accuracy of phenological parameter extraction. Fifty-three cloud-free, atmospherically corrected Landsat TM surface reflectance images (1995–2008) were acquired from the U.S. Geological Survey Earth Resources Observation and Science Center via the ESPA platform (https: / / espa.cr.usgs.gov / ). All images were converted to NDVI and projected onto a common UTM / WGS-84 coordinate system. Land cover in the study area was classified into five main types: cultivated land, woodland, grassland, built-up land, and others (e.g., water bodies). Simultaneously, the 2008 MOD09Q18 surface reflectance composite product was converted to an NDVI time series and resampled to a spatial resolution of 240 meters, corresponding to an aggregation of 8 × 8 Landsat 30-meter pixels, ensuring spatial compatibility of the spatiotemporal data fusion. Figure 4The day-series distribution of selected Landsat NDVI imagery from 1995 to 2007 is presented. For the spatiotemporal fusion experiment, Landsat-MODIS image pairs for four key phenological periods (DOY 81, 209, 241, and 353) in 2008 were selected as base images. These dates correspond to key growth stages: DOY 81 and 353 correspond to the jointing and emergence stages of winter wheat; DOY 209 and 241 capture the jointing and tasseling stages of summer maize. To evaluate prediction performance, three additional Landsat images (DOY 65, 129, and 289) were used as independent test scenes. These test dates reflect key seasonal transitions for winter wheat, including greening, tasseling, and post-harvest sowing, making them ideal for evaluating the temporal accuracy and phenological sensitivity of the fused NDVI output.

[0055] To comprehensively evaluate FSATS, it was benchmarked against three widely used fusion algorithms under consistent input conditions: ESTARFM (see X. Zhu, J. Chen, F. Gao, X. Chen, and JG Masek, "An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions," Remote Sensing of Environment, vol. 114, no. 11, pp. 2610-2623, 2010), LMGM (see Y. Rao, X. Zhu, J. Chen, and J. Wang, "An Improved Method for Producing High Spatial-Resolution NDVI Time Series Datasets with Multi-Temporal MODIS NDVI Data and Landsat TM / ETM plus Images," Remote Sensing, vol. 7, no. 6, pp. 7865-7891, 2015), and PK-STEM (see Y. Wang, G. Yan, ...). D. Xie, R. Hu, and H. Zhang, "Generating Long Time Series of HighSpatiotemporal Resolution FPAR Images in the Remote Sensing Trend SurfaceFramework," IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp.1-15, 2021), all algorithms were used to generate 30-meter NDVI time series at 8-day intervals throughout 2008, using the same base image pairs.For ESTARFM and LMGM, the nearest-date base image strategy was adopted, that is, using the base images of DOY 081 and 209 to predict the image of DOY 129, as described by Wang et al. (see Y. Wang, G. Yan, D. Xie, R. Hu, and H. Zhang, "Generating Long Time Series of High SpatiotemporalResolution FPAR Images in the Remote Sensing Trend Surface Framework," IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp. 1-15, 2021). The spatial window size of ESTARFM was set to 5×5 Landsat pixels, and all other parameters followed the default configuration specified by Zhu et al. (see X. Zhu, J. Chen, F. Gao, X. Chen, and JG Masek, "An enhanced spatial and temporal adaptive reflectance fusion model for complex heterogeneous regions," Remote Sensing of Environment), while LMGM followed the parameters specified by Rao et al. (see Y. Rao, X. Zhu, J. Chen, and J. Wang, "An Improved Method for Producing High Spatial-Resolution NDVI Time Series Datasets with Multi-Temporal MODIS NDVI Data and Landsat TM / ETM plus Images," Remote Sensing, vol. 7, no. 6, pp. 7865-7891). For PK-STEM and FSATS, the prior NDVI time series (RSTS) data were constructed from 46 historical Landsat images (1995–2007) using DOY synthesis. Temporal filtering was then applied to smooth the synthesized NDVI trajectory. Specifically, PK-STEM employed a DL function to preserve key phenological inflection points, while FSATS used a Hants filter, resulting in average RMSEs of 0.090 and 0.084, respectively.The implementation of PK-STEM strictly follows the original configuration described by Wang et al. (see Y. Wang, G. Yan, D. Xie, R. Hu, and H. Zhang, "Generating LongTime Series of High Spatiotemporal Resolution FPAR Images in the RemoteSensing Trend Surface Framework," IEEE Transactions on Geoscience and RemoteSensing, vol. 60, pp. 1-15, 2021). FSATS requires the computation of spatially explicit initial state errors and observation error statistics, with the perturbation matrix derived using dynamic model estimation and the standard deviation of observations, respectively. The key innovation of FSATS lies in its spatiotemporally adaptive uncertainty modeling. Unlike PK-STEM, which assumes constant spatial invariance uncertainty, FSATS estimates dynamically propagated uncertainty based on local pixel-level errors, achieving fine-grained control over prediction reliability. Figure 5 (a) to (d) illustrate the heterogeneity of the uncertainty distributions at DOY=001, 065, 129 and 289.

[0056] Both ensemble-based methods (PK-STEM and FSATS) used a diagonal error covariance matrix with an ensemble size of N = 100. FSATS was tested in three observation modes: Mode I using only MODIS NDVI observations, Mode II using only Landsat NDVI observations, and Mode III integrating MODIS and Landsat NDVI observations.

[0057] To rigorously evaluate performance, this invention conducted three experiments: 1) intra-FSATS comparison (modes I-III) using MODIS, Landsat, and joint observations; 2) inter-algorithm comparison between FSATS (mode III), ESTARFM, LMGM, and PK-STEM.

[0058] All experiments used the same input: four Landsat base images (DOY = 81, 209, 241, 353) and three independent test images (DOY = 65, 129, 289).

[0059] Performance metrics include: correlation coefficient (R), root mean square error (RMSE), and mean absolute error (MAE).

[0060] The correlation coefficient (R) is calculated as follows: (17) in: The correlation coefficient is... The total number of pixels, Number of pixels , The actual NDVI value. This represents the true mean of the normalized vegetation index. The predicted normalized vegetation index value. The mean of the predicted normalized vegetation index values. The standard deviation of the true normalized vegetation index value. This represents the standard deviation of the predicted normalized vegetation index value.

[0061] The root mean square error is calculated as follows: (18) in: This is the root mean square error; The Mean Absolute Error (MAE) is calculated as follows: (19) 1) Prediction results of FSATS under different modes Figure 6 (a)–(j) use spatial visualization, accuracy metrics, and scatter plots to comprehensively compare FSATS predictions across three operating modes. Figure 6 In Model III shown in (d), FSATS reconstructed the heterogeneous landscape with high fidelity, closely approximating the actual Landsat NDVI of DOY 289. This included clearly depicted field boundaries, roads, and varying vegetation structure. In contrast, Figure 6 (b) The pattern I shown lacks sufficient spatial detail, obscuring the transition between cultivated and non-cultivated land, while Figure 6 (c) shows that Pattern II has a significant overestimation of NDVI, especially in high biomass areas.

[0062] These differences stem from the input constraints of each mode: Mode I is limited to coarse-resolution MODIS input, Mode II relies solely on sparse Landsat data, while Mode III combines two data sources, leveraging their respective strengths.

[0063] From a quantitative perspective, Model III consistently achieved the best performance across the three validation dates (DOY 065, 129, and 289), with RMSE values ​​ranging from 0.053 to 0.116, MAE values ​​from 0.039 to 0.098, and Pearson correlation coefficients (R) from 0.806 to 0.915. Model I performed moderately, while Model II performed the worst, particularly at DOY 289 (R = 0.426), reflecting its inability to track rapid phenological changes during aging.

[0064] Overall, the prediction accuracy was ranked as Mode III > Mode I > Mode II. Evaluated across all test imagery, Mode III achieved the highest R-value (0.838) and the lowest RMSE (0.092) and MAE (0.071), reinforcing the advantages of combining MODIS and Landsat data. Figure 6 (h)-(j) show that in Mode III, the data points are closely clustered near the 1:1 reference line, indicating minimal deviation and strong consistency.

[0065] Figure 7 (a)–(c) show the NDVI temporal curves for representative cultivated land, woodland, and grassland. In all cases, FSATS produces smooth NDVI tracks due to internal temporal filtering. Mode I closely approximates the MODIS trend but deviates from the Landsat peak; Mode II better matches the Landsat peak but fails in woodland areas due to limited Landsat NDVI data; Mode III captures both temporal dynamics and fine-resolution details, showing strong consistency with both MODIS and Landsat NDVI imagery.

[0066] 2) Comparison of FSATS with ESTARFM, LMGM and PK-STEM Figure 8 (a)–(j) compare the predictions of DOY289 using the base image pairs DOY 241 and DOY 353 by FSATS, ESTARFM, LMGM, and PK-STEM. Visually, FSATS predictions are closest to the actual Landsat NDVI, especially in areas of rapidly aging farmland (highlighted by black circles). In contrast, ESTARFM and LMGM fail to capture such abrupt temporal changes due to their linear dynamics assumptions. PK-STEM also performs poorly at longer temporal distances from the base image because it relies solely on statistics derived from MODIS. FSATS shows a significant advantage in preserving fine-scale spatial heterogeneity and tracking abrupt drops in farmland NDVI (e.g., Figure 8(The black circle in the image highlights aging fields). In contrast, ESTARFM and LMGM, assuming linear or piecewise linear transitions, underestimate NDVI variations and obscure spatial details, failing to capture abrupt temporal changes in aging fields. PK-STEM, relying solely on MODIS-derived statistics and temporally distant baseline imagery, also failed to recover short-term dynamics.

[0067] Figure 9 (a)–(p) show scatter plots of predicted NDVI versus actual NDVI for all three validation dates. FSATS achieved the highest overall accuracy, with R-values ​​consistently above 0.806 and RMSE below 0.116. In DOY 289, FSATS recorded R = 0.839 and RMSE = 0.116, outperforming all other methods. In DOY 129, both FSATS and PK-STEM performed well, benefiting from the relative homogeneity of vegetation during peak growth. In DOY 065, results from all methods were comparable due to low vegetation variability and close proximity of baseline image pairs.

[0068] Figure 10 The mean and standard deviation of NDVI generated by the four methods across the study area were evaluated. MODIS and prior NDVI signals showed a typical bimodal pattern typical of the winter wheat-summer maize system, while Landsat observations were sparse and temporally distant. Nevertheless, FSATS, ESTARFM, and LMGM recovered the double-season structure. PK-STEM failed to clearly distinguish the crop cycle. LMGM underestimated the NDVI of the first season, likely due to unmixing errors. Notably, FSATS produced a smoother and more accurate curve than ESTARFM, indicating its superior performance in capturing complex temporal patterns.

[0069] Figure 11 The pixel-level time series in (a)–(c) show that FSATS consistently aligns with Landsat observations and MODIS temporal curves across all land cover types. ESTARFM and LMGM exhibit greater temporal fluctuations, while PK-STEM and FSATS produce clearer, more interpretable curves. FSATS performs best in distinguishing seasonal boundaries for bi-seasonal vegetation.

[0070] This invention compares the performance of four representative spatiotemporal fusion methods in reconstructing vegetation dynamics in complex agricultural landscapes: ESTARFM, LMGM, PK-STEM, and FSATS. These methods represent three methodological categories within the STF field: (i) weight-based fusion models (e.g., ESTARFM), (ii) spectral unmixing-based models (e.g., LMGM), and (iii) trend surface-based models (RSTS) (e.g., PK-STEM and FSATS). The comparison highlights their comparative advantages and limitations in addressing spatial heterogeneity, phenological complexity, data scarcity, and computational efficiency.

[0071] As a representative weighted method, ESTARFM exhibits relatively stable performance in the spatial dimension and produces reasonable temporal profiles when vegetation dynamics follow a gradual and predictable trend. However, its reliance on the assumption of linear variation between observation dates limits its ability to capture rapid phenological transitions, such as abrupt aging or crop rotation. Furthermore, the lack of an explicit temporal smoothing mechanism often results in temporally noisy or discontinuous NDVI trajectories, especially under sparse, high-resolution observation conditions.

[0072] LMGM, as a linear unmixing-based method, exhibits moderate performance in terms of spatial accuracy and temporal consistency, with its main advantage being computational efficiency. However, the inherent uncertainties in linear unmixing within mixed-cell environments lead to unstable predictions and amplified temporal noise. Although subsequent smoothing can partially mitigate these effects, initial temporal biases often propagate into phenological parameter estimations, reducing overall reliability.

[0073] Both PK-STEM and FSATS belong to the RSTS class of methods, which integrate long-term high-resolution imagery to guide fusion under data-scarce conditions. PK-STEM assumes a single-season phenological structure and heavily relies on MODIS-derived statistics to modulate high-resolution temporal trajectories. While this strategy is effective in homogeneous single-season planting areas, its performance deteriorates in heterogeneous agricultural landscapes where single-season and double-season planting systems coexist.

[0074] FSATS represents a significant generalization and advancement of the RSTS concept. By introducing an adaptive trend surface that uses historical NDVI signals to distinguish between single-season and double-season vegetation, FSATS relaxes the restrictive phenological assumptions embedded in earlier RSTS-based methods. More importantly, FSATS incorporates a unified ensemble Kalman smoother with a model of spatiotemporal uncertainty, enabling robust assimilation of both high-resolution and coarse-resolution observations. This design allows FSATS to effectively handle asynchronous sub-pixel phenological dynamics under sparse and irregular observation conditions.

[0075] From a performance perspective, FSATS consistently achieved the highest accuracy in spatial, temporal, and phenological indices (RMSE = 0.092, R = 0.838). The results indicate that FSATS is particularly effective in capturing complex phenological transitions that are often oversimplified or missed by existing STF methods.

[0076] Because it relies on long-term, high-resolution prior knowledge and adaptive harmonic analysis, FSATS holds promise for applicability to a wide range of agricultural regions dominated by single- and double-cropping systems, which constitute the majority of the world's arable land. In regions with more complex phenological patterns, such as triple-cropping systems, the current framework remains scalable, although additional harmonic components or region-specific parameterization may be required to adequately represent multiple growing seasons.

[0077] In terms of computational efficiency, FSATS achieves a favorable balance between model complexity and operational feasibility. Similar to PK-STEM, the most time-consuming step in FSATS is pixel-by-pixel harmonic fitting to construct the prior trend surface (which takes approximately one minute in the experiments of this invention). Subsequent fusion and update steps are performed on the full high-resolution image using matrix-based EnKS computation, which is computationally efficient, taking approximately 30 seconds for a one-year forecast with an 8-day interval. While FSATS is less efficient than unmixing-based methods (such as LMGM) that operate only at coarse resolution, it is far more efficient than weight-based methods (such as ESTARFM) for dense time series reconstruction.

[0078] FSATS is particularly well-suited for long-term vegetation reconstruction tasks, where the cost of historical data preparation can be amortized over extended time series. For short-term studies or those relying on single image pairs, methods like ESTARFM or LMGM may be more practical due to their minimal preprocessing requirements. This underscores that FSATS is a robust and scalable solution for long-term, regional-scale phenological monitoring under conditions of persistent data scarcity.

[0079] In summary, FSATS is not limited to specific vegetation indices or sensor configurations. Although this invention focuses on NDVI-based fusion, the framework is generally applicable to other vegetation parameters as long as appropriate observation operators and uncertainty models are defined. Due to its sensor-independent design, FSATS can be easily extended to integrate Sentinel-2 and Sentinel-3 observations, including red-edge-based vegetation indices, which have great potential for enhancing phenological monitoring.

[0080] The various embodiments in this specification are described in a related manner. Similar or identical parts between embodiments can be referred to mutually. Each embodiment focuses on describing the differences from other embodiments. In particular, the system embodiments are basically similar to the method embodiments, so the description is relatively simple; relevant parts can be referred to the descriptions of the method embodiments.

[0081] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention are included within the scope of protection of the present invention.

Claims

1. A high-resolution normalized differential vegetation index time-series reconstruction method based on adaptive uncertainty, characterized in that, Includes the following steps: S1: Collect historical high-resolution normalized difference vegetation index (NDVI) images and coarse-resolution NDVI time series, preprocess them respectively to obtain a priori NDVI time series with fixed time intervals and a coarse-resolution NDVI time series. S2: Establish a dynamic state-space model based on the time series of prior normalized vegetation index with fixed time intervals; S3: Estimate the uncertainty of high-resolution observations, the uncertainty of dynamic state-space models, and the uncertainty of coarse-resolution observations, respectively; S4: Construct an ensemble Kalman smoother to perform multi-scale data assimilation and obtain the reconstructed high-resolution normalized differential vegetation index time series.

2. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 1, characterized in that, The specific steps for S1 are as follows: S1.1: Collect historical high-resolution normalized vegetation index images, compile images of the same day each year in chronological order to obtain a historical time series; S1.2: Harmonic analysis is performed on the historical time series to obtain the prior normalized differential vegetation index time series. As shown in equation (1); （1） in: For the prior normalized differential vegetation index time series. For the first The time points of each observation sample Here are the indices of the observed time points in the time series, j=1,2,...,N. The coefficient at zero frequency, The number of periodic terms in the series. For harmonic component indices, n=1, 2, ..., N f N is the number of observation samples in the time observation series. Harmonic frequency The corresponding coefficients of the first triangular component, For harmonic frequencies, Harmonic frequency The corresponding coefficients of the second triangular component, For frequency; S1.3: Resample the prior normalized differential vegetation index time series to obtain a prior normalized differential vegetation index time series with fixed time intervals; S1.4: Preprocess the coarse-resolution normalized differential vegetation index time series to obtain the coarse-resolution prior normalized differential vegetation index time series.

3. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 1, characterized in that, The specific steps for S2 are as follows: S2.1: Define the state equations for the dynamic state-space model, expressed as: （2） in: For time t, the high-resolution normalized differential vegetation index field is defined. For the high-resolution normalized differential vegetation index field at time t-1, Here is the transition matrix; S2.2: Based on the prior normalized differential vegetation index time series, the transition matrix is ​​calculated using the prior values ​​at time t and time t-1. Establish a dynamic state space model, as shown in equation (3), and obtain the dynamic state space model; （3） in: Let be the prior normalized differential vegetation index at time t. The normalized differential vegetation index at time t-1. It is the first constant. This is the dynamic transition matrix.

4. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 1, characterized in that, The specific steps for S3 are as follows: S3.1: Addressing the uncertainty of high-resolution observations Let it be a constant; S3.2: Calculate the uncertainty of the dynamic state-space model ; S3.3: Calculating the uncertainty of coarse-resolution observations ; S3.4: respectively for and Apply a lower bound constraint so that neither is less than .

5. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 4, characterized in that, The specific steps in S3.2 are as follows: S3.2.1: Obtain observations of various high-resolution normalized differential vegetation indices ; S3.2.2: Using the transition matrix The observation of each high-resolution normalized differential vegetation index By projecting backwards to the initial time, we obtain the implicit initial state. , represented as: （4） in: For high-resolution normalized differential vegetation index observation The implicit initial state For high-resolution normalized differential vegetation index observation, Let be the inverse matrix of the normalized differential vegetation index at time t; S3.2.3: Calculate all implicit initial states The standard deviation of is taken as the uncertainty of the dynamic state-space model and expressed as: （5） in: For high-resolution normalized differential vegetation index observation Quantity, The data collection time for high-resolution normalized difference vegetation index observations. The mean of the initial time-normalized differential vegetation index values ​​was reconstructed for all high-resolution normalized differential vegetation index observations. This refers to a specific observation time within the data collection period for high-resolution normalized difference vegetation index (NDVI) observations.

6. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 4, characterized in that, The specific steps in S3.3 are as follows: S3.3.1: Construct a coarse-resolution time transition matrix based on the time series of the coarse-resolution normalized differential vegetation index. , represented as: （6） in: This is a coarse-resolution time transition matrix. Let be the coarse-resolution normalized differential vegetation index value at time t. This represents the coarse-resolution normalized differential vegetation index value at time t-1. It is the second constant; S3.3.2: Utilizing a coarse-resolution transition matrix Observations of various high-resolution normalized differential vegetation indices Back-projecting to the initial time, we obtain the normalized differential vegetation index observations derived from the coarse-resolution time series. Implicit initial state , represented as: （7） in: For high-resolution normalized differential vegetation index observation The initial state, Let be the inverse of the coarse-resolution time transition matrix at time t. It is the inverse of the coarse-resolution time transition matrix at time t-1; S3.3.3: Calculate the rate-normalized differential vegetation index observations derived from all coarse-resolution time series. Implicit initial state The standard deviation, as a coarse-resolution observation uncertainty. , represented as: （8） in: The mean of the coarse-resolution normalized differential vegetation index values ​​at the initial time is used to dynamically reconstruct the dynamics of all high-resolution normalized differential vegetation index observations. For high-resolution normalized differential vegetation index observation Quantity, This refers to the data collection time for high-resolution normalized differential vegetation index observations.

7. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 1, characterized in that, The specific steps for S4 are as follows: S4.1: Generate the normalized differential vegetation index value corresponding to the initial time using a dynamic state-space model, and add an unbiased Gaussian distribution perturbation. Generate a set and the set As the initial time analysis set ; S4.2: For each subsequent time step, utilize the dynamic transition matrix. The analysis set from the previous moment Forward propagation yields the prediction set for the current time step. ; S4.3: Determine whether there are any new observations at the current time. If it is a high-resolution observation or a coarse-resolution observation, calculate the Kalman gain using the corresponding observation uncertainty and update the prediction set to obtain the current analysis set. If there are no new observations, recursively obtain the normalized differential vegetation index value at the current time based on the predicted value at the previous time and the dynamic transition matrix. Integrate the results of each time to obtain the high spatial resolution normalized differential vegetation index time series.

8. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 7, characterized in that, When performing high-resolution observations, the specific steps to obtain a high spatial resolution normalized differential vegetation index time series are as follows: S4.3.1: Selecting High-Resolution Observation Uncertainties And calculate its square. As observation variance; S4.3.2: Adding variance to high-resolution observations Gaussian perturbation, generating a set of perturbation observations. ; S4.3.3: Using measurement operators The prediction set at the current moment Mapping to a high-resolution observation space yields the predicted observation set. ; S4.3.4: Calculate the prediction error covariance matrix And construct the diagonal observation error covariance matrix. ; S4.3.5: Calculate the Kalman gain matrix; S4.3.6: Update the prediction set at the current time step according to the Kalman filter. Updated to forward analytic set ; S4.3.7: Using the backward propagation algorithm, perform backward smoothing updates on all historical times less than t to generate backward analysis sets for each historical time. ; S4.3.8: Calculate the mean of the forward analysis set and the mean of the backward analysis set for each historical time point; S4.3.9: The forward analysis set mean and the backward analysis set mean are weighted and fused to obtain a high-resolution normalized differential vegetation index time series.

9. The method for time-series reconstruction of high-resolution normalized differential vegetation index based on adaptive uncertainty according to claim 7, characterized in that, When using coarse-resolution observations, the specific method for obtaining high spatial resolution normalized differential vegetation index time series is as follows: When coarse-resolution observations are acquired, the current prediction set is updated based on the uncertainty of the coarse-resolution observations to obtain the analysis set. The mean of the analysis set is taken as the current high spatial resolution normalized differential vegetation index value. No backward smoothing is performed. The analysis set is then integrated step by step to form the high spatial resolution normalized differential vegetation index time series.

Images