30.8 meter day space-time seamless normalized vegetation difference index efficient generation method suitable for large area scale

Abstract

The invention discloses a 30.8 meter day space-time seamless normalized vegetation difference index efficient generation method suitable for a large area scale, and the method comprises the steps: selecting a representative sample point, and preferably selecting a homogeneity representative sample point; aggregating cloudless 30-meter Landsat surface reflectance data and angle data of representative sample points to 250 meters as model input, extracting 250-meter NDVI data of GLASS as model output, and constructing a GPR regression model to estimate a 30-meter resolution clear sky NDVI value; for a target date cloud shielding Landsat pixel, searching a similar pixel by evaluating the similarity of a time sequence GLASS NDVI, constructing a time window by taking the date where the similar pixel is located as a center, and searching a reference image; establishing a regression model between a similar pixel GLASS NDVI value on the reference image and a similar pixel GLASS NDVI value on the target image, reconstructing a missing NDVI value on the target image by using the regression model, and weighting a plurality of reconstruction results to obtain a space-time seamless 30.8-meter-day Landsat NDVI; the method can provide important data support for regional scale vegetation dynamic monitoring and ecological environment assessment.

Description

Technical Field

[0001] The invention relates to the technical field of remote sensing product production, and in particular to a method for efficiently generating a 30-meter 8-day spatiotemporal seamless normalized vegetation difference index suitable for a large regional scale. Background Art

[0002] High spatiotemporal resolution NDVI products have important application value in the fields of vegetation dynamic monitoring, ecosystem assessment and agricultural management. With the development of remote sensing technology, a large amount of satellite data has provided the basis for the production of NDVI products, but how to efficiently generate 30-meter resolution NDVI time series on a large regional scale is still an important scientific issue. At present, the generation of NDVI products on a large regional scale mainly faces three challenges: first, the amount of Landsat data is huge, and the traditional pixel-by-pixel NDVI calculation and reconstruction methods are inefficient and difficult to meet the needs of operational production; second, the surface heterogeneity is significant, and simple random sampling or uniform sampling strategies are difficult to ensure the representativeness of the samples, affecting the accuracy of the NDVI inversion model; third, the cloud pollution is serious, resulting in the spatiotemporal discontinuity of NDVI observations, making it difficult to form a continuous and stable time series.

[0003] To address the above issues, existing NDVI generation methods mainly include: (1) NDVI generation method based on block processing. The study area is divided into smaller sub-areas and processed separately, but the spatial continuity of NDVI between sub-areas is ignored, resulting in obvious splicing marks in the product; (2) NDVI reconstruction method based on post-classification processing. NDVI reconstruction models are established according to different vegetation types, but the continuous change characteristics of vegetation growth conditions are not considered, making it difficult to accurately characterize the spatial gradient of NDVI; (3) NDVI generation method based on multi-source data fusion. Reconstruction is performed by combining NDVI products of different spatial resolutions, but 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 impact of terrain 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 presence of mixed pixels affects the quality of model training. In terms of model construction, current studies mostly use complex machine learning models, such as deep neural networks, which have high computational overhead and difficult parameter optimization when applied in large areas. Although some simple empirical models have high computational efficiency, their accuracy often cannot meet the requirements. In terms of cloud pollution data reconstruction, existing methods mainly rely on time series interpolation or spatial neighborhood interpolation, lacking the systematic use of multi-phase and multi-scale NDVI data. At the same time, the adaptability of the reconstruction strategy is insufficient, making it difficult to cope with different degrees of cloud pollution. In recent years, although some studies have begun to focus on the efficiency of large-area NDVI product production, most methods are still limited to the simple extension of small-area algorithms, and have failed to fundamentally solve the efficiency bottleneck of large-area processing. In addition, these methods often ignore efficiency while pursuing accuracy, resulting in huge consumption of computing resources and difficulty in realizing business applications.

[0005] Therefore, it is urgent 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 on a large regional scale. This requires systematic optimization of all links from sample selection, model construction to cloud pollution reconstruction to achieve the goals of high efficiency and high accuracy. Summary of the invention

[0006] In order to solve the problems existing in the prior art, the purpose of the present invention is to provide a method for efficiently generating a 30-meter 8-day spatiotemporal seamless normalized vegetation difference index suitable for large regional scales. The present invention can quickly generate long-term series of high-resolution NDVI products, 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: a method for efficiently generating a 30-meter 8-day spatiotemporal seamless normalized vegetation difference index suitable for large regional scales, comprising the following steps:

[0008] S1: Taking terrain factors and vegetation factors as constraints, a constrained sampling method based on representativeness level is used to select representative sample points, and homogeneous representative sample points are further selected based on the surface classification product;

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

[0010] S3: For the cloud-covered Landsat pixels on the target date, a spatial window is constructed with its location as the center, similar pixels are found by evaluating the similarity of the time series GLASS NDVI, a time window is constructed with its date as the center, and a reference image is found by calculating the structural similarity index;

[0011] S4: A regression model was established between the GLASS NDVI values ​​of similar pixels on the reference image and the GLASS NDVI values ​​of similar pixels on the target image. The regression model was applied to reconstruct the missing NDVI values ​​on the target image from the Landsat clear sky 30-m NDVI corresponding to the reference image. Multiple reconstruction results were weighted to obtain the final spatiotemporal seamless 30-m 8-day Landsat NDVI.

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

[0013] Topography and vegetation are selected as environmental factors to reflect the spatial distribution differences of surface characteristics; in terms of topographic factors, the elevation, slope, aspect, and curvature along contour lines in the digital elevation model (DEM) are selected to reflect the control effect of topography; in terms of vegetation factors, leaf area index and vegetation coverage parameters are selected to represent the constraints of vegetation;

[0014] The fuzzy c-means clustering method is used to perform fuzzy clustering on environmental factors to determine the main distribution areas of different surface features. The frequency distribution map of the center position of the environmental factor combination is generated through superposition analysis. The representative level is determined according to the frequency of the center position. The environmental factor combination appearing on each pixel is expressed as an environmental factor combination chain, and the sample points are designed according to the average membership value of the corresponding pixel.

[0015] The 30-meter surface classification product corresponding to the sample point is obtained. If at one sample point, that is, all 30-meter pixels corresponding to 250 meters have a preset proportion or more belonging to the same category, the sample point is considered to be a homogeneous sample point.

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

[0017] Aggregating the cloud-free 30-meter Landsat surface reflectance and angle data of representative sample points involves data from multiple periods. Only when all 30-meter Landsat pixels corresponding to the sample point are marked as clear sky in a certain period, the 30-meter Landsat surface reflectance and angle data are aggregated to 250 meters as model input;

[0018] When constructing the GPR regression model, the ten-fold cross-validation method was used to evaluate the generalization ability of the model, the kernel function was used to construct the covariance matrix, and the log-likelihood function was maximized 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] A 50km×50km spatial window was constructed with the GLASS 250m target pixel corresponding to the target cloud-covered Landsat pixel as the center. The spatial window corresponded to 200×200 GLASS pixels. For each pixel in the spatial window, the GLASS NDVI time series within 5 years was averaged to obtain the NDVI reference sequence of the pixel. The correlation coefficient between the NDVI reference sequence of the target pixel and the NDVI reference sequence of other pixels was calculated. Pixels with a correlation coefficient higher than 0.8 were identified as similar pixels.

[0021] A 5-year time window was constructed with the date corresponding to the target cloud-obscured Landsat pixel as the center. The GLASS NDVI image corresponding to the target date was called the target NDVI, and the GLASS NDVI images corresponding to other dates in the time window were called reference NDVI. The structural similarity index of all reference NDVIs and the target NDVI was calculated, and the images were sorted from high to low according to the structural similarity index to select the benchmark image.

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

[0023] For Landsat NDVI images that are only contaminated by clouds, we first select the benchmark image with the highest structural similarity index, and establish a regression model between the GLASS NDVI values ​​of similar pixels on the benchmark image and the GLASS NDVI values ​​of similar pixels on the target image. The regression model is used to reconstruct the missing NDVI values ​​on the target image from the Landsat clear sky 30-meter NDVI corresponding to the benchmark image. If the missing NDVI values ​​on the target image are not completely reconstructed, the benchmark image with the second highest structural similarity index is used to continue reconstruction until all missing NDVI values ​​on the target image are reconstructed.

[0024] As a further improvement of the present invention, when multiple reference images are used to reconstruct the missing NDVI values ​​of the same target image, if there are multiple reconstruction results, a weight is calculated for each reconstruction result based on the GLASS NDVI, and all reconstruction results are weighted to obtain the final spatiotemporal seamless 30-meter 8-day Landsat NDVI.

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

[0026] 1. The present invention establishes a complete sample selection system based on environmental factors, which significantly improves the representativeness and reliability of samples and provides 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 accuracy, making NDVI estimation on a large regional scale possible.

[0028] 3. The present invention designs a reconstruction strategy based on multiple similarities, which makes full use of spatiotemporal information and effectively solves the problem of data missing 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 regional-scale vegetation dynamic monitoring and ecological environment assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 is a flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

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

[0032] Example 1

[0033] like Figure 1 As shown in the figure, an efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index suitable for large regional scales includes:

[0034] S1: Taking terrain factors and vegetation factors as constraints, a constrained sampling method based on representativeness level is used to select representative sample points, and homogeneous representative sample points are further selected based on the surface classification product;

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

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

[0037] S4: A regression model was established between the GLASS NDVI values ​​of similar pixels on the reference image and the GLASS NDVI values ​​of similar pixels on the target image. The regression model was applied to reconstruct the missing NDVI values ​​on the target image from the Landsat clear sky 30-m NDVI corresponding to the reference image. Multiple reconstruction results were weighted to obtain the final spatiotemporal seamless 30-m 8-day Landsat NDVI.

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

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

[0040] In order to identify representative sample points, the fuzzy c-means clustering method was used to systematically analyze environmental factors. The maximum number of iterations was set to 1000, the fuzzy factor was set to 2, and the convergence threshold was set to 0.001. These parameter settings were verified by repeated experiments and can better reflect the spatial differentiation characteristics of environmental factors. The frequency distribution map of the central position of the environmental factor combination was generated by superposition analysis, the frequency of each position was calculated, and the representative level of different regions was determined according to the frequency. On this basis, the homogeneity evaluation was carried out in combination with the Globa lLand30 surface classification product, and more than 90% of the 30-meter pixels within the 250-meter pixel belonged to the same type of sample points strictly selected. This multi-level screening strategy ensures the representativeness and reliability of the sample, and finally selected about 200,000 high-quality training samples in the study area.

[0041] In one embodiment, step S2 includes:

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

[0043] The GPR model was constructed using a systematic optimization strategy. The generalization ability of the model was evaluated by a ten-fold cross-validation method, the covariance matrix was constructed using a kernel function, and the hyperparameters were optimized by maximizing the log-likelihood function. In order to improve the processing efficiency of large areas, a partitioned parallel processing strategy based on the UTM projection band was adopted to divide the study area into sub-areas of moderate size for simultaneous processing. After the model training was completed, the boundaries of each sub-area were smoothed to ensure the spatial continuity of the final result.

[0044] In one embodiment, step S3 includes:

[0045] Through a large number of experimental verifications, it is determined that it is most appropriate to construct a 50km×50km spatial window centered on the target cloud-covered 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, which not only meets the needs of statistical analysis, but also facilitates parallel processing. For each pixel, GLASS NDVI observations within 5 years are collected, and the NDVI reference series is obtained by averaging the time series. This multi-year averaging method can effectively reduce the impact of outliers and noise.

[0046] In the time dimension, a 5-year time window is constructed with the target date as the center. The basis for choosing 5 years as the time window is that the time span is long enough to include enough potential benchmark images, but not too long to avoid significant changes in surface cover. The GLASS NDVI corresponding to the target date is defined as the target NDVI, and the GLASS NDVI of other dates in the time window are used as candidate benchmark NDVI. The similarity between them is evaluated by calculating the structural similarity index, which takes into account the brightness, contrast and structural characteristics of the image at the same time, and can fully reflect the similarity between images.

[0047] In one embodiment, step S4 includes:

[0048] For each cloud-contaminated Landsat image in the study area, the optimal reference image is first selected based on the sorting results of the structural similarity index. On this reference image, a regression model between the GLASS NDVI value and the GLASS NDVI value of the target image is established using the similar pixels identified in the previous step. This regression relationship based on GLASS NDVI can capture the systematic differences in NDVI values ​​between different phases. Using the established regression model, the missing values ​​on the target image are reconstructed from the Landsat clear sky 30-meter NDVI corresponding to the reference image.

[0049] Considering that a single reference image may not be able to fully cover all cloud-contaminated areas, an iterative reconstruction strategy is adopted. When the optimal reference image cannot be fully reconstructed, the reference images with the highest structural similarity index are used in turn to continue reconstruction until all missing areas are filled. During the reconstruction process, the same location may come from the reconstruction results of multiple reference images. At this time, the weight coefficient of each reconstruction result is calculated 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, Landsat NDVI products with a spatial resolution of 30 meters and a temporal resolution of 8 days are finally generated in the Chinese region from 2000 to 2023.

[0051] This embodiment successfully solves the computational efficiency and accuracy problems of generating high spatiotemporal resolution NDVI products at large regional scales through effective sample selection strategies and efficient model building methods, generates high-precision, spatiotemporally continuous 30-meter resolution NDVI products, and provides important data support for regional-scale vegetation dynamic monitoring and ecological environment assessment. It integrates constrained sampling and efficient reconstruction strategies, fully utilizes the advantages of multi-source remote sensing data and advanced algorithms, considers the spatial heterogeneity and temporal continuity characteristics of NDVI, and realizes the rapid generation of high spatiotemporal resolution NDVI products at large regional scales.

[0052] Example 2

[0053] A method for efficiently generating a 30-meter 8-day spatiotemporal seamless normalized vegetation difference index comprises the following steps:

[0054] Step 1: Taking terrain factors and vegetation factors as constraints, a constrained sampling method based on representativeness level is used to select representative sample points, and homogeneity representative sample points are further selected based on the surface classification product;

[0055] Taking China as an example, we first construct a complete set of environmental factor systems. In terms of terrain factors, basic parameters such as elevation, slope, and aspect are extracted from the 30-meter resolution SRTM DEM, and composite indicators such as terrain relief and curvature along contour lines are calculated. Taking into account the regional differences in China's terrain, different parameter thresholds are used for different geographical divisions: in the plains, the focus is on micro-topography changes, while in the mountainous areas, more attention is paid to the expression of large-scale terrain characteristics. In terms of vegetation factors, in addition to the basic leaf area index and vegetation coverage, auxiliary indicators such as the normalized difference water index (NDWI) and the enhanced vegetation index (EVI) are also introduced. These indicators characterize the growth status and moisture characteristics of vegetation from different angles to improve the representativeness of the samples. For regions with obvious seasonality, such as the Northeast, special attention is also paid to the spatial differentiation of phenological characteristics.

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

[0057] In the homogeneity evaluation stage, a double verification strategy was adopted by combining the latest Global Land 30 surface classification product and China's land use data. 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, it will be determined as the final homogeneity sample point. Through this rigorous screening, about 1 million high-quality training samples were finally selected nationwide.

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

[0059] In the implementation process in China, we first designed data preprocessing and aggregation in a targeted manner. Considering the vast territory of my country, the diverse climate types, and the significant differences in atmospheric conditions and observation geometry, a regional and time-divided processing strategy is adopted for data aggregation. Based on climate and surface characteristics, the country is divided into six major regions: Northeast China, North China, South China, Northwest China, Qinghai-Tibet Plateau, and Southwest China. Each region conducts data aggregation and model training independently.

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

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

[0062] To ensure the generalization performance of the model, a stratified cross-validation strategy was adopted. In the spatial dimension, the samples of each partition were divided into 10 sub-regions according to the latitude and longitude grid; in the temporal dimension, considering the significant seasonal differences in China, the samples were grouped according to the phenological period (rather than simply the month). Through this stratified cross-validation, the model was ensured to have stable performance in different geographical locations and different growth stages.

[0063] Step 3: Establish a similarity analysis method based on space-time windows. For cloud-contaminated areas, 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] In the implementation process in China, the construction of the spatial window needs to consider the complex and diverse geographical environment characteristics. With the target cloud pollution pixel as the center, a dynamic spatial window is set based on the partition strategy: in the eastern plain area, considering the gentle terrain and relatively uniform 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 transition zone of the Qinghai-Tibet Plateau, an irregular spatial window is set according to the altitude gradient to ensure the relative consistency of surface conditions within the search range. For each pixel in the spatial window, its 5-year GLASS NDVI time series is extracted. Considering the vegetation growth characteristics of different climate zones in my country, the processing of the time series adopts a partition strategy: in temperate regions, the focus is on analyzing the dynamic characteristics of NDVI in the growing season; in subtropical and tropical regions, more attention is paid to the interannual variation of NDVI; in arid and semi-arid areas, special attention is paid to 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, the correlation coefficient of the NDVI time series is calculated, and a basic threshold of 0.8 is set. However, in different regions, the threshold is dynamically adjusted according to the degree of vegetation heterogeneity: in the farmland concentration area, due to the influence of agricultural activities, the threshold is appropriately raised to 0.85; in the natural vegetation area, considering the gradual characteristics of community succession, the threshold can be appropriately lowered to 0.75. In addition to the correlation coefficient, the main phenological characteristics of the time series (such as the start and end time of the growing season, peak value, etc.) are also calculated, requiring that the difference in the key phenological period of similar pixels does not exceed 8 days.

[0065] In terms of time window construction, the target date is used as the center and the time window is extended forward and backward for 5 years. The actual utilization strategy of the time window varies according to local conditions: in the northern region, historical data of the same season is used; in the evergreen region in the south, data of different seasons can be used more flexibly. The structural similarity index is calculated for each date image in the time window, which comprehensively considers the brightness, contrast and structural characteristics of the image. The calculation parameters of the structural similarity index are adjusted according to regional characteristics. In plain agricultural areas, more attention is paid to the consistency of local structural characteristics; in forest areas, the weights of brightness and contrast are appropriately increased; in farmland-urban mixed areas around cities, boundary preservation items are introduced to avoid mismatching between different land types. Based on the calculated structural similarity index, all candidate images in the time window are sorted to provide priority reference for subsequent NDVI reconstruction. To improve computational efficiency, the entire similarity analysis process adopts a hierarchical parallel computing strategy. First, a coarse-scale similarity assessment is performed on the GLASS NDVI resolution to screen out potential candidate areas, and then a refined analysis is performed in these areas. At the same time, the boundary information of the geographical division is used to avoid similarity searches across significant geographical boundaries, such as the areas on both sides of the Hengduan Mountains. In actual processing, in order to deal with possible data missing, an alternative plan is also established. When there is insufficient valid data within the 5-year time window (such as areas with poor observation conditions such as mountains and deserts), sufficient reference data is obtained by expanding the time window or adjusting the similarity threshold. At the same time, data quality indicators are recorded to provide a basis for uncertainty assessment for subsequent reconstruction results.

[0066] Step 4: Implement an efficient NDVI reconstruction strategy based on multiple reference images. By establishing a regression relationship between the reference image and the target image, using the structural similarity index to guide the reconstruction process, and using a weighted fusion method based on GLASS NDVI, a high-resolution NDVI product that is continuous in time and space is finally generated.

[0067] The regression model was established using a multi-level strategy. First, the benchmark image with the highest structural similarity index was used. Considering the diversity of vegetation types in China, the surface type factor was introduced when establishing the regression relationship between GLASS NDVI values. For farmland areas, piecewise linear regression was used to capture the stage characteristics of crop growth; for forest areas, polynomial regression was used to better describe continuous changes; for grassland areas, precipitation lag effects were introduced as auxiliary variables. Strict quality control was implemented during the reconstruction process. When the first benchmark image was used for reconstruction, an accuracy threshold based on vegetation type was set: the relative error between the reconstruction result and GLASS NDVI for farmland areas must be less than 10%, 15% for forest areas, and 20% for grassland areas. For areas that failed to meet the accuracy requirements or were not fully reconstructed, an iterative reconstruction process based on structural similarity index sorting was initiated.

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

[0069] The quality assessment of the final product was cross-validated using multiple methods: direct comparison with ground flux station measurement data; up-scale comparison with existing medium-resolution NDVI products such as MODI S; and verification of typical sample areas using high-resolution drone observation data. The verification results show that the overall accuracy of the product has reached R 2 >0.85, more than 95% of the pixels in the temporal continuity evaluation showed reasonable seasonal variation characteristics.

[0070] The above-mentioned embodiments only express the specific implementation of the present invention, and the description thereof is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. An efficient method for generating a 30-meter 8-day spatiotemporal seamless normalized vegetation difference index suitable for large regional scales, characterized in that: The following steps are involved: S1: Taking terrain factors and vegetation factors as constraints, a constrained sampling method based on representativeness level is used to select representative sample points, and homogeneous representative sample points are further selected based on the surface classification product; S2: Aggregate the cloud-free 30-meter Landsat surface reflectance and angle data of representative sample points to 250 meters as model input, extract GLASS 250-meter NDVI data as model output, and build a GPR regression model to estimate the 30-meter resolution clear sky NDVI value from the Landsat surface reflectance and angle data; S3: For the cloud-covered Landsat pixels on the target date, a spatial window is constructed with its location as the center, similar pixels are found by evaluating the similarity of the time series GLASS NDVI, a time window is constructed with its date as the center, and a reference image is found by calculating the structural similarity index; S4: A regression model was established between the GLASS NDVI values ​​of similar pixels on the reference image and the GLASS NDVI values ​​of similar pixels on the target image. The regression model was applied to reconstruct the missing NDVI values ​​on the target image from the Landsat clear sky 30-m NDVI corresponding to the reference image. Multiple reconstruction results were weighted to obtain the final spatiotemporal seamless 30-m 8-day Landsat NDVI.

2. The efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index applicable to large regional scale according to claim 1 is characterized in that: Step S1 is specifically as follows: Topography and vegetation are selected as environmental factors to reflect the spatial distribution differences of surface characteristics; in terms of topographic factors, the elevation, slope, aspect, and curvature along contour lines in the digital elevation model (DEM) are selected to reflect the control effect of topography; in terms of vegetation factors, leaf area index and vegetation coverage parameters are selected to represent the constraints of vegetation; The fuzzy c-means clustering method is used to perform fuzzy clustering on environmental factors to determine the main distribution areas of different surface features. The frequency distribution map of the center position of the environmental factor combination is generated through superposition analysis. The representative level is determined according to the frequency of the center position. The environmental factor combination appearing on each pixel is expressed as an environmental factor combination chain, and the sample points are designed according to the average membership value of the corresponding pixel. The 30-meter surface classification product corresponding to the sample point is obtained. If at one sample point, that is, all 30-meter pixels corresponding to 250 meters have a preset proportion or more belonging to the same category, the sample point is considered to be a homogeneous sample point.

3. The efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index applicable to large regional scale according to claim 1 is characterized in that: Step S2 is specifically as follows: Aggregating the cloud-free 30-meter Landsat surface reflectance and angle data of representative sample points involves data from multiple periods. Only when all 30-meter Landsat pixels corresponding to the sample point are marked as clear sky in a certain period, the 30-meter Landsat surface reflectance and angle data are aggregated to 250 meters as model input; When constructing the GPR regression model, the ten-fold cross-validation method was used to evaluate the generalization ability of the model, the kernel function was used to construct the covariance matrix, and the log-likelihood function was maximized by optimizing the hyperparameters to finally obtain the optimal GPR model parameter combination.

4. The efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index applicable to large regional scale according to claim 1 is characterized in that: Step S3 is as follows: A 50km×50km spatial window was constructed with the GLASS 250m target pixel corresponding to the target cloud-covered Landsat pixel as the center. The spatial window corresponded to 200×200 GLASS pixels. For each pixel in the spatial window, the GLASS NDVI time series within 5 years was averaged to obtain the NDVI reference sequence of the pixel. The correlation coefficient between the NDVI reference sequence of the target pixel and the NDVI reference sequence of other pixels was calculated. Pixels with a correlation coefficient higher than 0.8 were identified as similar pixels. A 5-year time window was constructed with the date corresponding to the target cloud-obscured Landsat pixel as the center. The GLASS NDVI image corresponding to the target date was called the target NDVI, and the GLASS NDVI images corresponding to other dates in the time window were called reference NDVI. The structural similarity index of all reference NDVIs and the target NDVI was calculated, and the images were sorted from high to low according to the structural similarity index to select the benchmark image.

5. The efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index applicable to large regional scale according to claim 1 is characterized in that: Step S4 is specifically as follows: For Landsat NDVI images that are only contaminated by clouds, we first select the benchmark image with the highest structural similarity index, and establish a regression model between the GLASS NDVI values ​​of similar pixels on the benchmark image and the GLASS NDVI values ​​of similar pixels on the target image. The regression model is used to reconstruct the missing NDVI values ​​on the target image from the Landsat clear sky 30-meter NDVI corresponding to the benchmark image. If the missing NDVI values ​​on the target image are not completely reconstructed, the benchmark image with the second highest structural similarity index is used to continue reconstruction until all missing NDVI values ​​on the target image are reconstructed.

6. The efficient generation method of 30-meter 8-day spatiotemporal seamless normalized vegetation difference index applicable to large regional scale according to claim 5 is characterized in that: When multiple reference images are used to reconstruct the missing NDVI values ​​of the same target image, if there are multiple reconstruction results, a weight is calculated for each reconstruction result based on the GLASS NDVI, and all reconstruction results are weighted to obtain the final spatiotemporal seamless 30-meter 8-day Landsat NDVI.

Images