A remote sensing vegetation index time series missing value reconstruction method adaptive to smoothing and trend constraint cooperation

Abstract

The application discloses a kind of remote sensing vegetation index time series missing value reconstruction methods of adaptive smoothing and trend constraint coordination, comprising: obtaining the remote sensing vegetation index time series data of pre-processing, extract effective observation value and identify missing value to construct one-dimensional time series;The observation density of time series is calculated, and the smoothing parameter and trend fusion weight are determined accordingly;Based on smoothing parameter, an adaptive weighted Whittaker smoothing model is constructed to obtain a preliminary smoothing sequence by fitting;According to the number of effective observation value, the trend reference value at the missing value is calculated by adaptively selecting the trend estimation method;Based on the trend fusion weight, the preliminary smoothing sequence and the trend reference value are weighted and fused to generate the filling value, and only the missing value is replaced to obtain the preliminary reconstruction sequence;Multi-level physical constraint verification is applied to the filling value, and the final reconstruction result is output.The application overcomes the defects of excessive smoothing and abnormal interpolation of traditional methods, and realizes the unity of data authenticity, continuity and physical rationality.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing data processing, and in particular, it is a method for reconstructing time-series missing values ​​of remote sensing vegetation indices through the synergy of adaptive smoothing and trend constraints. Background Technology

[0002] Remote sensing vegetation index time series are key physical quantities characterizing the growth status of surface vegetation, playing an irreplaceable fundamental supporting role in continuous dynamic tasks such as agricultural monitoring, drought assessment, and habitat quality analysis. However, due to limitations in sensor observation conditions and the frequent influence of uncontrollable factors such as cloud cover and atmospheric interference, the actual acquired time series data generally contain a large amount of missing data and abrupt changes, disrupting the spatiotemporal continuity of the data and limiting the effectiveness of large-scale applications of high-precision dynamic surface monitoring.

[0003] Currently, the mainstream techniques for reconstructing such time-series data can be divided into two main categories. One category is interpolation methods based on local functions, such as linear interpolation and spline interpolation. These methods heavily rely on the data status of a few neighboring observation points when performing local calculations, making it difficult to accurately characterize the long-term nonlinear evolution of vegetation growth processes under complex natural environments. The other category is smoothing methods based on global or local filtering, such as Savitzky-Golay filtering and the standard Whittaker smoothing algorithm. These algorithms apply uniform filtering to the entire sequence using constant parameters. However, when faced with severely uneven spatial pixel-level data loss distributions, these statistical algorithms with fixed processing strength are prone to over-smoothing in densely populated peak and trough regions, thus weakening the true phenological inflection point signals. Conversely, in extremely sparse frequency bands, over-extrapolation can cause results to deviate from the true physical change trajectory.

[0004] In summary, existing processing methods still have certain limitations when dealing with large-scale, irregular, and inconsistent observation data sequences. Their ability to balance filtering environmental noise with preserving accurate phenological details needs further improvement, and their ability to synergistically characterize the overall continuity and local evolutionary features of the reconstructed results also needs enhancement. Therefore, it is necessary to further investigate a remote sensing vegetation index time series reconstruction method that can balance data accuracy, temporal smoothness, and physical plausibility. Summary of the Invention

[0005] The purpose of this invention is to provide a method for reconstructing time-series missing values ​​of remote sensing vegetation indices by combining adaptive smoothing and trend constraints, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solution: A method for reconstructing time-series missing values ​​of remote sensing vegetation indices through adaptive smoothing and trend constraint coordination, comprising:

[0007] Acquire preprocessed remote sensing vegetation index time series data, extract valid observations and identify missing values, and construct a one-dimensional time series;

[0008] Calculate the observation density of a one-dimensional time series and determine the smoothing parameter and trend fusion weights accordingly;

[0009] An adaptive weighted Whittaker smoothing model is constructed based on smoothing parameters, and a one-dimensional time series is fitted accordingly to obtain a preliminary smoothed sequence.

[0010] Based on the number of valid observations in the one-dimensional time series, conditional branching is determined, and the corresponding trend estimation method is adaptively selected to calculate the trend reference value at the missing value.

[0011] Based on trend fusion weights, the preliminary smoothed sequence and trend reference values ​​are weighted and fused to generate imputation values. The imputation values ​​are used to replace the missing values ​​in the one-dimensional time series to obtain the preliminary reconstructed sequence, while keeping the effective observation values ​​unchanged.

[0012] Apply multi-level physical constraint verification to the imputation values ​​in the initial reconstructed sequence, and output the reconstruction results of the temporal missing values.

[0013] Beneficial effects: This invention overcomes the defects of excessive smoothing and abnormal interpolation in traditional methods, and achieves a unity of data authenticity, continuity and physical rationality. Attached Figure Description

[0014] Figure 1 A flowchart illustrating the steps of a remote sensing vegetation index time-series missing value reconstruction method that combines adaptive smoothing and trend constraint, provided in an embodiment of this application.

[0015] Figure 2 A flowchart illustrating the steps for obtaining a preliminary smooth sequence, provided in an embodiment of this application.

[0016] Figure 3 A flowchart illustrating the steps for determining smoothing parameters and trend fusion weights in an embodiment of this application.

[0017] Figure 4 A flowchart illustrating the steps of calculating the trend reference value at the missing value using the adaptive selection of the corresponding trend estimation method provided in this application embodiment.

[0018] Figure 5 This is a schematic diagram of the technical process for a remote sensing vegetation index time-series missing value reconstruction method that combines adaptive smoothing and trend constraint, provided in an embodiment of this application.

[0019] Figure 6 Example time series diagram of remote sensing EVI reconstruction results at the pixel scale provided in this application embodiment. Detailed Implementation

[0020] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. 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 should fall within the scope of protection of the present invention.

[0021] It should be noted that the terms "first," "second," etc., used in this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0022] like Figure 1 As shown, a method for reconstructing time-series missing values ​​of remote sensing vegetation indices using a combination of adaptive smoothing and trend constraints includes the following steps:

[0023] We acquire preprocessed remote sensing vegetation index time series data, extract valid observations and identify missing values, and construct a one-dimensional time series.

[0024] In this embodiment, the remote sensing vegetation index time-series data can specifically be enhanced vegetation index data acquired by a medium-resolution imaging spectrometer. In practical applications, due to the complexity of the atmospheric environment, the acquired raw remote sensing images are often contaminated to varying degrees. Therefore, it is necessary to perform cloud removal processing on the raw images using the quality control bands inherent in the data, removing low-quality pixels affected by cloud cover, aerosol contamination, etc. The uncontaminated pixels retained after removal are the valid observation values. The locations of the removed pixels and locations where the sensor itself did not acquire observation data are uniformly marked as invalid missing values.

[0025] Furthermore, after completing the pixel screening in the spatial dimension, for each pixel coordinate with valid observations within the study area, its series of observation data in the temporal dimension are arranged in chronological order. This transforms the two-dimensional spatial image data into a pixel-by-pixel one-dimensional time series. The one-dimensional time series contains discretely distributed valid and missing values, providing a standardized underlying data structure for subsequent temporal reconstruction.

[0026] Calculate the observation density of a one-dimensional time series and determine the smoothing parameter and trend fusion weights based on the observation density.

[0027] Specifically, observation density is used to quantitatively characterize the completeness of the current pixel time series data. A higher observation density indicates better data quality, as does a larger number of valid observations in a time series; conversely, a lower observation density indicates better data quality. Based on the calculated observation density, the system dynamically determines smoothing parameters to control the subsequent smoothing intensity, as well as trend fusion weights to balance global structure and local features. This density-based dynamic adjustment mechanism enables the system to adaptively address data quality differences between different pixels caused by varying cloud coverage frequencies, avoiding the problems of insufficient smoothing in dense data areas or over-smoothing in sparse data areas caused by using fixed parameters.

[0028] An adaptive weighted Whittaker smoothing model is constructed based on smoothing parameters. The adaptive weighted Whittaker smoothing model is then used to fit a one-dimensional time series to obtain a preliminary smoothed sequence.

[0029] Specifically, after determining the smoothing parameters, they are substituted into a pre-built model framework to form an adaptive weighted Whittaker smoothing model. During the fitting process, the adaptive weighted Whittaker smoothing model can flexibly adjust the penalty for changes in data curvature based on the overall data quality of the pixel, i.e., the smoothing parameters.

[0030] Furthermore, the one-dimensional time series containing missing values ​​is input into the adaptive weighted Whittaker smoothing model for global solution, outputting a continuous fitted curve. The values ​​corresponding to each time node on the fitted curve are extracted to obtain the preliminary smoothed sequence. The preliminary smoothed sequence exhibits good overall smoothness and reflects the global phenological cycle characteristics of vegetation growth.

[0031] Conditional branching is determined based on the number of valid observations in the one-dimensional time series, and the corresponding trend estimation method is adaptively selected to calculate the trend reference value at the missing value.

[0032] In this embodiment, the abundance of valid data varies greatly within the local time period where the missing values ​​are located. Blindly using the same method to extract local trends can easily lead to overfitting or underfitting. Therefore, conditional branching logic is introduced.

[0033] Specifically, the system counts the total number of valid observations throughout the one-dimensional time series. If the total number of valid observations is sufficient, indicating that the data contains enough information to support complex global statistical inferences, the system will select a trend estimation method that can utilize global covariance information to calculate the potential baseline value for the missing location. If the number of valid observations is extremely small, indicating severe data scarcity, the system will downgrade to a simpler local neighborhood interpolation method, relying only on a few valid data points before and after the missing point to infer the local trend. The calculated potential baseline value is the trend reference value, mainly used to characterize the local true change trend of vegetation within a specific missing period.

[0034] Based on trend fusion weights, the preliminary smoothed sequence and trend reference values ​​are weighted and fused to generate imputed values. The imputed values ​​are then used to replace the missing values ​​in the one-dimensional time series to obtain a preliminary reconstructed sequence, while keeping the effective observation values ​​unchanged.

[0035] In other words, based on the preliminary smoothed sequence and trend reference values, imputation values ​​are generated at the missing values ​​by combining trend fusion weights. The imputation values ​​are used to replace the missing values ​​in the one-dimensional time series to obtain a preliminary reconstructed sequence while keeping the effective observation values ​​unchanged.

[0036] Specifically, after obtaining the preliminary smoothed sequence reflecting the global structure and the trend reference value reflecting local characteristics, the two are weighted and summed using trend fusion weights. Through the fusion mechanism, the generated imputed value retains the overall smooth trend while reasonably approximating the actual direction of local changes.

[0037] It should be noted that this embodiment follows a data fidelity strategy. After generating imputation values, they are only used to replace missing values ​​at corresponding positions in the original one-dimensional time series. Valid and reliable observations in the original data are not modified or overwritten. The new sequence formed after the replacement is the preliminary reconstructed sequence, at which point the sequence has continuity in the time dimension.

[0038] Apply multi-level physical constraint verification to the imputation values ​​in the initial reconstructed sequence, and output the reconstruction results of the temporal missing values.

[0039] In this embodiment, although the initial reconstructed sequence achieves temporal continuity, due to the limitations of purely mathematical algorithms, some fill-in values ​​may deviate from the physical principles of vegetation growth, such as producing negative anomalies or abnormal peaks and troughs. Therefore, it is necessary to introduce prior physical knowledge of vegetation indices and statistical patterns from current pixel historical observations to construct multi-level constraint boundaries.

[0040] Specifically, each imputed value in the initial reconstructed sequence is fed into the constraint boundary for verification. If a imputed value is found to exceed physical norms or statistical limits, it is forcibly restricted to a reasonable boundary range, or it is gently morphologically repaired based on the morphology of nearby real mountain peaks. After all post-processing verification and correction, the final output is a high-quality temporal missing value reconstruction result, which can be directly used for downstream tasks such as agricultural monitoring or ecological analysis.

[0041] like Figure 2 As shown, in one possible implementation, an adaptive weighted Whittaker smoothing model is constructed based on smoothing parameters. This model is then used to fit a one-dimensional time series to obtain a preliminary smoothed sequence, including:

[0042] Different weighting coefficients are assigned to valid observations and missing values ​​in a one-dimensional time series, with the weighting coefficient for valid observations being greater than that for missing values.

[0043] Alternatively, different weighting coefficients are assigned to the valid observations and missing values ​​in a one-dimensional time series, with a larger weighting coefficient assigned to the valid observations and a zero weighting coefficient assigned to the missing values.

[0044] In a preferred implementation, when performing differential weighting, different weight coefficients are assigned to valid observations and missing values ​​in the one-dimensional time series. Specifically, valid observations are assigned a weight coefficient greater than one, and missing values ​​are assigned a weight coefficient equal to zero, so that missing values ​​do not produce a fidelity penalty in smooth fitting.

[0045] In other words, assign a larger weight coefficient to valid observations and a weight coefficient equal to zero to missing values.

[0046] Specifically, in a time series, valid observations represent the actual physical state acquired by remote sensing sensors, while missing values ​​indicate cloud cover or data failure at that time point. By constructing a diagonal weight matrix of the same length as the time series, the diagonal elements corresponding to valid observations are set to preset values ​​greater than 1. For example, the preset value can be specifically set to 10, mathematically constraining and forcibly improving the absolute fidelity of the computational model to the original real observation data. Simultaneously, the diagonal elements corresponding to missing values ​​are set to 0. When the weight coefficient equals 0, missing values ​​are cleared from the residual term of the objective function, thus preventing invalid numerical placeholders from interfering with the overall trend of the smooth curve.

[0047] Optionally, in some specific implementations, the weighting coefficients assigned to valid observations can be dynamically mapped based on the confidence score of the high-quality control band corresponding to that pixel. Specifically, a value of 10 is assigned to the observation with the best quality, and a weight between 2 and 5 is assigned to observations with slight quality interference, further refining the control granularity of data fidelity.

[0048] Construct a difference matrix to constrain temporal curvature changes.

[0049] Alternatively, based on the length of a one-dimensional time series, a difference matrix is ​​constructed to constrain the change in temporal curvature.

[0050] In this embodiment, a difference matrix is ​​constructed to constrain the curvature changes of the time series, specifically a second-order difference matrix. The smoothness of a time series signal can be objectively measured spatially by the curvature of the line connecting adjacent data points. In the specific implementation of the algorithm, the system constructs a first-order difference matrix to extract the linear slope between adjacent time nodes. Based on the first-order difference matrix, a second-order difference matrix is ​​further constructed to calculate and extract the rate of change of the slope, which is the discrete curvature of the time series. Each row of the second-order difference matrix corresponds to a local curvature term in the one-dimensional time series, and its non-zero elements on the main diagonal and adjacent diagonals are closely distributed according to a specific difference coefficient rule.

[0051] Furthermore, if the monitored physical object exhibits extremely high-frequency abrupt changes, as an alternative, the system can also construct a third-order difference matrix to forcibly constrain higher-order volatility. However, in specific application scenarios targeting remotely sensed vegetation phenological cycles, constructing a second-order difference matrix achieves the optimal technical balance between filtering out high-frequency atmospheric environmental noise and preserving the true inflection points of crop growth periods.

[0052] A smoothing objective function is constructed by combining weighting coefficients, difference matrices, and smoothing parameters, and then the smoothing objective function is transformed into a sparse symmetric positive definite linear equation system.

[0053] In other words, we construct a smoothing objective function by combining weighting coefficients, difference matrices, and smoothing parameters, and transform the optimization problem into a standard quadratic form.

[0054] Specifically, the principle for constructing the smoothing objective function lies in establishing an algebraic balance between the data fidelity term and the curvature penalty term. The system extracts the difference between the initial smoothed sequence and the one-dimensional time series to be solved, and performs matrix multiplication with the diagonal weight matrix to calculate the weighted sum of squares, thus forming the fidelity term. At the same time, it uses the second-order difference matrix to perform matrix multiplication on the initial smoothed sequence to extract the continuous curvature at time nodes, and combines it with preset smoothing parameters to calculate the sum of squares of curvature, thus forming the penalty term.

[0055] Etotal =E fidelity +λ*E smooth ;

[0056] Among them, E total To smooth the total computational error of the objective function, E fidelity E represents the fidelity error term strongly constrained by the diagonal weight matrix, λ is the smoothing parameter driven by the observation density, and E smooth This is the time-series curvature penalty error term constrained by the second-order difference matrix.

[0057] By taking the partial derivative of the smoothing objective function and setting the derivative to zero, the system transforms the unconstrained optimization problem of finding the minimum total error into a standard system of linear algebraic equations. Combining the diagonal weight matrix formed by the weighting coefficients with the smoothing parameters, the coefficient matrix of the transformed sparse symmetric positive definite linear equations exhibits a banded sparse distribution. This banded sparse distribution ensures that the non-zero elements in the coefficient matrix are highly concentrated on a few diagonals near the main diagonal, enabling direct and rapid solutions to the vast spatial dimensions of large-scale remote sensing images, while reducing memory usage during computation at the physical hardware level.

[0058] A direct solution method is used to numerically solve the coefficient matrix of a sparse symmetric positive definite linear equation system, and a preliminary smooth sequence is output.

[0059] Alternatively, the direct solution method can be used to numerically solve the coefficient matrix involved in the standard quadratic form and output a preliminary smooth sequence.

[0060] In this embodiment, for linear equation systems exhibiting a sparse, banded distribution, the computation module directly calls the sparse matrix direct solver tool in the numerical algorithm library to obtain accurate numerical solutions. For example, the system can use a banded solver based on Cholsky decomposition to perform rapid downward decomposition and upward substitution operations on sparse symmetric positive definite matrices. Compared to iterative solver algorithms such as the conjugate gradient method that rely on convergence thresholds, the direct solver method avoids, from its underlying mathematical principles, the abnormal states of slow convergence or divergent non-convergence that may occur with highly incomplete and complex time-series data, ensuring the operational stability of the platform's pixel-by-pixel processing mechanism. The final result vector output by the numerical solution is the preliminary smoothed sequence. The preliminary smoothed sequence not only achieves basic filling of missing positions in the original time series but also presents a smooth and coherent waveform that conforms to the natural growth and change patterns of surface vegetation over the entire annual time span.

[0061] For example, a difference matrix is ​​constructed to constrain the temporal curvature variation. Specifically, a second-order difference matrix is ​​constructed and combined with weighting coefficients and smoothing parameters so that the coefficient matrix of the transformed sparse symmetric positive definite linear equation system exhibits a banded sparse distribution characteristic, in order to adapt to the direct solution of large-scale spatial dimensions.

[0062] In some possible embodiments, constructing a difference matrix to constrain the temporal curvature variation can also be: constructing a second-order difference matrix, and combining weight coefficients and smoothing parameters to transform the optimization problem into a standard quadratic form, and employing a sparse matrix solution strategy to adapt to the computational efficiency of large-scale spatial dimensions of time series data.

[0063] like Figure 3 As shown, in an exemplary embodiment, the observation density of a one-dimensional time series is calculated, and the smoothing parameter and trend fusion weight are determined based on the observation density, including:

[0064] Extract the total length of the one-dimensional time series, calculate the ratio of the number of valid observations to the total length, and use the ratio as the observation density.

[0065] Specifically, when processing remote sensing vegetation index time-series data, the inconsistent frequency of cloud and rain weather impacts on different pixels leads to differences in the amount of effective data distributed along the time axis for time series corresponding to different pixels. To objectively quantify this data completeness, the system extracts the total length of the currently processed one-dimensional time series in the time dimension. It then iterates through the one-dimensional time series, counting the absolute number of valid observations not marked as missing values. Based on this, it calculates the ratio of the number of valid observations to the total sequence length using a division operation, defining this ratio as the observation density. The specific numerical calculation relationship is expressed as follows:

[0066] ρ = N / T';

[0067] Where ρ is the observation density, N is the number of valid observations in the one-dimensional time series, and T' is the total length of the one-dimensional time series.

[0068] In some optional implementations, for time series with extremely long spans, the system can also use a sliding time window approach to calculate local observation density. Specifically, a local time window of fixed length is extracted centered on the currently missing value to be filled, and the proportion of valid observations within the local time window is calculated as the local observation density, thereby more precisely quantifying the data missing distribution characteristics of the time series in different seasons or months.

[0069] Construct a continuous adjustment function with observation density as the independent variable.

[0070] In a preferred implementation, a continuous adjustment function with observation density as the independent variable is constructed, specifically as follows:

[0071] By introducing a preset adjustment index, a power function calculation model is constructed with the observation density as the base and the adjustment index as the exponent.

[0072] When calculating the smoothing parameters and trend fusion weights, different preset adjustment exponent constraint power function calculation models are used to control the degree of attenuation of the influence of observation density on smoothing intensity and local trend dependence.

[0073] In this embodiment, after acquiring the observation density, the system uses the observation density as an independent driving variable to construct a nonlinear mapping model for subsequent parameter calculations. Specifically, a preset adjustment exponent is introduced to construct a power function calculation model with the observation density as the base and the adjustment exponent as the exponent. Using a power function model as a continuous adjustment function leverages the curvature characteristics of the nonlinear function to smoothly transition the parameter response when the observation density changes slightly, avoiding the abrupt artifacts at different pixel boundaries caused by using a step-type piecewise function.

[0074] Furthermore, the continuous adjustment function forms a monotonically increasing mapping relationship from 0 to 1 in its value range. Based on this physical property, the higher the observation density, the closer the output value of the adjustment function is to 1; the lower the observation density, the closer the output value of the adjustment function is to 0. The positive monotonic mathematical characteristic provides a unified driving signal for the subsequent two-parameter coordinated decay.

[0075] The smoothing parameter is calculated by using a continuous adjustment function in combination with a pre-configured basic smoothing parameter, and the trend fusion weight is calculated by using a continuous adjustment function in combination with a pre-configured basic trend fusion weight, so that the smoothing parameter and the trend fusion weight can be positively coordinated to adjust with the change of observation density.

[0076] Specifically, when calculating the smoothing parameters, the system uses a pre-set global constant as the basic smoothing parameter and adopts a first preset adjustment exponent constraint power function calculation model. The specific numerical calculation relationship is expressed as follows:

[0077] λ=λ0*ρ γ ;

[0078] Where λ is the smoothing parameter actually applied to the smoothing model, λ0 is the pre-configured basic smoothing parameter, ρ is the observation density, and γ is the first preset adjustment index, preferably set to 0.5.

[0079] It should be noted that the pre-configured basic smoothing parameter λ0 is preferably set to 100, and the specific value of the basic smoothing parameter can be adaptively adjusted according to the sensor time resolution and time series length.

[0080] Based on this, the range of the first preset adjustment index is limited to greater than or equal to 0 and less than or equal to 1. When the observation density approaches 0, it indicates that the data observation is sparse, and the smoothing parameter calculated by the power function also approaches 0, thereby weakening the smoothing constraint of the system's underlying equations and avoiding over-smoothing that could lead to the loss of true phenological characteristics due to insufficient information. Conversely, when the observation density approaches 1, it indicates that the data observation is dense, and the calculated smoothing parameter tends to the basic smoothing parameter. In this case, the system enhances its noise resistance by increasing the smoothing strength.

[0081] Similarly, when calculating the trend fusion weights, the system uses an independent second preset adjustment exponent constrained power function calculation model, and the specific numerical calculation relationship is expressed as follows:

[0082] α=α0*ρ q ;

[0083] Wherein, α is the actual trend fusion weight applied to interpolation fusion; α0 is the pre-configured basic trend fusion weight, which can be set to 0.4; ρ is the observation density; and q is the second preset adjustment index, preferably set to 0.5.

[0084] The pre-configured base trend fusion weights are set to a range greater than or equal to 0 and less than or equal to 1. A second preset adjustment index is used to individually control the degree of attenuation of the influence of observation density on the trend fusion weights. By performing trend fusion weight calculations, when the observation density is high, the system outputs a larger trend fusion weight value, making the final reconstruction result more dependent on the local change trends generated by the branch extraction model. When the observation density is low, the trend fusion weight value output by the system decreases non-linearly according to a power-law, reducing the dependence on local interpolation results and causing the interpolated values ​​to converge more towards the global smoothing result to maintain overall structural stability. Accordingly, the system utilizes a unified underlying logic based on observation density to achieve an adaptive positive collaborative adjustment mechanism between smoothing intensity and local trend dependence at the spatial pixel scale.

[0085] like Figure 4 As shown, according to one aspect of this application, adaptively selecting a corresponding trend estimation method to calculate the trend reference value at the missing value includes:

[0086] The number of valid observations is compared with a preset observation threshold.

[0087] Specifically, in time series reconstruction tasks, a single mathematical model often cannot adapt to the drastic spatial fluctuations in data missing rates. Therefore, when processing each pixel sequence, the system extracts the total number of valid observations and compares it with a pre-set baseline parameter. The preset observation threshold is used to determine whether the current time series has the minimum information capacity to support a globally complex statistical algorithm. In specific parameter configurations, the preset observation threshold can be set to 8. By using comparison operations, the system introduces the computational process into different algorithmic branches, avoiding divergence problems caused by forcibly fitting complex models to extremely sparse data.

[0088] When the number of valid observations is greater than or equal to the observation threshold, Gaussian process regression is selected to fit the global trend, and trend reference values ​​at the missing values ​​are calculated.

[0089] Alternatively, when the number of valid observations is greater than or equal to the observation threshold, Gaussian process regression is used to fit the global trend using the valid observations, and trend reference values ​​are calculated at the missing values.

[0090] In this embodiment, when there are relatively sufficient valid observations retained in the sequence, the system determines that the current data has global statistical value, and at this time, the Gaussian process regression algorithm module is triggered. Gaussian process regression can utilize the temporal distribution characteristics of all valid points to provide more reasonable nonlinear trend extrapolation results. Through global fitting, the system predicts the potential optimal expected value for the current missing position to be filled, which serves as the reference coordinate to guide subsequent fusion calculations.

[0091] In a preferred embodiment, a Gaussian process regression is selected to fit the global trend, and trend reference values ​​at missing values ​​are calculated, including:

[0092] Construct a covariance matrix that reflects the correlation between each valid observation, and calculate the covariance vector between the current missing value to be filled and each valid observation.

[0093] For example, constructing a covariance matrix that reflects the correlation among the valid observations includes:

[0094] A squared exponential kernel function is introduced as the covariance function.

[0095] Specifically, before performing Gaussian process regression, the system needs to establish a mathematical tool to measure the correlation between different time points. This embodiment abandons linear or periodic kernels and directly introduces the squared exponential kernel function based on the characteristics of Gaussian distribution. The squared exponential kernel function assumes that the closer two time points are, the higher the correlation of their vegetation index characteristics, and that the correlation decays smoothly exponentially with time distance, which is consistent with the continuous physical evolution law of remotely sensed vegetation growth.

[0096] By combining the preset length scale parameter and the preset noise variance parameter, the time distance decay relationship between any two effective observation points in the time node sequence is calculated using the squared exponential kernel function, and the matrix elements in the covariance matrix are obtained. The length scale parameter is used to control the range of time correlation.

[0097] In this embodiment, the system extracts the time indices of any two valid observation points and substitutes them into the squared exponential kernel function for calculation. During this process, the preset length scale parameter determines the smooth span of the time correlation; a larger value indicates a longer memory effect of vegetation growth. The preset noise variance parameter is used to absorb the inherent instrument background noise in the observation data, preventing the model from overfitting to noisy points. In actual operational deployment, the preset length scale parameter can be set to 3, which physically corresponds to an influence period of approximately 48 days on the time axis; the preset noise variance parameter can be set to 0.01. The specific numerical calculation relationship is expressed as follows:

[0098] k val =exp(-1*(tt prime ) 2 / (2*l scale 2 ))+σ n 2 *δ val ;

[0099] Where, k val The covariance element values ​​between two time points, t and t0. prime These are the time node indices for two valid observation points, l scale σ is a preset length scale parameter. n σ is the noise standard deviation. n 2 δ is the preset noise variance parameter. val Let be the Kronecker function variable, when t equals t prime δ val The value is 1 if it is not 1, otherwise the value is 0.

[0100] Furthermore, the system iterates through all combinations of valid observation points and assembles a square matrix using the covariance element values, which is the covariance matrix. The covariance matrix characterizes the strength of cross-correlation within the valid observation data network on a global scale. Simultaneously, for the currently specified missing location to be filled, the system calculates the time distance decay between that location and all known valid observation points, forming a one-dimensional covariance vector. The covariance vector represents the weighted distribution network of the unknown point's influence by the globally known data.

[0101] To avoid directly inverting the covariance matrix, matrix decomposition is used to decompose the covariance matrix.

[0102] Specifically, in traditional statistical solutions, obtaining the predicted mean typically requires calculating the inverse of the covariance matrix. However, when the remote sensing time series is extremely long, directly calculating the inverse of a high-dimensional matrix is ​​not only computationally very complex but also prone to producing numerically unstable and anomalous floating-point results when the matrix condition number is poor. Therefore, the system uses a matrix factorization algorithm to decompose the constructed covariance matrix into a downward triangular matrix. The specific decomposition operation is expressed as follows:

[0103] K=L*L T ;

[0104] Where K is the covariance matrix, L is the lower triangular matrix obtained from the decomposition, and L T It is the transpose of the lower triangular matrix.

[0105] The first intermediate vector is obtained by combining the decomposed covariance matrix with the effective observations, and the second intermediate vector is obtained by combining the decomposed covariance matrix with the covariance vector.

[0106] In this embodiment, after obtaining the lower triangular matrix and its transpose, the system replaces the originally complex inversion operation with two efficient forward and backward substitution processes for solving the linear equation system. The system uses the lower triangular matrix and a column vector containing all original observation values ​​to extract the first intermediate vector, and uses the lower triangular matrix and the covariance vector to extract the second intermediate vector. The specific numerical solution relationship is expressed as follows:

[0107] L*α vec =Y vec ;

[0108] L*β vec =K star_vec ;

[0109] Where L is a lower triangular matrix, α vec Y is the first intermediate vector. vec Let β be a column vector consisting of all valid observations. vec K is the second intermediate vector. star_vec This is the covariance vector between the missing values ​​to be filled and each valid observation.

[0110] The predicted mean is obtained by multiplying the first intermediate vector and the second intermediate vector, and is used as the trend reference value at the missing values.

[0111] Specifically, the system extracts the independently solved first and second intermediate vectors and performs a multiplication calculation based on the vector inner product. The resulting scalar result is equivalent to the standard predicted mean of Gaussian process regression, and the system directly outputs it as a local trend reference value for the current missing value location. The specific numerical calculation relationship is expressed as follows:

[0112] T i =β vec T *α vec ;

[0113] Among them, T i β is the calculated predicted mean, which is the trend reference value at the missing values. vec T Let α be the transpose of the second intermediate vector. vec This is the first intermediate vector.

[0114] Alternatively, when dealing with specific scenarios where computing hardware resources are extremely abundant, the covariance matrix has a small dimension, and the internal condition number of the matrix is ​​well-performed, as an alternative to matrix factorization, the system can skip the decomposition step, directly calculate the inverse of the covariance matrix, and multiply it sequentially with the covariance vector and the effective observation vector to directly output the predicted mean.

[0115] When the number of valid observations is less than the observation threshold, local interpolation is selected to fit the local trend, and the trend reference value at the missing value is calculated.

[0116] Alternatively, when the number of valid observations is less than the observation threshold, the valid observations are used to select local interpolation to fit the local trend, and the trend reference value at the missing value is calculated.

[0117] Specifically, when the number of valid observations is insufficient to support the Gaussian model, the system automatically cuts off the global computation link and triggers the local interpolation module. The local interpolation module abandons the dependence on distant time series data and only uses the very few valid observations closest to the missing location to infer the local geometric trend, thereby minimizing the risk of overfitting while ensuring the robustness of the algorithm.

[0118] In a preferred implementation, local interpolation is selected to fit the local trend, and trend reference values ​​at the missing values ​​are calculated, including:

[0119] Extract valid observation points within a preset range before and after the time node of the missing value to be filled, form a neighborhood point set, and count the number of valid neighborhood points in the neighborhood point set.

[0120] In this embodiment, the system uses the current missing value as the central pointer and scans bidirectionally along the timeline to the past and future, extracting known data within a preset range to form a computational set. The preset range can be set to two valid observation points before and after the missing position. The total number of valid observation points actually acquired in the neighborhood point set is then counted.

[0121] When the number of valid neighborhood points is greater than or equal to a preset threshold, a quadratic polynomial is used to fit the neighborhood point set to calculate the trend reference value at the missing value.

[0122] Specifically, the system performs a secondary judgment by comparing the number of valid neighborhood points obtained from statistics with a preset threshold. The preset threshold can be set to 4. When the number of valid data points collected in the neighborhood is sufficient to support the calculation of the parabolic shape, the system calls the least squares algorithm engine. The engine uses the collected time index and observation values ​​to construct an overdetermined system of equations, solves for the polynomial coefficients, and uses the quadratic polynomial curve to calculate the geometric ordinate of the missing positions as the trend reference value for the final output.

[0123] When the number of valid neighborhood points is less than the threshold, linear interpolation is used to calculate the trend reference value at the missing value in the neighborhood point set.

[0124] Specifically, when effective data points in the neighborhood are extremely scarce, the system adopts the most conservative geometric estimation strategy. It extracts the two nearest effective observation points located directly in front of and behind the missing location, performs interpolation calculations using the geometric principle that two points determine a straight line, and uses the projected value of the corresponding missing time node on the straight line as the trend reference value.

[0125] In some alternative implementations, when faced with extreme conditions of extremely sparse data, as an alternative to linear interpolation, the system can also call the nearest neighbor interpolation algorithm to calculate the trend reference value, that is, directly copy the nearest valid observation as the filling benchmark.

[0126] In one embodiment of this application, based on trend fusion weights, a preliminary smoothed sequence and a trend reference value are weighted and fused to generate an imputed value. The imputed value is then used to replace the missing values ​​in the one-dimensional time series to obtain a preliminary reconstructed sequence, including:

[0127] Determine whether there are valid observations on both sides of the time point where the missing value to be filled exists.

[0128] In this embodiment, after calculating the smoothed sequence and local trend reference values, the system enters the data reconstruction and stitching stage. The index node of the currently processed missing value on the time axis is extracted, and all historical records before and after that index node are scanned. If at least one valid observation can be retrieved in both the time periods before and after this time node, it is determined that valid observations exist on both sides. If the missing value is located at the beginning or end of the entire one-dimensional time series, resulting in no valid observation records before or after its time node, it is determined that a valid observation exists only on one side.

[0129] In some alternative implementations, the system can accelerate the above-mentioned determination process by constructing a global Boolean mask matrix, using the non-zero element index of the matrix to quickly locate the valid boundaries of the time series, avoiding repeated sequence scanning operations when processing missing values ​​one by one.

[0130] When there are valid observations on both sides of the time node, the smoothed value and the trend reference value at the corresponding position in the preliminary smoothed sequence are weighted and summed based on the trend fusion weight to generate the fill value.

[0131] In this embodiment, for missing positions located within the time series, the system employs a trend soft-constraint fusion mechanism. Specifically, the preliminary smoothed sequence reflects the global structural characteristics of vegetation growth, while the trend reference value mainly reflects the local change trend around the missing period. The system uses trend fusion weights dynamically calculated based on observation density to perform linear fusion calculations on these two different dimensions of feature data. The specific numerical calculation relationship is expressed as follows:

[0132] z fill =(1-α)*z smooth +α*T ref ;

[0133] Among them, z fill The generated imputation value, α is the trend fusion weight, z smooth T represents the smoothed value at the corresponding time node in the initial smoothing sequence. ref These are the trend reference values ​​calculated at the corresponding time points.

[0134] By employing a weighted summation mechanism, the overall trend is preserved while the weighting coefficients guide missing values ​​to reasonably approximate local trends. This soft constraint strategy effectively avoids the problem of abnormally flattened peaks caused by directly using trend reference values ​​for hard replacement, ensuring the smoothness of data transition.

[0135] When there are valid observations on only one side of a time node, the standard deviation of all valid observations in the one-dimensional time series is calculated. A symmetric constraint interval is constructed by combining the preset scaling factor with the valid observation closest to the current missing value. The smoothed value at the corresponding position in the preliminary smoothed sequence is then restricted to the symmetric constraint interval to generate the imputed value.

[0136] Specifically, when missing values ​​are located at the beginning or end of the sequence and a trend cannot be constructed on both sides, any complex fitting algorithm is prone to generating divergent and outlier extrapolated values ​​due to the lack of constraint from one side. Therefore, the system switches to a one-sided symmetric interval constraint logic. The standard deviation of all valid observations is calculated by traversing the entire time series. The standard deviation reflects the overall discrete fluctuation level of the current pixel throughout the entire observation period. A preset scaling factor is introduced to calculate the half-width parameter of the constraint interval. The specific numerical calculation relationship is expressed as follows:

[0137] δ=u*σ;

[0138] Where δ is the half-width parameter of the constraint interval, u is the preset scaling factor, and σ is the standard deviation of all valid observations. The preset scaling factor is preferably set to 0.15, which is used to control the physical width of the constraint interval.

[0139] The nearest valid observation to the current missing position is extracted as the baseline anchor point. Symmetrical constraint intervals are constructed upwards and downwards using the half-width parameter, and a truncation function is called to limit the smoothed values. The specific numerical calculation relationship is expressed as follows:

[0140] z fill =clip(z smooth v obs -δ, v obs +δ);

[0141] Among them, z fill The generated padding value, clip is the truncation function, z smooth v is the smoothing value corresponding to the initial smoothing sequence. obs The nearest valid observation to the current missing value is represented by δ, which is the half-width parameter of the constraint interval. If the smoothed value exceeds the symmetric constraint interval, the truncation function will force the boundary extreme values ​​of the output interval to be used as filler values, thus cutting off the risk of divergence of endpoint data through hard constraints.

[0142] In another specific implementation, if the time series spans several years, the data source used to calculate the standard deviation can be replaced by all valid observations with locally valid observations from the year in which the current missing value is located, in order to eliminate the interference of long-term climate drift on the data standard deviation.

[0143] The generated imputation values ​​are used to replace the corresponding missing values, keeping the original valid observations in the one-dimensional time series unchanged, to obtain a preliminary reconstructed sequence.

[0144] In this embodiment, after completing the numerical calculations for all missing positions, the system executes a data fidelity strategy that only fills in the missing values. The filled values ​​are written one by one into the invalid placeholders of the original one-dimensional time series to construct the reconstructed sequence. Simultaneously, the existing valid observations in the original sequence are locked, and no smoothing filtering or numerical correction is applied to them.

[0145] The specific replacement logic is as follows: if the original input is a valid observation, the output maintains the original valid observation; if the original input is a missing value, the output is the generated imputation value. This underlying mechanism ensures the authenticity of the remote sensing inversion data and avoids secondary contamination of high-quality observation signals by traditional full-sequence smoothing methods. The one-dimensional array generated by splicing is the preliminary reconstructed sequence, awaiting further verification by physical rules.

[0146] According to one aspect of this application, multi-level physical constraint verification is applied to the imputation values ​​in the preliminary reconstructed sequence, including:

[0147] The first-level range constraint verification is performed on the fill value using the preset physical value range of the remote sensing vegetation index, and the fill value that exceeds the physical value range of the remote sensing vegetation index is restricted to the boundary extreme value.

[0148] Specifically, after obtaining the imputed values, due to the purely mathematical nature of various fitting extrapolation algorithms, the calculation results may exceed the domain of a specific physical quantity in the real world. Therefore, the system performs the most basic official physical range truncation. Taking the enhanced vegetation index acquired by the medium-resolution imaging spectrometer as an example, its officially recommended effective value lower limit is -0.2 and upper limit is 1.0. The system uses this range as the preset physical value range for the remote sensing vegetation index and calls the truncation algorithm to filter all generated imputed values. The specific numerical calculation relationship is expressed as follows:

[0149] z clip1 =clip(z fill (, -0.2, 1.0);

[0150] Among them, z clip1 The values ​​are the fill values ​​after the first-level range constraint verification, where clip is the truncation constraint function, and z is the fill value. fill The input filler values ​​are -0.2 and 1.0, which are the physical lower limit and physical upper limit of the enhanced vegetation index, respectively.

[0151] Extract the maximum and minimum values ​​from all valid observations, and calculate the range between the maximum and minimum values.

[0152] In this embodiment, after completing the basic physical boundary verification, due to the significant differences in the fluctuation range of vegetation indices corresponding to different land cover types, the system needs to formulate personalized data statistical boundaries for each individual pixel. The system extracts the set of valid observations from the entire one-dimensional time series being processed, iterates through the set of valid observations, and extracts the global maximum and global minimum values. The system performs a subtraction calculation, subtracting the minimum value from the obtained maximum value to obtain the range parameter representing the historical fluctuation span of that pixel.

[0153] Using the range and preset boundary coefficients, the maximum and minimum values ​​are calculated to obtain the upper limit and lower limit of the dynamic amplitude.

[0154] Specifically, considering that the missing period corresponding to the filled value may represent the extreme growth condition of that pixel in a specific year, directly using the historical maximum and minimum values ​​as absolute hard boundaries would obliterate normal interannual fluctuations. Therefore, the system introduces preset boundary coefficients to perform outward expansion calculations. The preset boundary coefficients include a first small positive number for the lower boundary and a second small positive number for the upper boundary. The first small positive number can be set to 0.08, and the second small positive number can be set to 0.1. The system multiplies the range by the two small positive numbers respectively, and adds or subtracts the product results to the extreme values ​​respectively, thereby constructing a shrinking container with a certain fault tolerance margin above the maximum value and below the minimum value.

[0155] By using the dynamic amplitude upper limit and dynamic amplitude lower limit, the fill value after the first level of range constraint verification is clipped and restricted to complete the data-driven amplitude constraint of the fill value.

[0156] In this embodiment, the system uses a dynamic amplitude upper limit and a dynamic amplitude lower limit to perform secondary filtering and interception on the filler values ​​that have passed the first level of constraints. To ensure that the final output does not violate physical common sense, the system further stipulates that the final effective dynamic lower limit must be the larger value between the physical lower limit and the calculated lower limit, and the final effective dynamic upper limit must be the smaller value between the physical upper limit and the calculated upper limit. The specific numerical calculation relationship is expressed as follows:

[0157] z final =clip(z clip1 y min -c1*R,y max +c2*R);

[0158] Among them, z final To complete the fill value after data-driven amplitude constraints, clip is the truncation constraint function, z clip1 The filled value after the first-level range constraint verification, y min The minimum value of the effective observations, y max R is the maximum effective observation value, C1 is the range parameter, and C2 is the preset boundary coefficient for controlling the expansion of the lower and upper limits, respectively.

[0159] In a further embodiment, after completing the data-driven amplitude constraint on the filled value, applying multi-level physical constraint verification also includes performing local peak protection and repair on the filled value, specifically including:

[0160] For each valid observation in a one-dimensional time series, determine whether the current valid observation is a local maximum within a preset local window, and determine whether the local maximum is greater than a preset peak threshold.

[0161] Specifically, to prevent the smoothing model from excessively smoothing out the peak values ​​of real alpine vegetation or crops during their peak growing season, the system adds a peak protection mechanism in the final stage of constraint processing. Using each valid observation as a central pointer, a preset local window is expanded to both sides. The radius of the preset local window can be set to two time steps. The system determines whether the observation pointed to by the pointer is a local maximum within the window; if so, it further determines whether the local maximum exceeds a preset peak threshold. The preset peak threshold can be set to 0.25 to filter out minor fluctuation noise in bare land or dormant backgrounds.

[0162] When the above judgment conditions are met, further verification is made to determine whether the local maximum value satisfies the mountain morphology characteristics that the effective observation value adjacent to the left is less than or equal to the local maximum value and the effective observation value adjacent to the right is less than the local maximum value.

[0163] In this embodiment, when a value is sufficiently large locally, the system further verifies from a geometric morphological perspective whether it constitutes a real isolated peak. The system extracts the first valid observation value adjacent to the left and the first valid observation value adjacent to the right of the local maximum value. The system performs a logical judgment operation to verify whether the value on the left is less than or equal to the local maximum value and whether the value on the right is less than the local maximum value. If the combined conditional logic is satisfied, it is confirmed that there is a typical peak shape with high confidence at that location, and the local maximum point is qualified for targeted protection by the system.

[0164] If the mountain peak morphology characteristics are met, the data-driven amplitude constraint-bound fill value is extracted from the neighborhood of the local maximum value, and it is determined whether the extracted fill value is lower than the product of the local maximum value and the preset judgment coefficient.

[0165] Specifically, the system retrieves and locates the imputed values ​​near the center point of the typical mountain peak shape due to missing data. The system extracts these imputed values ​​and introduces a preset judgment coefficient to quantify the severity of the algorithm's flattening of the imputed values. The preset judgment coefficient can be set to 0.75. The system multiplies the maximum center value with the preset judgment coefficient. If the extracted imputed value is lower than the product, the system determines that the original algorithm has caused significant and unacceptable signal attenuation in the temporal reconstruction of this area, and then triggers a mild repair module.

[0166] If so, the filled value is linearly weighted and repaired by combining the preset repair weight coefficient and the preset peak recovery ratio to compensate for the real peak features that are weakened during the smoothing process.

[0167] In this embodiment, the system invokes a linear interpolation engine to inject the original true peak signal component into the over-attenuated fill value. The specific numerical calculation relationship is expressed as follows:

[0168] z corrected =(1-η)*z final +η*p*r;

[0169] Among them, z corrected The final fill value after repair is achieved, η is the preset repair weight coefficient, and z final The current fill value is after data-driven amplitude constraint, p is the maximum value of the effective observation points that are judged to be typical mountain peak shape, and r is the preset peak recovery ratio.

[0170] For ease of understanding, a normalized example without a specific physical context is provided. Assume that, after morphological verification, the effective observation peak p within the local window is 0.8, and the filled value z after range constraints... final The value is 0.5. The system calculates the judgment standard as 0.8 * 0.75 = 0.6 based on the preset judgment coefficient of 0.75. Since 0.5 is lower than 0.6, it meets the conditions for triggering repair. The system extracts the preset repair weight coefficient η as 0.3 and the preset peak recovery ratio r as 0.88. Substituting the above parameters into the repair calculation relationship, the calculation process is as follows: z corrected =(1-0.3)*0.5+0.3*0.8*0.88=0.35+0.2112=0.5612. It can be seen that the system, through gentle repair calculations without changing the original valid observations and with a small impact, stably raises the data points that were originally excessively flattened to 0.5 to 0.5612. This approach approximates the true phenological peak while avoiding the overcompensation effect caused by artificially introduced spurious fluctuations. The repaired and filled values ​​are finally merged with other data to output a data file set for downstream tasks.

[0171] In a detailed embodiment, to verify the objective effectiveness of the multi-level constraint and branch fusion algorithm, the system performed an artificial missing data simulation experiment. Specifically, time-series data of the medium-resolution imaging spectrometer-enhanced vegetation index for a specific year were acquired as the experimental baseline dataset, for example, data from January 1, 2024 to December 31, 2024. An image dataset missing 30% of pixel values ​​was artificially constructed using a random masking missing data method to simulate cloud cover pollution in a real-world environment. As a comparative test case, data missing on specific dates was simulated, such as image data from May 8, 2024 and August 28, 2024.

[0172] The system invokes the methods provided in this application, traditional local smoothing filtering methods, global regularization smoothing methods, and basic geometric interpolation methods to perform parallel reconstruction calculations on the aforementioned mask dataset. For specific comparison, traditional methods include Savitzky-Golay filtering, traditional Whittaker smoothing, and linear interpolation. After calculation, the system extracts the reconstructed values ​​at each missing location and performs error statistical comparison with the original ground truth images without missing data for the corresponding dates.

[0173] The system calculates multiple statistical indicators to quantify reconstruction accuracy, including root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (CCD). Statistical comparison results show that the method provided in this application achieves optimal performance across all error metrics. Taking the mask reconstruction on May 8, 2024, as an example, the method in this application calculates an RMSE of 0.0624, an MAE of 0.0476, and a CCD of 0.7832. Comparative calculations show that traditional Whittaker smoothing, under the same conditions, has an RMSE as high as 0.0973 and a CCD of only 0.4720. It can be seen that the method in this application, by introducing an observation density-driven parameter adjustment engine and multi-level physical morphology constraints, avoids the local divergence and peak flattening caused by fixed empirical parameters in traditional methods, achieving a unity of data authenticity, temporal continuity, and physical rationality.

[0174] like Figure 5 As shown, according to another aspect of this application, a method for reconstructing time-series missing values ​​of remote sensing vegetation indices using adaptive smoothing and trend constraint synergy, taking the MODIS EVI time series as an example, mainly includes the following steps:

[0175] Step 1: Perform remote sensing time series data preprocessing.

[0176] MODIS EVI (Enhanced Vegetation Index of Medium Resolution Imaging Spectroradiometer) product data were acquired from the remote sensing cloud computing platform GEE. Low-quality pixels, such as those affected by cloud contamination, were identified and removed using the quality control (QA) band, resulting in high-quality Enhanced Vegetation Index (EVI) time series data for the study area. Based on this, a time series data cube was constructed: Stack(t, x, y); where t represents the time dimension and (x, y) represents the spatial pixel coordinates. Furthermore, invalid values ​​in the EVI time series were uniformly converted to missing values ​​(NaN) for subsequent standardized modeling and processing. This standardized organization of time series remote sensing data provides a foundation for pixel-by-pixel time series analysis and effectively avoids interference from invalid data in subsequent calculations.

[0177] Step 2: Perform effective pixel identification and time series construction.

[0178] An effective pixel mask is constructed, processing only pixel sequences with valid observations to avoid invalid computations on pixel sequences with completely missing values. The time series for each effective pixel (x, y) can be represented as:

[0179] Y = [y1, y2, ..., y] T＇ ];

[0180] Among them, y i Let represent the EVI value at time i. If there is no valid observation at that time, it is denoted as NaN.

[0181] Transforming the spatial problem into a pixel-by-pixel one-dimensional time series reconstruction problem provides a foundation for subsequent detailed modeling while ensuring computational efficiency.

[0182] Step 3: Perform adaptive weighted Whittaker smoothing modeling.

[0183] To achieve smooth reconstruction of the remotely sensed vegetation index time series while ensuring the accuracy of valid observations and suppressing noise interference, an adaptive weighted Whittaker smoothing model is constructed to fit the pixel-by-pixel time series. Specifically, based on a pixel's time series Y, the reconstructed sequence is calculated as Z = [z1, z2, ..., z...]. T＇ This allows it to maintain good smoothness while preserving the original observation characteristics.

[0184] (1) Construct the optimization objective function.

[0185] The objective function is as follows:

[0186] ;

[0187] Among them, y i The original observation value; z i The smoothed time series values; w i These are weighting coefficients used to adjust the contribution of different observations to the fitting process; when y i When y is a valid observation, it is assigned a large value, such as a constant of 10, to strengthen its constraint on the objective function; when y i When a value is missing, the value is set to 0, so that the corresponding item does not participate in the data consistency constraint; λ is the smoothing parameter; Δ d This represents the d-th order difference operator, typically with a value of 2; T' is the length of the time series. By increasing the weight of the effective observations, the fidelity of the smoothing curve to the original observations can be enhanced.

[0188] (2) Construct the matrix form of the difference operator.

[0189] Represent the time series as a vector form: Z = [z1, z2, ..., z2] T＇ ]T ,in T To transpose; construct a first-order difference matrix D1∈R (T＇-1)×T＇ :

[0190] ;

[0191] Based on this, a second-order difference matrix D2∈R is further constructed. (T＇-2)×T＇ This is used to constrain the curvature changes of a time series. Each row of the second-order difference matrix corresponds to a curvature term, and its specific form is as follows:

[0192] .

[0193] (3) Perform matrix form transformation.

[0194] The above optimization problem can be transformed into a standard quadratic form:

[0195] ;

[0196] Where W=diag(w i )∈R T＇×T＇ Let W be the diagonal weight matrix; taking the derivative of the objective function and setting it to zero, we obtain the following system of linear equations: (W+λD2) T D2)Z=WY, the coefficient matrix corresponding to this system of equations is A=(W+λD2) T D2) is a banded sparse symmetric matrix that exhibits symmetric positive definite properties under the condition that the weights of the effective observations are sufficiently constrained. Thus, the optimization problem is transformed into solving a system of sparse symmetric positive definite linear equations.

[0197] (4) Perform numerical solutions.

[0198] Considering that remote sensing time series data typically have large spatial dimensions, a sparse matrix solution strategy is adopted to improve computational efficiency: coefficient matrix A = W + λD² T D2 is a banded sparse symmetric matrix; it is preferable to use a direct solution method based on sparse matrices, which ensures numerical stability while taking into account computational efficiency.

[0199] Step 4: Perform adaptive smoothing parameter adjustment based on observation density.

[0200] To address the issue that traditional Whittaker smoothing has fixed smoothing parameters and is difficult to adapt to different data qualities, an adaptive adjustment mechanism based on observation density is proposed. Furthermore, trend fusion weights are also incorporated into the unified driving framework to achieve coordinated adjustment of smoothing intensity and trend dependence.

[0201] (1) Observation density.

[0202] For a time series of length T', the observation density is defined as:

[0203] ρ = N / T';

[0204] Where N is the number of valid observations in the time series, and T' represents the total length of the time series. The observation density ρ reflects the completeness of the time series data.

[0205] (2) Smoothing parameter.

[0206] Based on the observation density, a continuous smoothing parameter function is constructed:

[0207] λ=λ0ρ γ ;

[0208] Where λ0 is the basic smoothing parameter and γ∈[0,1] is the adjustment exponent. When ρ→0 (sparse observations), λ→0, thus weakening the smoothing constraint and avoiding oversmoothing in the case of insufficient information; when ρ→1 (dense observations), λ→λ0, which enhances the smoothing ability and improves noise resistance.

[0209] (3) Trend fusion weight.

[0210] Furthermore, the trend fusion weight α is also constructed as a function of the observation density:

[0211] α=α0·ρ q ;

[0212] Where α0∈[0,1] is the basic trend fusion weight; q is the adjustment index, used to control the degree of influence of observation density on trend weight; and ρ is the observation density. When the observation density is high, α takes a larger value, making the reconstruction result more dependent on the local trend; when the observation density is low, α takes a smaller value, making the interpolation value more dependent on the smoothing result.

[0213] Step 5: Perform missing value imputation and trend soft constraint.

[0214] After obtaining the smooth sequence Z, optimization and reconstruction are performed by filling in only missing values ​​and using soft trend constraints. Specifically, when there are valid observations on both sides of the missing value position, the following procedures (2) and (3) are executed; when there are only one-sided observations at the missing value position, the following procedure (4) is executed; when the missing value is located at the beginning or end of the sequence and a trend cannot be constructed, z is directly used. i .

[0215] (1) Initial filling.

[0216] Retaining the original observations, the missing values ​​are filled in only with smoothed values ​​to construct the reconstructed sequence z. i ':

[0217] z i '={yi If y i For valid observations; z i If y i For missing values};

[0218] Among them, y i The original observation value, z i These are the smoothed time series values.

[0219] (2) Construct trend reference values.

[0220] For a missing position i, if there are valid observations before and after it, an adaptive trend estimation strategy is adopted, and different trend reference value generation methods are selected according to the number of valid observations N.

[0221] 1) Gaussian process regression (GPR) fitting.

[0222] When the number of valid observations N≥N min This fitting method is used when N is selected. min =8. Let the set of effective observation time points be {t1, t2, ..., t...} N The corresponding set of observations is {y1, y2, ..., y}. N A zero-mean Gaussian process model is adopted, and the covariance function uses the squared exponential kernel (RBF kernel):

[0223] k(t, t') = exp(-(t-t') 2 / (2ζ 2 ))+σ n 2 ·δ t，t' ;

[0224] Where ζ is the length scale, and in this embodiment, ζ=3, corresponding to approximately 48 days; σ n 2 For the noise variance, σ is taken in this embodiment. n 2 =0.01; δ t，t' Let be the Kronecker function; t and t' are two different time points.

[0225] Construct the covariance matrix:

[0226] ;

[0227] For the time index i of the missing point, calculate its covariance vector with that of the valid observation points:

[0228] k * =[k(i,t1) k(i,t2) … k(i,t N )] T ;

[0229] Then the predicted mean of GPR, T i (i.e., trend reference value) is:

[0230] T i =k * T K -1 y;

[0231] Where K is the covariance matrix, k(t) N , t N Let k be the covariance between the Nth time point and itself. * Let y = [y1, y2, ..., y3] be the covariance vector between missing points and valid observations. N ] T In practical solutions, to avoid direct inversion, the Cholesky decomposition is used: Calculate LL... T =K, solve for Lα=y and Lβ=k * Where α and β are intermediate vectors, then T i =β T α.

[0232] 2) Local interpolation.

[0233] When the number of valid observations N < N min This fitting method is used. P valid observation points are taken before and after the missing point, preferably p=2, to form a neighborhood point set {t_left1, ..., t_left}. p ,t_right1,…,t_right p} and its corresponding value.

[0234] If the total number of valid neighborhood points is M (2≤M≤2p), then the following steps are taken:

[0235] If M≥4, a quadratic polynomial fitting is used to solve for the coefficients a, b, and c, minimizing the following expression:

[0236] ;

[0237] Then the trend reference value T i =a i 2 +b i +c;

[0238] If M < 4, use linear interpolation:

[0239] T i =y b +((y a -y b ) / (ab))·(ib);

[0240] Wherein, the previous observation point is (b, y) b ), and the subsequent observation point is (a, y a In particular, when a=b, T i =y b ;y a The vegetation index observation value is the closest valid observation point preceding the missing point; y b is the vegetation index observation value of the nearest valid observation point after the missing point; i is the time index of the missing position that needs to be filled.

[0241] Through an adaptive strategy, when there are sufficient effective observations, GPR can use global information from all effective points to provide a more reasonable trend extrapolation; when there are insufficient observations, local interpolation avoids the risk of overfitting, thereby significantly improving the reconstruction accuracy of sparse periods while ensuring robustness.

[0242] (3) Trend soft constraint fusion mechanism.

[0243] A weighted fusion method is adopted:

[0244] z i ''=(1-α)·z i '+α·T i ;

[0245] Among them, z i '' represents the final reconstructed value after trend soft constraint fusion; trend fusion weight α = α0·ρ q Trend reference value T i It mainly reflects the local trend of change, while the smoothing result z i It reflects the global structural characteristics. Through a weighted fusion mechanism, while preserving the smooth overall trend, it guides missing values ​​to reasonably approximate local change trends, avoiding the peak flattening problem caused by hard constraints.

[0246] (4) One-sided symmetric interval constraint.

[0247] When the missing location has only one valid observation (either front or back only), a symmetric interval constraint is used:

[0248] z i ''=clip(z i ', v obs -δ, v obs +δ);

[0249] Among them, v obs This represents the most recent valid observation; δ = u·σ, where σ represents the standard deviation of the valid observations in the time series; u is the proportionality coefficient, with a default value of 0.15, used to control the width of the constraint interval.

[0250] Step 6: Apply multi-level physical constraints.

[0251] To further enhance the physical plausibility and stability of the results, multi-level constraints were applied to the missing values ​​filled in the reconstructed time series, while the original valid observations remained unchanged. The multi-level constraints followed a progressive principle, specifically in the following order: first, range constraints were applied, forcibly limiting the filled values ​​to the official physical range of the vegetation index; second, data-driven amplitude constraints were applied, further shrinking outliers using the statistical distribution of valid observations; and third, peak protection was implemented, gently correcting peak reduction issues in typical mountain morphologies. Specifically, this included:

[0252] (1) Value range constraints.

[0253] According to the official MODIS EVI value range, which limits EVI to [-0.2, 1.0], the constraint can be expressed as: z i ''' = clip(z i '', -0.2, 1.0).

[0254] (2) Data-driven amplitude constraints.

[0255] Let the set of valid observations be: Y valid Its extreme value is y min =min(Y valid ), y max =max(Y valid If the range R* = y, then the range R* = y. max -y min Regarding the reconstruction result z i Apply constraints:

[0256] z final =clip(z i ''' ,y min -c1·R*,y max +c2·R*);

[0257] Where c1 and c2 are small positive numbers, and in this embodiment, they are taken as 0.08 and 0.1, respectively. y is also specified... min The intersection of -c1·R* and [-0.2, 1.0] is the final lower bound, which is max(y). min -c1·R*,-0.2), with an upper limit of min(y) max +c2·R*,1).

[0258] (3) Peak protection.

[0259] To further prevent peak clipping caused by smoothing, a peak protection mechanism is added after physical constraints. For each valid observation point, if its value is the maximum value within a local window (radius 2) and greater than a threshold (e.g., 0.25), and satisfies the condition that the valid value on the left ≤ peak value > valid value on the right (i.e., typical peak shape), then it is checked whether the filling value of the missing points in its neighborhood is significantly lower than the peak value. If it is lower than 0.75 times the peak value, then gentle repair is performed.

[0260] z corrected =(1-η)·z final +η·p·r;

[0261] Where η is the repair weight coefficient, preferably 0.3; r is the peak recovery ratio, preferably 0.88; p is the peak value of the effective observation point, z final This is the current filled value after applying range constraints and data-driven amplitude constraints. This repair only applies to the originally missing points, and the repair magnitude is small, preserving the true peak value while avoiding the artificial introduction of false fluctuations.

[0262] It should be noted that steps 5 and 6 are executed sequentially. Step 5 can integrate global smoothing structures and local trend features, while step 6, as a post-processing step, only verifies the physical rationality and statistical consistency of the filled values ​​and does not interfere with the trend fusion process in reverse. A single forward processing strategy is used, without iterative updates.

[0263] Step 7: Output Results and Verify Accuracy.

[0264] Finally, the reconstructed result is restored to a time-series raster image and output in formats such as GeoTIFF. To verify the effectiveness of this embodiment, an artificial missing pixel simulation experiment is used for evaluation. MODIS EVI time-series data from 2024, covering the period from January 1, 2024 to December 31, 2024, is selected as the experimental data. Two EVI image data sets with 30% missing pixel values ​​are constructed using a random masking missing pixel method, corresponding to May 8, 2024 and August 28, 2024, respectively. The original EVI image corresponding to the corresponding date is used as the reference ground truth. Based on this, the method proposed in this embodiment and the comparison method (Savitzky-Golay filtering, linear interpolation, Whittaker smoothing) are used for reconstruction. The error between the reconstructed result and the true value is calculated at the missing location. Evaluation indicators include root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²). 2 The specific results are shown below:

[0265] In the method proposed in this embodiment, the RMSE of image MOD13Q1_EVI_20240508 is 0.0624, the MAE is 0.0476, and the R...2 The RMSE is 0.7832; in the method proposed in this embodiment, the RMSE of image MOD13Q1_EVI_20240828 is 0.0762, the MAE is 0.0592, and the R... 2 The RMSE is 0.6547; in the Savitzky-Golay filtering method, the RMSE of image MOD13Q1_EVI_20240508 is 0.0636, MAE is 0.0486, and R... 2 The RMSE is 0.7687; in the Savitzky-Golay filtering method, the RMSE of image MOD13Q1_EVI_20240828 is 0.0780, MAE is 0.0605, and R... 2 The RMSE is 0.6240; in the traditional Whittaker method, the RMSE of image MOD13Q1_EVI_20240508 is 0.0973, MAE is 0.0765, and R... 2 The RMSE is 0.4720; in the traditional Whittaker method, the RMSE of image MOD13Q1_EVI_20240828 is 0.1515, MAE is 0.1289, and R... 2 The value is -0.3648; in the linear interpolation method, the RMSE of image MOD13Q1_EVI_20240508 is 0.0635, the MAE is 0.0487, and the R... 2 The RMSE is 0.7764; in the linear interpolation method, the RMSE of image MOD13Q1_EVI_20240828 is 0.0826, the MAE is 0.0647, and the R... 2 It is 0.5982.

[0266] Experimental results show that the RMSE, MAE, and R of the method proposed in this embodiment are satisfactory. 2 The experimental results achieved values ​​of 0.0624, 0.0476, and 0.7832, respectively, demonstrating good reconstruction accuracy and stability while effectively maintaining the temporal series trend. An example of the time series of remote sensing EVI reconstruction results at the pixel scale proposed in this embodiment is shown below. Figure 6 As shown.

[0267] This embodiment constructs an adaptive reconstruction mechanism and a multi-constraint collaborative optimization framework for remote sensing vegetation index time series. On the one hand, by introducing a data fidelity strategy that only fills in missing values ​​while keeping the original observations unchanged, the authenticity of the reconstruction results is guaranteed from the source. At the same time, an adaptive function of smoothing parameter λ and trend fusion weight α is constructed based on the effective observation density of the time series, realizing the collaborative dynamic matching of smoothing intensity and trend dependence, breaking through the limitation of traditional methods that rely on fixed parameters. On the other hand, for missing locations, the trend estimation method is adaptively selected according to the number of effective observations: when observations are sufficient, Gaussian process regression (GPR) is used to fit the global trend, making full use of the overall temporal structure; when observations are insufficient, local quadratic or linear interpolation is used to take into account local variation characteristics, effectively improving the reconstruction accuracy of sparse periods. On this basis, the smoothing results and trend reference values ​​are weighted and fused, with the smoothing results carrying the global structure and the trend reference values ​​guiding local changes, which avoids peak flattening caused by hard constraints and maintains the continuity of temporal changes. By superimposing value range constraints and data-driven amplitude constraints, and adding a peak protection mechanism, the peak-shaving problem of typical mountain peak morphology is gently repaired to ensure that the reconstruction result is physically reasonable and consistent with the statistical characteristics of the original data. Through synergy, this method systematically solves the problems of excessive smoothing, insufficient adaptability, lack of physical constraints, and weakening of key phenological features in existing methods, achieving a unity of data authenticity, temporal continuity, and physical reasonableness.

[0268] This invention addresses the specific shortcomings of fixed parameters leading to excessive smoothing or extrapolation failure by introducing a dual-parameter collaborative adjustment mechanism based on observation density. This mechanism dynamically drives smoothing strength and trend dependence through data completeness, reducing constraints to prevent divergence during reconstruction when data is sparse and enhancing filtering to resist noise during dense data, thus improving the adaptive processing capability for complex cloud and rain interference conditions. A data-driven conditional branch trend estimation and soft-constraint fusion framework is constructed. While preserving the overall temporal structure of the smooth curve, it guides missing values ​​to smoothly approximate the local true change trend, effectively eliminating the step-like artifacts caused by rigid stitching. By superimposing multi-level physical range constraints and a gentle restoration mechanism for typical mountain peak morphology, key phenological extreme signals attenuated by the algorithm are identified and quantitatively compensated. In summary, this invention not only ensures the fidelity of the original valid observation signals from the underlying mechanism but also reduces the average error of temporal filling in large-scale missing scenarios from the high level of traditional empirical algorithms, achieving a unity of spatiotemporal continuity, physical authenticity, and statistical rationality in remote sensing monitoring data.

[0269] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for reconstructing time-series missing values ​​of remote sensing vegetation indices using adaptive smoothing and trend constraint synergy, characterized in that, include: Acquire preprocessed remote sensing vegetation index time series data, extract valid observations and identify missing values, and construct a one-dimensional time series; Calculate the observation density of a one-dimensional time series and determine the smoothing parameter and trend fusion weights accordingly; An adaptive weighted Whittaker smoothing model is constructed based on smoothing parameters, and a one-dimensional time series is fitted accordingly to obtain a preliminary smoothed sequence. Based on the number of valid observations in the one-dimensional time series, conditional branching is determined, and the corresponding trend estimation method is adaptively selected to calculate the trend reference value at the missing value. Based on trend fusion weights, the preliminary smoothed sequence and trend reference values ​​are weighted and fused to generate imputation values. The imputation values ​​are used to replace the missing values ​​in the one-dimensional time series to obtain the preliminary reconstructed sequence, while keeping the effective observation values ​​unchanged. Apply multi-level physical constraint verification to the imputation values ​​in the initial reconstructed sequence, and output the reconstruction results of the temporal missing values.

2. The method according to claim 1, characterized in that, The preliminary smoothed sequence is obtained, including: Different weighting coefficients are assigned to valid observations and missing values ​​in a one-dimensional time series, with the weighting coefficient for valid observations being greater than that for missing values. Construct a difference matrix to constrain temporal curvature variations; A smoothing objective function is constructed by combining weighting coefficients, difference matrices, and smoothing parameters, and then the smoothing objective function is transformed into a sparse symmetric positive definite linear equation system. A direct solution method is used to numerically solve the coefficient matrix of a sparse symmetric positive definite linear equation system, and a preliminary smooth sequence is output.

3. The method according to claim 1, characterized in that, Determine the smoothing parameters and trend fusion weights, including: Extract the total length of the one-dimensional time series, calculate the ratio of the number of valid observations to the total length, and use the ratio as the observation density; Construct a continuous adjustment function with observation density as the independent variable; The smoothing parameters are calculated using a continuous adjustment function combined with pre-configured basic smoothing parameters, and the trend fusion weights are calculated using a continuous adjustment function combined with pre-configured basic trend fusion weights.

4. The method according to claim 1, characterized in that, Adaptively select the corresponding trend estimation method to calculate the trend reference value at the missing value, including: The number of valid observations is compared with a preset observation threshold; When the number of valid observations is greater than or equal to the observation threshold, Gaussian process regression is selected to fit the global trend, and the trend reference value at the missing value is calculated. When the number of valid observations is less than the observation threshold, local interpolation is selected to fit the local trend, and the trend reference value at the missing value is calculated.

5. The method according to claim 4, characterized in that, Gaussian process regression is selected to fit the global trend, and trend reference values ​​at missing values ​​are calculated, including: Construct a covariance matrix that reflects the correlation between each valid observation, and calculate the covariance vector between the current missing value to be filled and each valid observation; The covariance matrix is ​​decomposed using matrix decomposition. The first intermediate vector is obtained by combining the decomposed covariance matrix and the effective observations, and the second intermediate vector is obtained by combining the decomposed covariance matrix and the covariance vector. The predicted mean is obtained by multiplying the first intermediate vector and the second intermediate vector, and is used as the trend reference value at the missing values.

6. The method according to claim 4, characterized in that, Select local interpolation to fit the local trend, and calculate the trend reference value at the missing value location, including: Extract valid observation points within a preset range before and after the time node of the missing value to be filled, form a neighborhood point set, and count the number of valid neighborhood points in the neighborhood point set; When the number of valid neighborhood points is greater than or equal to the preset threshold, a quadratic polynomial is used to fit the neighborhood point set to calculate the trend reference value at the missing value. When the number of valid neighborhood points is less than the threshold, linear interpolation is used to calculate the trend reference value at the missing value in the neighborhood point set.

7. The method according to claim 1, characterized in that, The preliminary reconstructed sequence was obtained, including: Determine whether there are valid observations on both sides of the time point where the missing value to be filled exists; When there are valid observations on both sides of the time node, the smoothed value and the trend reference value at the corresponding position in the preliminary smoothed sequence are weighted and summed based on the trend fusion weight to generate the fill value; When there are only valid observations on one side of a time node, the standard deviation of all valid observations in the one-dimensional time series is calculated. A symmetric constraint interval is constructed by combining the preset scaling factor with the valid observation closest to the current missing value. The smoothed value at the corresponding position in the preliminary smoothed sequence is restricted to the symmetric constraint interval to generate the imputed value. The generated imputation values ​​are used to replace the corresponding missing values, keeping the original valid observations in the one-dimensional time series unchanged, to obtain a preliminary reconstructed sequence.

8. The method according to claim 1, characterized in that, Apply multi-level physical constraint verification to the imputation values ​​in the initial reconstructed sequence, including: The first-level range constraint verification of the fill value is performed using the preset physical value range of the remote sensing vegetation index. Extract the maximum and minimum values ​​from all valid observations, and calculate the range between the maximum and minimum values; Using the range and preset boundary coefficients, the maximum and minimum values ​​are calculated by expansion to obtain the upper limit and lower limit of the dynamic amplitude. By using the dynamic amplitude upper limit and dynamic amplitude lower limit, the fill value after the first level of range constraint verification is clipped and restricted to complete the data-driven amplitude constraint of the fill value.

9. The method according to claim 3, characterized in that, Construct a continuous adjustment function with observation density as the independent variable, specifically as follows: By introducing a preset adjustment index, a power function calculation model is constructed with the observation density as the base and the adjustment index as the exponent. When calculating the smoothing parameters and trend fusion weights, different preset adjustment exponent constraint power function calculation models are used to control the degree of attenuation of the influence of observation density on smoothing intensity and local trend dependence.

10. The method according to claim 5, characterized in that, Construct a covariance matrix that reflects the correlation among the various valid observations, including: Introduce the squared exponential kernel function as the covariance function; By combining the preset length scale parameter and the preset noise variance parameter, the time distance decay relationship between any two effective observation points in the time node sequence is calculated using the squared exponential kernel function, and the matrix elements in the covariance matrix are obtained.

Images