Photosynthetic phenological inversion system combining meteorological suitability and vegetation index

Abstract

This invention provides a photosynthetic phenological inversion system combining meteorological suitability and vegetation indices, belonging to the field of remote sensing and ecological environment monitoring technology. The system includes a data input module, a growing season index calculation module, an improved phenological index synthesis module, a time-series reconstruction and phenological extraction module, and a result output and verification module, connected sequentially. It acquires and preprocesses multi-source data; calculates the diurnal growing season index based on temperature, precipitation, photoperiod, and shortwave radiation factor; normalizes the index and couples it with the remote-sensed vegetation index to generate an improved photosynthetic phenological index; extracts the photosynthetic start and end dates using curve fitting and a dynamic threshold method; and finally verifies and determines the optimal technical solution using ground observation data to generate the final product. This invention effectively solves the problems of inaccurate physiological process characterization, high susceptibility to noise interference, and lack of reliable verification in traditional remote-sensing phenological extraction, achieving higher accuracy and more physiologically meaningful large-scale vegetation photosynthetic phenological information extraction.

Description

Technical Field

[0001] This invention belongs to the field of remote sensing and ecological environment monitoring technology, specifically relating to a photosynthetic inversion system that combines meteorological suitability and vegetation index. Background Technology

[0002] Vegetation photosynthetic phenology is a core indicator characterizing the seasonal changes in photosynthesis in ecosystems, and its accurate monitoring is of great significance for studying global carbon cycles, climate change, and ecosystem management. Currently, the extraction of large-scale phenological information mainly relies on remote sensing vegetation index time series, and the start and end dates of the growing season are determined by analyzing its seasonal fluctuations.

[0003] However, existing technologies have significant shortcomings. First, commonly used vegetation indices primarily reflect changes in canopy greenness or structure, which are temporally offset from the actual photosynthetic physiological activities of vegetation, making it difficult for the extracted "greenness phenology" to accurately represent "photosynthetic phenology." Second, remote sensing data is susceptible to interference from clouds, atmosphere, and soil background, resulting in significant time-series noise that affects the stability and accuracy of phenological parameter extraction. Furthermore, existing phenological products typically lack systematic validation with ground-measured photosynthetic data, making it difficult to meet the reliability requirements of high-precision ecological models.

[0004] Although existing studies have attempted to improve the situation through multi-index fusion or the introduction of meteorological factors, most have failed to couple the mechanisms of meteorological conditions with those of remote sensing observations, nor have they fundamentally solved the problem of the disconnect between vegetation indices and photosynthetic processes. Therefore, there is an urgent need for a new technical solution that can integrate meteorological driving mechanisms with remote sensing observations to provide more accurate and reliable large-scale vegetation photosynthetic phenological information. Summary of the Invention

[0005] The purpose of this invention is to provide a photosynthetic inversion system that combines meteorological suitability and vegetation index, which can effectively solve the problems in the background art.

[0006] To achieve the above objectives, the technical solution adopted by this invention is: a photosynthetic phenological inversion system combining meteorological suitability and vegetation index, comprising: The data input module is used to synchronously acquire flux tower data, reanalysis meteorological data and multi-source remote sensing vegetation index data, and perform unified spatiotemporal alignment preprocessing on the data to construct a multidimensional input dataset covering multiple vegetation types and multiple years. The growing season index calculation module is used to select four key meteorological factors that are highly correlated with the total primary productivity of vegetation, namely temperature, precipitation, photoperiod and shortwave radiation, based on a multidimensional input dataset. It calculates the growth scalar corresponding to each factor according to preset physiological threshold parameters corresponding to the vegetation type, and synthesizes the daily meteorological suitability growing season index through product operation. An improved phenological index synthesis module is used to perform annual minimum-maximum normalization on the growing season index to obtain a normalized growing season index, which is then coupled with a multi-source remote sensing vegetation index to generate an improved phenological index candidate sequence, and the candidate sequence is subjected to time-series smoothing. The time series reconstruction and phenology extraction module is used to perform nonlinear annual curve fitting on the smoothed improved phenological index candidate sequence, and extract the photosynthesis start date and photosynthesis end date based on the fitted curve using the dynamic threshold method. The results output and verification module is used to verify and evaluate the accuracy of candidate phenological products generated based on different combinations of remote sensing vegetation indices and different phenological extraction methods, using the phenological periods extracted from the total primary productivity data of the flux tower as a benchmark. It determines the optimal index-method combination and generates standardized photosynthetic phenological products based on this combination.

[0007] Furthermore, in the data input module: The flux tower data comes from the FLUXNET2015 dataset. The selected sites cover multiple vegetation types and exclude evergreen broad-leaved forest and farmland sites. The meteorological data was then analyzed from the ERA5-Land dataset. Preprocessing included resampling the original spatial resolution to the target resolution using bilinear interpolation, and aggregating hourly data into daily-scale variables. The multi-source remote sensing vegetation index data includes normalized vegetation index, enhanced vegetation index, nuclear normalized vegetation index, near-infrared reflectance vegetation index, and leaf area index. Preprocessing includes interpolating the time series to generate continuous diurnal series.

[0008] Furthermore, the growing season index calculation module includes: Temperature scalar quantum unit is used to calculate the temperature growth scalar based on the daily average air temperature using a three-segment nonlinear function that includes a lower limit threshold, an optimal temperature threshold, and an upper limit threshold. The precipitation scalar unit is used to calculate the precipitation growth scalar based on the cumulative average precipitation over the previous few days using a linear piecewise function that includes a minimum threshold and a maximum threshold. A photoperiodic scalar unit is used to calculate the photoperiodic growth scalar based on the length of day using a linear piecewise function. The shortwave radiation scalar unit is used to calculate the shortwave radiation growth scalar based on the daily average shortwave radiation intensity using a linear piecewise function.

[0009] Furthermore, the growing season index calculation module also includes a parameter optimization unit, which uses an optimization algorithm to optimize and determine the various physiological threshold parameters involved in the growth scalars of temperature, precipitation, photoperiod, and shortwave radiation, respectively, under the cross-validation framework, based on the phenological period extracted from the flux tower site as the baseline truth.

[0010] Furthermore, in the improved phenological index synthesis module, the coupling method is to perform daily, pixel-by-pixel multiplication of the normalized growing season index and the multi-source remote sensing vegetation index, and the temporal smoothing process adopts the sliding window moving average filtering algorithm.

[0011] Furthermore, in the time series reconstruction and phenology extraction module, the nonlinear annual curve fitting adopts a double logistic function model, which is a function containing six parameters: baseline offset, amplitude, spring half-saturation point, spring growth rate, autumn half-saturation point, and autumn senescence rate. The curve fitting process uses the multi-year average improved phenological index curve of the target pixel as the initial parameter value.

[0012] Furthermore, in the dynamic threshold method, the threshold is set as a fixed percentage of the annual variation of the improved phenological index; the photosynthesis start date is defined as the date on which the fitted curve value first exceeds the threshold during the spring curve rise; and the photosynthesis end date is defined as the date on which the fitted curve value last exceeds the threshold during the autumn curve decline.

[0013] Furthermore, in the results output and verification module, the process of extracting the phenological baseline from the total primary productivity of the flux tower includes: smoothing the daily total primary productivity sequence with a Gaussian weighted moving average, calculating its annual maximum value, and using a fixed proportion of this maximum value as a threshold to extract the photosynthetic start date and photosynthetic end date.

[0014] Furthermore, the results output and verification module uses at least one statistical indicator among correlation coefficient, root mean square error, and Kling-Gupta efficiency to quantitatively evaluate the accuracy of candidate phenological products. Based on the evaluation results, the optimal solution is determined to be a combination of coupling the normalized vegetation index as the base index and the normalized growing season index, and using a combination of dual logistic function fitting and dynamic thresholding to extract phenological parameters.

[0015] Furthermore, the system is configured to run in a high-performance computing cluster environment, capable of processing global-scale, daily temporal resolution remote sensing and meteorological data, and ultimately generating a global photosynthetic phenology raster product with a spatial resolution of 0.05 degrees and its metadata.

[0016] In summary, this application includes the following beneficial technical effects: This invention innovatively constructs a meteorological suitability growing season index, dynamically coupling four key meteorological factors—temperature, precipitation, photoperiod, and radiation—to form a core indicator directly representing the climate-driven force of photosynthesis. Subsequently, this index is used as a weight to fuse with a vegetation index, generating an improved phenological index. This allows remote sensing observation signals to be constrained and corrected to the dimension of photosynthetic physiological processes at the data source. This fundamental improvement ensures that the final extracted phenological dates more closely reflect the actual physiological activity sequence of vegetation, significantly reducing misjudgments of growing season length and ecosystem productivity estimation biases caused by physiological disconnects in traditional methods. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the architecture of a photosynthetic phenology inversion system that combines meteorological suitability and vegetation index; Figure 2 This is a schematic diagram illustrating the coupling mechanism between the growing season index (GSI) and the improved photosynthetic phenological index (MPI). Figure 3 This is a logical flowchart of the time sequence reconstruction and phenological extraction stages. Detailed Implementation

[0018] This embodiment is applied to a high-precision extraction task of vegetation photosynthetic phenology on a global scale. It aims to construct a phenological index that accurately reflects the photosynthetic physiological state of vegetation by integrating remote sensing observations and meteorological driving mechanisms, along with an optimized time-series reconstruction and parameter extraction process. The system is deployed in a high-performance computing cluster environment and is capable of processing global time-series data with a spatial resolution of 0.05 degrees and a daily scale from 2001 to 2022.

[0019] Example 1 First, the hardware and software collaborative platform of this invention is constructed from a system architecture perspective. This system comprises five core functional modules: a data input module, a growing season index (GSI) calculation module, an improved phenological index (MPI) synthesis module, a time series reconstruction and phenological extraction module, and a result output and verification module. Each module is interconnected with a distributed file system via a high-speed internal bus and runs in a parallel computing environment based on the Linux kernel.

[0020] First, there's the data input module, which serves as the system's multi-source heterogeneous data access subsystem. Its hardware configuration includes dual 10 Gigabit Ethernet interface cards, specifically Intel X710-DA2 models, for high-speed network data transmission. Local storage utilizes an NVMe solid-state storage array with a capacity of at least 20TB and a read / write bandwidth of at least 3GB / s to ensure efficient caching and reading / writing of massive amounts of input data. In the data preprocessing stage, the system can employ GPU acceleration units for parallel computing acceleration, such as using an NVIDIA A100 80GB GPU to perform computationally intensive preprocessing tasks like interpolation and resampling.

[0021] Those skilled in the art should understand that the above hardware configuration is a recommended high-performance implementation. In actual deployments, a general-purpose CPU combined with a corresponding scientific computing library can also be used to complete the same function, depending on the available computing resources.

[0022] This module is responsible for synchronously acquiring and integrating three types of key input data: The first category consists of flux tower ground validation data, sourced from the internationally recognized FLUXNET2015 dataset. The data was downloaded from the dataset's official FTP server via HTTPS, in NetCDF4 format. The dataset includes the following essential variables: total primary productivity, daily mean temperature, precipitation, shortwave radiation, and sunshine duration. These variables will serve as the ground validation benchmark for subsequent phenological extraction results.

[0023] The second category consists of meteorological reanalysis data, sourced from the ERA5-Land dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). The data was acquired in batches via an application programming interface (API) provided by the ECMWF Climate Data Store, with the original file format being GRIB2. The raw spatial resolution of this data is 0.1 degrees multiplied by 0.1 degrees, and the temporal resolution is hourly. The meteorological variables to be extracted include hourly near-surface air temperature, downward shortwave radiation intensity, and total precipitation.

[0024] The third category is remote sensing vegetation index data, specifically including Normalized Difference Vegetation Index (NDV), Enhanced Vegetation Index (EDI), Nuclear Normalized Difference Vegetation Index (NND), Near-Infrared Reflectance Vegetation Index (NIRS), and Leaf Area Index (LAI). The NND and EDI are derived from the MOD13C1 V061 product, while the LAI is derived from the GLASS dataset. This data was exported via the Google Earth Engine cloud platform's application programming interface (API) in GeoTIFF format. The spatial resolution is 0.05 degrees x 0.05 degrees, and the temporal resolution varies depending on the product, ranging from 8 to 16 days.

[0025] The system has a built-in subroutine for calculating vegetation indices: Calculation of Kernel Normalized Difference Vegetation Index (kNDVI): The kNDVI photosynthetic phenological inversion system, combining meteorological suitability and vegetation index, transforms the traditional linear relationship into a nonlinear one using the kernel function method to better capture the saturation effect under high vegetation cover. Its calculation is based on the reflectance of the red light band (e.g., Band 1, approximately 620-670 nm) of MODIS or similar sensors, combined with the photosynthetic phenological inversion system of meteorological suitability and vegetation index. A photosynthetic phenology retrieval system combining meteorological suitability and vegetation index, and a photosynthetic phenology retrieval system combining near-infrared band (e.g., Band 2, approximately 841-876 nm) reflectance with meteorological suitability and vegetation index. First, a photosynthetic inversion system combining kernel similarity, meteorological suitability, and vegetation index is calculated. : Among them, the photosynthetic phenology inversion system combines meteorological suitability and vegetation index. The photosynthetic phenology inversion system combining meteorological suitability and vegetation index is typically used as the scale parameter. The formula for calculating kNDVI is: Calculation of Near-Infrared Reflectance Vegetation Index (NIRv): The photosynthetic phenology inversion system combining meteorological suitability and vegetation index, NIRv is the product of near-infrared reflectance and NDVI, designed to enhance sensitivity to the proportion of photosynthetically active radiation absorbed. Its calculation is based on acquired NDVI data and the photosynthetic phenology inversion system combining near-infrared reflectance, meteorological suitability, and vegetation index. : The above calculations are performed after the spatiotemporal alignment step of data preprocessing and before time interpolation to ensure that all vegetation indices are generated and have a consistent spatiotemporal reference before entering subsequent modules.

[0026] Before entering the core processing flow, all input data must undergo uniform spatiotemporal alignment preprocessing to form a spatiotemporally consistent input dataset. This series of preprocessing operations can be implemented using scientific computing libraries on general-purpose computing platforms, such as GDAL, NumPy, and SciPy libraries in Python, or hardware acceleration using FPGA accelerator cards such as the Xilinx Alveo U280 to improve efficiency.

[0027] The preprocessing of meteorological data involves two steps. First, spatial resampling is performed: using bilinear interpolation, the original 0.1-degree resolution raster data is resampled to a target 0.05-degree grid. Bilinear interpolation is a common method based on a weighted average of the four nearest neighbor pixel values. Second, temporal aggregation is performed: the hourly data is aggregated to a daily scale. For temperature and shortwave radiation intensity, the arithmetic mean of all hourly data for the day is calculated to obtain the daily average. For precipitation, the sum of all hourly data for the day is calculated to obtain the total daily precipitation.

[0028] For remote sensing data preprocessing, the core is time-dimension interpolation. Due to irregular gaps in the original data (mainly caused by cloud cover), interpolation is needed for the time series of each pixel to generate a continuous daily series. A cubic spline interpolation method is employed. This method constructs a piecewise cubic polynomial to ensure that the curve has continuous first and second derivatives at the interpolation points, thereby generating a smooth daily-scale vegetation index series.

[0029] The preprocessing of flux tower site data mainly involves spatial matching. Based on the latitude and longitude coordinates of the flux tower sites, the nearest geographic raster cell is found in the preprocessed 0.05-degree spatial resolution remote sensing and meteorological raster data. The value of this cell is then assigned to the site, achieving spatial alignment between the site data and the raster data. This matching method is called the nearest neighbor method.

[0030] In summary, the data input module, through dedicated hardware and standardized preprocessing procedures, enables the synchronous acquisition, format unification, and spatiotemporal alignment of multi-source, heterogeneous, and multi-temporal resolution data, forming a high-quality, consistent input dataset. This dataset provides a reliable data foundation for subsequent growing season index calculations, improved phenological index synthesis, and phenological parameter extraction, ensuring high consistency and verifiability of the entire system from the data source, thereby supporting the generation of subsequent high-precision photosynthetic phenological products.

[0031] The GSI calculation module, as the core algorithm execution unit of the system, is responsible for constructing the diurnal growing season index. This module can be implemented on general-purpose computing nodes equipped with multi-core CPUs and parallel computing capabilities, such as servers using Intel Xeon series processors and more than 128GB of memory. Logically, the module is divided into four parallel sub-computing units, each processing one of the four meteorological factors: temperature, precipitation, photoperiod, and shortwave radiation.

[0032] Each sub-computation unit needs to load corresponding physiological threshold parameters based on the vegetation type of the target area during execution. The system supports seven vegetation types: deciduous broad-leaved forest, evergreen coniferous forest, grassland, mixed forest, shrubland, wetland, and savanna. For each vegetation type, there is an independent set of meteorological factor threshold parameters, which are stored in local files in common formats such as JSON for reading and calling during calculation.

[0033] Temperature growth scalar calculation: Based on the daily average temperature data T, a three-segment nonlinear response function is used to calculate the temperature growth scalar Ts. This function sets three key temperature thresholds: Minimum temperature threshold; Optimal temperature threshold; Maximum temperature threshold; The formula for calculating the temperature growth scalar Ts is as follows: When T ≤ or T ≥ hour: when < T < hour: The coefficients a, b, and c are obtained by solving the following system of linear equations: Calculation of precipitation growth scalar; Based on the cumulative average precipitation MP of the previous 10 days, the precipitation growth scalar Ps is calculated. First, the total precipitation of the previous 10 days is calculated, then divided by 10 to obtain the daily average precipitation MP. Based on preset minimum and maximum precipitation thresholds P_min and P_max, Ps is calculated using a linear piecewise function. When MP ≤ hour: When MP ≥ hour: when < MP < hour: Scalar calculation of photoperiodic growth; Based on the day length data PH, the photoperiodic growth scalar PHs is calculated. Day length PH is calculated using a standard astronomical formula, with the input parameters being geographic latitude and Julian day. This is based on a preset minimum day length threshold. and maximum threshold PHs is calculated using a linear piecewise function: When pH ≤ hour: When pH ≥ hour: when < PH < hour: Scalar calculation of shortwave radiation growth; The radiation growth scalar RADs are calculated based on the daily average shortwave radiation intensity RAD. This is done according to a preset minimum radiation intensity threshold. and maximum threshold RADs are calculated using a linear piecewise function: When RAD ≤ hour: When RAD ≥ hour: when < RAD < hour: Synthesis of the growing season index; After the above four growth scalars are calculated, the daily growth seasonal index (GSI) is synthesized through element-wise multiplication: This multiplication operation can be implemented using array operations in general-purpose programming languages, or it can be accelerated using the CPU's SIMD instruction set or GPU parallel computing.

[0034] Parameter optimization process; The threshold parameters in each of the above growth scalar functions are determined through optimization algorithms, including: Temperature parameters: , , ; Precipitation parameters: n、 ; Optical period parameters: , ; Shortwave radiation parameters: , ; The optimization process uses the phenological periods extracted from flux tower sites as the baseline truth and employs a genetic algorithm as the optimization engine, conducted within a 5-fold cross-validation framework. The specific optimization steps are as follows: 1. Set the search range for each parameter, for example: ∈ [-15, 15]℃ Other parameters are set within a reasonable range based on their physical meaning. 2. Initialize the genetic algorithm parameters: Population size: 50-200; Maximum number of iterations: 50-200; Crossover probability: 0.6-0.9; Mutation probability: 0.01-0.1; 3. Define the fitness function as Kling-Gupta efficiency (KGE), and the optimization objective is to maximize the KGE value.

[0035] 4. In each generation, evaluate the model performance corresponding to each parameter combination, and generate a new generation of population through selection, crossover, and mutation operations.

[0036] 5. Stop optimization when the maximum number of iterations is reached or the fitness value converges, and select the parameter combination with the highest KGE value as the optimal solution.

[0037] 6. Store the optimized parameters in the parameter library according to vegetation type for direct use in subsequent phenological extraction tasks.

[0038] In summary, through the structured computational process described above, the GSI calculation module transforms raw meteorological data into four independent scalars reflecting vegetation physiological limitations, ultimately synthesizing a comprehensive diurnal meteorological suitability index. The calculation of each scalar is based on explicit mathematical functions and data-driven optimization threshold parameters, ensuring the model's scientific validity and interpretability. The core innovation of the GSI calculation module lies in dynamically coupling discrete meteorological factors into continuous growth suitability indicators, providing precise meteorological constraints for subsequent fusion with remotely sensed vegetation indices. This effectively removes non-climate noise and significantly improves the accuracy and reliability of photosynthetic phenology extraction. The introduction of optimization further enhances the model's adaptability to different vegetation types and geographical regions, giving the entire system good scalability and practicality.

[0039] The core function of the MPI synthesis module is to receive the diurnal growth season index (GSI) sequence generated by the preceding module, and through three key steps—normalization, coupling, and smoothing—generate a set of improved phenological index candidate sequences that are of higher quality and more suitable for extracting photosynthetic physiological phenology.

[0040] This module first normalizes the input GSI sequence to eliminate the differences in the absolute magnitude of GSI between different regions and years, and unifies its numerical range to between 0 and 1, so that it can be used as a scalar weight.

[0041] For each individual cell, its daily-scale GSI time series for the entire year (365 or 366 days) is extracted. The original GSI value of each cell is converted into the Normalized Growing Season Index (NGSI) using a min-maximum normalization method.

[0042] The specific calculation formula is as follows: in: This indicates that the cell is in the th order. The normalized growth season index value for the day.

[0043] This indicates that the cell is in the th order. The original growing season index value for the day.

[0044] This represents the minimum value in the annual GSI sequence for that pixel.

[0045] This represents the maximum value of the cell in the annual GSI sequence.

[0046] After this step, in the NGSI sequence, the values ​​of dates representing completely unsuitable climate conditions for growth approach 0, and the values ​​of dates representing the most suitable climate conditions approach 1, thus forming a scalar field that can measure the potential stress of climate on photosynthesis on a daily basis.

[0047] After normalization, the module couples the meteorological suitability scalar NGSI with vegetation indices acquired through remote sensing, aiming to constrain and correct vegetation indices that purely reflect canopy greenness or structure using climatic conditions. The system uses five widely used vegetation indices as candidate bases, including the Normalized Difference Vegetation Index (NDVI), the Enhanced Vegetation Index (EVI), the Nuclear Normalized Difference Vegetation Index (kNDVI), the Near-Infrared Reflectance Vegetation Index (NIRv), and the Leaf Area Index (LAI).

[0048] The coupling operation is performed synchronously in the spatiotemporal dimensions, specifically in the form of daily, pixel-by-pixel multiplication. For each candidate vegetation index... ,in Representing the index type, calculate the daily value series of its corresponding Improved Phenological Index (MPI): in: Indicates based on the first Vegetation index, in the The improved phenological index value for the day.

[0049] Indicates the first Vegetation index in the 1st The original value of the day.

[0050] Through this multiplication operation, the original vegetation index signal is modulated by NGSI: on days with low meteorological suitability, even if the vegetation itself has a certain degree of greenness, its MPI value will be significantly suppressed; conversely, on days with high meteorological suitability, the vegetation index signal is preserved or even relatively enhanced. This process dynamically embeds meteorological information into phenological observations, generating five preliminary MPI candidate sequences.

[0051] Since the original remote sensing vegetation index and the initially synthesized MPI sequence usually contain high-frequency noise and outliers caused by clouds, aerosols, sensor noise, etc., the module needs to perform time-series smoothing on the MPI sequence in order to extract stable phenological change trends.

[0052] The smoothing process employs a sliding window moving average filtering algorithm. For each MPI candidate sequence, a fixed-length sliding window is defined. The window length is set to 21 days, with the center of the window being the date to be smoothed; that is, the window covers the 10 days before the current date, the current date, and the 10 days after the current date.

[0053] For the first in the sequence The smoothed MPI value is obtained by calculating the arithmetic mean of all valid MPI values ​​within the sliding window. The calculation formula is: in: Indicates the smoothed first... The index in the 1st The value of a day.

[0054] The summation symbol applies to the window from... arrive All days Accumulate.

[0055] This represents the number of valid data points within the current sliding window. For dates at the ends of the time series, when the window exceeds the series boundary, The value will be less than 21. In this case, only the average of the valid points within the window is taken, which is the conventional method for handling boundary cases.

[0056] This smoothing process effectively filters out random fluctuations shorter than the window scale (such as several days), while retaining low-frequency signals that reflect the seasonal process of vegetation growth and senescence, resulting in five smoothed, daily-scale MPI final candidate sequences.

[0057] Throughout the processing, the system needs to organize and manage intermediate and final data. Typically, five sets of MPI candidate sequences (before and after smoothing), the corresponding NGSI sequences, and the vegetation index identifier used for each pixel are stored together. This data can be organized into a multidimensional array or stored in a database by index, thus providing a complete and traceable data foundation for subsequent modules to perform unified accuracy verification and optimal scheme selection for multiple schemes.

[0058] The MPI synthesis module systematically incorporates meteorological suitability knowledge into remote sensing phenological analysis through a standard process of "normalization-coupling-smoothing". Normalization ensures a standardized expression of climate constraints, multiplicative coupling enables dynamic, pixel-level modulation of vegetation indices by meteorological factors, and moving average improves the signal-to-noise ratio of the data.

[0059] The module outputs five high-quality phenological index time series that have undergone meteorological correction and noise suppression. These series serve directly as input for the next stage, the "Time Series Reconstruction and Phenological Extraction Module." Their integration of meteorological and remote sensing information provides fundamental data support for the subsequent accurate extraction of the start and end dates of photosynthesis and the transformation from "greenness phenology" to "photosynthesis phenology."

[0060] The core function of the time series reconstruction and phenology extraction module is to model the input improved phenological index (MPI) time series and extract the photosynthetic start date and photosynthetic end date from it.

[0061] The module first performs a pre-smoothing of the original, uncorrected vegetation index sequence. The purpose of this step is to provide a high-quality reference sequence for subsequent comparative analyses, while also initially removing extreme noise caused by cloud pollution and other factors.

[0062] The Savitzky-Golay filter is used for processing. This is a filtering method based on local polynomial least squares fitting in the time domain. In specific implementation, the width of the filtering window is set to 15 days, meaning that each fitting uses data from the current date plus 7 days before and after it (a total of 15 points). The polynomial order used for fitting is set to 2. The algorithm is implemented through convolution operations. For each point in the sequence, a set of fixed convolution coefficients (obtainable through standard scientific computing library functions such as scipy.signal.savgol_coeffs) is used to perform a weighted average with the data within the window, directly outputting the smoothed value. For points at both ends of the time series, conventional methods such as mirror filling or window truncation are used to handle the boundaries.

[0063] The core of the module is to perform annual curve fitting on each MPI candidate sequence to characterize its complete seasonal growth trajectory. A bilogistic function is used as the fitting model, which is a six-parameter nonlinear function with the following standard mathematical expression: in: : Independent variable, representing the date sequence in a year, with a value range of 1 to 365 or 366.

[0064] The baseline offset parameter determines the baseline level of the curve throughout the year.

[0065] Amplitude parameter: The range of change of the control curve from the baseline level before (or after) the growing season to the peak of the growing season.

[0066] The half-saturation point parameter of the spring logistic function is related to the center position of the rising segment of the spring growth curve.

[0067] Spring growth rate parameter, which controls the steepness of the spring growth curve; the larger the value, the faster the rise.

[0068] The half-saturation point parameter of the autumn logistic function is related to the center position of the descending segment of the autumn aging curve.

[0069] The autumn aging rate parameter controls the steepness of the autumn aging curve decline; the higher the value, the faster the decline.

[0070] The goal of fitting is to find a set of optimal parameters. This makes the model curve Compared with the actual MPI observation sequence The goal is to minimize the sum of squared errors between the possible values. This is a nonlinear least squares optimization problem.

[0071] The Levenberg-Marquardt algorithm is used for solving the problem. This algorithm combines the Gauss-Newton method with gradient descent, and can effectively handle nonlinear fitting. The algorithm requires a set of initial parameter values. The method for setting initial values ​​is as follows: 1. Calculate the multi-year average MPI curve: For a target pixel, average its MPI sequence over many years (e.g., 5-10 years) daily to obtain a smooth, representative annual curve.

[0072] 2. Automatic Inflection Point Identification: On the multi-year average curve, identify the two points with the largest changes in curvature, and use them as the inflection points for the spring rise and autumn fall, respectively, and estimate their corresponding dates. and .

[0073] 3. Calculate the initial parameters: Let , .make , (An empirical starting value) , These values ​​will be used as Input optimization algorithm.

[0074] The optimization process is iterative until the parameter change or error decreases to less than the preset tolerance (e.g., ...). Stop when the optimal parameters are output. Thus, the fitted curve is obtained.

[0075] After obtaining the annual MPI fitting curve, the specific phenological dates are extracted using the dynamic thresholding method. First, the maximum value for that year is found from the fitting curve. and minimum value Calculate the annual variation range. Set dynamic threshold for: That is, the threshold is located above the minimum value, at a position of 20% of the range.

[0076] The steps for extracting the photosynthetic initiation date are as follows: 1. Determine the search range: usually from day 1 to the middle of the year (e.g., day 180).

[0077] 2. Forward scan: Starting from the beginning of the search interval, check the values ​​of the fitted curve day by day. .

[0078] 3. Judgment condition: Record the first condition that is met. Date .

[0079] 4. Output: This date This date was determined to be the start of photosynthesis.

[0080] The extraction steps for the end of photosynthesis are as follows: 1. Determine the search range: usually from the middle of the year (e.g., day 181) to the last day (day 365 or 366).

[0081] 2. Reverse Scan: Starting from the end date of the search interval, check the values ​​of the fitted curve day by day. .

[0082] 3. Judgment condition: Record the last one that meets the condition. Date .

[0083] 4. Output: This date It was determined to be the end date of photosynthesis.

[0084] In addition to the core method combination mentioned above, the module also runs several other phenological extraction methods in parallel as benchmarks or alternatives for performance comparison. These methods include, but are not limited to: 1. Derivative method: Calculate the first derivative of the fitted curve or the smoothed original sequence, and take the date corresponding to the maximum value in spring as SOS and the date corresponding to the minimum value in autumn as EOS.

[0085] 2. Polynomial fitting method: Use a high-order polynomial (such as 4th or 5th order) to directly fit the MPI sequence, and then find the extreme points by taking the derivative to determine the phenological date.

[0086] The system maintains the computation flow of these multiple methods (a total of 8). For each MPI candidate sequence, the SOS and EOS results extracted by all methods, along with the method identifier and the MPI index type used, are stored in a structured manner (e.g., stored in a Redis in-memory database or written to a temporary file), forming a complete set of candidate results for subsequent verification modules to perform system evaluation.

[0087] The time series reconstruction and phenological extraction module transforms continuous, noisy MPI index sequences into precise, physiologically meaningful phenological dates through a coherent process of "pre-smoothing reference sequence - nonlinear curve fitting - adaptive threshold extraction." The dual logistic function model effectively describes the seasonal initiation, peak, and decline of photosynthetic activity, while the dynamic threshold rule provides a stable and objective extraction criterion.

[0088] This module produces a series of phenological date "candidate products" based on different vegetation indices and extraction methods. Their diversity and traceability lay a solid foundation for the next stage of results output and verification, which involves rigorous accuracy evaluation and optimal solution selection. This module is the core link in realizing the crucial transformation from "index" to "date".

[0089] Parallel implementation of multiple phenological extraction methods and generation of candidate results. To determine the optimal phenological extraction strategy for the Improved Phenological Index (MPI), this module runs multiple combinations of time-series reconstruction and phenological extraction methods in parallel, generating a series of candidate phenological products for subsequent validation and screening. These combinations include, but are not limited to, the following eight: 1. Double Logistic Function and First Derivative Method (DLF-FD)After fitting the MPI sequence with a double logistic function, the first derivative of the fitted curve is calculated. The date corresponding to the maximum value of the spring derivative is taken as the start date of photosynthesis, and the date corresponding to the minimum value of the autumn derivative is taken as the end date of photosynthesis.

[0090] 2. Double Logistic Function and Second Derivative Method (DLF-SD) Calculate the second derivative of the fitted curve, and take the date corresponding to the maximum value of the second derivative in spring as the start date of photosynthesis and the date corresponding to the maximum value of the second derivative in autumn as the end date of photosynthesis.

[0091] 3. Fixed threshold and polynomial fitting method (FT-PF) A fixed threshold is determined based on the multi-year average MPI change rate. Higher-order polynomials are used to fit the MPI sequences for the first and second halves of the year, and the date when the fitted value equals the threshold is taken as the phenological period.

[0092] 4. Dual Logistic Function and Multi-Scale Dynamic Thresholding (DLF-DT) This is the core methodology combination for system evaluation. After fitting the MPI sequence with a bilogistic function, the fitted curve values ​​are normalized to obtain the MPI_ratio. Subsequently, a set of preset methods are applied in parallel. Dynamic threshold ratio Extraction is performed at percentages of 10%, 20%, 30%, 40%, and 50%. For each threshold percentage r, the dynamic threshold Th_r is calculated using the following formula: Wherein, MPI_max and MPI_min are the annual maximum and minimum values ​​of the fitted curve, respectively. Photosynthesis start date Defined as the date on which the MPI_ratio value first exceeds Th_r during the spring curve's upward trend; Photosynthesis End Date Defined as the date on which the MPI_ratio value last exceeds Th_r during the autumn phenological curve's decline. This process, using only one method combination, can generate five sets (corresponding to five threshold ratios) of candidate phenological dates.

[0093] The system manages the above and other alternative method combinations in a unified manner. For each MPI candidate sequence (such as NDVI-based MPI), the phenological date results extracted by all methods, along with the method identifiers and the threshold parameters used, are all stored in a structured manner, forming a complete... Candidate Phenological Product Set This product set serves as direct input for the next stage of accuracy verification and optimal solution selection.

[0094] The results output and verification module is responsible for quality control and scheme optimization before the final product is generated. This module first needs to establish a reliable ground phenology benchmark to objectively evaluate the accuracy of different remote sensing extraction schemes, and then automatically select the optimal scheme and generate standardized global phenology products.

[0095] To obtain accurate phenological information from the ground, this module processes total primary productivity (GPP) time series data from networks such as FLUXNET. To improve data comparability and highlight seasonal trends, the original daily GPP series needs to be smoothed.

[0096] A Gaussian weighted moving average algorithm is used for smoothing. Specifically, a sliding window of length 30 days is defined. For the center day t of the window, the weight of day k within the window is... Calculated using the Gaussian function: The standard deviation σ was set to 5. Then, a weighted average was calculated for all GPP values ​​within the window to obtain the smoothed GPP values.

[0097] Find the annual maximum value from the smoothed multi-year average GPP series, denoted as . .

[0098] Set a fixed ratio coefficient, such as 20%, and calculate the phenological extraction threshold accordingly: .

[0099] Photosynthesis Initiation Date Reference The extraction process is as follows: Starting from the beginning of the year (day 1), the smoothed GPP sequence is checked daily, and the first sequence that satisfies GPP(t) ≥ P's date t is recorded as .

[0100] Photosynthesis End Date Baseline The extraction process is as follows: starting from the middle of the year (e.g., day 180), the smoothed GPP sequence is checked day by day in reverse order, and the last one that satisfies GPP(t) ≥ The date t is recorded as .

[0101] The module receives 40 candidate phenological products from the preceding module, which are generated by combining 5 MPI indices with 8 extraction methods.

[0102] Using all available flux tower sites as validation points, the accuracy of each candidate product was verified. For each site, the remote sensing extraction results (SOS, EOS) corresponding to that site location were compared with the aforementioned ground baseline (…). , A year-by-year comparison was conducted.

[0103] Statistical indicators such as correlation coefficient (R), root mean square error (RMSE), and Kling-Gupta efficiency (KGE) were used to quantitatively evaluate the accuracy of each candidate product.

[0104] The specific calculation method is as follows: 1. Correlation coefficient R: The Pearson correlation coefficient between the remote sensing extracted sequence and the ground reference sequence is calculated to assess the consistency of their interannual variations. in, and These represent the remote sensing extraction date and ground reference date for year i, respectively. and It is their mean. This represents the total number of years included in the comparison.

[0105] 2. Root Mean Square Error (RMSE): Measures the average absolute deviation between the remote sensing extraction date and the ground reference date.

[0106] 3. Kling-Gupta efficiency (KGE): A comprehensive evaluation index that considers correlation, bias ratio, and variability ratio.

[0107] in, It is the correlation coefficient. It is the variability ratio. It is the deviation ratio. and These are the standard deviations of the remote sensing extracted sequence and the ground reference sequence, respectively. and It is their mean.

[0108] The system sorts and filters all 40 candidate solutions based on the comprehensive accuracy metrics of all verification sites. The operation process is as follows: First, calculate the median or average of the SOS and EOS accuracy metrics (R, RMSE, KGE) for each scheme across all validation sites.

[0109] Next, set the filtering rules. For example, prioritize the option with the highest average KGE value. Simultaneously, constraints can be set, such as requiring the average RMSE to be less than an acceptable threshold, such as 10 days.

[0110] The system uses automated scripts or decision logic to rank the performance of all schemes according to preset rules, and finally selects one or more schemes that perform best and most stably on most indicators. For example, through large-scale verification, a typical optimal combination has been determined as follows: an improved phenological index based on NDVI, curve fitting using a double logistic function, and extraction of phenological parameters using a 20% dynamic threshold.

[0111] After determining the optimal solution, the system processes the input data of all pixels globally based on that solution to generate the final global-scale photosynthetic phenological products.

[0112] The product is output in a standard geospatial raster data format, such as GeoTIFF. Each file corresponds to the global distribution of a phenological parameter (SOS or EOS) for a given year, with the cell value being the Julian Day.

[0113] Simultaneously, a detailed metadata file is generated and output. Metadata typically includes: product name, spatiotemporal range, spatial resolution, data format, production date, a detailed description of the optimal solution used, explanation of missing data values, a spatial uncertainty layer based on validation results, referenced vegetation type classification information, and data referencing methods, etc.

[0114] The results output and verification module serves as the final quality controller and product finalizer for the entire technical process. It establishes a rigorous ground-based benchmark verification system, quantitatively evaluates multiple candidate solutions, and automatically selects the optimal technical path using a data-driven approach.

[0115] This module not only ensures the scientific validity and reliability of the final phenological products, but also reflects the objectivity and repeatability of the methodology itself. Its standardized global photosynthetic phenological products directly serve fields such as global change ecology and carbon cycle simulation, effectively addressing the core issue of the disconnect between traditional remote sensing phenological products and real photosynthetic physiological processes.

[0116] With the support of the above system architecture, the workflow of this embodiment unfolds as follows according to time sequence logic: After the system starts, the data input module first concurrently downloads global ERA5-Land meteorological data, MOD13C1 / GLASS remote sensing data, and FLUXNET2015 flux tower data from 2001 to 2022. All data undergo unified spatiotemporal alignment preprocessing (including bilinear interpolation resampling of meteorological data to 0.05°, temporal interpolation of remote sensing data to generate continuous daily series, and spatial nearest neighbor matching between flux tower stations and raster data) to form a multidimensional data cube that is updated daily with a 0.05° grid unit.

[0117] Subsequently, the GSI calculation module traverses each grid cell, loading corresponding optimized threshold parameters (physiological thresholds for temperature, precipitation, photoperiod, and shortwave radiation) based on its vegetation type label (determined by the IGBP classification map). It then calculates the temperature growth scalar (a three-segment nonlinear function), precipitation growth scalar (a linear piecewise function based on the cumulative average precipitation of the previous 10 days), photoperiod growth scalar (a linear piecewise function based on the length of day), and shortwave radiation growth scalar (a linear piecewise function based on the average daily radiation intensity) in parallel on a daily basis. Finally, it synthesizes the daily-scale growing season index (GSI) through element-wise multiplication. This process is executed in parallel across approximately 250 million land grids globally, taking approximately 48 hours (on a 128-node cluster).

[0118] The generated GSI sequences are then fed into the MPI synthesis module. This module first performs year-round minimum-maximum normalization on the GSI to obtain the Normalized Growing Season Index (NGSI), ensuring it approaches zero during the non-growing season, thus effectively suppressing soil background and snow cover interference. Next, the system performs daily, pixel-by-pixel multiplication coupling of the NGSI with five vegetation indices: NDVI, EVI, kNDVI, NIRv, and LAI, generating five sets of improved phenological indices (MPIs). Each MPI set is filtered using a 21-day sliding window moving average to eliminate signal fluctuations caused by short-term weather events (such as single-day heavy rainfall) and preserve the true phenological trend.

[0119] The smoothed MPI sequence is then passed to the time-series reconstruction and phenological extraction module. This module first pre-smooths the original vegetation index time series using Savitzky-Golay filtering, serving as a reference sequence. Subsequently, each MPI candidate sequence is fitted with an annual curve using a dual logistic function (six-parameter model). During the fitting process, the algorithm uses the multi-year average curve as initial values ​​and adjusts parameters through Levenberg-Marquardt nonlinear least squares optimization to ensure the fitted curve closely matches the physiological trajectory of MPI changes. After fitting, the system calculates the annual MPI variation amplitude from the fitted curve and sets a 20% dynamic threshold (threshold = minimum value + 0.2 × variation amplitude). In spring, the algorithm scans the fitted curve forward day by day from day 1, recording the date the first time it exceeds the threshold as the photosynthetic start date (SOS). In autumn, it scans backward day by day from the last day, recording the last date the threshold is exceeded as the photosynthetic end date (EOS). This process runs simultaneously with a combination of multiple phenological extraction methods, generating a total of 40 candidate phenological products.

[0120] Finally, the results output and validation module calls the GPP phenological benchmarks from 84 flux tower sites (extracted by applying a 20% threshold method to the smoothed GPP sequences) to cross-validate 40 candidate products. Statistical results show that the combination of NDVI-MPI with dual logistic function fitting and 20% dynamic threshold extraction significantly outperforms the original NDVI (R=0.75, RMSE=18.3 days, KGE=0.54) in SOS extraction (R=0.82, RMSE=8.2 days, KGE=0.75); the improvement is even more significant in EOS (R increased from 0.61 to 0.75, RMSE decreased from 21.5 days to 7.1 days). Based on this, the system automatically locks the optimal solution and processes all global image data according to this solution to generate the final global-scale photosynthetic phenological products (SOS / EOS), which are written to a distributed storage system in GeoTIFF format for use by downstream ecosystem models.

[0121] This embodiment fully realizes the automated processing of the entire chain from multi-source data fusion, meteorological suitability index construction, improved phenological index synthesis to high-precision phenological parameter extraction, which fully demonstrates the high degree of synergy and technological advancement of the present invention in system architecture and method process.

[0122] Example 2 This second embodiment, based on the first embodiment, provides alternative implementation schemes for certain technical aspects, aiming to further illustrate the flexibility and scalability of the invention and enhance the robustness of the technical solution in different application scenarios. The application scenario also focuses on the extraction of global vegetation photosynthetic phenology, but a different strategy is adopted for key algorithm components compared to the first embodiment.

[0123] At the overall system architecture level, this embodiment retains the basic structure and functions of the data input module and the result output and verification module. However, alternative algorithm units are integrated into the GSI calculation module, MPI synthesis module, and time series reconstruction and phenological extraction module, as detailed below: In this embodiment, the optimization algorithm used in the GSI calculation module to determine the physiological threshold parameters of each growth scalar (including the thresholds of temperature, precipitation, photoperiod, and shortwave radiation) no longer uses the genetic algorithm (GA) selected in Embodiment 1, but instead integrates the Particle Swarm Optimization (PSO) engine as the optimization core.

[0124] The PSO engine can be deployed on dedicated AI acceleration hardware (such as the Graphcore IPU-M2000 accelerator card), utilizing its numerous internal parallel processing cores (e.g., 1472 independent IPU-Cores) to perform parallel evaluation of thousands of parameter particles, thereby accelerating the optimization process. The PSO engine's inputs include: the true values ​​of phenological periods extracted from flux tower data (as the optimization target), and preset reasonable search ranges for the thresholds of each meteorological factor. Its output is a set of optimal threshold parameters optimized for each of the seven vegetation types (deciduous broad-leaved forest, evergreen coniferous forest, grassland, mixed forest, shrubland, wetland, and savanna).

[0125] The optimization process employs a standard particle swarm optimization algorithm with iterative update rules. For a given set of particles... A population of particles, of which the first... The particles are in the iteration number of The state of time is determined by its position in the parameter space. and speed Indicates. Location. For a specific set of threshold parameters, the speed It determines the direction and distance of movement in the next iteration.

[0126] Particle state updates follow the formula below: in: and : respectively represent the first The particle in the first Second and third The velocity vector at the next iteration.

[0127] and : respectively represent the first The particle in the first Second and third The position vector at the next iteration.

[0128] Inertia weight, used to balance global exploration and local development capabilities; in this embodiment, it is taken as... .

[0129] , The learning factor adjusts the step size of the particle's flight towards the individual's historical best and the group's historical best, respectively. In this embodiment, it is taken as... .

[0130] , : A random number uniformly distributed within the interval [0,1], used to introduce randomness.

[0131] : No. The position with the highest fitness found by an individual particle in its historical iterations (the individual optimal solution).

[0132] The position with the highest fitness found by the entire particle swarm in the historical iterations (the global optimal solution).

[0133] The objective function is the same as in Example 1, which is to maximize Kling-Gupta efficiency (KGE). Iterations are performed within a cross-validation framework until the preset maximum number of iterations or fitness convergence is reached. Experiments show that the PSO algorithm, with the same computational resources, can achieve parameter optimization results comparable to the genetic algorithm (GA), and typically has a faster convergence speed (an improvement of approximately 15%).

[0134] In the MPI synthesis module of this embodiment, the coupling method between the Normalized Growing Season Index (NGSI) and the Remote Sensing Vegetation Index (VI) has been expanded. It is no longer limited to the simple element-wise product described in Embodiment 1. Instead, it introduces an exponentially weighted coupling form, the mathematical expression of which is: in: Based on the first Improved phenological indices are calculated using vegetation indices (such as NDVI, EVI, etc.).

[0135] : No. The raw value of the vegetation index.

[0136] Normalized growth season index at the same time and space location.

[0137] : The weighting index of the vegetation index item, which is an adjustable parameter.

[0138] The weighting index of the National Weather Suitability Index (NGSI) is an adjustable parameter.

[0139] and These are learnable weight parameters whose values ​​determine the relative contribution strength of remote sensing observation signals and meteorological driving signals in the final MPI. The system incorporates a weight optimization submodule, which uses a subset of flux tower station GPP phenological baseline data as a validation set and employs a grid search method within a predefined parameter space (e.g., , Within ) , traverse different Parameter combinations. For each set of parameters, the corresponding MPI sequence is calculated and phenological periods are extracted, then compared with a ground baseline. Finally, the set that minimizes the root mean square error (RMSE) of the phenological extraction results is selected. As the optimal weight for this region or vegetation type.

[0140] For example, the optimization results in the temperate deciduous forest region show that the optimal parameters are: , This indicates that in this region, the weighting of the meteorological suitability index (NGSI) contribution to the final phenological index signal is slightly higher than that of the original vegetation greenness index (VI) itself. This flexible coupling form allows the system to adaptively adjust the balance between meteorological and remote sensing information according to different ecological environments (such as arid and semi-arid savanna regions with severe water stress), and more accurately reflect the actual control effect of environmental stress on the initiation of photosynthetic phenology.

[0141] In this embodiment, different curve fitting models and threshold extraction methods were explored in the time series reconstruction and phenology extraction module: 1. Asymmetric Gaussian Function (AGF) Fitting Model: As an alternative to the bilogistic function, this embodiment uses AGF to fit the annual time series of the MPI. The AGF model can better characterize the asymmetric phenological curves during the rising and falling phases of the growing season. Its expression is: in: : The date sequence within a year.

[0142] : The amplitude of the fitted curve, i.e., the intensity of the peak during the growing season relative to the baseline.

[0143] The fitted curve reached its peak value. The date (peak day).

[0144] Parameters that control the width of the rising segment of the spring growth curve. The larger the value, the more gradual the rise.

[0145] : Parameters that control the width of the downward segment of the autumn aging curve The larger the value, the more gradual the decline.

[0146] The fitting process uses nonlinear least squares methods (such as the Levenberg-Marquardt algorithm) to solve for the parameters. To obtain better initial values ​​and avoid local optima, Bayesian optimization and other strategies can be used for parameter initialization. The phenological date extraction still employs a dynamic thresholding method (e.g., a 20% threshold), but the threshold is set based on the maximum, minimum, and magnitude of change of the AGF-fitted curve. Cross-validation shows that in regions such as the Mediterranean climate zone, where summer drought leads to significant asymmetry in the growing season curve, the accuracy of the photosynthetic end date (EOS) extracted using AGF fitting is approximately 3.2% higher on average than that using the dual logistic function.

[0147] 2. Fixed Threshold Extraction Method: As a simplified alternative to the dynamic thresholding method, this embodiment also presets multiple sets of fixed MPI absolute value thresholds (e.g., 0.2, 0.3, 0.4), allowing the system to adaptively select the applicable fixed threshold based on the vegetation type of the pixel. For example, a threshold of 0.25 can be set for grassland type, and a threshold of 0.35 can be set for forest type. The extraction rule is similar to the dynamic thresholding method: on the fitted curve or smoothed MPI sequence, find the date that first exceeds (for SOS) or last exceeds (for EOS) the fixed threshold. Although this method may sacrifice accuracy slightly (e.g., RMSE may increase by about 2 days), it greatly reduces computational complexity and is suitable for edge computing or real-time monitoring scenarios with limited computing resources (e.g., lightweight phenological monitoring systems mounted on UAVs).

[0148] In summary, Example 2, by introducing alternative technical components such as the Particle Swarm Optimization (PSO) algorithm, the exponentially weighted coupling model, the asymmetric Gaussian function (AGF) fitting, and the fixed threshold method, fully verifies the high flexibility, configurability, and robustness of the technical framework described in this invention. These alternative solutions provide feasible technical options for dealing with different data types, regional characteristics, accuracy requirements, and computational constraints, further expanding the scope of application of this invention and consolidating its technical advantages.

[0149] In summary, Example 2, by replacing the core algorithm components, verifies the flexibility and scalability of the technical framework of the present invention, provides diverse implementation schemes for different computing environments and accuracy requirements, and further consolidates the technical advantages and application breadth of the present invention.

[0150] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the present invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the present invention. Therefore, the embodiments should be regarded as exemplary and non-limiting in all respects.

[0151] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A photosynthetic phenology inversion system combining meteorological suitability and vegetation index, characterized in that, include: The data input module is used to synchronously acquire flux tower data, reanalysis meteorological data and multi-source remote sensing vegetation index data, and perform unified spatiotemporal alignment preprocessing on the data to construct a multidimensional input dataset covering multiple vegetation types and multiple years. The growing season index calculation module is used to select four key meteorological factors that are highly correlated with the total primary productivity of vegetation, namely temperature, precipitation, photoperiod and shortwave radiation, based on a multidimensional input dataset. It calculates the growth scalar corresponding to each factor according to preset physiological threshold parameters corresponding to the vegetation type, and synthesizes the daily meteorological suitability growing season index through product operation. An improved phenological index synthesis module is used to perform annual minimum-maximum normalization on the growing season index to obtain a normalized growing season index, which is then coupled with a multi-source remote sensing vegetation index to generate an improved phenological index candidate sequence, and the candidate sequence is subjected to time-series smoothing. The time series reconstruction and phenology extraction module is used to perform nonlinear annual curve fitting on the smoothed improved phenological index candidate sequence, and extract the photosynthesis start date and photosynthesis end date based on the fitted curve using the dynamic threshold method. The results output and verification module is used to verify and evaluate the accuracy of candidate phenological products generated based on different combinations of remote sensing vegetation indices and different phenological extraction methods, using the phenological periods extracted from the total primary productivity data of the flux tower as a benchmark. It determines the optimal index-method combination and generates standardized photosynthetic phenological products based on this combination.

2. The photosynthetic phenology inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, In the data input module: The flux tower data comes from the FLUXNET2015 dataset. The selected sites cover multiple vegetation types and exclude evergreen broad-leaved forest and farmland sites. The meteorological data was then analyzed from the ERA5-Land dataset. Preprocessing included resampling the original spatial resolution to the target resolution using bilinear interpolation, and aggregating hourly data into daily-scale variables. Multi-source remote sensing vegetation index data include directly acquired normalized vegetation index, enhanced vegetation index, leaf area index, as well as calculated kernel normalized vegetation index and near-infrared reflectance vegetation index. Among them, the nuclear normalized vegetation index is calculated based on the reflectance of red light and near-infrared bands using the kernel function method, and the near-infrared reflectance vegetation index is obtained by multiplying the near-infrared band reflectance and the normalized vegetation index. Preprocessing includes interpolating the time series to generate continuous diurnal scale series, and performing corresponding band operations and synthesis on the indices to be calculated.

3. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, The growing season index calculation module includes: Temperature scalar quantum unit is used to calculate the temperature growth scalar based on the daily average air temperature using a three-segment nonlinear function that includes a lower limit threshold, an optimal temperature threshold, and an upper limit threshold. The precipitation scalar unit is used to calculate the precipitation growth scalar based on the cumulative average precipitation over the previous few days using a linear piecewise function that includes a minimum threshold and a maximum threshold. A photoperiodic scalar unit is used to calculate the photoperiodic growth scalar based on the length of day using a linear piecewise function. The shortwave radiation scalar unit is used to calculate the shortwave radiation growth scalar based on the daily average shortwave radiation intensity using a linear piecewise function.

4. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, The growing season index calculation module also includes a parameter optimization unit, which uses optimization algorithms to optimize and determine various physiological threshold parameters related to temperature, precipitation, photoperiod, and shortwave radiation growth scalars for different vegetation types, based on the phenological period extracted from flux tower sites as the baseline true value and under a cross-validation framework.

5. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, In the improved phenological index synthesis module, the coupling method is to perform daily, pixel-by-pixel multiplication of the normalized growing season index and the multi-source remote sensing vegetation index, and the time-series smoothing process adopts the sliding window moving average filtering algorithm.

6. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, In the time series reconstruction and phenology extraction module, the nonlinear annual curve fitting adopts a double logistic function model, which is a function containing six parameters: baseline offset, amplitude, spring half-saturation point, spring growth rate, autumn half-saturation point, and autumn senescence rate. The curve fitting process uses the multi-year average improved phenological index curve of the target pixel as the initial parameter value.

7. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index as described in claim 1 or 6, characterized in that, In the dynamic thresholding method, multiple candidate threshold ratios are set; the result output and verification module determines the optimal threshold ratio for extracting the photosynthetic start day and photosynthetic end day from the multiple candidate threshold ratios through accuracy verification.

8. The photosynthetic phenology inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, In the results output and verification module, the process of extracting the phenological baseline from the total primary productivity of the flux tower includes: smoothing the daily total primary productivity series with a Gaussian weighted moving average, calculating its annual maximum value, and using a fixed proportion of this maximum value as a threshold to extract the photosynthetic start date and photosynthetic end date.

9. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, The results output and verification module uses at least one statistical indicator among correlation coefficient, root mean square error, and Kling-Gupta efficiency to quantitatively evaluate the accuracy of candidate phenological products. Based on the evaluation results, the optimal solution is determined to be a combination of coupling the normalized vegetation index as the base index and the normalized growing season index, and using a combination of dual logistic function fitting and dynamic thresholding to extract phenological parameters.

10. The photosynthetic phenological inversion system combining meteorological suitability and vegetation index according to claim 1, characterized in that, The system is configured to run in a high-performance computing cluster environment, capable of processing global-scale, daily temporal resolution remote sensing and meteorological data, and ultimately generating a global photosynthetic phenology raster product with a spatial resolution of 0.05 degrees and its metadata.

Images