把 46 分钟压到 9 分钟：AERMAP 临界山高计算（CALCHC）的纬度带剪枝

目录

1. 摘要
2. 问题现象：不是卡死，是真的要算 47 分钟
3. AERMAP 在算什么：临界山高与 10% 坡度判据
4. 剪枝的数学依据
5. 实现：只跳过“不可能分支”
6. 另一处改动：NED 边界的 off-by-one
7. 实测与逐位回归
8. 被排除的假设
9. 合规性：为什么可以改 AERMAP，但不能碰扩散内核
10. 复现步骤
11. 局限与后续
12. 附录

摘要

AERSCREEN 在“点源 + 复杂地形”工况下，一次筛选要跑 46.7 分钟；用户侧表现为界面长时间无输出， 旧版 900 s 超时误杀父进程，还留下孤儿进程，看起来像“卡死”。

真正的耗时在 AERMAP 的临界山高例程 CALCHC：对每个受体、每个 DEM 剖面、整列约 1900 个节点， 逐个做坐标更新、取高程、并调用大地线反算 GPNHRI。1 万个受体约 17 亿次节点级运算， 占 AERMAP 全部运行时间约 75%。

改进没有动判据、没有动数据、没有动任何默认参数，只利用判据本身的单调性：

节点入选条件是 ZNODE − AZELEV ≥ 0.1 × RDIST，而 RDIST ≥ |南北向距离|。 因此当 |Δ南北向| > 10 × (剖面最大高程 − 受体高程) 时，该节点必然不入选，可整段跳过。

一致性不是“基本一致”，而是逐字节相同：AERMAP.KIL、.OUT、 MAPDETAIL、aermap_srcloc.out、aermap_rcpt.out 全部逐字节一致； 下游 AERMOD 输入文本与 .PLT 浓度逐字节一致。

一、问题现象：不是卡死，是真的要算 47 分钟

AERSCREEN 的典型调用链是：

AERSCREEN → AERMAP（受体高程 / 临界山高）→ AERMOD → 结果

一次筛选会调用 AERMAP 三次（源高程一次、受体高程一次、山高/详细地形一次）， 每次都独立跑 DEM 检索。当项目位于丘陵或山地、又使用 1 弧秒 DEM 时：

• 单个 AERMAP 调用可达十几分钟；
• 三次叠加后，单次筛选 40–80 分钟；
• AERSCREEN 自身把 AERMAP/AERMOD 作为子进程调用，控制台只有极少量输出，从外部看就是“无响应”。

旧版前端设了 900 s 总超时，超时后只杀掉 AERSCREEN 父进程，AERMAP 子进程仍在跑， 于是出现“进程已停止、CPU 仍在满负荷”的孤儿进程现象。

定位结论：不是死锁，也不是 DEM 下载失败，而是 AERMAP 的算法常数太大。

二、AERMAP 在算什么：临界山高与 10% 坡度判据

2.1 临界山高的用途

复杂地形下，AERMOD 需要一个“地形影响尺度”来约束初始扩散。AERMAP 为每个受体计算 临界山高（hill height scale） AZHILL(IREC)，其物理含义是：

受体周围、相对受体仰角不低于 10%（坡度 0.1）的地形点中，最高的那个高程。

也就是说，只有“从受体看过去足够陡、足够高”的地形才会抬高山高尺度；平缓的远山不影响结果。

2.2 代码级判据

在 sub_calchc.f 的节点循环里（EPA 原始逻辑）：

C*              Determine the Lateral Distance (R) Between
C*              Receptor and that of this Point
                IF (IPLAN(IDEM) .EQ. 0) THEN
                   PLAT(2) = (YNODE / 3600.0D0) * dtorad
                   PLON(2) = (XNODE / 3600.0D0) * dtorad
                   RDIST = DISTLL(PLAT,PLON)          ! 大地线反算
                ELSE IF (IPLAN(IDEM) .EQ. 1) THEN
                   XARG = XNODE
                   YARG = YNODE
                   RDIST = DSQRT((XRSHFT-XARG)**2
     &                         + (YRSHFT-YARG)**2)
                END IF

                BP = RDIST * 0.1D0

                IF (ZNODE - AZELEV(IREC) .GE. BP) THEN
                   HEFF = ZNODE
                   IF (HEFF .GT. HEFMAX) THEN
                      HEFMAX = HEFF
                      AZHILL(IREC) = HEFF
                   END IF
                END IF

两个关键事实：

1. 判据是 10% 坡度：ZNODE − AZELEV ≥ 0.1 × RDIST，等价于仰角 ≥ 5.71°。
2. 距离只在判据里出现一次，而且是以 0.1 × RDIST 的下界形式出现。

2.3 计算量

以 1 弧秒 DEM、约 25 km 探测半径为例：

• 每个剖面（profile）是一列南北向节点，约 1900 个；
• 单受体涉及的剖面数 × 节点数 ≈ 16.5 万节点；
• 每个节点至少一次 GETNOD（坐标累加）、一次 GETDEM（直接访问文件随机读）、 一次 DISTLL → GPNHRI（Helmert–Rainsford 大地线反算，含多次三角函数与迭代）。

1 万个受体 → 约 17 亿次节点级运算，其中大地线反算是绝对大头。

单受体微基准也印证了这一点：打桩进入 CALCHC 约 3.86 s；把节点循环截断到 120 个后为 4.82 s； 不截断时 15–18 s。耗时与节点数近似线性，而与“取到多少高程”关系不大—— 说明瓶颈在“整列遍历 + 每点算距离”，不在 I/O。

三、剪枝的数学依据

3.1 定理

记：

• Z_r = AZELEV(IREC)：受体高程；
• M_p = MAXEL(IPROF)：该剖面所有节点高程的最大值（已在索引文件中预先算好）；
• y_k：节点 k 的南北向坐标（经纬度 DEM 为纬度弧秒，UTM DEM 为北向坐标）；
• y_r：受体的南北向坐标（RLATSHFT / YRSHFT）；
• R_k：节点到受体的大地线（或平面）距离。

若

|y_k − y_r|  >  10 × (M_p − Z_r)

则该节点不可能满足 10% 坡度判据。

3.2 证明

节点入选需要

Z_k − Z_r  ≥  0.1 × R_k

而距离必不小于南北向分离：

R_k  ≥  |y_k − y_r|

又因为 Z_k ≤ M_p，所以

0.1 × |y_k − y_r|  ≤  0.1 × R_k  ≤  Z_k − Z_r  ≤  M_p − Z_r

整理即得

|y_k − y_r|  ≤  10 × (M_p − Z_r)

逆否命题成立：带外节点必然不入选。

这条不等式只用到“距离 ≥ 南北向分离”和“节点高程 ≤ 剖面最大值”两个平凡事实， 与 DEM 分辨率、投影方式、地形复杂度无关，因此适用范围很宽。

3.3 两种投影下的带宽

UTM 型 DEM（IPLAN = 1）

R_k = sqrt(ΔE² + ΔN²) ≥ |ΔN|

带宽直接取米：

B_p = 10 × (M_p − Z_r)        单位：m

经纬度型 DEM（IPLAN = 0）

节点纬度以弧秒存储，距离由 DISTLL（Helmert–Rainsford）给出。 真实子午线每度长度最小约 110574 m（赤道），向两极增大到约 111694 m。取

110000 m/度  <  110574 m/度

折算：

B_p = 10 × (M_p − Z_r) × 3600 / 110000     单位：弧秒

因为分母取小了，带宽算大了——只会多算、不会漏算，这是刻意的保守处理。

3.4 为什么只剪南北向

剖面（profile）是一列等东西向坐标、变南北向坐标的节点。循环沿列方向推进， 每个节点的 x 相同、y 递增。所以：

• 南北向坐标 y_k 与循环下标 INOD 一一对应，可以直接换算成“起止下标”；
• 东西向坐标在整列内是常数，剪不掉任何节点。

换句话说，剪枝不是“少算一点”，而是把“逐点判断”提升为“区间判断”： 整段落在带外的节点，一律不进内层。

四、实现：只跳过“不可能分支”

补丁位于 patch_aermap_calchc_band.py，对 sub_calchc.f 做三处插入。

4.1 循环前算带边界

C*           --- AERSCREEN-ZH PATCH: limit nodes to the lat band ---
             IF (IPLAN(IDEM) .EQ. 0) THEN
                CHC_BAND = 10.0D0*(MAXEL(IPROF)-AZELEV(IREC))
     &                   * (3600.0D0/110000.0D0)
                CHC_YR   = RLATSHFT
             ELSE
                CHC_BAND = 10.0D0*(MAXEL(IPROF)-AZELEV(IREC))
                CHC_YR   = YRSHFT
             END IF
             IF (DYM(IDEM) .GT. 0.0D0 .AND. CHC_BAND .GT. 0.0D0) THEN
                CHC_I1 = MAX(1, INT((CHC_YR-CHC_BAND-YBASE(IPROF))
     &                 / DYM(IDEM)) - 1)
                CHC_I2 = MIN(NODES(IPROF), INT((CHC_YR+CHC_BAND
     &                 -YBASE(IPROF)) / DYM(IDEM)) + 2)
             ELSE
                CHC_I1 = 1
                CHC_I2 = NODES(IPROF)
             END IF

三个细节：

• DYM 是行间距（弧秒或米），YBASE 是剖面首行坐标；两者都来自 DEM 索引， 不需要额外读取高程数据。
• 下标两侧各留 −1 / +2 的余量：INT() 是向零截断，且行列基址可能差 1， 余量用于吸收浮点求值与下标基准的不确定度。余量只影响“多算几个”，不影响结果。
• 若 M_p ≤ Z_r（受体高于整条剖面），带宽非正，代码退回全循环—— 这种情况下本来一个节点也不可能入选，属于可再优化的空间，但不改变正确性。

4.2 循环内只跳过重活

               CALL GETNOD(IPROF,INOD)

                JRECRD = JRECRD + 1
C**             --- AERSCREEN-ZH PATCH: skip out-of-band nodes ---
                IF (INOD .LT. CHC_I1 .OR. INOD .GT. CHC_I2)
     &              CYCLE NODE_LOOP

                ZNODE = GETDEM(IDEM,JRECRD)
                ...  DISTLL / 10% 判据 ...

注意保留了 GETNOD 与 JRECRD = JRECRD + 1：

• GETNOD 负责把节点坐标累加进 YNODE / XNODE，跳过它会改变后续节点的坐标值；
• JRECRD 是直接访问文件的记录号，跳过自增会错位读取。

所以被剪掉的只有三件事：GETDEM 随机读、DISTLL 大地线反算、判据比较。 这三件事对带外节点本来就不可能产生“入选”结果。

4.3 为什么输出能逐位不变

CALCHC 的输出只有 HEFF / HEFMAX / AZHILL(IREC)，而它们只在 ZNODE − AZELEV ≥ 0.1 × RDIST 成立时被更新。带外节点被证明恒不成立，因此：

• 更新序列完全相同；
• HEFMAX 的浮点比较顺序完全相同（只有入选节点参与比较）；
• HILLDBG 调试输出也完全相同（调试行只在 HEFF > HEFMAX 分支内写）。

换言之，这不是“结果接近”，而是同一段程序在等价输入上的同一条执行路径。

4.4 固定格式的坑

AERMAP 源码是 Fortran 固定格式（fixed-form），正文必须落在第 7–72 列， 超出部分会被静默截断。补丁脚本在写入前逐行校验 72 列限制， 避免“编译通过但语句被截断”的隐性错误。

五、另一处改动：NED 边界的 off-by-one（真实缺陷）

这是修 bug，不是优化，性质必须区分清楚。

5.1 症状

源点或受体恰好落在 DEM 瓦片边缘时，AERMAP 报：

NOT Inside DEM File

随后 AERSCREEN 静默回退为不考虑地形继续计算——用户看不到报错， 但结果已经变成“平坦地形”，评价等级可能被低估。

5.2 根因

sub_nedchk.f 判断点是否落在 NED/GeoTIFF 范围内时，行列上界用了

(dble(nRows(idem)-1))
(dble(nCols(idem)-1))

而 GeoTIFF 的像元中心坐标范围本应覆盖 nRows / nCols 个像元。 少算一行一列，导致最外圈像元被判为“域外”。

修正：

(dble(nRows(idem)))
(dble(nCols(idem)))

5.3 影响面

只影响落在瓦片最外一列/一行的点；域内点的判定与取值完全不变。 这也解释了为什么该问题长期只在“源点贴着瓦片边界”的项目里偶发。

六、实测与逐位回归

6.1 性能

测试条件：点源 + 复杂地形，25 km 探测半径，1 弧秒 DEM，10765 受体（粗 DEM 组为完整受体）。

全分辨率 DEM 下单剖面节点多（约 1900 个），剪枝收益稳定在约 5×； 粗分辨率 DEM 下单剖面节点数少，固定开销占比上升，收益回落到 2.9×（但仍与节点数近似成正比）。

6.2 一致性

以下文件在两版之间逐字节相同（cmp 无差异）：

• AERMAP.KIL（关键中间量）
• AERMAP.OUT
• MAPDETAIL（逐受体地形明细）
• aermap_srcloc.out（源高程）
• aermap_rcpt.out（受体高程）

下游：

• AERMOD 输入文件（CO/SO/RE/ME/OU 全部 pathway）逐字节相同；
• AERMOD .PLT 浓度逐值相同。

6.3 为什么可以断言“逐位”

输出是由 HEFF/HEFMAX/AZHILL 经固定格式写出的确定性文本。 只要这三个量的更新序列一致，格式化结果就一致； 而更新序列只取决于“哪些节点入选”，这一点已由第三节的不等式严格保证。

七、被排除的假设（工程方法）

性能问题最容易误判在 I/O 或编译器上。以下假设都做了实测排除：

排除这些之后，问题才收敛到“算法复杂度”这一层。

八、合规性：为什么可以改 AERMAP，但不能碰扩散内核

项目的原则是：合规性优先，不动扩散内核。

• AERMAP 是地形预处理器，负责把 DEM 转成受体高程与临界山高， 不参与浓度计算；AERMOD / AERSCREEN 的扩散方程、参数、默认值一律未改。
• 剪枝补丁不改变任何输出数值（逐字节验证），在数学上等价于“跳过不可能成立的分支”。
• 边界 off-by-one 修的是明确的实现错误（少算一行一列），修后域内取值不变。
• 所有改动以补丁脚本形式提供，源码随安装包分发，可自行 --check / --revert / 重新编译复算。

这套组合（源码 + 可复现补丁 + 逐位回归报告）是改动可以被审计的前提。

九、复现步骤

# 1) 查看补丁状态
python3 patch_aermap_calchc_band.py  <aermap源码目录> --check
python3 fix_aermap_ned_corners.py    <aermap源码目录> --check

# 2) 一键应用（幂等）
bash apply_all.sh <aermap源码目录>

# 3) 编译（示例：Windows 交叉编译）
bash build_aermap_one.sh x86_64-w64-mingw32-gfortran \
     aermap.exe <已打补丁的源码目录>

# 4) 逐位比对两次运行的输出
cmp AERMAP.KIL  AERMAP_ref.KIL
cmp AERMAP.OUT  AERMAP_ref.OUT
cmp aermap_rcpt.out aermap_rcpt_ref.out

十、局限与后续

• 复杂度仍是 O(受体 × 剖面 × 节点)，改进的是常数因子；受体规模再上一个量级，仍需从数据结构入手。
• 当地形极陡、M_p − Z_r 很大时，带宽变宽，剪枝收益下降； 但最坏情况只退化到原始全循环，不会变慢（多出的只是几条整数运算）。
• 曾尝试的两条路线已实测撤回：节点索引缓存 +16%、目录缓存 +13%， 均不如纬度带剪枝，代码保留在仓库中供参考。
• 可继续探索：剖面级预筛（先按 M_p 排序）、受体级多线程、SIMD 化 GPNHRI。

附录

附录 A：符号表

附录 B：补丁清单

两脚本均支持 --check / --revert，幂等； CALCHC 补丁写入前校验 Fortran 固定格式 72 列。