An efficient generation method for 30-meter 8-day spatio-temporal seamless Normalized Difference Vegetation Index applicable to large regional scales

Abstract

The present invention discloses a method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales, including: selecting representative sample points, preferably homogeneous representative sample points; aggregating cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points to 250 meters as model inputs, extracting GLASS 250-meter NDVI data as model outputs, and constructing a GPR regression model to estimate the 30-meter resolution clear-sky NDVI value; for Landsat pixels with cloud cover on the target date, finding similar pixels by evaluating the similarity of time-series GLASS NDVI, constructing a time window centered on the date where they are located, and finding a reference image; establishing a regression model between the GLASS NDVI values of similar pixels on the reference image and the GLASS NDVI of similar pixels on the target image, applying the regression model to reconstruct the missing NDVI values on the target image, and weighting multiple reconstruction results to obtain a spatio-temporal seamless 30-meter 8-day Landsat NDVI; the present invention can provide important data support for regional-scale vegetation dynamic monitoring and ecological environment assessment.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing product production, in particular to an efficient generation method for a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales. Background Art

[0002] High spatio-temporal resolution NDVI products have important application values in the fields of vegetation dynamic monitoring, ecosystem assessment, agricultural management, etc. With the development of remote sensing technology, a large amount of satellite data provides a basis for the production of NDVI products. However, how to efficiently generate a 30-meter resolution NDVI time series on a large regional scale remains an important scientific issue. Currently, generating NDVI products on a large regional scale mainly faces three challenges: First, the amount of Landsat data is huge, and the traditional per-pixel NDVI calculation and reconstruction methods are inefficient and difficult to meet the requirements of operational production; second, the surface heterogeneity is significant, and simple random sampling or uniform sampling strategies are difficult to ensure the representativeness of samples, affecting the accuracy of the NDVI inversion model; third, the cloud contamination is serious, resulting in spatio-temporal discontinuity of NDVI observations and making it difficult to form a continuous and stable time series.

[0003] In response to the above problems, the existing NDVI generation methods mainly include: (1) The NDVI generation method based on block processing. The study area is divided into smaller sub-regions for separate processing, but the spatial continuity of NDVI between sub-regions is ignored, resulting in obvious stitching marks in the product; (2) The NDVI reconstruction method based on post-classification processing. NDVI reconstruction models are established separately according to different vegetation types, but the continuous change characteristics of vegetation growth conditions are not considered, and it is difficult to accurately depict the spatial gradual change of NDVI; (3) The NDVI generation method based on multi-source data fusion. By combining NDVI products with different spatial resolutions for reconstruction, but the existing methods are often computationally complex and lack physical constraints on the scale conversion process.

[0004] In terms of sample selection, existing methods fail to fully consider the influence of topographic and vegetation factors on the spatial distribution of NDVI, resulting in insufficient sample representativeness. At the same time, the evaluation of sample homogeneity is not strict enough, and the existence of mixed pixels affects the quality of model training. In terms of model construction, current research mostly uses complex machine learning models, such as deep neural networks, etc. These models have large computational overhead and difficult parameter optimization when applied to large areas. While some simple empirical models have high computational efficiency, their accuracy often fails to meet the requirements. In terms of cloud-contaminated data reconstruction, existing methods mainly rely on time-series interpolation or spatial neighborhood interpolation, lacking systematic utilization of multi-temporal and multi-scale NDVI data. At the same time, the adaptability of the reconstruction strategy is insufficient and it is difficult to cope with different degrees of cloud contamination. In recent years, although some research has begun to focus on the efficiency of large-area NDVI product production, most methods are still limited to simple extensions of small-area algorithms and have not fundamentally solved the efficiency bottleneck of large-area processing. In addition, these methods often neglect efficiency while pursuing accuracy, resulting in huge consumption of computing resources and difficulty in realizing operational applications.

[0005] Therefore, there is an urgent need to develop a new method to significantly improve the processing efficiency while ensuring the accuracy of NDVI products, and to achieve the rapid production of 30-meter resolution NDVI products at the large-area scale. This requires systematic optimization of each link from sample selection, model construction to cloud contamination reconstruction to achieve the goals of high efficiency and high accuracy. Summary of the Invention

[0006] To solve the problems existing in the prior art, the purpose of the present invention is to provide an efficient generation method for 30-meter 8-day spatio-temporally seamless normalized difference vegetation index applicable to large-area scales. The present invention can quickly generate high-resolution NDVI products with long time series, providing important data support for regional-scale vegetation dynamic monitoring and ecological environment assessment.

[0007] To achieve the above purpose, the technical solution adopted by the present invention is: an efficient generation method for 30-meter 8-day spatio-temporally seamless normalized difference vegetation index applicable to large-area scales, including the following steps:

[0008] S1: Using topographic factors and vegetation factors as constraints, a constraint-based sampling method based on representativeness levels is used to select representative sample points, and homogeneous representative sample points are further optimized based on land surface classification products;

[0009] S2: Aggregate cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points to 250 meters as model inputs, extract GLASS 250-meter NDVI data as model outputs, and construct a GPR regression model to estimate 30-meter resolution clear-sky NDVI values from Landsat surface reflectance and angular data;

[0010] S3: For the Landsat pixels with cloud cover on the target date, construct a spatial window centered on its location, find similar pixels by evaluating the similarity of the time-series GLASS NDVI, construct a time window centered on its date, and find the reference image by calculating the structural similarity index;

[0011] S4: Establish a regression model between the GLASS NDVI values of the similar pixels on the reference image and the GLASS NDVI of the similar pixels on the target image. Apply the regression model to reconstruct the missing NDVI values on the target image from the 30-meter Landsat clear-sky NDVI corresponding to the reference image, and weight multiple reconstruction results to obtain the final 30-meter 8-day Landsat NDVI with seamless spatio-temporal coverage.

[0012] As a further improvement of the present invention, step S1 is specifically as follows:

[0013] Select terrain and vegetation as environmental factors to reflect the spatial distribution differences of surface features; in terms of terrain factors, select elevation, slope, aspect, and curvature along contour lines in the digital elevation model (DEM) to reflect the controlling effect of terrain; in terms of vegetation factors, select leaf area index and vegetation coverage parameters to represent the constraint conditions of vegetation;

[0014] Use the fuzzy c-means clustering method to perform fuzzy clustering on environmental factors, determine the main distribution areas of different surface features, generate a frequency distribution map of the central positions of environmental factor combinations through overlay analysis, and determine the level of representativeness according to the frequency of the central positions; express the environmental factor combinations appearing on each pixel as an environmental factor combination chain, and design sample points according to the level of the average membership value of the corresponding pixels;

[0015] Obtain the 30-meter surface classification product corresponding to the sample points. If, under a sample point, that is, among all the 30-meter pixels corresponding to 250 meters, more than a preset proportion belong to the same class, then this sample point is considered a homogeneous sample point.

[0016] As a further improvement of the present invention, step S2 is specifically as follows:

[0017] When aggregating the cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points, multiple periods of data are involved. Only when all the 30-meter Landsat pixels corresponding to a sample point in a certain period are marked as clear sky, aggregate its 30-meter Landsat surface reflectance data and angular data to 250 meters as the model input;

[0018] When constructing the GPR regression model, use the ten-fold cross-validation method to evaluate the generalization ability of the model, construct the covariance matrix using the kernel function, and maximize the log-likelihood function by optimizing the hyperparameters to finally obtain the optimal GPR model parameter combination.

[0019] As a further improvement of the present invention, step S3 is specifically as follows:

[0020] Construct a 50km×50km spatial window with the GLASS 250m target pixel corresponding to the Landsat pixel obscured by the target cloud as the center. This spatial window corresponds to 200×200 GLASS pixels. For each pixel within the spatial window, average the GLASS NDVI time series over 5 years to obtain the NDVI reference sequence of this pixel, calculate the correlation coefficient between the NDVI reference sequence of the target pixel and the NDVI reference sequences of other pixels, and the pixels with a correlation coefficient higher than 0.8 are determined as similar pixels;

[0021] Construct a 5-year time window with the date corresponding to the Landsat pixel obscured by the target cloud as the center. The GLASS NDVI image corresponding to the target date is called the target NDVI, and the GLASS NDVI images corresponding to other dates within the time window are called reference NDVI. Calculate the structural similarity index between all reference NDVI and the target NDVI, and sort them in descending order of the structural similarity index to select the reference image accordingly.

[0022] As a further improvement of the present invention, step S4 is specifically as follows:

[0023] For the Landsat NDVI image only contaminated by clouds, first select the reference image with the highest structural similarity index, establish a regression model between the GLASS NDVI values of the similar pixels on the reference image and the GLASS NDVI of the similar pixels on the target image, and apply this regression model to reconstruct the missing NDVI values on the target image from the Landsat clear-sky 30m NDVI corresponding to the reference image; if the missing NDVI values on the target image are not fully reconstructed, continue to reconstruct using the reference image with the second highest structural similarity index until all the missing NDVI values on the target image are reconstructed.

[0024] As a further improvement of the present invention, when using multiple reference images to reconstruct the missing NDVI values of the same target image, if there are multiple reconstruction results, calculate the weights for each reconstruction result based on GLASS NDVI, and weight all the reconstruction results to obtain the final spatio-temporal seamless 30m 8-day Landsat NDVI.

[0025] The beneficial effects of the present invention are as follows:

[0026] 1. The present invention establishes a complete sample selection system based on environmental factors, significantly improving the representativeness and reliability of the samples and providing high-quality training data for subsequent modeling.

[0027] 2. The present invention adopts an efficient GPR regression model, which greatly improves the calculation efficiency while ensuring the accuracy, making it possible to estimate NDVI at a large regional scale.

[0028] 3. The present invention designs a reconstruction strategy based on multiple similarities, makes full use of spatio-temporal information, and effectively solves the problem of data loss caused by cloud pollution.

[0029] 4. The present invention can efficiently generate 30-meter 8-day NDVI products at a large regional scale, providing important data support for vegetation dynamic monitoring and ecological environment assessment at the regional scale. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is a flowchart of an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0031] The embodiments of the present invention will be described in detail below with reference to the drawings.

[0032] Embodiment 1

[0033] As Figure 1 shown, a method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized vegetation difference index applicable to a large regional scale includes:

[0034] S1: Using topographic factors and vegetation factors as constraints, a representative sampling method based on representative levels is adopted to select representative sample points, and homogeneous representative sample points are further optimized according to the surface classification product;

[0035] S2: Aggregate cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points to 250 meters as the model input, extract GLASS 250-meter NDVI data as the model output, and construct a GPR regression model to estimate the 30-meter resolution clear-sky NDVI value from Landsat surface reflectance and angular data;

[0036] S3: For Landsat pixels obscured by clouds on the target date, a 50 km × 50 km spatial window is constructed centered on its location, similar pixels are found by evaluating the similarity of the time series GLASS NDVI, and a 5-year time window is constructed centered on its date, and a reference image is found by calculating the structural similarity index;

[0037] S4: Establish a regression model between the GLASS NDVI values of similar pixels on the reference image and the GLASS NDVI of similar pixels on the target image, apply the regression model to reconstruct the missing NDVI value on the target image from the Landsat clear-sky 30-meter NDVI corresponding to the reference image, and weight multiple reconstruction results to obtain the final spatio-temporal seamless 30-meter 8-day Landsat NDVI.

[0038] In one implementation, step S1 further includes:

[0039] Taking China as an example, first, a comprehensive environmental factor database is established. Topographic factors such as elevation, slope, aspect, and curvature along contour lines are extracted from SRTM DEM. These factors can comprehensively reflect the surface undulation characteristics and the spatial differentiation law of hydrothermal conditions. At the same time, vegetation factors such as leaf area index and vegetation coverage are obtained from GLASS products to characterize the spatial distribution characteristics of vegetation growth status.

[0040] To identify representative sample points, the fuzzy c-means clustering method is used to systematically analyze the environmental factors. The maximum number of iterations is set to 1000, the fuzzy factor is set to 2, and the convergence threshold is set to 0.001. These parameter settings have been verified through repeated experiments and can better reflect the spatial differentiation characteristics of environmental factors. A frequency distribution map of the central positions of environmental factor combinations is generated through overlay analysis, the frequency of each position is calculated, and the representative levels of different regions are determined according to the frequency. On this basis, combined with the GlobalLand30 land cover product, a homogeneity evaluation is carried out, and sample points where more than 90% of the 30-meter pixels within a 250-meter pixel belong to the same type are strictly screened. This multi-level screening strategy ensures the representativeness and reliability of the samples, and finally about 200,000 high-quality training samples are selected in the study area.

[0041] In one implementation, step S2 includes:

[0042] In the data aggregation stage, a strict quality control strategy is adopted. First, the Fmask4.6 algorithm is used to detect clouds in Landsat data. This algorithm has high accuracy in identifying clouds, cloud shadows, and snow cover. Data aggregation at 250-meter resolution is only performed when all 30-meter pixels corresponding to the sample points in a certain period are all marked as clear sky. For the samples that pass the screening, their surface reflectance data and geometric information such as solar zenith angle are aggregated to 250-meter resolution, and the corresponding GLASS NDVI data is extracted as the model output.

[0043] The construction of the GPR model adopts a systematic optimization strategy. The generalization ability of the model is evaluated through the ten-fold cross-validation method, the covariance matrix is constructed using the kernel function, and the hyperparameters are optimized by maximizing the log-likelihood function. To improve the processing efficiency of large areas, a partitioned parallel processing strategy based on UTM projection zones is adopted, and the study area is divided into sub-regions of appropriate size for simultaneous processing. After the model training is completed, the boundaries of each sub-region are smoothed to ensure the spatial continuity of the final result.

[0044] In one implementation, step S3 includes:

[0045] Through a large number of experiments, it is determined that it is most appropriate to construct a 50km×50km spatial window centered on the target cloud-obscured pixel. This range can not only contain enough potential similar pixels but also ensure the relative consistency of surface conditions. The spatial window corresponds to 200×200 GLASS pixels. This scale not only meets the needs of statistical analysis but also facilitates parallel processing. For each pixel, collect the GLASS NDVI observation values within 5 years and obtain the NDVI reference sequence through time series averaging. This method of multi-year averaging can effectively reduce the influence of outliers and noise.

[0046] In the time dimension, construct a 5-year time window centered on the target date. The reason for choosing 5 years as the time window is that this time span is long enough to contain enough potential reference images, and at the same time, it is not too long to avoid significant changes in land cover. The GLASS NDVI corresponding to the target date is defined as the target NDVI, and the GLASS NDVI of other dates within the time window is used as the candidate reference NDVI. Calculate the structural similarity index to evaluate the similarity between them. This index takes into account the brightness, contrast, and structural features of the image and can comprehensively reflect the similarity between images.

[0047] In one implementation, step S4 includes:

[0048] For each Landsat image contaminated only by clouds in the study area, first select the optimal reference image based on the ranking results of the structural similarity index. On this reference image, establish a regression model between the GLASS NDVI value and the GLASS NDVI value of the target image using the previously identified similar pixels. This regression relationship based on GLASS NDVI can capture the systematic differences in NDVI values between different time phases. Use the established regression model to reconstruct the missing values on the target image from the Landsat clear-sky 30-meter NDVI corresponding to the reference image.

[0049] Considering that a single reference image may not fully cover all cloud-contaminated areas, an iterative reconstruction strategy is adopted. When the optimal reference image cannot be fully reconstructed, successively use the reference images with higher rankings of the structural similarity index to continue the reconstruction until all missing areas are filled. During the reconstruction process, the reconstruction results at the same location may come from multiple reference images. At this time, calculate the weight coefficient of each reconstruction result based on GLASS NDVI. The weight calculation takes into account factors such as the structural similarity index and time distance of the reference image to ensure the reliability of the final result.

[0050] Through the above steps, finally generate a Landsat NDVI product with a 30-meter spatial resolution and an 8-day temporal resolution in the Chinese region from 2000 to 2023.

[0051] In this embodiment, through an effective sample selection strategy and an efficient model construction method, the computational efficiency and accuracy problems in the generation of high spatio-temporal resolution NDVI products at a large regional scale are successfully solved, and high-precision, spatio-temporally continuous 30-meter resolution NDVI products are generated, providing important data support for regional-scale vegetation dynamic monitoring and ecological environment assessment; by integrating a constrained sampling and an efficient reconstruction strategy, taking full advantage of multi-source remote sensing data and advanced algorithms, considering the spatial heterogeneity and temporal continuity characteristics of NDVI, the rapid generation of high spatio-temporal resolution NDVI products at a large regional scale is realized.

[0052] Example 2

[0053] An efficient method for generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index includes the following steps:

[0054] Step 1: Using topographic factors and vegetation factors as constraints, a constrained sampling method based on representative levels is adopted to select representative sample points, and homogeneous representative sample points are further optimized according to the surface classification product.

[0055] Taking the Chinese region as an example, a complete set of environmental factor systems is first constructed. In terms of topographic factors, basic parameters such as elevation, slope, and aspect are extracted from the 30-meter resolution SRTM DEM, and composite indices such as terrain relief and curvature along contour lines are calculated at the same time. Considering the regional differences in the terrain of China, different parameter thresholds are adopted for different geographical sub-regions: in the plain area, micro-topographic changes are focused on, while in the mountainous area, the expression of large-scale topographic features is more emphasized. In terms of vegetation factors, in addition to the basic leaf area index and vegetation coverage, auxiliary indices such as the normalized difference water index (NDWI) and the enhanced vegetation index (EVI) are introduced. These indices characterize the growth status and water characteristics of vegetation from different angles, improving the representativeness of the samples. For regions with obvious seasonality, such as Northeast China, the spatial differentiation law of phenological characteristics is also particularly concerned.

[0056] When using the fuzzy c-means clustering method for analysis, different parameters are set for different geographical sub-regions. For example, in the southern region with high vegetation coverage, the fuzzy factor is increased to better characterize the gradual change characteristics; in the northern region with relatively single vegetation types, the fuzzy factor can be appropriately reduced. The clustering results are verified through spatial continuity analysis and geographical expert knowledge to ensure their rationality. When generating the frequency distribution map of the central positions of the environmental factor combinations through overlay analysis, a zonal weighting strategy is adopted. For regions with complex terrain, such as the Qinghai-Tibet Plateau, the weight of the topographic factors is increased; for regions with significant vegetation changes, such as the Northeast Plain, the influence of the vegetation factors is strengthened. Finally, according to the frequency of occurrence of the central positions of the combinations, the representative levels are divided into five levels, and a minimum sample size requirement is set for each level.

[0057] In the homogeneity evaluation stage, a dual-verification strategy is adopted in combination with the latest Global Land 30 land surface classification product and Chinese land use data. Only when more than 90% of the 30-meter pixels within a 250-meter sample point not only belong to the same type but also remain stable in the past five years, can it be determined as the final homogeneous sample point. Through such strict screening, approximately 1 million high-quality training samples are finally selected nationwide.

[0058] Step 2: Construct a high-resolution NDVI inversion strategy based on the GPR regression model. By aggregating the processed Landsat surface reflectance data, angular information, and GLASS NDVI data, establish the mapping relationship from multi-source remote sensing observations to high-resolution NDVI, and achieve the efficient estimation of 30-meter resolution NDVI at the large regional scale;

[0059] During the implementation in the Chinese region, targeted designs are first carried out for data preprocessing and aggregation. Considering the vast territory, diverse climate types, and significant differences in atmospheric conditions and observation geometric relationships in China, a processing strategy of dividing regions and time periods is adopted for data aggregation. Based on climate and surface characteristics, the whole country is divided into six major regions: Northeast Region, North China Region, South China Region, Northwest Region, Qinghai-Tibet Plateau Region, and Southwest Region, and data aggregation and model training are carried out independently for each region.

[0060] A strict quality control strategy is adopted in the data aggregation stage. For the Landsat pixels corresponding to each 250-meter sample point, the Fmask4.6 algorithm is first used for cloud detection. This algorithm is optimized for the atmospheric and surface conditions in different regions of China: in the cloudy areas along the southeast coast, the recognition ability of thin clouds and cirrus clouds is enhanced; in the western plateau and mountainous areas, the discrimination rules for cloud-snow confusion are improved; in the northern regions, the recognition threshold for dust weather is optimized. Only when all 30-meter pixels corresponding to a sample point in a certain period are marked as clear sky, can the data aggregation at 250-meter resolution be carried out.

[0061] When constructing the GPR regression model, a kernel function design strategy based on physical knowledge is adopted. Considering the non-linear relationship between NDVI and surface reflectance, the RBF kernel function is selected as the basis, and a seasonal cycle term is introduced to capture the phenological change characteristics of vegetation. The model parameter optimization adopts a two-stage strategy: first, use grid search in each region to determine the initial range of parameters, and then finely adjust through the Bayesian optimization method to maximize the log-likelihood function.

[0062] To ensure the generalization performance of the model, a hierarchical cross-validation strategy is adopted. In the spatial dimension, the samples in each partition are divided into 10 sub-regions according to the longitude and latitude grids; in the time dimension, considering the significant seasonal differences in China, the samples are grouped according to the phenological period (instead of simply by month). Through this hierarchical cross-validation, it is ensured that the model has stable performance in different geographical locations and different growth stages.

[0063] Step 3, establish a similarity analysis method based on spatio-temporal windows. For the areas affected by cloud contamination, construct large-scale spatial windows and long-time series windows, and identify the best reference data by evaluating the similarity characteristics of GLASS NDVI, providing a reliable benchmark for subsequent NDVI reconstruction;

[0064] During the implementation in the Chinese region, the construction of spatial windows needs to consider the complex and diverse geographical environment characteristics. Taking the target cloud-contaminated pixel as the center, set dynamic spatial windows based on the partitioning strategy: in the eastern plain region, considering the flat terrain and relatively homogeneous vegetation types, the spatial window can be expanded to 70km×70km; in the western mountainous area, due to the fragmented terrain and diverse vegetation types, the spatial window is reduced to 30km×30km; in the transitional zone of the Qinghai-Tibet Plateau, set irregularly shaped spatial windows according to the altitude gradient to ensure the relative consistency of surface conditions within the search range. For each pixel within the spatial window, extract its 5-year GLASS NDVI time series. Considering the vegetation growth characteristics in different climate regions of our country, the processing of time series adopts a partitioning strategy: in the temperate region, focus on analyzing the dynamic characteristics of NDVI during the growing season; in the subtropical and tropical regions, pay more attention to the interannual variation law of NDVI; in the arid and semi-arid regions, especially note the response characteristics of vegetation to precipitation. Through this partitioned time series analysis, the pertinence of similarity evaluation is improved. The similarity evaluation adopts a multi-index comprehensive analysis method. First, calculate the correlation coefficient of the NDVI time series and set a basic threshold of 0.8. However, in different regions, dynamically adjust the threshold according to the degree of vegetation heterogeneity: in the concentrated farmland area, due to the influence of farming activities, appropriately increase the threshold to 0.85; in the natural vegetation area, considering the gradual change characteristics of community succession, the threshold can be appropriately reduced to 0.75. In addition to the correlation coefficient, also calculate the main phenological characteristics of the time series (such as the start and end times of the growing season, peak values, etc.), and require that the difference in key phenological periods between similar pixels does not exceed 8 days.

[0065] In terms of time window construction, expand 5 years forward and backward centered on the target date. The actual utilization strategy of the time window varies according to local conditions: in the northern regions, focus on using historical data of the same season; in the evergreen regions of the south, data from different seasons can be used more flexibly. Calculate the structural similarity index for each date image within the time window, which comprehensively considers the brightness, contrast, and structural features of the image. The calculation parameters of the structural similarity index are adjusted according to regional characteristics. In the plain agricultural areas, more attention is paid to the consistency of local structural features; in forest areas, the weights of brightness and contrast are appropriately increased; in the farmland-urban mixed areas around cities, a boundary preservation term is introduced to avoid incorrect matching between different land use types. Based on the calculated structural similarity index, sort all candidate images within the time window to provide a priority reference for subsequent NDVI reconstruction. To improve the calculation efficiency, a hierarchical parallel computing strategy is adopted for the entire similarity analysis process. First, conduct a coarse-scale similarity assessment at the GLASS NDVI resolution to screen out potential candidate regions, and then perform refined analysis within these regions. At the same time, use the boundary information of geographical partitions to avoid similarity searches across significant geographical boundaries, such as regions on both sides of the Hengduan Mountains. In actual processing, to cope with possible data missing situations, alternative solutions are also established. When the available data within the 5-year time window is insufficient (such as in areas with poor observation conditions like high mountains and deserts), ensure sufficient reference data by expanding the time window or adjusting the similarity threshold, etc. At the same time, record the data quality indicators to provide a basis for uncertainty assessment of subsequent reconstruction results.

[0066] Step 4: Implement an efficient NDVI reconstruction strategy based on multi-reference images. By establishing the regression relationship between the reference image and the target image, using the structural similarity index to guide the reconstruction process, and adopting a weighted fusion method based on GLASS NDVI, finally generate a high-resolution NDVI product with spatio-temporal continuity.

[0067] The establishment of the regression model adopts a multi-level strategy. First, use the reference image with the highest structural similarity index. Considering the diversity of vegetation types in the Chinese region, introduce the surface type factor when establishing the regression relationship between GLASS NDVI values. For farmland areas, use piecewise linear regression to capture the phased characteristics of crop growth; for forest areas, adopt polynomial regression to better describe the continuous changes; for grassland areas, introduce the precipitation lag effect as an auxiliary variable. Implement strict quality control during the reconstruction process. When using the first reference image for reconstruction, set the accuracy threshold based on vegetation types: the relative error between the reconstruction result and GLASS NDVI in the farmland area should be less than 10%, 15% in the forest area, and 20% in the grassland area. For areas that fail to meet the accuracy requirements or are not fully reconstructed, start the iterative reconstruction process based on the sorting of the structural similarity index.

[0068] For pixels with multiple reconstruction results, an intelligent weight assignment system was developed. The weight calculation not only considers the structural similarity index of the reference image, but also introduces a time-distance decay factor and an observation quality index. In areas with obvious seasons, the weight of the reference image in the same phenological period is increased; in areas with drastic changes, the weight of the reconstruction result with a closer time distance is increased. To improve the processing efficiency, a distributed computing framework was established. The processing tasks across the country were organized according to the UTM projection zones, and the size of each block was set to 100 km × 100 km, ensuring a certain overlap to eliminate the boundary effect. Within each block, parallel reconstruction was carried out based on the parameters designed according to the geographical partition, and the reconstruction results were seamlessly mosaicked after passing the spatial consistency test.

[0069] The quality assessment of the final product is cross-validated by multiple methods: direct comparison with the measurement data of ground flux stations; upscaling comparison with existing medium-resolution NDVI products such as MODIS; verification of typical sample areas using high-resolution observation data of unmanned aerial vehicles. The verification results show that the overall accuracy of the product reaches R 2 > 0.85 nationwide, and more than 95% of the pixels show reasonable seasonal variation characteristics in the time continuity evaluation.

[0070] The above-described embodiments merely represent the specific implementation manners of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the invention patent. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications and improvements can still be made, and these all fall within the protection scope of the present invention.

Claims

1. An efficient generation method for 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales, characterized in that, Including the following steps: S1: Taking the terrain factor and vegetation factor as constraints, a constrained sampling method based on representative levels is used to select representative sample points, and homogeneous representative sample points are further optimized according to the surface classification product; S2: Aggregate the cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points to 250 meters as the model input, extract the GLASS 250-meter NDVI data as the model output, and construct a GPR regression model to estimate the clear-sky NDVI value at 30-meter resolution from the Landsat surface reflectance and angular data; S3: For Landsat pixels obscured by clouds on the target date, construct a spatial window centered on its location, find similar pixels by evaluating the similarity of time-series GLASS NDVI, construct a time window centered on its date, and find the reference image by calculating the structural similarity index; S4: Establish a regression model between the GLASS NDVI values of similar pixels on the reference image and the GLASS NDVI of similar pixels on the target image, apply the regression model to reconstruct the missing NDVI value on the target image from the Landsat clear-sky 30-meter NDVI corresponding to the reference image, and weight multiple reconstruction results to obtain the final spatio-temporally seamless 30-meter 8-day Landsat NDVI.

2. The method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales according to claim 1, characterized in that, The specific steps of S1 are as follows: Select terrain and vegetation as environmental factors to reflect the spatial distribution differences of surface features; in terms of terrain factors, select elevation, slope, aspect, and curvature along contour lines in the digital elevation model (DEM) to reflect the control of terrain; in terms of vegetation factors, select leaf area index and vegetation coverage parameters to represent the constraint conditions of vegetation; Use the fuzzy c-means clustering method to perform fuzzy clustering on environmental factors, determine the main distribution areas of different surface features, generate a frequency distribution map of the central positions of environmental factor combinations through overlay analysis, and determine the level of representativeness according to the frequency of the central positions; express the environmental factor combinations appearing on each pixel as an environmental factor combination chain, and design sample points according to the level of the average membership value of its corresponding pixel; Obtain the 30-meter surface classification product corresponding to the sample points. If, under a sample point, that is, more than a preset proportion of all 30-meter pixels corresponding to 250 meters belong to the same class, then this sample point is considered a homogeneous sample point.

3. The method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales according to claim 1, wherein The specific steps of S2 are as follows: When aggregating the cloud-free 30-meter Landsat surface reflectance data and angular data of representative sample points, multiple periods of data are involved. Only when all 30-meter Landsat pixels corresponding to this sample point in a certain period are marked as clear sky, aggregate its 30-meter Landsat surface reflectance data and angular data to 250 meters as the model input; When constructing the GPR regression model, use the ten-fold cross-validation method to evaluate the generalization ability of the model, construct the covariance matrix using the kernel function, and maximize the log-likelihood function by optimizing the hyperparameters to finally obtain the optimal GPR model parameter combination.

4. The method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales according to claim 1, characterized in that, The specific steps of S3 are as follows: Construct a 50 km × 50 km spatial window centered on the GLASS 250 - meter target pixel corresponding to the Landsat pixel obscured by the target cloud. This spatial window corresponds to 200 × 200 GLASS pixels. For each pixel within the spatial window, the GLASS NDVI time series over 5 years is averaged to obtain the NDVI reference series for this pixel. Calculate the correlation coefficient between the NDVI reference series of the target pixel and the NDVI reference series of other pixels. Pixels with a correlation coefficient higher than 0.8 are determined as similar pixels; Construct a 5 - year time window centered on the date corresponding to the Landsat pixel obscured by the target cloud. The GLASS NDVI image corresponding to the target date is called the target NDVI, and the GLASS NDVI images corresponding to other dates within the time window are called reference NDVI. Calculate the structural similarity index between all reference NDVI and the target NDVI, and sort them from high to low according to the structural similarity index to select the reference image accordingly.

5. The method for efficiently generating the 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales according to claim 1, wherein, Step S4 is specifically as follows: For Landsat NDVI images contaminated only by clouds, first select the reference image with the highest structural similarity index, establish a regression model between the GLASS NDVI values of similar pixels on the reference image and the GLASS NDVI of similar pixels on the target image, and apply this regression model to reconstruct the missing NDVI values on the target image from the Landsat clear - sky 30 - meter NDVI corresponding to the reference image. If the missing NDVI values on the target image are not fully reconstructed, use the reference image with the second - highest structural similarity index to continue the reconstruction until all missing NDVI values on the target image are reconstructed.

6. The method for efficiently generating a 30-meter 8-day spatio-temporal seamless normalized difference vegetation index applicable to large regional scales according to claim 5, characterized in that, When using multiple reference images to reconstruct the missing NDVI values of the same target image, if there are multiple reconstruction results, calculate weights for each reconstruction result based on GLASS NDVI, and weight all reconstruction results to obtain the final spatio - temporally seamless 30 - meter 8 - day Landsat NDVI.

Images