邢貞相 ，傅 爽，孫明新，紀 毅 ，付 強 *，李 衡

（1.東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，哈爾濱 150030；2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)水資源高效利用重點實驗室，哈爾濱 150030）

蒸散發(fā)（ET）是水文循環(huán)中重要環(huán)節(jié)，影響流域產(chǎn)流量計算結(jié)果。在濕潤地區(qū)，年蒸散發(fā)量約占年雨量總量50%；在干旱地區(qū)，則高達90%[1]。因此，蒸散發(fā)估算為影響水文模型模擬精度重要因素?，F(xiàn)有蒸發(fā)觀測設(shè)備多為E601 型蒸發(fā)器和20 cm 口徑蒸發(fā)器，通過觀測水面蒸發(fā)量估算陸面蒸發(fā)為常用方法。但蒸發(fā)觀測站布置相對雨量站布置更稀疏，較大流域分布極不均勻，難以監(jiān)測流域蒸發(fā)情況。在相對較小流域無蒸發(fā)觀測站，需間接計算方法推求，如蒸發(fā)量等值線或鄰近測站估算數(shù)據(jù)。因此，傳統(tǒng)地面觀測手段難以獲得準確大尺度實際蒸散發(fā)及其空間分布特征，利用遙感反演可獲得大尺度區(qū)域蒸散發(fā)數(shù)據(jù)[2]。遙感蒸散數(shù)據(jù)經(jīng)同化后，綜合遙測數(shù)據(jù)和地面數(shù)據(jù)信息，其精度更高。因此，遙感蒸散發(fā)數(shù)據(jù)同化問題已成為遙感數(shù)據(jù)應(yīng)用研究熱點，當前與遙感蒸散發(fā)相關(guān)同化研究已取得一定成果。Kalma 等綜述利用遙感估算蒸散發(fā)量方法，并指出該方法局限性[3]。Schuurmans等利用陸面能量平衡模型（SEBAL）反演ET，改善水文模型對高海拔地區(qū)ET 估算[4]。Jang等利用MODIS（MODerate resolution Imaging spectroradiometer）和MM5（5th generation meso-scale meteorological model）同化蒸散量，發(fā)現(xiàn)可估算ET[5]。Dong等對比分析ET同化常用水文模型（即SWAT、VIC、MIKE SHE、VIP、TOPMODEL、IHDM、Xinanjiang 和DTVGM），結(jié)果表明，分布式時變增益模型（DTVGM）更適合構(gòu)造同化方案[6]。Xu等基于集合卡爾曼濾波，首次將新一代中國氣象衛(wèi)星（FY3AVIRR）數(shù)據(jù)同化，結(jié)果表明陸地模型（CoLM）低估蒸發(fā)率[7]。Xie等研究表明基于集合卡爾曼濾波同化后蒸散發(fā)數(shù)據(jù)可有效用于SWAT模型參數(shù)估算[8]。同化后蒸散發(fā)數(shù)據(jù)應(yīng)用于徑流模擬研究取得成果。Qin 等利用擴展卡爾曼濾波對SEBS 反演ET 和分布式水文模型WEP-L模擬ET作數(shù)據(jù)同化，但更新ET未對徑流模擬產(chǎn)生影響[9]。Yin 等利用集合卡爾曼濾波算法（EnKF）基于雙層遙感模型與分布式時變增益水文模型（DTVGM）作同化，同化后流量過程模擬精度略降低[10]。Zou 等采用DTVGM 模型，分別以MOD16 和FLUXNETMET 蒸散發(fā)產(chǎn)品為ET 觀測值，通過EnKF 同化提高徑流模擬精度[11]。

目前，蒸散發(fā)同化數(shù)據(jù)在徑流模擬中應(yīng)用較少，次洪模擬應(yīng)用更少，且模擬效果不穩(wěn)定。本文利用集合卡爾曼濾波算法將蒸散發(fā)數(shù)據(jù)同化后，應(yīng)用于洪水模擬，為無地面蒸發(fā)資料地區(qū)洪水模擬提供參考。為更好地在水文模型中考慮蒸散發(fā)空間分布不均勻性，本研究采用半分布式水文模型，利用高分辨率DEM（Digital Elevation Model）結(jié)合遙感反演區(qū)域蒸散發(fā)數(shù)據(jù)模擬洪水。

1 研究區(qū)及數(shù)據(jù)資料

七星河流域位于中國黑龍江省三江平原東北部，東經(jīng)131°9′～131°44′、北緯46°12′～46°38′?？偯娣e為1 344 km2，其中山地占總面積38.3%，主要分布于流域西南部和南部；平原占總面積61.7%，主要分布于流域北部和中部。七星河流域?qū)俸疁貛Т箨懶詽駶櫦撅L(fēng)氣候區(qū)，多年平均降雨量534.5 mm，年內(nèi)分配不均，6至9月份降雨量占全年72%。流域多年平均水面蒸發(fā)量486～703 mm。夏季最高氣溫38.3 ℃，冬季最低氣溫-37.2 ℃，年平均氣溫2.3～3.4 ℃。

七星河流域內(nèi)共設(shè)4個雨量站、1個水文站（同時監(jiān)測降雨、徑流），平均每268.8 km2布設(shè)一個站點。以上各站可提供2006～2013 年10 場實測場次降水數(shù)據(jù)。采用反距離權(quán)重法（IDW）[12]對雨量數(shù)據(jù)作空間插值，得到七星河流域30′′×30′分辨率小時降水插值場。流量數(shù)據(jù)采用該流域出口斷面保安水文站2006～2013年間實測徑流數(shù)據(jù)。由于流域內(nèi)無蒸發(fā)觀測站，本文移用與七星河流域具有相近緯度和氣象條件寶清站（132°15′E 46°20′N）實測蒸發(fā)資料作為七星河流域蒸發(fā)數(shù)據(jù)。遙感蒸散發(fā)數(shù)據(jù)選取SSEBop 月尺度蒸散發(fā)產(chǎn)品[13]。該蒸散發(fā)產(chǎn)品通過operational Simplified Surface Energy Balance（SSEBop）模型產(chǎn)生[14]。SSEBop 使用基于簡化表面能量平衡（SSEB）方法，結(jié)合遙感MODIS 熱圖像生成Evapotranspiration（ET）（每10 d 采集1 次）和參考ET 熱 指 數(shù)[15-16]。最 初SSEB 公 式 基 于SEBAL 和METRIC 模型冷與熱像素原理[17-18]。該蒸散發(fā)數(shù)據(jù)空間分辨率為1 km。此外，本研究中還采用美國USGS提供30〃×30〃分辨率DEM（見圖1）。

2 研究方法

本文采用集合卡爾曼濾波算法同化SSEBop 蒸散發(fā)產(chǎn)品和地面蒸發(fā)器獲得估測值。分別將地面估測值和同化后蒸散發(fā)值應(yīng)用到基于蓄滿產(chǎn)流半分布水文模型中，以檢驗同化后遙感蒸散發(fā)數(shù)據(jù)有效性。以實測流量過程為參照標準，選擇納什效率系數(shù)（ENS）、洪峰誤差、洪峰滯時為評價指標，評價水文模型徑流模擬精度。

2.1 集合卡爾曼濾波方法

卡爾曼濾波理論由Kalman 于1960 年首先提出，是一種利用線性系統(tǒng)狀態(tài)方程，根據(jù)系統(tǒng)輸入輸出觀測數(shù)據(jù)，對系統(tǒng)狀態(tài)作最優(yōu)估計算法，應(yīng)用廣泛[19]。為克服傳統(tǒng)卡爾曼濾波僅限于處理線性問題不足，使其可應(yīng)用于非線性問題，Evensen首次提出集合卡爾曼濾波[20]。集合卡爾曼濾波方法采用蒙特卡洛隨機采樣，依靠隨機生成一系列樣本，通過計算樣本隨機誤差，得到模型誤差分布，用于處理非線性序列同化問題。

本研究中，對當前時刻蒸發(fā)皿實測月蒸發(fā)量（ETm）加入均值為0，標準差為10 mm符合高斯分布隨機擾動，樣本數(shù)取為30，作為背景場。選擇該組標準差和樣本數(shù)量獲得較好同化效果。將遙感蒸散發(fā)產(chǎn)品作為觀測值更新背景場，得蒸散發(fā)分析值及其相應(yīng)誤差協(xié)方差矩陣，將分析值誤差協(xié)方差作為下一時刻背景場誤差，如此循環(huán)，可得到一組新蒸散發(fā)序列。詳細步驟可參見文獻[11]。

2.2 基于蓄滿產(chǎn)流半分布式水文模型

七星河流域?qū)儆诎霛駶櫟貐^(qū)，夏季降雨較多，基本符合蓄滿產(chǎn)流應(yīng)用條件。模型以DEM 柵格為計算單元，并假設(shè)柵格單元內(nèi)氣象、地形、土壤、植被等條件空間分布均勻，僅考慮以上各要素在不同柵格間差異性。模型將土壤分為上層、下層、深層3 層，采用3 層蒸散發(fā)模型[1]，將徑流劃分為地表徑流、壤中流和地下徑流。由于降雨初期流域下滲能力通常大于降雨強度[21]，降水先補充深層地下水，后補充下層和上層土壤水，待蓄滿后產(chǎn)生地表徑流。地下水和壤中流均采用線性出流模式，以此求出每個柵格出流量。地下徑流和壤中流出流之后，先補給下一個柵格土壤水，蓄滿后產(chǎn)生地表徑流。地表徑流則根據(jù)DEM提取出柵格間流向，求得每個柵格到流域出口距離，水流在坡面上運動速度用公式（1）估算：

式中，vi為網(wǎng)格單元i 處流速，si為網(wǎng)格單元i處坡面坡度，可由DEM求得，a為經(jīng)驗系數(shù)。

水流在河道中運動速度由Manning 公式（2）計算：

式中，vi為網(wǎng)格單元i 處流速（m·s-1）；ni為網(wǎng)格單元i 粗糙度；Sli為網(wǎng)格單元格i 處縱向斜率（m·m-1）；Ri是網(wǎng)格單元i處水力半徑。

將每個柵格出流過程累加后，可得流域出口徑流過程。具體細節(jié)參見文獻[22]。

2.3 參數(shù)率定方法

遺傳算法模擬生物自然選擇和群體遺傳機制，有效解決復(fù)雜非線性組合問題及多目標函數(shù)優(yōu)化問題。遺傳算法具體步驟參見文獻[23]。

3 結(jié)果與分析

3.1 水文模型率定

本文所構(gòu)建分布式水文模型共有9 個參數(shù)（見表1）。在模型參數(shù)率定過程中，先根據(jù)模型參數(shù)物理意義初步確定每個參數(shù)取值范圍（見表1），然后采用遺傳算法同時優(yōu)化所有參數(shù)。參數(shù)率定期選取七星河流域2006～2011 年間6 場次洪水資料，洪水編號分別為No.060614、No.070805、No.090714、No.090818、No.100804、No.110603。率定期蒸發(fā)資料移用寶清站E601 型陸面蒸發(fā)器實測數(shù)據(jù)，時間尺度為1 h。參數(shù)率定結(jié)果見表1，率定期洪水模擬結(jié)果見圖2，率定期各場洪水模擬精度見表2。

由表2可知，模型率定期各場洪水模擬納什效率系數(shù)均在0.7 以上，平均納什效率系數(shù)達0.84；單場洪水總徑流量相對誤差均在20%以內(nèi)，平均總徑流量誤差為10.1%；單場洪峰流量相對誤差在10%以內(nèi)，平均洪峰流量相對誤差為4.55%；單場洪水峰現(xiàn)時差最大為4 h、最小為1 h。率定期090714 場次和100804 場次效率系數(shù)較低，原因在于這兩場洪水均有復(fù)峰，說明本文構(gòu)建水文模型復(fù)峰洪水模擬精度略低于單峰洪水模擬精度，但總體精度符合場次洪水模擬要求。

表1 七星河流域參數(shù)率定結(jié)果Table 1 Calibration results of parameters of Qixing river basin

3.2 土壤初始自由水含水量確定

土壤初始自由水含水量反映土壤初始水量情況，其對水文模型徑流模擬結(jié)果有較大影響，在模型中需先估算土壤初始自由水含水量。

由于土壤初始自由水含水量難以適時測量，本文根據(jù)地下水和河道水互相補給機制，依據(jù)河道初始徑流量估算土壤初始自由水含水量。即事先建立七星河流域不同場次降雨開始時土壤初始自由水含水量和相應(yīng)河道初始徑流量散點圖（見圖3），并采用非線性擬合技術(shù)獲取兩者間相關(guān)函數(shù)式（3）。利用式（3）計算各場降雨開始時土壤初始自由水含水量用于水文模型徑流計算。

土壤初始含水量與河道徑流量擬合曲線方程為式（3）：

式中，F(xiàn)M為土壤初始自由水含水量（mm），IR為河道初始徑流量（m3·s-1）。式（3）反映七星河河道初始徑流量和流域土壤初始自由水含水量之間數(shù)量關(guān)系。

3.3 蒸散發(fā)數(shù)據(jù)同化

由于七星河流域年內(nèi)降雨主要集中在夏季。因此，選取2006～2013年間6～9月汛期蒸散發(fā)數(shù)據(jù)同化處理。先將1 km 分辨率遙感蒸散發(fā)數(shù)據(jù)重采樣為30′′分辨率，將臨近蒸發(fā)站寶清站E601 蒸發(fā)器測量值乘以蒸散發(fā)折算系數(shù)1（該系數(shù)基于以往經(jīng)驗值得出），得出實際地面蒸散發(fā)。在月尺度和每個30′′×30′柵格上，將臨近蒸發(fā)站寶清站測量值作為模型算子，以SSEBop 遙感蒸散發(fā)模型反演蒸散發(fā)作為觀測數(shù)據(jù)，輸入集合卡爾曼濾波算法中同化計算，流域出口寶安站所在柵格計算結(jié)果如圖4所示。由圖4可知，同化后蒸散發(fā)量更接近于地面蒸發(fā)站實測數(shù)據(jù)，說明EnKF可有效同化遙感蒸散發(fā)數(shù)據(jù)。

同化后蒸散發(fā)數(shù)據(jù)時間尺度單位（月），洪水模擬時段長度單位（h），為使同化后蒸散發(fā)數(shù)據(jù)時間尺度適應(yīng)洪水模擬單位時段長度，驗證期內(nèi)各月同化后蒸散發(fā)量根據(jù)式（4）轉(zhuǎn)化為每小時蒸散發(fā)量：

式中，ETh為小時蒸散發(fā)量（mm），ETEnKF為同化后月蒸散發(fā)量，ETm為蒸發(fā)皿實測月蒸發(fā)量，ETd為蒸發(fā)皿實測日蒸發(fā)量。

3.4 模型驗證及應(yīng)用

為驗證率定后水文模型和同化后蒸散發(fā)數(shù)據(jù)有效性，現(xiàn)選取2012～2013 年間4 場洪水資料（No.120728，No.120912，No.130629，No.130722）驗證，驗證期各場次洪水模擬結(jié)果見圖5。圖5 中同時給出利用蒸發(fā)皿實測數(shù)據(jù)計算洪水模擬結(jié)果（Flow-calculated）和利用EnKF 同化后蒸散發(fā)數(shù)據(jù)計算洪水模擬結(jié)果（EnKF-Calculated）。驗證期各場洪水模擬精度評價見表2。

由圖5可知，驗證期各場洪水模擬過程線與實測洪水過程線擬合精度均較高。圖5a、5b 中肉眼僅見2 條過程線，實際圖5 中有3 條過程線。因為No.120728和No.120912兩場洪水分別利用蒸發(fā)站實測蒸散發(fā)數(shù)據(jù)和同化遙感蒸散發(fā)數(shù)據(jù)模擬，兩條洪水過程曲線幾乎重合而難以分辨。

由表2可知，無論蒸散發(fā)數(shù)據(jù)采用地面站實測還是遙感蒸散發(fā)同化，驗證期水文模型對各場洪水模擬效率系數(shù)均在0.8以上。各場洪水總徑流量誤差均在20%以內(nèi)，洪峰流量誤差在15%以內(nèi)，峰現(xiàn)時間在5%以內(nèi)，且≤3 個時間段，可滿足模擬精度要求。說明所構(gòu)建水文模型模擬效果較好，所率定水文模型參數(shù)可較好描述七星河流域產(chǎn)匯流機制。因此，本文構(gòu)建基于蓄滿產(chǎn)流半分布式水文模型可用于次洪模擬。

對比表2中利用地面蒸發(fā)皿估算蒸散發(fā)數(shù)據(jù)洪水模擬效果和利用同化后遙感蒸散發(fā)數(shù)據(jù)洪水模擬精度發(fā)現(xiàn)，利用同化后遙感蒸散發(fā)數(shù)據(jù)替代地面蒸發(fā)皿實測估算蒸散發(fā)數(shù)據(jù)，水文模型對洪水No.120728總徑流模擬精度提高0.3%，洪峰模擬精度提高0.3%；120812 洪水總徑流模擬精度提高0.2%，洪峰模擬精度提高0.1%，130629洪水效率系數(shù)提高0.01，總徑流模擬精度提高4%，洪峰模擬精度提高1.6%，130722洪水效率系數(shù)提高0.01，洪峰模擬精度提高0.8%。上述分析發(fā)現(xiàn)：①同化蒸散發(fā)數(shù)據(jù)可較好代表蒸發(fā)站實測數(shù)據(jù)；②與傳統(tǒng)利用地面蒸發(fā)皿估算蒸散發(fā)數(shù)據(jù)相比，利用集合卡爾曼濾波同化后蒸散發(fā)數(shù)據(jù)作洪水模擬，所有模擬精度指標均有提升，尤其是總徑流模擬精度提高最為明顯。

4 結(jié) 論

本研究利用遙感蒸散發(fā)產(chǎn)品SSEBop 同化地面測站蒸發(fā)數(shù)據(jù)，結(jié)合基于蓄滿產(chǎn)流半分布式水文模型研究遙感蒸散發(fā)產(chǎn)品在次洪模擬中應(yīng)用效果。半分布式水文模型可較好同時考慮大尺度流域蒸散發(fā)和下墊面空間不均勻性，對降雨徑流過程有效模擬。遙感蒸散發(fā)產(chǎn)品可獲得大空間尺度面蒸散發(fā)數(shù)據(jù)，可克服傳統(tǒng)地面蒸發(fā)數(shù)據(jù)受地面蒸發(fā)站點密度影響缺陷，并有效考慮蒸散發(fā)空間不均勻性。

由于降雨等輸入、水文模型結(jié)構(gòu)存在不確定性，本文提出次洪模擬方法有待進一步改進。例如，利用土壤含水量遙感反演替代本文的自由初始徑流量估算，提高水文模型初始條件精度。此外，可嘗試利用遙感降水產(chǎn)品替代本文面雨量數(shù)據(jù)，更好契合半分布式水文模型中下墊面不均勻情景設(shè)計，充分考慮降雨空間不均勻性。