Method for adaptively reconstructing NDVI missing value

Abstract

The invention discloses a method for adaptively reconstructing an NDVI (normalized difference vegetation index) missing value, which belongs to the technical field of data missing value reconstruction, and is characterized in that an NDVI time sequence is decomposed into trend, season and abnormal components through dynamic Kalman filtering to identify abnormal observation caused by factors such as cloud and aerosol, and reconstruction is completed by using a space-time low-rank structure. Compared with a traditional method, the method can generate a high-quality NDVI time sequence without external quality information, is suitable for observation products from different satellite sources, and further provides support for continuous monitoring and accurate description of long-term vegetation dynamics.

Description

Technical Field

[0001] This invention relates to the field of data missing value reconstruction technology, and in particular to an adaptive method for reconstructing NDVI missing values. Background Technology

[0002] Since the 1970s, optical remote sensing satellites have provided continuous vegetation observation data at a global scale, offering long-term evidence for revealing the response of ecosystems to global change. Among these, the Normalized Difference Vegetation Index (NDVI), calculated by normalizing the difference in reflectance between the near-infrared and red light bands, is one of the most commonly used indicators characterizing vegetation productivity, phenological rhythms, and ecological status. However, due to the influence of optical imaging mechanisms and atmospheric radiative transfer characteristics, NDVI observations are easily affected by atmospheric disturbances such as clouds and aerosols, often exhibiting low values ​​or abnormal fluctuations, leading to reduced continuity and stability of the time series. Therefore, constructing high-quality NDVI time series is fundamental for accurately monitoring vegetation dynamics. Existing NDVI reconstruction methods mainly fall into two categories: those based on temporal information (time domain) and those based on spatiotemporal information (spatial domain).

[0003] In daily practice, the existing technical solutions have been found to have the following problems: Time-domain-based methods are based on the assumption that NDVI changes smoothly over time. They process pixels individually and rely solely on time series data for reconstruction. Typical methods include those based on smooth curve fitting (such as Savitzky-Golay, asymmetric Gaussian, and bilogistic methods), regularization methods that suppress anomalous fluctuations by balancing smoothness and observation fidelity (such as Whittaker filtering and L1 trend filtering), and frequency domain analysis methods that utilize the high-frequency characteristics of noise (such as Fourier and wavelet transforms). Because they ignore spatial domain information, these methods perform poorly under continuous or large-scale cloud cover conditions.

[0004] In contrast, spatial domain-based methods utilize auxiliary observations of neighboring pixels to construct spatial correlations, thereby improving the robustness of reconstructions under atmospheric disturbance conditions. Examples include the spatiotemporal Savitzky-Golay method, which estimates pollution observations by fitting inter-pixel mapping relationships, and tensor methods, which construct NDVI as a three-dimensional tensor in spatial, temporal, and annual cycle directions and utilize its low-rank properties for completion. These methods can suppress outliers while preserving valid observations, resulting in more stable reconstructions. However, they typically rely on standardized data quality labels, thus primarily applicable to data products with quality information, such as MODIS. Although datasets like GIMMS-3g+ offer advantages such as long observation time spans and wide coverage, the lack of reliable quality labels makes existing quality-information-dependent methods difficult to apply directly, limiting their potential application in long-term ecological research.

[0005] Therefore, it is necessary to provide a new technical solution to solve the above problems. Summary of the Invention

[0006] To address the aforementioned technical issues, this application provides an adaptive method for reconstructing missing NDVI values. This method can generate high-quality NDVI time series without the need for external quality information and is applicable to observation products from different satellite sources, thereby providing support for continuous monitoring and accurate depiction of long-term vegetation dynamics.

[0007] An adaptive method for reconstructing missing NDVI values ​​includes: S1. Kalman filtering is used to dynamically fit the original NDVI time series. By using the residual characteristics between the filtered predicted values ​​and the actual observed values, abnormal observation points caused by cloud and aerosol interference are adaptively identified. S2. Construct a three-dimensional tensor from the NDVI observation data within the region, and rearrange the tensor dimensions based on spatial similarity; S3. The low-rank tensor completion algorithm is used to reconstruct missing observation points or identified outliers to obtain a continuous, smooth and physically consistent NDVI sequence.

[0008] Optionally, in step S1, the model expression for NDVI is: ; In the formula, This represents the NDVI observation value at time t; This represents the regression variable, consisting of the seasonal and trend terms of NDVI; Sensitivity parameters corresponding to the regression variables; Corresponding observation noise.

[0009] Optionally, the prior distribution of the sensitivity parameter of the regression variable is expressed as: ; In the formula, Represents the state transition matrix; express Sensitivity parameters of the regression variable at any given time; This indicates process noise.

[0010] Optionally, in step S1, given the prior distribution and NDVI observations, the posterior distribution of the sensitivity parameters of the optimal regression variable can be solved by Kalman filtering; The expression for the posterior distribution of the sensitivity parameter of the regression variable is as follows: ; in, Sensitivity parameters representing regression variables Given observation data The prior distribution is as follows; Indicates that given prior knowledge Distribution The distribution of .

[0011] Optionally, the sensitivity parameter of the regression variable Given observation data The prior distribution of the following is represented as a Gaussian distribution with mean A and variance A; The expression for the mean A is: ; The expression for variance A is: ; In the formula, Let A represent the mean. Indicates the variance B; From the previous moment The mean of the posterior distribution is given; From the previous moment The variance of the posterior distribution is given.

[0012] Optionally, given prior knowledge Distribution The distribution of B can be determined by a Gaussian distribution of mean B and variance B. The expression for the mean B is: ; The expression for variance B is: ; In the formula, Indicates the mean B; Let B represent the variance.

[0013] Based on the sensitivity parameters of the regression variables Prior distribution and The sensitivity parameters of the regression variables can be used. The posterior distribution is represented as a Gaussian distribution with mean C and variance C. The expression for the mean C is: ; The expression for variance C is: ; In the formula, This represents the mean C; Indicates the variance C; Indicate prior The difference between the distribution and the observed distribution ; This represents the Kalman gain matrix, which determines the weights of prior information. .

[0014] Optionally, the low-rank tensor completion algorithm is used to reconstruct missing observation points or identified outliers. The low-rank tensor reconstruction algorithm treats noisy observations in remote sensing as missing values ​​and uses the optimization criterion of the lowest tensor rank after reconstruction to reconstruct noisy observations using observation information with good observation quality. The noise observation model reconstructed using observation information of good observation quality is as follows: ; In the formula, This indicates a good observation in the original observations; This indicates the noise-free observation after reconstruction; Tensor The weights of the nth mode.

[0015] Compared with the prior art, this application has at least the following beneficial effects: The adaptive NDVI missing value reconstruction method proposed in this invention, based on Kalman filtering and low-rank tensor reconstruction, can automatically identify anomalous observations caused by clouds and aerosols without the need for external quality control information, and achieve high-precision reconstruction by utilizing the spatiotemporal low-rank characteristics. This method can significantly improve the continuity and spatial coverage of NDVI time series, providing more reliable observational data support for vegetation carbon cycle assessment, atmospheric CO2 inversion, and climate change research. Attached Figure Description

[0016] The following sections will describe some specific embodiments of the invention in detail by way of example and not limitation, with reference to the accompanying drawings. The same reference numerals in the drawings denote the same or similar parts or portions. Those skilled in the art should understand that these drawings are not necessarily drawn to scale. In the drawings: Figure 1 This is a schematic diagram of the overall process of the present invention; Figure 2 A diagram illustrating the process of constructing a three-dimensional tensor from NDVI observation data within a region and rearranging the tensor dimensions based on spatial similarity; Figure 3 The average error comparison chart shows the inversion effect of the adaptive reconstruction method for missing NDVI values ​​of the present invention compared with other traditional inversion methods. Figure 4 This is a comparison chart of the mean absolute error between the adaptive reconstruction method for missing NDVI values ​​of the present invention and other traditional inversion methods. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this application clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. Based on the embodiments in this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0018] like Figure 1 As shown, an adaptive method for reconstructing missing NDVI values ​​includes: S1. Kalman filtering is used to dynamically fit the original NDVI time series. By using the residual characteristics between the filtered predicted values ​​and the actual observed values, abnormal observation points caused by cloud and aerosol interference are adaptively identified.

[0019] S2. Construct a three-dimensional tensor from the NDVI observation data within the region, and rearrange the tensor dimensions based on spatial similarity.

[0020] S3. The low-rank tensor completion algorithm is used to reconstruct missing observation points or identified outliers to obtain a continuous, smooth and physically consistent NDVI sequence.

[0021] The adaptive reconstruction method for missing NDVI values ​​of the present invention can generate high-quality NDVI time series without the need for external quality information and is applicable to observation products from different satellite sources, thereby providing support for continuous monitoring and accurate depiction of long-term vegetation dynamics.

[0022] The method for adaptively reconstructing missing NDVI values ​​according to embodiments of the present invention includes the following steps: Step 1: Model the NDVI observations of each pixel in the GIMMS product using Kalman filtering, transforming them into observations composed of trend and seasonality. The model expression for NDVI is: (1) In the formula, This represents the NDVI observation value at time t; This represents the regression variable, consisting of the seasonal and trend terms of NDVI; Sensitivity parameters corresponding to the regression variables; Corresponding observation noise.

[0023] using the state transition matrix supply prior distribution Given a prior distribution and observations In this case, the optimal solution can be obtained through Kalman filtering. posterior distribution The formula is as follows: (2) in, express Given observation data The prior distribution under the given condition can be represented as following a mean of . The variance is The Gaussian distribution.

[0024] and This can be expressed as the following formula: ; (3) In the formula, From the previous moment The mean of the posterior distribution is given; From the previous moment The variance of the posterior distribution is given.

[0025] Indicates that given prior knowledge Under the distribution, The distribution can be represented by the mean being The variance is The Gaussian distribution determines this.

[0026] and This can be expressed as the following formula: ; (4) Therefore, according to Prior distribution and Can The posterior distribution is represented by the mean being The variance is The Gaussian distribution.

[0027] and This can be expressed as the following formula: ; (5) In the formula, Indicate prior The difference between the distribution and the observed distribution; This represents the Kalman gain matrix, which determines the weights of prior information. ; .

[0028] Each optimal posterior can be solved using formulas 1-5. .

[0029] For each optimal posterior obtained by solving the problem Substituting into Formula 1 yields the optimal posterior. .calculate and The difference between the two values ​​is considered to be a noisy observation if the difference is less than -0.05.

[0030] Step 2: Construct a three-dimensional tensor from the NDVI observation data within the acquisition area, and rearrange the NDVI observation sequences marked as noise or missing in Step 1 according to spatial similarity.

[0031] Specifically, such as Figure 2 As shown, by rearranging the NDVI observation sequences marked as noise or missing in step 1 according to spatial neighborhood, phenological cycle and long-term trend, the rearranged tensor can represent spatial similarity, seasonal phenological similarity and long-term trend similarity in three dimensions, so as to construct an NDVI three-dimensional tensor with higher spatiotemporal similarity.

[0032] Step 3: For the NDVI tensor obtained after rearrangement, the low-rank tensor reconstruction algorithm is used to complete the missing observations, so as to recover the NDVI values ​​weakened by atmospheric interference using spatiotemporal structure information and generate a continuous and physically consistent NDVI time series.

[0033] The low-rank tensor reconstruction algorithm treats noisy observations in remote sensing as missing values ​​and uses the lowest rank of the reconstructed tensor as the optimization criterion to reconstruct noisy observations using observation information of good observation quality. Its expression is as follows: (6) In the formula, This represents a good observation from the original observations. This indicates the noise-free observation after reconstruction. Tensor The weights of the nth mode.

[0034] In low-rank tensor reconstruction, this method employs the Alternating Direction Multiplier Method (ADMM) for optimization.

[0035] Specifically, the original optimization problem (Equation 6) is transformed into an equivalent form (Equation 7), and the optimal solution is gradually approximated through iterative loops.

[0036] (7) In the k-th iteration, intermediate variables The optimal estimate can be expressed as Equation 8.

[0037] (8)

[0038] in, For soft singular value thresholding, the truncation threshold is... This updated formula is the optimal solution to the optimization problem (Formula 7).

[0039] In obtaining This can then be used to estimate the optimal reconstruction tensor for the k-th iteration. (Formula 9).

[0040]

[0041] When iterating tensor The difference between this update and the previous one is below a set threshold. When the iteration stops, the tensor that has finally converged is... This serves as the NDVI tensor after missing value reconstruction. Subsequently, the tensor is reverse-rearranged according to the spatial structure before rearrangement to obtain high-quality NDVI observation data with noise interference removed.

[0042] Based on the above method, missing NDVI data generated under widespread cloud and fog conditions in southern China in 2020 were reconstructed. The reconstruction results are attached. Figure 3 As shown in the figure, compared to the abnormally low value areas in the original observations, the reconstructed NDVI exhibits a more continuous spatiotemporal distribution, and the values ​​better match the regional vegetation growth characteristics. The results of simulation experiments comparing the results of different reconstruction methods are shown in the figure. Figure 4 As shown, the method of this invention can effectively reconstruct missing NDVI values ​​even without quality indicators.

[0043] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments according to this application. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0044] It should be noted that the terms "first," "second," etc., used in the specification, claims, and accompanying drawings of this application 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 this application described herein can be implemented in sequences other than those illustrated or described herein.

[0045] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for adaptively reconstructing missing NDVI values, characterized in that, include: S1. Kalman filtering is used to dynamically fit the original NDVI time series. By using the residual characteristics between the filtered predicted values ​​and the actual observed values, abnormal observation points caused by cloud and aerosol interference are adaptively identified. S2. Construct a three-dimensional tensor from the NDVI observation data within the region, and rearrange the tensor dimensions based on spatial similarity; S3. The low-rank tensor completion algorithm is used to reconstruct missing observation points or identified outliers to obtain a continuous, smooth and physically consistent NDVI sequence.

2. The method for adaptively reconstructing missing NDVI values ​​as described in claim 1, characterized in that, In step S1, the model expression for NDVI is: ； In the formula, This represents the NDVI observation value at time t; This represents the regression variable, consisting of the seasonal and trend terms of NDVI; Sensitivity parameters corresponding to the regression variables; Corresponding observation noise.

3. The method for adaptively reconstructing missing NDVI values ​​as described in claim 2, characterized in that, The expression for the prior distribution of the sensitivity parameter of the regression variable is: ； In the formula, Represents the state transition matrix; express Sensitivity parameters of the regression variable at any given time; This indicates process noise.

4. The method for adaptively reconstructing missing NDVI values ​​as described in claim 3, characterized in that, In step S1, given the prior distribution and NDVI observations, the posterior distribution of the sensitivity parameters of the optimal regression variable can be solved by Kalman filtering. The expression for the posterior distribution of the sensitivity parameter of the regression variable is as follows: ； in, Sensitivity parameters representing regression variables Given observation data The prior distribution is as follows; Indicates that given prior knowledge Distribution The distribution of .

5. The method for adaptively reconstructing missing NDVI values ​​as described in claim 4, characterized in that, Sensitivity parameters of regression variables Given observation data The prior distribution of the following is represented as a Gaussian distribution with mean A and variance A; The expression for the mean A is: ； The expression for variance A is: ； In the formula, Let A represent the mean. Indicates the variance B; From the previous moment The mean of the posterior distribution is given; From the previous moment The variance of the posterior distribution is given.

6. The method for adaptively reconstructing missing NDVI values ​​as described in claim 5, characterized in that, Given prior Distribution The distribution of can be determined by a Gaussian distribution of mean B and variance B; The expression for the mean B is: ； The expression for variance B is: ； In the formula, Indicates the mean B; Let B represent the variance.

7. The method for adaptively reconstructing missing NDVI values ​​as described in claim 6, characterized in that, Based on the sensitivity parameters of the regression variables Prior distribution and The sensitivity parameters of the regression variables can be used. The posterior distribution is represented as a Gaussian distribution with mean C and variance C. The expression for the mean C is: ； The expression for variance C is: ； In the formula, This represents the mean C; Indicates the variance C; Indicates prior The difference between the distribution and the observed distribution ; This represents the Kalman gain matrix, which determines the weights of prior information. .

8. The method for adaptively reconstructing missing NDVI values ​​as described in claim 7, characterized in that, In the process of reconstructing missing observation points or identified outliers using the low-rank tensor completion algorithm, the low-rank tensor reconstruction algorithm treats noisy observations in remote sensing as missing values ​​and uses the optimization criterion of the lowest tensor rank after reconstruction to reconstruct noisy observations using observation information with good observation quality. The noise observation model reconstructed using observation information of good observation quality is as follows: ； In the formula, This indicates a good observation in the original observations; This indicates the noise-free observation after reconstruction; Tensor The weights of the nth mode.

Images