简介这是一份基于Watterson模型的短波航空移动信道建模与仿真文档面向短波通信、航空通信及信道仿真相关研究人员、工程师或高年级学生用于理解飞行器相对运动与电离层反射共同作用下的信道衰落特性。文档从Watterson模型入手推导了抽头延迟线、时变频响与抽头增益函数表达式并结合多普勒效应构建适用于远距离航空移动通信的短波信道模型还针对不同飞行器参数进行差异化仿真及典型民航航迹定制化仿真。资源包含1个docx文件约526KB内容涵盖引言、模型原理、建模方法与仿真实现便于快速掌握关键公式与应用思路。已有159人学习下载适合作为课程设计、课题预研或技术文档写作的参考。1. 短波航空移动信道为什么非用 Watterson 模型不可短波航空移动通信有个很磨人的特点飞机在巡航速度下飞电离层反射路径还在不断抖动信道同时是“快衰落”和“大频移”双难场景。做这套信道的建模与仿真Watterson模型是绕不开的标准起点——ITU-R F.520把它列为窄带HF信道仿真模型结构上是一条多径抽头延迟线每条路径只需要四个参数时延、增益、多普勒频移和多普勒扩展。这篇笔记按“模型拆解 → 参数计算 → MATLAB实现 → 场景复跑 → 踩坑排查 → 统计验证”的顺序把短波航空移动信道从物理概念变成可运行的代码。照参数表和脚本走新手能搭起自己的仿真链路熟手能直接拿去评估HF数据链波形和误码率。2. Watterson 模型拆解与航空场景参数计算2.1 把 Watterson 模型拆成一张抽头延迟线图Watterson模型的数学形式并不复杂本质上就是把接收信号写成多条路径的加权叠加y(t) Σ G_i · g_i(t) · x(t - τ_i) · exp(j2πf_d,i t)下标 i 表示第 i 条路径τ_i 是这条路径的传播时延G_i 是平均幅度增益f_d,i 是多普勒频移g_i(t) 是一个单位功率的复高斯随机过程。信号先延迟 τ_i再乘上时变的复增益 G_i · g_i(t)最后搬移到 f_d,i 附近所有路径叠加就是接收信号。这里的 g_i(t) 是关键它由两个正交的高斯白噪声分别通过一个高斯形多普勒滤波器得到实部和虚部都服从零均值高斯分布所以合成包络的模服从瑞利分布相位在 0~2π 内均匀分布。这正好和电离层散射传播的统计特性对上——散射体数量多、彼此独立中心极限定理保证叠加结果趋向复高斯。模型假设每条路径的频谱都是高斯形中心在 f_d,i 附近标准差 σ_i 就是多普勒扩展。这套结构的工程价值在于把“频率选择性衰落”和“时间选择性衰落”解耦了。频率选择性由多径时延 τ_i 的不同体现时间选择性由 g_i(t) 的时变特性体现。航空场景里电离层多径数量少典型只有 2~3 条强路径带宽又窄恰好满足 Watterson 模型“窄带、路径可分离”的适用前提。做传播模型仿真时别一上来就堆十几条路径先考虑物理上是否存在这么多独立反射体。2.2 航空移动场景的三个关键参数时延、频移、谱扩展航空移动场景和地面固定站最大的区别是飞机速度带来的多普勒效应。多普勒频移用几何关系直接算f_d v · cos(θ) · f_c / cv 是飞机地速θ 是航向与来波方向的夹角f_c 是工作载频c 是光速。举例v240 米/秒θ60°f_c12 MHz算出来 f_d 240 × 0.5 × 12e6 / 3e8 4.8 Hz。这个数值在 HF 频段不算小对 3 kHz 窄带波形来说已经会影响解调性能尤其是多普勒扩展再叠加之后信道相干时间会短到几十毫秒量级。多普勒频移不是恒定值。飞机转向、爬升时 θ 变化f_d 可能几秒钟内漂移好几个赫兹天线方向图不对称时不同路径的多普勒频移也会有差异。所以仿真里我一般把 f_d 设计成按帧更新的参数每帧几十到几百毫秒重新计算一次模拟机动过程。千万不要把 f_d 写成固定常数那只能覆盖平直飞行且航向不变的理想情况。多普勒扩展 σ_i 描述的是路径谱被展宽的程度主要来源是电离层电子密度随机扰动、飞机姿态晃动和反射点不稳定。工程上典型取值地波 0.1~0.3 HzE 层反射 0.5~1 HzF 层反射 1~2 Hz强扰动条件下可以到 5 Hz 以上。σ 直接决定信道相干时间工程近似可以用 T_c ≈ 1 / (2σ) 估σ1 Hz 时相干时间约 0.5 秒σ5 Hz 时只有 0.1 秒。交织深度和均衡器更新速率都要按这个尺度设计。时延 τ_i 由传播路径长度差决定。航空 HF 场景典型取三条路径地波、E 层反射、F 层反射。地波时延接近 0E 层反射路径短时延约 0.3~0.8 msF 层反射路径长时延约 0.8~3 ms。时延差决定了频率选择性衰落的位置频域上表现为周期性的衰落谷点谷点间隔约 1/Δτ。比如 Δτ1 ms第一个深衰落谷就在 1 kHz 处这对 3 kHz 波形的子载波布点影响很大选频时要把导频放在谷点间隔之外。2.3 典型航空场景参数表可直接抄下面这张表是一组典型的短波航空移动信道参数载频 12 MHz、飞机速度 240 m/s、航向与来波夹角 60°三条路径都能给出明确取值。实际使用时把多普勒频移按自己的载频和速度等比换算时延和扩展可以在表内范围里取。路径典型来源时延 τ (ms)幅度增益 G (dB)多普勒频移 f_d (Hz)多普勒扩展 σ (Hz)说明径1地波0-10≈00.3距离近时增益更高扩展最小径2E 层反射0.50f_d 0.20.8主径时延短谱较窄径3F 层反射1.2-5f_d - 0.31.5时延长扩展大附加频移有偏差这套参数的直接含义是接收信号在时域上会看到三个能量峰值径2 是主径径1 是微弱地波径3 是延迟了 1.2 ms 的 F 层反射。如果只评估恶劣场景把 F 层增益抬到和 E 层相当或者把 σ 调到 3 Hz 以上就能模拟电离层扰动较强的状态。要注意 Watterson 模型的适用边界它要求信号带宽远小于信道相干带宽。HF 窄带数据链带宽一般是 3 kHz 以内相干带宽在几百 Hz 到几 kHz 量级满足条件。如果把几十 kHz 以上的宽带波形硬塞进来路径内频谱不再平坦模型失真这种情况得换宽带电离层信道模型。选型时先确认这一点否则后面所有仿真结果都没有参考意义。3. 用 MATLAB/Octave 实现 Watterson 信道抽头延迟线与多普勒滤波器3.1 先定采样率和帧长窄带模型不要盲目升采样率动手写代码前先把两个基础参数定下来采样率 fs 和单次处理帧长 N。HF 窄带波形带宽一般取 3 kHz采样率用 9600 Hz 就很稳妥既满足奈奎斯特又给滤波和延迟量化留了余量。没有必要为了“更精细”把采样率提到 48 kHz 或更高——Watterson 模型本身是窄带假设带宽外的频谱细节没有物理意义采样率翻倍只会让每条路径的延迟变成一长串零计算量和内存都白白浪费。帧长 N 的选择要兼顾频率分辨率和仿真时长。N 取 2 的幂便于 FFT我一般用 2^16也就是 65536 点。在 fs9600 Hz 下这对应约 6.8 秒数据频率分辨率 Δf fs/N ≈ 0.146 Hz足够分辨 0.3 Hz 量级的多普勒谱结构。如果 N 太小多普勒扩展在频域只有两三个点滤波输出形状会明显失真如果 N 太大单次 FFT 的延迟和内存占用都会上来而且一段数据内飞机航向可能已经变了f_d 的时变假设就不成立。另外要注意整个仿真链路里的信号都应该是基带复信号。发射端在 12 MHz 载频上做的调制进信道前要下变频成 I/Q 复包络信道输出再搬回载频做解调。Watterson 模型的乘法运算发生在基带这样多普勒频移才是相对基带的偏移量直接用实信号会算出双边的错误谱。3.2 频域法生成多普勒谱白噪声滤波怎么做生成 g_i(t) 的常用做法是时域滤波白噪声通过一个高斯形带通滤波器再加一个混频器搬移中心频率。但时域滤波有三个麻烦FIR 滤波器阶数不够时谱形会带肩IIR 滤波器会引入相位失真还要额外设计混频器。我更喜欢直接在频域构造高斯谱一次 FFT 乘法就得到任意形状的色噪声谱形和控制量完全一致。频域法分四步。第一步生成复白噪声 n_cplx实部和虚部各自是单位方差高斯白噪声整体功率为 1第二步对 n_cplx 做 FFT第三步在频域乘上以 f_d 为中心的高斯谱 H第四步 IFFT 回时域。这样就得到一条功率为 1、中心频率在 f_d、标准差为 σ 的高斯形随机过程。f_fft (0:N-1) * fs / N; f_fft(f_fft fs/2) f_fft(f_fft fs/2) - fs; % 映射到 [-fs/2, fs/2) H exp(-(f_fft - fd).^2 / (2 * sigma^2)); H H / sqrt(mean(abs(H).^2)); % 能量归一化 n_cplx (randn(N, 1) 1j * randn(N, 1)) / sqrt(2); g ifft(fft(n_cplx) .* H(:));这里 H 的构造有几个细节要注意。频率轴 f_fft 是按 FFT 输出顺序排列的先用 (0:N-1) 生成再把大于 fs/2 的频点映射到负半轴这样 H 和 FFT 输出自然对齐不需要做 fftshift/ifftshift 的配对。能量归一化用 sqrt(mean(|H|^2)) 而不是 max(H)保证滤波后噪声总功率保持不变——这样后级每径的增益 G_i 才能真正按 dB 值控制功率。sigma 的单位是 Hz它决定高斯谱的宽度半功率带宽约为 2.355σ工程上可以直接把 σ 当多普勒扩展用。3.3 watterson_channel 函数代码与参数说明把上面的思路封装成一个可复用的 MATLAB/Octave 函数输入是基带复信号、采样率、路径参数结构体和随机种子输出是经过信道的信号。随机种子参数很关键仿真可复现完全靠它。function y watterson_channel(x, fs, taps, n_rand) % x 输入基带复信号列向量 % fs 采样率单位 Hz % taps 路径参数结构体数组字段 % delay_s 时延单位 s % gain_db 平均幅度增益单位 dB % doppler_shift_hz 多普勒频移单位 Hz % doppler_spread_hz 多普勒扩展标准差单位 Hz % n_rand 随机种子索引保证结果可复现 rng(n_rand, twister); N length(x); y zeros(N, 1); for k 1:length(taps) % 频域高斯谱 f_fft (0:N-1) * fs / N; f_fft(f_fft fs/2) f_fft(f_fft fs/2) - fs; sigma taps(k).doppler_spread_hz; fd taps(k).doppler_shift_hz; H exp(-(f_fft - fd).^2 / (2 * sigma^2)); H H / sqrt(mean(abs(H).^2)); % 两个正交高斯白噪声源生成复包络 g(t) n_cplx (randn(N, 1) 1j * randn(N, 1)) / sqrt(2); g ifft(fft(n_cplx) .* H(:)); % 多径时延整数倍采样点 dly round(taps(k).delay_s * fs); if dly 0 x_delay [zeros(dly, 1); x(1:end-dly)]; else x_delay x; end % 该径增益并叠加 a 10^(taps(k).gain_db / 20); y y a * g .* x_delay; end end逐段解释一下。第一部分构造频域高斯谱原理就是 3.2 节说的频域滤波sigma 大于 0 是硬性要求设为 0 会在归一化时除零这是后面避坑章要重点讲的。第二部分生成复包络 g利用随机种子保证每次调用相同 n_rand 得到相同 g这对蒙特卡洛仿真里“每个信噪比点用同一组信道实现”至关重要。第三部分做时延先用 round 取整到采样点本文第 5 章会说明分数时延怎么补。第四部分把每径乘增益叠加a 用 10^(dB/20) 是因为增益定义在幅度域。注意这个函数内部没有加 AWGN。信噪比处理我习惯放在调用侧做原因很简单不同接收机的噪声基底定义不同信道模型只负责输出多径叠加信号加噪是接收链路的事。如果每径预算、滤波器能量归一化有改动信道输出功率会变在调用侧统一按实测功率加噪信噪比才可控。4. 跑通一组航空场景仿真参数表与完整脚本4.1 把航空场景翻译成参数结构体仿真第一步是把 2.3 节那张航空场景参数表翻译成 MATLAB 结构体。载频用 12 MHz飞机速度 240 m/s航向与来波夹角 60 度先算出几何多普勒频移 4.8 Hz再分别设置三条路径自己的频移和扩展。fs 9600; N 2^16; fc 12e6; v 240; theta deg2rad(60); fd v * fc * cos(theta) / 3e8; % 约 4.8 Hz taps(1).delay_s 0; taps(1).gain_db -10; taps(1).doppler_shift_hz 0; taps(1).doppler_spread_hz 0.3; taps(2).delay_s 0.5e-3; taps(2).gain_db 0; taps(2).doppler_shift_hz fd 0.2; taps(2).doppler_spread_hz 0.8; taps(3).delay_s 1.2e-3; taps(3).gain_db -5; taps(3).doppler_shift_hz fd - 0.3; taps(3).doppler_spread_hz 1.5;这里第二径和第三径的多普勒频移故意在几何频移上加了 ±0.2 Hz 和 -0.3 Hz 的偏差模拟不同反射点对飞机相对速度的差异。幅度增益按表取值径1 是弱地波径2 是主径径3 是较强的 F 层反射。如果做最恶劣场景评估把 taps(3).gain_db 改成 -2 dBdoppler_spread_hz 改成 3就能看到明显更深的衰落。4.2 用单音验证信道冲激响应参数结构体建好后先用单音信号跑一遍。单音的好处是输出频谱直接反映信道的多普勒结构多径时延体现在频域干涉纹样里。代码里我生成一个 100 Hz 的复单音过信道后画包络幅度和多普勒谱。t (0:N-1). / fs; x exp(1j * 2 * pi * 100 * t); % 100 Hz 复单音 y watterson_channel(x, fs, taps, 1); subplot(2, 1, 1); plot(t, abs(y)); xlabel(时间 (s)); ylabel(包络幅度); title(信道输出包络); subplot(2, 1, 2); Y fftshift(fft(y)); f_plot (-N/2:N/2-1) * fs / N; plot(f_plot, 10*log10(abs(Y).^2 / N)); xlim([-20 20]); xlabel(频率 (Hz)); ylabel(功率谱 (dB)); title(输出多普勒谱);跑完看两个现象。第一包络幅度随时间剧烈起伏这就是瑞利衰落径2 主径的 g(t) 在起作用第二条轨迹的包络因为叠加了三路独立随机过程看起来不像标准瑞利那么规则正常。第二多普勒谱上能量集中在 0~6 Hz 范围其中径1 的峰在 0 Hz 附近径2 的峰在 5 Hz 附近径3 的峰在 4.5 Hz 附近谱峰位置和参数表设定吻合。如果谱峰明显不对先检查 fc、v、cos(theta) 的单位和数值这类问题多半出在几何频移换算上。4.3 换成 QPSK 做误码率摸底单音验证通过后把信号换成 QPSK开始看通信性能。我这里用最简单的直接映射 QPSK比特流按两比特映射到 ±1±j过信道加高斯白噪声然后硬判决解调。这个脚本没有均衡和信道估计BER 会偏高但正好用来对比信道参数变化对链路的影响。Nsym 2^15; bits randi([0 1], Nsym, 1); symbols ((1 - 2*bits(1:2:end)) 1j * (1 - 2*bits(2:2:end))) / sqrt(2); rx watterson_channel(symbols, fs, taps, 5); snr_db 15; Ps mean(abs(rx).^2); Pn Ps / 10^(snr_db / 10); noise (randn(length(rx), 1) 1j * randn(length(rx), 1)) / sqrt(2); rx_noisy rx sqrt(Pn) * noise; demod sign(real(rx_noisy)) 1j * sign(imag(rx_noisy)); tx_ref sign(real(symbols)) 1j * sign(imag(symbols)); ber mean(demod ~ tx_ref); fprintf(SNR %d dB, BER %.4f\n, snr_db, ber);这段脚本里有个容易疏忽的地方SNR 定义在信道输出信号功率上。Ps 是在 rx 上加噪前实测的Pn Ps / 10^(snr_db/10) 保证加到 rx 上的噪声功率严格按设定 SNR 来。这样无论信道增益怎么归一化只要每次按实测功率算噪声BER 曲线就不会因为信道能量变化漂移几个 dB。扫 SNR 从 5 dB 到 25 dB每个点固定同一个随机种子跑同一组信道实现再取不同随机种子做多帧平均就能得到一条可复现的 BER 曲线。加的噪声用复数高斯白噪声实部虚部各占一半功率除以 sqrt(2) 保证总功率为 Pn。5. Watterson 仿真翻车现场5 个高频坑与排查方法5.1 随机种子没固定误码率曲线像心电图现象同一套参数、同一个 SNR 点跑两次 BER 差好几倍扫点画出的曲线上下乱跳完全看不出趋势。原因信道每条路径的 g(t) 是随机过程每次调用都重新生成单帧 BER 的统计方差本来就大如果没有固定种子等效于每次用不同的信道实现做一次试验结果当然不稳定。解决在蒙特卡洛仿真里显式控制随机种子信道函数入口传入 n_rand 参数每个 SNR 点先用一组固定种子生成信道对一组不同种子统计后求平均 BER。不要靠 rng shuffle 碰运气。5.2 多普勒扩展设成 0滤波器直接“仿真发散”现象把 doppler_spread_hz 设成 0 想模拟静态径结果函数运行报告 NaN 或者功率谱出现巨大尖峰整条 BER 曲线失效。原因sigma0 时 H 退化为一个冲激能量归一化里 mean(|H|^2) 可能为 0除法产生 NaN即便绕开归一化一个频点上的增益无限窄时域信号也会变成数值爆炸。解决任何随机径的 sigma 不要低于 0.01 Hz真正想模拟完全静态的地波单独加一条确定增益支路不走高斯随机滤波器。这是仿真发散最常见的原因先查参数表里有没有 sigma 被误设成 0。5.3 时延不是采样周期的整数倍现象仿真输出的衰落周期和理论对不上波形上出现不该有的高频抖动BER 曲线有下限地板。原因函数里时延用 round(delay_s * fs) 取整量化误差最大达到 ±0.5/fs也就是约 52 μs。对于 3 kHz 波形这个误差相当于符号周期的 15% 以上足以改变多径干涉的相位关系。解决常见做法是把时延拆成整数部分加分数部分整数部分用延迟线分数部分用 sinc 插值或者 Farrow 结构。工程上如果只关心统计性能也可以用插值后的过采样信号来跑但要付出 4 倍计算量。我一般会优先做分数时延滤波器而不是盲目提高采样率。5.4 把宽带波形硬塞进窄带模型现象带宽 20 kHz 以上的波形过信道后频谱出现奇怪的翘曲包络统计也不符合预期和实测数据差距很大。原因Watterson 模型假设每条路径内信道是频率平坦的也就是路径时延差造成的相干带宽大于信号带宽。宽带信号跨过多个衰落谷每条路径的增益随频率变化模型不再成立。解决先算相干带宽粗略用 1/最大时延差估算带宽超过相干带宽就要换宽带电离层信道模型或者把信号分成多个窄带子载波每个子载波分别过独立 Watterson 信道。这个坑在 OFDM 类波形上特别容易踩因为总带宽看着不大子载波一多实际占用带宽就上去了。5.5 盲调 SNR 导致能量不守恒现象调大某径增益后BER 曲线整体平移好几个 dB甚至出现 SNR 更高误码率反而更差的悖论。原因接收端噪声功率是按某个固定参考功率加的但信道输出功率已经因为增益调整变了信噪比实际值和标注不一致。解决把加噪环节做成“先实测信道输出功率再按目标 SNR 计算噪声功率”每次修改信道参数后重新标定一次。这个坑不只在 Watterson 仿真里出现凡是链路级仿真都要养成“实测功率加噪”的习惯不要用理论功率估算因为多径叠加、滤波归一化都可能改变实际输出功率。6. 验证与进阶用统计特性确认信道模型靠不靠谱6.1 三种验证手段包络分布、自相关、BER 对照实现完信道别急着往链路里接先做三个验证确认模型本身没有错。第一种是包络分布验证针对单条瑞利径统计输出包络的 CDF 和理论瑞利 CDF 对比。第二种是自相关验证计算 g(t) 的自相关函数确认相干时间和多普勒扩展的换算关系。第三种是 BER 对照在无多径、只保留一条静态径的条件下BER 曲线应当接近理想 AWGN 情况如果明显偏移说明增益或加噪环节有 bug。包络分布的验证代码可以复用 watterson_channel输入用全 1 信号提取单径时变增益跑 20 个独立实现采样统计tap.delay_s 0; tap.gain_db 0; tap.doppler_shift_hz 0; tap.doppler_spread_hz 1; x ones(N, 1); samples []; for k 1:20 y watterson_channel(x, fs, tap, k); samples [samples; abs(y)]; end [F, X] ecdf(samples); plot(X, F); hold on; plot(X, 1 - exp(-X.^2 / 2), r--);如果蓝色 CDF 和红色理论瑞利曲线贴合良好说明复高斯随机包络和能量归一化都没问题。如果有明显偏差检查 H 的归一化分母以及滤波输出功率是否真的为 1。这个脚本每次只跑 20 帧任何数值异常都会立刻暴露比直接看 BER 曲线定位问题快得多。6.2 进阶洛伦兹谱替换、分数时延与时变频移如果场景不是电离层反射而是对流层散射多普勒谱往往更接近洛伦兹形而不是高斯形。替换方法很简单把 H 的表达式从高斯改成洛伦兹形分母变成 sigma^2 加上频率差平方同样做能量归一化。Watterson 模型本身不限制谱形只要谱函数能积分归一化整个函数框架不用动。分数时延的工程实现我常用频域相位法在 FFT 域乘上 exp(-j2πf τ_frac)等价于一次理想分数延迟滤波。这个方案实现只有一行但要注意时延不要超过帧长。最后把多普勒频移做成按帧更新的变量模拟飞机转弯和机动帧长取 100 ms 到 1 sfd 按运动模型每帧重算再代入信道函数。我现在的习惯是每写一次信道代码先跑单音验证和包络 CDF再进链路联调。以前偷懒跳过验证直接上 QPSK 才发现多径系数符号反了来回排查浪费两天。这套流程本身不复杂但每一步都值得先确认再往前走。希望帮到你。本文还有配套的精品资源点击获取 短波航空移动信道仿真：Watterson模型参数解析与MATLAB实现
简介这是一份基于Watterson模型的短波航空移动信道建模与仿真文档面向短波通信、航空通信及信道仿真相关研究人员、工程师或高年级学生用于理解飞行器相对运动与电离层反射共同作用下的信道衰落特性。文档从Watterson模型入手推导了抽头延迟线、时变频响与抽头增益函数表达式并结合多普勒效应构建适用于远距离航空移动通信的短波信道模型还针对不同飞行器参数进行差异化仿真及典型民航航迹定制化仿真。资源包含1个docx文件约526KB内容涵盖引言、模型原理、建模方法与仿真实现便于快速掌握关键公式与应用思路。已有159人学习下载适合作为课程设计、课题预研或技术文档写作的参考。1. 短波航空移动信道为什么非用 Watterson 模型不可短波航空移动通信有个很磨人的特点飞机在巡航速度下飞电离层反射路径还在不断抖动信道同时是“快衰落”和“大频移”双难场景。做这套信道的建模与仿真Watterson模型是绕不开的标准起点——ITU-R F.520把它列为窄带HF信道仿真模型结构上是一条多径抽头延迟线每条路径只需要四个参数时延、增益、多普勒频移和多普勒扩展。这篇笔记按“模型拆解 → 参数计算 → MATLAB实现 → 场景复跑 → 踩坑排查 → 统计验证”的顺序把短波航空移动信道从物理概念变成可运行的代码。照参数表和脚本走新手能搭起自己的仿真链路熟手能直接拿去评估HF数据链波形和误码率。2. Watterson 模型拆解与航空场景参数计算2.1 把 Watterson 模型拆成一张抽头延迟线图Watterson模型的数学形式并不复杂本质上就是把接收信号写成多条路径的加权叠加y(t) Σ G_i · g_i(t) · x(t - τ_i) · exp(j2πf_d,i t)下标 i 表示第 i 条路径τ_i 是这条路径的传播时延G_i 是平均幅度增益f_d,i 是多普勒频移g_i(t) 是一个单位功率的复高斯随机过程。信号先延迟 τ_i再乘上时变的复增益 G_i · g_i(t)最后搬移到 f_d,i 附近所有路径叠加就是接收信号。这里的 g_i(t) 是关键它由两个正交的高斯白噪声分别通过一个高斯形多普勒滤波器得到实部和虚部都服从零均值高斯分布所以合成包络的模服从瑞利分布相位在 0~2π 内均匀分布。这正好和电离层散射传播的统计特性对上——散射体数量多、彼此独立中心极限定理保证叠加结果趋向复高斯。模型假设每条路径的频谱都是高斯形中心在 f_d,i 附近标准差 σ_i 就是多普勒扩展。这套结构的工程价值在于把“频率选择性衰落”和“时间选择性衰落”解耦了。频率选择性由多径时延 τ_i 的不同体现时间选择性由 g_i(t) 的时变特性体现。航空场景里电离层多径数量少典型只有 2~3 条强路径带宽又窄恰好满足 Watterson 模型“窄带、路径可分离”的适用前提。做传播模型仿真时别一上来就堆十几条路径先考虑物理上是否存在这么多独立反射体。2.2 航空移动场景的三个关键参数时延、频移、谱扩展航空移动场景和地面固定站最大的区别是飞机速度带来的多普勒效应。多普勒频移用几何关系直接算f_d v · cos(θ) · f_c / cv 是飞机地速θ 是航向与来波方向的夹角f_c 是工作载频c 是光速。举例v240 米/秒θ60°f_c12 MHz算出来 f_d 240 × 0.5 × 12e6 / 3e8 4.8 Hz。这个数值在 HF 频段不算小对 3 kHz 窄带波形来说已经会影响解调性能尤其是多普勒扩展再叠加之后信道相干时间会短到几十毫秒量级。多普勒频移不是恒定值。飞机转向、爬升时 θ 变化f_d 可能几秒钟内漂移好几个赫兹天线方向图不对称时不同路径的多普勒频移也会有差异。所以仿真里我一般把 f_d 设计成按帧更新的参数每帧几十到几百毫秒重新计算一次模拟机动过程。千万不要把 f_d 写成固定常数那只能覆盖平直飞行且航向不变的理想情况。多普勒扩展 σ_i 描述的是路径谱被展宽的程度主要来源是电离层电子密度随机扰动、飞机姿态晃动和反射点不稳定。工程上典型取值地波 0.1~0.3 HzE 层反射 0.5~1 HzF 层反射 1~2 Hz强扰动条件下可以到 5 Hz 以上。σ 直接决定信道相干时间工程近似可以用 T_c ≈ 1 / (2σ) 估σ1 Hz 时相干时间约 0.5 秒σ5 Hz 时只有 0.1 秒。交织深度和均衡器更新速率都要按这个尺度设计。时延 τ_i 由传播路径长度差决定。航空 HF 场景典型取三条路径地波、E 层反射、F 层反射。地波时延接近 0E 层反射路径短时延约 0.3~0.8 msF 层反射路径长时延约 0.8~3 ms。时延差决定了频率选择性衰落的位置频域上表现为周期性的衰落谷点谷点间隔约 1/Δτ。比如 Δτ1 ms第一个深衰落谷就在 1 kHz 处这对 3 kHz 波形的子载波布点影响很大选频时要把导频放在谷点间隔之外。2.3 典型航空场景参数表可直接抄下面这张表是一组典型的短波航空移动信道参数载频 12 MHz、飞机速度 240 m/s、航向与来波夹角 60°三条路径都能给出明确取值。实际使用时把多普勒频移按自己的载频和速度等比换算时延和扩展可以在表内范围里取。路径典型来源时延 τ (ms)幅度增益 G (dB)多普勒频移 f_d (Hz)多普勒扩展 σ (Hz)说明径1地波0-10≈00.3距离近时增益更高扩展最小径2E 层反射0.50f_d 0.20.8主径时延短谱较窄径3F 层反射1.2-5f_d - 0.31.5时延长扩展大附加频移有偏差这套参数的直接含义是接收信号在时域上会看到三个能量峰值径2 是主径径1 是微弱地波径3 是延迟了 1.2 ms 的 F 层反射。如果只评估恶劣场景把 F 层增益抬到和 E 层相当或者把 σ 调到 3 Hz 以上就能模拟电离层扰动较强的状态。要注意 Watterson 模型的适用边界它要求信号带宽远小于信道相干带宽。HF 窄带数据链带宽一般是 3 kHz 以内相干带宽在几百 Hz 到几 kHz 量级满足条件。如果把几十 kHz 以上的宽带波形硬塞进来路径内频谱不再平坦模型失真这种情况得换宽带电离层信道模型。选型时先确认这一点否则后面所有仿真结果都没有参考意义。3. 用 MATLAB/Octave 实现 Watterson 信道抽头延迟线与多普勒滤波器3.1 先定采样率和帧长窄带模型不要盲目升采样率动手写代码前先把两个基础参数定下来采样率 fs 和单次处理帧长 N。HF 窄带波形带宽一般取 3 kHz采样率用 9600 Hz 就很稳妥既满足奈奎斯特又给滤波和延迟量化留了余量。没有必要为了“更精细”把采样率提到 48 kHz 或更高——Watterson 模型本身是窄带假设带宽外的频谱细节没有物理意义采样率翻倍只会让每条路径的延迟变成一长串零计算量和内存都白白浪费。帧长 N 的选择要兼顾频率分辨率和仿真时长。N 取 2 的幂便于 FFT我一般用 2^16也就是 65536 点。在 fs9600 Hz 下这对应约 6.8 秒数据频率分辨率 Δf fs/N ≈ 0.146 Hz足够分辨 0.3 Hz 量级的多普勒谱结构。如果 N 太小多普勒扩展在频域只有两三个点滤波输出形状会明显失真如果 N 太大单次 FFT 的延迟和内存占用都会上来而且一段数据内飞机航向可能已经变了f_d 的时变假设就不成立。另外要注意整个仿真链路里的信号都应该是基带复信号。发射端在 12 MHz 载频上做的调制进信道前要下变频成 I/Q 复包络信道输出再搬回载频做解调。Watterson 模型的乘法运算发生在基带这样多普勒频移才是相对基带的偏移量直接用实信号会算出双边的错误谱。3.2 频域法生成多普勒谱白噪声滤波怎么做生成 g_i(t) 的常用做法是时域滤波白噪声通过一个高斯形带通滤波器再加一个混频器搬移中心频率。但时域滤波有三个麻烦FIR 滤波器阶数不够时谱形会带肩IIR 滤波器会引入相位失真还要额外设计混频器。我更喜欢直接在频域构造高斯谱一次 FFT 乘法就得到任意形状的色噪声谱形和控制量完全一致。频域法分四步。第一步生成复白噪声 n_cplx实部和虚部各自是单位方差高斯白噪声整体功率为 1第二步对 n_cplx 做 FFT第三步在频域乘上以 f_d 为中心的高斯谱 H第四步 IFFT 回时域。这样就得到一条功率为 1、中心频率在 f_d、标准差为 σ 的高斯形随机过程。f_fft (0:N-1) * fs / N; f_fft(f_fft fs/2) f_fft(f_fft fs/2) - fs; % 映射到 [-fs/2, fs/2) H exp(-(f_fft - fd).^2 / (2 * sigma^2)); H H / sqrt(mean(abs(H).^2)); % 能量归一化 n_cplx (randn(N, 1) 1j * randn(N, 1)) / sqrt(2); g ifft(fft(n_cplx) .* H(:));这里 H 的构造有几个细节要注意。频率轴 f_fft 是按 FFT 输出顺序排列的先用 (0:N-1) 生成再把大于 fs/2 的频点映射到负半轴这样 H 和 FFT 输出自然对齐不需要做 fftshift/ifftshift 的配对。能量归一化用 sqrt(mean(|H|^2)) 而不是 max(H)保证滤波后噪声总功率保持不变——这样后级每径的增益 G_i 才能真正按 dB 值控制功率。sigma 的单位是 Hz它决定高斯谱的宽度半功率带宽约为 2.355σ工程上可以直接把 σ 当多普勒扩展用。3.3 watterson_channel 函数代码与参数说明把上面的思路封装成一个可复用的 MATLAB/Octave 函数输入是基带复信号、采样率、路径参数结构体和随机种子输出是经过信道的信号。随机种子参数很关键仿真可复现完全靠它。function y watterson_channel(x, fs, taps, n_rand) % x 输入基带复信号列向量 % fs 采样率单位 Hz % taps 路径参数结构体数组字段 % delay_s 时延单位 s % gain_db 平均幅度增益单位 dB % doppler_shift_hz 多普勒频移单位 Hz % doppler_spread_hz 多普勒扩展标准差单位 Hz % n_rand 随机种子索引保证结果可复现 rng(n_rand, twister); N length(x); y zeros(N, 1); for k 1:length(taps) % 频域高斯谱 f_fft (0:N-1) * fs / N; f_fft(f_fft fs/2) f_fft(f_fft fs/2) - fs; sigma taps(k).doppler_spread_hz; fd taps(k).doppler_shift_hz; H exp(-(f_fft - fd).^2 / (2 * sigma^2)); H H / sqrt(mean(abs(H).^2)); % 两个正交高斯白噪声源生成复包络 g(t) n_cplx (randn(N, 1) 1j * randn(N, 1)) / sqrt(2); g ifft(fft(n_cplx) .* H(:)); % 多径时延整数倍采样点 dly round(taps(k).delay_s * fs); if dly 0 x_delay [zeros(dly, 1); x(1:end-dly)]; else x_delay x; end % 该径增益并叠加 a 10^(taps(k).gain_db / 20); y y a * g .* x_delay; end end逐段解释一下。第一部分构造频域高斯谱原理就是 3.2 节说的频域滤波sigma 大于 0 是硬性要求设为 0 会在归一化时除零这是后面避坑章要重点讲的。第二部分生成复包络 g利用随机种子保证每次调用相同 n_rand 得到相同 g这对蒙特卡洛仿真里“每个信噪比点用同一组信道实现”至关重要。第三部分做时延先用 round 取整到采样点本文第 5 章会说明分数时延怎么补。第四部分把每径乘增益叠加a 用 10^(dB/20) 是因为增益定义在幅度域。注意这个函数内部没有加 AWGN。信噪比处理我习惯放在调用侧做原因很简单不同接收机的噪声基底定义不同信道模型只负责输出多径叠加信号加噪是接收链路的事。如果每径预算、滤波器能量归一化有改动信道输出功率会变在调用侧统一按实测功率加噪信噪比才可控。4. 跑通一组航空场景仿真参数表与完整脚本4.1 把航空场景翻译成参数结构体仿真第一步是把 2.3 节那张航空场景参数表翻译成 MATLAB 结构体。载频用 12 MHz飞机速度 240 m/s航向与来波夹角 60 度先算出几何多普勒频移 4.8 Hz再分别设置三条路径自己的频移和扩展。fs 9600; N 2^16; fc 12e6; v 240; theta deg2rad(60); fd v * fc * cos(theta) / 3e8; % 约 4.8 Hz taps(1).delay_s 0; taps(1).gain_db -10; taps(1).doppler_shift_hz 0; taps(1).doppler_spread_hz 0.3; taps(2).delay_s 0.5e-3; taps(2).gain_db 0; taps(2).doppler_shift_hz fd 0.2; taps(2).doppler_spread_hz 0.8; taps(3).delay_s 1.2e-3; taps(3).gain_db -5; taps(3).doppler_shift_hz fd - 0.3; taps(3).doppler_spread_hz 1.5;这里第二径和第三径的多普勒频移故意在几何频移上加了 ±0.2 Hz 和 -0.3 Hz 的偏差模拟不同反射点对飞机相对速度的差异。幅度增益按表取值径1 是弱地波径2 是主径径3 是较强的 F 层反射。如果做最恶劣场景评估把 taps(3).gain_db 改成 -2 dBdoppler_spread_hz 改成 3就能看到明显更深的衰落。4.2 用单音验证信道冲激响应参数结构体建好后先用单音信号跑一遍。单音的好处是输出频谱直接反映信道的多普勒结构多径时延体现在频域干涉纹样里。代码里我生成一个 100 Hz 的复单音过信道后画包络幅度和多普勒谱。t (0:N-1). / fs; x exp(1j * 2 * pi * 100 * t); % 100 Hz 复单音 y watterson_channel(x, fs, taps, 1); subplot(2, 1, 1); plot(t, abs(y)); xlabel(时间 (s)); ylabel(包络幅度); title(信道输出包络); subplot(2, 1, 2); Y fftshift(fft(y)); f_plot (-N/2:N/2-1) * fs / N; plot(f_plot, 10*log10(abs(Y).^2 / N)); xlim([-20 20]); xlabel(频率 (Hz)); ylabel(功率谱 (dB)); title(输出多普勒谱);跑完看两个现象。第一包络幅度随时间剧烈起伏这就是瑞利衰落径2 主径的 g(t) 在起作用第二条轨迹的包络因为叠加了三路独立随机过程看起来不像标准瑞利那么规则正常。第二多普勒谱上能量集中在 0~6 Hz 范围其中径1 的峰在 0 Hz 附近径2 的峰在 5 Hz 附近径3 的峰在 4.5 Hz 附近谱峰位置和参数表设定吻合。如果谱峰明显不对先检查 fc、v、cos(theta) 的单位和数值这类问题多半出在几何频移换算上。4.3 换成 QPSK 做误码率摸底单音验证通过后把信号换成 QPSK开始看通信性能。我这里用最简单的直接映射 QPSK比特流按两比特映射到 ±1±j过信道加高斯白噪声然后硬判决解调。这个脚本没有均衡和信道估计BER 会偏高但正好用来对比信道参数变化对链路的影响。Nsym 2^15; bits randi([0 1], Nsym, 1); symbols ((1 - 2*bits(1:2:end)) 1j * (1 - 2*bits(2:2:end))) / sqrt(2); rx watterson_channel(symbols, fs, taps, 5); snr_db 15; Ps mean(abs(rx).^2); Pn Ps / 10^(snr_db / 10); noise (randn(length(rx), 1) 1j * randn(length(rx), 1)) / sqrt(2); rx_noisy rx sqrt(Pn) * noise; demod sign(real(rx_noisy)) 1j * sign(imag(rx_noisy)); tx_ref sign(real(symbols)) 1j * sign(imag(symbols)); ber mean(demod ~ tx_ref); fprintf(SNR %d dB, BER %.4f\n, snr_db, ber);这段脚本里有个容易疏忽的地方SNR 定义在信道输出信号功率上。Ps 是在 rx 上加噪前实测的Pn Ps / 10^(snr_db/10) 保证加到 rx 上的噪声功率严格按设定 SNR 来。这样无论信道增益怎么归一化只要每次按实测功率算噪声BER 曲线就不会因为信道能量变化漂移几个 dB。扫 SNR 从 5 dB 到 25 dB每个点固定同一个随机种子跑同一组信道实现再取不同随机种子做多帧平均就能得到一条可复现的 BER 曲线。加的噪声用复数高斯白噪声实部虚部各占一半功率除以 sqrt(2) 保证总功率为 Pn。5. Watterson 仿真翻车现场5 个高频坑与排查方法5.1 随机种子没固定误码率曲线像心电图现象同一套参数、同一个 SNR 点跑两次 BER 差好几倍扫点画出的曲线上下乱跳完全看不出趋势。原因信道每条路径的 g(t) 是随机过程每次调用都重新生成单帧 BER 的统计方差本来就大如果没有固定种子等效于每次用不同的信道实现做一次试验结果当然不稳定。解决在蒙特卡洛仿真里显式控制随机种子信道函数入口传入 n_rand 参数每个 SNR 点先用一组固定种子生成信道对一组不同种子统计后求平均 BER。不要靠 rng shuffle 碰运气。5.2 多普勒扩展设成 0滤波器直接“仿真发散”现象把 doppler_spread_hz 设成 0 想模拟静态径结果函数运行报告 NaN 或者功率谱出现巨大尖峰整条 BER 曲线失效。原因sigma0 时 H 退化为一个冲激能量归一化里 mean(|H|^2) 可能为 0除法产生 NaN即便绕开归一化一个频点上的增益无限窄时域信号也会变成数值爆炸。解决任何随机径的 sigma 不要低于 0.01 Hz真正想模拟完全静态的地波单独加一条确定增益支路不走高斯随机滤波器。这是仿真发散最常见的原因先查参数表里有没有 sigma 被误设成 0。5.3 时延不是采样周期的整数倍现象仿真输出的衰落周期和理论对不上波形上出现不该有的高频抖动BER 曲线有下限地板。原因函数里时延用 round(delay_s * fs) 取整量化误差最大达到 ±0.5/fs也就是约 52 μs。对于 3 kHz 波形这个误差相当于符号周期的 15% 以上足以改变多径干涉的相位关系。解决常见做法是把时延拆成整数部分加分数部分整数部分用延迟线分数部分用 sinc 插值或者 Farrow 结构。工程上如果只关心统计性能也可以用插值后的过采样信号来跑但要付出 4 倍计算量。我一般会优先做分数时延滤波器而不是盲目提高采样率。5.4 把宽带波形硬塞进窄带模型现象带宽 20 kHz 以上的波形过信道后频谱出现奇怪的翘曲包络统计也不符合预期和实测数据差距很大。原因Watterson 模型假设每条路径内信道是频率平坦的也就是路径时延差造成的相干带宽大于信号带宽。宽带信号跨过多个衰落谷每条路径的增益随频率变化模型不再成立。解决先算相干带宽粗略用 1/最大时延差估算带宽超过相干带宽就要换宽带电离层信道模型或者把信号分成多个窄带子载波每个子载波分别过独立 Watterson 信道。这个坑在 OFDM 类波形上特别容易踩因为总带宽看着不大子载波一多实际占用带宽就上去了。5.5 盲调 SNR 导致能量不守恒现象调大某径增益后BER 曲线整体平移好几个 dB甚至出现 SNR 更高误码率反而更差的悖论。原因接收端噪声功率是按某个固定参考功率加的但信道输出功率已经因为增益调整变了信噪比实际值和标注不一致。解决把加噪环节做成“先实测信道输出功率再按目标 SNR 计算噪声功率”每次修改信道参数后重新标定一次。这个坑不只在 Watterson 仿真里出现凡是链路级仿真都要养成“实测功率加噪”的习惯不要用理论功率估算因为多径叠加、滤波归一化都可能改变实际输出功率。6. 验证与进阶用统计特性确认信道模型靠不靠谱6.1 三种验证手段包络分布、自相关、BER 对照实现完信道别急着往链路里接先做三个验证确认模型本身没有错。第一种是包络分布验证针对单条瑞利径统计输出包络的 CDF 和理论瑞利 CDF 对比。第二种是自相关验证计算 g(t) 的自相关函数确认相干时间和多普勒扩展的换算关系。第三种是 BER 对照在无多径、只保留一条静态径的条件下BER 曲线应当接近理想 AWGN 情况如果明显偏移说明增益或加噪环节有 bug。包络分布的验证代码可以复用 watterson_channel输入用全 1 信号提取单径时变增益跑 20 个独立实现采样统计tap.delay_s 0; tap.gain_db 0; tap.doppler_shift_hz 0; tap.doppler_spread_hz 1; x ones(N, 1); samples []; for k 1:20 y watterson_channel(x, fs, tap, k); samples [samples; abs(y)]; end [F, X] ecdf(samples); plot(X, F); hold on; plot(X, 1 - exp(-X.^2 / 2), r--);如果蓝色 CDF 和红色理论瑞利曲线贴合良好说明复高斯随机包络和能量归一化都没问题。如果有明显偏差检查 H 的归一化分母以及滤波输出功率是否真的为 1。这个脚本每次只跑 20 帧任何数值异常都会立刻暴露比直接看 BER 曲线定位问题快得多。6.2 进阶洛伦兹谱替换、分数时延与时变频移如果场景不是电离层反射而是对流层散射多普勒谱往往更接近洛伦兹形而不是高斯形。替换方法很简单把 H 的表达式从高斯改成洛伦兹形分母变成 sigma^2 加上频率差平方同样做能量归一化。Watterson 模型本身不限制谱形只要谱函数能积分归一化整个函数框架不用动。分数时延的工程实现我常用频域相位法在 FFT 域乘上 exp(-j2πf τ_frac)等价于一次理想分数延迟滤波。这个方案实现只有一行但要注意时延不要超过帧长。最后把多普勒频移做成按帧更新的变量模拟飞机转弯和机动帧长取 100 ms 到 1 sfd 按运动模型每帧重算再代入信道函数。我现在的习惯是每写一次信道代码先跑单音验证和包络 CDF再进链路联调。以前偷懒跳过验证直接上 QPSK 才发现多径系数符号反了来回排查浪费两天。这套流程本身不复杂但每一步都值得先确认再往前走。希望帮到你。本文还有配套的精品资源点击获取