60 GHz氧气吸收与毫米波大气衰减Matlab复现

简介：面向通信工程、电磁波传播与毫米波技术课程学习及科研入门需求，这份代码围绕ITU-R P.676-11建议书复现无线电波在大气气体中的衰减曲线，帮助读者理解氧气、水汽等大气吸收对传播损耗的影响，并可在MATLAB环境中直接运行得到对应图像。压缩包共6个文件，以5个.m脚本为主，另有1个txt说明文件，整体约3KB，体积轻量但结构集中，便于快速查阅和改参。代码每行均有注释，main入口可输出复现图像，适合用于课程作业、论文示意图、算法验证及二次开发前的原理梳理。对初次接触大气衰减模型的读者而言，既能省去从公式到程序的重复搭建，也可通过逐行阅读掌握吸收系数计算、参数设置与曲线绘制之间的对应关系。目前已有418人学习下载，说明其在相关学习和入门场景中具有一定参考价值。

1. 60 GHz 附近的那条“断崖”：毫米波大气衰减到底被谁吃掉了

同样一条 1 km 的水平路径，架在 35 GHz，大气吸收损耗只有零点几 dB；挪到 60 GHz，直接跳到十几 dB。做 matlab 复现大气损耗曲线时，这个落差最容易让人误判成代码写错了——反复查系数、换单位、重画图，结果发现代码是对的，是氧气在吃能量。电磁波的大气衰减里，毫米波大气衰减占绝对主导的就是氧气和水汽的吸收谱线，它决定了 5G 毫米波、卫星链路、射电天文窗口频率的选型边界。把这条曲线算对，等于拿到了一张“哪些频点能通、能通多远”的地图。这份基于 P.676-11 建议书的 matlab 代码下载资源，做的正是这件事：把建议书里那份几十条谱线的表格，变成一条可以直接查、直接用的损耗曲线。

真正值得搞懂的不是“怎么画图”，而是“为什么曲线长这样”。为什么 60 GHz 有个又宽又深的坑，为什么 22 GHz 和 183 GHz 各有一个尖峰，为什么同一频点下雨天和晴天的差距能到一个数量级，为什么有人在 118 GHz 附近算出负衰减。这些问题都指向同一套机制：复折射率虚部、谱线形状因子、连续谱修正、以及沿路径的分层积分。下面按“物理量怎么来 → 单条谱线怎么展开 → 多条谱线怎么叠加 → 路径上怎么积分 → 复现时怎么翻车”的顺序拆。

2. 逐线求和：0.1820f·N'' 这个式子为什么能把吸收算准

大气气体衰减的逐线法（line-by-line）看上去只有一个核心式子：

γ = 0.1820 f · N''(f) （dB/km，f 单位 GHz）

所有复杂度都藏在 N'' 里。搞懂这个 0.1820 是怎么来的，比背公式有用得多，因为它是后面所有单位换算错误的源头，也是判断“代码算出来量级对不对”的唯一标尺。

2.1 复折射率虚部：吸收损耗的物理出处

电磁波在介质中传播，折射率写成复数 n = n' − j n''。n' 决定相速度（折射、弯曲、时延），n'' 决定衰减。电场沿 z 传播时幅度按 exp(−k₀ n'' z) 衰减，k₀ = 2πf/c 是自由空间波数。功率正比于场强平方，所以功率衰减指数是 2k₀n''z。

把功率衰减换算成分贝：A = 10·log₁₀(e^(2k₀n''z)) = 8.686·k₀·n''·z。这里有个容易漏的细节：无线电工程里习惯用折射率 N = (n−1)×10⁶ 表示，所以 n'' 对应的是 N''×10⁻⁶，不是 N'' 本身。代入 f 用 GHz、z 用 km：

γ = 8.686 × (2π f×10⁹ / 3×10⁸) × N''×10⁻⁶ × 10³ = 8.686 × (2π/3) × 10⁻² × f × N'' ≈ 0.1819 f N''

四舍五入就是建议书里的 0.1820。这个推导值得自己手推一遍，因为一旦你把 N'' 当成无量纲小数、或者把 f 用 MHz 代进去，曲线形状对、数值整体偏 10³，而 semilogy 画出来依然“像那么回事”，这种错误在图上几乎看不出来。

顺带解释了另一个常见困惑：为什么很多资料里 N'' 小到 10⁻⁴ 量级，衰减却有几十 dB/km。因为 f 是 GHz 量级，0.1820×60×10⁻⁴ 只有 0.001 dB/km——等等，这说明 10⁻⁴ 量级的 N'' 对应的是极小的衰减。真正的强吸收线附近 N'' 能到 0.1 以上，乘上 0.1820×60 就是 1 dB/km 起步。所以判断代码对不对，先看 60 GHz 附近 N'' 有没有到 0.1 量级，这是最快的体检指标。

2.2 谱线强度 S_i：θ³ 和 exp(a₂(1−θ)) 分别管什么

N'' 由两部分组成：一堆离散谱线的贡献叠加，加上一个缓变的连续谱底。每条谱线的强度写成（氧气）：

S_i = a₁ × 10⁻⁷ × p × θ³ × exp[a₂(1 − θ)]

水汽的形式类似，只是把 p 换成水汽分压 e，10⁻⁷ 换成 10⁻¹，θ 的幂次从 3 变成 3.5。其中 θ = 300/T 是归一化温度，p、e 单位都是 hPa。

这两个温度因子来源不同，理解它们能避免一类很隐蔽的错误。exp[a₂(1−θ)] 来自玻尔兹曼分布：谱线强度正比于处于低能级的分子数，能级差除以 kT 就得到指数项，a₂ 是负值（氧气谱线大约在 −6 到 −8 之间），温度升高、θ 变小、(1−θ) 变正、指数项趋近 1 甚至更大——高温让更多分子爬到低能级，吸收增强。θ³（水汽 θ^3.5）则来自跃迁偶极矩和配分函数的温度定标，是量子力学推导的结果，不是拟合出来的经验指数。

压力 p 出现在正比位置，这一条决定了气压随高度下降时吸收会同步变弱，也是后面做高度积分的基础。注意这里的 p 是干空气压力，不是总压；如果直接把总压代进去，在湿度大的低空会引入百分之几的偏差，链路预算里这类偏差会被误认为“模型不准”。

2.3 向量化实现：把 44 条谱线一次性铺成矩阵

建议书里氧气有 44 条谱线（延伸到 800 GHz 以上），水汽有 35 条。如果写双重 for 循环，频率扫描 1000 个点就是 79000 次内层计算，在 matlab 里慢得让人想砸键盘。正确做法是广播成矩阵，一次算完。

function [g_o, g_w] = gas_specific_atten(f, p, e, T) % f : 频率(GHz)，列向量 N×1 % p : 干空气压力(hPa)，标量或与 f 同长 % e : 水汽分压(hPa) % T : 绝对温度(K) % g_o / g_w : 氧气 / 水汽比衰减，dB/km theta = 300 ./ T; [fo, a1, a2, a3, a4, a5, a6, a7] = oxygen_table(); % 44 行谱线参数 [fw, b1, b2, b3, b4, b5] = vapour_table(); % 35 行谱线参数 F = f(:); % N×1 fo = fo(:).'; % 1×44 fw = fw(:).'; % 1×35 % ---- 氧气谱线 ---- So = a1 .* 1e-7 .* p .* theta.^3 .* exp(a2 .* (1 - theta)); dfo = p .* theta.^(0.8 - a4) .* (a3 + a5 .* theta); % 压力展宽宽度 Dfo = 0.535 .* dfo + sqrt(0.217 .* dfo.^2 + 2.1316e-12 .* fo.^2 ./ theta); dlt = (a6 + a7 .* theta) .* 1e-4 .* (p + e) .* theta.^0.8; % 干涉/混合项 No_line = sum(So .* shape(F, fo, Dfo, dlt), 2); % 沿谱线维求和 % ---- 干空气连续谱 ---- d = 5.6e-4 .* (p + e) .* theta.^0.8; No_cont = F .* p .* theta.^2 .* ( ... 6.14e-5 ./ (d .* (1 + (F ./ d).^2)) + ... 1.40e-12 .* F ./ (1 + 1.9e-5 .* F.^1.5) ); g_o = 0.1820 .* F .* (No_line + No_cont);

shape是谱线形状因子的独立函数，第 3 章展开。这里要强调sum(..., 2)的第二维参数：矩阵是 N×44，谱线在列方向，必须沿第 2 维求和；写成sum(x)会得到 1×44 的向量，画图时 matlab 不报错，只是曲线完全不对——这是复现过程中最高频的静默失败之一。

dfo那一行里的theta.^(0.8 - a4)也值得停一下：a₄ 是每条谱线自己的温度指数，不是全局常数，所以必须按列广播，写成theta^(0.8 - a4)（去掉点）在某些 matlab 版本会直接报错，在另一些版本会悄悄只取第一个元素。

3. 谱线形状因子：为什么每条线要在 ±f 两个位置各算一次

谱线不是一根针。真实气体里分子不断碰撞、不断热运动，能级宽度被展宽，吸收峰变成一个有一定宽度的钟形。形状因子 F_i 描述的就是这个钟形的具体形状，它决定了 60 GHz 那个坑到底有多宽、22 GHz 的尖峰到底有多尖。

3.1 洛伦兹线型与镜像项

建议书给的形状因子是：

F_i = (f / f_i) × [ (Δf − δ(f_i − f)) / ((f_i − f)² + Δf²) + (Δf − δ(f_i + f)) / ((f_i + f)² + Δf²) ]

第一项是标准的洛伦兹线型，中心在 f_i，半宽 Δf。第二项中心在 −f_i，看着像是“负频率上的幽灵”，其实是 Van Vleck–Weisskopf 修正的必然产物：真实的谱线响应函数必须满足某种对称性，否则在远离共振点的地方会算出不符合物理的结果。少了这一项，低频段（比如 1–10 GHz）的衰减会明显偏小，而这个频段恰好是很多微波链路的工作区，偏差会直接体现在链路预算里。

注意前面的 f/f_i 因子。它让 F_i 在 f 远离 f_i 时不会衰减得太快，保证远翼贡献被正确计入。这也是为什么 183 GHz 的水汽线会影响到 300 GHz 以上的衰减——远翼不是零。

3.2 Δf 的合成：压力展宽与多普勒展宽怎么合并

线宽 Δf 由两种展宽机制合成：

Δf = 0.535·Δf_p + sqrt(0.217·Δf_p² + 2.1316×10⁻¹²·f_i²/θ)

Δf_p 是压力（碰撞）展宽宽度，正比于气压，低空占主导，线型是洛伦兹型；后面那一项对应多普勒展宽，正比于频率、反比于温度的平方根，高空占主导，线型是高斯型。两者量级相当时，用这种半经验的平方和开方形式做 Voigt 近似，比精确算 Voigt 卷积快几个数量级，精度损失在工程上可以接受。

这个设计取舍很典型：建议书要的是“全球任何站点、任何高度都能算”，所以宁可牺牲一点极端条件下的精度，也要保证公式在 p 从 1013 hPa 降到 0.1 hPa 的整个范围内数值稳定。理解这一点，就不会在“为什么不用精确 Voigt”上钻牛角尖。

系数 0.535 和 0.217 是拟合出来的，不是推导的。这类常数在复现时最忌讳手打，抄错一位，整条曲线会整体抬升或压低几十个百分点，而且形状看不出异常。

3.3 干涉项 δ：不是所有谱线都对称

δ 是谱线混合（mixing）参数，描述碰撞过程中谱线之间的相互干扰，量级在 10⁻⁴ 乘气压。它让谱线形状略微不对称——一侧比另一侧宽一点。氧气谱线用 a₆、a₇ 两个参数算 δ，水汽谱线在这个版本里取零。

写代码时要注意符号：δ 出现在分子上的形式是 (Δf − δ(f_i − f))，δ 为正时，(f_i − f) 在 f < f_i 一侧为正，会减小有效宽度；f > f_i 一侧 (f_i − f) 为负，会增宽。符号抄反了，谱线两侧的宽度关系就颠倒，在单条线上看不出来，在 60 GHz 那种多条线紧密重叠的区域会让整个吸收带形状明显走样。

function F = shape(f, fi, Df, delta) % f: N×1；fi / Df / delta: 1×M，广播成 N×M df1 = fi - f; % 共振点左侧距离 df2 = fi + f; % 镜像项 F = (f ./ fi) .* ( ... (Df - delta .* df1) ./ (df1.^2 + Df.^2) + ... (Df - delta .* df2) ./ (df2.^2 + Df.^2) ); % 返回 N×M，调用方按第 2 维求和 end

这段的数值安全性来自分母：Df 最小也有 MHz 量级（0.1 hPa 高空的氧气线宽约 0.6 MHz），永远不为零，所以 f 恰好等于 f_i 时也只是分母取最小值，不会出现 Inf。真正会产生 NaN 的地方在别处，第 5 章会讲。

4. 连续谱与分层积分：为什么不能拿单点比衰减乘距离

算出一条 γ(f) 曲线，很多人的第一反应是“衰减 = γ × 距离”。在 1 km 的实验台上这么算没问题，在一条几十公里的斜路径上，这个假设会带来百分之几十的误差。原因有两个：连续谱项的存在，以及大气参数沿高度剧烈变化。

4.1 干空气连续谱与水汽连续谱：谱线之间的“地板”

谱线之间不可能是真空。氮气、氧气的碰撞会产生一个缓变的宽带吸收底，这就是干空气连续谱项 N''_D。它的结构里有两个分式：前一项在低频段占主导，用 d 这个“参考宽度”做平滑，d 正比于 (p+e)·θ^0.8；后一项在 100 GHz 以上逐渐抬头，形如 f/(1+1.9×10⁻⁵f^1.5)，是一种经验拟合的远翼叠加。

水汽连续谱更麻烦，它包含三项：自加宽（ρ² 正比，湿空气里水汽分子互相碰撞）、水汽与干空气的碰撞加宽（ρ·e 正比）、以及一个远翼项。三项都是温度相关的幂函数，系数直接抄建议书。

这两项的实际意义在于：在“窗口频点”——比如 35 GHz、94 GHz、140 GHz——衰减几乎全部来自连续谱，谱线贡献很小。如果只算了谱线、漏了连续谱，窗口频点的衰减会算成接近零，曲线在窗口处会塌到坐标轴上去，形状非常假。反过来，如果连续谱系数用错，窗口频点的误差能到 30% 以上，而峰值频点反而看起来正常，这类错误最难通过肉眼发现。

顺带提一句水汽密度与水汽分压的换算：ρ(g/m³) = 216.7·e/T，e = ρT/216.7。建议书不同小节用的输入变量不一样，逐线法用 e，简化式常用 ρ，中间切错一次，湿度直接差一个数量级。

4.2 沿路径积分：γ 是高度的强函数

γ 依赖 p、e、T，这三个量都随高度变化。标准大气下，压力按 p(h) = 1013.25(1−2.25577×10⁻⁵h)^5.2559 衰减，温度按 6.5 K/km 的递减率下降（到 11 km 后等温），水汽密度按指数衰减，标高只有约 2 km。

所以真实的斜路径衰减是：

A = ∫₀^L γ(h(l)) dl

直接数值积分当然可以，取 100 层、用梯形法，几毫秒就算完。但建议书提供了一个更省事的近似：因为 γ 随高度的形状跟频率有关，可以把它拆成“地面比衰减 × 等效高度”的形式：

A = (γ_o·h_o + γ_w·h_w) / sin θ （仰角 θ ≥ 5°）

h_o 和 h_w 分别是氧气和水汽的等效高度。氧气等效高度在 6 km 左右，但在 60 GHz 附近会被 t₁ 这一项显著抬高——因为 60 GHz 的吸收在低空就已经饱和，再往上的贡献占比变大，等效高度自然上移。水汽等效高度只有 1.7 km 上下，σ_w 这个 sigmoid 项描述的是水汽在垂直方向上的“塌缩”程度，气压越低、水汽越集中。

这个近似的设计动机很清楚：链路预算工程师要的是“输入地面温湿压，输出总衰减”，而不是一条积分曲线。把积分预先解析化，用几个闭合公式代替数值积分，代价是仰角低于 5° 时精度下降（那时地球曲率和路径长度修正都要另外处理）。复现代码时，两条路都留着，用数值积分去校验等效高度近似，是很好的自检手段。

4.3 从单点到全图：一条完整曲线的生成流程

把上面的部件串起来，一条完整的大气损耗曲线生成流程是：

1. 定义频率向量，1–350 GHz，步长 0.1 GHz，用列向量避免转置错误；
2. 给定地面气象参数：p = 1013.25 hPa、T = 288.15 K、ρ = 7.5 g/m³，换算成 e；
3. 调逐线函数，得到 γ_o、γ_w，相加得总比衰减；
4. semilogy 画图，横轴频率、纵轴 dB/km；
5. 换一组参数（比如 ρ = 1 g/m³ 的干燥天气）重画，两条线叠在一起对比。

f = (1:0.1:350).'; p = 1013.25; T = 288.15; rho = 7.5; e = rho .* T / 216.7; % 水汽密度 -> 水汽分压 [g_o, g_w] = gas_specific_atten(f, p, e, T); gamma = g_o + g_w; semilogy(f, gamma, 'LineWidth', 1.2); hold on; grid on; xlabel('Frequency (GHz)'); ylabel('Specific attenuation (dB/km)');

画出来应该看到：22 GHz 一个尖峰（约 0.2 dB/km），60 GHz 一片宽平台（十几 dB/km），118 GHz 一个次峰，183 GHz 一个更高的尖峰（几十 dB/km），中间 35、94、140、220 GHz 是明显的低谷。如果 183 GHz 的峰没有 60 GHz 的高，八成是水汽谱线强度那一项的指数用错了。

5. 复现时最容易翻车的几处：单位、发散与验证手段

同一份建议书，不同人复现出来曲线能差一倍，绝大多数不是公式理解错，而是踩在这几个坑里。下面几条是从反复对比中攒下来的血泪经验。

5.1 单位与量纲的三处高频翻车

现象：曲线整体平移，形状完全正确，峰值位置也对，就是数值差一个常数倍。

原因之一是把 p 和 e 的单位搞混。逐线法里 p、e 全部用 hPa，谱线强度式里的 10⁻⁷ 和 10⁻¹ 这两个因子就是为了配合 hPa 标定的。如果换成 Pa（乘 100），衰减会整体放大 100 倍。

原因之二是水汽密度和水汽分压混用。简化式那张表用的是 ρ(g/m³)，逐线法用的是 e(hPa)，两者差一个 216.7/T ≈ 0.75 的系数，在 288 K 时约等于 0.75，看似不大，但湿度项在 22 GHz 和 183 GHz 附近的贡献是主导项，0.75 倍的偏差足以让峰值对不上。

原因之三是温度用了摄氏度。θ = 300/T 里的 T 必须是开尔文。用 15 代替 288.15，θ 变成 20，θ³ 就是 8000 倍，所有谱线强度直接爆掉，曲线会变成一片几十万 dB/km 的怪物。

解决：函数入口处写三段断言，把量纲钉死。

assert(T > 100 && T < 400, '温度必须是开尔文'); assert(p > 0 && p < 1200, '压力单位应为 hPa'); assert(e >= 0 && e < 100, '水汽分压单位应为 hPa');

这三行看着土，但能省掉大量“对着图发呆”的时间。

5.2 频率网格与数值发散

现象：曲线在某些频点出现 NaN 或者突然下探到 10⁻³⁰。

原因：频率向量从 0 开始。f = 0 时，γ = 0.1820×0×N'' = 0，semilogy 会把 0 画到坐标轴外面去，看起来像断线。更隐蔽的是 f 极小但非零时，谱线形状因子里的 f/f_i 趋于 0，整个 N'' 也趋于 0，log 坐标下同样是断崖。

另一个来源是压力或温度传成了向量且长度不匹配。matlab 的隐式扩展在维度兼容时会悄悄广播，比如 f 是 1000×1、p 是 1×1000，结果会变成 1000×1000 的矩阵，内存瞬间吃满，然后返回一个维度错误——报错信息指向g_o，跟真正的错误位置差得很远。

解决：频率向量从 1 GHz 起步，不要从 0 起；所有输入在函数内部统一(:)化成列向量；关键中间量用assert(all(isfinite(x)))卡一道。

5.3 结果验证：怎么确认曲线没有 bug

现象：代码跑通了，图画出来了，但不确定对不对。

原因：大气衰减没有“标准答案”可以直接比对，建议书本身也只给了有限几个校验点。

解决：用三个物理约束做交叉验证。第一，60 GHz 附近在标准地面气象条件下，比衰减应该落在 10–15 dB/km；第二，35 GHz 和 94 GHz 这两个窗口频点，应该落在 0.1–0.3 dB/km 量级；第三，把水汽密度设成 0，22 GHz 和 183 GHz 的峰应该基本消失，而 60 GHz 的峰几乎不变——因为那是氧气的。这三条只要有一条不满足，就说明某一部分的系数或者单位出了问题，比逐行读代码快得多。

第四个验证是数值积分与等效高度近似互相校验，两者在仰角 30° 以上应该相差不到 5%。这个交叉验证尤其有价值，因为两条实现路径几乎不共享代码，能同时错的概率很低。

6. 逐线法还是简化式：精度、算力与场景的取舍

建议书里其实给了两套算法。逐线法精度高，但要读两张几十行的系数表，还要处理展宽、干涉项、连续谱；简化式只有几个闭合公式，输入 f、p、T、ρ 就能出结果，代价是适用范围被限制在 1–350 GHz，且在某些频点的误差能到 10% 以上。

这个取舍值得展开想一下。逐线法的“贵”不在于算力——44 条谱线在 matlab 里向量化之后，1000 个频点也就毫秒级。真正贵的是“可移植性”：那两张表有几百个数字，任何一个抄错都可能导致局部失真，而且失真藏在曲线里，肉眼极难发现。简化式把几百个数字压成十几个拟合参数，抄错的风险低得多，代价是精度。这也是为什么很多商用链路预算工具默认跑简化式，只在做频谱规划、干扰分析这类需要高保真的场景才切到逐线法。

复现的时候我建议两条都实现，用逐线法当基准去标定简化式。这样既能拿到高精度曲线，又能在需要实时计算的场合有一份轻量实现，还能顺便验证自己对谱线参数的理解——如果两套算法在你的输入下差得离谱，那一定有一边写错了，而它们几乎不可能同时错成一致的结果。

还有一个容易被忽略的取舍是频率分辨率。0.1 GHz 步长在 60 GHz 那片宽平台上够用，但在 22.235 GHz 的尖峰附近，峰顶可能被跳过，画出来的峰高偏低。正确做法是在已知谱线中心附近做局部加密，或者在画图前先按谱线表生成一组“必采频点”并进频率向量。这个细节不会影响链路预算的结论，但会影响你对“峰到底有多尖”的判断。

从那以后我每次做这类谱线复现，都会强制走一遍“三频点体检 + 两算法互校 + 量纲断言”的流程，哪怕只是改一个系数，也先跑一遍再画图。吃过一次亏：某个高湿条件下 183 GHz 的峰值莫名比 60 GHz 低，查了两天才发现是水汽谱线表少录了一行，而那一行恰好是 183 GHz 的主线——曲线整体看起来“挺像”，只有跟基准对比才暴露。这套流程连同谱线表、绘图脚本和校验函数，我整理成了一份可直接运行的 matlab 代码，放在下载包里，拿到之后先跑check_points.m看三个体检频点是否在预期区间，再去改参数探索不同温湿条件，会省掉很多对着图猜原因的夜晚。希望帮到你。

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