标准实验报告——DSP与卫星导航教研室

1、电子科技大学通信与信息工程学院本科教学卫星与导航系列实验标准实验报告课程名称： l 定位与导航原理与应用 l 定位与导航工程 电子科技大学教务处制表电子科技大学通信与信息工程学院标准实验报告实验名称：导航信号传输模型仿真电子科技大学教务处制表电子科技大学实 验 报 告学生姓名： 学号： 提交日期： 得分 项目名称：实时卫星位置、速度和时间解算（PVT解算）及结果分析【实验目的】1) 理解实时卫星位置解算在卫星导航解算过程中所起的作用，了解为完成卫星位置解算所需的条件； 2) 了解 GPS 时间、卫星的额定轨道周期的含义，了解星历的构成、周期及应用条件； 3) 了解 Doppler 频移的成因、

2、作用以及根据已知条件预测 Doppler 频移的方法； 4) 了解 Doppler 频移的变化范围及其与卫星仰角之间的关系； 5) 能够根据实验数据编写求解 Doppler 频移的相关程序。【实验原理】实时卫星位置解算在整个导航解算过程中具有举足轻重的作用，通常我们为了获得接收机的地理位置，需要对卫星发射导航电文时的时间及运行速度有所了解，所以可以说，卫星的实时速度和时间是解算卫星实时位置的基础，而卫星的实时位置又是解算接收机三维位置坐标的基础。可见卫星实时位置、时间及速度在整个定位过程中的重要地位。一般来说要确定接收机的三维位置，需要同时解算出至少四颗卫星的实时位置。 卫星某一时刻发出的信号

3、可以分为三部分：载波（L1）、测距码（CA）、导航电文。对GPS某颗卫星进行实时位置的解算，需要已知这颗卫星的星历和周内时，这些信息都包含在速率为 50bps 的导航电文中（图3.1中的数据码）。导航电文通过测距码（C/A 码）进行扩频，然后用扩频的信号去调制频率为 L1的正弦波载波，然后卫星将调制后的载波信号播发出去。其模型可以用如下公式表示： 数据码测距码 载波 (3.1)其中和是调制幅度，是精码和粗码，它们都是对数据码的扩频码，数据码经过扩频后分为两路进行调制。在本地接收机收到卫星信号后，通过剥离载波L1，还原其扩频之后的信号，然后按照导航电文的格式最终将数据码编译成导航电文（数据码）。

4、它可分为5个部分：l 遥测字 导航电文前导+授权用户信息l 交接字 由于P码数据长，不易捕获，需要CA码的辅助来捕获l 数据块1 钟修正参数、期龄、星期编号等l 数据块2 星历（在子帧2、3中），用于计算卫星位置l 数据块3 历书（子帧4、5中）提供卫星布局，健康状况等信息本实验的一个重点在于通过导航电文来获取其卫星发射时间和星历从而得到卫星的实时位置。GPS卫星在空间中的位置是时间t的函数，要计算卫星的位置首先要收集齐时钟和星历参数，然后需要确定卫星的发射时刻。导航信号的发射时刻可以通过导航电文在每一子帧的的交接字中的周内时计数器（Z计数）得知，通过该计数器可以得到估算的发射时刻。同时在子帧

5、1中包含钟改正参数来对估算的时钟进行修正。导航电文中的数据块2是卫星星历信息，数据块3是卫星的历书信息。星历主要向用户提供有关计算卫星运行位置的信息，而历书主要向用户提供GPS卫星的概略星历及卫星的工作状态等。在数据块2（子帧2和3）中包含有许多重要星历参数。星历数据参数如下表所示（一颗卫星）: 表3.1 星历数据参数 参数（子帧2中） 意义 参数（子帧3中） 意义 轨道半径的正弦调和修正幅度 倾斜角的余弦调和修正项幅度 n计算值的平均移动误差 e在每星期历元轨道平面上升点的经度 参考时间的平近点角 倾斜角的正弦调和修正项的幅度 纬度辅角的余弦调和修正项幅度 参考时间的倾斜角 离心率 轨道半径

6、的正弦调和修正项幅度 纬度辅角的正弦调和修正项幅度 近地点的幅度 长半轴平方根 赤径的速度 参考时间星历 IODE数据与星历发布号 IDOT 倾斜角的速度后续导航解算单元根据导航电文中相应的参数进行卫星位置解算、各种实时误差的消除、本地接收机位置解算以及定位精度因子（DOP）的计算等工作。也就是说，根据收到的导航电文，接收端就可以通过相关公式计算出发送电文时刻卫星的大致位置，这对于解算出接收机的地理位置尤为重要。卫星的角速度和切向速度可以通过卫星轨道模型来进行估计。GPS卫星的额定轨道周期是半个恒星日，或者说 11 小时 58 分钟 2.05秒；各轨道接近于圆形，轨道半径（即从地球质心到卫星的

7、额定距离）大约为26560km。由此可以计算得到卫星的运行的角速度和切向速度（如图3.3）： =2/(11*3600+58*60+2.05)0.0001458 rad/s (3.2)然后通过角速度和已知的轨道半径rs（26560km）计算切向速度 vs=rs*26560km*0.00014583874m/s （3.3）本实验的另一个重要内容是 Doppler 频移的预测，即通过对卫星进行相隔时间为 1s 的多点测量（本实验给出了三点），进而估计 Doppler 频移。l Doppler 频移产生原理：由于卫星与接收机存在相对径向运动，因此信号接收频率会随这种相对运动而发生偏移，称为 Doppl

8、er频移，其直接表现是接收机接收到的卫星信号不在 L1 频点（1575.42MHz）上，而是在其基础上叠加了一个大约为-5KHz到+5KHz 左右的频率偏移。Doppler频移的存在给前端相关器进行频域搜索和卫星信号捕获带来了诸多困难。如果事先能够估算出Doppler频偏，就会降低捕获难度，缩短捕获时间，进而缩短接收机的启动时间。接收机的启动时间是衡量接收机性能好坏的重要参数之一，而实现卫星信号的快速捕获，缩短接收机的启动时间也是目前GNSS领域的热点问题。l Doppler 频移的计算：已知卫星位置和本地接收机的初始位置，根据空间两点间的距离公式，可以得出卫星到接收机的距离d。设某卫星在短时

9、间 t 内经过点 S1、S2，接收机到点 S1、S2 的距离分别为d1、d2,。在 t 相对较小（本实验取 t=1s）的情况下，我们可近似认是卫星与接收机在 t 时间内的平均相对径向运动速度，再将此速度转换为频率的形式就可以得到 Doppler 频移的估计值。设本地接收机的初始位置为 R（xr，yr，zr），卫星所经历的空间两点的坐标分别为 S1 （x1，y1，z1）、S2（x2，y2，z2），间隔时间为 t，卫星与接收机平均相对径向运动速度为 vd，光速为 c，Doppler 频移为 fd，则 Doppler 频移预测的具体公式如下所示：d1d222L1vd22 图 3.1其中，为一颗卫星不

10、同时刻到接收机的伪距，为两个时刻之间的径向速度。Doppler 频移与卫星的仰角有关，卫星仰角越大，Doppler 频移越小。当卫星的仰角为 90 度（即卫星在接收机正上方的天顶上）时，理论上 Doppler 频移为零。本实验根据卫星位置和本地接收机的初始位置算出卫星的仰角，来验证 Doppler 频移同卫星仰角之间的关系.【实验步骤】1) 查阅资料建立相应模型，在C/C+或者Matlab平台上根据星历数据及其定义实现对卫星实时位置的解算；2) 运行主程序以取得可视卫星的实时导航数据（如 GPS 时间、各颗卫星的星历等）。将实验平台仪器的USB端口接入电脑，待驱动安装成功后，打开实验一程序；3

11、) 在“选择 GPS 时刻”列表框的下拉菜单中，任意选择一个GPS 时刻。（注：北斗和 GPS 系统由于存在系统时差而具有不同的周内时。这里的 GPS 时刻，对于GPS卫星指其系统周内时，对于北斗卫星则表示将北斗的周内时加上系统时差换算之后的 GPS 系统周内时）；4) 在“所选时刻可视卫星星历”列表框中出现所选时刻天空中所有可视卫星的星历信息，如图 3.6 所示。选定一颗卫星，将“所选时刻可视卫星星历”中该卫星对应的参数输入到1）中的解算代码中，计算卫星位置。5) 在“选择卫星号”列表框的下拉菜单中，出现所选时刻天空中所有可视卫星的序号。北斗卫星的编号从 101 开始，即北斗 1 号星的编号

12、为 101。选择与4）中对应的卫星序号，在“卫星位置信息”中会列出所选时刻该卫星的实时位置如图3.7。对比该位置与之前代码解算的结果。并将其记录在表格中（表格一）；6) 在“卫星位置信息”列表框中同时会出现所选卫星在所选 GPS 时刻一秒和两秒后的所对应的 ECEF 坐标系下的三维坐标以及接收机在 ECEF 坐标系下的初始位置坐标，这些数据用于求解 Doppler 频移，根据附表记录其值（表格一）；7) 在“卫星位置信息”列表框中还会出现该卫星在 11 小时 58 分后的ECEF 位置坐标，这是根据卫星在所选 GPS 时刻的星历数据推算出来的，用以验证卫星的额定轨道周期。根据附表记录其值（表格

13、一）；8) 根据步骤 6）记录的数据，在C/C+或Matlab环境下编写代码，实现对Doppler 频移估值的求解，将所得数据记录在附表中（表格一）；9) 重复前面实验，记录并解算出所选时刻天空中所有可视卫星的相关数据，按附表格式将所得数据记录下来（表格二）；10) 重复前面实验，比较并分析不同时刻同一卫星的仰角、ECEF 坐标系下的坐标以及Doppler频移的差异和同一时刻不同卫星仰角、坐标及Doppler频移差异；11) 重复步骤 2 到步骤 11，选择不同时间段进行记录、求解、分析。【核心程序代码】% Find satellites position - %Restore semi-ma

14、jor axis a = eph(prn).sqrtA * eph(prn).sqrtA; %Time correction tk = check_t(time - eph(prn).t_oe); %Initial mean motion n0 = sqrt(GM / a3); %Mean motion n = n0 + eph(prn).deltan; %Mean anomaly M = eph(prn).M_0 + n * tk; %Reduce mean anomaly to between 0 and 360 deg M = rem(M + 2*gpsPi, 2*gpsPi); %In

15、itial guess of eccentric anomaly E = M; %- Iteratively compute eccentric anomaly - for ii = 1:10 E_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 loop break; end end %Reduce eccentric anomaly to between 0 and 360

16、deg E = rem(E + 2*gpsPi, 2*gpsPi); %Compute relativistic correction term dtr = F * eph(prn).e * eph(prn).sqrtA * sin(E); % %Calculate the true anomaly nu = atan2(sqrt(1 - eph(prn).e2) * sin(E), cos(E)-eph(prn).e); %Compute angle phi phi = nu + eph(prn).omega; %Reduce phi to between 0 and 360 deg phi

17、 = 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 radius r = a * (1 - eph(prn).e*cos(E) + . eph(prn).C_rc * cos(2*phi) + . eph(prn).C_rs * sin(2*phi); %Correct inclination i = eph(prn).i_0 + eph(prn).iDot * tk + . eph(

18、prn).C_ic * cos(2*phi) + . eph(prn).C_is * sin(2*phi); %Compute the angle between the ascending node and the Greenwich meridian Omega = eph(prn).omega_0 + (eph(prn).omegaDot - Omegae_dot)*tk - . Omegae_dot * eph(prn).t_oe; %Reduce to between 0 and 360 deg Omega = rem(Omega + 2*gpsPi, 2*gpsPi); %- Co

19、mpute 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、 按附表格式整理实验数据，并整理所编程序。（以一颗卫星为例）附表一卫星时间 卫星序号:10仰角Doppler 频移解

20、算坐标ECEF坐标308640-6898114.365083-6898114.84023169.3668914.836185542725942 KHz10231973.05697610231973.40023111426817.71196811426816.812402308641-6899460.241363-6899460.72645369. 3668914.836197123663138 KHz10232948.98024310232949.32603511424268.23045111424267.987032308642-6890805.435930-6890805.87403269

21、. 3668914.836208705610145 KHz10233924.76035410233924.78534311421718.24329911421718.20765311小时58 分后的坐标ECEF 位置坐标X:-6898413.749678Y:10229384.563445Z:11423157.6302542、对同一时刻不同卫星的 Doppler 频移进行比较，根据实际数据得出卫星仰角Doppler频移之间的关系。 附表二GPS 时间可视卫星序号仰角Doppler 频移3086411069. 3668914.836197123663138 KHz3086412016.607849

22、4.893789169568091 KHz3086413244.3657424.869897124663638 KHz由上表可得: 卫星的仰角越大，其多普勒频移越小。3、比较并分析不同时刻同一卫星的仰角、ECEF 坐标以及 Doppler 频移的差异。 附表三GPS世间卫星序号ECEF坐标仰角Doppler 频移30864110-6898413.74967869. 3668914.836197123663138 KHz10229384.56344511423157.63025430864210-6898413.20938469. 3668914.836208705610145 KHz 1022

23、9384.04893211423157.103970由上表分析得出什么结论？答：在时间相差很近的两个时刻里，卫星的坐标，仰角以及多普勒频移变化都很小，几乎不变化。仰角越大，多普勒频移越小。4、比较当前时刻卫星的 ECEF 坐标及由当前星历推算出的该卫星在 11 小时 58 分后的 ECEF 坐标，思考二者为何只是大致相同而不是绝对一致。 当前时刻与 11 小时 58 分后 2 号卫星的位置信息由下表给出：GPS世间卫星序号ECEF坐标仰角30864110X:-6898413.74967813.058204Y:10229384.563445Z:11423157.63025435172110X:-

24、6898413.03417813.058204Y:10229384.984645Z:11423157.204354分析其原因可能是： 可能是卫星的严格周期为11 小时 58 分2秒，而不是简单的11 小时 58 分；可能是因为接收机的计时误差；可能因为电离层、对流层误差导致。【意见与建议】对本实验过程及方法、手段的改进建议：希望可以把那款软件拷下来，这样方便学生在寝室调整代码。电子科技大学实 验 报 告学生姓名： 学号： 提交日期： 得分 项目名称：接收机位置解算及结果分析【实验目的】1）理解接收机位置解算的原理及基本公式；2）学习接收机位置解算的方法；3）分析将本地接收机钟差作为未知量进行导

25、航解算的原因；4）能够根据实验数据编写程序解算接收机位置。【实验原理】GPS 接收机位置解算指的是解出本地接收机的三维位置坐标，这是接收机的核心部分。为了得到GPS接收机的位置，需要得到卫星与接收机的距离，也就是伪距。伪距由传播时间和光速的乘积得到：(4.1)传播时间可以由发射时刻与接收时刻相减得到。由于卫星位置已知，这样就可以把接收机定位在以卫星为球心半径为的球面上。求出本地接收机的位置也就是算出多个以卫星为球心的球面的相交点。因为需要解出本地接收机的三维坐标，存在三个未知量，所以理论上至少需要三颗卫星的位置等信息来解算接收机位置。实际上，上述求解过程是在假定卫星时钟与接收机时钟同步的前提下

26、进行的。然而二者并不同步，彼此都存在这钟差。GPS 卫星时钟与世界协调时（UTC，即格林尼治时间）的内在系统时间标度是同步的，且卫星钟差可根据导航电文中给出的有关钟差参数加以修正前面提到，由于GPS卫星属于中轨道卫星，卫星的信号在传输过程中会产生诸如对流层误差，电离层误差，相对论效应的等各种实时传输误差。因此，之前求得的距离实际上市包含上述各种误差的伪距，直接带入上述方程求解误差会很大，需要通过实验一所介绍的方法进行误差消除和伪距的矫正，即有：(4.2)其中为电离层延迟误差，为对流层延迟误差，为卫星的钟差。当然上式仅对电离层延迟、对流层延迟和卫星钟差对伪距造成的影响进行了矫正，而未考虑其他如多

27、径效应，接收机噪声等方面的误差影响。修正误差后其基准频率的稳定度为10-13左右。而本地接收机的钟差在解算时包含在接收机位置中，严重影响定位精度。为解决这一问题，通常情况下，我们将接收机钟差和本地接收机的三维位置坐标一起作为未知量进行求解。因此，至少需要4颗卫星，4个观测量，才能进行方程组的求解。而解出的接收机钟差可以用来校正本地接收机的时钟，这也使得接收机具有了授时功能。求解卫星位置的基本方程组为：(4.3)(4.4)(4.5)(4.6)其中，卫星的三维坐标为，接收机的三维坐标位置为，为接收机钟差，c为光速，为各颗卫星的距离。根据这四个方程，接收机可以开始对卫星位置进行解算。因为这个方程组是

28、非线性的联立方程，所以直接求解是很难的。要解出这4个解，有泰勒级数法、最小二乘法等，其中最常用的是最小二乘法。l 利用最小二乘法求解接收机位置坐标：首先，设置方程初始解。迭代前给方程组的4个未知数设定一个初始值。初始值的设置分为两种情况，若是首次定位，则全部设为0。若已经成功定位，则将上一次的结果设置为本次迭代的初始值。然后将方程组线性化：(4.7)对方程微分，结果为(4.8)将（4.8）写成矩阵形式可得：(4.9)其中，(4.10)(4.11)(4.12)在这里，和表示第k-1次迭代求出的钟差和接收机与对应卫星的距离。k=1时表示第1步中设置的初始值。之后利用最小二乘法求解方程组（4.9）：

29、(4.13)将得出的结果进行更新：(4.14)最后，重复以上步骤，反复迭代。一般来说，迭代7次左右即可得到接收机的坐标及时钟钟差。l 本地接收机钟差作为未知量进行解算的原因：之所以对接收机钟差进行求解，其重要原因在于接收机钟差对定位精度的影响不可忽略。而通过实验发现，同一时刻不同卫星所对应的接收机钟差大致相同，因此，这也为我们将接收机钟差作为一个未知量参与运算提供了依据。假定某一时刻，实际测得的伪距用Pseudorange表示，由两点间距离公式得出的真实距离用Range表示，则(4.15)表示该时刻的所有实时误差之和。（提示：这里的真实距离可用本实验的验证 1）中得到的近似接收机真实位置坐标来

30、求取）根据实验一的介绍，我们已经了解了电离层延迟、对流层误差等主要实时传输误差的误差范围，可结合实验数据进行比较可知，这只占R中很小的一部分，而其中绝大部分误差来源于接收机的钟差。这也从侧面反映出接收机钟差对定位精度产生的影响远远大于其他实时传输误差所造成的影响。同时，比较同一时刻各可视卫星的R会发现，其间差异不大。可以得出一个一般结论，即一般情况下，同一时刻不同卫星所对应的接收机钟差大致相同。因此，可以将接收机钟差作为一个未知变量进行求解。【实验步骤】1）运行主程序以取得求解接收机位置所需的实时参数。包括可视卫星的实时导航数据、实时传输误差、伪距等；2）运行实验四程序，步骤1中记录的所有GP

31、S时间就会出现在“选择 GPS 时刻”列表框的下拉菜单中，选择任意一个GPS时刻；理解程序中列出来的各参数的含义，并按照表一进行记录；3）利用程序进行本地接收机位置解算。本地接收机位置解算需要至少 4 颗卫星的信息，如果所选GPS时刻可视卫星数小于4颗，则需要重新选择其它GPS 时刻进行解算。4）得到接收机的三维坐标信息以及其经纬度和高度，按照表一对数据进行记录；5）选择一个任意GPS时刻，将当前时刻可视卫星的数据信息、各卫星的ECEF坐标、伪距以及计算所得的接收机位置坐标分别记录在表二的相应位置；6）结合星历数据和表二中的数据在Matlab环境下编写程序，求取接收机位置坐标，按表二记录结果；

32、7）选择任意一个 GPS 时刻，计算R，数据记录在表三。并比较同一时刻不同卫星所对应的R之间的关系。【数据记录与分析】1) 按附表格式整理实验数据和程序。表一：GPS 时间：_95377_可视卫星数目为_10_卫星序号伪距（m）电离层延迟（t）对流层误差（m）卫星位置（m）接收机位置（m）1421120718.53-0.0012985.472487-22118920.702871 23058042.153582 14429021.492768 -1690547.538724 1621120719.53-0.0043883.840251644544.586168 21911537.650155

33、26103781.609453 3565396.441737 2021120720.53-0.0010765.79253914201414.036339 23234877.974117 18409787.843196 12714826.561875 2321120721.53-0.0008676.18513413064535.114457 23569400.999140 23059866.086854 21909917.459505 23235325.782198 23567500.605946 9022469.272299 21439646.169100 1421120722.53-0.00

34、12655.476648-22118792.450246 14427952.924797 -1700138.905040 1621120723.53-0.0043953.83841643720.862923 26102479.027585 3575047.100148 2021120724.53-0.0010755.79347514201656.478025 18414817.143634 12707399.910391 2321120725.53-0.0008696.18018413058337.438534 9027668.223318 21441117.181021 表二：GPS 时间：

35、_106285_卫星位置（m）伪距（m）接收机位置（m）-22118920.702871 21120718.5323058042.153582 14429021.492768 -1690547.538724 644544.586168 21120719.5321911537.650155 26103781.609453 3565396.441737 14201414.036339 21120720.5323234877.974117 18409787.843196 12714826.561875 13064535.114457 21120721.5323569400.999140 230598

36、66.086854 21909917.459505 23235325.782198 23567500.605946 9022469.272299 21439646.169100 -22118792.450246 21120722.5314427952.924797 -1700138.905040 643720.862923 21120723.5326102479.027585 3575047.100148 14201656.478025 21120724.5318414817.143634 12707399.910391 13058337.438534 21120725.539027668.223318 21441117.181021 Matlab程序代码：表三：GPS 时间：_106288_卫星序号卫星位置（m）接收机位置(m)Range（m）伪距（m）R（m）14-22118920.702871 -2182286.6