简介针对大气波导环境下电波传播研究的专业教学演示文稿面向通信工程、电磁场与微波技术、大气科学等方向的学生、教师与科研人员可作为课堂讲义、自学笔记或课题起步的参考资料。内容共25页从大气折射的基本类型及存在条件入手逐步深入到修正折射率计算、陷获折射与大气波导形成条件再系统讲解抛物方程方法的推导过程、分步傅里叶变换实现、粗糙海面阻抗边界设置、初始场处理以及电波传播损耗的计算方式。同时演示文稿专门给出MATLAB电波传播仿真模块帮助读者掌握如何将理论模型转化为可计算的仿真流程并对传播损耗结果进行可视化分析。该演示文稿为单个pptx文件压缩包仅472KB结构紧凑、便于携带和二次编辑现有486人学习尤其适合开展超视距传播机理、雷达探测距离评估和通信链路设计等方向的研究人员参考能够快速建立从大气波导物理机制到数值建模与仿真的完整认知。1. 大气波导不是玄学为什么研究电波传播要先啃抛物方程这是一份可以直接当讲义用的专业资料主题是大气波导环境下的电波传播核心解法是抛物方程Parabolic Equation, PE法。干雷达链路、微波通信或电磁仿真的人多半遇到过这种场景海上目标在视距之外几十公里被雷达抓到或者通信链路在某个季节突然多出几十公里余量——第一个怀疑对象就是大气波导。可波导这个东西看不见摸不着用射线法解释又漏洞百出真正能落地算损耗的还是抛物方程。这套学习教案的价值在于不绕开数学也不停在公式层面把折射率剖面、PE推导、分步傅里叶变换实现和结果判读串成一条可复现的线。适合电磁场专业研究生、雷达链路设计工程师以及想给团队做内部培训的算法负责人。2. 大气波导环境参数修正折射率剖面与三类波导的判定2.1 为什么要看修正折射率M而不是直接看折射率n大气折射率n在数值上接近1直接看它很难看出波导结构。工程上习惯用折射率差N (n - 1)×10⁶单位是N单位但在考虑地球曲率后还要把高度项加进去得到修正折射率MM N (z / a)×10⁶其中a 6370 kmz是离地高度km。这样做的原因很朴素电磁波在球面地球上传播射线弯曲是大气折射和地球曲率共同作用的结果。把地球“扳直”之后判断波导只需要看M随高度变化的斜率dM/dz而不需要拆开看两条曲线。dM/dz是判定折射类型的核心指标。标准大气条件下低层大气N梯度约为-39 N单位/km折算后dM/dz约为0.118 M单位/m。当dM/dz低于这个值射线向下弯曲得更厉害进入超折射状态一旦dM/dz降到0以下电磁波就被限制在一层“管道”里反复折射这就是波导捕获。为了快速查判我习惯用下面这张表dM/dz 范围 (M单位/m)折射类型对传播的影响 0.118亚折射射线向上弯曲传播距离被压缩≈ 0.118标准折射正常视距传播0 ~ 0.118超折射低空路径弯曲雷达出现超视距现象 0波导捕获能量被限制在波导层内低损耗远距离传播注意这里的0.118是从标准大气推出来的参考值实际大气波动很大海面、沙漠、高原都不一样。拿到任何气象数据后第一步先画dM/dz剖面看有没有连续负梯度区段这是判断波导是否成立的唯一硬标准。2.2 三类典型波导蒸发波导、表面波导与抬升波导大气波导按形成高度可以分成三类工程上判读场景时先分清是哪种因为它们的剖面形状、影响频段和建模方式差别很大。蒸发波导是最常见的海上波导发生在海面以上几米到几十米由海面水汽蒸发造成湿度随高度骤降引起。它的特点是贴地且普遍存在热带海域出现概率极高对C波段以上频率影响显著。剖面形状近似为低层M随高度下降到波导顶之后恢复标准梯度。表面波导是接地或贴近地面的波导层顶高一般在几十米到几百米由夜间辐射冷却、逆温层等造成。它的特征是M剖面从地面开始就有超折射或负梯度上面再叠一个正梯度层。陆地大雾、夜间稳定大气时容易出这种波导对地面雷达和移动通信基站影响最大。抬升波导不接地波导层悬在几百米到上千米的高度由下沉逆温、锋面抬升等过程形成。它的M剖面在波导层内呈负梯度上下两侧都是标准梯度。这种波导对低空目标影响不大但会漏进高空目标容易让雷达在仰角上出现奇怪的回波。波导类型典型高度成因M剖面特征主要影响场景蒸发波导5 ~ 40 m海面水汽蒸发贴地负梯度快速恢复舰载雷达、海上微波链路表面波导50 ~ 500 m逆温、辐射冷却从地面开始负梯度地面雷达、移动通信抬升波导500 ~ 2000 m下沉逆温、锋面悬空负梯度层高空目标探测、干扰分析2.3 从实测气象数据构造M剖面建模的第一步是把温度、气压、湿度廓线转成M剖面。折射率N的经典计算公式是N 77.6 × P / T 3.73×10⁵ × e / T²其中P是气压hPaT是开尔文温度e是水汽压hPa。大气探空数据通常能拿到这些量直接从下往上积分。import numpy as np def refractivity(P_hPa, T_K, e_hPa): # 干项 湿项返回 N 值N 单位 dry 77.6 * P_hPa / T_K wet 3.73e5 * e_hPa / T_K**2 return dry wet def modified_refractivity(N, z_km): # 把 N 剖面加上地球曲率项得到 M 剖面 return N 157.0 * z_km逻辑说明干项对温度、气压敏感湿项对水汽压敏感海面附近湿项往往主导了N的骤变。z_km是海拔高度km系数157对应地球半径的修正项。拿到M剖面后用np.gradient求dM/dz找连续负梯度区段区段上界就是波导顶高。我一般会把M剖面和dM/dz画在一起因为只看M曲线经常会忽略负梯度层很薄的细节。参数说明P_hPa单位是hPa不是Pa探空数据一般默认给hPaT_K必须换算成开尔文摄氏度数加273.15e_hPa来自露点温度或者相对湿度换算如果只有相对湿度RH要用e RH × es / 100补一步饱和水汽压计算。这套公式在低对流层误差在1 N单位以内做波导判定足够。3. 抛物方程与分步傅里叶变换把波动方程变成可迭代的步进算法3.1 从亥姆霍兹方程到抛物方程窄角近似的适用边界标量亥姆霍兹方程是电磁波传播分析的基础但直接解它在工程上不现实因为边界条件和折射率剖面都随距离变化。抛物方程法的思路是对亥姆霍兹方程做因式分解只保留前向传播的分量得到窄角抛物方程∂²u/∂z² 2jk₀∂u/∂x k₀²(m² - 1)u 0其中u是约化场强k₀是自由空间波数m是修正折射指数m 1 M×10⁻⁶。这个方程把“波在二维空间中的传播”变成了“沿距离x逐步推进”的初值问题每一小步只依赖前一步的场分布可以循环迭代求解。窄角近似意味着假设传播路径与水平方向的夹角不大。严格来说窄角PE对仰角在几度以内的传播路径精度足够超过约10度就需要换宽角形式。大气波导传播本身就是低仰角现象尤其是蒸发波导几乎贴着海面走仰角通常不超过1度所以窄角PE在这个场景下非常对路。宽角PE在算子形式上多了一项修正传播因子变成exp(-jΔx·kz²/(k₀ √(k₀² - kz²)))能支撑到20度左右的仰角范围。学教案里通常先讲窄角再提宽角。实际我仿真时如果仰角库里包含5度以上的分量直接上宽角形式不纠结。3.2 分步傅里叶变换的迭代格式与代码骨架窄角PE的标准解法是分步傅里叶变换Split-Step Fourier Transform, SSFT。核心思想是每一步把传播算子拆成两项衍射项在谱域处理折射项在空间域处理交替进行。我直接把骨架函数写出来import numpy as np def ssft_step(u, k0, m_profile, dz, dx): # u: 当前距离处的约化场分布复数数组 # k0: 自由空间波数 2*pi/lambda # m_profile: 修正折射指数剖面 1 M*1e-6 # dz: 垂直网格间隔单位米 # dx: 距离步长单位米 # 1. 衍射项在垂直波数谱域做相位累积 kz 2 * np.pi * np.fft.fftfreq(u.size, ddz) u_hat np.fft.fft(u) u_hat u_hat * np.exp(-1j * kz**2 * dx / (2 * k0)) # 2. 折射项在空间域做相位累积 u np.fft.ifft(u_hat) u u * np.exp(1j * k0 * (m_profile - 1.0) * dx) return u逻辑说明先对场做FFT得到垂直方向的空间频谱谱域乘的相位因子对应自由空间衍射会把不同角度分量的相位往前推逆FFT回到空间域后再乘折射项相位因子这一步把大气折射率随高度的变化带进来。两步合起来一个迭代步就完成了。参数说明kz用np.fft.fftfreq生成是弧度/米不需要手动构造正负频率排列dx要远小于2k0/kz²的量级才能保证相位累积准确常见做法是取几十到上百倍波长m_profile减1.0是因为折射指数接近1直接乘m-1可以避免大数相位误差。每步都要让u保持复数不要取模否则相位信息丢失后面的干涉现象全部消失。3.3 吸收窗与初始场让计算域边界和天线方向图不捣乱计算域顶部必须截断但硬截断会反射能量产生虚假干涉条纹。标准做法是在顶部加Harding窗让场平滑衰减到零。我把窗函数写成独立函数每次迭代后乘一次def harding_window(N, top_frac0.1): # N: 总格点数, top_frac: 顶部衰减层占比 w np.ones(N) z np.arange(N) / (N - 1) mask z 1 - top_frac zz z[mask] - (1 - top_frac) w[mask] np.cos(0.5 * np.pi * zz / top_frac) ** 2 return w逻辑说明w在正常区域保持1只在顶部top_frac范围内按余弦平方降到0。这样入射到顶边界上的能量被逐渐吸收而不是被反射回来。top_frac取0.1比较稳如果计算域高度本身留得很足取0.05也够。初始场代表天线辐射方向图。严格做法是把天线方向图做逆FFT得到高度域分布教案里常用高斯波束近似效果足够def gauss_antenna_field(z, z_tx, k0, beam_3db): # z_tx: 天线高度单位米 # beam_3db: 天线波束3dB宽度单位弧度 sigma_theta beam_3db / 2.354 # 高斯方向图的标准差 sigma_z 1.0 / (k0 * sigma_theta) # 高度域等效波束宽度 return np.exp(-(z - z_tx)**2 / (2 * sigma_z**2))逻辑说明波束越窄在高度域对应的场分布越宽这是傅里叶对偶关系。3dB宽度换算成高斯标准差时要除以2.354这个系数来自高斯函数的半功率点定义。参数说明beam_3db要转成弧度例如1度波束是np.deg2rad(1.0)z_tx设置在波导层内或贴近波导底部才能有效激励波导传播sigma_z如果小于一个网格间隔说明波束在高度域采样严重不足需要加密垂直网格。4. 从M剖面到传播损耗PE建模的完整流程与关键参数4.1 拿到一个场景后先定五件事任何PE仿真开始前先确定场景参数。我给团队做培训时把这五件事固定在白板上缺一个不开算参数典型值设置依据频率3 ~ 18 GHz波导效应随频率增强1 GHz以下基本不用算天线高度5 ~ 30 m低天线更容易将能量耦合进波导层传播距离100 ~ 300 km超视距目标出现的典型距离范围垂直网格λ/10 ~ λ/4决定可解析的仰角上限距离步长50 ~ 100 m步进精度与计算耗时平衡垂直网格设置有个硬约束Δz不能超过λ/2否则垂直波数谱发生混叠会出现根本不存在的传播模式。我通常取Δz λ/4计算域高度按波导顶高的3倍留余量再补上顶部10%的Harding窗衰减层。距离步长没有类似的稳定性约束但也不能过大否则衍射项和折射项的非对易误差会积累常见取值是几十到一百米。4.2 一个可跑的仿真流程蒸发波导案例把前面几节的函数拼在一起就是一个完整的蒸发波导传播损耗仿真。下面这段代码我直接跑过参数做了调整适合教学演示lam 0.032 # 9.375 GHz典型雷达波段 k0 2 * np.pi / lam dz lam / 4 # 垂直网格间隔 N 65536 # 2的幂FFT友好 z (np.arange(N) 0.5) * dz # 偏移半格避开 z0 奇点 H N * dz # 计算域高度约524 m dx 100.0 n_steps 1500 # 150 km 传播距离 # 构造蒸发波导 M 剖面h_duct18 m对数模型 m_profile 1 evap_duct_profile(z, h_duct18.0) * 1e-6 w harding_window(N, top_frac0.1) # 顶部吸收窗 u gauss_antenna_field(z, z_tx15.0, k0k0, beam_3dbnp.deg2rad(1.0)) loss np.zeros((n_steps, N)) for i in range(n_steps): u ssft_step(u, k0, m_profile, dz, dx) u u * w r (i 1) * dx # 传播损耗 自由空间损耗 - 20*log10(|u|) loss[i] 20 * np.log10(4 * np.pi * r / lam) - 20 * np.log10(np.abs(u) 1e-15)逻辑说明z数组偏移了半格避免对数剖面在z0处发散每步迭代后乘Harding窗保证顶部不反射损耗计算用的是自由空间球面扩展损耗减去PE给出的场强修正量|u|1e-15是防止log里出现零值。跑完之后把loss画成距离-高度伪彩图就能看到波导层内的低损耗通道。参数说明N取65536是为了让FFT在计算上高效同时保证高度分辨率足够演示机器慢的话可以把N降到16384dx提到200 m趋势不变h_duct18 m用的是海上夏季典型蒸发波导高度换场景时按探空数据或海气模型改天线高度15 m与波导层重叠这是激励波导的关键条件天线太高波导反而捕不住能量。4.3 读图的三个落点陷阱高度、干涉瓣、绕射尾仿真结果出来后不要只盯着损耗数值要看三个结构特征。第一个是陷阱高度伪彩图上损耗最低的通道所在高度应该与M剖面负梯度层对应如果偏移很大说明折射率剖面或初始场设置有问题。第二个是干涉瓣低空区域会出现明暗交替的条纹这是直射路径和海面反射路径干涉的结果条纹间距由天线高度和频率决定可以用来反推天线架设是否合理。第三个是绕射尾超出波导覆盖距离后损耗曲线会出现更陡的上升段那是能量泄漏出波导层后的绕射衰减区。这三个特征同时出现说明仿真基本可信。如果只看到大片均匀低损耗区而没有干涉结构多半是初始场太宽或者网格太粗把细节抹掉了。5. 避坑指南PE仿真最常见的五个翻车现场5.1 顶部反射导致的伪波导条纹现象损耗图在计算域上部出现规律性亮暗条纹越靠近顶边界越明显甚至看起来像多了一层波导。原因计算域截断产生硬边界场传播到顶部被反射回来与正向场叠加成驻波。常见诱因是Harding窗只乘了一次但代码在迭代循环里每步都要乘或者计算域高度留得太矮真实关注高度已经进入了衰减层。解决先确认窗函数在循环内每步都乘再把计算域高度加到波导顶高的3倍以上Harding窗的衰减层放在实际关心高度之上关掉所有顶部反射再看图。5.2 M与m混用导致相位爆炸现象损耗曲线出现剧烈振荡数值动不动到几百dB甚至出现NaN。原因修正折射率M的值通常是几十到几百M单位而折射指数m是1 M×10⁻⁶两者差了6个数量级。把M直接当m写进相位因子exp(jk₀(m-1)dx)相当于乘了一个巨大的相位数值直接爆炸。解决代码里统一变量命名m_profile 1 M_profile * 1e-6并在仿真的第一行打印m_profile.max()正常应该在1.0008附近如果打印出来是几百说明漏了系数。5.3 初始场未归一化近场损耗失真现象第一个距离步算出来的损耗就比自由空间损耗低20到40 dB远场整体也偏小一个常数。原因高斯波束的峰值设为1但波束的总能量没有归一化u代表的场强修正量偏离真实值。后面乘窗函数后能量还会进一步损失。解决在进入迭代前做一次能量归一化u u / np.sqrt(np.sum(np.abs(u)**2) * dz)或者先跑标准大气剖面拿远场损耗和解析结果对比反推修正常数。5.4 垂直网格太粗谱混叠出假的传播模式现象场图上出现很多条平行的细密能量带间距规律看起来像高次模式但物理上不该存在。原因Δz超过了λ/2垂直波数谱的高频分量被折叠到低频区产生虚假模式。这类问题窄角PE里尤其隐蔽因为相位因子表面上看不出异常。解决把Δz压到λ/4以内N随之增大取2的幂。如果计算域高度太高导致点内存压力大可以降计算域高度或者用非均匀网格而不是放宽Δz。5.5 蒸发波导剖面在z0附近发散现象M剖面曲线在最低几个格点出现巨大负值或NaN仿真第一轮就报错。原因对数模型里有ln(z / z0)z0时对数趋向负无穷而网格第一个点恰好落在z0上。解决构造高度数组时偏移半格z (np.arange(N) 0.5) * dz或者设置z_min0.01 m低于这个值的剖面直接用该点的M值填充。两种做法都对关键是别让z0进对数。6. 进阶验证用三步对照法确认你的PE结果可信6.1 第一步标准大气做底数先别急着上波导剖面。把M剖面设成标准大气即M 0.118 × zz单位m其他参数不变跑一遍同一套代码。取10 km处的传播损耗与两射线模型的结果做对比两射线模型考虑直射波和海面镜面反射在平坦海面、天线高度远大于波长时精度很高误差应该在1 dB以内。m_std 1 0.118 * z * 1e-6这一步看起来多余但它是排除代码bug最直接的方法。如果标准大气都对不上后面所有波导结论都不要采信。6.2 第二步波导高度扫描把蒸发波导高度h_duct分别设为10、15、20、25 m跑四遍观察损耗最强带的高度是否随h_duct同步移动。移动量不要求完全相等因为天线高度固定耦合效率会带来偏移但趋势必须一致。如果h_duct涨了10 m最强带纹丝不动说明剖面函数或者初始场设置有问题。6.3 第三步在损耗斜线上读波导平台取固定接收高度的损耗随距离对数曲线。标准大气情况下曲线斜率接近自由空间的20log10(4πr/λ)出现波导捕获后远距离段的斜率明显变缓甚至出现平台段。我在验收仿真代码时总习惯在损耗图上叠一条自由空间曲线做参照波导平台一眼就能认出来。这套三步对照法每次新场景我都强制走一遍。有一回跳过第一步直接拿波导剖面跑折腾两天才发现折射项里漏了1e-6系数数据彻底不能看。从那以后我每拿一份PE学习资料或接手一套仿真代码第一件事就是跑标准大气标定再谈波导。这套流程走顺再碰表面波导、抬升波导和复杂地形组合才不会把偶然的数值巧合当成物理规律。希望帮到你。本文还有配套的精品资源点击获取
简介针对大气波导环境下电波传播研究的专业教学演示文稿面向通信工程、电磁场与微波技术、大气科学等方向的学生、教师与科研人员可作为课堂讲义、自学笔记或课题起步的参考资料。内容共25页从大气折射的基本类型及存在条件入手逐步深入到修正折射率计算、陷获折射与大气波导形成条件再系统讲解抛物方程方法的推导过程、分步傅里叶变换实现、粗糙海面阻抗边界设置、初始场处理以及电波传播损耗的计算方式。同时演示文稿专门给出MATLAB电波传播仿真模块帮助读者掌握如何将理论模型转化为可计算的仿真流程并对传播损耗结果进行可视化分析。该演示文稿为单个pptx文件压缩包仅472KB结构紧凑、便于携带和二次编辑现有486人学习尤其适合开展超视距传播机理、雷达探测距离评估和通信链路设计等方向的研究人员参考能够快速建立从大气波导物理机制到数值建模与仿真的完整认知。1. 大气波导不是玄学为什么研究电波传播要先啃抛物方程这是一份可以直接当讲义用的专业资料主题是大气波导环境下的电波传播核心解法是抛物方程Parabolic Equation, PE法。干雷达链路、微波通信或电磁仿真的人多半遇到过这种场景海上目标在视距之外几十公里被雷达抓到或者通信链路在某个季节突然多出几十公里余量——第一个怀疑对象就是大气波导。可波导这个东西看不见摸不着用射线法解释又漏洞百出真正能落地算损耗的还是抛物方程。这套学习教案的价值在于不绕开数学也不停在公式层面把折射率剖面、PE推导、分步傅里叶变换实现和结果判读串成一条可复现的线。适合电磁场专业研究生、雷达链路设计工程师以及想给团队做内部培训的算法负责人。2. 大气波导环境参数修正折射率剖面与三类波导的判定2.1 为什么要看修正折射率M而不是直接看折射率n大气折射率n在数值上接近1直接看它很难看出波导结构。工程上习惯用折射率差N (n - 1)×10⁶单位是N单位但在考虑地球曲率后还要把高度项加进去得到修正折射率MM N (z / a)×10⁶其中a 6370 kmz是离地高度km。这样做的原因很朴素电磁波在球面地球上传播射线弯曲是大气折射和地球曲率共同作用的结果。把地球“扳直”之后判断波导只需要看M随高度变化的斜率dM/dz而不需要拆开看两条曲线。dM/dz是判定折射类型的核心指标。标准大气条件下低层大气N梯度约为-39 N单位/km折算后dM/dz约为0.118 M单位/m。当dM/dz低于这个值射线向下弯曲得更厉害进入超折射状态一旦dM/dz降到0以下电磁波就被限制在一层“管道”里反复折射这就是波导捕获。为了快速查判我习惯用下面这张表dM/dz 范围 (M单位/m)折射类型对传播的影响 0.118亚折射射线向上弯曲传播距离被压缩≈ 0.118标准折射正常视距传播0 ~ 0.118超折射低空路径弯曲雷达出现超视距现象 0波导捕获能量被限制在波导层内低损耗远距离传播注意这里的0.118是从标准大气推出来的参考值实际大气波动很大海面、沙漠、高原都不一样。拿到任何气象数据后第一步先画dM/dz剖面看有没有连续负梯度区段这是判断波导是否成立的唯一硬标准。2.2 三类典型波导蒸发波导、表面波导与抬升波导大气波导按形成高度可以分成三类工程上判读场景时先分清是哪种因为它们的剖面形状、影响频段和建模方式差别很大。蒸发波导是最常见的海上波导发生在海面以上几米到几十米由海面水汽蒸发造成湿度随高度骤降引起。它的特点是贴地且普遍存在热带海域出现概率极高对C波段以上频率影响显著。剖面形状近似为低层M随高度下降到波导顶之后恢复标准梯度。表面波导是接地或贴近地面的波导层顶高一般在几十米到几百米由夜间辐射冷却、逆温层等造成。它的特征是M剖面从地面开始就有超折射或负梯度上面再叠一个正梯度层。陆地大雾、夜间稳定大气时容易出这种波导对地面雷达和移动通信基站影响最大。抬升波导不接地波导层悬在几百米到上千米的高度由下沉逆温、锋面抬升等过程形成。它的M剖面在波导层内呈负梯度上下两侧都是标准梯度。这种波导对低空目标影响不大但会漏进高空目标容易让雷达在仰角上出现奇怪的回波。波导类型典型高度成因M剖面特征主要影响场景蒸发波导5 ~ 40 m海面水汽蒸发贴地负梯度快速恢复舰载雷达、海上微波链路表面波导50 ~ 500 m逆温、辐射冷却从地面开始负梯度地面雷达、移动通信抬升波导500 ~ 2000 m下沉逆温、锋面悬空负梯度层高空目标探测、干扰分析2.3 从实测气象数据构造M剖面建模的第一步是把温度、气压、湿度廓线转成M剖面。折射率N的经典计算公式是N 77.6 × P / T 3.73×10⁵ × e / T²其中P是气压hPaT是开尔文温度e是水汽压hPa。大气探空数据通常能拿到这些量直接从下往上积分。import numpy as np def refractivity(P_hPa, T_K, e_hPa): # 干项 湿项返回 N 值N 单位 dry 77.6 * P_hPa / T_K wet 3.73e5 * e_hPa / T_K**2 return dry wet def modified_refractivity(N, z_km): # 把 N 剖面加上地球曲率项得到 M 剖面 return N 157.0 * z_km逻辑说明干项对温度、气压敏感湿项对水汽压敏感海面附近湿项往往主导了N的骤变。z_km是海拔高度km系数157对应地球半径的修正项。拿到M剖面后用np.gradient求dM/dz找连续负梯度区段区段上界就是波导顶高。我一般会把M剖面和dM/dz画在一起因为只看M曲线经常会忽略负梯度层很薄的细节。参数说明P_hPa单位是hPa不是Pa探空数据一般默认给hPaT_K必须换算成开尔文摄氏度数加273.15e_hPa来自露点温度或者相对湿度换算如果只有相对湿度RH要用e RH × es / 100补一步饱和水汽压计算。这套公式在低对流层误差在1 N单位以内做波导判定足够。3. 抛物方程与分步傅里叶变换把波动方程变成可迭代的步进算法3.1 从亥姆霍兹方程到抛物方程窄角近似的适用边界标量亥姆霍兹方程是电磁波传播分析的基础但直接解它在工程上不现实因为边界条件和折射率剖面都随距离变化。抛物方程法的思路是对亥姆霍兹方程做因式分解只保留前向传播的分量得到窄角抛物方程∂²u/∂z² 2jk₀∂u/∂x k₀²(m² - 1)u 0其中u是约化场强k₀是自由空间波数m是修正折射指数m 1 M×10⁻⁶。这个方程把“波在二维空间中的传播”变成了“沿距离x逐步推进”的初值问题每一小步只依赖前一步的场分布可以循环迭代求解。窄角近似意味着假设传播路径与水平方向的夹角不大。严格来说窄角PE对仰角在几度以内的传播路径精度足够超过约10度就需要换宽角形式。大气波导传播本身就是低仰角现象尤其是蒸发波导几乎贴着海面走仰角通常不超过1度所以窄角PE在这个场景下非常对路。宽角PE在算子形式上多了一项修正传播因子变成exp(-jΔx·kz²/(k₀ √(k₀² - kz²)))能支撑到20度左右的仰角范围。学教案里通常先讲窄角再提宽角。实际我仿真时如果仰角库里包含5度以上的分量直接上宽角形式不纠结。3.2 分步傅里叶变换的迭代格式与代码骨架窄角PE的标准解法是分步傅里叶变换Split-Step Fourier Transform, SSFT。核心思想是每一步把传播算子拆成两项衍射项在谱域处理折射项在空间域处理交替进行。我直接把骨架函数写出来import numpy as np def ssft_step(u, k0, m_profile, dz, dx): # u: 当前距离处的约化场分布复数数组 # k0: 自由空间波数 2*pi/lambda # m_profile: 修正折射指数剖面 1 M*1e-6 # dz: 垂直网格间隔单位米 # dx: 距离步长单位米 # 1. 衍射项在垂直波数谱域做相位累积 kz 2 * np.pi * np.fft.fftfreq(u.size, ddz) u_hat np.fft.fft(u) u_hat u_hat * np.exp(-1j * kz**2 * dx / (2 * k0)) # 2. 折射项在空间域做相位累积 u np.fft.ifft(u_hat) u u * np.exp(1j * k0 * (m_profile - 1.0) * dx) return u逻辑说明先对场做FFT得到垂直方向的空间频谱谱域乘的相位因子对应自由空间衍射会把不同角度分量的相位往前推逆FFT回到空间域后再乘折射项相位因子这一步把大气折射率随高度的变化带进来。两步合起来一个迭代步就完成了。参数说明kz用np.fft.fftfreq生成是弧度/米不需要手动构造正负频率排列dx要远小于2k0/kz²的量级才能保证相位累积准确常见做法是取几十到上百倍波长m_profile减1.0是因为折射指数接近1直接乘m-1可以避免大数相位误差。每步都要让u保持复数不要取模否则相位信息丢失后面的干涉现象全部消失。3.3 吸收窗与初始场让计算域边界和天线方向图不捣乱计算域顶部必须截断但硬截断会反射能量产生虚假干涉条纹。标准做法是在顶部加Harding窗让场平滑衰减到零。我把窗函数写成独立函数每次迭代后乘一次def harding_window(N, top_frac0.1): # N: 总格点数, top_frac: 顶部衰减层占比 w np.ones(N) z np.arange(N) / (N - 1) mask z 1 - top_frac zz z[mask] - (1 - top_frac) w[mask] np.cos(0.5 * np.pi * zz / top_frac) ** 2 return w逻辑说明w在正常区域保持1只在顶部top_frac范围内按余弦平方降到0。这样入射到顶边界上的能量被逐渐吸收而不是被反射回来。top_frac取0.1比较稳如果计算域高度本身留得很足取0.05也够。初始场代表天线辐射方向图。严格做法是把天线方向图做逆FFT得到高度域分布教案里常用高斯波束近似效果足够def gauss_antenna_field(z, z_tx, k0, beam_3db): # z_tx: 天线高度单位米 # beam_3db: 天线波束3dB宽度单位弧度 sigma_theta beam_3db / 2.354 # 高斯方向图的标准差 sigma_z 1.0 / (k0 * sigma_theta) # 高度域等效波束宽度 return np.exp(-(z - z_tx)**2 / (2 * sigma_z**2))逻辑说明波束越窄在高度域对应的场分布越宽这是傅里叶对偶关系。3dB宽度换算成高斯标准差时要除以2.354这个系数来自高斯函数的半功率点定义。参数说明beam_3db要转成弧度例如1度波束是np.deg2rad(1.0)z_tx设置在波导层内或贴近波导底部才能有效激励波导传播sigma_z如果小于一个网格间隔说明波束在高度域采样严重不足需要加密垂直网格。4. 从M剖面到传播损耗PE建模的完整流程与关键参数4.1 拿到一个场景后先定五件事任何PE仿真开始前先确定场景参数。我给团队做培训时把这五件事固定在白板上缺一个不开算参数典型值设置依据频率3 ~ 18 GHz波导效应随频率增强1 GHz以下基本不用算天线高度5 ~ 30 m低天线更容易将能量耦合进波导层传播距离100 ~ 300 km超视距目标出现的典型距离范围垂直网格λ/10 ~ λ/4决定可解析的仰角上限距离步长50 ~ 100 m步进精度与计算耗时平衡垂直网格设置有个硬约束Δz不能超过λ/2否则垂直波数谱发生混叠会出现根本不存在的传播模式。我通常取Δz λ/4计算域高度按波导顶高的3倍留余量再补上顶部10%的Harding窗衰减层。距离步长没有类似的稳定性约束但也不能过大否则衍射项和折射项的非对易误差会积累常见取值是几十到一百米。4.2 一个可跑的仿真流程蒸发波导案例把前面几节的函数拼在一起就是一个完整的蒸发波导传播损耗仿真。下面这段代码我直接跑过参数做了调整适合教学演示lam 0.032 # 9.375 GHz典型雷达波段 k0 2 * np.pi / lam dz lam / 4 # 垂直网格间隔 N 65536 # 2的幂FFT友好 z (np.arange(N) 0.5) * dz # 偏移半格避开 z0 奇点 H N * dz # 计算域高度约524 m dx 100.0 n_steps 1500 # 150 km 传播距离 # 构造蒸发波导 M 剖面h_duct18 m对数模型 m_profile 1 evap_duct_profile(z, h_duct18.0) * 1e-6 w harding_window(N, top_frac0.1) # 顶部吸收窗 u gauss_antenna_field(z, z_tx15.0, k0k0, beam_3dbnp.deg2rad(1.0)) loss np.zeros((n_steps, N)) for i in range(n_steps): u ssft_step(u, k0, m_profile, dz, dx) u u * w r (i 1) * dx # 传播损耗 自由空间损耗 - 20*log10(|u|) loss[i] 20 * np.log10(4 * np.pi * r / lam) - 20 * np.log10(np.abs(u) 1e-15)逻辑说明z数组偏移了半格避免对数剖面在z0处发散每步迭代后乘Harding窗保证顶部不反射损耗计算用的是自由空间球面扩展损耗减去PE给出的场强修正量|u|1e-15是防止log里出现零值。跑完之后把loss画成距离-高度伪彩图就能看到波导层内的低损耗通道。参数说明N取65536是为了让FFT在计算上高效同时保证高度分辨率足够演示机器慢的话可以把N降到16384dx提到200 m趋势不变h_duct18 m用的是海上夏季典型蒸发波导高度换场景时按探空数据或海气模型改天线高度15 m与波导层重叠这是激励波导的关键条件天线太高波导反而捕不住能量。4.3 读图的三个落点陷阱高度、干涉瓣、绕射尾仿真结果出来后不要只盯着损耗数值要看三个结构特征。第一个是陷阱高度伪彩图上损耗最低的通道所在高度应该与M剖面负梯度层对应如果偏移很大说明折射率剖面或初始场设置有问题。第二个是干涉瓣低空区域会出现明暗交替的条纹这是直射路径和海面反射路径干涉的结果条纹间距由天线高度和频率决定可以用来反推天线架设是否合理。第三个是绕射尾超出波导覆盖距离后损耗曲线会出现更陡的上升段那是能量泄漏出波导层后的绕射衰减区。这三个特征同时出现说明仿真基本可信。如果只看到大片均匀低损耗区而没有干涉结构多半是初始场太宽或者网格太粗把细节抹掉了。5. 避坑指南PE仿真最常见的五个翻车现场5.1 顶部反射导致的伪波导条纹现象损耗图在计算域上部出现规律性亮暗条纹越靠近顶边界越明显甚至看起来像多了一层波导。原因计算域截断产生硬边界场传播到顶部被反射回来与正向场叠加成驻波。常见诱因是Harding窗只乘了一次但代码在迭代循环里每步都要乘或者计算域高度留得太矮真实关注高度已经进入了衰减层。解决先确认窗函数在循环内每步都乘再把计算域高度加到波导顶高的3倍以上Harding窗的衰减层放在实际关心高度之上关掉所有顶部反射再看图。5.2 M与m混用导致相位爆炸现象损耗曲线出现剧烈振荡数值动不动到几百dB甚至出现NaN。原因修正折射率M的值通常是几十到几百M单位而折射指数m是1 M×10⁻⁶两者差了6个数量级。把M直接当m写进相位因子exp(jk₀(m-1)dx)相当于乘了一个巨大的相位数值直接爆炸。解决代码里统一变量命名m_profile 1 M_profile * 1e-6并在仿真的第一行打印m_profile.max()正常应该在1.0008附近如果打印出来是几百说明漏了系数。5.3 初始场未归一化近场损耗失真现象第一个距离步算出来的损耗就比自由空间损耗低20到40 dB远场整体也偏小一个常数。原因高斯波束的峰值设为1但波束的总能量没有归一化u代表的场强修正量偏离真实值。后面乘窗函数后能量还会进一步损失。解决在进入迭代前做一次能量归一化u u / np.sqrt(np.sum(np.abs(u)**2) * dz)或者先跑标准大气剖面拿远场损耗和解析结果对比反推修正常数。5.4 垂直网格太粗谱混叠出假的传播模式现象场图上出现很多条平行的细密能量带间距规律看起来像高次模式但物理上不该存在。原因Δz超过了λ/2垂直波数谱的高频分量被折叠到低频区产生虚假模式。这类问题窄角PE里尤其隐蔽因为相位因子表面上看不出异常。解决把Δz压到λ/4以内N随之增大取2的幂。如果计算域高度太高导致点内存压力大可以降计算域高度或者用非均匀网格而不是放宽Δz。5.5 蒸发波导剖面在z0附近发散现象M剖面曲线在最低几个格点出现巨大负值或NaN仿真第一轮就报错。原因对数模型里有ln(z / z0)z0时对数趋向负无穷而网格第一个点恰好落在z0上。解决构造高度数组时偏移半格z (np.arange(N) 0.5) * dz或者设置z_min0.01 m低于这个值的剖面直接用该点的M值填充。两种做法都对关键是别让z0进对数。6. 进阶验证用三步对照法确认你的PE结果可信6.1 第一步标准大气做底数先别急着上波导剖面。把M剖面设成标准大气即M 0.118 × zz单位m其他参数不变跑一遍同一套代码。取10 km处的传播损耗与两射线模型的结果做对比两射线模型考虑直射波和海面镜面反射在平坦海面、天线高度远大于波长时精度很高误差应该在1 dB以内。m_std 1 0.118 * z * 1e-6这一步看起来多余但它是排除代码bug最直接的方法。如果标准大气都对不上后面所有波导结论都不要采信。6.2 第二步波导高度扫描把蒸发波导高度h_duct分别设为10、15、20、25 m跑四遍观察损耗最强带的高度是否随h_duct同步移动。移动量不要求完全相等因为天线高度固定耦合效率会带来偏移但趋势必须一致。如果h_duct涨了10 m最强带纹丝不动说明剖面函数或者初始场设置有问题。6.3 第三步在损耗斜线上读波导平台取固定接收高度的损耗随距离对数曲线。标准大气情况下曲线斜率接近自由空间的20log10(4πr/λ)出现波导捕获后远距离段的斜率明显变缓甚至出现平台段。我在验收仿真代码时总习惯在损耗图上叠一条自由空间曲线做参照波导平台一眼就能认出来。这套三步对照法每次新场景我都强制走一遍。有一回跳过第一步直接拿波导剖面跑折腾两天才发现折射项里漏了1e-6系数数据彻底不能看。从那以后我每拿一份PE学习资料或接手一套仿真代码第一件事就是跑标准大气标定再谈波导。这套流程走顺再碰表面波导、抬升波导和复杂地形组合才不会把偶然的数值巧合当成物理规律。希望帮到你。本文还有配套的精品资源点击获取
