HISR Method for Parameter Inversion of Solute Transport through Unsaturated Soils
-
摘要: 以非饱和土中溶质迁移参数反演问题为背景, 依据正则化方法的思路, 以Itakura Saito距离作为同伦函数中的平凡问题, 将同伦方法引入非线性参数反演问题的求解, 进而提出一种求解非线性参数反演问题的大范围收敛(HomotoyItakura SaitoRegularization, HISR) 方法.为保证迭代稳定性, 并同时削弱观测噪声的影响, 同伦参数的修正采用了连续化修正方法.本文将HISR方法应用于求解带有平衡及非平衡吸附效应的一维非饱和土中溶质迁移参数反演问题, 计算结果表明HISR方法具有大范围收敛性及计算稳健性, 同时有较强的抵抗观测噪声的能力.Abstract: Taken the parameter inversion of solute transport through unsaturated soils as the application background and homotopy method imported for solving the nonlinear parameter inversion problem, this paper presents a global convergence method, HISR (homotopy Itakura-Saito regularization) method. The Itakura-Saito distance between the iterative parameter value and the estimated parameter value was adopted as the trivial problem in the homotopy function based on the idea of regularization. Moreover, a continuous method is employed to modify the homotopy parameter during iteration in order to ensure the ability of the iteration and the capability of the noise resistance. Problems of parameter inversion of solute transport coupled with equilibrium and non-equilibrium effects through one-dimensional unsaturated soils are carried out as numerical examples, and the computational results clearly demonstrated the features of global convergence and stability of the HISR method. Besides, a favorable solution is obtained even though experimental quantities were contaminated heavily by noise by HISR method.
-
研究溶质在土壤和地下水中的迁移规律时, 参数反演是重要内容之一.溶质在非饱和土中的迁移转化方程可以表达为耦合不同效应的对流扩散方程形式(Nielsen et al., 1986), 但该方程中的某些参数很难直接测得, 通常需要借助于某些实验手段.依据计算结果与实验结果最佳拟合的原则, 借助于数值反演方法来确定其数值, 这就是参数反演问题.
非饱和土中溶质迁移参数反演问题可以转化为非线性无约束优化问题进行求解(Sun, 1994).基于Levenberg-Marquardt迭代格式的非线性最小二乘方法是较早用于求解该类问题, 且目前仍广泛使用的一种方法. Parker and Van Genuchten (1984)所开发的一维溶质迁移过程的水动力学参数及溶质迁移参数估计软件包CXTFIT即基于这种方法.该软件包在实际中已得到广泛使用(Torsten, 1998; Pang and Murray, 1999; Kamra et al., 2001), 并由Toride等(1995) 进行完善和改进.Levenberg-Marquardt方法是对传统的Gauss-Newton法的改进, 该方法具有较好的计算效率, 相对于Gauss-Newton法而言, 其收敛范围也大了很多.但在实际计算中, 若迭代初值选取不好, 仍可能造成迭代的发散(刘钦圣, 1989).启发式算法(如模拟退火(Barry and Morris, 1999) 等) 可以对此进行改善, 并可能寻求到全局最优解.但该类方法本质上属于随即搜索方法, 相对而言, 其计算效率较低.
同伦方法(王宇, 1996) 是求解非线性问题的一种有效方法.该方法对于求解初值的选取没有严格限制, 能够做到全局收敛.本文以带有平衡及非平衡吸附效应的一维非饱和土中溶质迁移规律研究为背景, 将同伦方法引入非线性参数反演问题的求解; 同时在分析了同伦方法正则化效应的基础上, 以Itakura-Saito距离(El-Jaroudi and Makhoul, 1991) 作为同伦函数中的平凡问题, 进而提出求解非饱和土中溶质迁移参数反演的一种方法, 我们称之为HISR方法.为保证迭代稳定性, 并同时削弱观测噪声的影响, 采用了连续化同伦修正方法.数值算例依据数值模拟一维实验室土柱实验产生的模拟实验数据, 对纵向弥散度α、吸附系数K、沉淀速度KA、裂解速度KB和相对最大沉淀量γ进行反演计算.计算结果表明了该方法的大范围收敛性及计算的稳健性.
1. 溶质在非饱和土中的迁移规律模型
溶质在非饱和土中的迁移转化是一个十分复杂的物理化学过程.这里我们考虑平衡吸附及非平衡吸附效应.平衡吸附一般可以使用Henry公式、Freundlich公式或Langmuir公式表达(Nielsen et al., 1986), 若使用Henry公式, 则平衡吸附方程可以表达为:
(1) 其中: R为迟滞系数, 即R=θ+ρK; C为待求的溶质浓度; K为吸附系数; ρ为干土样的密度(g/cm3); θ为饱和度(体积含水率); α为纵向弥散度; ν为液相流速(cm/h); q为达西流速(cm/h); z为空间坐标(cm).
非平衡吸附包括物理非平衡吸附及化学非平衡吸附, 相应的可以采用2区域(two-region) 模型及2位置(two-site) 模型描述.虽然物理非平衡吸附和化学非平衡吸附基于不同的概念, 但可归于相同的表达形式(Pang and Murray, 1999), 对于一维非饱和土中的溶质迁移过程, 相应的方程可以表达为:
(2) 
(3) 其中: S为吸附浓度; β为交换系数; D为分配系数, S0为固相中最大吸附浓度.
方程(2)、(3) 的初始及边界条件根据相应的实验室土柱实验过程确定.所谓实验室土柱实验(图 1a), 就是将实验用土样封装于一维土柱中, 在土柱入流端分时不断供入一定浓度的溶液或蒸馏水, 在土柱出流端定时收集出流溶液, 进而测得出流溶液浓度随时间变化的曲线, 即穿透曲线(图 1b).穿透曲线一般作为反演的依据.
图 1 实验室土柱实验(a) 及穿透曲线(b) 示意图Fig. 1. Sketch map of laboratory soils column experiment (a) and breakthrough curve (BTC) (b)为计算方便, 常以入流溶液的原始浓度C0做为浓度C的特征量, 以最大吸附浓度S0作为吸附浓度S的特征量, 将方程(2) 和(3) 做无量纲化处理(Li et al., 1999), 无量纲的吸附浓度以S2表示, 无量纲的浓度仍以C表示, 整理后可得:
(4) 
(5) 其中: KA为沉淀速度, KB为裂解速度, γ为相对最大沉淀量, 均为无量纲参数; L为土柱长度(cm).
当入流采用方波输入时(即首先注入溶液, 模拟实际中的污染过程, 一段时间后, 再注入蒸馏水, 以模拟实际中的清洗过程), 入流端边界条件给定如下:
(6) t1和t2分别为入流端加入溶液及加蒸馏水的时间.出流端边界条件可以给定为第二类边界条件(王超, 1996; Inoue, 2000) :
(7) 
(8) 式(4) ~ (8) 即为考虑平衡及非平衡吸附效应时溶质在非饱和土中迁移的控制方程.当相关参数给定时, 通过数值方法(如有限元法) 可以解得出流浓度C.
实际中, 纵向弥散度α、吸附系数K、沉淀速度KA、裂解速度KB和相对最大沉淀量γ很难通过实验直接测得.可根据土柱实验所获得的穿透曲线来反演得到.将实验所测得的出流端溶液浓度分布记为
(9) 式(9) 即为土壤中溶质迁移的参数反演问题的一般表达式, 亦即本文所主要关心的问题.
2. HISR方法
依据同伦方法的思想(王则柯和高堂安, 1990;王宇, 1996), 构造式(9) 的线性同伦函数:
(10) 进而将原问题转化为其同伦函数的零点求解问题, 即:
(11) 式(10) 中, λ为同伦参数.
式(9) 所示的非线性参数反演问题一般为不确定的, 即解的存在性、唯一性很难保证, 同时解不连续依赖于可观测量.对于不适定问题的求解, 正则化方法(Tikhonov and Arsenin, 1977; Engl and Hanke, 1996) 是一种很有效的方法.该方法主要通过引入正则化泛函及正则化参数, 使用一个连续的问题来逼近原不适定问题, 同时通过选择适当的正则化泛函为解提供先验信息, 从而有效地抵抗原问题的不适定性.
依据正则化方法的思想, 将式(10) 中的平凡问题
(12) 式中, m为参数的个数, pi0为参数真实值的预估值, 一般选取迭代中上一步所求得的参数值.
选取式(12) 的原因在于: ①该函数是一个连续的光滑函数; ②随着迭代的进行, 该函数的值将不断减小, 从而可以使得同伦函数(10) 不断地向原问题(9) 逼近.
实际问题中不可避免地包含有观测噪声, 因此方程(11) 一般不能精确满足.由于同伦函数,
(13) 依据预估-校正算法(王宇, 1996) 的思路, 对于问题(13) 采用如下思路求解: 给定同伦参数λ及
按照上述思路, 在每迭代步中, 对问题(13) 做线性化.假定迭代进行至第n步, 将第n+1步计算得到的溶质边界浓度值在第n步计算得到的溶质边界浓度值附近展开为一阶Taylor级数, 并略去高阶小量, 得:
(14) 将(14) 式代入(13) 式, 并取(13) 式中的
(15) 问题(15) 中, H的雅戈比矩阵和海塞矩阵分别为:
(16) 
(17) 同伦参数λ依据下式调整:
(18) 式(18) 中, n为迭代步数, 相应的λn为第n个迭代步时同伦参数的取值.β为修正步长调整参数, 若采用较大的β值, 可以增大同伦参数的修正步长, 加速求解的收敛速度, 但过大的β值可能会引起迭代的发散.N0为同伦参数预估初值的调整参数, 该参数越大, 则同伦参数的初始预估值越接近1.对于本文的问题, 建议β的取值介于0.5~0.8之间, N0可选取为1~5之间的整数.
依据式(18) 调整同伦参数, 其原因如下: ①同伦参数的变化是连续的, 这样可以使得同伦参数稳定地下降, 有利于求解过程稳定地进行; ②随着n的减小, 相应的同伦参数会趋近于零, 而不会等于零, 从正则化的角度考虑, 这也可以在一定程度上抑制观测噪声的影响.
综上分析, 可以最终得到HISR方法的迭代公式如下:
其中ᐁH及ΔH分别根据式(16) 及(17) 求得.
3. 数值算例
反演计算时, 在方程(4) ~ (8) 中可直接测得的实验参数(模拟实验参数) 给定的前提下, 首先给定一组待反演参数的数值(参数真实值), 即:
表 1 模拟实验参数数值Table Supplementary Table Value of simulated experimental parameters
算例 1 选取参数真实值
表 2 不同初始值时的参数反演结果Table Supplementary Table Numerical inversion results for different initial values of parameters
表 2中, IT为迭代次数.从表 2可以看出, 当待反演参数的任一参数求解初始值选取距离真值较远时(结果No.2~No.5), 迭代仍能准确地收敛到正确的结果.同时, 在求解初始值全部选取为距离真值较远的一般估计值[0.1, 0.1, 0.1, 0.1, 0.1]的条件下(结果7), 计算也能精确地收敛到参数真值.这说明了HISR方法具有较大的收敛范围, 对于求解初始值较难估计的实际问题也具有实用价值.该算例同时表明, 选取不同的求解初始值时, 迭代次数基本保持不变, 这说明该方法具有很好的计算稳健性.
算例 2 选取参数真实值
表 3 不同噪声强度下的反演结果Table Supplementary Table Numerical inversion value for different noise levels
表 3中, NA为模拟观测噪声的相对强度, IT为迭代次数, Err为反演所得参数与参数真实值间的2范数距离.
从算例2的计算结果可以看出, 在较小观测噪声情况下, HISR方法及LM方法均能收敛到很好的解点, 且其计算误差基本相同; 而在存在较大观测噪声的情况下, 采用HISR方法进行反演计算仍能收敛到一个较好的解点, 但LM方法可能发散(如误差强度为10%时), 同时, HISR方法的反演结果略好于LM方法.此外, 在不同的观测噪声强度下, 从同一求解初始点出发, 至迭代收敛时所用的迭代步数基本相同, 且明显少于LM方法.这说明在存在观测噪声的情况下, 该方法仍具有很好的计算效率及计算稳健性.
较大模拟观测噪声强度下(表 3的后4组结果, 相对噪声强度分别为4%, 8%, 10%, 15%), 根据反演结果计算得到的穿透曲线与模拟试验结果穿透曲线及真实穿透曲线的拟合情况如图 2所示.
图 2 模拟实验穿透曲线(○)、反演穿透曲线(实线) 与真实穿透曲线(虚线) 的拟合a.噪声强度为4%;b.噪声强度为8%;c.噪声强度为10%;d.噪声强度为15%Fig. 2. BTCs fitting of simulated experimental results, inverse results and true results4. 结论
从计算结果可以看出, 即使求解初值与参数真实值距离较远, 采用HISR方法求解仍能收敛到正确的解点, 同时迭代次数基本不变, 这说明HISR方法具有很大的收敛范围及较好的计算稳健性, 采用HISR方法求解非饱和土中溶质迁移参数反演问题是可行而有效的.此外, 应该说明的是, HISR方法要求待求参数取值应为非负, 这与本文问题中对待求参数的要求是一致的.本文的方法也可以拓展至其他待求参数为非负的非线性参数反演问题.
由于考虑了正则化效应, 使得HISR方法也具有一定的抵抗观测噪声的能力, 同时在存在观测噪声的情况下该方法仍具有很好的计算稳健性.这说明对于存在观测噪声的实际问题, HISR方法也是一种很有效的求解方法.但随着观测噪声的增大, 反演结果的误差也逐渐增大.对于正则化方法而言, 合理地选择正则化参数可以在相当程度上抵抗观测噪声(Engl and Hanke, 1996).因此, 对HISR方法中的同伦参数修正方法作进一步的探讨, 使得在保证计算稳定性的前提下, 进一步提高HISR方法抵抗强观测噪声的能力将是值得继续研究的问题.
-
图 1 实验室土柱实验(a) 及穿透曲线(b) 示意图
Fig. 1. Sketch map of laboratory soils column experiment (a) and breakthrough curve (BTC) (b)
图 1 实验室土柱实验(a) 及穿透曲线(b) 示意图
Fig. 1. Sketch map of laboratory soils column experiment (a) and breakthrough curve (BTC) (b)
图 2 模拟实验穿透曲线(○)、反演穿透曲线(实线) 与真实穿透曲线(虚线) 的拟合
a.噪声强度为4%;b.噪声强度为8%;c.噪声强度为10%;d.噪声强度为15%
Fig. 2. BTCs fitting of simulated experimental results, inverse results and true results
图 2 模拟实验穿透曲线(○)、反演穿透曲线(实线) 与真实穿透曲线(虚线) 的拟合
a.噪声强度为4%;b.噪声强度为8%;c.噪声强度为10%;d.噪声强度为15%
Fig. 2. BTCs fitting of simulated experimental results, inverse results and true results
表 1 模拟实验参数数值
Table 1. Value of simulated experimental parameters
表 1 模拟实验参数数值
Table 1. Value of simulated experimental parameters
表 2 不同初始值时的参数反演结果
Table 2. Numerical inversion results for different initial values of parameters
表 2 不同初始值时的参数反演结果
Table 2. Numerical inversion results for different initial values of parameters
表 3 不同噪声强度下的反演结果
Table 3. Numerical inversion value for different noise levels
表 3 不同噪声强度下的反演结果
Table 3. Numerical inversion value for different noise levels
-
El-Jaroudi, A., Makhoul, J., 1991. Discrete all-pole modeling. IEEE Trans. Acoustics, Speech, Signal Processing, 39(2): 411-423. doi: 10.1109/78.80824 El-Jaroudi, A., Makhoul, J., 1991. Discrete all-pole modeling. IEEE Trans. Acoustics, Speech, Signal Processing, 39(2): 411-423. doi: 10.1109/78.80824 Engl, H.W., Hanke, M., 1996. Neubauer, regularization of inverse problems. Kluwer, Dordrecht. Engl, H.W., Hanke, M., 1996. Neubauer, regularization of inverse problems. Kluwer, Dordrecht. Inoue, M., 2000. Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments. Advances in Water Resources, 23: 677-688. doi: 10.1016/S0309-1708(00)00011-7 Inoue, M., 2000. Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments. Advances in Water Resources, 23: 677-688. doi: 10.1016/S0309-1708(00)00011-7 Kamra, S.K., Lennartz, B., Van Genuchten, M.T., et al., 2001. Evaluating non-equilibrium solute transport in small soil columns. Journal of Contaminant Hydrology, 48: 189-212. doi: 10.1016/S0169-7722(00)00156-X Kamra, S.K., Lennartz, B., Van Genuchten, M.T., et al., 2001. Evaluating non-equilibrium solute transport in small soil columns. Journal of Contaminant Hydrology, 48: 189-212. doi: 10.1016/S0169-7722(00)00156-X Li, B.Y., Zhu, Z.H., Zhang, A.X. et al., 1998. The experimental study on the transport of Na+ and Ca++ in unsaturated soil. In: Proceeding of second international conference on unsaturated soils. International Academic Publishers, Beijing, 575-580. Li, B.Y., Zhu, Z.H., Zhang, A.X. et al., 1998. The experimental study on the transport of Na+ and Ca++ in unsaturated soil. In: Proceeding of second international conference on unsaturated soils. International Academic Publishers, Beijing, 575-580. Li, L., Barry, D.A., Morris, J.C., et al., 1999. CXTANNEAL: An improved program for estimating solute transport parameters. Environmental Modelling & Software, 14: 607-611. Li, L., Barry, D.A., Morris, J.C., et al., 1999. CXTANNEAL: An improved program for estimating solute transport parameters. Environmental Modelling & Software, 14: 607-611. Liu, Q.S., 1989. Numerical methods for least square problems. Beijing University of Technology Press, Beijing(in Chinese). Liu, Q.S., 1989. Numerical methods for least square problems. Beijing University of Technology Press, Beijing(in Chinese). Nielsen, D.R., Van Genuchten, M.T., Biggar, J.W., 1986. Water flow and solute transport process in the unsaturated zone. Water Resources Research, 22(9): 895-1085. Nielsen, D.R., Van Genuchten, M.T., Biggar, J.W., 1986. Water flow and solute transport process in the unsaturated zone. Water Resources Research, 22(9): 895-1085. Pang, L.P., Murray, E.C., 1999. Non-equilibrium transport of Cd in alluvial gravels. Journal of Contaminant Hydrology, 36: 185-206. doi: 10.1016/S0169-7722(98)00110-7 Pang, L.P., Murray, E.C., 1999. Non-equilibrium transport of Cd in alluvial gravels. Journal of Contaminant Hydrology, 36: 185-206. doi: 10.1016/S0169-7722(98)00110-7 Parker, J.C., Van Genuchten, M.T., 1984. Determining transport parameters from laboratory and field tracer experiments. Vir. Agric. Ex. Stat. Bull., Blacksburg. Parker, J.C., Van Genuchten, M.T., 1984. Determining transport parameters from laboratory and field tracer experiments. Vir. Agric. Ex. Stat. Bull., Blacksburg. Sun, N.Z., 1994. Inverse problems in groundwater modeling. Kluwer Academic Publisher, Boston. Sun, N.Z., 1994. Inverse problems in groundwater modeling. Kluwer Academic Publisher, Boston. Tikhonov, A., Arsenin, V., 1977. Solutions of ill-posed problems. John Wiley & Sons, Inc., Washington. Tikhonov, A., Arsenin, V., 1977. Solutions of ill-posed problems. John Wiley & Sons, Inc., Washington. Torsten, Z., 1998. Capability of convection-dispersion transport models to predict transient water and solute movement in undisturbed soil columns. Journal of Contaminant Hydrology, 30: 101-128. doi: 10.1016/S0169-7722(97)00034-X Torsten, Z., 1998. Capability of convection-dispersion transport models to predict transient water and solute movement in undisturbed soil columns. Journal of Contaminant Hydrology, 30: 101-128. doi: 10.1016/S0169-7722(97)00034-X Wang, C., 1996. Numerical simulation for water and non-conservative contaminant transport in soils. Journal of Hohai University, 24(3): 5-11(in Chinese with English abstract). Wang, C., 1996. Numerical simulation for water and non-conservative contaminant transport in soils. Journal of Hohai University, 24(3): 5-11(in Chinese with English abstract). Wang, Y., 1996. Homotopy optimization algorithm via computer. Dalian Marine University Press, Dalian(in Chinese). Wang, Y., 1996. Homotopy optimization algorithm via computer. Dalian Marine University Press, Dalian(in Chinese). 刘钦圣, 1989. 最小二乘问题的计算方法. 北京: 北京工业大学出版社. 刘钦圣, 1989. 最小二乘问题的计算方法. 北京: 北京工业大学出版社. 王超, 1996. 水分及非保守污染物在土壤中运移数值模拟. 河海大学学报, 24(3): 5-11. doi: 10.3321/j.issn:1000-1980.1996.03.002 王超, 1996. 水分及非保守污染物在土壤中运移数值模拟. 河海大学学报, 24(3): 5-11. doi: 10.3321/j.issn:1000-1980.1996.03.002 王宇, 1996. 计算机优化同伦算法. 大连: 大连海洋大学出版社. 王宇, 1996. 计算机优化同伦算法. 大连: 大连海洋大学出版社. -
下载:
下载:
点击查看大图
计量
- 文章访问数: 3978
- HTML全文浏览量: 426
- PDF下载量: 9
- 被引次数: 0
