Potential evapotranspiration estimation method and system

Abstract

The invention provides a potential evapotranspiration estimation method and system, and the method comprises the steps: substituting corrected potential evapotranspiration observation data into a potential evapotranspiration Penman model based on a wind function, and carrying out the inversion of an equivalent wind function label. And training a machine learning model by using the equivalent wind function label, the meteorological data and the earth surface data, and outputting a nonlinear wind function by using an estimation model based on the meteorological data and the earth surface data of the target station. And replacing a wind function in the potential evapotranspiration model with a nonlinear wind function, and outputting a potential evapotranspiration estimation result based on the potential evapotranspiration model after replacement. By learning an equivalent wind function, dependence on a complex and non-migratable resistance network is partially avoided. By learning a nonlinear wind function, the regional limitation of a traditional linear coefficient is overcome, and robust generalization is guaranteed. A data-driven wind function is embedded into a potential evapotranspiration model in a physical item replacement mode instead of end-to-end black box regression, so that the traceability of model output is ensured, and the model can be butted with an existing service system.

Description

Technical Field

[0001] This invention relates to the fields of land surface processes and hydrometeorology, and more specifically, to a method and system for estimating potential evapotranspiration. Background Technology

[0002] Evapotranspiration (ET) is a crucial component of the land surface hydrological cycle, and its generation and evolution are closely linked to climate change and underlying surface cover evolution. Evapotranspiration plays a vital role not only in global and regional climate but also influencing watershed groundwater storage and surface runoff. Global climate change has exacerbated the spatiotemporal variability of extreme hydrological events and processes, leading to regional and even global water resource unevenness. Therefore, research on evapotranspiration estimation methods has become a hot topic in the field of hydrological science.

[0003] Potential evapotranspiration (PET) refers to the maximum possible evapotranspiration under adequate water supply conditions, determined by meteorological conditions and underlying surface characteristics. It is a core driver and constraint of hydrological models, and a fundamental input for drought indices (such as the SPEI) and actual evapotranspiration estimation. Classical PET models include the Penman method and its variants, which are based on the principles of energy balance and mass conservation. The wind function in the aerodynamic term of the PET model characterizes the promoting effect of atmospheric dynamics on water vapor transport and is an empirical approximation connecting wind speed and turbulent exchange intensity. Traditional wind functions often employ linear or piecewise linear forms.

[0004] Surface turbulent transport processes are nonlinearly modulated by multiple factors, including wind speed, vegetation morphology (LAI, canopy height, roughness), climate zones (temperature, humidity gradients), and atmospheric stability. Within the Penman framework, drag parameterization is the primary source of uncertainty. Existing studies, by comparing one-source and two-source PM models, have found that complex drag networks (such as soil-canopy separation and diurnal stratification) do not necessarily reduce errors (the average relative error MAPE still reaches 32–53%), and the transferability of parameter calibration is poor. This indicates that the nonlinearity and spatial heterogeneity of drag are difficult to characterize with fixed formulas. Based on global eddy covariance station data, canopy drag (rc) under non-water-constrained conditions was inverted using a modified PM equation, revealing a significant upward trend in global rc (0.43 s·m). - ¹·yr -¹), primarily driven by increased CO2 concentration, increased LAI, and changes in VPD. Using static or linear parameterization would systematically underestimate the long-term trend of PET. Existing studies, through thermal response matrix (TRM) attribution analysis, have confirmed that in the cooling effect of land greening on land surface temperature (LST), 93% of the regions show a negative correlation, and in 82% of the regions, the dominant pathway is change in aerodynamic drag (ra) rather than land surface albedo. This highlights the central role of the wind function (or equivalent aerodynamic transport term) in land-atmosphere coupling and its sensitivity to vegetation dynamics.

[0005] However, traditional linear wind functions cannot capture the aforementioned nonlinear processes, resulting in insufficient accuracy and robustness in cross-scale and cross-regional applications.

[0006] In recent years, machine learning (ML) has demonstrated feasibility in estimating ET / PET related parameters. However, existing ML work largely focuses on end-to-end ET / PET regression or single-parameter learning. End-to-end black-box regression lacks interpretability. While directly fitting PET with neural networks or ensemble models can achieve high accuracy, the internal logic of the model is opaque, making it difficult to trace the source of errors, and it cannot be integrated with existing business systems (such as Penman-based irrigation decision models). Furthermore, its generalization ability in regions with scarce data is unknown. Summary of the Invention

[0007] The purpose of this invention is to provide a potential evapotranspiration estimation method and system to overcome the regional limitations of linear systems, ensure robust generalization, and ensure that the model output is traceable and can be integrated with existing business systems.

[0008] In a first aspect, the present invention provides a method for estimating potential evapotranspiration, the method comprising: Meteorological and surface data from multiple stations are acquired, and the meteorological and surface data are preprocessed. Obtain potential state samples of the multiple sites during the water-free period; An energy closure correction is performed on the latent heat flux in the potential state sample. Based on the corrected latent heat flux, potential evapotranspiration observation data is obtained. The potential evapotranspiration observation data is substituted into the potential evapotranspiration model to obtain the equivalent wind function label. A sample set is constructed using the equivalent wind function label, meteorological data, and surface data. The constructed machine learning model is trained using the sample set to obtain the trained estimation model. Obtain meteorological and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological and surface data of the target station and the estimation model. The nonlinear wind function replaces the wind function in the aerodynamic term of the potential evapotranspiration model, and the potential evapotranspiration estimation result is output based on the replaced potential evapotranspiration model.

[0009] In a second aspect, the present invention provides a potential evapotranspiration estimation system, the system comprising: The preprocessing module is used to acquire meteorological data and surface data from multiple stations and to preprocess the meteorological data and surface data. The acquisition module is used to obtain potential state samples of the multiple sites during the water-free period; The inversion module is used to perform energy closure correction on the latent heat flux in the potential state sample, obtain potential evapotranspiration observation data based on the corrected latent heat flux, substitute the potential evapotranspiration observation data into the potential evapotranspiration model, and invert to obtain the equivalent wind function label. The training module is used to construct a sample set using the equivalent wind function label, meteorological data, and surface data, and to train the constructed machine learning model using the sample set to obtain the trained estimation model. The processing output module is used to obtain meteorological data and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological data and surface data of the target station and the estimation model. The estimation module is used to replace the wind function in the aerodynamic term of the potential evapotranspiration model with the nonlinear wind function, and output the potential evapotranspiration estimation result based on the replaced potential evapotranspiration model.

[0010] This invention provides a method and system for estimating potential evapotranspiration. It preprocesses meteorological and surface data from multiple stations to obtain potential state samples during periods of no water constraint. Energy closure correction is applied to the latent heat flux in these samples, and potential evapotranspiration observation data is obtained based on the corrected latent heat flux. The potential evapotranspiration observation data is substituted into a potential evapotranspiration model to obtain equivalent wind function labels. A sample set is constructed using the equivalent wind function labels, meteorological data, and surface data. This sample set is used to train a machine learning model, resulting in a trained estimation model. Meteorological and surface data from the target station are obtained. Based on this data, the corresponding nonlinear wind function is output using the estimation model. The nonlinear wind function replaces the wind function in the aerodynamic term of the potential evapotranspiration model, and the potential evapotranspiration estimation result is output based on the replaced potential evapotranspiration model.

[0011] This solution partially avoids reliance on complex, non-transferable drag networks by learning an equivalent wind function, thus mitigating the main source of error in the PM framework at its source. By learning a nonlinear wind function, it overcomes the regional limitations of traditional linear coefficients, achieving robust generalization across multiple regions and different PFTs. The data-driven wind function is embedded into the potential evapotranspiration model through physical term substitution, rather than end-to-end black-box regression, ensuring the model output is traceable and compatible with existing business systems. Attached Figure Description

[0012] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments of the present invention will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0013] Figure 1 A flowchart of a potential evapotranspiration estimation method provided in an embodiment of the present invention; Figure 2 This is a hierarchical diagram of the system architecture for potential evapotranspiration estimation in an embodiment of the present invention; Figure 3 A scatter plot comparison and evaluation diagram of the potential evapotranspiration model provided in the embodiments of the present invention and various existing models; Figure 4 A schematic diagram comparing the PET spatial distribution of the potential evapotranspiration model provided in this embodiment of the invention with existing models; Figure 5 This is a schematic diagram illustrating the characteristic importance of the wind function provided in an embodiment of the present invention. Detailed Implementation

[0014] The technical solutions of the present invention will now be described with reference to the accompanying drawings in the embodiments of the present invention.

[0015] Please see Figure 1 The following is a flowchart of a potential evapotranspiration estimation method provided in an embodiment of the present invention. The detailed steps of the potential evapotranspiration estimation method are described below.

[0016] S11: Acquire meteorological and surface data from multiple stations and preprocess the meteorological and surface data; S12, obtain potential state samples of multiple sites during the water-free period; S13, perform energy closure correction on the latent heat flux in the potential state sample, obtain potential evapotranspiration observation data based on the corrected latent heat flux, substitute the potential evapotranspiration observation data into the potential evapotranspiration model, and invert to obtain the equivalent wind function label; S14. The equivalent wind function label, meteorological data and surface data constitute a sample set. The machine learning model is trained using the sample set to obtain the trained estimation model. S15: Obtain meteorological and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological and surface data of the target station and the estimation model. S16 replaces the wind function in the aerodynamic term of the potential evapotranspiration model with a nonlinear wind function, and outputs the potential evapotranspiration estimation result based on the replaced potential evapotranspiration model.

[0017] In this embodiment, meteorological and surface data from multiple stations can be acquired, and the data can be obtained from various public or commercial data sources. Meteorological data includes at least, for example, wind speed at a height of 2 meters. Temperature Air humidity (or water vapor pressure) Data include atmospheric pressure (Pa), net radiation (Rn), soil heat flux (G), and evapotranspiration; surface data should include at least leaf area index (LAI) and land cover type (IGBP).

[0018] Specifically, for example, the following data can be obtained through ERA5-Land (ECMWF fifth-generation reanalysis land surface product): spatial resolution 0.1°×0.1°, temporal resolution hourly, variables including 2 m temperature, 2 m dew point temperature, 10 m wind speed (converted to 2 m), net surface radiation, soil heat flux, etc. Similar data can be obtained through GLDAS (Global Land Data Assimilation System): spatial resolution 0.25°×0.25°, temporal resolution 3 hours or day, providing similar variables.

[0019] In addition, the required data can also be obtained through other data sources, such as MERRA-2, NCEP / NCAR reanalysis, and regional meteorological station observation interpolation products.

[0020] Surface data can be obtained in the following ways: MODIS LAI (MOD15A2H / MYD15A2H): Spatial resolution 500 m, temporal resolution 8 days, used after quality control and spatiotemporal interpolation.

[0021] GLASS LAI: Spatial resolution 0.05°, temporal resolution 8 days, generated based on multi-source fusion algorithm, can be used as a replacement or supplement to MODIS.

[0022] MCD12Q1 (MODIS Land Cover Type): Annual product, spatial resolution 500 m, offering multiple classification schemes such as IGBP.

[0023] Preprocessing of acquired meteorological and surface data can be achieved in the following ways: Spatiotemporal alignment of meteorological and surface data from different data sources is performed, and gridding is applied based on the corresponding station information. Quality assessment of meteorological and surface data is conducted, and quality optimization is performed based on the assessment results. Unit conversion and standardization are performed on meteorological and surface data.

[0024] The steps of conducting quality assessments on meteorological and surface data, and then optimizing the quality of these data based on the assessment results, can be implemented in the following ways: Missing items in meteorological and surface data are interpolated and filled using a preset interpolation method; data exceeding a preset range in meteorological and surface data are filtered out or replaced based on preceding and following data; closure rate is calculated for meteorological and surface data, and meteorological or surface data with closure rates exceeding a preset threshold are marked as suspicious data. Suspicious data are removed or have their weight reduced during inversion processing.

[0025] Specifically, in this embodiment, the spatial resolution of meteorological and surface data from different data sources can be unified to a target resolution (such as a 0.1° or 0.25° grid), achieved using bilinear interpolation or nearest neighbor resampling. Example code is as follows: Python from rasterio.warp import reproject, Resampling reproject(source_array, destination_array, src_transform=src_transform, dst_transform=dst_transform, src_crs=src_crs, dst_crs=dst_crs, resampling=Resampling.bilinear) Then, time scale alignment is performed. If the target scale is daily PET, the hourly or 3-hourly data is aggregated into daily averages (temperature, humidity) or daily cumulatives (radiation, flux). If the target scale is monthly PET, it is further aggregated into monthly averages or monthly cumulatives.

[0026] When performing quality control on meteorological and surface data, for example, missing items in the meteorological and surface data can be filled in using interpolation methods such as temporal interpolation (linear interpolation, cubic spline interpolation) or spatial neighborhood interpolation. Furthermore, if the proportion of defective items is too high, for example, above 30%, the time period or grid is judged as low quality, triggering a fallback strategy during subsequent inference.

[0027] In outlier removal, reasonable ranges can be set for each data point, as shown below: ; Data that exceeds the preset range can be filtered out (set as missing data), or it can be interpolated and replaced based on the data before and after it.

[0028] In the energy closure check, the closure rate can be calculated. If the closure rate exceeds a preset threshold, the half-hour / hour data will be marked as suspicious data and will be removed or have its weight reduced when the wind function label is retrieved in the future.

[0029] The conversion and standardization process mainly involves the following aspects: 1. Wind speed to altitude conversion: ERA5-Land provides 10 m wind speed Converted to 2 m;

[0030] in, The roughness length can be found in a table based on the land cover type (e.g., farmland). m, forest m) or dynamic calculation (related to LAI).

[0031] 2. Radiation unit: Convert W / m² to MJ·m - ²·d - ¹: .

[0032] 3. Feature standardization (for ML training): Continuous variables (e.g.) Z-score standardization is performed. Or Min-Max normalization: ; Categorical variables (IGBP) use One-Hot encoding or embedding.

[0033] In addition, potential state samples during the water-limited period of multiple sites are obtained. In this embodiment, the water-limited period of each site is determined based on energy balance and soil moisture information. Specifically, this can be achieved in the following ways: For each site, observation days with net radiation flux greater than soil heat flux were identified, and observation days with precipitation were removed. For the remaining observation days, water-free stress periods were screened based on evaporation fraction as an energy balance index, and water-free limitation periods were determined based on soil moisture within the water-free stress periods. Potential state samples within the water-free limitation periods were obtained using a flux tower.

[0034] In this embodiment, firstly, at each site, only observation days with net radiation flux (Rn) greater than soil heat flux (G) are retained, and all precipitation days are removed to avoid the impact of rainwater interception on latent heat flux (LE). Secondly, the evaporation fraction (EF = LE / (LE+H)) is used as the energy balance index. When the EF value of a certain day exceeds the 95th percentile of the EF distribution of that site, it is determined to be a water-free stress day. If the number of days meeting this condition does not exceed 15 days, the 15 days with the highest EF values ​​are taken as the water-free restriction period. Finally, during the water-free stress period, soil moisture must be greater than 50% of the maximum soil moisture (98th percentile of the soil moisture sequence) at that site. In this way, the water-free restriction period is determined, and a flux tower is used to obtain a sample of the potential state during the water-free restriction period.

[0035] The potential evapotranspiration state corresponds to a situation where "moisture supply is sufficient and evapotranspiration is mainly controlled by energy and aerodynamic conditions." Specifically, in this embodiment, it is achieved through a multi-criteria joint determination method. The specific criteria and their physical meanings are as follows: (Criterion 1) Upper quantile of evapotranspiration fraction (EF):

[0036] in, The upper quantile threshold of EF, preferred range In the example, take = 0.8, meaning that the top 20% of samples with the highest EF are retained.

[0037] Physical meaning: A high EF indicates that the latent heat accounts for a large proportion of the available energy, and that moisture has a smaller limitation.

[0038] (Guideline 2) Daytime conditions:

[0039] in, For incident shortwave radiation, Preferred concentrations are 200 W / m² (half-hour average) or 5 MJ·m. - ²·d -¹(Daily cumulative).

[0040] Physical significance: Eliminating low radiation periods at night and during dawn and dusk ensures active transpiration of vegetation.

[0041] (Guideline 3) Sufficient precipitation or soil moisture: Option A (based on precipitation):

[0042] in, For the past The cumulative precipitation over 7 days (e.g., d=7 days) is the median of regional precipitation climatology or 10 mm.

[0043] Option B (based on soil moisture):

[0044] in, These are the wilting point and field water holding capacity, respectively.

[0045] Physical significance: Ensures sufficient soil moisture supply to support potential evapotranspiration.

[0046] (Guideline 4) Exclude extreme drought and freezing:

[0047] Preferred kPa, °C, avoid extreme drought or permafrost conditions.

[0048] Example code is as follows: Python def identify_potential_state(data_df, ef_quantile=0.8, sw_threshold=200, precip_days=7, precip_threshold=10): """ Potential evapotranspiration state identification Parameters: ----------- data_df : pd.DataFrame Includes columns: LE, H, SWin, precip, Ta, VPD, etc. ef_quantile : float EF upper quantile threshold sw_threshold : float Shortwave radiation threshold (W / m²) precip_days : int cumulative number of precipitation days precip_threshold : float Cumulative precipitation threshold (mm) Returns: -------- potential_mask : pd.Series (bool) True indicates a potential state. """ # Calculate EF data_df['EF'] = data_df['LE'] / (data_df['LE'] + data_df['H']) # Rule 1: Upper quantile of EF ef_thresh = data_df['EF'].quantile(ef_quantile) mask1 = data_df['EF']>= ef_thresh # Guideline 2: Daytime mask2 = data_df['SWin']>sw_threshold # Guideline 3: Sufficient Rainfall data_df['precip_cumsum'] = data_df['precip'].rolling(window=precip_days).sum() mask3 = data_df['precip_cumsum']>precip_threshold # Rule 4: Eliminate extremes mask4 = (data_df['VPD']<3.0)&(data_df['Ta']>0) # comprehensive potential_mask = mask1&mask2&mask3&mask4 return potential_mask Based on the above, an energy closure correction is performed on the latent heat flux in the latent state sample. Specifically, this can be achieved in the following way: An energy balance relationship is constructed based on the net radiation, soil heat flux, latent heat flux, and sensible heat flux in the potential state sample. Based on the energy balance relationship, the unclosed energy is proportionally allocated to the latent heat flux and sensible heat flux to obtain the energy closure corrected latent heat flux.

[0049] Eddywise observations commonly suffer from energy non-closure issues, i.e., measured... The common closure ratio is between 0.7 and 0.9. To obtain a more accurate estimate of potential evapotranspiration, a closure correction needs to be performed on the latent heat flux (LE).

[0050] In this embodiment, various closure correction methods can be used to perform closure correction. The various closure correction methods and their physical meanings are as follows: (Method 1) Bowen ratio closure correction (β method):

[0051]

[0052]

[0053] Physical meaning: Assuming that the Bowen ratio is unaffected by closure error, it distributes unclosed energy proportionally.

[0054] (Method 2) Forced Closure Method:

[0055] Physical meaning: Directly scale up LE proportionally to satisfy energy balance.

[0056] (Method 3) Dual-path uncertainty: Simultaneously using method 1 and method 2, we obtain and , forming intervals Subsequent inversion When the interval estimate is obtained It is used for uncertainty quantification.

[0057] Based on this, the corrected potential evapotranspiration observation data are substituted into the potential evapotranspiration model to obtain the equivalent wind function label.

[0058] Based on the Penman formula in the potential evapotranspiration model, the following inversion formula can be obtained:

[0059] in, This represents the slope of the saturated water vapor curve. Represents the wet / dry constant. Indicates net radiation. Indicates soil heat flux. Indicates the latent heat of vaporization. Represents the wind function. Indicates saturated water vapor pressure. Indicates the actual water vapor pressure. Represents the coefficient.

[0060] The equivalent wind function label obtained by inversion Outliers may occur (e.g., a very small denominator when VPD is close to 0, or an excessively large energy closure correction). The following strategies can be used to perform robust regression or quantile pruning on the equivalent wind function labels obtained from the inversion to remove outliers and ensure that 0 < VPD < 0. obs<10: 1. Physical boundary check: like Discard the sample (the physical wind function should be positive); like (e.g., a coefficient of 10 times FAO24) is considered abnormal and discarded.

[0061] 2. Reject if VPD is too small: when When the denominator is unstable at kPa, the sample is discarded.

[0062] 3. Segmentation and cropping: reserve

[0063] For samples within the range, extreme values ​​are removed.

[0064] Example code is as follows: Python def invert_fobs(data_df, lambda_=2.45): """ The equivalent wind function f_obs is retrieved from observations. Parameters: ----------- data_df : pd.DataFrame Includes columns such as: LE_corrected, Rn, G, Delta, gamma, VPD, etc. lambda_ : float Latent heat of vaporization (MJ / kg) Returns: -------- f_obs : pd.Series """ # PET_obs (mm / day) PET_obs = data_df['LE_corrected'] / lambda_ # LE units must be MJ / m² / d # PET_rad (mm / day) PET_rad = (data_df['Delta'] / (data_df['Delta']+ data_df['gamma'])) \\ (data_df['Rn'] - data_df['G']) / lambda_ # Inversion of f_obs numerator = (PET_obs - PET_rad) (data_df['Delta'] + data_df['gamma']) denominator = data_df['gamma'] data_df['VPD'] f_obs = numerator / denominator # Outlier Handling f_obs = f_obs[(f_obs>0)&(f_obs<f_obs.quantile(0.99))&(data_df['VPD']> 0.2)] return f_obs Based on the above, an equivalent wind function label, meteorological data, and surface data are combined to form a sample set, which is then used to train the constructed machine learning model. This sample set can be applied to global FLUXNET sites, with hundreds to thousands of valid samples obtained from each site, which are then aggregated to form the sample set.

[0065] Specifically, in this embodiment, the machine learning model can be trained to obtain the estimation model in the following manner: The sample set is divided into a training set and a test set, which include wind speed at a specific height, leaf area index, and land cover type codes. The machine learning model is trained using the training set under set constraints, including a monotonically increasing constraint on wind speed at a specific height, an output boundary constraint on the output of the machine learning model, and an extrapolation distance monitoring constraint on each training sample in the training set. The machine learning model is evaluated using the test set according to set evaluation indicators until it meets the preset requirements, thus obtaining the estimation model trained by the machine learning model.

[0066] In this embodiment, { LAI, IGBP, [region / seasonal indication, latitude / longitude or climate zone, stability surrogate (these parameters are additional, unnecessary redundancy) (e.g., Richardson number)] Using )]} as input features, and the equivalent wind function labels obtained above... obs are the labels, and the training pairs are... A monotonically increasing machine learning model constrained by physical boundaries outputs a nonlinear wind function. .

[0067] First, the input features are processed. For example, for IGBP, the following methods can be used: Scheme A (One-Hot encoding: expanding IGBP=1,...,17 into a 17-dimensional binary vector; Option B (Embedding Layer): Map the neural network to a low-dimensional continuous space (e.g., 8-dimensional) using an embedding layer; Option C (Target Encoding): Using f for each IGBP category obs Replace the IGBP values ​​with the mean (to prevent overfitting).

[0068] For different climate zones, the following methods are used: Similar to One-Hot or ordinal encoding (A=1, B=2, ...).

[0069] Machine learning models include, but are not limited to, the following categories: Random Forest (RF): Captures nonlinearity and interactions through ensemble learning; Gradient Boosting Decision Trees (GBDTs): such as XGBoost and LightGBM, can be subject to monotonic constraints. Deep Neural Networks (DNNs): Employ monotonic activation functions or force monotonicity through post-processing.

[0070] Constraint mechanisms are set during model training, including monotonically increasing constraints on wind speed at specific heights, output boundary constraints, and extrapolation distance monitoring constraints: right Monotonically increasing constraints: Set `monotone_constraints={'u2': 1}` in GBDT, or use positive weights and ReLU activation in DNN to ensure... ; Output boundary limits: setting ,in Determined based on a physically reasonable range (e.g., referring to the coefficient range of FAO24 / KP1 / KP2) or the quantiles of the training data; Boundary extrapolation distance monitoring constraint: For inputs outside the training set (such as extreme high wind speeds or abnormal LAI), limit the extrapolation range or trigger a fallback strategy.

[0071] During model training, feature interactions can be set, for example, tree models (RF / GBDT) can automatically learn interactions (such as... No manual construction is required. If a linear model or shallow DNN is used, interaction terms can be added manually.

[0072] The following are exemplary code examples of training machine learning models under different types of machine learning models and different constraint mechanisms: (Preferred Solution A) Gradient Boosting Tree (GBDT) + Monotonic Constraints: Example code (Python + LightGBM): Python import lightgbm as lgb # Features: ['u2', 'LAI', 'IGBP', 'lat', 'lon', 'month_sin', 'month_cos',...] # Tag: f_obs # Building the training set train_data = lgb.Dataset(X_train, label=y_train, feature_name=['u2', 'LAI', 'IGBP', 'lat', 'lon', ...], categorical_feature=['IGBP']) # Parameter Settings params = { 'objective': 'regression', 'metric': 'rmse', 'boosting_type': 'gbdt', 'num_leaves': 63, 'learning_rate': 0.05, 'feature_fraction': 0.8, 'bagging_fraction': 0.8, 'bagging_freq': 5, 'monotone_constraints': [1, 0, 0, 0, 0, ...], # u2 is monotonically increasing (1), the rest are unconstrained (0) 'max_bin': 255 } # train model = lgb.train(params, train_data, num_boost_round=500, valid_sets=[valid_data], early_stopping_rounds=50) Key parameters: `monotone_constraints`: The list length is the same as the number of features; 1 indicates monotonically increasing, -1 indicates monotonically decreasing, and 0 indicates no constraints. Here, `monotone_constraints` is used to specify the constraints. Set it to 1.

[0073] (Preferred Option B) Random Forest (RF): The RF model is implemented using the scikit-learn library in Python. The main hyperparameters are set as follows: The more decision trees (n_estimators), the higher the model's stability, but the higher the computational cost; in this embodiment, it is set to 100. Storing out-of-bag (OOB) data can be used for subsequent evaluation; in this embodiment, out-of-bag scoring is enabled (oob_score=True). The maximum feature (max_features) determines the number of features randomly selected to find the best split when splitting a node in each decision tree; it affects the model's diversity, training time, and generalization ability, and is set to 0.5. The criterion parameter defines the evaluation criterion for splitting nodes; it determines how to measure the "quality" of each split, i.e., the standard for selecting the optimal split, and is set to "squared_error". The random state (random_state) ensures the reproducibility of the randomness of the sampling process and is set to 42.

[0074] RF itself does not directly support monotonic constraints, but post-processing can be used: 1. Train a RandomForest (RF) model (such as `RandomForestRegressor` from sklearn); 2. During prediction, for each sample, keeping other features constant, scan... Examine the predicted values ​​from minimum to maximum. Is it monotonically increasing? 3. If monotonicity is violated, use isotonic regression.

[0075] Monotonicization of the curve Alternatively, use monotonic RF variants (such as the `monot` and `monmlp` packages).

[0076] (Option C) Deep Neural Network (DNN) + Monotonic Activation: Example code is as follows (PyTorch): Python import torch import torch.nn as nn class MonotonicDNN(nn.Module): def __init__(self, input_dim): super().__init__() self.fc1 = nn.Linear(input_dim, 64) self.fc2 = nn.Linear(64, 32) self.fc3 = nn.Linear(32, 1) def forward(self, x): # x: [batch, input_dim], assuming x[:, 0] is u2 u2 = x[:, 0:1] other_features = x[:, 1:] # Force the weights of u2 to be positive to ensure monotonicity. h1 = torch.relu(self.fc1(other_features) + torch.abs(self.fc1.weight[:, 0:1]) @ u2.T) h2 = torch.relu(self.fc2(h1)) out = self.fc3(h2) return out In practice, GBDT is more commonly used because it is easy to implement monotonic constraints and has strong interpretability.

[0077] In the above output boundary constraints, after training is completed, during the prediction phase... Perform the following cropping:

[0078] The boundary values ​​are set as follows: Avoid zero or negative values; Reference FAO24 coefficient The limit has been slightly relaxed to 10.

[0079] In the above extrapolation distance monitoring constraints, for each inference sample Calculate its distance from the sample set:

[0080] like (e.g., quantiles of a sample set) Mark it as "extrapolation" to lower the confidence level or trigger a fallback.

[0081] In this embodiment, the training set in the sample set is used to train the model in the manner described above, and the training results are tested and evaluated using the test set. The model is evaluated under multiple evaluation indicators, including: (1) the coefficient of determination (R²). 2 This measures whether the model can reduce errors and capture changing trends. R 2 (1) The higher the value, the better the model performance. (2) Root Mean Square Error (RMSE) indicates how much error the model will produce in the simulation. The larger the error, the larger the RMSE. (3) Mean Bias (MB) describes the difference between the expected value of the predicted value (estimated value) and the true value. The larger the MB, the more it deviates from the true data. (4) Taylor Skill Score (TSS) combines the correlation coefficient (R) and standard deviation to evaluate the linear relationship and variability between two variables. TSS ranges from 0 to 1, and the higher the value, the better the model performance. (5) Since it is more complicated to compare models using multiple individual indicators, the Global Performance Index (GPI) provides a unified comprehensive scoring method. R 2 TSS can assess the linearity of the model, while RMSE and MB assess the bias of the model.

[0082] During the training and testing process, each PFT and climate zone was evaluated in groups, and indicators such as R², RMSE, MB, and GPI were calculated for each group.

[0083] In addition, ablation experiments can be performed, i.e., training control models without LAI or IGBP, to quantify the contribution and necessity of each feature.

[0084] During training, you can also use the sample set of one station as the test set each time, and perform training on the sample sets of the remaining stations, looping through all stations to evaluate the spatial generalization ability. An example code is as follows: Python from sklearn.model_selection import LeaveOneGroupOut logo = LeaveOneGroupOut() for train_idx, test_idx in logo.split(X, y, groups=site_ids): X_train, X_test = X[train_idx], X[test_idx] y_train, y_test = y[train_idx], y[test_idx] # Training the model model.fit(X_train, y_train) # Evaluate y_pred = model.predict(X_test) score = evaluate(y_test, y_pred) In addition, after training, R², RMSE, MB, and GPI are calculated by IGBP or Köppen grouping, as shown in the example code below: Python for pft in [1, 2, ..., 17]: # IGBP type mask = (X_test[:, igbp_idx] == pft) y_true_pft = y_test[mask] y_pred_pft = y_pred[mask] r2_pft = r2_score(y_true_pft, y_pred_pft) rmse_pft = np.sqrt(mean_squared_error(y_true_pft, y_pred_pft)) print(f"PFT {pft}: R²={r2_pft:.3f}, RMSE={rmse_pft:.3f}") In ablation experiments, this can be implemented using the following exemplary methods: Baseline: Used only Train a linear model (reproducing FAO24); +LAI: Add LAI features; +IGBP: Further add IGBP; +Monotonic constraints: Enable monotonic constraints in GBDT; Full model: All features + monotonic constraints.

[0085] Finally, the performance of each model is compared to quantify the contribution of each component.

[0086] In this embodiment, after training the machine learning model to obtain the estimation model, the estimation model can also be validated, which can be achieved in the following ways: Based on the estimation model, partial dependency plots were plotted on wind speed at a specific height, leaf area index, and land cover type coding to verify whether the modulation of wind speed at a specific height, leaf area index, and land cover type coding by the estimation model conforms to physical expectations; the SHAP framework was called to calculate the contribution of changes in wind speed at a specific height, leaf area index, and land cover type coding to the output of the estimation model.

[0087] In this embodiment, feature importance verification can be performed, as illustrated in the following example code: Python importance = model.feature_importance(importance_type='gain') feature_names = ['u2', 'LAI', 'IGBP', 'lat', 'lon', ...] for name, imp in zip(feature_names, importance): print(f"{name}: {imp}") Results: u_2, LAI, and IGBP are expected to be important features.

[0088] In addition, partial dependency graph (PDP) drawing is performed, as shown in the following example code: Python from sklearn.inspection import partial_dependence, plot_partial_dependence fig, ax = plt.subplots(1, 3, figsize=(15, 4)) plot_partial_dependence(model, X_train, features=['u2', 'LAI', 'IGBP'], ax=ax) plt.show() The following methods were used to verify whether the modulation of the estimation model for wind speed, leaf area index, and land cover type coding at a specific height conforms to the following physical expectations: The PDP should be monotonically increasing (verifying the effectiveness of the monotonic constraint); the PDP of LAI should reflect the nonlinear effect of vegetation on roughness / aerodynamic transmission; different IGBTs The mean should vary (e.g., forest > grassland > farmland).

[0089] In addition, the SHAP framework is called to calculate the contribution of each feature change to the output of the estimation model, as shown in the following example code: Python import shap explainer = shap.TreeExplainer(model) shap_values ​​= explainer.shap_values(X_test) shap.summary_plot(shap_values, X_test, feature_names=feature_names) SHAP plots show the contribution of each feature to each prediction, helping to diagnose whether the model has learned reasonable physical relationships.

[0090] Based on the estimation model obtained through training, the meteorological and surface data collected from the target site are processed using the estimation model to output the corresponding nonlinear wind function.

[0091] Specifically, the target site can collect global grid data (such as ERA5-Land 0.1° grid), and perform the following processing on a pixel-by-pixel and time-by-time basis to obtain the nonlinear wind function: 1. Input preparation: Reading ; 2. Calculate auxiliary variables: ; 3. Model Prediction: ; The obtained nonlinear wind function is used to replace the wind function in the aerodynamic term of the potential evapotranspiration model, and potential evapotranspiration estimation results are generated according to the daily or monthly scale.

[0092] The potential evapotranspiration model after the replacement is as follows:

[0093] Finally, the output potential evapotranspiration estimation results include potential evapotranspiration raster data (main product, in mm / day or mm / month), nonlinear wind function value raster data (intermediate product, for easy diagnosis and comparison), quality assurance indicators (QA indicators, recording input integrity (whether there are missing measurements), extrapolation distance (the distance between the input and the feature space of the training samples), fallback path (whether backoff is triggered)) and model uncertainty indicators (such as the standard deviation of the ensemble model, or the quantile range (e.g., 10%–90%)).

[0094] Example code is as follows (Python + xarray): Python import xarray as xr import numpy as np # Assume that 2D arrays of PET, f_hat, qa, and uncertainty have been calculated. ds = xr.Dataset({ 'PET': (['lat', 'lon'], PET), 'f_hat': (['lat', 'lon'], f_hat), 'QA': (['lat', 'lon'], qa), 'uncertainty': (['lat', 'lon'], uncertainty) }, coords={'lat': lats, 'lon': lons}) ds['PET'].attrs = {'units': 'mm / day', 'long_name': 'PotentialEvapotranspiration'} ds['f_hat'].attrs = {'units': 'm / s', 'long_name': 'Learned windfunction'} ds['QA'].attrs = {'long_name': 'Quality flag', 'flag_values': [0,1,2,3,4], 'flag_meanings': 'normal KP_fallback PT_fallback missing_inputextrapolation'} ds.to_netcdf('PETws_daily_20200101.nc') Based on the above, the potential evapotranspiration estimation method provided in this embodiment further includes the following steps: Determine whether the preset backoff mechanism is met. The preset backoff mechanism is determined when the meteorological or surface data input to the estimation model is missing, or the meteorological or surface data input to the estimation model exceeds a preset range, or the distance between the meteorological or surface data input to the estimation model and the feature space of the training samples of the estimation model exceeds a preset threshold, or the estimation scenario is a preset scenario. If the preset backoff mechanism is met, the potential evapotranspiration estimation result is output using the original potential evapotranspiration model containing the wind function.

[0095] In this embodiment, the QA flag indicates the fallback level for each site, for example: 0: Normal. Valid; 1: First-level catch-all (KP1 / KP2 / FAO24); 2: Second-level catch-all (PT); 3: Input missing; 4: Extrapolation distance exceeds threshold.

[0096] Among them, input missing, input exceeding the preset range, input extrapolation distance exceeding the threshold, and estimated scenario being a preset scenario all meet the preset rollback mechanism.

[0097] (Triggering Condition 1) Input Missing: If If any variable in LAI or IGBP is missing, prediction is impossible. This triggers a first-level safety net.

[0098] (Triggering Condition 2) Input Out of Bounds: If (e.g., 20 m / s) or (e.g., 10 m² / m²), exceeding the training range, triggers an extrapolation warning or fallback.

[0099] (Triggering condition 3) Extrapolation distance exceeds threshold: Calculate input The minimum distance d from the training set, if Low confidence level triggers a fallback mechanism.

[0100] (Triggering Condition 4) Extreme drought with high VPD: If VPD > 3 kPa and The aerodynamic term may be unstable, and a rollback of PT is allowed.

[0101] Level 1 safety net: Revert to Wright KP1 / KP2 or FAO24, and calculate PET using a fixed coefficient; Secondary safety net: In extremely arid and high VPD scenarios (e.g., VPD > 3 kPa and...) Allows a fallback to the Priestley-Taylor method. This avoids the instability of the aerodynamic terms.

[0102] Example code is as follows: Python def compute_PET_with_fallback(u2, LAI, IGBP, Rn, G, Delta, gamma,VPD, lambda_, SM=None): qa = 0 # Initialize QA to normal. # Check input integrity if np.isnan(u2) or np.isnan(LAI) or np.isnan(IGBP): qa = 3 # Input missing PET = compute_PET_KP1(u2_filled, VPD, Rn, G, Delta, gamma, lambda_) f_hat = np.nan return PET, f_hat, qa # Check extrapolation if u2>20 or LAI>10: qa = 4 # Extrapolation PET = compute_PET_FAO24(u2, VPD, Rn, G, Delta, gamma, lambda_) f_hat = np.nan return PET, f_hat, qa # Normal Prediction f_hat = model.predict([[u2, LAI, IGBP, ...]])[0] # Check extreme scenarios if VPD>3.0 and SM is not None and SM<0.1: qa = 2 # Secondary catch-all PT PET = compute_PET_PT(Rn,G,Delta,gamma,lambda_,alpha=1.26); return PET , f_hat , q # New Penman dress PET_rad = (Delta / (Delta + gamma)) (Rn - G) / lambda_ PET_aero = (gamma / (Delta + gamma)) f_hat VPD PET = PET_rad + PET_aero return PET , f_hat , q def compute_PET_KP1(u2, VPD, Rn, G, Delta, gamma, lambda_, season='summer'): # Wright KP1 snowflake(s) if season == 'summer': yes, bw = 0.5, 0.6 else: yes, bw = 1.0, 0.54 f_KP1 = yes + bw u2 PET_rad = (Delta / (Delta + gamma)) (Rn - G) / lambda_ PET_aero = (gamma / (Delta + gamma)) f_KP1 VPD return PET_rad + PET_aero def compute_PET_FAO24(u2, VPD, Rn, G, Delta, gamma, lambda_): f_FAO24 = 1.0 + 0.54 u2 PET_rad = (Delta / (Delta + gamma)) (Rn - G) / lambda_ PET_aero = (gamma / (Delta + gamma)) f_FAO24 VPD return PET_rad + PET_aero def compute_PET_PT(Rn, G, Delta, gamma, lambda_, alpha=1.26): # Priestley-Taylor return alpha (Delta / (Delta + gamma)) (Rn - G) / lambda_ In this embodiment, the model uncertainty indicators in the potential evapotranspiration estimation results mainly include the standard deviation of the ensemble model, quantile regression, and the interval of the two-path closure correction. The specific information is as follows: Ensemble model standard deviation: When training multiple models (such as RF or GBDT with different random seeds), the mean and standard deviation of the output are calculated during prediction.

[0103]

[0104] Will Transmission to PET:

[0105] Quantile regression: Train a quantile regression model (such as Quantile GBDT) to output the 10th, 50th, and 90th percentiles, forming the prediction interval.

[0106] The interval for the two-path closure correction: using the equivalent wind function labels obtained above. interval During training The range is defined, and the upper and lower bounds of PET are output during inference.

[0107] In this embodiment, to implement the aforementioned potential evapotranspiration estimation method, a corresponding system architecture is provided. This system architecture mainly includes four layers, such as... Figure 2 As shown, it includes a data layer, a feature and label layer, a model layer, and an application layer. The key modules and their functional settings in each layer are as follows: Data layer: (Module 1) Data Fetcher: Functionality: Automatically downloads required data from public data sources (CDS, NASA EarthData, FLUXNET portal). Technology Stack: Python + cdsapi (ERA5) / wget / requests. Output: Raw NetCDF / HDF5 files.

[0108] (Module 2) Regridding Engine: Function: Unifies data of different resolutions to a target grid (e.g., 0.1° or 0.25°). Technology Stack: CDO (Climate Data Operators) / xESMF (Python xarray-based) / rasterio. Methods: Conservative resampling, bilinear resampling, nearest neighbor resampling.

[0109] (Module 3) Quality Control Module: Functions: Outlier detection, missing value imputation, energy closure check. Algorithms: Range check, temporal interpolation (linear / cubic spline), spatial interpolation (inverse distance weighting / Kriging).

[0110] (Module 4) Caching and Indexing: Function: Stores preprocessed data in an efficient format, supporting fast slice reading. Technology Stack: HDF5 (with chunking) / Zarr / Parquet (for tabular data). Advantages: Avoids redundant preprocessing, accelerating training and inference.

[0111] Feature and label layers: (Module 5) Potential State Identifier: Input: FLUXNET site data (half-hour / hourly LE, H, SWin, precip, Ta, VPD, SM, etc.). Output: Boolean mask, marking potential state samples. Implementation: Joint decision based on multiple criteria.

[0112] (Module 6) Energy Closure Corrector: Input: Original LE, H, Rn, G. Output: Corrected LE_corrected (two output methods are optional). Implementation: Bowen's ratio method, forced closure method.

[0113] (Module 7) f_obs Inverter: Input: LE_corrected, Rn, G, Delta, gamma, VPD, etc. Output: f_obs label. Implementation: The inversion formula shown above, including outlier filtering.

[0114] (Module 8) Sample Builder: Function: Aggregates (feature, f_obs) pairs from each site and divides them into training, validation, and test sets. Output: Training set DataFrame or NumPy array.

[0115] Model layer: (Module 9) Model Trainer: Input: Training set (X, y). Output: Trained model object (e.g., LightGBM Booster, RF pkl).

[0116] Features: Hyperparameter search (Grid Search / Bayesian Optimization); Cross-validation (LOSO / K-Fold); Monotonic constraint setting; Model saving (pickle / joblib / ONNX).

[0117] (Module 10) Model Inference Engine: Input: Global grid data (u2, LAI, IGBP, ...). Output: f_hat raster. Mode: Batch inference: Processing a whole month / year of data at once; Online inference: Real-time API service. Optimization: Multiprocessing / GPU acceleration (if DNN).

[0118] (Module 11) QA & Fallback Controller: Input: Raw input data + f_hat prediction. Output: QA flag grid. Logic: Check input integrity, extrapolation distance, extreme scenarios, and trigger fallback.

[0119] (Module 12) Uncertainty Quantifier: Input: Multiple predictions from the ensemble model. Output: Uncertainty raster (standard deviation / quantile range). Method: Standard deviation or quantile range of the ensemble model.

[0120] (Module 13) Model Versioning: Functionality: Manages different versions of a model (v1.0, v1.1, ...), recording training configurations, data sources, and performance metrics. Technology Stack: MLflow / DVC / Weights & Biases. Advantages: Reproducibility, A / B testing, and rollback capabilities.

[0121] Application layer: (Module 14) Product Publisher: Input: PET, QA, uncertainty raster. Output: NetCDF / GeoTIFF file with complete metadata (coordinates, units, timestamps, attributes).

[0122] Metadata example (CF-compliant): PET:units = "mm day-1" PET:long_name = "Potential Evapotranspiration estimated by PETws" PET:standard_name = "water_potential_evapotranspiration_flux" PET:grid_mapping = "crs" PET:valid_range = 0.0, 20.0 (Module 15) Visualization Service: Functionality: The web front-end displays the spatiotemporal distribution, comparison charts, and time series of PET (Peak Emissions). Technology Stack: Flask / Django (backend) + Leaflet / OpenLayers (maps) + Plotly (charts).

[0123] (Module 16) API Service: Function: Provides a RESTful or gRPC interface for external systems to query PET.

[0124] Example endpoint: `GET / api / v1 / pet?lat=40.0&lon=116.0&date=2020-01-01` → Returns the PET value for that point. POST / api / v1 / pet / batch → Batch query of multiple points / time periods Technology stack: FastAPI / Flask-RESTful.

[0125] (Module 17) Monitoring & Logging: Functionality: Real-time monitoring of system operation status, recording of error logs, and statistical analysis of QA flag distribution. Technology Stack: Prometheus (metrics) + Grafana (dashboard) + ELK Stack (logging).

[0126] The system framework provided in this embodiment is deployed using the following deployment scheme: Offline batch processing mode: Suitable for generating global PET products for historical periods.

[0127] Workflow: 1. Data preparation: Download and preprocess ERA5-Land, MODIS LAI, and MCD12Q1 data for the specified period (e.g., 2000–2020).

[0128] 2. Batch inference: Multiple processes run in parallel, with each process handling a number of time steps or space blocks.

[0129] 3. Product Generation: The output is summarized as annual / monthly NetCDF files.

[0130] 4. Quality Inspection: Automated scripts check product integrity and statistically analyze QA distribution.

[0131] Deployment environment: Computing clusters (such as HPC) or cloud virtual machines (AWS EC2, Google Cloud Compute).

[0132] Resource requirements: A typical configuration is a 64-core CPU + 128 GB RAM, which takes several hours to one day to process global 0.1° diurnal scale data.

[0133] Online service mode: Applicable to real-time or near real-time PET queries.

[0134] Architecture: Front-end: Web application, where users input spatiotemporal coordinates.

[0135] Backend: API server, which calls the model inference engine.

[0136] Caching: Pre-calculate PET for frequently used regions / time periods and store it in Redis or a database to speed up queries.

[0137] Model deployment: The model is packaged as ONNX or TensorFlow SavedModel and deployed on TensorFlowServing / ONNX Runtime.

[0138] Response time: Single query <100 ms; Batch query (hundreds of points) ~1 s.

[0139] Containerization and Automation: Containerizing modules using Docker facilitates cross-platform deployment: Dockerfile # Dockerfile example FROM python:3.9-slim WORKDIR / app Copy requirements.txt. RUN pip install -r requirements.txt COPY. CMD ["python", "inference.py"] Automated CI / CD process (taking GitLab CI as an example): yaml stages: - test - train - deploy test: script: - pytest tests / train: script: - python train_model.py - mlflow log-model model.pkl deploy: script: - docker build -t petws:latest . - docker push registry.example.com / petws:latest - kubectl apply -f k8s / deployment.yaml To verify the effectiveness of the potential evapotranspiration estimation method provided in this embodiment, three specific embodiments are described below.

[0140] Example 1: Global Model Training and Validation Based on FLUXNET2015 Data preparation: 1. Site selection: 152 sites were selected from the FLUXNET2015 dataset, covering 17 IGBP types (ENF, DBF, MF, WSA, SAV, GRA, WET, CRO, URB, etc.) and the main Köppen climate zones (A–E).

[0141] 2. Time range: 2001–2014, a total of 14 years of data.

[0142] 3. Variable extraction: Eddy observations: LE_F_MDS (gap-filled latent heat flux), H_F_MDS (gap-filled sensible heat flux), NETRAD (net radiation), G_F_MDS (soil heat flux), TA_F (air temperature), VPD_F (VPD), WS_F (wind speed), P_F (precipitation), SW_IN_F (incident shortwave).

[0143] Auxiliary data: Site metadata (latitude, longitude, IGBP, altitude).

[0144] 4. LAI data: LAI time series of the pixels where the station is located extracted from MODIS MOD15A2H (8-day composite, interpolated to daily scale).

[0145] Latent state identification and label inversion: For each site, perform latent state identification in step 4: EF threshold: SW threshold: 200 W / m²; Precipitation: >10 mm cumulatively over the past 7 days; VPD <3 kPa, Ta >0°C; Approximately 58,000 half-hour samples were obtained (accounting for ~15% of the total sample).

[0146] Energy closure correction: Using the Bowen ratio method, the closure rate was improved from 0.75 to 1.0.

[0147] Inverting f_obs: Samples with VPD < 0.2 kPa or f_obs < 0 were removed, resulting in 52,000 valid labels.

[0148] Model training: Features: u2, LAI, IGBP (17-dimensional One-Hot encoding), lat, lon, month_sin, month_cos, a total of 23 dimensions.

[0149] Model: LightGBM, parameters are as follows: Python params = { 'objective': 'regression', 'metric': 'rmse', 'num_leaves': 63, 'learning_rate': 0.05, 'feature_fraction': 0.8, 'bagging_fraction': 0.8, 'bagging_freq': 5, 'monotone_constraints': [1, 0, 0,...,0], # Only u2 is monotonic 'max_depth': 10, 'num_iterations': 500, 'early_stopping_round': 50 } Cross-validation: LOSO (152-fold), leaving one site out for testing each time.

[0150] Training time: Approximately 2 hours on a 64-core CPU server.

[0151] Verification results: Site-wide aggregate metrics: R² = 0.88; RMSE = 0.31 (after normalization, the original unit is approximately 0.9 m / s); MB = 0.02.

[0152] Partial PFT assessment (results): ENF (Evergreen Coniferous Forest): R²=0.85, RMSE=0.29; GRA (Grassland): R²=0.90, RMSE=0.27; CRO (Farmland): R²=0.89, RMSE=0.28; DBF (Deciduous Broadleaf Forest): R²=0.87, RMSE=0.32.

[0153] Monotonicity check: For all test samples, with other features fixed, scan u2 from 0 to 15 m / s to confirm that f_hat is monotonically increasing (no violation).

[0154] Feature importance (Gini): 1. u2: 42%; 2. LAI: 28%; 3. IGBP (merged): 18%; 4. lat: 7%; 5. month_sin / cos: 5%.

[0155] Example 2: Global 0.1° Daily Scale PET Product Generation Inference settings: Objective: To generate global 0.1° daily-scale PET products in 2020.

[0156] Input data: ERA5-Land: 2 m temperature, 2 m dew point, 10 m wind speed, net radiation, soil heat flux; MODIS LAI: MOD15A2H v006 (8-day synthesis, interpolated to day scale); MCD12Q1: 2020 IGBP classification.

[0157] Grid count: Approximately 3600 × 1800 = 6.48 million pixels globally.

[0158] Time steps: 366 days (2020 is a leap year).

[0159] Total sample size: 6.48M × 366 ≈ 237 million.

[0160] Reasoning process: Block-based parallel processing: The global network is divided into 10°×10° sub-blocks (a total of 36×18=648 blocks); each block is inferred independently, and a temporary NetCDF file is output; finally, they are merged into a global NetCDF.

[0161] Single-block inference time: Approximately 5 minutes (64 cores in parallel). Total inference time: Approximately 5 hours (including I / O).

[0162] QA Statistics: QA distribution of global pixel-to-daily samples: QA=0 (Normal): 92.3%; QA=1 (First-level fallback, KP1): 5.2% (mainly in high-latitude winters where LAI data is missing); QA=2 (Second-level fallback, PT): 1.8% (mainly in extremely arid regions such as the Sahara and Arabian Deserts); QA=3 (Input missing): 0.5%; QA=4 (Extrapolation): 0.2%.

[0163] Product verification: Comparison with FLUXNET site observations (32 sites available in 2020): PETws: R²=0.84, RMSE=0.58 mm / day, MB=0.03 mm / day; FAO24: R²=0.71, RMSE=1.12 mm / day, MB=0.52 mm / day; PETKP1: R²=0.68, RMSE=1.38 mm / day, MB=0.78 mm / day.

[0164] Comparison with FAO56 Penman-Monteith reference crop evapotranspiration (ETo) product (farmland area): Correlation: R² = 0.79; PETws is slightly higher than ETo (MB = 0.15 mm / day), which is in line with expectations (PET is usually ≥ ETo).

[0165] Example 3: Effect of Regional System Bias Correction Comparison of arid regions (Northwest China NWC): Regional definition: Xinjiang, Gansu, Ningxia, and western Inner Mongolia (35–50°N, 75–110°E).

[0166] Time period: Average from 2010 to 2020.

[0167] result: FAO56 annual average PET: 1320 mm / year; PETws annual average PET: 1180 mm / year (a decrease of 10.6%); Measured ET (water balance method): approximately 250 mm / year; Effective precipitation: approximately 200 mm / year.

[0168] Analysis: FAO56 overestimated PET under high wind speed (annual average 3.2 m / s) and high VPD (>2 kPa) conditions; PETws corrected this bias by using a nonlinear wind function, making PET closer to the reasonable upper limit of energy and moisture constraints.

[0169] Comparison of humid regions (South China SC): Regional definition: Guangdong, Guangxi, Fujian, and southern Jiangxi (22–28°N, 105–120°E).

[0170] Time period: Average from 2010 to 2020.

[0171] result: FAO56 annual average PET: 1050 mm / year; PETws annual average PET: 1120 mm / year (an increase of 6.7%); measured ET (MODIS ET product): approximately 850 mm / year; annual precipitation: approximately 1600 mm / year.

[0172] Analysis: Under conditions of sufficient moisture and high LAI (>4), PETws captured the enhancing effect of vegetation on aerodynamic transport, making PET estimates more sensitive to energy-driven effects; the fixation factor of FAO56 underestimated this effect.

[0173] In summary, among them, Figure 3 The diagram illustrates a scatter plot comparing the solution provided in this embodiment with several existing models under different vegetation types. Figure 3 The data is arranged in a 6x5 matrix. Rows represent different vegetation types (ENF, GRA, CRO, SAV, DBF, EBF, etc.), and columns represent different models (PETow, FAO24, PETKP1, PETKP2, PETws). Each subplot has the following dimensions: horizontal axis = simulated PET, vertical axis = observed PET, black dashed line = 1:1 line, red dashed line = fitted line, and labeled GPI values.

[0174] Depend on Figure 3 It can be seen that the GPI of the model (PETws) column (far right) in this scheme is all >1.1, significantly higher than other models. The scatter points are closely clustered around the 1:1 line, and R² > 0.83. The GPI of the FAO24 / PETKP1 / PETKP2 columns are mostly negative or close to zero, and the scatter points deviate significantly from the 1:1 line.

[0175] also, Figure 4 The model (PETws) provided in this embodiment schematically illustrates the annual average distribution and interannual variation of the national PET spatial distribution between 1985 and 2015. This includes (a) the annual average and (b) the interannual trend bar distributions for seven climate zones.

[0176] Potential evapotranspiration shows an increasing trend across different climatic regions. Spatially, the average annual PET gradually increases from northwest to southeast. Figure 4(a) The average annual PET (petrol oxide) levels are relatively low in the Qinghai-Tibet Plateau (QTP) and the Northwest Arid Region (NWC), at 550 mm / yr and 707 mm / yr, respectively. In contrast, the average annual PET levels are highest in Central China (CC) and South China (SC), at 975 mm / yr and 1101 mm / yr, respectively. Looking at the trend, PET levels in most regions showed an upward trend between 1985 and 2015, with the national average PET growth rate being approximately 1.8 mm / yr. 2 ( Figure 4 (b) The increase was most significant in North China (NC) and Central China (CC), with some areas experiencing growth rates exceeding 6 mm / yr. 2 The growth rate of PET in QTP and NWC regions was relatively low, at 0.8 mm / yr. 2 and 1.3 mm / yr 2 .

[0177] Figure 5 The diagram illustrates the dominant factors influencing the wind function. LAI is the most dominant factor, followed by wind speed (WS_F / u2) and IGBP. The wind function increases with increasing LAI and decreases with increasing wind speed.

[0178] Based on the results of Examples 1-3 above, it is evident that this scheme possesses significant advantages. Specifically, in Example 1: the quantitative verification results of R²=0.88 and RMSE=0.31 are significantly superior to benchmarks that do not employ this invention (such as linear wind function R²~0.6). In Example 2: the global product-site comparison GPI>1.1, while the GPI of FAO24 / PETKP1 / PETKP2 <0.5. In Example 3: the regional system bias correction is 10.6% (arid region) and 6.7% (humid region), addressing key deficiencies of existing technologies.

[0179] Furthermore, all input data in this solution is publicly available (ERA5-Land, MODIS, FLUXNET). The algorithm and code framework are open and transparent (Python + LightGBM / scikit-learn and other open-source libraries). The computational resource requirements are reasonable (a 64-core server can complete global product generation within hours). Therefore, this solution is feasible.

[0180] Furthermore, Table 1 shows the comparison between the model (PETws) provided in this embodiment and the existing models (FAO24 / FAO56) in several comparison items. While maintaining the Penman physical framework, the present invention overcomes the regional limitations of the fixed coefficients of FAO through data-driven optimization, and the measured RMSE is reduced by more than 50% (from 1.2 to 0.6 mm / day).

[0181] Table 1. Comparison of the present invention model (PETws) with existing FAO24 / FAO56

[0182] Furthermore, compared to the existing Wright KP1 / KP2 methods, the Wright KP series introduces seasonal parameters (such as summer parameters). ;winter However, it is still piecewise linear in nature, and the coefficients are derived from temperate grassland experiments.

[0183] In comparison, this scheme has the following advantages: Continuous nonlinearity: It learns the complex interaction of u2-LAI-IGBP through machine learning, rather than discrete piecewise computation. Global training: KP coefficients are only applicable to the temperate zone of North America, while this invention covers major PFT and climate zones worldwide. Quantitative validation: As shown in the attached figure, PET... KP1 / PET KP2 In most PFTs, GPI < 0, while PETws > 1.1.

[0184] In addition, compared with existing end-to-end ML PET estimation methods, existing deep learning direct regression PET methods (such as LSTM, CNN-LSTM) can achieve high accuracy, but at least have the following problems: black box problem: the internal logic is not transparent and is difficult to physically interpret; data dependence: a large number of PET labels are required (while PET itself is difficult to observe directly), and training data is limited; poor engineering adaptability: it is difficult to interface with existing Penman-based hydrological models.

[0185] In comparison, this solution has the following advantages: Physical embedding: It retains the Penman structure, only replaces the wind function submodule, and outputs intermediate quantities. It can be analyzed independently. Innovative label construction: f-labels are cleverly constructed using eddy covariance (LE) observations through latent state identification and FOBS inversion, instead of directly using PET. Interpretability: SHAP / PDP analysis confirms the physical relationships learned by the model (u2 monotonicity, LAI-enhanced roughness, etc.).

[0186] Compared to existing PT α-learning and rc-learning methods, which are both data-driven explorations of the physics submodule and learn α or rc respectively, this invention learns the wind function f (aerodynamic coefficients in the Penman framework). α / rc is often post-processed (i.e., first estimated using PM, then adjusted α / rc), while this invention directly inverts f_obs from observations as a supervisory signal. α is used in the PT method (a simplified Penman method, without aerodynamic terms), and rc is used in PM (requiring simultaneous estimation of ra and rc, resulting in high complexity). In this invention, f is directly used in Penman, achieving moderate complexity and high accuracy.

[0187] In summary, the potential evapotranspiration estimation method provided in this embodiment has at least the following beneficial effects: 1. Improve the accuracy and consistency of PET estimation at a global scale, across vegetation types and climate zones: by learning nonlinear wind functions It overcomes the regional limitations of traditional linear coefficients and achieves robust generalization in arid / semi-arid regions, humid / semi-humid regions, and under different PFTs.

[0188] 2. Maintain physical interpretability and traceability: Embed the data-driven wind function into the Penman framework by replacing physical terms, rather than using end-to-end black-box regression, to ensure that the model output is traceable and can be integrated with existing business systems.

[0189] 3. Introducing physical constraints and engineering quality control: for wind speed Apply monotonically increasing constraints to limit the output boundary, and provide QA flags and fallback strategies (fallback to KP1 / KP2 / FAO24 or PT method) to ensure stable operation of the model when the input data is abnormal or missing.

[0190] 4. Reduce the uncertainty of drag parameterization: By learning the equivalent wind function, the dependence on complex and non-transferable drag networks is partially avoided, thus mitigating the main source of error in the PM framework from the source.

[0191] This solution can be applied to multiple scenarios, including but not limited to the following: Hydrological simulation: PET serves as input to hydrological models such as VIC, SWAT, and HBV, improving the accuracy of runoff simulation. In watershed water resource assessment, accurate PET is fundamental for estimating actual ET and groundwater recharge.

[0192] Agricultural Irrigation: Replaces FAO56 ETo as a reference for calculating irrigation water demand, reducing over-irrigation in arid areas and improving water use efficiency. Supports precision agriculture decision-making systems, optimizing irrigation scheduling by combining soil moisture and crop coefficients.

[0193] Drought Monitoring and Early Warning: PET input to the SPEI (Standardized Precipitation Evapotranspiration Index) improves the accuracy of drought identification. In the context of global change, it dynamically tracks PET trends to assess the impact of changes in evaporation demand on water resources.

[0194] Based on the same inventive concept, this embodiment of the invention also provides a potential evapotranspiration estimation system. This embodiment can divide the potential evapotranspiration estimation system into functional modules according to the above method embodiments: The preprocessing module is used to acquire meteorological data and surface data from multiple stations and to preprocess the meteorological data and surface data. The acquisition module is used to obtain potential state samples of the multiple sites during the water-free period; The inversion module is used to perform energy closure correction on the latent heat flux in the potential state sample, obtain potential evapotranspiration observation data based on the corrected latent heat flux, substitute the potential evapotranspiration observation data into the potential evapotranspiration model, and invert to obtain the equivalent wind function label. The training module is used to construct a sample set using the equivalent wind function label, meteorological data, and surface data, and to train the constructed machine learning model using the sample set to obtain the trained estimation model. The processing output module is used to obtain meteorological data and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological data and surface data of the target station and the estimation model. The estimation module is used to replace the wind function in the aerodynamic term of the potential evapotranspiration model with the nonlinear wind function, and output the potential evapotranspiration estimation result based on the replaced potential evapotranspiration model.

[0195] The potential evapotranspiration estimation system provided in this embodiment can be used to execute the potential evapotranspiration estimation method under any of the above embodiments. For details not covered in this embodiment, please refer to the corresponding descriptions in the above embodiments. This embodiment will not elaborate further here.

[0196] In the embodiments provided by this invention, it should be understood that the disclosed apparatus and method can be implemented in other ways. The apparatus embodiments described above are merely illustrative. For example, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. Furthermore, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Additionally, the displayed or discussed mutual couplings, direct couplings, or communication connections may be through some communication interfaces; indirect couplings or communication connections between devices or units may be electrical, mechanical, or other forms.

[0197] Furthermore, the units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0198] Furthermore, the functional modules in the various embodiments of the present invention can be integrated together to form an independent part, or each module can exist independently, or two or more modules can be integrated to form an independent part.

[0199] It should be noted that if the functionality is implemented as a software module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0200] The above description is merely an embodiment of the present invention and is not intended to limit the scope of protection of the present invention. For those skilled in the art, the present invention can have various modifications and variations. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for estimating potential evapotranspiration, characterized in that, The method includes: Meteorological and surface data from multiple stations are acquired, and the meteorological and surface data are preprocessed. Obtain potential state samples of the multiple sites during the water-free period; An energy closure correction is performed on the latent heat flux in the potential state sample. Based on the corrected latent heat flux, potential evapotranspiration observation data is obtained. The potential evapotranspiration observation data is substituted into the potential evapotranspiration model to obtain the equivalent wind function label. A sample set is constructed using the equivalent wind function label, meteorological data, and surface data. The constructed machine learning model is trained using the sample set to obtain the trained estimation model. Obtain meteorological and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological and surface data of the target station and the estimation model. The nonlinear wind function replaces the wind function in the aerodynamic term of the potential evapotranspiration model, and the potential evapotranspiration estimation result is output based on the replaced potential evapotranspiration model.

2. The potential evapotranspiration estimation method according to claim 1, characterized in that, The method further includes: Determine whether a preset rollback mechanism is met. The preset rollback mechanism is defined as follows: the meteorological data or surface data input to the estimation model is missing; the meteorological data or surface data input to the estimation model exceeds a preset range; the distance between the meteorological data or surface data input to the estimation model and the feature space of the training samples of the estimation model exceeds a preset threshold; or the estimation scenario is a preset scenario. If the preset backoff mechanism is met, the potential evapotranspiration estimation result is output using the original potential evapotranspiration model that includes the wind function.

3. The potential evapotranspiration estimation method according to claim 1, characterized in that, The steps for preprocessing the meteorological and surface data include: Meteorological and surface data from different data sources are spatiotemporally aligned and then gridded based on the corresponding station information. The meteorological data and surface data are subjected to quality assessment, and the meteorological data and surface data are optimized based on the quality assessment results; The meteorological and surface data are subjected to unit conversion and standardization processing.

4. The potential evapotranspiration estimation method according to claim 3, characterized in that, The steps of conducting quality assessments on the meteorological and surface data, and optimizing the meteorological and surface data based on the quality assessment results, include: Missing items in the meteorological and surface data are interpolated and filled in using a preset interpolation method. Data exceeding a preset range in the meteorological and surface data are filtered out or interpolated and replaced based on the preceding and following data. Closure rate is calculated for the meteorological data and surface data. Meteorological data or surface data with closure rate exceeding a preset threshold are marked as suspicious data. Suspicious data are removed or have their weight reduced during inversion processing.

5. The potential evapotranspiration estimation method according to claim 1, characterized in that, The step of obtaining potential state samples of the multiple sites during the water-free period includes: For each of the aforementioned sites, observation days with net radiation flux greater than soil heat flux were identified, and observation days with precipitation were excluded. For the retained observation days, the water-free stress period in the observation days is screened based on the evaporation fraction as an energy balance index, and the water-free limitation period is determined based on soil moisture within the water-free stress period; A flux column was used to obtain a sample of the potential state during the water-free period.

6. The potential evapotranspiration estimation method according to claim 1, characterized in that, The step of performing energy closure correction on the latent heat flux in the latent state sample includes: An energy balance relationship is constructed based on the net radiation, soil heat flux, latent heat flux, and sensible heat flux in the potential state sample. Based on the energy balance relationship, the unclosed energy is proportionally allocated to the latent heat flux and the sensible heat flux to obtain the latent heat flux after energy closure correction.

7. The potential evapotranspiration estimation method according to claim 1, characterized in that, The step of training the constructed machine learning model using the sample set to obtain the trained estimation model includes: The sample set is divided into a training set and a test set, and the training set and test set include wind speed at a specific height, leaf area index and land cover type coding; The machine learning model is trained using the training set and under the set constraint mechanism, which includes a monotonically increasing constraint on the wind speed at the specific height, an output boundary constraint on the output result of the machine learning model, and an extrapolation distance monitoring constraint on each training sample in the training set. The machine learning model is evaluated using the test set according to the set evaluation indicators until the preset requirements are met, thus obtaining the estimation model trained by the machine learning model.

8. The potential evapotranspiration estimation method according to claim 7, characterized in that, The method further includes a step of validating the estimation model, which includes: Based on the estimation model, partial dependency plots were plotted on the wind speed, leaf area index, and land cover type coding at the specific height, respectively, to verify whether the modulation of the estimation model on the wind speed, leaf area index, and land cover type coding at the specific height conforms to physical expectations. The SHAP framework is invoked to calculate the contribution of changes in wind speed at the specified altitude, leaf area index, and land cover type coding to the output of the estimation model.

9. The potential evapotranspiration estimation method according to claim 1, characterized in that, The potential evapotranspiration estimation results include potential evapotranspiration raster data, nonlinear wind function value raster data, quality assurance indicators, and model uncertainty indices.

10. A potential evapotranspiration estimation system, characterized in that, The system includes: The preprocessing module is used to acquire meteorological data and surface data from multiple stations and to preprocess the meteorological data and surface data. The acquisition module is used to obtain potential state samples of the multiple sites during the water-free period; The inversion module is used to perform energy closure correction on the latent heat flux in the potential state sample, obtain potential evapotranspiration observation data based on the corrected latent heat flux, substitute the potential evapotranspiration observation data into the potential evapotranspiration model, and invert to obtain the equivalent wind function label. The training module is used to construct a sample set using the equivalent wind function label, meteorological data, and surface data, and to train the constructed machine learning model using the sample set to obtain the trained estimation model. The processing output module is used to obtain meteorological data and surface data of the target station, and output the corresponding nonlinear wind function based on the meteorological data and surface data of the target station and the estimation model. The estimation module is used to replace the wind function in the aerodynamic term of the potential evapotranspiration model with the nonlinear wind function, and output the potential evapotranspiration estimation result based on the replaced potential evapotranspiration model.

Images