張玉昶,榮富強(qiáng),王明閃,王建國,萬 芳
(1.河南省地質(zhì)礦產(chǎn)勘查開發(fā)局第二地質(zhì)環(huán)境調(diào)查院,河南鄭州 450000;2.華北水利水電大學(xué)水資源學(xué)院,河南鄭州 450011)
近年來,極端氣候頻發(fā),導(dǎo)致山洪災(zāi)害越來越頻繁發(fā)生,山洪災(zāi)害已經(jīng)成為世界各國自然災(zāi)害中的主要災(zāi)種之一,造成了大量的人員傷亡和社會(huì)經(jīng)濟(jì)損失[1],為了防御山洪災(zāi)害,減輕災(zāi)害損失,國內(nèi)外許多專家學(xué)者對(duì)山洪災(zāi)害預(yù)警進(jìn)行了研究。王協(xié)康等[2]采用TRI-GRS 模型分析了小流域的滑坡易發(fā)性,將暴雨山洪洪水劃分為預(yù)警區(qū)、中、高風(fēng)險(xiǎn)區(qū)。我國下墊面等條件復(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ù)雜的線性、非線性的相關(guān)關(guān)系,本文引入Copula 函數(shù)來描述降雨量與峰值雨強(qiáng)變量間的相關(guān)性,Copula 函數(shù)可將聯(lián)合分布的函數(shù)和它各自的邊緣分布函數(shù)聯(lián)系在一起,基于變量之間的非線性相關(guān)關(guān)系而建立的,可以描述變量間非線性、非對(duì)稱和對(duì)稱的相關(guān)關(guān)系,目前在水文科學(xué)中得到廣泛的應(yīng)用[9-11]。
西峽縣位于河南省西南邊陲,與陜西省接界,地處秦嶺山系東南余脈伏牛山南。東南與本市的內(nèi)鄉(xiāng)、淅川相連,北西與欒川、盧氏和陜西省的商南縣毗鄰。地處北亞熱帶向暖溫帶的過度地區(qū),屬于亞熱帶季風(fēng)型大陸性氣候,季風(fēng)影響明顯,四季分明,春季溫度升降劇烈,雨水分布不均,夏季熱雨集中,暴雨頻發(fā),地質(zhì)地貌復(fù)雜,加上受人類活動(dòng)的影響,導(dǎo)致山洪災(zāi)害發(fā)生頻繁。山洪災(zāi)害在縣境內(nèi)分布廣,具有來勢猛、流速快、破壞力大、突發(fā)性強(qiáng)的特點(diǎn),預(yù)報(bào)、預(yù)測、預(yù)防難度較大,不僅對(duì)山區(qū)的基礎(chǔ)設(shè)施造成毀滅性破壞,而且對(duì)人民群眾的生命安全構(gòu)成極大的威脅,已經(jīng)成為當(dāng)前防災(zāi)減災(zāi)工作中的突出問題,是山區(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ù)西峽縣近年來的統(tǒng)計(jì)資料表明,每次山洪災(zāi)害的發(fā)生,都是在降雨量較大這一基本前提條件下形成的;高強(qiáng)度的降雨是引起山洪災(zāi)害的最主要原因之一。
在相同的條件下,降雨歷時(shí)越長,降雨量越大,產(chǎn)生的徑流量越大,山洪災(zāi)害損失也越嚴(yán)重。在溪河洪水、泥石流、滑坡三種山洪災(zāi)害類型中,溪河洪水災(zāi)害和泥石流災(zāi)害特點(diǎn)相似,在下墊面條件滿足的情況下,只要有足夠大的降雨強(qiáng)度和降雨量,溪河洪水和泥石流災(zāi)害就可發(fā)生,降雨歷時(shí)越長,所產(chǎn)生的洪量越大,災(zāi)害損失越大。滑坡災(zāi)害與降雨歷時(shí)的關(guān)系更為明顯,一般而言,滑坡和降雨并不是同時(shí)發(fā)生,滑坡一般滯后于降雨,由溪河洪水誘發(fā)的滑坡災(zāi)害會(huì)隨著降雨歷時(shí)的增長,洪流對(duì)坡腳的掏蝕,滑坡也會(huì)越來越多。
1959年Copula函數(shù)理論提出,Sklar提出可把一個(gè)聯(lián)合分布函數(shù)分解為一個(gè)Copula 函數(shù)和K個(gè)邊緣分布函數(shù),而分解出來的Copula 函數(shù)可用來描述變量之間相關(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 的三類常用函數(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ì)角線上的元素為1的對(duì)稱正定矩陣;|ρ|表示與矩陣ρ相對(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á)為:
(1)不同類型的Copula 函數(shù)比較分析。一般情況下,多元正態(tài)Copula 函數(shù)通常描述變量之間的相關(guān)關(guān)系,但是由于此函數(shù)有對(duì)稱性的特點(diǎn),所以對(duì)于變量的非對(duì)稱相關(guān)關(guān)系很難擬合。
而Gumbel Copula 函數(shù)和Clayton Copula 函數(shù)都具有非對(duì)稱性的特點(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ì)稱性的特點(diǎn),所以很難描述水文變量之間的非對(duì)稱關(guān)系,F(xiàn)rank Copula 函數(shù)只適合描述具有對(duì)稱相關(guān)結(jié)構(gòu)的變量之間的相關(guān)關(guān)系,即各個(gè)水文變量間極大值相關(guān)性與極小值相關(guān)性是對(duì)稱增長的,但是其尾部的分布變量是比較獨(dú)立的,所以Frank Copula 函數(shù)無論在描述上尾部還是下尾部的相關(guān)性中都是不敏感的,故無法描述尾部變化的相關(guān)性。
(2)混合Copula 函數(shù)的構(gòu)造與相關(guān)性分析。Gumbel Copu‐la 函數(shù)、Clayton Copula 函數(shù)、Frank Copula 函數(shù)三類常用的阿基米德函數(shù)能夠捕捉尾部相關(guān)的情形:上尾部的相關(guān)、下尾部的相關(guān)、上尾部和下尾部的對(duì)稱相關(guān)。這些Copula 函數(shù)具有描述各種水文變量模式之間的關(guān)系,特別是尾部相關(guān)關(guān)系,水文系統(tǒng)中各變量之間的關(guān)系是復(fù)雜多變的,不是拘泥于某種特定關(guān)系,它們只能反映水文變量間相關(guān)性的某個(gè)側(cè)面,因此很難用一個(gè)簡單的Copula 函數(shù)來全面的刻畫水文系統(tǒng)中各變量之間的相關(guān)模式,所以一個(gè)更加靈活的Copula 函數(shù)需要構(gòu)造,來描述各種水文變量模式之間的關(guān)系。應(yīng)用不同Copula 函數(shù)的優(yōu)點(diǎn),本文選用Gumbel Copula 函數(shù)、Clayton Copula 函數(shù)、Frank Copula 函數(shù)的線性組合來構(gòu)造混合的Copula 函數(shù),可以更加靈活的刻畫水文系統(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ù)向量(α,λ,θ)用來描述變量間相關(guān)的程度,變量間相關(guān)的模式由線性權(quán)重參數(shù)向量(wG,wF,wcl)表示。
(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)則來計(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為樣本中滿足條件x≤xk的個(gè)數(shù),即滿足:x1≤x1k,x2≤x2k,x3≤x3k。
設(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)合分布概率,本文主要考慮超過概率和條件概率,即當(dāng)總雨量X超過某特定值的條件下,峰值雨強(qiáng)超過某特定值時(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ì)于超過概率,重現(xiàn)期計(jì)算公式如下式所示。
分別采用Gumbel、Clayton、Frank 及混合Copula(M-Copula)函數(shù)構(gòu)建12 h 場次降雨量和峰值雨強(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ù)均通過K-S 檢驗(yàn),能較好的擬合臨界雨強(qiáng)的邊緣分布,根據(jù)均方根誤差RSME最小準(zhǔn)則,選取M-Copu‐la函數(shù)為連接函數(shù)來擬合降雨量與峰值雨強(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
綜合考慮防災(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)的暴雨洪水頻率差異過大。
通常情況下,由于臨界雨量是從成災(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í)間及洪水過程線綜合進(jìn)行“折減”處理。
以軍馬河鄉(xiāng)白果村回龍灣組為例,在不同的土壤含水量的情況下,利用風(fēng)險(xiǎn)概率雨型分別計(jì)算1、3、6、24 h 臨界雨量,其臨界雨量預(yù)警閾值見表2 所示;水位預(yù)警閾值見圖1 所示;準(zhǔn)備轉(zhuǎn)移預(yù)警和立即轉(zhuǎn)移預(yù)警雨量臨界線圖見圖2所示。
圖1 水位預(yù)警閾值圖Fig.1 Threshold map of water level warning
圖2 準(zhǔn)備轉(zhuǎn)移預(yù)警和立即轉(zhuǎn)移預(yù)警雨量臨界線圖Fig.2 Threshold maps for the preparation and immediate diversion of early warning rainfall
表2 臨界雨量預(yù)警閾值表Tab.2 Table of critical rainfall warning thresholds
對(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
通過對(duì)西峽縣的山洪災(zāi)害預(yù)警分析可知:混合Copula 函數(shù)可以更加靈活的刻畫水文系統(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ù)來擬合臨界雨量聯(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ù)警等均具有重要作用。