SPI指数(Standardized Precipitation Index,标准化降水指数)是反映干湿状况的一个指标,主要计算步骤如下:

收集研究区域过去30年或以上时间尺度(一般选取30年)的月降水量资料。

对月降水量资料进行统计分析,拟合出最适合的概率分布函数。常用的有Pearson III 分布、Gamma分布等。

根据所选取的概率分布函数,估计出各个时间尺度下的平均值和标准差。

对于任意一个时间尺度,用其降水量减去该时间尺度下的平均值,再除以标准差,即可计算出该时间的SPI值。

根据SPI的值确定干湿状况。一般来说,SPI>0表示湿润,SPI<0表示干旱。干旱和湿润的强度根据SPI绝对值的大小判断。

绘制不同时间尺度下的SPI变化曲线,分析各个时间尺度的干湿状况。

综上,SPI指数借助长期历史资料,能很好地反映不同时间尺度下的干湿状况,是评估干旱的重要指标之一。

MATLAB代码:

%% SPI指数计算

clc;close all;clear all;%清除空间

%% 载入数据

data=xlsread('降水.xls');

x=data(:,3);

y=data(:,1);

% x=data(:,2);

% y=data(:,1);

%% 计算伽马分布的参数

% %生成降雨数据

% N=1000;%天数

% x=randi([0,100],N,1);%降雨量

n=length(x);

H1= x==0;

m=sum(H1);%0的项数

% A=sum(log(x(x~=0)))/m-log(mean(x));

% A=log(mean(x))-sum(log(x(x~=0)))/m;

x2=x(x~=0);

% [alpha,beta] = gamma_fit(x2);

[p,ci] = gamfit(x2);

alpha=p(1)

beta=p(2)

% alpha=(1+sqrt(1+4*A/3))/(4*A);

% beta=mean(x)/alpha;

q=m/n;

%% 参数设定

c0=2.515517;

c1=0.802853;

c2=0.010328;

d1=1.432788;

d2=0.189269;

d3=0.001308;

% 计算SPI指数

SPI=SPIfun(q,alpha,beta,c0,c1,c2,d1,d2,d3,x);

figure;

plot(y,SPI,'b-');

hold on;

plot(y,zeros(length(y),1),'r-');

xlabel('时间');

ylabel('SPI');

title('降雨量SPI');

function [a,b] = gamma_fit(x,s)

% GAMMA_FIT     Maximum-likelihood gamma distribution.

%

% GAMMA_FIT(x) returns the MLE (a,b) for the data in vector x.

%

% GAMMA_FIT(m,s) returns the MLE (a,b) for data with sufficient statistics

% given by

%   m = mean(x)

%   s = log(m) - mean(log(x))

%

% The gamma distribution is parameterized as

%   p(x) = x^(a-1)/(Gamma(a) b^a) exp(-x/b)

%   E[x] = ab

%

% The algorithm is a generalized Newton iteration, described in

% "Estimating a Gamma distribution" by T. Minka.

% Written by Tom Minka

if nargin == 1

    m = mean(x);

    s = log(m) - mean(log(x));

else

    % suff stats given

    m = x;

end

a = 0.5/s;

if 0

    % lower bound

    for iter = 1:1000

        old_a = a;

        a = inv_digamma(log(a) - s);

        if(abs(a - old_a) < 1e-8)

            break;

        end

    end

end

% gen Newton

for iter = 1:100

    old_a = a;

    g = log(a)-s-digamma(a);

    h = 1/a - trigamma(a);

    a = 1/(1/a + g/(a^2*h));

    if(abs(a - old_a) < 1e-8)

        break;

    end

end

b = m/a;

function H=Hfun(q,alpha,beta,x)

%% 计算累积概率

% gamcdf(x,a,b)

G=gamcdf(x,alpha,beta);

H=q+(1-q)*G;

function SPI=SPIfun(q,alpha,beta,c0,c1,c2,d1,d2,d3,x)

%% 计算SPI

H=Hfun(q,alpha,beta,x);

SPI=zeros(1,length(H));

for i=1:length(H)

    if H(i)<=0.5

        k=sqrt(log(1/(H(i).^2)));

        SPI(i)=-(k-(c0+c1*k+c2*k^2)/(1+d1*k+d2*k^2+d3*k^3));

    else

        k=sqrt(log(1/(1-H(i))^2));

        SPI(i)=k-(c0+c1*k+c2*k^2)/(1+d1*k+d2*k^2+d3*k^3);

    end

end

function k=kfun(q,alpha,beta,c0,c1,c2,d1,d2,d3,x)

%% 计算k

H=Hfun(q,alpha,beta,x);

for i=1:length(H)

    if H(i)<=0.5

       

    else

       

    end

end

程序结果:

alpha =

    0.6792

beta =

   10.4584

>>