孙 强，顾艳玲，张 磊，姓海涛，张 敏

(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).