电子科技大学通信与信息工程学院本科教学卫星与导航系列实验标准实验报告课程名称：●定位与导航原理与应用●定位与导航工程电子科技大学教务处制表电子科技大学通信与信息工程学院标准实验报告实验名称：导航信号传输模型仿真电子科技大学教务处制表电子科技大学实验报告学生姓名：侯玉皓学号：2012019030016 提交日期：2015.6.24项目名称：实时卫星位置、速度和时间解算（PVT解算）及结果分析【实验目的】1)理解实时卫星位置解算在卫星导航解算过程中所起的作用，了解为完成卫星位置解算所需的条件；2)了解 GPS 时间、卫星的额定轨道周期的含义，了解星历的构成、周期及应用条件；3)了解 Doppler 频移的成因、作用以及根据已知条件预测 Doppler 频移的方法；4)了解 Doppler 频移的变化范围及其与卫星仰角之间的关系；5)能够根据实验数据编写求解 Doppler 频移的相关程序。
【实验原理】实时卫星位置解算在整个导航解算过程中具有举足轻重的作用，通常我们为了获得接收机的地理位置，需要对卫星发射导航电文时的时间及运行速度有所了解，所以可以说，卫星的实时速度和时间是解算卫星实时位置的基础，而卫星的实时位置又是解算接收机三维位置坐标的基础。
可见卫星实时位置、时间及速度在整个定位过程中的重要地位。
一般来说要确定接收机的三维位置，需要同时解算出至少四颗卫星的实时位置。
卫星某一时刻发出的信号可以分为三部分：载波（L1）、测距码（C\A ）、导航电文。
对GPS 某颗卫星进行实时位置的解算，需要已知这颗卫星的星历和周内时，这些信息都包含在速率为 50bps 的导航电文中（图3.1中的数据码）。
导航电文通过测距码（C/A 码）进行扩频，然后用扩频的信号去调制频率为 L1的正弦波载波，然后卫星将调制后的载波信号播发出去。
其模型可以用如下公式表示：11111()(()())cos()(()())sin()L p i i L c i i L S t A P t D t t A C t D t t =?+?w y w y (3.1)其中p A 和c A 是调制幅度，,i i P C 是精码和粗码，它们都是对数据码i D 的扩频码，数据码经过扩频后分为两路进行调制。
在本地接收机收到卫星信号后，通过剥离载波L1，还原其扩频之后的信号，然后按照导航电文的格式最终将数据码编译成导航电文（数据码）。
它可分为5个部分：● 遥测字 导航电文前导+授权用户信息● 交接字 由于P 码数据长，不易捕获，需要C\A 码的辅助来捕获 ● 数据块1 钟修正参数、期龄、星期编号等● 数据块2 星历（在子帧2、3中），用于计算卫星位置● 数据块3 历书（子帧4、5中）提供卫星布局，健康状况等信息本实验的一个重点在于通过导航电文来获取其卫星发射时间和星历从而得到卫星的实时位置。
GPS 卫星在空间中的位置是时间t 的函数，要计算卫星的位置首先要收集齐时钟和星历参数，然后需要确定卫星的发射时刻。
导航信号的发射时刻可以通过导航电文在每一子帧的的交接字中的周内时计数器（Z 计数）得知，通过该计数器可以得到估算的发射时刻c t 。
同时在子帧1中包含钟改正参数012,,a a a 来对估算的时钟进行修正。
导航电文中的数据块2是卫星星历信息，数据块3是卫星的历书信息。
星历主要向用户提供有关计算卫星运行位置的信息，而历书主要向用户提供GPS 卫星的概略星历及卫星的工作状态等。
数据码 测距码 载波在数据块2（子帧2和3）中包含有许多重要星历参数。
星历数据参数如下表所示（一颗卫星）:表3.1 星历数据参数后续导航解算单元根据导航电文中相应的参数进行卫星位置解算、各种实时误差的消除、本地接收机位置解算以及定位精度因子（DOP）的计算等工作。
也就是说，根据收到的导航电文，接收端就可以通过相关公式计算出发送电文时刻卫星的大致位置，这对于解算出接收机的地理位置尤为重要。
卫星的角速度和切向速度可以通过卫星轨道模型来进行估计。
GPS卫星的额定轨道周期是半个恒星日，或者说 11 小时 58 分钟 2.05秒；各轨道接近于圆形，轨道半径（即从地球质心到卫星的额定距离）大约为26560km。
由此可以计算得到卫星的运行的角速度和切向速度（如图3.3）：ω=2π/(11*3600+58*60+2.05)≈0.0001458 rad/s (3.2) 然后通过角速度ω和已知的轨道半径r s （26560km ）计算切向速度 v s =r s *ω≈26560km*0.0001458≈3874m/s（3.3）本实验的另一个重要内容是 Doppler 频移的预测，即通过对卫星进行相隔时间为 1s 的多点测量（本实验给出了三点），进而估计 Doppler 频移。
● Doppler 频移产生原理：由于卫星与接收机存在相对径向运动，因此信号接收频率会随这种相对运动而发生偏移，称为 Doppler 频移，其直接表现是接收机接收到的卫星信号不在 L1 频点（1575.42MHz ）上，而是在其基础上叠加了一个大约为-5KHz 到+5KHz 左右的频率偏移。
Doppler 频移的存在给前端相关器进行频域搜索和卫星信号捕获带来了诸多困难。
如果事先能够估算出Doppler 频偏，就会降低捕获难度，缩短捕获时间，进而缩短接收机的启动时间。
接收机的启动时间是衡量接收机性能好坏的重要参数之一，而实现卫星信号的快速捕获，缩短接收机的启动时间也是目前GNSS 领域的热点问题。
● Doppler 频移的计算：已知卫星位置和本地接收机的初始位置，根据空间两点间的距离公式，可以得出卫星到接收机的距离d 。
设某卫星在短时间 t 内经过点 S 1、S 2，接收机到点 S 1、S 2 的距离分别为d 1、d 2,。
在 t 相对较小（本实验取 t=1s ）的情况下，我们可近似认21d d t-是卫星与接收机在 t 时间内的平均相对径向运动速度，再将此速度转换为频率的形式就可以得到 Doppler 频移的估计值。
设本地接收机的初始位置为 R （xr ，yr ，zr ），卫星所经历的空间两点的坐标分别为 S1 （x1，y1，z1）、S2（x2，y2，z2），间隔时间为 t ，卫星与接收机平均相对径向运动速度为 vd ，光速为 c ，Doppler 频移为 fd ，则 Doppler 频移预测的具体公式如下所示：1d 2d =21d d d v t -=1d d L f v c =?d1vd22图 3.1其中，123,,d d d 为一颗卫星不同时刻到接收机的伪距，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 semi-major axisa = eph(prn).sqrtA * eph(prn).sqrtA;%Time correctiontk = check_t(time - eph(prn).t_oe);%Initial mean motionn0 = 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 degM = rem(M + 2*gpsPi, 2*gpsPi);%Initial guess of eccentric anomalyE = M;%--- Iteratively compute eccentric anomaly ----------------------------for ii = 1:10E_old = E;E = M + eph(prn).e * sin(E);dE = rem(E - E_old, 2*gpsPi);if abs(dE) < 1.e-12% Necessary precision is reached, exit from the loopbreak;endend%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); %Ïà¶ÔÂÛÐÞÕýÏî%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 degphi = rem(phi, 2*gpsPi);%Correct argument of latitudeu = 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)*tk - ...Omegae_dot * eph(prn).t_oe;%Reduce to between 0 and 360 degOmega = rem(Omega + 2*gpsPi, 2*gpsPi);%--- Compute satellite coordinates ------------------------------------satPositions(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、按附表格式整理实验数据，并整理所编程序。
