孫 強，顧艷玲，張 磊，姓海濤，張 敏

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

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

1 譜風險測度模型

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 時，Gξ,δ(x)分別對應于 ParetoⅡ分布，指數(shù)分布和普通帕累托分布。

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

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

(2)

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

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

(3)

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

(4)

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

1.2 風險厭惡函數(shù)

1.2.1 譜風險測度數(shù)學描述

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

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

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

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

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

φ(p)(p)dp

(5)

樣本數(shù)據(jù)分布函數(shù)為離散型時，則譜風險定義為：

φi

(6)

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

1.2.2 指數(shù)風險厭惡函數(shù)

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

(1)指數(shù)風險厭惡函數(shù)。

(7)

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

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

(2)冪風險厭惡函數(shù).

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

(8)

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

(3)雙曲型風險厭惡函數(shù)。

(9)

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

2 極值譜風險度量模型

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

(10)

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

3 建模流程

建模流程見圖2。

圖2 建模流程Fig.2 Modeling process

4 應用實例

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

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

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

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

對于數(shù)據(jù)分布的判斷，可以從數(shù)據(jù)統(tǒng)計的峰度、偏度，還有Q-Q圖中進行說明。從表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)計表Tab.1 Descriptive statistics of daily variation of seepage flow

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

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

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

4.3 用譜風險測度計量實證

選取極值譜風險測度方法進行大壩滲流量每日變化量的風險測量時，應注意此方法應假設每日滲流量變化序列的尾部分布服從廣義累托分布。首先，我們要采用滲流量每日變化的正值分布的Hill估計值來確定閾值。

圖5 Hill圖Fig.5 Hill plot

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

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

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

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

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

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

5 結(jié) 語

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

□

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

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

[3] 任 杰， 蘇懷智， 陳 蘭，等.基于POT模型的大壩位移預警指標實時估計[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] 益 智,楊敏敏. 基于極值譜風險測度的金融市場風險度量[J].商業(yè)經(jīng)濟與管理，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] 呂世超,田 凱,楊永愉.基于譜風險度量的風險譜函數(shù)的研究[J].北京化工大學學報(自然科學版)，2011,(6).

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