在遥感反演获得地表温度、反照率和植被覆盖度等地表特征参数基础上^{[ 22]},选用基于张仁华等^{[ 23]}的像元排序对比法和分层能量切割算法结合蒸散发日尺度扩展模块,开发出模拟日尺度的双层ET遥感模型^{[ 24]},获取1999—2007年流域典型季节代表日的ET。

像元排序对比法主要根据混合像元理论,在8~14 μm波段范围内建立混合地表温度与组分温度间的混合地表温度分解的基本方程。

(1)

其中, T_m, T_s, T_v分别为混合像元、土壤表面和植被冠层温度。 ε_m、 ε_s和 ε_v分别为混合像元、土壤表面和植被冠层的比辐射率, σ为玻尔兹曼常数, f为混合像元的植被覆盖率。本研究中, T_m使用单窗算法^{[ 25]}反演而得, f使用其与植被指数的经验关系获得^{[ 22]}, ε_m则是利用线性混合理论计算^{[ 22]}。

分层能量切割算法的本质是求算土壤波文比 β_s和植被波文比 β_v:

(2)

其中, T_SH和 T_SL表示研究区内等覆盖率像元的最高土壤温度和最低土壤温度, T_VH和 T_VL表示研究区内等覆盖率像元的最高植被温度和最低植被温度,其计算来自于梯形散点图干、湿边的确定^{[ 23]}。在波文比估算基础上,利用波文比能量平衡法切割地表可利用能量获得土壤蒸发和植被蒸腾过程,进而求出蒸发比,再根据蒸发比不变法实现ET的日尺度扩展。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 流域蒸散同化技术路线图Fig.2 The schematic framework of watershed evapotranspiration data assimilation

分布式时变增益模型(Distributed Time Variant Gain Model,DTVGM)是在水量平衡机理基础上,基于水文非线性理论和分布式水文模拟技术建立的流域水循环分析模型^{[ 21]}。该模型应用在黄河、海河、淮河和黑河流域,得到较好的反馈^{[ 26]}。在DTVGM中,ET基于潜在蒸散发和土壤湿度的关系获得,DTVGM将土壤湿度对ET的影响考虑在内,构建了ET与土壤湿度时间上的响应关系,实现了诊断变量和状态变量的传递,具备展开同化的基本条件。 t和 t+1连续时段之间的计算公式可以概括为:

(3)

(4)

式中, ET为实际蒸散发, K_r为土壤水出流系数, AW是土壤湿度, ETp表示潜在蒸散发, WM为饱和土壤湿度, P为实测降雨量, RS为地表径流, KAW为蒸散发折算系数。其中,土壤湿度为模型状态变量,初始时刻的 AW₀基于基于沙河流域自然条件,根据多年统计平均给定;潜在蒸散发ET p实测降雨量 P作为输入变量,实时驱动模型;其余参数通过率定获得。从公式(3)、(4)可以看出 t时刻土壤湿度 AW_t是模型模拟获得ET的状态变量, AW_t的改变将直接影响 t+1时刻的 AW_t+1和ET。

本研究将主要基于EnKF算法构建水文模型顺序同化系统。EnKF针对状态变量加入一组高斯噪声(通常为高斯白噪声)扰动,获得一组符合正态分布的初始状态,并将此作为初始时刻的状态,驱动一组模型模拟,计算下一时刻的状态集合。EnKF利用了集合预报的观点成功的将非线性模型的状态关系转化至线性空间,并且借助集合的形式实现了误差的反馈,具有良好的可操作性^{[ 27]}。主要计算步骤为:

首先,初始化背景场,即对初始时刻的先验状态变量(算例中是土壤湿度) x_t^b加入 m个服从期望为0、方差为 Q的标准正态分布的误差扰动,生成初始背景场,驱动模型并行运行。 t时刻(t≥0)的状态背景场记为 X_t^b。相应的t时刻的状态变量分析场 X_t^a记为:

(5)

式中,分析场 X_t^a是与背景场与观测场对应的后验状态的集合,用来替代水文模型中该时刻的背景场,实现状态变量的更新。 Z_t表示 t时刻的观测场,本研究为遥感蒸散发加一组期望为0、方差为 R的高斯白噪声扰动组成的 m维向量; Z_t^b为观测的预报集合,即水文模型并行模拟的t时刻的蒸散发集合。 K_t为对应时刻的更新算子卡尔曼增益矩阵。

当第 t时刻不存在遥感观测时, K_t为0,状态变量背景场不更新,模型继续向前模拟;当第t时刻存在遥感观测时,进行同化更新,并替换背景场,此时, K_t计算如下:

(6)

式中, R为观测误差协方差, P_t^b为背景场误差协方差, H为观测算子,表示状态变量和观测值之间的计算关系。根据Pipunic等^{[ 20]}的研究, P_t^b和 H并不需要显示计算,公式(6)可化简为:

(7)

其中:

(8)

式中,上划线变量表示由对应变量平均值组成的向量。从上式可以看出,当集合成员数 m为无穷大时, 就是背景场的期望值,而相应的分析场的平均值 就是经过同化计算得到的期望值,在实际操作中选择一定大小的集合空间可以近似的表示这种关系。此外,通过集合同化计算卡尔曼增益时,间接将协方差的更新返回,得到一个包含了误差更新分析场集合。对于EnKF而言,分析场与背景场、分析场与观测场的方差之和最小,实现了概率最优。

采用EnKF数据同化对水文模型校正前,需要确定集合成员数 m、模型误差 Q和观测误差 R^{[ 23]}。理论上,样本空间越大越能代表真实情况,然而过大的集合取值会产生巨大的计算负荷。根据研究经验,随着集合数的增大,同化的性能会逐渐提高,但是提高的速率存在递减趋势,也就是说当集合数较大时,样本空间足可以表达状态整体^{[ 20]}。研究选取1 000作为集合数上限,先固定观测误差和模型误差,基于增长步距法^{[ 28]}以10为步距进行参数优选,最终确定集合数为100。参照遥感模型和水文模型的蒸散发估算误差,设定模型误差方差上限和观测误差上限分别为0.2和0.02,同理采用增长步距的方法确定模型误差方差 Q和观测误差方差 R分别取0.05和0.006。

在划分28个子流域的基础上(图1),通过SUE-CA全局优化算法^{[ 29]}对DTVGM模型进行参数率定,选取1990—1998年作为预热、率定期,1999—2007年作为验证期,计算步长为1天。选择水量平衡系数(Water Balance Coefficient, WB)和纳西效率系数(Nash-Sutcliffe Efficiency Coefficient, NS)2个目标函数评价模拟,具体公式为:

(9)

(10)

式中, Q_s和 Q_o是模拟流量和观测流量,上划线表示是相应流量的均值, n为时间序列长度。在优化过程中,为了综合考虑多个目标函数的影响,再构建一个综合目标函数 OBJ:

(11)

上述指标中, WB和 NS趋向1表示模拟可靠, OBJ趋向0表示模拟可靠。经过SCE-UA参数优化,验证期(1999—2007年)DTVGM模型模拟流量的评价函数 WB、 NS和 OBJ达到1.100、0.778和0.161,表明DTVGM模型对径流的模拟效果较好。观察日ET过程线轨迹(图3)发现,水文模型模拟ET过程时,出现了一些与事实不符的情况:如在2000年4月、2001年5月和2002年5月份日ET平均值分别只有0.65 mm,0.78 mm和0.52 mm;而根据吴炳方等^{[ 30]}和Qin^{[ 19]}等的研究,该流域春季4、5月份流域日均ET均略大于1 mm。此外,水文模型模拟结果与遥感反演结果对照发现(图4),春季水文模拟ET普遍小于遥感反演ET;夏季模拟ET大于遥感反演的ET,冬季小于遥感反演ET,模型模拟可能夸大了ET的波动情况。

针对遥感反演的1999—2007年流域典型季节代表日的ET,利用地面通量站实测数据验证。结果显示:遥感反演与实测相对误差介于0~12%之间,平均绝对百分比误差(Mean Absolute Percentage Error, MAPE)约为10%,相关系数为0.97,有较高的反演精度,具体验证过程参见文献[24]。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 DTVGM模拟流域月蒸散量(1999—2007年)Fig.3 The monthly ET simulated by DTVGM with out data assimilation (1999-2007)

选取沙河流域内的小汤山通量观测站的实测数据来比较双层ET遥感模型、DTVGM水文模型以及同化结果。观测站位于沙河流域1号水文单元内(图1),从2002年开始,利用涡动相关仪和大孔径闪烁仪实时观测了一系列的地表通量数据,并进行了时间尺度转化,生成了日尺度ET数据。该子流域地势平坦,土地利用类型较为单一(农用地占95%以上),土壤类型均为潮土,因此,可以认为该水文单元平均日ET与单元内一定面积内的平均日ET相差不大,可以利用观测站获得的ET作为验证参照。基于此,图4比较了1号水文单元内,对应日的各类计算结果。DTVGM、遥感反演、数据同化的日蒸散发结果与实测数据之间的MAPE分别为25.8%、7.3%和8.2%,表明同化使模拟结果向遥感观测方向偏移。这一结论同时印证了遥感反演精度高于水文模型模拟的结论,也符合模型误差扰动大于观测误差扰动的设定。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同方式获得的区域日蒸散发比较Fig.4 The ET of typical day by different method

图5比较了1999—2007年同化前、同化后的流域平均ET以及对应的遥感反演日ET,为了便于观察,对模拟轨迹取5日移动平均曲线。从图5可以看出,同化前DTVGM模型模拟的ET在夏季大于数据同化后的结果,在冬季小于同化结果,整体的波动范围要大于同化结果。经过数据同化,DTVGM水文模型在2000年4月、2001年和2002年的5月出现ET明显偏小的现象得到了修正。上述3个月的日均ET由小于1 mm分别升至1.39 mm,1.55 mm和2.01 mm,与同类研究结果近似^{[ 19, 30]}。同时,从图5可以直观的看出,遥感观测输入的当日蒸散发改变最为明显;从2001年开始到2003年上半年,遥感观测较多,相应时段内数据同化前后的ET轨迹差别较大;在遥感数据缺乏的2004年、2005年和2006年3年中,ET的同化轨迹与同化前模拟轨迹逐渐接近,2006年和2007年二者几乎重合,在2007年又出现了一次遥感观测,此后蒸散发过程的同化轨迹开始偏移同化前模拟的轨迹。上述现象反映了遥感ET的输入是决定同化的关键,当存在遥感观测时同化效果好;当没有遥感ET输入时,随着模拟的进行,同化轨迹逐渐偏向同化前的模拟轨迹,直到下一次遥感观测出现,同化轨迹又出现一定幅度的改变;观测的输入频率决定着数据同化对模拟的改变程度。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 1999—2007年沙河流域逐日流域平均蒸散发(5日移动均线)Fig.5 Average daily ET processes in the Shahe river basin in 1999-2007 (5 days moving average curve)

选取沙河流域内通量观测站的ET实测数据来比较DTVGM水文模型同化前、后的输出结果。由于地面通量站观测并不连续,因此分别选取全年观测效果较好的2002年和2007年的实测数据进行验证,比较对应日的DTVGM模拟、数据同化与实测ET,如图6所示。

图6

Fig. 6

	Figure Option ViewDownloadNew Window
	图6 日蒸散发的DTVGM模拟和同化结果比较(2002年、2007年)Fig. 6 The comparation of ET in typical days through DTVGM and data assimilation (2002, 2007)

由图6结果得出,同化前DTVGM模型模拟与地面实测相差较大,经过数据同化,2002年ET模拟得到较为显著的优化,与实测值更加接近,且相关系数较高( R²=0.8352)。相比较而言,2007年由于仅有一天的观测输入,较少的观测输入导致数据同化对模拟结果的校正不明显。

为进一步论证同化效果,选择验证效果较好的ET-Watch模型^{[ 30]}的年蒸散发结果进行比较(表2)。同化结果与ET-Watch的年尺度流域ET差别在10%以内,较为合理,而DTVGM直接模拟结果与ET-Watch的差别较大。同化系统除了借助EnKF算法外,还有一个关键环节就是误差反馈,是通过更新状态变量实现水循环过程模拟的整体更新,由于流域产流等水文变量均是从状态变量转化而来,因此沙河流域地表径流等变量也会发生一定程度的变化。经验证,水文模型经同化后的 WB, NS和 OBJ分别为1.025,0.758和0.134,与同化前结果(1.100、0.778和0.161)相比,从水量平衡的角度来看模拟精度有所提升,流量过程的模拟精度略有降低,综合评价指标表明模拟效果提高,整体水循环模拟是改善的。

表2

Table 2

表2(Table 2)

表2 DTVGM模拟、数据同化输出与ET-Watch模型估算的年蒸散发比较 Table 2 The Comparation of yearly ET simulated by DTVGM, DA and ET-Watch 年份 同化前年ET(mm) 模拟相对误差(%) 同化后年ET(mm) 同化相对误差(%) ET-Watch结果(mm) 2002 452.18 12.74 564.78 1.90 575.50 2003 454.78 -4.52 450.56 -3.42 435.12 2004 543.82 -12.09 506.31 -4.18 485.15 2005 501.09 10.85 520.25 8.04 562.05

表2 DTVGM模拟、数据同化输出与ET-Watch模型估算的年蒸散发比较 Table 2 The Comparation of yearly ET simulated by DTVGM, DA and ET-Watch