遙感蒸散發(fā)在無(wú)測(cè)站資料地區(qū)洪水模擬中的應(yīng)用

邢貞相 ，傅 爽，孫明新，紀(jì) 毅 ，付 強(qiáng) *，李 衡

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

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

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

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

七星河流域位于中國(guó)黑龍江省三江平原東北部，東經(jīng)131°9′～131°44′、北緯46°12′～46°38′。總面積為1 344 km2，其中山地占總面積38.3%，主要分布于流域西南部和南部；平原占總面積61.7%，主要分布于流域北部和中部。七星河流域?qū)俸疁貛Т箨懶詽駶?rùn)季風(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個(gè)雨量站、1個(gè)水文站（同時(shí)監(jiān)測(cè)降雨、徑流），平均每268.8 km2布設(shè)一個(gè)站點(diǎn)。以上各站可提供2006～2013 年10 場(chǎng)實(shí)測(cè)場(chǎng)次降水?dāng)?shù)據(jù)。采用反距離權(quán)重法（IDW）[12]對(duì)雨量數(shù)據(jù)作空間插值，得到七星河流域30′′×30′分辨率小時(shí)降水插值場(chǎng)。流量數(shù)據(jù)采用該流域出口斷面保安水文站2006～2013年間實(shí)測(cè)徑流數(shù)據(jù)。由于流域內(nèi)無(wú)蒸發(fā)觀測(cè)站，本文移用與七星河流域具有相近緯度和氣象條件寶清站（132°15′E 46°20′N(xiāo)）實(shí)測(cè)蒸發(fā)資料作為七星河流域蒸發(fā)數(shù)據(jù)。遙感蒸散發(fā)數(shù)據(jù)選取SSEBop 月尺度蒸散發(fā)產(chǎn)品[13]。該蒸散發(fā)產(chǎn)品通過(guò)operational Simplified Surface Energy Balance（SSEBop）模型產(chǎn)生[14]。SSEBop 使用基于簡(jiǎn)化表面能量平衡（SSEB）方法，結(jié)合遙感MODIS 熱圖像生成Evapotranspiration（ET）（每10 d 采集1 次）和參考ET 熱 指 數(shù)[15-16]。最 初SSEB 公 式 基 于SEBAL 和METRIC 模型冷與熱像素原理[17-18]。該蒸散發(fā)數(shù)據(jù)空間分辨率為1 km。此外，本研究中還采用美國(guó)USGS提供30〃×30〃分辨率DEM（見(jiàn)圖1）。

2 研究方法

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

2.1 集合卡爾曼濾波方法

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

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

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

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

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

水流在河道中運(yùn)動(dòng)速度由Manning 公式（2）計(jì)算：

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

將每個(gè)柵格出流過(guò)程累加后，可得流域出口徑流過(guò)程。具體細(xì)節(jié)參見(jiàn)文獻(xiàn)[22]。

2.3 參數(shù)率定方法

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

3 結(jié)果與分析

3.1 水文模型率定

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

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

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

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

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

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

土壤初始含水量與河道徑流量擬合曲線方程為式（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ā)器測(cè)量值乘以蒸散發(fā)折算系數(shù)1（該系數(shù)基于以往經(jīng)驗(yàn)值得出），得出實(shí)際地面蒸散發(fā)。在月尺度和每個(gè)30′′×30′柵格上，將臨近蒸發(fā)站寶清站測(cè)量值作為模型算子，以SSEBop 遙感蒸散發(fā)模型反演蒸散發(fā)作為觀測(cè)數(shù)據(jù)，輸入集合卡爾曼濾波算法中同化計(jì)算，流域出口寶安站所在柵格計(jì)算結(jié)果如圖4所示。由圖4可知，同化后蒸散發(fā)量更接近于地面蒸發(fā)站實(shí)測(cè)數(shù)據(jù)，說(shuō)明EnKF可有效同化遙感蒸散發(fā)數(shù)據(jù)。

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

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

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

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

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

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

對(duì)比表2中利用地面蒸發(fā)皿估算蒸散發(fā)數(shù)據(jù)洪水模擬效果和利用同化后遙感蒸散發(fā)數(shù)據(jù)洪水模擬精度發(fā)現(xiàn)，利用同化后遙感蒸散發(fā)數(shù)據(jù)替代地面蒸發(fā)皿實(shí)測(cè)估算蒸散發(fā)數(shù)據(jù)，水文模型對(duì)洪水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í)測(cè)數(shù)據(jù)；②與傳統(tǒng)利用地面蒸發(fā)皿估算蒸散發(fā)數(shù)據(jù)相比，利用集合卡爾曼濾波同化后蒸散發(fā)數(shù)據(jù)作洪水模擬，所有模擬精度指標(biāo)均有提升，尤其是總徑流模擬精度提高最為明顯。

4 結(jié) 論

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

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