A cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints

Abstract

The present invention discloses a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, belonging to the field of remote sensing technology. The method comprises extracting surface reflectance data from acquired remote sensing data and performing band calculation to obtain full-time series Landsat NDVI and full-time series MODIS NDVI, which are used as reconstruction data and reference data, respectively; extracting PIXEL QA data as quality mark data; obtaining a reconstruction result based on the quality mark data, reconstruction data, reference data, and a pre-trained reconstruction model; constructing a time series arrangement diagram based on the reconstruction result to obtain a reconstructed image; wherein, a third-order tensor matrix is ​​established using the quality mark data, reconstruction data, and reference data, and the reconstruction model is constructed based on the third-order tensor matrix. The present invention effectively compensates for the insufficient spatial or temporal resolution of high-resolution images and enhances the accuracy and integrity of data reconstruction.

Description

Technical Field

[0001] The invention relates to a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, and belongs 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 changes in vegetation indices over long time spans, they can not only monitor vegetation cover dynamics, assess ecosystem health, and identify trends in land degradation and ecological restoration, but also have important application value in crop growth forecasting, yield estimation, precision agriculture practices, and meteorological disaster assessment. For example, the Normalized Difference Vegetation Index (NDVI), as the most widely used vegetation index, is often used to analyze seasonal and interannual variations in vegetation, helping to understand the relationship between vegetation growth and climatic factors such as precipitation and temperature, revealing the impact of climate change on vegetation growth and the carbon cycle, and contributing 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 are widely used for vegetation index monitoring due to their high temporal resolution, effectively reflecting temporal changes in surface vegetation. Although MODIS data are also affected by cloud cover, various methods have been developed to address this issue, such as the Maximum Composition (MVC) method, which synthesizes relatively high-quality products on 16-day and even monthly timescales. Furthermore, researchers domestically and internationally have developed a series of time series reconstruction methods for such data, including Savitzky-Golay filters and statistical interpolation methods. However, the low spatial resolution of MODIS data limits its application in fine-scale monitoring of surface feature changes. Compared to MODIS, Landsat data offers a higher spatial resolution of 30 meters, enabling detailed capture of surface features and finding widespread application in industrial and agricultural production, resource exploration, and environmental disaster assessment. However, due to its low temporal resolution, Landsat data is highly susceptible to cloud cover, often resulting in significant gaps in the time series data, making it difficult to construct a complete time series. The gaps in Landsat data are particularly severe in areas with high cloud cover, limiting its potential for fine-scale dynamic monitoring. Although some methods have been developed for the reconstruction of Landsat time series, they have not been able to achieve stable results due to the lack of reference information. Summary of the Invention

[0004] The purpose of the present invention is to provide a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints. By making full use of the high frequency and full coverage of low-resolution images, it solves the insufficient spatial or temporal resolution of high-resolution images in the existing technology, enhances the accuracy and completeness of data reconstruction, and thus implements efficient, robust and high-precision image reconstruction.

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

[0006] The present invention provides a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, comprising:

[0007] Based on the acquired remote sensing data, the surface reflectance data was extracted and band calculation was performed to obtain the full time series Landsat NDVI and the full time series MODIS NDVI, which were used as reconstruction data and reference data respectively;

[0008] According to the Landsat QA data in the Landsat8 C2L2 data, PIXEL QA data is extracted as quality mark data;

[0009] Obtain reconstruction results based on quality labeled data, reconstruction data, reference data, and pre-trained reconstruction models;

[0010] Construct a time series arrangement diagram based on the reconstruction results to obtain a reconstructed image;

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

[0012] Furthermore, based on the acquired remote sensing data, the surface reflectance data is extracted to perform 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, wherein the expression for band calculation is:

[0013] ;

[0014] 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.

[0015] Furthermore, a third-order tensor matrix is ​​established using the quality mark data, the reconstruction data, and the reference data, and the reconstruction model is constructed based on the third-order tensor matrix, including:

[0016] Using the quality mark data, the reconstruction data and the reference data, respectively, a third-order tensor matrix is ​​generated, wherein the third-order tensor matrix includes a mark tensor, a reconstruction tensor and a reference tensor;

[0017] The reconstructed tensor is divided into multiple small areas, and the reconstructed tensor is completed according to the abnormal points in the marked tensor and the reference tensor to obtain all the block results;

[0018] All the block results are processed in parallel using a sliding window, and all the block results are weighted averaged to obtain a complete image tensor as the reconstruction model.

[0019] Furthermore, the method further 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 reconstruction tensor of each round through the tensor update function;

[0021] Checking the reconstructed tensor of each round according to the convergence condition to minimize the loss function, and continuously optimizing the value of the reconstructed tensor during the iteration process;

[0022] The difference between the completed reconstructed tensor of each round and the reference tensor is measured according to the singular value distribution and the cumulative ratio of the singular values, and the normalized weight matrix of the loss function is updated to optimize the reconstruction model.

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

[0024] ;

[0025] Where, To reconstruct the tensor The loss value, is the normalized weight matrix, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage Has low-rank properties, represents the total number of dimensions of the reference variable matrix, is a penalty parameter used to control the importance of the quality label data consistency constraint in the loss function. As the iteration proceeds, it increases by 1.2 times the initial value each time. Consistency constraints for quality labeling data to control the reconstruction tensor Deviation quality marker tensor the degree, is a regularization parameter used to control the influence of the reference data consistency constraint on the recovery process. is the reference data consistency constraint used to control the reconstruction tensor Deviation from the reference tensor The degree of is the square of the Frobenius norm.

[0026] Furthermore, the difference between the completed reconstructed tensor of each round and the reference tensor is measured according to the singular value distribution and the cumulative ratio of the singular values, and the normalized weight matrix of the loss function is updated 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 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 the preset threshold, and record the total number of singular values ​​and the number of singular values ​​when the preset threshold is reached;

[0030] Calculate the reconstruction tensor based on the normalized singular value distribution vector and the number of singular values ​​when the preset threshold is reached No. The normalized weight matrix of the loss function to be updated is obtained by calculating the ratio of the number of singular values ​​of the dimension when the singular values ​​reach a preset threshold to the total number of singular values, and the reconstruction model is optimized according to the normalized weight matrix.

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

[0032] ;

[0033] Where, is the normalized singular value distribution vector, containing the reconstruction tensor All singular values ​​of For each round of reconstructed tensors after convergence check Perform singular value decomposition, The reconstructed tensor for each round after the convergence check The sum of all singular values ​​of ;

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

[0035] ;

[0036] Where, is the cumulative proportion of singular values, is the number of singular values ​​when the preset threshold z is reached, represents the first in the normalized singular value distribution vector singular values;

[0037] Reconstructing the tensor No. The ratio of the number of singular values ​​in the dimension that actually participate in the calculation of the cumulative proportion of singular values ​​to the total number of singular values ​​in the normalized singular value distribution vector is expressed as:

[0038] ;

[0039] Where, To reconstruct the tensor No. The ratio of the number of singular values ​​actually involved in calculating the cumulative proportion of singular values ​​in the dimension to the total number of singular values ​​in the normalized singular value distribution vector. is the number of singular values ​​that actually participate 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] Where, is the normalized weight matrix, Reconstruct the tensor by The vector consisting of the k(i) values ​​corresponding to the i-th dimension, Represents the reconstruction tensor The total number of dimensions.

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

[0044] < ;

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

[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, and the tensor update function is expressed as:

[0047] ;

[0048] Where, is the reconstruction tensor of the current iteration, is the reconstructed tensor in the previous iteration, For the dimensional reference variable matrix, represents the total number of dimensions of the reference variable matrix, are all Lagrange multipliers obtained by the Lagrange multiplier matrix update function of and, is a regularization parameter used to balance the reconstruction tensor in the previous iteration The weight of is the reference tensor, which is updated by the reference variable update function; The total number of dimensions of the reconstructed tensor, used to normalize the weights.

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

[0050] ;

[0051] Where, For the The first iteration dimensional Lagrange multiplier matrix, is the first dimensional Lagrange multiplier matrix, For the The reconstruction tensor updated in iterations , For the In the iteration dimensional reference variable matrix ;

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

[0053] ;

[0054] Where, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage Has low rank, Represents the objective function to be solved , get the minimum , represents the total number of dimensions of the reference variable matrix, Used for constraints and Proximity, balance The low-rank representation of the current reconstruction tensor consistency, is the Lagrange multiplier matrix, used as the dual variable to constrain and Updates;

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

[0056] ;

[0057] Where, is a left singular matrix, is the singular value vector, is a right singular matrix, 、 and Together they constitute dimensional reference variable matrix The singular value decomposition result of To construct a diagonal matrix, To take the larger of the two values, It is a soft threshold operation, which shrinks the singular value to ensure dimensional reference variable matrix Keep rank low, is the nuclear norm regularization parameter, , is the diagonal matrix after shrinking the singular values.

[0058] Compared with the prior art, the present invention has the following beneficial effects:

[0059] 1. This method combines multi-source data from 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 repair missing data more precisely in local areas, significantly improving the detail expression of the reconstruction process. At the same time, by combining the outlier point analysis of the marked tensor and the reference tensor, it avoids the problem of excessive computational complexity caused by global processing and improves computational efficiency and stability.

[0061] 3. 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, weighted averaging processing smoothes the boundary transitions between block results, making the reconstructed image more continuous and natural.

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

[0063] Figure 1 FIG2 is a flow chart of a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints provided by an embodiment of the present invention;

[0064] Figure 2 FIG2 is a timing diagram of reconstruction results provided by an embodiment of the present invention;

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

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

[0067] Figure 5 Shown is an explanatory diagram of the pixel quality label PIXELQA provided by an embodiment of the present invention. DETAILED DESCRIPTION

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

[0069] The term "and / or" simply describes a relationship between related objects, indicating that three possible relationships exist. For example, "A and / or B" can mean: A exists alone, A and B exist simultaneously, or B exists alone. Additionally, the character " / " generally indicates an "or" relationship between the related objects.

[0070] Example 1

[0071] like Figure 1 As shown, this embodiment introduces a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, including:

[0072] Step 1: Based on the acquired remote sensing data, the surface reflectance data is extracted 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.

[0073] In this embodiment, Landsat NDVI is the Landsat Normalized Difference Vegetation Index (NDVI), and MODIS NDVI is the Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI). This method utilizes NDVI data calculated from surface reflectance, leveraging the high spatial resolution of Landsat and the high temporal resolution of MODIS to construct reconstruction input data that combines temporal continuity with spatial detail, laying the foundation for subsequent image reconstruction in the reconstruction model.

[0074] Step 2: Based on the Landsat QA data in the Landsat8 C2L2 data, extract the PIXEL QA data as quality mark data.

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

[0076] Step 3: Obtain the reconstruction result based on the quality labeled data, reconstruction data, reference data, and the pre-trained reconstruction model.

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

[0078] Step 4: Construct a time series arrangement diagram based on the reconstruction results to obtain the reconstructed image.

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

[0080] Wherein, a third-order tensor matrix is ​​established using the quality mark data, the reconstruction data and the reference data, and the reconstruction model is trained based on the third-order tensor matrix.

[0081] A reconstruction model was established by integrating Landsat8 C2L2 data with MODIS surface reflectance data, Landsat QA quality-labeled data, and reference MODIS NDVI data. This model leverages the strengths 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 imagery exhibits significant improvements in both temporal continuity and spatial resolution, ultimately generating a high-quality remote sensing image product with a high degree of detail.

[0082] Example 2

[0083] Based on the same inventive concept as Example 1, this example introduces a specific implementation step of a cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, including:

[0084] The study area was selected as an area in northwestern Hubei Province, intercepted from the Landsat satellite row number (125, 38). The area is a forested area with relatively abundant vegetation. Due to the iteration of Landsat satellites and the alternating time of different satellites, this embodiment selects Landsat8 images from each quarter from 2014 to 2020 as the reconstruction data. At the same time, MOD13Q1 data is selected as reference data, also covering the period from 2014 to 2020, with one scene every 16 days. By weighted averaging the data, the corresponding MODIS NDVI quarterly data is obtained. To facilitate comparison of reconstruction effects, in this embodiment, the reconstruction method implemented in this patent is abbreviated as ST-M. The reconstruction method that does not reference MODIS NDVI data is called ST.

[0085] Step 1: Based on the remote sensing data obtained in the study area, the surface reflectance data is extracted and band calculation is performed 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: Based on the Landsat QA data in the Landsat8 C2L2 data, extract the PIXEL QA data as quality mark data.

[0090] Step 3: Obtain the reconstruction result based on the quality labeled data, reconstruction data, reference data, and the pre-trained reconstruction model.

[0091] The method of establishing a third-order tensor matrix using the quality mark data, the reconstruction data, and the reference data, and constructing the reconstruction model based on the third-order tensor matrix, includes:

[0092] Using the quality mark data, the reconstruction data and the reference data, respectively, a third-order tensor matrix is ​​generated, wherein the third-order tensor matrix includes a mark tensor, a reconstruction tensor and a reference tensor;

[0093] The reconstructed tensor is divided into multiple small areas, and the reconstructed tensor is completed according to the abnormal points in the marked tensor and the reference tensor to obtain all the block results;

[0094] All the block results are processed in parallel using a sliding window, and all the block results are weighted averaged 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. The weight matrix is ​​initialized and normalized so that the weight sum is 1 to obtain the normalized weight matrix to balance the contribution of multi-dimensional data in the tensor completion process. The penalty parameter is set. In this embodiment, the initial value is 50, which is increased to 1.5 times of the original value after each iteration. It is used to control the importance of the consistency constraint term of the quality label data in the loss function and set the regularization parameter In this embodiment, the value is 200, which is used to control the influence 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 is set to 4, which is used to define the distance 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 reconstruction tensor of each round through the tensor update function;

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

[0099] The difference between the completed reconstructed tensor of each round and the reference tensor is measured according to the singular value distribution and the cumulative ratio of the singular values, and the normalized weight matrix of the loss function is updated to optimize the reconstruction model.

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

[0101] ;

[0102] Where, To reconstruct the tensor The loss value, is the normalized weight matrix, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage Has low-rank properties, represents the total number of dimensions of the reference variable matrix, is a penalty parameter used to control the importance of the quality label data consistency constraint in the loss function. As the iteration proceeds, it increases by 1.2 times the initial value each time. Consistency constraints for quality labeling data to control the reconstruction tensor Deviation quality marker tensor the degree, is a regularization parameter used to control the influence of the reference data consistency constraint on the recovery process. is the reference data consistency constraint used to control the reconstruction tensor Deviation from the reference tensor The degree of is the square of the Frobenius norm.

[0103] In this embodiment, the difference between the completed reconstructed tensor of each round and the reference tensor is measured according to the singular value distribution and the cumulative ratio of singular values, 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 normalize to obtain the normalized singular value distribution;

[0106] Calculate the cumulative proportion of singular values ​​starting from the first singular value until the cumulative proportion 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] Calculate the reconstruction tensor based on the normalized singular value distribution vector and the number of singular values ​​when the preset threshold is reached No. The normalized weight matrix is ​​obtained by calculating the ratio of the number of singular values ​​of the dimension that reaches a preset threshold to the total number of singular values; and the reconstruction model is optimized according to the normalized weight matrix.

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

[0109] ;

[0110] Where, is the normalized singular value distribution vector, containing the reconstruction tensor All singular values ​​of For each round of reconstructed tensors after convergence check Perform singular value decomposition, The reconstructed tensor for each round after the convergence check The sum of all singular values ​​of ;

[0111] Express the cumulative proportion of singular values ​​as:

[0112] ;

[0113] Where, is the cumulative proportion of singular values, To reach the preset threshold The number of singular values ​​when , represents the first in the normalized singular value distribution vector singular values, wherein, in this embodiment, the preset threshold =0.85;

[0114] Will reconstruct the tensor No. The ratio of the number of singular values ​​in the dimension that actually participate in the calculation of the cumulative proportion of singular values ​​to the total number of singular values ​​in the normalized singular value distribution vector is expressed as:

[0115] ;

[0116] Where, To reconstruct the tensor No. The ratio of the number of singular values ​​actually involved in calculating the cumulative proportion of singular values ​​in the dimension to the total number of singular values ​​in the normalized singular value distribution vector. is the number of singular values ​​that actually participate in calculating the cumulative proportion of singular values, is the total number of singular values ​​in the normalized singular value distribution vector;

[0117] The normalized weight matrix is ​​expressed as:

[0118] ;

[0119] Where, is the normalized weight matrix, Reconstruct the tensor by The vector consisting of the k(i) values ​​corresponding to the i-th dimension, Represents the reconstruction tensor The total number of dimensions.

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

[0121] < ;

[0122] Where, is the relative residual of the current iteration, , is the reconstruction tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the L2 norm, used to calculate and 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. The tensor update function is expressed as:

[0124] ;

[0125] Where, is the reconstruction tensor of the current iteration, For the dimensional reference variable matrix, n represents the total number of dimensions of the reference variable matrix, are all Lagrange multipliers obtained by the Lagrange multiplier matrix update function of and, is a regularization parameter used to balance the reconstruction tensor in the previous iteration The weight of Indicates the last iteration, is the reference tensor, which is updated by the reference variable update function; The number of dimensions of the reconstructed tensor, used to normalize the weights.

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

[0127] ;

[0128] Where, For the The first iteration dimensional Lagrange multiplier matrix, is the first dimensional Lagrange multiplier matrix, For the The reconstruction tensor updated in iterations , For the In the iteration dimensional reference variable matrix ;

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

[0130] ;

[0131] Where, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage Has low rank, Indicates the total number of dimensions of the reference variable matrix; Represents the objective function to be solved , get the minimum , Used for constraints and Proximity, balance The low-rank representation of the current reconstruction tensor consistency, is the Lagrange multiplier matrix, used as the dual variable to constrain and Updates;

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

[0133] ;

[0134] Where, is a left singular matrix, is the singular value vector, is a right singular matrix, 、 and Together they constitute dimensional reference variable matrix The singular value decomposition result of To construct a diagonal matrix, To take the larger of the two values, It is a soft threshold operation, which shrinks the singular value to ensure dimensional reference variable matrix Keep rank low, 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 it according to the time series, and obtain the time series diagram.

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

[0137] = ;

[0138] Where, is a normalized pixel weight matrix used to weight each pixel that will constitute the final restored image. It is the pixel coverage matrix, which indicates the number of times each pixel is covered during the sliding window block process.

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

[0140] ;

[0141] Where, is the final restored image, is the normalized pixel weight matrix used to weight each pixel, Calculate the accumulation of all small block grayscale results obtained from the block results in steps 1 to 3.

[0142] like Figure 2 As shown in the figure, after reconstruction using the reconstruction method (ST-M) of this embodiment on the complete time series, 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 actual situation of ground vegetation in a nearly cloudless situation; Landsat NDVI represents the normalized vegetation index of land satellite images, which shows the observed ground vegetation affected by cloud and fog; ST NDVI represents the normalized vegetation index obtained by the reconstruction method without reference 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 with reference 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 reference to MODIS data.

[0143] like Figure 3As shown, data reconstruction for missing points in the complete time series is effective. Using the reconstruction method of this embodiment, the normalized vegetation index (NDVI) at the missing points is greater than the corresponding value in the original data. Compared to the reconstruction method (ST) that does not use MODIS NDVI as a reference, this reconstruction method (ST-M) produces higher NDVI values. Landsat NDVI represents the normalized vegetation index (NDVI) from landsat satellite imagery, demonstrating that the NDVI in the original image is significantly affected by cloud obscuration, resulting in lower values. ST NDVI represents the NDVI obtained using the reconstruction method without reference to MODIS NDVI. This method achieves preliminary image reconstruction and significantly improves NDVI values. ST-M NDVI represents the NDVI obtained using the reconstruction method with reference to MODIS NDVI. This method, based on the preliminary reconstruction, further references MODIS NDVI data, further improving the NDVI value of the reconstructed image.

[0144] like Figure 4 As shown, in the fourth quarter of 2014, Figure 4 The bottom image is the original image from the fourth quarter of 2014. The original data was affected by the weather in the study area, with heavy cloud cover and poor data quality. Figure 4 The image on the upper right is the reconstructed image obtained using the reconstruction method (ST) without reference to MODIS NDVI. After reconstruction using the reconstruction model without reference to MODIS data, a large amount of clouds is removed while retaining as much detail as possible; Figure 4 The upper left image is a reconstructed image obtained using the reconstruction method (ST-M) based on MODIS NDVI. After reconstruction using the reconstruction model based on MODIS data, the remaining clouds were further removed while retaining a large amount of details. The image is clear and the reconstruction effect is good.

[0145] like Figure 5 As shown, these are the values ​​corresponding to pixels of different qualities in the PIXELQA band. In this embodiment, points with a value of 21824 are set as non-missing points, and the remaining points are set as missing points, thereby obtaining quality mark data.

[0146] In summary, the present invention introduces low-resolution data MODIS NDVI and Normalized Difference Vegetation Index as reconstruction reference data, combined with the spatiotemporal characteristics of high-resolution data Landsat NDVI, to achieve the complementary advantages of multi-source data and reconstruct a high-precision high-resolution vegetation index sequence. The low-resolution data, with its excellent temporal continuity, provides temporal information constraints for the missing areas of the high-resolution data, significantly enhancing the temporal dynamic consistency of the reconstruction results. The present invention particularly adopts a variational tensor completion method. By constructing a reconstruction model, the Landsat time series vegetation index is expanded into a third-order tensor. The pixel quality marker PIXELQA in the Landsat QA band is combined to locate the missing point data. The temporal dynamic characteristics of the low-resolution MODIS NDVI data and the spatial detail characteristics of the high-resolution Landsat data are used to perform multiple iterative optimization and completion. Through the low-rank constraint of the reconstruction model, the spatial detail expression of the image is enhanced, and the reconstruction accuracy and reliability are significantly improved. Ultimately, the reconstructed image generated by the present invention has excellent continuity and consistency in the spatiotemporal dimension, achieving a high spatial resolution of 10 meters, and can fully preserve the dynamic change information of the surface vegetation. This method not only provides a new technical path for the reconstruction of high-resolution time series data, but also has broad application prospects. It can provide high-quality remote sensing data support for environmental monitoring, ecological protection, agricultural and forestry management and other fields.

[0147] Those skilled in the art will appreciate that embodiments of the present invention may be provided as methods, systems, or computer program products. Thus, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

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

[0149] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0150] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operational steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions executed on the computer or other programmable device for implementing the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0151] The embodiments of the present invention are described above in conjunction with the accompanying drawings, but the present invention is not limited to the above-mentioned specific implementation methods. The above-mentioned specific implementation methods are merely illustrative and not restrictive. Under the guidance of the present invention, ordinary technicians in this field can also make many forms without departing from the scope of protection of the purpose of the present invention and the claims, which are all protected by the present invention.

Claims

1. A cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints, characterized by: include: Based on the acquired remote sensing data, the surface reflectance data was extracted for band calculation to obtain the full time series Landsat NDVI and the full time series MODIS NDVI, which were used as reconstruction data and reference data respectively; According to the Landsat QA data in the Landsat8 C2L2 data, PIXEL QA data is extracted as quality mark data; Obtain reconstruction results based on quality labeled data, reconstruction data, reference data, and pre-trained reconstruction models; Construct a time series arrangement diagram based on the reconstruction results to obtain a reconstructed image; wherein, a third-order tensor matrix is ​​established using the quality mark data, the reconstruction data and the reference data, and the reconstruction model is constructed based on the third-order tensor matrix; Establishing a third-order tensor matrix using the quality mark data, the reconstruction data, and the reference data, and constructing the reconstruction model based on the third-order tensor matrix, including: Generating a third-order tensor matrix using the quality mark data, the reconstruction data and the reference data, respectively, wherein the third-order tensor matrix includes a mark tensor, a reconstruction tensor and a reference tensor; The reconstructed tensor is divided into multiple small areas, and the reconstructed tensor is completed according to the abnormal points in the marked tensor and the reference tensor to obtain all the block results; All the block results are processed in parallel using a sliding window, and all the block results are weighted averaged to obtain a complete image tensor as the reconstruction model.

2. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 1 is characterized in that: Based on the acquired remote sensing data, the surface reflectance data is extracted and band calculation is performed to obtain the normalized vegetation index; the normalized vegetation index includes the full time series Landsat NDVI and the full time series MODIS NDVI, where the band calculation expression is: ； 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.

3. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 1 is characterized in that: The method further 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; Checking the reconstructed tensor of each round according to the convergence condition to minimize the loss function, and continuously optimizing the value of the reconstructed tensor during the iteration process; The difference between the completed reconstructed tensor of each round and the reference tensor is measured according to the singular value distribution and the cumulative ratio of the singular values, and the normalized weight matrix of the loss function is updated to optimize the reconstruction model.

4. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 3 is characterized in that: The loss function is expressed as: ； In the formula, To reconstruct the tensor The loss value, is the normalized weight matrix, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage It has low-rank properties, represents the total number of dimensions of the reference variable matrix, is a penalty parameter used to control the importance of the quality label data consistency constraint in the loss function. As the iteration proceeds, it increases by 1.2 times the initial value each time. Consistency constraints for quality labeling data to control the reconstruction tensor Deviation quality marker tensor the degree, is a regularization parameter used to control the influence of the reference data consistency constraint on the recovery process. is the reference data consistency constraint used to control the reconstruction tensor Deviation from the reference tensor The degree of is the square of the Frobenius norm.

5. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 4 is characterized in that: Measuring the difference between the completed reconstructed tensor and the reference tensor in each round according to the singular value distribution and the cumulative ratio of the singular values, and updating the normalized weight matrix of the loss function to optimize the reconstruction model, including: 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 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 the preset threshold, and record the total number of singular values ​​and the number of singular values ​​when the preset threshold is reached; Calculate the reconstruction tensor based on the normalized singular value distribution vector and the number of singular values ​​when the preset threshold is reached No. The normalized weight matrix of the loss function to be updated is obtained by calculating the ratio of the number of singular values ​​in the singular values ​​of the dimension that reaches a preset threshold to the total number of singular values, and the reconstruction model is optimized according to the normalized weight matrix.

6. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 5 is characterized in that: The normalized singular value distribution vector is expressed as: ； Where, is the normalized singular value distribution vector, containing the reconstruction tensor All singular values ​​of For each round of reconstructed tensors after convergence check Perform singular value decomposition, The reconstructed tensor for each round after the convergence check The sum of all singular values ​​of ; The cumulative proportion of singular values ​​is expressed as: ； Where, is the cumulative proportion of singular values, is the number of singular values ​​when the preset threshold z is reached, represents the first in the normalized singular value distribution vector singular values; Reconstructing the tensor No. The ratio of the number of singular values ​​in the dimension that actually participate in the calculation of the cumulative proportion of singular values ​​to the total number of singular values ​​in the normalized singular value distribution vector is expressed as: ； Where, To reconstruct the tensor No. The ratio of the number of singular values ​​actually involved in calculating the cumulative proportion of singular values ​​in the dimension to the total number of singular values ​​in the normalized singular value distribution vector. is the number of singular values ​​that actually participate in calculating the cumulative proportion of singular values, is the total number of singular values ​​in the normalized singular value distribution vector; The normalized weight matrix is ​​expressed as: ； Where, is the normalized weight matrix, Reconstruct the tensor by The vector consisting of the k(i) values ​​corresponding to the i-th dimension, Represents the reconstruction tensor The total number of dimensions.

7. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 4 is characterized in that: The convergence condition is expressed as: < ； Where, is the relative residual of the current iteration, , is the reconstruction tensor of the current iteration, is the reconstructed tensor in the previous iteration, is the L2 norm, used to calculate and Euclidean distance of the difference.

8. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 4 is 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. The tensor update function is expressed as: ； Where, is the reconstruction tensor of the current iteration, is the reconstructed tensor in the previous iteration, For the dimensional reference variable matrix, represents the total number of dimensions of the reference variable matrix, are all Lagrange multipliers obtained by the Lagrange multiplier matrix update function of and, is a regularization parameter used to balance the reconstruction tensor in the previous iteration The weight of is the reference tensor, which is updated by the reference variable update function; The total number of dimensions of the reconstructed tensor, used to normalize the weights.

9. The cross-scale remote sensing vegetation index reconstruction method based on time series tensor constraints according to claim 8 is characterized in that: The Lagrange multiplier matrix update function is expressed as: ； Where, For the The first iteration dimensional Lagrange multiplier matrix, is the first dimensional Lagrange multiplier matrix, For the The reconstruction tensor updated in iterations , For the In the iteration dimensional reference variable matrix ; The reference variable matrix update function is expressed as: ； Where, For the dimensional reference variable matrix, For the dimensional reference variable matrix The nuclear norm is used to encourage Has low rank, Represents the objective function to be solved , get the minimum , represents the total number of dimensions of the reference variable matrix, Used for constraints and Proximity, balance The low-rank representation of the current reconstruction tensor consistency, is the Lagrange multiplier matrix; The singular value decomposition update function is expressed as: ； Where, is a left singular matrix, is the singular value vector, is a right singular matrix, 、 and Together they constitute dimensional reference variable matrix The singular value decomposition result of To construct a diagonal matrix, To take the larger of the two values, is the nuclear norm regularization parameter.

Images