基于灾损曲线的综合风险与指标体系综合风险验证

1. 计算基于灾损曲线的综合风险

（1）首先得到不同灾害的经济损失、人口损失曲线（表达式）

（2）计算空间上灾害的重现期，重现期越大，代表灾害程度越严重，发生可能性越小。计算过程：得到灾害各危险性要素37年间的内蒙古全域平均值，构建单变量边际分布，选择最优copula，三种灾害均为Clayton Copula函数，计算联合重现期，计算联合概率。

% 加载内蒙古全域频率数据和强度数据
clc;clear;
fname = "G:\Data\Drought\disaster-loss-curve\frequency_neimenggu_mean.csv";
freq = load(fname); % 频率数据
mname = "G:\Data\Drought\disaster-loss-curve\magnitude_neimenggu_mean.csv";
mag = load(mname); % 强度数据

X = freq;
Y = abs(mag);
alpha=0.05;
% 得出结果：freq服从gev分布(p=0.9696)，mag服从gev分布(p=0.8333)
%通常选择通过检验的分布，以gev和ev拟合和kstest和ADtest为例
[para0,pci0] = gevfit(X);
U=gevcdf(X,para0(1),para0(2),para0(3)); % 返回数组X在gev分布的形状参数、范围参数以及位置参数处的cdf概率值
[hx0,px0]=kstest(X,[X,U],alpha);  % h = 0 接受原假设，服从假定的分布

[para1,pci1] = gevfit(Y);
V=gevcdf(Y,para1(1),para1(2),para1(3)); % 返回数组X在gev分布的形状参数、范围参数以及位置参数处的cdf概率值
[hx1,px1]=kstest(Y,[Y,V],alpha);  % h = 0 接受原假设，服从假定的分布

% 读取nc文件
frequency_filename = "G:\Data\Drought\drought_indices\by_year\drought_frequency_yearmean_1985_2021.nc";
magnitude_filename = "G:\Data\Drought\drought_indices\by_year\drought_magnitude_yearmean_1985_2021.nc";
frequency = ncread(frequency_filename, 'spei_01_frequency');
magnitude = ncread(magnitude_filename, 'spei_01_magnitude');
latitude = ncread(magnitude_filename, 'lat');
longitude = ncread(magnitude_filename, 'lon');

%%
orperiod = zeros(301,171);
for i=1:length(longitude)
    for j=1:length(latitude)
        freq_series = frequency(i, j); % 加载volume变量37个时间序列数据
        magnitude_series = magnitude(i, j); % 加载magnitude变量37个时间序列数据
        magnitude_series = abs(magnitude_series);

        U0 = gevcdf(freq_series,para0(1),para0(2),para0(3));
        V0 = gevcdf(magnitude_series,para1(1),para1(2),para1(3));
        if U0>0.99
           U0=0.99;
        end
        if V0>0.99
            V0=0.99;
        end
        
        % 尝试 Clayton Copula 拟合 若拟合失败，给参数赋予9999；
        paramhat2 = copulafit('Clayton',[U(:), V(:)]);
        Ccla = copulacdf('Clayton',[U(:), V(:)],paramhat2);
       
        CClayor = copulacdf('Clayton',[U0, V0],paramhat2);
        orperiod(i,j) = 1./(1-CClayor);      
    end
end

%%

% 将结果写入nc
ncname = 'G:\Data\Drought\drought_indices\by_year\drought_dangerous1.nc';
% 构造变量名和变量尺寸
nccreate(ncname,'orperiod',...
    'Dimensions', {'longitude',length(longitude),'latitude',length(latitude)},...
    'FillValue','disable');
% nccreate(ncname,'longitude','Dimensions',{'latitude',length(latitude),'longitude',length(longitude)});
% nccreate(ncname,'latitude','Dimensions',{'latitude',length(latitude),'longitude',length(longitude)});
% 这里Fillvalue设置为disable是不允许matlab用默认值填充你的矩阵
nccreate(ncname,'longitude','Dimensions',{'longitude',length(longitude)});
nccreate(ncname,'latitude','Dimensions',{'latitude',length(latitude)});

ncwrite(ncname,'longitude',longitude);
ncwrite(ncname,'latitude',latitude);
ncwrite(ncname,'orperiod',orperiod); % 将数据储存到该变量中

ncwriteatt(ncname,'orperiod','units','y'); % 给该变量添加描述
ncwriteatt(ncname,'orperiod','missing_value','nan'); % 给该变量添加描述：nan为copula构建失败的值

（3）得到空间上的联合重现期后，在python中根据灾损曲线分别计算经济损失（风险）、人口损失（风险）。

（4）加载到arcgis中，在工具箱中打开Band Collection Statistics工具，计算空间相关系数。

 2. 得到相关性系数，结果：