CN114463202B - Vegetation index time series reconstruction method combining matrix completion and trend filtering - Google Patents

Abstract

The present invention provides a vegetation index time series reconstruction method combining matrix completion and trend filtering, comprising the following steps: firstly, determining the missing position in the time series data through the quality mark data in the vegetation index product, then, for the vegetation index sequence of a single pixel, transforming from a one-dimensional vector to a two-dimensional matrix through matrix change, then, for each pixel after the transformation, establishing the optimal energy equation for low-rank matrix completion, and implementing matrix completion through the inexact augmented Lagrangian algorithm to obtain a preliminary time series completion matrix without data missing. Finally, the completion matrix is vectorized, and an energy optimization equation for weighted trend filtering is established on the basis of the one-dimensional vector, and the model is solved through the alternating direction multiplier method, thereby further filtering out the residual noise, obtaining smooth and clean high-quality vegetation index time series data, and realizing high-precision reconstruction of long-term remote sensing vegetation index series.

Description

Vegetation index time series reconstruction method combining matrix completion and trend filtering

Technical Field

The invention relates to the field of remote sensing technology and vegetation ecological technology, and in particular to a vegetation index time series reconstruction method combining matrix completion and trend filtering.

Background Art

Remote sensing vegetation index is an important indicator for studying the status and changes of vegetation in terrestrial ecosystems, and can reflect the long-term dynamic changes of vegetation structure and function in large-scale regions. In the past few decades, with the launch of a large number of optical remote sensing satellites and sensor imagers, such as Landsat, AVHRR, MODIS, SPOT, etc., a series of remote sensing vegetation index time series with a duration of more than ten years or even decades have been formed, including the normalized vegetation index NDVI and the enhanced vegetation index EVI. At present, these standard NDVI time series products have been widely used in the fields of vegetation dynamic monitoring, agricultural yield estimation, land cover change, and terrestrial carbon cycle, and have played an important role in supporting ecological environment construction, global change response and cognition. However, in the process of optical sensor imaging, it is inevitable that it will be affected by sensor failure, clouds, aerosols, and atmospheric conditions, resulting in a lot of noise and information missing in the vegetation index time series. Although the current standard vegetation index products such as MODIS, AVHRR, and SPOT use the maximum synthesis method to reduce the impact of noise and missing information, there are still many deviations, which greatly affects their subsequent applications.

To this end, scholars at home and abroad have developed a series of time domain filtering methods to remove noise and missing problems in remote sensing vegetation index time series, which can be mainly divided into the following categories: First, based on the function fitting method, the entire vegetation index time series is fitted through a pre-set function form to remove the influence of noise. Typical methods include asymmetric Gaussian function (AG) and double logic curve (DL). This type of method has unique advantages in characterizing the seasonal regular change pattern of vegetation growth, and is therefore widely used in the extraction of phenological information, but has limited application capabilities in areas with irregular seasonal changes and complex vegetation. The second category is the filtering method based on the local sliding window, which smoothes the target time series through local windowing, such as IDR iterative interpolation and Savitzky-Golay (SG filtering). These two methods have higher computational efficiency, but the selection of different parameters has a greater impact on the results. The third category is the frequency domain method, including harmonic analysis (HANTS), wavelet transform and other methods, and their parameter determination is also a big problem. In addition, there are other methods that cannot be classified into the above three categories, such as methods based on data assimilation and Whittaker filtering.

These time-domain filtering methods have been widely used in vegetation-related research due to their simple models and high efficiency, but they usually have the problem of difficulty in handling continuous missing data in time, making it impossible to obtain high-quality NDVI time series data.

Summary of the invention

The present invention proposes a vegetation index time series reconstruction method combining matrix completion and trend filtering to solve the technical problem that high-quality NDVI time series data cannot be obtained in the prior art. Since the method proposed in the present invention takes into account the phenological periodicity and smoothness characteristics of NDVI time series data, high-quality NDVI time series data can be obtained.

The technical solution adopted by the present invention is: a vegetation index time series reconstruction method combining matrix completion and trend filtering, the method comprising the following steps:

S1: Input the original vegetation index time series data y and the corresponding quality mark data f, the quality mark data is used to characterize the quality of each pixel in the original vegetation index time series data;

S2: For each pixel's vegetation index sequence, the one-dimensional vector is transformed into a two-dimensional matrix Y through matrix transformation, where the two dimensions represent the annual coverage length of the time series and the monitoring frequency within a year, respectively, to characterize the seasonal variation characteristics of vegetation and the similarity of inter-annual periodic variation. The constructed two-dimensional matrix is the NDVI time series data matrix; then, based on the quality marker data, the missing data in the two-dimensional matrix is found, and the corresponding position is recorded as missing;

S3: For each pixel transformed two-dimensional matrix, the optimal energy equation for low-rank matrix completion is established, and the matrix completion is realized by the inexact augmented Lagrangian algorithm to obtain the preliminary time series completion matrix X without data missing;

S4: vectorize the completed matrix X to obtain the NDVI time series data x without missing items;

S5. Based on the NDVI time series data without missing items and the quality marker data f, an optimization equation for weighted trend filtering is established. Noise filtering is achieved through the alternating direction multiplier method to obtain smooth, clean, high-quality NDVI time series data z as the reconstructed NDVI time series data.

In one embodiment, in step S1, the input vegetation index time series raw data y has n years, each year has k time series points, and the total length of y and f is kn.

In one embodiment, when the vegetation index time series raw data y is data in the MODIS vegetation index product MOD13, the NDVI data is divided into four categories: clean pixels, possibly clean pixels, cloud-covered pixels, and snow-covered pixels, and the quality marks are 0, 1, 2, and 3 respectively.

In one embodiment, step S2 includes:

S2.1: Raw vegetation index time series data Matrixization, forming a matrix with k rows and n columns Each column is one year's NDVI time series data, including k time series point data, with a total of n columns;

S2.2: According to the quality marking data f, find out the pixels covered by clouds and snow, mark the found pixels as missing pixels, and record the corresponding elements in the matrix Y as missing elements. Mark other pixels as clean pixels, and record the corresponding elements in the matrix Y as clean elements.

In one embodiment, step S3 includes:

S3.1: For each pixel transformed two-dimensional matrix, establish the optimal energy equation for low-rank matrix completion:

In the formula, is the required completion matrix, i represents the i-th row of the matrix, j represents the j-th column, the subscript (i, j) represents the position of the element, and the set of the positions (i, j) of the clean elements in Y is denoted as Ω;

S3.2: We use the minimization of the matrix nuclear norm instead of the minimization of the matrix rank to transform equation (1) into a convex optimization problem, namely:

Where ||X|| _* is the nuclear norm of the matrix;

S3.3: Solve using the inexact augmented Lagrangian algorithm. First, by introducing the auxiliary variable E, equation (2) is transformed into an equivalent form:

Reconstruct the augmented Lagrangian function of equation (3), that is:

Where M is the Lagrange multiplier, μ is the penalty parameter, and equation (4) obtains the optimal solution through the following alternating iterative steps:

Wherein, the superscript k represents the number of iterations, ρ represents a constant greater than 1, which is used to update the penalty parameter μ, P _Ω (·) is the matrix projection operator, S ₁ /μ _k (·) is the singular value threshold operator, and the optimal solution X obtained by the above iterative solution is the NDVI time series data completion matrix.

In one embodiment, step S4 includes:

Complete the matrix with NDVI time series data Convert to time series vector The NDVI time series data completion matrix does not contain missing data.

In one embodiment, step S5 includes:

S5.1: According to the quality mark data f, the weight of the clean pixel is set to 1, the weight of the missing reconstructed pixel is set to 0.8, and the optimization equation of the weighted trend filter is established as:

In the formula, is the desired denoising result, λ is the regularization parameter, W is a diagonal matrix, the diagonal elements are the weights of the corresponding pixels, and D is the second-order difference matrix;

S5.2: Using the alternating direction multiplier method, first transform equation (6) into a constrained model:

Reconstruct the augmented Lagrangian function of equation (7), that is:

Where m is the Lagrange multiplier, β is the penalty parameter, and equation (6) obtains the optimal solution through the following alternating iterative steps:

Among them, θ represents a constant greater than 1, which is used to update the penalty parameter β.

is a soft threshold function, and the superscript k represents the number of iterations. The optimal solution z obtained through the above iterative solution is the final reconstruction result of NDVI time series data.

The above one or more technical solutions in the embodiments of the present application have at least one or more of the following technical effects:

(1) NDVI time series data has obvious phenological periodicity characteristics, so the matrix constructed from NDVI time series data has significant low-rank characteristics, so that missing values can be filled through low-rank matrix completion. The present invention can effectively fill the missing values in the NDVI time series, especially for the continuous missing values in the time series.

(2) In addition, considering that the NDVI time series data has the characteristic of smoothness, the noise is further filtered out through weighted trend filtering. The present invention can effectively remove the noise in the NDVI time series and obtain a smooth NDVI time series result while maintaining key points such as the maximum point and the minimum point.

In summary, the NDVI time series data reconstruction method combining low-rank matrix completion and weighted trend filtering proposed in the present invention can effectively improve the quality of NDVI time series data and has strong stability and versatility.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

FIG1 is a flow chart of a method for reconstructing time series NDVI based on low-rank matrix completion and weighted trend filtering in an embodiment of the present invention;

FIG2 is a diagram showing the variation of singular values of a time series NDVI matrix according to an embodiment of the present invention;

FIG. 3 is a comparison diagram of the original data and the reconstruction results of the time series NDVI data in an embodiment of the present invention.

DETAILED DESCRIPTION

The inventors of this application have found through extensive research and practice that the time domain filtering methods in the prior art have been widely used in vegetation-related research due to their simple models and high efficiency, but they usually have the problem of difficulty in handling continuous missing time. This is because these methods process the vegetation index sequence of each pixel separately. Therefore, when there are continuous missing time in the time series, it is difficult to have sufficient reference for effective reconstruction of the target pixel, making it difficult to obtain ideal results. In view of this, it is necessary to develop a new time domain filtering method that can effectively guarantee the processing accuracy of different missing situations while ensuring processing efficiency. By combining low-rank matrix completion and weighted trend filtering, the vegetation index time series is modeled as a whole, and the reference information of other years in the time series is effectively utilized through matrix changes, so that high-quality and efficient processing of the vegetation index time series can be achieved.

The present invention provides a time series NDVI data reconstruction method combining low-rank matrix completion and weighted trend filtering, which takes into account the phenological periodicity and smoothness characteristics of NDVI time series data, can effectively fill in the missing data in the NDVI time series data and filter out noise, and obtain high-quality NDVI time series data.

In order to achieve the above object, the main concepts of the present invention are as follows:

For the long-term remote sensing vegetation index series, the missing position in the time series data is first determined by the quality mark data in the vegetation index product. Then, for the vegetation index series of a single pixel, the matrix is transformed from a one-dimensional vector to a two-dimensional matrix. The two dimensions represent the annual coverage length of the time series and the monitoring frequency within a year, which characterizes the seasonal variation characteristics of vegetation and the similarity of inter-annual periodic variation. In addition to the time neighborhood, more reference information can be provided for the missing pixels, and high-precision reconstruction can be achieved under different missing scenarios. For each pixel after the transformation, the optimal energy equation for low-rank matrix completion is established (see formula (1)). The matrix completion is achieved by the inexact augmented Lagrangian algorithm to obtain a preliminary time series completion matrix without missing data. Finally, the completion matrix is vectorized, and the energy optimization equation for weighted trend filtering is established on the basis of the one-dimensional vector (see formula (6)). The model is solved by the alternating direction multiplier method, thereby further filtering out the residual noise and obtaining smooth and clean high-quality vegetation index time series data, realizing high-precision reconstruction of long-term remote sensing vegetation index series.

In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

The embodiment of the present invention provides a vegetation index time series reconstruction method combining matrix completion and trend filtering, including:

S1: Input the original vegetation index time series data y and the corresponding quality mark data f, the quality mark data is used to characterize the quality of each pixel in the original vegetation index time series data;

S2: For each pixel's vegetation index sequence, the one-dimensional vector is transformed into a two-dimensional matrix Y through matrix transformation, where the two dimensions represent the annual coverage length of the time series and the monitoring frequency within a year, respectively, to characterize the seasonal variation characteristics of vegetation and the similarity of inter-annual periodic variation. The constructed two-dimensional matrix is the NDVI time series data matrix; then, based on the quality marker data, the missing data in the two-dimensional matrix is found, and the corresponding position is recorded as missing;

S3: For each pixel transformed two-dimensional matrix, the optimal energy equation for low-rank matrix completion is established, and the matrix completion is realized by the inexact augmented Lagrangian algorithm to obtain the preliminary time series completion matrix X without data missing;

S4: vectorize the completed matrix X to obtain the NDVI time series data x without missing items;

S5. Based on the NDVI time series data without missing items and the quality marker data f, an optimization equation for weighted trend filtering is established. Noise filtering is achieved through the alternating direction multiplier method to obtain smooth, clean, high-quality NDVI time series data z as the reconstructed NDVI time series data.

As shown in FIG1 , it is a flow chart of a method for reconstructing NDVI time series data combining low-rank matrix completion and weighted trend filtering provided by the present invention.

In the specific implementation process, the original vegetation index time series data y is the NDVI time series data, and the quality mark data f marks the quality of each pixel in y.

S2: The missing position in the time series data is determined by the quality marker data in the vegetation index product. Then, for the vegetation index sequence of a single pixel, the one-dimensional vector is transformed into a two-dimensional matrix (NDVI time series data matrix) through matrix transformation. The two dimensions represent the annual coverage length of the time series and the monitoring frequency within a year, which characterizes the seasonal variation characteristics of vegetation and the similarity of interannual periodic changes. In this way, more reference information can be provided for the missing pixels in addition to the time neighborhood, and high-precision reconstruction can be achieved under different missing scenarios.

In step S3, for each pixel transformed matrix, the optimal energy equation for low-rank matrix completion (see formula (1)) is established, and the matrix completion is realized by the inexact augmented Lagrangian algorithm to obtain a preliminary time series completion matrix without data missing. Then, the completion matrix is vectorized by step S4, and finally, the energy optimization equation of weighted trend filtering is established on the basis of the one-dimensional vector by step S5 (see formula (6)). The model is solved by the alternating direction multiplier method, thereby further filtering out the residual noise, obtaining smooth and clean high-quality vegetation index time series data, and realizing high-precision reconstruction of long-term remote sensing vegetation index series.

Please refer to FIG2 , which is a diagram showing the variation of singular values of the time series NDVI matrix in an embodiment of the present invention.

As can be seen from Figure 2, the singular value curve of matrix Y decreases very quickly, so matrix Y has a significant low-rank characteristic, so the missing data can be reconstructed by low-rank matrix completion. The optimization equation for low-rank matrix completion is established, and the matrix completion is realized by the inexact augmented Lagrangian algorithm to obtain the NDVI time series data completion matrix X.

In one embodiment, in step S1, the input vegetation index time series raw data y has n years, each year has k time series points, and the total length of y and f is kn.

The NDVI time series data show periodic similarity characteristics with an annual cycle.

In one embodiment, when the vegetation index time series raw data y is data in the MODIS vegetation index product MOD13, the NDVI data is divided into four categories: clean pixels, possibly clean pixels, cloud-covered pixels, and snow-covered pixels, and the quality marks are 0, 1, 2, and 3 respectively.

In one embodiment, step S2 includes:

S2.1: Raw vegetation index time series data Matrixization, forming a matrix with k rows and n columns Each column is one year's NDVI time series data, including k time series point data, with a total of n columns;

S2.2: According to the quality marking data f, find out the pixels covered by clouds and snow, mark the found pixels as missing pixels, and record the corresponding elements in the matrix Y as missing elements. Mark other pixels as clean pixels, and record the corresponding elements in the matrix Y as clean elements.

In one embodiment, step S3 includes:

S3.1: For each pixel transformed two-dimensional matrix, establish the optimal energy equation for low-rank matrix completion:

In the formula, is the required completion matrix, i represents the i-th row of the matrix, j represents the j-th column, the subscript (i, j) represents the position of the element, and the set of the positions (i, j) of the clean elements in Y is denoted as Ω;

S3.2: We use the minimization of the matrix nuclear norm instead of the minimization of the matrix rank to transform equation (1) into a convex optimization problem, namely:

Where ||X|| _* is the nuclear norm of the matrix;

S3.3: Solve using the inexact augmented Lagrangian algorithm. First, by introducing the auxiliary variable E, equation (2) is transformed into an equivalent form:

Reconstruct the augmented Lagrangian function of equation (3), that is:

Where M is the Lagrange multiplier, μ is the penalty parameter, and equation (4) obtains the optimal solution through the following alternating iterative steps:

Wherein, the superscript k represents the number of iterations, ρ represents a constant greater than 1, which is used to update the penalty parameter μ, P _Ω (·) is the matrix projection operator, S ₁ /μ _k (·) is the singular value threshold operator, and the optimal solution X obtained by the above iterative solution is the NDVI time series data completion matrix.

Specifically, since the optimization equation established in step 3.1 is non-convex, the rank minimization of the matrix is replaced by the nuclear norm minimization of the matrix through S3.2, and equation (1) is transformed into a convex optimization problem.

Steps 1 to 6 in S3.3 are the algorithm for solving X.

In one embodiment, step S4 includes:

Complete the matrix with NDVI time series data Convert to time series vector The NDVI time series data completion matrix does not contain missing data.

In one embodiment, step S5 includes:

S5.1: According to the quality mark data f, the weight of the clean pixel is set to 1, the weight of the missing reconstructed pixel is set to 0.8, and the optimization equation of the weighted trend filter is established as:

In the formula, is the desired denoising result, λ is the regularization parameter, W is a diagonal matrix, the diagonal elements are the weights of the corresponding pixels, and D is the second-order difference matrix;

S5.2: Using the alternating direction multiplier method, first transform equation (6) into a constrained model:

Reconstruct the augmented Lagrangian function of equation (7), that is:

Where m is the Lagrange multiplier, β is the penalty parameter, and equation (6) obtains the optimal solution through the following alternating iterative steps:

Among them, θ represents a constant greater than 1, which is used to update the penalty parameter β.

is a soft threshold function, a and b represent the two variables of the function, and the superscript k represents the number of iterations. The optimal solution z obtained by the above iterative solution is the final reconstruction result of the NDVI time series data.

Since the optimization equation (6) cannot be directly solved analytically, the alternating direction multiplier method is used to solve it. First, equation (6) is transformed into a constrained model:

Please refer to FIG3 , which is a comparison diagram of the original data and the reconstruction result of the time series NDVI data in an embodiment of the present invention.

It should be understood that the above description of the preferred embodiment is relatively detailed and cannot be regarded as limiting the scope of patent protection of the present invention. Under the enlightenment of the present invention, ordinary technicians in this field can also make substitutions or modifications without departing from the scope of protection of the claims of the present invention, which all fall within the scope of protection of the present invention. The scope of protection requested for the present invention shall be based on the attached claims.

Claims (7)

1. The vegetation index time sequence reconstruction method combining matrix complement and trend filtering is characterized by comprising the following steps:

S1: inputting vegetation index time series original data y and quality marking data f corresponding to the vegetation index time series original data y, wherein the quality marking data are used for representing the quality of each pixel in the vegetation index time series original data;

S2: changing a vegetation index sequence of each pixel from a one-dimensional vector to a two-dimensional matrix Y through matrix change, wherein two dimensions respectively represent the annual coverage length of a time sequence and the monitoring frequency in one year, and are used for representing seasonal change characteristics and annual periodic change similarity of vegetation, and the constructed two-dimensional matrix is an NDVI time sequence data matrix; then, based on the quality mark data, finding out missing data in the two-dimensional matrix, and recording the corresponding position as missing;

S3: aiming at the two-dimensional matrix after each pixel transformation, establishing an optimized energy equation of low-rank matrix completion, and realizing matrix completion through a non-precise augmentation Lagrangian algorithm to obtain a preliminary time sequence completion matrix X without data loss;

S4: vectorizing the complement matrix X to obtain NDVI time sequence data X without deletion;

s5, based on the NDVI time sequence data without the deletion and the quality marking data f, an optimization equation of weighted trend filtering is established, noise filtering is achieved through an alternate direction multiplier method, and smooth and clean high-quality NDVI time sequence data z is obtained and used as reconstructed NDVI time sequence data.

2. The vegetation index time series reconstruction method of claim 1, wherein in step S1, the input vegetation index time series raw data y has n years, k time series points each year, and the total length of y and f is kn.

3. The method for reconstructing a vegetation index time series according to claim 1, wherein when the vegetation index time series raw data y is data in MOD is vegetation index product MOD13, NDVI data is divided into four types, i.e., clean pixel, possible clean pixel, cloud cover pixel and snow cover pixel, and the quality marks are respectively 0,1,2 and 3.

4. A method for reconstructing a time series of vegetation indices as claimed in claim 3 wherein step S2 comprises:

s2.1: time-series raw data of vegetation index Matrixing to form a matrix of k rows and n columnsWherein each column is NDVI time series data of one year, and comprises k time series point data, and n columns are all arranged;

S2.2: and according to the quality marking data f, finding out pixels covered by the cloud and snow, marking the found pixels as missing pixels, recording corresponding elements in the matrix Y as missing elements, marking other pixels as clean pixels, and recording the corresponding elements in the matrix Y as clean elements.

5. The vegetation index time series reconstruction method of claim 1, wherein step S3 comprises:

s3.1: for each two-dimensional matrix after pixel transformation, establishing an optimized energy equation of low-rank matrix completion:

In the formula, For the complement matrix to be solved, i represents the ith row of the matrix, j represents the jth column, the position of the element is represented by the subscript (i, j), and the set of the positions (i, j) of the clean elements in Y is represented as omega;

s3.2: the use of a matrix's kernel norm minimization instead of a matrix's rank minimization translates equation (1) into a convex optimization problem, namely:

Wherein X _* is the kernel norm of the matrix;

S3.3: solving by using a non-precise augmented Lagrange algorithm, firstly, converting the formula (2) into an equivalent form by introducing an auxiliary variable E:

reconstructing an augmented lagrangian function of formula (3), namely:

Wherein M is Lagrangian multiplier, mu is penalty parameter, and formula (4) obtains optimal solution through the following alternate iteration steps:

Wherein the upper index k represents the iteration number, ρ represents a constant greater than 1, used to update the punishment parameter μ, P _Ω (·) is the matrix projection operator, And (3) obtaining an optimal solution X for the singular value threshold operator through the iterative solution, and complementing the matrix for the NDVI time sequence data.

6. The vegetation index time series reconstruction method of claim 1, wherein step S4 comprises:

complementing NDVI time series data into matrix Conversion to time series vectorsThe NDVI time series data complement matrix does not contain missing data.

7. The vegetation index time series reconstruction method of claim 1, wherein step S5 comprises:

S5.1: according to the quality marking data f, the weight of the clean pixel is set to be 1, the weight of the missing reconstructed pixel is set to be 0.8, and an optimization equation for establishing weighted trend filtering is as follows:

In the formula, For the denoising result to be solved, lambda is a regularization parameter, W is a diagonal matrix, diagonal elements are weights of corresponding pixels, and D is a second-order differential matrix;

s5.2: solving by using an alternate direction multiplier method, firstly converting the formula (6) into a model with constraint:

Reconstructing an augmented lagrangian function of formula (7), namely:

wherein m is Lagrangian multiplier, beta is penalty parameter, formula (6) obtains optimal solution through the following alternate iteration steps,

Wherein θ represents a constant greater than 1, for updating the punishment parameter β,

And as a soft threshold function, the superscript k represents an optimal solution z obtained by the iterative solution of the iteration times, and the optimal solution z is a final NDVI time sequence data reconstruction result.