全文下载链接： tecdat.cn/?p=15301
在本文中，结合POT模型的洪水风险评估能够从有限的实测资料中获取更多的洪水风险信息,得到更贴近事实的风险评估结果,能为决策者提供更多的依据,从而使决策结果更加可靠实用。
对于这些同样面临挑战的人，我希望这个博客将有助于简化工作。
案例POT序列在47年的记录期内提供了高于74 m 3 / s 阈值的47个峰值。
我们的目标是将概率模型拟合到这些数据并估算洪水分位数。
我从获取了每次洪水的日期，并将其包含在文件中。有趣的是，最早的洪水流量是1943年，而最后一次是1985年，是43年的记录，而不是47年。这是因为1939年至1943年的洪水都小于74 m 3 / s的阈值。
首先计算这些数据点的绘制位置。
请注意，这是记录的年数，而不是峰值数。
同样，重要的是要认识到，方程式1对POT系列的作用与年度系列不同。让我们看一个显示这种差异的示例。考虑以下情况：我们根据47年的数据分析了POT系列的94个峰。在这种情况下，最小的峰的等级为94。重复间隔为：
这大约是半年或6个月，这似乎是合理的（47年中有94个高峰，因此平均每年有2个高峰，平均相隔约6个月）。
将绘图位置解释为年度超出概率将得出以下结果：
也就是说，概率大于1，这没有意义。因此，我们不能使用绘图位置公式来计算阈值峰值序列中的数据的AEP。取而代之的是，方程式1的逆可以解释为EY，即每年的预期超出次数。
ARR示例将指数分布拟合为概率模型。
为了计算L2，我们使用QJ Wang（Wang，1996）的公式
L2 <- function(q){ q <- sort(q) n <- length(q) 0.5*(1/choose(n,2))*sum((0:(n-1) - (n-1):0)*q)}
qi从最小到最大的顺序是流量（POT）
n是流的数量
L2 = 79.12
指数分布的参数可以用L矩表示。我们使用的是广义帕累托（GP）公式。
对于指数分布：
这些参数估计值的置信区间可以使用bootstrapping计算得出。
Beta的95％置信区间是（37.4，89.4）和（120.6，244.7）。参数之间的相关性约为-0.5。参数的不确定性如图1所示。
param_errors_df %>%ggplot(aes(x = V1, y = V2)) +geom_point(size = 0.1) +scale_x_continuous(name = 'beta') +scale_y_continuous(name = bquote('q'['*'])) +stat_ellipse(colour = 'red') + # 95% 置信区间theme_gray(base_size = 7)
图1：参数的不确定性。椭圆显示置信限度为95％
点击标题查阅往期内容
R语言极值推断：广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法
左右滑动查看更多
01
02
03
04
指数分布将超出概率与流的大小相关。在这种情况下，在任何POT事件中，峰值流量超过某个值的概率
为：
这是针对超额概率的。在水文学中，我们通常使用超出概率（洪水大于特定值的概率），因此所需方程式为一个减去所示方程式。
通过将每年超过阈值的洪峰平均数乘以POT概率，我们可以将POT概率转换为每年的预期超标次数。
74 m 3 / s阈值，POT系列中有47个值，并且有47年的数据，因此每年的平均峰值数为1。因此，EY可以表示为：
其中，q是每年POT洪水的平均数量，
w也可以与EY以下内容相关 ：
我们还可以EY通过以下公式与年度超出概率相关：
因此，通过以下等式，洪水幅度可以与AEP相关，而AEP可以与洪水幅度相关。
这些方程式可用于估计标准EY值的分位数。使用bootstrap自举法估计了置信区间（95％）（表1）。
res
| EY | AEP | AEP (1 IN X) | ARI | Q | LWR | UPR |
|---|---|---|---|---|---|---|
| 1.00 | 0.63 | 1.6 | 1.0 | 68 | 39 | 90 |
| 0.69 | 0.50 | 2.0 | 1.4 | 127 | 106 | 150 |
| 0.50 | 0.39 | 2.5 | 2.0 | 178 | 151 | 215 |
| 0.22 | 0.20 | 5.1 | 4.5 | 308 | 253 | 404 |
| 0.20 | 0.18 | 5.5 | 5.0 | 323 | 267 | 427 |
| 0.11 | 0.10 | 9.6 | 9.1 | 417 | 335 | 565 |
| 0.05 | 0.05 | 20.5 | 20.0 | 542 | 434 | 754 |
| 0.02 | 0.02 | 50.5 | 50.0 | 687 | 548 | 965 |
| 0.01 | 0.01 | 100.5 | 100.0 | 797 | 627 | 1167 |
现在绘制数据和拟合度（图2）。x值是根据等式1的逆计算的EY；y值是流量。拟合基于等式6。使用bootstrap自举法计算分位数的置信区间。
图2：河流的部分序列显示契合度和置信区间
我个人更希望该图向右增加，这通常是洪水频率曲线的绘制方式。这仅涉及使用ARI作为纵坐标（图3）。
图3：河流部分序列显示契合度和置信区间
点击文末 “阅读原文”
获取全文完整代码数据资料。
本文选自《R语言POT超阈值模型在洪水风险频率极值分析中的应用研究》。
点击标题查阅往期内容
【视频】R语言极值理论EVT：基于GPD模型的火灾损失分布分析|数据分享
POT超阈值模型和极值理论EVT分析
R语言极值推断：广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法
R语言极值理论EVT：基于GPD模型的火灾损失分布分析
R语言有极值（EVT）依赖结构的马尔可夫链(MC)对洪水极值分析
R语言POT超阈值模型和极值理论EVT分析
R语言混合正态分布极大似然估计和EM算法
R语言多项式线性模型：最大似然估计二次曲线
R语言Wald检验 vs 似然比检验
R语言GARCH-DCC模型和DCC（MVT）建模估计
R语言非参数方法：使用核回归平滑估计和K-NN(K近邻算法)分类预测心脏病数据
matlab实现MCMC的马尔可夫转换ARMA - GARCH模型估计
R语言基于Bootstrap的线性回归预测置信区间估计方法
R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型
Matlab马尔可夫链蒙特卡罗法（MCMC）估计随机波动率（SV，Stochastic Volatility） 模型
Matlab马尔可夫区制转换动态回归模型估计GDP增长率R语言极值推断：广义帕累托分布GPD使用极大似然估计、轮廓似然估计、Delta法
