Soil hydraulic parameter robust calculation method based on one-dimensional upward infiltration test

Abstract

The invention discloses a soil hydraulic parameter robust calculation method based on a one-dimensional upward infiltration test. According to the method, a piston type moisture profile model (PASCP) is constructed, and the application range of the model is effectively expanded through an innovative method for converting heterogeneous moisture content distribution of an unsaturated region into an equivalent homogeneous region. Based on the model, an analytic infiltration model considering saturation time is established, and an interval constraint optimization objective function is introduced, so that the inversion precision and robustness of the key hydraulic parameters of the Brooks-Corey model are remarkably improved. Experiments prove that the method shows excellent applicability under different soil texture conditions, and particularly shows higher result reliability under the condition that measurement errors exist. According to the parameter estimation method based on the standard cutting ring, an efficient and practical technical means is provided for rapid and accurate representation of the hydraulic characteristics of the soil, and the method has wide application prospects in the fields of agricultural irrigation, hydrological simulation and the like.

Description

Technical Field

[0001] The invention relates to a robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test, and belongs to the technical field of soil physics and soil hydrology. Background Art

[0002] Soil hydraulic properties, including the soil moisture characteristic curve θ(h) and the unsaturated hydraulic conductivity curve K(θ), are key input parameters for hydrological models based on the Richards equation. Together, they determine the transport characteristics of soil water in the unsaturated zone. Accurately measuring these soil hydraulic parameters is crucial not only for field-scale water and solute transport studies but also for global-scale hydrological and energy cycle simulations. Traditional laboratory methods for determining θ(h) and K(θ) are the pressure film method and the transient profile method, respectively. While theoretically simple, the practical procedures are cumbersome and time-consuming. To improve efficiency, numerical inversion techniques for simultaneously estimating soil θ(h) and K(θ) based on transient flow experiments have rapidly developed, including single- and multi-step drainage experiments, evaporation methods, and upward infiltration methods. Among these, upward infiltration inversion techniques based on ring-cut soil sampling have achieved significant breakthroughs in the past decade. This method utilizes capillary-driven infiltration against gravity to slow water movement, making it particularly suitable for accurately measuring water flux through 5-cm-high ring cutters, thus providing an effective path for large-scale soil hydraulic investigations. However, existing numerical inversion techniques still face a fundamental challenge: the lack of a universal and efficient solution algorithm for the highly nonlinear Richards equation. Current numerical schemes often require customized spatial discretization strategies based on different soil types and boundary conditions. This strong scenario-dependence severely restricts their wider practical application.

[0003] In order to address the universality and complexity issues faced by the above numerical methods, researchers have gradually turned their attention to analytical and semi-analytical physical models in order to more efficiently predict soil hydraulic parameters. Physical models have also been applied to parameter inversion under upward infiltration conditions. Moret-Fernandez and Latorre (2016) combined the saturated hydraulic conductivity K measured based on Darcy's law with the Haverkamp model. s , successfully inverted the infiltration rate S and shape parameter β, and further estimated θ(h). Chinese patent publication number CN114397427A discloses a soil hydraulic property prediction method based on the infiltration process of ring cutter soil samples, which provides a new algorithm to simultaneously estimate θ(h) and the infiltration rate S. However, both of the above methods require additional measurement of the saturated hydraulic conductivity K of the soil. sAs the input parameter, θ(h) is estimated. As pointed out by van Genuchten (1991), the originator of the classic soil hydraulics model (van Genuchten model), K s Parameters are often difficult to measure accurately. To address this issue, Wu et al. published a paper in the Journal of Hydrology in 2022, based on patent CN114397427A, and proposed a constant pressure percolation physical model (ASCP). This model can simultaneously estimate the key parameters describing θ(h) and K(θ), including the pore distribution index n, the inverse of the intake suction 1 / h, and the relative humidity. d and saturated hydraulic conductivity K s However, the accurate estimation of parameters by the ASCP method relies on strict equality constraints, that is, it is necessary to accurately measure the single-point cumulative infiltration (I * ) and the corresponding soil column height (L * However, in actual laboratory measurements, due to the limited accuracy of the pressure sensor and the interference of bubbles in the Malvern flask on the indoor air pressure (Latorre et al., 2015), the I * It is often difficult to measure accurately, which leads to the failure to meet strict equality constraints, thereby weakening the performance of the ASCP method in real scenarios. Therefore, how to reduce the measurement error is crucial to the I * The impact of infiltration and the robustness of the method have become key issues that need to be addressed urgently. Furthermore, experimental evidence from Moret-Fernandez and Latorre (2019) indicates that after the wetting front reaches the soil surface, the soil water absorption process will continue slowly until it reaches a final stable state. Given the significant differences in the time required for soils of different textures to reach a stable state, infiltration data at this stage may contain more critical information about soil hydraulic properties. Incorporating this information into the ASCP method may further improve the accuracy of hydraulic parameter estimates. Summary of the Invention

[0004] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test, so as to further improve the accuracy of hydraulic parameter estimation.

[0005] The technical solution adopted in the present invention is:

[0006] A robust calculation method for soil hydraulic parameters based on a one-dimensional upward infiltration test comprises the following steps:

[0007] Step A: Apply a constant lower boundary water head h to the ring cutter soil sample p, carry out a one-dimensional upward infiltration test on homogeneous soil, record the relationship between the cumulative infiltration volume I of the ring cut soil sample and the time t, the time t* when the wetting front first reaches the surface of the ring cut soil sample, and the cumulative infiltration volume I when the infiltration curve tends to be stable end ;

[0008] Step B: Based on the initial volumetric water content θ of the ring cut soil sample i , the ring height L* and the cumulative infiltration volume I when the soil reaches full saturation end , calculate the saturated volumetric water content of the soil θ s , the calculation formula is:

[0009]

[0010] Step C: Based on the relationship between the cumulative infiltration volume I and time t obtained in step A, determine the single-point cumulative infiltration volume I* corresponding to the time t* when the wetting front first reaches the surface of the ring cutter soil sample;

[0011] Step D: Combine the saturated volumetric water content θ of the ring cut soil sample s , initial volume water content θ i , residual volume water content θ r , the single-point cumulative infiltration volume I* when the wetting front first reaches the soil surface, and the observed data of cumulative infiltration volume changing with time I obs (t i ), the parameter optimization inversion method based on interval constraints is used to obtain the estimated parameters n, 1 / h d and K s .

[0012] Preferably, in step D, the objective function of the interval-constrained parameter optimization inversion method is constructed according to formula (14):

[0013]

[0014] Preferably, the first term of formula (14), I obs (t i ) is the measured value obtained in step A, representing the observation time point t i The cumulative infiltration volume corresponding to the time I, I sim (t i ; n,1 / h d ,K s ) is the observation time point t i The simulated cumulative infiltration volume is calculated according to the following formula:

[0015] ① Time for soil to reach saturation (t s ) is:

[0016]

[0017] Where z s * is the equivalent wetting front z fe Function expression, when the equivalent wet front z fe Reaching the soil surface (i.e. fe =L * ), the soil reaches saturation, at this time z s * It can be clearly expressed as:

[0018]

[0019] ②For infiltration time t <t s When , the functional relationship between soil cumulative infiltration I and time t is:

[0020]

[0021] ③When t≥t s When the soil reaches full saturation, the cumulative infiltration remains unchanged, and the calculation formula is as follows:

[0022] I=(θ s -θ i )L * (10);

[0023] The parameters to be estimated n, 1 / h d and K s Substitute it into the measured value I obs (t i ) and analog value I sim (t i ; n,1 / h d ,K s ) tends to be the smallest;

[0024] The second term of formula (14), I sim (L * ; n,1 / h d ,K s ) is the cumulative infiltration volume of the simulation, which is calculated according to the following formula:

[0025]

[0026] where z f =L*, σ is the penalty factor [–], ∈ is the constraint threshold [mm], and its value depends on the upper limit of the maximum measurement error in actual observation.

[0027] The remaining parameters are expressed according to the following relations:

[0028]

[0029] Preferably, formula (14) is solved using the genetic algorithm (GA package) in R language: the algorithm generates 1000 sets of estimated parameters n, 1 / h in each generation. d and K s ; When the value of φ in formula (14) does not increase within 500 consecutive generations, or the cumulative number of iterations reaches 5000, the algorithm stops running; the local search strategy is enabled; the parameter search range is K s =[10 -6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Logarithmic transformation is performed; the constraint threshold ε is set to 0.4 mm; σ is set to Here, .Machine$double.xmax represents the maximum double-precision floating-point number representable in the R language.

[0030] The present invention proposes a novel robust calculation method for soil hydraulic parameters based on one-dimensional upward infiltration test (PASCP method). Unlike the ASCP method, which can only describe the infiltration process before the wetting front reaches the soil surface, PASCP further expands the scope of application and incorporates the soil hydraulic characteristic information contained in the steady-state stage (the infiltration stage that continues after the wetting front reaches the soil surface) into the parameter estimation, thereby achieving a complete simulation of the entire upward infiltration process of the soil from the initial wet state to full saturation, and significantly improving the accuracy of hydraulic parameter estimation. In addition, in order to reduce the measurement error, the single-point cumulative infiltration amount (I * ) values, a parameter optimization method based on interval constraints is proposed based on the PASCP method to robustly estimate the soil hydraulic parameters (n, 1 / h d and K s ).

[0031] The beneficial effects of the present invention are as follows:

[0032] 1. The PASCP method effectively addresses the limited applicability of the ASCP method after the wetting front reaches the soil surface, providing an important new approach for robust estimation of soil hydraulic parameters.

[0033] 2. By introducing the interval-constrained parameter optimization inversion method, the measurement error is reduced to the single-point cumulative infiltration when the wetting front first reaches the soil surface (I * ) values, significantly improving the robustness of soil hydraulic parameter estimation;

[0034] 3. The parameter estimation method based on standard ring cutter provides an effective technical means for quickly and reliably characterizing soil hydraulic properties. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 Schematic diagram of the evolution of the water profile during soil infiltration. (a) The profile shape when the wetting front first reaches the soil surface; (b) The profile changes as the wetting front reaches the surface and the infiltration process continues; (c) The profile shape when the soil reaches a fully saturated steady state.

[0036] Figure 2 It is a calculation flow chart of the present invention.

[0037] Figure 3 This is an example of a parameter inversion experiment for a clay soil sample. (a) Schematic diagram of the one-dimensional infiltration experimental setup and measurement process for a laboratory ring-knife soil sample; (b) The measured cumulative infiltration curve and key time points, including the time when the wetting front first reaches the soil surface (t * ) and the moment when the soil reaches full saturation (t s ); (c) Iterative trajectories of PASCP and ASCP methods in the parameter estimation process; (d) Comparison of soil moisture characteristic curves inverted by PASCP and ASCP methods with the measured values.

[0038] Figure 4 This is a parameter estimation chart of five soils with significantly different textures measured in the laboratory.

[0039] Figure 5 Figure 1 shows parameter estimates for seven soils with distinct textures based on numerical simulations. (a) The PASCP method includes infiltration information during the steady-state phase; (b) The PASCP method excludes infiltration information during the steady-state phase.

[0040] The present invention will be further described below with reference to specific embodiments. DETAILED DESCRIPTION

[0041] The present invention will be further described below in conjunction with a soil infiltration test, but this should not be construed as limiting the present invention. Modifications and substitutions made to the present method, steps, or conditions without departing from the spirit and substance of the present invention are intended to fall within the scope of the present invention. Unless otherwise specified, the technical means used in the following examples are conventional means well known to those skilled in the art.

[0042] Example 1 Deduction process of PASCP method

[0043] 1. Piston-type moisture profile construction method

[0044] The soil hydraulic parameter inversion method based on the piston-type moisture profile model (PASCP) of the present invention is specifically an improvement and extension of the ASCP method previously proposed by the applicant (Wu et al., 2022, Journal of Hydrology). In the ASCP method previously proposed by the applicant, the soil moisture profile consists of a saturated zone and an unsaturated zone, where the moisture distribution in the saturated zone is uniform, while the unsaturated zone shows a significant gradient change ( Figure 1 a), its cross-section mathematical expression is:

[0045]

[0046] Where z represents the soil depth [cm], its zero point is at the bottom of the soil, and the coordinate axis is directed upward; f is the distance the wetting front advances [cm]; z s is the saturation zone length [cm]; S ei corresponds to the initial water content θ i The effective saturation [–]. n is the pore distribution index. The parameters a and b are defined as follows:

[0047]

[0048] However, the ASCP method can only describe the infiltration process before the wetting front reaches the soil surface. When the wetting front reaches the surface, the unsaturated zone gradually becomes saturated ( Figure 1 ), the moisture profile deviates from the description of the ASCP formula (Formula 1). To address this limitation, the present invention proposes a new piston-type moisture profile construction method (Modified ASCP with piston-type water profile, referred to as PASCP) based on the piston flow hypothesis of Green and Ampt (1911). Its innovation lies in avoiding the explicit simulation of the gradual saturation process of the unsaturated zone, while accurately reflecting the cumulative infiltration process of this stage, effectively expanding the scope of application of the ASCP method. Specifically, the integral method is used to transform the heterogeneous moisture content distribution of the unsaturated zone into an equivalent homogeneous zone, and the equivalent length z of the unsaturated zone is defined. ufe , which is calculated as follows:

[0049]

[0050] Among them, z ufe represents the equivalent length of the unsaturated zone, z s is the saturation region length. s With z ufe Add together to construct the equivalent piston-type wetting front depth z fe , whose expression is:

[0051] zfe =z s +z ufe (5)

[0052] After adopting the piston-shaped moisture profile model, the PASCP method can effectively simulate the upward advancement of the moist front in the form of a piston-shaped moisture profile. fe When it reaches the soil surface (e.g. Figure 1 c), the soil column enters a fully saturated state and reaches a steady state, which is defined as the soil saturation time (t s )[min].

[0053] 2. Analytical model derivation

[0054] Using the piston-shaped moisture profile mentioned above, a one-dimensional moisture infiltration analytical model suitable for constant head boundary conditions was further derived based on the principle of mass conservation.

[0055] ④ Time for soil to reach saturation (t s ) is:

[0056]

[0057] Where z s * is the equivalent wetting front z fe Function expression of the equivalent wet front z fe Reaching the soil surface (i.e. fe =L * ), the soil reaches saturation, at this time z s * It can be clearly expressed as:

[0058]

[0059] It is worth noting that t s Depends on parameters n, 1 / h d and K s , which shows its application value in parameter inversion.

[0060] ⑤For infiltration time t <t s When , the functional relationship between soil cumulative infiltration I and time t is:

[0061]

[0062] ⑥When t≥t s When the soil reaches full saturation, the cumulative infiltration remains unchanged, and the calculation formula is as follows:

[0063] I=(θ s -θ i )L* (10)

[0064] ⑦ Cumulative infiltration volume I and the distance of the wetting front z f The functional relationship is:

[0065]

[0066] The unknown parameters in the above formula are defined as follows: θ i is the initial volumetric water content of soil [cm -3 cm -3 ],θ s is the saturated volumetric water content of soil [cm -3 cm -3 ], h p Water head at the lower boundary of soil [cm], L * is the soil height [cm], and the remaining parameters are expressed according to the following relationships:

[0067]

[0068] Soil residual volume water content θ r The determination method is: the sieved soil is placed in a dry air environment with a relative humidity of less than 15% for several months until its mass remains constant. The water content of the soil at this time is the soil residual mass water content, and then multiplied by the bulk density of the ring cut soil sample to obtain the soil residual volume water content; the initial soil volume water content θ i The determination method is the drying method; the saturated volumetric water content of soil θ s The determination method is: by inverting formula (10), the expression is:

[0069]

[0070] Among them, I end It represents the stable value [cm] when the cumulative infiltration curve reaches a plateau, representing the total amount of water absorbed by the soil sample.

[0071] 3. Interval-constrained parameter optimization inversion method

[0072] The present invention further proposes an interval-constrained optimization inversion method based on the PASCP method to improve the soil hydraulic parameters (n, h d and K s ) estimation accuracy and robustness. The objective function expression of the specific algorithm is:

[0073]

[0074] Among them, the first item represents the measured data I obs (t i ) and simulated data Isim (t i ; n,1 / h d ,K s ), the root mean square error (RMSE) between the simulated value I sim (t i ; n,1 / h d ,K s ) is calculated by formula (6)-(10), where t i is the observation time point, I sim is the cumulative infiltration of soil at the corresponding moment. The second term represents the introduction of interval constraints through the penalty function to ensure that the simulated cumulative infiltration amount I sim (L * ; n,1 / h d ,K s ) and the single-point cumulative infiltration when the wetting front first reaches the soil surface I * The deviation between them is maintained within the range of ±ε to consider the influence of experimental measurement error on the inversion accuracy. sim (L * ; n,1 / h d ,K s ) is calculated by formula (11) (i.e. z f =L*); σ is the penalty factor [–], ∈ is the constraint threshold [mm], and its value depends on the upper limit of the maximum measurement error in actual observation.

[0075] The above measured data (cumulative infiltration volume I i Over time t i The relationship between the change of the soil hydraulic properties and the field characteristics of the soil can be obtained by using a portable rapid determination device for soil hydraulic properties disclosed by the applicant (patent publication number: CN216209116U). A constant water head h is applied to the bottom of the ring cutter soil sample. p The measurement frequency is once per second. A camera is used to determine the exact time (t*) when the wetting front reaches the soil surface. From the measurement data (I~t), according to t * The cumulative infiltration amount I when the wetting front reaches the soil surface is obtained * . Record the cumulative infiltration volume I when the infiltration curve tends to be stable end . Measure the ring knife height L * .

[0076] Formula (14) is solved using the genetic algorithm (GA package) in R language: each generation of the algorithm generates 1000 sets of parameters to be estimated, n, 1 / h d and K s ; When the value of φ in formula (14) does not increase within 500 consecutive generations, or the cumulative number of iterations reaches 5000, the algorithm stops running; the local search strategy is enabled; the parameter search range is K s =[10-6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Logarithmic transformation is performed. The constraint threshold ε is set to 0.4 mm. σ is set to Here, .Machine$double.xmax represents the maximum double-precision floating-point number representable in the R language.

[0077] Example 2

[0078] This example conducts parameter inversion experiments based on the PASCP method ( Figure 2 ), the specific operations are as follows:

[0079] Step A: Infiltration test measurement phase:

[0080] (1) Apply a constant water head h at the bottom of the ring cutter soil sample p ;

[0081] (2) Record the relationship between the cumulative infiltration volume I of the soil sample and time t;

[0082] (3) Use a camera to record the time t when the wetting front first reaches the soil surface * .

[0083] Step B: Experimental parameter acquisition stage:

[0084] (1) From the measured data, according to t * The cumulative infiltration amount I when the wetting front reaches the soil surface is obtained * ;

[0085] (2) Determine the initial moisture content of the soil sample θ i and residual moisture content θ r ;

[0086] (3) Record the cumulative infiltration volume I when the infiltration curve tends to be stable end , calculate the saturated water content θ s ;

[0087] (4) Measure the ring cutter height L * .

[0088] Step C: Model parameter inversion stage:

[0089] (1) Constructing the optimization objective function based on the R language genetic algorithm (GA), including the RMSE term and the interval constraint term;

[0090] (2) Set the population size to 1000, the maximum number of iterations to 5000 (or 500 generations), and enable the local search strategy;

[0091] (3) The parameter search range is K s =[10 -6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Perform logarithmic transformation;

[0092] (4) Obtain soil hydraulic parameters n and 1 / h through optimization calculation d and K s The best estimate of .

[0093] The measured data in step A (the relationship between the cumulative infiltration volume I and the cumulative time t, I obs (t i )) was obtained using a portable field soil hydraulic properties rapid measurement device disclosed by the applicant (patent publication number: CN216209116U), with a measurement frequency of once per second;

[0094] The residual volumetric water content of the soil in step B is θ r The determination method is: the sieved soil is placed in a dry air environment with a relative humidity of less than 15% for several months until its mass remains constant. The water content of the soil at this time is the soil residual mass water content, and then multiplied by the bulk density of the ring cut soil sample to obtain the soil residual volume water content; the initial soil volume water content θ i The determination method is the drying method; the saturated volumetric water content of soil θ s Determined using formula (13).

[0095] The objective function in step C is specifically shown in formula (14). It is solved by using the genetic algorithm (GA) based on R language. The population size is set to 1000, and the algorithm is stopped after a maximum of 5000 iterations or 500 generations, provided that the fitness is not improved. The local search strategy is enabled. The parameter search range is K s =[10 -6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Logarithmic transformation is performed. The constraint threshold ε is set to 0.4 mm. σ is set to Where .Machine$double.xmax represents the maximum double-precision floating point number that can be represented in R language. The simulated infiltration amount I in formula (14) sim (L * ; n,1 / h d ,K s ) is calculated using formula (11); Simulated data I sim (t i ; n,1 / h d ,K s ) is calculated using formulas (6)-(10).

[0096] Example 3

[0097] This example aims to verify the ability of the PASCP method to handle measurement errors under conditions where actual measurement errors exist, and to demonstrate its improved robustness in parameter estimation compared to the ASCP method. Following the steps of Example 2, the following specific experiments were performed:

[0098] The collected soil samples were first air-dried and passed through a 2 mm sieve. The sieved soil was stored in a dry air environment with a relative humidity below 15% for several months until its mass remained constant. The residual mass water content of the soil was then determined gravimetrically. Based on a predetermined bulk density of 1.4, the residual water content θ was calculated. r =0.025cm 3 cm -3 The soil samples were then slowly packed layer by layer into a ring cutter with an inner diameter of 50.46 mm, a height of 50 mm, and a nylon mesh at the bottom. During the packing process, the ring cutter was gently tapped and roughened between layers to achieve a uniform bulk density. After the ring cutter was filled, the initial water content of the soil in the ring cutter was θ i =0.025cm 3 cm -3 .

[0099] For the filled ring cutter, the cumulative infiltration curve (I~t) is measured using a soil hydraulic properties rapid measurement device. Specifically, a constant water head h is used. p = 0cm. The infiltration data per second was collected by a portable infiltration device, and the time t when the wetting front reached the surface was recorded using a camera. * =3356s, corresponding to the cumulative infiltration amount I when the wetting front reaches the soil surface * =1.848cm. Record the cumulative infiltration volume I when the infiltration curve tends to be stable end =1.984cm, calculate the saturated water content θ s =0.422cm 3 cm -3 . Measure the ring knife height L * =5cm.

[0100] The optimization objective function (Formula (14)) was constructed using the genetic algorithm (GA) based on R language. Specifically, the population size was set to 1000, the maximum number of iterations was 5000 (or 500 generations), and the local search strategy was enabled; the parameter search range was K s =[10 -6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Perform logarithmic transformation; s ,θ i ,θ r , I * , L * and h p As a known parameter, by combining formulas (6)-(11) to calculate I in formula (14) sim (L * ; n,1 / h d ,K s ) and I sim (t i ; n,1 / h d ,K s ), the constraint threshold ε is set to 0.4 mm. Finally, the soil hydraulic parameters n and 1 / h are obtained by optimizing the calculation formula (14). d and K s The best estimate of .

[0101] Figure 3 The parameter iterative estimation processes of the PASCP and ASCP methods are demonstrated respectively. In order to evaluate the improvement effect of the PASCP method in parameter estimation robustness, the ASCP method is used as a comparative reference. Figure 3 Figure 1 shows that the PASCP method has a wider parameter search range, while the ASCP method, due to its strict equality constraints, has a limited search space. Interestingly, the parameter values ultimately determined by the ASCP method overlap with some of the parameter values found during the PASCP iteration. This phenomenon indicates that the PASCP method's parameter search space encompasses the ASCP method's search range, highlighting the PASCP method's robustness in the presence of large measurement errors. Figure 3 b further shows that the SWRC predicted by the PASCP method is more consistent with the measured value (RETC), while the prediction by the ASCP method shows obvious deviations (RMSE_ASCP=0.030, RMSE_PASCP=0.015, RMSE_RETC=0.010).

[0102] Example 4

[0103] This example aims to extensively verify the accuracy and parameter identification ability of the present invention in estimating soil hydraulic parameters. Five soils with significantly different textures are used as research objects (as shown in Table 1. The properties of the five soils include: soil texture composition (percentage of sand, silt, and clay), soil organic carbon content (SOC), bulk density (ρ b ), soil residual moisture content (θ r ), soil saturated water content (θ s ), initial soil moisture content (θ i ) and the single-point cumulative infiltration when the wetting front reaches the soil surface (I * The data in the table are actual measurements using a standard method (see Klute, 1986). These data are compared with the ASCP method to highlight the advantages of the present invention. Since the specific measurement procedures for the five soils are the same as in the previous examples, only the final results are presented here.

[0104] Figure 4 The PASCP and ASCP methods were used to show the degree of dispersion between the estimated parameters and the measured parameters for the five soil types mentioned above. The purpose was to verify whether the PASCP method can accurately capture the parameter variation trends between different soil textures. By comparing with the ASCP method, the PASCP method is more robust in parameter estimation. The regression analysis results show that the parameters estimated by the PASCP method (n,1 / h d , and log 10 (K s There is a significant linear correlation between the measured parameters (R 2 ≥0.804, P<0.05), indicating that the PASCP method can effectively reflect the parameter changes caused by soil texture differences. Among them, the estimated values of parameters n and log10(Ks) are closer to the measured values (Slope_n=0.866, Slope_log 10 (K s )=0.803); and the parameter 1 / h d The estimated values of the PASCP method exhibit some deviation, which may be related to the hysteresis effect in the measured soil moisture characteristic curve. In contrast, the ASCP method exhibits a significant decrease in parameter estimation accuracy under measurement error conditions, and the correlation and statistical significance between the parameter estimates and the measured values also decrease significantly. This further highlights the superior stability and reliability of the PASCP method under conditions of measurement error interference.

[0105] In summary, the PASCP method can not only effectively capture the parameter variation patterns among different soil textures, but also has excellent robustness and reliability in the presence of actual measurement errors.

[0106] Table 1 Properties of five soils

[0107]

[0108] Example 5

[0109] This example has two main purposes: First, to evaluate the soil saturation time (t s ) on the accuracy of parameter estimation; secondly, to compensate for the hysteresis effect in the actual measured moisture characteristic curve, which leads to the error of parameter 1 / h d Therefore, the advantage of the present invention in parameter estimation accuracy may be underestimated. Therefore, this embodiment uses HYDRUS-1D-based numerical simulation to select seven virtual soils with significant texture differences as research objects (see Table 2, the properties of the seven soils include: soil residual moisture content (θ r ), soil saturated water content (θ s ), initial soil moisture content (θ i ), pore distribution index n, inverse of intake suction force 1 / h d and saturated hydraulic conductivity K s ), each soil simulation was repeated ten times to evaluate the discreteness of the estimated parameters. The measurement steps used in the simulation are consistent with the previous embodiment, and only the simulation results are shown here.

[0110] Figure 5 The parameter inversion accuracy of the PASCP method is demonstrated on seven virtual soils with different textures. s ) as additional inversion constraint information ( Figure 5 a), PASCP estimated parameters (n, 1 / h d and K s ) shows a highly consistent linear relationship with the theoretical value (R 2 >0.99, P<0.001), where the estimated value of parameter n is almost consistent with the theoretical value. d and K s It is slightly overestimated, but its deviation shows a systematic trend, and PASCP can still accurately reflect the 1 / h difference between different soil textures. d and K s The changing law of R 2 ≥0.998). In addition, the standard deviations of all inversion parameters are low, indicating that this method can uniquely estimate the hydraulic parameters. In contrast, when the saturation time (t s) as the inversion constraint ( Figure 5 b), parameter n is significantly underestimated (Slope = 0.691), the estimation accuracy decreases and the standard deviation increases significantly. At the same time, the parameter 1 / h d and K s The prediction accuracy of the saturation time (t s ) as inversion constraint information, plays an important role in improving the accuracy and uniqueness of parameter estimation. In summary, the algorithm proposed in the present invention can accurately and uniquely estimate soil hydraulic parameters.

[0111] Table 2 Properties of seven soils

[0112]

Claims

1. A robust calculation method for soil hydraulic parameters based on one-dimensional upward infiltration test, characterized by The steps include: Step A: Apply a constant lower boundary water head h to the ring cutter soil sample p , carry out a one-dimensional upward infiltration test on homogeneous soil, record the relationship between the cumulative infiltration volume I of the ring cut soil sample and the time t, the time t* when the wetting front first reaches the surface of the ring cut soil sample, and the cumulative infiltration volume I when the infiltration curve tends to be stable end ; Step B: Based on the initial volumetric water content θ of the ring cut soil sample i , the ring knife height L* and the cumulative infiltration amount I when the infiltration curve tends to be stable end , calculate the saturated volumetric water content of the soil θ s , the calculation formula is: Step C: Based on the relationship between the cumulative infiltration volume I and time t obtained in step A, determine the single-point cumulative infiltration volume I* corresponding to the time t* when the wetting front first reaches the surface of the ring cutter soil sample; Step D: Combine the saturated volumetric water content θ of the ring cut soil sample s , initial volume water content θ i , residual volume water content θ r , the single-point cumulative infiltration volume I* when the wetting front first reaches the soil surface, and the observed data of cumulative infiltration volume changing with time I obs (t i ), the parameter optimization inversion method based on interval constraints is used to obtain the estimated parameters n, 1 / h d and K s .

2. The robust calculation method for soil hydraulic parameters according to claim 1, characterized in that: In step D, the objective function of the interval-constrained parameter optimization inversion method is constructed according to formula (14):

3. The robust calculation method for soil hydraulic parameters according to claim 2, characterized in that: The first term of formula (14), I obs (t i ) is the measured value obtained in step A, representing the observation time point t i The cumulative infiltration volume corresponding to the time I, I sim (t i ; n,1 / h d ,K s ) is the observation time point t i The simulated cumulative infiltration volume is calculated according to the following formula: ① Time for soil to reach saturation (t s ) is: Where z s * is the equivalent wetting front z fe Function expression, when the equivalent wet front z fe Reaching the soil surface (i.e. fe =L * ), the soil reaches saturation, at this time z s * It can be clearly expressed as: ②For infiltration time t <t s When , the functional relationship between soil cumulative infiltration I and time t is: ③When t≥t s When the soil reaches full saturation, the cumulative infiltration remains unchanged, and the calculation formula is as follows: I=(θ s -θ i )L * (10); The parameters to be estimated n, 1 / h d and K s Substitute it into the measured value I obs (t i ) and analog value I sim (t i ; n,1 / h d ,K s ) tends to be the smallest; The second term of formula (14), I sim (L * ; n,1 / h d ,K s ) is the cumulative infiltration volume of the simulation, which is calculated according to the following formula: where z f =L*, σ is the penalty factor [–], ∈ is the constraint threshold [mm], and its value depends on the upper limit of the maximum measurement error in actual observation. The remaining parameters are expressed according to the following relations:

4. The robust calculation method for soil hydraulic parameters according to claim 3, characterized in that: Formula (14) is solved using the genetic algorithm in R language: each generation of the algorithm generates 1000 sets of estimated parameters n, 1 / h d and K s When the value of φ in formula (14) does not improve within 500 consecutive generations, or the cumulative number of iterations reaches 5000, the algorithm stops running; Enable local search strategy; parameter search range is K s =[10 -6 ,10 -2 ]cm·s -1 , 1 / h d =[0.01,1]cm -1 , n=[0.1,0.7], for K s and 1 / h d Logarithmic transformation is performed; the constraint threshold ε is set to 0.4 mm; σ is set to Here, .Machine$double.xmax represents the maximum double-precision floating-point number representable in the R language.

Images