ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

多路径效应全解析:从Matlab仿真到手机GNSS定位误差诊断

多路径效应全解析:从Matlab仿真到手机GNSS定位误差诊断 1. 多路径效应到底从哪里来一场“回声”如何干扰导航信号做高精度定位的人基本都遇到过这种怪事明明站点设在视野开阔的楼顶四周没有遮挡RTK固定解好好的某一天基线解算突然跳了几厘米过一会儿又自己回来CORS站的伪距残差出现规律性波浪怎么检查星历、有没有周跳都查不出明显毛病。后来绕着天线转了一圈发现楼下多了块施工钢板或者不远处停了辆金属车身的卡车甚至只是前几天下雨地面形成了一层均匀水膜。问题基本就出在一个词上多路径效应。卫星导航信号从两万多公里高空打下来到达地面时的功率非常低通常只有约-130 dBm比室内WiFi信号还要弱几十个数量级。接收机天线除了收到卫星直达信号还会收到经地面、墙面、水面、车顶等反射后的“迟到信号”。直达信号和反射信号混在一起进天线就像你在会议室里讲话除了原声还能听到墙壁反射的回声。人对回声还能忍受因为人耳分辨能力有限但GNSS接收机对信号到达时间的测量是按纳秒级别计算的一个迟到的“回声”压过来伪距、载波相位都会出现偏差。这就是开头那些诡异跳动的来源。理解多路径不能只把它当成一个小故障它是高精度GNSS定位中最顽固的误差源之一。它与测站周围环境强相关无法通过双频消电离层、差分方法彻底消除也很难用普通白噪声模型去描述。这篇文章不打算只讲理论我会把多路径从物理来源、数学模型、Matlab仿真到手机GNSS原始数据诊断再到工程上常用的抑制手段完整串一遍方便你需要做多路径分析或者排查定位异常时能直接照着操作。1.1 直达信号这么强为什么反射信号还能搅局很多第一次接触多路径的人会问反射信号本来就比直达信号弱功率上差了好几倍甚至十几倍为什么还能干扰接收机关键在于接收机测量的不是“信号里最强的那个分量”而是天线口面上所有信号的矢量合成结果。可以这样想直达信号是一个正弦波反射信号是另一个同频率但相位延迟的正弦波两者在天线端叠加后合成波的幅度和相位都会发生改变。相位变化多大取决于反射信号的相对幅度和相对延迟。如果反射信号与直达信号相位一致合成波变强如果相位相反合成波变弱甚至出现深衰落接收机跟踪环路会明显抽搐。载波相位测量依赖的是对合成信号相位的估计合成波相位已经偏离了直达信号相位测出来的距离自然就偏了。伪距测量同样会被带偏。伪距依靠码相关器去匹配接收到的码序列相关峰的形状本来是近似三角形的被反射信号污染后相关峰会变形、顶部变缓甚至出现多峰码环路的锁定点发生偏移伪距输出就不再是真实几何距离。由于C/A码一个码片约293米相关峰哪怕只是被“掰弯”一小部分反映到伪距上就是几米到几十米的偏差。所以在现实中伪距多路径误差往往比载波相位多路径误差大一个数量级但两者都需要认真对待。1.2 伪距和载波相位同样是多路径伤的深度不一样多路径对伪距和载波相位的影响不能用同一个尺度衡量。载波相位测量精度很高L1波长约19厘米L2波长约24.4厘米多路径引起的相位偏差通常为毫米到厘米级极限情况下最多到约四分之一波长也就是L1大约5厘米伪距就惨得多普通测地型接收机在开阔环境下伪距多路径可能只有几十厘米到几米但在反射严重的环境中到十几米也不罕见手机这类消费级设备受影响更大。观测值类型典型多路径误差范围内在原因L1/L2载波相位毫米级到约5厘米合成信号相位偏移受波长约束C/A码伪距0.5米到30米码相关峰畸变受码片长度影响P(Y)码伪距0.3米到10米码片更窄相关峰畸变稍小载噪比CN0变化3-10 dB-Hz反射信号与直达信号干涉产生起伏对比一下其他误差源电离层误差通过双频组合基本可以消掉一阶项对流层误差可以通过模型加参数估计削弱卫星钟差和接收机钟差在差分或精密单点定位里都有成熟解法。多路径不同它取决于你天线周围两米还是二十米内有什么东西同一颗卫星在不同测站上的多路径完全独立。这也是为什么高精度定位里常有这种说法多路径是“最后一块难啃的骨头”。2. 仿真之前先把数学模型吃透反射路径、幅度比和相位合成拿到Matlab源码直接跑确实很爽但如果不理解背后的数学模型改参数等于瞎猜。做多路径仿真本质上是在算三笔账反射信号比直达信号多走了多远、反射信号比直达信号弱了多少、两个信号叠加后接收机测出来的相位和伪距究竟偏了多少。这三笔账算清楚了源码里每一行都能对上号。2.1 反射路径差天线高度和仰角决定了误差的“节奏”最常见的多路径场景是地面反射假设地面是一个理想水平反射面天线相位中心到地面的垂直高度为 h卫星仰角为 el那么镜面反射点位置的几何关系很明确反射信号到达天线时相对直达信号多走的路程为$$\Delta r 2h\cdot\sin(el)$$这个公式是整个地面反射模型的地基。看出规律来了吗h 越小路径差越小仰角越低路径差也越小。路径差对应到载波相位上就是$$\Delta\varphi \frac{2\pi \Delta r}{\lambda}$$如果天线高2米卫星仰角30度路径差就是2米L1波长约0.19米一个路径差里塞了10多个波长相位差在360度范围内反复转圈。所以实际观测到的载波多路径误差看上去是一条不断振荡的曲线振荡的快慢受天线高度和卫星仰角共同控制。还有一个很实际的经验多路径误差不是高频随机噪声它的变化速度取决于卫星仰角的变化率低仰角卫星一天里仰角变化慢多路径误差往往是慢悠悠的波浪形这给滤波处理带来了很大麻烦。2.2 反射系数不是所有反射面都一样反射信号与直达信号的幅度比记为 α它由两部分决定反射面的菲涅尔反射系数以及天线在反射信号入射方向上的增益。反射系数不是一个固定数它与信号入射角、反射面材质、表面粗糙度、介电常数都有关。表里给的是工程上常用的近似参考值。反射面材质近似反射系数范围说明平静水面0.7-1.0反射最强雨后或水面反光是重灾区金属板/车顶0.8-0.9施工钢板、金属围挡极其难缠混凝土/沥青0.3-0.6常见地面影响中等裸土/草地0.2-0.5相对柔和但低仰角时仍不可忽视积雪表面0.3-0.8随含水和表面状态变化明显还有一个关键概念是镜面反射和漫反射的区别。表面凸起高度相对信号波长远小于某个阈值时信号才会形成规则的镜面反射水面、抛光金属、玻璃幕墙都是典型如果表面很粗糙信号散射到各个方向进入天线的反射能量被摊薄反而没那么危险。所以做站点环境调查时看到大面积平整表面就要格外警惕。2.3 叠加信号从相位偏移到距离误差的转换设直达信号幅度为 A反射信号幅度为 αA两者相位差为 θ合成信号可以写成$$S(t) A\cos(\omega t) \alpha A\cos(\omega t \theta)$$把合成信号展开后可以得到合成波相对直达波的相位偏移$$\psi \arctan\left(\frac{\alpha\sin\theta}{1\alpha\cos\theta}\right)$$载波相位多路径误差就是把这个相位偏移换算成距离$$\delta L \frac{\lambda}{2\pi}\psi$$当 α 接近1也就是反射信号几乎和直达信号等强时ψ 最大可以到 π/2对应 L1 载波约4.8厘米的偏差。伪距码信号的分析要复杂一些因为码相关器的鉴别器是非线性的反射信号会让相关峰变形伪距多路径误差不能用上面这个简单相位公式精确描述工程仿真里常用数值积分算相关函数或者用近似模型再乘一个经验上限。理解这一点很重要源码里算载波多路径可以直接用解析式伪距多路径想要贴近真实最好走数值仿真路线。3. Matlab源码拆解从卫星位置解算到误差曲线可视化这类“多路径效应分析”源码包网上流传的版本不少文件结构大同小异。我以一份常见结构为例拆解你自己拿到其他版本时也能快速定位关键函数。整个源码的工作思路是先设定接收机和反射面参数再计算卫星在某一时段内的位置并得到仰角、方位角序列然后代入反射几何模型生成多路径误差最后用曲线和图形展示结果。3.1 源码包里有什么目录结构与运行入口一份完整可用的源码包通常包含以下文件文件名作用main_mp_analysis.m主入口设置参数并调用各子函数satellite_position.m由星历或简化轨道参数计算卫星位置compute_azel.m由卫星位置和接收机位置计算仰角、方位角reflect_geometry.m根据反射面高度计算路径差和反射点位置mp_error_generation.m生成载波和伪距多路径误差序列plot_mp_results.m绘图输出误差曲线、半天球图等运行主函数后通常你会看到三类图多路径误差随卫星仰角的变化曲线、多路径误差随时间的变化序列、误差在方位角-仰角网格上的分布图。其中第三种最有用它直接把“哪个方向来的卫星被污染得最厉害”画在了平面上相当于给测站做了一张多路径CT。3.2 核心代码逐段解读从卫星位置到误差曲线下面这份代码是对源码核心逻辑的等价重写我用注释把关键步骤标清楚方便你做二次开发。运行环境建议R2019b以上不需要额外工具箱坐标转换部分自己写了几行简洁代码。% main_demo_mp.m % 地面反射多路径误差仿真主干脚本 clear; clc; close all; c 299792458.0; % 光速, m/s f1 1575.42e6; % GPS L1 lambda1 c / f1; % 约0.1903 m f2 1227.60e6; % GPS L2 lambda2 c / f2; % 约0.2442 m % 接收机参数 lat0 40.0; % 纬度, 度 lon0 116.0; % 经度, 度 h0 50.0; % 大地高, m h_ant 2.0; % 天线相位中心到反射面的垂直高度, m % 将大地坐标转换为地心直角坐标(系数简写适合仿真) latr lat0 * pi / 180; lonr lon0 * pi / 180; a_wgs 6378137.0; % 长半轴, m f_wgs 1.0 / 298.257223563; % 扁率 e2 f_wgs * (2 - f_wgs); N a_wgs / sqrt(1 - e2 * sin(latr)^2); rec_x (N h0) * cos(latr) * cos(lonr); rec_y (N h0) * cos(latr) * sin(lonr); rec_z (N * (1 - e2) h0) * sin(latr); % 简化卫星轨迹产生一个从5度到85度再回落的仰角序列 % 实际源码中由satellite_position.m读星历生成这里仅演示几何关系。 el_seq [5:1:85, 85:-1:5]; n length(el_seq); % 镜面反射路径差 delta_r 2 * h_ant * sind(el_seq); % 单位: m % 反射系数幅度比取0.5代表中等强度地面反射 alpha 0.5; % 载波相位多路径误差 theta 2 * pi .* delta_r / lambda1; % L1相位差 psi atan(alpha .* sin(theta) ./ (1 alpha .* cos(theta))); mp_carrier_l1 psi / (2 * pi) * lambda1; % 转为米 % 伪距多路径用相关峰变形的简化近似 mp_code alpha .* delta_r ./ (1 alpha); % 米实际需限制幅度 mp_code(mp_code 20) 20; % 工程经验上限 % 绘图 figure(Name, 多路径误差随仰角变化, Color, w); subplot(2,1,1); plot(el_seq, mp_carrier_l1*100, b-, LineWidth, 1.5); xlabel(卫星仰角 (deg)); ylabel(载波多路径误差 (cm)); title(L1载波相位多路径误差); grid on; subplot(2,1,2); plot(el_seq, mp_code, r-, LineWidth, 1.5); xlabel(卫星仰角 (deg)); ylabel(伪距多路径误差 (m)); title(伪距多路径误差(简化近似)); grid on;这份代码把核心链路压缩到了最短跑完不超过两秒。实际源码里components会改成调用子函数但数学内核就是上面这段。改参数时重点看三个量h_ant天线高度、alpha反射系数、el_seq仰角区间。把h_ant从2米改成5米你会发现振荡频率明显变快这就是路径差公式在起作用。3.3 结果图怎么读仰角误差曲线与半天球图仿真跑完后第一张仰角-误差图中能看到很典型的特征低仰角区域内误差幅度大且振荡剧烈仰角超过40度以后误差整体衰减但不会完全归零。原因是路径差公式里仰角越低反射信号越容易进入天线而且低仰角方向天线增益通常也不高反射信号相对权重更大。如果你拿到的是带方位角信息的源码建议把结果画成半天球图横轴为方位角纵轴为仰角颜色表示该方向卫星的多路径误差均值。半天球图上的“红色弧段”就是反射污染主要方向顺着那个方位出去找反射面十有八九能找到水面、金属屋顶或者玻璃幕墙。这招在野外站点验收时极其好用。4. 处理手机GNSS原始数据把多路径从观测值里抽出来很多人以为多路径分析是测量型接收机的专利其实现在Android手机也能参与。这几年“手机GNSS原始数据”越来越火小米、华为、三星等旗舰机型都开放了原始GNSS测量接口虽然手机天线性能比不上测量型天线多路径也严重得多但用来做信号质量诊断、验证算法思路完全够用。这里要区分一个概念手机能调用的原始观测值是系统底层输出的伪距、载波相位、多普勒和载噪比不是NMEA语句里的经纬度数据量要大得多。4.1 手机数据从哪来Android原始GNSS测量与日志Android 7.0之后系统提供GnssMeasurement API可以从底层拿到每一颗卫星的原始测量数据。小米手机通常在开发者选项里可以开启“GNSS日志”开启后会记录日志文件其他品牌可以用Geo RINEX Logger这类工具它能直接输出RINEX 3.03观测文件方便后续处理。日志里几个关键字段对应关系如下系统字段含义对应观测值ReceivedSvTimeNanos接收信号时间伪距的原始计数CarrierFrequencyHz信号频率L1/L5/E1/E5a等Cn0DbHz载噪比信号质量指示AccumulatedDeltaRangeMeters累积载波相位(ADR)载波相位观测值PseudorangeRateMetersPerSecond伪距率多普勒拿到这些字段后可以直接在Matlab里解析也可以转成RINEX后用成熟的观测值处理函数读取。需要注意手机载波相位不如测量型接收机稳定周跳频繁处理时要多留个心眼。4.2 MP组合把伪距多路径从观测值里剥离出来提取伪距多路径最常用的手段是MP组合Multipath Combination它利用双频伪距和载波相位的线性组合把几何距离、卫星钟差、接收机钟差、对流层甚至一阶电离层全部消掉剩下主要是伪距多路径、载波多路径和噪声。以GPS L1为例$$MP_1 P_1 - \frac{f_1^2f_2^2}{f_1^2-f_2^2} L_1 \frac{2f_2^2}{f_1^2-f_2^2} L_2$$其中 P1 是L1伪距米L1、L2是载波相位观测值换算成米。这个组合里还包含一个由模糊度和硬件延迟组成的常数项所以实际使用时要把整个序列减去它的均值剩下围绕0波动的部分就是伪距多路径加噪声。Matlab处理代码很简洁% 输入: P1, L1, L2 均为列向量, 单位: 米 % 说明: L1/L2 由 AccumulatedDeltaRangeMeters 提供 f1 1575.42e6; f2 1227.60e6; % MP组合 mp1 P1 - (f1^2 f2^2) / (f1^2 - f2^2) .* L1 ... (2 * f2^2) / (f1^2 - f2^2) .* L2; % 去掉常数偏差 mp1 mp1 - mean(mp1, omitnan); % 绘制MP1序列 plot(mp1);使用时有三个坑第一如果数据里有周跳MP序列会整体跳变一段去均值前必须先做周跳探测或分段处理否则整段偏差都会失真第二手机CN0跳变剧烈卫星刚捕获或信号衰落时伪距可能异常建议把CN0过低的数据点直接剔除第三MP组合依赖双频观测如果手机不支持第二频率可以退而求其次用单频伪距减去载波相位平滑值但那已经不是严格的MP组合了。4.3 CN0诊断信噪比是判断多路径的一把尺手机数据里另一个诊断利器是载噪比CN0。正常环境下卫星仰角越高信号越强CN0总体随仰角上升存在多路径反射时直达信号和反射信号干涉CN0会出现锯齿状起伏甚至在某几个仰角区间发生断崖式下跌这就是反射信号在某些相位条件下抵消了直达信号功率。实际操作中可以画两张图第一张是CN0随仰角的散点图每一个点代表一个观测历元第二张是MP1时间序列。两张图叠加看凡是MP1剧烈波动的地方CN0往往也在同步抖动两者对应关系非常明显。判断标准上我个人习惯把CN0低于30 dB-Hz的载波观测值标记为低质量低于25 dB-Hz直接剔除对手机来说这个阈值要更保守一些因为手机天线增益低CN0普遍比测量型接收机低不少。5. 天线、选址与数据处理三层防线的工程实践分析做得再好最终还是要回到工程问题怎么把多路径压下去。我的经验是三层防线——天线层、选址层、数据层一层比一层便宜但每一层都有独立价值。只靠其中某一层效果都有限。5.1 天线层扼流圈天线在低仰角处扼住多路径地面反射信号进入天线的路径通常是从低仰角方向来的所以天线本身要能压低低仰角增益。测量型天线里的“扼流圈天线”就是为了这个目的设计的天线底部有一圈圈同心金属环形成扼流结构对低仰角表面波产生衰减。装上扼流圈天线后低仰角方向的多路径增益能降几个dB到十几dB不等抗多路径能力显著提升。对固定站、监测站来说双频扼流圈天线属于标准配置对RTK流动站常规测量型天线已经自带较稳定的相位中心和一定的多路径抑制能力但如果你在港口、桥梁这类金属结构多的环境作业尽量选低仰角增益抑制更好的天线型号。手机天线没得选它内部是线极化或贴片天线方向图很宽多路径天生严重这也是手机定位到了厘米级之后反复出现环境误差的根源。除了天线本身还有一个低成本手段是抑径板或地面板在基准站天线下方铺一块金属圆盘人为抬高有效反射面高度使得地面反射信号被挡在天线主瓣之外。有些临时监测项目会用这个方法快速改善信号环境成本比换扼流圈天线低很多。5.2 测站选址先看反射区再埋天线天线高程和反射区距离之间存在直接几何关系对水平地面卫星仰角 el 时镜面反射点到天线底部的水平距离约等于$$d \frac{h_{ant}}{\tan(el)}$$天线高2米仰角10度的信号在约11米外的地面上发生镜面反射天线高0.5米同样仰角信号反射点只有约2.8米远。反射点离天线越远要清除的反射物范围就越大。所以天线不能架太低也不能选在大面积平整表面正中间。选址的经验规则可以整理成几条天线周围10米内尽量不要有水面、金属围挡、大面积玻璃幕墙20米内不要有大体量金属结构天线正下方地面最好是不平整的草地或碎石不要铺平整混凝土如果周围无法避开水面尽量让天线高度提高使反射路径差变大同时接受低仰角多路径增强的现实再通过数据处理压低。一天中哪些方向容易污染用半天球图看最直观。新站点埋好后先连续观测一天画出各颗卫星的MP1和残差半天球图把红色污染集中在哪个方向、哪个仰角区间记下来。之后所有数据策略都以这张图为基础。这一步看起来费时间实际上是最值钱的前期投入。5.3 数据处理层高度角定权和恒星日滤波天线和选址解决的是“物理层”多路径数据层还要处理漏进来的残余部分。最常见的做法是高度角定权卫星仰角越低观测值方差越大权重越低。经典模型是$$\sigma^2(el) a^2 \frac{b^2}{\sin^2(el)}$$a和b是经验常数大致取 a0.003-0.005mb0.003-0.01m具体值根据接收机和天线调整。这个模型把低仰角观测的权重压得很低间接削弱了低仰角多路径的影响。更高级的手段是恒星日滤波和半天球建模。GPS卫星的轨道周期约为11小时58分与恒星日高度相关也就是说同一颗卫星每天几乎在同一恒星时刻经过同一片天空如果测站周围环境没有变化多路径误差也具有重复性。把昨天提取的MP序列按23小时56分4秒平移叠加到今天的观测上就能把多路径当作系统误差校正掉一部分。这个方法的局限也很明显只要反射环境变了比如到场车辆停车位置变了、雨后水面状态变了昨天的误差模型今天就失效所以它适合长期固定站点不适合随时变化的流动站。6. 实战排查中的几个反直觉经验这部分是我在实际项目中踩过坑之后沉淀下来的体会未必写在论文里但排查问题时非常管用。如果你正在和一段诡异的定位漂移作斗争大概率能在这里找到灵感。6.1 多路径不是白噪声滤波器里最容易翻车的地方很多处理系统把观测噪声假定为高斯白噪声方差按经验设置这种近似对热噪声成立对多路径完全不成立。多路径是低频、缓慢变化、与几何构型相关的系统误差一段弧段内它可能持续十分钟朝着一个方向偏。如果滤波器把这段偏差当成噪声“平均”掉结果不是收敛到正确位置而是收敛到一个被多路径拉偏的位置并且由于滤波器过度自信固定解置信度甚至很高非常具有欺骗性。应对办法是让滤波器不要那么“相信”受污染观测。具体操作包括适当调大观测方差尤其是低仰角观测使用双频观测时优先保留载波相位质量好且CN0高的数据卡方检验和固定解Ratio检验要保留不要为了固定率好看而无限放宽条件。我在处理监测站数据时宁可让固定率从95%降到88%也不允许有系统性偏差混进最终解算。6.2 MP组合的去均值操作要小心常数项和周跳是两回事MP组合公式里的常数项包含模糊度组合和硬件延迟本来就不是零均值。所以做MP分析时去均值是常规操作但你必须清楚这个常数项只在没有周跳的连续弧段内成立。一旦发生周跳模糊度组合发生整数跳变MP序列的均值会整体移动如果不分段处理整条序列都会被这个跳变污染看起来像一条台阶状曲线。我处理手机数据时习惯先把载波相位序列做一次周跳检测最简单的方法是看相邻历元ADR的差值是否发生异常跳变。分段之后每一段各自去均值再拼接成完整的多路径时间序列。别嫌麻烦用手机数据做过一次你就会明白那个台阶跳变的坑太容易踩了。6.3 低仰角卫星砍不砍不同场景要分开决策很多人一看到多路径就想着抬高截止高度角把25度以下卫星全砍了。这个思路在静态变形监测里问题不大但动态RTK测量里风险很高。低仰角卫星对空间几何构型DOP值非常重要砍掉它们会放大其他方向误差反而降低定位精度。我的习惯是分场景处理固定站静态解算截止角可以放到15到20度配合CN0加权动态RTK和PPP-RTK不轻易砍卫星而是用CN0和高度角联合定权让低仰角信号参与解算但权重压低如果是手机数据做算法研究CN0阈值比仰角阈值更有效因为手机低仰角信号受城市楼宇反射影响太大严重时CN0已经低于可跟踪门限系统自己就丢弃了。归根结底多路径抑制的目标不是把所有低仰角卫星都赶走而是让每颗卫星的贡献与它的可信度匹配。回到开头那个站点异常的问题我现在的排查路径很固定先看半天球图找污染方向再看MP组合确认是伪距多路径还是载波相位多路径最后查现场有没有新增反射物。多路径这东西只要方向判断对了处理手段其实都是现成的最怕的是把它误判成钟差跳变或者卫星故障绕一大圈才回到原点。希望这篇文章能帮你少绕那一段路。
返回列表