变强从来不是因为社会的阳光面，而是为了当社会的阴暗面向我或我们在意的人伸手时，我们能够干净利索地绞杀它们。
游程理论识别灾害事件---基于MATLAB语言的编程实现
前言
1. 版本
1.1 2021年3月16日，Version 1
2. 摘要
- 游程理论，灾害事件，MATLAB
3. 微信公众号GISRSGeography
- 公众号GISRSGeography的内容涉及GIS，遥感和作物模型等的内容，会坚持更新，
欢迎大家关注，谢谢！。
一、游程理论简述
游程理论通常用于识别灾害事件，图1展示了基于游程理论识别事件事件的概念图。当灾害指数低于某一阈值时，且持续时间超过一定长度时则认为灾害事件发生。其中灾害事件的持续时间就是灾害事件持续时间，灾害事件持续时间内灾害指数之和就是灾害事件的严重性，灾害严重性除以灾害持续时间就是灾害强度。
图1. 游程理论识别灾害事件概念图
二、游程理论识别灾害事件的MATLAB程序实现
1. DisasterEventExtra.m函数-灾害事件提取函数实现
基于MATLAB2020编程实现游程理论，具体函数代码如下：
%{
1. 函数目的：基于游程理论提取灾害事件
2. 版本
2.1 山东青岛 2021年3月11日 Version 1
2.2 山东青岛 2021年3月15日 Version 2
(1) 针对干旱事件可能不存在的情况进行修改
3. 输入参数
(1) DS_Array：array,灾害数据时间序列，n行*3列【年，月，灾害指数】
注意：数组的前两列并不一定是年月，如果修改，函数的某些部分可能需要修改
(2) threshold_value：灾害事件识别阈值
4. 返回参数
(1) ds_event_properties_table：table, 灾害事件统计结果，包括：
1) 灾害事件数目 2) 灾害事件开始的index 3) 灾害事件结束的index
4) 灾害事件开始年 5) 灾害事件开始月 6)灾害事件结束年 7) 灾害事件结束月
8) 灾害事件持续时间 9) 灾害事件严重性 10) 灾害事件强度
%}
%%
function ds_event_properties_table = DisasterEventExtra(DS_Array,threshold_value)
ds_series = DS_Array(:,3); % 灾害指数时间序列
% 1. 获取小于某一特定灾害阈值的灾害序列的index
index_ds = find(ds_series <= threshold_value);
% 2. 剔除灾害事件长度为1的灾害事件
index_ds_array = nan(length(index_ds),5);
index_ds_array(:,1) = index_ds;
index_ds_array(1:length(index_ds)-1,2) = index_ds(2:length(index_ds))-index_ds(1:length(index_ds)-1);
index_ds_array(2:length(index_ds),3) = index_ds(2:length(index_ds))-index_ds(1:length(index_ds)-1);
index_ds_array(index_ds_array(:,3)==1 | index_ds_array(:,2)==1 ,4) = 1;
index_ds_array(index_ds_array(:,3)==1 | index_ds_array(:,2)==1 ,5) = index_ds_array(index_ds_array(:,3)==1 | index_ds_array(:,2)==1,1);
index_event_ge2 = index_ds_array(index_ds_array(:,3)==1 | index_ds_array(:,2)==1 ,5); % 灾害事件持续时间大于1的index
if ~isempty(index_event_ge2) % 如果灾害事件数目不为0，统计灾害事件相关信息
event_array = nan(nansum(index_ds_array(:,4)),3);
event_array(:,1) = index_event_ge2; %灾害事件持续时间大于1的事件的Index
% 3. 获取每个灾害事件的开始时间和结束时间的index
event_array(2:length(event_array),2) = event_array(2:length(event_array),1)-event_array(1:length(event_array)-1,1);
event_array(1:length(event_array)-1,3) = event_array(2:length(event_array),1)-event_array(1:length(event_array)-1,1);
% 3.1 灾害事件数目
event_count = sum(event_array(:,2)>1)+1;
% 3.2 每个灾害事件的开始时间index
start_index_array = nan(event_count,1);
start_index_array(1,1) = event_array(1,1);
start_index_array(2:end,1) = event_array(event_array(:,2)>1,1);
% 3.3 每个灾害事件的结束时间index
end_index_array = nan(event_count,1);
end_index_array(1:event_count-1,1) = event_array(event_array(:,3)>1,1);
end_index_array(end,1) = event_array(end,1);
% 3.4 灾害事件属性-开始时间，结束时间，持续时间
ds_event_properties = nan(event_count,10);
ds_event_properties(:,1) = event_count;
ds_event_properties(:,2) = start_index_array;
ds_event_properties(:,3) = end_index_array;
ds_event_properties(:,4) = DS_Array(start_index_array,1);
ds_event_properties(:,5) = DS_Array(start_index_array,2);
ds_event_properties(:,6) = DS_Array(end_index_array,1);
ds_event_properties(:,7) = DS_Array(end_index_array,2);
ds_event_properties(:,8) = (end_index_array - start_index_array) + 1;
% 3.5 灾害事件属性-灾害严重性，灾害强度
for ii_event = 1:event_count
st_index = ds_event_properties(ii_event,2);
ed_index = ds_event_properties(ii_event,3);
ds_event_properties(ii_event,9) = sum(ds_series(st_index:ed_index)); % 灾害严重性
ds_event_properties(ii_event,10) = mean(ds_series(st_index:ed_index)); % 灾害强度
end
% 3.6 灾害事件属性整理
table_name = {'EventCount','StIndex','EdIndex','StYr','StMt',...
'EdYr','EdMt','Duration','Severity','Intensity'};
ds_event_properties_table = array2table(ds_event_properties,...
'VariableNames',table_name);
else % 如果灾害事件为0，即没有发生灾害，灾害事件属性设置为-99
ds_event_properties = nan(1,10);
ds_event_properties(1,1) = 0;
ds_event_properties(1,2) = -99;
ds_event_properties(1,3) = -99;
ds_event_properties(1,4) = -99;
ds_event_properties(1,5) = -99;
ds_event_properties(1,6) = -99;
ds_event_properties(1,7) = -99;
ds_event_properties(1,8) = -99;
ds_event_properties(1,9) = -99;
ds_event_properties(1,10) = -99;
table_name = {'EventCount','StIndex','EdIndex','StYr','StMt',...
'EdYr','EdMt','Duration','Severity','Intensity'};
ds_event_properties_table = array2table(ds_event_properties,...
'VariableNames',table_name);
end
end
2.DisasterEventExtra.m函数测试
通过构建测试数据集，编程实现对DisasterEventExtra.m灾害事件提取函数的测试。函数测试代码如下：
%{
% 1. 程序目的
(1) 学习如何从灾害指数时间序列中提取灾害事件
% 2. 版本
2.1 山东青岛 2021年3月15日
Version 1
%}
%%
clc
clear all
%% 测试数据构建
ds_array = [2019,1,-0.700000000000000;2019,2,-1;2019,3,-0.700000000000000;2019,4,1.20000000000000;2019,5,-1;2019,6,-0.700000000000000;2019,7,1;2019,8,-1;2019,9,0.500000000000000;2019,10,1;2019,11,1.20000000000000;2019,12,-0.200000000000000;2020,1,-0.800000000000000;2020,2,1.20000000000000;2020,3,-0.900000000000000;2020,4,-0.950000000000000;2020,5,1.25000000000000;2020,6,-0.900000000000000;2020,7,-1.20000000000000;2020,8,-1.50000000000000;2020,9,-0.900000000000000;2020,10,0.800000000000000;2020,11,1.10000000000000;2020,12,-0.800000000000000];
% 灾害指数阈值
threshold = -0.5;
%% 测试灾害事件提取函数
ds_event_properties = DisasterEventExtra(ds_array,threshold); % 调用函数
disp(ds_event_properties) % 查看灾害事件提取结果
%%
disp('Finished')
测试程序的运行结果如下：
三、总结
- 上述基于游程理论识别灾害时间的函数针对的时间序列是年、月时间序列，如果是年、旬时间序列同样适用，但是对于多长持续时间算作一次灾害事件可能就需要修改。最简单的修改思路是在提取出灾害事件之后，从灾害事件中选择持续时间符合条件的灾害事件。
- 上述程序没有考虑当灾害事件间隔为1个时间单位的情况下事件的合并问题，后续需要改进。
- 游程理论的实现算法上是否有利于程序的修改和完善尚有待进一步考虑和优化。
版权声明：本文为EWBA_GIS_RS_ER原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接和本声明。