卫星导航实验报告

本科教学卫星与导航系列实验标准实验报告课程名称：•定位与导航原理与应用•定位与导航工程标准实验报告实验名称：导航信号传输模型仿真电子科技大学实验报告学生姓名：侯玉皓 学号：2012019030016提交日期：2015.6.24^分项目名称：实时卫星位置、速度和时间解算Pvt解算)及结果分析【实验目的】1) 理解实时卫星位置解算在卫星导航解算过程中所起的作用，了解为完成 卫星位置解算所需的条件；2) 了解GPS时间、卫星的额定轨道周期的含义，了解星历的构成、周期 及应用条件；3) 了解Doppler频移的成因、作用以及根据已知条件预测 Doppler频移 的方法；4) 了解Doppler频移的变化范围及其与卫星仰角之间的关系；5) 能够根据实验数据编写求解Doppler频移的相关程序实验原理】实时卫星位置解算在整个导航解算过程中具有举足轻重的作用，通常我们为 了获得接收机的地理位置，需要对卫星发射导航电文时的时间及运行速度有所 了解，所以可以说，卫星的实时速度和时间是解算卫星实时位置的基础，而卫 星的实时位置又是解算接收机三维位置坐标的基础可见卫星实时位置、时间 及速度在整个定位过程中的重要地位。

一般来说要确定接收机的三维位置，需 要同时解算出至少四颗卫星的实时位置卫星某一时刻发出的信号可以分为三部分：载波(L1)、测距码(C\A)、导 航电文对GPS某颗卫星进行实时位置的解算，需要已知这颗卫星的星历和周内 时，这些信息都包含在速率为50bps的导航电文中(图3.1中的数据码)导航 电文通过测距码C/A码)进行扩频，然后用扩频的信号去调制频率为的正弦 波载波，然后卫星将调制后的载波信号播发出去其模型可以用如下公式表示：S (t) = A (P (t)㊉ D (t)) cos(① t +w ) + A (C (t)㊉ D (t)) sin(① t + (3.1) L1 p jAi ■ L1 1 c i i L1w 1) \ T H测距码数据码载波其中^和Ac是调制幅度，P；C是精码和粗码，它们都是对数据码D的扩 频码，数据码经过扩频后分为两路进行调制在本地接收机收到卫星信号后， 通过剥离载波L1，还原其扩频之后的信号，然后按照导航电文的格式最终将数 据码编译成导航电文(数据码)它可分为5个部分：•遥测字导航电文前争授权用户信息•交接字 由于P码数据长，不易捕获，需要AA码的辅助来捕获•数据块1钟修正参数、期龄、星期编号等•数据块2星历(在子帧2、3中)，用于计算卫星位置•数据块3历书(子帧4、5中)提供卫星布局，健康状况等信息本实验的一个重点在于通过导航电文来获取其卫星发射时间和星历从而得到 卫星的实时位置。

GPS卫星在空间中的位置是时间t的函数，要计算卫星的位置首先要收集齐 时钟和星历参数，然后需要确定卫星的发射时刻导航信号的发射时刻可以通过 导航电文在每一子帧的的交接字中的周内时计数器(Z计数)得知，通过该计数 器可以得到估算的发射时刻［同时在子帧1中包含钟改正参数a0,%,a2来对 估算的时钟进行修正导航电文中的数据块2是卫星星历信息，数据块3是卫星的历书信息星历主要向用户提供有关计算卫星运行位置的信息，而历书主要向用户提供?、卫星的概略星历及卫星的工作状态等在数据块2（子帧2和3）中包含有许多重要星历参数星历数据参数如 下表所示（一颗卫星）：表3.1星历数据参数参数（子帧2中）意义参数（子帧 3中）意义crs轨道半径的正弦调和修正幅度c ic倾斜角的余弦调和修正项幅度An计算值的平均移动 误差Qe在每星期历元轨道 平面上升点的经度m 0参考时间的平近点 角c is倾斜角的正弦调和修正项的幅度cuc纬度辅角的余弦调和修正项幅度i0参考时间的倾斜角es离心率crc轨道半径的正弦调和修正项幅度cus纬度辅角的正弦调和修正项幅度⑦近地点的幅度as长半轴平方根Q赤径的速度toe参考时间星历IODE数据与星历发布号IDOT倾斜角的速度后续导航解算单元根据导航电文中相应的参数进行卫星位置解算、各种实时 误差的消除、本地接收机位置解算以及定位精度因子（DOP）的计算等工作。

也 就是说，根据收到的导航电文，接收端就可以通过相关公式计算出发送电文时刻 卫星的大致位置，这对于解算出接收机的地理位置尤为重要卫星的角速度和切向速度可以通过卫星轨道模型来进行估计GPS卫星的额 定轨道周期是半个恒星日，或者说1小时58分钟2.05秒；各轨道接近于圆形，轨道半径（即从地球质心到卫星的额定距离）大约为26560km由此可以计算得 到卫星的运行的角速度和切向速度（如图3）：s=2n/（11*3600+58*60+2.05）=0.0001458 rad/s （3.2）然后通过角速度3和已知的轨道半径rs（26560km）计算切向速度vs=rs* 3 2 26560km*0.0001458 ^ 3874m/s （ 3.3）本实验的另一个重要内容是Doppler频移的预测，即通过对卫星进行相隔时 间为1s的多点测量（本实验给出了三点进而估计Doppler频移• Doppler频移产生原理：由于卫星与接收机存在相对径向运动，因此信号接收频率会随这种相对运 动而发生偏移，称为Doppler频移，其直接表现是接收机接收到的卫星信号不 在L1频点（1575.42MHz）上，而是在其基础上叠加了一个大约为~5KHz到 +5KHz左右的频率偏移。

Doppler频移的存在给前端相关器进行频域搜索和卫星 信号捕获带来了诸多困难如果事先能够估算出Doppler频偏，就会降低捕获 难度，缩短捕获时间，进而缩短接收机的启动时间接收机的启动时间是衡量 接收机性能好坏的重要参数之一，而实现卫星信号的快速捕获，缩短接收机的 启动时间也是目前GNSS领域的热点问题• Doppler频移的计算：已知卫星位置和本地接收机的初始位置，根据空间两点间的距离公式，可 以得出卫星到接收机的距离d设某卫星在短时间t内经过点S「S2,接收机 到点S「S2的距离分别为d「d2在t相对较小（本实验取t=1s）的情况下，\d -d我们可近似认是卫星与接收机在t时间内的平均相对径向运动速度，再t将此速度转换为频率的形式就可以得到oppler频移的估计值设本地接收机的初始位置为R （xr，yr，zr），卫星所经历的空间两点的坐 标分别为S1 （x1，y1，z1）、S2（x2, y2, z2），间隔时间为t，卫星与接收机 平均相对径向运动速度为vd，光速为c， Doppler频移为fd，贝V Doppler频 移预测的具体公式如下所示：vd ]22d1d=Wi P +°1 -乃）2 +（% p d2=正2 -y +化-y）2 +§ -y|d2 - d 12 lLv 二r C fv=图3.1其中，d1 d四为一颗卫星不同时刻到接收机的伪距，V为两个时刻之间 的径向速度。

Doppler频移与卫星的仰角有关，卫星仰角越大，Doppler频移越小当卫 星的仰角为90度（即卫星在接收机正上方的天顶上）时，理论上Doppler频 移为零本实验根据卫星位置和本地接收机的初始位置算出卫星的仰角，来验 证Doppler频移同卫星仰角之间的关系.【实验步骤】1） 查阅资料建立相应模型，在C/C++或者Matlab平台上根据星历数据及其定 义实现对卫星实时位置的解算；2） 运行主程序以取得可视卫星的实时导航数据（如GPS时间、各颗卫星的星 历等）将实验平台仪器的USB端口接入电脑，待驱动安装成功后，打开 实验一程序；3） 在“选择GPS时刻”列表框的下拉菜单中，任意选择一个GPS时刻注：北斗和GPS系统由于存在系统时差而具有不同的周内时这里的GPS时刻，对于GPS卫星指其系统周内时，对于北斗卫星则表示将北斗的 周内时加上系统时差换算之后的GPS系统周内时）；4） 在“所选时刻可视卫星星历”列表框中出现所选时刻天空中所有可视卫星的 星历信息，如图3.6所示选定一颗卫星，将“所选时刻可视卫星星历”中 该卫星对应的参数输入到1）中的解算代码中，计算卫星位置5） 在“选择卫星号”列表框的下拉菜单中，出现所选时刻天空中所有可视卫 星的序号。

北斗卫星的编号从101开始，即北斗1号星的编号为101 选择与4）中对应的卫星序号，在“卫星位置信息”中会列出所选时刻该卫 星的实时位置如图3.7对比该位置与之前代码解算的结果并将其记录在 表格中（表格一;）6） 在“卫星位置信息”列表框中同时会出现所选卫星在所选GPS时刻一秒和 两秒后的所对应的ECEF坐标系下的三维坐标以及接收机在ECEF坐标系 下的初始位置坐标，这些数据用于求解Doppler频移，根据附表记录其值（表格一）；7） 在“卫星位置信息”列表框中还会出现该卫星在11小时58分后的ECEF 位置坐标，这是根据卫星在所选GPS时刻的星历数据推算出来的，用以验 证卫星的额定轨道周期根据附表记录其值（表格一）；8） 根据步骤6）记录的数据，在C/C++或Matlab环境下编写代码，实现对 Doppler频移估值的求解，将所得数据记录在附表中（表格一）；9） 重复前面实验，记录并解算出所选时刻天空中所有可视卫星的相关数据，按 附表格式将所得数据记录下来（表格二）10） 重复前面实验，比较并分析不同时刻同一卫星的仰角ECEF坐标系下的坐标 以及Doppler频移的差异和同一时刻不同卫星仰角、坐标及Doppler频移差 异；11） 重复步骤2到步骤11，选择不同时间段进行记录、求解、分析。

【核心程序代码】%% Find satellite's position%Restore semimajor axis a = eph（prn）.sqrtA 火 eph（prn）.sqrtA%Time correctiontk = check_t（time eph（prn）.t_oe）%Initial mean motion n0 = sqrt（GM / a八3） %Mean motionn = n0 + eph（prn）.deltan%Mean anomalyM = eph（prn）.M_0 + n * tk %Reduce mean anomaly to between 0 and 360 deg M = rem（M + 2火gpsPi, 2火gpsPi）%Initial guess of eccentric anomalyIteratively compute eccentric anomalyfor ii1:10EolddEM + eph(prn).e 火 sin(E) rem(E E_old, 2火gpsPi)if abs(dE) < 1.e12% Necessary precision is reached, exit from the loop breakendend%Reduce eccentric anomaly to between 0 and 360 degE = rem(E + 2火gpsPi, 2火gpsPi) %Compute relativistic correction termdtr = F * eph(prn).e 火 eph(prn).sqrtA 火 sin(E) %Ia^OAUBPOyIi%Calculate the true anomalynu = atan2(sqrt(1 eph(prn).e八2) * sin(E), cos(E)eph(prn).e)%Compute angle phiphi = nu + eph(prn).omega %Reduce phi to between 0 and 360 deg phi = rem(phi, 2火gpsPi)%Correct argument of latitude u = phi + ...eph(prn).C_uc 火 cos(2火phi) + ...eph(prn).C_us 火 sin(2火phi)%Correct radiusr = a * (1 eph(prn).e*cos(E)) + ...eph(prn).C_rc * cos(2*phi) + ...eph(prn).C_rs * sin(2*phi)%Correct inclinationi = eph(prn).i_0 + eph(prn).iDot * tk + ... eph(prn).C_ic * cos(2*phi) + ... eph(prn).C_is * sin(2*phi)%Compute the angle between the ascending node and the Greenwich meridianOmega = eph(prn).omega_0 + (eph(prn).omegaDot Omegae_dot)*tkOmegae_dot * eph(prn).t_oe%Reduce to between 0 and 360 degOmega = rem(Omega + 2*gpsPi, 2*gpsPi)Compute satellite coordinatessatPositions(1, satNr) = cos(u)*r 火 cos(Omega) sin(u)*r 火 cos(i)*sin(Omega)satPositions(2, satNr) = cos(u)*r 火 sin(Omega) + sin(u)*r 火 cos(i)*cos(Omega)satPositions(3, satNr) = sin(u)*r * sin(i);【结果记录与分析】1、按附表格式整理实验数据，并整理所编程序以一颗卫星为例)附表一卫星时间卫星序号:10仰角Doppler 频移解算坐标ECEF坐标308640-6898114.365083-6898114.84023169.366891°4.836185542725942 KHz10231973.05697610231973.40023111426817.71196811426816.812402308641-6899460.241363-6899460.72645369. 366891°4.836197123663138 KHz10232948.98024310232949.32603511424268.23045111424267.987032308642-6890805.435930-6890805.87403269. 366891°4.836208705610145 KHz10233924.76035410233924.78534311421718.24329911421718.20765311小时58X:-6898413.749678分后的坐标Y:10229384.563445ECEF 位Z:11423157.630254置坐标2、对同一时刻不同卫星的Doppler频移进行比较，根据实际数据得出卫星仰 角Doppler频移之间的关系。

附表二GPS时间可视卫星序号仰角Doppler 频移3086411069. 366891°4.836197123663138 KHz3086412016.607849°4.893789169568091 KHz3086413244.365742°4.869897124663638 KHz由上表可得：卫星的仰角越大，其多普勒频移越小3、比较并分析不同时刻同一卫星的仰角、ECEF坐标以及Doppler频移的差 异附表三GPS世间卫星序号ECEF坐标仰角Doppler 频移30864110-6898413.74967869. 366891°4.836197123663138 KHz10229384.56344511423157.63025430864210-6898413.20938469. 366891°4.836208705610145 KHz10229384.04893211423157.103970由上表分析得出什么结论？答：在时间相差很近的两个时刻里，卫星的坐标，仰角以及多普勒频移变化都很小，几乎不变化仰角越大，多普勒频移越小4、比较当前时刻卫星的ECEF坐标及由当前星历推算出的该卫星在11小时58分后的ECEF坐标，思考二者为何只是大致相同而不是绝对一致。

当前时刻与11小时58分后2号卫星的位置信息由下表给出：GPS世间卫星序号ECEF坐标仰角30864110X:-6898413.74967813.058204°Y:10229384.563445Z:11423157.63025435172110X:-6898413.03417813.058204°Y:10229384.984645Z:11423157.204354分析其原因可能是：可能是卫星的严格周期为11小时58分2秒，而不是简单的11小时58分;可能是因为接收机的计时误差；可能因为电离层、对流层误差导致意见与建议】对本实验过程及方法、手段的改进建议：希望可以把那款软件拷下来，这样方便学生在寝室调整代码。