劉永偉，王 文，劉元波，劉 慶，凌 哲

(1.中國科學(xué)院南京地理與湖泊研究所流域地理學(xué)重點實驗室，江蘇 南京 210008；2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；3.江西省水利規(guī)劃設(shè)計研究院，江西 南昌 330029； 4.江蘇省水利工程規(guī)劃辦公室，江蘇 南京 210029)

數(shù)據(jù)同化是將多源觀測信息應(yīng)用于模型模擬預(yù)報的一種有效手段，被公認(rèn)是提高模型模擬預(yù)報精度、量化不確定性方面最具前途的方法[1-2]。徑流是流域產(chǎn)匯流模擬中最為重要的輸出變量，可以反映整個流域或匯水區(qū)域綜合的水文信息，且地面實測徑流數(shù)據(jù)具有很高的可靠性。因此，將站點徑流實測信息同化到水文模型模擬預(yù)報過程中，對于改善流域狀態(tài)變量及模型參數(shù)的估計精度，進(jìn)而改進(jìn)流域的降雨徑流過程模擬及預(yù)測具有重要意義。

徑流的同化始于20世紀(jì)80年代，主要是在流域/區(qū)域的水文模擬預(yù)報中基于站點徑流觀測數(shù)據(jù)進(jìn)行模型輸出結(jié)果的修正，也稱為誤差校正[3]。當(dāng)時所使用的同化方法主要為卡爾曼濾波(KF)法和擴(kuò)展KF法。隨著數(shù)據(jù)同化技術(shù)與模型模擬技術(shù)的發(fā)展，徑流的同化不僅體現(xiàn)在線性系統(tǒng)或弱非線性系統(tǒng)對系統(tǒng)狀態(tài)變量以及輸出結(jié)果的校正上，更表現(xiàn)為在非線性系統(tǒng)中通過對系統(tǒng)狀態(tài)變量及參數(shù)的更新或優(yōu)化來改進(jìn)整個系統(tǒng)的過程模擬及預(yù)報精度[4-9]。用于同化的水文模型由簡單的集總式概念性模型(如新安江模型[10]、格林-安普特模型[11])逐漸發(fā)展為基于物理機(jī)制的復(fù)雜的分布式模型(如SWAT模型[12]、TopNet模型[13])，同化方法也由線性KF法逐漸發(fā)展為各種集合KF法(如EnKF法)、粒子濾波和變分法等?；诤唵蔚募偸礁拍钚阅Ｐ偷膹搅鲾?shù)據(jù)同化方面，大量研究表明，通過徑流數(shù)據(jù)同化能夠獲得較好的模型參數(shù)與狀態(tài)估計結(jié)果，從而有效地改進(jìn)流域水文模擬的預(yù)報精度和可靠性[8，10，14]。但是，基于物理機(jī)制的分布式水文模型，由于其復(fù)雜的模型結(jié)構(gòu)、大量的模型變量與參數(shù)，導(dǎo)致其徑流數(shù)據(jù)同化過程中存在高維度現(xiàn)象，大大增加了數(shù)據(jù)同化效果的不確定性及同化實施的難度[12-13，15-16]。分布式水文模型中，徑流數(shù)據(jù)的同化效果更容易受到模型的物理機(jī)制、模型結(jié)構(gòu)、模型對流域空間異質(zhì)性的概化能力以及同化系統(tǒng)規(guī)則等方面的影響。目前，分布式水文模型中徑流數(shù)據(jù)的同化研究，基于徑流觀測數(shù)據(jù)進(jìn)行模型狀態(tài)更新與參數(shù)估計方面都相對滯后，主要表現(xiàn)在對分布式水文模型模擬誤差與徑流觀測誤差的合理量化[7，13，17]、對模型參數(shù)的同化處理(包括參數(shù)在同化過程中的演變以及參數(shù)過擬合問題的處理等)[5，15]等方面。

因此，本文在分布式水文模型SWAT中開展基于集合卡爾曼濾波(EnKF)法的徑流數(shù)據(jù)同化研究，在合理量化模型模擬誤差和徑流觀測誤差、有效處理參數(shù)演變及過擬合問題的基礎(chǔ)上，構(gòu)建基于數(shù)據(jù)同化的水文模型參數(shù)優(yōu)化方案，綜合評估徑流數(shù)據(jù)同化對SWAT模型參數(shù)的優(yōu)化效果。

1 基于EnKF法和SWAT模型的數(shù)據(jù)同化方案

1.1 SWAT模型

SWAT模型是由美國農(nóng)業(yè)部農(nóng)業(yè)研究中心(ARG-USDA)開發(fā)的分布式流域尺度水文模型[18]，主要包含水文過程子模型、土壤侵蝕子模型和污染負(fù)荷子模型。SWAT模型水文過程子模型涉及的關(guān)鍵模塊有地表徑流模塊、蒸散發(fā)模塊、土壤水分運動模塊、地下徑流模塊和河道匯流模塊。地表徑流的模擬一般采用SCS曲線法；蒸散發(fā)的模擬是在潛在蒸散發(fā)量的基礎(chǔ)上計算實際蒸散發(fā)量，包括植被冠層截留蒸發(fā)、植被蒸騰和土壤蒸發(fā)；土壤水分運動采用垂直分層水量平衡模擬，考慮下滲與壤中流過程；地下徑流的模擬分為淺層和深層地下徑流，涉及潛水再蒸發(fā)、灌溉等過程，基于水量平衡方程實現(xiàn)；河道匯流基于變動蓄量法或馬斯經(jīng)跟法對進(jìn)入河道中的水流進(jìn)行匯流演算。

1.2 基于EnKF法的狀態(tài)-參數(shù)優(yōu)化更新

1.2.1 基于EnKF法的狀態(tài)更新

EnKF法是Evensen[19]于1994年基于隨機(jī)動態(tài)預(yù)報理論，將集合預(yù)報思想引入KF法，提出的一種順序數(shù)據(jù)同化算法。EnKF法主要分為兩步：預(yù)報(狀態(tài)轉(zhuǎn)移過程)和分析(狀態(tài)更新過程)?；贛onte Carlo方法，將模型狀態(tài)變量加入初始狀態(tài)誤差(通常為高斯白噪聲)，生成一組(數(shù)目N)隨機(jī)變量集合Xk，i(k=0，i=1，2，…，N)，每組隨機(jī)變量Xk，i在k時刻的一步預(yù)報方程為

(1)

當(dāng)k+1時刻有觀測存在時，模型狀態(tài)變量的預(yù)報值基于觀測信息更新如下：

(2)

1.2.2 狀態(tài)擴(kuò)展、參數(shù)演變與過擬合問題處理

順序數(shù)據(jù)同化中參數(shù)估計可以通過狀態(tài)擴(kuò)展方式[20]實現(xiàn)，即將待估計的參數(shù)在其物理區(qū)間內(nèi)均勻采樣后，獲得待估計的參數(shù)集合，將其納入狀態(tài)向量X中，并基于觀測信息與更新狀態(tài)變量一樣對估計參數(shù)進(jìn)行更新[21]。假設(shè)待優(yōu)化估計的參數(shù)向量為δ，則擴(kuò)展后的狀態(tài)向量可以表示為[XTδT]T，基于EnKF法，k時刻狀態(tài)與參數(shù)的更新方程如下：

(3)

H*=[H0]

(4)

式中0為零矩陣，其維數(shù)為nob×npar，nob為觀測向量Yt中觀測變量的個數(shù)，npar為δ中參數(shù)的個數(shù)。

在以上基于狀態(tài)擴(kuò)展方式的EnKF法參數(shù)估計框架下，隨著EnKF法更新過程的繼續(xù)，參數(shù)集合將快速收斂，使得分析場很快趨近于背景場，最終導(dǎo)致觀測信息對參數(shù)更新不起作用，即出現(xiàn)濾波發(fā)散。為避免濾波發(fā)散問題，模型參數(shù)在同化過程中需要與狀態(tài)變量一樣具有一個隨時間的演變過程(狀態(tài)變量隨時間演變過程通過模型本身的物理機(jī)制實現(xiàn))。但是，這與一般意義上對模型參數(shù)的理解(隨機(jī)動態(tài)模型中的參數(shù)通常被假定為常數(shù))有沖突，因此，通常期望參數(shù)的變化相比于狀態(tài)變量的變化要小得多[1]。采用條件協(xié)方差法，模擬參數(shù)演變：

(5)

另外，參數(shù)更新過程中，一個重要的問題是參數(shù)的過擬合(過度更新)問題，其有可能導(dǎo)致模型模擬/預(yù)報過程中斷或模型無法正常運行。參數(shù)的過度更新，可能是由于數(shù)據(jù)同化過程中異常信號的加入(如強(qiáng)降雨)或?qū)δＰ湍M與觀測誤差的量化不合理，造成數(shù)據(jù)同化過程中對模型輸出結(jié)果與實測數(shù)據(jù)的過度擬合導(dǎo)致。為抑制參數(shù)同化過程中的過度更新，采用參數(shù)平滑法，引入平滑因子β∈(0，1.0)進(jìn)行處理[22]：

(6)

當(dāng)β= 0時，表示參數(shù)更新過程中沒有平滑作用，β越大，表示參數(shù)更新過程中的平滑作用越強(qiáng)。參考已有研究[22]，取β為0.8。

1.2.3 模型誤差與徑流觀測誤差的概化

SWAT模型產(chǎn)匯流模擬中，模型輸入的不確定性主要考慮降雨和最高/最低氣溫輸入誤差。假定站點降雨輸入誤差服從對數(shù)正態(tài)分布[23]，且該降雨輸入誤差在時間和空間上獨立：

(7)

(8)

式中：T為實測最高/最低氣溫；Te為誤差擾動后最高/最低氣溫集合。

模型結(jié)構(gòu)的不確定性對模型模擬預(yù)報的影響主要通過擾動模型狀態(tài)變量來刻畫，即對模型狀態(tài)變量的預(yù)報集合加上一定的高斯誤差擾動[15]：

(9)

模型參數(shù)的不確定性通過對敏感參數(shù)加以高斯誤差擾動：

(10)

本文研究冗余機(jī)械臂液壓驅(qū)動系統(tǒng)能量消耗優(yōu)化方法,給出機(jī)械臂液壓驅(qū)動簡圖模型,推導(dǎo)末端執(zhí)行器運動方程式.根據(jù)機(jī)械臂末端執(zhí)行器運動路徑建立動力學(xué)模型,并且對運動參數(shù)進(jìn)行離散化,構(gòu)造多目標(biāo)函數(shù),添加約束條件.簡化機(jī)械臂液壓驅(qū)動動力學(xué),引用動態(tài)規(guī)劃算法優(yōu)化液壓驅(qū)動控制系統(tǒng),在Matlab軟件中對動態(tài)規(guī)劃算法優(yōu)化后的液壓驅(qū)動系統(tǒng)進(jìn)行能量仿真,并且與無動態(tài)規(guī)劃的能量消耗仿真結(jié)果形成對比.結(jié)果表明:冗余機(jī)械臂液壓系統(tǒng)采用動態(tài)規(guī)劃算法優(yōu)化后,能量消耗較小,機(jī)械臂運動相對穩(wěn)定,有利于節(jié)約資源.

徑流觀測誤差主要源于水位觀測誤差和水位流量轉(zhuǎn)換誤差。假定徑流觀測誤差服從正態(tài)分布：

(11)

1.3 基于EnKF法的徑流數(shù)據(jù)同化方案

數(shù)據(jù)同化方案主要包括模型算子、觀測算子和同化算法。模型算子為SWAT模型，同化算法為EnKF法。觀測算子用于連接需要優(yōu)化的模型狀態(tài)變量/參數(shù)并輔助同化的觀測數(shù)據(jù)。輔助同化的觀測數(shù)據(jù)即為水文站點的實測徑流量，而需要優(yōu)化的模型狀態(tài)變量與參數(shù)的選擇根據(jù)其對降雨徑流過程模擬的重要性和敏感性，選擇地表徑流滯留量Qstor、土壤含水量SWly和河道蓄水量Vstored三類狀態(tài)變量和控制地表徑流、壤中流、地下徑流以及河道匯流過程的7類參數(shù)，見表1。當(dāng)狀態(tài)變量與觀測變量為同一類型的變量(如都為徑流)時，觀測算子為一個權(quán)重矩陣，其權(quán)重值可以通過狀態(tài)變量到觀測變量的插值過程獲得。然而，由于SWAT模型的高度非線性，無法用一個矩陣顯式地表示出SWAT模型中狀態(tài)變量X到與觀測變量Y所對應(yīng)的模型輸出變量的轉(zhuǎn)化過程，為此，需要把Y對應(yīng)的模型輸出變量(即觀測站點上的模擬徑流)納入X中。但在EnKF法更新過程中不對觀測站點的模擬徑流進(jìn)行更新處理，因為只有對狀態(tài)變量的更新才有意義[13]。此時，X中包括nst個模型狀態(tài)變量和被更新的參數(shù)與nob個觀測站點的模型模擬徑流，這樣X與Y中存在對應(yīng)變量，二者的狀態(tài)轉(zhuǎn)移矩陣H*可以寫成：

表1 SWAT模型徑流數(shù)據(jù)同化中被更新的模型狀態(tài)變量與參數(shù)

[0]nob×nstInob×nob

(12)

式中狀態(tài)轉(zhuǎn)移矩陣H*由nob×nst維的0矩陣和nob×nob維的單位矩陣構(gòu)成。

基于EnKF法和SWAT模型的站點觀測徑流數(shù)據(jù)同化方案的實現(xiàn)是將EnKF法耦合到SWAT模型模擬過程中，該同化方案的流程如圖1所示。

圖1 基于EnKF法與SWAT模型耦合的站點觀測徑流數(shù)據(jù)同化方案流程Fig.1 Flow chart of the streamflow data assimilation scheme based on the coupling of EnKF and SWAT

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

研究區(qū)為淮河上游淮濱水文站以上流域(簡稱淮濱流域)，流域控制面積為16 005 km2，整個流域除西部和西南部地區(qū)為山區(qū)外，大部分地區(qū)為平地，高程范圍25～1 117 m。土地利用類型主要為農(nóng)業(yè)用地(包括旱地34.2%和水田37.6%)和林地(23.6%)。土壤類型主要分為水稻土、黃褐土、黃棕壤、沙姜黑土等七類土壤。該流域冬春干旱少雨，夏秋悶熱多雨，多年平均降水量約900 mm，其中有50%～80%發(fā)生在6—9月，多年平均氣溫約15℃。

SWAT模型的輸入數(shù)據(jù)主要包括下墊面數(shù)據(jù)和氣象數(shù)據(jù)兩類(見表2)。降雨數(shù)據(jù)主要由流域內(nèi)共106個雨量站提供。最高/最低氣溫、太陽輻射、風(fēng)速和相對濕度數(shù)據(jù)從研究區(qū)內(nèi)(信陽)及其附近(固始、廣水、駐馬店和棗陽)共5個氣象站獲得，來自中國氣象科學(xué)數(shù)據(jù)共享網(wǎng)。另外，流域內(nèi)有6個水文站，分別為大坡嶺、長臺關(guān)、竹竿輔、息縣、潢川和淮濱，提供流域的逐日徑流量數(shù)據(jù)(因潢川站流量數(shù)據(jù)質(zhì)量有問題，故未采納)。研究區(qū)氣象、水文站點分布見圖2，圖中數(shù)字1～37為子流域代碼。

表2 淮濱流域 SWAT模型構(gòu)建中采用的數(shù)據(jù)

圖2 淮濱流域氣象、水文站站點分布及子流域劃分 Fig.2 Map of the meteorological and hydrological stations，and the subbasin division in Huaibin River Basin

3 徑流數(shù)據(jù)同化實現(xiàn)與結(jié)果分析

3.1 模型模擬與徑流觀測誤差量化

EnKF法可以顯式地處理模型與觀測誤差，并假定其都服從高斯分布。由EnKF法可以看出，對模型誤差(模型輸入、參數(shù)及結(jié)構(gòu))和徑流觀測誤差的合理量化在一定程度上決定了數(shù)據(jù)同化的效果。本次基于徑流數(shù)據(jù)同化的SWAT模型參數(shù)優(yōu)化研究中，主要考慮降雨和最高/最低氣溫輸入誤差、模型狀態(tài)誤差、模型參數(shù)誤差以及徑流觀測誤差，誤差均方差分別記為σp、σT、σs、σpar和σobs項。其中，模型狀態(tài)誤差又由模型狀態(tài)變量誤差和與觀測變量對應(yīng)的模型輸出變量誤差組成，分別記為σstate和σsflux。由于最高/最低氣溫輸入誤差對降雨徑流模擬的影響較為敏感，借鑒以往研究[24]，設(shè)σT為1℃。而其他因子誤差均方差的設(shè)置是在模型參數(shù)敏感性分析的基礎(chǔ)上，采用自動率定與人工調(diào)整相結(jié)合的方法對模型參數(shù)進(jìn)行調(diào)整后，基于調(diào)整后的模型采用最大后驗估計(MAP)法并結(jié)合全局優(yōu)化SCE-UA法[25]對其進(jìn)行量化的。

在MAP法中，需要事先給定各個因子誤差均方差的先驗分布。借鑒已有研究[8]，各因子誤差均方差的先驗分布設(shè)置如下：①σobs的先驗分布設(shè)置為均值9.25%、標(biāo)準(zhǔn)差0.925%的高斯分布；②其他因子σp、σstate、σsflux和σpar的先驗分布設(shè)置為給定閾值內(nèi)的均勻分布，其給定的閾值分別為0～0.5、0～0.3、0～0.5和0～0.3，另外σobs給定的閾值為0～0.5。確定各個因子誤差均方差的先驗分布后，即可獲得抽樣得到的某組誤差均方差的先驗概率?；谠摻M誤差均方差對模型降雨輸入、狀態(tài)變量、參數(shù)以及輸出變量進(jìn)行擾動，得到率定期內(nèi)的模型模擬徑流輸出系列集合。以流域出口控制站(淮濱站)的實測徑流為驗證數(shù)據(jù)，計算整個率定期內(nèi)該組誤差參數(shù)的似然函數(shù)。根據(jù)以上誤差均方差的先驗概率與似然函數(shù)計算得到其后驗概率。通過SCE-UA方法搜索使后驗概率達(dá)到最大時的誤差參數(shù)組，作為σobs、σp、σstate、σsflux和σpar的最優(yōu)估計，即最大后驗估計。

對以上誤差均方差的估計，采用了淮濱站2002—2004年共3年的日徑流量數(shù)據(jù)。在SCE-UA方法中，模型截斷運行次數(shù)為60 000次，而實際運行過程中，當(dāng)模型運行到34 860次時已收斂，最終獲得σobs、σp、σstate、σsflux和σpar的最優(yōu)估計結(jié)果，分別為0.087、0.305、0.028、0.29、0.113。

圖3為上述均方差參數(shù)條件下，2005—2008年的模型模擬徑流(95%的置信區(qū)間)與站點實測徑流過程，圖中黑色實線為實測徑流，灰色區(qū)域為模擬徑流95%的置信區(qū)間。從圖3可以看出，4年內(nèi)絕大多數(shù)觀測徑流數(shù)據(jù)都落在95%的置信區(qū)間內(nèi)，說明該誤差參數(shù)估計具有一定的合理性。同時也注意到，低流量時(Q<100 m3/s)，徑流模擬存在一定程度的低估，這與觀測誤差參數(shù)估算前模型的參數(shù)率定有一定關(guān)系。

圖3 最優(yōu)誤差均方差條件下2005—2008年的模型模擬(95%置信區(qū)間)與站點實測徑流過程Fig.3 Modeling(with 95% confidence interval) and observed runoff process based on the optimal error parameters over 2005—2008

3.2 模型參數(shù)同化更新過程

根據(jù)流域內(nèi)5個站點(大坡嶺、長臺關(guān)、竹竿輔、息縣和淮濱站)2003—2005年的日徑流實測數(shù)據(jù)，基于EnKF法實時更新SWAT模型敏感參數(shù)與狀態(tài)變量，以期獲得模型參數(shù)的最優(yōu)估計。從被同化更新的7類共59個敏感參數(shù)中(表1)，隨機(jī)選取2個代表性參數(shù)，即水文條件Ⅱ下SCS曲線數(shù)和表征土壤有效持水能力參數(shù)(SOL_AWC)，給出其在數(shù)據(jù)同化過程中的演變過程，見圖4。從圖4可以看出，隨著觀測信息的不斷加入，這兩個參數(shù)的集合分布越來越窄，集合均值逐漸趨于穩(wěn)定。從參數(shù)集合在某一天的概率密度分布(圖4小圖)可以看出，參數(shù)集合由初始時的均勻分布，通過實時的同化更新，逐漸演變?yōu)楦咚狗植?；同時，集合的方差/標(biāo)準(zhǔn)差逐漸變小，說明參數(shù)估計的不確定性在逐漸減小。

圖4 基于EnKF法徑流數(shù)據(jù)同化的SWAT模型參數(shù)CN2和SOL_AWC的演變過程Fig.4 Evolution process of CN2 and SOL_AWC in EnKF based streamflow assimilation

圖5為SWAT模型中被優(yōu)化的七類共59個敏感參數(shù)其集合均方差隨徑流數(shù)據(jù)同化更新的變化過程。從圖5可以看出，七類參數(shù)的集合均方差隨同化更新過程的繼續(xù)逐漸減小，表明參數(shù)集合的分布越來越窄，即參數(shù)估計的不確定性逐漸減小。最終，參數(shù)集合的均方差都收斂到初始集合均方差的5%以內(nèi)，說明基于EnKF法的徑流數(shù)據(jù)同化對SWAT模型參數(shù)具有一定的優(yōu)化更新能力。

圖5 SWAT模型中七類參數(shù)的集合均方差隨徑流數(shù)據(jù)同化的變化過程 Fig.5 Variations of the ensemble mean square errors for 7 types of parameters in SWAT updated with the streamflow assimilation

3.3 模型參數(shù)優(yōu)化效果評估

基于實測徑流數(shù)據(jù)同化進(jìn)行SWAT模型參數(shù)優(yōu)化時，由于真實(或最優(yōu))的模型參數(shù)值未知且很多參數(shù)無法直接量測，很難直接對參數(shù)估計效果進(jìn)行評估。本文通過基于優(yōu)化后的參數(shù)模擬獲得的流域內(nèi)站點徑流過程間接地對參數(shù)優(yōu)化效果進(jìn)行評估。具體方案為：基于對淮濱流域內(nèi)5個站點2003—2005年(率定期)實測徑流數(shù)據(jù)同化獲得的參數(shù)估計結(jié)果運行模型，得到2006—2008年(驗證期)的徑流模擬過程(記為EnKF)；將該徑流模擬過程與站點實測徑流過程進(jìn)行比較，通過模擬徑流與站點實測徑流過程線的接近/貼合程度以及相應(yīng)的評估指標(biāo)，對基于EnKF法徑流數(shù)據(jù)同化獲得的模型參數(shù)優(yōu)化效果進(jìn)行評估。同時，將未經(jīng)率定的參數(shù)(即在給定模型輸入條件下，SWAT模型默認(rèn)的參數(shù))在驗證期的徑流模擬過程(記為SIM)作為參照。

表3為基于EnKF法優(yōu)化的參數(shù)及未經(jīng)率定的參數(shù)獲得的流域內(nèi)5個水文站點在驗證期的徑流模擬結(jié)果。從表3可以看出，基于EnKF法獲得的參數(shù)估計結(jié)果所對應(yīng)的流域內(nèi)5個站點在驗證期的徑流模擬精度相較于未經(jīng)率定的情況有大幅度的提高。EnKF法對應(yīng)5個站點徑流模擬的相對偏差(除竹竿輔站外)都保持在15%以內(nèi)，均方根誤差顯著減小，尤其是流域出口淮濱站的均方根誤差減小了34.5%，而納什效率系數(shù)和相關(guān)系數(shù)r均明顯增加，5個站點的納什效率系數(shù)都在0.6以上，其中大坡嶺、息縣、淮濱3站在0.8以上，尤其流域出口淮濱站達(dá)到了0.88，表明基于EnKF法在SWAT模型中同化站點徑流數(shù)據(jù)、更新模型參數(shù)與狀態(tài)變量可以得到較好的模型參數(shù)估計結(jié)果。

表3 基于EnKF法優(yōu)化及未經(jīng)率定的參數(shù)條件下流域內(nèi)5個站點2006—2008年的徑流模擬結(jié)果

圖6為基于EnKF法優(yōu)化的參數(shù)(EnKF)及未經(jīng)率定參數(shù)(SIM)條件下，流域出口淮濱站在驗證期(2006-01-01—2008-12-31)的徑流模擬過程。從圖6可以看出，基于EnKF法參數(shù)估計結(jié)果的驗證期徑流(EnKF)，其洪峰過程(圖6(a))和基流過程(圖6(b)，圖6(b)縱坐標(biāo)為以10為底的對數(shù)坐標(biāo))都明顯優(yōu)于未經(jīng)參數(shù)率定條件下的徑流模擬結(jié)果(SIM)。雖然基于EnKF法估計的參數(shù)所對應(yīng)的徑流過程相較于實測徑流過程還存在一定程度的低估，但其相對于未率定過程而言其低估程度大為減小，表明通過EnKF法同化站點實測徑流對SWAT模型參數(shù)具有一定的優(yōu)化能力，可以作為未來SWAT模型參數(shù)率定的手段。

圖6 基于EnKF法及未經(jīng)率定的參數(shù)條件下，流域出口淮濱站在驗證期的徑流過程Fig.6 Runoff process at the outlet of Huaibin River Basin based on the EnKF optimized parameters and uncalibrated model parameters

4 結(jié) 語

本文基于狀態(tài)擴(kuò)展方式的EnKF法，在對SWAT模型誤差(模型輸入、模型結(jié)構(gòu)、模型參數(shù))和站點徑流觀測誤差采用最大后驗估計方法合理量化的基礎(chǔ)上，以條件協(xié)方差法并引入平滑因子有效刻畫參數(shù)隨時間的演變過程，構(gòu)建了基于徑流數(shù)據(jù)同化的水文模型參數(shù)優(yōu)化方案。通過站點實測徑流數(shù)據(jù)的同化對SWAT模型參數(shù)及狀態(tài)變量進(jìn)行實時更新，以淮河上游淮濱水文站以上流域為研究區(qū)，對該方案中模型參數(shù)的優(yōu)化效果進(jìn)行了評估。

結(jié)果顯示，在站點實測徑流數(shù)據(jù)同化過程中，被優(yōu)化更新的參數(shù)集合逐漸收斂并趨于穩(wěn)定。通過對比分析基于穩(wěn)定后的參數(shù)獲得的模擬徑流和實測徑流過程，發(fā)現(xiàn)通過徑流數(shù)據(jù)同化獲得的參數(shù)具有較好的模擬預(yù)測精度，基于數(shù)據(jù)同化方法得到的參數(shù)進(jìn)行降雨徑流模擬，流域出口淮濱站徑流模擬的納什效率系數(shù)可達(dá)0.88，間接證明了基于EnKF法的徑流數(shù)據(jù)同化對SWAT模型參數(shù)具有一定的優(yōu)化估計能力，且采用數(shù)據(jù)同化方式進(jìn)行水文模型參數(shù)的率定具有一定的可行性。

相較于目前常用的參數(shù)率定方法(如全局優(yōu)化方法)將模型模擬誤差全部歸咎于模型參數(shù)誤差這一不足[26-27]，通過同化站點徑流進(jìn)行模型參數(shù)優(yōu)化，可綜合考慮模型輸入、模型結(jié)構(gòu)和模型參數(shù)的不確定性，理論上，基于數(shù)據(jù)同化進(jìn)行模型參數(shù)的優(yōu)化可獲得更優(yōu)的參數(shù)估計。從研究實例來看，當(dāng)同化進(jìn)行到600～800時間步長后，模型參數(shù)就可收斂且趨于穩(wěn)定，說明相較于一般的參數(shù)率定方法(一般需要一個豐、平、枯水文周期的實測數(shù)據(jù))，基于數(shù)據(jù)同化的參數(shù)優(yōu)化要求更短的實測數(shù)據(jù)系列，這對于資料匱乏地區(qū)(通常數(shù)據(jù)系列較短)降雨徑流模型中參數(shù)的優(yōu)化具有重要意義。然而，基于數(shù)據(jù)同化的模型參數(shù)優(yōu)化必須建立在對模型模擬誤差與徑流觀測誤差合理量化的基礎(chǔ)上，而模型模擬誤差與徑流觀測誤差的估計一直以來都是水文界的一個難題[8，28]，這也是當(dāng)前該方法在應(yīng)用中面臨的一大挑戰(zhàn)。另外，由于水文模型，尤其具有物理機(jī)制的分布式模型的變量與參數(shù)眾多[29-30]，導(dǎo)致徑流數(shù)據(jù)同化過程中不可避免地存在高維度現(xiàn)象，這樣同化效果的不確定性及同化實施的難度都將大大增加，因此，基于徑流數(shù)據(jù)同化進(jìn)行水文模型參數(shù)的優(yōu)化還需要發(fā)展新的、更優(yōu)的數(shù)據(jù)同化方法與策略。