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

王復生，田 娟

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

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

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

本文以巢湖流域桃溪水文站控制子流域為研究對象，建立MIKE NAM模型，綜合考慮NAM模型中各類參數的取值情況，運用Sobol方法對模型參數的敏感性進行定量分析，評估各參數對模型輸出結果方差的影響程度，將高敏感性參數作為主要率定參數，以期提高模型率定效率和模擬精度。

1 模型與方法

1.1 MIKE NAM模型

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

1.2 Sobol參數敏感性分析方法

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

Sobol方法假設模型結構可以表示為函數形式u=f(x)，其中模型參數x=x1，x2，…，xn可視為n維離散點，u為標量輸出。

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

(1)

式中，k=i1，…，is。

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

(2)

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

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

(3)

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

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

(4)

則模型的敏感性指數可以表示為：

式中，Si—參數xi單獨作用時的敏感度；Sij—參數xi和xj相互作用的敏感度；STi—參數xi和其它參數共同作用的敏感度；D～i—除參數xi外其它參數共同作用產生的方差。

2 研究區域與模型構建

2.1 研究區域概況

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

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

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

根據研究區域內實測降雨資料及雨量站的分布情況，選取9個代表站點進行降雨徑流計算，利用泰森多邊形法分析得到各個雨量站的計算權重值見表1。

表1 桃溪站控制流域面雨量計算各雨量站權重

2.2 NAM模型相關參數

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

表2 NAM模型主要參數

3 計算結果與分析

3.1 參數敏感性分析

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

圖2 參數敏感性分析成果圖

由敏感性分析成果可知：

以洪峰流量作為目標函數時，坡面流匯流系數CQOF和坡面流匯流時間常數CK1、2的一階敏感性指數、全局敏感性指數遠高于其它參數，說明這兩個參數對洪峰流量大小起決定性作用。CQOF的一階敏感性指數大于0，CK1、2的一階敏感性指數小于0，說明坡面流匯流系數與流量呈正相關，坡面流匯流時間常數與流量呈負相關，這與參數代表的物理意義相符。

以峰現時間作為目標函數時，坡面流匯流時間常數CK1、2的敏感性指數高于其它參數，說明該參數對洪峰出現的時間有較大影響。

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

3.2 模型率定

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

表3 桃溪站參數率定成果

圖3 桃溪站模擬洪水與實測洪水擬合效果

圖4 桃溪站累計洪量擬合效果

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

表4 計算誤差統計表

4 結論

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

(2)對CQOF和CK1、2兩個參數采用人工率定結果，其它參數采用模型自動率定結果，以豐樂河桃溪水文站控制流域為研究區，模擬了2020年汛期降雨作用下的區域降雨徑流過程，并根據實測流量數據對模擬效果進行擬合檢驗，結果表明模型計算誤差在合理范圍內。

(3)本文構建的NAM模型中未考慮流域內圩區的調蓄以及流域下游水位實時變化的作用，對模擬精度有一定的影響，并且參數敏感性分析成果缺少不同區域不同場次降雨的驗證。后續可以將NAM模型與MIKE11水動力模型進行耦合，針對不同降雨強度、不同流域特征進行建模和率定，進而實現巢湖全流域的降雨徑流模擬。本研究可為巢湖流域水文水動力模型的構建和驗證提供技術支撐和應用基礎。