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

简介：这份文档面向短波通信与航空信道建模方向的研究生、科研人员及工程技术人员，聚焦远距离航空移动通信中短波信道衰落特性复杂、缺乏标准模型的问题。内容以Watterson抽头延迟线模型为基础，引入飞行器相对运动产生的多普勒频移，推导时变频响与抽头增益函数，并针对不同种类、不同参数飞行器开展差异化仿真，在航线航迹已知时还可实现特定民用航空场景的定制化信道仿真。资源包为单个docx文档，约526KB，内含引言、Watterson信道模型、短波航空移动信道建模等章节，配有公式推导与模型结构说明，便于读者直接理解建模思路并复现仿真流程。目前已有159人学习，适合需要快速掌握短波航空信道建模方法、撰写相关论文或搭建仿真验证环境的研究者参考。

1. 短波航空移动信道仿真：从 Watterson 模型到定制化航迹的完整复现路径

短波频段（3～30 MHz）靠电离层反射实现超视距通信，在卫星链路不可用时，它往往是飞行器远程指挥控制的唯一手段。但飞行器一动，收发两端的相对运动就会叠加时变多普勒频移和扩展，再叠上电离层本身的多径与时延，信道衰落特性比固定台站复杂得多。这份文档给出的方案，是在经典 Watterson 模型基础上引入飞行器运动状态参数，把多普勒效应拆成电离层分量和相对运动分量，并针对民用航空、私人飞机、无人载具三类机动频率做了差异化仿真。它适合做短波链路预算、通信体制验证、半实物仿真前端设计的工程师，也适合需要快速搭出可调信道模型的研究生。下面按“模型原理 → 参数怎么定 → 代码怎么落 → 坑在哪 → 进阶怎么用”的顺序拆开讲。

2. Watterson 模型为什么能当基线：抽头延迟线结构与三个核心假设

2.1 抽头延迟线到底在模拟什么

Watterson 模型的结构并不复杂：输入信号进入一条抽头延迟线，每个抽头对应一种电离层传播模式或路径，抽头增益函数 (G_i(t)) 对延迟后的信号做调制，模拟该路径上的衰落，最后各路相加再叠加性噪声。它的时变频响写成：

[ H(f,t)=\sum_{i=1}^{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)^2+j2\pi v_{ia}\Delta t]+C_{ib}(0)\exp[-2\pi^2\sigma_{ib}^2(\Delta t)^2+j2\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.1～1 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_i=f_{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_0+a(t)\cos\alpha(t)\cdot t ]

(v_0) 是初速度，(a(t)) 是加速度，(\alpha(t)) 是平直平面内速度与加速度方向的夹角。(\alpha=0) 时做直线运动，其余情况做圆周运动。这个式子能覆盖变速直线、匀速圆周、变速圆周等常见机动。

但飞行器受发动机功率限制，不可能一直加速。发动机功率为：

[ P=Fv=(f+ma)v=(kv^2+ma)v ]

(f=kv^2) 是空气阻力，(k) 是比例系数，(m) 是质量。当 (|v/v_{\max}|\ge0.8) 时，假设发动机功率保持最大恒定：

[ P_{\max}=kv_{\max}^3=kv^3+mav ]

推出加速度衰减规律：

[ a=k_0\left(\frac{1}{x}-x^2\right),\quad x\in[0.8,1] ]

其中 (x=|v/v_{\max}|)，(k_0=v_{\max}^2k/m)。这一步很关键：如果不加这个约束，仿真里飞行器速度会无限增长，多普勒频移直接发散，频谱图会变成一条不断外扩的斜线，完全失真。

3.3 三类机动频率的参数参考

文档按机动频率等级给出了三类典型场景的参考值：

机动频率越高，飞行器速度和加速度变化越快，额外的多普勒频移和扩展越大。仿真时假设最大加速度 80 m/s²、最大飞行速度 600 m/s，随机生成 100 s 速度变化，可以看到低机动频率下速度曲线平滑，高机动频率下速度在短时间内剧烈波动。这个表是后续所有差异化仿真的参数入口，改一个数就能切换场景。

4. 仿真实现：希尔伯特滤波器、时变多普勒与扩展的代码落地

4.1 希尔伯特滤波器与 I/Q 两路系数生成

实际短波信号载频在 3～30 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, fs=fs) 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_q

numtaps取 64～128 之间通常够用，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.1～1 Hz。如果飞行器做圆周运动，theta要按运动方向实时更新，不能取常数。

4.3 多普勒扩展的高斯滤波实现

Watterson 模型假设多普勒扩展功率谱服从高斯分布。高斯白噪声功率谱均匀，用高斯滤波器滤一下就能得到高斯功率谱噪声序列，再与输入信号相乘，时域相乘对应频域卷积，实现频谱扩展。噪声采样率太高会导致高斯滤波器阶数过大，常见做法是用低采样率噪声序列插值滤波到系统采样率：

from scipy.signal import gaussian, lfilter from scipy.interpolate import interp1d def apply_doppler_spread(signal, fs, sigma, fs_noise=1000): """ 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)), std=fs_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, kind='linear', fill_value='extrapolate') noise_sys = interp(t_sys) return signal * noise_sys

sigma按路径分别取 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, dtype=complex) for tau, sigma in zip(path_delays, sigmas): delay_samples = int(tau * fs) delayed = np.concatenate([np.zeros(delay_samples, dtype=complex), 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至少取 64，f0不超过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 以上，中机动跑 10～20 s，高机动跑 1 s 左右。或者把多个机动片段拼接起来，模拟连续机动。

6. 定制化航迹仿真：已知航线时怎么把信道模型钉到具体场景

民用航空的航线航迹往往是确定的，这时候可以把运动参数直接写死，做定制化仿真。文档给了一个典型三段式航迹：0～5 s 匀加速直线，初速度 0，加速度 50 m/s²；5～10 s 匀速圆周，半径 1000 m，速度 250 m/s，加速度大小 (v^2/R) 方向时变；10～15 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，确认频谱退化成单根谱线。这三步走完，模型基本不会有大问题。

从那以后我每次搭信道仿真，都强制先跑一遍退化验证，再上真实参数。希望帮到你。

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