Global total primary productivity simulation method based on EVI

Abstract

The invention provides a global total primary productivity simulation method based on EVI. According to the method, high-spatial-resolution Sentinel-2 satellite data and ERA5-Land meteorological data are obtained, photosynthetically active radiation PAR and enhanced vegetation index EVI data are combined, a general primary productivity GPP simulation model suitable for a global terrestrial ecosystem is constructed, different sensitivities are presented in combination with normalized vegetation index NDVI and EVI data in different natural environments, and the global terrestrial ecosystem can be constructed. According to the method, two abnormal value thresholds are set from two dimensions of NDVI and EVI, EVI data are reconstructed together, GPP is simulated based on the reconstructed EVI and PAR, the accuracy of abnormal data screening is further improved, sudden change values caused by noise or abnormal events are effectively removed, and therefore the data quality is improved, and the precision and spatial resolution of GPP simulation are effectively improved.

Description

Technical Field

[0001] The present invention relates to the field of vegetation remote sensing. By utilizing high-resolution Sentinel-2 satellite data and site observations, a gross primary productivity (GPP) simulation method is constructed and applied to achieve a high spatial resolution GPP simulation method for global terrestrial ecosystems. Background Art

[0002] Gross Primary Production (GPP) describes the amount of atmospheric carbon dioxide fixed by vegetation through photosynthesis over a given period of time. It can reflect the dynamics of terrestrial ecosystems as carbon sinks and is crucial for monitoring the global carbon cycle and climate change. Numerous studies have effectively constructed GPP simulation models and successfully developed regional and global-scale GPP products, enhancing our understanding of vegetation photosynthesis and the global carbon cycle. Currently, there are five main GPP simulation algorithms: the eddy covariance method, which relies on flux towers and lacks spatially continuous observations; process-based models, which involve a large number of parameters and input data and are computationally difficult; data-driven statistical models, which rely on the quality of input data and lack interpretability; simulation methods based on chlorophyll fluorescence, which generally have low spatial resolution; and semi-empirical light energy utilization models, which are relatively simple in structure but also need to consider factors such as different ecosystem types.

[0003] Simulating gross primary productivity (GPP) typically requires multiple parameters, including meteorological, vegetation, soil, and other environmental parameters. Current GPP simulation models primarily use one or more of the following parameters: photosynthetically active radiation (PAR), leaf area index (LAI), normalized difference vegetation index (NDVI), and enhanced vegetation index (EVI). Different parameter choices are generally applicable to specific vegetation types and can affect the accuracy and effectiveness of GPP simulation models.

[0004] The current mainstream satellite remote sensing data used in the simulation process include AVHRR NDVI, MODIS NDVI / EVI, and GOME-2SIF. AVHRR NDVI data, acquired using the Advanced Very High Resolution Radiometer (AVHRR) sensor on NOAA satellites, is a normalized vegetation index (NDVI) widely used in vegetation monitoring and ecological and environmental research. The AVHRR sensor has a native spatial resolution of 1.1 km (corresponding to the instantaneous field of view at satellite altitude). However, after processing and resampling, NDVI data often has spatial resolutions of 1 km and 0.5 km, which vary between datasets and processing versions. MODIS data is primarily used to extract EVI and NDVI data, with a spatial resolution of 0.25 km to 1 km. GOME-2 SIF (Solar Induced Chlorophyll Fluorescence) primarily extracts photosynthetically active radiation (PAR) and chlorophyll fluorescence (SIF) from vegetation. SIF is closely related to vegetation photosynthesis and can be used to simulate vegetation primary productivity (GPP). By analyzing SIF changes, we can understand the photosynthetic capacity and health of vegetation. Its spatial resolution is approximately 80 km. Therefore, the spatial resolution provided by current mainstream data is relatively coarse, with accuracy typically ranging from 200 meters to tens of kilometers.

[0005] In summary, the current GPP simulation algorithms have various deficiencies. Due to reasons of data acquisition and parameter selection, existing GPP simulation products are difficult to simultaneously meet the requirements of high spatial resolution and global coverage, and cannot support fine-scale ecological research.

[0006] The European Space Agency's Sentinel-2 satellite provides 5-day revisits, multiple spatial resolutions of 10-60 meters, and observation content from visible light to shortwave infrared. Compared with data such as AVHRR NDVI, MODIS NDVI / EVI and GOME-2SIF, it has an advantage in improving spatial resolution and can clearly display surface details. It is suitable for fine-scale ecological research and environmental monitoring, becoming a new tool for detecting landscape heterogeneity and ecosystem status, and is expected to solve the current problem of coarse spatial resolution.

[0007] In summary, there is currently a lack of accurate and effective global high-spatial-resolution GPP simulation models. There is an urgent need to simulate global terrestrial ecosystems with the help of high-spatial-resolution satellite data and select better parameters to set up GPP simulation models to further enhance our understanding of the global carbon cycle and climate change. Summary of the Invention

[0008] To solve the above problems, this application provides a global gross primary productivity simulation method based on EVI. Through Sentinel-2 satellite data, high-spatial-resolution vegetation index information is obtained. Combined with factors such as parameter sensitivity and photosynthesis feedback, a GPP model suitable for simulating global high-spatial-resolution is constructed to achieve global-scale terrestrial ecosystem GPP simulation.

[0009] Simulating gross primary productivity (GPP) usually requires multiple parameters. In the parameter selection of the model, different parameters reflect different conditions of vegetation. The more parameters, the more accurate and detailed the response. However, too many parameters will lead to increased computational complexity, which is not conducive to practical application.

[0010] Considering that photosynthetically active radiation (PAR) is the energy driver of vegetation photosynthesis, it directly determines the rate of photosynthesis and plant productivity. It can better reflect the vegetation's conversion of light energy into chemical energy through photosynthesis, thereby fixing carbon. PAR is also closely related to vegetation photosynthetic efficiency, ecosystem carbon cycling, and vegetation growth status. It is an important foundational indicator for vegetation productivity models and ecosystem carbon balance research. PAR data is relatively easy to obtain and monitor, with multiple observation methods and data sources available globally, facilitating large-scale application. Therefore, PAR is selected as one of the key drivers of GPP simulations.

[0011] However, PAR only provides the energy input for photosynthesis and does not take into account the situation where vegetation is subjected to growth stress such as drought, nutrient deficiency, and environmental changes. These stress factors will affect the photosynthetic efficiency and productivity of vegetation. Therefore, there are many uncertainties in estimating GPP based solely on PAR, and it is necessary to conduct comprehensive analysis and correction in combination with other factors to improve the accuracy and reliability of GPP simulation.

[0012] Currently, the most commonly used parameter for simulating vegetation growth, health, and biomass is the Normalized Difference Vegetation Index (NDVI). This is based on vegetation's high reflectance in the near-infrared band and low reflectance in the red band. Higher NDVI values ​​indicate more lush vegetation. Compared to the NDVI, the Enhanced Vegetation Index (EVI) incorporates blue wavelengths to correct for atmospheric influences and soil background noise, thereby more accurately reflecting the actual growth, health, and biomass of vegetation. It is more sensitive to vegetation health and photosynthetic activity. The EVI is applicable to various ecosystem types, including forests, grasslands, farmlands, and wetlands. Whether in forest ecosystems with high vegetation cover or sparse grasslands or farmlands with low vegetation cover, the EVI effectively reflects vegetation growth and photosynthetic capacity, providing a reliable basis for GPP simulations. Therefore, a combination of photosynthetically active radiation (PAR) and EVI data is chosen to construct a GPP simulation model to better reflect the actual photosynthetic activity of vegetation.

[0013] However, due to the influence of clouds, aerosols, etc., vegetation indices usually have outliers. For example, compared with NDVI and EVI, NDVI is more sensitive to precipitation and radiation factors, while EVI is more affected by terrain factors (especially in areas with undulating terrain) than NDVI. Therefore, under different natural environments, the two data show different sensitivities. In order to improve the quality of the data and reduce the impact of unreliable outliers caused by cloud cover, atmospheric effects, precipitation, terrain and other factors, this application selects these two vegetation indices to complement each other, sets two outlier thresholds from different dimensions, and jointly reconstructs the EVI data.

[0014] Specifically, the GPP simulation model is as follows: ; in, represents the GPP in the case of CX vegetation, specifically, Represents the reconstructed enhanced vegetation index EVI, as well as The coefficients corresponding to different vegetation conditions are respectively, specifically, the corresponding neural network model can be trained by a sample set of multiple experimental data under multiple different vegetation conditions to obtain the corresponding optimal coefficients.

[0015] The present invention is achieved through the following technical solutions: Step S1, acquiring and preprocessing corresponding data; the data includes Sentinel-2L2A multi-band data and ERA5-Land incident shortwave radiation SSRD data; Specifically: Step S11: Acquire Sentinel-2 L2A multi-band data and ERA5-Land incident shortwave radiation (Surface Solar Radiation Downward, SSRD) data from Google Earth Engine (GEE); Step S12: Extract the QA60 band of the Sentinel-2 L2A multi-band data as a quality control band, read the cloud and cirrus flag bits therein through bitwise operation, and treat it as cloud removal processing; Step S2, calculating the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) using Sentinel-2 data; Step S3, by counting the abnormal characteristics of the NDVI time series within N days, the corresponding EVI time series is reconstructed for the first time based on the NDVI time series; Statistical analysis was performed on the NDVI time series, and the corresponding first outlier threshold was set based on its statistical characteristics. The data in the NDVI time series below the first outlier threshold was regarded as anomaly, and the corresponding data in the EVI time series was marked as outliers and removed. The EVI time series after the first reconstruction was obtained by linear interpolation of adjacent valid data.

[0016] Step S31: For each pixel Perform statistics on time series and set outlier thresholds based on their statistical characteristics : ; ; in, , Represented respectively The 65th and 35th percentiles of the time series are subtracted to obtain a threshold , and compare its size with 0.4, and select the minimum value as the threshold for the first outlier screening.

[0017] Step S32: traverse the NDVI time series of all pixels and take the NDVI time series of the same pixel with an increase greater than the first outlier threshold within N days. The value of is regarded as an outlier, where N is the number of days in a period of time, and the EVI value in the EVI time series that is on the same day as the outlier in the NDVI time series is marked as an outlier and removed, and the outlier is reconstructed.

[0018] Step S4: reconstructing the EVI time series after the first reconstruction for the second time by counting the abnormal features of the EVI time series after the first reconstruction; Further, the EVI time series obtained from the first reconstruction is smoothed, and the data obtained after smoothing is statistically analyzed. According to its statistical characteristics, a second outlier threshold is obtained. , and the EVI time series obtained from the first reconstruction is subjected to secondary outlier detection and reconstruction.

[0019] Considering that although the above steps correct based on the abnormal data, some outliers may not be detected. Further, the second outlier threshold of the EVI time series after the first reconstruction for each pixel is calculated : ; ; where represents the EVI time series after the first reconstruction, represents a new series filtered by a Savitzky-Golay filter (m = 5, d = 2, m is the half-width of the smoothing window, d is the degree of the polynomial); and respectively represent the 85th and 15th percentiles of the time series, is the second outlier threshold.

[0020] Points greater than are regarded as outliers, and the corresponding data values are replaced by linear interpolation between adjacent two points, and the time series that has undergone two outlier reconstructions is obtained. time series.

[0021] Step S5: Using the SG iterative smoothing algorithm, reconstruct the EVI time series after the secondary reconstruction, calculate the mean value every M days, where 1 < M < N, for example, M takes the value of 16, and obtain the EVI time series data of the third reconstruction with a time resolution of 16 days.

[0022] Step S6: Resample the obtained SSRD data of ERA5-Land, calculate the mean value every M days, and multiply by the conversion coefficient to obtain the photosynthetically active radiation PAR.

[0023] Step S7: Based on the EVI time series data and PAR of the third reconstruction, simulate the gross primary productivity GPP: ; where represents the GPP in the case of CX vegetation. Specifically, represents the enhanced vegetation index (EVI) after three reconstructions, as well as The coefficients corresponding to different vegetation conditions are as follows. Specifically, the corresponding neural network model can be trained by multiple sample sets of experimental data under different vegetation conditions to obtain the corresponding optimal coefficients and obtain the corresponding Simulation method.

[0024] Compared with the prior art, the advantages of the above method provided by the present invention include: 1. This application optimizes the GPP simulation model. By considering the performance of multiple parameters, it selects photosynthetically active radiation PAR and enhanced vegetation index EVI data to construct a GPP simulation model to better reflect the actual situation of vegetation photosynthesis. Among them, photosynthetically active radiation PAR directly determines the rate of photosynthesis and plant productivity, and PAR data usually has a high temporal and spatial resolution, which is convenient for application on a large scale. Compared with NDVI, the enhanced vegetation index EVI introduces a blue light band to correct for atmospheric effects and soil background noise, thereby more accurately reflecting the actual growth, health status and biomass of vegetation. EVI is applicable to a variety of terrestrial ecosystems around the world, including forests, grasslands, farmlands and wetlands. The two work together to provide a reliable basis for GPP simulation from multiple dimensions, better reflect the actual situation of vegetation photosynthesis, and provide a GPP simulation model with high spatial resolution applicable to a variety of ecosystems around the world.

[0025] 2. This application determines a first outlier threshold by traversing the NDVI time series of each pixel and combining it with the statistical characteristics of the NDVI time series (such as percentiles). The values ​​in the NVDI time series that exceed this threshold are marked as outliers. The data in the corresponding EVI time series is determined to be abnormal based on the outliers in the NVDI time series, and the EVI abnormal data is reconstructed for the first time. On this basis, the EVI time series after the first correction is further screened for secondary outliers. The second outlier threshold is determined by combining the Savitzky-Golay filter and the statistical characteristics of the data time series (such as percentiles). The second outlier threshold is used to further screen possible abnormal data values ​​and perform a secondary reconstruction. Compared with the traditional method of setting outlier thresholds for EVI time series data based on experience to eliminate abnormal data, this application takes into account that NDVI time series and EVI time series data show different sensitivities under different natural environments, and selects these two vegetation indices to complement each other. Two outlier thresholds are set from the two dimensions of NDVI time series and EVI time series to jointly reconstruct the EVI time series. When setting the threshold, the numerical statistical characteristics of each time series (such as percentiles) are taken into account to further improve the accuracy of abnormal data screening. While maintaining the overall trend of the time series, the mutation values ​​caused by noise or abnormal events can be effectively removed, thereby improving data quality and providing more reliable input for subsequent vegetation index analysis and application.

[0026] 3. In the selection of satellite remote sensing data, the observation data with multiple spatial resolutions of 10-60 meters, from visible light to short-wave infrared, provided by the European Space Agency's Sentinel-2 satellite with high spatial resolution are selected. Compared with the existing mainstream AVHRR NDVI, MODIS NDVI / EVI and GOME-2SIF data, it has great advantages in improving spatial resolution, can clearly display surface details, and is suitable for fine-scale ecological research and environmental monitoring. It can achieve 10-meter global land GPP simulation. Combined with the GPP simulation method based on secondary reconstruction EVI of this application, a GPP simulation method with high spatial resolution applicable to various ecosystems around the world is provided, which solves the problem of the existing GPP simulation method being applicable to a small range of ecosystems and coarse spatial resolution, and enhances the monitoring effect of the global carbon cycle and climate change. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following briefly introduces the drawings required for describing the embodiments of the present application. Obviously, the drawings described below are only some embodiments of the present application. Those skilled in the art can also derive other drawings based on these drawings without inventive effort. Figure 1 This is a schematic diagram of the GPP simulation algorithm of this application; Figure 2 It is an experimental data diagram of the embodiment of the present application. DETAILED DESCRIPTION

[0028] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are part of the embodiments of this application, not all of them. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0029] Please refer to Figure 1 , Figure 1 This is a schematic diagram of the GPP simulation algorithm proposed in Example 1 of this application. Figure 1 As shown, the present invention is implemented through the following technical solutions: Step S1, performing corresponding data acquisition and data preprocessing; the data includes Sentinel-2L2A multi-band data and ERA5-Land incident shortwave radiation SSRD data.

[0030] Specifically: Step S11: Acquire Sentinel-2 L2A multi-band data and ERA5-Land incident shortwave radiation (Surface Solar Radiation Downward, SSRD) data from Google Earth Engine (GEE); Table 1 shows the corresponding data information:

[0031] Table 1 Sentinel-2 L2A multi-band data and ERA5-Land incident shortwave radiation (SSRD) data.

[0032] Step S12: Extract the QA60 band of the Sentinel-2 L2A multi-band data as a quality control band, read the cloud and cirrus flag bits therein through bitwise operation, and treat it as cloud removal processing; Specifically: ; ; The above formula is used to perform bitwise operations on 1 to generate a bit mask. The two bit masks are then used to perform a bitwise AND (bitwiseAnd()) operation on QA60. A value of 0 indicates that there is no cloud or cirrus cloud at that pixel.

[0033] After obtaining the data, the GPP simulation algorithm is constructed. Specifically: Step S2, calculating the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) using Sentinel-2 data; Specifically: ; ; in, , , These are the blue, red, and near-infrared reflectance parameters from Sentinel-2 satellite data. By directly extracting the corresponding bands in GEE, the corresponding NDVI and EVI time series are calculated using the above formulas.

[0034] Step S3, by counting the abnormal characteristics of the NDVI time series within N days, the corresponding EVI time series is reconstructed for the first time based on the NDVI time series; Furthermore, a statistical analysis is performed on the NDVI time series, and the corresponding first outlier threshold is set based on its statistical characteristics. The data in the NDVI time series that is lower than the first outlier threshold is regarded as anomaly, and the corresponding data in the EVI time series is marked as anomaly and eliminated. The EVI time series after the first reconstruction is obtained by linear interpolation replacement of adjacent valid data.

[0035] Step S31: Perform statistics on the NDVI time series of each pixel and set the first outlier threshold according to its statistical characteristics : ; ; in, , Represented respectively The 65th and 35th percentiles of the time series are extracted using the GEE built-in function ee.Reducer.percentile(), and the difference is taken to obtain a threshold. , and compare its size with 0.4, and select the minimum value as the threshold for the first outlier screening.

[0036] The value of 0.4 is an empirical value. Analyzing a large amount of vegetation index data reveals that most valid vegetation index data are concentrated within a certain range, while values ​​outside this range are often considered outliers. The 0.4 threshold was determined based on experience and analysis of data distribution to distinguish between normal vegetation growth and possible anomalies. An increase in EVI or NDVI values ​​exceeding 0.4 is considered a possible indicator of data anomalies. Therefore, 0.4 is set as a comparison threshold to compare with the statistical characteristics of the NDVI time series, allowing for optimal identification of possible anomalous vegetation index data.

[0037] Step S32: traverse the NDVI time series of all pixels and take the NDVI time series of the same pixel with an increase greater than the first outlier threshold within N days. The value of is regarded as an outlier, where N is the number of days in a period of time, and the EVI value in the EVI time series that is on the same day as the outlier in the NDVI time series is marked as an outlier and removed, and the outlier is reconstructed.

[0038] Specifically, select the NDVI time series of the same pixel over a period of time: = }, where N is the number of days in a period of time. For example, N=20, which represents the NDVI time series within 20 days, and the data is screened for outliers.

[0039] Calculate the above NDVI time series The average value of For example, when N=20, any With the average The absolute value of the difference between the first outlier threshold Compare, if it is greater than Consider it as abnormal and mark it as an outlier. Further, the time series of the pixel within 20 days Perform difference operation on any two values ​​​​excluding abnormal values, take the absolute value of the difference obtained for statistics, and select the one with a value greater than The data is recorded as the preliminary outlier, and the preliminary outlier is analyzed. The preliminary outlier is compared with the data on both sides of its adjacent data, and the nearest normal value (that is, non-outlier and non-preliminary outlier) is selected from the data on both sides for difference comparison. If the absolute value of the difference between it and the adjacent two sides is not greater than , it is marked as a normal value, otherwise it is marked as an abnormal value. If the preliminary abnormal value is the endpoint data, the nearest normal value (that is, non-abnormal value and non-preliminary abnormal value) in the data on one side of the time series is compared with it. If the absolute value of the difference between it and the nearest normal value is not greater than , then it is marked as a normal value, otherwise it is marked as an abnormal value.

[0040] Filter out the corresponding EVI time series: = }, where N is the number of days in a period of time. For example, N=20, which represents the EVI time series within 20 days. The NDVI time series Outliers on the same day Values ​​are marked as outliers and removed based on the adjacent The values ​​are re-interpolated linearly: ; in , , and Represents the re-interpolated values ​​at the nth time point , the original value at the n-1th time point and the original at time point n+1 , and the re-linear difference is obtained Marked as normal value. So far, the new EVI time series after the first correction is obtained and named .

[0041] Furthermore, if any one or both of the adjacent EVI values ​​on both sides of the abnormal EVI value are also considered abnormal values, the adjacent normal values ​​on the left and right sides closest to it are selected to re-interpolate and perform linear interpolation. The value is an endpoint value (n is 1 or N), or the exception If there is no normal value on one side of the data value, the nearest normal value on the side with a normal value is selected as the replacement value.

[0042] Step S4: reconstruct the EVI time series after the first reconstruction for the second time by counting the abnormal features of the EVI time series after the first reconstruction; Considering that although the above steps correct the EVI based on the abnormal data of NDVI, some abnormal values ​​may not be detected. Further, the EVI time series obtained by the first reconstruction is smoothed and the data obtained after smoothing is statistically analyzed to obtain the second abnormal value threshold according to its statistical characteristics. , perform secondary outlier detection and reconstruction on the EVI time series obtained by the primary reconstruction.

[0043] Furthermore, for each pixel, Outlier threshold for time series: ; ; in, represents the EVI time series after the first reconstruction, represents a new sequence obtained by filtering through a Savitzky-Golay filter (m=5, d=2, m is the half-width of the smoothing window, d is the degree of the polynomial); and Respectively The 85th and 15th percentiles of the time series, is the second outlier threshold.

[0044] Will Greater than The point is considered as abnormal, and the corresponding The data value is replaced by linear interpolation between two adjacent points: ; in, , , and Represents the re-interpolated values ​​at the nth time point , at the n-1th time point and the n+1th time point , and replace the linear difference with Marked as normal value.

[0045] Furthermore, if the exception Adjacent values ​​on both sides If any one or two values ​​are also considered as outliers, the nearest normal values ​​on the left and right sides are selected to re-interpolate the values. The value is an endpoint value (n is 1 or N), or the exception If there is no normal value on one side of the data value, the nearest normal value on the side with a normal value is selected as the replacement value.

[0046] So far, we have obtained the Time series.

[0047] Step S5: Using SG iterative smoothing algorithm, the secondary reconstruction Reconstruct the time series and calculate the mean value every M days, where 1 < M < N. For example, when M is 16, time series data with a time resolution of 16 days is obtained. .

[0048] Step S6: Resample the obtained SSRD data of ERA5-Land, calculate the mean value every M days, and multiply by the conversion coefficient to obtain the photosynthetic active radiation PAR.

[0049] Specifically, it can be resampled to a resolution of 10 meters, and the conversion coefficient can be taken as 0.45, and M can be taken as 16 to obtain the photosynthetic active radiation PAR with a resolution of 16 days: ; Among them, the conversion coefficient 0.45 is obtained from long-term ground observations and measured data of the above data. Through synchronous measurements under different vegetation types and environmental conditions, at different locations and times, the ratio of PAR to SSRD is approximately between 0.4 and 0.5, and 0.45 is selected as the fixed conversion coefficient through approximate processing.

[0050] Step S7: Based on and PAR, simulate the gross primary productivity GPP: ; Among them, represents the GPP in the case of CX vegetation. Specifically, represents the enhanced vegetation index EVI after three reconstructions, and correspond to the coefficients in different vegetation cases respectively. Specifically, through multiple cases of different vegetation, regression analysis is performed based on the corresponding photosynthetic utilization efficiency models to obtain the corresponding coefficients and the corresponding simulation method.

[0051] Figure 2 is the experimental data shown in the embodiments of the present disclosure. Specifically, the corresponding simulation models in C3 vegetation and C4 vegetation are obtained through regression analysis of the corresponding photosynthetic utilization efficiency models: ; ; Among them, represents the GPP of C3 vegetation, represents the GPP of C4 vegetation.

[0052] Figure 2For the comparison of experimental data for the above vegetation, small figures (A1)-(A3) are experimental data in the original time scale, small figures (B1)-(B3) are experimental data after 8 days, and small figures (C1)-(C3) are experimental data after 16 days.

[0053] Among them, the small pictures (B1)-(B3) and (C1)-(C3) are smoothed by SG algorithm compared with the small picture (A1)-(A3). It can be seen that the determination coefficient in the small picture (C3) is The results show that the smoothing process can effectively reduce the noise in the data and improve the fit of the model.

[0054] In summary, the present invention solves the problems of coarse resolution and small scope of application in GPP simulation, and provides a GPP simulation method that can be used for 10-meter spatial resolution and is applicable to various terrestrial ecosystems on a global scale.

[0055] The various embodiments in this specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the various embodiments can be referenced to each other.

[0056] Although preferred embodiments of the present invention have been described, those skilled in the art may make additional changes and modifications to these embodiments once they become aware of the basic inventive concepts. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the embodiments of the present invention.

[0057] Finally, it should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or terminal device that includes a series of elements includes not only those elements, but also other elements not explicitly listed, or elements inherent to such process, method, article, or terminal device. In the absence of further limitations, an element defined by the phrase "comprising a ..." does not exclude the presence of additional identical elements in the process, method, article, or terminal device that includes the element.

Claims

1. A method for simulating global gross primary productivity based on EVI, characterized by: The steps include: Step S1, acquiring and preprocessing corresponding data, wherein the data includes Sentinel-2 L2A multi-band data and ERA5-Land incident shortwave radiation SSRD data; Step S2, calculating the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) using Sentinel-2 data; Step S3, by counting the abnormal characteristics of the NDVI time series within N days, the corresponding EVI time series is reconstructed for the first time based on the NDVI time series; Step S4: reconstruct the EVI time series after the first reconstruction for the second time by counting the abnormal features of the EVI time series after the first reconstruction; Step S5: Use the SG iterative smoothing algorithm to reconstruct the EVI time series after secondary reconstruction and calculate the mean value of each M days, where 1 <M<N; Step S6: Resample the ERA5-Land SSRD data, calculate the mean value for each M days, and multiply it by the conversion coefficient. , get photosynthetically active radiation PAR; Step S7: Based on the three-time reconstructed EVI time series data and PAR, simulate the gross primary productivity (GPP).

2. A global gross primary productivity simulation method based on EVI according to claim 1, characterized in that: The corresponding data acquisition and data preprocessing further include the following steps: Step S11: Obtain Sentinel-2 L2A multi-band data and ERA5-Land incident shortwave radiation SSRD data from Google Earth Engine (GEE); Step S12: Extract the QA60 band of the Sentinel-2 L2A multi-band data as a quality control band, and read the cloud and cirrus flag bits therein through bitwise operation, which is regarded as cloud removal processing.

3. The method for simulating global gross primary productivity based on EVI according to claim 1, characterized in that: The first reconstruction of the corresponding EVI time series based on the NDVI time series by counting abnormal characteristics of the NDVI time series within N days also includes: Statistical analysis was performed on the NDVI time series, and the corresponding first outlier threshold was set based on its statistical characteristics. The data in the NDVI time series below the first outlier threshold was regarded as anomaly, and the corresponding data in the EVI time series was marked as outliers and removed. The EVI time series after the first reconstruction was obtained by linear interpolation of adjacent valid data.

4. The method for simulating global gross primary productivity based on EVI according to claim 3, characterized in that: Obtaining the first reconstructed EVI time series includes the following steps: Step S31: Perform statistics on the NDVI time series of each pixel and set the first outlier threshold according to its statistical characteristics ; ； ； in, , Represent the 65th and 35th percentiles of the NDVI time series, respectively. Difference is taken to get a threshold value. , is the first outlier threshold; Step S32: traverse the NDVI time series of all pixels and take the NDVI time series of the same pixel with an increase greater than the first outlier threshold within N days. The value of is regarded as an outlier, where N is the number of days in a period of time, and the EVI value in the EVI time series that is on the same day as the outlier in the NDVI time series is marked as an outlier and removed, and the outlier is reconstructed.

5. The method for simulating global gross primary productivity based on EVI according to claim 1, characterized in that: The second reconstruction of the EVI time series after the first reconstruction by counting abnormal features of the EVI time series after the first reconstruction further includes: The second outlier threshold is obtained by smoothing the EVI time series obtained by the first reconstruction and performing statistics on the data obtained after smoothing. , perform secondary outlier detection and reconstruction on the EVI time series obtained by the primary reconstruction.

6. The method for simulating global gross primary productivity based on EVI according to claim 5, characterized in that: Performing secondary outlier detection and reconstruction on the EVI time series obtained by the first reconstruction further includes: Calculate the second outlier threshold of the EVI time series after the first reconstruction of each pixel ; ； ； in, represents the EVI time series after the first reconstruction, represents a new sequence obtained by filtering through a Savitzky-Golay filter (m=5, d=2, m is the half-width of the smoothing window, d is the degree of the polynomial); and Respectively The 85th and 15th percentiles of the time series, is the second outlier threshold.

7. The method for simulating global gross primary productivity based on EVI according to claim 6, characterized in that: Performing secondary outlier detection and reconstruction on the EVI time series obtained by the primary reconstruction further includes: Will Greater than The point is considered as abnormal, and the corresponding The data value is replaced by linear interpolation between two adjacent points, and the data value is reconstructed twice by the outlier value. Time series.

8. The method for simulating global gross primary productivity based on EVI according to claim 1, characterized in that: The SG iterative smoothing algorithm is used to reconstruct the EVI time series after the second reconstruction, and the calculation of the mean value of each M days also includes: where N is 20 and M is 16, and the EVI time series data of the third reconstruction with a time resolution of 16 days is obtained.

9. The method for simulating global gross primary productivity based on EVI according to claim 1, characterized in that: Resample the ERA5-Land SSRD data, calculate the mean value every M days, and multiply it by the conversion coefficient , the photosynthetically active radiation PAR also includes: resampling to 10-meter resolution, conversion coefficient It can be taken as 0.45, and M can be taken as 16 to obtain photosynthetically active radiation PAR with a resolution of 16 days.

10. The method for simulating global gross primary productivity based on EVI according to claim 1, characterized in that: Based on the three-time reconstructed EVI time series data and PAR, the simulated gross primary productivity (GPP) also includes: ； in, represents the GPP in the case of CX vegetation, specifically, Represents the enhanced vegetation index EVI after three reconstructions, as well as The coefficients corresponding to different vegetation conditions.

Images