基于Morris—Sobol的臨界雨量參數敏感性分析

原文林 高倩雨 張曉蕾 郝鵬

摘要：山洪災害臨界雨量計算參數眾多，部分取值受主觀因素影響較大。對此，應用敏感性分析理論，建立基于Morris-Sobol的臨界雨量計算參數敏感性分析二層模型。該模型首先利用修正的Morris篩選法對參數進行定性分析，以定性分析結果為基礎，選取中等靈敏度以上參數作為定量分析的參數，再通過Sobol指數法對選取參數進行定量敏感性分析。通過定性分析與定量計算相結合的分層遞進，分析計算各參數對臨界雨量的影響程度。以河南省安陽縣都里鄉都里前街小流域為研究對象，對其臨界雨量計算過程中的參數敏感性進行分析，結果表明：最大th點雨量均值H1的敏感性最大，說明在臨界雨量計算中，短歷時暴雨對計算結果的影響較大。不論是利用設計暴雨還是實測降雨進行臨界雨量計算，均應注意短歷時暴雨的資料精度。

關鍵詞：山洪災害；臨界雨量；參數敏感性分析；Moms-Sobol二層模型：水位流量反推法

中圖分類號：TV121+.1

文獻標志碼：A

doi： 10. 3969/j.issn.1000-1379.2018.07.008

1引言

山洪災害是具有成災迅速、破壞力強、危害嚴重等特點的一種自然災害，國內外一般選取臨界雨量作為山洪災害預警預報的重要指標。臨界雨量是指導致一個流域或區域發生山溪洪水可能致災時，降雨量達到或超過的量級和強度

。

目前，國外最具代表性的臨界雨量計算方法為美國的Flash Flood Guidance（ FFG）方法。國內山洪災害臨界雨量計算方法主要有數據驅動法和水文水力學法。我國山丘區面積占國土面積的近三分之二，保證山丘區小流域預警預報質量才能有效減少災害損失，而山丘區小流域雨量、水文等監測設施不完善，缺乏實測資料，臨界雨量計算通常采用水文水力學法。在應用水文水力學法時，部分計算參數取值受主觀因素影響較大，需通過參數敏感性及其相互關聯性分析，辨識參數的敏感程度，確定各參數對臨界雨量計算結果的影響，為臨界雨量計算及參數率定提供參考依據。

敏感性分析一般分為局部分析和全局分析，局部敏感性分析僅針對單個參數對計算結果的影響，全局敏感性分析考慮整個參數空間對計算結果的影響8。全局敏感性分析方法又可分為定性和定量兩類方法，定性分析法包括Morris篩選法、多元回歸法和LH -OAT方法等，定量分析法包括Sobol指數法、ExtendFAST法和GLUE方法等。修正的Morris篩選法利用微分計算各參數的敏感性，具有計算量小且能夠處理多個輸人參數的優點，但缺乏參數間交互效應的分析：Sobol指數法基于方差分解的原理，具有對單個參數的主效應、全效應及多個參數的交互效應進行分析的優點，在很多領域得到應用，但該方法計算量較大。

山洪災害臨界雨量計算涉及參數眾多，為了研究各參數對臨界雨量的影響程度以及各參數間的交互效應對臨界雨量的影響，本文利用修正的Moms篩選法計算量小的優點對各參數進行敏感性定性分析，結合Sobol指數法的定量計算，建立基于Morris-Sobol的臨界雨量計算參數敏感性分析二層模型。以河南省安陽縣都里鄉都里前街小流域的臨界雨量計算為例，結合臨界雨量計算模型中的參數，辨識計算參數對臨界雨量的影響程度及相互關聯性，為臨界雨量計算參數取值提供參考依據。

2臨界雨量計算參數分析

水位流量反推法是水文水力學法中最典型的方法，本文主要針對水位流量反推法計算臨界雨量中涉及的參數進行敏感性分析與討論。水位流量反推法假定降雨與洪水同頻率，根據河道控制斷面成災水位，由水位流量關系計算對應的成災流量，由流量頻率曲線確定成災洪水頻率，最后由降雨頻率曲線確定臨界雨量。水位流量反推法計算臨界雨量主要包括設計暴雨計算、設計洪水計算和水位流量關系推求三部分。

（1）設計暴雨計算參數。山丘區小流域設計暴雨主要是基于當地暴雨圖集等資料進行計算，即暴雨圖集法。根據研究對象所處地理位置，利用暴雨圖集確定相應參數取值并進行設計暴雨計算。涉及的參數主要有最大24h點雨量均值H24、最大1h點雨量均值H1、變差系數Cv和暴雨遞減指數n1、n2、n3。

（2）設計洪水計算參數。小流域設計洪水的計算方法主要有推理公式法、等流時線法和瞬時單位線法等，本文采用推理公式法進行設計洪水計算，其中涉及的參數為降雨損失參數μ。

（3）水位流量關系參數。水位流量關系推求常用方法為曼寧公式法。曼寧公式中涉及的參數主要有糙率n和水面比降i。

3基于Morris-Sobol的臨界雨量計算參數敏

感性分析二層模型

在進行參數敏感性分析時，為了提高計算效率，筆者提出了Morris篩選法定性分析和Sobol指數法定量計算的臨界雨量計算參數敏感性分析二層模型，其基本流程見圖1。

3.1定性分析

山洪災害臨界雨量計算中的參數敏感性定性分析采用修正的Morris篩選法，其基本原理為選取山洪災害臨界雨量計算模型中的某一參數xi，以固定步長C進行擾動變化，最大變幅為M，計算其靈敏度判別因子S，其計算公式為式中：Yo為參數xi初始取值對應的臨界雨量；Yi、Yi+1分別為參數xi第i次、第i+1次擾動變化后對應的臨界雨量；Pi、Pi+1分別為參數xi第i次、第i+1次擾動變化后的取值相對于初始值的變化率，%：Z為擾動變化總次數，即試驗次數，其值由固定步長C和最大變幅M決定。

根據各參數的敏感性判別因子進行參數敏感性分級，敏感性分級標準見表1。

根據各參數敏感性定性分析結果，即敏感性分級情況，去除不靈敏參數，并對剩余參數進行敏感性定量分析。

3.2定量分析

山洪災害臨界雨量計算中的參數敏感性定量分析采用Sobol指數法，其基本原理是將臨界雨量計算系統分解成單個參數及參數間相互組合的函數，計算單個參數和參數組合的方差，并計算分析各部分方差對總方差的影響，即分解方差，以此來分析參數的重要程度和參數間的相互影響程度

。

假設Y=f（x）=f（x1，x2，…，xi，…，xn），其中：Y 為臨界雨量，xi為敏感性定量分析選取的參數。根據方差分解理論有

主效應指標Sxi反映單個參數xi單獨對臨界雨量的影響，主效應指標越大，說明該參數對臨界雨量的影響越大。全效應指標ST包含了參數xi的主效應和該參數與其他參數的交互作用，若參數的全效應指標較大，說明該參數不僅對臨界雨量影響較大，而且與其他參數的交互效應也很大。

采用蒙特卡洛法模擬得到參數的一階、二階及更高階次的敏感度。在利用修正的Morris篩選法進行定性分析的基礎上，根據靈敏度判別因子S選取山洪災害臨界雨量計算中的k個參數進行定量分析，在k個參數的取值范圍內利用拉丁超立方抽樣t次，抽取A和B兩個矩陣：

矩陣的每一行表示k個參數的組合，每一列表示t次隨機抽取參數xi的值，將式（5）中矩陣A的第i列換成矩陣B的第i列，其余列保持不變記為Ci；再將矩陣B的第i列換成矩陣A的第i列，其余列保持不變，得到矩陣記為C-i。

利用A、B、Ci、C-i矩陣中的每一列參數分別進行臨界雨量計算，并按照式（7）至式（12）進行相應的方差估計和參數敏感性指標計算。式中：f為矩陣A和B均值的估計值；V（Y）為臨界雨量的方差估計值；Ui為矩陣A和Ci均值的估計值；U-i為A和C-i均值的估計值；Sxi為主效應指標的估計值；ST為全效應指標的估計值。

4實例分析

4.1研究區概況

以河南省安陽縣都里鄉都里前街小流域為研究對象，分析山洪災害臨界雨量計算中的參數敏感性。安陽縣地處河南省北部，北緯35°35-36°21'，東經113°35'-114°45，地勢由西到東呈階梯降低，有山區、丘陵、平原、洼地多種地表形態。安陽縣屬暖溫帶季風氣候區，全年氣溫變化較大，四季分明，因受季風影響，故降雨時空分布極不均勻，多年平均降雨量510mm。研究小流域形狀見圖2，小流域面積為69.11km^2，河長為21.07km。研究對象控制斷面成災水位為220.15m。

根據《河南省中小流域設計暴雨洪水圖集》，結合研究區自然地理、地形地貌、水文氣象等情況及其所在水文分區，確定各分析參數取值范圍及初始值，見表2。

4.2定性分析結果

根據修正的Morris篩選法，以步長C=0.5%對某一參數值進行擾動，擾動最大變幅M=±20%，且保證其余參數固定不變。計算各參數的靈敏度判別因子，結果見表3，敏感參數分布見圖3。

由表3可知，山洪災害臨界雨量計算中不靈敏參數有3個，分別是變差系數Cv，、暴雨遞減指數n1、暴雨遞減指數n2；中等靈敏參數有1個，為降雨損失參數μ；靈敏參數有5個，分別是最大24 h點雨量均值H24、最大1h點雨量均值H1、暴雨遞減指數n2、糙率n和水面比降i。

在模型第一層定性分析結果的基礎上，去除不靈敏參數，減少定量分析輸入參數，提高定量分析計算效率，同時對研究系統中的參數進行充分的分析計算，初步選取中等靈敏參數及以上等級共6個參數作為定量分析的參數，即降雨損失參數μ、最大24h點雨量均值H24、最大1 h點雨量均值H，、暴雨遞減指數n2、糙率n和水面比降i。

4.3定量分析結果

根據選定的6個敏感性定量分析參數及其取值范圍，進行拉丁超立方抽樣，得到2000個樣本，采用Sobol指數法計算6個參數的主效應指標和全效應指標估計值，結果見表4。

由表4可見，6個參數的敏感性由強到弱依次為最大1h點雨量均值H、糙率n、水面比降i、暴雨遞減指數n2降雨損失參數μ、最大24h點雨量均值H24。

最大1h點雨量均值H，不僅自身的變化對臨界雨量計算結果的影響遠大于其余參數，且其與其他參數的交互效應對臨界雨量的影響程度也最大，這說明在臨界雨量計算中，短歷時暴雨對計算結果的影響較大。因此，不論是利用設計暴雨還是實測降雨進行臨界雨量計算，均應注意短歷時暴雨的數據精度。

糙率n對臨界雨量的影響程度也較大。在實際工作中，其取值通常根據河道特征參照天然河道糙率表進行確定。為了提高臨界雨量計算精度，在確定山洪災害臨界雨量計算參數值的時候，應重點分析糙率n的取值。當實測水文資料充足時，應利用實測水文資料進行推算：當實測水文資料缺乏時，應綜合考慮河流的溝道特征、植被生長和床面粗糙情況等多方面因素，選取一個最佳值。

在確定山洪災害臨界雨量計算參數取值時，按照參數的敏感性大小排序情況，加強對敏感性較強參數的研究，適當減少對敏感性較弱參數的研究，提高臨界雨量計算精度，保證山洪災害預謦預報精準度。

此外，通過基于Morris-Sobol的臨界雨量計算參數敏感性分析二層模型對計算中涉及的設計暴雨參數、設計洪水參數和水位流量關系參數進行敏感性分析，綜合分析結果可知，水位流量關系參數糙率n和水面比降i的敏感性分別排在第2位和第3位，對臨界雨量的影響均較大。在確定山洪災害臨界雨量計算參數取值時，應加強水位流量關系參數的研究，提高山洪災害預警預報質量。

5結論

利用Morris篩選法和Sobol指數法建立敏感性分析模型，通過定性分析與定量計算相結合的分層遞進，分析計算參數對臨界雨量的影響程度，可為山洪災害臨界雨量計算時參數取值提供參考依據。以河南省安陽縣都里前街小流域為例的分析，結果表明：在山洪災害臨界雨量計算參數中，最大1h點雨量均值H，對山洪災害臨界雨量的影響程度最大，因此在山洪災害臨界雨量計算時應注意短歷時暴雨的資料精度。

此外，通過對各參數的敏感性進行綜合分析，水位流量關系參數的敏感性較強，因此應對水位流量關系參數進行重點研究，尤其是糙率n對山洪災害臨界雨量的影響程度較大，應利用實測水文資料進行推算，或綜合考慮河流的溝道特征、植被生長和床面粗糙情況等多方面因素，選取最佳值。

本研究僅針對最為常用的水位流量反推法進行了參數敏感性討論，在今后的分析研究中可針對多種方法分別研究不同方法下的敏感性參數，以便實際應用選取不同方法時作為參考。