简介这是一份面向通信工程、电磁场与电波传播方向研究者的MATLAB仿真资源聚焦大气波导环境下的射线描迹。大气波导使电波在对流层逆温层中反复折射形成超视距传播是高频通信与雷达探测中需要特别分析的现象这份代码能帮助初学者和工程师直观模拟多径射线路径理解直达波、地面反射波和大气层反射波的成因。包内共3个文件以MATLAB脚本.m为主配合2张结果图.jpg方便对照验证压缩包仅41KB轻量易用适合快速跑通模型。已有220人浏览学习可用于无线通信课程设计、科研预研或日常教学演示。运行后不仅能掌握射线追踪的基本流程还能观察不同大气结构对电波路径的影响为通信系统覆盖预测与干扰分析提供实用的仿真基础。1. 大气波导射线描迹在修正折射率剖面里追一条电磁波的路径雷达屏幕上多出几百公里的幽灵回波、微波链路越过地平线还能稳定收发大气波导常常是幕后推手。对流层中温度湿度分布异常时电波弯曲程度超过地球曲率能量被限制在几百米厚的薄层里形成超视距传播。要把这种现象算清楚最直接的工具就是射线描迹用几何光学近似把传播问题变成沿路径积分高度和仰角的初值问题。用 matlab 实现时只需要三样东西一条修正折射率剖面、一个射线微分方程、一段绘图脚本。标题里的 daqibodao.rar 就是这类问题常见的代码包。解压后多数情况下能整理成剖面构造、射线求解、可视化三个模块。文章按这个顺序展开先把修正折射率 M 和波导判据讲透再给出可直接复现的 matlab 脚本最后做仰角扫描和与抛物方程法的对标。适合做雷达覆盖评估、微波链路设计的工程师也适合刚接触电波传播仿真、想用 matlab 把射线描迹跑通的研究生。2. 修正折射率剖面与波导判据把 dM/dh 这条主线立起来2.1 为什么要用修正折射率 M 而不是直接用折射指数 N真实大气折射指数 n 非常接近 1垂直方向变化只有百万分之一量级直接画 n 的剖面基本是一条平线看不出任何结构。工程上先放大成折射率 N (n-1)×10^6单位是 N-unit典型地面值在 300 到 350 之间。但 N 剖面仍然不够直观因为地球是个球面即便大气完全均匀一条沿直线传播的射线相对地表也会“越来越高”——地表本身在往下弯。修正折射率 M 就是为解决这个问题提出的M N (h/a)×10^6 ≈ N 0.157×h。其中 h 是海拔米a 取地球半径 6371 km0.157 是 10^6 除以地球半径的结果。这个公式把“大气折射”和“地球曲率”合并成一个量等效于把球面大气拉平让所有分析都变成在平地上的折射问题。标准大气下 dM/dh 约等于 0.118 M/m此时射线恰好沿地表平行方向弯曲一旦某段高度范围内 dM/dh 小于 0说明电波向下弯得比地表还快这就是波导层。2.2 用双线性剖面构造表面波导与悬空波导实测剖面来源很多探空数据、再分析资料、海气耦合模式输出都能用但做机理分析时最常用的是双线性近似一段负梯度陷获层加一段正常梯度层。负梯度层是波导的核心正常层用标准 0.118 M/m 就行。下面这组参数覆盖了三种最常见的波导形态。波导类型陷获层厚度 (m)层内 dM/dh (M/m)层外 dM/dh (M/m)常见场景蒸发波导840-0.5-0.250.118海面、热带海域雷达接地表面波导50300-0.3-0.10.118辐射逆温、陆地早晨悬空波导200800-0.15-0.050.118上下两侧副热带高压、锋面附近注意这是做仿真捏合剖面用的经验范围不是气候统计值。梯度绝对值越大陷获能力越强层越厚能约束的波长范围越宽。悬空波导和表面波导的本质区别在于陷获层悬浮在空中射线在它的上边界和下边界都会被弯回来而表面波导的下边界直接是地面。在 matlab 里把这样一个剖面写成函数并不复杂关键是分段点要连续。function M m_profile_bi(hvec, h_trap, g1, g2, M0) % 双线性修正折射率剖面 % hvec 高度向量 (m) % h_trap 陷获层顶高 (m) % g1 陷获层内梯度 (M/m)应为负数 % g2 陷获层上方梯度 (M/m)通常取 0.118 % M0 地面处 M 值典型 330 左右 M zeros(size(hvec)); for k 1:numel(hvec) h hvec(k); if h h_trap M(k) M0 g1 * h; else M(k) M0 g1 * h_trap g2 * (h - h_trap); end end end调用方式很直接Mvec m_profile_bi((0:0.5:2000), 100, -0.3, 0.118, 330);这段代码在 h_trap 处保证 M 值连续因为第二段里补偿了 g1×h_trap 这一项。分段线性剖面最常见的错误就是忽略这个补偿项导致连接点出现一个台阶后续求梯度时会产生一个虚假的尖峰轨迹会被这个假尖峰弹一下。2.3 用 dM/dh 快速判断射线是穿透还是被陷获剖面构造完第一件事不是算射线而是把梯度画出来看。在 matlab 里计算 M 的垂直梯度用 gradient 函数就行。dh_step 0.5; dMdh gradient(Mvec, dh_step); figure(Color, w); plot(Mvec, Hvec, LineWidth, 1.5); hold on; yyaxis right; plot(dMdh, Hvec, --, LineWidth, 1); grid on; xlabel(M (M-unit) / dMdh (M-unit/m)); ylabel(高度 (m)); legend(M 剖面, dM/dh, Location, best); title(波导剖面梯度检查);运行后会看到 M 曲线在 100 米以下向下弯dM/dh 在负区间出现一个明显的负值平台。判据只有一条dM/dh 0 的连续区间就是潜在陷获层。但这里有个容易误判的点负梯度绝对值不够大、或者层厚不够时即使 dM/dh 为负也陷不住特定波长的电波。定量判断要交给第 5 章的临界仰角扫描梯度检查只是先把剖面中的结构错误筛掉比如拐点断裂、负层厚度只有两三个网格点这类一眼能看出的毛病。3. matlab 射线描迹实现从射线方程到可复现的 ode45 脚本3.1 球面分层大气下射线方程的简化形式射线描迹的理论起点是球面分层大气中的 Bouguer 不变量n(ah)cosθ 沿射线路径守恒其中 a 是地球半径h 是高度θ 是射线与当地水平面的夹角。对路径弧长 s 求导再用 n 与 M 的换算关系做近似最终会得到一个非常干净的常微分方程组dh/ds sinθdx/ds cosθdθ/ds 1e-6 × (dM/dh) × cosθ推导过程中用到了两个近似n 接近于 1以及 h 远小于 a。这两个条件在对流层低层完全成立。结果就是地球曲率项和折射梯度项在 M 坐标下互相抵消方程里只剩 dM/dh 这一个剖面量。对写代码的人来说这是最舒服的形式不用再单独引入地球半径、不用区分“几何高度”和“折射高度”输入剖面算梯度剩下的交给积分器。这个方程组的状态量是高度 h 和仰角 θ自变量是弧长 s。初始条件给出发射点高度和初始仰角比如海面雷达架高 20 米、波束仰角 0.05°对应 h020、th00.05×pi/180。积分到预设最大距离或者地面事件触发时停止。3.2 最小可跑脚本在 matlab 中定义微分方程并求解解这个方程组用 ode45 就够了光照强度不高自由度也只有两个。关键是在 matlab 中定义微分方程的右侧函数把剖面梯度按当前高度插值出来。function dyds ray_rhs(s, y, Hvec, Mvec, dMdh_vec) % 射线描迹微分方程状态 y [h; theta] % s 为弧长本函数不使用但 ode45 要求位置一致 h y(1); th y(2); if h Hvec(1) || h Hvec(end) dMdh 0; % 超出剖面范围按自由空间处理 else dMdh interp1(Hvec, dMdh_vec, h, linear); end dyds [sin(th); 1e-6 * dMdh * cos(th)]; end直线飞行段用 dMdh0 是个工程折中剖面顶以上没有数据按折射率梯度消失处理射线走直线不会影响波导内的轨迹形态。地面截断通过 odeset 的事件函数实现。function [value, isterminal, direction] ground_evt(~, y) value y(1); % 高度降到 0 时触发 isterminal 1; direction -1; % 只捕获从正到负的穿越 end主脚本把剖面、初始条件和积分器串起来。Hmax 2000; dh_step 0.5; Hvec (0:dh_step:Hmax); Mvec m_profile_bi(Hvec, 100, -0.3, 0.118, 330); dMdh_vec gradient(Mvec, dh_step); h0 20; th0 0.05 * pi/180; opt odeset(Events, ground_evt, RelTol, 1e-6); [s, y] ode45((s,y) ray_rhs(s, y, Hvec, Mvec, dMdh_vec), ... [0 300e3], [h0; th0], opt); x cumtrapz(s, cos(y(:,2))); % 弧长积分出水平距离 figure(Color, w); plot(x/1000, y(:,1), b, LineWidth, 1.2); xlabel(水平距离 (km)); ylabel(高度 (m)); grid on;代码逻辑说明y 矩阵第一列是高度第二列是弧度制仰角。水平距离不能用 s 直接代替因为射线有仰角水平投影是 s×cosθ 的积分所以用 cumtrapz 累积。RelTol 取 1e-6 与 1e-9 相比100 公里处的高度差通常小于 1 米不需要更紧。3.3 地面截断、剖面插值与数值步长的三个细节第一个细节是事件函数的 direction。如果省略 directionRay 落地穿到负高度后ode45 会在射线重新从地下穿回时再次触发事件一条轨迹可能被截成两段。明确写 direction-1 后事件只在高度从正变负的瞬间触发一次。第二个细节是剖面拐点的网格对齐。线性插值在拐点处会把尖角抹成平滑过渡带过渡带宽等于一个网格间距。陷获层只有 50 米厚时0.5 米网格的过渡带占比 1%影响不大但网格放到 5 米过渡带就占了 10%临界仰角会偏移明显。让 h_trap 落在网格节点上是零成本的修正。第三个细节是网格加密验证。把 dh_step 减半重跑同一条射线对比 100 公里处的高度差异超过几十米就说明原网格偏粗。网格问题在表面波导里尤其容易被忽视因为负梯度层本身薄插值误差会被 dM/dh 的假抖动放大。4. 仰角扫描与射线簇可视化把单条轨迹扩展成覆盖图4.1 批量算一组初始仰角并保存轨迹数据单条射线只能回答“这一个方向会不会进波导”。实际雷达波束有一定垂直张角天线方向图主瓣覆盖 -0.2° 到 0.8° 是常有的事。把所有仰角的轨迹都算一遍才能回答波导把能量送到了哪些距离、哪些区域被打出盲区。把主脚本包进 for 循环用结构数组存结果即可。th0_list (-0.2:0.02:0.8) * pi/180; traj struct([]); for k 1:numel(th0_list) [sk, yk] ode45((s,y) ray_rhs(s,y,Hvec,Mvec,dMdh_vec), ... [0 300e3], [h0; th0_list(k)], opt); traj(k).x cumtrapz(sk, cos(yk(:,2))); traj(k).h yk(:,1); traj(k).th0 th0_list(k) * 180/pi; end循环内每次都调用 ode4530 个仰角、300 公里传播距离在 0.5 米网格下通常几秒内完成不需要用 parfor 提前优化。之所以把轨迹存下来而不是边算边画是因为后面找交点、统计射线密度都要反复访问这些数据。th0 同时存弧度和度数是为了画图时图例直接显示度数。4.2 把射线簇和波导层边界画到同一张图上射线簇单独画出来只是几十条曲线没有高度参考看不出波导层的约束效果。常见做法是上下双子图上图画射线簇和陷获层顶下图画 M 剖面并标出负梯度区间。figure(Color, w, Position, [100 100 760 680]); subplot(2,1,1); hold on; for k 1:numel(traj) plot(traj(k).x/1000, traj(k).h, LineWidth, 0.7); end plot([0 300], [100 100], r--, LineWidth, 1.5); plot([0 300], [0 0], k, LineWidth, 2); xlabel(水平距离 (km)); ylabel(高度 (m)); ylim([0 1200]); grid on; subplot(2,1,2); hold on; plot(Mvec, Hvec, LineWidth, 1.5); neg_idx dMdh_vec 0; area(Hvec(neg_idx), Mvec(neg_idx), FaceAlpha, 0.2); xlabel(M (M-unit)); ylabel(高度 (m)); grid on;从这张图能直接读出三件事负仰角射线快速砸向地面接近 0° 的射线被陷获层顶压住在 100 米以下来回弯曲初始仰角超过某个值的射线直接穿出波导层高度一路抬升。area 函数把 dM/dh 的负区间填成半透明色块和上图的红色虚线完全对应剖面设置错误在两级对照下很容易发现。4.3 从射线交点识别聚焦区与通信盲区几何光学里能量沿射线管流动射线汇聚的地方能量密度高射线稀疏的地方就是覆盖盲区。虽然射线描迹不含衍射和干涉信息但用射线密度判断聚焦区位置和边界一直是工程上最常用的快速方法。密度统计可以在 matlab 里用距离门实现每隔若干公里统计该断面附近 20 米高度窗内的射线数量。xq 20:2:300; hit zeros(1, numel(xq)); for k 1:numel(traj) xk traj(k).x / 1000; for i 1:numel(xq) xc xq(i); if xc xk(1) || xc xk(end) continue; end hk interp1(xk, traj(k).h, xc); hit(i) hit(i) double(abs(hk - traj(k).h(end)) 20); end end这段代码为了可读性用了双层循环射线数在 100 条以内时没有必要向量化。把 hit 画成阶梯图后峰值位置通常对应“越障超视距”回波最明显的距离段。值得强调的是射线密度只是能量分布的粗略代理天线方向图加权和距离扩散损耗都没有计入它回答的是“哪一段大概率有信号”而不是“信号有多少 dBm”。5. 进阶陷获角数值判据与抛物方程法交叉验证5.1 用二分法找临界陷获仰角给定剖面后最大能陷获的初始仰角称为临界陷获角。解析计算需要解超越方程数值上二分法更省事。先定义一个逻辑函数判断某仰角下射线最大高度是否超过陷获层顶。function flag is_trapped(th0, h0, Hvec, Mvec, dMdh_vec, Htop) [~, y] ode45((s,y) ray_rhs(s,y,Hvec,Mvec,dMdh_vec), ... [0 300e3], [h0; th0]); flag max(y(:,1)) Htop 50; end从 0° 开始以 0.1° 步长向上探测找到第一个穿透仰角后在它和前一个仰角之间二分收敛到 0.001° 即可。Htop50 的余量是为了容忍波导顶部的轻微溢出振荡。注意 is_trapped 用的是 300 公里最大距离这意味着只关心射线在长距离内不逃逸的情况。5.2 与抛物方程法对比时边界条件怎么对齐射线描迹是几何光学近似频率偏低或距离偏远时衍射和干涉效应会让它高估聚焦区能量。抛物方程法是目前公认更完整的标量波解法matlab 里分步傅里叶实现的公开代码很多。与 PE 对标时三件事必须一致剖面换算一致PE 用折射指数实部射线这边用 M 剖面两边通过 MN0.157h 互相转换初始场一致PE 用窗函数加方向图射线这边按方向图采样出的仰角逐条计算边界条件一致PE 顶部要有吸收层底部在海上用阻抗边界射线这边只有地面截断。两者数值不可能完全重合拿“100 公里处损耗主峰的位置”来比对偏差小于一个波束宽度就认为自洽。5.3 检查 M 剖面垂直分辨率的快速脚本最后提供一个开工前必跑的检查把剖面网格加密一倍重新求临界陷获角和 100 公里处落点两次结果差异大就说明原始剖面分辨率不足。用 matlab 写就是重新采样后调用同一个二分函数临界角差超过 0.01° 时给出告警。这个检查不依赖任何解析解却能筛掉一批实测剖面——比如原始数据只有 50 米间隔插值到 0.5 米后临界角看似没问题但落点能差好几公里。现在用 codex 这类工具改 matlab 脚本已经很快它能像操作 python 任务一样把循环、事件函数的样板补好但剖面分辨率够不够仍然只有这段检查脚本能回答。本文还有配套的精品资源点击获取 
长河君 · 2026/9/16 1:18:02 · 阅读 1.2k 大气波导射线描迹与MATLAB仿真：从修正折射率剖面到电波路径
简介这是一份面向通信工程、电磁场与电波传播方向研究者的MATLAB仿真资源聚焦大气波导环境下的射线描迹。大气波导使电波在对流层逆温层中反复折射形成超视距传播是高频通信与雷达探测中需要特别分析的现象这份代码能帮助初学者和工程师直观模拟多径射线路径理解直达波、地面反射波和大气层反射波的成因。包内共3个文件以MATLAB脚本.m为主配合2张结果图.jpg方便对照验证压缩包仅41KB轻量易用适合快速跑通模型。已有220人浏览学习可用于无线通信课程设计、科研预研或日常教学演示。运行后不仅能掌握射线追踪的基本流程还能观察不同大气结构对电波路径的影响为通信系统覆盖预测与干扰分析提供实用的仿真基础。1. 大气波导射线描迹在修正折射率剖面里追一条电磁波的路径雷达屏幕上多出几百公里的幽灵回波、微波链路越过地平线还能稳定收发大气波导常常是幕后推手。对流层中温度湿度分布异常时电波弯曲程度超过地球曲率能量被限制在几百米厚的薄层里形成超视距传播。要把这种现象算清楚最直接的工具就是射线描迹用几何光学近似把传播问题变成沿路径积分高度和仰角的初值问题。用 matlab 实现时只需要三样东西一条修正折射率剖面、一个射线微分方程、一段绘图脚本。标题里的 daqibodao.rar 就是这类问题常见的代码包。解压后多数情况下能整理成剖面构造、射线求解、可视化三个模块。文章按这个顺序展开先把修正折射率 M 和波导判据讲透再给出可直接复现的 matlab 脚本最后做仰角扫描和与抛物方程法的对标。适合做雷达覆盖评估、微波链路设计的工程师也适合刚接触电波传播仿真、想用 matlab 把射线描迹跑通的研究生。2. 修正折射率剖面与波导判据把 dM/dh 这条主线立起来2.1 为什么要用修正折射率 M 而不是直接用折射指数 N真实大气折射指数 n 非常接近 1垂直方向变化只有百万分之一量级直接画 n 的剖面基本是一条平线看不出任何结构。工程上先放大成折射率 N (n-1)×10^6单位是 N-unit典型地面值在 300 到 350 之间。但 N 剖面仍然不够直观因为地球是个球面即便大气完全均匀一条沿直线传播的射线相对地表也会“越来越高”——地表本身在往下弯。修正折射率 M 就是为解决这个问题提出的M N (h/a)×10^6 ≈ N 0.157×h。其中 h 是海拔米a 取地球半径 6371 km0.157 是 10^6 除以地球半径的结果。这个公式把“大气折射”和“地球曲率”合并成一个量等效于把球面大气拉平让所有分析都变成在平地上的折射问题。标准大气下 dM/dh 约等于 0.118 M/m此时射线恰好沿地表平行方向弯曲一旦某段高度范围内 dM/dh 小于 0说明电波向下弯得比地表还快这就是波导层。2.2 用双线性剖面构造表面波导与悬空波导实测剖面来源很多探空数据、再分析资料、海气耦合模式输出都能用但做机理分析时最常用的是双线性近似一段负梯度陷获层加一段正常梯度层。负梯度层是波导的核心正常层用标准 0.118 M/m 就行。下面这组参数覆盖了三种最常见的波导形态。波导类型陷获层厚度 (m)层内 dM/dh (M/m)层外 dM/dh (M/m)常见场景蒸发波导840-0.5-0.250.118海面、热带海域雷达接地表面波导50300-0.3-0.10.118辐射逆温、陆地早晨悬空波导200800-0.15-0.050.118上下两侧副热带高压、锋面附近注意这是做仿真捏合剖面用的经验范围不是气候统计值。梯度绝对值越大陷获能力越强层越厚能约束的波长范围越宽。悬空波导和表面波导的本质区别在于陷获层悬浮在空中射线在它的上边界和下边界都会被弯回来而表面波导的下边界直接是地面。在 matlab 里把这样一个剖面写成函数并不复杂关键是分段点要连续。function M m_profile_bi(hvec, h_trap, g1, g2, M0) % 双线性修正折射率剖面 % hvec 高度向量 (m) % h_trap 陷获层顶高 (m) % g1 陷获层内梯度 (M/m)应为负数 % g2 陷获层上方梯度 (M/m)通常取 0.118 % M0 地面处 M 值典型 330 左右 M zeros(size(hvec)); for k 1:numel(hvec) h hvec(k); if h h_trap M(k) M0 g1 * h; else M(k) M0 g1 * h_trap g2 * (h - h_trap); end end end调用方式很直接Mvec m_profile_bi((0:0.5:2000), 100, -0.3, 0.118, 330);这段代码在 h_trap 处保证 M 值连续因为第二段里补偿了 g1×h_trap 这一项。分段线性剖面最常见的错误就是忽略这个补偿项导致连接点出现一个台阶后续求梯度时会产生一个虚假的尖峰轨迹会被这个假尖峰弹一下。2.3 用 dM/dh 快速判断射线是穿透还是被陷获剖面构造完第一件事不是算射线而是把梯度画出来看。在 matlab 里计算 M 的垂直梯度用 gradient 函数就行。dh_step 0.5; dMdh gradient(Mvec, dh_step); figure(Color, w); plot(Mvec, Hvec, LineWidth, 1.5); hold on; yyaxis right; plot(dMdh, Hvec, --, LineWidth, 1); grid on; xlabel(M (M-unit) / dMdh (M-unit/m)); ylabel(高度 (m)); legend(M 剖面, dM/dh, Location, best); title(波导剖面梯度检查);运行后会看到 M 曲线在 100 米以下向下弯dM/dh 在负区间出现一个明显的负值平台。判据只有一条dM/dh 0 的连续区间就是潜在陷获层。但这里有个容易误判的点负梯度绝对值不够大、或者层厚不够时即使 dM/dh 为负也陷不住特定波长的电波。定量判断要交给第 5 章的临界仰角扫描梯度检查只是先把剖面中的结构错误筛掉比如拐点断裂、负层厚度只有两三个网格点这类一眼能看出的毛病。3. matlab 射线描迹实现从射线方程到可复现的 ode45 脚本3.1 球面分层大气下射线方程的简化形式射线描迹的理论起点是球面分层大气中的 Bouguer 不变量n(ah)cosθ 沿射线路径守恒其中 a 是地球半径h 是高度θ 是射线与当地水平面的夹角。对路径弧长 s 求导再用 n 与 M 的换算关系做近似最终会得到一个非常干净的常微分方程组dh/ds sinθdx/ds cosθdθ/ds 1e-6 × (dM/dh) × cosθ推导过程中用到了两个近似n 接近于 1以及 h 远小于 a。这两个条件在对流层低层完全成立。结果就是地球曲率项和折射梯度项在 M 坐标下互相抵消方程里只剩 dM/dh 这一个剖面量。对写代码的人来说这是最舒服的形式不用再单独引入地球半径、不用区分“几何高度”和“折射高度”输入剖面算梯度剩下的交给积分器。这个方程组的状态量是高度 h 和仰角 θ自变量是弧长 s。初始条件给出发射点高度和初始仰角比如海面雷达架高 20 米、波束仰角 0.05°对应 h020、th00.05×pi/180。积分到预设最大距离或者地面事件触发时停止。3.2 最小可跑脚本在 matlab 中定义微分方程并求解解这个方程组用 ode45 就够了光照强度不高自由度也只有两个。关键是在 matlab 中定义微分方程的右侧函数把剖面梯度按当前高度插值出来。function dyds ray_rhs(s, y, Hvec, Mvec, dMdh_vec) % 射线描迹微分方程状态 y [h; theta] % s 为弧长本函数不使用但 ode45 要求位置一致 h y(1); th y(2); if h Hvec(1) || h Hvec(end) dMdh 0; % 超出剖面范围按自由空间处理 else dMdh interp1(Hvec, dMdh_vec, h, linear); end dyds [sin(th); 1e-6 * dMdh * cos(th)]; end直线飞行段用 dMdh0 是个工程折中剖面顶以上没有数据按折射率梯度消失处理射线走直线不会影响波导内的轨迹形态。地面截断通过 odeset 的事件函数实现。function [value, isterminal, direction] ground_evt(~, y) value y(1); % 高度降到 0 时触发 isterminal 1; direction -1; % 只捕获从正到负的穿越 end主脚本把剖面、初始条件和积分器串起来。Hmax 2000; dh_step 0.5; Hvec (0:dh_step:Hmax); Mvec m_profile_bi(Hvec, 100, -0.3, 0.118, 330); dMdh_vec gradient(Mvec, dh_step); h0 20; th0 0.05 * pi/180; opt odeset(Events, ground_evt, RelTol, 1e-6); [s, y] ode45((s,y) ray_rhs(s, y, Hvec, Mvec, dMdh_vec), ... [0 300e3], [h0; th0], opt); x cumtrapz(s, cos(y(:,2))); % 弧长积分出水平距离 figure(Color, w); plot(x/1000, y(:,1), b, LineWidth, 1.2); xlabel(水平距离 (km)); ylabel(高度 (m)); grid on;代码逻辑说明y 矩阵第一列是高度第二列是弧度制仰角。水平距离不能用 s 直接代替因为射线有仰角水平投影是 s×cosθ 的积分所以用 cumtrapz 累积。RelTol 取 1e-6 与 1e-9 相比100 公里处的高度差通常小于 1 米不需要更紧。3.3 地面截断、剖面插值与数值步长的三个细节第一个细节是事件函数的 direction。如果省略 directionRay 落地穿到负高度后ode45 会在射线重新从地下穿回时再次触发事件一条轨迹可能被截成两段。明确写 direction-1 后事件只在高度从正变负的瞬间触发一次。第二个细节是剖面拐点的网格对齐。线性插值在拐点处会把尖角抹成平滑过渡带过渡带宽等于一个网格间距。陷获层只有 50 米厚时0.5 米网格的过渡带占比 1%影响不大但网格放到 5 米过渡带就占了 10%临界仰角会偏移明显。让 h_trap 落在网格节点上是零成本的修正。第三个细节是网格加密验证。把 dh_step 减半重跑同一条射线对比 100 公里处的高度差异超过几十米就说明原网格偏粗。网格问题在表面波导里尤其容易被忽视因为负梯度层本身薄插值误差会被 dM/dh 的假抖动放大。4. 仰角扫描与射线簇可视化把单条轨迹扩展成覆盖图4.1 批量算一组初始仰角并保存轨迹数据单条射线只能回答“这一个方向会不会进波导”。实际雷达波束有一定垂直张角天线方向图主瓣覆盖 -0.2° 到 0.8° 是常有的事。把所有仰角的轨迹都算一遍才能回答波导把能量送到了哪些距离、哪些区域被打出盲区。把主脚本包进 for 循环用结构数组存结果即可。th0_list (-0.2:0.02:0.8) * pi/180; traj struct([]); for k 1:numel(th0_list) [sk, yk] ode45((s,y) ray_rhs(s,y,Hvec,Mvec,dMdh_vec), ... [0 300e3], [h0; th0_list(k)], opt); traj(k).x cumtrapz(sk, cos(yk(:,2))); traj(k).h yk(:,1); traj(k).th0 th0_list(k) * 180/pi; end循环内每次都调用 ode4530 个仰角、300 公里传播距离在 0.5 米网格下通常几秒内完成不需要用 parfor 提前优化。之所以把轨迹存下来而不是边算边画是因为后面找交点、统计射线密度都要反复访问这些数据。th0 同时存弧度和度数是为了画图时图例直接显示度数。4.2 把射线簇和波导层边界画到同一张图上射线簇单独画出来只是几十条曲线没有高度参考看不出波导层的约束效果。常见做法是上下双子图上图画射线簇和陷获层顶下图画 M 剖面并标出负梯度区间。figure(Color, w, Position, [100 100 760 680]); subplot(2,1,1); hold on; for k 1:numel(traj) plot(traj(k).x/1000, traj(k).h, LineWidth, 0.7); end plot([0 300], [100 100], r--, LineWidth, 1.5); plot([0 300], [0 0], k, LineWidth, 2); xlabel(水平距离 (km)); ylabel(高度 (m)); ylim([0 1200]); grid on; subplot(2,1,2); hold on; plot(Mvec, Hvec, LineWidth, 1.5); neg_idx dMdh_vec 0; area(Hvec(neg_idx), Mvec(neg_idx), FaceAlpha, 0.2); xlabel(M (M-unit)); ylabel(高度 (m)); grid on;从这张图能直接读出三件事负仰角射线快速砸向地面接近 0° 的射线被陷获层顶压住在 100 米以下来回弯曲初始仰角超过某个值的射线直接穿出波导层高度一路抬升。area 函数把 dM/dh 的负区间填成半透明色块和上图的红色虚线完全对应剖面设置错误在两级对照下很容易发现。4.3 从射线交点识别聚焦区与通信盲区几何光学里能量沿射线管流动射线汇聚的地方能量密度高射线稀疏的地方就是覆盖盲区。虽然射线描迹不含衍射和干涉信息但用射线密度判断聚焦区位置和边界一直是工程上最常用的快速方法。密度统计可以在 matlab 里用距离门实现每隔若干公里统计该断面附近 20 米高度窗内的射线数量。xq 20:2:300; hit zeros(1, numel(xq)); for k 1:numel(traj) xk traj(k).x / 1000; for i 1:numel(xq) xc xq(i); if xc xk(1) || xc xk(end) continue; end hk interp1(xk, traj(k).h, xc); hit(i) hit(i) double(abs(hk - traj(k).h(end)) 20); end end这段代码为了可读性用了双层循环射线数在 100 条以内时没有必要向量化。把 hit 画成阶梯图后峰值位置通常对应“越障超视距”回波最明显的距离段。值得强调的是射线密度只是能量分布的粗略代理天线方向图加权和距离扩散损耗都没有计入它回答的是“哪一段大概率有信号”而不是“信号有多少 dBm”。5. 进阶陷获角数值判据与抛物方程法交叉验证5.1 用二分法找临界陷获仰角给定剖面后最大能陷获的初始仰角称为临界陷获角。解析计算需要解超越方程数值上二分法更省事。先定义一个逻辑函数判断某仰角下射线最大高度是否超过陷获层顶。function flag is_trapped(th0, h0, Hvec, Mvec, dMdh_vec, Htop) [~, y] ode45((s,y) ray_rhs(s,y,Hvec,Mvec,dMdh_vec), ... [0 300e3], [h0; th0]); flag max(y(:,1)) Htop 50; end从 0° 开始以 0.1° 步长向上探测找到第一个穿透仰角后在它和前一个仰角之间二分收敛到 0.001° 即可。Htop50 的余量是为了容忍波导顶部的轻微溢出振荡。注意 is_trapped 用的是 300 公里最大距离这意味着只关心射线在长距离内不逃逸的情况。5.2 与抛物方程法对比时边界条件怎么对齐射线描迹是几何光学近似频率偏低或距离偏远时衍射和干涉效应会让它高估聚焦区能量。抛物方程法是目前公认更完整的标量波解法matlab 里分步傅里叶实现的公开代码很多。与 PE 对标时三件事必须一致剖面换算一致PE 用折射指数实部射线这边用 M 剖面两边通过 MN0.157h 互相转换初始场一致PE 用窗函数加方向图射线这边按方向图采样出的仰角逐条计算边界条件一致PE 顶部要有吸收层底部在海上用阻抗边界射线这边只有地面截断。两者数值不可能完全重合拿“100 公里处损耗主峰的位置”来比对偏差小于一个波束宽度就认为自洽。5.3 检查 M 剖面垂直分辨率的快速脚本最后提供一个开工前必跑的检查把剖面网格加密一倍重新求临界陷获角和 100 公里处落点两次结果差异大就说明原始剖面分辨率不足。用 matlab 写就是重新采样后调用同一个二分函数临界角差超过 0.01° 时给出告警。这个检查不依赖任何解析解却能筛掉一批实测剖面——比如原始数据只有 50 米间隔插值到 0.5 米后临界角看似没问题但落点能差好几公里。现在用 codex 这类工具改 matlab 脚本已经很快它能像操作 python 任务一样把循环、事件函数的样板补好但剖面分辨率够不够仍然只有这段检查脚本能回答。本文还有配套的精品资源点击获取