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

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

(武漢大學(xué)水資源與水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

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

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

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

1 方法與試驗(yàn)

1.1 SWAT模型

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

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

在進(jìn)行降雨-徑流計(jì)算時(shí)徑流流量由3部分組成[12]：

Q=Qsurf+Qlat+Qgw

(1)

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

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

(2)

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

(3)

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

HRU計(jì)算單元的水量平衡總方程為：

(4)

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

1.2 集合卡爾曼濾波

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

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

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

(5)

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

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

yt=[mTt,uTt,dTt]

(6)

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

其中觀測(cè)值可表示為：

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

(7)

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

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

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

(8)

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

Kt=CmHT[HCmHT+CD]-1

(9)

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

2 算例研究

2.1 算例說明

為了驗(yàn)證徑流、土壤水含量及地下水位等觀測(cè)信息的作用，本文構(gòu)造了一個(gè)虛擬的流域水文試驗(yàn)，使用漳河灌區(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個(gè)子流域(見圖1)，14個(gè)水文響應(yīng)單元(HRU)。徑流、土壤水含量及地下水位觀測(cè)數(shù)據(jù)由已知參數(shù)的參照模型運(yùn)行得到。

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

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

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

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

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

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

為了評(píng)價(jià)同化效果，使用均方根誤差RMSE來量化模擬值與真實(shí)值之間的差異：

(10)

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

3 分析及討論

3.1 徑流估計(jì)

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

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

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

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

圖3 加入不同觀測(cè)時(shí)流域徑流的RMSEFig.3 RMSE series of streamflow when adding different observations

3.2 土壤水估計(jì)

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

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

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

3.3 CN2估計(jì)

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

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

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

4 結(jié) 語

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

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

□

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

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

[3] 芮孝芳. 流域水文模型研究中的若干問題[J]. 水科學(xué)進(jìn)展, 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模型在我國(guó)的研究和應(yīng)用進(jìn)展[J]. 中國(guó)農(nóng)村水利水電，2015,(5):109-113.

[14] 王中根, 劉昌明, 黃友波. SWAT 模型的原理, 結(jié)構(gòu)及應(yīng)用研究[J]. 地理科學(xué)進(jìn)展, 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.