短波航空移动信道仿真：Watterson模型改进与定制化航迹实现

简介这份文档面向短波通信与航空信道建模方向的研究生、科研人员及工程技术人员聚焦远距离航空移动通信中短波信道衰落特性复杂、缺乏标准模型的问题。内容以Watterson抽头延迟线模型为基础引入飞行器相对运动产生的多普勒频移推导时变频响与抽头增益函数并针对不同种类、不同参数飞行器开展差异化仿真在航线航迹已知时还可实现特定民用航空场景的定制化信道仿真。资源包为单个docx文档约526KB内含引言、Watterson信道模型、短波航空移动信道建模等章节配有公式推导与模型结构说明便于读者直接理解建模思路并复现仿真流程。目前已有159人学习适合需要快速掌握短波航空信道建模方法、撰写相关论文或搭建仿真验证环境的研究者参考。1. 短波航空移动信道仿真从 Watterson 模型到定制化航迹的完整复现路径短波频段330 MHz靠电离层反射实现超视距通信在卫星链路不可用时它往往是飞行器远程指挥控制的唯一手段。但飞行器一动收发两端的相对运动就会叠加时变多普勒频移和扩展再叠上电离层本身的多径与时延信道衰落特性比固定台站复杂得多。这份文档给出的方案是在经典 Watterson 模型基础上引入飞行器运动状态参数把多普勒效应拆成电离层分量和相对运动分量并针对民用航空、私人飞机、无人载具三类机动频率做了差异化仿真。它适合做短波链路预算、通信体制验证、半实物仿真前端设计的工程师也适合需要快速搭出可调信道模型的研究生。下面按“模型原理 → 参数怎么定 → 代码怎么落 → 坑在哪 → 进阶怎么用”的顺序拆开讲。2. Watterson 模型为什么能当基线抽头延迟线结构与三个核心假设2.1 抽头延迟线到底在模拟什么Watterson 模型的结构并不复杂输入信号进入一条抽头延迟线每个抽头对应一种电离层传播模式或路径抽头增益函数 (G_i(t)) 对延迟后的信号做调制模拟该路径上的衰落最后各路相加再叠加性噪声。它的时变频响写成[ H(f,t)\sum_{i1}^{n}\exp(-j2\pi\tau_i f)G_i(t) ]其中 (i) 是路径标号(n) 是路径总数(\tau_i) 是第 (i) 条路径的时延(G_i(t)) 是第 (i) 条路径的抽头增益函数。这个式子把“多径”和“时变”两件事同时装进了一个线性系统里实现上就是延迟线加复数乘法器复杂度可控。Watterson 等人当年基于长期实测给出了三条关键假设每条路径增益是复高斯随机过程各路径增益函数相互独立每条路径增益对应频谱是两个高斯谱的叠加。这三条假设决定了后续所有参数的计算方式也决定了模型的适用边界——有限时间小于 10 分钟和有限带宽小于 12 kHz内信道可视为稳定。超出这个范围模型就不再是“准确”的只能算“近似”。2.2 抽头增益函数的数学形式与物理含义路径增益的表达式为[ G_i(t)G_{ia}(t)\exp(j2\pi f_{ia}t)G_{ib}(t)\exp(j2\pi f_{ib}t) ](G_{ia}) 和 (G_{ib}) 是相互独立的复高斯随机过程均值为零包络服从瑞利分布相位服从均匀分布。(f_{ia}) 和 (f_{ib}) 是第 (i) 条路径的两个多普勒频移分量对应电离层中两个磁离子分量。其自相关函数为[ C_i(\Delta t)C_{ia}(0)\exp[-2\pi^2\sigma_{ia}^2(\Delta t)^2j2\pi v_{ia}\Delta t]C_{ib}(0)\exp[-2\pi^2\sigma_{ib}^2(\Delta t)^2j2\pi v_{ib}\Delta t] ]对应的抽头增益谱函数是两个高斯谱的叠加[ f_i(f)\frac{C_{ia}(0)}{\sqrt{2\pi}\sigma_{ia}}\exp\left[-\frac{(f-f_{ia})^2}{2\sigma_{ia}^2}\right]\frac{C_{ib}(0)}{\sqrt{2\pi}\sigma_{ib}}\exp\left[-\frac{(f-f_{ib})^2}{2\sigma_{ib}^2}\right] ](\sigma_{ia})、(\sigma_{ib}) 是多普勒扩展(C_{ia}(0))、(C_{ib}(0)) 是自相关函数在零时延处的值代表功率。载频较低时两个磁离子分量的多普勒频移和扩展几乎一致功率谱也几乎重合仿真时通常简化成一个高斯谱。这一步简化在 6 MHz 以下基本不会带来可观测误差但在 15 MHz 以上要谨慎两个分量的分离度会变大。2.3 为什么航空场景不能直接套用固定台站参数固定短波台站的多普勒频移主要来自电离层本身的快速运动和反射层高度变化量级通常在 0.11 Hz。而飞行器以 250 m/s 速度飞行、载频 15 MHz 时相对运动产生的多普勒频移可达[ f_A\frac{f_c}{c}v\cos\theta_i\frac{15\times10^6}{3\times10^8}\times250\times\cos\theta_i\approx12.5\cos\theta_i\ \text{Hz} ]这个量级已经远超电离层分量不能再忽略。更麻烦的是飞行器做变速或圆周运动时(\theta_i) 和 (v) 都在变多普勒频移在每条路径内快速变化导致多普勒扩展在短时间内急剧变化。固定台站那套“扩展取常数”的做法在航空场景下会直接把信道模型的时变特性抹平。3. 航空移动信道建模把飞行器运动参数塞进多普勒频移3.1 多普勒频移的两分量拆分航空移动信道的每条路径多普勒频移由两部分组成[ f_if_{iA}f_{iB} ](f_{iB}) 是电离层产生的分量(f_{iA}) 是相对运动产生的分量计算式为[ f_{iA}\frac{f_c}{c}v\cos\theta_i ](f_c) 是载波频率(c) 是光速(v) 是飞行器运动速率(\theta_i) 是接收端入射电波与其运动方向的夹角。不同路径下的 (\theta_i) 往往不同常见做法是假设它在一定范围内服从均匀分布。这一步是航空模型和固定模型的分水岭固定模型里 (f_{iA}0)航空模型里它是主项。3.2 飞行器速度模型与最大速度约束飞行器速度写成[ v(t)v_0a(t)\cos\alpha(t)\cdot t ](v_0) 是初速度(a(t)) 是加速度(\alpha(t)) 是平直平面内速度与加速度方向的夹角。(\alpha0) 时做直线运动其余情况做圆周运动。这个式子能覆盖变速直线、匀速圆周、变速圆周等常见机动。但飞行器受发动机功率限制不可能一直加速。发动机功率为[ PFv(fma)v(kv^2ma)v ](fkv^2) 是空气阻力(k) 是比例系数(m) 是质量。当 (|v/v_{\max}|\ge0.8) 时假设发动机功率保持最大恒定[ P_{\max}kv_{\max}^3kv^3mav ]推出加速度衰减规律[ ak_0\left(\frac{1}{x}-x^2\right),\quad x\in[0.8,1] ]其中 (x|v/v_{\max}|)(k_0v_{\max}^2k/m)。这一步很关键如果不加这个约束仿真里飞行器速度会无限增长多普勒频移直接发散频谱图会变成一条不断外扩的斜线完全失真。3.3 三类机动频率的参数参考文档按机动频率等级给出了三类典型场景的参考值机动频率等级典型应用场景机动频率参考值Hz飞行状态持续时间s低民用航空0.01≥200中私人飞机等0.11020高无人载具等1≤1机动频率越高飞行器速度和加速度变化越快额外的多普勒频移和扩展越大。仿真时假设最大加速度 80 m/s²、最大飞行速度 600 m/s随机生成 100 s 速度变化可以看到低机动频率下速度曲线平滑高机动频率下速度在短时间内剧烈波动。这个表是后续所有差异化仿真的参数入口改一个数就能切换场景。4. 仿真实现希尔伯特滤波器、时变多普勒与扩展的代码落地4.1 希尔伯特滤波器与 I/Q 两路系数生成实际短波信号载频在 330 MHz直接处理不现实常见做法是数字下变频分离出 I、Q 两路基带信号。带通滤波器的设计思路是先设计一个 FIR 低通滤波器通带为所需带通滤波器通带的 1/2再转换成 I、Q 两路系数[ h_{IBP}(n)2h_{LP}(n)\cos\left(2\pi f_0\left[n-\frac{N-1}{2}\right]T\right) ][ h_{QBP}(n)2h_{LP}(n)\sin\left(2\pi f_0\left[n-\frac{N-1}{2}\right]T\right) ](h_{LP}(n)) 是低通滤波器系数(f_0) 是通带中心频率(N) 是滤波器阶数(T) 是采样周期。用 Python 实现时低通滤波器可以用scipy.signal.firwin生成再按上式转成 I、Q 系数import numpy as np from scipy.signal import firwin def design_iq_bandpass(numtaps, cutoff, fs, f0): numtaps: 滤波器阶数 cutoff: 低通截止频率Hz为带通带宽的一半 fs: 采样率Hz f0: 带通中心频率Hz h_lp firwin(numtaps, cutoff, fsfs) n np.arange(numtaps) delay (numtaps - 1) / 2 h_i 2 * h_lp * np.cos(2 * np.pi * f0 * (n - delay) / fs) h_q 2 * h_lp * np.sin(2 * np.pi * f0 * (n - delay) / fs) return h_i, h_qnumtaps取 64128 之间通常够用cutoff按信号带宽的一半设f0是下变频后的中心频率。注意f0不能超过fs/2否则余弦项会混叠。这一步的常见翻车点是滤波器阶数取太小导致 I、Q 两路幅度不一致后续频谱会出现镜像分量。4.2 时变多普勒频移的实现时变多普勒频移由电离层分量和相对运动分量叠加而成。让一个复信号频移 (\Delta f)等价于时域乘以 (\exp(j2\pi\Delta f t))。假设输入信号为 (\exp(j2\pi ft))频移后为 (\exp(j2\pi(f\Delta f)t))。代码实现def apply_doppler(signal, fs, fc, v, theta, f_ion): signal: 输入复基带信号 fs: 采样率 fc: 载波频率 v: 飞行器速度序列m/s长度与 signal 一致 theta: 入射角序列rad长度与 signal 一致 f_ion: 电离层多普勒频移Hz c 3e8 t np.arange(len(signal)) / fs f_rel fc / c * v * np.cos(theta) f_total f_rel f_ion phase 2 * np.pi * np.cumsum(f_total) / fs return signal * np.exp(1j * phase)这里用np.cumsum(f_total)/fs而不是f_total * t是因为 (f_total) 本身是时变的直接乘时间会引入相位误差。v和theta按 3.2 节的速度模型生成f_ion按经验取 0.11 Hz。如果飞行器做圆周运动theta要按运动方向实时更新不能取常数。4.3 多普勒扩展的高斯滤波实现Watterson 模型假设多普勒扩展功率谱服从高斯分布。高斯白噪声功率谱均匀用高斯滤波器滤一下就能得到高斯功率谱噪声序列再与输入信号相乘时域相乘对应频域卷积实现频谱扩展。噪声采样率太高会导致高斯滤波器阶数过大常见做法是用低采样率噪声序列插值滤波到系统采样率from scipy.signal import gaussian, lfilter from scipy.interpolate import interp1d def apply_doppler_spread(signal, fs, sigma, fs_noise1000): signal: 输入复基带信号 fs: 系统采样率 sigma: 多普勒扩展Hz fs_noise: 噪声生成采样率低于 fs 以降低滤波器阶数 n_noise int(len(signal) / fs * fs_noise) noise (np.random.randn(n_noise) 1j * np.random.randn(n_noise)) / np.sqrt(2) # 高斯滤波器标准差按 sigma 映射到噪声采样率 taps gaussian(int(6 * fs_noise / (2 * np.pi * sigma)), stdfs_noise / (2 * np.pi * sigma)) taps / np.sum(taps) noise_filt lfilter(taps, 1.0, noise) # 插值到系统采样率 t_noise np.arange(n_noise) / fs_noise t_sys np.arange(len(signal)) / fs interp interp1d(t_noise, noise_filt, kindlinear, fill_valueextrapolate) noise_sys interp(t_sys) return signal * noise_syssigma按路径分别取 0.5 Hz、1 Hz、1.5 Hz对应三条路径。fs_noise取 1000 Hz 通常够太低会丢失扩展细节太高滤波器阶数爆炸。插值用线性即可高阶插值在噪声序列上容易过冲。4.4 完整仿真链路与参数设置把上面几块串起来仿真参数参考 ITU 对典型短波电离层反射信道的建议载波频率 6 MHz 和 15 MHz 各跑一遍输入信号为从高频下变频到 200 Hz 的单音三条路径时延分别为 0、1 ms、2 ms多普勒扩展 0.5 Hz、1 Hz、1.5 Hz信噪比 10 dB。完整链路def simulate_hf_aviation_channel(fc, fs, duration, path_delays, sigmas, snr_db, v_seq, theta_seq, f_ion): t np.arange(0, duration, 1/fs) sig np.exp(1j * 2 * np.pi * 200 * t) # 200 Hz 单音 out np.zeros_like(sig, dtypecomplex) for tau, sigma in zip(path_delays, sigmas): delay_samples int(tau * fs) delayed np.concatenate([np.zeros(delay_samples, dtypecomplex), sig[:-delay_samples]]) if delay_samples 0 else sig path_sig apply_doppler(delayed, fs, fc, v_seq, theta_seq, f_ion) path_sig apply_doppler_spread(path_sig, fs, sigma) out path_sig # 加噪声 sig_power np.mean(np.abs(out)**2) noise_power sig_power / (10**(snr_db/10)) noise np.sqrt(noise_power/2) * (np.random.randn(len(out)) 1j*np.random.randn(len(out))) return out noise跑完后画冲激响应和频谱能看到同一机动频率下不同时刻衰落不同机动频率越高信号失真越明显15 MHz 下频谱扩展比 6 MHz 更显著。这些现象和理论一致说明链路搭对了。5. 避坑与排查仿真跑飞、频谱异常、参数不收敛的常见原因5.1 速度发散导致频谱无限外扩现象频谱图上信号能量不断向两侧扩散时间越长越宽完全看不出高斯谱形状。原因速度模型里没有加最大速度约束飞行器一直加速多普勒频移随时间线性增长。解决按 3.2 节的功率约束当 (|v/v_{\max}|\ge0.8) 时切换到加速度衰减公式确保速度收敛到 (v_{\max})。我一般会在速度生成后加一句断言assert np.max(np.abs(v)) v_max * 1.01跑飞了直接报错。5.2 多普勒扩展滤波器阶数过大导致内存爆掉现象程序卡在lfilter或直接 MemoryError。原因fs_noise设得和fs一样高高斯滤波器阶数按fs_noise/sigma算出来几千阶卷积时内存扛不住。解决把fs_noise降到 1000 Hz 左右滤波后再插值回系统采样率。插值带来的误差在sigma大于 0.1 Hz 时可忽略。5.3 I/Q 两路幅度不一致导致镜像分量现象频谱上除了有用信号还出现一个对称的镜像峰。原因希尔伯特滤波器设计时numtaps取太小或者f0接近fs/2余弦和正弦项的正交性被破坏。解决numtaps至少取 64f0不超过fs/2的 0.4 倍。生成 I、Q 系数后检查np.sum(h_i**2)和np.sum(h_q**2)是否接近相等差超过 5% 就加大阶数。5.4 相位累积用错导致频移方向反了现象多普勒频移应该是正的频谱却往负方向偏。原因apply_doppler里用了f_total * t而不是np.cumsum(f_total)/fs当f_total时变时相位计算错误。解决统一用累积和计算相位。另外检查theta的定义cos(theta)为正表示接近为负表示远离符号搞反了频移方向就反。5.5 机动频率切换后飞行状态持续时间不匹配现象高机动场景下仿真 100 s速度曲线在 1 s 内就完成了所有变化后面 99 s 几乎静止。原因表 1 里高机动频率对应的飞行状态持续时间是 ≤1 s仿真时长设 100 s 会导致大部分时间无机动。解决按场景调整仿真时长低机动跑 200 s 以上中机动跑 1020 s高机动跑 1 s 左右。或者把多个机动片段拼接起来模拟连续机动。6. 定制化航迹仿真已知航线时怎么把信道模型钉到具体场景民用航空的航线航迹往往是确定的这时候可以把运动参数直接写死做定制化仿真。文档给了一个典型三段式航迹05 s 匀加速直线初速度 0加速度 50 m/s²510 s 匀速圆周半径 1000 m速度 250 m/s加速度大小 (v^2/R) 方向时变1015 s 匀速直线速度 250 m/s加速度 0。仿真参数和前面一致取 Watterson 模型的一条径和航空移动信道的一条径在 5 s 内做频谱对比。实现时把v_seq和theta_seq按时间段分段生成def generate_trajectory(fs, duration): t np.arange(0, duration, 1/fs) v np.zeros_like(t) theta np.zeros_like(t) # 0-5s 匀加速直线 mask1 t 5 v[mask1] 50 * t[mask1] theta[mask1] 0 # 5-10s 匀速圆周 mask2 (t 5) (t 10) v[mask2] 250 omega 250 / 1000 # 角速度 theta[mask2] omega * (t[mask2] - 5) # 10-15s 匀速直线 mask3 t 10 v[mask3] 250 theta[mask3] theta[mask2][-1] if np.any(mask2) else 0 return v, theta跑完对比频谱能看到直线运动下信号整体频偏比 Watterson 模型大匀加速状态下频谱扩展比匀速直线大匀速圆周运动伴随方向改变出现正负频偏交替。这些差异在无噪声时非常明显加噪声后需要做多次平均才能看清。验证方法上我一般会做三件事一是把f_ion设为 0只看相对运动分量确认频偏方向和大小符合 (f_c v \cos\theta/c) 的手算值二是把v设为 0退化成 Watterson 模型确认频谱和标准 Watterson 输出一致三是把sigma设为 0确认频谱退化成单根谱线。这三步走完模型基本不会有大问题。从那以后我每次搭信道仿真都强制先跑一遍退化验证再上真实参数。希望帮到你。本文还有配套的精品资源点击获取