导航技术综合实验报告

导航技术综合实验报告实验号：实验九、实验十、实验十一、实验十二2012年4月23日星期一一、实验目的

理论：基本掌握解码星历表（数制转换）、求伪距观测量、求卫星位置和卫星钟差、最小二乘法求接收机位置的概念、方法和运算的过程。实践：通过对解码星历表（数制转换）、求伪距观测量、求卫星位置和卫星钟差、最小二乘法求接收机位置的程序仿真实际掌握方法。二、实验的方法和步骤实验方法：结合课本的理论知识使用matlab仿真软件完成实验。实验步骤：1、仔细阅读课本基本了解理论知识，掌握基本问题的解决方法与方式。2、认真研究老师所给实验题目，明确实验要求和目的，确立完成实验的基本思路。3、对照实验题目回顾书本理论知识，结合理论知识确定解决问题的方式方法。结算出理论值。4、绘制原理框图，方案流程图。5、确立仿真环境。6、绘制流程图‘7、编程调试，得到结果与理论推导结果做对比。三、实验结果实验9解码星历表（数制转换）一．函数模块解释函数功能：根据导航数据位信息，将从子帧1，子帧2和子帧3中得到星历数据转换成十进制格式，并求第一子帧的时间周TOW。一个子帧包含300数据位，导航信息共5个子帧。函数名称：function[eph,TOW]=ephemeris(bits,D30Star)1.实验要求：实验代码见实验9中的ephemeris.m文件，理解实现星历解算的代码，并把函数代码和流程图中的各个环节对应起来。获得输入的导航数据信息function[eph,TOW]=ephemeris(bits,D30Star)利用function形式，调用子函数，输入参数。检查导航数据的bits长度是否满足导航数据5个子帧的长度要求iflength(bits)<1500error('TheparameterBITSmustcontain1500bits!');end检查参数bits和D30Star是否是字符串if~ischar(bits)error('TheparameterBITSmustbeacharacterarray!');endif~ischar(D30Star)error('TheparameterD30Starmustbeachar!');end循环操作，一一对5个子帧进行解码。fori=1:5…end即外层循环。完成4.1提取第一个子帧的位；subframe=bits(300*(i-1)+1:300*i)4.2检查所提供的30个位的奇偶forj=1:10[subframe(30*(j-1)+1:30*j)]=...checkPhase(subframe(30*(j-1)+1:30*j),D30Star);D30Star=subframe(30*j);End4.3用subframeID=bin2dec(subframe(50:52))提取子帧号4.4基于子帧的ID号来进行解码，将位信息转换为十进制数。switchsubframeIDcase1eph.weekNumber=bin2dec(subframe(61:70))+1024;eph.accuracy=bin2dec(subframe(73:76));%用户测距精度eph.health=bin2dec(subframe(77:82));%卫星的健康状况eph.T_GD=twosComp2dec(subframe(197:204))*2^(-31);%群的延迟差别评估该数据为二进制补码，twosComp2dec将补码转换为对应的十进制数eph.IODC=bin2dec([subframe(83:84)subframe(197:204)]);%数据、时钟的发布号eph.t_oc=bin2dec(subframe(219:234))*2^4;%卫星时钟修正参数eph.a_f2=twosComp2dec(subframe(241:248))*2^(-55);%卫星时钟修正参数eph.a_f1=twosComp2dec(subframe(249:264))*2^(-43);%卫星时钟修正参数eph.a_f0=twosComp2dec(subframe(271:292))*2^(-31);%卫星时钟修正参数case2eph.IODE_sf2=bin2dec(subframe(61:68));%数据与星历发布号eph.C_rs=twosComp2dec(subframe(69:84))*2^(-5);%轨道半径的正弦调和修正项的幅度eph.deltan=...twosComp2dec(subframe(91:106))*2^(-43)*gpsPi;%计算值的平均移动误差eph.M_0=...twosComp2dec([subframe(107:114)subframe(121:144)])...*2^(-31)*gpsPi;eph.C_uc=twosComp2dec(subframe(151:166))*2^(-29);%纬度辐角的余弦调和修正项的幅度eph.e=...bin2dec([subframe(167:174)subframe(181:204)])...%离心率*2^(-33);eph.C_us=twosComp2dec(subframe(211:226))*2^(-29);%纬度辐角的正弦调和修正项的幅度eph.sqrtA=...bin2dec([subframe(227:234)subframe(241:264)])...%长半轴的平方根eph.t_oe=bin2dec(subframe(271:286))*2^4;%参考时间星历case3eph.C_ic=twosComp2dec(subframe(61:76))*2^(-29);%倾斜角的余弦调和修正项的幅度eph.omega_0=...twosComp2dec([subframe(77:84)subframe(91:114)])...%在每星期历元轨道平面上的升点经度*2^(-31)*gpsPi;eph.C_is=twosComp2dec(subframe(121:136))*2^(-29);%倾斜角的正弦调和修正项的幅度eph.i_0=...twosComp2dec([subframe(137:144)subframe(151:174)])...%参考时间的倾斜角*2^(-31)*gpsPi;eph.C_rc=twosComp2dec(subframe(181:196))*2^(-5);%轨道半径的正弦调和修正项的幅度eph.omega=...twosComp2dec([subframe(197:204)subframe(211:234)])...%近地点的辐角*2^(-31)*gpsPi;eph.omegaDot=twosComp2dec(subframe(241:264))*2^(-43)*gpsPi;%赤经的速率eph.IODE_sf3=bin2dec(subframe(271:278));%数据与星历发布号eph.iDot=twosComp2dec(subframe(279:292))*2^(-43)*gpsPi;%倾斜角的速率case4%Almanac,ionosphericmodel,UTCparameters.%SVhealth(PRN:25-32).%Notdecodedatthemoment.%年历，电离层模型，UTC参数。SV正常（PRN：25-32）。%此时不解码。case5%Itissubframe5这是子帧5%SValmanacandhealth(PRN:1-24).%Almanacreferenceweeknumberandtime.%Notdecodedatthemoment.%SV年历和正常值（PRN：1-24）。年历参考的周数和时间。%此刻不解码。end%switchsubframeID...转换子帧的IDend%forall5sub-frames...对所有的5子帧纠正第一子帧中的周时间TOW，当前子帧的时间要在当前的结果下减6sTOW=bin2dec(subframe(31:47))*6-30;%结论：通过对实验九的仿真，明确里在实际中结算星历的方法与步骤。对导航星历的求解，可以求解出相关信息。实验十求伪距观测量函数功能：根据跟踪的结果，以及计时的起始索引点，来求伪距观测值函数名称：function[pseudoranges]=calculatePseudoranges(absoluteSample,settings)输入相应的跟踪变量和每颗卫星的c/a码起点及初值实验要求：（1）补充完整实验代码，并设计程序，根据所给输入数据和结果验证代码的正确性。初始化接收机接收卫星信号的时间为无穷大根据通道c/a码记录结果，求每个卫星的传输时间然后求所有传输时间相对最先到达的卫星的时间差用相对时延差加上设置的一般到达时间作为该卫星到达接收机的时间Pseudoranges:用到达时间乘以光速获得卫星接收机的观测伪距，作为输出三．实验操作1.实验代码补充：function[pseudoranges]=calculatePseudoranges(absoluteSample,settings)%absoluteSample是导航C/A码起始索引点找到对应的数字点，除以每码采样点得到传输时间%absoluteSample是一个四维的数组，包含四颗星的数据%初始化各卫星传输时间为无穷大travelTime=inf(1,length(absoluteSample));%求每码采样数据点（samplesPerCode=16368）samplesPerCode=round(settings.samplingFreq/...(settings.codeFreqBasis/settings.codeLength));%补充代码%求每一个卫星的数据点对应的时间travelTime=absoluteSample/samplesPerCode;minimum=floor(min(travelTime));%求相对时延dt=travelTime-minimum;%求估计的传输时间=时延+设定的初值（到达时间）travelTime=dt+settings.startOffset;%求伪距pseudoranges=travelTime*(settings.c/1000);end2.实验数据说明：settings=msToProcess:1000numberOfChannels:4skipNumberOfBytes:16000000fileName:'usbdata.bin'dataType:'int8'IF:4092000settings=msToProcess:1000numberOfChannels:4skipNumberOfBytes:16000000fileName:'usbdata.bin'dataType:'int8'IF:4092000samplingFreq:16368000pseudoranges=1.0e+007*2.09262.15842.17912.0647结论：测量伪距的方法是采用相对方法测量伪距。选择最靠前一个子帧1的起始点作为参考点，先在不考虑所有卫星的时钟修正的问题时，可以认为所有卫星对应子帧的起始点是在同一时刻发射的，而不同卫星对应子帧的起始点是在不同时刻接收到的，其它卫星相对参考点的时间差求出来，那么这个时间差就可以认为是不同卫星到接收机的信号传输时间相对差值。在相对时延的基础上加上给定的参考点的卫星到接收的传输时间就是每个卫星的传输时间。程序对应的伪距提取的主要过程：（1）在对应伪距提取时刻，找到各个通道子帧的起始位置，根据bit同步，找到子帧起点对应于跟踪结果中的起点（精度为C／A码周期，即1ms)；（2）寻找对应的起点的数据点，然后用其除以C／A码周期的采样点数，得到相对应的传输时间；（3）求所有通道中最小的传输时间，把它作为参考点；（4）求各个卫星相对参考点的时间差∆ti，然后再时间差的基础上加上给定卫星到接收机传输时间的初值(ti=t0+∆ti)，那么最后ti是对应第i颗卫星的传输时间，乘以光速就获得该卫星的伪距观测值。每次提取伪距，耗时是500ms，第n次循环提取时都要把前面用去的解算时间加上去，所以程序的最后一步，每次循环后，在传输时间上加上了一个解算周期时间。实验十一、求卫星位置和卫星钟差一．函数模块解释函数功能：利用测量的伪距求卫星在ECEF坐标系中的坐标位置和卫星钟差，其中的求解原理见GPS简化版文件。函数名称：function[satPositions,satClkCorr]=satpos(transmitTime,prnList,...eph,settings)实验要求：（1）根据原理将上述代码的原理在书上找到相应的公式，并加以注释和理解；（2）根据satpos.m把实验代码与流程图联系起来理解，并结合书上原理，理解求卫星位置的函数过程。u=phi+...eph(prn).C_uc*cos(2*phi)+...eph(prn).C_us*sin(2*phi);所对应公式：%校正后的卫星地心相径r=a*(1-eph(prn).e*cos(E))+...eph(prn).C_rc*cos(2*phi)+...eph(prn).C_rs*sin(2*phi);所对应公式：%校正后的轨道倾角i=eph(prn).i_0+eph(prn).iDot*tk+...eph(prn).C_ic*cos(2*phi)+...eph(prn).C_is*sin(2*phi);所对应公式：% 计算升交点的经度Omega=eph(prn).omega_0+(eph(prn).omegaDot-Omegae_dot)*tk-...Omegae_dot*eph(prn).t_oe;Omega=rem(Omega+2*gpsPi,2*gpsPi);所对应公式：%求卫星ECEF坐标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);所对应公式：%本块作用：求该卫星的时钟修正项satClkCorr(satNr)=(eph(prn).a_f2*dt+eph(prn).a_f1)*dt+...eph(prn).a_f0-...eph(prn).T_GD+dtr;所对应公式：结论：卫星的无摄运动对应平滑的椭圆轨道，一般可通过一组适宜的参数来描述。但是，这组参数的选择并不是惟一的。在定位的过程中，接收机的位置一般是用ECEF坐标描述的，为了后续计算卫星和接收机之间的距离计算，就要统一两者的坐标参考系。一般是把卫星由其轨道旋转到ECEF参照系。实验十二、最小二乘法求接收机位置一．函数模块解释函数功能：根据卫星的估计位置和钟差，以及测量的伪距来求接收机位置和时钟误差，用的算法是最小二乘法，原理见书。函数名称：function[pos,el,az,dop]=leastSquarePos(satpos,obs,settings)实验要求：（1）结合流程图和求解原理补充完整函数代码，然后设计实验，将实验文件中所给的数据载入作为上输入；（2）验证实验结果是否和文件说明中的结果是一致的。2.实验数据说明：（1）伪距.mat和卫星坐标.mat是用于测试的卫星文件和相应的卫星对应的观测伪距（2）结果在实验12的说明文档中。3.实验要求：（1）结合流程图和求解原理补充完整函数代码，然后设计实验，将实验文件中所给的数据载入作为上输入；3.实验要求：（1）结合流程图和求解原理补充完整函数代码，然后设计实验，将实验文件中所给的数据载入作为上输入；（2）验证实验结果是否和文件说明中的结果是一致的。function[pos,el,az,dop]=leastSquarePos(satpos,obs,settings)%用最小二乘的迭代算法求接收机的位置，原理参照课本%迭代次数设置nmbOfIterations=7;dtr=pi/180;pos=zeros(4,1);X=satpos;nmbOfSatellites=size(satpos,2);A=zeros(nmbOfSatellites,4);omc=zeros(nmbOfSatellites,1);az=zeros(1,nmbOfSatellites);el=az;%开始迭代运算foriter=1:nmbOfIterationsfori=1:nmbOfSatellitesifiter==1%InitializevariablesatthefirstiterationRot_X=X(:,i);trop=2;else%Updateequationsrho2=(X(1,i)-pos(1))^2+(X(2,i)-pos(2))^2+...(X(3,i)-pos(3))^2;traveltime=sqrt(rho2)/settings.c;%补充code%Correctsatelliteposition(dotoearthrotation)Rot_X=e_r_corr(traveltime,X(:,i));%修正第i颗卫星的位置%求仰角，方位角（有函数代码）[az(i),el(i),dist]=topocent(pos(1:3,:),Rot_X-pos(1:3,:));if(settings.useTropCorr==1)%Calculatetroposphericcorrectiontrop=tropo(sin(el(i)*dtr),...0.0,1013.0,293.0,50.0,0.0,0.0,0.0);else%Donotcalculateorapplythetroposphericcorrectionstrop=0;endend%ifiter==1......else%伪距的修正补充omc(i)=(obs(i)-norm(Rot_X-pos(1:3),'fro')-pos(4)-trop);%建立A设计矩阵补充A(i,:)=[(-(Rot_X(1)-pos(1)))/obs(i)...(-(Rot_X(2)-pos(2)))/obs(i)...(-(Rot_X(3)-pos(3)))/obs(i)...1];end%fori=1:nmbOfSatellites%判断A的稚是否为4，否则方程无解补充ifrank(A)~=4pos=zeros(1,4);returnendx=A\omc;pos=pos+x;%位置迭代变换更新--补充endpos=pos';%计算精度因子--见原理ifnargout==4dop=zeros(1,5);Q=inv(A'*A);dop(1)=sqrt(trace(Q));%GDOPdop(2)=sqrt(Q(1,1)+Q(2,2)+Q(3,3));%PDOPdop(3)=sqrt(Q(1,1)+Q(2,2));%HDOPdop(4)=sqrt(Q(3,3));%VDOPdop(5)=sqrt(Q(4,4));%TDOPend（2）验证实验结果是否和文件说明中的结果