A multi-source remote sensing ecosystem carbon flux estimation method based on phenology

Abstract

The present invention discloses a multi-source remote sensing ecosystem carbon flux estimation method based on phenology, comprising the following steps: using multi-source remote sensing data to calculate the direct and diffuse photosynthetically active radiation of the target area; obtaining regional biological characteristic data of the same time series, using double S-curves to perform nonlinear least squares fitting to extract vegetation phenological characteristics, setting thresholds according to phenological characteristics, and constructing phenological constraint factors; based on the two-leaf light energy utilization efficiency framework, combining phenological constraint factors, establishing an ecosystem carbon flux model; and optimizing the model parameters using the Markov chain Monte Carlo method. Using the model after parameter optimization, global gross primary productivity data is generated. The method of the present invention is superior to the traditional light energy utilization efficiency model, and has advantages in simulating the temporal and spatial changes of gross primary productivity, and improves the simulation accuracy of global carbon flux.

Description

Technical Field

[0001] The present invention belongs to the field of carbon flux estimation of earth system spheres, and specifically relates to a multi-source remote sensing ecosystem carbon flux estimation method based on phenology. Background Art

[0002] Vegetation productivity is the beginning of the ecosystem carbon cycle and one of the key issues in the global carbon cycle. The gross primary productivity (GPP) of vegetation is defined as the total amount of organic carbon fixed by vegetation through photosynthesis per unit space and per unit time.

[0003] There are two common methods for measuring GPP in terrestrial ecosystems. Accurate measurements of GPP are usually made at the laboratory scale, and the development of eddy covariance technology has enabled local GPP measurements. However, it is complicated to scale up ecosystem processes from the laboratory scale or local scale to the regional or even global scale due to the influence of multiple environments on the flux. The development of quantitative remote sensing technology and computer science has provided a basis for monitoring the spatiotemporal changes of GPP on a global scale. Current GPP models can be divided into physical process-based, machine learning, and light energy utilization efficiency models. Process-based land-atmosphere models explicitly describe the physical processes involved in the carbon, water, and energy cycles, such as products such as BESS and TRENDY. These models have a good ability to monitor the interannual changes of GPP. However, this type of model usually requires a large number of input variables and parameter data sets, some of which may be difficult to obtain or have limited accuracy. Machine learning models are data-driven models based on mathematical and statistical principles. EC (eddy covariance) fluxes measured at the field scale are scaled up to regional, national, and global scales, such as products such as FluxCom and MTE_GPP, based on the unique patterns found in the input data. These models have demonstrated the effectiveness of solving complex classification and estimation problems through nonlinear regression fitting. However, the data mining process is often implemented through black-box models, which may limit the analysis and interpretation of data trends. The light use efficiency model follows the basic principles of biogeochemistry and considers light use efficiency to estimate GPP, assuming that the light use efficiency of plants may be reduced below its theoretical potential value due to environmental stress such as low temperature or water shortage. Currently, the light use efficiency model has become a very popular global GPP estimation method due to its high computational efficiency and few parameters, such as MODIS, GLASS and VPM.

[0004] At present, there is no two-leaf light energy utilization model that takes phenological factors into consideration. Photosynthesis is closely related to the phenological status of vegetation, and phenology is one of the most important controlling variables of the interannual variation of total primary productivity. The light energy utilization efficiency of vegetation is significantly affected by internal physiological traits such as carotenoid content and chlorophyll content. Studies have shown that the maximum light energy utilization efficiency of vegetation can significantly distinguish between green vegetation and senescent vegetation. The total amount of chlorophyll in the canopy during the growth stage of vegetation may be more than twice that of the senescent stage. Parameterizing the maximum light energy utilization efficiency of vegetation leaves in the mature and senescent stages, respectively, can significantly improve the accuracy of the light energy utilization efficiency model. Summary of the invention

[0005] In view of the problems existing in the prior art, the present invention provides a phenological-based multi-source remote sensing ecosystem carbon flux estimation method, which is based on a two-leaf light energy utilization efficiency framework and can verify its advantages in simulating the spatiotemporal changes of GPP.

[0006] To solve the above technical problems, the present invention provides the following technical solution: a multi-source remote sensing ecosystem carbon flux estimation method based on phenology, comprising the following steps:

[0007] S1. Based on multi-source satellite remote sensing data and meteorological data, the satellite instantaneous shortwave radiation is calculated, and then the daily solar radiation flux is calculated; according to the solar altitude angle and the clear sky index, the proportion of scattered photosynthetically active radiation is calculated, and then the high temporal and spatial resolution photosynthetically active radiation data is inverted, including scattered photosynthetically active radiation, direct photosynthetically active radiation, and scattered photosynthetically active radiation under the canopy.

[0008] S2. Collect biological characteristic data of the same time series in the target area, calculate vegetation phenological characteristics, and provide state variables for the establishment of phenological constraint factors;

[0009] S3, according to step S1, obtaining photosynthetically active radiation data with high temporal and spatial resolution, calculating the photosynthetically active radiation absorbed by the sun leaves and shade leaves of the vegetation, defining thresholds according to the state variables of step S2, establishing phenological constraint factors based on the surface water index, and then calculating and obtaining temperature constraint factors, humidity constraint factors, and carbon dioxide constraint factors, and further constructing an environmental constraint function; comprising the following sub-steps:

[0010] S3.1. Calculate the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation APAR su , APAR sh , as follows:

[0011] ,

[0012] ,

[0013] In the formula, LAI su and LAI sh are the leaf area index of sun leaves and shade leaves respectively, LAI is the leaf area index, is the clustering index, is the mean leaf sun angle, is the solar altitude angle, C is the multiple scattering contribution of direct radiation in the vegetation canopy, PAR dif It is diffuse photosynthetically active radiation. PAR dir Direct photosynthetically active radiation PAR dif，u is the diffuse photosynthetically active radiation below the canopy;

[0014] Area index of sun leaf LAI su and leaf area index of shade leaves LAI sh , calculated as follows:

[0015] ,

[0016] ,

[0017] S3.2. Establish phenological constraint factors based on surface water index as follows:

[0018] ,

[0019] ,

[0020] Where LSWI is the vegetation moisture index, and Represent the near infrared and short wave infrared bands of satellite sensors, SIF se and SIF ma are the SIF values ​​corresponding to the maturity and fall of vegetation leaves, respectively;

[0021] S3.3. Calculate the temperature limit factor as follows:

[0022] ,

[0023] The humidity limiting factor is calculated as follows:

[0024] ,

[0025] Calculate the carbon dioxide limitation factor as follows:

[0026] ,

[0027] The construction environment restriction function is as follows:

[0028] ,

[0029] In the formula, f T is the temperature limiting factor, f w is the humidity limiting factor, f p is the phenological limiting factor, f c is the CO2 limiting factor, T opt is the optimal growth temperature for vegetation. T a is the air temperature, LSWI max It is the maximum vegetation moisture content index during the vegetation growing season. is the sensitivity coefficient of GPP to the change of atmospheric CO2, CO2 is the atmospheric carbon dioxide concentration, CO 2base is the baseline CO2 concentration.

[0030] S4. Based on the framework of double-leaf light energy utilization, the ecosystem carbon flux model was established according to the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation and the environmental limiting function to estimate the total primary productivity; the ecosystem carbon flux model is as follows:

[0031] ,

[0032] In the formula, , are the maximum light energy utilization rates of the sun leaves and shade leaves of vegetation, APAR su , APAR sh The photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation respectively. is the environmental limiting factor.

[0033] S5. Optimize the parameters of sun leaf light energy utilization efficiency and shade leaf light energy utilization efficiency in the ecosystem carbon flux model using the Markov chain Monte Carlo algorithm to obtain an optimized ecosystem carbon flux model;

[0034] S6. Using the collected multi-source satellite remote sensing data as input, the global gross primary productivity with high temporal and spatial resolution is calculated using the ecosystem carbon flux model with optimized parameters.

[0035] Furthermore, the aforementioned step S1 calculates the satellite T 0 Instantaneous shortwave radiationI s As follows:

[0036] ,

[0037] In the formula, S 0 represents the solar constant at the top of the atmosphere, is the empirical coefficient, is the solar zenith angle parameter, is the surface water pressure;

[0038] Calculate the daily solar radiation flux based on the sinusoidal model of the daily variation of surface temperature and the length of daytime I g As follows:

[0039] ,

[0040] ,

[0041] In the formula, I s ( T ) is any pixel in the image T The solar radiation value at the time, T r It's sunrise time. T s It is the sunset time. The sunrise time and sunset time are determined by the date and the latitude of the geographical location. T Any time of day.

[0042] Furthermore, in the aforementioned step S1, the diffuse photosynthetically active radiation ratio TD is calculated according to the solar altitude angle and the clear sky index. f , as follows:

[0043] ,

[0044] ,

[0045] ,

[0046] ,

[0047] In the formula, CI represents the influence of the atmosphere on the extraterrestrial radiation, I g is solar radiation, I0 is external radiation, I sc is the solar constant, t d Julian Day, θ is the solar altitude angle, D f is the scattering ratio, CI d is the daily clear sky index, is the total radiation persistence factor.

[0048] Furthermore, in the aforementioned step S1, the scattered photosynthetically active radiation PAR dif , direct photosynthetically active radiation PAR dir , diffuse photosynthetically active radiation under the canopy PAR dif，u The calculation is as follows:

[0049] ,

[0050] ,

[0051] ,

[0052] ,

[0053] In the formula, PAR is photosynthetically active radiation, 𝛼 is the proportional factor, is the clumping index, which depends on the vegetation type, and LAI is the leaf area index.

[0054] Furthermore, the aforementioned step S2 specifically includes: using a double S-shaped curve to perform nonlinear least square fitting to calculate vegetation phenological characteristics, which is used to evaluate and calculate the growth state of vegetation leaves and provide state variables for the establishment of phenological constraint factors, as shown in the following formula:

[0055] ,

[0056] Wherein, t represents time, a, b, c, d, and f and g represent fitting parameters, respectively. a and b determine the minimum and maximum values ​​of the curve, respectively. c and d control the time and slope of the rising phase of the curve, respectively, corresponding to the leaf growth process. g and f control the time and slope of the falling phase of the curve, respectively, corresponding to the leaf falling process.

[0057] Furthermore, the aforementioned step S5 specifically refers to two parameters of the ecosystem carbon flux model: maximum light energy utilization rate and maximum light energy utilization rate of shade leaves The optimization includes the following sub-steps:

[0058] S5.1. Using Bayesian theory, estimate parameters The posterior distribution of :

[0059] Remember the given The conditional density of the random variable X is , and 𝑥=(𝑥1,…𝑥𝑛) is a sample taken from the population X, and its likelihood function is:

[0060] ,

[0061] Sample 𝑥 and parameters The joint probability density of , obtained by Bayesian formula The posterior distribution of is as follows:

[0062] ,

[0063] Where 𝜋 is the probability distribution of the parameter, is the prior probability of the parameter, is the posterior probability of the parameter, is the likelihood function, i.e. the distribution of observed samples; integral solution interval Refers to the parameters All possible value ranges;

[0064] S5.2. Samples are extracted from the probability density using the Markov chain Monte Carlo algorithm, and the approximate distribution is obtained using kernel density estimation. By continuously updating sample information, the Markov chain model can fully search the model parameter space and eventually converge to the high probability density area. Specifically, the Metropolis–Hastings algorithm is used for sample sampling to complete parameter optimization: ,

[0065] In the formula, is the acceptance rate, is the optimal parameter, The proposed parameters generated by the proposal function conform to the normal distribution. Through the iterative process, when the acceptance rate If it is greater than the random number generated from a 0-1 uniform distribution during the iteration, the proposed parameter is the optimal parameter.

[0066] Furthermore, the aforementioned step S6 includes the following sub-steps:

[0067] S6.1 Collect multi-source satellite data and meteorological data, use the ecosystem carbon flux model with optimized parameters to invert the global total primary productivity, and verify its simulation accuracy based on the measured data of flux stations;

[0068] S6.2. Collect land surface process models based on physical mechanisms and verify the ability of the parameter-optimized ecosystem carbon flux model in simulating the spatiotemporal changes of total primary productivity based on their ensemble average values.

[0069] Compared with the prior art, the beneficial technical effects of the above technical solutions adopted in the present invention are as follows: a multi-source remote sensing ecosystem carbon flux estimation method based on phenology disclosed in the present invention uses multi-source remote sensing data to calculate the direct photosynthetically active radiation and scattered photosynthetically active radiation of the study area; obtains regional biological characteristic data of the same time series, such as induced chlorophyll fluorescence (SIF), etc., and uses a double S-shaped curve for nonlinear least squares fitting to extract vegetation phenological characteristics for evaluating the growth status of vegetation leaves; sets thresholds according to phenological characteristics and constructs phenological constraint factors; based on the double-leaf light energy utilization efficiency framework and combined with phenological constraint factors, establishes an ecosystem carbon flux model under the influence of phenology with high computational efficiency and few parameters; based on the measured data of flux sites, uses Markov chain Monte Carlo chains to optimize model parameters; based on multi-source remote sensing data, uses a parameter-optimized ecosystem carbon flux model to generate global high temporal and spatial resolution total primary productivity data and verify its ability to simulate temporal and spatial trend changes. The present invention proposes a phenology-based multi-source remote sensing ecosystem carbon flux estimation method that is superior to the traditional light energy utilization efficiency model and has significant advantages in simulating the spatiotemporal changes in total primary productivity, thereby improving the simulation accuracy of global carbon flux. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 It is a schematic diagram of the method provided by the present invention.

[0071] Figure 2 It is a vegetation phenology node map.

[0072] Figure 3 It is a parameter optimization diagram based on the Markov chain Monte Carlo chain; in the figure, (a) is the maximum light energy utilization optimization diagram of the sun leaves, and (b) is the maximum light energy utilization optimization diagram of the shade leaves.

[0073] Figure 4 This is a scatter plot comparison of the PTEC simulation values ​​of the ecosystem carbon flux model after parameter optimization and the flux observation values.

[0074] Figure 5 This is a trend comparison chart of the ecosystem carbon flux model PTEC after parameter optimization and the physics-based land-atmosphere-surface model. DETAILED DESCRIPTION

[0075] In order to better understand the technical content of the present invention, specific embodiments are described below in conjunction with the accompanying drawings.

[0076] Various aspects of the invention are described herein with reference to the accompanying drawings, in which many illustrative embodiments are shown. The embodiments of the invention are not limited to those described in the accompanying drawings. It should be understood that the invention is implemented by any of the various concepts and embodiments described above, as well as the concepts and embodiments described in detail below, because the concepts and embodiments disclosed in the invention are not limited to any implementation. In addition, some aspects disclosed in the invention may be used alone or in any appropriate combination with other aspects disclosed in the invention.

[0077] like Figure 1 As shown, the present invention provides a multi-source remote sensing ecosystem carbon flux estimation method based on phenology, comprising the following steps:

[0078] S1. Based on multi-source satellite remote sensing data and meteorological data, the satellite instantaneous shortwave radiation is calculated, and then the daily solar radiation flux is calculated; according to the solar altitude angle and the clear sky index, the proportion of scattered photosynthetically active radiation is calculated, and then the high temporal and spatial resolution photosynthetically active radiation data is inverted, including scattered photosynthetically active radiation, direct photosynthetically active radiation, and scattered photosynthetically active radiation under the canopy.

[0079] S2. Collect biological characteristic data of the same time series in the target area, calculate vegetation phenological characteristics, and provide state variables for the establishment of phenological constraint factors;

[0080] S3, according to step S1, obtaining photosynthetically active radiation data with high temporal and spatial resolution, calculating the photosynthetically active radiation absorbed by the sun leaves and shade leaves of the vegetation, defining the threshold according to the state variable of step S2, establishing a phenological constraint factor based on the surface water index, and then calculating and obtaining the temperature constraint factor, the humidity constraint factor, and the carbon dioxide constraint factor, and further constructing the environmental constraint function;

[0081] S4. Based on the two-leaf light energy utilization framework, the ecosystem carbon flux model is established according to the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation and the environmental limiting function to estimate the total primary productivity;

[0082] S5. Optimize the parameters of sun leaf light energy utilization efficiency and shade leaf light energy utilization efficiency in the ecosystem carbon flux model using the Markov chain Monte Carlo algorithm to obtain an optimized ecosystem carbon flux model;

[0083] S6. Using the collected multi-source satellite remote sensing data as input, the global gross primary productivity with high temporal and spatial resolution is calculated using the ecosystem carbon flux model with optimized parameters.

[0084] As a preferred embodiment of the present invention, step S1 uses multi-source satellite remote sensing data and reanalysis data to calculate satellite T 0 Instantaneous shortwave radiationI s As follows:

[0085] ,

[0086] In the formula, S 0 represents the solar constant at the top of the atmosphere, which can be obtained from MOD03 satellite data and is 1367W ; is the empirical coefficient, with a value of 0.0455; is the solar zenith angle parameter, the unit is rad; is the surface water pressure; the unit is hPa, which can be calculated using the Clausius-Clapeyron formula combined with the MOD07_L2 waypoint temperature; I s The unit is W·m -2 .

[0087] Calculate the daily solar radiation flux based on the sinusoidal model of the daily variation of surface temperature and the length of daytime I g As follows:

[0088] ,

[0089] ,

[0090] In the formula, I s ( T ) is any pixel in the image T The solar radiation value at the time, T r It's sunrise time. T s It is the sunset time. The sunrise time and sunset time are determined by the date and the latitude of the geographical location. T Any time of the day, I g The unit is W·m -2 .

[0091] As a preferred embodiment of the present invention, in step S1, the diffuse photosynthetically active radiation ratio TD is calculated according to the solar altitude angle and the clear sky index. f (%), as follows:

[0092] ,

[0093] ,

[0094] ,

[0095] ,

[0096] Where CI represents the influence of the atmosphere on the extraterrestrial radiation; I g is the solar radiation, the unit is W·m -2 ; I0 is the extraterrestrial radiation, the unit is W·m -2 ;I sc is the solar constant, with a value of 1370 W·m -2 ;t d Julian Day, θ is the solar altitude angle, D f is the scattering ratio, CI d is the daily clear sky index, is the total radiation persistence factor.

[0097] As a preferred embodiment of the present invention, in step S1, the scattered photosynthetically active radiation PAR dif , direct photosynthetically active radiation PAR dir , diffuse photosynthetically active radiation under the canopy PAR dif，u The calculation is as follows:

[0098] ,

[0099] ,

[0100] ,

[0101] ,

[0102] In the formula, PAR is photosynthetically active radiation, 𝛼 is the proportional factor, is the clumping index, which depends on the vegetation type, and LAI is the leaf area index, which can be obtained from the GLASS_LAI data product.

[0103] As a preferred embodiment of the present invention, refer to Figure 2 , step S2 specifically includes: using multi-source satellite data to collect regional biometric data of the same time series, such as induced chlorophyll fluorescence (SIF), and using double S-shaped curves for nonlinear least squares fitting to calculate vegetation phenological characteristics, which are used to evaluate and calculate the growth status of vegetation leaves and provide state variables for the establishment of phenological constraint factors, as shown in the following formula:

[0104] ,

[0105] Wherein, t represents time, a, b, c, d, f and g represent fitting parameters respectively, a and b determine the minimum and maximum values ​​of the curve respectively, c and d control the time and slope of the rising stage of the curve respectively, corresponding to the leaf growth process, g and f control the time and slope of the falling stage of the curve respectively, corresponding to the leaf fall process. SIF data can be obtained from GOSIF data products.

[0106] As a preferred embodiment of the present invention, the aforementioned step S3 includes the following sub-steps:

[0107] S3.1. Calculate the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation APAR su , APAR sh , as follows:

[0108] ,

[0109] ,

[0110] In the formula, LAI su and LAI sh are the leaf area index of sun leaves and shade leaves respectively, LAI is the leaf area index, is the clustering index, is the mean leaf sun angle, set to 60 degrees; is the solar altitude angle, C is the multiple scattering contribution of direct radiation in the vegetation canopy, PAR dif It is diffuse photosynthetically active radiation. PAR dir Direct photosynthetically active radiation PAR dif，u is the diffuse photosynthetically active radiation below the canopy;

[0111] Area index of sun leaf LAI su Leaf area index of leaves and shade leaves LAI sh , calculated as follows:

[0112] ,

[0113] ,

[0114] S3.2. Establish phenological constraint factors based on surface water index as follows:

[0115] ,

[0116] ,

[0117] Where LSWI is the vegetation moisture index, and Represent the near infrared and short wave infrared bands of satellite sensors, SIF se and SIF ma are the SIF values ​​corresponding to the maturity and fall of vegetation leaves, respectively;

[0118] S3.3. Calculate the temperature limit factor as follows:

[0119] ,

[0120] The humidity limiting factor is calculated as follows:

[0121] ,

[0122] Calculate the carbon dioxide limitation factor as follows:

[0123] ,

[0124] The construction environment restriction function is as follows:

[0125] ,

[0126] In the formula, f T is the temperature limiting factor, f w is the humidity limiting factor, f p is the phenological limiting factor, f c is the CO2 limiting factor, T opt is the optimal growth temperature for vegetation, which can be obtained directly from the global vegetation optimum temperature map; T a is the air temperature, obtained from the ERA5 data-driven product; LSWI max is the maximum vegetation moisture content index during the vegetation growing season. is the sensitivity coefficient of GPP to atmospheric CO2 changes, and in the embodiment, the value is 0.129; CO2 is the atmospheric carbon dioxide concentration, which can be obtained from the Global Monitoring Laboratory of the National Oceanic and Atmospheric Administration of the United States; CO 2base It is the baseline concentration of CO2, with a value of 370.69 ppm, based on the CO2 concentration in the atmosphere in 2001.

[0127] As a preferred embodiment of the present invention, the ecosystem carbon flux model is established as follows:

[0128] ,

[0129] In the formula, , are the maximum light energy utilization rates of the sun leaves and shade leaves of vegetation, APAR su , APAR sh The photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation respectively. is the environmental limiting factor.

[0130] As a preferred embodiment of the present invention, a multi-source remote sensing ecosystem carbon flux model based on phenology is constructed. The model has only two parameters, namely, the maximum light energy utilization rate of the sun leaf and maximum light energy utilization rate of shade leaves , for these two parameters, such as Figure 3 As shown in the figure, (a) is the maximum light energy utilization optimization diagram of the sun leaves, and (b) is the maximum light energy utilization optimization diagram of the shade leaves.

[0131] The Markov Chain Monte Carlo algorithm (MCMC) method is used for parameter optimization, which includes the following sub-steps:

[0132] S5.1. Using Bayesian theory, estimate parameters The posterior distribution of :

[0133] Remember the given The conditional density of the random variable X is , and 𝑥=(𝑥1,…𝑥𝑛) is a sample taken from the population X, and its likelihood function is:

[0134] ,

[0135] Sample 𝑥 and parameters The joint probability density of , obtained by Bayesian formula The posterior distribution of is as follows:

[0136] ,

[0137] Where 𝜋 is the probability distribution of the parameter, is the prior probability of the parameter, is the posterior probability of the parameter, is the likelihood function, i.e. the distribution of observed samples; integral solution interval Refers to the parameters All possible value ranges;

[0138] S5.2. Samples are extracted from the probability density using the Markov chain Monte Carlo algorithm, and the approximate distribution is obtained using kernel density estimation. By continuously updating sample information, the Markov chain model can fully search the model parameter space and eventually converge to the high probability density area. Specifically, the Metropolis–Hastings algorithm is used for sample sampling to complete parameter optimization: ,

[0139] In the formula, is the acceptance rate, is the optimal parameter, The proposed parameters generated by the proposal function conform to the normal distribution. Through the iterative process, when the acceptance rate If it is greater than the random number generated from a 0-1 uniform distribution during the iteration, the proposed parameter is the optimal parameter.

[0140] Furthermore, the aforementioned step S6 includes the following sub-steps:

[0141] S6.1 collects multi-source satellite data and meteorological data, uses an ecosystem carbon flux model with optimized parameters to invert the global gross primary productivity, and verifies its simulation accuracy based on measured data from flux stations. Data from 142 flux stations around the world can be obtained free of charge from the FLUNET website.

[0142] S6.2. Collect land surface process models based on physical mechanisms (TREDNY_V8 products) and verify the ability of the parameter-optimized ecosystem carbon flux model in simulating the spatiotemporal changes of total primary productivity based on their ensemble average values.

[0143] Figure 4 This is a scatter plot of the 8-day mean GPP value simulated by the ecosystem carbon flux model PTEC after parameter optimization and the 8-day mean GPP value observed at 142 flux sites around the world. The figure shows that there is a high consistency between the GPP simulation value and the observed value, and the determination coefficient r 2 The root mean square error is 1.99gCm -2 d -1 The global GPP annual average values ​​from 2001 to 2021 simulated by the ecosystem carbon flux model PTEC after parameter optimization show that high GPP value areas are mainly distributed in tropical regions, while low GPP value areas are mainly distributed in polar regions. Figure 5This is a comparison chart of the global GPP sum simulated by the ecosystem carbon flux model PTEC after parameter optimization and the GPP product ensemble mean of the physical-based land-atmosphere-surface model (TRENDY). TRENDY products are the main reference products for the global carbon budget report. They are composed of 18 global dynamic vegetation model assessments, and their ensemble mean has a good response relationship to CO2 trend changes. It can be seen from the figure that the interannual change trend of the PTEC model and the ensemble mean of the TRENDY product are highly consistent, with a correlation coefficient of up to 0.98. This shows that the PTEC model has the same ability to identify the spatiotemporal changes of GPP as the dynamic vegetation model.

[0144] Although the present invention has been described above with preferred embodiments, it is not intended to limit the present invention. A person skilled in the art of the present invention may make various modifications and improvements without departing from the spirit and scope of the present invention. Therefore, the scope of protection of the present invention shall be determined by the definition of the claims.

Claims

1. A multi-source remote sensing ecosystem carbon flux estimation method based on phenology, characterized in that: The following steps are involved: S1. Based on multi-source satellite remote sensing data and meteorological data, the satellite instantaneous shortwave radiation is calculated, and then the daily solar radiation flux is calculated; According to the solar altitude angle and the clear sky index, the proportion of diffuse photosynthetically active radiation is calculated, and then the high temporal and spatial resolution photosynthetically active radiation data is inverted, including diffuse photosynthetically active radiation, direct photosynthetically active radiation, and diffuse photosynthetically active radiation under the canopy. S2. Collect biological characteristic data of the same time series in the target area, calculate vegetation phenological characteristics, and provide state variables for the establishment of phenological constraint factors; S3, according to step S1, obtaining photosynthetically active radiation data with high temporal and spatial resolution, calculating the photosynthetically active radiation absorbed by the sun leaves and shade leaves of the vegetation, defining thresholds according to the state variables of step S2, establishing phenological constraint factors based on the surface water index, and then calculating and obtaining temperature constraint factors, humidity constraint factors, and carbon dioxide constraint factors, and further constructing an environmental constraint function; specifically comprising the following sub-steps: S3.

1. Calculate the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation APAR su , APAR sh , as follows: , , In the formula, LAI su and LAI sh are the leaf area index of sun leaves and shade leaves respectively, LAI is the leaf area index, is the clustering index, is the mean leaf sun angle, is the solar altitude angle, C is the multiple scattering contribution of direct radiation in the vegetation canopy, PAR dif It is diffuse photosynthetically active radiation. PAR dir Direct photosynthetically active radiation PAR dif，u is the diffuse photosynthetically active radiation below the canopy; Area index of sun leaf LAI su and leaf area index of shade leaves LAI sh , calculated as follows: , , S3.

2. Establish phenological constraint factors based on surface water index as follows: , , Where LSWI is the vegetation moisture index, and Represent the near infrared and short wave infrared bands of satellite sensors, SIF se and SIF ma are the SIF values ​​corresponding to the maturity and fall of vegetation leaves, respectively; S3.

3. Calculate the temperature limit factor as follows: , The humidity limiting factor is calculated as follows: , Calculate the carbon dioxide limitation factor as follows: , The construction environment restriction function is as follows: , In the formula, f T is the temperature limiting factor, f w is the humidity limiting factor, f p is the phenological limiting factor, f c is the CO2 limiting factor, T opt is the optimal growth temperature for vegetation. T a is the air temperature, LSWI max It is the maximum vegetation moisture content index during the vegetation growing season. is the sensitivity coefficient of GPP to the change of atmospheric CO2, CO2 is the atmospheric carbon dioxide concentration, CO 2base is the baseline CO2 concentration; S4. Based on the framework of double-leaf light energy utilization, the ecosystem carbon flux model was established according to the photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation and the environmental limiting function to estimate the total primary productivity; the ecosystem carbon flux model is as follows: , In the formula, , are the maximum light energy utilization rates of the sun leaves and shade leaves of vegetation, APAR su , APAR sh The photosynthetically active radiation absorbed by the sun leaves and shade leaves of vegetation respectively. It is the environmental limiting factor; S5. Optimize the parameters of sun leaf light energy utilization efficiency and shade leaf light energy utilization efficiency in the ecosystem carbon flux model using the Markov chain Monte Carlo algorithm to obtain an optimized ecosystem carbon flux model; S6. Using the collected multi-source satellite remote sensing data as input, the global gross primary productivity with high temporal and spatial resolution is calculated using the ecosystem carbon flux model with optimized parameters.

2. The multi-source remote sensing ecosystem carbon flux estimation method based on phenology according to claim 1 is characterized in that: Step S1: Calculate satellites T 0 Instantaneous shortwave radiation I s As follows: ， In the formula, S 0 represents the solar constant at the top of the atmosphere, is the empirical coefficient, is the solar zenith angle parameter, is the surface water pressure; Calculate the daily solar radiation flux based on the sinusoidal model of the daily variation of surface temperature and the length of daytime I g As follows: ， ， In the formula, I s ( T ) is any pixel in the image T The solar radiation value at the time, T r It's sunrise time. T s It is the sunset time. The sunrise time and sunset time are determined by the date and the latitude of the geographical location. T Any time of day.

3. The multi-source remote sensing ecosystem carbon flux estimation method based on phenology according to claim 1 is characterized in that: In step S1, the diffuse photosynthetically active radiation ratio TD is calculated according to the solar altitude angle and the clear sky index. f , as follows: , , ， ， In the formula, CI represents the influence of the atmosphere on the extraterrestrial radiation, I g is solar radiation, I0 is external radiation, I sc is the solar constant, t d Julian Day, θ is the solar altitude angle, D f is the scattering ratio, CI d is the daily clear sky index, is the total radiation persistence factor.

4. The method for estimating carbon flux of a multi-source remote sensing ecosystem based on phenology according to claim 3 is characterized in that: In step S1, scattered photosynthetically active radiation PAR dif , direct photosynthetically active radiation PAR dir , diffuse photosynthetically active radiation under the canopy PAR dif，u The calculation is as follows: , , , , In the formula, PAR is photosynthetically active radiation, 𝛼 is the proportional factor, is the clumping index, which depends on the vegetation type, and LAI is the leaf area index.

5. The method for estimating carbon flux of a multi-source remote sensing ecosystem based on phenology according to claim 4, characterized in that: Step S2 specifically includes: using a double S-shaped curve to perform nonlinear least squares fitting to calculate vegetation phenological characteristics, which is used to evaluate and calculate the growth state of vegetation leaves and provide state variables for the establishment of phenological constraint factors, as shown in the following formula: , Wherein, t represents time, a, b, c, d, and f and g represent fitting parameters, respectively. a and b determine the minimum and maximum values ​​of the curve, respectively. c and d control the time and slope of the rising phase of the curve, respectively, corresponding to the leaf growth process. g and f control the time and slope of the falling phase of the curve, respectively, corresponding to the leaf falling process.

6. The method for estimating carbon flux of a multi-source remote sensing ecosystem based on phenology according to claim 5, characterized in that: Step S5 specifically includes two parameters of the ecosystem carbon flux model: maximum light energy utilization rate and maximum light energy utilization rate of shade leaves The optimization includes the following sub-steps: S5.

1. Using Bayesian theory, estimate parameters The posterior distribution of : Remember the given The conditional density of the random variable X is , and 𝑥=(𝑥1,…𝑥𝑛) is a sample taken from the population X, and its likelihood function is: , Sample 𝑥 and parameters The joint probability density of , obtained by Bayesian formula The posterior distribution of is as follows: , Where 𝜋 is the probability distribution of the parameter, is the prior probability of the parameter, is the posterior probability of the parameter, is the likelihood function, i.e. the distribution of observed samples; integral solution interval Refers to the parameters All possible value ranges; S5.

2. Samples are extracted from the probability density using the Markov chain Monte Carlo algorithm, and the approximate distribution is obtained using kernel density estimation. By continuously updating sample information, the Markov chain model can fully search the model parameter space and eventually converge to the high probability density area. Specifically, the Metropolis–Hastings algorithm is used for sample sampling to complete parameter optimization: , In the formula, is the acceptance rate, is the optimal parameter, The proposed parameters generated by the proposal function conform to the normal distribution. Through the iterative process, when the acceptance rate If it is greater than the random number generated from a 0-1 uniform distribution during the iteration, the proposed parameter is the optimal parameter.

7. The method for estimating carbon flux of a multi-source remote sensing ecosystem based on phenology according to claim 6, characterized in that: Step S6 includes the following sub-steps: S6.1 Collect multi-source satellite data and meteorological data, use the ecosystem carbon flux model with optimized parameters to invert the global total primary productivity, and verify its simulation accuracy based on the measured data of flux stations; S6.

2. Collect land surface process models based on physical mechanisms and verify the ability of the parameter-optimized ecosystem carbon flux model in simulating the spatiotemporal changes of total primary productivity based on their ensemble average values.

Images