企业官网建设流程全解析

简介：这是一份面向通信工程、电磁场与电波传播方向研究者的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 km，0.157 是 10^6 除以地球半径的结果。这个公式把“大气折射”和“地球曲率”合并成一个量，等效于把球面大气拉平，让所有分析都变成在平地上的折射问题。标准大气下 dM/dh 约等于 0.118 M/m，此时射线恰好沿地表平行方向弯曲；一旦某段高度范围内 dM/dh 小于 0，说明电波向下弯得比地表还快，这就是波导层。

2.2 用双线性剖面构造表面波导与悬空波导

实测剖面来源很多，探空数据、再分析资料、海气耦合模式输出都能用，但做机理分析时最常用的是双线性近似：一段负梯度陷获层加一段正常梯度层。负梯度层是波导的核心，正常层用标准 0.118 M/m 就行。下面这组参数覆盖了三种最常见的波导形态。

注意这是做仿真捏合剖面用的经验范围，不是气候统计值。梯度绝对值越大，陷获能力越强；层越厚，能约束的波长范围越宽。悬空波导和表面波导的本质区别在于陷获层悬浮在空中，射线在它的上边界和下边界都会被弯回来，而表面波导的下边界直接是地面。

在 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(a+h)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°，对应 h0=20、th0=0.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

直线飞行段用 dMdh=0 是个工程折中：剖面顶以上没有数据，按折射率梯度消失处理，射线走直线，不会影响波导内的轨迹形态。地面截断通过 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。如果省略 direction，Ray 落地穿到负高度后，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

循环内每次都调用 ode45，30 个仰角、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° 即可。Htop+50 的余量是为了容忍波导顶部的轻微溢出振荡。注意 is_trapped 用的是 300 公里最大距离，这意味着只关心射线在长距离内不逃逸的情况。

5.2 与抛物方程法对比时边界条件怎么对齐

射线描迹是几何光学近似，频率偏低或距离偏远时，衍射和干涉效应会让它高估聚焦区能量。抛物方程法是目前公认更完整的标量波解法，matlab 里分步傅里叶实现的公开代码很多。与 PE 对标时三件事必须一致：剖面换算一致，PE 用折射指数实部，射线这边用 M 剖面，两边通过 M=N+0.157h 互相转换；初始场一致，PE 用窗函数加方向图，射线这边按方向图采样出的仰角逐条计算；边界条件一致，PE 顶部要有吸收层，底部在海上用阻抗边界，射线这边只有地面截断。两者数值不可能完全重合，拿“100 公里处损耗主峰的位置”来比对，偏差小于一个波束宽度就认为自洽。

5.3 检查 M 剖面垂直分辨率的快速脚本

最后提供一个开工前必跑的检查：把剖面网格加密一倍，重新求临界陷获角和 100 公里处落点，两次结果差异大就说明原始剖面分辨率不足。用 matlab 写就是重新采样后调用同一个二分函数，临界角差超过 0.01° 时给出告警。这个检查不依赖任何解析解，却能筛掉一批实测剖面——比如原始数据只有 50 米间隔，插值到 0.5 米后临界角看似没问题，但落点能差好几公里。现在用 codex 这类工具改 matlab 脚本已经很快，它能像操作 python 任务一样把循环、事件函数的样板补好，但剖面分辨率够不够，仍然只有这段检查脚本能回答。

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