简介一篇围绕GPS卫星观测数据仿真的学术论文PDF原载《北京测绘》2018年第32卷第6期作者崔立鲁、熊旭等人从卫星数据特性出发系统阐述GPS观测数据仿真的基本原理面向GNSS导航定位、系统开发与算法验证领域的研究人员及测绘专业学习者。论文内容涵盖卫星轨道参数构建、信号传播时延计算、载波相位模拟等关键环节并利用仿真观测数据解算测站点坐标与真实伪距值及坐标值进行对比验证了仿真原理的正确性与可靠性为后续GPS数据处理软件检验提供了方法参考。资源为单个PDF文档整体大小约1.27MB学术文献格式完整包含摘要、关键词、引言、仿真模型、数据比较及结论等要素。当前已有103人学习浏览适合需要系统了解GPS观测数据仿真流程、开展算法验证工作或作为教学培训参考资料的读者使用。1. GPS卫星数据仿真一份技术文档背后的三个现实问题GPS卫星数据仿真这个词核心是把真实卫星的星历、钟差、信号传播延迟变成一组可重复的伪距和载波相位观测值供接收机算法在实验室里做测试。之前有个项目接收机算法改了一版去外场跑一次要一天而且天气和卫星几何没法复现后来用仿真数据三分钟就能把静态定位和动态轨迹全跑一遍。这份技术文档讲的就是这条流水线从广播星历到观测文件再到算法验证。适合正在做接收机基带、导航解算、组合导航或者自动化测试的工程师也适合想给课程实验搭一套真实数据的同学。2. 先搞清楚仿什么从广播星历到伪距/载波相位的完整链路GPS卫星数据仿真第一步不是写代码是选星历输入。常见做法是选广播星历理由是文件小、格式公开、字段齐全。广播星历每两个小时更新一次轨道精度在米级对接收机算法验证够用了。如果是做厘米级定位的仿真比如验证载波相位差分解算那就要用精密星历采样间隔一般15分钟或30分钟轨道精度能到厘米级。我一般会先在项目目录里建两个目录brdc和precise把两种输入分开后面换场景方便。2.1 广播星历与精密星历仿真输入选哪个广播星历里真正用的参数是开普勒轨道根数加摄动修正。RINEX导航文件里把这些参数按行排布解析时注意单位必须转成弧度或秒很多仿真做出来位置误差几十公里就是单位没转对。下面这张表是我每次写解析器都要对照的清单参数含义单位toe星历参考时刻秒周内秒sqrt(A)轨道长半轴平方根根号米e偏心率无M0参考时刻平近点角弧度omega近地点幅角弧度Omega0参考时刻升交点赤经弧度i0参考时刻轨道倾角弧度OmegaDot升交点赤经变化率弧度/秒IDOT轨道倾角变化率弧度/秒C_uc, C_us升交距角余弦/正弦摄动修正弧度C_rc, C_rs轨道半径余弦/正弦摄动修正米C_ic, C_is轨道倾角余弦/正弦摄动修正弧度选广播星历还是精密星历取决于你要验证什么。做单点定位、接收机自主算法、多径抑制这类测试广播星历就够了因为广播星历本身包含轨道误差和钟差刚好能验证算法在真实精度下的表现。做精密单点定位、RTK模糊度固定、或者评估电离层模型时我会换精密星历。精密星历通常用SP3格式存储每15分钟给一个卫星坐标需要插值到任意时刻。插值方法常用拉格朗日插值阶数取9到12插值边缘误差会明显变大所以取时间窗口时要留余量。还有一种做法是直接使用精密星历附带的卫星钟差它的时间分辨率比广播星历高可以到30秒甚至5秒对高动态接收机仿真很重要。实际项目里我不会在同一个仿真器里混用两种星历。广播星历和精密星历算出来的卫星位置差异可能达到米级如果基准站和流动站一个用广播星历一个用精密星历差分残差里会引入系统性偏差。这个坑在后面避坑章节还会提到。2.2 三个观测量模型伪距、载波相位、多普勒的数学表达式有了卫星位置再往下就是生成接收机看到的观测量。核心是三个伪距、载波相位、多普勒。伪距的模型可以写成这样P ρ c·(dt_r − dt_s) I T MP ερ是卫星到接收机的几何距离dt_r是接收机钟差dt_s是卫星钟差I是电离层延迟T是对流层延迟MP是多径误差ε是测量噪声。载波相位的模型要特别注意电离层符号φ ρ c·(dt_r − dt_s) − I T λ·N ελ是载波波长N是整周模糊度。伪距和载波相位的电离层项符号相反这是仿真里最容易漏掉的一条。如果两个都加成正号双频无电离层组合会算出异常大的电离层残差算法再聪明也救不回来。整周模糊度在静态仿真里通常设成一个固定整数动态环境下如果出现周跳需要单独生成周跳标志而不是简单地把模糊度重置。多普勒观测量在动态场景里价值最高D −(f/c)·dρ/dt f·(d dt_r/dt − d dt_s/dt) εf是载波频率c是光速dρ/dt是几何距离变化率。很多入门仿真直接用位置差分来近似距离变化率导致输出频率抖动。正确做法是从轨道速度计算卫星速度再结合接收机速度和单位视线向量点乘。对于接收机钟漂如果仿真里不需要模拟高稳定度晶振可以设为零如果要做组合导航建议加上随机游走钟漂否则多普勒残差太干净卡尔曼滤波器容易过度信任观测值。2.3 最小仿真链路从星历输入到观测输出的标准步骤一个可复现的最小链路分五步第一步读星历第二步算卫星位置和速度第三步给定接收机位置和时间第四步计算几何距离并注入各项误差第五步写RINEX观测文件。这五步串起来就是一套仿真器。不要一开始就追求实时星座和射频信号先把观测文件跑通用软件接收机解算位置能对上再考虑加场景和误差。我这边给一个判断链路是否完整的检查点导航文件里至少要包含当前时刻的星历观测文件里每颗卫星的伪距、载波相位、多普勒、信噪比四个量都在并且观测时刻能对应上。只要有这些单点定位就能跑起来。下面是链路各环节的输入输出清单步骤输入输出解析星历RINEX导航文件卫星轨道参数、钟差参数卫星位置解算轨道参数、观测时刻ECEF坐标、速度观测量生成接收机坐标、卫星坐标、误差参数伪距、载波相位、多普勒文件输出观测量、时间、系统信息RINEX观测文件闭环验证观测文件、导航文件接收机位置、残差时间系统是这条链路里最麻烦的地方。RINEX导航文件里的toe通常用GPS周内秒表示观测文件的TIME OF FIRST OBS一般用UTC表示。如果仿真器里直接拿UTC时间算卫星位置结果会差好几个微秒乘上光速就是几百米。我习惯在代码开头定义一个时间转换函数统一转成GPST只在写文件时转回UTC。这样既方便计算又不会在输出格式上出错。3. 搭一套能跑的仿真环境卫星位置解算与RINEX文件生成3.1 卫星位置解算用广播星历参数算ECEF坐标含Python代码GPS卫星数据仿真里最核心的算法是卫星位置解算。我用广播星历参数按经典算法算代码不长但每一步都有讲究。先看代码再解释为什么这样写。import math def satellite_position(ep, t): 根据广播星历参数 ep 计算 t 时刻卫星在地心地固系中的位置。 ep 是从 RINEX 导航文件解析出来的字典键与 2.1 节表格一致。 t 是接收机观测时刻单位是秒GPS 周内秒。 # 时间差从星历参考时刻 toe 到当前时刻 t tk t - ep[toe] if tk 302400: tk - 604800 if tk -302400: tk 604800 # 平近点角 M ep[M0] math.sqrt(ep[mu] / ep[sqrtA]**6) * tk # 偏近点角迭代求解 E M for _ in range(10): E M ep[e] * math.sin(E) # 真近点角和升交距角 v math.atan2(math.sqrt(1 - ep[e]**2) * math.sin(E), math.cos(E) - ep[e]) phi v ep[omega] # 摄动修正 delta_u ep[Cuc] * math.cos(2*phi) ep[Cus] * math.sin(2*phi) delta_r ep[Crc] * math.cos(2*phi) ep[Crs] * math.sin(2*phi) delta_i ep[Cic] * math.cos(2*phi) ep[Cis] * math.sin(2*phi) u phi delta_u r ep[sqrtA]**2 * (1 - ep[e] * math.cos(E)) delta_r i ep[i0] delta_i ep[IDOT] * tk # 升交点经度考虑地球自转 omega ep[Omega0] (ep[OmegaDot] - 7.2921151467e-5) * tk \ - 7.2921151467e-5 * ep[toe] # 轨道平面坐标 x_orb r * math.cos(u) y_orb r * math.sin(u) # 转到 ECEF x x_orb * math.cos(omega) - y_orb * math.cos(i) * math.sin(omega) y x_orb * math.sin(omega) y_orb * math.cos(i) * math.cos(omega) z y_orb * math.sin(i) return x, y, z这个函数有几个必须注意的地方。第一tk 要处理跨周问题超过半周要加减一周。第二偏近点角迭代必须初始化成平近点角不能从0开始。第三升交点经度要考虑地球自转项否则到中午和午夜的位置偏差会差一个方向。第四RINEX导航文件里 sqrt(A) 的单位是根号米M0、omega、Omega0 都是弧度如果按角度传进来位置直接偏出去几百公里。我还会把 mu、地球自转速率这些常量集中放一个配置里避免散落在代码里。3.2 生成RINEX 3.04观测文件头与数据块的必填字段卫星位置算完下一步是把观测量写到RINEX观测文件里。RINEX 3.04版本和2.11相比文件命名规则变了观测码也变成两个字符。我一般用3.04兼容性好很多新算法直接读3.04。文件头至少要有下面几个字段字段说明RINEX VERSION / TYPE必须是 3.04观测文件为 OTIME OF FIRST OBS第一个观测时刻UTCINTERVAL采样间隔单位秒# / TYPES OF OBSERV观测量类型数量及类型码END OF HEADER头文件结束标记因为仿真数据是自己生成的最容易出问题的是观测量类型码。RINEX 3.04里伪距是C1C、C1W载波相位是L1C、L1W多普勒是D1C、D1W信噪比是S1C、S1W。不要把2.11里的C1、L1、D1混进去否则很多解析器会直接报错。观测数据块按卫星编号顺序排列每行要写对应的伪距、载波相位带LLI标志位、多普勒和信号强度。载波相位如果带半周期修正要用0和1标记LLI位。实际仿真中我建议先不给整周模糊度赋值写成0后面做差分验证时再单独设置。3.3 用单点定位做闭环验证误差从哪看写好RINEX文件后一定要做闭环验证。常见做法是用开源的单点定位解算器输入导航文件和观测文件输出接收机位置。把解算位置和仿真时设定的接收机真值比较看误差在哪个量级。如果用的是理想观测值没有任何误差注入单点定位结果应该在毫米级或更小。这时候误差来源主要是星历本身广播星历的轨道误差会映射到定位结果里这是正常的。如果解算出来偏差几十米甚至几百公里先检查卫星位置代码用一段已知的星历数据和某个标准算法交叉对比而不是直接怀疑观测量。我验证时的习惯是把同一份RINEX导航文件喂给两个独立实现的解析器比较同一颗卫星在同一个时刻的XYZ坐标。最大差异小于1米说明卫星位置没写错差异在厘米级说明解析流程基本一致。然后再去看观测文件能不能被解算器正常解析。上面这些步骤全部走通这套仿真环境才算真的立住了。另一个值得做的验证是给接收机设置一个已知轨迹仿真完解算出来把逐历元位置和真值画在一起看残差是否有周期性。周期性通常来自卫星几何重复而不是算法错误判断的时候要区分开。4. 参数决定可信度误差注入与场景化仿真设置4.1 卫星钟差、群延迟与码偏差漏一个就是米级误差仿真里最容易被忽略的是卫星钟差。真实的GPS卫星钟差可以用广播星历里的af0、af1、af2三个多项式系数表达误差在纳秒量级乘上光速后就是分米到米级。很多仿真直接忽略卫星钟差伪距只给几何距离解算结果会带上一个随时间漂移的公共误差。看起来位置还能解但残差曲线会暴露问题。群延迟TGD也要单独建模。TGD是卫星硬件对不同频率信号造成的延迟差会在不同观测码之间产生系统偏差。如果仿真同时输出C1C和C1W两种伪距或者做双频组合TGD不设置组合结果里会出现常值偏差。我一般把TGD从导航文件中读出来在生成伪距时单独加上同时把对应的群延迟差分项也写进观测文件保证后处理能恢复。另外还有一个容易被忽略的是接收机硬件延迟它在不同频点之间也会有偏差。真实接收机里这部分偏差被时钟吸收但仿真里如果接收机钟差模型太简单硬件延迟就会以常值残差形式出现在定位结果里。4.2 电离层与对流层延迟模型选择与适用边界电离层延迟在不同频点上不同最粗糙的做法是用广播电离层模型系数输入8个参数输出斜向延迟改正值。这个模型在中纬度地区白天能消掉50%到60%的电离层误差但在低纬度和太阳活动高峰期会明显不足。对仿真来说如果只是验证单点定位用这个模型就够了如果要验证双频消电离层算法就直接在伪距上加相反的符号项让算法自己去消。对流层延迟分干项和湿项。干项可以用气压和温度经验公式计算湿项依赖水汽分布很难精确建模。在仿真里我习惯把对流层延迟拆成天顶干延迟、天顶湿延迟和映射函数三部分分别配置。静态场景用标准大气模型动态场景最好让高度角参与映射计算因为低仰角卫星的对流层延迟会被放大很多倍。映射函数如果直接用1/sin(elevation)在5度以下会趋向无穷大仿真时要给仰角设置下限否则会生成极其离谱的观测量。4.3 多径和遮挡场景给观测值加多少噪声算合理市区和隧道场景的多径是最难仿真的。最简单的做法是给伪距加上慢变的正弦误差幅度按环境区分开阔地0.1米到0.3米城市峡谷0.5米到1米树下信号遮挡场景可以到几米。载波相位的多径要小得多一般只有四分之一波长的量级。我通常先在理想观测值上跑一遍算法再逐步加多径和噪声对比两次结果才能看出算法在恶劣场景下损失了多少性能。遮挡仿真的关键是卫星仰角阈值。开阔场景阈值设10度城市峡谷设25度到30度楼宇遮挡还可以按方位角做扇形掩膜。不要用同一个阈值跑所有场景否则隧道定位测试完全没有参考价值。还有一点多径误差不是白噪声它有相关性如果直接加高斯白噪声很多抗多径算法会测不出来。至少要生成几个主反射路径让误差存在明显的时变趋势。4.4 误差预算把噪声参数组织成配置文件仿真的参数多了之后最好做一个误差预算表把每项误差的标准差和相关性写清楚。下面是我在一个模拟项目里用过的参数模板误差项分布类型标准差/范围备注卫星钟差多项式随机游走等效距离1米需要从星历读取接收机钟差随机游走等效距离10米动态场景按晶振模型电离层残差与高度角相关2米到10米广播模型只能消一部分对流层湿延迟高斯分布0.1米到0.3米低仰角按映射函数放大多径有色噪声0.3米到1米不要用白噪声测量噪声高斯分布0.1米到0.5米伪距比载波相位大很多把这些参数单独放一个JSON或配置文件里场景切换时不用改代码只改配置。这个习惯能省下很多返工时间。5. 避坑手册GPS卫星数据仿真的5个高频问题和排查路径5.1 卫星仰角全为负坐标系旋转顺序错了现象解算卫星位置后发现所有卫星仰角都在地平线以下接收机一颗星都看不到。 原因把ECEF坐标转换到接收机站心坐标系时旋转矩阵的先后顺序写反了。常见做法是先转经度再转纬度有人把两个旋转矩阵相乘时顺序弄反。 解决检查旋转矩阵确认是先绕Z轴旋转经度再绕Y轴旋转负纬度。或者用标准的地固系到站心系的变换矩阵直接替换。我给自己的排查顺序是先算地心距离和经纬度确认卫星位置本身在轨道上再看仰角。5.2 伪距输出不连续接收机时间基准没对齐现象仿真得到的伪距在同一颗卫星上相邻历元间出现几十公里跳变看起来像整周跳变。 原因接收机钟差在仿真里可能被设成随时间变化的量但没有和观测时刻对齐。常见的是在生成伪距时用了不同的时间系统一会儿GPST一会儿UTC导致计算几何距离时卫星位置和接收机位置不在同一个时刻。 解决统一用GPST。如果需要UTC输出在最后转换时加固定偏移。我还会在代码里单独定义一个时间对象把周内秒、GPS周和UTC偏移放在一起避免手动换算。5.3 RINEX文件被解算软件拒绝头信息与观测量不匹配现象自己的RINEX文件用文本编辑器看没问题但喂给解算程序直接报错或者只解出几颗星。 原因头文件里写的观测量数量比实际数据块少或者观测类型码与RINEX 3.04标准不一致。最常见的是在头文件写了C1C、L1C但数据块里还混着D1C、S1C之外的旧格式。 解决用解析器输出错误行号或者自己写一个RINEX校验函数统计每个历元的卫星数和观测量个数。我常写的校验逻辑是遍历文件检查每个数据块的卫星PRN是否在导航文件里存在不存在就跳过。5.4 多普勒频移和运动速度对不上动态场景运动模型不一致现象接收机匀速运动时某颗卫星的多普勒变化率和理论值差很多速度解算异常。 原因仿真时接收机位置更新率低比如每100毫秒更新一次但多普勒是按一个历元50毫秒算的时间基准不一致。或者是接收机速度方向用错了导致视线方向投影算错。 解决多普勒计算时统一使用接收机在相邻两个历元间的平均速度不要用瞬时速度近似。如果仿真支持动态轨迹轨迹生成和高频观测值生成使用同一个运动模型避免位置、速度、加速度互相矛盾。5.5 仿真结果好得像假数据忘记注入偏差与噪声现象单点定位误差一直在毫米级换了场景也一样看起来算法完美。 原因没有注入任何误差项或者只加了高斯白噪声忽略了钟差、电离层、对流层这些系统性偏差。真实接收机输出从不会有那么干净的残差。 解决至少加入卫星钟差、接收机钟差、电离层延迟和一组有色噪声。我常用的做法是先生成一组理想观测值再按误差预算分开加每个误差项单独存一份这样定位算法调试时可以根据误差类型回溯。6. 进阶把静态仿真改成动态差分验证的场景设计静态仿真跑通之后真正有价值的是把同一套仿真数据用在差分定位验证上。差分定位需要两台接收机同时观测同一颗卫星仿真器可以很容易做到这一点给定同一个星座分别在基准站和流动站位置生成观测值。关键是不要让两路观测值共用同一套噪声否则做站间单差后误差会被完全抵消算法收敛速度看起来不真实。我常用的做法是先确定基准站坐标和流动站运动轨迹对每个历元分别计算两颗接收机的伪距和载波相位。卫星钟差和电离层误差在两台接收机之间有强相关性单差后能消掉大部分但残余项要按基线长度重新生成不能直接设成零。这样生成的双差观测值才能真正检验RTK整周模糊度固定算法。载波相位整周模糊度在基准站和流动站之间没有固定约束仿真时可以随机给一组整数但要在文件中记录真值方便算法比对。动态场景里还有一个习惯值得保留每次仿真都记录一份真值文件内容包括接收机真实位置、速度、卫星位置、钟差和所有误差项。后处理算法如果解算失败拿真值做一个逐历元残差分析能快速定位是模糊度固定出错还是周跳探测阈值太高。这套做法后来被我用到很多项目里比单纯看定位精度指标有效得多。希望帮到你。本文还有配套的精品资源点击获取