基于EnKF法的徑流數(shù)據(jù)同化對(duì)SWAT模型參數(shù)優(yōu)化效果評(píng)估

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

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

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

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

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

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

1.1 SWAT模型

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

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

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

EnKF法是Evensen[19]于1994年基于隨機(jī)動(dòng)態(tài)預(yù)報(bào)理論，將集合預(yù)報(bào)思想引入KF法，提出的一種順序數(shù)據(jù)同化算法。EnKF法主要分為兩步：預(yù)報(bào)(狀態(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時(shí)刻的一步預(yù)報(bào)方程為

(1)

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

(2)

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

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

(3)

H*=[H0]

(4)

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

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

(5)

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

(6)

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

1.2.3 模型誤差與徑流觀測(cè)誤差的概化

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

(7)

(8)

式中：T為實(shí)測(cè)最高/最低氣溫；Te為誤差擾動(dòng)后最高/最低氣溫集合。

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

(9)

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

(10)

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

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

(11)

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

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

圖1 基于EnKF法與SWAT模型耦合的站點(diǎn)觀測(cè)徑流數(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ū)為淮河上游淮濱水文站以上流域(簡(jiǎn)稱淮濱流域)，流域控制面積為16 005 km2，整個(gè)流域除西部和西南部地區(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個(gè)雨量站提供。最高/最低氣溫、太陽輻射、風(fēng)速和相對(duì)濕度數(shù)據(jù)從研究區(qū)內(nèi)(信陽)及其附近(固始、廣水、駐馬店和棗陽)共5個(gè)氣象站獲得，來自中國(guó)氣象科學(xué)數(shù)據(jù)共享網(wǎng)。另外，流域內(nèi)有6個(gè)水文站，分別為大坡嶺、長(zhǎng)臺(tái)關(guān)、竹竿輔、息縣、潢川和淮濱，提供流域的逐日徑流量數(shù)據(jù)(因潢川站流量數(shù)據(jù)質(zhì)量有問題，故未采納)。研究區(qū)氣象、水文站點(diǎn)分布見圖2，圖中數(shù)字1～37為子流域代碼。

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

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

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

3.1 模型模擬與徑流觀測(cè)誤差量化

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

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

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

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

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

圖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個(gè)敏感參數(shù)其集合均方差隨徑流數(shù)據(jù)同化更新的變化過程。從圖5可以看出，七類參數(shù)的集合均方差隨同化更新過程的繼續(xù)逐漸減小，表明參數(shù)集合的分布越來越窄，即參數(shù)估計(jì)的不確定性逐漸減小。最終，參數(shù)集合的均方差都收斂到初始集合均方差的5%以內(nèi)，說明基于EnKF法的徑流數(shù)據(jù)同化對(duì)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)化效果評(píng)估

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

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

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

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

圖6 基于EnKF法及未經(jīng)率定的參數(shù)條件下，流域出口淮濱站在驗(yàn)證期的徑流過程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法，在對(duì)SWAT模型誤差(模型輸入、模型結(jié)構(gòu)、模型參數(shù))和站點(diǎn)徑流觀測(cè)誤差采用最大后驗(yàn)估計(jì)方法合理量化的基礎(chǔ)上，以條件協(xié)方差法并引入平滑因子有效刻畫參數(shù)隨時(shí)間的演變過程，構(gòu)建了基于徑流數(shù)據(jù)同化的水文模型參數(shù)優(yōu)化方案。通過站點(diǎn)實(shí)測(cè)徑流數(shù)據(jù)的同化對(duì)SWAT模型參數(shù)及狀態(tài)變量進(jìn)行實(shí)時(shí)更新，以淮河上游淮濱水文站以上流域?yàn)檠芯繀^(qū)，對(duì)該方案中模型參數(shù)的優(yōu)化效果進(jìn)行了評(píng)估。

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

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