尹 劍，邱遠宏，張 斌

(1. 貴州財經(jīng)大學(xué)西部現(xiàn)代化研究中心，貴州 貴陽 550025；2. 貴州財經(jīng)大學(xué)大數(shù)據(jù)應(yīng)用與經(jīng)濟學(xué)院，貴州 貴陽 550025； 3.東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，黑龍江 哈爾濱 150030)

蒸散發(fā)(Evapotranspiration，ET)是流域水文過程中關(guān)鍵且較難估算的一個分量[1]。遙感模型具有高時效特征，可得出較為準(zhǔn)確的蒸散發(fā)結(jié)果[2]。但是可見光-熱紅外測量易受天氣影響，難以獲取長時間序列連續(xù)的地表蒸散發(fā)[3-4]。水文模型以閉合流域的整個水循環(huán)系統(tǒng)為研究對象，可以模擬時間上連續(xù)的產(chǎn)匯流和地表蒸散發(fā)等水文變量[5]。水文模型等存在的問題是精度會隨著模型的向前模擬產(chǎn)生誤差積累[6]。如果把精度較高的遙感ET和水文模型結(jié)合起來實現(xiàn)優(yōu)勢互補，對獲得高精度且時間連續(xù)的蒸散發(fā)有重要的意義[7]。數(shù)據(jù)同化作為一種可以將多元數(shù)據(jù)融合而達到概率最優(yōu)解的理想算法，為這一目的提供了思路[8]。它將精度較高的蒸散發(fā)遙感反演結(jié)果作為觀測信息，調(diào)整水文模型的運行，使蒸散發(fā)水文模擬結(jié)果不斷地向觀測值靠攏，以減小誤差積累，進而校正水文模型[9]。本研究采用數(shù)據(jù)同化算法，構(gòu)建一個基于遙感蒸散發(fā)和水文模型的同化系統(tǒng)，并從流量精度和蒸散發(fā)時空格局分析同化效應(yīng)，以驗證數(shù)據(jù)同化對水文過程特別是對蒸散發(fā)過程模擬的影響。

1 研究方法

同化的流程有5個部分：①準(zhǔn)備精度較高的模型驅(qū)動數(shù)據(jù)；②利用觀測數(shù)據(jù)生成經(jīng)過真實性檢驗的觀測產(chǎn)品；③將驅(qū)動數(shù)據(jù)輸入水文模型，生成當(dāng)前時刻的水文狀態(tài)參數(shù)集；④利用觀測產(chǎn)品同化當(dāng)前時刻的水文模型，優(yōu)化狀態(tài)變量；⑤水文模型向前積分，直至下一個觀測時再次同化。

數(shù)據(jù)同化系統(tǒng)中關(guān)鍵的3個部分是模型算子、觀測以及同化算法。模型算子選擇日尺度的DTVGM水文模型[10-11]。觀測選擇遙感日蒸散發(fā)。對于同化算法，要考慮計算效率和同化性能這對矛盾，以實現(xiàn)效率和精度的統(tǒng)一，為避免濾波發(fā)散，選擇確定性集合卡爾曼濾波(Deterministic Ensemble Kalman Filter，DEnKF)算法[12]。DEnKF由Sakov和OKE[13]提出，是一種不需要加入觀測擾動的確定性同化方法，魯棒性較強，可在相對較小的集合下表現(xiàn)出良好的同化效果[14]，其更新過程為：

(1)

(2)

(3)

(4)

DEnKF同化方法，通過在分析擾動更新階段將卡爾曼增益矩陣減半的方式維持集合傳播：

(5)

2 研究區(qū)及數(shù)據(jù)

研究選擇北京沙河流域為實驗區(qū)，流域隸屬海河流域，分布在北緯40°00′～40°30′和東經(jīng)115°50′～116°20′。流域內(nèi)海拔落差1 309 m，面積約為1 138 km2。沙河流域是北京市重要的林業(yè)區(qū)和農(nóng)作物區(qū)，林地覆蓋率超過50%，耕地約占35%。流域內(nèi)植被覆蓋不均勻，四季分明，下墊面在時空上具有一定的分異性。流域共28個子流域，有5個主要水庫參與水資源調(diào)控，見圖1。近年來流域降雨量持續(xù)減少，出現(xiàn)干旱趨勢。準(zhǔn)確估算流域地表蒸散發(fā)變化過程，對流域水循環(huán)研究、應(yīng)對地區(qū)旱情、提高水資源利用效率具有重要的應(yīng)用價值。

圖1 研究區(qū)概況

考慮到研究區(qū)空間尺度，遙感數(shù)據(jù)采用Landsat影像。遙感蒸散發(fā)模型選擇一種雙層ET遙感模型反演日尺度ET[15]。遙感ET基于Landsat衛(wèi)星影像反演獲得，選取了1999—2007年云遮蓋小于10%的晴好日的TM/ETM數(shù)據(jù)，軌道號為123/32。從1999年7月1日起開始有成像質(zhì)量較好的TM影像，因此從該日遙感計算ET并導(dǎo)入同化系統(tǒng)。水文過程模擬通過沙河閘實測流量做驗證。選取納西效率系數(shù)(NES)和水量平衡系數(shù)(WB)作為精度評定指標(biāo)。DTVGM的參數(shù)采用SCE-UA算法率定[16]。

3 結(jié)果分析

3.1 流量模擬驗證

DTVGM和同化2種情況下，1999—2007年不同年份的WB和NES值見表1。為了更清楚地比較，圖2繪制了基于觀測和模擬的日均流量和降雨過程。可見，大部分年份模擬與觀測擬合效果較好。2種情況下NES均值大于0.75，一般滿足水文模擬的精度要求[17]。同化后NES從0.76提升到0.79。DTVGM直接模擬中NES大于0.75的年份比為55.6%，同化結(jié)果占比為77.8%。2種情景下的WB差異是明顯的。同化后的WB是1.09，比DTVGM同化前的1.20要好。也就是說，直接模擬的相對誤差約為20%，同化的相對誤差小于10%。DTVGM直接模擬，水量平衡相對誤差小于20%的年份的百分比為33.3%，而在同化后，該比例上升到88.9%，說明同化改善了徑流過程的模擬。

表1 同化前后模擬評價指標(biāo)比較

圖2 流量過程模擬結(jié)果比較

整體上來看，除了2003年優(yōu)化得到的模擬流量與實測流量差異比較明顯，多數(shù)年份的擬合效果滿足日尺度模擬精度的要求。從圖2可以看出，2003年春、夏季降水較少，而春、夏季又是農(nóng)作物需水量最大的季節(jié)，為了保證灌溉，流域內(nèi)水庫關(guān)閘截流。人工調(diào)控較為顯著地影響了徑流過程，導(dǎo)致實測流量與自然狀態(tài)下的降雨徑流關(guān)系不一致。雖然水文模型一定程度上整合了人類用水和地下水模塊，但由于數(shù)據(jù)的不完善，在降雨變異性較大的年份，模擬精度受到影響。但數(shù)據(jù)同化借助實時遙感信息，在一定程度上對該影響進行了校準(zhǔn)，2003年的WB和NES分別由同化前的1.52和0.44提升為同化后的1.23和0.66。

3.2 流域蒸散發(fā)過程分析

3.2.1蒸散發(fā)同化結(jié)果驗證

研究區(qū)1號子流域(圖1)在2002—2007年開展了連續(xù)的水熱通量觀測[7,15]。因此采用該觀測站擬合的逐日蒸散發(fā)與1號子流域的模擬日蒸散發(fā)進行比較，時間區(qū)段為2002年1月1日至2007年12月31日，分析同化前后的精度。圖3表示了同化前使用DTVGM模擬輸出的蒸散發(fā)與實測值的對比，圖4表示了同化系統(tǒng)的蒸散發(fā)輸出值與實測值的對比。同化前后決定系數(shù)R2分別為0.814 9、0.937 4，且同化后趨勢線斜率更接近1，表明經(jīng)過數(shù)據(jù)同化ET模擬精度更加接近實測值。

圖3 蒸散發(fā)的DTVGM模擬值與實測值對比

圖4 蒸散發(fā)同化系統(tǒng)輸出值與實測值對比

3.2.2蒸散發(fā)過程分析

選擇1999—2007年的逐月蒸散發(fā)進行對比，圖5為同化系統(tǒng)最終生成的1999—2007年的月蒸散量、潛在蒸散量(ETp)、降雨量變化軌跡圖以及各年匯總值。沙河流域ET年內(nèi)變化比較明顯，在不同年份的相同月份也呈現(xiàn)出差異性。各年月蒸散發(fā)過程均呈現(xiàn)出從1—12月先增大后減小的趨勢，年內(nèi)變化幅度小于潛在蒸散發(fā)。月蒸散發(fā)過程線峰值主要出現(xiàn)在7、8月份，且同一年內(nèi)這兩個月數(shù)值差異不大；月蒸散發(fā)過程線谷值出現(xiàn)在1月或12月。降雨過程線的峰值各年不同，主要落在6—8月，其峰值略提前于月蒸散發(fā)過程線的峰值，谷值與蒸散發(fā)過程線谷值一致，在降雨量增大的同月或次月蒸散發(fā)也隨之出現(xiàn)增大。潛在蒸散發(fā)過程線各年差異不是很明顯，谷值主要出現(xiàn)在1月或12月，峰值主要出現(xiàn)在5、6月。各年月蒸散發(fā)波動較大的情況主要發(fā)生在春季。潛在蒸散量和降雨量的變化過程可以反映氣候變化的影響，在兩者年際差異不是很大的情況下，如果蒸散發(fā)過程年際差異較大則說明人類活動對蒸散發(fā)產(chǎn)生了重要影響。

a) 1999年

b) 2000年

c) 2001年

d) 2002年

e) 2003年

f) 2004年圖5 1999—2007年蒸散發(fā)、潛在蒸散發(fā)和降水量變化過程

g) 2005年

h) 2006年

i) 2007年

j) 1999—2007年續(xù)圖5 1999—2007年蒸散發(fā)、潛在蒸散發(fā)和降水量變化過程

與其他年份不同，2003年月蒸散發(fā)峰值出現(xiàn)在6月，且全年月蒸散發(fā)波動幅度較小。對比降雨過程線發(fā)現(xiàn)，2003年春季，除了3月份有小幅降水外，4、5月的降雨量均較其他年同期偏小，尤其是4月份降雨幾乎為0，而流域蒸散發(fā)并沒有相應(yīng)地發(fā)生較大的偏差。4月份是農(nóng)業(yè)灌溉的主要月份，在降雨量低的情況下，為了緩解區(qū)域旱情，沙河流域各水庫提高了蓄水量，進行增加灌溉補充。這一結(jié)論與流域在2003年全年降水與其他年相近而流量卻較小有著一致性。

圖6對比了月蒸散發(fā)過程的峰谷值。各年月蒸散發(fā)峰值變化不明顯，對應(yīng)月份的降雨量波動較大，流域出口流量與降雨量有著一定的正相關(guān)關(guān)系，截流保證了流域需水，進而表現(xiàn)出各年蒸散發(fā)峰值未出現(xiàn)較大差異。相比其他年月蒸散發(fā)峰值都大于80 mm的情況，2003年的蒸散發(fā)峰值最小，僅有69.23 mm，且時間提前到6月，與該年春夏季降雨量較其他年同期偏小有關(guān)。經(jīng)計算，各年沙河流域月蒸散發(fā)峰值與降雨峰值之間的相關(guān)系數(shù)為0.186，與流量相關(guān)系數(shù)為0.534。因此得出結(jié)論，流域月蒸散發(fā)的峰值受流量的影響較大，受降雨的影響較小，加之沙河流域各年蒸發(fā)能力變化不明顯，從而蒸散發(fā)受潛在蒸發(fā)影響也較??；由于流量是受流域水庫閘壩控制的變量，沙河流域各年月蒸散發(fā)峰值受流域內(nèi)人類活動對水資源的調(diào)節(jié)影響大于氣候因子的影響。

圖6 1999—2007年沙河流域月蒸散發(fā)峰值、谷值及峰值出現(xiàn)月的降雨、流量對比

3.2.3典型年蒸散發(fā)分析

選擇2003年為代表年，分析該年各子流域蒸散發(fā)變化特征。2003年各月子流域平均蒸散發(fā)見圖7。

圖7、8分別顯示了2003年流域蒸散發(fā)的空間和時間的變化過程。前文分析2003年的春夏進行了水庫蓄水，補給水庫上游水資源需求。從圖7可以看出，蒸散發(fā)的空間分布在多數(shù)月份均呈現(xiàn)西北部上游山區(qū)大于下游平原區(qū)的特點。在1—3月上下游蒸散發(fā)空間差異明顯，3月蒸散發(fā)上升很快，可能與3月有一次比往年偏大的降雨有關(guān)。2003年4月流域降雨量幾乎為零，而4月又是農(nóng)業(yè)需水量大的月份，流域水庫開始截流蓄水用于補給灌溉，比較明顯的有上莊水庫、沙河閘以及十三陵水庫。水庫周邊子流域以及相鄰春種作物面積較大的子流域月蒸散量大于周圍區(qū)域。進入5月，流域降雨量依然較少，從圖中可以看出王家園水庫相關(guān)子流域蒸散發(fā)大于周邊地區(qū)，可能是由于該水庫調(diào)水補充周邊區(qū)域?qū)е碌?。同理也可分析推?—8月流域內(nèi)夏種作物分布區(qū)蒸散發(fā)較大的原因是區(qū)域灌溉導(dǎo)致的。9月降雨量開始同比回升，沙河閘正常開閘放水。整個秋季，流域蒸散發(fā)上游山區(qū)大于下游平原區(qū)的空間分布特征又得以出現(xiàn)。

a) 1月

d) 4月

g) 7月

b) 2月

e) 5月

h) 8月

c) 3月

f) 6月

i) 9月

j) 10月

k) 11月

圖8 2003年數(shù)據(jù)蒸散發(fā)計算結(jié)果(5日移動平均)

考察月蒸散發(fā)過程發(fā)現(xiàn)，在遙感觀測輸入后，數(shù)據(jù)同化修正了水文模型模擬春季蒸散發(fā)出現(xiàn)的偏差。對比蒸散發(fā)同化結(jié)果和氣象、水文實測數(shù)據(jù)發(fā)現(xiàn)，水文模型對蒸散發(fā)的模擬主要是沒有準(zhǔn)確反映人類活動，特別是在春旱出現(xiàn)的2003年。數(shù)據(jù)同化一定程度上校正了人類活動對水文模型模擬精度的影響，提高了模型對流域情景的反映，保證了流域蒸散發(fā)過程的高精度模擬。

4 結(jié)語

研究構(gòu)建了基于確定性數(shù)據(jù)同化的沙河流域水文模型，從流量模擬驗證和蒸散發(fā)模擬驗證分析同化效應(yīng)。結(jié)果發(fā)現(xiàn)數(shù)據(jù)同化系統(tǒng)改善了徑流過程和蒸散發(fā)過程的模擬精度。通過分析典型年蒸散發(fā)時空，看出同化可以改善水文模型對蒸散發(fā)時空分布的刻畫，同時可以在一定程度上削弱資料不完備地區(qū)的模擬偏差，如灌溉用水資料缺乏對水文模擬的影響。數(shù)據(jù)同化可提高水文模型對流域情景的實時反映能力，從而實現(xiàn)模擬校正。