土壤水力性质是水分运动方程的重要参数，主要包括水力传导率K（θ）、土壤水分特征曲线θ（h）等，是建立土壤水分运动数学模型的基础。参数的可靠性直接影响水分运动模型预测的准确性。尽管有许多实验室和田间方法可以用来确定包气带中以土壤水分特征曲线和水力传导率为代表的土壤水力性质，但大多数方法需要静态或稳态流条件来满足相应分析方法的假设条件或需要较精密的仪器，这使得测量耗时耗力，可行性不高，特别是当土壤的空间变异性很强且区域范围较大时，更增加了测量的难度。因此，人们开始发展间接方法来估计土壤水力性质。间接方法有土壤转换函数法、物理-经验法、分形几何法、土壤形态学网络模型和数值反演方法等^[1]。其中数值反演方法由于具有以下优点受到研究者的关注^[2-4]：（1）无需精密的测量仪器；（2）不需要达到稳态水流，在初始和边界条件上较直接方法更加灵活；（3）可以从单个瞬态流实验同时估算土壤水分特征曲线和水力传导率；（4）为优化参数提供置信区间。但数值反演方法在确定土壤水力性质时通常存在解的非唯一性问题。

关于如何减少数值反演方法解的非唯一性和提高参数估计的准确性，国内外学者做了许多相关研究^[5-7]。Toorman等^[8]借助响应曲面方法对目标函数进行分析，使用出流量、压力水头、含水量的不同组合作为目标函数的测量变量，结果表明，如果在目标函数中包含压力水头值，对参数估计的敏感性将一定程度有所提高。此外，压力水头的测量位置对参数估计也有影响，可以据此确定需要测定的变量及位置。van Dam等^[9]通过出流实验模拟结果，说明仅采用累计出流量数据容易导致解的非唯一性，而加入θ（h）信息可以得到可靠的参数解。Eching等^[10-11]在多步出流实验的累积出流量数据基础上添加测定的土壤压力水头数据，以三次不同初始参数反演结果的变异系数表征解的唯一性，在反演土壤水分特征曲线和非饱和水力传导率时得到了较好的拟合效果，验证了增加压力水头信息可以提高解的唯一性。Šimunek和van Genuchten^[12]应用了张力圆盘入渗数据进行数值反演，求解水力特征参数，通过参数与目标函数的响应曲面分析使用哪些入渗数据可以提高解的唯一性，认为可以适当增加剖面的压力水头信息和含水量信息来提高计算准确性，而加入累积入渗量信息不能改善拟合结果。Zijlstra和Dane^[13]将反演方法应用于层状土壤，土壤剖面由一个、两个或三个不同的层位组成，设定水分含量测量值为时间和深度的函数，下边界设为压力水头关于时间的函数，以最小化目标函数法优化参数，结果认为随着优化参数数量的增加，反演问题可能变得不适定，反演的不适定性取决于输入数据集的大小。张俊等^[14-15]利用数值反演方法耦合瞬时出流实验推求土壤水力学参数，结果显示，单步出流实验在目标函数中增加压力水头信息、多步出流实验联合累积出流量和压力水头信息，均可以减小反演问题的非唯一性。

综上所述，利用数值反演方法推求土壤水力性质方法虽然有多种，前人也做了关于降低数值反演方法解的非唯一性和提高参数估计准确性的相关研究，但通常只考虑单层土壤，而对层状土壤反演其土壤水力性质的研究还较少，需要深入开展这方面的工作。本研究结合室内层状土柱入渗试验，运用Hydrus-1D^[16]中的反演（inverse solution）模块求解层状土壤水力性质，并分析目标函数中包含累积入渗量、土壤压力水头及两者组合对解的唯一性及准确性的影响。以此探究获得可靠层状土壤水力性质需要测量的试验变量。

1 材料与方法 1.1 供试土壤

土样取自青岛市即墨区移风店镇农田，土壤类型为棕壤。经测定土壤的基本物理性质如表 1所示。

由表 1结果可知，采样点的土体构型相对属于“上粗下细型”。土壤质地分为两种：0~50 cm为质地较粗的砂壤，其砂粒含量明显高于下层，渗透性较好；50 cm以下为粉砂壤，下部土壤黏粒含量较高，渗透性较差。

1.2 试验方法

试验采用室内层状土柱定水头积水入渗方法，有机玻璃柱填装土柱，马氏瓶供水，使马氏瓶内管底端与液面相平以维持2 cm恒定水头（图 1）。有机玻璃柱规格：内径19 cm，高80 cm，一侧开六孔安装张力计，上端设进水口，下端设出水口。马氏瓶规格：内径9 cm。土柱分三层按田间土层排列顺序填装，从上至下，土层厚度分别为20 cm、20 cm和25 cm。由于田间原状土容重较大，室内土柱无法还原田间状态，故试验按照设计容重填装，从上至下容重分别为1.57 g·cm^–3、1.45 g·cm^–3、1.40 g·cm^–3。土柱表层和底层分别铺设2 cm、5 cm石英砂，由于石英砂质地与表层土相差较小，与底层相差较大，为了减小模拟误差，上部土层厚度概化为22 cm，下部土层厚度概化为27 cm，土柱总厚度为69 cm。6个张力计从上至下分别放置在9、19、29、39、49和59 cm深度处。土柱底端张力计读数开始变化时停止通水。试验全程监测各张力计读数和马氏瓶读数。

1.3 土壤水分运动模型的建立

Hydrus-1D模型是由国际地下水模型中心公布的数值模拟软件，用来模拟一维饱和-非饱和多孔介质中水、热、溶质的运移过程。

土壤水分以垂向运动为主，因此模型采用一维垂直入渗Richards方程描述土壤水分运动：

$C\left( h \right)\frac{{\partial h}}{{\partial t}} = \frac{\partial }{{\partial z}}\left[ {K\left( h \right)\frac{{\partial h}}{{\partial z}}} \right] - \frac{{\partial K\left( h \right)}}{{\partial z}}$

(1)

土壤水力函数用van Genuchten-Mualem模型表示：

$ \theta \left(h\right)=\left\{ \begin{array}{l}{\theta }_{r}+\frac{{\theta }_{s}-{\theta }_{r}}{{\left[1+{\left|\alpha h\right|}^{n}\right]}^{m}}h<0\\ {\theta }_{s}{\rm{\;\;\;\;\;\;\;\;\;\;}}h≥0\end{array} \right.$

(2)

$K\left( h \right) = {K_s}{S_e}^l{\left[ {1 - {{(1 - S_e^{\frac{1}{m}})}^m}} \right]^2}$

(3)

${S_e} = \frac{{\theta - {\theta _r}}}{{{\theta _s} - {\theta _r}}}$

(4)

式中，C（h）为比水容量（cm^–1），K_s为饱和水力传导率（cm·min^–1），K（h）为非饱和水力传导率（cm·min^–1），S_e为饱和度，θ_s为土壤饱和含水量（cm³·cm^–3），θ_r为残余含水量（cm³·cm^–3），α为进气吸力相关的参数（cm^–1），m和n为形状系数，h为压力水头（cm）；z为土壤深度（cm）；t为时间（min）。其中，m=1–1/n，l=0.5。

以上方程包含了5个独立参数，θ_r，θ_s，α，n和K_s，未知参数通过最小化目标函数获取^[15]。参数的优化计算采用Hydrus-1D，version 4.16（US Salinity Laboratory）。Hydrus-1D中目标函数最小化的方法采用LM法，该方法结合了牛顿法和最速下降法，为优化参数提供置信区间，并作为一种标准方法在水力学中求非线性拟合最小平方和而被普遍采用^[17-19]。

初始条件：

$ h={h}_{0}(z)\;\;\;\;0<z<L \;\;\;\;t=0$

(5)

上边界：

$ h\left(0, t\right)={H}_{0}\;\;\;\;z=0\;\;\;\;t>0$

(6)

下边界：

$\frac{{\partial h}}{{\partial Z}} = 0 \;\;\;\; L = 69{\rm{cm}} \;\;\;\; t > 0$

(7)

式中，h₀（z）为模拟开始时的土壤压力水头值，监测点处为实测值，其余节点采用插值法计算得到，H₀为土壤表面积水深度（2 cm）。

时间和空间离散：对试验条件下的一维水分运动进行数值模拟，模拟时长为T=1 220 min，对应于59 cm深度处水势传感器读数开始变化的时间。计算过程采用可变时间步长。空间步长为0.5 cm，依据土柱总厚度69 cm将土壤剖分为138个单元，共139个节点。

初始参数的选取：Hydrus-1D模型各层土壤水力学参数的选取至关重要，虽然最小化目标函数方法计算复杂度低，但对参数初始值敏感，且得到的优化参数可能是局部最优，而非全局最优，因此需要尝试不同的初始参数值。对于van Genuchten- Mualem中的六个参数θ_r，θ_s，α，n，K_s，l，其中l通常设为0.5。θ_r采用风干土的质量含水量与体积含水量的关系计算得到，θ_s通过试验时土柱内各层土壤容重计算出的孔隙度获得。为降低反演的复杂程度，反演计算时固定参数θ_r和θ_s（表 2），主要对3层土壤的α，n，K_s共9个参数进行优化。Carsel和Parrish^[20]根据概率分布的方法分析大量数据得到了不同土壤类型的van Genuchten非饱和水力参数，选取相应土质的α，n，K_s参数值作为反演的初值，并根据反演效果调整反演时α，n，K_s的初值，以达到更好反演效果，调整后的参数见表 3。

1.4 评价指标

为了探究层状土壤一维水分运动中测量变量的类型和数量对数值反演参数唯一性和准确性的影响，分别将累积入渗量（Q）、不同位置压力水头（H）及两者组合（Q+H）作为目标函数中的测量变量，对土壤水力学参数进行优化求解。每次反演过程使用3组不同初始值（表 3）进行运算。对目标函数为Q、H和Q+H的9次水分运动模拟过程分别用情景（Ⅰ，Ⅱ，Ⅲ）、（Ⅳ，Ⅴ，Ⅵ）和（Ⅶ，Ⅷ，Ⅸ）表示。采用R²，RMSE指标对拟合效果进行定量化评价。R²也称为决定系数，代表了回归方程对实测值的拟合度，值越接近1，拟合效果越好。RMSE为均方根误差，值越小代表拟合效果越好。对于解的非唯一性用变异系数（CV_s）进行评价，并结合优化参数的相关性矩阵以及优化参数的标准误差对反演效果进行分析。

2 结果与讨论 2.1 不同目标函数的拟合结果

图 2表示目标函数中测量变量为累积入渗量时拟合结果，图 3表示目标函数中测量变量为压力水头时拟合结果，图 4和图 5为目标函数中测量变量为累积入渗量与压力水头组合时拟合结果。由图 2~图 5可以看出模拟值与实测值符合程度较好，但19 cm处压力水头的模拟较差。导致该结果的可能原因如下：（1）19 cm处张力计读数受气泡影响读数上升较慢，未能及时反映出压力水头的变化，即由于仪器引起的测量误差；（2）19 cm处张力计位置靠近中层土壤上界面，土体构型为上粗下细，进水受空气阻力影响，水分在19 cm处运移较慢，土壤达到饱和需要较长时间，而模型未考虑进气阻力致使模型默认水分到19 cm处后很快达到饱和，因此模拟值与测量值吻合度较差。

2.2 不同目标函数的模拟精度

表 4为9种情景的反演效果定量化评价，总体拟合效果良好，R²均在0.855~0.999之间，RMSE值在0.005 50 cm~0.038 4 cm之间。三种目标函数下R²大小顺序为Q > Q+H > H，RMSE值大小顺序为Q < Q+H < H。结果表明：以R²和RMSE值为评价依据，且目标函数测量变量为累积入渗量时，拟合效果最优；压力水头与累积入渗量组合时效果次之；只用压力水头数据反演拟合效果较差。这可能是由于19 cm处压力水头模拟值和实测值相差较大导致的。

2.3 不同目标函数反演结果的变异系数

由表 5可以看出，目标函数仅有累积入渗量时，变异系数总体最小，目标函数仅有压力水头时变异系数总体最大。这可能由于累积入渗量实验误差较小，而19 cm处压力水头模拟值与实测值相差较大，导致目标函数为压力水头时变异系数较大。当目标函数为压力水头和累积入渗量的组合时，变异系数总体较目标函数只有压力水头时小。变异系数越小，解的非唯一性越小。说明在压力水头数据基础上加入累积入渗量数据可以减小解的非唯一性。

2.4 不同目标函数反演参数的相关性

参数相关矩阵反映了两个参数值之间的相关性。值为–1或+1表示完全线性相关，而0表示不相关。参数之间的高度相关性会导致低估参数的不确定性，降低收敛速度并增加非唯一性^[21]。根据Hydrus-1D软件模拟结果反演信息（Inverse solution information）中的相关性矩阵（correlation matrix），分析不同目标函数在9种情景下数值反演的优化参数之间的相关系数，相关系数绝对值大于0.9的数值表示相关性显著。通过对比不同目标函数下参数间相关性显著的数量，分析参数的不确定性和非唯一性。

三种目标函数情况下反演参数相关性显著的总数量，目标函数为Q、H和Q+H时分别为24、6和7个。由此可见，当目标函数中只有Q时，反演参数相关性显著的总数明显多于其他两种情况。说明目标函数只包含累积入渗量时，参数之间相关性较强。此时模型建议固定相关性较强的两个参数中的一个，以便于对另一个参数更准确求解，否则在此情况下得出的参数优化结果不确定性和非唯一性很强，易出现异参同效现象。这主要是由于层状土壤水力性质参数较多和水分运动模型的高度非线性造成的。而目标函数为压力水头及累计入渗量和压力水头组合时，参数相关性显著的数量基本相等，且大幅减少，这是由于目标函数中含压力水头时，提供了不同深度处土壤压力水头随时间变化的信息，限定了更多约束条件，减少了异参同效现象。因此，目标函数中包含压力水头时对降低反演过程参数非唯一性有明显效果。

2.5 不同目标函数反演参数的标准误差

Hydrus-1D提供了关于拟合参数的统计信息，其中参数的标准误差是度量拟合结果精确度的指标，由最优值附近的95%置信区间确定，标准误差越小，参数预测的95%置信区间越窄，拟合结果确定性越高。对三种目标函数下优化参数的标准误差进行分析，9种情景下优化参数的标准误差列在表 6中，9个参数的标准误差中，除α₁之外，其余参数的标准误差大小顺序均为Q+H < H < Q。表明累积入渗量与压力水头组合的参数拟合精确度最高，压力水头次之，累积入渗量最差。

2.6 验证反演参数的准确性

为进一步比较目标函数中不同测量变量对数值反演求得的参数准确度的影响，对累积入渗量反演优化的参数结果（表 7）进行正演，预测各观测点压力水头值随时间变化，并与实测值进行对比（图 6），验证参数的准确性^[15]。同样方法对压力水头反演优化的参数结果用累计入渗量预测值与实测值进行对比，得出结果如图 7所示。

由图 6可以看出目标函数只包含累积入渗量时所求参数对压力水头的预测值与实测值相差较大，三次模拟R²分别为0.334、0.338、0.332，RMSE分别为0.096 1 cm、0.095 6 cm、0.096 5 cm。由图 7可知目标函数值包含压力水头时所求参数对累积入渗量的预测值与实测值符合程度较高，三次模拟R²分别为0.996、0.994、0.993，RMSE分别为0.028 76 cm、0.034 02 cm、0.038 56 cm。由此可见运用累积入渗量数据对土壤水力学参数进行数值反演的准确性较差，而目标函数中包含压力水头时反演结果准确性较好，基本能反映累计入渗量和压力水头随时间的变化关系。

3 结论

目标函数仅含有累积入渗量时反演模型模拟值与实测值符合度虽然较高，决定系数R²达到0.999，均方根误差RMSE小于0.005 62 cm；但用压力水头实测值验证参数准确性的效果较差，R²最高仅为0.338，RMSE小于0.096 5 cm。使用压力水头作为目标函数进行参数优化求解，拟合度R²达到0.855以上，RMSE小于0.038 4 cm；对反演得到的参数用累计入渗量实测值进行验证，R²可达到0.99以上，RMSE小于0.038 56 cm，效果较好。累积入渗量和压力水头联合反演求土壤水力性质拟合度较高，R²达到0.905以上，RMSE小于0.035 4 cm。虽然累积入渗量和压力水头联合反演的R²和RMSE指标次于累积入渗量，但累积入渗量和压力水头联合反演可以一定程度降低优化参数的相关性和标准误差，从而降低解的非唯一性，提高解的准确性。综上所述，由于层状土壤水力性质的参数较多，数值反演方法求参数的不确定性和非唯一性较强，若单独使用一种测量数据进行反演计算，很难确定参数的准确性和唯一性。此外，若这一测量变量误差较大，将严重降低结果准确性。因此数值反演过程应尽量使用不同类型和较多数量的测量变量以增加约束条件，提高数值反演解的唯一性、确定性和准确性。