简介面向通信工程、电磁波传播与毫米波技术课程学习及科研入门需求这份代码围绕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/kmf 单位 GHz所有复杂度都藏在 N 里。搞懂这个 0.1820 是怎么来的比背公式有用得多因为它是后面所有单位换算错误的源头也是判断“代码算出来量级对不对”的唯一标尺。2.1 复折射率虚部吸收损耗的物理出处电磁波在介质中传播折射率写成复数 n n − j n。n 决定相速度折射、弯曲、时延n 决定衰减。电场沿 z 传播时幅度按 exp(−k₀ n z) 衰减k₀ 2πf/c 是自由空间波数。功率正比于场强平方所以功率衰减指数是 2k₀nz。把功率衰减换算成分贝A 10·log₁₀(e^(2k₀nz)) 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 换成水汽分压 e10⁻⁷ 换成 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×1fi / 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 正比于 (pe)·θ^0.8后一项在 100 GHz 以上逐渐抬头形如 f/(11.9×10⁻⁵f^1.5)是一种经验拟合的远翼叠加。水汽连续谱更麻烦它包含三项自加宽ρ² 正比湿空气里水汽分子互相碰撞、水汽与干空气的碰撞加宽ρ·e 正比、以及一个远翼项。三项都是温度相关的幂函数系数直接抄建议书。这两项的实际意义在于在“窗口频点”——比如 35 GHz、94 GHz、140 GHz——衰减几乎全部来自连续谱谱线贡献很小。如果只算了谱线、漏了连续谱窗口频点的衰减会算成接近零曲线在窗口处会塌到坐标轴上去形状非常假。反过来如果连续谱系数用错窗口频点的误差能到 30% 以上而峰值频点反而看起来正常这类错误最难通过肉眼发现。顺带提一句水汽密度与水汽分压的换算ρ(g/m³) 216.7·e/Te ρ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–350 GHz步长 0.1 GHz用列向量避免转置错误给定地面气象参数p 1013.25 hPa、T 288.15 K、ρ 7.5 g/m³换算成 e调逐线函数得到 γ_o、γ_w相加得总比衰减semilogy 画图横轴频率、纵轴 dB/km换一组参数比如 ρ 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/km60 GHz 一片宽平台十几 dB/km118 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 0semilogy 会把 0 画到坐标轴外面去看起来像断线。更隐蔽的是 f 极小但非零时谱线形状因子里的 f/f_i 趋于 0整个 N 也趋于 0log 坐标下同样是断崖。另一个来源是压力或温度传成了向量且长度不匹配。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 量级第三把水汽密度设成 022 GHz 和 183 GHz 的峰应该基本消失而 60 GHz 的峰几乎不变——因为那是氧气的。这三条只要有一条不满足就说明某一部分的系数或者单位出了问题比逐行读代码快得多。第四个验证是数值积分与等效高度近似互相校验两者在仰角 30° 以上应该相差不到 5%。这个交叉验证尤其有价值因为两条实现路径几乎不共享代码能同时错的概率很低。6. 逐线法还是简化式精度、算力与场景的取舍建议书里其实给了两套算法。逐线法精度高但要读两张几十行的系数表还要处理展宽、干涉项、连续谱简化式只有几个闭合公式输入 f、p、T、ρ 就能出结果代价是适用范围被限制在 1–350 GHz且在某些频点的误差能到 10% 以上。这个取舍值得展开想一下。逐线法的“贵”不在于算力——44 条谱线在 matlab 里向量化之后1000 个频点也就毫秒级。真正贵的是“可移植性”那两张表有几百个数字任何一个抄错都可能导致局部失真而且失真藏在曲线里肉眼极难发现。简化式把几百个数字压成十几个拟合参数抄错的风险低得多代价是精度。这也是为什么很多商用链路预算工具默认跑简化式只在做频谱规划、干扰分析这类需要高保真的场景才切到逐线法。维度逐线法简化式适用频率1–1000 GHz 以上1–350 GHz输入参数f、p、e、Tf、p、T、ρ系数规模约 500 个约 30 个典型误差参考基准窗口频点 5%–15%计算开销毫秒级向量化后微秒级主要风险系数抄录、展宽公式边界频点外推复现的时候我建议两条都实现用逐线法当基准去标定简化式。这样既能拿到高精度曲线又能在需要实时计算的场合有一份轻量实现还能顺便验证自己对谱线参数的理解——如果两套算法在你的输入下差得离谱那一定有一边写错了而它们几乎不可能同时错成一致的结果。还有一个容易被忽略的取舍是频率分辨率。0.1 GHz 步长在 60 GHz 那片宽平台上够用但在 22.235 GHz 的尖峰附近峰顶可能被跳过画出来的峰高偏低。正确做法是在已知谱线中心附近做局部加密或者在画图前先按谱线表生成一组“必采频点”并进频率向量。这个细节不会影响链路预算的结论但会影响你对“峰到底有多尖”的判断。从那以后我每次做这类谱线复现都会强制走一遍“三频点体检 两算法互校 量纲断言”的流程哪怕只是改一个系数也先跑一遍再画图。吃过一次亏某个高湿条件下 183 GHz 的峰值莫名比 60 GHz 低查了两天才发现是水汽谱线表少录了一行而那一行恰好是 183 GHz 的主线——曲线整体看起来“挺像”只有跟基准对比才暴露。这套流程连同谱线表、绘图脚本和校验函数我整理成了一份可直接运行的 matlab 代码放在下载包里拿到之后先跑check_points.m看三个体检频点是否在预期区间再去改参数探索不同温湿条件会省掉很多对着图猜原因的夜晚。希望帮到你。本文还有配套的精品资源点击获取