白河上游流域洪水模擬及參數(shù)敏感性分析

胡昊,李磊,王文川

(1.黃河水利職業(yè)技術(shù)學(xué)院,河南 開封 475004;2.華北水利水電大學(xué) 水資源學(xué)院,河南 鄭州 450046;3.中水珠江規(guī)劃勘測設(shè)計有限公司,廣州 510610)

近年來,受氣候變化及人類活動的影響,山洪災(zāi)害等自然災(zāi)害頻發(fā)[1]。這不僅對人民的生命財產(chǎn)安全帶來威脅,還影響社會和諧發(fā)展進(jìn)程。如何準(zhǔn)確預(yù)報山洪災(zāi)害,是水文專家及學(xué)者迫切要解決的難題。應(yīng)用水文模型模擬洪水,是一種切實有效的非工程防災(zāi)手段。水文模型是通過數(shù)學(xué)物理公式來描述水文過程的高維復(fù)雜耦合函數(shù),因此，在水文模型的構(gòu)建中會存在很大的不確定性,而這種不確定性會影響模型的模擬精度和結(jié)果的準(zhǔn)確度[2-3]。這種不確定性有一部分源自模型參數(shù)的不確定性,由于水文模型的參數(shù)較多,而不同的水文模型的參數(shù)不同,并且不是所有參數(shù)都具有實際物理含義。因此,很難做到同時提高所有參數(shù)的精度,且模型中一般也只有少數(shù)參數(shù)會對結(jié)果起到?jīng)Q定性作用。對于某一個特定的流域,對模型輸出有重要影響以及有率定需要的參數(shù)是有限的[4]。所以要對參數(shù)的影響進(jìn)行定量評估,并對參數(shù)進(jìn)行敏感性分析,這樣才能為模型實現(xiàn)更加高效的優(yōu)化與率定提供支撐。

為了更好地應(yīng)對山洪災(zāi)害,中國水科院提出了時空變源混合產(chǎn)流模型[5]。該模型本質(zhì)上是一種能夠在超滲、蓄滿兩種產(chǎn)流機(jī)制之間發(fā)生時空轉(zhuǎn)變的水文模型。本文以白河上游流域為研究對象,利用時空變源混合產(chǎn)流模型對流域進(jìn)行洪水模擬,運(yùn)用GLUE法和Sobol法對參數(shù)敏感性進(jìn)行分析,并判斷流域的產(chǎn)流方式。此研究對于類似白河上游山丘區(qū)等中小流域開展時空變源混合產(chǎn)流模型參數(shù)的敏感性分析及參數(shù)選取具有一定參考價值。

1 方法介紹

1.1 時空變源混合產(chǎn)流模型

時空變源混合產(chǎn)流模型是由劉昌軍等[6]提出的一種可為山洪災(zāi)害預(yù)報和預(yù)警工作提供便利的水文模型,適用于短歷時、強(qiáng)降雨的中小流域。該模型主要包括兩個方面：時空變源和混合產(chǎn)流。在外因(降雨入滲和蒸發(fā))和內(nèi)因(重力和基質(zhì)吸力)共同作用下土壤含水量的時空變化過程稱為時空變源;在土壤含水量時空變化下流域產(chǎn)流呈現(xiàn)出超滲/蓄滿時空動態(tài)的組合過程稱為混合產(chǎn)流。時空變源混合產(chǎn)流模型根據(jù)每一個計算時段內(nèi)的水文響應(yīng)單元的含水量以及累積下滲量,計算出超滲產(chǎn)流以及蓄滿產(chǎn)流的面積變化,同時判別出雨強(qiáng)和下墊面滲流能力之間的關(guān)系,實現(xiàn)超滲/蓄滿產(chǎn)流在每一個地貌水文響應(yīng)單元上的時空轉(zhuǎn)化。以子流域為計算單元的時空變源混合產(chǎn)流模型,其混合產(chǎn)流理論包含的產(chǎn)流過程有:流域蒸發(fā)、土壤蒸發(fā)、截流填洼、下滲、超滲、蓄滿、優(yōu)先流、壤中流、地下徑流。其中，優(yōu)先流是水在土壤中的某一部分的運(yùn)動要快于其他部分的產(chǎn)流，產(chǎn)生這部分徑流的區(qū)域稱為優(yōu)先流區(qū)域。混合產(chǎn)流成分由超滲地表產(chǎn)流、蓄滿產(chǎn)流、優(yōu)先流、壤中流和地下徑流組成。蓄滿產(chǎn)流由上層土壤優(yōu)先蓄滿得到的蓄滿產(chǎn)流和深層土壤蓄滿后向上層土壤補(bǔ)給產(chǎn)生的蓄滿徑流兩部分組成;超滲地表產(chǎn)流包括不透水面積的直接超滲地表產(chǎn)流和透水面積在雨強(qiáng)大于下滲能力時的超滲產(chǎn)流;優(yōu)先流為上層土壤出流;壤中流為下層土壤出流;地下徑流為地下水出流。時空變源混合產(chǎn)流模型產(chǎn)流過程如圖1所示[7]。

圖1 時空變源混合產(chǎn)流過程概化示意圖

時空變源混合產(chǎn)流模型是一種基于物理機(jī)制建立的水文模型,可分為兩種計算模型:日模型與次洪模型。模型是以水文循環(huán)模塊化為基礎(chǔ)搭建的,主要有降雨以及徑流兩種時間序列值。模型及其相關(guān)計算的主要參數(shù)及其取值范圍見表1[8]。

表1 時空變源混合產(chǎn)流模型主要產(chǎn)流參數(shù)及取值范圍

1.2 GLUE法

普適似然函數(shù)不確定性評估(Generalized Likelihood Uncertainty Estimation,GLUE)方法[9-10],起源于區(qū)域敏感性分析(Regionalized Sensitivity Analysis,RSA)方法。水文模型的模擬性能通過多個參數(shù)的互相作用來體現(xiàn),單個參數(shù)并不能完全反映水文模型的預(yù)報情況。通常,根據(jù)水文模型參數(shù)的物理意義或依據(jù)經(jīng)驗設(shè)定參數(shù)的取值范圍,然后使用蒙特卡洛隨機(jī)抽取參數(shù)組合樣本,并將隨機(jī)模擬參數(shù)組代入模型進(jìn)行運(yùn)算,根據(jù)模型模擬的洪水過程與實際洪水過程之間的似然函數(shù)值,來推求出各參數(shù)組合樣本的似然值。參數(shù)敏感性是根據(jù)參數(shù)與似然值之間的分布確定的,具體步驟如下:

步驟1確定GLUE似然判據(jù)。選擇反映實測洪水過程與模擬洪水過程擬合程度的確定性系數(shù)作為似然依據(jù),計算公式為:

(1)

步驟2確定模型參數(shù)的先驗分布,若不確定是哪種分布,以均勻分布表示。

步驟3根據(jù)在參數(shù)取值范圍內(nèi)生成的隨機(jī)參數(shù)組,利用式(1)計算似然值,并繪制參數(shù)與似然值點的函數(shù)分布圖。

1.3 Sobol法

Sobol法[11-12]作為一種基于方差的全局敏感性定量分析方法,能夠非常有效地分析高維非線性水文模型多個變量之間相互作用的敏感性,可根據(jù)Sobol法的分析結(jié)果決定每個模型輸入?yún)?shù)及其互相作用對整個模型輸出方差的貢獻(xiàn)度。

若用y=f(X)=f(x1,x2,…,xm)表示模型結(jié)構(gòu),x1、x2、…、xm表示模型參數(shù)(m為模型參數(shù)個數(shù)),則方差分解公式可表示為:

(2)

式中:V(y)為模型輸出y的總方差;Vi為第i個參數(shù)作用下的方差項;Vij為第i個和第j個參數(shù)共同作用下的方差項;V1,2,…,m為所有參數(shù)作用下的方差項。

定義一階敏感度Si、二階敏感度Sij與總敏感度STi指標(biāo)分別為:

(3)

(4)

(5)

式中V-i為不含第i個參數(shù)作用的方差項。一階敏感度Si表征單一模型參數(shù)組合對模型輸出的影響,二階敏感度Sij表征兩個模型參數(shù)組合對模型輸出的影響,總敏感度STi刻畫的是包括第i個模型參數(shù)在內(nèi)的所有參數(shù)組合對于模型輸出的影響。

2 實例分析

2.1 研究區(qū)概況

本次研究區(qū)域為白土崗水文站以上的白河上游流域,流域位于河南省境內(nèi),大部分面積在南陽市境內(nèi)和洛陽市欒川縣以及嵩縣的部分地區(qū),地理坐標(biāo)為東經(jīng)111°48′～112°25′、北緯33°19′～33°45′。白河上游流域地理位置如圖2所示。

圖2 白河上游流域地理位置圖

研究區(qū)地處北亞熱帶北部,是南北方氣候過渡地帶,冬季主要受亞洲大陸強(qiáng)冷高壓控制,寒冷高壓干燥少雨,夏季主要受西太平洋副熱帶高壓控制,西南季風(fēng)活動較為頻繁,溫度高,雨量充沛,為河南省3個典型的暴雨中心之一,極易形成暴雨災(zāi)害。研究區(qū)年降雨量為700～1 500 mm,多年平均降雨量為891.5 mm,降雨分布不均,主要集中在汛期,其降雨量約占全年降雨量的一半以上,且暴雨頻繁,經(jīng)常出現(xiàn)洪澇,無霜期長;年平均蒸發(fā)量為892～1 244 mm,平均值為1 000 mm,蒸發(fā)的年變化與溫度的年變化大致相同,通常12月至次年1月的蒸發(fā)量最小,而6—9月的蒸發(fā)量最大;河流徑流和降水的變化大體相似,具有明顯的季節(jié)性和區(qū)域性,同時在區(qū)域和時間上分布極不均勻。

2.2 數(shù)據(jù)來源

本文收集到的白河上游流域雨量站、水文站的數(shù)據(jù)情況見表2,雨量站與水文站的分布情況如圖2所示。

表2 研究區(qū)域雨量站、水文站資料情況表

2.3 流域洪水模擬

本次收集的降雨徑流資料相對齊全的年份為2012—2018年,因此截取此段時間內(nèi)白河上游流域22場洪水資料,將白土崗水文站的實測徑流資料以及周邊雨量站的降水資料利用線性插值程序插值為逐小時流量,再運(yùn)用時空變源混合產(chǎn)流模型進(jìn)行洪水模擬,模擬誤差見表3,部分洪水的模擬結(jié)果如圖3所示。

表3 白河上游流域洪水模擬結(jié)果指標(biāo)統(tǒng)計表

圖3 白河上游流域部分洪水模擬結(jié)果

對挑選的22場洪水的模擬結(jié)果進(jìn)行評判,將評判結(jié)果分為較好、可接受、較差3個等級。其中,較好等級的納什系數(shù)在0.8以上,且洪峰的相對誤差在10%以下;可接受等級的納什系數(shù)為0.6～0.8,洪峰相對誤差為10%～20%;較差等級的納什系數(shù)小于0.6,洪峰相對誤差大于20%。

經(jīng)統(tǒng)計:22場洪水中,8場洪水的模擬結(jié)果為較好,10場洪水的模擬結(jié)果為可接受,4場洪水的模擬結(jié)果為較差;平均納什系數(shù)為0.711,洪峰平均相對誤差為6.80%,可接受以上等級的洪水場次占所有洪水場次的81.8%。這說明模擬的場次洪水過程與實測洪水過程具有較好的一致性,可以認(rèn)為時空變源混合產(chǎn)流模型在由白土崗水文站控制的白河上游流域具有較好的適用性。

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

3.1 GLUE法分析結(jié)果

在時空變源混合產(chǎn)流模型參數(shù)的取值范圍內(nèi),運(yùn)用Matlab平臺編程,采用蒙特卡洛隨機(jī)取樣法生成2 000組隨機(jī)參數(shù)組,并將其輸入至?xí)r空變源混合產(chǎn)流模型中進(jìn)行模擬,確定以反映實測洪水與模擬洪水過程重合程度的確定性系數(shù)為似然值,并繪制參數(shù)與似然值點的分布圖,如圖4所示。圖4中,橫坐標(biāo)為各項參數(shù)的取值范圍,縱坐標(biāo)為各參數(shù)對應(yīng)的似然值(納什系數(shù)NS)。

圖4 GLUE法參數(shù)敏感度示意圖

從圖4中可以看出,可發(fā)生蓄滿產(chǎn)流的面積比pref_flow_den在(0,0.5)區(qū)間、表層土壤厚度layer_top_depth在(0,0.2)區(qū)間、快速壤中流線性出流系數(shù)fastcoef_lin在(0,0.3)區(qū)間、初始土壤含水量initial_water_content在整個區(qū)間都有明顯的變化,說明這4個參數(shù)在各自劃定的取值范圍內(nèi)非常敏感;而飽和水力傳導(dǎo)度conductivity、土壤厚度layer_depth、土壤向地下水庫最大排水系數(shù)soil2gw_max、慢速壤中流線性出流系數(shù)slowcoef_lin變化不明顯,說明取值范圍內(nèi)的參數(shù)組對洪水模擬影響并不大。由此,可以判斷出本次研究的敏感參數(shù)有pref_flow_den、layer_top_depth、fastcoef_lin、initial_water_content;非敏感參數(shù)有conductivity、layer_depth、soil2gw_max、slowcoef_lin。

3.2 Sobol法分析結(jié)果

同GLUE法,在時空變源混合產(chǎn)流模型的參數(shù)取值范圍內(nèi),使用蒙特卡洛隨機(jī)抽樣法,生成2 000組隨機(jī)參數(shù)組合,將其分別輸入至?xí)r空變源混合產(chǎn)流模型中,根據(jù)確定性系數(shù)來計算總敏感度STi(式(3)),選取敏感度大于0.01的參數(shù)為主要敏感參數(shù),計算結(jié)果如圖5所示。由圖5可知:fastcoef_lin、pref_flow_den、layer_top_depth、slowcoef_lin、initial_water_content、soil2gw_max為敏感參數(shù),其總敏感度分別為0.565、0.524、0.493、0.080、0.044、0.011,敏感度占比分別是33%、30%、29%、5%、2%、1%;非敏感參數(shù)為conductivity、layer_depth。

圖5 Sobol法參數(shù)敏感性分析圖

3.3 參數(shù)敏感性結(jié)果分析

根據(jù)GLUE法和Sobol法對應(yīng)的參數(shù)敏感性分析結(jié)果,可以判斷在白河上游流域,模型參數(shù)fastcoef_lin、pref_flow_den、layer_top_depth、initial_water_content為主要敏感參數(shù),其數(shù)值變動對結(jié)果有較大影響。由此說明:由白土崗水文站控制的白河上游流域的產(chǎn)流方式符合蓄滿產(chǎn)流機(jī)制,可劃歸為濕潤地區(qū),其主要徑流成分是蓄滿產(chǎn)流和壤中流,其次是地下徑流,超滲產(chǎn)流的比例最小,土壤特征對洪水過程影響顯著。

4 結(jié)論

本文以白河上游流域作為研究區(qū)域,運(yùn)用時空變源混合產(chǎn)流模型模擬了研究區(qū)2012—2018年間的22場洪水,并采用GLUE法和Sobol法對模型參數(shù)的敏感性進(jìn)行了分析。得到如下結(jié)論:

1)22場洪水的平均納什系數(shù)為0.711,洪峰平均相對誤差為6.80%,可接受以上標(biāo)準(zhǔn)的洪水場次占所有洪水場次的81.8%,模擬場次洪水過程與實測洪水過程吻合較好,說明用時空變源混合產(chǎn)流模型分析白河上游流域時具有較好的適用性。

2)使用GLUE法得出的敏感參數(shù)有pref_flow_den、layer_top_depth、fastcoef_lin、initial_water_content;使用Sobol法得出的敏感參數(shù)為astcoef_lin、pref_flow_den、layer_top_depth、slowcoef_lin、initial_water_content、soil2gw_max。結(jié)合實際計算過程和兩種參數(shù)敏感性方法分析結(jié)果認(rèn)為，白河上游流域的敏感參數(shù)為fastcoef_lin、pref_flow_den、layer_top_depth、initial_water_content。

3)白河上游流域以蓄滿產(chǎn)流方式為主。