企业官网建设流程全解析

做水旱灾害风险分析的朋友，大概都有这种体会：数据本身不会缺，逐日降水、逐日径流、水位过程线、标准化降水指数SPI，一拉就是几十年甚至上百年，真正让人头大的是怎么把“灾害事件”一个不落地挑出来，还要把发生的起止时间、持续多久、累积量多少、峰值多极端这些特征准确地算清楚。游程理论（Run Theory）就是专门干这件事的经典方法，而这套基于MATLAB语言实现的V2版本，是我在实际项目中反复打磨后沉淀下来的一套完整方案。它会告诉你如何从原始时间序列出发，自动识别出每一个灾害事件，输出一张可直接用于后续统计和绘图的事件特征表，并且把“识别—合并—计算—可视化—导出”整条链路整合成一个可复用工具箱，而不是一次性脚本。

1. 项目背景与V2版本设计思路

1.1 游程理论的核心思想与适用场景

游程理论在工程水文和气象灾害分析中用了很多年，原理其实一句话就能说清楚：给定一个阈值，把时间序列的每一个样本都标记成“超过阈值”或“低于阈值”，连续处于同一种状态的样本段就是一个“游程”（Run）。游程就是候选事件，接下来要做的只是对每个游程去计算持续时间、累积量、峰值等一系列特征。

这个方法真正厉害的地方在于，它把肉眼判断变成了一个完全可复现的算法流程。以前我看历史降水资料时，习惯在Excel里拿眼睛一段一段扫，数据短还行，一旦面对100年逐日序列，人工标注根本不可行，而且不同人标出来的结果还不一样。用游程理论做事件提取后，只要阈值确定，事件清单就是唯一确定的，后面的重现期分析、Mann-Kendall趋势检验、空间风险区划，所有统计结果都能在这个清单基础上稳定复现。

V2这套实现里，我默认支持两类最常见的目标场景：

从表里也能看出来，游程理论本质上是个通用的事件分割器，不挑领域，关键是阈值选得合理、后续特征统计符合业务需要。这套代码在设计时就没有绑死“干旱”或“洪水”，而是把“低值事件”和“高值事件”都做进了同一个函数，切换场景只用改一个参数。

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(end+1, 1) = startIdx(k); newEnd(end+1, 1) = endIdx(k); end end end

这个合并函数写起来简单，但设计时要想清楚一个细节：合并只会影响起止索引，不能在中途把两个游程的简单长度相加，因为中间那段超过阈值的间隔其实也属于事件过程的“复发间隔”。正确的做法是合并完成后，再用新的起止索引去原始序列里重新截取数据，重新计算累积量和峰值。这个逻辑我在核心函数里做了严格分离，避免出现“只合并区间但不重算特征”的隐性bug。

mergeGap取多少完全看业务需求。识别干旱事件时，逐日序列建议先试3天到5天；识别暴雨过程时，因为降水过程连续性更强，一般可以放宽到7天左右。后面第5章我会专门讲参数敏感性怎么检查。

2.3 特征指标怎么选：历时、累积量、峰值和平均强度

事件识别出来之后，特征指标体系就决定了你的分析能上升到什么层次。V2版本每个事件最终都输出六个字段，前四个是核心特征，后两个是派生特征。核心特征直接决定事件的“几何形状”，派生特征则用于事件严重程度的横向比较。

为什么要单独区分低值和高值的计算方法？因为干旱的“强度”体现在缺水量，也就是阈值减去实际值的累积；而洪水的“强度”体现在超限量，也就是实际值减去阈值的累积。如果统一用原始值累加，两个方向的物理意义都会乱套。

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(end+1, 1) = startIdx(k); newEnd(end+1, 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或者fill？area在x轴上有基线，不适合直接标注一段背景；fill会被多条曲线干扰图层顺序；patch配合FaceAlpha是控制背景标注最灵活的方式。

图片导出方面，我现在基本弃用print，全部改用exportgraphics，因为它在分辨率控制和字体渲染上更稳。要出论文图就300dpi起步，要矢量图就导出eps：

exportgraphics(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));

运行后的事件特征表示例长这样（具体数值随随机种子不同会有差异）：

从表格可以直观看到，mergeGap=2把很多原本零散的小游程合并成了持续几个月的事件。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这套代码之所以比第一版顺手，不是因为原理变了，而是因为每个开关都做了明确的输入输出约束，不再需要每次复制脚本后去改逻辑关系。如果你也要在项目中应用游程理论，我的建议是先用模拟数据把函数跑通，再用真实数据做参数敏感性检查，最后才进入正式统计环节。编码本身只是把成熟的思路落地，真正花心力的是理解和验证你设定的每一个参数。