游程理论这名字听起来像数学教材里的冷门概念但在灾害特征提取里它是最朴素也最耐用的一把刀。搞过干旱评估、洪峰识别、极端温度过程分析的人都绕不开它——把连续时间序列按阈值切成一段段低于阈值和高于阈值的区间每一段游程就是一个潜在的事件。就这么个思路从1967年Yevjevich提出到现在水文气候领域一直在用连发高水平论文都够用。但问题在于理论好懂写代码是另一回事。你手上有十年逐日降水、径流或SPI指数序列怎么自动把所有干旱事件捞出来计算出每场事件的起止时间、持续时间、累计烈度、峰值强度还要处理合并、滤噪、缺测一套流程跑完不出错——这个工程量比你想象中大得多。这篇博文分享的是我用MATLAB实现的第二版框架主要解决V1版里参数写死、阈值单一、事件合并依赖人工干预、大数据量循环效率低这几块痛点完整代码结构、核心函数、参数敏感性分析经验都会摊开讲适合做水文分析、气候诊断、灾害风险评估的学生和同行直接参考。1. 游程理论的数学含义与灾害事件的定义约定1.1 正负游程的本质阈值切割后的连续时间片段游程理论不复杂。给一个时间序列 \(x_t\) 和某个截断水平 \(x_0\)把所有时间点标记为两类\(x_t x_0\) 和 \(x_t \ge x_0\)。连续处于同一状态的时段就是一个游程。低于阈值的那一段称为负游程或亏水量游程高于阈值称为正游程。在干旱研究里我们关心负游程在洪水、极端降水研究里通常以超过阈值的正游程作为事件主体。关键是连续两个字。序列里哪怕只隔一个点高于阈值前后两段负游程就被切开了。这个特性直接决定了后续要不要做事件合并——真实世界不会因为你中间有两天下雨就把前后两场干旱硬拆成两段所以算法层面必须有合并机制。数学表达也不绕状态序列\(s_t I(x_t x_0)\)游程长度\(L \sum_{ta}^{b} s_t\)其中 \([a,b]\) 是连续状态区间游程强度累计亏水量\(D \sum_{ta}^{b} (x_0 - x_t)\)游程峰值\(P \max(x_0 - x_t)\)MATLAB里实现状态切换检测我习惯用diff配合逻辑数组不必循环逐点判断。比如below data threshold; change_idx find(diff([0; below(:); 0]) ~ 0);这样change_idx里就存了所有状态切换位置奇数位置是游程起点偶数位置是终点一次向量化操作拿到全部游程区间。这个思路是整套代码的地基V1版当时用的是for循环逐点判断数据一长十几万步就明显卡顿V2全面改成这种边界点定位法。1.2 阈值是常数还是动态的这决定了事件定义是否合理很多人做游程分析时只想着找个固定阈值一刀切但灾害事件的阈值往往不该是常数。以干旱为例不同季节、不同地区的降水基准差异巨大一套固定阈值要么在湿润区什么都识别不出来要么在干旱区把所有时段都判定为事件。V2版支持两种阈值模式绝对阈值直接指定一个常数例如径流小于5 m³/s视为枯水。动态阈值按时间窗口滑动计算百分位数比如逐日SPI序列低于-0.8或者径流低于过去30年同期30%分位值。动态阈值的实现核心在于窗口对齐。比如逐月数据要算每月的多年平均或分位数不能错位用movmedian或自定义分组函数都比较方便。代码上我封装了一个getThreshold函数输入数据和模式参数输出与数据等长的阈值数组function th getThreshold(data, time, mode, win) % mode: absolute, moving_percentile, monthly_clim switch mode case absolute th repmat(win, size(data)); case moving_percentile th movpercentile(data, win, 20); % 20%分位 case monthly_clim th monthlyPercentile(data, time, 20); end end动态阈值有一个容易踩的坑滑动窗口的边界段会出现NaN或样本不足的问题。MATLAB的movmedian、movprctile需要统计或金融工具箱在窗口不满时会默认用Endpoints参数填充我一般设为fill并手动补最小值避免边界处误判出大片事件。1.3 事件特征向量除了起止时间你还需要算什么游程识别只是第一步提取灾害事件特征才是用户真正关心的结果。做灾害分析时事件特征向量通常包括这些量特征名称计算公式物理含义事件编号按时间顺序计数唯一标识起始时间游程起点对应的日期事件开始结束时间游程终点对应的日期事件结束持续时间长度L单位为月/天事件持久性峰值强度min(x_t)或max(threshold - x_t)最严重时刻累计强度sum(threshold - x_t)总亏水量/总超阈值量平均强度累计强度/持续时间事件平均严重程度峰值出现时间最小值所在时刻事件演变的关键节点这些特征量组合起来足以支撑后续的趋势分析、概率分布拟合、时空聚类、事件频次变化等研究。V2版在输出特征矩阵的同时也会生成一张事件时间进度条图纵轴每行是一个事件横向是发生时间范围这对快速查看识别的合理性特别有用。2. V2版程序框架设计从一次性脚本到可配置工具箱2.1 V1版的痛点一条产品线里三个脚本改了又改我第一版的时候是接到一个任务——分析某个流域近六十年逐日径流的枯水事件。当时图省事写了一个主脚本阈值固定、事件合并规则硬编码每次换数据都要打开.m文件改参数换一个研究区域就复制一个脚本最后桌面上多出十几个不同命名的文件。更麻烦的是一旦上游数据更新了时间段整个过程重跑一遍输出表格的格式还不统一下游画图脚本经常对不上列名。V1的另一个致命弱点是事件合并太粗暴。我当时简单地规定两个事件间隔小于3天就合并但这个3天是拍脑袋定的而且没有考虑事件本身持续时间长短——间隔3天对于持续20天的长事件可能该合并对于只持续2天的短事件则可能不该合并。这种一刀切的规则直接导致识别结果在长时间尺度上波动很大。V2的核心目标就是规范化这几个问题参数不写死在代码里事件合并规则可配置输出统一结构代码模块化到负责人能直接跑的程度。2.2 模块划分与主控流程我把整个程序拆成六个模块run_analysis.m % 主控脚本配置参数、调用各模块 load_input.m % 数据读取与预处理缺测插值、NaN处理 get_threshold.m % 阈值生成绝对/动态 detect_runs.m % 游程检测返回游程起止索引 postprocess_events.m % 事件合并、最小持续时间过滤、边界处理 extract_features.m % 特征向量计算与输出 plot_events.m % 可视化检视主控脚本长这样参数集中在开头几乎不需要翻后面代码%% 用户配置区 data_file precipitation_daily_1960_2020.csv; time_col 1; value_col 2; threshold_mode moving_percentile; window_years 30; percentile_val 20; min_duration 3; % 最小事件持续天数 merge_gap 5; % 事件合并最大间隔天数 %% 固定流程 [data, time] load_input(data_file, time_col, value_col); th get_threshold(data, time, threshold_mode, window_years, percentile_val); runs detect_runs(data, th); events postprocess_events(runs, merge_gap, min_duration); feat extract_features(events, data, th, time); plot_events(time, data, th, events);这种设计下后续换数据、换阈值方式只需要动第一段配置区其他部分完全不用碰。我自己的体会是写科研分析代码抽象程度不必太高但配置与逻辑分离这一条绝对值得坚持因为你三个月后回来看自己的脚本唯一还能快速上手的入口就是清晰的配置区。2.3 面向过程还是封装类MATLAB场景下的选择有人会问要不要用MATLAB的classdef面向对象重写我试过但最终放弃了。原因是灾害特征提取虽然逻辑不复杂但研究人员通常需要在不同数据、不同阈值、不同合并规则之间快速试出合理结果函数式模块比类封装更轻量也更符合大多数做水文气候分析的MATLAB用户的使用习惯。但如果你的项目要复用性很高比如要打包给团队其他人调用、要集成到更大的分析流程里那建议把核心函数做成独立文件加上help注释做成工具箱的形式。V2版的妥协方案是所有函数以独立.m文件存储互相之间只依赖输入输出结构体不共享全局变量。这样既保留函数式调用的灵活又能调用addpath方式打包共享。以结构体events为例字段定为events struct(... start_idx, {}, ... % 游程起点索引 end_idx, {}, ... % 游程终点索引 start_time,{}, ... end_time, {}, ... duration, {}, ... % 持续时间 peak_idx, {}, ... % 峰值所在索引 peak_value,{}, ... % 峰值强度 intensity, {}); % 累计强度模块之间的数据通过这种结构体传递后续扩展新特征量的时候不需要改函数签名只要在extract_features里加字段即可。3. 核心函数实现游程检测、事件合并与边界处理3.1 detect_runs的向量化实现逻辑detect_runs是整套代码里最核心也最容易被写慢的地方。前面提到用diff找状态切换点完整实现如下function runs detect_runs(data, threshold) % DETECT_RUNS 基于游程理论识别低于阈值/高于阈值的连续区间 % 输入: % data - 列向量时间序列 % threshold - 列向量与data等长的阈值序列 % 输出: % runs - N x 2 矩阵每行是[起始索引, 结束索引] below data threshold; % 状态切换点前后状态不同之处 change diff([0; below(:); 0]); start_idx find(change 1); end_idx find(change -1) - 1; % 如果最后一个切换是状态开始但没有结束说明序列末尾仍在游程中 if length(start_idx) length(end_idx) end_idx(end1) length(data); end runs [start_idx, end_idx]; end你可能会问为什么change用diff([0; below; 0])而不是diff(below)因为这样能把序列首尾的游程边界也抓出来——第一个点如果是belowtrue前面补0后diff会得到1表示游程从这里开始最后一个点如果仍是belowtrue末尾补0后diff会得到-1表示游程在这里结束。否则首尾的事件会被吃掉这个错误在V1版踩过V2直接杜绝。这个函数比逐点if判断快一个数量级以上。十万个时间点循环写法需要十几秒这版向量化写法在零点几秒内完成。如果数据量更大后面可以考虑用find(diff(below))的C加速路径但一般来说这个效率已经够了。3.2 事件合并为什么单纯按间隔阈值不可靠前面说V1版拍脑袋定了间隔3天就合并V2做了改进把合并改成两级判定第一级间隔小于等于merge_gap的事件先被合并。第二级如果合并后的事件持续时间仍小于min_duration则丢弃。这个顺序是刻意的。先合并再滤短可以避免一个长事件被短暂间歇切成两段后每段都过短被误删的问题。实现上合并逻辑就是循环判断两个相邻事件之间的索引间隔。由于游程数量一般几百到几千循环在这里性能损失不大所以用循环反而更直观function events mergeRuns(runs, merge_gap) % MERGERUNS 将间隔过短的相邻游程合并 if isempty(runs) events runs; return; end merged []; cur_start runs(1, 1); cur_end runs(1, 2); for i 2:size(runs, 1) gap runs(i, 1) - cur_end - 1; if gap merge_gap % 合并延伸当前事件结束位置 cur_end runs(i, 2); else merged(end1, :) [cur_start, cur_end]; %#okAGROW cur_start runs(i, 1); cur_end runs(i, 2); end end merged(end1, :) [cur_start, cur_end]; events merged; end我在这里踩过一个真实的坑这个gap的计算公式是runs(i,1) - cur_end - 1还是runs(i,1) - cur_end这要看start_idx和end_idx是否包含边界。我定义的是闭区间索引[start, end]内所有点都是游程的一部分所以中间隔着的天数等于下个起点减当前终点再减1。比如第一个事件在5号结束第二个事件从8号开始实际间隔是6号、7号两天gap 8 - 5 - 1 2。如果写成8 - 5 3就会多算一天可能触发不该发生的合并。这种小细节看代码不觉得跑数据发现事件数量不对劲排查老半天才会意识到。3.3 最小持续时间过滤与缺测段的特殊处理min_duration过滤相当简单但缺测段NaN的处理就有讲究。游程检测时数据里不能有NaN否则data threshold得到false会把一个事件从中间劈开V1版在这个问题上出过很诡异的Bug。V2版采用的做法是数据加载后先做缺测插值。日值降水如果用线性插值可能会把有雨和无雨之间造出一些虚假的中间值这在干旱事件识别里不能接受。更稳妥的做法是对缺测天数极少小于总长度的1%的序列用临近值填充对缺测段过长的数据序列应直接设置thresholdNaN并显式把这些区段排除在事件识别之外同时在结果里标记该时段数据缺失事件边界不可靠。实际实现中我在load_input里增加了一个返回字段quality_flag缺失部分标记为0在detect_runs里遇到缺测段就强制将状态置为非事件并在postprocess_events里把横跨缺测段的事件拆开。function below markBelow(data, threshold, missed) below data threshold; below(missed) false; % 缺测段强制为非事件 end这样处理之后即使数据有间断也不会污染事件识别结果。4. 关键参数对事件识别结果的影响最小持续时间与合并窗口4.1 参数敏感性同一个序列不同参数得到不同事件数我用自己的流域径流数据1960-2020逐日做了个快速参数扫描这里给一组典型结果最小持续时间(天)合并间隔(天)识别事件数平均持续时间(天)0018711.20514314.63010318.2357824.5755233.17104438.7这组数字非常直观参数变一下事件数量几乎翻倍变化。不存在标准答案参数必须根据研究目的来确定。如果关注的是极端长期干旱min_duration就该设大一些如果关注的是所有短暂缺水过程对农业的影响min_duration应设小或为0。但不管怎么设有一点必须谨记在报告或论文里必须明确写出参数值而且最好做敏感性分析说明事件特征对参数的依赖程度。否则审稿人或合作者很容易质疑结果的可靠性。4.2 周期性与游程合并的关系季节性因子不可忽视游程合并还有个季节性陷阱。如果研究区域是典型的夏季风气候干季和湿季分明那么干季的负游程天然很长不需要合并但如果在春秋过渡期降水序列会频繁跨越阈值导致一场秋季缺水过程被切成好几段。这种情况用单一的merge_gap参数很难兼顾。V2版在配置里增加了一个可选参数seasonal_merge允许按月份设置不同的合并间隔比如夏季设7天、冬季设3天。实现上就是把前面mergeRuns函数里的固定merge_gap换成一个与事件位置相关的数组gaps getMonthlyMergeGap(merge_gap_seasonal, time, runs);这个扩展我用过之后觉得很值——它避免了一个窘境参数设大了把独立事件错误拼在一起设小了又把真实的长事件拦腰截断。如果你处理的是有明显季节性的数据强烈建议加这个机制。4.3 阈值分位数选择20%和10%差异巨大阈值分位数是另一个敏感参数。同样是动态阈值用20%分位值和10%分位值识别出的干旱事件强度特征完全不同。10%分位数识别出的事件频率低、但平均强度大更适合做极端干旱研究20%分位数则更适合做整体缺水状况评估。这也带出一个方法论的提醒阈值本质上是事件定义的一部分而不是单纯的算法参数。在动手写代码之前先想清楚你的灾害类型对应什么样的重现期或分位数标准。比如很多国内干旱研究采用径流距平百分率低于-20%或SPI值低于-1.0作为阈值这些是领域内约定俗成的标准直接选用能让结果更好被同行比较。5. 实测案例全流程从逐日径流序列到干旱事件特征表5.1 案例数据与预处理为展示完整流程我用一段模拟生成的逐日径流数据2000天约5.5年来走一遍。数据包含明显的季节性周期和几段人为设置的枯水期。读者可以把这部分替换成自己的实测数据。% 模拟数据 rng(42); n 2000; time datetime(2000, 1, 1) days(0:n-1); seasonal 30 20 * sin(2 * pi * day(time, dayofyear) / 365.25); noise 8 * randn(n, 1); flow seasonal noise; flow(flow 0) 0; % 人为制造两段持久的枯水期 flow(400:460) flow(400:460) * 0.2; flow(1200:1300) flow(1200:1300) * 0.1;预处理阶段我先用fillmissing线性插值补齐可能的缺失值再设定动态阈值为过去一年滑动窗口的20%分位数。窗口设1年既保证样本量足够又能跟上季节变化。5.2 游程识别与合并的可视化检视跑完detect_runs和postprocess_events后第一件事不是直接看表格而是先可视化因为人眼是判断算法合理性的第一道防线。我画了这样一张图径流曲线灰色显示动态阈值用红色虚线识别出的枯水事件用浅蓝色阴影标出。结果中发现开头第100天附近识别出一个持续6天的事件但我肉眼看着那一段其实只是正常的季节低值并没有真正枯到异常。为什么被识别因为动态阈值恰好在高值季节那条低值段的径流仍然低于滑动20%分位。这说明动态阈值在高变异季节容易产生误报。解决方案之一是抬高阈值分位数或者改用超过特定距平的人为标准。我在代码里预留了threshold_bias参数允许用户对动态阈值增加一个绝对偏移量这种情况直接手动调2就能压制误报。5.3 特征表输出与下游分析最终特征表长这样节选事件编号开始日期结束日期持续时间(天)峰值亏水量(m³/s)累计亏水量(m³·s⁻¹·天)平均强度12000-06-122000-06-26158.2351.43.4322001-02-022001-03-053212.67168.95.2832002-12-102003-01-184018.02296.37.4142003-02-132003-02-2085.9621.82.73输出格式上我用了writetable直接写出CSV列名与特征结构体字段保持一致。下游做频率分析时直接读表即可不用再改代码。为了避免列名歧义CSV里附带一行注释说明各列的物理单位这一招在处理合作项目时帮我避免了很多沟通成本。6. 调试过程避坑记录与性能优化经验6.1 首尾游程丢失问题diff边界补齐的教训V1版我用diff(below)查找切换点时序列开头就是负游程的情况下diff结果不会出现1导致第一段事件永远不会被识别。当时给我的分析结果带来很大的偏误——最干旱的区域反而被识别出最少的事件。V2版在状态序列前后各补一个0再diff正是为了解决这个边界问题。这里再强调一次change diff([0; below(:); 0]);这行代码是整个程序的隐形功臣。你如果发现自己跑出来的事件数量异常偏少优先检查这行事件从第2个点开始的情况往往就是这里出了问题。6.2 峰值搜索在等值平台段的处理clip还是中心点计算峰值强度时如果序列中出现连续多个点数值完全相同这种情况在低分辨率数据里很常见find会返回第一个最小值这没问题。但如果你用min(data(start:end))找到的是数值还想要峰值出现的时刻就要小心等值平台段里第一个最小点可能明显偏早或偏晚从而影响你对峰值出现时间的估计。我采用的策略是取平台段的中间索引作为峰值出现时间% 找到最小值的全部位置取中间点 [min_val, ~] min(segment); min_positions find(segment min_val); if length(min_positions) 1 peak_idx round(mean(min_positions)); else peak_idx min_positions(1); end6.3 性能对比向量化 vs 循环以及memory的取舍我用一个60万步的日值序列约1643年做了测试detect_runs循环版耗时约85秒向量化版本0.3秒。差距巨大所以向量化是必选。但如果序列极度长例如5分钟间隔水文数据一年就超过10万步构建below逻辑数组和diff结果本身也会消耗内存。可以用uint8类型省空间MATLAB逻辑数组已经只占1字节所以below本身开销不大diff结果的double类型相对占内存但序列长度百万量级也就几MB可以忽略。真正的大内存大户是动态阈值数组的构建过程movpercentile在长序列上可能要反复排序内存峰值会高。这种场景下把序列分段处理段间保留一个merge_gap的overlap是更稳妥的路线。6.4 批量处理多个站点的数据组织技巧做区域分析时经常要遍历几十上百个站点的数据。不要在主脚本里叠一个巨大for循环把站点的参数数据文件路径、阈值模式、合并间隔、最小持续时间放在一个表格里用parfor并行处理每个站点独立输出特征表和诊断图。我实际跑过一个小流域42个站点parfor在四核机器上从42分钟降到14分钟节省的时间足够用来多调几轮参数。需要注意parfor里头的函数必须全部是显式输入、显式输出不能依赖全局变量和脚本环境变量这也是我把所有逻辑都独立成函数文件的原因之一。6.5 关于MATLAB版本与工具箱依赖这套代码尽量只用基础MATLAB函数核心用到diff、find、min、datetime、movmedian此函数在R2016a之后才可用。如果手头版本太老movmedian可以用circular滑动窗口自写替代。画图部分用到patch和datetime轴这些也都在基础版里。有些人问要不要用工具箱里的findpeaks或信号处理函数。我的建议是不要过度依赖工具箱尤其当你需要把代码分享给没有相应工具箱的合作者时一个工具箱函数可能就是阻碍复现的门槛。游程理论本身足够简单自己实现更有掌控感也能随时按需修改。7. 从V2到后续还有哪些扩展切入点一套游程识别程序跑通以后你会发现它适合做的事情远不止干旱识别。给径流序列加一个正阈值它就变成洪水事件提取器把输入从径流换成气温低于阈值的游程就是寒潮过程换成空气质量指数高于阈值的就是污染过程。本质都是阈值切割连续区间识别特征量计算这一套逻辑。在这个框架基础上我建议你继续往这几个方向扩展一是把识别结果直接接入概率分布拟合如GEV、Gamma分布去算重现期很多灾害风险评估项目都需要这步二是增加空间模块多个站点的事件特征可以汇合到GIS环境中做时空制图三是考虑事件强度与时间的关系加入趋势检验Mann-Kendall函数直接输出特征序列的趋势显著性结果。就个人经验而言代码库的维护最关键的不是写得多炫而是命名一致、结构稳定、配置清晰。V2版做完后我后续给同一流域做极端降水分析时只写了两个新模块一个处理正游程一个把事件过程从原始序列中截取出来做事件合成分析。核心的detect_runs和postprocess_events一行没改。这让我对通用工具四个字有了新的理解——不是写出无所不能的框架而是把最基本的环节做扎实让上层变化有稳固的落点。 基于游程理论的MATLAB干旱事件自动识别与特征提取
游程理论这名字听起来像数学教材里的冷门概念但在灾害特征提取里它是最朴素也最耐用的一把刀。搞过干旱评估、洪峰识别、极端温度过程分析的人都绕不开它——把连续时间序列按阈值切成一段段低于阈值和高于阈值的区间每一段游程就是一个潜在的事件。就这么个思路从1967年Yevjevich提出到现在水文气候领域一直在用连发高水平论文都够用。但问题在于理论好懂写代码是另一回事。你手上有十年逐日降水、径流或SPI指数序列怎么自动把所有干旱事件捞出来计算出每场事件的起止时间、持续时间、累计烈度、峰值强度还要处理合并、滤噪、缺测一套流程跑完不出错——这个工程量比你想象中大得多。这篇博文分享的是我用MATLAB实现的第二版框架主要解决V1版里参数写死、阈值单一、事件合并依赖人工干预、大数据量循环效率低这几块痛点完整代码结构、核心函数、参数敏感性分析经验都会摊开讲适合做水文分析、气候诊断、灾害风险评估的学生和同行直接参考。1. 游程理论的数学含义与灾害事件的定义约定1.1 正负游程的本质阈值切割后的连续时间片段游程理论不复杂。给一个时间序列 \(x_t\) 和某个截断水平 \(x_0\)把所有时间点标记为两类\(x_t x_0\) 和 \(x_t \ge x_0\)。连续处于同一状态的时段就是一个游程。低于阈值的那一段称为负游程或亏水量游程高于阈值称为正游程。在干旱研究里我们关心负游程在洪水、极端降水研究里通常以超过阈值的正游程作为事件主体。关键是连续两个字。序列里哪怕只隔一个点高于阈值前后两段负游程就被切开了。这个特性直接决定了后续要不要做事件合并——真实世界不会因为你中间有两天下雨就把前后两场干旱硬拆成两段所以算法层面必须有合并机制。数学表达也不绕状态序列\(s_t I(x_t x_0)\)游程长度\(L \sum_{ta}^{b} s_t\)其中 \([a,b]\) 是连续状态区间游程强度累计亏水量\(D \sum_{ta}^{b} (x_0 - x_t)\)游程峰值\(P \max(x_0 - x_t)\)MATLAB里实现状态切换检测我习惯用diff配合逻辑数组不必循环逐点判断。比如below data threshold; change_idx find(diff([0; below(:); 0]) ~ 0);这样change_idx里就存了所有状态切换位置奇数位置是游程起点偶数位置是终点一次向量化操作拿到全部游程区间。这个思路是整套代码的地基V1版当时用的是for循环逐点判断数据一长十几万步就明显卡顿V2全面改成这种边界点定位法。1.2 阈值是常数还是动态的这决定了事件定义是否合理很多人做游程分析时只想着找个固定阈值一刀切但灾害事件的阈值往往不该是常数。以干旱为例不同季节、不同地区的降水基准差异巨大一套固定阈值要么在湿润区什么都识别不出来要么在干旱区把所有时段都判定为事件。V2版支持两种阈值模式绝对阈值直接指定一个常数例如径流小于5 m³/s视为枯水。动态阈值按时间窗口滑动计算百分位数比如逐日SPI序列低于-0.8或者径流低于过去30年同期30%分位值。动态阈值的实现核心在于窗口对齐。比如逐月数据要算每月的多年平均或分位数不能错位用movmedian或自定义分组函数都比较方便。代码上我封装了一个getThreshold函数输入数据和模式参数输出与数据等长的阈值数组function th getThreshold(data, time, mode, win) % mode: absolute, moving_percentile, monthly_clim switch mode case absolute th repmat(win, size(data)); case moving_percentile th movpercentile(data, win, 20); % 20%分位 case monthly_clim th monthlyPercentile(data, time, 20); end end动态阈值有一个容易踩的坑滑动窗口的边界段会出现NaN或样本不足的问题。MATLAB的movmedian、movprctile需要统计或金融工具箱在窗口不满时会默认用Endpoints参数填充我一般设为fill并手动补最小值避免边界处误判出大片事件。1.3 事件特征向量除了起止时间你还需要算什么游程识别只是第一步提取灾害事件特征才是用户真正关心的结果。做灾害分析时事件特征向量通常包括这些量特征名称计算公式物理含义事件编号按时间顺序计数唯一标识起始时间游程起点对应的日期事件开始结束时间游程终点对应的日期事件结束持续时间长度L单位为月/天事件持久性峰值强度min(x_t)或max(threshold - x_t)最严重时刻累计强度sum(threshold - x_t)总亏水量/总超阈值量平均强度累计强度/持续时间事件平均严重程度峰值出现时间最小值所在时刻事件演变的关键节点这些特征量组合起来足以支撑后续的趋势分析、概率分布拟合、时空聚类、事件频次变化等研究。V2版在输出特征矩阵的同时也会生成一张事件时间进度条图纵轴每行是一个事件横向是发生时间范围这对快速查看识别的合理性特别有用。2. V2版程序框架设计从一次性脚本到可配置工具箱2.1 V1版的痛点一条产品线里三个脚本改了又改我第一版的时候是接到一个任务——分析某个流域近六十年逐日径流的枯水事件。当时图省事写了一个主脚本阈值固定、事件合并规则硬编码每次换数据都要打开.m文件改参数换一个研究区域就复制一个脚本最后桌面上多出十几个不同命名的文件。更麻烦的是一旦上游数据更新了时间段整个过程重跑一遍输出表格的格式还不统一下游画图脚本经常对不上列名。V1的另一个致命弱点是事件合并太粗暴。我当时简单地规定两个事件间隔小于3天就合并但这个3天是拍脑袋定的而且没有考虑事件本身持续时间长短——间隔3天对于持续20天的长事件可能该合并对于只持续2天的短事件则可能不该合并。这种一刀切的规则直接导致识别结果在长时间尺度上波动很大。V2的核心目标就是规范化这几个问题参数不写死在代码里事件合并规则可配置输出统一结构代码模块化到负责人能直接跑的程度。2.2 模块划分与主控流程我把整个程序拆成六个模块run_analysis.m % 主控脚本配置参数、调用各模块 load_input.m % 数据读取与预处理缺测插值、NaN处理 get_threshold.m % 阈值生成绝对/动态 detect_runs.m % 游程检测返回游程起止索引 postprocess_events.m % 事件合并、最小持续时间过滤、边界处理 extract_features.m % 特征向量计算与输出 plot_events.m % 可视化检视主控脚本长这样参数集中在开头几乎不需要翻后面代码%% 用户配置区 data_file precipitation_daily_1960_2020.csv; time_col 1; value_col 2; threshold_mode moving_percentile; window_years 30; percentile_val 20; min_duration 3; % 最小事件持续天数 merge_gap 5; % 事件合并最大间隔天数 %% 固定流程 [data, time] load_input(data_file, time_col, value_col); th get_threshold(data, time, threshold_mode, window_years, percentile_val); runs detect_runs(data, th); events postprocess_events(runs, merge_gap, min_duration); feat extract_features(events, data, th, time); plot_events(time, data, th, events);这种设计下后续换数据、换阈值方式只需要动第一段配置区其他部分完全不用碰。我自己的体会是写科研分析代码抽象程度不必太高但配置与逻辑分离这一条绝对值得坚持因为你三个月后回来看自己的脚本唯一还能快速上手的入口就是清晰的配置区。2.3 面向过程还是封装类MATLAB场景下的选择有人会问要不要用MATLAB的classdef面向对象重写我试过但最终放弃了。原因是灾害特征提取虽然逻辑不复杂但研究人员通常需要在不同数据、不同阈值、不同合并规则之间快速试出合理结果函数式模块比类封装更轻量也更符合大多数做水文气候分析的MATLAB用户的使用习惯。但如果你的项目要复用性很高比如要打包给团队其他人调用、要集成到更大的分析流程里那建议把核心函数做成独立文件加上help注释做成工具箱的形式。V2版的妥协方案是所有函数以独立.m文件存储互相之间只依赖输入输出结构体不共享全局变量。这样既保留函数式调用的灵活又能调用addpath方式打包共享。以结构体events为例字段定为events struct(... start_idx, {}, ... % 游程起点索引 end_idx, {}, ... % 游程终点索引 start_time,{}, ... end_time, {}, ... duration, {}, ... % 持续时间 peak_idx, {}, ... % 峰值所在索引 peak_value,{}, ... % 峰值强度 intensity, {}); % 累计强度模块之间的数据通过这种结构体传递后续扩展新特征量的时候不需要改函数签名只要在extract_features里加字段即可。3. 核心函数实现游程检测、事件合并与边界处理3.1 detect_runs的向量化实现逻辑detect_runs是整套代码里最核心也最容易被写慢的地方。前面提到用diff找状态切换点完整实现如下function runs detect_runs(data, threshold) % DETECT_RUNS 基于游程理论识别低于阈值/高于阈值的连续区间 % 输入: % data - 列向量时间序列 % threshold - 列向量与data等长的阈值序列 % 输出: % runs - N x 2 矩阵每行是[起始索引, 结束索引] below data threshold; % 状态切换点前后状态不同之处 change diff([0; below(:); 0]); start_idx find(change 1); end_idx find(change -1) - 1; % 如果最后一个切换是状态开始但没有结束说明序列末尾仍在游程中 if length(start_idx) length(end_idx) end_idx(end1) length(data); end runs [start_idx, end_idx]; end你可能会问为什么change用diff([0; below; 0])而不是diff(below)因为这样能把序列首尾的游程边界也抓出来——第一个点如果是belowtrue前面补0后diff会得到1表示游程从这里开始最后一个点如果仍是belowtrue末尾补0后diff会得到-1表示游程在这里结束。否则首尾的事件会被吃掉这个错误在V1版踩过V2直接杜绝。这个函数比逐点if判断快一个数量级以上。十万个时间点循环写法需要十几秒这版向量化写法在零点几秒内完成。如果数据量更大后面可以考虑用find(diff(below))的C加速路径但一般来说这个效率已经够了。3.2 事件合并为什么单纯按间隔阈值不可靠前面说V1版拍脑袋定了间隔3天就合并V2做了改进把合并改成两级判定第一级间隔小于等于merge_gap的事件先被合并。第二级如果合并后的事件持续时间仍小于min_duration则丢弃。这个顺序是刻意的。先合并再滤短可以避免一个长事件被短暂间歇切成两段后每段都过短被误删的问题。实现上合并逻辑就是循环判断两个相邻事件之间的索引间隔。由于游程数量一般几百到几千循环在这里性能损失不大所以用循环反而更直观function events mergeRuns(runs, merge_gap) % MERGERUNS 将间隔过短的相邻游程合并 if isempty(runs) events runs; return; end merged []; cur_start runs(1, 1); cur_end runs(1, 2); for i 2:size(runs, 1) gap runs(i, 1) - cur_end - 1; if gap merge_gap % 合并延伸当前事件结束位置 cur_end runs(i, 2); else merged(end1, :) [cur_start, cur_end]; %#okAGROW cur_start runs(i, 1); cur_end runs(i, 2); end end merged(end1, :) [cur_start, cur_end]; events merged; end我在这里踩过一个真实的坑这个gap的计算公式是runs(i,1) - cur_end - 1还是runs(i,1) - cur_end这要看start_idx和end_idx是否包含边界。我定义的是闭区间索引[start, end]内所有点都是游程的一部分所以中间隔着的天数等于下个起点减当前终点再减1。比如第一个事件在5号结束第二个事件从8号开始实际间隔是6号、7号两天gap 8 - 5 - 1 2。如果写成8 - 5 3就会多算一天可能触发不该发生的合并。这种小细节看代码不觉得跑数据发现事件数量不对劲排查老半天才会意识到。3.3 最小持续时间过滤与缺测段的特殊处理min_duration过滤相当简单但缺测段NaN的处理就有讲究。游程检测时数据里不能有NaN否则data threshold得到false会把一个事件从中间劈开V1版在这个问题上出过很诡异的Bug。V2版采用的做法是数据加载后先做缺测插值。日值降水如果用线性插值可能会把有雨和无雨之间造出一些虚假的中间值这在干旱事件识别里不能接受。更稳妥的做法是对缺测天数极少小于总长度的1%的序列用临近值填充对缺测段过长的数据序列应直接设置thresholdNaN并显式把这些区段排除在事件识别之外同时在结果里标记该时段数据缺失事件边界不可靠。实际实现中我在load_input里增加了一个返回字段quality_flag缺失部分标记为0在detect_runs里遇到缺测段就强制将状态置为非事件并在postprocess_events里把横跨缺测段的事件拆开。function below markBelow(data, threshold, missed) below data threshold; below(missed) false; % 缺测段强制为非事件 end这样处理之后即使数据有间断也不会污染事件识别结果。4. 关键参数对事件识别结果的影响最小持续时间与合并窗口4.1 参数敏感性同一个序列不同参数得到不同事件数我用自己的流域径流数据1960-2020逐日做了个快速参数扫描这里给一组典型结果最小持续时间(天)合并间隔(天)识别事件数平均持续时间(天)0018711.20514314.63010318.2357824.5755233.17104438.7这组数字非常直观参数变一下事件数量几乎翻倍变化。不存在标准答案参数必须根据研究目的来确定。如果关注的是极端长期干旱min_duration就该设大一些如果关注的是所有短暂缺水过程对农业的影响min_duration应设小或为0。但不管怎么设有一点必须谨记在报告或论文里必须明确写出参数值而且最好做敏感性分析说明事件特征对参数的依赖程度。否则审稿人或合作者很容易质疑结果的可靠性。4.2 周期性与游程合并的关系季节性因子不可忽视游程合并还有个季节性陷阱。如果研究区域是典型的夏季风气候干季和湿季分明那么干季的负游程天然很长不需要合并但如果在春秋过渡期降水序列会频繁跨越阈值导致一场秋季缺水过程被切成好几段。这种情况用单一的merge_gap参数很难兼顾。V2版在配置里增加了一个可选参数seasonal_merge允许按月份设置不同的合并间隔比如夏季设7天、冬季设3天。实现上就是把前面mergeRuns函数里的固定merge_gap换成一个与事件位置相关的数组gaps getMonthlyMergeGap(merge_gap_seasonal, time, runs);这个扩展我用过之后觉得很值——它避免了一个窘境参数设大了把独立事件错误拼在一起设小了又把真实的长事件拦腰截断。如果你处理的是有明显季节性的数据强烈建议加这个机制。4.3 阈值分位数选择20%和10%差异巨大阈值分位数是另一个敏感参数。同样是动态阈值用20%分位值和10%分位值识别出的干旱事件强度特征完全不同。10%分位数识别出的事件频率低、但平均强度大更适合做极端干旱研究20%分位数则更适合做整体缺水状况评估。这也带出一个方法论的提醒阈值本质上是事件定义的一部分而不是单纯的算法参数。在动手写代码之前先想清楚你的灾害类型对应什么样的重现期或分位数标准。比如很多国内干旱研究采用径流距平百分率低于-20%或SPI值低于-1.0作为阈值这些是领域内约定俗成的标准直接选用能让结果更好被同行比较。5. 实测案例全流程从逐日径流序列到干旱事件特征表5.1 案例数据与预处理为展示完整流程我用一段模拟生成的逐日径流数据2000天约5.5年来走一遍。数据包含明显的季节性周期和几段人为设置的枯水期。读者可以把这部分替换成自己的实测数据。% 模拟数据 rng(42); n 2000; time datetime(2000, 1, 1) days(0:n-1); seasonal 30 20 * sin(2 * pi * day(time, dayofyear) / 365.25); noise 8 * randn(n, 1); flow seasonal noise; flow(flow 0) 0; % 人为制造两段持久的枯水期 flow(400:460) flow(400:460) * 0.2; flow(1200:1300) flow(1200:1300) * 0.1;预处理阶段我先用fillmissing线性插值补齐可能的缺失值再设定动态阈值为过去一年滑动窗口的20%分位数。窗口设1年既保证样本量足够又能跟上季节变化。5.2 游程识别与合并的可视化检视跑完detect_runs和postprocess_events后第一件事不是直接看表格而是先可视化因为人眼是判断算法合理性的第一道防线。我画了这样一张图径流曲线灰色显示动态阈值用红色虚线识别出的枯水事件用浅蓝色阴影标出。结果中发现开头第100天附近识别出一个持续6天的事件但我肉眼看着那一段其实只是正常的季节低值并没有真正枯到异常。为什么被识别因为动态阈值恰好在高值季节那条低值段的径流仍然低于滑动20%分位。这说明动态阈值在高变异季节容易产生误报。解决方案之一是抬高阈值分位数或者改用超过特定距平的人为标准。我在代码里预留了threshold_bias参数允许用户对动态阈值增加一个绝对偏移量这种情况直接手动调2就能压制误报。5.3 特征表输出与下游分析最终特征表长这样节选事件编号开始日期结束日期持续时间(天)峰值亏水量(m³/s)累计亏水量(m³·s⁻¹·天)平均强度12000-06-122000-06-26158.2351.43.4322001-02-022001-03-053212.67168.95.2832002-12-102003-01-184018.02296.37.4142003-02-132003-02-2085.9621.82.73输出格式上我用了writetable直接写出CSV列名与特征结构体字段保持一致。下游做频率分析时直接读表即可不用再改代码。为了避免列名歧义CSV里附带一行注释说明各列的物理单位这一招在处理合作项目时帮我避免了很多沟通成本。6. 调试过程避坑记录与性能优化经验6.1 首尾游程丢失问题diff边界补齐的教训V1版我用diff(below)查找切换点时序列开头就是负游程的情况下diff结果不会出现1导致第一段事件永远不会被识别。当时给我的分析结果带来很大的偏误——最干旱的区域反而被识别出最少的事件。V2版在状态序列前后各补一个0再diff正是为了解决这个边界问题。这里再强调一次change diff([0; below(:); 0]);这行代码是整个程序的隐形功臣。你如果发现自己跑出来的事件数量异常偏少优先检查这行事件从第2个点开始的情况往往就是这里出了问题。6.2 峰值搜索在等值平台段的处理clip还是中心点计算峰值强度时如果序列中出现连续多个点数值完全相同这种情况在低分辨率数据里很常见find会返回第一个最小值这没问题。但如果你用min(data(start:end))找到的是数值还想要峰值出现的时刻就要小心等值平台段里第一个最小点可能明显偏早或偏晚从而影响你对峰值出现时间的估计。我采用的策略是取平台段的中间索引作为峰值出现时间% 找到最小值的全部位置取中间点 [min_val, ~] min(segment); min_positions find(segment min_val); if length(min_positions) 1 peak_idx round(mean(min_positions)); else peak_idx min_positions(1); end6.3 性能对比向量化 vs 循环以及memory的取舍我用一个60万步的日值序列约1643年做了测试detect_runs循环版耗时约85秒向量化版本0.3秒。差距巨大所以向量化是必选。但如果序列极度长例如5分钟间隔水文数据一年就超过10万步构建below逻辑数组和diff结果本身也会消耗内存。可以用uint8类型省空间MATLAB逻辑数组已经只占1字节所以below本身开销不大diff结果的double类型相对占内存但序列长度百万量级也就几MB可以忽略。真正的大内存大户是动态阈值数组的构建过程movpercentile在长序列上可能要反复排序内存峰值会高。这种场景下把序列分段处理段间保留一个merge_gap的overlap是更稳妥的路线。6.4 批量处理多个站点的数据组织技巧做区域分析时经常要遍历几十上百个站点的数据。不要在主脚本里叠一个巨大for循环把站点的参数数据文件路径、阈值模式、合并间隔、最小持续时间放在一个表格里用parfor并行处理每个站点独立输出特征表和诊断图。我实际跑过一个小流域42个站点parfor在四核机器上从42分钟降到14分钟节省的时间足够用来多调几轮参数。需要注意parfor里头的函数必须全部是显式输入、显式输出不能依赖全局变量和脚本环境变量这也是我把所有逻辑都独立成函数文件的原因之一。6.5 关于MATLAB版本与工具箱依赖这套代码尽量只用基础MATLAB函数核心用到diff、find、min、datetime、movmedian此函数在R2016a之后才可用。如果手头版本太老movmedian可以用circular滑动窗口自写替代。画图部分用到patch和datetime轴这些也都在基础版里。有些人问要不要用工具箱里的findpeaks或信号处理函数。我的建议是不要过度依赖工具箱尤其当你需要把代码分享给没有相应工具箱的合作者时一个工具箱函数可能就是阻碍复现的门槛。游程理论本身足够简单自己实现更有掌控感也能随时按需修改。7. 从V2到后续还有哪些扩展切入点一套游程识别程序跑通以后你会发现它适合做的事情远不止干旱识别。给径流序列加一个正阈值它就变成洪水事件提取器把输入从径流换成气温低于阈值的游程就是寒潮过程换成空气质量指数高于阈值的就是污染过程。本质都是阈值切割连续区间识别特征量计算这一套逻辑。在这个框架基础上我建议你继续往这几个方向扩展一是把识别结果直接接入概率分布拟合如GEV、Gamma分布去算重现期很多灾害风险评估项目都需要这步二是增加空间模块多个站点的事件特征可以汇合到GIS环境中做时空制图三是考虑事件强度与时间的关系加入趋势检验Mann-Kendall函数直接输出特征序列的趋势显著性结果。就个人经验而言代码库的维护最关键的不是写得多炫而是命名一致、结构稳定、配置清晰。V2版做完后我后续给同一流域做极端降水分析时只写了两个新模块一个处理正游程一个把事件过程从原始序列中截取出来做事件合成分析。核心的detect_runs和postprocess_events一行没改。这让我对通用工具四个字有了新的理解——不是写出无所不能的框架而是把最基本的环节做扎实让上层变化有稳固的落点。