A permeability coefficient testing device and testing method

Abstract

The application discloses a permeability coefficient testing device and testing method, and belongs to the technical field of geotechnical engineering and hydrogeological testing. The method comprises the following steps: constructing a hydrological rhythm curve benchmark library through empirical wavelet transform and a variational self-encoder, and extracting a target water supply rhythm curve; generating an executable water supply curve and identifying an initial value by using a multi-scale layered Fourier feature physical information neural network; performing multi-mode flow closed-loop tracking of rolling time domain optimization and adaptive sliding mode feedback; performing multi-physical field collaborative loading by using a three-layer hierarchical fuzzy rule base; obtaining an unsaturated permeability coefficient function and a confidence interval thereof by means of deep operator network transfer learning and Bayesian inversion; and automatically generating a multi-dimensional test report. The application realizes high hydrological rhythm reproduction, high control tracking accuracy, multi-field collaborative reality and inversion result uncertainty quantification.

Description

Technical Field

[0001] This invention relates to the fields of geotechnical engineering and hydrogeological testing technology, specifically to a permeability coefficient testing device and testing method. Background Technology

[0002] The unsaturated permeability coefficient is a key hydraulic parameter characterizing the ability of water to migrate in unsaturated soils. It directly affects the analytical accuracy of engineering and scientific problems such as slope stability assessment, groundwater resource evaluation, pollutant migration prediction, and agricultural irrigation system design. However, accurately obtaining this parameter has always been a challenge in the field of geotechnical testing.

[0003] Traditional unsaturated permeability testing devices often employ constant head or constant flow water supply methods, resulting in relatively simple boundary conditions that fail to simulate the real hydrological rhythms of drastic fluctuations in rainfall intensity and frequent alternations between wet and dry conditions in nature. When a soil sample undergoes a complete cycle from a dry state to heavy rain infiltration, and then from slow evaporation to re-wetting, the internal water transport process exhibits significant hysteresis effects and nonlinear characteristics. A constant water supply mode often fails to elicit these dynamic responses, leading to measured permeability coefficients that are difficult to reflect the hydraulic behavior of the soil under actual field conditions.

[0004] Furthermore, existing testing methods also have shortcomings in terms of control accuracy and testing efficiency. Accurate tracking control of the water supply flow rate is a prerequisite for ensuring the reliability of the test boundary conditions. However, current methods mostly employ conventional proportional-integral-derivative (PID) regulation, which is prone to overshoot, oscillation, or tracking lag when dealing with large-scale, rapidly changing reference flow trajectories. This is especially true at the steep rising and falling edges of the hydrological rhythm curve, where the flow execution error increases significantly. Simultaneously, the coordinated loading control of the temperature and evaporation fields is also a complex process with multiple inputs and multiple outputs and strong coupling. Existing independent loop control methods struggle to effectively decouple the interactions between temperature, evaporation, and water supply, easily causing setpoint deviations in one physical field to be amplified and transmitted to other loops, affecting the accuracy of the spatiotemporal distribution of soil sample moisture content.

[0005] In the identification of permeability coefficients, traditional methods are typically based on simplified analytical models and parameter estimation through manual fitting or local optimization algorithms. These methods are not only sensitive to initial values ​​and prone to getting trapped in local optima, but also struggle to provide confidence intervals and uncertainty measures for the inversion results, making it difficult for engineers to objectively assess the reliability of the test results. While deep learning surrogate models offer the possibility of efficient inversion, systematic technical solutions are still lacking for migration fine-tuning strategies for specific soil samples and Bayesian inversion architectures that incorporate physical constraints.

[0006] Therefore, there is an urgent need for a permeability coefficient testing device and method that can highly reproduce natural hydrological rhythms, achieve high-precision coordinated control of water supply flow and multiple environmental physical fields, and provide reliable uncertainty quantification inversion results. Summary of the Invention

[0007] The purpose of this invention is to provide a permeability coefficient testing device and method, which achieves high reproduction of hydrological rhythms, high control tracking accuracy, real multi-field coordination, and quantification of uncertainty in inversion results.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A permeability testing device, comprising: The sample loading and sealing unit is used to contain the sample and provide a one-dimensional vertical permeation path and temperature control boundary. A water supply flow generation and closed-loop control unit is used to provide a controllable flow rate of water to the bottom of the sample; The environmental simulation loading unit is used to apply radiation, airflow, and temperature environments to the sample. The multi-parameter sensing and measurement unit includes a moisture content-temperature composite sensor array, a sample top weighing sensor, an air temperature and humidity sensor, and a pipeline pressure sensor. Piping switching and drainage collection unit; An integrated control and data acquisition unit, comprising an industrial computer and a high-speed data acquisition board and a programmable logic controller connected thereto, wherein the industrial computer is configured with: The hydrological rhythm digitization extraction module is used to perform empirical wavelet transform eigenmode decomposition and variational autoencoder reconstruction, build a digital hydrological rhythm curve benchmark library, and extract the target water supply rhythm curve from it; The physical information neural network pre-modeling module is used to generate an executable target water supply curve using a multi-scale hierarchical Fourier feature physical information neural network, and to identify the initial values ​​of soil sample hydraulic parameters during the training process. The multi-mode closed-loop flow tracking control module is used to generate composite control commands based on the executable target water supply curve and the measured flow rate fed back by the Coriolis mass flow meter, through rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so as to drive the variable frequency booster pump group and the electric proportional regulating valve. The multi-physics field collaborative environment loading module is used to generate four control commands for the environment simulation loading unit by using a three-layer hierarchical fuzzy rule base to perform multi-objective decoupling coordination based on the evaporation capacity setting and temperature setting extracted from the digital hydrological rhythm curve. The transfer learning permeability coefficient inversion module is used to utilize the measured flow rate and the measurement data of the multi-parameter sensing and measurement unit to perform transfer learning and Bayesian uncertainty inversion based on the pre-trained cross-attention deep operator network proxy model, and obtain the unsaturated permeability coefficient function of the soil sample and its confidence interval. The multidimensional test report generation module is used to summarize the process data of the above modules and generate a test report that includes hydrological rhythm adaptability analysis, water supply execution accuracy assessment, environmental collaborative control quality assessment, and permeability coefficient inversion reliability assessment.

[0009] The present invention also provides a testing method for a permeability testing device, applied to the aforementioned permeability testing device, comprising: Step 1: Obtain hydrological data of the target sampling area, construct a digital hydrological rhythm curve benchmark library through empirical wavelet transform and variational autoencoder, and extract the target water supply rhythm curve from it; Step 2: Based on the physical properties of the soil sample, a multi-scale hierarchical Fourier feature physical information neural network is used for course learning and training to generate an executable target water supply curve. At the same time, the initial values ​​of the hydraulic parameters of the soil sample are identified during the training process. Step 3: Based on the executable target water supply curve and the real-time measured water supply flow, perform multi-mode closed-loop flow tracking control, including rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so that the flow entering the bottom of the sample tracks the executable target water supply curve. Step 4: Based on the evaporation capacity and temperature settings extracted from the digital hydrological rhythm curve, and the real-time measured sample moisture content and temperature, multi-physics field collaborative loading control is carried out through a three-layer hierarchical fuzzy rule base to adjust the environmental simulation loading unit to achieve collaborative simulation of temperature, evaporation and water supply. Step 5: Using the measured water supply flow rate, internal volumetric water content and temperature of the sample, and environmental actuator status data throughout the entire time period, the pre-trained deep operator network proxy model is fine-tuned through transfer learning, and the unsaturated permeability coefficient function of the soil sample and its uncertainty quantification results are obtained by using the Bayesian approximate inversion method. Step 6: Summarize the process data from all steps and automatically generate a multi-dimensional test report that includes hydrological rhythm adaptability, water supply execution accuracy, environmental collaborative control quality, and permeability coefficient inversion reliability.

[0010] Further, in step 1, hydrological data of the target sampling area is acquired, and a digital hydrological rhythm curve benchmark library is constructed using empirical wavelet transform and variational autoencoder. The target water supply rhythm curve is then extracted from this library. Specifically: After normalizing the sequences of rainfall intensity, evaporation and surface runoff depth, empirical wavelet transform is performed on each sequence to extract the two intrinsic mode functions with the highest energy proportion. These functions are then input into a bidirectional long short-term memory variational autoencoder for reconstruction. The reconstructed modes are then added together to obtain the hydrological rhythm fundamental wave. The residuals between the original sequence and the fundamental wave are used to generate random perturbation terms using a second-order Markov chain. The fundamental wave and the perturbation terms are added together to form a digital hydrological rhythm curve and stored in the reference library.

[0011] Furthermore, in step 2, the course learning and training is divided into three stages: The first phase uses constant flow boundary data for training; The second phase uses variable flow boundary data for training; The third stage uses the rainfall component in the digitized hydrological rhythm curve as unlabeled boundary data for unsupervised training. The training loss function includes a data fitting term, a physical residual term, a boundary condition term, and an initial condition term, and uses a residual adaptive weighting strategy to dynamically adjust the weights of each loss term.

[0012] Furthermore, in step 3, the multi-mode closed-loop flow tracking control specifically includes: Before each test, the liquid capacity and liquid resistance of the water supply pipeline were identified online through short impulse response experiments, and a lumped parameter hydraulic transmission model was established. In each control cycle, a lumped parameter hydraulic transfer model is connected in series with a physical information neural network to form a differentiable mapping, and a feedforward control sequence is solved to minimize the control increment penalty while ensuring that the predicted flow in the future multi-step flow tracks the reference trajectory. A non-singular terminal sliding surface is constructed using the flow tracking error, and combined with the disturbance estimation based on one-step delay observation, a feedback compensation quantity containing equivalent control and super-spiking switching control is generated. The sum of feedforward and feedback is decomposed into low-frequency and high-frequency components according to the frequency band, which are used to adjust the speed of the variable frequency booster pump and the opening of the electric proportional regulating valve, respectively.

[0013] Furthermore, in step 4, the three-layer hierarchical fuzzy rule base includes: The upper-level operating condition identification layer divides the system operating status into six macro-operating conditions based on the current water supply flow, moisture content and environmental settings: heavy rain infiltration period, moderate rainfall infiltration period, light rain wetting period, post-rain evaporation period, long-term drought period and intermittent recovery period. The middle-layer multi-objective decoupling coordination layer takes temperature deviation, moisture content deviation and evaporation rate deviation as inputs and outputs temperature correction, moisture content correction and evaporation capacity correction. Its fuzzy rule consequent is corrected online by the physical information neural network based on the prediction deviation. The underlying actuator allocation layer maps the target correction amount to the execution instructions of the circulating water bath jacket, infrared radiation plate, fan and radiant heat source through a preset allocation matrix; The membership function of the fuzzy rule base adopts an asymmetric Gaussian type and is wrapped with an interval type II fuzzy set to handle sensor noise and model uncertainty.

[0014] Furthermore, the online gradient correction of the fuzzy rule consequent of the middle-layer multi-objective decoupling coordination layer is specifically as follows: on the set of active rules where the excitation intensity exceeds a preset threshold, a loss function is constructed based on the deviation between the moisture content and temperature of the next control cycle predicted by the physical information neural network and the set value. The gradient of this loss with respect to the rule conclusion value is calculated through backpropagation, and the rule conclusion value is updated along the gradient descent direction.

[0015] Furthermore, in step 5, the deep operator network proxy model consists of a parameter encoding branch, an environment context encoder, a spatiotemporal coordinate backbone network, and a cross-attention fusion module; during the transfer learning fine-tuning, the residual block weights and Fourier mapping frequency coefficients of the spatiotemporal coordinate backbone network are frozen, and the partially fully connected layers of the parameter encoding branch, the cross-attention module, the output head, and the environment context encoder are trained with different learning rates.

[0016] In summary, the present invention has at least one of the following beneficial technical effects: Firstly, by constructing a digital hydrological rhythm extraction mechanism, the measured hydrological sequence is transformed into a hydrological rhythm benchmark library using empirical wavelet transform and variational autoencoder. This can generate a target water supply rhythm curve that closely approximates the natural rainfall-evaporation alternation characteristics, making the test conditions more realistically reflect the in-situ hydrological dynamics, thereby stimulating the complete hydraulic response of soil samples under different wet and dry states.

[0017] Secondly, by using a multi-scale hierarchical Fourier feature physical information neural network as a digital twin, the saturation change and seepage flow of the sample under the target water supply curve can be predicted before the test. It can automatically determine whether the sample is saturated too early or exceeds the physical limits of the equipment, and perform constraint optimization and correction on the water supply curve. This avoids sample damage or invalid data caused by unreasonable boundary conditions in traditional tests, and improves the success rate and safety of the test.

[0018] Third, a multi-mode closed-loop flow tracking control strategy combining rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control is adopted, supplemented by a meta-learning regulator based on gated cyclic units to dynamically adjust the controller parameters. This strategy can maintain high-precision flow reproduction capability under a wide range of variable flow conditions and effectively reduce the impact of disturbances such as atomization and vibration on flow tracking.

[0019] Fourth, temperature, evaporation and water supply are loaded in a coordinated manner through a three-layer hierarchical fuzzy rule base. In the middle layer, online correction of rule consequents based on physical information neural networks and uncertainty processing of interval type II fuzzy sets are introduced, which effectively decouples and dynamically compensates the coupling effects between multiple physical fields. The evaporation intermittent enhancement mechanism further improves the simulation fidelity of the evaporation process.

[0020] Fifth, based on the pre-trained cross-attention deep operator network proxy model, the model is fine-tuned to fit specific soil samples through transfer learning and combined with the Bayesian approximate inversion method. This not only efficiently solves the posterior mean of the parameters of the unsaturated permeability coefficient function, but also outputs the confidence interval of each parameter and the uncertainty heatmap of the whole time and space, providing a quantitative credibility basis for engineering decision-making.

[0021] Sixth, the aforementioned functional modules are highly integrated into a unified device, forming a fully automated closed loop from hydrological rhythm extraction, pre-modeling verification, multi-mode flow control, multi-physics field collaborative loading to parameter inversion and report generation, which greatly reduces manual intervention and improves testing efficiency and data consistency. Attached Figure Description

[0022] Figure 1 This is a schematic flowchart of a test method for a permeability testing device according to an embodiment of the present invention; Figure 2 Comparison of unsaturated permeability coefficient function curves and uncertainty envelopes for different soil samples; Figure 3 A schematic diagram comparing the target water supply flow rate curve and the measured flow rate curve for a silty clay sample over all time periods. Figure 4 The absolute tracking error distribution of silty clay samples over the entire time period is shown. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention. Furthermore, the technical features involved in the various embodiments of this invention described below can be combined with each other as long as they do not conflict with each other.

[0024] This invention provides a permeability coefficient testing device, which executes the methods described in steps 1 to 6, realizing fully automated testing from digital hydrological rhythm extraction, physical information neural network pre-modeling, multi-mode flow closed-loop tracking, multi-physics field collaborative environmental loading, to permeability coefficient inversion based on transfer learning and multi-dimensional report output. The device consists of six units: a sample loading and sealing unit, a water supply flow generation and closed-loop control unit, an environmental simulation loading unit, a multi-parameter sensing and measurement unit, a comprehensive control and data acquisition unit, and a pipeline switching and drainage collection unit. These units will be described in detail below: 1. Sample loading and sealing unit The unit includes a vertically arranged stainless steel sample tube. Inside the tube, from bottom to top, are arranged a bottom porous permeable plate, a sample, a top porous permeable plate, and a top cap. Sensor mounting holes are located at different heights on the side wall of the sample tube for inserting a moisture content-temperature composite sensor array. A spiral circulating water bath jacket surrounds the outer wall of the sample tube. This jacket is connected to an external constant-temperature circulating water bath via piping to control the temperature of the sample's side surface. The bottom of the sample tube is connected downstream of a Coriolis mass flow meter via piping. An infrared radiation plate and a controllable-speed fan are positioned directly above the top cap, exposing the upper surface of the sample to controlled radiation and airflow to simulate evaporation. Other structures can be used depending on specific requirements, as long as the functionality is achieved; this application does not impose any limitations on this.

[0025] 2. Water supply flow generation and closed-loop control unit This unit includes a variable frequency booster pump set, an electric proportional control valve, a buffer pressure stabilizing tank, and a Coriolis mass flow meter. The inlet of the variable frequency booster pump set is connected to an external distilled water source, and the outlet is connected to the buffer pressure stabilizing tank via a pipeline connected in series with the electric proportional control valve. The outlet of the buffer pressure stabilizing tank is connected to the Coriolis mass flow meter, and then the water flows through the downstream pipeline of the flow meter into the bottom inlet of the sample loading and sealing unit. The Coriolis mass flow meter measures the instantaneous mass flow rate through its pipeline in real time and outputs a frequency... The flow signal is sent to the integrated control and data acquisition unit. The integrated control and data acquisition unit operates a multi-mode closed-loop flow tracking control module, generating a composite control command that is decomposed into a pump speed signal and a regulating valve opening signal, which drive the variable frequency booster pump and the electric proportional regulating valve, respectively. Other structures can also be used depending on specific requirements, as long as the function is achieved; this application does not impose any limitations on this.

[0026] 3. Environment simulation loading unit This unit includes an infrared radiation plate, a controllable-speed fan, a controllable-power radiant heat source, and a circulating water bath jacket. The infrared radiation plate is installed directly above the top cap of the sample, radiating heat downwards onto the sample surface. The controllable-speed fan is positioned above and to the side of the sample, generating airflow that sweeps across the sample surface to enhance evaporation. The controllable-power radiant heat source is located to the side of the sample cylinder, providing auxiliary heat load. The circulating water bath jacket provides or absorbs heat to regulate the temperature boundary of the sample sidewall. The control commands for the above actuators are generated by the multi-physics collaborative environment loading module executed by the integrated control and data acquisition unit. Other structures can also be used depending on specific requirements, as long as the functionality is achieved; this application does not impose any limitations on this.

[0027] 4. Multi-parameter sensing and measurement unit This unit includes: a moisture content-temperature composite sensor embedded in the middle and at multiple depths above and below the sample sidewall, whose measurement principle combines time-domain reflectometry with a thermistor to simultaneously output volumetric moisture content and temperature signals; a weighing sensor installed on the top of the sample to monitor the total mass change of the sample in real time and calculate the evaporation rate from the rate of weight loss; an air temperature and humidity sensor placed near the infrared radiation plate to acquire the ambient air temperature and humidity to calculate the evaporation driving potential and long-wave radiation correction; and pressure sensors installed in the water supply and drainage pipelines to monitor the inlet and outlet water pressures. All sensor signals are synchronously sent to the integrated control and data acquisition unit after analog-to-digital conversion, with a sampling frequency of not less than [missing information]. .

[0028] 5. Integrated Control and Data Acquisition Unit This unit serves as the control center and data aggregation node of the device, primarily comprising an industrial computer, high-speed data acquisition boards, a programmable logic controller (PLC), and communication interface modules. The industrial computer has six built-in software functional modules, corresponding to steps 1 through 6 respectively: The hydrological rhythm digitization extraction module is used to execute step 1, realizing empirical wavelet transform eigenmode decomposition, bidirectional long short-term memory variational autoencoder reconstruction, Markov chain random perturbation fitting, and the construction and retrieval of a digital hydrological rhythm curve benchmark library. The Physical Information Neural Network Pre-modeling and Digital Twin Module is used to execute step 2, realizing the construction of a multi-scale hierarchical Fourier feature physical information neural network, course learning and training, digital twin prediction and executability verification, and curve correction. The multi-mode closed-loop flow tracking control module is used to execute step 3, realizing rolling time-domain feedforward optimization, adaptive non-singular terminal sliding mode feedback control, meta-learning adjustment based on gated loop unit, and frequency band allocation and closed-loop execution of composite control commands. The multi-physics field collaborative environment loading module is used to execute step 4 to realize online reasoning of the three-layer hierarchical fuzzy rule base, online correction of rule consequents based on physical information neural network, uncertainty processing of interval type II fuzzy set, and enhanced control of evaporation interval; The transfer learning penetration coefficient inversion module is used to execute step 5, realizing the transfer learning fine-tuning of the cross-attention deep operator network proxy model, the two-layer optimization of the Bayesian uncertainty-aware inversion loss function, and the reconstruction of the penetration coefficient function. The multidimensional test report generation module is used to execute step 6, and realize the automatic generation of hydrological rhythm adaptability analysis, water supply execution accuracy assessment, environmental collaborative control quality assessment, permeability coefficient inversion reliability assessment, and comprehensive conclusions and data archiving index.

[0029] The high-speed data acquisition board synchronously acquires analog and digital signals from flow meters and various sensors, and stores them in a real-time database after signal conditioning. The programmable logic controller (PLC) converts control commands output from the industrial computer into industry-standard signals, such as... or The data is output to the variable frequency pump unit, regulating valve, and various environmental load actuators. The industrial computer is also equipped with a monitor to display key test data trends and system status in real time, and automatically generates a multi-dimensional test report in PDF format after the test is completed.

[0030] 6. Piping switching and drainage collection unit This unit, located downstream of the outlet at the bottom of the sample, includes an electrically operated three-way switching valve, a graduated cylinder, and a waste liquid tank. During normal testing, the outlet water is switched to the waste liquid tank. During specific calibration phases requiring measurement of permeate flow rate, the outlet water can be switched to the graduated cylinder for volume verification. This unit also includes an vent valve and a filter to prevent air locks and particulate matter blockage in the pipeline. Other structures may be used depending on specific needs, as long as the functionality is achieved; this application does not impose any limitations on this.

[0031] Finally, the working process of the device is as follows: During the test preparation phase, the samples are loaded into the sample cylinder, and pipeline connections and sensor calibration are completed. The integrated control and data acquisition unit executes the hydrological rhythm digitization extraction module to import or generate hydrological rhythm curves; the physical information neural network pre-modeling and digital twin module generates executable water supply curves based on input soil parameters using a pre-trained physical information neural network; the multi-mode closed-loop flow tracking control module initiates closed-loop tracking control of water supply flow to achieve high-precision reproduction of water supply rhythms; simultaneously, the multi-physics field collaborative environment loading module controls the environmental loading of various actuators to achieve coordinated changes in temperature, evaporation, and water supply, and continuously collects multi-physical quantity data; after the test, the transfer learning permeability coefficient inversion module uses the collected data and transfer learning methods to invert soil hydraulic parameters; the multi-dimensional test report generation module automatically summarizes and generates a multi-dimensional test report.

[0032] like Figure 1 As shown, the present invention also provides a testing method for a permeability testing device, applied to the aforementioned permeability testing device, comprising: Step 1: Obtain hydrological data of the target sampling area, construct a digital hydrological rhythm curve benchmark library through empirical wavelet transform and variational autoencoder, and extract the target water supply rhythm curve from it; Step 2: Based on the physical properties of the soil sample, a multi-scale hierarchical Fourier feature physical information neural network is used for course learning and training to generate an executable target water supply curve. At the same time, the initial values ​​of the hydraulic parameters of the soil sample are identified during the training process. Step 3: Based on the executable target water supply curve and the real-time measured water supply flow, perform multi-mode closed-loop flow tracking control, including rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so that the flow entering the bottom of the sample tracks the executable target water supply curve. Step 4: Based on the evaporation capacity and temperature settings extracted from the digital hydrological rhythm curve, and the real-time measured sample moisture content and temperature, multi-physics field collaborative loading control is carried out through a three-layer hierarchical fuzzy rule base to adjust the environmental simulation loading unit to achieve collaborative simulation of temperature, evaporation and water supply. Step 5: Using the measured water supply flow rate, internal volumetric water content and temperature of the sample, and environmental actuator status data throughout the entire time period, the pre-trained deep operator network proxy model is fine-tuned through transfer learning, and the unsaturated permeability coefficient function of the soil sample and its uncertainty quantification results are obtained by using the Bayesian approximate inversion method. Step 6: Summarize the process data from all steps and automatically generate a multi-dimensional test report that includes hydrological rhythm adaptability, water supply execution accuracy, environmental collaborative control quality, and permeability coefficient inversion reliability.

[0033] Next, the above method will be explained in detail with reference to specific embodiments and parameter settings: In step 1, hydrological data of the target sampling area is acquired, and a digital hydrological rhythm curve benchmark library is constructed using empirical wavelet transform and variational autoencoder. The target water supply rhythm curve is then extracted from this library. Specifically: The input is the measured hourly rainfall intensity sequence for a continuous hydrological year in the target sampling area. Measured sequence of evaporation amount in evaporating dish and measured sequence of surface runoff depth Time Index ,in The output is a digital hydrological rhythm curve. The benchmark library and the target water supply rhythm curve selected from it .

[0034] Step 1 is achieved through the following sub-steps: Step 101: Data Preprocessing Rainfall intensity sequence Evaporation amount sequence of evaporating dish and surface runoff depth sequence Missing value imputation and outlier removal were performed separately. Missing value imputation used linear interpolation between adjacent time points; if a consecutive missing time window exceeded 6 hours, it was filled with the historical average value for the same period over many years. Outlier removal followed the three-standard-deviation criterion; observations exceeding three standard deviations of the mean were considered outliers and replaced with the linear interpolation result. The processed sequences are denoted as follows: , , Then, perform max-min normalization on each dimension of the sequence to map it to the interval. The normalized three-dimensional sequence is denoted as .

[0035] Step 102: Empirical wavelet transform eigenmode decomposition right The three dimensions of the sequence are subjected to empirical wavelet transforms, adaptively decomposing the normalized sequence into a set of eigenmode functions. Taking the rainfall dimension as an example, let the normalized rainfall sequence be denoted as... First, calculate its Fourier spectrum. And an adaptive spectrum segmentation method is used to determine Frequency band boundary ,in , The segmentation is based on the detection of local minima in the Fourier spectrum. The spectral interval between two adjacent local minima is defined as a frequency band, and the total number of frequency bands is [number missing]. .

[0036] For each frequency band Constructing empirical wavelet functions For the lowest frequency band Additional construction of empirical scaling function The frequency domain expression of the empirical wavelet function is: ; The frequency domain expression of the empirical scaling function is: ; in, An auxiliary function to satisfy the compact support condition, For transition bandwidth, Pick For normalized rainfall sequences Perform empirical wavelet transform to obtain the first... Modal functions The result is the inverse Fourier transform of the detail coefficients. The low-frequency approximate component is denoted as Energy proportion of each mode The calculation formula is: ; Repeat the above process for the evaporation and runoff dimensions to obtain the mode sets for each dimension. , These correspond to the dimensions of rainfall, evaporation, and runoff, respectively.

[0037] Step 103: Principal Mode Extraction and Bidirectional Long Short-Term Memory Variational Autoencoder Reconstruction For each dimension, select the energy percentage. The two highest-order modes are used as the principal modes for subsequent reconstruction, while the remaining modes are discarded as high-frequency noise. The selected principal mode waveforms are then time-aligned and concatenated into a two-dimensional time series. , which serves as the input to a pre-constructed bidirectional long short-term memory variational autoencoder.

[0038] The bidirectional long short-term memory variational autoencoder consists of three parts: an encoder, a latent space sampling layer, and a decoder. The encoder's first layer is a bidirectional long short-term memory layer with an input size of 2 and 64 hidden units, while the second layer is a bidirectional long short-term memory layer with 32 hidden units. The second layer outputs the forward and backward hidden states at the final time step, which are mapped to the mean vector of the latent space through two independent fully connected layers. Sum of logarithmic variance vector Latent vector Sampling via reparameterization techniques: ,in , This represents element-wise multiplication. The decoder's first layer is a bidirectional long short-term memory layer with 32 hidden units, the second layer is a bidirectional long short-term memory layer with 64 hidden units, and finally, a fully connected layer with a temporal distribution maps the output back to two-dimensional space to obtain the reconstructed modal sequence. .

[0039] The loss function of an autoencoder is defined as a weighted combination of the reconstruction mean square error, the Pearson correlation coefficient penalty, and the KL divergence regularization: ; in For the index set of the two selected primary modes, The weights for the correlation coefficient loss, This is the KL divergence regularization coefficient. The Pearson correlation coefficient is calculated as follows: ; KL divergence term The calculation method is as follows ; in and They are respectively and The There are several components. Training uses an adaptive moment estimation optimizer with an initial learning rate of [missing information]. The batch size is 128, the maximum number of iterations is 200, and an early stop mechanism is set up to terminate training if the loss does not decrease for 10 consecutive iterations.

[0040] Step 104: Hydrological Rhythm Fundamental Wave Synthesis After training convergence, the corresponding reconstructed modes are obtained by forward inference of the two principal modes in each dimension through a variational autoencoder. During decoding, the latent vector is taken. To obtain a deterministic output, the reconstructed modes are directly added in the time domain to obtain the fundamental component of the hydrological rhythm for that dimension. Taking the precipitation dimension as an example, the fundamental wave... Represented as: ; Evaporation fundamental wave and runoff fundamental wave Generated in the same manner. The time alignment vectors of the three constitute the fundamental wave of the hydrological rhythm. .

[0041] Step 105: Fitting the random perturbation term The original normalized sequence Subtract the fundamental wave of the hydrological rhythm Obtain the residual sequence Residual sequences for the rainfall dimension A second-order Markov chain rainfall generator is used for random perturbation modeling. Firstly, a second-order Markov chain rainfall generator is used for modeling. The value range is divided into equal-frequency binning methods as follows: There are several states, each corresponding to a residual amplitude interval. The second-order transition frequencies from two consecutive states to the next state are counted, and then normalized to obtain the second-order state transition probability tensor. tensor elements The states at the first two time points are respectively: and The state at the next moment is The probability. Also, in each state... Empirical cumulative distribution function of the magnitude of the inner-fit residual The generation rule for the random perturbation term is as follows: randomly initialize the states at the first two time steps, and then, based on the second-order transition probability tensor... The state sequence is generated sequentially. Then, at each time step, an inverse transform is performed on the empirical cumulative distribution function corresponding to the current state to sample and generate random perturbation values ​​for the residuals. The residuals for the evaporation and runoff dimensions are processed in the same way to finally obtain the sequence of random perturbation terms. .

[0042] Step 106: Construction of Digital Hydrological Rhythm Curves and Extraction of Target Water Supply Rhythm Curves By adding the fundamental wave of the hydrological rhythm to the random disturbance term in the time dimension, a digital hydrological rhythm curve is obtained. The target sampling area corresponds to a continuous hydrological year. The system is divided into 365 daily segments based on natural days, with each daily segment stored as an independent entry in the benchmark database. Each record in the benchmark database is marked with a timestamp and a meteorological feature label, which includes four fields: total rainfall for that day, maximum daily rainfall intensity, number of consecutive dry days, and total evaporation.

[0043] Based on the testing objectives, operators retrieved typical time segments from the benchmark database using meteorological feature tags. The selected time segments were then extracted. The rainfall dimension component is then analyzed, and an inverse normalization operation is performed on this component to restore it to the instantaneous flow rate dimension with physical units, i.e., the target water supply rhythm curve. The unit is The inverse normalization process uses the maximum and minimum values ​​of rainfall dimensions saved in step 101 for inverse mapping.

[0044] Step 1 output target water supply rhythm curve This will be used as the upper boundary flow condition input to the physical information neural network in step 2 to generate an executable target water supply curve.

[0045] In step 2, based on the physical properties of the soil samples, a multi-scale hierarchical Fourier feature physical information neural network is used for course learning and training to generate an executable target water supply curve. Simultaneously, the initial values ​​of the hydraulic parameters of the soil samples are identified during the training process. Specifically: The inputs for step 2 include the particle size distribution curve and dry density of the sample. Organic matter content Preset ambient temperature Preset evaporation capacity and initial volumetric moisture content and the target water supply rhythm curve output in step 1. The output is an executable target water supply curve. and the initial values ​​of the identified hydraulic parameters .

[0046] Step 2 is achieved through the following sub-steps: Step 201: Vectorization and Normalization of Material Property Parameters The particle size distribution curve of the sample was discretized into 16 characteristic values, corresponding to the mass percentage of 16 sieve aperture ranges. , with dry density and organic matter content Combined into a soil parameter vector Set the preset ambient temperature Preset evaporation capacity and initial volumetric moisture content Combined into an environmental parameter vector .right and Perform min-max normalization on their respective dimensions and map to The intervals, after normalization, are concatenated into a vector of physical properties and environmental parameters. .

[0047] Step 202: Construct a multi-scale hierarchical Fourier feature physical information neural network The input to a physical information neural network is spatial coordinates. Time coordinates and parameter vector ,in The height of the sample. The total test duration is given. The network consists of a Fourier feature mapping layer, four parallel encoding branches, a multi-head cross-attention fusion layer, a two-stream decoder, and two output heads.

[0048] The Fourier feature mapping layer takes the input spatial coordinates and time coordinates These are respectively mapped to high-dimensional Fourier eigenvectors. Spatial coordinates are mapped as follows: ; Time coordinate mapping ; in and The frequency parameters are randomly sampled, respectively from Gaussian distributions. and Mid-sampling, , , and These are trainable hyperparameters. The two feature vectors have dimensions of 128 and 64 respectively, and are then compared with the normalized parameter vector. Concatenated into total input features .

[0049] The structure of the four parallel coding branches is as follows: Branch 1 is a high-frequency fine-grained branch, receiving... The input is then fed into a residual network consisting of five fully connected layers, each with 128 neurons. Batch normalization and the Swish activation function are used between layers. A skip connection is added between every two fully connected layers, and the output of the first layer is added to the input of the fourth layer before being fed into the fifth layer. This branch uses a larger standard deviation for the frequency parameter. and It is used to capture the steep gradient characteristics of the water transport front during infiltration.

[0050] Branch 2 is a low-frequency coarse-grained branch, with the same structure as Branch 1, but using a smaller standard deviation of the frequency parameters. and It is used to characterize the slowly changing matrix suction field and the overall water content evolution trend.

[0051] Branch 3 is a direct mapping branch of physical properties, which only accepts parameter vectors. The process is handled by four fully connected layers with the number of neurons being 128, 64, 64, and 32 respectively, and the output is a material property embedding vector.

[0052] Branch four is the boundary condition coding branch, which receives the normalized upper boundary flow value. and lower boundary pressure head value Boundary embedding vectors are generated through two fully connected layers to explicitly inject the time-varying features of the boundary conditions into the network.

[0053] The output dimensions of the four branches are uniformly mapped to 128 dimensions through linear projection, resulting in four feature vectors. Stack the four feature vectors as follows: The data is fed into a multi-head cross-attention fusion layer. The multi-head cross-attention fusion layer has 4 heads, and each head has a dimension of 32. For the... Size, embedded vector with physical properties Generate query by linear projection stacked feature matrices Bond generation via linear projection Sum ,in , , is a trainable projection matrix. The first... The attention output of the size is ; The outputs of the four heads are concatenated and then linearly projected to obtain the fused features. The data is then processed through residual connections and a feedforward network. The feedforward network consists of two fully connected layers, with the Swish activation function used in between.

[0054] The dual-stream decoder receives the fused feature vector and simultaneously feeds it into the volumetric water content decoding stream and the pressure head decoding stream. Each decoding stream consists of three fully connected layers with 128, 64, and 1 neuron respectively. The last layer has no activation function and directly outputs the result. or .

[0055] Step 203: Constructing the total loss function The training objective of the network is to simultaneously satisfy data constraints, physical constraints, boundary conditions, and initial conditions. The total loss function is defined as follows: ; in For data fitting loss, For physical residual loss, For boundary condition loss, Loss due to initial conditions , , , These are the weighting coefficients for each loss term.

[0056] Data fitting loss Enabled only when measured data is available, calculates the mean square error between network predictions and measured values. Physical residual loss. The network output must satisfy the one-dimensional Richards equation: ; in For the number of configuration points, and These are the volumetric water content and pressure head output by the network, respectively. The output value is calculated using the van Genuchten-Mualem formula. The van Genuchten-Mualem model consists of the following three equations combined: ; ; ; in For effective saturation, This represents the residual volumetric moisture content. This represents the saturated volumetric water content. It is the reciprocal of the intake value. For aperture distribution parameters, Satisfying Mualem constraints , The saturated permeability coefficient, This refers to the curvature factor. These six hydraulic parameters... Instead of being directly input, it is used as an additional trainable variable in the network and optimized synchronously with the network weights. This allows the network to not only approximate the solution of the partial differential equation during training, but also automatically identify the hydraulic parameters most suitable for the soil sample.

[0057] Boundary condition loss Including upper boundary flow conditions and lower boundary pressure conditions: ; in This is the upper boundary flow used in the current training phase. This represents the lower boundary pressure head.

[0058] Initial condition loss Constraint networks in The output at each time point is consistent with the preset initial conditions: ; The initial values ​​of the weights of each loss term are set to... , , , During training, a residual adaptive weighting strategy is used for dynamic adjustment. Every 100 iterations, the mean of each loss term is calculated. For terms whose mean exceeds 1.5 times the average loss, their weights are multiplied by [a factor not specified in the original text]. To accelerate the convergence of this term, the upper limit of the weight is set to .

[0059] Step 204: Course Learning and Training Strategies and Implementation The training dataset was constructed using a course-based learning strategy, with the difficulty gradually increasing in three stages.

[0060] The first phase uses constant upper boundary flow boundary condition data. It generates 1000 samples of constant head or constant flow conditions with different flow rates, covering a wide range of flow rates. Up to the upper limit of the Coriolis mass flow metering range The duration of each working condition is For each operating condition, the Richards equations were solved using finite element software to generate the volumetric water content. and pressure head The reference solution is used as the supervision label. Only data fitting loss and physical residual loss are enabled in this phase, and training is performed for 10,000 steps.

[0061] The second stage introduces variable flow boundary conditions. 5000 stepped or sinusoidal variable flow sequences with different amplitudes and rates of change are generated, and a reference solution is also generated through finite element analysis. All four loss terms are enabled in this stage, and the training run is performed for 50,000 steps.

[0062] The third stage uses the rainfall dimension components of multiple digitized hydrological rhythm curves randomly selected from the benchmark library in step 1 as the upper boundary flow. No reference solution is generated in this stage; unsupervised training is performed entirely based on physical residual loss and boundary loss, iterating for 30,000 steps.

[0063] Training uses the Adam optimizer with an initial learning rate of Phase Two was reduced to Phase 3 downgraded to Each batch is randomly sampled from the temporal and spatial domains. Each configuration point is used to calculate the physical residual. Each point is used for boundary condition loss. A number of points are used for initial conditional loss. During training, the model with the smallest validation physical residual is saved as the final pre-trained model.

[0064] Step 205: Digital Twin Prediction and Feasibility Verification The pre-trained physical information neural network was used as a digital twin of the soil sample for this test. Before the formal test began, the normalized physical property and environmental parameter vectors were... and target water supply rhythm curve It is fed into the network as input. The network moves along the time axis... Perform forward reasoning for the step size, and compute all spatial discrete points at each time step. and The distribution of unsaturated permeability coefficients across the entire field is then output using the van Genuchten-Mualem formula. Extract the effective saturation at the bottom boundary of the sample. If there are 10 consecutive time steps where the effective saturation exceeds If so, the sample is determined to have prematurely saturated.

[0065] Simultaneously, the bottom seepage flow rate is calculated based on the pressure head gradient output by the network: ; If the seepage flow rate is predicted The value at any given time exceeds the upper limit of the Coriolis mass flow metering range. If the target water supply rhythm curve exceeds the physical limit of the device, it is determined that the target water supply rhythm curve exceeds the physical limit of the device.

[0066] Step 206: Curve Correction and Generation of Executable Target Water Supply Curve If premature saturation or excessive flow occurs during the verification in step 205, then 95% of the maximum flow limit will be used as a constraint. A correction is made. The correction employs the augmented Lagrangian method, transforming the constrained optimization problem into an unconstrained saddle point problem. The corrected curve is defined as follows. Increment function The parameterization is performed by a radial basis function network, which contains 128 Gausky functions with their centers uniformly distributed in the time domain. The width parameter is set to... The optimization objective is to minimize The L2 norm satisfies two constraints simultaneously: all time steps And the effective saturation of the sample does not exceed The augmented Lagrangian function is ; in Penalties for violating the constraints For Lagrange multipliers, The penalty factor is the initial value. Each round of the outer cycle is multiplied by Iterate until the constraints are satisfied or the maximum number of outer loop iterations is reached, then output the corrected curve. If no constraint violation is found during verification, output directly. .

[0067] The executable target water supply curve output in step 2 This will serve as the reference trajectory for the flow tracking control module in step 3. Simultaneously, the hydraulic parameters identified by the pre-trained network... This will be used as the initial prior estimate for the inversion in step 5.

[0068] In step 3, based on the executable target water supply curve and the real-time measured water supply flow rate, multi-mode closed-loop flow tracking control is executed, including rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so that the flow rate entering the bottom of the sample tracks the executable target water supply curve, specifically: The input for step 3 is the executable target water supply curve generated in step 2. And the measured flow rate fed back in real time by the Coriolis mass flow meter. The output is a sequence of composite control commands sent to the variable frequency booster pump unit and the electric proportional regulating valve. and the full-time measured traffic record sequence. The measured sequence is then used in step 5 for permeability coefficient inversion.

[0069] Step 3 is achieved through the following sub-steps: Step 301: Discretization of the flow reference trajectory and online feedforward sequence generation Executable target water supply curve By time step Discretization into flow reference trajectory , ,in The total number of test steps is calculated. The first-order forward difference is also calculated. As prior information about flow changes.

[0070] To generate high-quality feedforward control inputs, a differentiable digital replica of the water supply system based on the pre-trained physical information neural network from step 2 is constructed. A lumped-parameter hydraulic transfer model is connected in series at the end of the physical information neural network to transmit the control inputs. Converted into flow rate input to the sample chamber The discrete-time state-space expression of the lumped-parameter hydraulic transfer model is as follows: ; in The equivalent pressure of the water supply pipeline is expressed in units of... , For normalized control commands, the value range is... , The fluid volume in the pipeline is expressed in units of... , The unit is the hydraulic resistance of the pipeline. . and The short-impulse response was identified online before each test. The identification method was as follows: applying a continuous short-impulse response... The step control command is received, the flow response curve is recorded, and the parameters of the first-order dynamic model are fitted using the least squares method. The physical information neural network receives the flow... As an upper boundary condition, the predicted effluent flow rate at the bottom of the output sample is calculated. This constitutes the control command To the desired outflow rate A complete differentiable mapping.

[0071] In each control cycle Starting with the measured state obtained from the previous cycle feedback, this differentiable mapping is used to predict the future. The rolling prediction is performed step by step, and the feedforward control sequence is obtained by solving the following trajectory tracking optimization problem with penalty: ; in and The positive weighting coefficients are... To control the incremental penalty term, optimization employs a mini-batch gradient descent method for online iteration. The initial value is taken as the time shift value of the sequence from the previous time step, and the number of iterations is fixed at 5. Only the first control variable of the optimized sequence is used. As the current time-of-flight instruction output.

[0072] Step 302: Adaptive Non-Singular Terminal Sliding Mode Feedback Controller Enhanced by Perturbation Observation Feedback compensation is introduced to eliminate parameter perturbations and unmodeled disturbances. Flow tracking error is defined. Constructing non-singular terminal sliding surfaces ; Operators , , , And satisfy . Depend on calculate.

[0073] Lumped disturbances exist in real systems This includes changes in pipeline resistance, actuator zero drift, and external vibration. The dynamics of the device's input and output are described using a nominal model with added disturbance terms: ; Among them, the nominal model The discrete output expression of the hydraulic transmission model in step 301 is taken. A perturbation estimator based on a one-step delayed observation is designed: ; ; This is the observer gain. This is the estimated disturbance value at the current moment.

[0074] Based on the principles of disturbance compensation and sliding mode equivalence, the feedback control law is designed as follows: ; Equivalent control items According to the order Furthermore, ignoring the conditions of disturbance changes, the derivation is based on online identification. The parameters are used to resolve the steady-state control components required for the flow command. The switching control term adopts a super-spiral form. ; ; in and For adaptive gain, its update law is: ; ; , , This adaptive mechanism automatically increases the gain to accelerate convergence when the tracking error is large, and gradually decreases the gain after the error converges.

[0075] Step 303: Meta-learning regulator based on gated recurrent unit Shape parameters in non-singular terminal sliding mode feedback controllers This has a significant impact on control adaptability. A meta-learning regulator based on a gated cyclic unit is designed to... The adjustment increment of shape parameters is dynamically output for the decision cycle.

[0076] The gated recurrent unit regulator consists of a two-layer gated recurrent unit network and two fully connected prediction heads. The first layer of the gated recurrent unit network has 64 hidden units, and the second layer has 32. The input state vector... Designed for ; in The reference flow rate change within the next second is calculated as the difference between the maximum and minimum values ​​within that window. The two output heads output shape parameters respectively. and The multiplicative adjustment factor, and The additive offset, the output range is scaled and limited to: , , .

[0077] The gated cyclic unit regulator was trained using an offline strategy based on model-independent meta-learning. A task distribution was constructed that included different soil sample permeability characteristics and different water supply rhythm patterns. Each task corresponds to a set of soil parameters and a water supply rhythm curve. During the meta-training phase, for each task, the regulator first performs 5 internal gradient updates, with the loss function defined as the root mean square of the tracking error. ; Next, calculate the external meta-loss and update the meta-parameters using the second-order gradient. The meta-learning rate is set to... The internal learning rate is set to After meta-training is complete, the regulator obtains an initialization that is highly sensitive to changes in the task.

[0078] During online operation, the gated cyclic unit regulator performs forward inference every 10 control cycles, outputting updated shape parameters. These updated values ​​are then smoothed using an exponential moving average before being fed into the feedback controller. The smoothing coefficient is [value missing]. .

[0079] Step 304: Assignment and Closed-Loop Execution of Composite Control Commands The composite control command at the current moment is ; This command needs to be assigned to the speed of the variable frequency booster pump unit. and electric proportional control valve opening Adopting the principle of frequency band division of labor: […] It is decomposed into low-frequency components by a first-order low-pass filter. and high frequency components The cutoff frequency is set to The distributive law is ; ; and Physical limit for pump speed, The balance valve position corresponds to the current pump speed. This refers to the valve position sensitivity.

[0080] Coriolis mass flow meter Frequency will be measured flow rate Feedback is provided to each link in the closed loop, forming a complete closed loop. Tracking performance is measured by the normalized root mean square error over all time periods. ; Design specifications are .

[0081] Step 3 outputs the full-time measured flow sequence This data will be used as one of the inputs to the hydraulic parameter inversion module in step 5, together with the sensor array time series, to form the multi-physical quantity dataset required for the inversion. The output parameter variation trajectory of the gated loop unit regulator in the feedback controller is recorded in the test report to illustrate the adaptive process of the control system.

[0082] In step 4, based on the evaporation capacity and temperature settings extracted from the digital hydrological rhythm curve, and the real-time measured sample moisture content and temperature, multi-physics field collaborative loading control is performed through a three-layer hierarchical fuzzy rule base to adjust the environmental simulation loading unit and achieve collaborative simulation of temperature, evaporation, and water supply. Specifically: The input for step 4 includes: the evaporation capacity setting curve obtained by inverse normalization of the evaporation dimension component in the digitized hydrological rhythm curve from step 1. and temperature set curve Step 2: The pre-trained physical information neural network model and its identified hydraulic parameters. Step 3: Real-time feedback of measured water supply flow rate and the water content of the sample's mid-volume, collected in real time by the sensor array. and temperature ,in The measurement location is at the middle of the sample. The height is the sample height. The output consists of a sequence of four control commands sent to the circulating water bath jacket, infrared radiation plate, controllable speed fan, and controllable power radiant heat source, as well as the measured environmental parameters throughout the entire time period. These records are ultimately incorporated into the inversion dataset of step 5 and the test report of step 6.

[0083] Step 4 is achieved through the following sub-steps.

[0084] Step 401: Construction of the water-thermal-vapor coupled dynamic model and establishment of the state space Based on the one-dimensional Richards equation described in step 2, heat conduction and convection mass transfer terms are introduced to construct a water-heat-vapor coupled kinetic model. The coupled governing equations are composed of the water transport equation and the heat transport equation.

[0085] The moisture transport equation is based on the isothermal Richards equation with the addition of a temperature gradient-driven water vapor diffusion flux: ; in This is the temperature-corrected unsaturated permeability coefficient. This is the effective diffusion coefficient of water vapor in soil pores, expressed in units of... , The density of water vapor in the pores, in units of . In step 2, the vanGenuchten-Mualem model Multiply by a temperature correction factor: ; For temperature Dynamic viscosity of wastewater, in units of , Reference temperature The ratio of the two dynamic viscosities was calculated using the Vogel-Fulcher-Tammann equation: ; in , .

[0086] The heat transport equation includes three terms: heat conduction, latent heat of vaporization, and convective heat transport. ; in Soil volumetric heat capacity, in units of The volume-weighted average of each component is used: ; Porosity The volumetric heat capacity of soil particles. This is the volumetric heat capacity of liquid water. This is the volumetric heat capacity of air. Soil thermal conductivity, in units of The Johansen model is adopted: ; ; and These are the thermal conductivity values ​​under saturated and dry conditions, respectively, estimated from soil properties. The latent heat of vaporization of water, in units of Using linear approximation . The term represents the heat flux carried by liquid water convection, of which Darcy velocity, in units of .

[0087] The upper boundary condition considers the surface energy balance: ; in For surface temperature, For air temperature, and These are the surface emissivity and the air emissivity, respectively. The Stefan-Boltzmann constant is... air density, The specific heat of air at constant pressure. For aerodynamic impedance, This refers to the surface evaporation rate. This refers to the heat flux entering the soil. Driven by the gradient between surface moisture content and relative humidity: ; For temperature The saturated water vapor density at that point, The density of water vapor in the air. It is the surface impedance, which is inversely proportional to the surface water content.

[0088] Step 402: Overall Architecture Design of a Three-Layer Hierarchical Fuzzy Rule Base Multiphysics field collaborative loading faces strong nonlinear coupling between temperature, evaporation, and water supply. The three-layer hierarchical fuzzy rule base constructed in this step decomposes the multi-loop coordination problem into three levels: the upper-level working condition identification and mode switching layer, the middle-level multi-objective decoupling and coordination layer, and the lower-level actuator allocation and fine-tuning layer.

[0089] The upper layer is based on the current water supply flow status. Moisture content Based on the external environment settings, the system operation is divided into six macroscopic operating conditions: heavy rainfall infiltration period, moderate rainfall infiltration period, light rain wetting period, post-rain evaporation period, long-term drought period, and intermittent recovery period. Each operating condition corresponds to a different control priority weight vector and an activation subset of the intermediate-level fuzzy rule base.

[0090] The middle layer receives operating condition tags and priority weights from the upper layer, based on temperature deviation. Moisture content deviation and evaporation rate deviation As input to fuzzy inference, the outputs are the target correction quantities for each decoupled physical quantity: temperature correction quantity. Evaporation capacity correction amount and moisture content correction .in The water flow rate is predicted online by the physical information neural network in step 2 based on the current water supply flow curve and initial conditions. The rate of weight loss recorded by the weighing sensor at the top of the sample is converted to units of... .

[0091] The bottom layer receives the target correction amount from the middle layer and distributes it to four actuators: the temperature increment of the circulating water bath jacket. Incremental power of infrared radiating plate Controllable fan speed increment and controllable power radiant heat source power increment .

[0092] Step 403: Parameterized Design of Membership Function for Fuzzy Rule Base The input and output variables of the fuzzy rule base are both fuzzified using asymmetric Gaussian membership functions across their respective universes of discourse. The analytical form of the asymmetric Gaussian membership function is as follows: ; in The membership function center value, and These are the left and right widthening parameters, respectively.

[0093] Each input variable is divided into 7 fuzzy subsets, with the following linguistic labels: NB, NM, NS, ZO, PS, PM, and PB. Temperature deviation is used as the criterion for... For example, its domain of discourse is normalized to The central value of each subset In order The corresponding left side is widened. In order widening on the right side In order Moisture content deviation and evaporation rate deviation The membership function parameters are designed according to the same principle, and the broadening parameters are independently calibrated according to their respective physical dimensions and response characteristics.

[0094] The membership function of the output variable adopts a symmetric Gaussian form, i.e. Similarly, it is divided into 7 subsets, with the central value evenly distributed.

[0095] Step 404: Rule structure, inference mechanism, and online adaptive correction of the mid-level fuzzy rule base The input to the mid-level fuzzy rule base is a three-dimensional deviation vector. The output is a three-dimensional correction vector. The complete rule base contains... There are 10 rules, each with the following structure: ; in For the first Rule No. Fuzzy subset labels of the dimensional input. To output the conclusion value, a single-valued consequent is used.

[0096] The reasoning process employs Mamdani minimal reasoning and weighted average defuzzification. For the current input... , No. Activation strength of the rule Determined by the minimum value of the membership degree of each dimension: ; Defuzzing output is ; An online correction mechanism for rule consequents based on a physical information neural network is introduced. Let the current time be... The physical information neural network uses the current measured state as initial conditions and the environmental boundary conditions corresponding to the current mid-layer output correction as input to predict... Moisture content at time and temperature Define prediction bias: ; Calculate the gradient using the differentiability of physical information neural networks. The gradient is backpropagated through the physical information neural network to its input boundary conditions, and then backpropagated from the boundary conditions to the consequents of the rules. This applies only to rules that satisfy... The active rule performs gradient descent updates: ; Learning rate Set as The update frequency is once every 10 control cycles.

[0097] Step 405: Handling rule uncertainty based on interval type-II fuzzy sets Wrap the interval type II fuzzy set outside the membership function. Stretch parameter. and Expand to interval and The lower bound is taken as the original widened portion. The upper bound is taken as the original expansion. This results in an uncertain coverage area, which is multiplied by a factor of 1.

[0098] For each rule, the activation strength becomes a range value: ; ; The lower bound membership degree Expand the upper limit The corresponding membership value, upper bound membership value Expand the lower boundary The corresponding membership value.

[0099] The defuzzy output is processed by the Karnik-Mendel order reduction algorithm to obtain an interval-type output. Finally, take their arithmetic mean. As a deterministic correction amount for the middle layer output. When the interval width When the preset threshold is exceeded, the system determines that the current uncertainty is too high, automatically reduces the control gain and triggers a warning flag, which is recorded in the test report in step 6.

[0100] Step 406: Low-level actuator allocation and constraint processing The bottom layer distributes the three target correction values ​​output from the middle layer to the four actuators. (Assignment matrix) Derived from physical relationships and corrected through calibration experiments: ; Assign matrix elements The physical meaning is: and These are the contribution coefficients of the circulating water bath jacket and the infrared radiation plate to the temperature correction, respectively. and These are the contribution coefficients of the circulating water bath jacket and the radiant heat source to the moisture content correction, respectively. , and These are the contribution coefficients of the infrared radiation plate, fan, and radiant heat source to the correction of evaporation capacity, respectively.

[0101] After allocation, physical limiting constraints are applied to the outputs of each actuator. Circulating water bath jacket temperature range. Temperature change rate limit Infrared radiating plate power range Fan speed range Radiant heat source power range If the calculation instruction exceeds the limit, a scaling rollback strategy is used to scale the out-of-bounds instruction back to the allowed range.

[0102] Step 407: Execution of Evaporation Intermittent Enhancement and Multiphysics Co-loading According to the rainfall-evaporation alternation rhythm defined by the digital hydrological rhythm curve in step 1, after the water supply flow execution module completes water supply for the current rainfall period and the flow rate returns to zero to enter the intermittent period in step 3, the underlying evaporation-related actuators automatically switch to enhanced mode. The evaporation setpoint is temporarily adjusted to the preset evaporation capacity curve during the first 30 minutes of the intermittent period. After 30 minutes, the evaporation setpoint will exponentially decay back to the preset curve. ; in The time constant is the end time of the current rainfall period. correspond The enhanced decay period, in units of .

[0103] Throughout the testing process, the three-layer hierarchical fuzzy rule base was used as... The control cycle executes the above reasoning and instruction generation repeatedly. The real-time status of the circulating water bath jacket, infrared radiation plate, fan, and radiant heat source is synchronously recorded by the integrated control and data acquisition module. Combined with the water supply flow data and sensor array data from step 3, a complete multi-physical quantity time series dataset is formed for the inversion required in step 5. .

[0104] In step 5, the pre-trained deep operator network proxy model is fine-tuned using real-time measured water flow rate, internal volumetric water content and temperature of the sample, and environmental actuator status data. A Bayesian approximate inversion method is then used to obtain the unsaturated permeability coefficient function of the soil sample and its uncertainty quantification results. Specifically: The inputs for step 5 include: the physical information neural network model pre-trained in step 2 and the initial values ​​of the trainable hydraulic parameters identified internally. Step 3: Measured water supply flow sequence recorded by full-time closed-loop tracking control. Step 4: The spatiotemporal sequence dataset acquired through multiphysics collaborative loading includes the sample interior. depth measurement points Volumetric water content at ,temperature Evaporation rate sequence obtained by converting the mass loss of distilled water from the top weighing sensor of the sample. Circulating water bath jacket set temperature sequence Infrared radiating plate power sequence Fan speed sequence and radiant heat source power sequence ; and the soil parameter vector used in step 2 and environmental parameter vector The output is the unsaturated permeability coefficient function of the tested soil sample. Optimization parameter vector Together with the 95% confidence interval of each component of the parameter, and the uncertain heatmap data of the full-time and space-time predicted water content and temperature field.

[0105] Step 5 is achieved through the following sub-steps.

[0106] Step 501: Reconstruction and Spatiotemporal Alignment of Multiphysics Observation Dataset The time series of physical quantities collected in steps 3 and 4 are arranged at a uniform sampling interval. Downsampling and synchronization alignment are performed to eliminate clock skew between data acquisition modules of each sensor. For each time step... Construct observation vectors ; Finally, the observation matrix is ​​formed. ,in This represents the inversion time step. Simultaneously, the input sequences of each environmental actuator in step 4 are also downsampled and recorded as a control input matrix. The four columns are as follows: , , and .

[0107] Step 502: Construction of a water-thermal coupled proxy model based on cross-attention deep operator networks Construct a surrogate model based on a deep operator network. The input to the surrogate model is a vector of hydraulic parameters. and spacetime coordinates The output is the volumetric water content predicted by the coupled model under these parameters. and temperature .

[0108] The deep operator network consists of a parameter encoding branch network, a spatiotemporal coordinate backbone network, an environment context encoder, and a cross-attention fusion module. The parameter encoding branch receives the normalized hydraulic parameter vector. The parameter latent representation is obtained after passing through a three-layer multilayer perceptron. The structure of a multilayer perceptron is: a fully connected layer. Swish activation; fully connected layer Swish activation; fully connected layer Linear output.

[0109] The environment context encoder receives a downsampled sequence of environment actuators. It is fed into a two-layer gated loop unit with 64 hidden units, and outputs an environment context vector. After linear projection to 128 dimensions, and with Element-by-element addition constitutes the joint latent representation of parameters and environment. ,in is a trainable projection matrix.

[0110] The spatiotemporal coordinate backbone network receives normalized spatial coordinates. and time coordinates First of all, and Apply Fourier feature maps respectively: ; ; frequency and From Gaussian distribution The spatiotemporal feature vectors are sampled and fixed. The concatenated spatiotemporal feature vectors have a dimension of 192 and are fed into the backbone network. The backbone network is a 6-layer fully connected residual network with 256 hidden units per layer, using Swish activation. Skip connections are added between every two layers to output the spatiotemporal features. .

[0111] The cross-attention fusion module combines parameter-environment joint latent representations. As a query, spatiotemporal features As keys and values, multi-head cross-attention is computed: ; in , , , The number of heads is 4. The cross-attention output is then processed by a feedforward network and a residual connection before being connected to... The layers are spliced ​​together and fed into two parallel output heads: the moisture content output head is a fully connected layer. Using Softplus activation ensures that the output is non-negative and limited to the physical range. Internally; the temperature output head is a fully connected layer. Add an offset after linear output The final proxy model output and .

[0112] Step 503: Generation of pre-training data for the surrogate model and multi-task pre-training Using the physical information neural network pre-trained in step 2 and the water-thermal coupled dynamics model constructed in step 4, a simulation dataset covering a wide parameter space is generated. In the parameter vector... Within the 6-dimensional space, the identification value in step 2 Centered on, Latin hypercube sampling is used to generate Group parameter samples, with the sampling range being the initial values ​​of each parameter. Within a physically reasonable range, the parameters are truncated. For each set of parameter samples, the water-thermal coupling control equations are solved using the finite element method, with boundary conditions consistent with those in step 4, and a time step of [missing information]. The space is discretized into 128 nodes, and the simulation duration is the total experiment duration. Randomly selected from simulation results a spacetime point Record the corresponding and The labels are used to form the pre-training dataset.

[0113] The pre-training loss function is ; in As the weight of the temperature term, For gradient consistency regularization: ; , This represents 10% of the spatiotemporal points sampled randomly. Training uses the AdamW optimizer with an initial learning rate of... Weight decay Cosine annealing scheduling, 200 training epochs, batch size 4096. After training, the surrogate model weights are saved. .

[0114] Step 504: Design of Transfer Learning Fine-tuning Strategy Proxy model parameter set The dataset is divided into three groups: a spatiotemporal feature encoding group, containing Fourier mapping frequency coefficients and backbone network residual block weights; a parametric encoding group, containing branch network multilayer perceptron weights; and a cross-attention and output head group, containing attention weight matrices and output head weights. The spatiotemporal feature encoding group is responsible for capturing common physical response patterns among soil samples, while the parametric encoding group and attention group are related to specific soil properties. During transfer learning fine-tuning, the spatiotemporal feature encoding group is completely frozen; the third fully connected layer in the parametric encoding group is unfrozen, and the learning rate is... The cross-attention group and output head group are all unfrozen, and the learning rate is [missing information]. The environment context encoder is fully unfrozen, with a learning rate of [missing information]. It is trained together with the thawed parameter set.

[0115] Step 505: Uncertainty-aware inversion loss function based on Bayesian approximation A Monte Carlo dropout layer is connected to the end of the output header of the proxy model, with dropout probability... Monte Carlo dropout maintains activation during both training and inference—it is used as a regular stochastic regularization method during training, and during inference, the output mean and variance are obtained through multiple forward propagations. The stochastic masking effect of dropout is used to construct a Bayesian approximate inference, thereby quantifying prediction uncertainty without adding additional model parameters.

[0116] Define the inversion loss function It consists of three parts.

[0117] The first part consists of multiphysics data fitting terms. For moisture content and temperature observations, a negative log-likelihood loss is used, assuming that the observation noise follows a zero-mean Gaussian distribution, and the model predicts the mean. and Prediction variance and Sample variance estimation using Monte Carlo dropout with multiple forward propagations: ; The second part is the physical constraint residual term. Random sampling within the solution domain. For each internal configuration point, calculate the mean square residuals of the Richards equation and the heat transport equation. The residual formula is the same as in step 4. ; The third part is the parameter prior regularization term. Construct a Gaussian prior using the initial parameter values ​​identified in step 2: ; in This is a diagonal prior covariance matrix, with diagonal elements set to the initial values ​​of each parameter. set up.

[0118] The total inversion loss function is ; , .

[0119] Step 506: Two-level optimization and posterior distribution estimation of parameters The parameter inversion employs a two-layer optimization strategy, with the outer layer optimizing the parameter vector. The inner layer estimates the variance after adjusting the Monte Carlo dropout layer based on the current parameters.

[0120] The outer layer uses an adaptive covariance matrix evolution strategy for global optimization, with a population size of [missing information]. Initial step size In each round of evolution, the parameter vector of each individual... Feed into the proxy model and propagate forward. Next, Monte Carlo dropout is enabled to obtain the predicted mean and variance and calculate... The average of the 50 losses is taken as the individual fitness. After 100 rounds of evolutionary iteration, the 10 individuals with the lowest fitness are selected as candidate solutions.

[0121] Initiate local refinement near candidate solutions, switching to a stochastic gradient Langevin dynamic sampler: ; Step length .run Step, abandon the past Preheating period, use the remaining The posterior mean of the estimated parameters for each sample The sum and covariance matrix. The parameters... The confidence interval is determined by the marginal distribution of its... and Quantiles are determined.

[0122] Step 507: Permeability Coefficient Function Reconstruction and Inversion Data Output Based on the posterior mean obtained from the inversion Reconstructing the unsaturated permeability coefficient function : ; in Simultaneously, it generates the uncertainty envelope for permeability coefficient prediction: using each of the posterior sample sets... Calculate separately Curves, take each The value of the curve Confidence intervals serve as a visual output of the inversion uncertainty.

[0123] The final output of step 5 includes: the optimal parameter vector. and the confidence intervals of its components, permeability coefficient function The curves and uncertainty envelope, the mean field and variance field of the full-time and spacetime moisture content and temperature predictions. These outputs serve as the core data source for the test report in step 6. The information will be entered directly into the soil sample hydraulic properties table, and the uncertainty envelope diagram will be attached to the quality control section of the report.

[0124] In step 6, the process data from all steps are summarized, and a multi-dimensional test report is automatically generated, including hydrological rhythm adaptability, water supply execution accuracy, environmental collaborative control quality, and permeability coefficient inversion reliability. Specifically: The input for step 6 is all the core outputs of the preceding five steps, including: the target water supply rhythm curve obtained in step 1. Retrieve records from the digital hydrological rhythm curve benchmark database; Step 2: Convergence curve of physical information neural network training, and identification of initial hydraulic parameters. and executable target water supply curve Correction records; Step 3: Measured flow sequence for full-time closed-loop flow tracking. Normalized root mean square error And control mode switching logs; Step 4: Full-time and spacetime temperature, water content, and evaporation rate sequences and environmental actuator status records obtained from multi-physics field collaborative loading; Step 5: Optimal hydraulic parameter vector obtained from inversion. Its 95% confidence interval and permeability coefficient function curve The system includes an uncertainty envelope, and predicts the mean and variance fields of water content and temperature across all time and space. The output is a structured, multi-dimensional test report, comprising hydrological rhythm adaptability analysis, water supply execution accuracy assessment, environmental collaborative control quality, soil sample hydraulic parameter inversion results and reliability evaluation, and comprehensive test conclusions and data archiving index.

[0125] Step 6 is achieved through the following sub-steps: Step 601: Overall Report Architecture and Data Aggregation The report is logically organized into five layers, from input rhythm to execution and loading, and then to inversion analysis. The first layer contains basic test information, including soil sample number and dry density. Organic matter content Initial volumetric moisture content Preset ambient temperature and preset evaporation capacity The first layer contains the test start and end timestamps. The second layer is the hydrological rhythm adaptability analysis. The third layer is the water supply flow tracking control accuracy assessment. The fourth layer is the environmental simulation collaborative loading quality assessment. The fifth layer is the permeability coefficient inversion results and uncertainty quantification. Each layer extracts key features from the original data and generates standardized charts and tables. The charts are embedded in vector format, and the tables are labeled in plain text.

[0126] Step 602: Hydrological Rhythm Adaptability Analysis Starting from the target water supply rhythm curve in step 1, calculate its similarity to similar rainfall pattern curves in the benchmark reservoir. Extract... The four characteristic values: total water supply volume Peak flow Flow variation coefficient and percentage of zero-traffic time The system searches the benchmark database for the five records with the smallest Euclidean distance to these four eigenvalues, calculates the average distance, and determines whether the target curve is located in a dense region of the benchmark database. If the average distance is below a threshold, it is marked as a typical rhythm; otherwise, it is marked as an atypical rhythm, and a warning is given in the report, indicating that the extrapolation of subsequent inversion results should be approached with caution.

[0127] Step 603: Evaluation of the accuracy of water supply flow tracking control Derive the normalized root mean square error from the closed-loop tracking control results in step 3. As an overall accuracy indicator, the total test time was divided into segments according to the six macroscopic working conditions identified in step 4, and the statistics for each segment were calculated separately. Maximum absolute deviation And the overshoot percentage. Plot a comparison chart of the target flow rate curve and the actual flow rate curve for the entire time period, marking different operating conditions with different background colors. Simultaneously, plot a histogram of the tracking error probability density distribution and use the Shapiro-Wilke test to check the normality of the error. Record the online update parameters of the gated cyclic unit regulator in step 3. , , The trajectory was analyzed, and its variation at the moment of operating condition switching was examined. Finally, the report provided the water supply execution accuracy level, when the global... And each working condition section The rating is A if it is applied at the time, otherwise it will be downgraded step by step.

[0128] Step 604: Environmental Simulation Co-loading Quality Assessment Based on the full-time and space-time sensor data from step 4, the evaluation is conducted from three aspects: temperature field uniformity, evaporation rate tracking, and consistency of coupling between different physical fields. For the temperature field, the data from the middle of the sample is calculated. The mean square error between the measured temperature and the set temperature is calculated, and the gradient analysis of temperature measurement points at both ends of the sample is used to determine whether there is abnormal heat accumulation or uneven heat dissipation. For the evaporation rate, a comparison graph of the preset evaporation capacity curve and the measured evaporation rate curve is plotted, and the lag time and transient error peak value are statistically analyzed in segments. The coupling consistency assessment uses the water-heat-vapor coupling model in step 4 as a reference, and the measured environmental actuator input sequence is used. Input the model, and normalize the theoretical temperature field and water content field obtained by the solution with the measured values ​​point-to-point using the coefficients of determination. Calculate. If If the co-loading quality is good, it is determined that the quality is excellent; otherwise, the source of the deviation is traced back to the subsystem analysis, such as fan efficiency decline or changes in jacket contact thermal resistance, and recorded in the anomaly and maintenance recommendation section of the report. Simultaneously, the deviation between the actual and theoretical values ​​of the enhancement factor during the intermittent evaporation enhancement period in step 4 is statistically analyzed.

[0129] Step 605: Permeability Coefficient Inversion Results and Reliability Assessment Directly referencing step 5 The residual volumetric water content and its 95% confidence interval are listed in tabular form. saturated volumetric water content Intake value reciprocal Aperture distribution parameters saturated permeability coefficient and curvature factor The initial identification values ​​from step 2 are appended to the table for comparison, and the relative deviation percentage is calculated. The curve, along with the uncertainty envelope generated in step 5, is plotted as a semi-logarithmic graph, and key saturation values ​​are analyzed. , , The permeability coefficient at a given location provides a 95% confidence interval. Referring to step 5, the spatial distribution of the relative variance of the full-time and spatial predictions is displayed using contour maps to assess the reliability of the inversion in different spatiotemporal regions. If the prediction variance of a certain region exceeds twice the global average, it is marked as a high-uncertainty area, and the possible cause is noted in the report analysis as insufficient sensor coverage or excessive nonlinearity of the physical processes in that region.

[0130] Step 606: Overall Conclusion and Data Archiving Index The system automatically generates comprehensive textual conclusions, including a test validity statement, an overview of the soil sample's hydraulic properties, and an inversion reliability rating. The test validity statement determines whether the water supply and loading errors are within the allowable tolerance zone. The overview of the soil sample's hydraulic properties provides... The report details the magnitude and typical unsaturated permeability coefficient variation range. The reliability rating is based on the confidence interval width and the proportion of high uncertainty areas, resulting in three levels: reliable, generally reliable, or requiring repeated testing. A test data archive index table is then generated, listing the storage paths and verification hash values ​​of the original data files for each item. The report concludes with an input / output summary diagram of the five-step main process, illustrating the data flow relationships between modules in a matrix diagram format.

[0131] The present invention also provides several embodiments for detailed description of the above-described apparatus and method: Example 1: Permeability test of silty clay samples under typical rainstorm-evaporation alternation rhythm This embodiment uses silty clay taken from a watershed in South China as the test object, and performs the entire process test according to the method and apparatus of this invention.

[0132] The basic physical properties of the sample are: dry density Organic matter content Initial volumetric water content Porosity Sample height Inner diameter of the sample cylinder .

[0133] In step 1, starting from hourly hydrological data of a continuous hydrological year from the target sampling area, a digital hydrological rhythm curve benchmark library containing 365 daily segments was constructed after empirical wavelet transform decomposition and reconstruction using a bidirectional long short-term memory variational autoencoder. During the retrieval, a representative continuous 72-hour period was selected from this area, containing a single intensity... The torrential rain event and the subsequent 48-hour evaporation and drying process. The extracted target water supply rhythm curve. The total water supply volume is Peak traffic .

[0134] In step 2, soil parameters and the target water supply rhythm curve are input into a pre-trained multi-scale hierarchical Fourier feature physical information neural network. The network converges after three stages of training (10,000 steps in the first stage, 50,000 steps in the second stage, and 30,000 steps in the third stage), with the total training loss decreasing from the initial value. Descending to Digital twin predictions show that at the peak of the rainstorm... At this location, the predicted saturation value at the bottom of the sample is The value did not exceed the 0.95 threshold; the predicted maximum infiltration flow rate was... The flow rate did not exceed the upper limit of the metering range. Therefore, there is no need to trigger curve correction; output can be done directly. Meanwhile, the initial values ​​of the hydraulic parameters identified by the physical information neural network are: , , , , , .

[0135] In step 3, the water supply flow generation and closed-loop control unit are activated. The short-pulse response experiment identifies the liquid capacity in the pipeline. and pipeline hydraulic resistance Prediction time domain of rolling time-domain feedforward optimization Mini-batch gradient descent iterates 5 times per step. The initial shape parameters of the adaptive non-singular terminal sliding mode feedback controller are set to... , , The gated cyclic unit learner dynamically adjusts the shape parameters with a 200ms period. After a full-time 72-hour test, the calculated full-time normalized root mean square error is... Meets design specifications The statistical results for each working condition are shown in Table 1.

[0136] In step 4, the three-layer hierarchical fuzzy rule base runs with a control cycle of 0.5 seconds. The activation rules of the 343 fuzzy rules in the middle layer are mainly concentrated in the NB (negative large) and NM (negative medium) subsets during the rainstorm infiltration period, and concentrated in the PS (positive small) and PM (positive medium) subsets during the post-rain evaporation period. Rule consequent correction based on the physical information neural network is executed once every 10 control cycles. A total of [number missing] rules were executed throughout the test. 1 rule (representing the total number of active rules) The value was corrected, with the average correction magnitude being equal to the initial value. During the 30-minute intermittent evaporation enhancement period, the measured enhancement factor of the evaporation rate was [missing value]. , compared with theoretical value The deviation is .

[0137] In step 5, with The size of the downsampled observation matrix is (5 depth measurement points). After 200 rounds of pre-training, the root mean square error of water content prediction on the validation set is: The root mean square error of temperature prediction is In the transfer learning fine-tuning phase, the spatiotemporal feature encoding group was frozen, while the third layer of the parameter encoding group, the cross-attention group, the output head group, and the context encoder were unfrozen. After 50 rounds of fine-tuning, the root mean square error of water content prediction was reduced to [missing value]. After 100 rounds of the adaptive covariance matrix evolution strategy, the stochastic gradient Langevin dynamic sampler was run for 5000 steps. The posterior mean and 95% confidence intervals of the hydraulic parameters are shown in Table 2. Permeability coefficient function. The curve at saturation The width of the 95% confidence interval at is This indicates that the inversion results have a high degree of reliability.

[0138] The automatically generated test report in step 6 shows: the water supply execution accuracy level is A; the coupling consistency of environmental co-loading quality is... The reliability rating is reliable.

[0139] Example 2: Permeability test of sandy loam soil samples under long-term intermittent light rain rhythm This embodiment uses sandy loam soil taken from a test site in North China as the test object. The basic physical properties of the sample are: dry density Organic matter content Initial volumetric water content Porosity The sample size was the same as in Example 1.

[0140] In step 1, a representative 120-hour period from a continuous hydrological year in this region is selected. This period is mainly characterized by intermittent light rain and includes 7 rainfall events, none of which have an intensity exceeding [a certain value]. The intervals between rainfalls last from 8 to 24 hours. The coefficient of variation in flow rate of the target water supply rhythm curve. Percentage of time with zero traffic It belongs to the typical hydrological rhythm of arid-semi-arid regions.

[0141] In step 5, the physical information neural network identifies the initial values ​​of the hydraulic parameters as follows: , , , , , Due to the high saturation permeability coefficient of sandy loam, digital twin predictions showed that the saturation at the bottom of the sample was only [value missing] at the peak of rainfall. It is far from reaching the premature saturation threshold, so no curve correction is needed.

[0142] In step 3, the gated recurrent unit learning regulator exhibited significant shape parameter adjustments at the start time of all seven rainfall events: The parameters averaged at the onset of rainfall from Upgraded to , Parameters are from Upgraded to , Parameters are from Downgraded to This indicates that the regulator automatically increased the aggressiveness of the controller's response during the transient phase of rainfall initiation. (Full-time) The results for each working condition are shown in Table 1.

[0143] In step 4, the three-layer hierarchical fuzzy rule base frequently switches between long-term drought periods and intermittent recovery periods. The output interval width of the interval type II fuzzy set increases on average at the time of operating condition switching. The system automatically reduced the control gain accordingly. No warning flags were triggered due to excessive uncertainty during the entire 120-hour test.

[0144] In step 5, the inversion result after fine-tuning through transfer learning is as follows: , , , , , The relative deviations from the initial identification values ​​in step two are respectively... , , , , and None exceeded This indicates that the initial identification in step two has achieved high accuracy, and the transfer learning inversion in step five further refines the parameter estimation and provides uncertainty quantification. The curve at saturation The width of the 95% confidence interval at is .

[0145] Step 6 Test Report Rating: Water supply execution accuracy level is A; Environmental co-loading quality The reliability rating is reliable.

[0146] Table 1. Statistics on water supply flow tracking accuracy of two soil samples under different macroscopic conditions.

[0147] In Table 1, the sandy loam soil examples did not include the heavy rain infiltration period, indicated by "—". Table 1 presents the statistical breakdown of water supply flow tracking accuracy for the two representative soil samples under six macroscopic conditions. As shown in the table, the NRMSE was highest during the heavy rain infiltration period due to the most drastic flow changes. m Both the overshoot percentage and accuracy were the highest; the highest accuracy was observed during prolonged drought periods due to the gradual change in flow rate. NRMSE for both soil samples under various operating conditions. m The NRMSE values ​​were all below 2%, and the NRMSE values ​​were all below 1% for all time periods. The water supply execution accuracy level was A, indicating that the multi-mode closed-loop flow tracking control method of the present invention has good adaptability to different soil types and different hydrological rhythms.

[0148] Table 2 Comparison of hydraulic parameter inversion results for silty clay samples

[0149] Table 2 presents a complete comparison of the hydraulic parameters of silty clay samples from the initial identification in step two to the Bayesian inversion in step five. As shown in the table, the relative deviations between the posterior mean of the inversion in step five and the initial identification values ​​in step two are all within [a certain range]. Within this range, it is evident that the physical information neural network has demonstrated good parameter identification capabilities during training. Step five further provides the 95% confidence intervals for each parameter, including the saturation permeability coefficient. The confidence interval width is Approximately 1 / 3 of the posterior mean This quantifies the uncertainty level of the inversion results, providing a reliable reference for engineering applications. Curvature factor It remains unchanged during inversion because it is less sensitive to the objective function, and the posterior distribution degenerates into the prior distribution.

[0150] Figure 2 The unsaturated permeability coefficient functions of two samples, silty clay and sandy loam, are shown. The curves and their 95% confidence interval envelopes are shown in the figure. It can be seen from the figure that the saturated permeability coefficient of sandy loam is about an order of magnitude higher than that of silty clay, and the slope of the unsaturated segment of its curve is steeper. This is consistent with the physical characteristics of sandy loam, where large pores account for a higher proportion and water is rapidly expelled when near saturation. The confidence interval envelope of silty clay is relatively narrow, especially at moderate saturation (…). to The presence of this region indicates a high degree of reliability in the soil sample inversion results.

[0151] Figure 3 To compare the target and measured flow rates of the silty clay samples over the entire time period, different background colors are used to distinguish different macroscopic operating conditions. As shown in the figure, the measured flow rate curve closely tracks the target curve even during the rapid flow changes in the infiltration period of heavy rainfall (0 to 16 hours), with the maximum deviation occurring at the peak of the heavy rainfall. Nearby, the absolute tracking error is . Figure 4 The absolute tracking error distribution over all time periods has an average error of [missing information]. More than 99.7% of the error data points fall within the upper limit of 3σ, indicating that the error distribution is stable and there is no systematic shift.

[0152] Embodiments of the present invention may be provided as methods, systems, or computer program products. Therefore, the present invention may take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0153] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0154] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0155] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0156] Contents not described in detail in this specification are prior art known to those skilled in the art. It is hereby indicated that the above description is intended to help those skilled in the art understand this invention, but does not limit the scope of protection of this invention. Any equivalent substitutions, modifications, improvements, or simplifications of the above descriptions that do not depart from the essential content of this invention fall within the scope of protection of this invention.

Claims

1. A permeability coefficient testing device, characterized in that, include: The sample loading and sealing unit is used to contain the sample and provide a one-dimensional vertical permeation path and temperature control boundary. A water supply flow generation and closed-loop control unit is used to provide a controllable flow rate of water to the bottom of the sample; The environmental simulation loading unit is used to apply radiation, airflow, and temperature environments to the sample. The multi-parameter sensing and measurement unit includes a moisture content-temperature composite sensor array, a sample top weighing sensor, an air temperature and humidity sensor, and a pipeline pressure sensor. Piping switching and drainage collection unit; The integrated control and data acquisition unit includes an industrial computer and a high-speed data acquisition board and a programmable logic controller connected thereto. The industrial computer is configured with: The hydrological rhythm digitization extraction module is used to perform empirical wavelet transform eigenmode decomposition and variational autoencoder reconstruction, build a digital hydrological rhythm curve benchmark library, and extract the target water supply rhythm curve from it; The physical information neural network pre-modeling module is used to generate an executable target water supply curve using a multi-scale hierarchical Fourier feature physical information neural network, and to identify the initial values ​​of soil sample hydraulic parameters during the training process. The multi-mode closed-loop flow tracking control module is used to generate composite control commands based on the executable target water supply curve and the measured flow rate fed back by the Coriolis mass flow meter, through rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so as to drive the variable frequency booster pump group and the electric proportional regulating valve. The multi-physics collaborative environment loading module is used to generate four control commands for the environment simulation loading unit by using a three-layer hierarchical fuzzy rule base to decouple and coordinate multiple objectives based on the evaporation capacity setting and temperature setting extracted from the digital hydrological rhythm curve. The transfer learning permeability coefficient inversion module is used to obtain the unsaturated permeability coefficient function and its confidence interval of the soil sample by using measured flow and measurement data from multi-parameter sensing and measurement units, based on a pre-trained cross-attention deep operator network proxy model for transfer learning and Bayesian uncertainty inversion. The multidimensional test report generation module is used to summarize the process data of the above modules and generate a test report that includes hydrological rhythm adaptability analysis, water supply execution accuracy assessment, environmental collaborative control quality assessment, and permeability coefficient inversion reliability assessment.

2. A testing method for a permeability testing device, applied to the permeability testing device according to claim 1, characterized in that, include: Step 1: Obtain hydrological data of the target sampling area, construct a digital hydrological rhythm curve benchmark library through empirical wavelet transform and variational autoencoder, and extract the target water supply rhythm curve from it; Step 2: Based on the physical properties of the soil sample, a multi-scale hierarchical Fourier feature physical information neural network is used for course learning and training to generate an executable target water supply curve. At the same time, the initial values ​​of the hydraulic parameters of the soil sample are identified during the training process. Step 3: Based on the executable target water supply curve and the real-time measured water supply flow, perform multi-mode closed-loop flow tracking control, including rolling time-domain feedforward optimization and adaptive non-singular terminal sliding mode feedback control, so that the flow entering the bottom of the sample tracks the executable target water supply curve. Step 4: Based on the evaporation capacity and temperature settings extracted from the digital hydrological rhythm curve, and the real-time measured sample moisture content and temperature, multi-physics field collaborative loading control is carried out through a three-layer hierarchical fuzzy rule base to adjust the environmental simulation loading unit to achieve collaborative simulation of temperature, evaporation and water supply. Step 5: Using the measured water supply flow rate, internal volumetric water content and temperature of the sample, and environmental actuator status data throughout the entire time period, the pre-trained deep operator network proxy model is fine-tuned through transfer learning, and the unsaturated permeability coefficient function of the soil sample and its uncertainty quantification results are obtained by using the Bayesian approximate inversion method. Step 6: Summarize the process data from all steps and automatically generate a multi-dimensional test report that includes hydrological rhythm adaptability, water supply execution accuracy, environmental collaborative control quality, and permeability coefficient inversion reliability.

3. The testing method of the permeability testing device according to claim 2, characterized in that, In step 1, hydrological data of the target sampling area is acquired, and a digital hydrological rhythm curve benchmark library is constructed using empirical wavelet transform and variational autoencoder. The target water supply rhythm curve is then extracted from this library. Specifically: After normalizing the sequences of rainfall intensity, evaporation and surface runoff depth, empirical wavelet transform is performed on each sequence to extract the two intrinsic mode functions with the highest energy proportion. These functions are then input into a bidirectional long short-term memory variational autoencoder for reconstruction. The reconstructed modes are then added together to obtain the hydrological rhythm fundamental wave. The residuals between the original sequence and the fundamental wave are used to generate random perturbation terms using a second-order Markov chain. The fundamental wave and the perturbation terms are added together to form a digital hydrological rhythm curve and stored in the reference library.

4. The testing method of the permeability testing device according to claim 3, characterized in that, In step 2, the course learning and training is divided into three stages: The first phase uses constant flow boundary data for training; The second phase uses variable flow boundary data for training; The third stage uses the rainfall component in the digitized hydrological rhythm curve as unlabeled boundary data for unsupervised training. The training loss function includes a data fitting term, a physical residual term, a boundary condition term, and an initial condition term, and uses a residual adaptive weighting strategy to dynamically adjust the weights of each loss term.

5. The testing method of the permeability testing device according to claim 4, characterized in that, In step 3, the multi-mode closed-loop flow tracking control specifically includes: Before each test, the liquid capacity and liquid resistance of the water supply pipeline were identified online through short impulse response experiments, and a lumped parameter hydraulic transmission model was established. In each control cycle, a lumped parameter hydraulic transfer model is connected in series with a physical information neural network to form a differentiable mapping, and a feedforward control sequence is solved to minimize the control increment penalty while ensuring that the predicted flow in the future multi-step flow tracks the reference trajectory. A non-singular terminal sliding surface is constructed using the flow tracking error, and combined with the disturbance estimation based on one-step delay observation, a feedback compensation quantity containing equivalent control and super-spiking switching control is generated. The sum of feedforward and feedback is decomposed into low-frequency and high-frequency components according to the frequency band, which are used to adjust the speed of the variable frequency booster pump and the opening of the electric proportional regulating valve, respectively.

6. The testing method of the permeability testing device according to claim 5, characterized in that, In step 4, the three-layer hierarchical fuzzy rule base includes: The upper-level operating condition identification layer divides the system operating status into six macro-operating conditions based on the current water supply flow, moisture content and environmental settings: heavy rain infiltration period, moderate rainfall infiltration period, light rain wetting period, post-rain evaporation period, long-term drought period and intermittent recovery period. The middle-layer multi-objective decoupling coordination layer takes temperature deviation, moisture content deviation and evaporation rate deviation as inputs and outputs temperature correction, moisture content correction and evaporation capacity correction. Its fuzzy rule consequent is corrected online by the physical information neural network based on the prediction deviation. The underlying actuator allocation layer maps the target correction amount to the execution instructions of the circulating water bath jacket, infrared radiation plate, fan and radiant heat source through a preset allocation matrix; The membership function of the fuzzy rule base adopts an asymmetric Gaussian type and is wrapped with an interval type II fuzzy set to handle sensor noise and model uncertainty.

7. The testing method of the permeability testing device according to claim 6, characterized in that, The online gradient correction of the fuzzy rule consequent of the middle-layer multi-objective decoupling coordination layer is specifically as follows: on the set of active rules where the excitation intensity exceeds a preset threshold, a loss function is constructed based on the deviation between the moisture content and temperature of the next control cycle predicted by the physical information neural network and the set value. The gradient of this loss with respect to the rule conclusion value is calculated through backpropagation, and the rule conclusion value is updated along the gradient descent direction.

8. The testing method of the permeability testing device according to claim 7, characterized in that, In step 5, the deep operator network proxy model consists of a parameter encoding branch, an environment context encoder, a spatiotemporal coordinate backbone network, and a cross-attention fusion module. During the transfer learning fine-tuning, the residual block weights and Fourier mapping frequency coefficients of the spatiotemporal coordinate backbone network are frozen, and some fully connected layers of the parameter encoding branch, the cross-attention module, the output head, and the environment context encoder are trained with different learning rates.

Images