极值分析介绍(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()计算出不同持续时间的降雨量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
1
2
3
4
5
6
7
8
9
10

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

output_4_1.png

<ggplot: (-9223371877569100372)>

针对每个duration进行极值分析

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

1
2
from scipy.stats import genextreme, gumbel_r
from statsmodels.distributions.empirical_distribution import ECDF
1
2
3
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。

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

对于每个(回归期,持续时间)的组合,计算极值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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') )

output_12_0.png

<ggplot: (-9223371877568797608)>

拟合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)R(T,D)为持续时间为D、回归期为T的降雨量,则R(2,D)R(2,D)为持续时间为D的降雨量的总体中值,因为回归期为2年的降雨量对应的累积概率F为0.5,选择R(2,D)R(2,D)为指标降雨量,作为DDF模型的缩放系数。假设在纵向上服从log-logistic分布,那么增长曲线的形式为

R(T,D)R(2,D)=(T1)c,T=1(1F),T>1,(1)\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)Ds(2)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)Ds(T1)c,T>1(3)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),两边取对数,得到:

logR(T,D)=logR(2,1)+slog(D)+clog(T1)(4)\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 \tag{5}\\ s_D &= e + f\log D \tag{6} \end{aligned}

其中指数项c和d均为D而非T的函数,当D=1时,cD=a=c24c_D=a=c_{24},即持续时间为1天的形状参数(公式(1)中c的对数)。

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

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

其中,当持续时间为12小时时,D=0.5, 以此类推。 指数项c仍然为D的函数,但是s为常数。因此,一旦通过公式(5)(6)确定c_1,即公式(6)c_24的值,则只需要确定s和h即可,反之亦然。

示例分析

最小二乘参数估计

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
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.2264D0.2963(T1)0.1580R(T,D) = 86.2264 \cdot D^{0.2963} \cdot (T-1)^{0.1580}

结果可视化

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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')
)

output_18_0.png

<ggplot: (-9223371877568792162)>

小结

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

  1. 极值分析的时候(拟合极值分布),如果我们对不同的Duration选择分布的话,是选择在所有Duration综合表现最好的同一个极值分布,还是针对每一个分布选择最好的极值分布呢?
  2. 我们在为DDF拟合曲线时,采用的是相同的分布(log-logistic),而这个分布也可以选择其他形式的。类似问题1,如果我们针对不同的Duration选择分布的话,又该如何选择呢? (曲线的拟合可以通过scipy.optimize.curve_fit实现。)
张da统帅 wechat
扫码订阅我的微信公众号