气象干旱综合指数MCI

气象干旱综合指数MCI

 

% MCI = Ka * (a*SPIW60 + b*MI30 + c*SPI90 +d*SPI150)
% a    SPIW   权重系数，北方及西部地区取 0.3,南方地区取 0.55
% b    MI     权重系数，北方及西部地区取 0.5,南方地区取 0.6；
% c    SPI90  权重系数，北方及西部地区取 0.3,南方地区取 0.2；
% d    SPI150 权重系数，北方及西部地区取 0.2,南方地区取 0.1；
% Ka   为季节调节系数，根据不同季节各地主要农作物生长发育阶段对土壤水分的敏感程度

clc;clear;
data = xlsread('三门峡1957.1.1-2021.3.1.xlsx');
data(any((data(:,3)*100+data(:,4)==229),2),:)=[];
P = data(:,5);
daysOfYear=365;
% a=0.3;b=0.5;c=0.3;d=0.2;%北方
a=0.55;b=0.6;c=0.2;d=0.1;%南方
Ka = [0.4 0.8 1.0 1.2 1.2 1.2 1.2 1 1 0.8 0.6 0.4];
% 计算SWAP60
a=0.85; N=60;
omega=zeros(N+1,1);
for n=0:N
     omega(n+1)= a^n;
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);
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      

% 计算SPI90
scale90= 90;
P90 = zeros(length(P)-scale90+1,1);
for j=1:length(P90)
    P90(j+scale90-1)=sum(P(j:j+scale90-1));
end
% 将第一年的其他值删去
P90=P90(366 :end ,1);
SPI90=zeros( length( P90 ),1);
for is=1:daysOfYear
    tind=is:daysOfYear:length(P90);
    Xn=P90(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));
    SPI90(tind)=norminv(Gam_xs);
end      

% 计算SPI150
scale150= 150;
P150 = zeros(length(P)-scale150+1,1);
for j=1:length(P150)
    P150(j+scale150-1)=sum(P(j:j+scale150-1));
end
% 将第一年的其他值删去
P150=P150(366 :end ,1);
SPI150=zeros( length( P150 ),1);
for is=1:daysOfYear
    tind=is:daysOfYear:length(P150);
    Xn=P150(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));
    SPI150(tind)=norminv(Gam_xs);
end 

%%%计算PET= E0  注意表格数据气压单位百帕 风俗 两米风速
ELE = data(1,14); LAT = data(1,13);as = 0.25;bs = 0.5;
M = data(:,3);
Ri = data(:,4); J = (M-1)*30.4 + Ri;%日
R = data(:,5);            %降水量
Pre = data(:,6);            %气压
Tmean = data(:,7);        %平均温度
Ws = data(:,10);           %风速
Hum = data(:,11);         %湿度
Sun = data(:,12);         %日照
u2=Ws;%u2为2米高处的风速
%∆为饱和水汽压曲线的斜率
delta=4098*0.6108*exp(17.27*Tmean./(Tmean+237.3))./(Tmean+237.3).^2;
eS=0.6108*exp(17.27*Tmean./(Tmean+237.3));  %饱和水汽压
eA=eS.*Hum/100;                    %ea为实际水汽压
gamma1 = 0.665/10^3 *Pre/10;%r为湿度计常数，kPa/℃
Phi=LAT*pi/180;     %纬度 由度数化成弧度
%黄赤交角，是指地球公转轨道面（黄道面）与赤道面（天赤道面）的交角
Delt=0.408*sin(2*pi*J/365-1.39);
ws=acos(-tan(Phi)*tan(Delt));  %日落角
dr=1+0.033*cos(2*pi*J/365); %地 日 距离
Gsc=0.0820;  %太阳常数
%Ra为用相对应的蒸发单位表达的地球辐射值，mm/d
Ra=((24*60)/pi)*Gsc*dr.*(ws.*sin(Phi).*sin(Delt)+cos(Phi).*cos(Delt).*sin(ws));
Rso=(as+bs+2*ELE/100000).*Ra;%晴空辐射
N=24*ws/pi;                  %可能日照时数
Rs=(as+bs*(Sun./N)).*Ra;
alpha=0.23;          %作物反射率
Rns=(1-alpha)*Rs;    %短波辐射
%计算净长波辐射RnL
Tkmax=Tmean+273.16;Tkmin=Tmean+273.16;
fcd=1.35*Rs./Rso-0.35;
Rnl=(4.903*10^(-9))*fcd.*(0.34-0.14*sqrt(eA)).*((Tkmax.^4+Tkmin.^4)/2);
Rn=Rns-Rnl;
E0=(0.408 .* delta .*Rn+gamma1 *900 ./(Tmean+273) .*u2 .*(eS-eA))./(delta+gamma1 .*(1+0.34*u2));

% 计算MI30
scale30= 30;
P30 = zeros(length(P)-scale30+1,1);
E30 = zeros(length(P)-scale30+1,1);
for j=1:length(P30)
    P30(j+scale30-1)=sum(P(j:j+scale30-1));
    E30(j+scale30-1)=sum(E0(j:j+scale30-1));
end
MI30 = (P30-E30) ./E30;
MI30 = MI30(366:end,1);

MCI0 = a*SWAP + b*MI30 + c*SPI90 +d*SPI150;
MCI0 = [M(366:end) Ri(366:end) MCI0];
X = find(MCI0(:,1)==1);MCI(X,1)=Ka(1)*MCI0(X,3);
X = find(MCI0(:,1)==2);MCI(X,1)=Ka(2)*MCI0(X,3);
X = find(MCI0(:,1)==3);MCI(X,1)=Ka(3)*MCI0(X,3);
X = find(MCI0(:,1)==4);MCI(X,1)=Ka(4)*MCI0(X,3);
X = find(MCI0(:,1)==5);MCI(X,1)=Ka(5)*MCI0(X,3);
X = find(MCI0(:,1)==6);MCI(X,1)=Ka(6)*MCI0(X,3);
X = find(MCI0(:,1)==7);MCI(X,1)=Ka(7)*MCI0(X,3);
X = find(MCI0(:,1)==8);MCI(X,1)=Ka(8)*MCI0(X,3);
X = find(MCI0(:,1)==9);MCI(X,1)=Ka(9)*MCI0(X,3);
X = find(MCI0(:,1)==10);MCI(X,1)=Ka(10)*MCI0(X,3);
X = find(MCI0(:,1)==11);MCI(X,1)=Ka(11)*MCI0(X,3);
X = find(MCI0(:,1)==12);MCI(X,1)=Ka(12)*MCI0(X,3);
MCI = [data(366:end,2:4),MCI];