一个「看起来像修 bug」的补丁，为什么最后被撤掉？

一、事情的起点

AERMAP 在给受体/源取地形高程时，先判断这个点落在哪片 DEM 瓦片内。判据用的是 DMCNR——从 GeoTIFF 的 tiepoint 推算出来的瓦片四角坐标，写在 sub_nedchk.f 里：

DMCNR(idem,2,1) = DMCNR(idem,2,2) - (DYM(idem) * (dble(nRows(idem)-1)))  ! 西南角纬度
DMCNR(idem,1,3) = DMCNR(idem,1,2) + (DXM(idem) * (dble(nCols(idem)-1)))  ! 东北角经度

现场出现过这样的现象：源点紧贴瓦片边界时，AERMAP 报：

OU E331 Gap Receptor NOT Inside DEM File
OU W330 DEMREC: Receptor NOT Inside DEM File - First Pass

整个 AERMAP 运行作废；AERSCREEN 更麻烦——它会静默回退为不考虑地形继续算，报告里看不出异常，但评价等级可能被低估。于是有了一个很自然的推断：

当时的推断

这两处减 1 是「差一」错误（off-by-one）：瓦片有 nRows 行，就该按 nRows 算，减 1 等于少算最外一行/一列，贴着边界的点因此被判在外。

顺着这个推断，补丁把两处改成 dble(nRows(idem)) / dble(nCols(idem))。改完确实“不再报 NOT Inside”（在一些算例上），于是它被当成边界修复合并进了随包二进制。

二、验收方式的问题

回过头看，当时的验收标准是两条：

1. 不再出现 NOT Inside DEM File；
2. 换成新二进制后，既有算例“能跑通、结果差不多”。

这两条都不足以证明“只修了边界”。边界判断和坐标取样共用同一组变量，动了这组变量，最先变化的往往不是边界，而是所有点的取样位置——而这在“能跑通、结果差不多”的验收里是看不见的。

三、用真值重新验收

2026-11-12，为把三端（演流 / AERSCREEN-ZH / Web）的地形口径统一，我们做了一次更严格的对照。环境统一为：

• DEM：Copernicus GLO-30，1°×1° 瓦片（N23_00_E113_00），30 m，NAD83 GeoTIFF
• 算例：项目中心 23.217°N / 113.483°E，受体 21×21、250 m 间距、5 km 范围（441 个）
• 气象：AERMET 生成的 2024 年 .SFC/.PFL；源：1 g/s SO₂、烟囱 50 m

3.1 高程：441 个点全部变化

只看一处数据可能还像“个别点不同”，做全量比对就很清楚了：441 个受体，870 行输出，没有一行相同。这不是边界修复，这是整幅 DEM 的取样网格平移。

3.2 与 GeoTIFF 真值比：改后更远

用 GeoTIFF 自身的 tiepoint 与 pixel scale 做双线性插值，作为独立真值，对 36 个内部点比较：

在丘陵地形上，一个像元（30 m）的水平位移会带来几米到几十米的高程差。NED 版不但没有更准，反而更偏离真值——因为它把取样原点整体挪了一格。

3.3 下游浓度：差 25.88%

把两套受体高程分别喂给同一颗 AERMOD（其余完全一致），比较 1 小时最大地面浓度：

• 266.24935：官方取样口径 · µg/m³（峰值位于 -500, 1250）
• 335.16366：NED 修正版 · µg/m³（峰值位于 -500, 1000）
• 同一工程、同一气象、同一源强下的差异 +25.88%

这不是舍入误差，也不在“差不多”的范围内。一个补丁让某个项目的地形影响被系统性放大四分之一，这才是必须撤掉它的理由。

四、更要命的是：它并没有修好边界

把点正好放在瓦片下边界（23.0000°N），只给一片瓦片，四种构建的表现是：

为什么？把 nRows-1 改成 nRows 后，西南角理论上正好落在 23.0，但浮点计算（24.0 − 3600 × (1/3600)）得到的值可能略高于或略低于 23.0，而判据是 .GE. / .LE. 的精确比较——压线的点在数学上就没有“内”可言。换句话说：

这个补丁付出的是确定的代价（全图取样平移一格），换来的却是不确定的收益（压线点仍然按浮点运气决定）。

五、真正解决边界的是什么

问题不在“边界怎么算”，而在“只给了一片瓦片”。相邻瓦片本身是接壤的，只要把覆盖范围向外扩一点、把相邻瓦片一并交给 AERMAP，落在缝隙里的点自然有另一片接住：

包围盒 → 各方向外扩 1 km → 取整度瓦片 → 与已有瓦片求差 → 缺哪片补哪片

实测同一颗内核、同一个 23.0000°N 的点：

于是三端都改成“覆盖检查 + 缺片补齐”：

• 演流（Tk 版）：_run_aermap() 拿上游已有的“项目中心经纬度 + 最远受体/源半径”构造包围盒，缺瓦片就自动补下并合并。
• AERSCREEN-ZH：AERMAP 由内核内部调用，改为下载半径 +1 km；并在写 DEMlist.txt 前用源点经纬度检查覆盖，缺片直接报错说明，不再静默回退平地（老工程没存经纬度时才退回用 UTM 反算）。
• 演流（Web 版）：本来就按“评价范围 + 1 km”取瓦片。

共用逻辑在 dem_cover.py（纯标准库：经纬度 ⇄ UTM 正反算、瓦片名解析、必需瓦片计算；正算用于把界面输入转 UTM，反算仅在老工程缺经纬度时兜底）。

六、剪枝为什么可以留

同一轮对照里，CALCHC 纬度带剪枝的表现正好相反：

• 60 个受体：剪枝版 9 s，官方版 45 s；
• 受体高程逐字节相同；AERMAP.OUT 去掉时间戳后无差异。

原因是它在数学上可以证明等价：节点入选条件是 ZNODE − AZELEV ≥ 0.1×RDIST，而 RDIST ≥ |Δ南北向|，所以当 |Δ南北向| > 10×(MAXEL − AZELEV) 时该节点必然不入选；剪枝只是把“必然不入选”的节点整段跳过，浮点路径与节点顺序都不变。

判据很简单

改内核前先问一句：这个改动会不会改变任何一次浮点运算的输入？会 → 必须逐字节回归；不会 → 才谈得上“纯性能优化”。

改内核前先问一句：这个改动会不会改变任何一次浮点运算的输入？会 → 必须逐字节回归；不会 → 才谈得上“纯性能优化”。

七、最终处置

因为取样口径回到官方，同一工程的进一步预测结果会变（本算例 −20.6%）。

• 演流 v1.4.158 起、AERSCREEN-ZH v1.2.48 起采用官方取样口径；
• 需要复现旧结果，请继续使用 v1.4.157 / v1.2.47。

八、可以带走的四条

1. “边界 bug”要先确认坐标系约定。 DMCNR 是像元中心口径，nRows−1 本身是对的；把角点口径套到中心口径上，就变成了整体平移。
2. 验收要用真值，不要用“能跑通”。 这次能翻案，靠的是 GeoTIFF 双线性插值 + 441 点全量逐字节比对。
3. 静默回退比报错更危险。 AERSCREEN 曾在地形丢失后继续算；现在缺片会明确报错并给处理建议。
4. 性能优化和数值修正要分开管理。 前者可以证明等价、可以长期保留；后者必须逐案回归，不能靠“看起来一样”。