极值分析介绍（4）——降水DDF分析

前面主要是以气温作为对象进行极值分析，接下来我们以降水量作为对象，介绍水文学中常用的DDF和IDF曲线分析。

IDF/DDF介绍

强度-时间-频率（intensity-duration-frequency,IDF）或深度持续时间频率（depth-duration-frequency,DDF）曲线, 顾名思义，是将降雨强度/深度，与其持续时间和发生频率联系起来的数学函数。IDF和DDF曲线是水资源管理中最常用的工具之一，通常用于水文学的洪水预报和土木工程中的城市排水设计。此外，处于对降水的时间集中或者时间结构分析的需要，IDF曲线也在水文气象学分析中得到应用。

无论是IDF还是DDF曲线，都来源于平稳假设下的历史降雨记录。 降雨IDF/DDF曲线信息描述了各种持续时间和强度的极端降雨事件的频率。根据观测降雨数据的适用性，我们可以采用理论分布、经验分布等不同的数学表达形式来描述IDF/DDF曲线。对于每一个持续时间（duration,例如5、10、60、180……分钟）,可以设置其ECDF、经验分布函数以及确定的频率或者回归期。

因此，经验的IDF/DDF曲线实际上是发生频率相同但是持续时间和强度/深度不同的点的集合。类似地，理论或者半经验的IDF/DDF曲线是物理上合理（具有实际的物理意义）但是其参数是通过经验拟合估计的数据表达式。

因为政府关门的原因，NOAA的网站暂时上不去,看了一下FTP上的原始数据，处理起来需要相当的时间。所以，只好先用示例数据凑活一下。

首先，导入逐小时的降水量数据，如果是分钟级的降雨量可以通过numpy.rolling()计算出不同持续时间的降雨量。

import pandas as pd
import numpy as np
import matplotlib.pylab as plt
from matplotlib.legend_handler import HandlerLine2D
from matplotlib.pylab import rcParams
from collections import OrderedDict
rcParams['figure.figsize'] = 15, 9

import seaborn as sns

import lmoments3
from lmoments3 import distr

from plotnine.ggplot import *
from plotnine.qplot import *
from plotnine.geoms import *
from plotnine.coords import *
from plotnine.labels import *
from plotnine.facets import *
from plotnine.scales import *
from plotnine.themes import *

%matplotlib inline

data = pd.read_csv('precipitation.csv')
data.head()


df = data.melt(id_vars=['anno'])
data.set_index('anno', inplace=True)
(qplot(x='anno', y='value', data=df, color='variable') +
 geom_line() + ggtitle('Annual Maximum PRCP(mm)') +
 labs(x = 'Year',y='Precipitation(mm)', color='duration'))

C:\ProgramData\Anaconda3\lib\site-packages\plotnine\layer.py:449: UserWarning: geom_point : Removed 20 rows containing missing values.
  self.data = self.geom.handle_na(self.data)
C:\ProgramData\Anaconda3\lib\site-packages\plotnine\geoms\geom_path.py:74: UserWarning: geom_path: Removed 1 rows containing missing values.
  warn(msg.format(n1-n2))

针对每个duration进行极值分析

采用前面介绍的方法，在这里我们采用Gumbel(GUM)分布进行分析

from scipy.stats import genextreme, gumbel_r
from statsmodels.distributions.empirical_distribution import ECDF

gumfit = OrderedDict()
for name in data.columns.values.tolist():
    gumfit[name] = distr.gum.lmom_fit(data[name].dropna())

对于不同的回归期，分别拟合DDF曲线

首先选定一些回归期，例如2,5,10,20,50,100,200。

T = np.array([2, 5, 10, 20, 50, 100, 200])
P = 1 - 1.0/T

对于每个（回归期，持续时间）的组合，计算极值

samprobs = OrderedDict()

for name in gumfit.keys():
    fitted_gum = distr.gum(**gumfit[name])
    samprobs[name] = fitted_gum.ppf(P)    

pts = pd.DataFrame.from_dict(samprobs, orient='index')
pts.columns=[2, 5, 10, 20, 50, 100, 200]
pts['duration'] = [1,3,6,12,24]
pts.head()
df = pts.melt(id_vars='duration')

(ggplot(df) +
 aes(x='duration', y='value', color='variable') +
 geom_line() + geom_point() +
 labs(x = 'Duration',y='Precipitation(mm)', color='Return Period') )

拟合DDF曲线

从上面的图上我们可以看到，DDF是一组Duration-PCRP的曲线，因为有三个变量，如果数据非常全面的话，展现的应当是三维空间里的曲面。在这样的情况下，去拟合DDF曲线，也就是估计DDF曲线的参数需要考虑“横向”和“纵向”两个方面的因素。

• 一方面是，考虑固定回归期（Return Period相同）情况下，PCRP随Duration变化的关系，我称之为“横向”（比较不同的Duration）；
• 另一方面是，考虑固定持续时间（Duration相同）情况下，PCRP随Return Period变化的关系，我称之为“纵向”（比较不同的Return Period，正如前面对最高温度的极值分析所做的那样）。

在考虑“横向”和“纵向”两方面约束的条件下，学者提出一个DDF模型应当包括：

1. 一个指标降雨量（index rainfall value），或者称之为“基准降雨量”更为合适，作为模型的缩放系；
2. 由指标降雨量和乘子表示的增长曲线。
   记$R(T,D)$为持续时间为D、回归期为T的降雨量，则$R(2,D)$为持续时间为D的降雨量的总体中值，因为回归期为2年的降雨量对应的累积概率F为0.5，选择$R(2,D)$为指标降雨量，作为DDF模型的缩放系数。假设在纵向上服从log-logistic分布，那么增长曲线的形式为

$$\frac{R(T,D)}{R(2,D)}=(T-1)^c, T=\frac{1}{(1-F)}, T>1,\tag{1}$$

其中F为累积概率密度函数（CDF）。

增长曲线考虑了纵向的关系，横向上则通过考虑不同Duration的R(2,D)的关系加以表征，具体而言, 一个用来表示R(2,D)随着D变化的合理形式为：

$$R(2,D)=R(2,1)D^s \tag{2}$$

其中，D=1也就是持续时间为24小时（1天）作为单位持续时间（unit duration），之所以这样选择是因D=1对于1d~25d（D>=1）的降水和持续时间小于24小时（D<1）的降水都较为合适。因此，持续时间一天的降雨量中值R(2,1)对于DDF模型而言至关重要。

综合公式(1)和（2）,可以得到完整的DDF模型为：

$$R(T,D)=R(2,1)D^s(T-1)^c,T>1 \tag{3}$$

其中指数s，即公式（2）中的s，确定得到对于不同持续时间的R(2,D)需要R(2,1)乘以的倍数；指数c，即公式（1）中的c，是log-logistic增长曲线的形状参数，作为R(2,D)的乘数得到回归期为T时，持续时间D的降雨量R(T,D)。

根据公式（3），两边取对数，得到：

$$\log R(T,D) = \log R(2,1) + s\cdot \log(D) + c\cdot \log(T-1) \tag{4}$$

进一步改写为

• 持续时间1~25天

$$\begin{aligned} c_D &= a + b\log D, 1\leq D\leq 25 \\ s_D &= e + f\log D \tag{5} \end{aligned}$$

其中指数项c和d均为D而非T的函数，当D=1时，$c_D=a=c_{24}$，即持续时间为1天的形状参数（公式（1）中c的对数）。

• 持续时间15分钟~24小时

$$\begin{aligned} c_D &= c_{24} + h \log{D} \\ s_D &= s \tag{6} \end{aligned}$$

其中，当持续时间为12小时时，D=0.5, 以此类推。 指数项c仍然为D的函数，但是s为常数。因此，一旦通过公式（5）确定$c_{D=1}$，即公式（5）$c_{24}$的值，则只需要确定s和h即可，反之亦然。

示例分析

最小二乘参数估计

以上是DDF模型的一些理论基础（具体的数学细节没有展开），那么回到示例。由于我们只有持续时间1、3、6、12和24小时的降水极值，那么应该通过公式（5）估计参数，回到公式（4）可知，方程即为$y= b + a_1 \cdot x_1 + a_2 \cdot x_2$形式，将其增广为线性方程组，采用最小二乘法求解即可，具体代码如下：

to_solve = pts.melt(id_vars=['duration'])
to_solve.head()
to_solve['y'] = np.log(to_solve.value)
to_solve['x1'] = np.log(to_solve.duration/24)
to_solve['x2'] = np.log((to_solve.variable-1.0).tolist())
to_solve['x3'] = 1
X = to_solve[['x1','x2','x3']]
Y = to_solve['y']
solutions = np.dot(np.dot(np.linalg.inv(np.dot(X.T,X)), X.T), Y)
params = {'R21': np.exp(solutions[2]),
          's': solutions[0],
          'c': solutions[1]
}
params

{'R21': 86.22639630504528, 's': 0.29631257152502616, 'c': 0.1580346980212747}

所以，求解的結果为

$$R(T,D) = 86.2264 \cdot D^{0.2963} \cdot (T-1)^{0.1580}$$

结果可视化

根据拟合的表达式作图，虽然公式适用于24小时以内，为了美观期间，做到30小时

samdurations = np.linspace(0.5, 26, 100)
est_prcp = pd.DataFrame()
for t in T:
    est_prcp[str(t)] = (params['R21'] * ((samdurations/24) ** params['s']) *
                      ((t-1.0) ** params['c']))
# est_prcp.columns = ["Tr = 2","Tr = 10","Tr = 5",
#                     "Tr = 20","Tr = 50","Tr = 100", "Tr = 200"]
est_prcp['duration']  = samdurations
est = est_prcp.melt(id_vars='duration')
est['variable'] = est.variable.astype(int)

(ggplot() +
 geom_line(aes(x='duration', y='value', color = pd.Categorical(est.variable)), data=est) +
 geom_point(aes(x='duration', y='value', color = pd.Categorical(df.variable)), data=df) +
 scale_color_discrete(name="Return Period", labels= ["Tr = 2","Tr = 5","Tr = 10",
                                                     "Tr = 20","Tr = 50","Tr = 100", "Tr = 200"])
 + scale_x_log10() + scale_y_log10() +
 labs(x='Duration(Hour) in Log', y = 'Depth of PRCP(mm) in Log')
 )

小结

本节重点介绍了DDF模型的基本形式，以及如何采用极值分析得到DDF曲线。需要注意的是，我们在选择极值分布时，直接采用了Gumbel分布，而没有根据拟合优度度量选择，而在选择DDF曲线的形式时则是通过最小二乘法得到的。那么，有两个问题：

1. 极值分析的时候（拟合极值分布），如果我们对不同的Duration选择分布的话，是选择在所有Duration综合表现最好的同一个极值分布，还是针对每一个分布选择最好的极值分布呢？
2. 我们在为DDF拟合曲线时，采用的是相同的分布（log-logistic），而这个分布也可以选择其他形式的。类似问题1，如果我们针对不同的Duration选择分布的话，又该如何选择呢？ （曲线的拟合可以通过scipy.optimize.curve_fit实现。）