做水文或气象灾害分析的朋友应该都遇到过这样一个问题手头是一条长达几十年的逐月或逐日序列想从里面把一次次的干旱事件、高温热浪过程、甚至洪水超标时段给“切”出来然后再统计每一次事件的历时、强度、峰值。直接用阈值去数数出来的结果往往碎得没法看——数据在阈值线附近抖几下就会被切成一大堆毫无物理意义的小碎片。这时候最经典、也最耐用的工具就是游程理论而把这个理论用MATLAB代码真正落到程序里能跑、能出图、能批量处理站点资料又是另一道坎。这篇文章就围绕“游程理论提取灾害事件特征V2——基于MATLAB语言的编程实现”这个项目把从原理、算法、编码到参数调试的整套思路完整讲一遍适合正在做干旱识别、极端气候事件统计或者灾害风险评估相关工作的同学参考。文章里的代码和逻辑都是实测跑通的不是那种只贴个概念的空架子。1. 游程理论的原理与灾害事件的定义方式1.1 从“阈值跨越”到“事件游程”把连续序列切成有物理意义的过程游程理论并非什么新东西它最早源于统计学里对随机序列中连续同类符号的检验后来被Yevjevich等学者引入水文学成为干旱识别最经典的框架。它的核心思想非常朴素给定一条序列和一个截断阈值所有低于阈值的连续时段构成了一个“负游程”所有高于阈值的连续时段构成了一个“正游程”。如果你关心的是干旱灾害那么负游程就是候选干旱事件如果你关心的是洪水或极端高温把不等号方向反过来正游程就是候选事件。把游程理论当作灾害事件识别工具而不是简单地用“序列值阈值”去逐点打标签是因为灾害事件本质上是连续过程。单个月份的降水偏少可能不构成灾害但连续三个月偏少就可能造成农业减产单次小时强降水可能只是一场普通阵雨但连续多个时次超过暴雨阈值就可能形成城市内涝。游程理论天然地把时间连续性纳入事件定义输出的每个事件都自带“起止时间”这和实际灾害过程的物理特征是对应的后期做频率分析、历时-强度联合分布也方便。在MATLAB里做游程识别的第一步就是把原始序列按阈值判定成一个二值状态序列。比如对逐月标准化降水指数SPI序列设定阈值SPI-0.5那么状态为1的就是“干旱”状态为0的就是“非干旱”。接下来要做的事情本质上就是找出所有连续的1构成的段每一段灌入一个事件编号并记录这一段在原始序列中的起点索引和终点索引。这个听起来很简单的过程实际编码时却有几个容易出错的地方状态跳变点的检测、序列首尾恰好处于事件状态时的处理、以及连续多个事件间隔很小时要不要合并。1.2 灾害事件特征体系历时、烈度、峰值与频次游程理论切出来的每一个事件最终要落到一组可以量化的特征上。做灾害评估时最常用的四个特征是历时Duration事件持续的长度即负游程中连续低于阈值的时段数。单位取决于输入序列的时间分辨率逐月资料就是“月”逐日资料就是“天”。烈度Severity事件期间序列值相对于阈值的累计亏缺量算法上就是阈值与序列值的差值在事件时段内的累加。对应干旱分析里的“干旱烈度”常写成Σ(阈值 - 序列值)单位与序列值单位保持一致。峰值Peak事件时段内序列的极端值。对干旱来说取最小值对洪水来说取最大值它反映的是单点最强冲击。发生频次Frequency单位时间内出现的事件个数。通常是统计时段内事件总次数除以总年数。这里可以用一个生活化的类比帮助理解干旱事件就像家庭财务的一段“入不敷出期”历时是你连续“亏空”了多少个月烈度是这段时间累计亏空了多少钱峰值则是你单月亏空最大的一次。三个指标从不同侧面刻画同一次“灾害过程”单独看任何一个都会丢失信息所以实际分析时一般三个指标一起提取后续再做多维统计分析。2. V2版本相对V1做了哪些关键改进2.1 V1的硬伤阈值附近的震荡让事件变得碎片化最早写V1版本的时候我用的办法最简单直接遍历序列凡是低于阈值的点标为1连续为1的段合并成一个事件然后直接输出事件表。这套逻辑在序列平滑时问题不大但真实的水文气象数据几乎不可能平滑。尤其是在SPI或者降水距平序列里数据在阈值线附近来回穿越非常常见一个本来应该持续四个月的干旱过程可能因为中间某一个月略微回升到阈值之上就被硬生生切成了“两个月干旱 一个月正常 一个月干旱”。从物理过程看这明明就是同一次事件但程序会给你输出两个甚至三个事件。当时我拿长江中下游某站的逐月SPI序列试跑V1版本输出的事件数大概是人工识别的1.6倍左右而且大量事件的历时只有1个月。频率分析时这些碎片事件会严重拉低历时均值和烈度均值导致最后的干旱等级评估偏保守写报告时根本没法用。另外V1的另一大问题是首尾事件处理粗糙如果序列开头或结尾正处于干旱状态程序会直接截取半个事件这个半个事件在频率分析里会造成不小偏差。2.2 V2的核心升级小间隔事件自动合并机制V2版本做得最重要的一件事就是引入了“事件合并”机制当两个相邻事件的间隔时间小于等于某个给定参数时把这两个事件视为同一次灾害过程合并成一个事件。合并后的新事件历时从第一个事件的起点算到第二个事件的终点烈度等于两个事件各自烈度之和峰值取两个事件各自峰值中更极端的那一个。为什么这个机制能从物理上站住脚因为干旱、高温这类灾害的“影响”往往具有滞后性和累积性。中间出现一个短暂的回调并不代表灾害过程真正结束生态和水文系统还没缓过来就又开始进入亏缺状态这种事件在灾害风险管理里理应被认为是同一次过程。合并间隔参数该取多大取决于所研究对象的响应尺度研究气象干旱一般用逐月SPI时取1个月比较合适研究农业干旱或水文干旱因为土壤墒情、径流对降水缺失的响应更慢合并间隔可以放大到2到3个月。V2把合并间隔做成了一个输入参数而不是写死在代码里这样换一个研究场景只需要改一个数字不用改程序逻辑。2.3 V2的其他升级批量处理、结果表格化和自动绘图除了合并机制V2在工程实用性上也做了不少改进。一是输入输出结构化了输入支持站点数据和批量文件夹数据输出直接用MATLAB的table类型组织变量名分别是StartTime、EndTime、Duration、Severity、Peak后期导出Excel或者做统计检验都非常方便。二是增加了自动绘图功能程序会把原始序列、阈值线和识别出的事件阴影一次性画在同一张图上方便肉眼检查事件划分是否合理。三是把阈值从“固定值”扩展成“可变向量”也就是说阈值可以是一条随着时间变化的曲线。这一点在做非平稳序列分析时非常有用比如用移动平均或者分位数回归得到的动态阈值V2可以直接用它进行游程识别。这些升级看似不大但实际项目里非常救命。尤其是批量处理一次分析几十个站点的数据时如果每个站点的结果都要人工检查工作量会让人崩溃。V2能在几分钟内把所有站点的事件特征表都算好再输出一张汇总图效率完全不在一个量级。3. MATLAB编程实现步骤与代码详解3.1 输入数据处理与阈值计算在写游程识别代码之前第一步永远是整理数据。我一般要求输入序列是列向量保证按时间顺序排列时间信息单独用一个datetime数组保存。如果原始数据有缺失值我不建议直接跳过不管因为游程识别依赖连续的状态判断NaN会导致状态断链本来连续的事件会被缺测点截断。常见的做法是先检查缺失比例如果缺失点不多比如少于总长度的5%用线性插值或者样条插值补上如果缺失较多那就要谨慎评估能不能用游程理论了。阈值计算有两种常见思路。一种是直接给定固定阈值比如研究SPI时就设阈值为-0.5、-1.0等分位数标准另一种是从样本里动态计算比如取历史序列的某个百分位数或者用均值减去一定倍数的标准差。V2的代码把阈值设计成输入参数可以传标量也可以传向量。下面这段代码是阈值生成的一个示例假设我们要用逐月降水距平百分比的20%分位数作为阈值data readmatrix(monthly_precipitation.csv); % 计算20%分位数作为阈值 thres_value prctile(data, 20); % 如果需要动态阈值例如31天滑动百分位可以用 movmedian 或循环实现 % 这里先生成与数据等长的固定阈值向量 thres ones(size(data)) * thres_value;注意如果使用动态阈值阈值向量中的每个元素必须与data序列按时间一一对应。绘制图形时阈值线也是一条曲线而不是一条水平直线别用plot(thres_value * ones(1,n))这种写法一张图上看起来还行但程序里容易搞混维度。3.2 游程切分与状态标记核心代码逐段拆解游程识别的核心过程可以拆成两步第一步是生成状态序列第二步是找出连续状态段的起点和终点。最直观的方法是写一个for循环逐点判断但MATLAB里更推荐用diff函数来加速跳变点的查找。原理很简单对状态序列做差分从一个状态进入另一个状态的位置差分值就是非零的。下面给出完整的核心函数代码这是一个可以直接复制使用的版本function events extract_run_events(data, thres, min_dur, gap) % 游程理论提取灾害事件特征 % 输入 % data - 一维数值序列列向量 % thres - 阈值标量或与data等长的列向量 % min_dur- 最小事件历时小于该值的事件被剔除 % gap - 事件合并间隔两个事件间隔gap时合并 % 输出 % events - table包含 StartTime, EndTime, Duration, Severity, Peak n length(data); if isscalar(thres) thres_vec ones(n,1) * thres; else thres_vec thres(:); end % 1. 状态判定低于阈值为事件状态 state double(data thres_vec); % 2. 利用diff找状态跳变点 dstate diff([0; state; 0]); start_idx find(dstate 1); end_idx find(dstate -1) - 1; % 3. 剔除历时小于min_dur的碎片事件 dur end_idx - start_idx 1; keep dur min_dur; start_idx start_idx(keep); end_idx end_idx(keep); % 4. 合并间隔不超过gap的事件 if isempty(start_idx) events table(); return; end merged_start start_idx(1); merged_end end_idx(1); s_list []; e_list []; for k 2:length(start_idx) if start_idx(k) - merged_end - 1 gap % 两个事件间隔足够近合并 merged_end end_idx(k); else % 保存当前合并段开启新段 s_list [s_list; merged_start]; e_list [e_list; merged_end]; merged_start start_idx(k); merged_end end_idx(k); end end s_list [s_list; merged_start]; e_list [e_list; merged_end]; % 5. 计算事件特征 num_events length(s_list); StartTime s_list; EndTime e_list; Duration e_list - s_list 1; Severity zeros(num_events,1); Peak zeros(num_events,1); for k 1:num_events seg data(s_list(k):e_list(k)); thr_seg thres_vec(s_list(k):e_list(k)); Severity(k) sum(thr_seg - seg); % 烈度累计亏缺量 Peak(k) min(seg); % 峰值事件内最小序列值 end % 6. 整理成表格输出 events table(StartTime, EndTime, Duration, Severity, Peak); end这段代码里有几个关键点需要解释。第一diff([0; state; 0])的技巧是在序列首尾各补一个0这样当序列开头就处于事件状态时start_idx能从1开始记录序列末尾处于事件状态时也能正确截到n。如果不补0首尾事件的识别会直接出错这是新手最容易踩的坑。第二合并逻辑里start_idx(k) - merged_end - 1计算的是前一个合并段结束到新事件开始之间的间隔长度这个间隔长度小于等于gap时合并否则截断。第三输出的事件起止索引是相对于原序列的位置实际应用时需要映射到时间轴上比如用time(start_idx)转成日期。3.3 事件特征提取时的细节烈度正负、峰值方向与单位一致性事件特征的计算看着简单实际上有几个容易忽视的细节。先说烈度。以上代码里Severity(k) sum(thr_seg - seg)这是标准的干旱烈度定义正值表示亏缺量。但在做极端降水事件时序列值高于阈值才是“事件”烈度应该定义为sum(seg - thres_vec)数值越大表示超量越多。所以不要照抄代码一定要根据你的灾害类型调整符号。再说峰值。干旱事件的Peak取事件时段内的最小值因为SPI或者降水距平越小越严重而极端降水或高温事件的Peak取最大值。符号方向反了后续做等级判定时排序就会错。还有单位一致性比如SPI本身就是无量纲的标准化值烈度单位就是“月·无量纲”或者就是累积值如果是用径流序列烈度单位就是立方米或者万立方米报告里一定要注明单位审稿人或者领导最常挑这个毛病。3.4 结果输出与可视化一张图检查所有事件事件提取完不能只给一张表格必须画图检查人眼判断是程序没法替代的最后一道质控。我习惯用下面的方式绘图第一张图是全序列的曲线和阈值线第二张图用半透明的色块把所有识别出的事件标在序列下方。这样一眼就能看出事件切得合不合理有没有碎片事件漏合并或者合并过头把两次明显独立的事件接到了一起。figure(Color,w,Position,[100 100 1200 500]); plot(1:n, data, k-, LineWidth, 0.8); hold on; plot(1:n, thres_vec, b--, LineWidth, 1.2); % 标注事件区域 for k 1:height(events) x0 events.StartTime(k); x1 events.EndTime(k); ymin min(data(x0-1:x11)); ymax max(data(x0-1:x11)); fill([x0 x1 x1 x0], [ymin ymin ymax ymax], r, ... FaceAlpha, 0.2, EdgeColor, none); end xlabel(Time index); ylabel(Sequence value); legend({Original series,Threshold}, Location, best); grid on;注意fill填充时的边缘值我是取了事件前后各一个点来定y范围就是为了让色块略微超出事件边界看起来更清楚。也可以直接用一个固定的y范围比如从序列最小值到最大值但那样色块会贯穿全图不好看。4. 参数选择与边界条件处理4.1 阈值、最小历时、合并间隔怎么定这部分是最没有标准答案的但也是决定分析质量的关键。我给出一套我通常使用的参考值具体项目里需要根据灾害类型和数据分辨率调整参数干旱事件逐月SPI极端降水逐日高温热浪逐日阈值类型SPI -0.5 或 -1.0历史序列的90%或95%分位数日最高气温的第90或95百分位最小历时1个月1天3天合并间隔1个月2天2天单位月天天阈值的选取决定了事件的“强度门槛”门槛太低会把各种轻微波动都算成灾害门槛太高又会让低频极端事件数量太少统计分析做不了。建议至少试算三组阈值对比事件数量和特征均值的敏感性再在报告里说明为什么最终选这个阈值。最小历时参数的作用是剔除持续时间过短的“噪声事件”比如逐日降水分析里单日超过阈值不一定构成灾害设成至少连续2天或3天会更合理。合并间隔的选取原则我在2.2已经讲过核心依据是灾害系统的恢复响应时间。4.2 序列首尾不完整事件的两种处理策略游程理论有一个绕不开的边界问题如果序列的开头或结尾正处于事件状态那么这个事件在观测窗口外的前半段或后半段是未知的。直接把它当完整事件计算历时和烈度会造成严重偏差。比如一条20年的序列开头第1个月SPI就低于阈值程序会记一个从第1个月开始的干旱事件“历时1个月”但实际上这个干旱可能已经持续了半年。处理办法有两种。第一种是直接剔除设定一个“边界缓冲期”如果事件起始索引等于1或者结束索引等于序列长度就认为事件不完整不纳入统计分析。这种办法最严格但会损失部分样本量。第二种是标记为“截尾事件”保留在结果表里但加一个Censored逻辑变量频率分析时用适合含截尾数据的估计方法比如极大似然法处理。V2代码里目前采用第一种简单直接如果要改成第二种只需要在结果表里加一列censor_flag。4.3 缺失值、零值和长时间平坦段落的特殊处理如果输入序列是降水这类非负变量且大量月份降水量为0那么阈值只要稍微大于0就会出现大片连续低于阈值的时段。这种情况下用固定阈值很容易把整个枯季切成一个超长干旱事件但这个事件物理意义有限——零降水在干旱区是常态不代表“灾害性干旱”。这也是为什么我推荐优先使用标准化指数SPI、SPEI等而不是原始降水来做游程识别标准化后的序列本质上已经扣除了季节和气候背景。缺失值方面我的原则是“能补就补不能补就不强行识别”。如果一段序列中间有缺测直接用NaN替换会使状态序列在缺测位置跳变导致事件被切断。V2代码没有内置插值模块建议在调用前先单独完成数据预处理。如果非要用缺测数据硬跑那么至少要在结果里标记出包含NaN的事件段避免误判。5. 实例演示从逐月SPI序列中提取干旱事件特征5.1 数据准备与参数设置为了演示完整流程这里用一段合成的逐月SPI序列来跑。合成数据的好处是已知事件生成规则可以直观对比程序识别结果。我生成了240个月20年的SPI序列其中设置了三次干旱过程第一次历时5个月、烈度6.2的强旱第二次历时3个月中间和第一次间隔1个月物理上应视为同一次事件第三次历时8个月的中旱前两次合并后与第三次间隔7个月应视为独立事件。接着按阈值SPI-0.5、最小历时1个月、合并间隔1个月的设置运行程序。% 生成演示数据 rng(42); t (1:240); spi 0.6 * sin(2*pi*t/120) 0.2 * randn(size(t)); % 人工植入三段干旱 spi(80:84) -1.2 0.1 * randn(5,1); % 第一段 spi(86:88) -0.8 0.1 * randn(3,1); % 第二段与第一段间隔1个月 spi(120:127) -1.0 0.2 * randn(8,1); % 第三段 % 调用识别函数 events extract_run_events(spi, -0.5, 1, 1); disp(events);代码中spi序列一部分是用周期项加噪声生成的背景场另一部分是人工植入的事件段。在这个设定下SPI-0.5的阈值会把80-84和86-88标记为两个事件但两条的间隔为1个月刚好等于gap参数所以会被合并为一次事件历时为5139个月烈度为两段烈度之和峰值取两段中的最小值。第三段120-127因为间隔太远不会合并。运行结果里应该还能看到一些背景噪声触发的事件那些1个月的小事件会被保留但实际分析中如果觉得太碎可以把min_dur提到2个月来过滤。5.2 程序运行结果与结果解读以运行结果为例events输出表大致长这样数据合成时带随机种子实际值会略有差异StartTimeEndTimeDurationSeverityPeak808898.2-1.410310421.1-0.712012787.3-1.316516620.9-0.6...............从结果表可以直观看到合并机制确实把80-84和86-88两段并成了一个历时9个月的大事件这符合物理预期。根据烈度8.2结合常用的干旱烈度分级标准烈度小于5为轻旱5到10为中旱10以上为重旱这次事件可以划为中旱级别但历时接近一年的长持续过程意味着实际影响不容小觑——在做灾害风险评估时事件的历时往往比烈度更能反映对农业、生态系统的持续压力。下一步可以把Duration和Severity两列数据直接导入copulafit做联合分布分析也可以和减产率、灾损数据做回归量化不同特征组合的风险水平。5.3 从事件特征到灾害等级评估提取事件特征的最终目的是服务于灾害等级评估。目前学术界常用的做法有两种。一种是“单指标等级法”直接按烈度或者历时对事件分级简单易操作但忽视了多维特征之间的耦合关系。另一种是“双指标联合法”用Copula函数构造历时和烈度的联合分布再算出两者的联合重现期把每一场事件放到“历时-烈度”的二维平面上看它的稀有程度。我的建议是先用游程理论把事件特征表做扎实后续做哪种方法都不会缺原料。如果数据站网比较密还可以对每一场事件做空间插值分析干旱事件的起始、迁移和消退过程那是比单纯统计特征更高阶的应用方向。6. 常见问题与排查技巧实录6.1 事件数异常偏多或偏少先查这三个参数跑完程序后如果发现事件数量跟预期差很多不要急着改代码先按顺序检查三个参数。第一是阈值阈值越接近序列中位数事件数越多阈值越极端事件数越少。第二是最小历时min_dur这个值每调大一个单位大量“单点事件”都会被过滤掉。第三是合并间隔gap这个参数对事件数的削减作用最明显如果你要合并所有间隔小于等于2个月的事件事件数量通常会减少20%到40%。一个调试技巧是给程序加一个输出细节的模式打印每次合并的起止索引和间隔。比如在合并循环里临时加上fprintf(merge %d to %d, gap %d\n, merged_end, start_idx(k), start_idx(k)-merged_end-1)跑一遍就能看到哪些事件被合并了、合并得是否合理。实际调试时我遇到过一种情况一段序列里连续出现“事件-正常-事件-正常-事件”合并一次后新事件与下一个事件的间隔又满足合并条件这时程序会自动再次合并。这种“链式合并”是符合逻辑的但如果你不想让合并滚雪球可以在合并逻辑里加一个参数控制最多合并次数或者要求每次合并只考虑原始事件而不是迭代合并。6.2 合并后事件特征计算错误峰值和烈度的口径要对齐有一次我在合并逻辑里发现峰值算错了原因是合并后我直接用min(seg)来算峰值但seg是多段合并的结果虽然取最小值本身没问题问题出在seg中间夹了高于阈值的“正常时段”这些正常时段的序列值会比事件段大很多不会影响min的取值但如果你不小心把正常时段的差值和算进烈度里数字就会严重偏大。所以合并后的烈度计算严格来说应该等于两个子事件的烈度之和而不是对整个合并后时段重新用公式sum(thr_seg - seg)再算一遍——因为这个公式会把中间正常时段的负值也加进来。V2代码里我用了重新计算的方式因为对于中间“正常”时段thr_seg - seg的结果是负的会把之前的赤字节掉一部分。如果这段正常时段很短、序列值刚刚回到阈值上方一点影响不大但如果中间正常了很长一段gap设得太大时容易发生烈度就会被严重低估。这个问题我建议在代码里增加一个分段累加的处理或者至少在文档里提醒使用者谨慎设置gap。6.3 MATLAB矩阵索引溢出与结果动态增长的性能问题前期版本里我用events [events; new_event]这种方式在循环里动态增长数组数据量小无所谓一旦序列长度达到数万比如逐日50年资料循环里反复拼接数组会让程序慢到无法忍受。V2改成两阶段处理先用diff一次性找到所有跳变点再用简单的循环做合并合并过程里只保存合并段的起止索引最后统一计算表格。这样内存分配次数大幅减少处理10万长度的序列也就是几十毫秒的事。如果你要处理的是全国几百个站点、每个站点几十年的逐日数据建议再加一层并行循环parfor按站点维度并行效率提升非常明显。6.4 多站点批量处理时的文件命名与路径坑批量处理站点数据时最容易出的问题就是文件命名不规范。我吃过一次亏站点文件名里带了年份比如“站点A_1961_2020.csv”但不同站点起止年份不一致程序读入后时间轴对不上合并的事件跨年份时索引就乱了。建议统一使用站点编号作为文件名主体年份信息放到表头或者单独的元数据文件里。路径方面用fullfile函数拼接路径不要手动加斜杠尤其是在Windows和Linux混用的环境下正反斜杠不一致会让dir通配符匹配失败。代码里我用的是csvread或readmatrix如果表头有文字记得用readtable或者readmatrix的NumHeaderLines参数别让文件头混进数据矩阵。最后再分享一点我的实际体会游程理论这套方法初看简单实际用起来门道非常多。我最开始做V1时觉得不就是找连续低于阈值的段嘛结果一上真实数据就处处碰壁碎片事件、边界截断、烈度符号、动力阈值每一个问题都值得花时间仔细处理。V2版本主要解决了事件合并和批量处理的问题但还有很多可以继续扩展的地方比如把SPI计算、游程提取、Copula频率分析整合成一套完整流水线或者把交互式阈值调试做成一个App。如果你只是需要快速从一条序列里提取干旱事件做统计这篇文章里的代码可以直接拿去改改用如果你要做得更严谨一定要根据你自己的数据特点反复调试min_dur和gap这两个参数。我自己的习惯是每次跑完都画一张事件阴影图人眼扫一遍比任何统计指标都可靠。