基于EnKF法的径流数据同化对SWAT模型参数优化效果评估

刘永伟，王 文，刘元波，刘 庆，凌 哲

(1.中国科学院南京地理与湖泊研究所流域地理学重点实验室，江苏 南京 210008；2.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；3.江西省水利规划设计研究院，江西 南昌 330029； 4.江苏省水利工程规划办公室，江苏 南京 210029)

数据同化是将多源观测信息应用于模型模拟预报的一种有效手段，被公认是提高模型模拟预报精度、量化不确定性方面最具前途的方法[1-2]。径流是流域产汇流模拟中最为重要的输出变量，可以反映整个流域或汇水区域综合的水文信息，且地面实测径流数据具有很高的可靠性。因此，将站点径流实测信息同化到水文模型模拟预报过程中，对于改善流域状态变量及模型参数的估计精度，进而改进流域的降雨径流过程模拟及预测具有重要意义。

径流的同化始于20世纪80年代，主要是在流域/区域的水文模拟预报中基于站点径流观测数据进行模型输出结果的修正，也称为误差校正[3]。当时所使用的同化方法主要为卡尔曼滤波(KF)法和扩展KF法。随着数据同化技术与模型模拟技术的发展，径流的同化不仅体现在线性系统或弱非线性系统对系统状态变量以及输出结果的校正上，更表现为在非线性系统中通过对系统状态变量及参数的更新或优化来改进整个系统的过程模拟及预报精度[4-9]。用于同化的水文模型由简单的集总式概念性模型(如新安江模型[10]、格林-安普特模型[11])逐渐发展为基于物理机制的复杂的分布式模型(如SWAT模型[12]、TopNet模型[13])，同化方法也由线性KF法逐渐发展为各种集合KF法(如EnKF法)、粒子滤波和变分法等。基于简单的集总式概念性模型的径流数据同化方面，大量研究表明，通过径流数据同化能够获得较好的模型参数与状态估计结果，从而有效地改进流域水文模拟的预报精度和可靠性[8，10，14]。但是，基于物理机制的分布式水文模型，由于其复杂的模型结构、大量的模型变量与参数，导致其径流数据同化过程中存在高维度现象，大大增加了数据同化效果的不确定性及同化实施的难度[12-13，15-16]。分布式水文模型中，径流数据的同化效果更容易受到模型的物理机制、模型结构、模型对流域空间异质性的概化能力以及同化系统规则等方面的影响。目前，分布式水文模型中径流数据的同化研究，基于径流观测数据进行模型状态更新与参数估计方面都相对滞后，主要表现在对分布式水文模型模拟误差与径流观测误差的合理量化[7，13，17]、对模型参数的同化处理(包括参数在同化过程中的演变以及参数过拟合问题的处理等)[5，15]等方面。

因此，本文在分布式水文模型SWAT中开展基于集合卡尔曼滤波(EnKF)法的径流数据同化研究，在合理量化模型模拟误差和径流观测误差、有效处理参数演变及过拟合问题的基础上，构建基于数据同化的水文模型参数优化方案，综合评估径流数据同化对SWAT模型参数的优化效果。

1 基于EnKF法和SWAT模型的数据同化方案

1.1 SWAT模型

SWAT模型是由美国农业部农业研究中心(ARG-USDA)开发的分布式流域尺度水文模型[18]，主要包含水文过程子模型、土壤侵蚀子模型和污染负荷子模型。SWAT模型水文过程子模型涉及的关键模块有地表径流模块、蒸散发模块、土壤水分运动模块、地下径流模块和河道汇流模块。地表径流的模拟一般采用SCS曲线法；蒸散发的模拟是在潜在蒸散发量的基础上计算实际蒸散发量，包括植被冠层截留蒸发、植被蒸腾和土壤蒸发；土壤水分运动采用垂直分层水量平衡模拟，考虑下渗与壤中流过程；地下径流的模拟分为浅层和深层地下径流，涉及潜水再蒸发、灌溉等过程，基于水量平衡方程实现；河道汇流基于变动蓄量法或马斯经跟法对进入河道中的水流进行汇流演算。

1.2 基于EnKF法的状态-参数优化更新

1.2.1 基于EnKF法的状态更新

EnKF法是Evensen[19]于1994年基于随机动态预报理论，将集合预报思想引入KF法，提出的一种顺序数据同化算法。EnKF法主要分为两步：预报(状态转移过程)和分析(状态更新过程)。基于Monte Carlo方法，将模型状态变量加入初始状态误差(通常为高斯白噪声)，生成一组(数目N)随机变量集合Xk，i(k=0，i=1，2，…，N)，每组随机变量Xk，i在k时刻的一步预报方程为

(1)

当k+1时刻有观测存在时，模型状态变量的预报值基于观测信息更新如下：

(2)

1.2.2 状态扩展、参数演变与过拟合问题处理

顺序数据同化中参数估计可以通过状态扩展方式[20]实现，即将待估计的参数在其物理区间内均匀采样后，获得待估计的参数集合，将其纳入状态向量X中，并基于观测信息与更新状态变量一样对估计参数进行更新[21]。假设待优化估计的参数向量为δ，则扩展后的状态向量可以表示为[XTδT]T，基于EnKF法，k时刻状态与参数的更新方程如下：

(3)

H*=[H0]

(4)

式中0为零矩阵，其维数为nob×npar，nob为观测向量Yt中观测变量的个数，npar为δ中参数的个数。

在以上基于状态扩展方式的EnKF法参数估计框架下，随着EnKF法更新过程的继续，参数集合将快速收敛，使得分析场很快趋近于背景场，最终导致观测信息对参数更新不起作用，即出现滤波发散。为避免滤波发散问题，模型参数在同化过程中需要与状态变量一样具有一个随时间的演变过程(状态变量随时间演变过程通过模型本身的物理机制实现)。但是，这与一般意义上对模型参数的理解(随机动态模型中的参数通常被假定为常数)有冲突，因此，通常期望参数的变化相比于状态变量的变化要小得多[1]。采用条件协方差法，模拟参数演变：

(5)

另外，参数更新过程中，一个重要的问题是参数的过拟合(过度更新)问题，其有可能导致模型模拟/预报过程中断或模型无法正常运行。参数的过度更新，可能是由于数据同化过程中异常信号的加入(如强降雨)或对模型模拟与观测误差的量化不合理，造成数据同化过程中对模型输出结果与实测数据的过度拟合导致。为抑制参数同化过程中的过度更新，采用参数平滑法，引入平滑因子β∈(0，1.0)进行处理[22]：

(6)

当β= 0时，表示参数更新过程中没有平滑作用，β越大，表示参数更新过程中的平滑作用越强。参考已有研究[22]，取β为0.8。

1.2.3 模型误差与径流观测误差的概化

SWAT模型产汇流模拟中，模型输入的不确定性主要考虑降雨和最高/最低气温输入误差。假定站点降雨输入误差服从对数正态分布[23]，且该降雨输入误差在时间和空间上独立：

(7)

(8)

式中：T为实测最高/最低气温；Te为误差扰动后最高/最低气温集合。

模型结构的不确定性对模型模拟预报的影响主要通过扰动模型状态变量来刻画，即对模型状态变量的预报集合加上一定的高斯误差扰动[15]：

(9)

模型参数的不确定性通过对敏感参数加以高斯误差扰动：

(10)

本文研究冗余机械臂液压驱动系统能量消耗优化方法,给出机械臂液压驱动简图模型,推导末端执行器运动方程式.根据机械臂末端执行器运动路径建立动力学模型,并且对运动参数进行离散化,构造多目标函数,添加约束条件.简化机械臂液压驱动动力学,引用动态规划算法优化液压驱动控制系统,在Matlab软件中对动态规划算法优化后的液压驱动系统进行能量仿真,并且与无动态规划的能量消耗仿真结果形成对比.结果表明:冗余机械臂液压系统采用动态规划算法优化后,能量消耗较小,机械臂运动相对稳定,有利于节约资源.

径流观测误差主要源于水位观测误差和水位流量转换误差。假定径流观测误差服从正态分布：

(11)

1.3 基于EnKF法的径流数据同化方案

数据同化方案主要包括模型算子、观测算子和同化算法。模型算子为SWAT模型，同化算法为EnKF法。观测算子用于连接需要优化的模型状态变量/参数并辅助同化的观测数据。辅助同化的观测数据即为水文站点的实测径流量，而需要优化的模型状态变量与参数的选择根据其对降雨径流过程模拟的重要性和敏感性，选择地表径流滞留量Qstor、土壤含水量SWly和河道蓄水量Vstored三类状态变量和控制地表径流、壤中流、地下径流以及河道汇流过程的7类参数，见表1。当状态变量与观测变量为同一类型的变量(如都为径流)时，观测算子为一个权重矩阵，其权重值可以通过状态变量到观测变量的插值过程获得。然而，由于SWAT模型的高度非线性，无法用一个矩阵显式地表示出SWAT模型中状态变量X到与观测变量Y所对应的模型输出变量的转化过程，为此，需要把Y对应的模型输出变量(即观测站点上的模拟径流)纳入X中。但在EnKF法更新过程中不对观测站点的模拟径流进行更新处理，因为只有对状态变量的更新才有意义[13]。此时，X中包括nst个模型状态变量和被更新的参数与nob个观测站点的模型模拟径流，这样X与Y中存在对应变量，二者的状态转移矩阵H*可以写成：

表1 SWAT模型径流数据同化中被更新的模型状态变量与参数

[0]nob×nstInob×nob

(12)

式中状态转移矩阵H*由nob×nst维的0矩阵和nob×nob维的单位矩阵构成。

基于EnKF法和SWAT模型的站点观测径流数据同化方案的实现是将EnKF法耦合到SWAT模型模拟过程中，该同化方案的流程如图1所示。

图1 基于EnKF法与SWAT模型耦合的站点观测径流数据同化方案流程Fig.1 Flow chart of the streamflow data assimilation scheme based on the coupling of EnKF and SWAT

2 研究区与数据

研究区为淮河上游淮滨水文站以上流域(简称淮滨流域)，流域控制面积为16 005 km2，整个流域除西部和西南部地区为山区外，大部分地区为平地，高程范围25～1 117 m。土地利用类型主要为农业用地(包括旱地34.2%和水田37.6%)和林地(23.6%)。土壤类型主要分为水稻土、黄褐土、黄棕壤、沙姜黑土等七类土壤。该流域冬春干旱少雨，夏秋闷热多雨，多年平均降水量约900 mm，其中有50%～80%发生在6—9月，多年平均气温约15℃。

SWAT模型的输入数据主要包括下垫面数据和气象数据两类(见表2)。降雨数据主要由流域内共106个雨量站提供。最高/最低气温、太阳辐射、风速和相对湿度数据从研究区内(信阳)及其附近(固始、广水、驻马店和枣阳)共5个气象站获得，来自中国气象科学数据共享网。另外，流域内有6个水文站，分别为大坡岭、长台关、竹竿辅、息县、潢川和淮滨，提供流域的逐日径流量数据(因潢川站流量数据质量有问题，故未采纳)。研究区气象、水文站点分布见图2，图中数字1～37为子流域代码。

表2 淮滨流域 SWAT模型构建中采用的数据

图2 淮滨流域气象、水文站站点分布及子流域划分 Fig.2 Map of the meteorological and hydrological stations，and the subbasin division in Huaibin River Basin

3 径流数据同化实现与结果分析

3.1 模型模拟与径流观测误差量化

EnKF法可以显式地处理模型与观测误差，并假定其都服从高斯分布。由EnKF法可以看出，对模型误差(模型输入、参数及结构)和径流观测误差的合理量化在一定程度上决定了数据同化的效果。本次基于径流数据同化的SWAT模型参数优化研究中，主要考虑降雨和最高/最低气温输入误差、模型状态误差、模型参数误差以及径流观测误差，误差均方差分别记为σp、σT、σs、σpar和σobs项。其中，模型状态误差又由模型状态变量误差和与观测变量对应的模型输出变量误差组成，分别记为σstate和σsflux。由于最高/最低气温输入误差对降雨径流模拟的影响较为敏感，借鉴以往研究[24]，设σT为1℃。而其他因子误差均方差的设置是在模型参数敏感性分析的基础上，采用自动率定与人工调整相结合的方法对模型参数进行调整后，基于调整后的模型采用最大后验估计(MAP)法并结合全局优化SCE-UA法[25]对其进行量化的。

在MAP法中，需要事先给定各个因子误差均方差的先验分布。借鉴已有研究[8]，各因子误差均方差的先验分布设置如下：①σobs的先验分布设置为均值9.25%、标准差0.925%的高斯分布；②其他因子σp、σstate、σsflux和σpar的先验分布设置为给定阈值内的均匀分布，其给定的阈值分别为0～0.5、0～0.3、0～0.5和0～0.3，另外σobs给定的阈值为0～0.5。确定各个因子误差均方差的先验分布后，即可获得抽样得到的某组误差均方差的先验概率。基于该组误差均方差对模型降雨输入、状态变量、参数以及输出变量进行扰动，得到率定期内的模型模拟径流输出系列集合。以流域出口控制站(淮滨站)的实测径流为验证数据，计算整个率定期内该组误差参数的似然函数。根据以上误差均方差的先验概率与似然函数计算得到其后验概率。通过SCE-UA方法搜索使后验概率达到最大时的误差参数组，作为σobs、σp、σstate、σsflux和σpar的最优估计，即最大后验估计。

对以上误差均方差的估计，采用了淮滨站2002—2004年共3年的日径流量数据。在SCE-UA方法中，模型截断运行次数为60 000次，而实际运行过程中，当模型运行到34 860次时已收敛，最终获得σobs、σp、σstate、σsflux和σpar的最优估计结果，分别为0.087、0.305、0.028、0.29、0.113。

图3为上述均方差参数条件下，2005—2008年的模型模拟径流(95%的置信区间)与站点实测径流过程，图中黑色实线为实测径流，灰色区域为模拟径流95%的置信区间。从图3可以看出，4年内绝大多数观测径流数据都落在95%的置信区间内，说明该误差参数估计具有一定的合理性。同时也注意到，低流量时(Q<100 m3/s)，径流模拟存在一定程度的低估，这与观测误差参数估算前模型的参数率定有一定关系。

图3 最优误差均方差条件下2005—2008年的模型模拟(95%置信区间)与站点实测径流过程Fig.3 Modeling(with 95% confidence interval) and observed runoff process based on the optimal error parameters over 2005—2008

3.2 模型参数同化更新过程

根据流域内5个站点(大坡岭、长台关、竹竿辅、息县和淮滨站)2003—2005年的日径流实测数据，基于EnKF法实时更新SWAT模型敏感参数与状态变量，以期获得模型参数的最优估计。从被同化更新的7类共59个敏感参数中(表1)，随机选取2个代表性参数，即水文条件Ⅱ下SCS曲线数和表征土壤有效持水能力参数(SOL_AWC)，给出其在数据同化过程中的演变过程，见图4。从图4可以看出，随着观测信息的不断加入，这两个参数的集合分布越来越窄，集合均值逐渐趋于稳定。从参数集合在某一天的概率密度分布(图4小图)可以看出，参数集合由初始时的均匀分布，通过实时的同化更新，逐渐演变为高斯分布；同时，集合的方差/标准差逐渐变小，说明参数估计的不确定性在逐渐减小。

图4 基于EnKF法径流数据同化的SWAT模型参数CN2和SOL_AWC的演变过程Fig.4 Evolution process of CN2 and SOL_AWC in EnKF based streamflow assimilation

图5为SWAT模型中被优化的七类共59个敏感参数其集合均方差随径流数据同化更新的变化过程。从图5可以看出，七类参数的集合均方差随同化更新过程的继续逐渐减小，表明参数集合的分布越来越窄，即参数估计的不确定性逐渐减小。最终，参数集合的均方差都收敛到初始集合均方差的5%以内，说明基于EnKF法的径流数据同化对SWAT模型参数具有一定的优化更新能力。

图5 SWAT模型中七类参数的集合均方差随径流数据同化的变化过程 Fig.5 Variations of the ensemble mean square errors for 7 types of parameters in SWAT updated with the streamflow assimilation

3.3 模型参数优化效果评估

基于实测径流数据同化进行SWAT模型参数优化时，由于真实(或最优)的模型参数值未知且很多参数无法直接量测，很难直接对参数估计效果进行评估。本文通过基于优化后的参数模拟获得的流域内站点径流过程间接地对参数优化效果进行评估。具体方案为：基于对淮滨流域内5个站点2003—2005年(率定期)实测径流数据同化获得的参数估计结果运行模型，得到2006—2008年(验证期)的径流模拟过程(记为EnKF)；将该径流模拟过程与站点实测径流过程进行比较，通过模拟径流与站点实测径流过程线的接近/贴合程度以及相应的评估指标，对基于EnKF法径流数据同化获得的模型参数优化效果进行评估。同时，将未经率定的参数(即在给定模型输入条件下，SWAT模型默认的参数)在验证期的径流模拟过程(记为SIM)作为参照。

表3为基于EnKF法优化的参数及未经率定的参数获得的流域内5个水文站点在验证期的径流模拟结果。从表3可以看出，基于EnKF法获得的参数估计结果所对应的流域内5个站点在验证期的径流模拟精度相较于未经率定的情况有大幅度的提高。EnKF法对应5个站点径流模拟的相对偏差(除竹竿辅站外)都保持在15%以内，均方根误差显著减小，尤其是流域出口淮滨站的均方根误差减小了34.5%，而纳什效率系数和相关系数r均明显增加，5个站点的纳什效率系数都在0.6以上，其中大坡岭、息县、淮滨3站在0.8以上，尤其流域出口淮滨站达到了0.88，表明基于EnKF法在SWAT模型中同化站点径流数据、更新模型参数与状态变量可以得到较好的模型参数估计结果。

表3 基于EnKF法优化及未经率定的参数条件下流域内5个站点2006—2008年的径流模拟结果

图6为基于EnKF法优化的参数(EnKF)及未经率定参数(SIM)条件下，流域出口淮滨站在验证期(2006-01-01—2008-12-31)的径流模拟过程。从图6可以看出，基于EnKF法参数估计结果的验证期径流(EnKF)，其洪峰过程(图6(a))和基流过程(图6(b)，图6(b)纵坐标为以10为底的对数坐标)都明显优于未经参数率定条件下的径流模拟结果(SIM)。虽然基于EnKF法估计的参数所对应的径流过程相较于实测径流过程还存在一定程度的低估，但其相对于未率定过程而言其低估程度大为减小，表明通过EnKF法同化站点实测径流对SWAT模型参数具有一定的优化能力，可以作为未来SWAT模型参数率定的手段。

图6 基于EnKF法及未经率定的参数条件下，流域出口淮滨站在验证期的径流过程Fig.6 Runoff process at the outlet of Huaibin River Basin based on the EnKF optimized parameters and uncalibrated model parameters

4 结 语

本文基于状态扩展方式的EnKF法，在对SWAT模型误差(模型输入、模型结构、模型参数)和站点径流观测误差采用最大后验估计方法合理量化的基础上，以条件协方差法并引入平滑因子有效刻画参数随时间的演变过程，构建了基于径流数据同化的水文模型参数优化方案。通过站点实测径流数据的同化对SWAT模型参数及状态变量进行实时更新，以淮河上游淮滨水文站以上流域为研究区，对该方案中模型参数的优化效果进行了评估。

结果显示，在站点实测径流数据同化过程中，被优化更新的参数集合逐渐收敛并趋于稳定。通过对比分析基于稳定后的参数获得的模拟径流和实测径流过程，发现通过径流数据同化获得的参数具有较好的模拟预测精度，基于数据同化方法得到的参数进行降雨径流模拟，流域出口淮滨站径流模拟的纳什效率系数可达0.88，间接证明了基于EnKF法的径流数据同化对SWAT模型参数具有一定的优化估计能力，且采用数据同化方式进行水文模型参数的率定具有一定的可行性。

相较于目前常用的参数率定方法(如全局优化方法)将模型模拟误差全部归咎于模型参数误差这一不足[26-27]，通过同化站点径流进行模型参数优化，可综合考虑模型输入、模型结构和模型参数的不确定性，理论上，基于数据同化进行模型参数的优化可获得更优的参数估计。从研究实例来看，当同化进行到600～800时间步长后，模型参数就可收敛且趋于稳定，说明相较于一般的参数率定方法(一般需要一个丰、平、枯水文周期的实测数据)，基于数据同化的参数优化要求更短的实测数据系列，这对于资料匮乏地区(通常数据系列较短)降雨径流模型中参数的优化具有重要意义。然而，基于数据同化的模型参数优化必须建立在对模型模拟误差与径流观测误差合理量化的基础上，而模型模拟误差与径流观测误差的估计一直以来都是水文界的一个难题[8，28]，这也是当前该方法在应用中面临的一大挑战。另外，由于水文模型，尤其具有物理机制的分布式模型的变量与参数众多[29-30]，导致径流数据同化过程中不可避免地存在高维度现象，这样同化效果的不确定性及同化实施的难度都将大大增加，因此，基于径流数据同化进行水文模型参数的优化还需要发展新的、更优的数据同化方法与策略。