多源觀測在水文數(shù)據(jù)同化中的應(yīng)用

楊楚慧，林 琳，張秋汝，史良勝

(武漢大學(xué)水資源與水電工程國家重點實驗室，武漢 430072)

流域水文模型被廣泛應(yīng)用于灌區(qū)管理和洪澇災(zāi)害防治等方面[1,2]。然而由于模型誤差和觀測誤差等不確定性的影響，水文模型在進行徑流等水文狀態(tài)預(yù)報時存在不可避免的偏差[3]。因此，多目標函數(shù)自動調(diào)參、數(shù)據(jù)同化等方法[4,5]被用于修正誤差來改善水文狀態(tài)的估計。其中數(shù)據(jù)同化(data assimilation)由于綜合考慮了各種誤差來源，而被廣泛應(yīng)用于水文模型的參數(shù)校正和狀態(tài)預(yù)測。水文數(shù)據(jù)同化通過向模型中連續(xù)加入觀測數(shù)據(jù)來降低不確定性，最終提高模型預(yù)測能力。隨著觀測手段的迅速發(fā)展，數(shù)據(jù)同化技術(shù)在近年來得到了深入研究。

在水文數(shù)據(jù)同化模型中，研究者大多采用徑流和土壤水等觀測數(shù)據(jù)來進行流域水文狀態(tài)的更新和參數(shù)估計。Xie和Zhang[6]利用徑流觀測準確地預(yù)測了徑流及相關(guān)參數(shù)。Chen等[7]利用表層土壤含水量觀測進行同化，發(fā)現(xiàn)深層土壤水和徑流的估計值與實測有明顯的偏差。Sun和Lee等[8,9]聯(lián)合利用徑流和土壤水觀測進行同化，發(fā)現(xiàn)徑流和土壤水狀態(tài)的估計結(jié)果優(yōu)于同化其中任一觀測。此外，Shi等[10,11]用數(shù)值和實例驗證了同時加入徑流、地下水位、土壤水、以及蒸散發(fā)等多種觀測對多參數(shù)估計的有效性，但不同的觀測組合對于水文狀態(tài)預(yù)測的影響仍有待討論。其中地下水位作為一種簡單經(jīng)濟的觀測手段，受到上層土壤水補給和蒸發(fā)等因素的影響，能在一定程度上反映非飽和土壤水的變化情況，因此研究同化地下水位對于徑流和土壤水的影響過程具有一定價值。

為了進一步研究多源觀測對于水文數(shù)據(jù)同化的影響，本文使用集合卡爾曼濾波(EnKF)方法，與分布式水文模型SWAT結(jié)合，研究徑流、土壤水和地下水位等觀測在參數(shù)估計和狀態(tài)預(yù)測方面的性能。

1 方法與試驗

1.1 SWAT模型

SWAT(Soil and Water Assessment Tool)模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心(USDA)開發(fā)的半分布式水文模型。模型可以模擬多種水文循環(huán)物理過程，進行在長時間序列下不同土壤類型、土地利用和管理條件的復(fù)雜流域水文計算[12]，在我國被廣泛應(yīng)用于水文模擬和非點源模擬等方面[13]。

在建立模型過程中，基于DEM數(shù)據(jù)，流域被劃分成若干個子流域，并根據(jù)不同的下墊面條件進一步歸類成若干個水文響應(yīng)單元(HRU)[14]。SWAT模擬的逐日水文循環(huán)基于每個HRU的水量平衡方程進行，HRU是SWAT最基本的計算單元。

在進行降雨-徑流計算時徑流流量由3部分組成[12]：

Q=Qsurf+Qlat+Qgw

(1)

式中：Q為當日總徑流量；Qsurf為地表徑流量；Qlat為側(cè)向徑流量，(mmH2O)；Qgw為當日流入河道的地下水徑流量。

其中的地表徑流量一般采用SCS曲線數(shù)法來計算，其計算方程如下：

(2)

式中：P為第i天的降水量；S為滯留參數(shù)，隨土壤類型和土地利用等下墊面條件產(chǎn)生空間差異，隨土壤含水量的變化產(chǎn)生時間差異，其定義公式為：

(3)

式中：CN為日曲線數(shù)，取決于CN2和前期土壤水分條件；CN2為土壤水分條件為一般濕潤時的曲線數(shù)。

HRU計算單元的水量平衡總方程為：

(4)

式中：SWt為第i天的土壤含水量；SW0為土壤初始含水量；t為時間，d；Qsurf為地表徑流量；ETi為第i天的實際蒸散發(fā)量；Wseep為從土壤剖面進入包氣帶的水量。

1.2 集合卡爾曼濾波

數(shù)據(jù)同化通過向系統(tǒng)中加入連續(xù)的觀測來修正模型的變量狀態(tài)[15]，集合卡爾曼濾波(EnKF)是一種順序數(shù)據(jù)同化算法，由于其較低的計算成本和簡易靈活性，被廣泛應(yīng)用于水文領(lǐng)域。

水文模型在氣象因子的驅(qū)動下狀態(tài)和參數(shù)發(fā)生動態(tài)變化：

ut+1=f(ut,At+1)+wt+1

(5)

式中：u是水文狀態(tài)值如徑流、土壤水含量等；A是驅(qū)動因子如降水等；t是時間；w是誤差項，代表模型誤差、參數(shù)誤差等的影響。

進行同化的狀態(tài)向量由模型參數(shù)、狀態(tài)變量和觀測共同組成：

yt=[mTt,uTt,dTt]

(6)

式中：y是狀態(tài)向量；m是估計的模型參數(shù)，本算例中為CN2；u是模型包含的狀態(tài)變量，如徑流量Q和土壤水含量SW等；d是觀測值，包括徑流、剖面土壤水含量和地下水位；t是時間序列。

其中觀測值可表示為：

dt+1=G(mt,ut)+wt+1

(7)

式中：G(∶)表示模型狀態(tài)和觀測值之間的轉(zhuǎn)換關(guān)系；w為觀測誤差。

對于t+1時刻的樣本j，集合卡爾曼濾波的更新向量可表示為：

yat+1,j=ymt+1,j+Kj(dt+1,j-Hymt+1,j)

(8)

式中：y是狀態(tài)向量；a和m分別代表最優(yōu)估計值和預(yù)報值；d是觀測值；H是觀測矩陣，代表狀態(tài)向量與觀測向量之間的轉(zhuǎn)換關(guān)系；K是卡爾曼增益，表示在同化過程中觀測和狀態(tài)值的權(quán)重：

Kt=CmHT[HCmHT+CD]-1

(9)

式中：Cm是狀態(tài)向量的協(xié)方差矩陣；CD是觀測誤差矩陣；CmHT表示狀態(tài)向量與觀測向量的互協(xié)方差；HCmHT表示觀測向量的自協(xié)方差，可以直接利用樣本統(tǒng)計的方法得到。

2 算例研究

2.1 算例說明

為了驗證徑流、土壤水含量及地下水位等觀測信息的作用，本文構(gòu)造了一個虛擬的流域水文試驗，使用漳河灌區(qū)的下墊面和氣象資料，但調(diào)整了土壤類型及其分布。研究區(qū)總面積 1 089.7 km2，主要土地利用方式為落葉闊葉林(FRSD)，主要土壤類型為黏土(56.57%)和黏壤土(43.43%)。建立模型的主要輸入數(shù)據(jù)包括數(shù)字高程地形圖(DEM)、土地利用及土壤類型圖、日降水資料和其他氣象數(shù)據(jù)。

在SWAT模型中輸入以上數(shù)據(jù)后，根據(jù)DEM信息、土壤類型和土地利用方式，將流域劃分成7個子流域(見圖1)，14個水文響應(yīng)單元(HRU)。徑流、土壤水含量及地下水位觀測數(shù)據(jù)由已知參數(shù)的參照模型運行得到。

圖1 研究區(qū)域DEM和子流域劃分Fig.1 DEM and watershed delineation of research area

2.2 數(shù)據(jù)同化過程

整個模擬過程分為3個步驟：①預(yù)熱期，在長期氣象驅(qū)動下使模型達到合理的水文狀態(tài)，運行時間為2001年1月1日至2004年1月1日；②擾動期，通過正態(tài)分布對參數(shù)進行擾動生成若干個樣本，每一個樣本獨立向前演算，運行時間為2004年1月1日至2004年4月9日；③同化期，連續(xù)性地加入觀測值來更新水文狀態(tài)(徑流和表層土壤含水量等)和模型參數(shù)CN2，從而得到更為準確的估計值，運行時間為2004年4月10日至2004年11月1日。

本文設(shè)置了4個算例。為了考慮影響徑流估計的因素，設(shè)置算例1：僅加入流域出口處徑流觀測(見圖1)；算例2：僅加入流域表層5 cm土壤水觀測；算例3：綜合考慮以上2個算例，聯(lián)合利用徑流和表層含水量觀測，來討論多源觀測對于提高模型狀態(tài)描述的效果；算例4：同時加入徑流和地下水位觀測，來探討地下水位對于反演土壤水狀態(tài)的作用。

在上述算例中，所有算例估計的水文狀態(tài)量都相同，包括徑流Q和表層含水量SW等，不確定參數(shù)僅考慮CN2，其初始分布服從G(60,52)，所有觀測的頻率為每日一次，徑流觀測相對誤差為10%，土壤水觀測相對誤差為3%，地下水觀測絕對誤差為2 cm，土壤水和地下水模型誤差為5%，計算樣本共300個。

本文給出了參照解、EnKF解和非條件解來驗證數(shù)據(jù)同化的效果。參照解是指在正確的參數(shù)和邊界條件下的水文狀態(tài)；EnKF解是在錯誤的初始參數(shù)條件下利用數(shù)據(jù)同化模型計算得到的結(jié)果；非條件解的初始參數(shù)和邊界條件與EnKF解相同，但不進行參數(shù)更新和數(shù)據(jù)同化時得到的結(jié)果。

為了評價同化效果，使用均方根誤差RMSE來量化模擬值與真實值之間的差異：

(10)

式中：n為計算單元數(shù)目，等于子流域或者HRU的個數(shù)；Y(xi)為參照值；Yt(xi)為同化后的樣本均值。

3 分析及討論

3.1 徑流估計

圖2描繪了加入不同觀測時子流域7徑流的時間變化，同時給出了參照解和非條件解。由圖2可以看出，當加入的觀測包含徑流時，徑流估計值與參照值十分接近；僅使用土壤水觀測時，僅加入表層土壤水觀測在雨季能較好地預(yù)測徑流值，旱季時預(yù)測的徑流與參照值有較大的偏差。

圖2 子流域7加入不同觀測時徑流的時間變化 Fig.2 Estimated streamflow process in subbasin 7 when adding different observations

圖3給出了整個流域徑流RMSE的時間變化。圖3中在加入徑流觀測后，無論是否加入其他類型觀測，流域徑流估計的RMSE都能迅速減小到接近0值；僅加入淺層土壤水時雨季的RMSE迅速降低，在旱季則會造成徑流估計的惡化，RMSE值甚至大于非條件解。因此，僅需加入徑流觀測就能迅速校正徑流預(yù)測值，而表層土壤含水量只能反映該子流域地表產(chǎn)流的過程，而無法反映基流過程，因此無法改善徑流的估計效果。

由于流域出口處徑流包含了所在子流域的產(chǎn)流信息和整個流域的匯流信息，可以利用該點的徑流值反演上游各子流域的徑流過程。因此，在估計流域徑流時，對出口點的徑流進行觀測十分必要。算例1中雖然僅加入了流域出口處的徑流觀測，但數(shù)據(jù)同化后整個流域的徑流預(yù)測都獲得了改善，因此流域出口處徑流對于整個流域的徑流估計具有重要的價值，Xie[6]等人在研究中也提出了相似的觀點。

圖3 加入不同觀測時流域徑流的RMSEFig.3 RMSE series of streamflow when adding different observations

3.2 土壤水估計

圖4為加入不同觀測時子流域7中HRU 14的各層土壤含水量的RMSE變化情況。如圖4(a)所示，利用徑流觀測可以較好地估計出表層土壤水，這是由于表層土壤含水量對徑流產(chǎn)流過程有著強烈的影響，此時2者相關(guān)性較高。但隨著土壤深度的增加，相關(guān)性的減弱，加入徑流觀測對深層土壤水估計效果的產(chǎn)生明顯的偏差，且變化較為劇烈[見圖4(b)和4(c)]。總體而言，加入徑流觀測對土壤剖面整體的含水量估計作用較弱。

僅利用表層土壤水觀測時，深層土壤水估計會產(chǎn)生偏差，使得土壤水估計會發(fā)生劇烈變化[見圖4(b)和4(c)]，此時的RMSE也變化十分迅速。同時加入徑流和表層土壤水觀測能從一定程度上調(diào)整深層土壤水的偏離情況，圖4中該算例的RMSE值在大部分情況下最小，說明同時加入這2種觀測能調(diào)整模型偏差，使狀態(tài)估計更加準確。

圖4(c)中加入地下水和徑流觀測后的深層土壤水量RMSE也略小于僅加入徑流觀測，因此額外的地下水位數(shù)據(jù)能在一定程度上提高土壤含水量的估計效果。在SWAT中，補給和蒸發(fā)會影響地下水位的變化：在降雨較為強烈和持續(xù)的時候，土壤含水量的連續(xù)增加會導(dǎo)致地下水位的上升；土壤含水量的持續(xù)減小則會引起毛細作用的增強，并進一步帶來地下水位的降低。因此，地下水位作為一種廉價的觀測信息，可以用于改善土壤含水量估計。Shi等[11]研究發(fā)現(xiàn)，地下水位數(shù)據(jù)還可為土壤水參數(shù)的估計提供有用的信息。

3.3 CN2估計

圖5分別使用估計值和參照值作為橫縱坐標，反映所有HRU的CN2在同化期結(jié)束時的估計情況，若散點越接近45°線，則效果越好。從圖5(b)、5(d)和5(e)中可以發(fā)現(xiàn)，加入徑流的算例的CN2估計效果都比較均勻的分布在45°線左右，這是由于地表產(chǎn)流和CN2有強烈的相關(guān)性，因此僅加入流域出口點徑流能較為準確的描述整個流域的徑流過程，并獲取合理的CN2參數(shù)。而本算例中沒有考慮誤差和觀測點數(shù)目等因素的影響，可以在后續(xù)工作中進一步討論。

而圖5(c)使用僅表層5 cm土壤水觀測來同化CN2值時，參數(shù)估計值普遍小于參照值，因此僅利用表層土壤水不能準確的估計流域參數(shù)。

圖5 CN2模擬值與參照值的比較Fig.5 Final results of CN2 final assimilation results compared to reference values

4 結(jié) 語

本文建立了基于分布式水文模型SWAT的數(shù)據(jù)同化方法，構(gòu)造虛擬的流域水文模型，討論了不同的觀測數(shù)據(jù)對于水文變量預(yù)測和參數(shù)估計的價值。研究結(jié)果表明：①僅加入徑流觀測能獲得較為準確的徑流和CN2參數(shù)估計，但是深層土壤水估計會發(fā)生偏差；②加入表層含水量觀測能模擬雨季的地表徑流變化情況，但是利用表層土壤含水量估計深層的水分變化情況會產(chǎn)生較大的偏差，從而導(dǎo)致地下水補給和基流估計十分不準確，因此無法反映整個徑流過程；③同時加入表層含水量和徑流模型對狀態(tài)的估計效果要略優(yōu)于加入其中任何一種觀測；④由于地下水位變化受到補給和蒸發(fā)的影響，地下水觀測在一定程度上能反映土壤剖面含水量的變化趨勢，從而提高土壤水的估計效果。

對于水文數(shù)據(jù)同化問題，若只需要了解徑流過程，可以只加入徑流觀測進行同化，僅提供表層土壤水分信息無法獲取根區(qū)含水量情況。成本低、易獲取的地下水觀測可以從一定程度上反映深層土壤水變化趨勢。

□

[1] 芮孝芳, 蔣成煜, 張金存. 流域水文模型的發(fā)展[J]. 水文, 2006,26(3):22-26.

[2] 代俊峰, 崔遠來. SWAT模型及其在灌區(qū)管理中的應(yīng)用前景[J]. 中國農(nóng)村水利水電, 2006,(6):34-36.

[3] 芮孝芳. 流域水文模型研究中的若干問題[J]. 水科學(xué)進展, 1997,8(1):94-98.

[4] Muleta M K, Nicklow J W. Sensitivity and uncertainty analysis coupled with automatic calibration for a distributed watershed model[J]. Journal of Hydrology, 2005,306(1):127-145.

[5] Lei F, Huang C, Shen H, et al. Improving the estimation of hydrological states in the SWAT model via the ensemble Kalman smoother: synthetic experiments for the Heihe River Basin in northwest China[J]. Advances in Water Resources, 2014,67:32-45.

[6] Xie X, Zhang D. Data assimilation for distributed hydrological catchment modeling via ensemble Kalman filter[J]. Advances in Water Resources, 2010,33(6):678-690.

[7] Chen F, Crow W T, Starks P J, et al. Improving hydrologic predictions of a catchment model via assimilation of surface soil moisture[J]. Advances in Water Resources, 2011,34(4):526-536.

[8] Sun L, Seidou O, Nistor I, et al. Simultaneous assimilation of in situ soil moisture and streamflow in the SWAT model using the extended Kalman filter[J]. Journal of Hydrology, 2016,543:671-685.

[9] Lee H, Seo D J, Koren V. Assimilation of streamflow and in situ soil moisture data into operational distributed hydrologic models: effects of uncertainties in the data and initial model soil moisture states[J]. Advances in Water Resources, 2011,34(12):1 597-1 615.

[10] Shi Y, Davis K J, Zhang F, et al. Parameter estimation of a physically based land surface hydrologic model using the ensemble Kalman filter: a synthetic experiment[J]. Water Resources Research, 2014,50(1):706-724.

[11] Shi Y, Davis K J, Zhang F, et al. Parameter estimation of a physically-based land surface hydrologic model using an ensemble Kalman filter: a multivariate real-data experiment[J]. Advances in Water Resources, 2015,83:421-427.

[12] Neitsch S, Arnold J, Kiniry J, et al. SWAT 2009 theoretical documentation[R]∥ Texas Water Resources Institute Technical Report, 2011:406.

[13] 王曉朋, 喬 飛, 雷 坤, 等.SWAT模型在我國的研究和應(yīng)用進展[J]. 中國農(nóng)村水利水電，2015,(5):109-113.

[14] 王中根, 劉昌明, 黃友波. SWAT 模型的原理, 結(jié)構(gòu)及應(yīng)用研究[J]. 地理科學(xué)進展, 2011,22(1):79-86.

[15] Evensen G. The ensemble Kalman filter: theoretical formulation and practical implementation[J]. Ocean Dynamics, 2003,53(4):343-367.