Evapotranspiration estimation method based on data fusion and physical guided neural networks

Abstract

This application relates to the field of remote sensing hydrological monitoring technology and discloses a method for estimating evapotranspiration based on data fusion and a physics-guided neural network. This method acquires multi-source satellite remote sensing and meteorological data, utilizes an enhanced spatiotemporal adaptive reflectance fusion model and a data mining sharing algorithm to fuse coarse and fine resolution images, generating daily seamless high-resolution surface reflectance and surface temperature data. A physics-guided neural network model is constructed, comprising a parametric subnetwork and a main network embedding the surface energy balance physical equations. The parametric subnetwork is used to predict the thermodynamic roughness length, driving the main network to iteratively calculate sensible and latent heat fluxes in conjunction with meteorological data. The model is trained based on a loss function containing physical constraints. This invention, through deep coupling of data-driven and physical mechanisms, resolves the spatiotemporal resolution contradiction of remote sensing data and improves the model's generalization ability and physical consistency in data-scarce regions.

Description

Technical Field

[0001] This invention relates to the field of remote sensing hydrological monitoring technology, specifically to an evapotranspiration estimation method based on data fusion and a physically guided neural network. Background Technology

[0002] Evapotranspiration, as a major pathway for farmland water consumption, necessitates accurate estimation of its spatiotemporal variability for assessing irrigation efficiency and water resource management. While satellite remote sensing-based estimation methods can provide regional-scale spatial results, overcoming the limited representativeness of site observations, existing technologies still have significant limitations in model construction and data sources. Regarding model construction, while physics-based models (such as the single-source energy balance model SEBS and the Penman-Monteith model) are theoretically sound, their practical application is often limited by the difficulty in parameterizing key parameters such as surface impedance and thermodynamic roughness length. The uncertainty of these parameters directly impacts the energy distribution process, leading to estimation errors. On the other hand, while purely data-driven machine learning methods can establish strong nonlinear mappings, the lack of explicit physical constraints makes them prone to overfitting or learning spurious driving relationships outside the training data coverage, resulting in insufficient model generalization ability and interpretability.

[0003] To address these issues, existing research attempts to employ a physics- and data-driven coupling strategy, such as Physics-Guided Neural Networks (PINNs). However, current PINN-like methods are mostly applied to relatively simplified physical models and primarily rely on vegetation indices provided by optical sensors as key inputs. Vegetation indices reflect changes in canopy structure and morphology, and their response to water stress typically lags behind changes in surface temperature. Surface temperature can more directly and rapidly reflect the energy and water exchange between the surface and the atmosphere, but existing coupled models do not fully utilize the advantage of surface temperature in capturing early water changes, often marginalizing its constraining role in the energy balance equation, thus limiting the model's estimation accuracy under dynamic environments.

[0004] Furthermore, achieving refined management at the farmland plot scale requires the support of high spatiotemporal resolution remote sensing data. Current satellite sensor designs struggle to simultaneously meet the demands for both high spatial resolution and high temporal revisit frequency: for example, while sensors like MODIS or VIIRS possess daily observation capabilities, their spatial resolution is relatively coarse (kilometer-level), failing to identify plot details; while high-resolution sensors like Landsat have long revisit periods (16 days), making it difficult to capture rapid evapotranspiration fluctuations caused by irrigation or rainfall. This inherent contradiction in spatiotemporal resolution makes it difficult for existing evapotranspiration estimations to generate seamless, plot-level detailed, spatiotemporally continuous data, failing to meet the practical needs of refined agricultural water resource management.

[0005] Therefore, this invention proposes an evapotranspiration estimation method based on data fusion and physically guided neural networks to address the shortcomings of existing technologies. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides an evapotranspiration estimation method based on data fusion and a physically guided neural network. This method solves the problems of existing technologies being unable to obtain seamless, high-resolution daily remote sensing data and the difficulty in balancing the consistency of physical mechanisms and the generalization ability of the model in evapotranspiration estimation.

[0007] To achieve the above objectives, the present invention provides the following technical solution: an evapotranspiration estimation method based on data fusion and physically guided neural networks, comprising the following steps: Acquire multi-source satellite remote sensing data, reanalyze meteorological forcing data, ground-aided data and flux observation data, and perform spatiotemporal matching preprocessing on the data; The multi-source satellite remote sensing data is fused using an enhanced spatiotemporal adaptive reflectance fusion model to generate daily seamless high-resolution surface reflectance data, and daily seamless high-resolution surface temperature data is generated based on the high-resolution surface reflectance data using a data mining sharing algorithm. A physical-guided neural network model is constructed, which includes a cascaded parametric subnetwork and a physical computation layer. The physical computation layer contains a set of physical equations for the Earth's surface energy balance system. The vegetation index, high-resolution surface temperature data and reanalysis meteorological forcing data corresponding to the high-resolution surface reflectance data are input into the parameter subnetwork, and the thermodynamic roughness length is predicted by the neural network. The physical calculation layer is driven by the predicted thermodynamic roughness length, and combined with the reanalysis meteorological forcing data and ground-aided data, the sensible heat flux and latent heat flux are output through iterative calculation. A physical constraint loss function is constructed, which includes a data error term and a physical constraint term. The physical guidance neural network model is trained using the flux observation data, wherein the physical constraint term is used to constrain the predicted value of thermodynamic roughness length.

[0008] Preferably, the multi-source satellite remote sensing data includes coarse-resolution high-frequency data and fine-resolution low-frequency data; the step of acquiring multi-source satellite remote sensing data adopts the following data configuration strategy: for surface reflectance data, daily revisited coarse-resolution surface reflectance products are selected as coarse-resolution high-frequency data, and periodically revisited fine-resolution surface reflectance products are selected as fine-resolution low-frequency data; for surface temperature data, daily revisited coarse-resolution surface temperature products are selected as coarse-resolution high-frequency data, and periodically revisited fine-resolution surface temperature products are selected as fine-resolution low-frequency data.

[0009] Preferably, the steps for generating seamless daily high-resolution land surface temperature data specifically include: aggregating the high-resolution land surface reflectance data to the same spatial scale as the coarse-resolution land surface temperature data; constructing a nonlinear regression model between the coarse-resolution land surface temperature data and the aggregated high-resolution land surface reflectance data; applying the nonlinear regression model to the original high-resolution land surface reflectance data to obtain a preliminary prediction value of the high-resolution land surface temperature; calculating the residual between the coarse-resolution land surface temperature data and the aggregated preliminary prediction value of the high-resolution land surface temperature, and interpolating and superimposing the residual onto the preliminary prediction value of the high-resolution land surface temperature to obtain the final high-resolution land surface temperature data.

[0010] Preferably, the parameter subnetwork adopts a fully connected neural network architecture; the input features of the parameter subnetwork include: the normalized vegetation index, enhanced vegetation index and leaf area index calculated using the high-resolution surface reflectance data, the high-resolution surface temperature data, and the air temperature, wind speed and saturated vapor pressure difference in the reanalysis meteorological forcing data; the output layer of the parameter subnetwork is configured with an activation function to output a non-negative predicted value of the thermodynamic roughness length.

[0011] Preferably, the physical computation layer solves the embedded physical equations using an iterative algorithm based on the Moning-Obukhov similarity theory, the iterative algorithm including the calculation and updating of the atmospheric stability correction function.

[0012] Preferably, the method for constructing the physical constraint term in the physical constraint loss function specifically includes: calculating the prior reference value of the thermodynamic roughness length using a semi-empirical formula for the surface energy balance system; calculating the difference between the predicted value of the thermodynamic roughness length output by the parameter sub-network and the prior reference value; and incorporating the difference as a regularization term to penalize the output of the parameter sub-network into the physical constraint loss function.

[0013] Preferably, the method further includes a daily-scale evapotranspiration extension step based on the evaporation ratio method: calculating the instantaneous evaporation ratio at the time of satellite passage based on the sensible heat flux and the latent heat flux output by the physical computing layer; calculating the daily net radiation based on the daily-scale reanalysis meteorological forcing data; and, based on the constant evaporation ratio assumption, using the instantaneous evaporation ratio, daily net radiation, and the latent heat of vaporization of water, extending the instantaneous latent heat flux at the time of satellite passage to the total daily evapotranspiration.

[0014] Preferably, the reanalysis meteorological forcing data includes saturated vapor pressure difference, and the calculation steps of the saturated vapor pressure difference include: calculating the saturated vapor pressure using the temperature data in the reanalysis meteorological forcing data based on the Magnus empirical formula; and calculating the saturated vapor pressure difference based on the saturated vapor pressure and the relative humidity data in the reanalysis meteorological forcing data.

[0015] Preferably, the preprocessing steps for the flux observation data include: outlier removal and density correction of the original eddy covariance observation data; forced energy closure correction of the observed original sensible heat flux and original latent heat flux using the Bowen ratio method, and obtaining the true values ​​of the closed sensible heat flux and latent heat flux as true value labels for model training.

[0016] Preferably, the method further includes a model-driven factor analysis step: establishing a background dataset using the Shapley additive interpretation method, calculating the Shapley value of each input feature of the parameter sub-network on the evapotranspiration estimation result; quantifying the marginal contribution of each input feature based on the Shapley value, and verifying whether the physical guided neural network model conforms to the physical mechanism of evapotranspiration driven by meteorological factors and vegetation factors.

[0017] This invention provides a method for estimating evapotranspiration based on data fusion and a physically guided neural network. It has the following beneficial effects: 1. This invention integrates an enhanced spatiotemporal adaptive reflectivity fusion model with a data mining sharing algorithm to fuse coarse-resolution images with high revisit cycles and fine-resolution images with high spatial texture detail. This data processing strategy effectively overcomes the physical limitations of single sensors in terms of temporal frequency and spatial clarity, generating seamless daily 20-meter resolution surface reflectivity and surface temperature data. This provides fundamental data support with high spatiotemporal continuity for field-scale farmland water resource management and crop water consumption monitoring.

[0018] 2. This invention constructs a physics-guided neural network architecture that couples a parametric subnetwork with a physical computation layer. It utilizes the neural network to fit the thermodynamic roughness length, which is difficult to parameterize accurately, and then performs iterative calculations through an embedded set of surface energy balance equations. This design leverages the advantages of deep learning in nonlinear feature extraction while forcing the model output to adhere to the law of energy conservation through physical equations and a loss function containing physical constraints. Compared to purely data-driven models, this method significantly reduces dependence on training data, avoids overfitting, and exhibits stronger adaptability in regions lacking observational data.

[0019] 3. To address the low transparency of the internal computational process in deep neural networks, this invention introduces the Shapley additive interpretation method for posterior analysis of the model. By quantifying the marginal contributions of input features such as net radiation, saturated vapor pressure difference, and vegetation index to the evapotranspiration estimation results, it is possible to intuitively determine whether the model has learned the true physical driving relationship between energy supply and atmospheric moisture deficit on evapotranspiration. This verification mechanism eliminates the possibility that the model merely fits spurious correlations in the data, improves the credibility of the estimation results, and provides a clear physical direction for further model optimization. Attached Figure Description

[0020] Figure 1 This is a schematic flowchart of the method of the present invention; Figure 2 This is a structural block diagram of the system of the present invention.

[0021] The module includes: 10. Data preprocessing module; 20. Parameter fusion and generation module; 30. Evapotranspiration estimation module; and 40. Model evaluation module. Detailed Implementation

[0022] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0023] See attached document Figure 1 , Figure 1 This is a flowchart illustrating an evapotranspiration estimation method based on data fusion and a physically guided neural network according to an embodiment of the present invention. The present invention provides an evapotranspiration estimation method based on data fusion and a physically guided neural network, comprising the following steps: S100 acquires multi-source satellite remote sensing data, reanalyzes surface temperature data, meteorological forcing data, auxiliary data, and flux observation data, and performs spatiotemporal matching preprocessing. S200 uses an enhanced spatiotemporal adaptive reflectance fusion model to fuse multi-source satellite remote sensing data, generating daily seamless 20-meter resolution surface reflectance data, and combines data mining and sharing algorithms to generate daily seamless 20-meter resolution surface temperature data. S300 constructs a physical guided neural network model that couples the physical mechanism of the surface energy balance system. It cascades the parameter subnetwork with the physical calculation layer of the surface energy balance and uses the thermodynamic roughness length prediction value output by the parameter subnetwork to drive the physical calculation layer to calculate the sensible heat flux and latent heat flux. S400 trains a physics-guided neural network model based on a loss function that includes physical constraints, and validates the model and performs driving factor analysis using flux observation data and the Shapley additive interpretation method.

[0024] See attached document Figure 2 , Figure 2 This is a block diagram of an evapotranspiration estimation system based on data fusion and a physically guided neural network according to an embodiment of the present invention. The system is used to perform the above-described method and includes: The data preprocessing module 10 is used to perform step S100. The input data acquired by the data preprocessing module 10 includes coarse-resolution high-frequency imagery, fine-resolution low-frequency imagery, reanalysis surface temperature data, meteorological forcing data, and digital elevation model data covering the target area; the meteorological forcing data includes air temperature, solar radiation, saturated vapor pressure difference, wind speed, and air pressure; the data preprocessing module 10 performs cloud removal, resampling, projection transformation, and cropping on the above data to ensure the consistency of all spatial data in the geographic coordinate system, and processes the flux observation data into daily-scale evapotranspiration true values.

[0025] The parameter fusion generation module 20 is used to execute step S200. This module utilizes spatiotemporal linear interpolation and an enhanced spatiotemporal adaptive reflectance fusion model to fuse coarse-resolution high-frequency images and fine-resolution low-frequency images, outputting daily seamless 1km resolution land surface temperature data and 20m resolution land surface reflectance data. Based on this, the module uses a data mining sharing algorithm to establish a regression model between low-resolution land surface temperature and aggregated shortwave band reflectance, and applies this model to the 20m resolution land surface reflectance data, downscaling the land surface temperature to 20m resolution. Furthermore, the module calculates the normalized vegetation index and enhanced vegetation index based on the 20m resolution land surface reflectance data.

[0026] The evapotranspiration estimation module 30 is used to execute step S300. The physically guided neural network model constructed by the evapotranspiration estimation module 30 includes a parameter subnetwork, a main network, and a physical constraint loss function. The parameter subnetwork receives meteorological forcing data, vegetation index, and surface temperature data at a resolution of 20 meters, and outputs an initial estimate of the thermodynamic roughness length through neural network layers. The main network embeds the physical equations of the surface energy balance system, receives the initial estimate of the thermodynamic roughness length, and iteratively calculates the sensible heat flux and latent heat flux in combination with meteorological and remote sensing parameters. The evapotranspiration estimation module 30 uses the physical constraint loss function to guide the update of network parameters. This loss function includes a data error term and a physical constraint term. The physical constraint term applies a penalty to the output of the parameter subnetwork based on the thermodynamic roughness length calculated by the prior formula of the surface energy balance system.

[0027] The model evaluation module 40 is used to execute step S400. This module calculates the coefficient of determination, root mean square error, and mean deviation index to quantify the consistency between the evapotranspiration estimation results and flux observation data. Simultaneously, the module employs the Shapley additive interpretation method to calculate the contribution values ​​of the input features, identify key driving factors affecting the 20-meter resolution evapotranspiration estimation, and verify the model's learning effect on the physical mechanisms.

[0028] See attached document Figure 1 In the evapotranspiration estimation method based on data fusion and physical guided neural network described in this embodiment of the invention, the multi-source spatiotemporal data acquisition and preprocessing stage mainly includes data acquisition configuration and standardized cleaning processing.

[0029] In step S100, the data set acquired by the system specifically includes multi-source satellite remote sensing data, reanalysis meteorological forcing data, ground-aided data, and flux station observation data used for model training and verification.

[0030] To support the subsequent generation of high spatiotemporal resolution parameters for multi-source satellite remote sensing data, this embodiment adopts a data configuration strategy that combines "coarse-resolution high-frequency" and "fine-resolution low-frequency" data. The specific configuration scheme is as follows: Surface reflectance data configuration: The daily 1 km resolution surface reflectance product (such as VNP09GA) provided by the Visible Infrared Imaging Radiometer Suite (VIIRS) is selected as the coarse resolution high-frequency data source to capture the daily variation trend of surface reflectance; the 20 m resolution surface reflectance product (such as Level-2A product) provided by the Multispectral Imager (MSI) on the Sentinel-2 satellite is selected as the fine resolution low-frequency data source to provide fine spatial texture information.

[0031] Land surface temperature (LST) data configuration: Select daily 1 km resolution VIIRS land surface temperature products (such as VNP21A1D) as high-frequency data sources; and select 100 m resolution land surface temperature products provided by Landsat series satellites or 70 m resolution land surface temperature products provided by the Ecosystem Spaceborne Experimental Thermal Radiometer as low-frequency, high spatial resolution data sources.

[0032] Data from different sources are matched using band mapping relationships to construct data pairs suitable for subsequent spatiotemporal fusion algorithms.

[0033] For reanalysis of meteorological forcing data, this embodiment selects a high spatiotemporal resolution land surface data assimilation system (such as CLDAS-V2.0); the acquired variables include 2-meter air temperature. Surface incident solar radiation 10-meter wind speed Atmospheric pressure , specific moisture and relative humidity Considering the difference between the original resolution (e.g., 0.0625°) and the target resolution (20 meters) of the reanalysis data, this step includes spatial downscaling or resampling. Furthermore, based on the acquired temperature and relative humidity data, the saturated vapor pressure difference is calculated. As a key meteorological factor driving crop evapotranspiration; saturated vapor pressure difference The calculation follows the physical formula: ; in, Temperature (unit: degrees Celsius) Relative humidity (unit: percentage).

[0034] For ground-based auxiliary data, digital elevation model (DEM) data (e.g., NASA-DEM, 30-meter resolution) and land use cover data (e.g., GLC_FCS30D) covering the target area were acquired. Based on the land use type data, the canopy height of different vegetation types within the study area was determined using a lookup table method. and zero-plane displacement This provides static surface feature inputs for subsequent aerodynamic parameter calculations.

[0035] Based on flux station observation data, obtain the raw observation values ​​of the eddy covariance system, including sensible heat flux. Latent heat flux Net radiation and soil heat flux .

[0036] After obtaining the above raw data, perform the data cleaning and spatiotemporal alignment preprocessing sub-steps.

[0037] First, quality control is performed on multi-source satellite remote sensing data; the built-in quality assessment band (QA-band) is used to identify and remove pixels contaminated by clouds, aerosols or shadows to generate a high-quality effective observation mask; missing or invalid pixel values ​​are marked as invalid values ​​to avoid error propagation.

[0038] Secondly, a unified spatiotemporal benchmark correction is performed; all spatial data (including remote sensing imagery, reanalysis raster data, DEM, and land use data) are projected and transformed to a unified geographic coordinate system or projected coordinate system (such as WGS84 or UTM projection) to ensure strict spatial overlap. Simultaneously, all raster data is spatially cropped to preserve the coverage of the target study area.

[0039] Next, multi-scale resampling processing is performed; for data sources with inconsistent resolutions, an adaptive resampling strategy is adopted; for categorical data (such as land use types), nearest neighbor interpolation is used to resample to a standard resolution of 20 meters to maintain the category attributes; for continuous variables (such as meteorological data and DEM), bilinear interpolation or cubic convolution interpolation is used to resample to a resolution of 20 meters to ensure the spatial smoothness of the data; for CLDAS land surface temperature data, it is resampled to a resolution of approximately 7 kilometers to facilitate subsequent fusion calculations with the 1-kilometer resolution VIIRS-LST product.

[0040] Finally, energy closure correction and time-domain aggregation were performed on the flux observation data; outlier removal, coordinate rotation, and WPL density correction were performed on the original half-hour scale flux data, and the Bowen ratio method was used to force energy closure, solving the common energy non-closure problem in eddy covariance observations; the closed sensible heat flux... and latent heat flux The calculation is as follows: , ; The processed instantaneous flux data were further aggregated into diurnal evapotranspiration values, which served as ground truth labels for model training and validation. Through these steps, a spatiotemporal matching dataset containing meteorological, remote sensing, topographic, and ground truth labels was constructed, laying the data foundation for subsequent parameter generation and model training.

[0041] See attached document Figure 1 The daily seamless high-resolution land surface parameter generation stage in this embodiment of the invention specifically includes spatiotemporal fusion of land surface reflectance, downscaling of land surface temperature, and calculation of vegetation index. This stage aims to address the incompatibility of temporal and spatial resolution between a single remote sensing data source, providing standardized 20-meter resolution input features for subsequent physics-guided neural network models.

[0042] For daily coarse-resolution imagery with a resolution of 1 km (e.g., VIIRS / MODIS) and periodic fine-resolution imagery with a resolution of 20 m (e.g., Sentinel-2), the system executes step S210 to generate daily seamless 20 m surface reflectance data based on the Enhanced Spatiotemporal Adaptive Reflectance Fusion Model (ESTARFM). This step assumes that the change in surface reflectance over time is linear and uses a pair of reference images from two different time points to capture this trend. Specifically, at the prediction time... The two most recent pairs of cloudless coarse-resolution images acquired before and after this period. With fine resolution images As the baseline data, it is denoted as time. and For the predicted time To determine the target fine-resolution pixel value, the system searches for spectrally similar pixels within a moving window and uses a weighting function to comprehensively calculate the predicted value. The calculation process is shown in the following formula: ; in, Indicates the position (x, y) and time. band The fine-resolution reflectance prediction value; Reference time ( Pick or ); This represents the number of similar pixels within the sliding window. For the first The normalized weights of similar pixels are determined by spectral distance, temporal distance and spatial distance. The conversion coefficients are used to correct systematic errors between coarse and fine resolution sensors. For the specific implementation details of weight calculation and similar pixel selection in the ESTARFM algorithm, those skilled in the art can refer to relevant existing technical literature, as these are well-known techniques in the field and will not be elaborated upon here. Through this step, the system fills the temporal gap in the Sentinel-2 data, outputting continuous 20-meter surface reflectance data daily.

[0043] After acquiring high-resolution reflectance data, the system executes step S220, which performs surface temperature downscaling based on the Data Mining Sharing (DMS) algorithm, transforming the 1-kilometer resolution surface temperature (LST) into a 20-meter resolution. This step is based on the physical premise that there is a strong correlation between surface temperature, surface shortwave reflectance, and vegetation indices. First, the 20-meter surface reflectance data generated in step S210 is aggregated to the same spatial scale as the coarse-resolution LST (1 kilometer). Second, at the coarse-resolution scale, a nonlinear regression model is constructed between surface temperature and the aggregated shortwave reflectance. In this embodiment, a random forest regression algorithm is used to fit the nonlinear relationship, which can effectively handle the nonlinear characteristics under complex land cover. After the model is built, it is applied to the original 20-meter resolution land reflectance data to obtain a preliminary prediction of the land surface temperature at a resolution of 20 meters. To ensure energy conservation and eliminate system bias, the system further calculates the residual term, i.e., the difference between the original coarse resolution LST and the aggregated predicted LST, and then interpolates and superimposes the residual term onto the preliminary prediction value to obtain the final daily seamless 20-meter resolution land surface temperature. The mathematical expression of this process is as follows: ; in, This represents a trained random forest regression model; These are multi-band reflectance feature vectors with a resolution of 20 meters. This is the residual correction term after bilinear interpolation.

[0044] Subsequently, the system executes step S230, calculating key vegetation indices as vegetation parameter inputs to the physical-guided neural network. Based on the 20-meter resolution surface reflectance data generated in step S210, the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) are calculated. These indices reflect the growth status, cover, and leaf area density of vegetation, and are key biophysical parameters determining plant transpiration rates. The formula for calculating NDVI is as follows: ; The formula for calculating EVI is as follows: ; in, , , These represent the surface reflectance values ​​for the near-infrared band, red band, and blue band, respectively. This is the gain factor (with a value of 2.5). and These are the aerosol impedance coefficients, used to correct the influence of atmospheric aerosols on the red light band. Take 6. Take 7.5); The canopy background adjustment factor is set to 1. Through steps S210 to S230 above, the system completes the transformation from multi-source raw data to high spatiotemporal resolution model input features, and constructs a complete feature space that includes thermal conditions and vegetation conditions.

[0045] See attached document Figure 2 The Physically Guided Neural Network (PINN) construction phase coupled with the SEBS mechanism described in this embodiment of the invention aims to deeply integrate the physical mechanism of the Earth Surface Energy Balance System (SEBS) with the data mining capabilities of deep learning. The model structure constructed in this phase includes a parametric subnetwork, a physical computation layer (i.e., the main network), and a loss function based on physical constraints.

[0046] In step S310, a parametric subnetwork is constructed to generate a predicted value for the thermodynamic roughness length.

[0047] The parametric subnetwork is configured to address the difficulty in accurately parameterizing key parameters in traditional SEBS models, particularly the thermodynamic roughness length. This parameter typically relies on empirical methods in traditional approaches. The model has poor universality under different vegetation covers. The parameter subnetwork designed in this embodiment adopts a fully connected neural network architecture, including an input layer, four hidden layers and an output layer.

[0048] Specifically, the input features of the parametric subnetwork Includes the 20-meter resolution feature generated in step S200: surface temperature. Surface reflectance and calculated leaf area index Normalized Difference Vegetation Index Enhanced vegetation index and meteorological forcing data (temperature) Wind speed saturated water vapor pressure difference ) and canopy height Each hidden layer has 64 neurons and uses ReLU (Rectified Linear Unit) as the activation function to introduce a non-linear feature transformation. To prevent overfitting, a Dropout strategy is introduced between hidden layers with a random deactivation ratio set to 0.2. The output layer outputs a non-negative thermodynamic roughness length prediction value through an activation function (such as Softplus or Sigmoid scaling). This is used for subsequent physical layer calculations.

[0049] In step S320, a physical computation layer (main network) with embedded SEBS physical equations is constructed to realize the physical derivation from parameters to flux.

[0050] The physics computation layer does not contain the trainable weights of a traditional neural network. Instead, it hard-codes the physical equations of the Earth's surface energy balance system, constructing them as differentiable computational graph nodes. This layer receives the outputs of the parametric subnetwork. In addition to relevant meteorological and surface auxiliary data, the sensible heat flux is output through iterative calculation. and latent heat flux .

[0051] The core logic of the physical calculation layer is based on the Monin-Obukhov similarity theory (MOST). First, according to the surface energy balance equation, the net surface radiation... With soil heat flux The difference is determined by the sensible heat flux. and latent heat flux Consumption; ; Among them, net radiation and soil heat flux Sensible heat flux is calculated directly based on input radiation data and surface parameters. The calculation is encoded as the following aerodynamic impedance formula: ; In the formula, air density, The specific heat capacity of air at constant pressure. For surface temperature, For reference altitude temperature, Aerodynamic impedance for heat transfer.

[0052] In order to solve The physical computation layer implements an iterative calculation process based on a stability correction function. Aerodynamic impedance The calculation formula is as follows: ; ; In the formula, is the von Kármán constant (usually taken as 0.4). For frictional wind speed, For observation altitude, Zero planar displacement The momentum roughness length is based on the canopy height. The calculations are performed using empirical formulas. ; and These are stability correction functions for heat and momentum, respectively. It is the length of the Moninobhof.

[0053] because It is itself and Functions: ; In the formula, The acceleration due to gravity is used; the physics computation layer resolves the coupling relationship between the above formulas through a pre-defined finite number of iterations or by solving the fixed-point equations, ultimately obtaining the convergent induced heat flux. And then calculate the latent heat flux through the residuals. The physical computing layer uses an iterative method to solve the problem. , and The coupling relationship between them. The specific calculation process is as follows: Initialization assumptions: First, it is assumed that the atmosphere is in a neutral and stable state. Under this state, the influence of atmospheric stability is ignored, and the momentum transport stability correction function is applied. and heat transfer stability correction function The initial values ​​are all set to zero.

[0054] Initiating iterative estimation: Based on the initial zero-stability correction value, the initial frictional wind speed is calculated using the input wind speed, air temperature, and surface parameters. and aerodynamic impedance Then, based on the difference between the surface air temperature and the air temperature at its reference altitude, the initial sensible heat flux is calculated. .

[0055] Loop update calculation: Enter the iterative loop and repeat the following sub-steps: Update the Moning-Obukhov length: Use the frictional wind speed obtained in the previous step (or initialization). and sensible heat flux Calculate the current Moning-Obukhov length. .

[0056] Calculate the stability correction: based on the calculated The value is used to determine the atmospheric stability state (stable or unstable), and the momentum stability correction function is calculated using the corresponding integral similarity function (such as the Paulson function or the Buntinger-Dyer formula). and thermal stability correction function .

[0057] Update frictional wind speed and impedance: [Update the settings] Substitute the wind speed profile equation and recalculate the friction wind speed. The updated version Substitute the temperature profile equations and recalculate the aerodynamic impedance. .

[0058] Updated sensible heat flux: Utilizing the updated aerodynamic impedance Recalculate the sensible heat flux .

[0059] Convergence Criteria and Output: Calculate the difference between the sensible heat flux obtained in this iteration and the result of the previous iteration. If the absolute value of this difference is less than a preset convergence threshold (e.g., 0.01 W / m²), then the convergence is considered complete. 2 If the number of iterations reaches a preset limit, the computation is considered converged, and the iteration stops; at this time, the output is... This is the final determined sensible heat flux, and then the latent heat flux is calculated based on the residuals of the surface energy balance equation. This process enables the direct expression of physical mechanisms in the forward propagation of neural networks.

[0060] In step S330, a physical constraint loss function based on prior physical knowledge is constructed to train and constrain the model.

[0061] To ensure the output of the parameter subnetwork To ensure compliance with physical laws and physical consistency in the final estimated evapotranspiration, this embodiment constructs a composite loss function that includes a data error term and a physical constraint term. : ; in, This is a data error term used to measure the latent heat flux output by the physical computing layer. True values ​​of flux station observations Mean squared error (MSE) between: ; This is a physical constraint term used as an intermediate output of the constraint parameter subnet. Specifically, the prior reference value of the thermodynamic roughness length is calculated using classic semi-empirical formulas in the SEBS model (such as the Massman or Su model). While this prior value may not be entirely precise, it provides a range and trend of values ​​that align with physical expectations; the physical constraint terms are used to calculate predicted values. , and prior reference value Differences between them: ; This is a balancing coefficient used to adjust the weight of physical constraints in the total loss. It is achieved by minimizing... By using backpropagation algorithms (such as the Adam optimizer) to update the weights of the parametric subnetwork, the model can both fit the observed data and follow the physical laws of surface energy balance, thus exhibiting better generalization ability in areas lacking observation data.

[0062] See attached document Figure 1 The training and application phase of the evapotranspiration estimation model described in this embodiment mainly involves the construction of input feature vectors, the implementation of model training strategies, and the calculation of daily evapotranspiration based on a time-scale extension algorithm. This phase aims to transform the neural network model coupled with physical mechanisms into a practically applicable regional evapotranspiration estimation tool, and to extend the instantaneous flux at the time of satellite transit into daily-scale evapotranspiration data with hydrological application value.

[0063] In step S410, a high-dimensional input feature vector of the model is constructed and a training strategy is configured.

[0064] The system aligns and reassembles the preprocessed multi-source data at the pixel scale to construct the input feature matrix used to drive the physically guided neural network. For each spatial pixel with a resolution of 20 meters, the input feature vector is... It includes three types of key parameters: The first category consists of meteorological forcing parameters, including shortwave radiation incident on the surface. 2 meters temperature 10-meter wind speed saturated water vapor pressure difference and atmospheric pressure These parameters are primarily derived from downscaled reanalysis data and reflect the atmosphere's driving force on evapotranspiration.

[0065] The second category consists of remote sensing observation parameters, including surface temperature at a resolution of 20 meters. Normalized Difference Vegetation Index Enhanced vegetation index Leaf area index and surface albedo These parameters are calculated by fusing multi-source remote sensing data and reflect the condition of the underlying surface of the land.

[0066] The third category is static auxiliary parameters, including digital elevation model data. Canopy height Zero plane displacement and observation altitude .

[0067] To eliminate the impact of differences in physical dimensions and their numerical ranges on the convergence speed of neural network gradient descent, the feature matrix is ​​processed before being fed into the model. Standardize the sample so that its mean is 0 and its variance is 1.

[0068] In terms of training strategy configuration, the Adam optimizer is used to update model parameters. An initial learning rate (e.g., 0.001) is set, and a learning rate decay strategy is introduced. When the validation set loss no longer decreases within a preset number of epochs, the learning rate is automatically reduced to help the model converge to the global optimum. Observational data from flux stations are used to divide the training, validation, and test sets. Early stopping is used to prevent overfitting, i.e., training is terminated early when the validation set error continues to rise, thus preserving the optimal weight parameters.

[0069] In step S420, the instantaneous flux is calculated using the trained model.

[0070] Large-scale raster data is input into a pre-trained physical guidance neural network model. The model predicts the thermodynamic roughness length through a parametric subnetwork. The data is then processed through an embedded SEBS physical computing layer for forward propagation calculations, outputting the instantaneous perceived heat flux at the satellite's transit time (e.g., 11:00 AM). and instantaneous latent heat flux This process maintains physical constraints, ensuring that the output flux components satisfy the surface energy balance equations. .

[0071] In step S430, a diurnal evapotranspiration extension calculation based on the evaporation ratio method is performed.

[0072] Since satellites can only provide information on the surface conditions at the moment of transit, this embodiment employs the constant evaporation ratio (EF) assumption for time-scale extension in order to obtain information on water consumption throughout the day. This method assumes that during the daytime, the proportion of latent heat flux at the surface in effective energy (net radiation minus soil heat flux) remains relatively constant.

[0073] First, calculate the instantaneous evaporation ratio at the moment of satellite transit. : ; in, The instantaneous latent heat flux output by the model. For instantaneous net radiation, This represents the instantaneous soil heat flux. For Its calculation is based on the Earth's surface radiation balance equation: ; In the formula, It is the surface albedo. For instantaneous incident shortwave radiation, The surface emissivity. The specific emissivity of the atmosphere. This is the Stefan-Boltzmann constant. For instantaneous soil heat flux... An empirical formula based on net radiation and vegetation cover was used to estimate: ; In the formula, For vegetation coverage, and The ratios of soil heat flux to net radiation under full vegetation cover and bare soil cover, respectively (typically...). Take 0.05, Take 0.315).

[0074] Subsequently, daily net radiation was calculated using diurnal meteorological radiation data. Based on the assumption of a constant evaporation ratio, daily evapotranspiration (Unit: mm / day) The calculation formula is as follows: ; in, To make the unit Convert to The coefficient; The latent heat of vaporization of water (unit: J / kg) varies with temperature. ); The density of water (approximately 1000 kg / m³) 3 ); This is the daily soil heat flux, which is usually assumed to be close to zero in calculations with a daily step size.

[0075] Through steps S410 to S430, the system achieves a spatiotemporal extension from instantaneous point-scale physical deduction to regional daily-scale expansion, generating 20-meter resolution daily evapotranspiration products with clear physical meaning and spatiotemporal continuity. This provides directly applicable data support for refined agricultural water resource management. The specific settings for the evaporation ratio method and its related parameters can be adjusted by those skilled in the art according to the actual geographical environment; its basic principles are well-known in the field and will not be elaborated upon here.

[0076] See attached document Figure 1 The model accuracy verification and physical mechanism interpretability analysis stage described in this embodiment of the invention aims to quantitatively evaluate the accuracy of the generated evapotranspiration data and verify whether the physical guidance neural network has truly learned the physical mechanism of surface energy balance. This stage mainly includes an accuracy evaluation step based on independent flux stations and an interpretability analysis step based on SHAP values.

[0077] In step S510, a model accuracy evaluation based on flux stations is performed.

[0078] To objectively evaluate the generalization ability of the physics-guided neural network model, the system was validated using independent test set data that was not used in model training. The test set data came from eddy covariance flux observation stations distributed across different climate zones and underlying surface types (such as Xinxiang Station, Huailai Station, Gucheng Station, and Luancheng Station), covering the complete growth period of key crops such as winter wheat and summer maize. The system extracted 20-meter resolution model estimates for the coordinates of the flux stations. And compare it with the true ground observation values ​​after energy closure correction. Perform spatiotemporal matching.

[0079] Three statistical indicators were used to quantitatively describe the model performance: coefficient of determination (COP). ), root mean square error (RMSE), and average deviation (Bias).

[0080] Coefficient of determination Used to measure the goodness of fit between model estimation results and observed data, it characterizes the model's ability to explain the variance of the observed data. The calculation formula is as follows: ; The root mean square error (RMSE) measures the average error between model estimates and observed values, reflecting the model's prediction accuracy. The calculation formula is as follows: ; The average bias is used to assess whether the model has a systematic tendency to overestimate or underestimate, and the calculation formula is as follows: ; In the above formula, To verify the total number of samples; For the first The model estimates the evapotranspiration value for each sample. For the first Ground-based observed evapotranspiration values ​​for each sample; This represents the arithmetic mean of the ground observations. The specific calculation process for the above indicators can be implemented by those skilled in the art using conventional statistical software or programming libraries, and is considered well-known technology in this field.

[0081] In step S520, the physical mechanism verification of the driving factor based on the SHAP value is performed.

[0082] To address the "black box" characteristic of deep neural networks—the difficulty in explaining how input features specifically affect output results—this embodiment introduces the SHAP method for interpretability analysis. The SHAP method originates from the Shapley value concept in cooperative game theory and is applied to machine learning interpretation. Its aim is to fairly distribute the deviation of the model's predicted values ​​to each input feature, thereby quantifying the marginal contribution of each feature to a specific prediction result.

[0083] Specifically, firstly, 100 samples are randomly selected from the model training set as the background dataset. Based on this background dataset, the system calculates the SHAP value of each sample in the validation set and constructs a feature importance ranking graph and a dependency graph. Specifically, for the input feature set of the physically guided neural network... (Including net radiation) saturated water vapor pressure difference Leaf area index Surface temperature (etc.), the model's prediction output for a specific sample. This can be interpreted as a linear sum of the SHAP values ​​of all features: ; in, The predicted evapotranspiration value output by the model; This is the baseline output value of the model (usually the mean of all sample predictions); For the first The SHAP value of a feature indicates whether that feature contributes positively or negatively to the current prediction result. The total number of features.

[0084] By analyzing the SHAP value distribution of each feature, we can verify whether the model follows physical laws. The specific verification logic is as follows: If not only do the model accuracy metrics meet the requirements, but the SHAP analysis also shows net radiation and saturated water vapor pressure difference It has the highest feature importance and shows a positive correlation trend (i.e., and The larger the value, the more positive and increasing the corresponding SHAP value, which drives up the ET prediction value. This proves that the model has successfully captured the dominant driving force of energy supply and atmospheric moisture deficit on evapotranspiration.

[0085] If leaf area index and vegetation index (such as) The fact that the SHAP value of the model increases significantly during the peak growth period of the crop proves that the model correctly reflects the biophysical control mechanism of vegetation transpiration.

[0086] Conversely, if noise features unrelated to physical mechanisms exhibit abnormally high SHAP values, or if key physical parameters (such as...) are... If the SHAP value of a model contradicts common sense in physics (e.g., increased radiation leads to a significant decrease in ET prediction), it indicates that the model may have only fitted spurious correlations in the data and has not been verified by physical mechanisms.

[0087] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. An evapotranspiration estimation method based on data fusion and physically guided neural networks, characterized in that, Includes the following steps: Acquire multi-source satellite remote sensing data, reanalyze meteorological forcing data, ground-aided data and flux observation data, and perform spatiotemporal matching preprocessing on the data; The multi-source satellite remote sensing data is fused using an enhanced spatiotemporal adaptive reflectance fusion model to generate daily seamless high-resolution surface reflectance data, and daily seamless high-resolution surface temperature data is generated based on the high-resolution surface reflectance data using a data mining sharing algorithm. A physical-guided neural network model is constructed, which includes a cascaded parametric subnetwork and a physical computation layer. The physical computation layer contains a set of physical equations for the Earth's surface energy balance system. The vegetation index, high-resolution surface temperature data and reanalysis meteorological forcing data corresponding to the high-resolution surface reflectance data are input into the parameter subnetwork, and the thermodynamic roughness length is predicted by the neural network. The physical calculation layer is driven by the predicted thermodynamic roughness length, and combined with the reanalysis meteorological forcing data and ground-aided data, the sensible heat flux and latent heat flux are output through iterative calculation. A physical constraint loss function is constructed, which includes a data error term and a physical constraint term. The physical guidance neural network model is trained using the flux observation data, wherein the physical constraint term is used to constrain the predicted value of thermodynamic roughness length.

2. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The multi-source satellite remote sensing data includes coarse-resolution high-frequency data and fine-resolution low-frequency data; the steps for acquiring multi-source satellite remote sensing data employ the following data configuration strategy: For surface reflectance data, coarse-resolution surface reflectance products revisited daily are selected as coarse-resolution high-frequency data, and fine-resolution surface reflectance products revisited periodically are selected as fine-resolution low-frequency data. For land surface temperature data, coarse-resolution land surface temperature products revisited daily are selected as coarse-resolution high-frequency data, and fine-resolution land surface temperature products revisited periodically are selected as fine-resolution low-frequency data.

3. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The specific steps for generating seamless, high-resolution daily land surface temperature data include: The high-resolution surface reflectance data is aggregated to the same spatial scale as the coarse-resolution surface temperature data; Construct a nonlinear regression model between the coarse-resolution surface temperature data and the aggregated high-resolution surface reflectance data; The nonlinear regression model was applied to the original high-resolution surface reflectance data to obtain a preliminary prediction of high-resolution surface temperature. The residual between the coarse-resolution land surface temperature data and the aggregated high-resolution preliminary land surface temperature prediction is calculated, and the residual is interpolated and superimposed on the high-resolution preliminary land surface temperature prediction to obtain the final high-resolution land surface temperature data.

4. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The parameter subnetwork adopts a fully connected neural network architecture; The input features of the parameter subnetwork include: the normalized vegetation index, enhanced vegetation index and leaf area index calculated using the high-resolution surface reflectance data, the high-resolution surface temperature data, and the air temperature, wind speed and saturated vapor pressure difference in the reanalysis meteorological forcing data. The output layer of the parameter subnetwork is configured with an activation function to output a non-negative predicted value of the thermodynamic roughness length.

5. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The physical computation layer solves the embedded physical equations using an iterative algorithm based on the Moning-Obukhov similarity theory. The iterative algorithm includes the calculation and updating of the atmospheric stability correction function.

6. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The method for constructing the physical constraint terms in the physical constraint loss function specifically includes: The prior reference value of thermodynamic roughness length is calculated using a semi-empirical formula for the Earth's surface energy balance system. Calculate the difference between the predicted thermodynamic roughness length output by the parameter subnetwork and the prior reference value; The difference is used as a regularization term to penalize the output of the parameter subnetwork and is included in the physical constraint loss function.

7. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The method also includes a diurnal evapotranspiration extension step based on the evaporation ratio method: Based on the sensible heat flux and latent heat flux at the moment of satellite passage output by the physical computing layer, the instantaneous evaporation ratio at the moment of satellite passage is calculated. Daily net radiation is calculated based on the reanalysis meteorological forcing data at the daily scale; Based on the assumption of a constant evaporation ratio, the instantaneous evaporation ratio, net daily radiation, and latent heat of vaporization of water are used to extend the instantaneous latent heat flux at the moment of satellite transit to the total daily evaporation.

8. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, Further analysis of meteorological forcing data includes saturated vapor pressure differential, the calculation steps of which include: Based on the Magnus empirical formula, the saturated vapor pressure is calculated using the temperature data in the reanalysis meteorological forcing data. The saturated vapor pressure difference is calculated based on the saturated vapor pressure and the relative humidity data in the reanalysis meteorological forcing data.

9. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The preprocessing steps for the flux observation data include: Outlier removal and density correction were performed on the original vorticity-related observation data; The Bowenbi method was used to perform forced energy closure correction on the observed raw sensible heat flux and raw latent heat flux, and the closed true values ​​of sensible heat flux and latent heat flux were obtained as true value labels for model training.

10. The evapotranspiration estimation method based on data fusion and physically guided neural networks according to claim 1, characterized in that, The method also includes a model-driven factor analysis step: A background dataset is established using the Shapley additive interpretation method, and the Shapley value of each input feature of the parameter subnetwork on the evapotranspiration estimation result is calculated. The marginal contribution of each input feature is quantified based on the Shapley value to verify whether the physical guided neural network model conforms to the physical mechanism of evapotranspiration driven by meteorological and vegetation factors.

Images