手机上跑GNSS原始数据做定位这个事放在五年前还是件挺折腾的事现在门槛低了不少但真正把链路跑通、再把几种算法放在一起对比出说服力强的结论依然有不少坑。这个项目做的就是典型的一条龙拿手机GNSS日志解析卫星观测值结合星历自己算定位结果然后对比WLS、EKF、MHE、RTS四种算法在相同数据下的表现。我最初接这个活的时候心里想的是不就解个方程嘛真正动手才发现从日志格式的坑到滤波器发散的坑每一个都能让人卡上小半天。这篇文章把我整个实现过程、参数整定思路、踩过的坑全部整理出来适合正在做GNSS数据处理课程设计、毕业设计或者想用低成本设备入门定位算法研究的朋友。你不需要买几十万的接收机一台支持原始GNSS输出的Android手机就够了。1. 项目概览从手机GNSS日志到定位结果的一整条链路1.1 这个项目到底在做什么展开说这个项目分四层第一层是数据获取用Android手机上的GNSS Logger类应用记录原始测量值包括伪距、载波相位、多普勒频移、信号强度以及接收机时间戳。第二层是数据解析把日志文件喂给MATLAB脚本提取出每个历元、每颗可见卫星的观测数据同时读入广播星历计算卫星在信号发射时刻的位置和钟差。第三层是定位解算核心部分是四套算法WLS加权最小二乘、EKF扩展卡尔曼滤波、MHE模型预测估计、RTS平滑。它们各自输出一条定位轨迹。第四层是对比评估把四条轨迹与参考真值比较统计水平误差、垂直误差、收敛时间、计算耗时等指标。这个项目最有价值的地方在于它把接收机-数据-算法-评估的完整链条都打通了不是那种只调现成库函数的演示而是每个环节都能看到、能改、能复现的完整实现。1.2 为什么用MATLAB而不是Python或C我承认Python在GNSS领域现在很流行尤其是配合georinex、gnsspy这类库解析RINEX确实方便。但这个项目我坚持用MATLAB有几个现实原因第一MATLAB的矩阵运算和滤波器设计工具箱非常成熟。EKF、卡尔曼平滑这类递推算法用矩阵表达几乎是翻译而非实现代码可读性极高后期调参非常顺手。第二MATLAB的绘图能力在论文和报告场景下无可替代。定位轨迹对比图、误差累积分布图、卫星天空图几行代码就能输出出版级质量的图。做算法对比研究图就是脸面。第三MATLAB自带一些GNSS相关的工具箱比如Navigation Toolbox里有卫星位置计算、定位解算的参考函数。虽然我没直接用但在验证自己写的星历解析代码是否正确时这些工具箱是不错的对照基准。当然MATLAB的缺点也明显最典型的就是循环效率。但在GNSS数据处理里历元数一般也就几千个MATLAB的向量化写法足够应付完全不会成为瓶颈。2. 数据链路构建日志读取与观测值解析2.1 手机GNSS原始数据的获取与格式要拿到原始观测值手机必须支持Android 7.0及以上版本的原始GNSS测量接口Android 10以后支持得更完整。市面上大多数中高端手机都支持但有些厂商的固件会限制输出。我在测试中主要用了两台设备一台是Google Pixel系列一台是某国产旗舰前者日志格式非常规范后者偶尔会出现伪距跳变建议有条件的话优先用Pixel系或者原生Android系统的手机。采集软件我用的是Google官方的GNSS Logger它可以直接输出一个TXT格式的日志文件。这个文件里面有几种类型的行#开头的是注释行记录接收机型号、固件版本、AGPS辅助数据等信息。Raw开头的行是原始测量值字段很多包括时间戳、卫星号、载波频率、伪距、伪距rate、载波相位、CN0信噪比等。Nav开头的行是导航电文每行对应一颗卫星的星历参数。Meas开头的行是GNSS测量汇总信息。值得说明的是GNSS Logger输出的是伪距和载波相位已经经过Android框架处理的测量值不是I/Q采样级别的原始中频数据所以不需要自己处理信号捕获跟踪的问题。这对做定位算法研究反而是好事能把精力集中在更高层的处理上。如果你需要把这些日志转成标准RINEX 3.03格式可以用GNSS Logger配套的转换脚本也可以直接用我的MATLAB解析脚本一步到位。我选择自己解析而不是转RINEX原因是跳过一步转换少一个可能出错的环节而且MATLAB脚本可以直接对接后续处理省去文件格式来回折腾。2.2 关键解析细节伪距与载波相位的提取伪距观测值的核心方程是ρ r c·δt_u - c·δt_s I T ε其中r是接收机到卫星的几何距离δt_u是接收机钟差δt_s是卫星钟差I是电离层延迟T是对流层延迟。在单频场景下电离层延迟只能用Klobuchar模型粗略改正如果手机支持双频L1L5就可以用无电离层组合消除一阶电离层项。解析时最容易踩的坑是时间同步伪距对应的测量时刻是接收机接收到信号的时刻即接收机时间戳。计算卫星位置时必须用信号的发射时刻也就是接收时刻减去伪距传播时间。传播时间约为伪距除以光速大约0.07秒左右。如果直接用接收时刻去查星历算卫星位置会产生几十米的误差这个在实现时必须严格处理。载波相位的使用门槛更高一些。手机上输出的载波相位通常有整数周模糊度而且手机的天线相位中心和时钟稳定性远不如测绘级接收机直接拿来做精密单点定位不现实我在项目里主要用它来做周跳探测和观测值质量评估辅助筛选数据不直接进定位方程。2.3 星历文件与卫星位置计算广播星历包含开普勒轨道参数和摄动改正项。每个参数的含义如下sqrtA轨道长半轴的平方根e轨道偏心率i0轨道倾角omega0升交点赤经omega近地点幅角M0平近点角DeltaN平均运动角速度改正iDot、omegaDot倾角和升交点赤经的变率Cuc, Cus, Crc, Crs, Cic, Cis谐波摄动改正系数toe星历参考时刻卫星位置计算的流程是标准的解Kepler方程得到偏近点角E M e·sin(E)迭代至收敛然后计算真近点角、升交角距加上摄动改正最后旋转到地心地固坐标系。MATLAB里实现这段逻辑大概100多行核心就是查表写方程。如果你手头有Navigation Toolbox可以用gnssconstellation和ephemeris相关函数做交叉验证。我第一次实现后随机抽了几颗GPS卫星与在线SP3精密星历对比位置误差在10米量级这个量级对单点定位完全够用。2.4 观测值质量筛查与预处理数据解析出来后不能直接拿去定位必须先做质量筛查。我按以下顺序处理一是高度角筛选。设置10度或15度的截止高度角。手机在城市峡谷环境下载波相位多径严重低高度角卫星观测噪声大且多径误差明显筛掉是划算的。二是信噪比筛选。CN0低于20dB-Hz的卫星直接剔除30dB-Hz以上算健康。手机上信噪比数据经常波动做滑动平均更稳。三是伪距合理性检查。计算伪距残差如果某颗星的残差超过100米大概率是周跳或者粗差剔除。四是卫星几何检查。保留的卫星少于4颗时直接放弃该历元另外还计算DOP值HDOP大于5的历元宁可丢点也不硬算。这一套预处理做完定位结果的稳定性会有质的提升。很多人算法写了半天效果不好问题往往不在算法本身而在数据质量。3. 四种定位算法原理拆解3.1 WLS加权最小二乘最直接的批处理定位WLS是GNSS单点定位的经典解法。核心思想是在每个历元用所有可见卫星的伪距观测值最小化加权残差平方和。对于n颗卫星、4个未知数三维位置加接收机钟差的情况观测方程线性化后可以写成Δρ H·Δx其中H矩阵的每一行是接收机到卫星方向的单位向量转置最后一列为1。WLS解为Δx (H^T·W·H)^{-1}·H^T·W·Δρ权重矩阵W通常取伪距噪声方差的倒数。我用的是基于高度角的权重模型σ² a² b² / sin²(el)高度角el越低sin(el)越小噪声方差越大权重越低。这是一种经典的仰角加权模型简单且实用。MATLAB实现时要注意迭代收敛判断。因为观测方程是线化后的第一次解算用的初始位置可能是地心或者上一次的粗略解需要迭代直到Δx的范数小于某个阈值一般1毫米或者位置变化小于0.001米就收敛了。我实测下来冷启动时通常3到5次迭代就能收敛。WLS的优点是简单、稳定、没有状态模型假设每历元独立解算误差不累积。缺点也很明显没有利用载体运动的物理规律单历元解在卫星几何差或者多径严重时跳变很厉害。3.2 EKF扩展卡尔曼滤波递推式动态定位EKF的思路是把定位问题建模为标准的状态估计问题。状态向量取x [X, Y, Z, Vx, Vy, Vz, cdt_u, cdt_dot]^T也就是三维位置、三维速度、接收机钟差、钟差变化率。有的实现不带速度状态位置直接用随机游走模型也可以但带上速度状态在中低速场景下会更平滑。状态方程是最简单的常速模型CV过程噪声主要来源于载体机动加速度。测量方程是伪距的完整观测模型和WLS里的观测方程一致。EKF的递推分两步时间更新用状态转移矩阵预测状态和协方差。预测协方差的公式是P_pred F·P·F^T QQ是过程噪声协方差矩阵反映你对运动模型不确定性的估计。测量更新计算卡尔曼增益K P_pred·H^T·(H·P_pred·H^T R)^{-1}然后更新状态和协方差。MATLAB里实现EKF最怕的是协方差矩阵不正定。原因通常是数值精度问题或者Q、R矩阵设置不当。我实际测试的结果是在手机这种观测噪声较大的场景下EKF的精度相比WLS提升有限但轨迹平滑度明显更好速度输出也很稳适合运动状态分析。还有一个关键点在滤波初值。EKF对初始位置偏差非常敏感。我的做法是先跑2秒WLS用WLS的结果作为EKF的初始状态这样能避免滤波初期发散。3.3 MHE模型预测估计有限时域优化MHE全称Moving Horizon Estimation它的核心思想是把状态估计表述为固定窗口内的优化问题。每个时刻只使用最近N个历元的观测数据估计窗口内的初始状态和噪声序列取窗口末尾的状态作为当前估计。数学上MHE在每步需要求解一个带约束的最小二乘问题minimize Σ_{k0}^{N-1} ||v_k||²_{R^{-1}} ||w_k||²_{Q^{-1}} ||x₀ - x̄₀||²_{P̄₀^{-1}}其中v_k是测量噪声w_k是过程噪声x̄₀和P̄₀是先验信息。在GNSS定位里约束条件可以包括载体速度的上下限等物理约束这是MHE相比EKF的优势——它能显式处理状态约束。但MHE的计算开销很可观。窗口长度N越大优化问题的规模越大。我在实现时窗口取10个历元用MATLAB的fmincon求解单步耗时约50到100毫秒实时性勉强能达到1Hz的实时处理要求但明显比EKF慢得多。MHE的位置精度与窗口大小正相关窗口长一些、数据质量好的时候精度比EKF略好尤其是在发生短时信号遮挡时MHE因为有历史窗口约束不容易像EKF那样被单个坏历元拉偏。但如果窗口内存在较多粗差MHE同样会受影响。3.4 RTS平滑离线后处理精度上限RTS全称Rauch-Tung-Striebel平滑器。它的思路是先用标准卡尔曼滤波正向跑一遍记录每一时刻的状态估计和协方差然后从末尾往前做反向平滑。反向平滑的核心公式是x_s,k x_k K_s,k·(x_s,k1 - F_k·x_k)其中K_s,k是平滑增益是协方差矩阵的函数。RTS的意义在于它同时利用了t时刻之前和之后的所有观测信息理论上在所有线性高斯估计器中精度是最优的。在GNSS场景下RTS把WLS/EKF那种只用过去数据的限制打破了。城市环境下卫星信号经常被遮挡几分钟遮挡期间的定位误差很大但RTS在信号恢复后能回头修正遮挡时段的位置估计。在静态或低速场景RTS的精度可以比EKF提升30%以上。这个方案的代价是只能在事后处理且需要保存整个时段的滤波器状态存储量不小。手机上采集一小时数据EKF正向和RTS反向各跑一遍总的处理时间大约几秒钟完全可接受。3.5 四种算法特性横向对比我把四套算法的核心特性做了一张对比表算法计算复杂度实时性相对精度抗粗差能力利用未来数据适用场景WLS低高基准弱否实时单点定位EKF低高略优于WLS中否实时动态定位MHE高中优于EKF中强部分(窗口内)有约束的实时估计RTS中低(事后)最优中是全部高精度后处理这里的相对精度结论是基于我的实测数据不是绝对结论。不同数据质量、不同场景排序可能会有变化但RTS在最末位平滑后精度最高这个结论在绝大多数场景都成立。4. 核心实现细节与参数调优4.1 坐标系统转换ECEF、LLA与ENU定位解算输出的坐标是地心地固系ECEF下的XYZ但要评估误差或者画轨迹图最好转成经纬高LLA再投影到以参考点为原点的东北天坐标系ENU。我最常用的转换流程是ECEF转LLA用标准的迭代公式或者MATLAB的ecef2lla函数DDM格式的经纬度注意不要混结算符ENU坐标系以真值轨迹的起点或者某个参考站的坐标为原点这样误差分析能直观看出水平方向和垂直方向的偏差。这里有个容易搞混的点手机GNSS输出的是WGS84坐标系的坐标而地图应用通常用的也是WGS84但如果你做的是国内的高精度定位评估可能需要了解CGCS2000与WGS84的差异。前者和后者在大多数民用场景下差异在厘米级对单点定位精度评估来说可以忽略但如果做高程方向的精密分析还是要注意参考框架。4.2 权重矩阵的构造策略WLS和EKF的测量噪声矩阵R的构造直接影响定位精度这是我调参过程中感受最深的地方。第一版脚本我用了等权模型所有卫星的伪距噪声方差都取1米²结果在城市环境下定位误差很大尤其是低高度角卫星的观测噪声被低估了导致结果被带偏。第二版改成基于高度角的模型效果好了很多但依然存在个别卫星异常拉偏的情况。第三版在高度角基础上加入信噪比校正CN0越高权重越大。模型形式是σ² a² b²/sin²(el) × (CN0_ref/CN0)^2。这个经验模型的物理含义是信噪比直接反映接收信号的质量低信噪比观测值往往伴随显著的多径误差压低权重是非常合理的处理。对于手机上经常出现的异常伪距跳变我还在测量更新前加了一步残差卡方检验。计算每个历元的标准化残差如果某颗星的残差超过3倍标准差就把该卫星的测量噪声方差乘以一个大因子等效于降低权重而不是直接粗暴剔除。这样处理保留了卫星的几何贡献同时抑制了粗差影响。4.3 EKF的过程噪声和测量噪声整定EKF里Q矩阵和R矩阵的比值决定了滤波器的响应速度和稳定性这是整定过程中最核心的平衡。Q矩阵物理含义是你对运动模型的信任程度Q越大代表你越相信测量值、越不信运动模型滤波器响应越快但越不平滑容易受噪声影响Q越小则相反轨迹越平滑但对真实机动响应越慢。我的初始设置参考了经验值位置过程噪声取0.1 m/√s量级对应的方差速度过程噪声取0.5 (m/s)/√s钟差过程噪声取10 m/√s对应的方差钟差变化率取1 (m/s)/√s。然后用前2分钟的静态数据做灵敏度分析把位置过程噪声从小到大扫一遍找到水平误差最小的量级。如果车载运动三轴加速度突变大位置过程噪声需要上调一个数量级否则转弯时滤波会明显滞后于真实轨迹。手机步行场景则相对温和可以维持较小值。还有一个容易被忽视的细节是R矩阵的时间一致性。伪距噪声方差和接收机跟踪环路带宽有关手机与专业接收机不同室内外切换时伪距噪声特性会突然变化建议在R矩阵里加一个基于CN0的时变因子。4.4 性能评估方法算法对比不能只看轨迹看起来对不对要有量化指标。我从三个维度评估水平误差2D位置误差的均方根RMS和95%分位误差。定义是ENU坐标系下北向和东向分量的平面距离误差。这是民用导航最重要的指标。垂直误差的RMS。手机单频伪距定位的垂直误差普遍比水平误差大1.5到2倍这主要是卫星几何特性决定的。我实测WLS的垂直误差RMS约8到12米水平误差RMS约5到8米参考真值是手机自带的加速度计行人航位推算结果。计算耗时。用MATLAB的tic/toc统计每个历元的平均处理时间在需要强调实时性的场景这个指标很重要。另外我还画了CDF误差累积分布图和逐历元时间序列误差曲线。CDF图能直观看出90%的时间误差落在什么范围内比单看RMS更有说服力。关于参考真值这是个关键问题。手机定位实验很难获得厘米级的真值我用了三种参考来源静态场景直接用已知坐标点动态场景用手机自带的GNSS融合定位输出做参考虽然有误差但量级远小于单频伪距定位作为对比基准可以接受如果有条件可以在开阔场地用RTK接收机同步观测得到真正的高精度参考轨迹。5. 实测结果与避坑实录5.1 静态场景下四种算法的表现在楼顶开阔环境做了15分钟的静态测试天线位置固定手机平放视野内可见卫星数平均10颗左右GPS北斗Galileo联合定位。从这个场景拿到的数据显示WLS水平误差RMS约4.6米95%误差9.2米。由于每历元独立解算位置序列呈点状分布能明显看到多径带来的随机跳变。EKF水平误差RMS约4.1米95%误差7.8米。轨迹平滑不少但静态下仍存在低频漂移这主要来自伪距多径的慢变分量。MHE窗口N10时水平误差RMS约3.8米95%误差6.5米略微优于EKF。RTS水平误差RMS约2.9米95%误差5.4米。平滑后静态精度提升最明显低频漂移被有效抑制。需要提醒的是这些数据只是单一场景的一次测试不代表算法在所有场景下的绝对水平。做对比研究时至少要在开阔地、半遮挡环境、城市峡谷各做一组综合看结论才靠谱。5.2 动态步行场景的对比拿着手机沿校园道路正常步行约600米途经一段两侧有建筑的半遮挡区域。速度约1.2m/s。动态场景下WLS单点解跳得非常厉害半遮挡区域甚至出现了连续多个历元定位结果偏离真实路径超过30米的情况。EKF因为有速度状态约束偏离幅度被压到10米内轨迹连续但出现过冲。MHE在遮挡区域的轨迹最稳定因为窗口内的历史信息形成了约束。RTS是整个时段处理完后输出在遮挡段落的轨迹被后续信号恢复后的信息拉回视觉效果明显更好实际误差也最小。动态和静态的结论一致RTS MHE ≥ EKF WLS精确的排序依数据质量略有浮动。5.3 常见问题与排查技巧速查表我把这个项目里实际遇到的高频问题整理成了一张速查表现象可能原因排查与解决解析日志后卫星数为0日志时间戳格式不匹配检查GNSS Logger输出的是GPS周秒还是Unix时间先统一时间基准伪距全部为负值或异常大未按发射时刻计算卫星位置卫星钟差未修正核实伪距对应的时间标签补上卫星钟差改正项定位结果偏向某个方向几十米电离层/对流层改正未加或符号错误单频场景至少加Klobuchar电离层改正检查对流层改正项符号WLS迭代不收敛初始位置距离真值太远或卫星几何退化先用粗略位置如手机网络定位做初值少于4颗星时放弃该历元EKF轨迹在信号遮挡后偏移很大过程噪声Q设置过小滤波器太相信运动模型适当增大Q或在遮挡期间提高位置先验噪声滤波协方差矩阵非正定数值精度问题或Q/R设置不合理改用平方根滤波或约瑟夫形式检查矩阵对称性同一历元双重定位结果差很多卫星观测值粗差未剔除加强残差卡门检验对超3倍标准差的观测值降权RTS结果比EKF还差平滑实现中矩阵求逆数值问题或状态转移矩阵时间步长不匹配检查平滑增益公式确认状态转移矩阵用的是均匀时间间隔必要时使用稳定化的矩阵求逆方法5.4 几个值得一说的实现技巧MATLAB循环效率不高解析日志文件时逐行处理几万行数据会有点慢。我一开始用的是fgetl循环300MB的日志要跑好几秒。后来改成文本擦操作一次读入后用strsplit和cellfun做向量化处理速度提升了一个量级。解析GNSS日志这类文本处理工作尽量用MATLAB的文本批处理函数别一行行读。另一个技巧是绘图。EKF和MHE的位置输出有协方差信息可以画误差椭圆。我在轨迹图上每隔固定历元画一个95%误差椭圆能直观看出定位置信度变化。卫星信号遮挡时误差椭圆明显变大这种可视化效果在报告里非常加分。还有一个容易被忽略的问题手机GNSS日志中的GNSS接收机时钟受钟差和漂移影响明显几十分钟内可能漂移数十毫秒。当间隔较长时间再处理时务必要考虑钟差的状态转移模型是否能捕捉这种漂移否则滤波会出现系统性偏差。EKF必须把钟差变化率纳入状态向量WLS则不需要担心因为每历元重新估计钟差。6. 扩展方向与个人经验做到这一步四种算法的对比已经完整了但这个项目的可扩展性很强给你几个我觉得很有价值的方向增加RTK或者PPP处理模块。虽然手机上行RTK的整数模糊度固定很难但利用双频观测值做PPP精度可以大幅提升。加入惯导数据。手机日志里通常同时有加速度计和陀螺仪数据和GNSS做松组合在城市峡谷场景下的连续性和精度都会上一个台阶。数据驱动的方式优化噪声参数。我在调参时做了灵敏度分析能不能做个自动化的网格搜索或Bayesian优化把Q/R参数自动整定出来效果大概率优于手工调参。把算法封装成可配置的框架。目前四套算法各自独立你可以设计一个统一的接口让新算法比如粒子滤波、无迹卡尔曼滤波插进来就能跑对比实验这样后续做研究写论文非常方便。最后分享一点个人体会。这个项目做完我最大的感触不是哪个算法精度最高而是数据质量决定上限算法决定逼近上限的程度。我在实验里用同一套数据反复测试结果差异最大的不是WLS和RTS之间的算法差距而是数据预处理做得好不好的差距。把观测值质量筛查、时间对齐、卫星位置计算这些基础细节做到位哪怕只用WLS也能有不错的效果反过来再先进的滤波算法喂进去一堆粗差数据输出的也只能是看起来很平滑的错误轨迹。做GNSS定位研究别急着上复杂算法先把数据链路搞扎实再谈算法创新。这个项目的价值也正在于此——它逼着你走完整个链路每一步都能看到数据是怎么流动、怎么被处理的。如果你也在做类似的事情建议一定保持这样的习惯每处理完一个环节就把中间结果用图画出来看着数据走一遍远比闷头调算法参数要有效得多。