复杂地形电磁波多径传输建模：基于相对余隙判据的仿真方案

简介：电磁波在真实山地、城市环境中并非沿直线传播，大气折射、吸收、多径干涉与绕射共同决定了最终路径损耗。工程上通常以自由空间损耗为基准，再叠加环境效应修正。复杂地形传播建模的核心是先判断目标处于视距多径还是绕射区，而这依赖于相对余隙与第一菲涅尔区半径的比值。该决策逻辑能够避免盲目套用单一公式，为雷达对抗仿真、链路预算、电子战效能评估提供可靠依据。文中进一步结合SEKE算法与数字高程模型剖面抽取，给出了完整可落地的仿真实现框架。通过六个仿真模式，可对频率扫描、DEM分辨率敏感性及动态轨迹进行参数化分析，帮助工程师和研究者快速掌握复杂地形电磁波多径传输建模与仿真的工程实践。

1. 复杂地形电磁波多径传输建模：一份能直接落地的仿真方案

做雷达电子战仿真的人，迟早会被同一件事卡住：明明自由空间损耗算得顺手，一到山地、丘陵、城市这类复杂地形，路径损耗怎么算都不对，因为电波根本不是按直线走的。这份《复杂地形环境中电磁波多径传输建模与仿真》把问题拆得比较干净——大气折射、大气吸收、多径干涉、绕射四种效应各有一套算法，再用“相对余隙”这两个判据决定当前场景到底该跑哪个模型。它最值钱的地方不是某个公式，而是这套“先判断用哪个模型，再计算衰减”的传播总模型决策逻辑。适合做雷达对抗仿真、链路预算、电子战效能评估的工程师，也适合刚入手电波传播建模、想找个完整参考框架的研究生。下面按我拆这篇论文的思路，把模型、参数和坑逐个过一遍。

2. 传播机理与四类环境效应：折射、吸收、多径与绕射的工程化处理

2.1 自由空间损耗是基准，所有环境效应都是修正项

论文把传播损耗拆成三层叠加：自由空间损耗是基准，大气吸收是媒质修正，多径干涉或绕射是地形修正。自由空间损耗公式是电波传播建模里最基础的那一个：

L_f = 32.45 + 20·lg d(km) + 20·lg f(MHz)

这个公式用 dB 表示，d 是收发距离（km），f 是工作频率（MHz）。注意它的对数结构：距离每翻一倍，损耗增加约 6dB；频率每翻一倍，同样增加约 6dB。这背后是能量球面扩散的物理本质，与大气和地形无关，所以它是整个模型的“底板”。

工程上我习惯先把这条曲线拉出来看一遍，再叠加其它效应。如果自由空间损耗本身就比系统链路预算大出 20dB 以上，后面那些环境效应修正的意义就不大了，可以直接判定链路不达标，省掉大量仿真时间。论文里没有明确写这个预判步骤，但做工程级仿真时这是必备的快速筛选手段。

2.2 大气折射：用射线追踪把弯曲路径算出来

电磁波在对流层里传播时，路径并非直线，而是向下弯曲的曲线。原因是大气折射指数随高度递减，波前上部速度略快于下部，波前逐步倾斜。这里有个反直觉点：频率对折射路径影响极小，决定因素是折射指数随高度的梯度。

论文给出了射线追踪的积分式，用来描述电波传播距离 R 与高度 h 的关系。工程实现时不直接解这个积分，而是把大气分层，每层视为折射指数恒定，逐层套用 Snell 定律计算折射角。这个“逐层折射”的过程就是射线追踪。

大气折射指数模型采用 GJB1655 给出的指数率模型：

n(h) = 338.5 · exp(-0.1404·h)

其中 h 是相对于海平面的高度，单位 km。这个模型是我国军用标准中的平均大气折射指数，数值 338.5 对应海平面附近的折射指数修正值，-0.1404 是指数衰减率。我这里提醒一点：这个模型是“平均”状态，实际试验场区如果温度和湿度分布异常，最好用探空气球实测数据替换指数率模型，否则折射修正量会有偏差。

射线追踪的典型做法是：从辐射源位置出发，以初始仰角 θ₀ 开始，按高度步长 Δh 逐层推进，每层计算一次折射角和水平位移，直到到达目标高度或射出对流层。高度步长一般取 50~100m 比较稳，步长太大会放大折射角误差，太小则计算量翻倍，实时仿真场景下压力比较大。

2.3 大气吸收：10GHz 是忽略与计入的分界线

大气吸收由对流层中的氧分子和水蒸气分子对电磁波的谐振吸收引起。氧分子在 60GHz 附近有强吸收峰，水蒸气在 22GHz 和 183GHz 附近有吸收峰。论文给出的吸收损耗计算公式是沿传播路径的积分：

L_A = ∫(γ_O2(r) + γ_H2O(r)) dr

其中 γ_O2 和 γ_H2O 分别是氧和水蒸气的吸收系数，单位 dB/km，r 是传播路径。

工程上这条最重要的经验值：工作频率低于 10GHz 时，大气吸收可以忽略，直接置零；在 10GHz~40GHz（毫米波段），必须计入。这条分界线在论文中有明确说明，但实现时要小心——所谓“忽略”不是真的零，而是损耗相对其它项小到不影响结论。做链路预算时我一般留 1dB 裕量来吸收这个近似误差。

吸收系数与大气压力 P、温度 T、水蒸气密度 ρ 三个参数相关。仿真时如果没有实测气象数据，可以直接从配置文件读取标准大气参数，或者通过界面人工输入。论文处理得比较务实：大气参数走配置文件，这样换试验场区时不用改代码。

2.4 多径干涉与绕射：此消彼长的两种地形效应

这两种效应都源自地球形状和地表特征，但物理机制不同。多径干涉是直射波与地面（或海面）反射波在接收点叠加形成干涉场；绕射是电波越过障碍物进入阴影区的现象。论文特别强调了一个工程判断：在一次传播过程中，目标要么处于多径干涉区，要么处于绕射区，要么处于过渡区，三者必居其一，所以这两类模型要放一起研究。

多径干涉建模时只考虑两条路径——直射路径和一次地面反射路径。两条路径距离不同导致相位差，反射过程还会产生额外相移，叠加后形成干涉。干涉的结果是接收功率与距离的关系不再遵守二次方反比定律，会出现一系列波峰和波谷。

绕射按障碍物形状分为球形地球绕射、单刃峰绕射、单圆顶形障碍物绕射、非刃峰绕射、多峰绕射。其中球形地球绕射由地球曲率引起，其余由山峰或高大建筑物引起。算法实现上，标准的绕射算法来自 Blake 和 Kerr 的经典文献；后来林肯实验室在此基础上发展出 SEKE 算法，对多径干涉、球形地球绕射和多刃峰绕射做了更精细的处理。论文选择把 SEKE 作为多径与绕射模型的实现基础。

这里有个很重要的工程选择：对于复杂地形，不要试图用一个通用模型适配所有场景，而是准备多个子模型，让程序根据地形剖面参数自动选择。这正是论文传播总模型的核心思想，下一章详细展开。

3. 传播总模型与相对余隙：两个判据决定跑哪个子模型

3.1 为什么需要一套传播总模型

如果只看单个效应模型，每个都有现成算法，但放到真实地形上就会出问题：山峰可能正好挡在收发路径中间，地面反射可能被树木遮挡，球面绕射和刃峰绕射的适用条件完全不同。直接套单一模型，算出来的损耗可能偏离实测值十几 dB。

论文给出的解法是建立一个“传播总模型”——先通过数字地图提取收发点之间的二维地形剖面，找到剖面最高点，计算相对余隙，再依据相对余隙值选择具体的传播模型。这个思路本质上是一个决策树，比盲目套公式靠谱得多。

3.2 第一判据：视距多径与绕射的边界

关键参数是相对余隙 p，定义为传播余隙 δ 与第一菲涅尔区半径 F₁ 的比值。第一菲涅尔区半径的计算公式为：

F₁ = sqrt(λ · d₁ · d₂ / (a · k))

其中 d₁、d₂ 分别是路径最高点到发射点和接收点的距离，λ 为波长，a 为地球半径（约 6370km），k 为等效地球半径因子（标准大气下通常取 4/3）。

判定规则如下：

• 若 p > 1：余隙充足，绕射影响可忽略，只采用视距直射加反射的多径模型；
• 若 1/2 < p < 1：多径干涉和绕射都不能忽略，传播损耗取两者的加权平均；
• 若 p ≤ 1/2：进入第二判据，在“多峰绕射”和“地表绕射”之间做选择。

这里我把 p 看作一个非线性开关，它是随频率变化的——频率越高，F₁ 越小，同样的地形余隙对应的 p 值越大，系统越倾向于多径模型。这解释了为什么同一片地形，L 波段可能判定为视距多径，Ka 波段却可能进入绕射区。

加权平均的权重也有讲究。过渡区的传播损耗表达式为：

L = a_MP · L_M + (1 - a_MP) · L_D

其中 L_M 是多路径损耗，L_D 是绕射损耗，加权因子 a_MP = 2p - 1。注意这个线性内插在 p=0.5 时权重为 0，在 p=1 时权重为 1，物理含义是：余隙越小，越偏向纯绕射；余隙越大，越偏向纯多径。

3.3 第二判据：多峰绕射还是地表绕射

当 p ≤ 1/2 时，路径最高点已经明显切入菲涅尔区，需要进一步判别障碍物形态。论文给出的第二判据基于 h_M/F₁ 的值——h_M 是路径最高点到地面的高度，F₁ 仍是第一菲涅尔区半径：

• 若 h_M/F₁ < 1/4：浅绕射区，采用地表绕射模型；
• 若 h_M/F₁ > 1/2：深绕射区，采用多峰绕射模型；
• 若 1/4 ≤ h_M/F₁ ≤ 1/2：过渡绕射区，采用多峰绕射与地表绕射加权平均。

第二判据的过渡区加权公式为：

L = a_KS · L_K + (1 - a_KS) · L_S

其中 L_K 为多峰绕射损耗，L_S 为地表绕射损耗，加权因子 a_KS = 4·h_M/F₁ - 1。当 h_M/F₁ = 1/4 时权重为 0，等于 1/2 时权重为 1，同样是线性内插。

用这套双判据决策树，任何一条传播路径都能落到一个具体的子模型上，不会出现“不知道该用哪个公式”的尴尬。我把两套判据整理成一个参数速查表：

3.4 总损耗的叠加方式

选定子模型后，总路径损耗按下式叠加：

L_total = L_f + L_atmos + L_prop

其中 L_f 是自由空间损耗，L_atmos 是大气吸收损耗——当工作频率 f < 10GHz 时直接取 0，L_prop 是多径干涉或绕射带来的附加损耗。

这里要强调一个实现细节：L_prop 的物理含义是“相对于自由空间的附加衰减”，它可能为正值（衰减增强），也可能为负值（多径增益）。多径干涉在某些距离上会形成建设性叠加，使接收功率高于自由空间预测值，即 L_prop 为负。做链路分析时如果发现某段距离上损耗曲线突然“凹陷”，别急着怀疑代码，先检查是不是进入了多径干涉增强区。

3.5 频率和余隙的耦合关系

用这套模型做频率扫描时，能明显看到一个规律：同一地形剖面，频率升高会让相对余隙 p 变大，因为第一菲涅尔区半径与波长正相关——频率越高，F₁ 越小，同样的物理余隙显得“更充裕”。这意味着频率升高，系统可能从绕射主导切换到多径主导。如果仿真中发现路径损耗随频率非单调变化，先检查模型切换点是不是就在扫描区间内。

4. 数字地图二维剖面抽取：从 DEM 到相对余隙的关键一步

4.1 为什么剖面的准确性直接决定判据可信度

传播总模型的两个判据都依赖两个关键量：剖面最高点位置和高度。如果剖面提取不准，后续所有计算都可疑。论文采用的数字地图是格网型（Grid）数字高程模型 DEM——以固定采样间隔按矩形或正方形网格排列，主体数据是每个网格节点的高程值。

收发点之间的路径是地球大圆的一段弧，不是平面直线。所以不能直接把 DEM 网格丢进笛卡尔坐标算，要先在球面上求出收发点间的大圆路径，再沿路径做高程插值。论文给出的剖面抽取流程是：

1. 输入收发点经纬度及采样点数；
2. 根据收发点经纬度求出大地线距离；
3. 求出发射点相对接收点的方位角；
4. 沿大圆路径计算各采样点经纬度；
5. 用插值法求各采样点高程，得到二维剖面。

这个流程里第四步是最关键的——采样点经纬度通常不会恰好落在 DEM 网格节点上，必须通过插值获得高程。

4.2 双线性多项式内插的工程实现

论文采用双线性多项式内插法，原理是用最靠近插值点的四个已知数据点组成一个四边形，确定一个双线性多项式来内插待插点高程。函数形式为：

z = a₀ + a₁·x + a₂·y + a₃·x·y

四个系数 a₀~a₃ 由四个参考点 P₁(x₁,y₁,z₁)、P₂(x₂,y₂,z₂)、P₃(x₃,y₃,z₃)、P₄(x₄,y₄,z₄) 联立求解，然后将待插点坐标代入即可得到高程值。

实际写代码时，我一般先判断采样点落在哪个 DEM 网格单元内，取该单元四个顶点做内插。Python 实现大致是这样：

def bilinear_interp(dem, x, y): """ dem: 2D numpy array, 行对应纬度方向, 列对应经度方向 x, y: 待插点坐标, 单位与 DEM 网格间距一致 """ x0, y0 = int(x), int(y) # 网格左下角索引 # 边界保护: 防止索引越界 x1, y1 = min(x0 + 1, dem.shape[1] - 1), min(y0 + 1, dem.shape[0] - 1) dx, dy = x - x0, y - y0 z00 = dem[y0, x0] z10 = dem[y0, x1] z01 = dem[y1, x0] z11 = dem[y1, x1] # 双线性内插: 先沿 x 方向, 再沿 y 方向 z0 = z00 * (1 - dx) + z10 * dx z1 = z01 * (1 - dx) + z11 * dx z = z0 * (1 - dy) + z1 * dy return z

这个实现与论文公式等价——a₀ = z00，a₁ = z10 - z00，a₂ = z01 - z00，a₃ = z00 - z10 - z01 + z11，展开后就是上面的双线性形式。注意 DEM 的坐标系方向，我这里约定 row 方向是纬度（y），col 方向是经度（x），如果你的数据是反的，内插公式要相应调整。

采样点数直接影响剖面分辨率。采样点太稀疏会漏掉真实最高点，把单刃峰判成平缓地形，导致误选模型；采样点太密则计算量大，实时仿真帧率掉得厉害。论文的仿真系统提供了一种模式专门研究“路径损耗随采样点数变化”，默认值可以由收发距离和 DEM 网格间距推算——一般取收发距离除以 DEM 网格间距再乘以 1.5，保证每个 DEM 格网至少被 1~2 个采样点覆盖。

4.3 从剖面到相对余隙的计算链

拿到剖面后，先找到剖面最高点 M，M 到收发直射线 AB 的垂直距离就是传播余隙 δ。这里有一个几何细节容易被忽略：直射线 AB 要考虑地球曲率，即等效地球半径修正，不然余隙会偏大。论文式(5)里引入了等效地球半径因子 k，标准大气下取 k=4/3。

余隙与第一菲涅尔区半径的比值 p = δ/F₁ 就是我们需要的相对余隙。我一般把这个计算封装成独立函数，输入为剖面数组和频点，输出为相对余隙和剖面最高点位置，方便在多个仿真模式里复用。

5. 仿真系统实现与经典避坑：从环境模型到加权因子的四个坑

5.1 系统架构：环境模型与效应模型分离

论文提出的系统设计分成两个独立模块，这个架构值得借鉴。环境模型封装试验场区的大气参数、地表形状和地表电特性数据，供查询使用；多径传输效应模型利用环境模型提供的数据计算衰减。大气参数可以通过界面直接输入或读取配置文件，地形特征参数则从数字地图（DTED 等）处理后获得。

开发环境是 Windows XP + Visual Studio.NET 2003 + Qt Commercial 4.3.3。Qt 负责 GUI 和交互，VS 负责核心计算模块。系统共支持六种仿真模式，这个扩展性和 Qt 的信号槽机制关系很大——每个模式的参数面板独立，仿真结果用 OpenGL 或 Qwt 类库做可视化。

5.2 四个典型的坑及解决方案

坑一：相对余隙恰好落在边界值附近，模型来回跳

现象：仿真时固定收发位置只修改频率，路径损耗曲线在某个频点附近出现明显的“台阶”或抖动，损耗值突变超过 10dB。

原因：频率连续变化导致相对余隙 p 从 1.0 附近穿越模型切换边界，系统在“多径模型”和“加权模型”之间来回切换，而两种模型计算出的损耗值本身存在跳变。

解决：在判据里加迟滞区间（hysteresis），比如 p > 1.05 才切换为纯多径模型，p < 0.95 才切换为加权模型，中间区域保持原模型不变。这个技巧在工程仿真里很常用，虽然论文没写，但实际做频率扫描时几乎必踩。

坑二：加权因子线性内插在极端地形下失真

现象：在地形剧烈起伏区域，加权模型算出的路径损耗明显偏离实测值，偏差集中在过渡区。

原因：a_MP = 2p - 1 是一个线性近似，前提假设是两种传播机制在过渡区的贡献随余隙线性变化。但真实场景中，绕射损耗随余隙的变化率远高于线性，尤其当障碍物形态接近尖峰时。

解决：将过渡区细分为多个子区间，每个子区间单独标定权重曲线，或者直接改用实测数据拟合权重函数。如果没有实测条件，至少把过渡区的仿真结果标注为“低置信度”，不做工程决策依据。

坑三：10GHz 以下大气吸收直接置零，但湿度极端时失效

现象：工作频率 8GHz，试验区湿度很高（水蒸气密度超过 20g/m³），实测损耗比仿真结果高 2~3dB。

原因：论文的规则是 f < 10GHz 时大气吸收损耗直接取 0，这个规则在标准大气条件下误差可以接受，但极端高湿条件下水蒸气吸收在 8GHz 已经不可忽略。

解决：程序里保留完整的吸收系数计算模块，只在界面层做“快速模式”和“精确模式”的切换。快速模式按论文规则忽略，精确模式关闭忽略开关，直接积分计算。做链路预算时用精确模式，做快速参数扫描时用快速模式。

坑四：DEM 边界区域的插值伪影导致剖面异常

现象：收发点接近 DEM 数据范围边缘时，剖面出现异常尖峰或凹陷，相对余隙判据被误导。

原因：双线性内插在数据边界外没有定义，采样点落在 DEM 边缘外时程序自动取边界值（clamp），导致剖面在边界处出现“悬崖”形状，被误判为山峰。

解决：剖面抽取前先检查所有采样点是否都在 DEM 覆盖范围内，如果有越界采样点，要么裁剪收发端距离，要么对边界外区域做地形外推（常见做法是填充平均高程值或最低高程值）。我在实现里更倾向直接裁剪，外推地形反而会引入不确定的人为地形。

5.3 六种仿真模式的实质区别

论文的六种模式本质上是对“参数扫描维度”的划分：固定收发端扫路径损耗、扫频率、扫采样点数，以及发射端固定接收端做三类运动——水平运动、垂直运动、沿预定轨迹运动。从实现角度，前三种是静态计算，后三种是动态计算——每帧更新接收端位置，重新抽取剖面并计算损耗。

动态模式下性能优化很关键。剖面抽取和余隙计算可以预计算固定部分（发射端到地形最高点的部分），只需要每帧更新接收端到最高点的那段路径。如果发射端固定而接收端沿轨迹运动，这个优化能把计算量降一半左右，在 Qt 的定时器驱动下可以保持流畅更新。

6. 六种仿真模式的工程运用：从链路预判到动态轨迹测算

6.1 模式选择的工作逻辑

我拿到一个新的试验场区，通常先跑两种模式摸底。先跑模式一，收发端固定的路径损耗值，确认基础链路余量；再跑模式二，路径损耗随频率变化，看频段内是否有模型切换点。如果切换点正好落在工作频率附近，说明这个试验场区对频率特别敏感，需要进一步细察。

模式三是比较特殊的参数敏感性分析——路径损耗随采样点数的变化。它直接反映数字地图分辨率对结果的影响程度。如果采样点数从 100 增加到 500 时损耗变化超过 3dB，说明当前 DEM 分辨率不足以支撑这个距离段的准确计算，需要换更高精度的地形数据或者加密采样。

动态模式的选择逻辑是：

6.2 动态轨迹模式的实现方式

模式六（预定轨迹运动）最有实战价值。实现时用插值把轨迹离散成一串位置点，每个位置点都做一次完整的剖面抽取和损耗计算。轨迹点数建议取 200~500 个，太少轨迹不够平滑，太多计算时间太长。输出结果是路径损耗随轨迹位置变化的曲线，配合等高线图一起看，能直观判断哪些航段处于绕射阴影区、哪些航段有多径增强。

6.3 一个值得记住的验证习惯

这套模型最后能否用于工程决策，取决于一件事：是否用实测数据标定过。论文里以实装对抗演练为背景做了仿真试验，但我建议你复现这套模型时，先用一小段简单地形做校验——找一个平整场地，实测自由空间加地面反射的损耗值，调多径模型参数直到与实测吻合，再去跑复杂地形。

从那以后我每次搭这类仿真链路，都强制自己先做一遍“单点标定”再上复杂场景，模型切换边界处的参数也每次单独核对。这个习惯帮我挡掉了不少后续的数据争议，希望帮到你。

本文还有配套的精品资源，点击获取