基于EnKS和SWAT模型的闽江流域径流数据同化

2023-09-11陈芸芝唐丽芳汪小钦

水资源与水工程学报 2023年4期

项勇, 陈芸芝, 唐丽芳, 汪小钦

(1.福州大学数字中国研究院(福建), 福建福州 350108; 2.福州大学空间数据挖掘和信息共享教育部重点实验室,福建福州 350108; 3.福建省水土保持实验站, 福建福州 350001)

1 研究背景

流域水文模型可以反映整个流域的综合水文信息,在水旱灾害防治、水资源管理和开发利用等方面应用广泛[1-3]。但受模型结构、观测数据、气象输入等不确定性因素的影响,水文模型在进行水文状态模拟时存在不可避免的偏差[4]。数据同化技术能够综合不同观测信息,以获得流域状态变量的最佳估计,从而减少水文模型径流模拟偏差[5]。已有研究表明,同化遥感数据,如土壤湿度、蒸散发、积雪等可以有效减少模型输出的不确定性[6-8]。然而由于径流观测数据与径流模拟数据直接相关,且地表实测径流数据具有较高可靠性,同化站点径流实测信息仍是径流预测的最优选择[9]。

数据同化方法可分为变分数据同化和顺序数据同化。集合卡尔曼滤波(ensemble Kalman filter, EnKF)是一种顺序同化方法,可以根据实测信息顺序更新模型状态。相较于普通卡尔曼滤波,EnKF解决了实际应用中背景误差协方差矩阵计算困难的问题,并且能够直接应用于非线性模型[10]。同时,相较于更复杂的粒子滤波等方法,EnKF的计算效率更高,在水文数据同化中得到了广泛的应用[11-13]。标准EnKF方法在小流域尺度以及状态变量与观测变量之间实时响应的情况下取得了良好的效果[14-15],但该方法未考虑土壤水等水文状态变量和河道出口流量之间的时间滞后问题,导致中大型流域的径流数据同化效果较差[16]。针对时间滞后问题,学者们提出了不同解决方法,如Weerts等[17]设置了固定滞后时间,通过当前观测值更新固定滞后时间前的模型状态变量;McMillan等[18]基于REnKF方法,利用当前观测值递归地更新最大滞后期前至当前观测时间的模型状态变量,结果表明该方法可以解决集水区土壤含水量和观测流量之间的响应滞后问题,但其使用同一观测数据迭代更新先前模型状态,存在过度校正的风险,且计算复杂。集合卡尔曼平滑器(ensemble Kalman smoother, EnKS)是EnKF的重要变体,相较于EnKF方法,EnKS综合考虑了时间窗口内的状态变量和观测值,使得预测结果更为平滑、准确,能够有效解决径流数据同化中的时间滞后问题。该方法在集总式水文模型中表现良好[10,19],但目前仍较少应用于分布式和半分布式水文模型。

本文以闽江流域为研究区,基于EnKS方法和半分布式水文模型SWAT(soil and water assessment tool),构建径流数据同化方案,在确定最佳时间窗口长度的前提下,评估EnKS方法对半分布式水文模型的优化效果,分析数据同化对于不同径流分量的影响,为闽江流域径流数据同化提供参考。

2 数据来源与研究方法

2.1 研究区概况

本文选取闽江流域为研究区,如图1所示。研究区流域面积为59 922 km2,地处东经116°23′～119°35′,北纬25°23′～28°16′,主要支流包括建溪、富屯溪和沙溪等。闽江流域属于亚热带季风气候,雨量充沛,但分布不均,降水主要集中在夏秋季节,春冬季节降水较少,多年平均降水量为1 710 mm。

图1 闽江流域水系及水文、气象、雨量站点分布

2.2 数据来源

闽江流域SWAT模型构建所使用的数据如表1所示。气象数据为中国气象数据网(http://data.cma.cn/)提供的泰宁、福州、浦城、建瓯等8个气象站点2017—2020年逐日气温、相对湿度、风速等数据,以及福建省水文水资源勘测局提供的宁化、大田、松溪、武夷山等10个雨量站点2017—2020年逐日降水数据。实测径流数据为福建省水文水资源勘测局提供的七里街水文站、沙县水文站、竹岐水文站2018—2020年逐日实测径流量数据。

表1 闽江流域SWAT模型构建中使用的数据

2.3 SWAT模型

SWAT模型是应用最广泛的半分布式水文模型之一[20],该模型根据地形将整个流域划分为多个子流域,并根据下垫面条件进一步划分水文响应单元(hydrological response unit,HRU),HRU是SWAT最基本的计算单元。

SWAT模拟的逐日水文循环基于每个HRU的水量平衡方程进行,其计算方程如下:

(1)

式中:SWt为第t天结束时的土壤含水量,mm;SW0为初始土壤含水量,mm;Ri为第i天的总降水量,mm;Qsurf,i为第i天的地表径流量,mm;Ei为第i天的蒸散量,mm;Wseep,i为第i天的下渗量和壤中流,mm;Qgw,i为第i天的地下水径流量,mm。

总径流Qi由第i天的地表径流Qsurf,i、壤中流Qlat,i和地下水径流Qgw,i组成:

Qi=Qsurf,i+Qlat,i+Qgw,i

(2)

各HRU产生的总径流在子流域尺度上汇集,继而进入河网汇总,并在河道出口断面处计算径流。

2.4 集合卡尔曼滤波

集合卡尔曼滤波基于Monte Carlo方法和标准卡尔曼滤波公式来估计非线性模型中的状态和参数[10]。EnKF主要分为预报和更新两步。

在预报步骤中,将初始化后的k时刻状态变量集合代入模型,得到k+1时刻的状态变量预测值:

(3)

在基于模型状态预报值获得模型状态预报协方差矩阵后,模型状态变量的预报值基于观测信息更新,其公式如下:

(4)

(5)

(6)

(vk+1～N(0,Rk+1))

(7)

2.5 集合卡尔曼平滑器

EnKS是EnKF的扩展,该方法将平滑窗口内存在的全部观测与模拟分别纳入观测变量值集合与模型状态变量集合,然后通过EnKF公式更新状态变量。EnKF方法即为EnKS时间窗口等于1天的特殊情况。

在初始化随机变量集合后,SWAT模拟运行至平滑窗口最终时刻,获取窗口内所有状态变量模拟值,因此扩展的状态变量X包括平滑窗口内所有L个状态值x,其可表示为:

X=[x1x2…xL]T

(8)

相应地,扩展的观测向量集合Y包括时空独立分布的所有m个观测值,即为:

Y=[y1y2…ym]T

(9)

而后使用EnKF公式计算卡尔曼增益矩阵,对模型状态变量进行更新。

2.6 评价指标

为评估模型预测与实测值的接近程度,采用纳什效率系数(Nash-Sutcliffe efficiency coefficient,NSE)、均方根误差(root mean square error,RMSE)共同评估模型模拟效果[21]。其中,NSE反映了模拟值和实测值的拟合度,NSE越趋近于1则模拟效果越好;RMSE反映了模拟值与实测值的偏差,RMSE越趋近于0,代表模拟效果越好。

3 径流数据同化实验

图2 径流观测数据同化方案流程

3.1 SWAT模型参数率定

选取闽江流域七里街、沙县、竹岐3个水文站实测径流量数据对SWAT模型进行率定和验证,以2018年1月1日—2019年12月31日为率定期,2020年1月1日—2020年12月31日为验证期,流域参数率定结果如表2所示。

表2 流域参数率定结果

结果表明,七里街、沙县、竹岐3个水文站的NSE在率定期分别为0.76、0.74、0.82,验证期分别为0.70、0.50、0.72。一般认为当NSE>0.5时,模型模拟结果可被接受[22]。因此该SWAT模型可用于后续闽江流域径流数据同化研究。

3.2 状态变量选择与更新

土壤水决定地表水文循环的诸多过程,是水文模型模拟中最重要的变量之一,在SWAT同化中可以将所有HRU的土壤含水量纳入状态向量中进行计算,但这增加了同化的计算量,因此本文定义了土壤水校正比率分别控制各个站点上游的土壤含水量,共设置3个比率,其控制范围如图3所示。HRU尺度上各层土壤含水量更新公式如下:

图3 研究区不同土壤水校正比率空间分布

sol_sta(c,b)=sol_st(c,b)·ratio

(10)

式中:sol_sta(c,b)为校正后第b个HRU中第c层土壤的含水量,sol_st(c,b)为未校正前第b个HRU中第c层土壤的含水量,ratio为土壤水校正比率。

为避免HRU尺度上土壤含水量的过度更新,将土壤水校正比率设置在0.8～1.2之间,将土壤含水量设置在1×10-6mm至饱和土壤含水量之间。

3.3 误差概化

降水是SWAT输入数据中最重要的数据之一。本文假设降水输入误差服从对数正态分布[23],且降雨输入误差在时间和空间上独立,即:

P′=Pεp(εp～lnN(μ,σ2))

(11)

式中:P′为降雨经误差扰动后的降雨集合;P为站点实测降水量;εp为降雨误差扰动项,服从对数正态分布,均值为μ,方差为σ2。参考以往研究设定μ=1,σ2=25%[24]。

采用比例因子法确定径流观测误差协方差矩阵[25]:

Rk=(αYk)2

(12)

式中:Rk为k时刻的观测误差协方差矩阵;Yk为k时刻的观测向量;α为比例因子,参考以往研究,设定α=0.1[25]。

4 结果与分析

4.1 滞后时间窗口确定

基于EnKS方法的同化模型性能受滞后时间窗口大小影响。随着时间窗口的增大,需要考虑更多与当前状态变量无关的误差信息,若低估窗口长度,则损失部分滞后径流信息。同时最佳时间窗口长度在枯水期和丰水期会发生变化。以七里街站为例,选取2018年1月1日—2018年4月1日为枯水期,2019年5月1日—2019年8月1日为丰水期,计算不同时期模型径流量随滞后时间窗口长度变化的模拟精度,结果见图4。

图4 不同时期和时间窗口长度下七里街站同化结果

由图4可以看出,七里街站点枯水期最佳时间窗口长度为5 d,丰水期最佳时间窗口长度为1 d。丰水期最佳时间窗口长度明显短于枯水期,这是因为丰水期河道径流量大,SWAT模型中河道径流流速随水位升高而快速增加,径流的滞后性降低,最佳时间窗口长度缩短。因此针对枯水期和丰水期的径流数据同化,应分别选取较长和较短的时间窗口。

2018年1月1日—2020年12月31日包含枯水期和丰水期,为确定该时期最佳时间窗口,本文对比了该时期不同站点在不同时间窗口长度下的逐日径流模拟精度,如图5所示。

图5 不同时间窗口长度下各水文站点同化结果

从图5中发现模型总体精度随时间窗口的增加呈现先增加后下降的趋势。当同化时间窗口较小时,随着时间窗口长度的增加,模型能够获取更多滞后径流信息,从而更好地改善结果。然而当同化时间窗口过长时,后续径流模拟误差并不与当前水文单元土壤水误差直接相关,主要由其他误差导致,而模型仍将后续误差直接归因于当前状态变量误差,使得同化模型模拟精度下降。综合丰、平、枯水文周期,最终确定七里街站点最佳时间窗口长度为2 d,沙县站点最佳时间窗口长度为5 d,竹岐站点最佳时间窗口长度为4 d。

4.2 径流量模拟

基于率定后的SWAT模型,利用EnKS方法和EnKF方法分别模拟闽江流域2018—2020年逐日径流量,并与原模型模拟结果进行对比,如图6所示,图6(d)、6(e)、6(f)中的纵坐标为以10为底数的对数坐标。由图6可见,原模型在枯水期径流量模拟中存在明显的低估现象,在丰水期则存在不同程度的低估和高估现象。而EnKS方法和EnKF方法模拟结果更为接近观测值,表明同化站点实测径流量能够有效校正模型状态误差。

图6 2018-2020年各水文站不同同化方法及原模型径流量模拟结果

相较于EnKS方法,EnKF方法存在明显的时间滞后现象。以2018年初为例,如图7所示,EnKF方法的径流峰值出现在观测峰值之后,形成虚假峰值。这是由于SWAT模型中土壤水和径流之间的映射存在滞后,在枯水期土壤水校正对当天径流的影响极小,导致EnKF方法低估了土壤水变化的影响,对土壤水的校正不足。此外,当土壤水校正对当天径流产生影响时,EnKF方法仅能获取当天径流变化信息,忽视了其余滞后径流信息,导致EnKF方法对土壤水过度校正,出现虚假峰值。EnKS方法综合考虑了滞后时间窗口内的模型模拟径流和观测径流,能在观测时间前对土壤含水量进行校正,模拟结果更为准确。

图7 2018年1—4月各水文站不同同化方法及原模型径流模拟结果

表3为不同同化方法及原模型对2018—2020年各水文站点径流模拟精度的评价。由表3中的评价参数可知,相较原模型模拟结果,EnKS方法和EnKF方法均显著改善了径流模拟结果。其中EnKF方法模拟的七里街、沙县、竹岐3个站点径流NSE分别提升了0.06、0.03、0.01,RMSE分别降低了10.29%、5.43%、8.15%;EnKS方法模拟的七里街、沙县、竹岐3个站点径流NSE分别提升了0.09、0.15、0.04,RMSE分别降低了16.95%、30.78%、12.05%。3个站点NSE均达到0.80以上,其中沙县站点效果最好,NSE达到0.86,表明EnKS方法在SWAT模型中同化站点径流数据能够显著改善径流估计精度。EnKF方法虽然同样改善了径流估计值,但改进程度不及EnKS方法。EnKS方法模拟结果的NSE在七里街、沙县、竹岐3个站点上比EnKF方法分别提升了0.03、0.12、0.03,RMSE分别降低了7.43%、26.81%、4.25%。

表3 2018—2020年各水文站不同方法径流模拟精度评价

4.3 校正土壤水对不同水文变量的影响

土壤水的校正影响地表径流、壤中流等其他水文组分,分析土壤水校正对不同水文组分的影响,能够为水文数据同化的状态变量选择提供一定参考。选择典型HRU分析土壤水校正对其他水文变量的影响,其下垫面性质如表4所示。

表4 所选典型HRU下垫面性质

图8为原SWAT模型和EnKS校正土壤水后模型变量基流(浅层地下水补给)、深层地下水补给、壤中流、地表径流在2018年1月1日—2020年12月31日的变化情况。

由图8可见,闽江流域SWAT模型中水文响应单元产流主要由基流、壤中流和地表径流组成,深层地下水补给在产流中占比极少,且受土壤水变化影响不大,因此深层地下水补给并不是水文数据同化的理想状态变量。而在坡度平缓地区,土壤水以下渗为主,较少产生侧向径流,壤中流在产流中的占比较少,因此对于坡度平缓地区,壤中流并不是理想的同化状态变量。

以2018年1月1日—2018年4月1日为枯水期,2019年5月1日—2019年8月1日为丰水期,分别计算同化前后HRU不同径流分量值在枯水期和丰水期之和,结果如表5所示。

表5 所选典型HRU在枯水期和丰水期同化前后不同径流分量 mm

表5中的计算结果表明,SWAT模型中土壤水的校正主要改变了模型中的基流、壤中流和地表径流。在枯水期,土壤水的校正主要改变了基流和壤中流;在丰水期,土壤水的校正主要改变了基流、壤中流和地表径流。

对比HRU-1167、HRU-1168、HRU-1169 3种不同坡度典型下垫面,发现壤中流随坡度的增加而增加,SWAT模型中壤中流的计算与坡度成正比,因此当流域地势较为平缓时,土壤水校正对壤中流的改进不明显。分别对比高坡度下垫面HRU-1166、HRU-1169和低坡度下垫面HRU-1205、HRU-1279,发现土壤性质差异同样影响土壤水校正对壤中流的改进程度。由于SWAT模型中壤中流计算与饱和水力传导系数SOL_K成正比,而不同类型土壤的饱和水力传导系数存在差异,导致土壤水校正在高渗透率土壤区域对壤中流的改进更显著。

5 讨论

EnKS方法和EnKF方法均可纠正SWAT模型土壤水误差,改善水文变量估计结果,提高SWAT模型对流量的预测准确度。实验结果显示基于EnKS方法的同化效果优于EnKF方法,这和Li等[9]在GR4H水文模型以及Meng等[19]在新安江模型中所进行的实验结果相一致,表明考虑水文模型的时间滞后性可以有效提升模型同化精度。

时间窗口的长度对EnKS方法的性能影响较大。为减少计算量,本文以比率的方式对土壤水进行校正,这需要综合不同HRU中土壤水至河道出口的时间,以确定最佳时间窗口长度。由于SWAT模型中径流流速与径流量相关,因此最佳时间窗口长度随径流量的变化而变化。此外,由于不同流域下垫面性质的差异,最佳时间窗口在不同流域也会发生变化,这为快速准确地确定最佳时间窗口长度带来了一定困难。本文综合了2018年1月1日—2020年12月31日径流数据同化总体精度,确定了闽江流域不同区域的最佳同化时间窗口长度。后续可以考虑采取自适应方法,自动确定不同时期不同区域最佳时间窗口长度,以更好地校正模型土壤含水量,改善径流估计结果。

由于水文模型机制和下垫面性质的空间异质性,土壤水校正对不同径流分量的改进程度在空间上存在差异。雷芳妮[26]的研究同样表明,土壤水校正在高渗透率土壤类型上对于壤中流等水文变量能取得更显著的改进。除土壤类型外,坡度差异也会导致土壤水校正对壤中流的改进程度不同。此外,数据同化算法对不同径流分量的改进程度存在时间异质性,在枯水期土壤水校正主要改变了基流和壤中流,而在丰水期主要改变了基流、壤中流和地表径流。因此,选择地表径流和壤中流等相关状态变量作为同化对象时,应考虑水文周期和区域下垫面性质,以避免与模型机制相冲突。