SPEI 算法及代码实现(Java)
前言
关于SPEI
SPEI:Standardized Precipitation Evapotranspiration Index 标准化降水蒸散发指数。该指数是使用降水和气温数据计算得到的表征干湿状态的指数,可以计算多种时间尺度(1,3,12,24月等等),SPEI值越大越湿润,值越小越干旱,在干旱研究中应用广泛。它是从SPI Standardized Precipitation Index 归一化降水指数进一步发展而来,同时纳入蒸散发的影响,在气温有显著变化趋势的地区更加适用,尤其适合长时间尺度的研究。
关于SPEI发展背景及详细描述,可以查阅文献[1]
SPEI算法
这里尽量详细讲解一下SPEI的计算过程。
首先,SPEI指数,使用一种简单的方法计算潜在蒸散发PET,计算公式如下:
P
E
T
=
16
K
(
10
T
I
)
m
PET=16K(\frac{10T}{I})^m
PET=16K(I10T)m
其中,
T
T
T是月均温,单位摄氏度;
I
I
I是热量指数,该指数通过12个月的月度指数相加得到,计算公式如下:
I
=
∑
m
o
n
=
1
12
i
=
∑
m
o
n
=
1
12
(
T
m
o
n
,
a
v
e
5
)
1.514
I=\sum_{mon=1}^{12}i = \sum_{mon=1}^{12}(\frac{T_{mon,ave}}{5})^{1.514}
I=mon=1∑12i=mon=1∑12(5Tmon,ave)1.514
T
m
o
n
,
a
v
e
T_{mon,ave}
Tmon,ave为某月研究时期内多年均值,如一月多年均值为
T
1
,
a
v
e
T_{1,ave}
T1,ave。
m
m
m是一个取决于热量指数的系数:
m
=
6.75
×
1
0
−
7
I
3
−
7.71
×
1
0
−
5
I
2
+
1.79
×
1
0
−
2
I
+
0.492
m=6.75×10^{-7}I^3-7.71×10^{-5}I^2+1.79×10^{-2}I+0.492
m=6.75×10−7I3−7.71×10−5I2+1.79×10−2I+0.492
系数
K
K
K 是通过 该月天数
N
D
M
NDM
NDM 和 最大日照时数
N
N
N 计算得到的:
K
=
(
N
12
)
(
N
D
M
30
)
K=(\frac{N}{12})(\frac{NDM}{30})
K=(12N)(30NDM)
最大日照实时则可以根据其所在地理位置(纬度
ϕ
\phi
ϕ,以弧度表示)和所在一年内的时间得到(地球公转规律)。
N
=
(
24
π
)
ω
s
=
(
24
π
)
a
r
c
c
o
s
(
−
t
a
n
ϕ
t
a
n
δ
)
N=(\frac{24}{\pi})\omega_s=(\frac{24}{\pi})arccos(-tan\phi tan\delta)
N=(π24)ωs=(π24)arccos(−tanϕtanδ)
上式中,
δ
\delta
δ为太阳赤纬(以弧度表示),即是通过一年内所在时间确定:
δ
=
0.4093
s
e
n
(
2
π
J
265
−
1.405
)
\delta=0.4093sen(\frac{2\pi J}{265}-1.405)
δ=0.4093sen(2652πJ−1.405)
其中,
J
J
J是该月平均Julian日,Junlian日就是以一年内日序表示的日期,其中,1月1日为第1 Julian日,12月31日为第365/366 Julian日。
至此就求出来了潜在蒸散发,进一步,根据 降水
P
P
P 和 潜在蒸散发
P
E
T
PET
PET 计算核心变量,即 降水与潜在蒸散发差
D
D
D,某月的
D
D
D计算公式如下:
D
i
=
P
i
−
P
E
T
i
D_i=P_i-PET_i
Di=Pi−PETi
不同时间尺度本质上是
D
i
D_i
Di在时间上的累计,对任意 时间尺度
k
k
k 的计算如下:
当
j
<
k
j<k
j<k
X
i
,
j
k
=
∑
l
=
13
−
k
+
j
12
D
i
−
1
,
l
+
∑
l
=
1
j
D
i
,
l
X_{i,j}^k=\sum_{l=13-k+j}^{12}D_{i-1,l}+\sum_{l=1}^{j}D_{i,l}
Xi,jk=l=13−k+j∑12Di−1,l+l=1∑jDi,l
当
j
≥
k
j\ge k
j≥k
X
i
,
j
k
=
∑
l
=
j
−
k
+
1
j
D
i
,
l
X_{i,j}^k=\sum_{l=j-k+1}^{j}D_{i,l}
Xi,jk=l=j−k+1∑jDi,l
其中,
D
i
,
l
D_{i,l}
Di,l 为第
i
i
i 年第一月。
D
i
,
j
k
D_{i,j}^k
Di,jk 为以
k
k
k为尺度的第
i
i
i年和
j
j
j月的降水蒸散发差。
之后就是对累计变量进行函数拟合,不同于SPI可以用两参数Gamma分布拟合,SPEI变量存在负值情况,需要三参数分布拟合。实验证明。三参数log-logistic分布可以很好拟合变量,该分布函数的表达形式如下:
F
(
x
)
=
[
1
+
(
α
x
−
γ
)
β
]
−
1
F(x)=\left[1+\left(\frac{\alpha}{x-\gamma}\right)^\beta\right]^{-1}
F(x)=[1+(x−γα)β]−1
参数通过下式估算(L-矩估计):
α
=
(
w
0
−
2
w
1
)
β
Γ
(
1
+
1
/
β
)
Γ
(
1
−
1
/
β
)
\alpha=\frac{(w_0-2w_1)\beta}{\Gamma(1+1/\beta)\Gamma(1-1/\beta)}
α=Γ(1+1/β)Γ(1−1/β)(w0−2w1)β
β
=
2
w
1
−
w
0
6
w
1
−
w
0
−
6
w
2
\beta=\frac{2w_1-w_0}{6w_1-w_0-6w_2}
β=6w1−w0−6w22w1−w0
γ
=
w
0
−
α
Γ
(
1
+
1
/
β
)
Γ
(
1
−
1
/
β
)
\gamma=w_0-\alpha\Gamma(1+1/\beta)\Gamma(1-1/\beta)
γ=w0−αΓ(1+1/β)Γ(1−1/β)
其中
w
s
w_s
ws通过下式计算:
w
s
=
1
M
∑
i
=
1
M
(
1
−
F
i
)
s
D
i
w_s=\frac{1}{M}\sum_{i=1}^{M}(1-F_i)^sD_i
ws=M1i=1∑M(1−Fi)sDi
F
i
=
i
−
0.35
M
F_i=\frac{i-0.35}{M}
Fi=Mi−0.35
标准化的分布函数值(概率)计算如下:
p
=
1
−
F
(
x
)
p=1-F(x)
p=1−F(x)
之后即可计算SPEI指数,当
p
≤
0.5
p\le0.5
p≤0.5时
w
=
−
2
l
n
(
p
)
w=\sqrt{-2ln(p)}
w=−2ln(p)
S
P
E
I
=
w
−
c
0
+
c
1
w
+
c
2
w
2
1
+
d
1
w
+
d
2
w
2
+
d
3
w
3
SPEI=w-\frac{c_0+c_1w+c_2w^2}{1+d_1w+d_2w^2+d_3w^3}
SPEI=w−1+d1w+d2w2+d3w3c0+c1w+c2w2
各常量系数为:c0=2.515517, c1=0.802853, c2=0.010328, d1=1.432788, d2=0.189269, and d3=0.001308。当
p
>
0.5
p\gt0.5
p>0.5时
w
=
−
2
l
n
(
1
−
p
)
w=\sqrt{-2ln(1-p)}
w=−2ln(1−p)
S
P
E
I
=
−
(
w
−
c
0
+
c
1
w
+
c
2
w
2
1
+
d
1
w
+
d
2
w
2
+
d
3
w
3
)
SPEI=-(w-\frac{c_0+c_1w+c_2w^2}{1+d_1w+d_2w^2+d_3w^3})
SPEI=−(w−1+d1w+d2w2+d3w3c0+c1w+c2w2)
代码实现(Java)
由于代码过长,就不在这里粘贴了,有需要可以点击这里自行下载(目前只能计算12个月及以内时间尺度):
本代码定义了一个SPEI计算类SPEII, 放在了SPEII.java,使用如下:
double[][] pre = new double[38][12];
double[][] temp = new double[38][12];
……
SPEII speii = new SPEII(temp, pre, lat,1,1981); // lat:纬度 1:scale 1981:起始年份
double[][] result = speii.getResult();
参考文献
[1]Vicente-Serrano S M , S Beguería, JI López-Moreno. A Multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index[J]. Journal of Climate, 2010, 23(7):1696-1718.