风险等级的计算 - 上方山地区森林火灾风险评估组

风险等级的计算

2021-07-01 00:00:00

计算上方山国家森林公园地区的火灾危险性、危害性和综合风险等级，需要用到大量的数据。由于计算量过大，所以我们采用MatLab编程进行计算，获得了研究区域内每个栅格的危险性、危害性以及综合风险等级。

以下是我们所使用的程序代码：

1.层次分析法算权重
function W = ccfx(filename,sheet,xlRange)
%% 读取判断矩阵
% disp('请输入判断矩阵A：')
% A = input('判断矩阵A= ');
A = xlsread(filename,sheet,xlRange);
[m,n] = size(A);
% %% 判断是否为正互反矩阵
% for i = 1:m
% for j = 1:n
% if (A(i,j)* A(j,i)~= 1)
% disp('矩阵不是正互反矩阵！')
% W = [];
% return
% end
% end
% end

%% 算术平均值求权值

%对列归一化
Sum_A1 = sum(A);
SUM_A1 = repmat(Sum_A1,m,1);
Standard_A = A./SUM_A1;

%将归一化的矩阵按行求和
Sum_A2 = sum(Standard_A,2);

%将相加后得到的向量除以n得到权重向量
W_mean = Sum_A2/n;
% disp('算数平均值求权重结果为：')
% disp(W_mean)

%% 几何平均法求权重

%将A元素按行相乘得到一个新的列向量
%prod函数和max函数类似，一个用于乘，一个用于加
Product_A1 = prod(A,2);

%将新的向量的每个分量开n次方
Product_A2 = Product_A1.^(1/n);

%对该列向量进行归一化即可得到权重向量
W_geometry = Product_A2/sum(Product_A2);
% disp('几何平均值求权重结果为：')
% disp(W_geometry)

%% 特征值求权重
%求特征值和特征向量
%V为特征值对应的特征向量，D为特征值构成的对角阵
[V,D] = eig(A);

%找出最大特征值及所对应的特征向量
Max_eig = max(D(:));
%找到D中第一个与最大特征值相等的元素的位置
[r,c] = find(D == Max_eig,1);
feature_vector = V(:,c);
%求权重
W_feature = feature_vector/sum(feature_vector);
% disp('特征值求权重结果为：')
% disp(W_feature)

%% 判断是否为一致
if n==2
W = (W_feature+W_geometry+W_mean)/3;
return;
end
CI = (Max_eig - n) / (n-1);
RI = [0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.46 1.49 1.52 1.54 1.56 1.58 1.59];
CR = CI / RI(n);
if CR < 0.1
% disp('CR<0.1，判断矩阵A可接受!');
W = (W_feature+W_geometry+W_mean)/3;
% disp('权重结果为：')
% disp(W)
else
disp(strcat('CR= ',num2str(CR)))
disp('CR>=0.1，判断矩阵A不可接受!需进行修改！');
return;
end
end

2.计算风险等级
clc,clear
%数据表地址
filename = 'C:\Users\xyn\Desktop\暑期实习\隶属矩阵数据.xlsx';
filename_A = 'C:\Users\xyn\Desktop\暑期实习\权重矩阵数据.xlsx';
sheet = [1 2 3];
xlRange_data = 'B2:M59488';
m = 8;

%获取总体栅格属性矩阵
R_RasterData = xlsread(filename,sheet(3),xlRange_data);
[Raster_num,Attribute_num] = size(R_RasterData);

%% 计算危险性下的各级指标权重
%获取权重向量
A__ = ccfx(filename_A,1,'B2:E5');%二级指标权重
% A__ = ccfx(filename_A,1,'B2:D4');%二级指标权重
A__1 = ccfx(filename_A,1,'B8:C9');%可燃物因子中三级指标权重
A__2= ccfx(filename_A,1,'B13:D15');%气象因子中三级指标权重
A__3=ccfx(filename_A,1,'B18:E21');%立地因子中三级指标权重
A__4=ccfx(filename_A,1,'B24:D26');%火源管理因子中三级指标权重

%% 计算危害性下的各级指标权重
%获取权重向量
A_ = ccfx(filename_A,2,'B2:E5');%二级指标权重
% A_ = ccfx(filename_A,2,'B2:D4');%二级指标权重
A_1 = ccfx(filename_A,2,'B8:C9');%可燃物因子中三级指标权重
A_2= ccfx(filename_A,2,'B13:D15');%气象因子中三级指标权重
A_3=ccfx(filename_A,2,'B18:E21');%立地因子中三级指标权重
A_4=ccfx(filename_A,2,'B24:D26');%火源管理因子中三级指标权重

%% 计算危险性下的各栅格等级评价
R_Vegetation_types = xlsread(filename,sheet(1),'B3:D8');%植被类型
% R_hierarchical_structure = xlsread(filename,sheet(1),xlRange);%层次结构
% R_Canopy_density = xlsread(filename,sheet(1),xlRange);%郁闭度
R_Shrub_layer_coverage = xlsread(filename,sheet(1),'B12:D16');%植被覆盖度
% R_herb_coverage = xlsread(filename,sheet(1),'B16:D18');%草本层盖度

R_average_temperature = xlsread(filename,sheet(1),'B20:D22');%月平均温度
R_average_precipitation = xlsread(filename,sheet(1),'B24:D26');%月平均降水
R_average_wind_speed = xlsread(filename,sheet(1),'B28:D29');%月平均相对湿度

R_altitude = xlsread(filename,sheet(1),'B32:D34');%海拔
R_Aspect = xlsread(filename,sheet(1),'B36:D39');%坡向
R_landforms = xlsread(filename,sheet(1),'B41:D43');%地貌
R_Slope = xlsread(filename,sheet(1),'B45:D47');%坡度

R_Artificial_intensity = xlsread(filename,sheet(1),'B49:D51');%人为强度
R_Road_network_density = xlsread(filename,sheet(1),'B53:D54');%防火宣传
R_Tourism_development = xlsread(filename,sheet(1),'B57:D58');%旅游开发

%计算危险性等级
for k = 1:Raster_num
%获取单个栅格属性矩阵
Rasterk = R_RasterData(k,:);
%获取隶属矩阵
% R1 = [R_Vegetation_types(Rasterk(1),:);R_Shrub_layer_coverage(Rasterk(2),:);R_herb_coverage(Rasterk(3),:)];
R1 = [R_Vegetation_types(Rasterk(1),:);R_Shrub_layer_coverage(Rasterk(2),:)];
R2 = [R_average_temperature(Rasterk(3),:);R_average_precipitation(Rasterk(4),:);R_average_wind_speed(Rasterk(5),:)];
R3 = [R_altitude(Rasterk(6),:);R_Aspect(Rasterk(7),:);R_landforms(Rasterk(8),:);R_Slope(Rasterk(9),:)];
R4 = [R_Artificial_intensity(Rasterk(10),:);R_Road_network_density(Rasterk(11),:);R_Tourism_development(Rasterk(12),:)];
B1 = A__1'*R1;
B2 = A__2'*R2;
B3 = A__3'*R3;
B4 = A__4'*R4;
B_DANGER(k,:) = A__'*[B1;B2;B3;B4];
% B_DANGER(k,:) = A__'*[B1;B2;B3];
C1(k,:) = B_DANGER(k,:)./sum(B_DANGER(k,:));
disp(strcat('栅格',num2str(k),'危险性等级已计算成功！'))
end

%% 计算危害性下的各栅格等级评价
R_Vegetation_types = xlsread(filename,sheet(2),'B3:D8');%植被类型
% R_hierarchical_structure = xlsread(filename,sheet(1),xlRange);%层次结构
% R_Canopy_density = xlsread(filename,sheet(1),xlRange);%郁闭度
R_Shrub_layer_coverage = xlsread(filename,sheet(2),'B12:D16');%植被覆盖率
% R_herb_coverage = xlsread(filename,sheet(2),'B16:D18');%草本层盖度

R_average_temperature = xlsread(filename,sheet(2),'B20:D22');%月平均温度
R_average_precipitation = xlsread(filename,sheet(2),'B24:D26');%月平均降水
R_average_wind_speed = xlsread(filename,sheet(2),'B28:D29');%月相对湿度

R_altitude = xlsread(filename,sheet(2),'B32:D34');%海拔
R_Aspect = xlsread(filename,sheet(2),'B36:D39');%坡向
R_landforms = xlsread(filename,sheet(2),'B41:D43');%地貌
R_Slope = xlsread(filename,sheet(2),'B45:D47');%坡度

R_Artificial_intensity = xlsread(filename,sheet(2),'B49:D51');%人为强度
R_Road_network_density = xlsread(filename,sheet(2),'B53:D54');%防火宣传
R_Tourism_development = xlsread(filename,sheet(2),'B57:D58');%旅游开发

%计算危害性等级
disp(' ')
for k = 1:Raster_num
%获取单个栅格属性矩阵
Rasterk = R_RasterData(k,:);
%获取隶属矩阵
% R1 = [R_Vegetation_types(Rasterk(1),:);R_Shrub_layer_coverage(Rasterk(2),:);R_herb_coverage(Rasterk(3),:)];
R1 = [R_Vegetation_types(Rasterk(1),:);R_Shrub_layer_coverage(Rasterk(2),:)];
R2 = [R_average_temperature(Rasterk(3),:);R_average_precipitation(Rasterk(4),:);R_average_wind_speed(Rasterk(5),:)];
R3 = [R_altitude(Rasterk(6),:);R_Aspect(Rasterk(7),:);R_landforms(Rasterk(8),:);R_Slope(Rasterk(9),:)];
R4 = [R_Artificial_intensity(Rasterk(10),:);R_Road_network_density(Rasterk(11),:);R_Tourism_development(Rasterk(12),:)];
B1 = A_1'*R1;
B2 = A_2'*R2;
B3 = A_3'*R3;
B4 = A_4'*R4;
B_HARM(k,:) = A_'*[B1;B2;B3;B4];
% B_HARM(k,:) = A_'*[B1;B2;B3];
C2(k,:) = B_HARM(k,:)./sum(B_HARM(k,:));
disp(strcat('栅格',num2str(k),'危害性等级已计算成功！'))
end

%% 计算风险等级
A = [0.5 0.5];
for k = 1:Raster_num
R = [B_DANGER(k,:) ; B_HARM(k,:)];
B = A*R;
C(k,:) = B./sum(B);%归一化
end

%% 将数据导出至excel表
outputfile = strcat('C:\Users\xyn\Desktop\暑期实习\导出数据\outputdata',num2str(m),'.xlsx');
num = 1:Raster_num;
% 将数据组集到data
data1 = [num', C];
data2 = [num', C1];
data3 = [num', C2];

% 添加变量名称
title1 = {'栅格', 'Ⅰ', 'Ⅱ','Ⅲ'};
title2 = {'栅格', 'Ⅰ', 'Ⅱ','Ⅲ'};
title3 = {'栅格', 'Ⅰ', 'Ⅱ','Ⅲ'};

% 将data写入到outputdata.xlsx文件中
xlRange = strcat('A1:D',num2str(Raster_num+1));
% s = xlswrite(outputfile, [title1;data1],1,xlRange_t);
s = xlswrite(outputfile, [title1;num2cell(data1)],'风险等级',xlRange);
% s = xlswrite(outputfile, title2,2,xlRange_t);
s = xlswrite(outputfile, [title2;num2cell(data2)],'危险性等级',xlRange);
% s = xlswrite(outputfile, title3,3,xlRange_t);
s = xlswrite(outputfile, [title3;num2cell(data3)],'危害性等级',xlRange);

disp('数据导出完成！')

Classify(outputfile,m,Raster_num);

3.导出数据分级
function Classify(filename,m,Raster_num)
%% 评定风险等级
% filename = 'C:\Users\xyn\Desktop\暑期实习\outputdata.xlsx';
xlRange_data = strcat('B2:D',num2str(Raster_num+1));
Data1 = xlsread(filename,'风险等级',xlRange_data);
Data2 = xlsread(filename,'危险性等级',xlRange_data);
Data3 = xlsread(filename,'危害性等级',xlRange_data);
[Raster_num,Attribute_num] = size(Data1);
%评定风险等级
for k = 1:Raster_num
Max = max(Data1(k,:));
ind = find(Data1(k,:) == Max);
C1(k,1) = ind;
end
%评定危险性等级
for k = 1:Raster_num
Max = max(Data2(k,:));
ind = find(Data2(k,:) == Max);
C2(k,1) = ind;
end
%评定危害性等级
for k = 1:Raster_num
Max = max(Data3(k,:));
ind = find(Data3(k,:) == Max);
C3(k,1) = ind;
end

%% 将数据导出至excel表
outputfile = strcat('C:\Users\xyn\Desktop\暑期实习\导出数据\ClassifiedData',num2str(m),'.xlsx');
num = 1:Raster_num;
% 将数据组集到data
data = [num', C1,C2,C3];

% 添加变量名称
title = {'栅格', '风险等级', '危险性等级','危害性等级'};

end

上方山地区森林火灾风险评估组

Addres

Email

Telephone

Contact