基于Sobol方法的MIKE NAM模型參數(shù)敏感性分析

王復(fù)生，田 娟

(安徽省水利水電勘測設(shè)計(jì)研究總院有限公司,安徽 合肥 230088)

降雨徑流模擬是防洪排澇計(jì)算、流域洪水演算以及水資源保護(hù)規(guī)劃的重要研究方法。隨著數(shù)字化技術(shù)的發(fā)展，基于模型對(duì)流域內(nèi)的降雨徑流過程進(jìn)行模擬計(jì)算已經(jīng)成為了一種主要趨勢。具有物理機(jī)制的降雨徑流水文模型可以綜合考慮降雨、土壤特性、下墊面條件等多方面因素，反演不同邊界條件下的降雨產(chǎn)匯流規(guī)律及空間分布特征[1- 2]，對(duì)洪水演進(jìn)理論發(fā)展和工程應(yīng)用均具有重要價(jià)值。

MIKE NAM模型是一種集總式水文模型，可以用于一個(gè)流域或由多個(gè)子分區(qū)組成的較大河流及由河流、渠道構(gòu)成的復(fù)雜河系的降雨徑流模擬過程，在國內(nèi)外得到了廣泛的研究與應(yīng)用，并取得了良好的模擬效果[3- 5]。NAM模型的參數(shù)較多，模擬過程中存在模型參數(shù)的不確定性與敏感性，參數(shù)和邊界條件之間的相互影響等問題[6- 7]。參數(shù)取值的合理性以及對(duì)參數(shù)含義的準(zhǔn)確理解對(duì)模型精度會(huì)產(chǎn)生很大影響，因此對(duì)模型參數(shù)進(jìn)行全局敏感性分析，篩選出高敏感性參數(shù)進(jìn)行重點(diǎn)率定十分有必要。

本文以巢湖流域桃溪水文站控制子流域?yàn)檠芯繉?duì)象，建立MIKE NAM模型，綜合考慮NAM模型中各類參數(shù)的取值情況，運(yùn)用Sobol方法對(duì)模型參數(shù)的敏感性進(jìn)行定量分析，評(píng)估各參數(shù)對(duì)模型輸出結(jié)果方差的影響程度，將高敏感性參數(shù)作為主要率定參數(shù)，以期提高模型率定效率和模擬精度。

1 模型與方法

1.1 MIKE NAM模型

MIKE是由丹麥水資源與水環(huán)境研究所(DHI)開發(fā)的專業(yè)工程軟件，適用于河流、灌溉渠道及其它水體的一二維水動(dòng)力、水質(zhì)、泥沙運(yùn)輸模擬等[8- 9]。MIKE NAM模型作為一種新興的水文模擬工具，在流域降雨產(chǎn)匯流模擬中具有良好的適用性，被廣泛應(yīng)用到世界各地不同氣象水文條件的流域，在國內(nèi)長江、黃河、松花江流域均有應(yīng)用案例，是一個(gè)經(jīng)過大量工程實(shí)踐驗(yàn)證的模型工具。NAM模型將土壤含水層分成積雪儲(chǔ)水層(Snow Storage)、地表儲(chǔ)水層(Surface Storage)、淺層儲(chǔ)水層(Lower Zone Storage)和地下水儲(chǔ)水層(Ground Water Storage)4個(gè)部分，模型分別進(jìn)行連續(xù)計(jì)算以模擬流域中各種相應(yīng)的水文過程。

1.2 Sobol參數(shù)敏感性分析方法

Sobol方法是一種基于方差分解的定量全局敏感性分析算法，相比于定性敏感性分析方法，該方法能夠通過敏感性量化指標(biāo)的計(jì)算直接給出模型參數(shù)的敏感性大小，其核心是將目標(biāo)函數(shù)的總方差分解為單個(gè)參數(shù)的方差和多個(gè)參數(shù)間相互作用的方差，目前已被廣泛應(yīng)用于經(jīng)濟(jì)、環(huán)境等多個(gè)領(lǐng)域[10]。

Sobol方法假設(shè)模型結(jié)構(gòu)可以表示為函數(shù)形式u=f(x)，其中模型參數(shù)x=x1，x2，…，xn可視為n維離散點(diǎn)，u為標(biāo)量輸出。

假設(shè)函數(shù)f(x)可積，xi在[0，1]上服從均勻分布，且f(x)滿足：

(1)

式中，k=i1，…，is。

則函數(shù)f(x)的方差分解表達(dá)式可以表示為：

(2)

式中，1≤i1<…is≤n(1≤s≤n)，則被加數(shù)共有2n項(xiàng)。

類似的，模型的總方差也可以分解為單個(gè)參數(shù)和多個(gè)參數(shù)相互作用的組合：

(3)

式中，D—模型的總方差；Di—參數(shù)xi產(chǎn)生的方差；Dij—參數(shù)xi和xj相互作用產(chǎn)生的方差，D1，2，…，n—n個(gè)參數(shù)共同作用產(chǎn)生的方差。

將上式歸一化即可得到各參數(shù)和參數(shù)之間相互作用的敏感性：

(4)

則模型的敏感性指數(shù)可以表示為：

式中，Si—參數(shù)xi單獨(dú)作用時(shí)的敏感度；Sij—參數(shù)xi和xj相互作用的敏感度；STi—參數(shù)xi和其它參數(shù)共同作用的敏感度；D～i—除參數(shù)xi外其它參數(shù)共同作用產(chǎn)生的方差。

2 研究區(qū)域與模型構(gòu)建

2.1 研究區(qū)域概況

本研究對(duì)巢湖流域上游豐樂河桃溪水文站控制流域進(jìn)行建模分析。

豐樂河發(fā)源于大別山余脈，流經(jīng)雙河、桃溪、豐樂至三河鎮(zhèn)下大潭灣與杭埠河匯合后流入巢湖，河道全長117.5km，流域面積2124km2，是杭埠河最大的支流。豐樂河干流上的桃溪水文站設(shè)立于1951年7月，控制集水面積1510km2，約占豐樂河流域總面積的75%，有長系列的實(shí)測流量資料，且精度可靠，是巢湖流域水文分析計(jì)算重要的參證站點(diǎn)。豐樂河桃溪站控制流域地形圖如圖1所示。

圖1 桃溪站控制流域地形圖

根據(jù)研究區(qū)域內(nèi)實(shí)測降雨資料及雨量站的分布情況，選取9個(gè)代表站點(diǎn)進(jìn)行降雨徑流計(jì)算，利用泰森多邊形法分析得到各個(gè)雨量站的計(jì)算權(quán)重值見表1。

表1 桃溪站控制流域面雨量計(jì)算各雨量站權(quán)重

2.2 NAM模型相關(guān)參數(shù)

NAM模型所需數(shù)據(jù)資料主要為流域降雨及蒸發(fā)數(shù)據(jù)，模型的初始條件包括開始時(shí)刻流域地表儲(chǔ)水層和根區(qū)儲(chǔ)水層的土壤含水率，以及坡面流、壤中流和基流的初始值[11]。NAM模型的輸出結(jié)果主要包括地表坡面徑流、壤中流和基流過程等。模型主要參數(shù)及取值范圍詳見表2。作為集總式水文模型，NAM模型中每個(gè)子流域的參數(shù)均為本區(qū)域的參數(shù)平均值。

表2 NAM模型主要參數(shù)

3 計(jì)算結(jié)果與分析

3.1 參數(shù)敏感性分析

由于Umax與Lmax相關(guān)性較高，可視為一個(gè)參數(shù)，本文以Umax、CQOF、CKIF、CK1、2、TOF、TIF、TG、CKBF共8個(gè)模型參數(shù)作為自變量函數(shù)，用蒙特卡羅法(Monte Carlo Simulation)對(duì)自變量進(jìn)行抽樣，每個(gè)參數(shù)得到10組樣本值，根據(jù)Sobol參數(shù)敏感性分析方法產(chǎn)生10×(8+2)=100個(gè)參數(shù)組合，將這100組參數(shù)分別帶入桃溪站控制流域NAM模型中，對(duì)2020年6月8日—8月8日的區(qū)域降雨進(jìn)行模擬。分別以洪峰流量、峰現(xiàn)時(shí)間和最大7d洪量作為NAM模型輸出結(jié)果的目標(biāo)函數(shù)，用一階敏感性指數(shù)Si和全局敏感性指數(shù)STi來量化輸出結(jié)果對(duì)各個(gè)參數(shù)的敏感度高低。模型參數(shù)敏感性分析成果如圖2所示。

圖2 參數(shù)敏感性分析成果圖

由敏感性分析成果可知：

以洪峰流量作為目標(biāo)函數(shù)時(shí)，坡面流匯流系數(shù)CQOF和坡面流匯流時(shí)間常數(shù)CK1、2的一階敏感性指數(shù)、全局敏感性指數(shù)遠(yuǎn)高于其它參數(shù)，說明這兩個(gè)參數(shù)對(duì)洪峰流量大小起決定性作用。CQOF的一階敏感性指數(shù)大于0，CK1、2的一階敏感性指數(shù)小于0，說明坡面流匯流系數(shù)與流量呈正相關(guān)，坡面流匯流時(shí)間常數(shù)與流量呈負(fù)相關(guān)，這與參數(shù)代表的物理意義相符。

以峰現(xiàn)時(shí)間作為目標(biāo)函數(shù)時(shí)，坡面流匯流時(shí)間常數(shù)CK1、2的敏感性指數(shù)高于其它參數(shù)，說明該參數(shù)對(duì)洪峰出現(xiàn)的時(shí)間有較大影響。

以最大7d洪量作為目標(biāo)函數(shù)時(shí)，坡面流匯流系數(shù)CQOF的敏感性指數(shù)高于其它參數(shù)，說明該參數(shù)對(duì)洪量有較大影響。

3.2 模型率定

根據(jù)敏感性分析結(jié)論，NAM模型模擬結(jié)果對(duì)不同參數(shù)變化的敏感性是有差異的，參數(shù)CQOF和CK1、2對(duì)NAM模型精度的影響比其它參數(shù)大，率定時(shí)需重點(diǎn)考慮。因此本研究在對(duì)2020年6月8日—8月8日降雨進(jìn)行模擬時(shí)，先采用NAM模型自帶的SCE全局尋優(yōu)算法對(duì)8個(gè)參數(shù)進(jìn)行自動(dòng)率定，迭代次數(shù)設(shè)置為1000次，在此基礎(chǔ)上對(duì)CQOF和CK1、2兩個(gè)重點(diǎn)參數(shù)進(jìn)行人工率定，最終得到的參數(shù)率定成果如表3所示。桃溪站模擬洪水與實(shí)測洪水?dāng)M合效果見圖3，洪量的擬合效果見圖4。

表3 桃溪站參數(shù)率定成果

圖3 桃溪站模擬洪水與實(shí)測洪水?dāng)M合效果

圖4 桃溪站累計(jì)洪量擬合效果

2020年桃溪站實(shí)測洪峰流量為1239m3/s，模型計(jì)算得到的洪峰流量為1304m3/s，差值為5.2%，原因在于2020年桃溪站上游柏林圩等圩口發(fā)生了潰破，對(duì)豐樂河洪峰流量起到了削峰作用。2020年桃溪站實(shí)測最大7d洪量為4.26億m3/s，模型計(jì)算最大7d洪量為4.41億m3，誤差為3.5%，桃溪站實(shí)測洪峰與模型模擬洪峰的出現(xiàn)時(shí)間相差10h，誤差在合理范圍內(nèi)?？紤]破圩還原后，誤差進(jìn)一步減小，計(jì)算誤差結(jié)果見表4。

表4 計(jì)算誤差統(tǒng)計(jì)表

4 結(jié)論

(1)基于Sobol方法對(duì)MIKE NAM模型中的主要參數(shù)進(jìn)行全局敏感性分析，分析結(jié)果表明：參數(shù)CQOF和CK1、2為高敏感度參數(shù)，對(duì)模型精度影響較大，建議作為重點(diǎn)率定參數(shù)。

(2)對(duì)CQOF和CK1、2兩個(gè)參數(shù)采用人工率定結(jié)果，其它參數(shù)采用模型自動(dòng)率定結(jié)果，以豐樂河桃溪水文站控制流域?yàn)檠芯繀^(qū)，模擬了2020年汛期降雨作用下的區(qū)域降雨徑流過程，并根據(jù)實(shí)測流量數(shù)據(jù)對(duì)模擬效果進(jìn)行擬合檢驗(yàn)，結(jié)果表明模型計(jì)算誤差在合理范圍內(nèi)。

(3)本文構(gòu)建的NAM模型中未考慮流域內(nèi)圩區(qū)的調(diào)蓄以及流域下游水位實(shí)時(shí)變化的作用，對(duì)模擬精度有一定的影響，并且參數(shù)敏感性分析成果缺少不同區(qū)域不同場次降雨的驗(yàn)證。后續(xù)可以將NAM模型與MIKE11水動(dòng)力模型進(jìn)行耦合，針對(duì)不同降雨強(qiáng)度、不同流域特征進(jìn)行建模和率定，進(jìn)而實(shí)現(xiàn)巢湖全流域的降雨徑流模擬。本研究可為巢湖流域水文水動(dòng)力模型的構(gòu)建和驗(yàn)證提供技術(shù)支撐和應(yīng)用基礎(chǔ)。