做GNSS数据处理的人，几乎都跟多路径效应打过照面。伪距误差从几米到几十米，载波相位也会跟着漂，而且最难缠的是——你换一台接收机、换一个环境，它的表现就完全不一样。最近在整理“GNSS多路径效应分析（含Matlab源码）”这类课题时，我把整个技术路线重新捋了一遍，从信号模型、误差机理到频谱分析和抑制算法，一条线全部走通。这篇文章就是写给正在做相关课程设计、期末project或者毕业设计的同学：我会先用大白话把多路径的来龙去脉讲清楚，再拆解这套Matlab源码背后的模块设计，然后一步步说明怎么仿真、怎么分析、怎么评价结果。文末还整理了我在实际跑数据时踩过的坑，尤其是那种“图全出来了，但结果一看就不对”的情况。
先说个生活化的类比。你在山谷里喊一声，除了直接听到的声音，还会听到山壁反射回来的回声，有时候甚至能听到好几声回声。GNSS多路径就是无线版的回声：卫星信号除了沿直线到达接收机天线，还被地面、水面、楼面、金属围栏等反射，形成延迟到达的“回声信号”。直射信号和反射信号在接收机内部叠加，导致伪距和载波相位测量值发生偏差。这个偏差不是随机白噪声，它有结构、有规律，跟卫星高度角、反射面几何、天线环境都有关系。所以要真正理解GNSS多路径效应，不能只停留在“知道有这回事”，得能建模、能仿真、能分析、能抑制。
1. 项目整体拆解：GNSS多路径效应分析到底在解决什么问题
1.1 原理铺垫：多路径信号是怎么“混”进接收机的
从信号层面看，接收机实际收到的信号可以简化为直射信号与若干反射信号的叠加：
直达信号幅度记为 A，反射信号幅度记为 a·A，其中 a 是反射衰减系数，通常在0.1到0.9之间变化。反射信号比直达信号多走了一段额外路程，对应一个时间延迟 τ，同时因为反射面的材质和入射角度不同，还会附加一个相位偏移。接收机用本地复现的伪随机码与接收信号做相关运算，得到的相关峰形状就变成了“直射相关峰”和“反射相关峰”的叠加结果。码环鉴别器根据相关峰的对称性来估计码相位，一旦相关峰被反射信号拖出“尾巴”，鉴别器的过零点就会偏移，伪距观测值因此产生偏差。
伪距多路径误差的量级相当可观。民用C/A码一个码片对应约293米的等效距离，反射信号只要延迟零点几个码片，就足以让伪距偏出数米到十几米。载波相位多路径稍微温和一点，典型条件下在L1载波波长的四分之一量级，也就是大约5厘米上下，但在反射信号组合极端的情况下也可能更大。别小看这5厘米，RTK和PPP这类高精度定位的指标都是以厘米级甚至毫米级衡量的，多路径一旦起来，浮点解直接掉进分米级甚至更难收敛。
多路径效应之所以“难搞”，主要有三个原因。第一，它不是高斯白噪声，而是有偏的、结构化的误差，常规的均值滤波和卡尔曼滤波对它效果有限。第二，它与接收机所在环境强相关，同一颗卫星的反射路径在不同站点完全不同，所以差分技术只能消除卫星钟差、对流层等公共误差，消除不掉站周围的多路径。第三，它与卫星几何构型有关，低高度角卫星的多路径尤其严重，而低高度角卫星在定位解算中又常常承担重要的几何约束作用。
1.2 源码架构：从信号生成到误差评估的五个模块
拿到“GNSS多路径效应分析”这个课题，第一反应不应该是直接翻代码，而是先把项目骨架弄清楚。这类Matlab项目通常会包含五个模块，它们的顺序恰好也是一条完整的分析流水线：
第一个模块是信号生成。根据指定卫星的PRN号生成C/A码序列，构造直射信号，再按设定的延迟码片数、反射衰减系数和相位生成反射信号，最后把两路信号叠加。第二个模块是接收机相关处理。仿真码环的本地码生成、超前/滞后相关运算、鉴别器和环路滤波，从相关结果中提取码相位误差，进而得到伪距多路径误差序列。第三个模块是误差特征提取。对伪距残差、载波相位残差或SNR序列进行处理，用双频MP组合把多路径分量从其他误差中分离出来。第四个模块是频谱分析。用FFT或Lomb-Scargle周期图分析误差序列的周期性，找出多路径对应的主频和周期。第五个模块是抑制算法与效果评估。在接收机环路、观测值和定位解算三个层面采用不同的抑制策略，对比抑制前后的RMS、峰峰值和CDF曲线。
为什么这个课题几乎总是跟Matlab绑定？因为上述流程里充满了矩阵运算和信号处理操作：相关积分本质上是序列的向量点乘，频谱分析直接调用FFT，滤波器设计一条命令就能完成，再加上Matlab自带的工具箱和绘图体系，非常适合做算法原型验证。如果你用的版本比较新，还可以留意一下Navigation Toolbox里提供的GNSS相关工具，它可以把卫星导航仿真和视觉、惯性仿真串起来。不过就多路径分析这个课题来说，核心代码并不依赖工具箱，基本的手写相关函数和FFT就能跑通。
2. 多路径效应分析的三个核心环节：信号、指纹与抑制
2.1 第一步：把直达信号和反射信号叠加成基带观测
在Matlab里做多路径仿真，不需要真的模拟1575.42MHz的L1载波中频信号，那样计算量太大而且没有必要。更聪明的做法是直接做基带仿真：生成C/A码序列，用采样率把码片离散化，然后叠加延迟的反射信号。这里有一个工程判断值得说一下：基带仿真的采样率不能太低，否则无法刻画零点几个码片的反射延迟。我习惯把采样率设在4MHz到10MHz之间，4MHz意味着每个码片约有4个采样点，对应约73米的距离分辨率，能分辨0.25码片量级的反射延迟；如果要研究更短延迟的多径，可以提高到10MHz以上。
反射信号的两个关键参数是延迟和衰减。延迟常用码片数来表示，一个C/A码片持续时间为1/1.023MHz，约977.5纳秒，对应293米等效路径差。仿真中常见的多径延迟范围是0.1到1.5个码片。举个例子，如果反射面比天线平面高5米，卫星高度角在30度附近，反射路径比直射路径多走大约2×5×sin30°=5米，对应延迟只有约16.7纳秒，折算成C/A码大约0.017个码片。实际反射场景中，由于楼面、地面、水面等反射体造成的路径延迟可能更大，所以取0.1到1.5码片能覆盖绝大多数典型情况。反射衰减系数则取决于反射面材质：水泥地、湿地面、金属围栏的反射较强，衰减系数可能到0.5到0.9；草地、灌木等漫反射明显，衰减系数往往在0.2以下。
下面是一段核心处理逻辑的示意，完整工程源码请以你拿到的项目包里的文件为准：
% 生成包含单路反射的基带信号（示意） function y = genMultipathBaseband(cfg) prn = generateCACode(cfg.sv); % 生成指定PRN的C/A码序列 t = (0:cfg.nSamples-1) / cfg.fs; % 基带采样时间轴 codeIndex = mod(floor(t * cfg.fsCode), 1023) + 1; direct = prn(codeIndex); % 直射信号 delaySamples = round(cfg.tau * cfg.fs); % 多径延迟转采样点数 refIndex = mod(floor((t - cfg.tau) * cfg.fsCode), 1023) + 1; reflected = prn(refIndex); % 延迟后的反射信号 y = cfg.amp * direct + cfg.refAmp * reflected * exp(1j * cfg.refPhase); end这里把反射信号做成了带相位偏移的复信号，是为了后续方便观察载波相位多路径效应。如果你只关心伪距多路径，也可以只用实数信号并配合衰减系数处理。代码里mod函数保证了延迟不会超出C/A码周期，但实际使用时要注意，如果延迟超过一个码片周期，反射码与直射码可能错开太远，相关峰分离过于明显，反而接近“两个卫星信号”的场景，和真实多路径行为不符。
2.2 第二步：用MP组合和SNR把多路径“指纹”提取出来
仿真跑完，得到了伪距误差时间序列，但真实场景中你拿不到“真值”，怎么从观测数据里剥离多路径分量？这是很多同学卡壳的地方。常用思路有两个。
第一个思路是构造双频MP组合（Multipath Combination）。伪距观测值里包含几何距离、卫星钟差、接收机钟差、电离层、对流层、多路径和噪声，而载波相位观测值里同样包含几何和钟差，但电离层影响与伪距相反且符号不同。通过伪距和载波相位的线性组合，可以把几何距离、卫星钟差和接收机钟差等公共项基本消掉，再通过双频组合消掉一阶电离层项，剩下的主要就是伪距多路径、接收机噪声和一个常数模糊度组合。这个组合序列的波动，就反映了多路径随时间的演化。实际代码里的MP组合系数取决于频率和伪距类型，你不需要死记公式，关键是理解它是在“做减法”，把不关心的公共误差减掉，把多路径暴露出来。
第二个思路是利用信噪比SNR或载噪比C/N0。反射信号与直射信号干涉时，合成信号的功率会随反射信号与直射信号的相位差周期性变化。随着卫星高度角变化，反射路径长度也在变化，相位差不断扫过2π的整数倍，SNR就会出现低频振荡。把SNR序列做频谱分析，主峰频率往往对应多路径反射的周期。这个思路对硬件依赖小，手机原始GNSS数据也能用，特别适合做定性判断。
拿到误差序列之后，接下来是频谱分析。这里有一个我踩过很多次的坑：多路径误差序列的频率非常低，周期通常是几分钟到几十分钟，而不是几赫兹。如果你只观测了5分钟，频率分辨率只有1/300≈0.00333Hz，对应的最小可分辨周期约300秒，基本看不出什么结构。建议静态观测时长至少30分钟以上，最好覆盖一颗卫星从低仰角升到高仰角再降下的完整过程。采样间隔1秒时，奈奎斯特频率为0.5Hz，多路径频率远低于这个值，所以采样率不是瓶颈，瓶颈是观测总时长。另外GNSS卫星高度角变化不是均匀的，时间序列严格来说不是等间隔周期信号，直接用FFT有时会糊掉主峰，此时可以试试Lomb-Scargle周期图，它对非均匀采样更友好。
2.3 第三步：窄相关、MEDLL、SNR加权哪种更值得用
多路径分析不只是为了看清误差，最终要落回到“怎么抑制它”。抑制手段分成几层。
天线层面的经典方案是扼流圈天线，通过同心圆环结构把低仰角反射信号挡在天线方向图之外。这种方案效果好，但体积大、成本高，适合基准站。接收机环路层面最常用的方案是窄相关技术，把超前/滞后相关器间隔从1个码片缩窄到0.1个码片左右。它的直观解释是“把视线收窄”，减少反射信号对鉴相器过零点的拖拽。窄相关对延迟小于1码片的短延迟多径效果很好，但对长延迟多径作用有限。再进阶一点是MEDLL（Multipath Estimating Delay Lock Loop），它同时估计多个反射信号的幅度、延迟和相位，把估计出的多径分量从相关函数中剥离，再用“干净”的相关函数解算直达信号。MEDLL效果更好，但算法复杂度明显上升，需要维护多维参数估计，实时性压力大。
观测值层面最简单的办法是SNR加权或高度角加权。多路径误差随高度角升高而快速衰减，因此在定位解算时降低低高度角观测值的权重，能在一定程度上压制多路径的影响。这个方法的缺点是它不消灭误差，只是限制误差的影响力，但在实际工程里很好用。后处理层面还可以用载波相位平滑伪距，或者对固定站建立多路径半天球图（Multipath Hemispherical Map），把每个卫星方向的历史多路径残差做成空间分布图，后续同一方向再出现卫星时直接用历史改正数去补偿。这个方法我在变形监测项目里实测过，对固定站点的多路径重复性利用非常有效。
3. 实操复盘：从数据到结果的一整套流程
3.1 数据准备与预处理：真实观测文件的三个必做步骤
想把这套多路径分析用到真实数据上，第一步是准备好观测数据。如果你手里有测量型接收机，可以输出RINEX格式的观测文件，或者直接导出伪距、载波相位、多普勒和信噪比的CSV。如果没有专业接收机，也可以用支持原始GNSS输出的Android手机，比如小米、华为、三星等机型配合GnssLogger应用，能拿到PseudoRange、CarrierPhase、CarrierFrequencyHz和Cn0DbHz等字段。手机数据做多路径分析定位是够用的，但你要有心理准备：手机天线抗多路径能力弱，伪距噪声比测量型接收机大一个量级，载波相位还经常发生周期跳变，建议只做定性分析，不要拿来做严谨的定量结论。
拿到数据后，预处理有三个必做步骤。第一是周跳探测。载波相位一旦发生周跳，MP组合序列会出现明显的阶跃跳变，整个频谱分析都会被污染。常用的探测组合有GF（Geometry-Free）组合和MW（Melbourne-Wübbena）组合，一旦检测到周跳，就把该弧段从分析中剔除。第二是粗差剔除。伪距残差突然跳几米甚至几十米的点通常不是多路径，而是接收机异常或非视距信号（NLOS），可以用中值滤波加3倍MAD（Median Absolute Deviation）的准则剔除。第三是高度角筛选。低于10度仰角的卫星信号穿过的对流层路径长、多路径极其严重，信号受遮挡和散射的概率也高，除非你专门研究低仰角多路径，否则建议只在分析时纳入，在定位解算时不纳入。
3.2 参数配置与核心函数调用：让源码按你的场景跑起来
跑通这套Matlab源码之前，先把参数配置梳理清楚。很多同学拿到代码直接点运行，图出来了觉得“好像有点东西”，但细看参数却是一团乱麻。我习惯用结构体把参数集中管理，比如这样：
cfg = struct( ... 'fs', 4e6, ... % 基带采样率，4MHz 'fsCode', 1.023e6, ... % C/A码速率，1.023MHz 'sv', 1, ... % 卫星PRN号 'tau', 0.5/1.023e6, ... % 多径延迟，0.5码片 'refAmp', 0.5, ... % 反射信号幅度，相对直达信号0.5 'refPhase', pi/4, ... % 反射信号附加相位 'dllBandwidth', 1, ... % 码环带宽，1Hz 'correlatorSpacing', 0.1, ... % 超前滞后相关器间隔，0.1码片 'duration', 3600); % 仿真时长，1小时这几个参数的作用都可以落实到物理含义上。采样率决定了延迟分辨率，码环带宽决定了环路对动态和噪声的权衡，带宽容则滤噪效果好但对动态响应慢，带太宽则噪声大。相关器间隔是窄相关设置的核心，0.1码片是经典取值。仿真时长直接影响频谱分辨率，1小时观测对应约0.000278Hz的频率分辨率，足够分辨周期为1小时的振荡；如果你的多路径周期是20分钟，那1小时观测也能看出主峰。跑代码的顺序建议是先用默认参数复现基线结果，确认没有多径时伪距误差接近零，再逐步把反射信号衰减系数调大、延迟调长，观察误差变化。
3.3 结果可视化与指标解读：判断分析有没有效
代码跑完，图会出不少，但最重要的是能读出有效信息。我认为这套分析至少要输出六类图：第一类是直射信号和叠加多径后的相关峰对比图，用来直观展示相关峰畸变；第二类是伪距多路径误差时间序列，反映误差随时间的振荡特性；第三类是SNR时间序列和对应的频谱图，用来寻找多路径主频；第四类是误差随高度角变化的散点图，用来体现“低仰角误差大”的规律；第五类是伪距误差的直方图和CDF累积分布曲线，用来做抑制前后的定量对比；第六类如果使用真实数据，建议画天空图，把每颗卫星位置的误差用颜色映射上去，能立刻看出多路径在哪个方位集中。
判断分析有没有效，有三个标志。第一，误差序列出现明显的低频平滑振荡，而不是高频白噪声；第二，频谱图上存在清晰的主峰，主峰频率对应的周期与卫星高度角变化周期、反射体几何位置能对上；第三，如果你测试了不同抑制算法，低仰角区域的误差改善要明显，高仰角区域改善较小，这是符合物理预期的。如果三条都不满足，那大概率是数据或参数有问题，赶紧回查。
4. 常见问题排查与避坑指南
4.1 现象：仿真曲线“平得像条直线”
这个问题的原因往往很朴实——反射信号根本没加进去，或者加得太弱被噪声淹没了。检查第一件事，看代码里有没有真的把refAmp和tau用到信号合成里；很多时候改参数只改了cfg结构体，但核心函数里还在用写死的旧值。第二件事，看反射衰减系数，如果设成0.01甚至0，相当于没有反射，曲线当然平。第三件事，检查相关器间隔和码环带宽，如果相关器间隔过大或者环路滤波太强，微小偏移会被平均掉。排查顺序建议是：先输出信号叠加前后的相关峰图形，用肉眼看有没有明显畸变，再检查伪距误差的统计量，最后才是调环路参数。
4.2 现象：频谱图上找不到预期峰值
频谱图空空如也，多半是观测时长不够或时间序列不均匀。前面说过，1分钟的数据只能看到周期1分钟以上的分量，多路径周期往往在几分钟以上，所以先把观测时长拉到30分钟以上。其次是时间序列不均匀的问题，卫星高度角变化率不是常数，直接用FFT经常让主峰变宽或者分裂；解决办法是先用插值把序列重采样到等间隔，或者直接用Lomb-Scargle周期图。第三是反射信号太弱，SNR振荡幅度被噪声淹没，可以提高反射衰减系数到0.3以上再试。还有一个容易被忽略的点，频谱分析的输入应该是去趋势后的残差序列，如果伪距误差序列里有一个很大的常数偏置，FFT的结果会在零频处出现巨峰，把其他周期成分都压得看不见，记得先去均值再去线性趋势。
4.3 现象：Matlab报错“未定义函数或变量”
这类报错在跑任何人给的Matlab源码时都容易出现。常见原因有几个。第一，缺少对应工具箱，比如xcorr、fft、filter这些函数分别属于Signal Processing Toolbox和Matlab基础，如果报错点落在这些函数上，大概率是工具箱没安装；到附加功能里补装就行。第二，函数名拼写或版本差异，一些较老或较新的Matlab版本对部分API有调整，比如某些导航函数在R2021a之后才进入Navigation Toolbox。第三，路径问题，源码目录没有加入Matlab搜索路径，运行前先cd到工程目录或者用addpath把子目录加进去。建议把报错信息、Matlab版本、缺失函数名一起搜索，基本都能定位。
4.4 现象：真实数据里多路径“遍地开花”
真实数据和仿真数据差异极大。真实环境不是理想镜面反射，建筑物、树木、车辆、地面坡度都会产生随机散射，多路径往往不是单一主频，而是连续谱。而且真实观测值里混着电离层残差、对流层残差、接收机钟差跳变、天线相位中心变化等问题，直接拿原始伪距做频谱分析会非常脏。我的建议是先用MP组合把公共误差消掉，再按卫星分组分析，每组只取连续无周跳弧段。同时要做好“低仰角大误差”的心理预期，10度仰角附近的伪距多路径可以到几十米，如果出现这种点先别急着当粗差剔除，它可能真是多路径的杰作。手机原始GNSS数据上更要小心，它的载波相位连续性差，用之前务必做周跳探测，否则频谱主峰全被阶跃跳变污染了。
结尾：一点个人体会
最后分享一个我自己的体会。多路径效应分析这个课题，从来不是一个“跑通源码就完事”的任务。真正有价值的是在使用这套分析流程时形成的对误差结构的直觉。当你能盯着伪距残差序列，一眼看出“这个卫星的高度角正在下降”“它刚好经过了一面金属幕墙”的时候，你对多路径的理解就超过了只会跑代码的人。所以我建议你拿到这个课题后，不要急着上复杂算法，先把仿真、频谱、仰角关系这一套跑熟，再拿真实数据去验证。多路径的克星不是某个神奇算法，而是系统化的分析和持续的环境感知，这两样东西，恰恰是Matlab这个平台能帮你养成的核心能力。
