基于Copula 函數(shù)風(fēng)險(xiǎn)組合的西峽縣山洪災(zāi)害預(yù)警分析

張玉昶，榮富強(qiáng)，王明閃，王建國(guó)，萬(wàn) 芳

（1.河南省地質(zhì)礦產(chǎn)勘查開(kāi)發(fā)局第二地質(zhì)環(huán)境調(diào)查院，河南鄭州 450000；2.華北水利水電大學(xué)水資源學(xué)院，河南鄭州 450011）

近年來(lái)，極端氣候頻發(fā)，導(dǎo)致山洪災(zāi)害越來(lái)越頻繁發(fā)生，山洪災(zāi)害已經(jīng)成為世界各國(guó)自然災(zāi)害中的主要災(zāi)種之一，造成了大量的人員傷亡和社會(huì)經(jīng)濟(jì)損失［1］，為了防御山洪災(zāi)害，減輕災(zāi)害損失，國(guó)內(nèi)外許多專(zhuān)家學(xué)者對(duì)山洪災(zāi)害預(yù)警進(jìn)行了研究。王協(xié)康等［2］采用TRI-GRS 模型分析了小流域的滑坡易發(fā)性，將暴雨山洪洪水劃分為預(yù)警區(qū)、中、高風(fēng)險(xiǎn)區(qū)。我國(guó)下墊面等條件復(fù)雜多樣，山洪災(zāi)害成災(zāi)機(jī)理復(fù)雜，而臨界雨量是山洪災(zāi)害的主要影響因素。目前，臨界雨量的計(jì)算方法主要有數(shù)據(jù)驅(qū)動(dòng)法及水文水力學(xué)法［3，4］。但是，臨界雨量受前期影響雨量影響較大［5］，江錦紅等［6］采用雙曲函數(shù)確定的暴雨臨界雨量最為山洪災(zāi)害的預(yù)警標(biāo)準(zhǔn)；劉媛媛等［7］針對(duì)小流域，利用臨界流量及水位計(jì)算相應(yīng)的臨界頻率；張珊珊等［8］基于HEC-HMS 模型，研究選用雨量預(yù)警指標(biāo)，綜合考慮土壤含水量、匯流時(shí)間等，反推臨界雨量。臨界雨量的研究成果豐碩，但是未考慮到幾種風(fēng)險(xiǎn)組合情況的臨界雨量及相應(yīng)的預(yù)警分析。水文科學(xué)領(lǐng)域的種隨機(jī)變量之間往往呈現(xiàn)出各種復(fù)雜的線(xiàn)性、非線(xiàn)性的相關(guān)關(guān)系，本文引入Copula 函數(shù)來(lái)描述降雨量與峰值雨強(qiáng)變量間的相關(guān)性，Copula 函數(shù)可將聯(lián)合分布的函數(shù)和它各自的邊緣分布函數(shù)聯(lián)系在一起，基于變量之間的非線(xiàn)性相關(guān)關(guān)系而建立的，可以描述變量間非線(xiàn)性、非對(duì)稱(chēng)和對(duì)稱(chēng)的相關(guān)關(guān)系，目前在水文科學(xué)中得到廣泛的應(yīng)用［9-11］。

1 研究區(qū)概況

西峽縣位于河南省西南邊陲，與陜西省接界，地處秦嶺山系東南余脈伏牛山南。東南與本市的內(nèi)鄉(xiāng)、淅川相連，北西與欒川、盧氏和陜西省的商南縣毗鄰。地處北亞熱帶向暖溫帶的過(guò)度地區(qū)，屬于亞熱帶季風(fēng)型大陸性氣候，季風(fēng)影響明顯，四季分明，春季溫度升降劇烈，雨水分布不均，夏季熱雨集中，暴雨頻發(fā)，地質(zhì)地貌復(fù)雜，加上受人類(lèi)活動(dòng)的影響，導(dǎo)致山洪災(zāi)害發(fā)生頻繁。山洪災(zāi)害在縣境內(nèi)分布廣，具有來(lái)勢(shì)猛、流速快、破壞力大、突發(fā)性強(qiáng)的特點(diǎn)，預(yù)報(bào)、預(yù)測(cè)、預(yù)防難度較大，不僅對(duì)山區(qū)的基礎(chǔ)設(shè)施造成毀滅性破壞，而且對(duì)人民群眾的生命安全構(gòu)成極大的威脅，已經(jīng)成為當(dāng)前防災(zāi)減災(zāi)工作中的突出問(wèn)題，是山區(qū)經(jīng)濟(jì)社會(huì)可持續(xù)發(fā)展的重要制約因素之一。

西峽縣山洪災(zāi)害致災(zāi)因素具有自然地理和經(jīng)濟(jì)社會(huì)的雙重屬性，降雨因素是誘發(fā)山洪災(zāi)害的直接因素和激發(fā)條件。降雨量大，多數(shù)情況下意味著雨強(qiáng)高、激發(fā)力強(qiáng)，在一定的下墊面條件下，易產(chǎn)生溪河洪水災(zāi)害、泥石流災(zāi)害和滑坡災(zāi)害。根據(jù)西峽縣近年來(lái)的統(tǒng)計(jì)資料表明，每次山洪災(zāi)害的發(fā)生，都是在降雨量較大這一基本前提條件下形成的；高強(qiáng)度的降雨是引起山洪災(zāi)害的最主要原因之一。

在相同的條件下，降雨歷時(shí)越長(zhǎng)，降雨量越大，產(chǎn)生的徑流量越大，山洪災(zāi)害損失也越嚴(yán)重。在溪河洪水、泥石流、滑坡三種山洪災(zāi)害類(lèi)型中，溪河洪水災(zāi)害和泥石流災(zāi)害特點(diǎn)相似，在下墊面條件滿(mǎn)足的情況下，只要有足夠大的降雨強(qiáng)度和降雨量，溪河洪水和泥石流災(zāi)害就可發(fā)生，降雨歷時(shí)越長(zhǎng)，所產(chǎn)生的洪量越大，災(zāi)害損失越大。滑坡災(zāi)害與降雨歷時(shí)的關(guān)系更為明顯，一般而言，滑坡和降雨并不是同時(shí)發(fā)生，滑坡一般滯后于降雨，由溪河洪水誘發(fā)的滑坡災(zāi)害會(huì)隨著降雨歷時(shí)的增長(zhǎng)，洪流對(duì)坡腳的掏蝕，滑坡也會(huì)越來(lái)越多。

2 混合Copula 函數(shù)的風(fēng)險(xiǎn)組合概率的臨界雨量確定

2.1 Copula理論及Copula函數(shù)分布模型

1959年Copula函數(shù)理論提出，Sklar提出可把一個(gè)聯(lián)合分布函數(shù)分解為一個(gè)Copula 函數(shù)和K個(gè)邊緣分布函數(shù)，而分解出來(lái)的Copula 函數(shù)可用來(lái)描述變量之間相關(guān)性，也就是可將聯(lián)合分布的函數(shù)和它各自的邊緣分布函數(shù)聯(lián)系在一起的函數(shù)。N元Copula函數(shù)C(u1，u2，…，uN)具有以下性質(zhì)：

常用的Copula 數(shù)有很多種，本文主要介紹常用的多元正態(tài)Copula 函數(shù)、阿基米德Copula 函數(shù)Archimedean 的三類(lèi)常用函數(shù)（Clayton 函數(shù)、Gumbel Copula 函數(shù)、Frank Copula 函數(shù)）和極值Copula函數(shù)。

（1）多元正態(tài)Copula函數(shù)。其中N元正態(tài)Copula函數(shù)的分布函數(shù)與密度函數(shù)可表達(dá)為：

式中：ρ為對(duì)角線(xiàn)上的元素為1的對(duì)稱(chēng)正定矩陣；|ρ|表示與矩陣ρ相對(duì)應(yīng)的行列式的值；Φρ(?)表示相關(guān)系數(shù)矩陣為ρ的標(biāo)準(zhǔn)正態(tài)分布函數(shù)；Φ-1(?)為標(biāo)準(zhǔn)多元正態(tài)函數(shù)Φ(?)的逆函數(shù)；?=(?1，?2，…，?N)′，?n=Φ-1(un)，n=1，2，…，N，I為單位矩陣。

（2）阿基米德Copula函數(shù)。其中阿基米德Copula函數(shù)可表示為：

Gumbel Copula 函數(shù)、Clayton Copula 函數(shù)、Frank Copula 函數(shù)是常用的二元阿基米德Copula 函數(shù)，分別由它們擴(kuò)展到N元阿基米德Copula函數(shù)可表示為：

①N元Gumbel Copula函數(shù)的表達(dá)式為：

②N元Clayton Copula函數(shù)的表達(dá)式為：

③N元Frank Copula函數(shù)的表達(dá)式為：

其中，α，θ，λ為相關(guān)參數(shù)。

（3）極值Copula 函數(shù)。極值Copula 函數(shù)（Extreme Value Copula，EVC）可表達(dá)為：

2.2 混合Copula函數(shù)的構(gòu)造與相關(guān)性分析

（1）不同類(lèi)型的Copula 函數(shù)比較分析。一般情況下，多元正態(tài)Copula 函數(shù)通常描述變量之間的相關(guān)關(guān)系，但是由于此函數(shù)有對(duì)稱(chēng)性的特點(diǎn)，所以對(duì)于變量的非對(duì)稱(chēng)相關(guān)關(guān)系很難擬合。

而Gumbel Copula 函數(shù)和Clayton Copula 函數(shù)都具有非對(duì)稱(chēng)性的特點(diǎn)，其中Gumbel Copula函數(shù)為“J”字形分布，下尾低而上尾高，對(duì)水文變量的下尾部變化不太敏感，而對(duì)上尾部的分布變化比較敏感，所以很難描述下尾部的相關(guān)變化情況，即當(dāng)一個(gè)水文變量出現(xiàn)極大值時(shí)，另外兩個(gè)水文變量也出現(xiàn)極大值的概率增大；

Clayton Copula 函數(shù)為“L”字形分布，下尾高而上尾低，對(duì)水文變量的下尾部變化比較敏感，而對(duì)上尾部的分布變化不太敏感，所以很難描述上尾部的相關(guān)變化情況，即當(dāng)一個(gè)水文變量出現(xiàn)極小值時(shí)，其他4個(gè)水文變量也出現(xiàn)極小值的概率增大。

Frank Copula 函數(shù)為“U”字形分布，它有對(duì)稱(chēng)性的特點(diǎn)，所以很難描述水文變量之間的非對(duì)稱(chēng)關(guān)系，F(xiàn)rank Copula 函數(shù)只適合描述具有對(duì)稱(chēng)相關(guān)結(jié)構(gòu)的變量之間的相關(guān)關(guān)系，即各個(gè)水文變量間極大值相關(guān)性與極小值相關(guān)性是對(duì)稱(chēng)增長(zhǎng)的，但是其尾部的分布變量是比較獨(dú)立的，所以Frank Copula 函數(shù)無(wú)論在描述上尾部還是下尾部的相關(guān)性中都是不敏感的，故無(wú)法描述尾部變化的相關(guān)性。

（2）混合Copula 函數(shù)的構(gòu)造與相關(guān)性分析。Gumbel Copu‐la 函數(shù)、Clayton Copula 函數(shù)、Frank Copula 函數(shù)三類(lèi)常用的阿基米德函數(shù)能夠捕捉尾部相關(guān)的情形：上尾部的相關(guān)、下尾部的相關(guān)、上尾部和下尾部的對(duì)稱(chēng)相關(guān)。這些Copula 函數(shù)具有描述各種水文變量模式之間的關(guān)系，特別是尾部相關(guān)關(guān)系，水文系統(tǒng)中各變量之間的關(guān)系是復(fù)雜多變的，不是拘泥于某種特定關(guān)系，它們只能反映水文變量間相關(guān)性的某個(gè)側(cè)面，因此很難用一個(gè)簡(jiǎn)單的Copula 函數(shù)來(lái)全面的刻畫(huà)水文系統(tǒng)中各變量之間的相關(guān)模式，所以一個(gè)更加靈活的Copula 函數(shù)需要構(gòu)造，來(lái)描述各種水文變量模式之間的關(guān)系。應(yīng)用不同Copula 函數(shù)的優(yōu)點(diǎn)，本文選用Gumbel Copula 函數(shù)、Clayton Copula 函數(shù)、Frank Copula 函數(shù)的線(xiàn)性組合來(lái)構(gòu)造混合的Copula 函數(shù)，可以更加靈活的刻畫(huà)水文系統(tǒng)中各變量之間的相關(guān)關(guān)系。

混合Copula函數(shù)M-Copula表達(dá)式為：

式中：CG，CF，Ccl分別為Gumbel Copula、Frank Copula、Clayton Copula 函數(shù)；wG，wF，wcl為相應(yīng)Copula 函數(shù)的權(quán)重系數(shù)。由式（8）～（10）可知，MC3中包括了6 個(gè)參數(shù)，參數(shù)向量(α，λ，θ)用來(lái)描述變量間相關(guān)的程度，變量間相關(guān)的模式由線(xiàn)性權(quán)重參數(shù)向量(wG，wF，wcl)表示。

2.3 Copula模型的參數(shù)估計(jì)、檢驗(yàn)與評(píng)價(jià)

（1）Copula 模型的參數(shù)估計(jì)。Copula 函數(shù)模型中的參數(shù)估計(jì)有很多種，一般采用極大似然和矩估計(jì)，而其中極大似然估計(jì)是較常用的Copula 模型參數(shù)估計(jì)方法，其聯(lián)合分布函數(shù)的密度函數(shù)為：

式中：θc為Copula 函數(shù)的1 ×mc維參數(shù)向量；Fn(xn；θn)為邊緣分布函數(shù)；θn為邊緣分布函數(shù)Fn(xn；θn)的1 ×mn維參數(shù)向量；θ=(θ1，θ2，…，θN，θc)′；n=1，2，…，N。

因此，可以得到樣本(x1t，x2t，…，xNt)，t=1，2，…，T的對(duì)數(shù)似然函數(shù)為：

使似然函數(shù)取的最大值的θ即是最大似然估計(jì)值。

（2）Copula 模型的檢驗(yàn)與評(píng)價(jià)。應(yīng)用的Copula 分布函數(shù)是否能夠很好的擬合變量之間的相關(guān)結(jié)構(gòu)及分布，所以Copula 函數(shù)的檢驗(yàn)與擬合優(yōu)度評(píng)價(jià)需要建立。K-S用于檢驗(yàn)樣本是否服從同一分布，故應(yīng)用其對(duì)Copula 分布函數(shù)進(jìn)行檢驗(yàn)；Copula 函數(shù)的擬合度評(píng)價(jià)采用均方根誤差RSME最小準(zhǔn)則來(lái)計(jì)算，其定義如式（14）：

式中：N是樣本容量；i為樣本序號(hào)；pc為模型計(jì)算的理論頻率；p0為聯(lián)合分布的經(jīng)驗(yàn)頻率。

其中K-S檢驗(yàn)的統(tǒng)計(jì)量D計(jì)算公式如（15）式所示：

式中：N為樣本容量；Ck為樣本xk=(x1k，x2k，x3k)的Copula 值；uk為樣本中滿(mǎn)足條件x≤xk的個(gè)數(shù)，即滿(mǎn)足：x1≤x1k，x2≤x2k，x3≤x3k。

2.4 風(fēng)險(xiǎn)組合概率的臨界雨量確定

設(shè)隨機(jī)變量X，Y分別服從邊緣分布FX(x)、FY(y)，其聯(lián)合分布函數(shù)為F(x，y)，其中X表示總雨量，Y表示峰值雨強(qiáng)，則存在一個(gè)Copula函數(shù)：

臨界雨量中總雨量和峰值雨強(qiáng)存在多種耦合情況，可以運(yùn)用上述的聯(lián)合概率分布模型計(jì)算出它們的聯(lián)合分布概率，本文主要考慮超過(guò)概率和條件概率，即當(dāng)總雨量X超過(guò)某特定值的條件下，峰值雨強(qiáng)超過(guò)某特定值時(shí)該雨型發(fā)生的條件概率，如式（17）所示。

式中：XA表示概率為AA的總雨量值；y表示某特定峰值雨強(qiáng)值。

由此，雨型風(fēng)險(xiǎn)組合概率(A，B)的定義為：

在山洪災(zāi)害預(yù)警預(yù)報(bào)中，通常利用重現(xiàn)期對(duì)暴雨量級(jí)進(jìn)行描述，對(duì)于超過(guò)概率，重現(xiàn)期計(jì)算公式如下式所示。

3 西峽縣基于臨界雨量的山洪災(zāi)害預(yù)警分析

3.1 Copula函數(shù)的選擇

分別采用Gumbel、Clayton、Frank 及混合Copula（M-Copula）函數(shù)構(gòu)建12 h 場(chǎng)次降雨量和峰值雨強(qiáng)聯(lián)合概率分布模型，K-S檢驗(yàn)統(tǒng)計(jì)量及擬合優(yōu)度評(píng)價(jià)指標(biāo)計(jì)算結(jié)果如表1所示。

表1 計(jì) 算 表 明，Gumbel-Copula、Clayton-Copula、Frank-Copula、M-Copula 函數(shù)均通過(guò)K-S 檢驗(yàn)，能較好的擬合臨界雨強(qiáng)的邊緣分布，根據(jù)均方根誤差RSME最小準(zhǔn)則，選取M-Copu‐la函數(shù)為連接函數(shù)來(lái)擬合降雨量與峰值雨強(qiáng)聯(lián)合分布情況。

表1 Copula概率分布函數(shù)的計(jì)算、檢驗(yàn)與評(píng)價(jià)結(jié)果Tab.1 Calculation results，inspection and evaluation of Copula cumulative distribution function

3.2 組合概率的臨界雨量的風(fēng)險(xiǎn)預(yù)警計(jì)算

綜合考慮防災(zāi)對(duì)象所處河段河谷形態(tài)、洪水上漲速率、預(yù)警響應(yīng)時(shí)間和站點(diǎn)位置等因素，在臨界雨量的基礎(chǔ)上綜合確定準(zhǔn)備轉(zhuǎn)移和立即轉(zhuǎn)移的預(yù)警指標(biāo)；并利用該預(yù)警指標(biāo)進(jìn)行暴雨洪水復(fù)核校正，以避免與成災(zāi)水位及相應(yīng)的暴雨洪水頻率差異過(guò)大。

通常情況下，由于臨界雨量是從成災(zāi)水位對(duì)應(yīng)流量的洪水推算得到，故在數(shù)值上認(rèn)為臨界雨量即立即轉(zhuǎn)移指標(biāo)；對(duì)于準(zhǔn)備轉(zhuǎn)移指標(biāo)，是在臨界雨量基礎(chǔ)上根據(jù)準(zhǔn)備轉(zhuǎn)移時(shí)間及洪水過(guò)程線(xiàn)綜合進(jìn)行“折減”處理。

以軍馬河鄉(xiāng)白果村回龍灣組為例，在不同的土壤含水量的情況下，利用風(fēng)險(xiǎn)概率雨型分別計(jì)算1、3、6、24 h 臨界雨量，其臨界雨量預(yù)警閾值見(jiàn)表2 所示；水位預(yù)警閾值見(jiàn)圖1 所示；準(zhǔn)備轉(zhuǎn)移預(yù)警和立即轉(zhuǎn)移預(yù)警雨量臨界線(xiàn)圖見(jiàn)圖2所示。

圖1 水位預(yù)警閾值圖Fig.1 Threshold map of water level warning

圖2 準(zhǔn)備轉(zhuǎn)移預(yù)警和立即轉(zhuǎn)移預(yù)警雨量臨界線(xiàn)圖Fig.2 Threshold maps for the preparation and immediate diversion of early warning rainfall

表2 臨界雨量預(yù)警閾值表Tab.2 Table of critical rainfall warning thresholds

3.3 結(jié)果合理性分析

對(duì)于山丘區(qū)河流洪峰模數(shù)具有一定的規(guī)律，流域面積越大，洪峰模數(shù)越小，面積越小，洪峰模數(shù)越大。由此可以根據(jù)河南省山丘區(qū)水文斷面的洪峰模數(shù)統(tǒng)計(jì)資料，間接評(píng)估設(shè)計(jì)洪峰流量的合理性。根據(jù)以往經(jīng)驗(yàn)資料，50年一遇的中小河流斷面的設(shè)計(jì)洪峰流量對(duì)應(yīng)的洪峰模數(shù)一般在10.0 左右，并且10.0 以下占多數(shù)?；佚垶辰M計(jì)算得到洪峰模數(shù)為8.86，因此設(shè)計(jì)流量計(jì)算成果是合理的。

表3 回龍灣組計(jì)算成果合理性分析Tab.3 Analysis on rationality of calculation results of Huilongwan

圖3 各級(jí)水位對(duì)比圖Fig.3 Comparison chart of water level at different stages

4 結(jié)論

通過(guò)對(duì)西峽縣的山洪災(zāi)害預(yù)警分析可知：混合Copula 函數(shù)可以更加靈活的刻畫(huà)水文系統(tǒng)中各變量之間的相關(guān)關(guān)系，Copula 分布函數(shù)能夠很好的擬合總雨量和峰值雨強(qiáng)的二維聯(lián)合分布函數(shù)變量之間的相關(guān)結(jié)構(gòu)及分布，根據(jù)均方根誤差RSME最小準(zhǔn)則，選取M-Copula 函數(shù)為連接函數(shù)來(lái)擬合臨界雨量聯(lián)合分布情況；基于對(duì)西峽縣軍馬河鄉(xiāng)白果村回龍灣組為例進(jìn)行的計(jì)算及分析，計(jì)算組合概率的臨界雨量，并對(duì)山洪災(zāi)害的水位、準(zhǔn)備轉(zhuǎn)移預(yù)警及立即轉(zhuǎn)移預(yù)警雨量臨界進(jìn)行預(yù)警及分析。結(jié)果表明：對(duì)于山丘區(qū)河流洪峰模數(shù)具有一定的規(guī)律，流域面積越大，洪峰模數(shù)越小，面積越小，洪峰模數(shù)越大，設(shè)計(jì)流量為639 m3∕s 的計(jì)算成果是合理的。因此，對(duì)山洪災(zāi)害預(yù)警分析要最大限度的保證人民群眾的安全，分析評(píng)價(jià)中以各沿河村落漲洪歷時(shí)為基礎(chǔ)，分析沿河村落暴雨洪水預(yù)警響應(yīng)時(shí)間及其分布情況，這對(duì)山洪災(zāi)害防治預(yù)案編制及預(yù)警等均具有重要作用。