Soil moisture dynamic inversion method and system based on mixed physical perception general differential equation

Abstract

This invention discloses a method and system for dynamic inversion of soil moisture based on a general differential equation using hybrid physical sensing, belonging to the field of artificial intelligence hydrological remote sensing technology. The invention first completes multi-source spatiotemporal data preprocessing and hydrological clustering, then generates physically constrained soil hydrological parameters through a static hypernetwork. Dynamic hydrological flux is modeled using a combination of spatiotemporal KAN and FiLM-Mamba networks. A general differential equation is constructed based on water balance to deduce soil moisture. A joint loss model with physical constraints is used to train the model, and finally, explicit hydrological analytical expressions are extracted through KAN symbolic regression to achieve prediction. This invention deeply integrates hydrological physical laws with neural networks, ensuring physical consistency of predictions, eliminating anti-physical anomalies, improving generalization ability in extreme scenarios, and possessing both high computational efficiency and strong interpretability.

Description

Technical Field

[0001] This invention belongs to the field of hydrological remote sensing inversion technology, specifically, it relates to a method and system for dynamic inversion of soil moisture based on a hybrid physical sensing general differential equation. Background Technology

[0002] Topsoil moisture is a core variable in the surface hydrological cycle and the exchange of energy between the land and the atmosphere. It directly affects the occurrence and evolution of hydrological processes such as precipitation infiltration, evapotranspiration, surface runoff, and deep seepage. It plays an irreplaceable role in agricultural drought monitoring, precision irrigation management of farmland, regional water resource allocation, and early warning of extreme climate events. With the development of hydrological remote sensing and land surface simulation technologies, the industry's demand for high spatiotemporal resolution and high precision regional soil moisture data continues to grow. Data-driven methods, represented by deep learning, are widely used in scenarios such as downscaling of multi-source remote sensing soil moisture data and time-series dynamic prediction due to their strong fitting ability to complex nonlinear relationships.

[0003] However, existing soil moisture inversion and prediction technologies still have certain shortcomings in practical engineering applications and scientific research:

[0004] Existing pure data-driven deep learning models based on architectures such as LSTM and Transformer essentially fit the statistical correlations between data points, without embedding hydrophysical laws into the model architecture, thus lacking strict physical constraints. In extreme climate scenarios not covered by the training data (such as torrential rainstorms or prolonged extreme droughts), the models are prone to outputting results that violate hydrophysical common sense, such as actual evapotranspiration exceeding the potential evapotranspiration limit, negative soil moisture values, or uncontrolled oversaturation. These anomalies may cause the model's extrapolation prediction accuracy to drop sharply in unseen scenarios, failing to meet the robustness requirements of actual business operations.

[0005] Traditional distributed hydrophysical models based on water balance equations, while strictly adhering to physical laws such as mass conservation, rely heavily on expert experience and massive amounts of in-situ measured data for manual calibration of key physical parameters such as saturated hydraulic conductivity, maximum storage capacity, and nonlinear porosity index. In large-scale study areas with significant spatial heterogeneity, the parameters vary greatly depending on different landforms and soil types. The time and data costs of parameter calibration are extremely high, and it is difficult to effectively extend from the station scale to the regional scale, which may not be suitable for the needs of large-scale, high-resolution soil moisture retrieval.

[0006] Long-term, high-resolution dynamic simulation of soil moisture requires models to simultaneously handle complex spatial-geographic interaction features and long-term time dependencies. Existing deep learning models are mostly based on multilayer perceptrons (MLP), LSTM, and Transformer as their core frameworks. Among them, MLP has limited ability to fit spatiotemporal nonlinear features and is prone to overfitting; LSTM is prone to gradient vanishing when processing long time series data and cannot effectively capture the hydrological evolution patterns across years; the computational complexity of the Transformer's self-attention mechanism increases quadratically with the sequence length, resulting in extremely high memory consumption and computation time; at the same time, existing Neural ODEs use the adjoint state method to solve long sequences, which may have serious memory leakage and computational efficiency bottlenecks, and cannot meet the needs of large-scale, long-term batch extrapolation; and the reasoning process of existing deep learning models is not transparent, and cannot provide hydrological researchers with intuitive physical evolution process analysis, nor can it extract universally applicable mathematical expressions of hydrological laws from the trained models. They can only achieve data fitting and prediction, which may not be able to support the research needs of scientific discovery such as exploring hydrological response mechanisms and identifying key control factors, thus limiting the application boundaries of the technology. Summary of the Invention

[0007] To address the aforementioned problems and technical deficiencies, this invention adopts the following technical solution: a soil moisture dynamic inversion method based on a hybrid physical sensing general differential equation. Step S1: Multi-source spatiotemporal data acquisition, preprocessing, and K-means clustering: acquire static geographic attribute data and dynamic meteorological time-series data of the study area. After spatiotemporal alignment and standardization, use the K-means clustering algorithm to divide the static geographic attribute data and obtain hydrological functional zones and their corresponding clustering zone labels.

[0008] Step S2: Generation of physical parameters based on static hypernet: Construct a static hypernet GeoHyperNet, taking the static geographic attribute data from step S1 and its corresponding clustering partition labels as input, and combining the physical prior bias mechanism and spatial nonlinear mapping to generate soil hydrological parameters with clear physical meaning. The soil hydrological parameters include at least saturated hydraulic conductivity Ksat, maximum water storage capacity Smax, and nonlinear porosity index c_exp.

[0009] Step S3: Dynamic flux modeling based on spatiotemporal network and feature modulation: The dynamic meteorological time series data is feature-encoded using a temporal Kolmogorov-Arnold network. The extracted time series features and the static features output by the static supernetwork are used as inputs to the dynamic flux network FiLMFluxNN. The dynamic flux parameters are output through the physical head PhysicsHead. The dynamic flux parameters include at least the actual evapotranspiration efficiency and drainage correction term.

[0010] Step S4: Dynamic estimation of soil moisture based on the general differential equation of hybrid physical sensing: Based on the principle of mass conservation, a general differential equation UDE based on the water balance equation is constructed, and the meteorological driving input, evapotranspiration consumption and drainage flux are used as the rate of change of the derivative of the equation; the soil hydrological parameters generated in step S2 and the dynamic flux parameters generated in step S3 are substituted into the general differential equation, and the solution is iteratively obtained along the time step by the explicit Euler integral algorithm to deduce the continuous predicted value sequence of soil moisture at the current time and the next time step.

[0011] Step S5: Model training based on physical constraint joint loss: Construct a joint loss function that integrates data-driven error, mass conservation constraint and flux boundary constraint, calculate the loss gradient using the soil moisture prediction value obtained in step S4, and perform backpropagation joint training optimization on the static supernetwork, dynamic flux network and general differential equation module.

[0012] Step S6: Explicit extraction and prediction loop of hydrological equations based on native symbolic regression: After the model training converges, the learned network weights are extracted from the Kolmogorov-Arnold network layer in the static supernetwork and spatiotemporal network. Through polynomial reconstruction and threshold pruning, an explicit analytical expression of hydrological evolution and the corresponding feature importance map are generated. The analytical expression is used to replace the neural network to realize the real-time prediction of soil moisture, forming a complete soil moisture prediction loop.

[0013] Preferably, the processing procedure for the static hypernetwork GeoHyperNet in step S2 includes:

[0014] The static geographic attribute data is divided into soil attribute branch, terrain feature branch, geographic coordinate branch, and clustering condition branch;

[0015] For the geographic coordinate branch, a multi-frequency sine and cosine basis function is used to encode the geographic coordinates. The encoding result is input into the spatial KAN encoder to learn the hidden spatial periodic features and nonlinear spatial interaction relationships.

[0016] For the clustering conditional branch, the clustering partition label is converted into a one-hot code, and a preset physical bias benchmark value is obtained through a matrix multiplication lookup table mechanism. A learnable bias modulation layer is then superimposed, and hydrophysical prior knowledge is injected.

[0017] For the topographic feature branch, a gated modulation signal is generated through the KAN layer, and this signal is used to adaptively weight and fuse the soil property branch features processed by the KAN layer.

[0018] The output features of the four branches are concatenated along the feature dimension and input into the deep KAN network to generate initial parameter logit values. The initial parameter logit values ​​are added to the physical prior benchmark values ​​and learnable bias terms output by the clustering conditional branches, and then mapped to absolute physical dimensions through a safety boundary constraint layer composed of Softplus and Sigmoid activation functions to obtain soil hydrological parameters with clear physical meaning.

[0019] Furthermore, the dynamic flux modeling in step S3 specifically includes:

[0020] Time feature coding: Extract time coding information from dynamic meteorological time series data, retain periodic sine and cosine variables as annual cycle benchmark coding, and fit the residual correction amount of seasonal asymmetric fluctuations through Chebyshev polynomial KAN layer;

[0021] Feature linear modulation: The encoded temporal features are input into the dynamic flux network FiLMFluxNN, which uses the state-space model Mamba as its backbone network and simultaneously receives scaling parameters generated by the projection of the static hypernetwork. Translation parameters The dynamic features are softened residual modulation through a characteristic linear modulation mechanism;

[0022] Flux parameter output: The modulated feature input physical head, PhysicsHead, employs a hierarchical fully connected network combined with an activation function to output the actual evapotranspiration efficiency constrained by the Sigmoid function. And the non-negative drainage correction term constrained by the Softplus function. .

[0023] Preferably, the calculation process of the general differential equation UDE in step S4 is as follows:

[0024] Construct the ordinary differential equation of water balance Where S is the soil moisture content, For precipitation replenishment, For actual evaporation, For drainage flux;

[0025] Among them, actual evaporation PET represents the potential evaporation of the input. The actual evapotranspiration efficiency generated by step S3;

[0026] Drainage flux ,in The current soil moisture level. For saturated hydraulic conductivity, For maximum water storage capacity, It is a nonlinear porosity index;

[0027] The ordinary differential equation is solved using the explicit Euler method. The soil moisture content at the next time step is: Apply rigid physical constraints to the updated soil moisture content. .

[0028] Furthermore, the joint loss function in step S5 is the sum of three parts: the rate-of-change weighted Huber loss, the stepwise mass conservation loss, and the flux boundary constraint loss, wherein:

[0029] Rate-of-change weighted Huber loss: Calculates the absolute change of the true observation over adjacent time intervals. Based on this, an adaptive rate of change weight is constructed. Multiply the weight by the Huber loss between the predicted and observed values;

[0030] Stepwise mass conservation loss: Calculating the mean square error of soil moisture change and net flux input. ;

[0031] Flux Boundary Constraint Loss: Minimizing the Mean Function The constraint is that the actual evapotranspiration is non-negative and not greater than the potential evapotranspiration.

[0032] Preferably, the explicit extraction of physical equations in step S6 specifically includes:

[0033] For a KAN layer that has converged during training, obtain its linear basis weight matrix. and Chebyshev polynomial weight tensor ;

[0034] Based on the function of the Kolmogorov-Arnold network, the mapping of each input feature is reconstructed into a univariate symbolic expression. ,in for Chebyshev polynomials of the first kind of order;

[0035] The symbolic expression is traversed and parsed, and only addition sub-items with absolute value coefficients greater than a preset pruning threshold are retained, while redundant nodes are removed.

[0036] Output explicit analytical expressions for hydrological evolution and feature importance maps.

[0037] A soil moisture dynamic inversion system based on a hybrid physical sensing general differential equation, used to execute the soil moisture dynamic inversion method according to any one of claims 1 to 6, characterized in that it includes:

[0038] The preloading and data management module is used to acquire static geographic data and dynamic time-series meteorological data, complete data preprocessing and preloading, and build a hierarchical cross-validation data loader.

[0039] The GeoHyperNet module, which incorporates a spatial KAN encoder and a hydrological clustering prior lookup table, is used to fuse static geographic features with hydrophysical prior biases and output soil hydrological parameters with clear physical meaning.

[0040] The dynamic flux neural state space module, which includes a time-based KAN encoder and a feature linear modulation network, is used to capture the long-range dependencies of meteorological time-series data and output physical process correction variables for evapotranspiration and drainage.

[0041] The explicit Euler general differential equation calculation module is used to solve the soil moisture change process based on the water balance differential equation and the explicit Euler integral method, and output the time series prediction value of soil moisture.

[0042] The multidimensional joint physical loss calculation module is used to construct a joint loss function that integrates adaptive weighted rate of change, mass conservation constraints, and flux boundary constraints, providing gradient descent directions for model training;

[0043] The symbolic regression and interpretability extraction module is used to parse the multinomial weights of the KAN layer, extract explicit mathematical expressions through symbolic regression, and output explicit formulas for hydrological evolution and feature importance maps.

[0044] Preferably, the GeoHyperNet module of the static physical constraint supernetwork adopts a four-branch parallel architecture, namely, soil attribute branch, terrain feature branch, geographic coordinate branch and clustering condition branch. After each branch completes the extraction and processing of the corresponding features, they are fused by the deep KAN network to output soil hydrological parameters that conform to the absolute physical dimension constraints.

[0045] Furthermore, the dynamic flux neural state space module uses the Mamba state space model as the backbone network and combines it with a feature linear modulation mechanism to achieve residual-type soft modulation of static geographic features on dynamic temporal features.

[0046] Furthermore, the symbolic regression and interpretability extraction module incorporates a SymPy algebraic computation engine to simplify the reconstructed symbolic expression, prune redundant terms, and ultimately output an explicit hydrological analytical expression and a feature importance ranking map in LaTeX format.

[0047] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0048] (1) This invention embeds the water balance equation deep into the core of the deep learning architecture and constructs a derivation module with the general differential equation as the core, ensuring that the evolution of soil moisture follows the principle of mass conservation from the model structure level; at the same time, it performs joint optimization through the joint loss function composed of stepwise mass conservation loss and flux boundary constraint loss, forcing the evapotranspiration and drainage flux output by the model to strictly conform to the common sense of hydrophysics; with the bidirectional physical hard constraint in the explicit Euler integral process, it effectively eliminates anti-physical phenomena such as negative soil moisture, uncontrolled oversaturation, and actual evapotranspiration exceeding the upper limit; so that the model can still maintain stable extrapolation prediction ability in extreme drought, torrential rain and other scenarios not covered by training data, significantly improving the generalization performance of the model and the robustness of actual business applications;

[0049] (2) By introducing K-means hydrological clustering partitioning and hydrological expert prior lookup table, this invention injects prior knowledge in the field of hydrology into the model initialization stage, which greatly reduces the model's dependence on massive in-situ measured calibration data. Combined with a multi-frequency sine and cosine basis spatial KAN encoder, it can effectively capture the nonlinear interaction characteristics of geographic space. Combined with the terrain gating modulation mechanism, it can accurately characterize the control effect of terrain on soil water retention capacity. It can adaptively map discrete static geographic attributes into a continuous spatial hydrological parameter field that conforms to the constraints of absolute physical dimensions, effectively alleviate the difficulty of parameter regionalization expansion caused by spatial heterogeneity, and significantly improve the applicability of the model in unmeasured areas and complex terrain.

[0050] (3) This invention uses Kolmogorov-Arnold network and Mamba state space model to completely replace the traditional LSTM and Transformer architecture. It not only overcomes the bottleneck of gradient vanishing and computational complexity explosion when the traditional architecture processes long time series data, but also achieves higher computational efficiency with extremely low memory consumption. Among them, the time KAN encoder adopts the design of linear periodic basis plus Chebyshev residual correction, which not only ensures the stable fitting of annual cycle hydrological fluctuations, but also accurately captures the nonlinear drift characteristics of seasonal changes in soil moisture in different regions. The Mamba backbone network and FiLM feature linear modulation mechanism are deeply coupled to efficiently capture the long-range time dependence of ultra-long sequences with linear time complexity, while realizing the progressive constraint control of static geographical environment on dynamic hydrological process, taking into account the fitting accuracy and computational efficiency of the model.

[0051] (4) Based on the characteristics of the KAN network, this invention has a built-in native symbolic regression engine. Without relying on complex external parsing libraries, it can directly reconstruct the symbolic expression of the feature map from the linear weights and Chebyshev polynomial weights of the KAN layer that has been trained and converged. After removing redundant terms through threshold pruning, it directly outputs the explicit algebraic analytical expression of hydrological evolution and the feature importance map. This allows hydrological researchers to directly extract interpretable hydrological response formulas from deep learning models and identify key control factors affecting soil moisture dynamics. It realizes a paradigm leap from data black-box fitting to the discovery of hydrological scientific laws. It can not only be used for operational soil moisture prediction, but also support scientific research such as hydrological mechanism exploration, which greatly expands the application boundaries of the technology. Attached Figure Description

[0052] In the attached diagram:

[0053] Figure 1 This is a flowchart of the present invention;

[0054] Figure 2 This is a system module diagram of the present invention. Detailed Implementation

[0055] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments. Generally, the components of the embodiments of the present invention described and shown in the accompanying drawings can be arranged and designed in various different configurations.

[0056] Example 1:

[0057] like Figure 1 As shown, the method for dynamic inversion of soil moisture based on the general differential equation of hybrid physical sensing has the following steps:

[0058] Step S1: Acquisition of multi-source spatiotemporal data, clustering, partitioning, and preloading processing:

[0059] S101. Data Acquisition and Composition: Acquire static geographic attribute data and dynamic meteorological time series data of the study area.

[0060] Static geographic attribute data includes: soil attributes (clay content, sand content, bulk density, organic carbon, porosity, a total of 5 dimensions), topographic factors (elevation, slope, a total of 2 dimensions), and geographic coordinates (longitude, latitude, a total of 2 dimensions).

[0061] The dynamic meteorological time series data includes precipitation, potential evapotranspiration, surface temperature, and leaf area index, generating lagged first-order variables for precipitation and potential evapotranspiration, which together with the original variables constitute a 6-dimensional dynamic meteorological driver.

[0062] S102, Hydrological Functional Zoning:

[0063] This step introduces the K-means clustering algorithm to perform hydrological functional zoning based on the static attributes of spatial grid points. For example, the number of clusters is set to K=4, and the study area is divided into transition zone, dry-hot / clay zone, and wet / sandy zone. A corresponding one-hot coding label (4 dimensions in total) for each spatial grid is generated for each cluster partition.

[0064] S103, Data Preprocessing and Efficient Preloading Mechanism:

[0065] All the above data undergoes spatiotemporal alignment and standardization. To overcome the I / O bottleneck of data loading in traditional deep learning frameworks, this invention introduces a pre-loading mechanism during the model training initialization phase: considering that dynamic meteorological data belongs to long-term high-dimensional tensors, all standardized tensor data is converted into bfloat16 or float32 mixed-precision format and loaded into the video memory of the computing device (such as GPU) in one complete operation. This mechanism completely eliminates the I / O transmission bottleneck caused by the repeated data transfer between the CPU and GPU by the multi-process DataLoader during training, realizing non-blocking batch slice reading and improving the training efficiency of long-term data.

[0066] Step S2: Generation of absolute dimensional physical parameters based on the static hypernet (GeoHyperNet):

[0067] Traditional spatial parameter calibration is extremely time-consuming. This invention constructs a static hypernetwork GeoHyperNet, which takes the 13-dimensional static features (5-dimensional soil + 2-dimensional topography + 2-dimensional latitude and longitude + 4-dimensional clustering one-hot code) output from step S1 as input. Through multi-branch parallel processing and physical constraint mapping, it generates soil hydrological parameters with absolute physical dimensions, as well as modulation parameters for subsequent dynamic modules.

[0068] S201, Multi-branch parallel feature extraction:

[0069] GeoHyperNet employs a four-branch parallel architecture to process different types of static input features separately:

[0070] Spatial KAN coding branch: Extracts latitude and longitude features, abandons traditional fixed-frequency location coding, and adopts a formula Multiply by an adaptive multi-frequency segment, then pass through The function is expanded to construct a spatial periodic basis containing low-frequency global trends and high-frequency global details. The expanded features are then concatenated with the original latitude and longitude and input into a two-layer Chebyshev multinomial Kolmogorov-Arnold network to learn highly nonlinear geospatial interactions.

[0071] Terrain-gated modulation branch: Elevation and slope features are input into an independent KAN module for feature transformation, followed by a Sigmoid activation function to generate a value range of... Gated modulation signal.

[0072] Soil property branch: The 5-dimensional soil properties are input into the KAN module to extract deep features, and then multiplied element-wise with the above-mentioned terrain-gated modulation signal. This step precisely describes, in terms of physical mechanisms, the effective weight of topography (such as slope) in determining the actual water retention capacity of the soil—a hydrological process.

[0073] Clustering Conditional Physical Prior Branch: A hydrological expert prior lookup table is constructed, mapping the 4D clustering one-heat encoding to preset physical bias benchmark values ​​through matrix multiplication. For example, for arid / clay regions, the lookup table directly assigns extremely low initial values ​​for saturated hydraulic conductivity and extremely high initial values ​​for water storage capacity, ensuring that the network initialization phase conforms to hydrological and physical common sense. Simultaneously, a linear layer initialized to zero is superimposed to provide learnable bias fine-tuning.

[0074] S202, Feature Fusion and Absolute Physical Dimension Mapping:

[0075] The output features of the four branches are concatenated and input into the deep fusion KAN network to generate a 3D initial parameter logit value. This logit value is then added to the physical bias output from the clustering prior branch in S201 and input into the absolute physical dimension mapping layer. Strict constraints are applied through an activation function with physical boundaries, ultimately outputting the following key hydrological parameters: saturated hydraulic conductivity (… ): Unit: mm / day;

[0076] Maximum water storage capacity ( ): The unit is mm, and the physical range is strictly limited to 10~90mm;

[0077] Brooks-Corey nonlinear exponent ( ): , dimensionless.

[0078] S203, Dynamic modulation parameter generation:

[0079] In addition to the physical parameters mentioned above, GeoHyperNet also generates scaling parameters through an additional linear header. Translation parameters These two parameters will be passed to the subsequent dynamic flux network to establish characteristic linear modulation control of the static geographical environment on the dynamic hydrological evolution process.

[0080] Step S3: Dynamic flux modeling based on spatiotemporal KAN and feature modulation (FiLM):

[0081] S301, Time Feature Coding:

[0082] Extract time variables (such as Julian Day DOY) and map them to periodic features using sine / cosine transforms. and Based on this, a time-based KAN encoder is used, preserving the basic... The encoding serves as a stable annual cycle basis for propagation, while the residual correction for seasonal fluctuations is learned using the Chebyshev KAN branch. The two are adaptively fused through a learnable gating mechanism, and the fusion formula is as follows:

[0083]

[0084] This design, in the early stages of training, uses basic... The encoding provides stable annual gradient guidance; in the later stages of training, the gradual activation of the KAN residual branch accurately captures nonlinear climate drift characteristics in different regions, such as "seasonal changes in soil moisture content caused by asymmetry in the start and end times of the dry season".

[0085] S302, Dynamic State Space Feature Extraction:

[0086] The encoded temporal features output from step S301 are concatenated with the dynamic meteorological features output from step S1 and then input into the dynamic flux backbone network (FiLMFluxNN). This network employs the Mamba state-space model (MambaBlock) with linear time complexity as the core temporal processor, replacing traditional LSTM or Transformer architectures. The Mamba model can efficiently maintain continuous state memory over time steps of hundreds of days, significantly improving the computational efficiency and modeling capability when handling long-term hydrological processes.

[0087] S303, FiLM residual softening modulation:

[0088] Dynamic temporal hidden features extracted from Mamba networks Introduce the scaling parameters generated by the static hypernet in step S2. Translation parameters It performs characteristic linear modulation to achieve implicit control of the dynamic process by the static environment. The modulation process uses a residual softening modulation formula:

[0089]

[0090]

[0091] in This is a learnable initial weak modulation intensity parameter used to control the degree of interference of static features on dynamic features. This structure deeply simulates the key hydrophysical mechanism of the gradual constraint and control of static soil and topographic environment on dynamic precipitation and evaporation processes, ensuring that dynamic processes remain physically plausible within a static environmental context.

[0092] S304, Physical process variable output:

[0093] The modulated features are input into the PhysicsHead network, which performs physical constraint mapping through a multi-layer fully connected network and a specific activation function, outputting the following two types of key process variables:

[0094] Actual evaporation efficiency Constrained by the Sigmoid activation function;

[0095] Drainage Neural Network Correction Term The range is limited to 0~50mm / day by the combined constraints of the Softplus activation function and the Clamp operation.

[0096] Step S4: Derivation of Euler integral based on the general differential equation:

[0097] This step follows the strict principle of stepwise water balance in terms of physical mechanism, constructs and solves a general differential equation system, and realizes the temporal evolution of soil moisture.

[0098] S401, Construction of Differential Equations:

[0099] Based on the law of conservation of mass, an ordinary differential equation for soil moisture change is constructed:

[0100] .

[0101] in, This represents the current surface soil moisture content. The effective precipitation supply is the input amount; The actual evapotranspiration consumption estimated by the model; This represents the combined drainage flux of deep drainage and surface runoff. The equation strictly requires that the rate of change of soil moisture over time must maintain a dynamic balance with the net flux (input minus output).

[0102] S402. Input Preprocessing and Flux Calculation:

[0103] First, the input lagged first-order precipitation and potential evapotranspiration are denormalized to restore them to their true physical dimensions. Then, a hard constraint is applied using the ReLU function to ensure that the meteorological driving term is non-negative. ).

[0104] Subsequently, based on the physical parameters generated in step S2 and the process variables output in step S3, each flux term is calculated:

[0105] Actual evaporation: ;

[0106] Drainage flux: ;

[0107] The drainage flux formula introduces a relative saturation S / S_max and allows for slight surface water accumulation (the upper limit of saturation is relaxed to 1.2) to avoid abnormal problems such as negative bases or exponential explosions in numerical calculations.

[0108] S403, Efficient Numerical Integration Solution:

[0109] To eliminate the computational overhead caused by the Global Interpreter Lock (GIL) in the Python environment, this invention preferably employs an explicit Euler integral algorithm optimized by static compilation based on PyTorch just-in-time compilation (@torch.jit.script), according to the soil moisture at the previous time step. Iterative solution along the time dimension:

[0110]

[0111] The integration process is executed entirely in parallel on the GPU, eliminating the need for frequent data exchanges between the CPU and GPU, which significantly improves the computational efficiency of long-term time-series simulations.

[0112] S404, Physical security boundary constraints:

[0113] To prevent the model from outputting outliers that violate physical laws under extreme climate gradients, a two-way safety hard constraint is forcibly applied after each integration step:

[0114] ,

[0115] This constraint ensures that the lower limit of soil moisture is... To avoid division by zero errors caused by negative values ​​in numerical calculations, the upper limit shall not exceed the maximum water storage capacity. times.

[0116] Step S5: Optimization of physical joint loss based on adaptive weighting:

[0117] To integrate the high accuracy of data-driven approaches with the conservation of physical equations, this invention constructs a three-in-one joint loss function for gradient updates of the model network:

[0118] S501. Composition of the joint loss function:

[0119] The joint loss function consists of three parts: rate-of-change weighted Huber loss, stepwise mass conservation loss, and flux boundary constraint loss. Its mathematical expression is as follows:

[0120] .

[0121] in, and These are the weight hyperparameters for the mass conservation loss and the flux boundary constraint loss, respectively, used to control the proportion of physical constraint penalties in the total loss.

[0122] S502, Rate of Change Weighted Huber Loss:

[0123] To address the issue that soil moisture sequences, where the majority of samples are from periods of relative stability, directly using the mean-gravitational osmosis (MSE) can easily lead to model degradation and difficulty in accurately capturing extreme events. This invention proposes an adaptive weighting mechanism:

[0124] First, calculate the absolute change in observations between adjacent time points:

[0125]

[0126] Construct an adaptive dynamic weight matrix based on this change:

[0127]

[0128] hyperparameters The weighted intensity coefficient is usually taken as... This weighting mechanism can dynamically amplify the loss gradient during periods of drastic changes in soil moisture (such as heavy rain infiltration or rapid drought), guiding the model to focus more on fitting key hydrological abrupt changes and effectively overcoming the characteristic inundation problem caused by long-term flat periods.

[0129] Multiplying the above weight matrix element-wise with the Huber loss between predicted and observed values ​​yields the rate-of-change weighted Huber loss:

[0130] .

[0131] This loss causes the model to automatically amplify the gradient penalty during "critical mutation frames" when rainfall infiltrates rapidly or drought worsens sharply.

[0132] S503, Stepwise mass conservation loss:

[0133] To ensure that the hydrological processes output by the model strictly follow the principle of water balance, this invention directly applies conservation constraints to the actual flux components calculated by the model's forward propagation:

[0134]

[0135] This constraint mandates that the change in soil moisture at each time step must be strictly equal to the sum of precipitation replenishment, evapotranspiration consumption, and drainage flux. This ensures that the evolution of the model over a long time span fully conforms to the water balance equation, completely eliminating the unbounded accumulation of errors and anti-physical water generation phenomena commonly found in pure black-box models.

[0136] S504, Flux Boundary Constraint Loss:

[0137] To impose physical boundary conditions and ensure the physical rationality of the output flux, this invention utilizes an asymmetric activation function to construct a flux boundary constraint loss:

[0138] .

[0139] Wherein, mean(ReLU(-E) act Penalizes cases where actual evapotranspiration is negative, forcing actual evapotranspiration to be non-negative; mean(ReLU(E) act -PET)) constrains actual evapotranspiration to not exceed the potential evapotranspiration limit, ensuring that the evapotranspiration process strictly adheres to the physical limitations of energy and moisture availability.

[0140] Step S6: Native KAN symbolic regression and explicit physics formula extraction:

[0141] The black-box nature of deep learning models is considered the biggest obstacle to their application in hydrological science. To overcome this technical bottleneck, this invention utilizes a native symbolic regression engine built into independent Chebyshev KAN layers to directly extract explicit physical formulas from a trained network; the steps are as follows:

[0142] S601, Weight Parameter Extraction:

[0143] After the model training converges, the linear basis weight matrix is ​​extracted from the target KAN layer to be analyzed. and Chebyshev polynomial weight matrix ;

[0144] S602, Symbolic Expression Reconstruction:

[0145] Using the Chebyshev polynomial expansion principle, the mapping relationship of each input feature term is reconstructed into a single-variable symbolic algebraic expression:

[0146] ;

[0147] in, The weight matrix is ​​a linear basis; For the corresponding number Weighting coefficients of the Chebyshev polynomial; It is the highest order; express A basis of first-order Chebyshev polynomials. Because Chebyshev polynomials strictly require the input domain to be within... Between, here is introduced The hyperbolic tangent function preprocesses the input variables, safely mapping unbounded latent features and converging them into the orthogonality domain.

[0148] S603, Redundancy pruning:

[0149] The symbolic expression reconstructed by the SymPy algebra engine is parsed and simplified term by term. Each additive term in the polynomial expression is traversed, and if the absolute value coefficient of that term is less than a preset pruning threshold (e.g., ...), it is pruned. When the condition is met, the sub-item is identified as redundant noise and removed.

[0150] S604, Analytical Formula Output and Feature Importance Analysis:

[0151] After redundancy pruning, a high-order analytical mathematical formula is generated and output in LaTeX format. This formula can serve as an empirical equation for explicit soil moisture content prediction in specific hydrological zones, providing an interpretable mathematical expression for hydrological research.

[0152] Furthermore, based on the absolute weights of the corresponding sub-items of each input feature, a feature importance ranking map is automatically generated. This map visually demonstrates the contribution of different geographical environmental factors (such as soil properties and topographic factors) to the hydrological process output, providing a quantitative basis for revealing hydrological response patterns and identifying key control factors.

[0153] S605, Final Soil Moisture Content Output:

[0154] After establishing the complete prediction loop described above, real-time or future static geographic attributes and dynamic meteorological time-series data are input into the system. The model generates the corresponding actual evapotranspiration and drainage fluxes based on the extracted explicit hydrological patterns, and then utilizes the solved generalized ordinary differential equation (UDE) and explicit Euler iterative algorithm (…). After applying two-way physical security constraints, it continuously outputs time-series prediction results of regional soil moisture content that conform to the water balance mechanism and have high accuracy.

[0155] Example 2: Model Engineering Deployment and Parallel Training Strategies

[0156] In practical engineering applications and large-scale data inversion scenarios, to further improve the training and inference efficiency of the HPA-UDE model described in Example 1, this example discloses a system deployment strategy:

[0157] First, regarding model structure compilation, we enable graph compilation optimization for deep learning frameworks (such as PyTorch), and perform modular graph compilation for the static hypernetwork GeoHyperNet and the dynamic flux network FiLMFluxNN, completely avoiding Python interpreter blocking during execution of complex dynamic control flows.

[0158] Secondly, regarding hardware computing power adaptation, a hierarchical cross-validation data loader is constructed, ensuring that the training, validation, and test sets contain data features of all typical terrains by using cluster IDs. During training, `bfloat16` automatic mixed precision is enabled in conjunction with gradient accumulation (e.g., the accumulation step count is set to 4) to achieve a large training batch size with relatively low GPU memory consumption, thereby ensuring the gradient stability of the complex physical system during the optimization process.

[0159] Finally, in training scheduling, a physical constraint progressive injection strategy is adopted. That is, in the first N epochs, only data-driven scheduling is used. Guide the model to learn the general data distribution, and then gradually scale it up. and Activating the physical constraint loss allows the model to smoothly transition between accuracy and physical rigor, avoiding severe gradient oscillations caused by the physical equations in the initial stage.

[0160] The above embodiments only illustrate preferred embodiments of the present invention, and their descriptions are relatively specific and detailed, but they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications, improvements, and substitutions without departing from the concept of the present invention, and these all fall within the protection scope of the present invention.

Claims

1. A soil moisture dynamic inversion method based on hybrid physical perception general differential equation, characterized in that, include: Step S1: Acquisition, preprocessing and K-means clustering of multi-source spatiotemporal data: Acquire static geographic attribute data and dynamic meteorological time series data of the study area. After spatiotemporal alignment and standardization, K-means clustering algorithm is used to divide the static geographic attribute data and obtain hydrological functional zones and their corresponding clustering zone labels. Step S2: Generation of physical parameters based on static hypernet: Construct a static hypernet GeoHyperNet, taking the static geographic attribute data from step S1 and its corresponding clustering partition labels as input, and combining the physical prior bias mechanism and spatial nonlinear mapping to generate soil hydrological parameters with clear physical meaning. The soil hydrological parameters include at least saturated hydraulic conductivity Ksat, maximum water storage capacity Smax, and nonlinear porosity index c_exp. Step S3: Dynamic flux modeling based on spatiotemporal network and feature modulation: The dynamic meteorological time series data is feature-encoded using a temporal Kolmogorov-Arnold network. The extracted time series features and the static features output by the static supernetwork are used as inputs to the dynamic flux network FiLMFluxNN. The dynamic flux parameters are output through the physical head PhysicsHead. The dynamic flux parameters include at least the actual evapotranspiration efficiency and drainage correction term. Step S4: Dynamic estimation of soil moisture based on the general differential equation of hybrid physical sensing: Based on the principle of mass conservation, a general differential equation UDE based on the water balance equation is constructed, and the meteorological driving input, evapotranspiration consumption and drainage flux are used as the rate of change of the derivative of the equation; the soil hydrological parameters generated in step S2 and the dynamic flux parameters generated in step S3 are substituted into the general differential equation, and the solution is iteratively obtained along the time step by the explicit Euler integral algorithm to deduce the continuous predicted value sequence of soil moisture at the current time and the next time step. Step S5: Model training based on physical constraint joint loss: Construct a joint loss function that integrates data-driven error, mass conservation constraint and flux boundary constraint, calculate the loss gradient using the soil moisture prediction value obtained in step S4, and perform backpropagation joint training optimization on the static supernetwork, dynamic flux network and general differential equation module. Step S6: Explicit extraction and prediction loop of hydrological equations based on native symbolic regression: After the model training converges, the learned network weights are extracted from the Kolmogorov-Arnold network layer in the static supernetwork and spatiotemporal network. Through polynomial reconstruction and threshold pruning, an explicit analytical expression of hydrological evolution and the corresponding feature importance map are generated. The analytical expression is used to replace the neural network to realize the real-time prediction of soil moisture, forming a complete soil moisture prediction loop.

2. The soil moisture dynamics inversion method based on the hybrid physical sensing general differential equation according to claim 1, characterized in that, The processing procedure for the static hypernetwork GeoHyperNet in step S2 includes: The static geographic attribute data is divided into soil attribute branch, terrain feature branch, geographic coordinate branch, and clustering condition branch; For the geographic coordinate branch, a multi-frequency sine and cosine basis function is used to encode the geographic coordinates. The encoding result is input into the spatial KAN encoder to learn the hidden spatial periodic features and nonlinear spatial interaction relationships. For the clustering conditional branch, the clustering partition label is converted into a one-hot code, and a preset physical bias benchmark value is obtained through a matrix multiplication lookup table mechanism. A learnable bias modulation layer is then superimposed, and hydrophysical prior knowledge is injected. For the topographic feature branch, a gated modulation signal is generated through the KAN layer, and this signal is used to adaptively weight and fuse the soil property branch features processed by the KAN layer. The output features of the four branches are concatenated along the feature dimension and input into the deep KAN network to generate initial parameter logit values. The initial parameter logit values ​​are added to the physical prior benchmark values ​​and learnable bias terms output by the clustering conditional branches, and then mapped to absolute physical dimensions through a safety boundary constraint layer composed of Softplus and Sigmoid activation functions to obtain soil hydrological parameters with clear physical meaning.

3. The soil moisture dynamics inversion method based on the hybrid physical sensing general differential equation according to claim 1, characterized in that, The dynamic flux modeling in step S3 is specifically... include: Time feature coding: Extract time coding information from dynamic meteorological time series data, retain periodic sine and cosine variables as annual cycle benchmark coding, and fit the residual correction amount of seasonal asymmetric fluctuations through Chebyshev polynomial KAN layer; Feature linear modulation: The encoded temporal features are input into the dynamic flux network FiLMFluxNN, which uses the state-space model Mamba as its backbone network and simultaneously receives scaling parameters generated by the projection of the static hypernetwork. Translation parameters The dynamic features are softened residual modulation through a characteristic linear modulation mechanism; Flux parameter output: The modulated feature input physical head, PhysicsHead, employs a hierarchical fully connected network combined with an activation function to output the actual evapotranspiration efficiency constrained by the Sigmoid function. And the non-negative drainage correction term constrained by the Softplus function. .

4. The method for soil moisture dynamic inversion based on the general differential equation of hybrid physical sensing according to claim 1, characterized in that, The calculation process of the general differential equation UDE in step S4 is as follows: Construct the ordinary differential equation of water balance Where S is the soil moisture content, For precipitation replenishment, For actual evaporation, For drainage flux; Among them, actual evaporation PET represents the potential evaporation of the input. The actual evapotranspiration efficiency generated by step S3; Drainage flux ,in The current soil moisture level. For saturated hydraulic conductivity, For maximum water storage capacity, The porosity index is nonlinear; the explicit Euler method is used along the time step. Solving the ordinary differential equation, the soil moisture content at the next time step is: By imposing hard physical constraints on the preliminary calculated values, the final soil moisture content for the next time step is obtained. The formula is:

5. The method for soil moisture dynamic inversion based on the general differential equation of hybrid physical sensing according to claim 1, characterized in that, The joint loss function in step S5 is the sum of three parts: rate-of-change weighted Huber loss, stepwise mass conservation loss, and flux boundary constraint loss, where: Rate of change weighted Huber loss: Calculating the true soil moisture observations The absolute change at adjacent time points Based on this, an adaptive rate of change weight is constructed. Multiply the weight by the Huber loss between the predicted and observed values; Stepwise mass conservation loss: Calculating the mean square error of soil moisture change and net flux input. ; Flux Boundary Constraint Loss: Minimizing the Mean Function The constraint is that the actual evapotranspiration is non-negative and not greater than the potential evapotranspiration.

6. The method for soil moisture dynamic inversion based on the general differential equation of hybrid physical sensing according to claim 1, characterized in that, The explicit extraction of physical equations in step S6 specifically includes: For a KAN layer that has converged during training, obtain its linear basis weight matrix. and Chebyshev polynomial weight tensor ; Based on the function of the Kolmogorov-Arnold network, the mapping of each input feature is reconstructed into a univariate symbolic expression. ,in for Chebyshev polynomials of the first kind of order; The symbolic expression is traversed and parsed, and only addition sub-items with absolute value coefficients greater than a preset pruning threshold are retained, while redundant nodes are removed. Output explicit analytical expressions for hydrological evolution and feature importance maps.

7. A soil moisture dynamic inversion system based on a hybrid physical sensing general differential equation, used to execute the soil moisture dynamic inversion method according to any one of claims 1 to 6, characterized in that, include: The preloading and data management module is used to acquire static geographic data and dynamic time-series meteorological data, complete data preprocessing and preloading, and build a hierarchical cross-validation data loader. The GeoHyperNet module, which incorporates a spatial KAN encoder and a hydrological clustering prior lookup table, is used to fuse static geographic features with hydrophysical prior biases and output soil hydrological parameters with clear physical meaning. The dynamic flux neural state space module, which includes a time-based KAN encoder and a feature linear modulation network, is used to capture the long-range dependencies of meteorological time-series data and output physical process correction variables for evapotranspiration and drainage. The explicit Euler general differential equation calculation module is used to solve the soil moisture change process based on the water balance differential equation and the explicit Euler integral method, and output the time series prediction value of soil moisture. The multidimensional joint physical loss calculation module is used to construct a joint loss function that integrates adaptive weighted rate of change, mass conservation constraints, and flux boundary constraints, providing gradient descent directions for model training; The symbolic regression and interpretability extraction module is used to parse the multinomial weights of the KAN layer, extract explicit mathematical expressions through symbolic regression, and output explicit formulas for hydrological evolution and feature importance maps.

8. The soil moisture dynamic inversion system based on the hybrid physical sensing general differential equation according to claim 7, characterized in that, The GeoHyperNet module, which is a static physical constraint hypernetwork, adopts a four-branch parallel architecture, namely, soil attribute branch, terrain feature branch, geographic coordinate branch, and clustering condition branch. After each branch completes the extraction and processing of its corresponding features, the features are fused by a deep KAN network to output soil hydrological parameters that conform to the constraints of absolute physical dimensions.

9. The soil moisture dynamic inversion system based on the hybrid physical sensing general differential equation according to claim 7, characterized in that, The dynamic flux neural state space module uses the Mamba state space model as its backbone network and combines it with a feature linear modulation mechanism to achieve residual-type soft modulation of dynamic temporal features by static geographic features.

10. The soil moisture dynamic inversion system based on the hybrid physical sensing general differential equation according to claim 7, characterized in that, The symbolic regression and interpretability extraction module has a built-in algebraic calculation engine, which is used to simplify the reconstructed symbolic expression, prune redundant terms, and finally output an explicit hydrological analytical expression and feature importance ranking map in standard mathematical typesetting format.

Images