Vegetation index image reconstruction method based on earth surface remote sensing technology

Abstract

The invention discloses a vegetation index image reconstruction method based on an earth surface remote sensing technology, and belongs to the technical field of remote sensing, and the method comprises the steps: extracting earth surface reflectance data according to obtained remote sensing data, carrying out the wave band calculation, obtaining a full time sequence Landsat NDVI and a full time sequence MODIS NDVI, and respectively taking the Landsat NDVI and the MODIS NDVI as reconstruction data and reference data; extracting PIXEL QA data as quality marking data; obtaining a reconstruction result according to the quality mark data, the reconstruction data, the reference data and a pre-trained reconstruction model; constructing a time sequence arrangement diagram according to a reconstruction result to obtain a reconstructed image; wherein a three-order tensor matrix is established by using the quality marker data, the reconstruction data and the reference data, and the reconstruction model is constructed based on the three-order tensor matrix. According to the method, the defect of insufficient spatial or temporal resolution of a high-resolution image is effectively overcome, and the accuracy and integrity of data reconstruction are enhanced.

Description

Technical Field

[0001] The present invention relates to a method for reconstructing vegetation index images based on surface remote sensing technology, belonging to the field of remote sensing technology. Background Art

[0002] Time-series vegetation indices are of great value in ecological monitoring, agricultural management, and climate change research. By analyzing the changes in vegetation indices over a long time span, not only can the dynamics of vegetation cover be monitored, the health of ecosystems be evaluated, land degradation and ecological restoration trends be identified, etc., but also it has important application values in crop growth prediction, yield estimation, precision agriculture practice, meteorological disaster assessment, etc. For example, NDVI, the normalized difference vegetation index, as the most widely used vegetation index, is often used to analyze the seasonal and interannual changes of vegetation, help understand the relationship between vegetation growth and climatic factors such as precipitation and temperature, reveal the impact of climate change on vegetation growth and carbon cycle, and contribute to global change research. In short, time-series vegetation indices, with their continuity and wide applicability, have irreplaceable value in global climate change and vegetation-related research.

[0003] MODIS data is widely used for vegetation index monitoring due to its high temporal resolution and can better reflect the temporal changes of surface vegetation. Although MODIS data is also affected by clouds, there are already various methods to handle these problems, such as using the maximum value composite method (MVC) to synthesize relatively high-quality products at 16-day or even monthly scales. In addition, domestic and foreign scholars have also developed a series of time-series reconstruction methods for this type of data, including the Savitzky-Golay filter, statistical interpolation methods, etc. However, the low spatial resolution of MODIS data limits its application in the fine monitoring of surface feature changes. Compared with MODIS, Landsat data provides a higher spatial resolution of 30 meters, which can finely capture surface features and has a wide range of applications in industrial and agricultural production, resource exploration, environmental disaster assessment, etc. However, due to its low temporal resolution, Landsat data is extremely vulnerable to cloud cover, often resulting in a large number of missing values in the time series data and making it difficult to form a complete time series. Especially in areas with high cloud cover, the missing problem of Landsat data is more serious, limiting its application potential in fine-scale dynamic monitoring. Although some methods have been developed for the reconstruction of Landsat time series, stable results have not been achieved due to the lack of reference information. Summary of the Invention

[0004] The object of the present invention is to provide a method for reconstructing vegetation index images based on surface remote sensing technology. By making full use of the high frequency and full coverage of low-resolution images, it solves the problem of insufficient spatial or temporal resolution of high-resolution images in the prior art, enhances the accuracy and integrity of data reconstruction, and then implements efficient, robust, and high-precision image reconstruction.

[0005] To solve the above technical problems, the present invention is implemented by adopting the following technical solutions.

[0006] The present invention provides a method for reconstructing vegetation index images based on surface remote sensing technology, including:

[0007] According to the acquired remote sensing data, extract the surface reflectance data for band calculation to obtain the full time series Landsat NDVI and the full time series MODIS NDVI, which are used as reconstruction data and reference data respectively;

[0008] According to the Landsat QA data in the Landsat8 C2L2 data, extract the PIXEL QA data as quality marking data;

[0009] According to the quality marking data, reconstruction data, reference data, and a pre-trained reconstruction model, obtain the reconstruction result;

[0010] According to the reconstruction result, construct a time series arrangement diagram to obtain the reconstructed image;

[0011] Wherein, a third-order tensor matrix is established using the quality marking data, reconstruction data, and reference data, and the reconstruction model is constructed based on the third-order tensor matrix.

[0012] Furthermore, according to the acquired remote sensing data, extract the surface reflectance data for band calculation to obtain the normalized difference vegetation index; the normalized difference vegetation index includes the full time series Landsat NDVI and the full time series MODIS NDVI. Among them, the expression for band calculation is:

[0013] ;

[0014] NDVI is the normalized difference vegetation index, NIR is the surface reflectance data of the near-infrared band, and Red is the surface reflectance data of the infrared band.

[0015] Furthermore, establishing a third-order tensor matrix using the quality marking data, reconstruction data, and reference data, and constructing the reconstruction model based on the third-order tensor matrix includes:

[0016] Using the quality marking data, reconstruction data, and reference data to generate third-order tensor matrices respectively, the third-order tensor matrices include a marking tensor, a reconstruction tensor, and a reference tensor;

[0017] Divide the reconstructed tensor into multiple small regions, and complete the reconstructed tensor according to the abnormal points and reference tensor in the marked tensor to obtain all block results;

[0018] Process all block results in parallel using a sliding window, and perform weighted averaging on all block results to obtain a complete image tensor as the reconstruction model.

[0019] Furthermore, it also includes defining a loss function and a tensor update function to train the reconstruction model through an iterative process to obtain a trained reconstruction model;

[0020] Among them, an iterative method is used to obtain the reconstructed tensor for each round through the tensor update function;

[0021] Check the reconstructed tensor for each round according to the convergence condition to minimize the loss function, and continuously optimize the value of the reconstructed tensor during the iterative process;

[0022] Measure the difference between the reconstructed tensor for each round and the reference tensor after completion according to the singular value distribution and singular value cumulative ratio, and update the normalized weight matrix of the loss function to optimize the reconstruction model.

[0023] Furthermore, the loss function is expressed as:

[0024] ;

[0025] In the formula, is the loss value of the reconstructed tensor , is the normalized weight matrix, is the -dimensional reference variable matrix, is the -dimensional reference variable matrix 's nuclear norm, which is used to encourage to have a low-rank property, represents the total number of dimensions of the reference variable matrix, is a penalty parameter, which is used to control the importance of the quality-labeled data consistency constraint term in the loss function, and increases by 1.2 times the initial value each time as the iteration progresses, is the quality-labeled data consistency constraint, which is used to control the deviation degree of the reconstructed tensor from the quality-labeled tensor , is a regularization parameter, which is used to control the influence degree of the reference data consistency constraint on the recovery process, is the reference data consistency constraint, which is used to control the deviation degree of the reconstructed tensor from the reference tensor , where is the square of the Frobenius norm.

[0026] Further, measure the difference between the reconstructed tensor of each round and the reference tensor after completion according to the singular value distribution and the cumulative proportion of singular values, and update the normalized weight matrix of the loss function to optimize the reconstruction model, including:

[0027] Perform singular value decomposition on the reconstructed tensor of each round to obtain a set of singular values;

[0028] Divide each singular value by the sum of all singular values and then normalize it to obtain the normalized singular value distribution;

[0029] Calculate the cumulative proportion of singular values starting from the first singular value until the cumulative proportion is less than a preset threshold, and record the total number of singular values and the number of singular values when the preset threshold is reached;

[0030] According to the normalized singular value distribution vector and the number of singular values when the preset threshold is reached, calculate the proportion of the number of singular values when the preset threshold is reached to the total number of singular values among the singular values of the -th dimension of the reconstructed tensor to obtain the normalized weight matrix of the loss function to be updated, and optimize the reconstruction model according to the normalized weight matrix.

[0031] Further, the normalized singular value distribution vector is expressed as:

[0032] ;

[0033] In the formula, is the normalized singular value distribution vector, including all singular values of the reconstructed tensor , is to perform singular value decomposition on the reconstructed tensor of each round after convergence check, is the reconstructed tensor of each round after convergence check, and

[0034] The cumulative proportion of singular values is expressed as:

[0035] ;

[0036] In the formula, is the cumulative proportion of singular values, is the number of singular values when the preset threshold z is reached, represents the -th singular value in the normalized singular value distribution vector;

[0037] Reconstructed tensor The ratio of the number of singular values actually participating in calculating the cumulative proportion of singular values in the -th dimension to the total number of singular values in the normalized singular value distribution vector is expressed as:

[0038] ;

[0039] In the formula, is the ratio of the number of singular values actually participating in calculating the cumulative proportion of singular values in the -th dimension of the reconstructed tensor to the total number of singular values in the normalized singular value distribution vector, is the number of singular values actually participating in calculating the cumulative proportion of singular values, is the total number of singular values in the normalized singular value distribution vector;

[0040] The normalized weight matrix is expressed as:

[0041] ;

[0042] In the formula, is the normalized weight matrix, is the vector composed of the k(i) values corresponding to the -th dimension of the reconstructed tensor represents the total number of dimensions of the reconstructed tensor .

[0043] Furthermore, the convergence condition is expressed as:

[0044] < ;

[0045] In the formula, is the relative residual of the current iteration, , is the reconstructed tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the L2 norm, used to calculate the Euclidean distance of the difference between and .

[0046] Furthermore, the tensor update function includes a Lagrange multiplier matrix update function, a reference variable matrix update function, and a singular value decomposition update function. The tensor update function is expressed as:

[0047] ;

[0048] In the formula, is the reconstructed tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the dimensional reference variable matrix, represents the total number of dimensions of the reference variable matrix, is the sum of all Lagrange multipliers obtained through the Lagrange multiplier matrix update function ; is the regularization parameter, used to balance the weight of the reconstructed tensor in the previous iteration ; is the reference tensor, updated through the reference variable update function; is the total number of dimensions of the reconstructed tensor, used to normalize the weight.

[0049] Furthermore, the Lagrange multiplier matrix update function is expressed as:

[0050] ;

[0051] In the formula, is the Lagrange multiplier matrix of the th dimension obtained in the th iteration, is the Lagrange multiplier matrix of the th dimension obtained in the t-th iteration, is the reconstructed tensor updated in the th iteration , is the th dimension of the reference variable matrix in the th iteration ;

[0052] The reference variable matrix update function is expressed as:

[0053] ;

[0054] In the formula, is the reference variable matrix of the th dimension, is the nuclear norm of the reference variable matrix of the th dimension, used to encourage to have low rankness, , represents solving the objective function to obtain the minimum , represents the total number of dimensions of the reference variable matrix, used to constrain to be close to , balancing the low-rank representation of and the consistency of the current reconstructed tensor ; is the Lagrange multiplier matrix, used as the dual variable for constraints and update;

[0055] The singular value decomposition update function is expressed as:

[0056] ;

[0057] In the formula, is the left singular matrix, is the singular value vector, is the right singular matrix, , and together constitute the singular value decomposition result of the -dimensional reference variable matrix , is the operation of constructing a diagonal matrix, is to take the larger value of two numerical values, where is the soft threshold operation, which shrinks the singular values to ensure that the -dimensional reference variable matrix remains low-rank, is the nuclear norm regularization parameter, , is the diagonal matrix after shrinking the singular values.

[0058] Compared with the prior art, the beneficial effects achieved by the present invention are:

[0059] 1. This method combines multi-source data Landsat and MODIS and their low-rank information, effectively solving the data missing problem caused by insufficient resolution of a single data source in traditional methods.

[0060] 2. The block completion method enables the algorithm to more finely repair missing data in a local area, significantly improving the detail expression ability of the reconstruction process. At the same time, by combining the outlier analysis of the marked tensor and the reference tensor, it avoids the problem of excessive computational complexity caused by global processing, improving the computational efficiency and stability.

[0061] 3. The sliding window technology combined with parallel processing greatly improves the computational efficiency of large-scale remote sensing data processing and reduces the time cost of traditional global calculation methods; at the same time, through weighted average processing, it smooths the boundary transition between block results, making the reconstructed image more continuous and natural.

[0062] 4. By updating the normalization weights of the loss function, the reconstruction model can adaptively adjust the attention degree to data of different dimensions during the iterative optimization process, improving the adaptability and robustness of the algorithm in complex data scenarios, thereby further enhancing the accuracy and reliability of the reconstructed tensor. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 The figure shows a schematic flowchart of a vegetation index image reconstruction method based on surface remote sensing technology provided by an embodiment of the present invention;

[0064] Figure 2 The figure shows a schematic diagram of the time series of the reconstruction result provided by an embodiment of the present invention;

[0065] Figure 3 The figure shows a schematic comparison diagram of the reconstruction results obtained by the present invention and other methods provided by an embodiment of the present invention;

[0066] Figure 4 The figure shows a schematic comparison diagram of the reconstructed images obtained by the present invention and other methods provided by an embodiment of the present invention;

[0067] Figure 5 The figure shows an explanatory diagram of the pixel quality flag PIXELQA provided by an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0068] The technical solution of the present invention will be described in detail below with reference to the drawings and specific embodiments. It should be understood that the specific features in the embodiments of the present invention are detailed descriptions of the technical solution of the present invention, rather than limitations on the technical solution of the present invention. Without conflict, the technical features in the embodiments of the present invention and the embodiments can be combined with each other.

[0069] The term "and / or" is merely a description of the association relationship of associated objects, indicating that there can be three relationships. For example, A and / or B can represent: A exists alone, A and B exist simultaneously, and B exists alone. In addition, the character " / " generally represents an "or" relationship between the associated objects before and after.

[0070] Embodiment 1

[0071] As Figure 1 shown, this embodiment introduces a vegetation index image reconstruction method based on surface remote sensing technology, including:

[0072] Step 1: According to the acquired remote sensing data, extract the surface reflectance data for band calculation to obtain the full-time series Landsat NDVI and the full-time series MODIS NDVI, which are used as the reconstruction data and the reference data respectively.

[0073] In this embodiment, Landsat NDVI is the normalized difference vegetation index of Landsat satellite images, and MODIS NDVI is the normalized difference vegetation index of Moderate Resolution Imaging Spectroradiometer (MODIS) images. By using the NDVI data calculated from surface reflectance, the present invention fully combines the characteristics of high spatial resolution of Landsat and high temporal resolution of MODIS to construct reconstructed input data with both temporal continuity and spatial details, laying a foundation for the reconstruction of subsequent images in the reconstruction model.

[0074] Step 2: Extract the PIXEL QA data as quality marking data according to the Landsat QA data in Landsat8 C2L2 data.

[0075] The present invention automatically extracts the position data of valid pixels and abnormal pixels by using the marking information in the Landsat QA quality control band, and generates a quality marking matrix, providing spatial positioning support for the data quality of the reconstruction model, effectively improving the sensitivity and robustness of the reconstruction model to abnormal data.

[0076] Step 3: Obtain the reconstruction result according to the quality marking data, the reconstructed data, the reference data, and the pre-trained reconstruction model.

[0077] The present invention constructs a third-order tensor matrix based on the quality marking data, the reconstructed data, and the reference data, combines the tensor completion algorithm, provides temporal dynamic constraints with low-resolution MODIS NDVI data, and the low-rank characteristics of high spatial resolution with Landsat data. By using the pre-trained reconstruction model, the missing points are complemented through multiple iterations, thereby generating a reconstruction result with higher accuracy.

[0078] Step 4: Construct a time series arrangement diagram according to the reconstruction result to obtain the reconstructed image.

[0079] The present invention conducts temporal analysis on the reconstruction result after tensor completion, constructs a time series arrangement diagram, and performs spatial smoothing processing and detail optimization on the reconstructed image, finally generating a high-resolution continuous vegetation index image, providing a high-quality data basis for subsequent monitoring of vegetation dynamic changes and ecological analysis.

[0080] Among them, a third-order tensor matrix is established by using the quality marking data, the reconstructed data, and the reference data, and the reconstruction model is trained based on the third-order tensor matrix.

[0081] The reconstruction model is established by integrating Landsat8 C2L2 data, MODIS surface reflectance data, Landsat QA quality tag data, and reference data MODIS NDVI, which fully combines the advantages of multi-source remote sensing data and achieves high-precision reconstruction through tensor iterative completion and low-rank optimization. After threshold verification, the completed high-resolution Landsat NDVI image has been significantly improved in both temporal continuity and spatial resolution, and finally generates a high-quality remote sensing image product with high detail expression.

[0082] Example 2

[0083] Based on the same inventive concept as Example 1, this example introduces specific implementation steps of a vegetation index image reconstruction method based on surface remote sensing technology, including:

[0084] The study area is selected as an area belonging to the northwestern Hubei region, which is intercepted from the area with the row number (125, 38) of the Landsat satellite. The area belongs to the forest area with a relatively large amount of vegetation. Due to the iteration of Landsat satellites and the alternating time of different satellites, this embodiment selects Landsat8 images of each quarter from 2014 to 2020 as the reconstructed data; at the same time, MOD13Q1 data is selected as the reference data, and the time is also 2014-2020, one scene every 16 days. By weighted averaging the data, the corresponding Modis NDVI quarterly data is obtained. To facilitate the comparison of reconstruction effects, in this embodiment, the reconstruction method implemented by this patent is referred to as ST-M, and the reconstruction method that does not refer to MODIS NDVI data in this reconstruction method is called ST.

[0085] Step 1: Based on the remote sensing data obtained in the study area, the surface reflectance data is extracted for band calculation to obtain the normalized vegetation index; the normalized vegetation index includes the full time series Landsat NDVI and the full time series MODIS NDVI, which are used as reconstruction data and reference data respectively.

[0086] In this embodiment, the expression for band calculation is:

[0087] ;

[0088] NDVI is the Normalized Difference Vegetation Index, NIR is the surface reflectance data in the near-infrared band, and Red is the surface reflectance data in the infrared band.

[0089] Step 2: Extract PIXEL QA data as quality marker data based on the Landsat QA data in Landsat8 C2L2 data.

[0090] Step 3: Obtain the reconstruction result according to the quality marked data, the reconstructed data, the reference data, and the pre-trained reconstruction model.

[0091] Among them, a third-order tensor matrix is established by using the quality marked data, the reconstructed data, and the reference data, and the reconstruction model is constructed based on the third-order tensor matrix, including:

[0092] Use the quality marked data, the reconstructed data, and the reference data to generate a third-order tensor matrix respectively. The third-order tensor matrix includes a marked tensor, a reconstructed tensor, and a reference tensor;

[0093] Divide the reconstructed tensor into multiple small regions, and complete the reconstructed tensor according to the abnormal points in the marked tensor and the reference tensor to obtain all the block results;

[0094] Process all the block results in parallel using a sliding window, and perform weighted averaging on all the block results to obtain a complete image tensor as the reconstruction model.

[0095] In this embodiment, the number of iterations is set to 100 to ensure the full optimization of the reconstruction model and the convergence of the final result. Initialize the weight matrix and perform normalization processing on it so that the sum of the weights is 1 to obtain a normalized weight matrix to balance the contributions of multi-dimensional data in the tensor completion process. Set the penalty parameter , the initial value in this embodiment is 50, and it is increased to 1.5 times the original after each iteration, which is used to control the importance of the quality marked data consistency constraint term in the loss function. Set the regularization parameter , the value in this embodiment is 200, which is used to control the influence degree of the reference data consistency constraint on the recovery process. The side length M of the window is set to 8, and the step size stepsize is set to 4, which is used to define the distance that the window slides each time.

[0096] In this embodiment, a loss function and a tensor update function are also defined to train the reconstruction model through an iterative process to obtain a trained reconstruction model.

[0097] Among them, an iterative method is used to obtain the reconstructed tensor of each round through the tensor update function;

[0098] Check the reconstructed tensor of each round according to the convergence condition to minimize the loss function, and continuously optimize the value of the reconstructed tensor during the iteration process;

[0099] Measure the difference between the reconstructed tensor of each round after completion and the reference tensor according to the singular value distribution and the cumulative proportion of singular values, and update the normalized weight matrix of the loss function to optimize the reconstruction model.

[0100] In this embodiment, the loss function is expressed as:

[0101] ;

[0102] Wherein, is the reconstructed tensor of the loss value, is the normalized weight matrix, is the dimensional reference variable matrix, is the dimensional reference variable matrix of the nuclear norm, used to encourage to have a low-rank property, represents the total number of dimensions of the reference variable matrix, is the penalty parameter, used to control the importance of the quality-labeled data consistency constraint term in the loss function, and increases by 1.2 times the initial value each time as the iteration progresses, is the quality-labeled data consistency constraint, used to control the reconstructed tensor deviates from the quality-labeled tensor degree, is the regularization parameter, used to control the influence degree of the reference data consistency constraint on the recovery process, is the reference data consistency constraint, used to control the reconstructed tensor deviates from the reference tensor degree, wherein, is the square of the Frobenius norm.

[0103] In this embodiment, according to the singular value distribution and the singular value cumulative ratio, the difference between the reconstructed tensor after completion and the reference tensor in each round is measured, and the normalized weight matrix of the loss function is updated to optimize the reconstruction model, including:

[0104] Perform singular value decomposition on the reconstructed tensor of each round to obtain a set of singular values;

[0105] Divide each singular value by the sum of all singular values, and then perform normalization to obtain the normalized singular value distribution;

[0106] Calculate the singular value cumulative ratio starting from the first singular value until the cumulative ratio is less than the preset threshold, and record the total number of singular values and the number of singular values when the preset threshold is reached;

[0107] According to the normalized singular value distribution vector and the number of singular values when the preset threshold is reached, calculate the ratio of the number of singular values when the preset threshold is reached to the total number of singular values among the singular values of the dimension of the reconstructed tensor to obtain the normalized weight matrix; optimize the reconstruction model according to the normalized weight matrix.

[0108] ​In this embodiment, the normalized singular value distribution vector is expressed as:

[0109] ;

[0110] In the formula, is the normalized singular value distribution vector, which contains all the singular values of the reconstruction tensor ; is the singular value decomposition of the reconstruction tensor in each round after convergence check, is the reconstruction tensor in each round after convergence check, and is the sum of all the singular values;

[0111] The singular value cumulative ratio is expressed as:

[0112] ;

[0113] In the formula, is the singular value cumulative ratio, is the number of singular values when reaching the preset threshold , represents the -th singular value in the normalized singular value distribution vector. Among them, in this embodiment, the preset threshold = 0.85;

[0114] The ratio of the number of singular values actually participating in the calculation of the singular value cumulative ratio in the -th dimension of the reconstruction tensor to the total number of singular values in the normalized singular value distribution vector is expressed as:

[0115] ;

[0116] In the formula, is the ratio of the number of singular values actually participating in the calculation of the singular value cumulative ratio in the -th dimension of the reconstruction tensor to the total number of singular values in the normalized singular value distribution vector, is the number of singular values actually participating in the calculation of the singular value cumulative ratio, is the total number of singular values in the normalized singular value distribution vector;

[0117] The normalized weight matrix is expressed as:

[0118] ;

[0119] In the formula, is the normalized weight matrix, is composed of the reconstruction tensor The vector composed of the k(i) values corresponding to the i-th dimension, represents the reconstructed tensor of the total number of dimensions.

[0120] In this embodiment, the convergence condition is expressed as:

[0121] < ;

[0122] In the formula, is the relative residual of the current iteration, , is the reconstructed tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the L2 norm, used to calculate and the Euclidean distance of the difference.

[0123] In this embodiment, the tensor update function includes a Lagrange multiplier matrix update function, a reference variable matrix update function, and a singular value decomposition update function, and the tensor update function is expressed as:

[0124] ;

[0125] In the formula, is the reconstructed tensor of the current iteration, is the reference variable matrix of the th dimension, n represents the total number of dimensions of the reference variable matrix, is the sum of all Lagrange multipliers obtained by the Lagrange multiplier matrix update function , is the regularization parameter, used to balance the weight of the reconstructed tensor in the previous iteration, represents the previous iteration, is the reference tensor, updated by the reference variable update function; is the number of dimensions of the reconstructed tensor, used to normalize the weight.

[0126] In this embodiment, the Lagrange multiplier matrix update function is expressed as:

[0127] ;

[0128] In the formula, is the Lagrange multiplier matrix of the th iteration obtained for the th dimension, is the Lagrange multiplier matrix of the th dimension obtained for the t-th iteration, is the reconstructed tensor updated in the th iteration, , is the -dimensional reference variable matrix in the th iteration; ;

[0129] The reference variable matrix update function is expressed as:

[0130] ;

[0131] In the formula, is the -dimensional reference variable matrix, is the nuclear norm of the -dimensional reference variable matrix , used to encourage to have low rank, represents the total number of dimensions of the reference variable matrix; represents solving the objective function to obtain the minimum , is used to constrain and to be close, balancing 's low-rank representation and the consistency of the current reconstructed tensor , is the Lagrange multiplier matrix, used as a dual variable to constrain and 's update;

[0132] The singular value decomposition update function is expressed as:

[0133] ;

[0134] In the formula, is the left singular matrix, is the singular value vector, is the right singular matrix, , and together constitute the singular value decomposition result of the -dimensional reference variable matrix , is the operation of constructing a diagonal matrix, is to take the larger value of two numbers, where is the soft threshold operation, which shrinks the singular values to ensure that the -dimensional reference variable matrix remains low rank, is the nuclear norm regularization parameter, , is the diagonal matrix after shrinking the singular values.

[0135] Step 4: Calculate the normalized pixel weight matrix, reconstruct the final restored image through the final restored image function, sort by time series, and obtain the time series diagram.

[0136] In this embodiment, the function for obtaining the normalized pixel weight matrix is expressed as:

[0137] = ;

[0138] In the formula, is the normalized pixel weight matrix, which is used to weight each pixel that will form the final restored image. is the pixel coverage times matrix, which represents the number of times each pixel is covered during the sliding window block division process.

[0139] In this embodiment, the final restored image function is expressed as:

[0140] ;

[0141] In the formula, is the final restored image. is the normalized pixel weight matrix, which is used to weight each pixel. is the accumulation of all small block gray results obtained by calculating the block division results in Steps 1 - 3.

[0142] As Figure 2 shown, on the complete time series, after reconstruction using the reconstruction method (ST - M) of this embodiment, the NDVI of the missing points is closer to the reference data. Among them, AUX NDVI represents the MODIS NDVI as the reference data, which shows the true situation of ground vegetation under near - cloud - free conditions; Landsat NDVI represents the normalized vegetation index of Landsat images, which shows the observed ground vegetation situation affected by cloud cover; ST NDVI represents the normalized vegetation index obtained by the reconstruction method without referring to MODIS NDVI, which shows the ground vegetation situation after preliminary reconstruction. Compared with before reconstruction, at most time points, the NDVI value is closer to the true value; ST - M NDVI represents the normalized vegetation index obtained by the reconstruction method referring to MODIS NDVI, which shows the ground vegetation situation after further reconstruction. At most time points, it is closer to the true value than the reconstruction method without referring to MODIS data.

[0143] As Figure 3As shown, the data reconstruction for the missing points is satisfactory over the complete time series. After adopting the reconstruction method of this embodiment, the normalized difference vegetation index (NDVI) at the missing points is greater than the corresponding values in the original data. Compared with the reconstruction method (ST) that does not use MODIS NDVI as a reference, the NDVI values obtained by this reconstruction method (ST-M) are higher. Among them, Landsat NDVI represents the normalized difference vegetation index of Landsat images, showing that the normalized difference vegetation index in the original images is significantly affected by cloud cover, resulting in lower values. ST NDVI represents the normalized difference vegetation index obtained by the reconstruction method that does not refer to MODIS NDVI. This method achieves preliminary image reconstruction, resulting in a significant increase in NDVI values. ST-M NDVI represents the normalized difference vegetation index obtained by the reconstruction method that refers to MODIS NDVI. On the basis of preliminary reconstruction, this method further refers to MODIS NDVI data, thus further increasing the normalized difference vegetation index values of the reconstructed images.

[0144] As Figure 4 shown, in the image of the fourth quarter of 2014, Figure 4 the bottommost picture is the original image of the fourth quarter of 2014. Due to weather reasons in the study area, there is a lot of cloud cover and the data quality is poor; Figure 4 the picture in the upper right is the reconstructed image obtained by the reconstruction method (ST) that does not refer to MODIS NDVI. After being reconstructed by the reconstruction model that does not refer to MODIS data, a large amount of cloud is removed while retaining as much detail as possible; Figure 4 the picture in the upper left is the reconstructed image obtained by the reconstruction method (ST-M) that refers to MODIS NDVI. After being reconstructed by the reconstruction model that refers to MODIS data, the remaining cloud is further removed while retaining a large amount of detail. The image is clear and the reconstruction effect is satisfactory.

[0145] As Figure 5 shown, the following are the values corresponding to pixel points of different qualities in the PIXELQA band. In this embodiment, except for the points with a value of 21824, which are set as non-missing points, the remaining points are all set as missing points to obtain the quality-marked data.

[0146] In summary of the above embodiments, the present invention introduces low-resolution data MODIS NDVI, the normalized difference vegetation index, as the reconstruction reference data, and combines the spatio-temporal characteristics of high-resolution data Landsat NDVI to achieve the complementary advantages of multi-source data, and can reconstruct a high-precision high-resolution vegetation index sequence. The low-resolution data provides temporal information constraints for the missing areas of the high-resolution data with its excellent temporal continuity, significantly enhancing the temporal dynamic consistency of the reconstruction results. The present invention particularly adopts the variational tensor completion method. By constructing a reconstruction model, the Landsat time-series vegetation index is unfolded into a third-order tensor, and the missing point data is located by combining the pixel quality flag PIXELQA in the Landsat QA band. Using the temporal dynamic characteristics of the low-resolution MODIS NDVI data and the spatial detail characteristics of the high-resolution Landsat data, multiple iterative optimizations are performed for completion. Through the low-rank constraint of the reconstruction model, the spatial detail expressiveness of the image is enhanced, and the reconstruction accuracy and reliability are significantly improved. Finally, the reconstructed image generated by the present invention has excellent continuity and consistency in the spatio-temporal dimension, reaches a high spatial resolution of 10 meters, and can comprehensively retain the dynamic change information of surface vegetation. This method not only provides a new technical path for the reconstruction of high-resolution time-series data, but also has a wide range of application prospects, and can provide high-quality remote sensing data support for fields such as environmental monitoring, ecological protection, and agricultural and forestry management.

[0147] Those skilled in the art should understand that the embodiments of the present invention can be provided as a method, a system, or a computer program product. Therefore, the present invention can be implemented in the form of a complete hardware embodiment, a complete software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention can be implemented in the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0148] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram can be implemented by computer program instructions, and the combination of the processes and / or blocks in the flowchart and / or block diagram can also be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate a device for realizing the functions specified in one Figure 1 one process or multiple processes and / or blocks Figure 1 one block or multiple blocks.

[0149] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to operate in a particular manner, such that the instructions stored in the computer-readable memory produce a manufacture including an instruction device that implements the functions specified in one or more of the procedures Figure 1 one or more procedures and / or blocks Figure 1 specified in the block or blocks.

[0150] These computer program instructions can also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer-implemented process, whereby the instructions executed on the computer or other programmable apparatus provide steps for implementing the functions specified in one or more of the procedures Figure 1 one or more procedures and / or blocks Figure 1 specified in the block or blocks.

[0151] The embodiments of the present invention have been described above in conjunction with the accompanying drawings. However, the present invention is not limited to the above specific embodiments. The above specific embodiments are merely illustrative and not restrictive. Under the inspiration of the present invention, those of ordinary skill in the art can also make many forms without departing from the spirit and scope protected by the claims of the present invention. All of these are within the protection scope of the present invention.

Claims

1. A method for reconstructing vegetation index images based on surface remote sensing technology, characterized in that, Including: According to the acquired remote sensing data, extract the surface reflectance data for band calculation to obtain the full-time series Landsat NDVI and the full-time series MODIS NDVI, which are used as reconstruction data and reference data respectively; Extract the PIXEL QA data from the Landsat QA data in the Landsat8 C2L2 data as quality marking data; Obtain the reconstruction result according to the quality marking data, the reconstruction data, the reference data, and the pre-trained reconstruction model; Construct a time-series arrangement diagram according to the reconstruction result to obtain the reconstructed image; Among them, a third-order tensor matrix is established using the quality marking data, the reconstruction data, and the reference data, and the reconstruction model is constructed based on the third-order tensor matrix.

2. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 1, wherein According to the acquired remote sensing data, extract the surface reflectance data for band calculation to obtain the normalized difference vegetation index; the normalized difference vegetation index includes the full-time series Landsat NDVI and the full-time series MODIS NDVI, where the expression for band calculation is: ； NDVI is the normalized difference vegetation index, NIR is the surface reflectance data of the near-infrared band, and Red is the surface reflectance data of the infrared band.

3. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 1, characterized in that Establishing a third-order tensor matrix using the quality marking data, the reconstruction data, and the reference data, and constructing the reconstruction model based on the third-order tensor matrix includes: Generate a third-order tensor matrix using the quality marking data, the reconstruction data, and the reference data respectively. The third-order tensor matrix includes a marking tensor, a reconstruction tensor, and a reference tensor; Divide the reconstruction tensor into multiple small regions, and complement the reconstruction tensor according to the abnormal points in the marking tensor and the reference tensor to obtain all block results; Adopt a sliding window to process all block results in parallel, and perform weighted averaging on all block results to obtain a complete image tensor as the reconstruction model.

4. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 3, wherein, It also includes defining a loss function and a tensor update function to train the reconstruction model through an iterative process to obtain a trained reconstruction model; Among them, an iterative method is used to obtain the reconstruction tensor of each round through the tensor update function; Check the reconstruction tensor of each round according to the convergence condition to minimize the loss function, and continuously optimize the value of the reconstruction tensor during the iterative process; Measure the difference between the reconstructed tensor of each round after complementation and the reference tensor according to the singular value distribution and the cumulative proportion of singular values, and update the normalized weight matrix of the loss function to optimize the reconstruction model.

5. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 4, characterized in that The loss function is expressed as: ； In the formula, is the reconstruction tensor 's loss value, is the normalized weight matrix, is the -dimensional reference variable matrix, is the -dimensional reference variable matrix 's nuclear norm, used to encourage to have a low-rank property, represents the total number of dimensions of the reference variable matrix, is the penalty parameter, used to control the importance of the quality-labeled data consistency constraint term in the loss function, and increases by 1.2 times the initial value each time as the iteration progresses, is the quality-labeled data consistency constraint, used to control the degree to which the reconstruction tensor deviates from the quality-labeled tensor , is the regularization parameter, used to control the influence degree of the reference data consistency constraint on the recovery process, is the reference data consistency constraint, used to control the degree to which the reconstruction tensor deviates from the reference tensor , where is the square of the Frobenius norm.

6. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 5, characterized in that, Measuring the difference between the reconstructed tensor of each round after complementation and the reference tensor according to the singular value distribution and the cumulative proportion of singular values, and updating the normalized weight matrix of the loss function to optimize the reconstruction model includes: Perform singular value decomposition on the reconstructed tensor of each round to obtain a set of singular values; Divide each singular value by the sum of all singular values and then normalize it to obtain the normalized singular value distribution; Calculate the cumulative proportion of singular values starting from the first singular value until the cumulative proportion is less than a preset threshold, and record the total number of singular values and the number of singular values when the preset threshold is reached; Calculate a reconstructed tensor according to the normalized singular value distribution vector and the number of singular values when a preset threshold is reached. For the -dimensional singular values, calculate the ratio of the number of singular values that reach the preset threshold to the total number of singular values, obtain the normalized weight matrix of the loss function to be updated, and optimize the reconstruction model according to the normalized weight matrix.

7. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 6, wherein: The normalized singular value distribution vector is expressed as: ； In the formula, is the normalized singular value distribution vector, containing all the singular values of the reconstructed tensor ; is the singular value decomposition of the reconstructed tensor for each round after convergence check, is the reconstructed tensor for each round after convergence check ; and The cumulative proportion of singular values is expressed as: ； In the formula, is the cumulative proportion of singular values, is the number of singular values when reaching the preset threshold z, represents the th singular value in the normalized singular value distribution vector; Reconstructed tensor The ratio of the number of singular values actually participating in the cumulative proportion calculation of singular values in the dimension to the total number of singular values in the normalized singular value distribution vector is expressed as: ； In the formula, is the ratio of the number of singular values actually participating in the calculation of the cumulative proportion of singular values in the -th dimension of the reconstructed tensor to the total number of singular values in the normalized singular value distribution vector, where is the number of singular values actually participating in the calculation of the cumulative proportion of singular values, is the number of singular values actually participating in the calculation of the cumulative proportion of singular values, and is the total number of singular values in the normalized singular value distribution vector; The normalized weight matrix is expressed as: ； In the formula, is the normalized weight matrix, is the vector composed of the k(i) values corresponding to the i-th dimension of the reconstruction tensor , and represents the total number of dimensions of the reconstruction tensor .

8. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 5, characterized in that The convergence condition is expressed as: < ； wherein, is the relative residual of the current iteration, , is the reconstructed tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the L2 norm used to calculate and the Euclidean distance of the difference.

9. The vegetation index image reconstruction method based on surface remote sensing technology according to claim 5, characterized in that The tensor update function includes a Lagrange multiplier matrix update function, a reference variable matrix update function, and a singular value decomposition update function, and the tensor update function is expressed as: ； Wherein, is the reconstructed tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the -dimensional reference variable matrix, represents the total number of dimensions of the reference variable matrix, is the sum of all Lagrange multipliers obtained by the Lagrange multiplier matrix update function , is the regularization parameter, used to balance the weight of the reconstructed tensor in the previous iteration, is the reference tensor, updated by the reference variable update function; is the total number of dimensions of the reconstructed tensor, used to normalize the weight.

10. The method for reconstructing vegetation index images based on surface remote sensing technology according to claim 9, wherein, The Lagrange multiplier matrix update function is expressed as: ； In the formula, is the Lagrange multiplier matrix of the -th iteration and the -th dimension, is the Lagrange multiplier matrix of the -th dimension obtained in the -th iteration, is the reconstructed tensor updated in the -th iteration, is the -th iteration and the -th dimension reference variable matrix ; The reference variable matrix update function is expressed as: ； In the formula, is the -dimensional reference variable matrix, is the -dimensional reference variable matrix 's nuclear norm, which is used to encourage to have low rank, represents solving the objective function to obtain the minimum , represents the total number of dimensions of the reference variable matrix, is used to constrain and 's proximity, and balance 's low-rank representation and the consistency of the current reconstructed tensor , is the Lagrange multiplier matrix; The singular value decomposition update function is expressed as: ； wherein, is the left singular matrix, is the singular value vector, is the right singular matrix, 、 and together constitute the singular value decomposition result of the -dimensional reference variable matrix , is the operation of constructing a diagonal matrix, is to take the larger value of two numerical values, is the nuclear norm regularization parameter.

Images