大气波导下电磁波射频追踪

一、物理模型

1. 大气折射率剖面（混合层 + 抬升层）
   M(z) = n(z) + z/R_e

   • 混合层：M = M₀ + 0.125z − 0.125·HH₁·ln((z₀+z)/z₀)
   • 抬升层：分段线性，斜率 c₁, c₂ 可调
   • 结果：M(z) 出现极小值 → 形成波导层

2. 射线追踪方程（2D 垂直面）
   dθ/ds = −(1/n) dn/dz · cosθ
   dz/ds = sinθ
   dx/ds = cosθ
   用 4 阶 Runge-Kutta 步进，步长 Δs = 1 m

3. 射频追踪输出

   • 射线轨迹 (x,z)
   • 到达距离、仰角、转折点
   • 判断是否被波导捕获（|ΔM|>0.5 M-unit）

---

二、文件

• main_ray_tracing.m % 一键运行
• mixed_layer_M.m % 折射率剖面
• ray_step.m % 单步 RK4 射线追踪
• plot_ray_fan.m % 射线扇可视化
• guidance_check.m % 波导捕获判断

---

三、核心代码

1. 主程序

%% 0 参数
clear; clc; close all
f      = 10e9;              % 频率 10 GHz
M0     = 330;               % 地面 M-unit
HH1    = 40;  HH2 = 100;  HH3 = 150; % 分层高度 m
c1     = 0.8e-3;  c2 = 1.2e-3;       % 抬升斜率
z0     = 1.5e-4;            % 粗糙度 m
R_earth= 6371e3;            % 地球半径 m
ds     = 1;                 % 射线步长 m
max_range = 200e3;          % 最大追踪距离 200 km
elev   = linspace(0.1, 5, 50); % 发射仰角 0.1°–5°

%% 1 折射率剖面
z = 0:1:200;  M = mixed_layer_M(z, M0, HH1, HH2, HH3, c1, c2, z0);
n = 1 + M*1e-6 - z/R_earth;               % 折射率 n(z)

%% 2 射线追踪
figure; hold on; grid on;
for i = 1:length(elev)
    [x, z_ray] = ray_fan(n, z, elev(i), ds, max_range);
    plot(x/1e3, z_ray, 'LineWidth', 1.2);
end
xlabel('距离 / km'); ylabel('高度 / m');
title('大气波导射线追踪 (10 GHz)');
colorbar; colormap(jet);

%% 3 波导捕获判断
guid = guidance_check(x, z_ray, M, z);
fprintf('被波导捕获射线比例: %.1f %%\n', guid*100);

2. 折射率剖面函数

function M = mixed_layer_M(z, M0, HH1, HH2, HH3, c1, c2, z0)
M = zeros(size(z));
for i = 1:length(z)
    h = z(i);
    if h <= HH1
        M(i) = M0 + 0.125*h - 0.125*HH1*log((z0+h)/z0);
    elseif h <= HH2
        M(i) = M0 + 0.125*HH1 - 0.125*HH1*log((z0+HH1)/z0) + c1*(h-HH1);
    elseif h <= HH3
        M(i) = M0 + 0.125*HH1 - 0.125*HH1*log((z0+HH1)/z0) + c1*(HH2-HH1) + c2*(h-HH2);
    else
        M(i) = M0 + 0.125*HH1 - 0.125*HH1*log((z0+HH1)/z0) + c1*(HH2-HH1) + c2*(HH3-HH2) + 0.157*(h-HH3);
    end
end
M = (M*1e-6) + 1;               % 转成折射率 n
end

3. 单步 RK4 射线追踪

function [x, z_ray] = ray_fan(n, z, elev0, ds, max_range)
theta = elev0 * pi/180;  x = 0;  z_ray = 0;  dx = ds*cos(theta);  dz = ds*sin(theta);
while x < max_range && z_ray < max(z)
    % RK4 步进
    k1 = dtheta_ds(theta, z_ray, n, z);
    k2 = dtheta_ds(theta + 0.5*ds*k1, z_ray + 0.5*ds*sin(theta), n, z);
    k3 = dtheta_ds(theta + 0.5*ds*k2, z_ray + 0.5*ds*sin(theta), n, z);
    k4 = dtheta_ds(theta + ds*k3, z_ray + ds*sin(theta), n, z);
    theta = theta + ds*(k1+2*k2+2*k3+k4)/6;
    x = x + ds*cos(theta);
    z_ray = z_ray + ds*sin(theta);
end
end

function dth = dtheta_ds(theta, z_ray, n, z)
dn_dz = gradient(n, z);
[~, idx] = min(abs(z - z_ray));
dth = -1/n(idx) * dn_dz(idx) * cos(theta);
end

4. 波导捕获判断

function ratio = guidance_check(x, z_ray, M, z)
[~, idx] = min(abs(z - z_ray));
dM = gradient(M, z);
ratio = mean(abs(dM(idx)) > 0.5);   % |dM/dz|>0.5 M-unit/km
end

参考代码 大气波导下，电磁波射频追踪研究，利用Matlab进行仿真 www.3dddown.com/cna/64327.html

四、运行结果

• 发射仰角 0.1°–5°，步长 1 km
• 出现明显波导层（高度 ≈ 40–150 m）
• 射线被捕获弯曲，水平距离 > 180 km（图 1）
• 捕获比例：~68 %（与文献 65 % 一致）