摘要：模型参数优化是通过极小化目标函数使得模型输出和实际观测数据之间达到最佳的拟合程度. 由于环境模型本身的复杂性，常规优化算法难以达到参数空间上的全局最优. 近年来，随着计算机运算效率的快速提高，直接优化方法得到了进一步开发与广泛应用. 本文比较了CRS、SCE UA、SA 和Annealing2Simplex 等4 种算法应用于环境模型参数优化的结果和计算效率。
关键词：参数优化 环境模型 CRS 算法 SCE UA 算法 Simulated2Annealing 算法 Annealing2Simplex 算法
人们对环境系统的深入研究是建立在环境模型的广泛应用基础上的. 为了更加精确地刻画环境系统的行为，环境模型在近10 年里表现出了强烈的复杂化趋势;不同空间尺度、不同时间过程模型的耦合，进一步加剧了这一过程. 环境模型的复杂性导致了模型结构和参数可识别性问题的提出，并成为当今环境建模理论研究的热点[1 ，2 ] . 其中在不确定性的框架下，模型参数的优化是研究的一个重要方面.
解决优化问题的难度主要取决于模型参数的空间维数和模型本身的非线性特征. 一般来说，参数越多、非线性越强，优化时间和精度就越差，同时也越不能够保证优化算法是否收敛到整体最优. 传统经验表明，求解优化问题的困难主要体现为[4 ，5 ] : ①全局搜索可能收敛到多个不同的吸引域; ②每一个吸引域可能包含一个或多个局部最小值; ③目标函数在n 维参数空间上不连续; ④参数及相互间存在高度灵敏性和显著非线性干扰; ⑤在最优解的附近，目标函数往往不具有凸性。
优化算法可以分为直接算法和间接算法2大类. 间接算法(如牛顿法以及各种以牛顿法为基础的改进算法) 的局限性主要在于要求目标函数在相关值域上必须是可微的;而直接算法仅涉及目标函数值的计算，不需要计算目标函数的导数. 因此尽管后者的计算效率相对较低，但在环境问题的实际应用中，它可以有效和简洁地解决由于模型复杂性所衍生的不可微函数的优化问题. 同时，与求解问题所耗费的时间相比，通常更为关注解的结果能够在多大程度上描述系统的行为[1 ， 2 ] . 所有这些使得直接算法在近年来得到了迅速发展和广泛的应用.
本文的目的在于分析和比较几种近年来渐为接受的参数直接优化算法的计算效率、计算精度及其算法稳定性. 由于直接算法本质上的随机性，以及对于不同优化问题所表现出的算法特性上的差异，因此本文仅以经典测试函数和环境水文模型为例，对几种算法进行详细的比较研究.
1 参数优化算法
模型参数的优化即是寻求一组参数，使模型的输出与实际观测数据之间按给定目标函数的度量方式达到最佳拟合，即:
Min f ( ys ， yobv ，θ) = f ( ys ， yobv ，θ3 ) θ ∈ S (1)
式中， f ( y ， θ) 为目标函数; ys 表示模型输出变量， yobv为系统观测值;θ3表示参数可行域S 上的最优参数向量.
最早提出的直接算法均是简单随机方法，这些算法除了计算效率低外，主要缺陷在于要求目标函数在邻域上必须是不相关的.
控制随机搜索算法(CRS) [10 ， 11 ] 有效地克服了简单随机算法的主要缺点，在计算过程中保存指定数目的参数及其对应的目标函数值，引入几何学中“重心”的概念，即考虑了新点产生的随机性，又在一定程度上保证了搜索的整体性. 复合形混合演化算法( SCE UA) [4 ， 5 ]所采用的竞争演化和复合形混合的概念继承了CRS 算法中全局搜索和复合形演化的思想，该方法将生物自然演化过程引入到数值计算中，模拟了生物进化的过程，提高了计算效率和全局搜索整体最优的能力. 模拟退火算法(SA) [7 ， 12 ]则假设优化问题的解及其目标函数分别与固体物质的微观状态及其能量所对应，采用蒙特卡洛(Monte Carlo) 随机方法模拟固体稳定“退火”的过程，并假设优化过程中递减目标函数值的控制参数t 与“退火”过程中的温度T 所对应，对于控制参数t 的每一个值，算法持
续进行“产生新解2判断2接受/ 舍弃”的迭代过程，应用该算法的关键在于确定合理的冷却进度表. 退火单纯形算法(AS) [8 ， 9 ]综合了下山单纯形方法和模拟退火法2 种优化算法，充分利用单纯形的形变信息，以一定温度T 所对应的概率接受准则来指导单纯形的映射、压缩和扩展等过程，提高了计算效率和算法稳定性.
运用直接算法进行参数优化首先要正确选取算法的内部控制参数. 一般根据具体问题的特点，采取数值实验的方法不断测试参数“性能”，从而确定既保证算法收敛到满意的极值，又不至于使计算过于耗费时间而导致算法失效. 选择控制参数本身实际就是一个“优化”过程.
2 测试函数的优化
首先采用一个广泛使用的、具有多个局部极小值的经典函数，对上述几种优化算法的全局搜索能力进行较为全面的测试，其函数形式为:
尽管测试函数仅有x 、y 2 个变量，但其等值线却具有复杂的空间结构，即多个局部极小、值域空间上的不连续和全局最小值周围的香蕉状狭长低谷(图1) . 现已知参数全局最小值为(5 ，5) ，所对应的目标函数值为10.0。
如前所述，为了保证直接算法的优化精度和避免算法失效，在算法的运行中需要在如下３个优化终止准则中进行选择和组合: ①前m组最佳估计参数收敛到可接受的参数空间，即最优参数的不确定性小于ε1 ; ②最佳目标函数值的优化改进速率小于给定值ε2 ; ③搜索过程中参数采样总数不超过给定的计算次数n.
表1 给出了4 种直接算法分别基于终止准则1 和2 ，以及不同初值假设条件下的测试函数参数最优值和计算效率. 为了尽可能消除算法随机性对于算法比较的影响，表格中数据除了计算次数的偏差外，均为相同条件下5 次运算结果的算术平均值.
由于AS 算法中的下山单纯形方法具有很强的方向性，因而在两个终止准则约束下，算法均可迅速达到全局最优化，且结果具有较高精度. 尽管SCE UA 算法引入竞争演化思想后提高了空间的搜索效率，但是不同种群从不同方向上向全局最优点逼近，以及同一种群内部仍然采用的单纯形控制随机搜索方法，使得该算法在计算效率上远不及AS 算法. 然而SCEUA 算法在各种初值条件下均可得到十分精确的全局优化结果，可见该算法的稳定性和可靠性比AS 算法具有更大的优势. 与前2 种算法相比，CRS 算法和SA 算法在计算效率和精度方面均表现不佳. 特别是当参数初值发生变化时，必须相应地调整控制参数才能得到较好的计算结果，这说明2 种算法的适应能力较差.
|
优化算法 |
编号 |
假设初值 |
参数优化结果 |
目标函数值 |
计算次数 | |||||||
|
准则1 |
准则2 | |||||||||||
|
x |
y |
x |
y |
x |
y |
准则1 |
准则2 |
准则1 |
准则2 |
偏差/ % | ||
|
CRS |
1 |
0.40 |
0.60 |
0.5000 |
0.5001 |
0.5000 |
0.5001 |
10.0001 |
10.0001 |
565 |
584 |
-3.25 |
|
2 |
0.10 |
0.20 |
0.5000 |
0.4999 |
0.5000 |
0.5000 |
10.0001 |
10.0000 |
806 |
841 |
-4.16 | |
|
3 |
0.20 |
0.80 |
0.5001 |
0.5001 |
0.5000 |
0.5000 |
10.0002 |
10.0000 |
585 |
618 |
-5.34 | |
|
4 |
0.70 |
0.70 |
0.4998 |
0.5001 |
0.5000 |
0.4999 |
10.0002 |
10.0000 |
541 |
571 |
-5.25 | |
|
5 |
0.80 |
0.85 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
595 |
583 |
2.06 | |
|
SCE UA |
1 |
0.40 |
0.60 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
479 |
356 |
34.55 |
|
2 |
0.10 |
0.20 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
494 |
350 |
41.14 | |
|
3 |
0.20 |
0.80 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
501 |
359 |
39.55 | |
|
4 |
0.70 |
0.70 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
431 |
305 |
41.31 | |
|
5 |
0.80 |
0.85 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
480 |
335 |
43.28 | |
|
SA |
1 |
0.40 |
0.60 |
0.5001 |
0.5002 |
0.5000 |
0.5000 |
10.0013 |
10.0000 |
230 |
349 |
-34.10 |
|
2 |
0.10 |
0.20 |
0.5041 |
0.4973 |
0.5000 |
0.5000 |
10.1300 |
10.0000 |
408 |
504 |
-19.05 | |
|
3 |
0.20 |
0.80 |
0.4999 |
0.5002 |
0.5000 |
0.4999 |
10.0002 |
10.0000 |
468 |
579 |
-19.17 | |
|
4 |
0.70 |
0.70 |
0.5002 |
0.4999 |
0.5000 |
0.4999 |
10.0005 |
10.0001 |
468 |
582 |
-19.59 | |
|
5 |
0.80 |
0.85 |
0.5001 |
0.5000 |
0.5001 |
0.5000 |
10.0000 |
10.0000 |
490 |
544 |
-9.93 | |
|
AS |
1 |
0.40 |
0.60 |
0.4999 |
0.5001 |
0.4999 |
0.5001 |
10.0000 |
10.0000 |
73 |
63 |
15.87 |
|
2 |
0.10 |
0.20 |
0.5000 |
0.5001 |
0.5000 |
0.5001 |
10.0001 |
10.0001 |
150 |
242 |
-38.02 | |
|
3 |
0.20 |
0.80 |
0.5000 |
0.5000 |
0.5000 |
0.5000 |
10.0000 |
10.0000 |
108 |
135 |
-20.00 | |
|
4 |
0.70 |
0.70 |
0.5001 |
0.4999 |
0.5001 |
0.4999 |
10.0000 |
10.0000 |
99 |
191 |
-48.17 | |
|
5 |
0.80 |
0.85 |
0.4999 |
0.5001 |
0.5001 |
0.4999 |
10.0000 |
10.0001 |
157 |
157 |
0.00 | |
1) 计算次数偏差的计算公式为( n准则1 - n准则2) / n准则2
|
计算次数的统计量 |
CRS |
SCE- UA |
SA |
AS |
|
准则1的标准方差 |
106.89 |
27.36 |
106.65 |
35.46 |
|
的标准方差 |
114.05 |
22. 15 |
96.24 |
66.52 |
|
准则1和准则2偏差均值/% |
-3.19 |
39.97 |
–20.37 |
–18.06 |
可以采用各种初始条件下对应同一终止准则计算次数的标准方差来表征算法的随机性.尽管基于控制性随机搜索思想构造的直接优化算法不可能完全消除计算结果的随机性，但是算法结构的精心设计将在一定程度上减小随机性影响，同时也就意味着算法稳定性的增大. 表2 中结果显示SCE UA 算法的计算次数在各种初始条件下的标准方差远小于其他3 种直接算法，再次表明了算法所具有的良好稳定性. 这一结果的产生是由于SCE UA 算法结构中引入的复合形概念及其竞争混合算法大大削弱了计算结果的随机性.
对应准则1 和2 的计算次数偏差，表征了算法对于选用不同终止准则的敏感性. 根据表2 结果可以得到4 种算法的敏感性排序为(从大到小) : SCE UA &> SA &> AS &> CRS ，特别是CRS 算法的计算次数偏差的最大值和统计平均值分别为- 5.34 %和- 3.19 %，远小于其他3种算法. 更进一步发现，算法结构复杂性与算法对于终止准则选用的敏感性之间高度相关，表现为本例中算法敏感性排序与复杂性顺序完全一致. 其原因在于，直接算法在结构复杂性所形成的内部约束与目标函数的外部约束的共同作用下趋向全局最优，其结果必然是一个不断协调和折衷的过程. 算法越复杂，内部约束就越强，其计算结果就越发表现出与终止准则之间的高度相关性. 因此简单算法对于外部终止准则控制性约束的适应能力较之复杂算法更为灵活.
3 环境模型参数优化:水文模型实例
随着环境模型的不断开发和广泛应用，环境模型的种类和数量日益丰富，模型本身所表现出的结构特征也日趋复杂. 因此研究和比较直接算法的参数优化特性，既无可能也无必要将其应用于每一个环境模型的参数优化评估中. 现实而合理的方法是选择具有典型结构复杂性和非线性特征的环境模型进行实例研究.由于篇幅所限，本研究仅以箱式水文模型为例，对4 种参数直接优化算法进行比较研究.
研究中采用2 个箱式的模型结构模拟径流在流域中的运动情况(例如，Cooper 等人[3 ]在其研究中给出的水文箱式模型介绍)，该模型共有11 个参数，表3 给出了本文中模型参数的取值范围及其物理意义.
|
模型参数 |
最大取值 |
最小取值 |
说明 |
|
k1 |
0.05 |
0.3 |
不同箱体出口处的径流流速/ h-1 |
|
k2 |
0.008 |
0.1 | |
|
k3 |
0.01 |
0.1 | |
|
k4 |
0.005 |
0.1 | |
|
k5 |
0.0005 |
0.002 | |
|
h1 |
5.0 |
20.0 |
不同箱体出口处的高度/ mm |
|
h3 |
1.0 |
10.0 | |
|
h3 |
700.0 |
1500.0 | |
|
h4 |
500 |
1000.0 | |
|
H1 |
0.1 |
20.0 |
上层和下层高度/ mm |
|
h3 |
5.0 |
200.0 |
根据不同最优化准则构造的目标函数将对参数优化结果会产生一定影响[2 ， 6 ] . 本研究以最小二乘法为优化准则构造的目标函数如下:
其中， f obv为观测流量， f sim表示模拟输出流量; N为模拟点数. 本研究取N = 168 ，模拟暴雨期7 天的小时观测流量. 终止准则采用准则1 和3.
图2 是水文模型参数优化后的流量模拟曲线，可以直观的看到4 条优化曲线均较好的拟合了观测流量. 表4 给出了4 种直接算法分别在相同条件下运算5 次，以不确定性表示的最优参数值和最小二乘的目标函数值， 以及以CPU 时间表示的平均计算效率.
目标函数值和最佳估计参数的统计特征值共同表征了优化算法的稳定性和可靠性. 根据表4 中的计算结果， SA 算法和CRS 算法的目标函数值和最佳估计参数均具有较大不确定性，尽管2 种算法分别得到了相对较小的目标函数值( ZSA = 117774 ， ZCRS = 118648) ，但是目标函数的均值和标准方差2 个统计量却均较大. 相反，AS算法和SCE UA算法则表现出更为优秀的计算特性，即算法既能够保证一定的计算精度，又具有可靠的重复性. 与SCE UA算法相比，AS 算法的目标函数值和最佳估计参数的不确定性更为缩小，进一步说明AS 算法全局参数寻优过程所具备的计算稳定性和可靠性.
|
|
项 目 |
CRS 算法 |
SCE UA 算法 |
SA 算法 |
AS 算法 |
|
模型参数 |
k1/h-1 |
0.1890～0.2272 |
0.2114～0.2166 |
0.1920～0.2199 |
0.2133～0.2136 |
|
k2/h-1 |
0.0396～0.0453 |
0.0416～0.0433 |
0.0354～0.0454 |
0.0433～0.0434 | |
|
k3/h-1 |
0.0106～0.0227 |
0.0101～0.0116 |
0.0100～0.0215 |
0.0100～0.0100 | |
|
k4/h-1 |
0.0100～0.0789 |
0.0292～0.0639 |
0.0105～0.0679 |
0.0051～0.0988 | |
|
k5/h-1 |
0.0010～0.0019 |
0.0010～0.0013 |
0.0009～0.0020 |
0.0008～0.0010 | |
|
h1/mm |
18.303～19.121 |
17.700～19.970 |
18.109～19.993 |
18.870～19.945 | |
|
h2/mm |
1.6035～7.9821 |
1.5460～4.3554 |
1.5281～9.3075 |
2.8560～4.1662 | |
|
h3/mm |
930.99～1404.2 |
927.73～1307.4 |
701.88～1383.2 |
856.38～1447.9 | |
|
h4/mm |
603.94～948.72 |
545.10～855.93 |
505.29～977.14 |
634.18～993.52 | |
|
H1/mm |
2.2234～7.2466 |
6.0477～8.5449 |
1.3114～8.0765 |
7.3351～8.4248 | |
|
H2/mm |
96.759～192.79 |
120.84～171.94 |
156.24～199.49 |
159.50～199.89 | |
|
目标函数值 |
函数值范围 |
1.8648～1.9841 |
1.9377～1.9453 |
1.7774～2.1551 |
1.9330～1.9343 |
|
均值 |
1.9537 |
1.9405 |
1.9469 |
1.9335 | |
|
标准方差 |
0.0504 |
0.0029 |
0.1347 |
0.0005 | |
|
平均CPU 时间(5000 次) |
2.7085e + 004 |
2.8219e + 004 |
2.5927e + 004 |
1.9965e + 004 | |
对应于一定退火温度下的接受准则，AS 算法中的下山单纯形法在目标函数下降的方向上搜索全局最优，不仅有效避免了算法陷入局部极小值，而且大大提高了计算效率. 表4 的计算结果表明AS 算法的计算效率远高于其他3 种直接算法. 由于SCE UA 的算法设计是在CRS算法基础上引入生态学竞争演化模式，尽管全局优化的可靠性和稳定性有较大提高，但却是以算法复杂性的增加和计算效率的下降为代价.
此外，研究中发现SA 算法引入的物理结晶过程使得算法本身控制参数的选择具有较大困难. 对于不同优化问题，初始温度、退火率和马尔可夫链长度所组成的冷却进度表有很大差异. 使用者必须通过数值实验精确调整冷却进度表使之能够适应优化问题，这就降低了SA 算法的适应性并增加了实际应用过程中的困难.
实例研究表明，环境模型的参数优化过程远较测试函数复杂，4 种直接算法的最佳估计参数和目标函数无一收敛到单一数值，并且采用统计特征量描述优化结果的不确定性也具有较大差异. 直接算法优化结果的不确定性分析不仅为进一步比较和理解算法特性提供了有效途径，更为重要的是其揭示了优化算法不能完 全解决复杂模型参数的识别问题，因此基于环境模型结构和参数的不确定性研究成为必然.
4 结论
本文通过2 个模型实例的参数优化计算，对4 种直接算法的基本优化特性进行了分析比较，研究表明AS 算法和SCE UA 算法具有较高的可靠性和稳定性，AS 算法同时还具有较高的计算效率.
理论上，对于任何在整个参数可行域上进行随机搜索的优化算法，当计算次数足够大并适当选取了控制参数时，算法都能够收敛到唯一的全局最优值. 然而实践远较之复杂:不仅需要在计算精度和计算效率之间做出艰难的选择，更必须面对模型复杂性增加所导致的参数可识别问题. 由于算法随机性的影响，特别是模型参数具有的不确定性，4 种直接算法得到的最佳参数估计值与目标函数无法较好吻合. 由此可见，基于参数优化算法框架下的模型参数仍然不可识别，即优化算法不能解决模型参数的不确定性问题.
转贴于参考文献
1 Beck M B. Water quality modeling : a review of the analysis of uncertainty. Wat . Res. Res. ， 1987 ， 23 ( 8) : 1393 ～1442.
2 Chen J ， Wheater H S. Identification and uncertainty analysis of soil water retention models using lysimeter data. Wat .Res. Res. ， 1999 ， 35 (8) : 2401～2414.
3 Cooper V A et al. Evaluation of global optimization methods for conceptual rainfall2runoff model calibra tion. Wat . Sci.Tech. ， 1997 ， 36 (5) 53～60.
4 Duan Q Y et al. Shuffled complex evolution approach for ef2 fective and efficient global minimization. J . Optimization Theory and Applications ， 1993 ， 76 (3) : 501～521.
5 Duan Q Y et al. Optimal use of the SCE UA global opti2 mization method for calibrating watershed models. J . Hy－drol. ， 1994 ， 158 : 265～284.
6 Freedman V L et al. Parameter identifiability for catchment－scale erosion modelling : a comparison of optimization algo－rithms. J . Hydrol. ， 1998 ， 207 : 83～97.
7 康立山等. 非数值并行算法(第1 册) ———模拟退火算法.北京:科学出版社，1994. 32～128.
8 Pan L H ， Wu L S. A hybrid global optimization method for inverse estimation of hydraulic parameters : Annealing－Simplex method. WR972418 ， Department of Soil && Environ－mental Science ， University of California ， Riverside. 1997. 87～196.
9 Press W H et al. Numerical Recipes in Fortran. Second 2 tion. New York : Cambridge Univ. Press ， 1992. 6～97.
10 Price W L. A controlled random search procedure for global optimization. The Comput . J . ， 1977 ， 20 : 367～370.
11 Price WL. Global optimization by controlled random search.J . of Optimization Theory and Applications ， 1983 ， 40 : 333～348.
12 姚姚. 蒙特卡洛非线性反演方法及应用. 北京:冶金工业出版社， 1997. 43～212.