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

电子科技大学通信与信息工程学院本科教学卫星与导航系列实验标准实验报告课程名称：定位与导航原理与应用定位与导航工程电子科技大学教务处制表电子科技大学通信与信息工程学院标准实验报告实验名称：导航信号传输模型仿真电子科技大学教务处制表电子科技大学实验报告学生姓名：学号：提交日期：项目名称：实时卫星位置、速度和时间解算（PVT解算）及结果分析【实验目的】1)理解实时卫星位置解算在卫星导航解算过程中所起的作用，了解为完成卫星位置解算所需的条件；2)了解GPS时间、卫星的额定轨道周期的含义，了解星历的构成、周期及应用条件；3)了解Doppler频移的成因、作用以及根据已知条件预测Doppler频移的方法；4)了解Doppler频移的变化范围及其与卫星仰角之间的关系；5)能够根据实验数据编写求解Doppler频移的相关程序。【实验原理】实时卫星位置解算在整个导航解算过程中具有举足轻重的作用，通常我们为了获得接收机的地理位置，需要对卫星发射导航电文时的时间及运行速度有所了解，所以可以说，卫星的实时速度和时间是解算卫星实时位置的基础，而卫星的实时位置又是解算接收机三维位置坐标的基础。可见卫星实时位置、时间及速度在整个定位过程中的重要地位。一般来说要确定接收机的三维位置，需要同时解得分算出至少四颗卫星的实时位置。卫星某一时刻发出的信号可以分为三部分：载波（L1）、测距码（C\A）、导航电文。对GPS某颗卫星进行实时位置的解算，需要已知这颗卫星的星历和周内时，这些信息都包含在速率为50bps的导航电文中（图3.1中的数据码）。导航电文通过测距码（C/A码）进行扩频，然后用扩频的信号去调制频率为L1的正弦波载波，然后卫星将调制后的载波信号播发出去。其模型可以用如下公式表示：11111()(()())cos()(()())sin()LpiiLciiLStAPtDttACtDtt=?+?wywy(3.1)其中pA和cA是调制幅度，,iiPC是精码和粗码，它们都是对数据码iD的扩频码，数据码经过扩频后分为两路进行调制。在本地接收机收到卫星信号后，通过剥离载波L1，还原其扩频之后的信号，然后按照导航电文的格式最终将数据码编译成导航电文（数据码）。它可分为5个部分：遥测字导航电文前导+授权用户信息交接字由于P码数据长，不易捕获，需要C\A码的辅助来捕获数据块1钟修正参数、期龄、星期编号等数据块2星历（在子帧2、3中），用于计算卫星位置数据块3历书（子帧4、5中）提供卫星布局，健康状况等信息本实验的一个重点在于通过导航电文来获取其卫星发射时间和星历从而得到卫星的实时位置。GPS卫星在空间中的位置是时间t的函数，要计算卫星的位置首先要收集齐时钟和星历参数，然后需要确定卫星的发射时刻。导航信号的发射时刻可以通过导航电文在每一子帧的的交接字中的周内时计数器（Z计数）得知，通过该计数器可以得到估算的发射时刻ct。同时在子帧1中包含钟改正参数012,,aaa来对估算的时钟进行修正。导航电文中的数据块2是卫星星历信息，数据块3是卫星的历书信息。星历主要向用户提供有关计算卫星运行位置的信息，而历书主要向用户提供GPS卫星的概略星历及卫星的工作状态等。数据码测距码载波在数据块2（子帧2和3）中包含有许多重要星历参数。星历数据参数如下表所示（一颗卫星）:表3.1星历数据参数后续导航解算单元根据导航电文中相应的参数进行卫星位置解算、各种实时误差的消除、本地接收机位置解算以及定位精度因子（DOP）的计算等工作。也就是说，根据收到的导航电文，接收端就可以通过相关公式计算出发送电文时刻卫星的大致位置，这对于解算出接收机的地理位置尤为重要。卫星的角速度和切向速度可以通过卫星轨道模型来进行估计。GPS卫星的额定轨道周期是半个恒星日，或者说11小时58分钟2.05秒；各轨道接近于圆形，轨道半径（即从地球质心到卫星的额定距离）大约为26560km。由此可以计算得到卫星的运行的角速度和切向速度（如图3.3）：参数（子帧2中）意义参数（子帧3中）意义rsc轨道半径的正弦调和修正幅度icc倾斜角的余弦调和修正项幅度Δn计算值的平均移动误差Ωe在每星期历元轨道平面上升点的经度0m参考时间的平近点角isc倾斜角的正弦调和修正项的幅度ucc纬度辅角的余弦调和修正项幅度0i参考时间的倾斜角se离心率rcc轨道半径的正弦调和修正项幅度usc纬度辅角的正弦调和修正项幅度ω近地点的幅度sa长半轴平方根Ω赤径的速度eot参考时间星历IODE数据与星历发布号IDOT倾斜角的速度ω=2π/(11*3600+58*60+2.05)≈0.0001458rad/s(3.2)然后通过角速度ω和已知的轨道半径rs（26560km）计算切向速度vs=rs*ω≈26560km*0.0001458≈3874m/s（3.3）本实验的另一个重要内容是Doppler频移的预测，即通过对卫星进行相隔时间为1s的多点测量（本实验给出了三点），进而估计Doppler频移。Doppler频移产生原理：由于卫星与接收机存在相对径向运动，因此信号接收频率会随这种相对运动而发生偏移，称为Doppler频移，其直接表现是接收机接收到的卫星信号不在L1频点（1575.42MHz）上，而是在其基础上叠加了一个大约为-5KHz到+5KHz左右的频率偏移。Doppler频移的存在给前端相关器进行频域搜索和卫星信号捕获带来了诸多困难。如果事先能够估算出Doppler频偏，就会降低捕获难度，缩短捕获时间，进而缩短接收机的启动时间。接收机的启动时间是衡量接收机性能好坏的重要参数之一，而实现卫星信号的快速捕获，缩短接收机的启动时间也是目前GNSS领域的热点问题。Doppler频移的计算：已知卫星位置和本地接收机的初始位置，根据空间两点间的距离公式，可以得出卫星到接收机的距离d。设某卫星在短时间t内经过点S1、S2，接收机到点S1、S2的距离分别为d1、d2,。在t相对较小（本实验取t=1s）的情况下，我们可近似认21ddt-是卫星与接收机在t时间内的平均相对径向运动速度，再将此速度转换为频率的形式就可以得到Doppler频移的估计值。设本地接收机的初始位置为R（xr，yr，zr），卫星所经历的空间两点的坐标分别为S1（x1，y1，z1）、S2（x2，y2，z2），间隔时间为t，卫星与接收机平均相对径向运动速度为vd，光速为c，Doppler频移为fd，则Doppler频移预测的具体公式如下所示：2221111()()()rrrdxxyyzx=-+-+-2222222()()()rrrdxxyyzx=-+-+-21dddvt-=1ddLfvc=?d1d222L1vd22图3.1其中，123,,ddd为一颗卫星不同时刻到接收机的伪距，dv为两个时刻之间的径向速度。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，选择不同时间段进行记录、求解、分析。【核心程序代码】%%Findsatellite'sposition----------------------------------------------%Restoresemi-majoraxisa=eph(prn).sqrtA*eph(prn).sqrtA;%Timecorrectiontk=check_t(time-eph(prn).t_oe);%Initialmeanmotionn0=sqrt(GM/a^3);%Meanmotionn=n0+eph(prn).deltan;%MeananomalyM=eph(prn).M_0+n*tk;%Reducemeananomalytobetween0and360degM=rem(M+2*gpsPi,2*gpsPi);%InitialguessofeccentricanomalyE=M;%---Iterativelycomputeeccentricanomaly----------------------------forii=1:10E_old=E;E=M+eph(prn).e*sin(E);dE=rem(E-E_old,2*gpsPi);ifabs(dE)1.e-12%Necessaryprecisionisreached,exitfromtheloopbreak;endend%Reduceeccentricanomalytobetween0and360degE=rem(E+2*gpsPi,2*gpsPi);%Computerelativisticcorrectiontermdtr=F*eph(prn).e*eph(prn).sqrtA*sin(E);%Ïà¶ÔÂÛÐÞÕýÏî%Calculatethetrueanomalynu=atan2(sqrt(1-eph(prn).e^2)*sin(E),cos(E)-eph(prn).e);%Computeanglephiphi=nu+eph(prn).omega;%Reducephitobetween0and360degphi=rem(phi,2*gpsPi);%Correctargumentoflatitudeu=phi+...eph(prn).C_uc*cos(2*phi)+...eph(prn).C_us*sin(2*phi);%Correctradiusr=a*(1-eph(prn).e*cos(E))+...eph(prn).C_rc*cos(2*phi)+...eph(prn).C_rs*sin(2*phi);%Correctinclinationi=eph(prn).i_0+