吸附性土壤溶质运移参数识别的粒子群-差分算法

K s 为饱和导水率 （ cm min － 1 ） ， D w 为自由水中 式中，
2 －1 7 /3 2 D0 、 m、 分子扩散系数 （ cm min ） ， = θ / θ s ， β为 经验常数。
（ mg L － 1 ） 。
2
2. 1
待求参数计算方法
优化模型
由于土壤含水率和溶质浓度单位量级不同， 故采 用归一化法将含水率和浓度转化为 0 ～ 1 区间上无量 记归一化后第 i 点第 j 时刻含水率和溶质浓 纲数值， ［18 ］ C sij ， 然后采用权重系数法 将式 （ 5 ） 度分别为 θ sij 、 多目标优化问题， 转化为式（ 6 ） 单目标优化模型。
Cxtfit 方法， 利用 Gm王康等
［9 ］
(
(
)
)
(
)
－
（ cJ w ） x
（ 1）
数值法求解此类问题具有一定的精度， 但也存在公 式推导过程繁琐复杂或者数值震荡等问题； Lou 和
* 国家自然科学基金项目（ 50979089 ， 51239009 ） 资助 通讯作者： 白
mail： baidan@ xaut. edu. cn 丹（ 1960 —） ， 男， 教授， 博士生导师， 主要从事节水灌溉理论与技术研究。E-
槡
∑ ∑ （ cij j =1 i =1
i 为采样位置点序号， n 为采样位置点总数， j 式中， h 为采样时间点总数， 为采样时间点序号， θ ij 为第 j
3 －3 c ij 为第 j 时刻 时刻第 i 点含水率实测值（ cm cm ） ， －1 D0 ， m） 第 i 点溶质浓度实测值 （ mg L ） ， θ ij （ K s ， β， 3 －3 c ij （ D0 ， 为第 j 时刻第 i 点含水率计算值（ cm cm ） ， m， Ks ， ， D ， k ） j i β w d 为第 时刻第 点溶质浓度计算值
第 50 卷 第 3 期 2013 年 5 月
土
壤
学
报
Vol. 50 ， No. 3 May， 2013
ACTA PEDOLOGICA SINICA
吸附性土壤溶质运移参数识别的粒子群 差分算法
任长江
1
*
白 丹
1
周蓓蓓
1
梁 伟
2
陈 燕
1
（ 1 西安理工大学水利水电学院， 西安 710048 ） （ 2 陕西师范大学旅游与环境学院， 西安 710062 ）
而土壤含水率 率以及浓度的实测值和计算值有关， D 、 m 、 K 、 、 D k d 这 6 个参 和溶质浓度计算值与 0 s β w、 数有关， 对此问题难以建立明确的函数关系式 ， 针 对这一复杂的优化问题采用有粒子群 有限差分 法进行求解， 粒子群算法参数识别具体计算过程见 图 1 所示。粒子群算法两个主要参数为粒子位置和 速度， 本例中是水分和溶质运移参数的函数 ， 粒子 位置和速度采用动态权重和异步时变学习因子进 行更新， 见式（ 7 ） ～ 式（ 13 ） 所示。 x i = rand（ D0 ， m， Ks ， Dw ， kd ） β， v i = rand（ D0 ， m， Ks ， Dw ， kd ） β， x i +1 = x i + v i （ 7） （ 8） （ 9）
487
3 －3 式中， θ 为体积含水率（ cm cm ） ， θa 为土壤初始含水
率（ cm cm ） ， θb 为土壤上边界含水率 （ cm cm 为入渗 溶 液 浓 度 （ mg L
－1
3
－3
3
－3
）， c
水率和溶质浓度实测值， 同时对式 （ 1 ） 水分和溶质 ， m、 Ks 、 Dw 、 kd 运移方程采用差分法 代入参数 D0 、 β、 便可获得和试验中对应时空分布条件下的土壤含 水率以及溶质浓度计算值。 分别以水分和溶质浓 度的实测值和计算值的标准差最小建立如下所示 多目标优化模型。
［11 ］ Li［10］、 Li 等［12］利用同伦正则化法和参数 崔凯等 、
土壤水分以及溶质运移方程参数的确定是土 壤水分和溶质运移数值模拟的基础。 目前对土壤 水分和溶 质 参 数 的 求 解 主 要 有： 借 助 现 代 测 量 设 模型 模 拟 软 件、 数 值 计 算、 数 值 反 演 等 方 法。 备、 利用高光谱影像技术通过对土壤光 Mohamed 和 Marwan［2］ 利 用 TDR （ Time 谱的 分 析， Domain Reflectometry） 对土壤介电常数的测定， Can［3 ］ demir 等 通过对土壤物理特性的测定等方法间接 地计算出了土壤水分运移相关参数。 现代测量设 备能够快速测量水分以及溶质运移参数 ， 但高昂的 设备 费 用 使 其 应 用 受 到 限 制； 刘 建 军 等 利 用 Hydrus1D 软 件， Zhang［5］ 采 用 Mike2D 软 件 包， 张 ［6 ］ 嘉和王明玉 利用 MT3DMS 软 件 分 别 对 一 维、 二 维、 三维等水分、 溶质运移参数进行了求解。 模拟 软件能够快速的对研究问题给出模拟结果 ， 但模拟 的精度 易 受 软 件 运 行 参 数 设 置 的 影 响。 Duan 和 Cui［7］基于小波函数用 WaveletGalerkin 法， 杨文涛 等 采用粒子平均移 动和分形理论等数值方法计算土壤溶质运移参数 。
［13 ］
Kandasamy 等
［1 ］
1
非 饱 和 水 分溶质运移方程及其 参数
定水头定浓度条件下一维垂直非饱和吸附性 ［15 ］ 水分溶质运移方程 ： θ θ K （ θ） t = x D（ θ） x － k r 1 + d b （ θc） = D Ih c θ t x x c = ca （ t = 0， x） θ = θa ， θ = θ ， （ t， x = 0） b c = cb θ = θ ， （ t = ∞， x） a c = ca
［8 ］ ［4 ］
， Sun H Y 等［14］ 采用遗传 算法对土壤水分参数进行了数值反演， 数值反演在 迭代优化算法， 杨坤等 解决此类问题中应用较广， 但也存在算法不收敛或 者陷入局部最优解等问题。 因此寻求设备简单、 方 计算精确、 算法收敛的计算水分和溶质运 法简便、 移参数方法是本文要研究的内容 。
mail： changjiang412032@ 163. com 作者简介： 任长江（ 1984 —） ， 男， 博士研究生， 主要从事节水灌溉理论与技术研究。E收稿日期： 2012 － 07 － 26 ； 收到修改稿日期： 2013 － 01 － 03
3期
任长江等： 吸附性土壤溶质运移参数识别的粒子群 差分算法
h h
）， ca 为 土 壤 初 始 浓 度
（ mg L － 1 ） ， c b 为土壤上边界浓度 （ mg L － 1 ） ， D（ θ） 为
2 －1 K （ θ ） 为非饱和导水率 非饱和扩散率 （ cm min ） ， －1 （ cm min ） ， DIh （ θ） 为水动力弥散系数（ cm2 min － 1 ） ，
k d 为吸附系数， v 为平均孔隙流速 （ cm min － 1 ） ， rb 为 J w 为水分通量 （ cm min 土壤容重 （ g cm ） ， t 为入渗时间（ min） 。 为入渗距离（ cm） ， 和非饱和导水率 ［17 ］ 力弥散系数 计算公式如下： D （ θ ） = D0 （ θ / θ b ） K （ θ） = K s e － β（ θb － θ） D Ih 因子， c1 i 、 c1 f 、 c2 i 、 c2 f 为 c1 式中， gbest i 为粒子历史最 和 c2 的初始值和最终值； pbest i 、 w i 为第 i 代惯性权重因子， 优位置和全局最优位置， w max 、 w min 为 最 大 和 最 小 惯 性 权 重 因 子， n 为迭代 次数。 待求参数优化计算过程实际上是联合运用一 维非饱和水分运移方程和溶质运移方程以及待求 参数计算的多目标数学模型进行的， 对一维非饱和 土壤水分 － 溶质运移方程采用差分法求解， 而待求 参数计算的多目标数学模型则采用粒子群算法求 解， 最终获得待求参数的最优值。
对于一维垂直吸附性非饱和水分 － 溶质运移 通过试验可获得不同时间不同位置上土壤含 问题，
minZ = （ 1 － α） minZ1 + αminZ2 h h h h 2 （ 6） = （ 1 － α） （ － （ K ， ， D ， m ） ） / n + （ csij － csij （ D0 ， m， Ks ， Dw ， kd ） ） 2 / n θ θ β α β， ∑ ∑ ∑ ∑ sij sij s 0 1 1 1 1 j = i = j = i = θ = （ θ － θ ） / （ θ － θ ） ； c = （ c － c ） / （ c － c ） sij max ij max min sij max ij max min vi +1 = wi vi + c1 rand（ pbest i － xi ） + c2 rand（ gbest i － x1 ） 式中， α 为权重系数（ 0 ＜ α ＜ 1 ） 。
槡
槡
2. 2
粒子群算法求解模型步骤 参数识别优化计算模型式 （ 6 ） 是一个含有 6 个 参数的非线性规划模型， 其目标函数值与各点含水
（ 10 ） w i = w min + （ w max － w min ） i / n c1 = c1 i + （ c1 f － c1 i ） i / n c2 = c2 i + （ c2 f － c2 i ） i / n （ 11 ） （ 12 ） （ 13 ）
488
土
壤
学
报