极值分析介绍（5）——降水IDF分析

发表于 分类于 阅读次数： Disqus：

接着前面对DDF分析的余温，这一部分简单介绍与其十分相像的IDF分析。

引言

降水强度是指单位时间的降水量。通常以毫米/日、毫米/小时、毫米/10分钟为单位，还可根据应用部门的专门需要而定。在桥梁、涵洞等工程中所取时段较短(几分钟至几小时)，在农业上有的选取时段较长(旬、月或作物生长发育期)。

前面介绍了DDF曲线，与之类似，强度-持续时间-频率就构成了IDF曲线。降雨IDF曲线信息描述了各种持续时间和强度的极端降雨事件的频率。从定义我们就可以发现，DDF与IDF曲线存在非常密切的关系，因为Intensity=Depth/DurationIntensity=Depth/Duration，实际上可以认为DDF实际上就是IDF一种不同的表达形式。

载入IDF数据

查阅文献发现，IDF曲线更为常用，与DDF类似，IDF值也可以通过分析不同持续时间的年度最大降水量获得。但是，网上很那查到高分辨率（例如5分钟、15分钟）的将于数据，所以我只好从NOAA的网站上直接找了IDF的数据来做示例分析。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
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 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
data = pd.read_csv('IDF_Miami.csv', sep='\t', header=0)
data.rename(columns={'by duration for AEP:': 'duration'}, inplace=True)
data = data.iloc[2:-7]
data = data.iloc[:,:-2]

data['duration'] = [15/60,0.5, 1, 2, 3, 6, 12, 24, 48, 72]
data

duration '1/2 '1/5 '1/10 '1/25 '1/50 '1/100 '1/200
2 0.25 110 141 165 196 220 244 269
3 0.50 86 111 130 155 174 194 214
4 1.00 58 74 88 107 122 139 156
5 2.00 36 47 55 68 79 90 103
6 3.00 27 35 42 52 61 70 81
7 6.00 16 21 25 32 38 45 52
8 12.00 9 12 15 20 23 27 32
9 24.00 5 7 9 12 14 16 18
10 48.00 3 4 5 7 8 9 10
11 72.00 2 3 4 5 5 6 7 
1
2
3
(ggplot(data.melt(id_vars='duration')) +
 aes(x='duration', y='value', color='variable') +
 geom_line())

output_4_0

拟合IDF的经验曲线

虽然我找到的IDF数据已经比较细致了，但是仍然有一些数据是缺失的，比如持续时间为45分钟的降雨强度在不同回归期的数据。为了得到这些数据，从而使得数据更加平滑，可以对IDF数据进行插值，那么在一定的假设下，广义的IDF关系可以表示为：

I=a⋅Tk(D+b)c,(1)I = \frac{a \cdot T^k}{(D+b)^c},\tag{1}

其中

  • I 为降水强度（mm/hour）；
  • D 为降水持续时间（hour）；
  • T 为回归年
  • a, b, c和k均为系数。

Python拟合函数库

DDF中，我们最后化简为线性方程组的求解问题，在IDF的拟合中，我们直接采用Python的非线性最小二成曲线拟合库（The Non-Linear Least-Squares Minimization and Curve-Fitting for Python， LMFIT)来拟合上述IDF曲线。

1
2
from lmfit import Model
from lmfit import Parameter

根据公式(1)，定义我们需要拟合的IDF函数：

1
2
def regidf(D, T, a, b, c, k):
    return a*(T**k)/((D+b)**c)

在定义的IDF函数的基础上，利用lmfit包创建一个回归模型，其中自变量为持续时间（D）和回归年（T），因变量为降水强度（I）：

1
model = Model(regidf, independent_vars=['D', 'T'])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import re

idfparams = OrderedDict()

T = np.array([2, 5, 10, 25, 50, 100, 200])*1.
for i in range(1, len(data.columns)):
    datax = data.iloc[:,i] * 1.0
    result = model.fit(datax,
                      D = data.duration,
                      T = T[i-1],
                      a = datax.mean(),
                      b = Parameter(value=1, min=0),
                      c = Parameter(value=1, min=0),
                      k = Parameter(value=1, min=0))

    idfparams[T[i-1]] = [result.params['a'].value,result.params['b'].value,
                       result.params['c'].value,result.params['k'].value]
fnl = pd.DataFrame.from_dict(idfparams, orient='index')
fnl.columns=["a","b","c","k"]
fnl
a b c k
2.0 38.257228 0.467853 0.890587 1.104794
5.0 35.178514 0.446387 0.867877 0.671393
10.0 24.728260 0.432919 0.846448 0.686858
25.0 37.761637 0.410247 0.801326 0.410080
50.0 148.194826 0.418142 0.784277 0.021104
100.0 154.817823 0.438750 0.772256 0.036860
200.0 183.639597 0.449728 0.753146 0.021370

绘制IDF曲线

与之前的工作类似，设定持续时间的采样区间为0.25小时(15分钟)~72小时，采样数为100：

1
samdurations = np.linspace(0.25,72.0,1000)

采用前面估计的回归系数，计算IDF对于采样样本的估计值：

1
2
3
4
5
6
7
8
9
fnlt = fnl.T
ddfvalues = OrderedDict()
for t in T:     
    ddfvalues[t] = regidf(samdurations,T=t, a= fnlt[t]["a"],
                             b=fnlt[t]["b"],c=fnlt[t]["c"],
                             k=fnlt[t]["k"],)

inh = pd.DataFrame.from_dict(ddfvalues, orient='index').T
inh['duration'] = samdurations
1
2
data.columns=['duration', 2, 5, 10, 25, 50, 100, 200]
df = data.melt(id_vars='duration')
1
2
3
4
5
6
7
8
9
10
11
est = inh.melt(id_vars='duration')

(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 = 25","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_19_0

小结

在DDF的基础上，我们进一步讨论了IDF的拟合问题，并采用LMFIT库完成了拟合。需要注意的是，对于LMFIT而言，参数的初始值选择非常重要，尤其是当方程存在欠定的情况下，IDF可能会收敛到不同的参数上。

此外，当Return Period到500年以上时，出现了无法求解的问题，说明此时的数据可能已经不适合用IDF模型假设的形式进行描述。尤其是，我们观察到，即使是对数表示的曲线，仍然存在比较明显的非线性，或者说Duration在较大和较小时服从不同的关系。这也提醒我们在拟合曲线时注意适用范围。

I’ve seen things you people wouldn’t believe.
Attack ships on fire off the shoulder of Orion.
I watched C-beams glitter in the dark near the Tannhäuser Gate.
All those moments will be lost in time, like tears in rain.
Time to die.

–Tears in Rain/C-Beams Speech (from Balde Runner)

blade_runner_2049