做水旱灾害风险分析的朋友大概都有这种体会数据本身不会缺逐日降水、逐日径流、水位过程线、标准化降水指数SPI一拉就是几十年甚至上百年真正让人头大的是怎么把“灾害事件”一个不落地挑出来还要把发生的起止时间、持续多久、累积量多少、峰值多极端这些特征准确地算清楚。游程理论Run Theory就是专门干这件事的经典方法而这套基于MATLAB语言实现的V2版本是我在实际项目中反复打磨后沉淀下来的一套完整方案。它会告诉你如何从原始时间序列出发自动识别出每一个灾害事件输出一张可直接用于后续统计和绘图的事件特征表并且把“识别—合并—计算—可视化—导出”整条链路整合成一个可复用工具箱而不是一次性脚本。1. 项目背景与V2版本设计思路1.1 游程理论的核心思想与适用场景游程理论在工程水文和气象灾害分析中用了很多年原理其实一句话就能说清楚给定一个阈值把时间序列的每一个样本都标记成“超过阈值”或“低于阈值”连续处于同一种状态的样本段就是一个“游程”Run。游程就是候选事件接下来要做的只是对每个游程去计算持续时间、累积量、峰值等一系列特征。这个方法真正厉害的地方在于它把肉眼判断变成了一个完全可复现的算法流程。以前我看历史降水资料时习惯在Excel里拿眼睛一段一段扫数据短还行一旦面对100年逐日序列人工标注根本不可行而且不同人标出来的结果还不一样。用游程理论做事件提取后只要阈值确定事件清单就是唯一确定的后面的重现期分析、Mann-Kendall趋势检验、空间风险区划所有统计结果都能在这个清单基础上稳定复现。V2这套实现里我默认支持两类最常见的目标场景应用场景典型数据游程方向常见阈值形式干旱识别逐日降水、SPI、径流低于阈值降水1mm/日、SPI-1、径流距平洪水/暴雨过程小时雨量、日径流、水位高于阈值设计频率P90、警戒水位、超阈值雨量极端高温/热浪逐日气温高于阈值历史95分位温度水质异常逐日浓度监测高于/低于阈值水质标准限值从表里也能看出来游程理论本质上是个通用的事件分割器不挑领域关键是阈值选得合理、后续特征统计符合业务需要。这套代码在设计时就没有绑死“干旱”或“洪水”而是把“低值事件”和“高值事件”都做进了同一个函数切换场景只用改一个参数。1.2 V1版踩过的坑与V2版重构目标最早的V1版本其实就是一段写死在脚本里的代码能跑但用起来非常难受。换一个研究对象就得复制脚本、改阈值、改列名数据前半段和后半段阈值不一样就完全无法处理相邻两个很短的小游程明明应该合并成一次完整的过程却被程序强行拆成两个事件导致最终统计里短过程数量虚高、长过程被严重低估。V1还有一个让印象比较深的问题它用循环逐点扫描整个序列几十万个样本跑起来要好几秒到几十秒。虽然不至于跑不动但调参时每改一次阈值都要等体验非常差。而且它只输出了起止索引和历时累积量、峰值都要自己再写一遍代码去算工程量被重复劳动摊薄了不少。所以V2版本在设计时定了几个明确目标。第一核心引擎独立成一个函数参数全部从外部传入不再写死任何业务逻辑。第二支持固定阈值、向量阈值、按月份或滑动窗口计算的阈值解决汛期和非汛期阈值不同、不同站点阈值不同的问题。第三增加事件合并规则允许把间隔小于指定天数的游程合并为一个完整事件。第四输出统一为MATLAB的table格式后续写CSV、画图、做透视统计都非常省事。第五可视化单独封装事件区间用半透明色带直接标在原始曲线上方便论文出图和项目汇报。这几个目标在V2里都落地了下面从数据结构和核心算法讲起。2. 数据结构设计与核心算法拆分2.1 用逻辑掩膜定位游程从边界diff到起止索引定位所有游程的起止位置是整个算法里最基础也最关键的一步。我最早写V1时用的是for循环逐点判断代码长而且慢。后来切成矩阵思维之后性能提升非常明显先通过一次比较运算得到逻辑掩膜mask再用diff在掩膜上找“从0变1”和“从1变0”的位置分别对应游程起点和终点。mask data threshold; % 低于阈值的位置记为true b diff([0; mask; 0]); % 首尾补0是为了处理边界 startIdx find(b 1); % 0-1的位置游程开始 endIdx find(b -1) - 1; % 1-0的位置前一位游程结束这里首尾补0非常关键。如果一个事件正好从序列第一个样本就开始那么mask的第一个值是1直接diff会丢掉这个起点同理如果事件一直持续到序列最后一个样本末尾也需要补一个0让它正常闭合。补0之后再diff就能保证每一个连续段都有明确的起点和终点。得到startIdx和endIdx之后一次事件的基本区间就有了。duration等于endIdx - startIdx 1这个计算本身就是向量化操作几十万个样本的数据也不会有性能压力。真正花时间的是后面的合并和特征计算但即使那样总耗时也远低于V1的循环扫描。2.2 事件合并规则间隔多大才能算同一个事件灾害事件在原始序列里往往不是干干净净的一段而是“主过程加若干小波动”的组合。以干旱为例可能连续60天降水低于阈值中间有两三天降了一点水超过阈值之后又进入低于阈值状态。从水文学角度讲这两段应该合并成一次完整的干旱事件不能因为中间下了两三天雨就把它拆成两次。所以V2引入了一个参数mergeGap含义是“相邻两个游程之间的最大间隔”。当间隔小于等于mergeGap时把两个游程合并成一个更大的事件区间从第一个游程的起点一直到第二个游程的终点。function [newStart, newEnd] merge_runs(startIdx, endIdx, gap) newStart startIdx(1); newEnd endIdx(1); for k 2:numel(startIdx) if startIdx(k) - newEnd(end) - 1 gap newEnd(end) endIdx(k); % 扩展当前事件终点 else newStart(end1, 1) startIdx(k); newEnd(end1, 1) endIdx(k); end end end这个合并函数写起来简单但设计时要想清楚一个细节合并只会影响起止索引不能在中途把两个游程的简单长度相加因为中间那段超过阈值的间隔其实也属于事件过程的“复发间隔”。正确的做法是合并完成后再用新的起止索引去原始序列里重新截取数据重新计算累积量和峰值。这个逻辑我在核心函数里做了严格分离避免出现“只合并区间但不重算特征”的隐性bug。mergeGap取多少完全看业务需求。识别干旱事件时逐日序列建议先试3天到5天识别暴雨过程时因为降水过程连续性更强一般可以放宽到7天左右。后面第5章我会专门讲参数敏感性怎么检查。2.3 特征指标怎么选历时、累积量、峰值和平均强度事件识别出来之后特征指标体系就决定了你的分析能上升到什么层次。V2版本每个事件最终都输出六个字段前四个是核心特征后两个是派生特征。核心特征直接决定事件的“几何形状”派生特征则用于事件严重程度的横向比较。特征字段含义低值事件干旱计算方式高值事件洪水计算方式StartIndex事件起点样本索引游程起点游程起点EndIndex事件终点样本索引游程终点游程终点Duration事件历时单位与数据一致EndIndex - StartIndex 1EndIndex - StartIndex 1PeakValue事件峰值/极值段内最小值段内最大值Accumulation事件累积量sum(阈值 - 原始值)sum(原始值 - 阈值)MeanIntensity平均强度Accumulation / DurationAccumulation / Duration为什么要单独区分低值和高值的计算方法因为干旱的“强度”体现在缺水量也就是阈值减去实际值的累积而洪水的“强度”体现在超限量也就是实际值减去阈值的累积。如果统一用原始值累加两个方向的物理意义都会乱套。MeanIntensity这个派生字段是我在V2里新加的。过去只报累计量和历时不好直接比较一个10天干旱和另一个30天干旱谁更严重因为历时和累计量互相干扰。平均强度把总量折算成每个样本的“缺多少/超多少”虽然只是简单除法但在做多个事件的严重程度排序时非常直观。3. MATLAB代码实现与部署细节3.1 核心函数run_events的完整实现整个V2工具箱的核心是一个单独的函数run_events.m它承担从原始序列到事件特征表的全部逻辑。我尽量把函数写得“参数够用但不臃肿”所有控制开关都通过opts结构体传入不依赖全局变量。function T run_events(data, threshold, opts) % RUN_EVENTS 基于游程理论提取灾害事件特征 % 输入 % data - 时间序列向量支持列向量或行向量 % threshold - 阈值标量或与data等长的向量 % opts.type - low提取低于阈值事件干旱high提取高于阈值事件洪水 % opts.minDur - 最小事件历时默认1小于该值的游程被过滤 % opts.mergeGap - 事件合并间隔默认0表示不合并 % 输出 % T - table类型事件特征表 arguments data (:,1) double threshold (:,1) double opts.type (1,:) char low opts.minDur (1,1) double 1 opts.mergeGap (1,1) double 0 end data data(:); threshold threshold(:); if isscalar(threshold) threshold repmat(threshold, length(data), 1); elseif length(threshold) ~ length(data) error(threshold长度必须为1或与data等长); end switch lower(opts.type) case low mask data threshold; case high mask data threshold; otherwise error(opts.type 必须为 low 或 high); end mask mask ~isnan(data); % 缺测值不参与事件识别 b diff([0; mask; 0]); startIdx find(b 1); endIdx find(b -1) - 1; if opts.mergeGap 0 ~isempty(startIdx) [startIdx, endIdx] merge_runs(startIdx, endIdx, opts.mergeGap); end n numel(startIdx); if n 0 T table(); return; end duration endIdx - startIdx 1; peak zeros(n, 1); acc zeros(n, 1); for k 1:n seg data(startIdx(k):endIdx(k)); thrSeg threshold(startIdx(k):endIdx(k)); if strcmpi(opts.type, low) acc(k) sum(thrSeg - seg); peak(k) min(seg); else acc(k) sum(seg - thrSeg); peak(k) max(seg); end end keep duration opts.minDur; startIdx startIdx(keep); endIdx endIdx(keep); duration duration(keep); peak peak(keep); acc acc(keep); T table(startIdx, endIdx, duration, peak, acc, ... VariableNames, {StartIndex,EndIndex,Duration, ... PeakValue,Accumulation}); T.MeanIntensity T.Accumulation ./ T.Duration; end function [newStart, newEnd] merge_runs(startIdx, endIdx, gap) newStart startIdx(1); newEnd endIdx(1); for k 2:numel(startIdx) if startIdx(k) - newEnd(end) - 1 gap newEnd(end) endIdx(k); else newStart(end1, 1) startIdx(k); newEnd(end1, 1) endIdx(k); end end end这里有几个实现细节值得展开说。第一threshold既支持标量也支持向量向量阈值意味着可以逐日使用不同阈值这是处理季节差异的关键。第二mask最后强制去掉NaN位置否则NaN会被当作低于阈值凭空生成一堆假事件。第三特征计算里的for循环不会成为性能瓶颈因为循环次数等于事件数量而不是样本数量即使100年逐日序列里有两三百个事件这个循环也是瞬间完成。第四arguments语法需要MATLAB R2019b及以上版本如果还在用R2018a及更早版本需要改成nargin和narginchk的写法我这里为了代码简洁直接用了新语法。3.2 阈值生成固定阈值、百分位阈值与按月阈值阈值是整个游程理论里最重要的参数选得好不好直接决定事件清单是否合理。V2把阈值计算和事件识别拆开了你可以先用任何方式算出阈值向量再传给run_events。我这里分享三种最常用的生成方式都可以用在实践中。固定阈值是最简单的一种直接一个标量传给函数就行。比如逐日降水干旱定义的阈值常取1mm/日SPI定义的干旱阈值常取-1。这种方式适合已有明确行业标准或设计标准的场景优点是结果可解释性强缺点是不同地区相同阈值可能会误判。百分位阈值适合没有明确标准、但数据长度足够的情况。比如逐日降水可以用历史90分位作为暴雨阈值用10分位作为干旱阈值。MATLAB里直接用prctile就能算thr repmat(prctile(data, 10), length(data), 1); % 固定10分位阈值按月阈值解决的是季节差异问题尤其适合中国这种降水高度集中、冬夏差异显著的地区。同一个绝对阈值在汛期和枯水期完全没有可比性因此按每个月份的多年分位值分别计算阈值更合理month month(t); % 根据时间轴提取月份编号 thrMonth accumarray(month, data, [], (x) prctile(x, 10)); thr thrMonth(month); % 展开为与data等长的向量如果数据量不够、按月阈值不稳定也可以用滑动窗口的百分位阈值。窗口长度通常取30天到90天但我个人建议慎重使用滑动窗口因为它会让阈值本身变得非常平滑事件识别的边界会被“磨”掉不少。当站点资料只有二三十年时按月阈值比滑动窗口更稳定、更接近气候态。3.3 可视化与结果导出事件区间标注和图片输出事件识别出来之后如果没有可视化的辅助检查很难让人完全放心。尤其是拿给别人看结果的时候一张带事件色带标注的过程线比一堆数字有说服力得多。V2里我习惯用patch画半透明色带把每个事件区间在原始曲线上框出来。figure(Color,w,Position,[100 100 1200 400]); plot(t, spi, k-, LineWidth, 1); hold on; plot(t, thr, r--, LineWidth, 1.2); yl ylim; for k 1:height(T) xStart t(T.StartIndex(k)); xEnd t(T.EndIndex(k)); patch([xStart xEnd xEnd xStart], [yl(1) yl(1) yl(2) yl(2)], ... [0.8 0.9 1.0], FaceAlpha, 0.4, EdgeColor, none); end xlabel(时间); ylabel(SPI); legend({SPI,阈值,事件区间}, Location,best); box on; grid on;这段代码里最值得说的是patch的用法。它的四个顶点坐标是左右边界和时间轴上下限配合FaceAlpha设置透明度就能做出“事件区间变亮”的效果同时不遮挡原始曲线。有人会问为什么不用area或者fillarea在x轴上有基线不适合直接标注一段背景fill会被多条曲线干扰图层顺序patch配合FaceAlpha是控制背景标注最灵活的方式。图片导出方面我现在基本弃用print全部改用exportgraphics因为它在分辨率控制和字体渲染上更稳。要出论文图就300dpi起步要矢量图就导出epsexportgraphics(gcf, drought_events.png, Resolution, 300); exportgraphics(gcf, drought_events.eps, ContentType, vector);结果表格导出CSV用writetable一行搞定writetable(T, drought_events.csv);如果需要把索引换算成实际时间只要提前把t(T.StartIndex)提取出来加成一列即可。4. 完整算例用模拟SPI序列识别干旱事件4.1 构造模拟序列并运行识别为了让整个流程看得见摸得着我用模拟的逐月SPI序列跑一遍完整算例。SPI本身就是标准化指数理论上近似标准正态分布所以这里用randn生成模拟序列是合理的。我生成80年逐月数据识别“SPI小于-1”的干旱事件要求最小历时3个月、间隔不超过2个月的相邻事件合并。rng(2024); n 12 * 80; t datetime(1944, 1, 1) calmonths(0:n-1); spi randn(n, 1); opts.type low; opts.minDur 3; opts.mergeGap 2; T run_events(spi, -1, opts); disp(head(T, 10));运行后的事件特征表示例长这样具体数值随随机种子不同会有差异StartIndexEndIndexDurationPeakValueAccumulationMeanIntensity31366-1.231.870.3152587-1.613.420.4989924-1.081.120.281211288-1.954.760.601771815-1.341.680.34从表格可以直观看到mergeGap2把很多原本零散的小游程合并成了持续几个月的事件。Duration列就是合并后的完整历时Accumulation列表示累计缺水量PeakValue列记录的是事件期间SPI最低值也就是最干旱的月份。4.2 识别结果解读与敏感性检查拿到事件表之后第一步不是急着做统计而是先做三件事。第一算一下所有事件的历时分布看是不是集中在3到8个月如果出现大量历时长达三四十个月的事件大概率是合并间隔给太大了。第二把事件区间画到原始SPI曲线上肉眼检查每个色带是否和实际低于阈值的片段吻合。第三画一张事件开始时间与月份的关系图验证干旱事件是否多发于某个季节。敏感性检查也很重要我建议固定其他参数只改变mergeGap分别取0、1、2、3、5观察事件数量的变化。一次合理的参数选择应该满足事件数量随mergeGap增大而先快速下降之后趋于平缓。如果mergeGap从2改到3时事件数量还在剧烈变化说明原来的参数太敏感需要谨慎取值。用模拟数据时有一个好处结果可以反复生成颗粒度很细。但换成真实站点数据后参数敏感性检查会直接决定研究结论会不会被审稿人质疑。我遇到过不少论文一上来就用某个默认间隔跑完全部统计完全不说明为什么取这个值这种漏洞在审查时很容易被抓住。V2把mergeGap、minDur都做成了显式参数就是希望每一步都能留痕、可解释。5. 常见问题排查与参数调优实录5.1 缺测值和边界效应对识别的干扰真实观测数据几乎不可能没有缺测尤其降水数据逐日序列里经常有NaN。代码里我对NaN的处理是直接不参与识别也就是mask里对应位置强制为false但这会带来一个副作用如果一段缺测正好发生在两个游程中间它会把本应连续的游程硬生生截断导致mergeGap必须设得很大才能把它们重新接上而mergeGap设大了又会误合并其他事件。这个问题没有完美解法我的习惯是分情况处理。如果缺测占比很小少于1%建议先用插值把序列补全再识别比如用spline插值或者相邻月份平均值。如果缺测占比大且集中最好不要强行插值而是把缺测段作为事件断点来处理在报告中明确说明识别结果会低估长事件。还有一种做法是把缺测段也视为低于阈值但这只在极端情况下使用因为它会让事件长度显著偏大必须谨慎。边界效应是所有序列分析方法都绕不开的问题。序列开头如果正处在一次事件中间那么这次事件从一开始就不完整序列结尾同理。V2代码通过首尾补0解决了“识别不闭合”的问题但它不能解决“事件被截断”的物理问题。在最终统计时我通常会把序列首尾各去掉半年或者单独标记事件是否与边界相邻避免边界事件干扰频率统计。5.2 合并间隔和最小历时怎么取这两个参数是用户最容易拍脑袋定的地方。minDur的物理含义是“多短的过程不算事件”。对SPI干旱分析小于3个月的事件一般不算气候意义上的干旱而更像短时异常所以minDur取3比较常见。对逐日降水如果只关注暴雨过程minDur可以取1只要单日超阈值就算一次事件如果关注洪水过程则建议minDur取2到3把短时孤立降水过滤掉。mergeGap的取值建议遵循“过程连续性”原则。降水过程里两场雨间隔超过7天基本可以认为是两次独立过程所以mergeGap最多取7。干旱事件里面中间偶尔一两天的降水不改变干旱本质所以逐日序列取3到5是合理区间。这些值不是拍脑袋出来的应该和业务专家确认“多少次小波动内仍然算同一次事件”再把确认结果翻译成参数。我做项目时还有个习惯就是做一张参数敏感性表放在分析文档里。列是mergeGap行是minDur中间的数值是事件总数这样评审人看到后心里会踏实很多。V2的代码跑一次只要几十毫秒遍历参数组合完全没压力。5.3 代码性能优化与版本兼容性性能方面V2相比V1最大的提升就是去掉了全序列循环。定位游程用逻辑索引和diff特征计算只在事件个数级别循环因此即使面对50000个样本的逐日百年序列一次完整运行也基本在0.01秒量级。如果序列更长比如测站很多、需要批量跑几千个格点还可以做进一步的批量并行化用parfor遍历站点每个站点内部还是调用同一个run_events函数。版本兼容性上面提过主要注意两点。第一arguments语法需要R2019b及以上如果你的环境是R2018a及以前需要把参数校验改成narginchk加属性判断。第二exportgraphics需要R2020a及以上老版本可以用print但分辨率控制不如exportgraphics直观。我实测过R2021a和R2023b环境代码都不需要改动核心逻辑只依赖MATLAB基础模块不调用任何工具箱所以即使只有最基础的MATLAB版本也能跑。我个人在实际使用中的体会是工具越通用越要把参数边界定义清楚。V2这套代码之所以比第一版顺手不是因为原理变了而是因为每个开关都做了明确的输入输出约束不再需要每次复制脚本后去改逻辑关系。如果你也要在项目中应用游程理论我的建议是先用模拟数据把函数跑通再用真实数据做参数敏感性检查最后才进入正式统计环节。编码本身只是把成熟的思路落地真正花心力的是理解和验证你设定的每一个参数。