做自然灾害类的数学建模，最容易掉进的一个误区，就是把每一问当成彼此独立的机器学习或偏微分方程习题：第一问跑个相关性排序，第二问算个安全系数，第三问拟合一条时间曲线，第四问再单独跑个流体仿真。但真实的泥石流灾害根本不是这样发生的——灾前气象背景如何改变孔隙水压？源区岩冰体在什么条件下失稳？崩塌物能不能转化为泥石流？进入沟道后又要多久才会冲到下游口岸？
在 2026 年第七届 MathorCup 大数据竞赛初赛（赛道 B“泥石流灾害预报及防灾策略”）中，我和队友围绕吉隆“岩冰崩塌—碎屑流—泥石流—山洪”灾害链做了一次完整的全流程探索。比起追求表面上好看的准确率，我们花最多精力的地方，其实是怎么把气象、遥感、地形与水文多源数据塞进同一套可回查的物理框架里，并且不掩盖任何一个失败算例或数据盲区。
一、先把四道赛题放进同一条物理链
拿到赛题后，我们做的第一件事不是急着写代码，而是先把四个子问题锚定到同一套空间坐标与时间戳上：固定共同源区、沟道入口断面 C 和口岸附近的参考断面 D，把整套模型拆解为四个首尾咬合的模块：
| 研究阶段 | 核心问题 | 我们采用的方法 |
|---|---|---|
| 1. 气象关联 | 灾前区域气象背景有哪些显著偏离？ | 112 个同季节共同有效背景窗对照 |
| 2. 失稳与形成 | 源区在什么力学与水源条件下失稳并成流？ | Mohr–Coulomb 条件力平衡与参数单调上下界 |
| 3. 触发时间 | 如何表达首次越界时刻并评估区间覆盖？ | 遥感前兆外推 + GFS 驱动条件预测（含右删失） |
| 4. 沟道传播 | 进入入口 C 后多久到达断面 D、规模多大？ | 一维有限流深守恒仿真与双网格收敛核验 |
这样做的好处是物理阶段非常清晰：源区失稳 → 泥石流启动 → 进入断面 C → 向下游断面 D 传播。 我们在代码里严格区分了“源区首次失稳时刻”、“C 入口起算时刻”和“D 断面强流检出时刻”，连地图桥点与模型断面之间的空间偏移也单独做了核对。只有把这些接口定义死，后面谈“传播历时”和“预警提前量”才不会变成一笔糊涂账。
二、动手建模前，先和脏数据对账
多源数据建模最怕“垃圾进、垃圾出”。我们在前期花了大量时间做原始资料对账：
- 气象表格逐格盘点：针对 24 份小时级 CSV（共 3,312,792 个测量单元），我们先排查缺测码、文本空值和质量标记，在保留原件的前提下生成严格质量视图。计算窗口均值时要求至少 80% 数据完整，而计算降水累积时则要求 100% 全窗完整——因为一旦把缺测当成 0 毫米混进累加器，就会直接伪造出“干旱无雨”的假象。
- 遥感影像的晴空交集：重建源区边界时，我们把灾前 14 景与灾后 6 景光学影像做了 84 种配对交叉检查，发现最佳共同晴空覆盖率其实只有 36.94%。这意味着任何边界判读都必须老老实实标出“未共同观测区”，而不能靠插值去硬画一条看似完整的闭合曲线。
- 物理意义的边界：多年 DEM 高程差可以作为体积代理，但里面混杂了期间演化与配准误差；InSAR 雷达相干性反映的是表面散射特征的变化，不能直接等同于滑坡位移；而河流水位如果缺乏实测断面评级曲线，也不能随意换算成洪峰流量。
为了保证论文里写的每一个数字都能查到源头，我们用 Git 把配置文件、清洗脚本、中间小表和报告生成管线全部锁在一起，做到任何一张表都可以从原始数据一键重算。
三、气象关联：为什么“没下雨”不能简单当成异常
第一问要找出灾前最异常的气象因素。我们在检查数据时发现，不同气象变量的缺测时段并不一致——如果每个变量各自挑能用的年份去算均值，最后排出来的名次很可能只是样本池不同造成的假象。因此，我们先把 11 项候选指标强制对齐到 112 个共同有效的历史同季节背景窗上，再比较事件值相对历史背景的偏离程度。
排出来的前五项依次是：气温、净辐射、相对湿度、垂直温度梯度和潜热通量。其中气温和净辐射非常稳健，而排在第五名的潜热通量与风切变同分，在年份重采样（Bootstrap）中的入选比例只有约 44.14%。与其在报告里只报一个光鲜的 Top 5，我们选择把这份对样本年份的敏感性直接画出来：
这里还有一个很有意思的细节：零降水的处理。 灾前 24 小时当地确实没有降雨，但如果直接按分位数硬算，很容易得出“极度干旱异常”的结论；可实际上，在 112 个历史同期窗口里，原本就有 64 个窗口滴雨未下。我们专门修正了并列零值（Ties）的统计处理，才避免了把高山地区司空见惯的无雨日误判成极端气象信号——这也进一步印证了吉隆这次事件更可能是受高温与辐射驱动的岩冰热力失稳，而非传统暴雨型泥石流。
四、源区失稳：用条件力平衡守住未知的边界
到了第二问，我们用 Mohr–Coulomb 极限平衡模型把坡角、内摩擦角、黏聚力和孔隙水压力装进安全系数 FS（当 FS ≤ 1 时达到临界失稳）。由于缺乏源区钻孔实测的岩土参数，直接拍脑袋定死一组参数是没有说服力的；所以我们先独立复算了 440 行临界面应力平衡，再利用安全系数对参数单调的数学性质，给出了 286 个条件参数组合下的 FS 上下界区间。
更重要的是：山体滑下来（失稳），不代表一定能变成泥石流并冲进沟道入口 C。 从失稳到形成泥石流，还必须同时满足松散物源充足、短时融水或水源汇入、以及地形连通三个条件。我们在代码里写了一条硬规则：只要任一形成条件未知或不满足，或者阶段转换时延缺失，后续状态就严格保持为空值（未知），绝不把“数据缺失”偷偷填成“零风险”。
五、触发时间：给“窗口内没滑下来”留出位置
第三问预测触发时间时，我们走了两条平行路线：一条用真实的遥感形变前兆序列对比常速、线性趋势和逆速度（Inverse Velocity）外推；另一条则固定起报批次的 GFS 数值天气预报温度输入，串接孔压响应与安全系数模型，推演未来 48 小时内何时首次跌破 FS = 1。
在 3,544 个初始稳定的构造样本中，有 1,638 个在 48 小时窗口内触发了失稳，而另外 1,906 个直到窗口结束依然稳定。做工程预测时，这 1,906 个样本千万不能直接丢掉——在生存分析里它们属于典型的右删失（Right-censored）数据，即“我们只知道它在 48 小时内没滑，但不知道它具体哪一刻会滑”。因此，我们在报告里同时给出了包含右删失的整体预测范围，以及仅针对窗内触发样本的条件预测区间。
为了检验名义 95% 预测区间的含金量，我们还特意做了“压力测试”：在 12 种假设偏离情景下测试覆盖率退化，甚至拿外地源区的两组保留物理算例去跑盲测——结果首次有限区间发布全部漏掉了真实目标时刻。我们把这个失败结果原封不动地写进了最终报告，因为诚实交代模型在跨区域外推时的失效边界，比编造一个“准确率 99%”的神话更有工程价值。
六、沟道传播：首达时间只差 10 秒，体积却差了一半？
第四问给定断面 C 的入口水沙过程后，我们搭建了一维有限流深守恒模型来算泥石流何时冲到下游口岸断面 D。在这个环节，我们踩到了整场比赛最值得庆幸的一个坑。
起初我们只盯着“强流到达 D 断面的时间”。在对比粗细两套空间网格时，有一组方案的首达时间仅仅相差 10 秒，看起来似乎已经完美收敛了；但当我们把入口、出口、断面通量和域内存量加起来做质量守恒对账时，却震惊地发现：600 秒内的净通过体积在两套网格之间竟然相差了 50.76%！
顺着这个异常往下查，我们定位到了数值格式上的病根：旧的一阶迎风分支在薄层陡坡流动下会严重低估重力加速度。于是我们重构了二阶线性空间重构格式、完整保留单元内部的床坡源项，并用静水保持、恒坡启动和溃坝解析解重新做了基准检验。改进之后，四组固定设计在 0.9375 米和 0.46875 米两套细网格上跑完 16 次仿真，核心指标的最大网格差异从 50% 降到了 0.52%：
即便数值收敛到了 0.52%，我们依然在文档里标明了它的物理局限：加密计算网格并不会凭空增加真实地形的采样精度，而等宽矩形沟道、固定床面和忽略沿程沟床卷吸（Entrainment）的假设，意味着这套一维模型适合用来做条件传播时间与通量量级评估，但不能直接拿来画下游村庄的二维淹没红线。
七、从仿真数字到真实预警的距离
很多建模论文算完“泥石流 15 分钟到达下游”就直接宣布预警成功了，但在实际防灾链路里，物理传播时间 ≠ 可用疏散时间。
我们在最后把预警拆成了“源区灾前预报”和“滑震/入沟后的在途告警”两本账：就算传感器在泥石流进入 C 断面的一瞬间捕捉到了波形，后面还要扣除数据回传、算法确认、人工研判、告警下发和居民响应的耗时。把这些工程延迟全部扣掉之后，才能看出下游究竟还剩多少真正的安全余量。
八、这次比赛留给我的工程复盘
为了防止最后写论文时出现“代码改了、图更新了，但正文里的数字还是旧版”的低级错误，我们把整个工程拆成了流程编排、科学计算、自动验证和报告生成四层：所有正文引用的关键数值都直接从通过校验的结果小表里读取，彻底告别手工复制粘贴。
回看这几天的密集建模，我最大的体会是：在复杂系统工程里，定义模块之间的接口与边界，和推导每个模块内部的公式同样重要。 只有把数据缺口、右删失样本、守恒误差和跨区域失效算例都坦诚地摆在台面上，我们才真正知道模型在哪一步是可靠的、在哪一步需要向现场实测数据低头。
资料来源
本文整理自我们团队的参赛论文《多源资料约束下的吉隆灾害链分析与条件预警》、全流程复盘笔记及配套验证代码库。文中两幅配图均来自论文原图，所有统计指标均保留其对应的有效样本分母与预设工况条件。