孫 強(qiáng)，顧艷玲，張 磊，姓海濤，張 敏

(1.河海大學(xué) 水利水電學(xué)院，南京 210098；2.中國水利水電科學(xué)研究院，北京 100038；3.中水北方勘測設(shè)計(jì)研究有限責(zé)任公司，天津 300222；4.南京市溧水區(qū)水務(wù)局，南京 211200)

水工建筑物的安全一直以來都是水利工程中十分關(guān)注的問題，所以大壩觀測量閾值的研究有很大的必要性。有關(guān)大壩的閾值研究成果已有諸多報(bào)道，如今大壩閾值一般由置信區(qū)間法、典型小概率法計(jì)算而獲得[1]；谷艷昌等[2]基于時(shí)效分量突變模型，結(jié)合遺傳尋優(yōu)算法，提出了水庫大壩閾值的確定方法；任杰等[3]提出基于POT模型的大壩預(yù)警指標(biāo)實(shí)時(shí)估計(jì)方法，說明基于POT模型來估計(jì)大壩預(yù)警指標(biāo)這種方法更加安全合理。由此可見對于閾值的研究方法都是以一定量的數(shù)學(xué)方法為基礎(chǔ)，通過實(shí)測資料得出時(shí)效模型，再結(jié)合相關(guān)計(jì)算方法，從而確定大壩的閾值，而現(xiàn)實(shí)中很少有針對滲流量分布提出閾值，所以本文提出用POT模型來刻畫大壩每日滲流值變化量的分布，并加入譜風(fēng)險(xiǎn)測度模型來計(jì)算出大壩的閾值指標(biāo)。本文中引入考慮專家風(fēng)險(xiǎn)態(tài)度的譜風(fēng)險(xiǎn)函數(shù)，來刻畫專家對待監(jiān)測數(shù)據(jù)的態(tài)度，從而從主觀和客觀兩個(gè)角度構(gòu)建極值譜風(fēng)險(xiǎn)度量模型來計(jì)算大壩運(yùn)行閾值指標(biāo)。

1 譜風(fēng)險(xiǎn)測度模型

1.1 極值理論的POT模型

在極值理論中，POT(Peaks Over Threshold)模型對于極端數(shù)據(jù)的處理是非常有效的。POT 模型主要基于廣義帕累托分布(Generalized Pareto Distribution，簡稱 GPD)，廣義帕累托分布由Pickands[4]，Balkema和de Haan[5]引入，公式如下：

(1)

式中：ξ為形狀參數(shù)；δ>0為尺度參數(shù)。ξ<0，ξ=0和ξ>0 時(shí)，Gξ,δ(x)分別對應(yīng)于 ParetoⅡ分布，指數(shù)分布和普通帕累托分布。

根據(jù)PDBH定理[6]，對足夠大的閾值μ，有尾部估計(jì)：

F(x)=[1-Fn(μ)]Gξ,δ(x)+Fn(μ)

(2)

式中：Fn(μ)是在閾值μ處的經(jīng)驗(yàn)分布函數(shù)；Gξ,δ(x)為廣義帕累托分布，其形狀參數(shù)ξ與(1)式相同。

將式子(2)中經(jīng)驗(yàn)分布函數(shù)Fn(μ)用n-Nμ/n估計(jì)并化簡，所以可以確定F(x)：

(3)

在給定的置信水平p的條件下，由(3)式求反函數(shù)得到F-1x(p)：

(4)

關(guān)于 GPD 分布參數(shù)的估計(jì)采用極大似然估計(jì)(Maximum Likelihood Estimate，即 MLE)，Hosking和 Wallis(1987年)證明在ξ≥-1時(shí)，MLE 方法在大樣本條件下比其他估計(jì)方法如概率加權(quán)矩估計(jì)更加有效[6]。

1.2 風(fēng)險(xiǎn)厭惡函數(shù)

1.2.1 譜風(fēng)險(xiǎn)測度數(shù)學(xué)描述

2002年，Acerbi提出了譜風(fēng)險(xiǎn)測度(Spectral Risk Measure)[7]，它考慮了專家的風(fēng)險(xiǎn)譜或風(fēng)險(xiǎn)厭惡函數(shù)，即對大壩風(fēng)險(xiǎn)的厭惡性。Acerbi提出譜風(fēng)險(xiǎn)測度概念時(shí)，也對風(fēng)險(xiǎn)厭惡函數(shù)提出了要求：

(1)φ(p) 是(0,1]上可積的實(shí)函數(shù)；

(2)規(guī)范性：φ(p)≥0,?p∈(0,1]；

(3)單調(diào)性：若p1≤p2?φ(p1)≤φ(p2)；

其中，φ(p)表示風(fēng)險(xiǎn)厭惡函數(shù)，代表對風(fēng)險(xiǎn)的厭惡程度。當(dāng)φ(p)滿足以上三種性質(zhì)時(shí)，φ(p)稱為可容許風(fēng)險(xiǎn)譜。Xp表示樣本數(shù)據(jù)分布函數(shù)的逆函數(shù)，可用F-1x(p)來表示，那么譜風(fēng)險(xiǎn)測度的計(jì)算公式也可以表示為：

φ(p)(p)dp

(5)

樣本數(shù)據(jù)分布函數(shù)為離散型時(shí)，則譜風(fēng)險(xiǎn)定義為：

φi

(6)

式中：φi表示風(fēng)險(xiǎn)厭惡函數(shù)，代表專家對風(fēng)險(xiǎn)的厭惡程度；X表示一個(gè)包含n個(gè)變化值的樣本；X[i]表示將樣本數(shù)據(jù)變化值按從大到小排列后排在第i位的變化值。

1.2.2 指數(shù)風(fēng)險(xiǎn)厭惡函數(shù)

常用的風(fēng)險(xiǎn)厭惡函數(shù)主要有以下三種[8]：

(1)指數(shù)風(fēng)險(xiǎn)厭惡函數(shù)。

(7)

式中：γ為風(fēng)險(xiǎn)厭惡參數(shù)，γ∈(0,∞)；p為累積概率。不同風(fēng)險(xiǎn)厭惡參數(shù)下的指數(shù)風(fēng)險(xiǎn)厭惡函數(shù)圖像如圖1。

圖1 不同風(fēng)險(xiǎn)厭惡系數(shù)下的指數(shù)風(fēng)險(xiǎn)函數(shù)fig.1 Index risk function under different risk aversion coefficient

(2)冪風(fēng)險(xiǎn)厭惡函數(shù).

φγ(p)=(1-γ)(1-p)-γ

(8)

式中：γ為風(fēng)險(xiǎn)厭惡參數(shù)，γ∈(0,1)；p為累積概率。

(3)雙曲型風(fēng)險(xiǎn)厭惡函數(shù)。

(9)

式中：γ為風(fēng)險(xiǎn)厭惡參數(shù)，γ∈[0,μ]；1-μ為置信水平，也就是說當(dāng)置信水平取0.99時(shí)，α=0.01；當(dāng)置信水平取0.95時(shí)，α=0.05。

2 極值譜風(fēng)險(xiǎn)度量模型

在風(fēng)險(xiǎn)厭惡函數(shù)選取中，我們選取指數(shù)風(fēng)險(xiǎn)厭惡函數(shù)，結(jié)合式子(4)在此基礎(chǔ)上構(gòu)造的極值譜風(fēng)險(xiǎn)測度模型為：

(10)

式中：μ為閾值；ξ為形狀參數(shù)；δ為尺寸參數(shù)；n為樣本總個(gè)數(shù)；Nμ為超過閾值的樣本個(gè)數(shù)；p為置信水平。

3 建模流程

建模流程見圖2。

圖2 建模流程Fig.2 Modeling process

4 應(yīng)用實(shí)例

4.1 數(shù)據(jù)的采集與分析

本文實(shí)證采用的是某壩的滲流量實(shí)測資料，本文開始時(shí)已把所得到的數(shù)據(jù)整理一遍，異常點(diǎn)已被刪除，并根據(jù)除險(xiǎn)加固的前與后，將滲流量的數(shù)據(jù)分成兩個(gè)時(shí)間段：從1974年4月27日到2009年12月31日和從1974年4月27日到2015年12月31日，第一時(shí)間段滲流觀測量共計(jì)918個(gè)數(shù)據(jù)，第二個(gè)時(shí)間段滲流觀測量共計(jì)1 062個(gè)數(shù)值。

令P為此壩的每日滲流量測值，并對Pt做對數(shù)化處理，以減少數(shù)據(jù)的波動(dòng)性。然后令R=lnpt-lnpt-1，得到：第一時(shí)間段滲流觀測值變化量共計(jì)917個(gè)數(shù)據(jù)，第二個(gè)時(shí)間段滲流觀測值變化量共計(jì)1 061個(gè)，最后得到此壩滲流量初始觀測值的變化量直方圖3和數(shù)據(jù)統(tǒng)計(jì)結(jié)果(表1)。

4.2 數(shù)據(jù)分布的判斷

對于數(shù)據(jù)分布的判斷，可以從數(shù)據(jù)統(tǒng)計(jì)的峰度、偏度，還有Q-Q圖中進(jìn)行說明。從表1可見，R值的偏度分別為0.373和0.432，說明此壩的滲流量測值出現(xiàn)了一定程度的右偏，而且出現(xiàn)較大的峰度，大于正態(tài)分布的峰度3，根據(jù)邊寬江、程波等人[9]對尖峰厚尾問題的研究，因此所得出的R值具有尖峰的特征，不服從正態(tài)分布，因而不能用正態(tài)分布來擬合R值的數(shù)列。

表1 滲流量每日變化值描述性資料統(tǒng)計(jì)表Tab.1 Descriptive statistics of daily variation of seepage flow

圖3 滲流量每日變化值直方圖Fig.3 histogram of daily variation of seepage flow

Q-Q圖是一種散點(diǎn)圖，是由標(biāo)準(zhǔn)正態(tài)分布的分位數(shù)為縱坐標(biāo)，樣本值為橫坐標(biāo)的散點(diǎn)圖，Q-Q圖利用了樣本數(shù)據(jù)的分位數(shù)與標(biāo)準(zhǔn)正態(tài)的分位數(shù)進(jìn)行對比，因此證明樣本數(shù)據(jù)的分布屬于哪種類型的分布。從圖4中可以看到Q-Q圖樣本點(diǎn)上凸，說明了該樣本屬于厚尾分布，所以滲流量測值的變化量分布屬于尖峰厚尾分布，滿足POT模型成立的前提條件。

圖4 Q-Q圖Fig.4 Q-Q Plots

4.3 用譜風(fēng)險(xiǎn)測度計(jì)量實(shí)證

選取極值譜風(fēng)險(xiǎn)測度方法進(jìn)行大壩滲流量每日變化量的風(fēng)險(xiǎn)測量時(shí)，應(yīng)注意此方法應(yīng)假設(shè)每日滲流量變化序列的尾部分布服從廣義累托分布。首先，我們要采用滲流量每日變化的正值分布的Hill估計(jì)值來確定閾值。

圖5 Hill圖Fig.5 Hill plot

由圖5可以看出，大壩滲流量變化量的變化值總數(shù)在兩個(gè)階段分別為916個(gè)和1060個(gè)，分布序列從40開始散點(diǎn)圖趨于穩(wěn)定，選取穩(wěn)定的標(biāo)準(zhǔn)是兩個(gè)相鄰變化量的差值小于0.5%，因此選取的閾值u為0.694 0和0.693 2，大于閾值的數(shù)據(jù)Nu分別有30個(gè)和33個(gè)，接下來根據(jù)以上粗略估計(jì)的閾值，依據(jù)超額X-u數(shù)據(jù)采用極大似然方法對廣義帕累托函數(shù)的參數(shù)進(jìn)行估計(jì)，得到參數(shù)δ、ξ的最大估計(jì)值為：

最后，根據(jù)上表的估計(jì)值，且選取風(fēng)險(xiǎn)厭惡參數(shù)γ為0.005，然后計(jì)算出該壩滲流變化量的閾值，如表3。

從表3可以看出，在1974-2009年極值譜風(fēng)險(xiǎn)模型算出來的數(shù)值比無風(fēng)險(xiǎn)厭惡函數(shù)的POT模型較小，但是相差不大，最大相差0.06%；但是在1974-2015年期間，兩者相差0.36%，說明在數(shù)據(jù)量比較大時(shí)，極值譜風(fēng)險(xiǎn)模型算出來的風(fēng)險(xiǎn)度量值更加能體現(xiàn)出此壩為安全狀態(tài)。

在實(shí)際上，由于此壩2010年進(jìn)行除險(xiǎn)加固，此壩正在良好的運(yùn)行當(dāng)中，說明2015年大壩的安全性與2009年之前相比更高。通過極度譜風(fēng)險(xiǎn)模型計(jì)算2015年時(shí)的閾值計(jì)算和2009年相比，極值譜風(fēng)險(xiǎn)模型更能說明大壩運(yùn)行狀況從較高風(fēng)險(xiǎn)降到較安全的狀態(tài)。

表2 參數(shù)估計(jì)值Tab.2 Parameter estimation

表3 結(jié)果對比Tab.3 Comparison results

5 結(jié) 語

POT模型通過合理地設(shè)定閾值，更全面精準(zhǔn)的考慮了所有變化較大的測值，因此更能客觀地表現(xiàn)出工程實(shí)際測值的分布情況；風(fēng)險(xiǎn)厭惡函數(shù)能考慮每個(gè)變化值對于大壩風(fēng)險(xiǎn)的重要性，因此能體現(xiàn)大壩的小概率風(fēng)險(xiǎn)的情況。雖然通過工程實(shí)例驗(yàn)證，依據(jù)極值譜風(fēng)險(xiǎn)測度模型測量出的閾值數(shù)值較小，但是極值譜風(fēng)險(xiǎn)測度具有更為嚴(yán)謹(jǐn)?shù)摹⑼晟频挠?jì)算過程以及全面的理論基礎(chǔ)，而且在現(xiàn)實(shí)中，大壩監(jiān)測數(shù)據(jù)變化值的分布總是存在尖峰厚尾特性，并可以為大壩安全評估提供了依據(jù)，所以極值譜風(fēng)險(xiǎn)測度的實(shí)用性比較高，是一個(gè)比較準(zhǔn)確的風(fēng)險(xiǎn)測度方法。

□

[1] 吳中如.水工建筑物安全監(jiān)控理論及其應(yīng)用[M].北京：高等教育出版社，2003.

[2] 谷艷昌,王士軍.水庫大壩結(jié)構(gòu)失穩(wěn)突發(fā)事件預(yù)警閾值研究[J].水利學(xué)報(bào)，2009,12.

[3] 任 杰， 蘇懷智， 陳 蘭，等.基于POT模型的大壩位移預(yù)警指標(biāo)實(shí)時(shí)估計(jì)[J].水力發(fā)電，2016,(4).

[4] Pickands J.Statistical inference using extreme order statistics [J].Ann Statist，1975，(3):119-131.

[5] Balkema A A，de Haan L.Residual life time at great age [J].Ann Probab,1974，2:792-804.

[6] 益 智,楊敏敏. 基于極值譜風(fēng)險(xiǎn)測度的金融市場風(fēng)險(xiǎn)度量[J].商業(yè)經(jīng)濟(jì)與管理，2009,(8).

[7] HOSKING J R M，WALLIS J R. Parameter and Quantile Estimation for the Generalized Pareto Distribution [J]. Technometrics，1987，29(3)：339-349.

[8] ACERBI C.Spectral Measures of Risk:A Coherent Representation of Subjective Risk Aversion[J].Journal of Banking and Finance，2002，26(7)：1 505-1 518.

[9] 呂世超,田 凱,楊永愉.基于譜風(fēng)險(xiǎn)度量的風(fēng)險(xiǎn)譜函數(shù)的研究[J].北京化工大學(xué)學(xué)報(bào)(自然科學(xué)版)，2011,(6).

[10] 邊寬江，程 波，王蕾蕾. 收益分布尖峰厚尾問題的統(tǒng)計(jì)檢驗(yàn)[J].統(tǒng)計(jì)與決策，2009,9(7).