在自然界，干旱一般有两种类型，一类是由气候、海陆分布、地形等相对稳定的因素在某个相对固定的地区常年形成的水分短缺现象，这类干旱也可称为干旱气候；另一类是各种气象因子（如降水、气温等）的年际变化形成的随机性异常水分短缺现象。

干旱指标计算

干旱分类

1 标准降水指数（Standardized Precipitation Index, SPI）

1993-The relationship of drought frequency and duration to time scale-Thomas B. McKee
数据资料： 月尺度降水资料
优点： 计算简单，且在时空上具有可比性；可计算不同尺度SPI用于表征各类不同可用水资源的干旱；可描绘干湿状况。
缺点： 不能满足任意尺度的应用，只能用在月尺度及以上；不考虑前期降雨的影响；忽略气温变化对干旱影响，在研究全球变暖背景下的干旱问题时适用性较差。
计算步骤：

具体来讲：

1.1 MATLAB代码实现

主函数：调用SPI函数

% scale 时间尺度
% nseas 季节数量
scale=[1,3,6,12];        
nseas=12;

% 一个月尺度SPI，其他时间尺度与此类似
Z = SPI(Data,scale(1),nseas);

输入的Data格式如下：为60年每个月的月均降雨量

计算SPI函数：

function [Z]=SPI(Data,scale,nseas)
%Standardized Precipitation Index （标准化降水指数）
% Input Data（输入数据）
% Data : Monthly Data vector not matrix (monthly or seasonal precipitation)
% 数据：月度数据向量而非矩阵(月度或季节性降水)
% scale （事件尺度）: 1,3,6,12,48  月尺度‘
% nseas : number of season (monthly=12) 季节数量 nseas=12
% Example
% Z=SPI(gamrnd(1,1,1000,1),3,12); 3-monthly scale, 
% Notice that  the rest of the months of the fist year are removed.
% eg. if scale =3, fist year data 3-12 SPI values are not estimated.

%if row vector then make coloumn vector 行/列向量转换
%if (sz==1) Data(:,1)=Data;end % 确保数据为列向量

erase_yr=ceil(scale/12);         % 向上取整 (若为1,3,6,12，erase_yr为1)

% Data setting to scaled dataset 数据设置为缩放数据集
A1=[];          % 初始化
for is=1:scale
    A1=[A1,Data(is:length(Data)-scale+is)];    % 按时间尺度列出数据
end
XS=sum(A1,2);        % 对X的每行分别求和

if (scale>1)        % n个月尺度则为自当前月延续n-1个月（年份取整数）
    XS(1:nseas*erase_yr-scale+1)=[];   
end
% q 降水量为0的概率
% Gam_xs 累积概率
% PARMHAT=gamfit(X);   伽玛分布的参数估计
% 返回拟合X中数据的gamma分布参数的最大似然估计
% PARMHAT（1）和PARMHAT（2）分别是形状和比例参数a和B的估计
% P=gamcdf（X，A，B）返回gamma累积分布函数
% X= NORMINV(P,MU,SIGMA)   其中，P取概率 MU取均值 SIGMA取方差
% 返回值为a    对连续函数，所有小于等于a的值，其出现概率的和为F(a)=P(x<a)
for is=1:nseas
    tind=is:nseas:length(XS);
    Xn=XS(tind);                   % 对应序数
    [zeroa]=find(Xn==0);    % Xn为0的序数
    Xn_nozero=Xn;
    Xn_nozero(zeroa)=[];
    q=length(zeroa)/length(Xn);
    parm=gamfit(Xn_nozero);
    Gam_xs=q+(1-q)*gamcdf(Xn,parm(1),parm(2));
    Z(tind)=norminv(Gam_xs);
end      

1.2 干旱等级说明

SPI实际是多年同一日历日的降水概率，是否可以与重现期等联系起来？

1.3 序列长短的补充说明

SPI在计算一个月时间尺度SPI值时，各月降水分析得到各月SPI，无所争议。但当时间尺度大于1个月时，如三个月时间尺度时，第一个SPI值是前3个月的降水和分析得到的，如此以来，前俩月无法得到相应的SPI3。

SPI函数中下列代码用于将年份取整：如60年数据共720个序列长度，有下列代码时求得708（即删除了前一年的12个数据）；将下列代码除去，得到718个序列长度的SPI3值（无前俩月的SPI3）。可根据需要进行相应修改。

if (scale>1)         % n个月尺度则为自当前月延续n-1个月（年份取整数）
    XS(1:nseas*erase_yr-scale+1)=[];   
end

思考： SPI中每一个月降水是平均日降水或者月降水量和，对结果是否都没影响？最终都是要化为标准正态分布。仔细想了想，月均降水应该更合理，将每一日降水程度都考虑到，这样也能和后面的日尺度SWAP指数联系起来。
but，如果输入为各个月的月均降水，那么在计算3个月及以上时间尺度的SPI时，将各个月的月均降水加起来是否并不合理？所以在使用时还是使用相应尺度上的累计降水？？？
两者的区别应该不大吧。

2 标准降水蒸散发指数（Standardize precipitation evaporation index, SPEI）

2010- A multiscalar Drought Index Sensitive to Global Warming: The Standardized Precipitation Evapotranspiration Index（提出SPEI）
数据资料： 月尺度降水资料及其它用于计算潜在蒸散发所需数据
优点： 考虑了温度因素，可用于气候变暖下的干旱分析。
缺点： 不能满足任意尺度的应用，只能用在月尺度及以上；
计算方法：

2.1 潜在蒸散发（Potential Evapotranspiration, PET）的计算

Thornthwaite法和FAO Penman-Monteith法常常被推荐使用。FAO Penman-Monteith法计算误差小，但需要的气象要素多，Thornthwaite法计算相对简单，需要的气象要素少，但有一定的局限性。
可参见另一篇文章MATLAB利用不同方法实现潜在蒸散发计算

2.2 MATLAB代码实现

3 标准径流指数（Standardized Runoff Index, SRI）

2008-Use of a standardized runoff index for characterizing hydrologic drought（提出SRI指数）
数据资料： 月径流资料
优点： 计算简单。
缺点： 类似SPI指标。
计算方法： 与SPI计算过程类似

function [Z]=SRI(Data,scale,nseas)
%Standardized Runoff Index （标准化径流指数）
% Input Data（输入数据）
% Data : Monthly Data vector not matrix (monthly or seasonal precipitation)
% 数据：月度数据向量而非矩阵
% scale （事件尺度）: 1,3,6,12,48  月尺度
% nseas : number of season (monthly=12) 季节数量 nseas=12
% Example
% Z=SRI(gamrnd(1,1,1000,1),3,12); 3-monthly scale, 
% Notice that  the rest of the months of the fist year are removed.
% eg. if scale =3, fist year data 3-12 SRI values are not estimated.

%if row vector then make coloumn vector 行/列向量转换
%if (sz==1) Data(:,1)=Data;end % 确保数据为列向量

erase_yr=ceil(scale/12);         % 向上取整 (若为1,3,6,12，erase_yr为1)

% Data setting to scaled dataset 数据设置为缩放数据集
A1=[];          % 初始化
for is=1:scale
    A1=[A1,Data(is:length(Data)-scale+is)];    % 按时间尺度列出数据
end
XS=sum(A1,2);        % 对X的每行分别求和

if (scale>1)         % 时间尺度为3、6、12时，进入判断
    XS(1:nseas*erase_yr-scale+1)=[];   
end
% q 径流量为0的概率
% Gam_xs 累积概率
% PARMHAT=gamfit(X);   伽玛分布的参数估计
% 返回拟合X中数据的gamma分布参数的最大似然估计
% PARMHAT（1）和PARMHAT（2）分别是形状和比例参数a和B的估计
% P=gamcdf（X，A，B）返回gamma累积分布函数
% X= NORMINV(P,MU,SIGMA)   其中，P取概率 MU取均值 SIGMA取方差
% 返回值为a    对连续函数，所有小于等于a的值，其出现概率的和为F(a)=P(x<a)
for is=1:nseas
    tind=is:nseas:length(XS);
    Xn=XS(tind);                   % 对应序数
    [zeroa]=find(Xn==0);    % Xn为0的序数
    Xn_nozero=Xn;
    Xn_nozero(zeroa)=[];
    q=length(zeroa)/length(Xn);
    parm=gamfit(Xn_nozero);
    Gam_xs=q+(1-q)*gamcdf(Xn,parm(1),parm(2));
    Z(tind)=norminv(Gam_xs);
end      

end

4 日尺度旱涝急转指数（Dry-wet abrupt alternation index，DWAAI）

指标出自论文：长江中下游流域旱涝急转事件特征分析及其与ENSO的关系-闪丽洁

function [ DWAAI, SPADrought ,SPAFlood ,Fast ,Transition]= computeDWAAI(P , a )
%computeDWAAI 函数用于计算DWAAI

% SPAbefore  前期标准化降水异常值
% SPAafter     后期标准化降水异常值
% SAPI   后期第i天的标准化降水异常值
% SAPI0 前期最后一天的标准化降水异常值
% a      权重系数
Ndrought=44;
Nflood=10;
Plength=length(P);             % 总天数

SPAbeforeN=computeSPAN(P ,Ndrought);
SPAafterN=computeSPAN(P ,Nflood );
API = computeAPI( P ); 
SAPIN = computeSAPI( API );     

n= Plength-Ndrought- Nflood+1;
SPADrought=zeros( n,1 );      % 旱期SPA
SPAFlood=zeros( n,1 );           % 涝期SPA
Fast=zeros( n,1 );                    % “急”的程度
Transition=zeros( n,1 );           % “转”的程度
K=zeros( n,1 );
DWAAI=zeros( n,1 );
for in=1: n
    SPAbefore=SPAbeforeN(in );
    SPAafter=SPAafterN( Ndrought+ in );
    SAP0=SAPIN(in);
    SAPI=SAPIN(  (in+1): (in+Nflood) ); 
    
    for ilength=1:Nflood
        K(in)=K(in)+( SAPI(ilength) -  SAP0 )/ilength;
    end
    DWAAI(in)=( K(in)+( SPAafter-SPAbefore )* (abs( SPAbefore )+ abs( SPAafter )  ) )*a^( -abs( SPAbefore+SPAafter ));
    SPADrought(in)=SPAbefore;      % 旱期SPA
    SPAFlood(in)=SPAafter;             % 涝期SPA
    Fast(in)=K(in) *a^( -abs( SPAbefore+SPAafter )) ;                    % “急”的程度
    Transition(in)=( SPAafter-SPAbefore )* (abs( SPAbefore )+ abs( SPAafter )  )*a^( -abs( SPAbefore+SPAafter ));           % “转”的程度
end
end

函数2：

function SPA = computeSPAN(P ,N)
% computeSPA函数用于计算SPA
% 确保P为列向量

Plength=length(P);
A1=[];          % 初始化
for is=1:N
    A1=[A1,P(is:Plength-N+is)];    % 按时间尺度列出数据
end
XS=sum(A1,2);        % 对X的每行分别求和


XSmean=mean( XS );          % mean precipitation
sigma=sqrt(var( XS ));       % standard deviation

SPA=(XS-XSmean)./sigma;

end

函数3：

function API = computeAPI( X )          % , k
% computeAPI函数用于计算API 
% X 为日降水序列

Xlength=length(X);             % 总天数
Ndays=44;
k=0.955;    % 衰减系数

API =zeros( Xlength-Ndays+1,1 );
for in=1: (Xlength-Ndays+1)
    for jn=1:Ndays
            API(in,1)= API(in,1)+ k^(jn-1)*X( in+Ndays - jn );
    end
end

end

函数4：

function SAPI = computeSAPI( API )     
% computeAPI函数用于计算SAPI 
% 采用Γ分布进行拟合


% PARMHAT=gamfit(X);   伽玛分布的参数估计
% 返回拟合X中数据的gamma分布参数的最大似然估计
% PARMHAT（1）和PARMHAT（2）分别是形状和比例参数a和B的估计
% P=gamcdf（X，A，B）返回gamma累积分布函数
% X= NORMINV(P,MU,SIGMA)   其中，P取概率 MU取均值 SIGMA取方差
% 返回值为a    对连续函数，所有小于等于a的值，其出现概率的和为F(a)=P(x<a)  

[zeroa]=find(API==0);    % Xn为0的序数
API_nozero=API;
API_nozero(zeroa)=[];
q=length(zeroa)/length(API);
parm=gamfit(API_nozero);
Gam_xs=q+(1-q)*gamcdf(API,parm(1),parm(2));
SAPI=norminv(Gam_xs); 

end

5 标准化加权平均降水量（Weighted Average Precipitation，SWAP）

J2009-Determining the start, duration, and strength of flood and drought with daily precipitation: Rational
J2014-The day-to-day monitoring of the 2011 severe drought in China
数据资料： 日尺度降水资料
优点： 计算简单；利用日尺度降水进行分析，。
缺点： 忽略气温变化对干旱影响

5.1 日尺度SWAP计算

MATLAB代码实现：

function SWAP = getSWAP(P)
% getSWAP函数用于计算SWAP
% a 权重随时间衰减参数 [0,1]
daysOfYear=365;
%% 计算WAP
% 1）a=0.9   N=44
% 2）a=0.95 N=60
a=0.9;
N=44;
omega=zeros(N+1,1);
for n=0:N
     omega(n+1)= (1-a)*a^n/(1-a^(N+1));
%     if  n==0
%         omega(n+1)= 1-a;
%     else
%         omega(n+1)= (1-a)*a^n/(1-a^(N+1));
%     end
end

ndays=length(P) - N ;
WAP=zeros( ndays , 1 );
for idays=1:ndays
    for n=0:N
        WAP(idays,1)=WAP(idays,1)+omega(n+1)*P(idays+N-n);
    end
end
% 将第一年的其他值删去
WAP=WAP(daysOfYear-N+1 :end ,1);
%% 进行标准化WAP->SWAP
% 按对应日序数逐年排列，根据多年同一日WAP值构建Gamma分布，并正态标准化，得到对应SWAP值
XS=WAP;
SWAP=zeros( length( XS ),1);
for is=1:daysOfYear
    tind=is:daysOfYear:length(XS);
    Xn=XS(tind);                   % 对应序数
    [zeroa]=find(Xn==0);    % Xn为0的序数
    Xn_nozero=Xn;
    Xn_nozero(zeroa)=[];
    q=length(zeroa)/length(Xn);
    parm=gamfit(Xn_nozero);
    Gam_xs=q+(1-q)*gamcdf(Xn,parm(1),parm(2));
    SWAP(tind)=norminv(Gam_xs);
end      

% Xn=WAP;
% [zeroa]=find(Xn==0);    % Xn为0的序数
% Xn_nozero=Xn;
% Xn_nozero(zeroa)=[];
% q=length(zeroa)/length(Xn);
% parm=gamfit(Xn_nozero);
% Gam_xs=q+(1-q)*gamcdf(Xn,parm(1),parm(2));
% SWAP=norminv(Gam_xs); 

end

5.2 日以上时间尺度SWAP计算

根据“Determining the start, duration, and strength of flood and drought with daily precipitation: Rational”一文中 “Once the daily values of the flood and drought extent are determined, the means of a month and longer periods can be calculated.”，推测日以上时间尺度SWAP为相应时段内的平均值。
但假若月尺度SWAP为一个月SWAP的平均值是否还服从标准的正态分布？
显然并不符合，统计学告诉我们，

也就是说，月尺度SWAP和日尺度SWAP的分级阈值应当是不同的。

6 有效干旱指数（Effective drought index，EDI）

J1999-Objective Quantification of Drought Severity and Duration-首次定义EDI
数据资料： 30年以上月降水数据
优点： 日时间尺度（daily time step）；使用可变时间段来计算降水量。
缺点：
计算步骤：

MATLAB代码实现：（暂时还没有学如何编程）

7 指数比较

7.1 SPI与SWAP比较

SWAP为WAP按年序排列，由Gamma函数拟合后转化为正态分布，可见，SWAP各月末日的计算结果与SPI的计算结果相同（在上述参数取值下）。
事实上，a值和N值对SWAP的范围并不会有太大影响。在转化为标准正态分布后是否取值范围十分接近？我迷惑了，求小伙伴解答…
将某站点月尺度SPI与月尺度SWAP序列进行比较，可以看出，SPI的取值范围更广。显然，运用相同的等级划分标准是不妥的，那么在使用SWAP值确定干旱的等级时应如何定义？

将月尺度SPI与日尺度SWAP序列绘制于图中，可见SWAP的取值范围会更广。

可以知道，SWAP在运用不同尺度进行分析时，相应的分级阈值应当是不同的。

来源：WW、forever