基于極值譜風險的大壩滲流量閾值問題的研究

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

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

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

1 譜風險測度模型

1.1 極值理論的POT模型

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

1.2 風險厭惡函數

1.2.1 譜風險測度數學描述

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

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

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

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

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

φ(p)(p)dp

(5)

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

φi

(6)

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

1.2.2 指數風險厭惡函數

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

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

(7)

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

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

(2)冪風險厭惡函數.

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

(8)

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

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

(9)

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

2 極值譜風險度量模型

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

(10)

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

3 建模流程

建模流程見圖2。

圖2 建模流程Fig.2 Modeling process

4 應用實例

4.1 數據的采集與分析

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

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

4.2 數據分布的判斷

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

表1 滲流量每日變化值描述性資料統計表Tab.1 Descriptive statistics of daily variation of seepage flow

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

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

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

4.3 用譜風險測度計量實證

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

圖5 Hill圖Fig.5 Hill plot

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

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

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

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

表2 參數估計值Tab.2 Parameter estimation

表3 結果對比Tab.3 Comparison results

5 結 語

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

□

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

[2] 谷艷昌,王士軍.水庫大壩結構失穩突發事件預警閾值研究[J].水利學報，2009,12.

[3] 任 杰， 蘇懷智， 陳 蘭，等.基于POT模型的大壩位移預警指標實時估計[J].水力發電，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].商業經濟與管理，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] 呂世超,田 凱,楊永愉.基于譜風險度量的風險譜函數的研究[J].北京化工大學學報(自然科學版)，2011,(6).

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