导航技术综合实验报告
实验号：实验九、实验十、实验十一、实验十二
20##年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文件，理解实现星历解算的代码，并把函数代码和流程图中的各个环节对应起来。
1、 获得输入的导航数据信息
function [eph, TOW] = ephemeris(bits, D30Star)
利用function形式，调用子函数，输入参数。
2、 检查导航数据的bits长度是否满足导航数据5个子帧的长度要求
if length(bits) < 1500
error('The parameter BITS must contain 1500 bits!');
end
3、 检查参数bits和D30Star是否是字符串
if ~ischar(bits)
error('The parameter BITS must be a character array!');
end
if ~ischar(D30Star)
error('The parameter D30Star must be a char!');
end
4、 循环操作，一一对5个子帧进行解码。
for i = 1:5
…end即外层循环。完成
4.1提取第一个子帧的位；
subframe = bits(300*(i-1)+1 : 300*i)
4.2检查所提供的30个位的奇偶
for j = 1:10
[subframe(30*(j-1)+1 : 30*j)] = ...
checkPhase(subframe(30*(j-1)+1 : 30*j), D30Star);
D30Star = subframe(30*j);
End
4.3用subframeID = bin2dec(subframe(50:52))提取子帧号
4.4基于子帧的ID号来进行解码，将位信息转换为十进制数。
switch subframeID
case 1
eph.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);%卫星时钟修正参数
case 2
eph.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;%参考时间星历
case 3
eph.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;%倾斜角的速率
case 4
% Almanac, ionospheric model, UTC parameters.
% SV health (PRN: 25-32).
% Not decoded at the moment.
%年历，电离层模型，UTC参数。SV正常（PRN：25-32）。
%此时不解码。
case 5 %--- It is subframe 5 这是子帧5-------------------------------------
% SV almanac and health (PRN: 1-24).
% Almanac reference week number and time.
% Not decoded at the moment.
%SV年历和正常值（PRN：1-24）。年历参考的周数和时间。
%此刻不解码。
end % switch subframeID ...转换子帧的ID
end % for all 5 sub-frames ...对所有的5子帧
5、 纠正第一子帧中的周时间TOW，当前子帧的时间要在当前的结果下减6s
TOW = bin2dec(subframe(31:47)) * 6 - 30;%
结论：通过对实验九的仿真，明确里在实际中结算星历的方法与步骤。对导航星历的求解，可以求解出相关信息。
实验十 求伪距观测量
函数功能： 根据跟踪的结果，以及计时的起始索引点，来求伪距观测值
函数名称：function [pseudoranges] = calculatePseudoranges(absoluteSample,settings) 输入相应的跟踪变量和每颗卫星的c/a码起点及初值
实验要求：
（1）补充完整实验代码，并设计程序，根据所给输入数据和结果验证代码的正确性。
1、 初始化接收机接收卫星信号的时间为无穷大
2、 根据通道c/a码记录结果，求每个卫星的传输时间
3、 然后求所有传输时间相对最先到达的卫星的时间差
4、 用相对时延差加上设置的一般到达时间作为该卫星到达接收机的时间
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);
end
2.实验数据说明：
settings =
msToProcess: 1000
numberOfChannels: 4
skipNumberOfBytes: 16000000
fileName: 'usbdata.bin'
dataType: 'int8'
IF: 4092000
settings =
msToProcess: 1000
numberOfChannels: 4
skipNumberOfBytes: 16000000
fileName: 'usbdata.bin'
dataType: 'int8'
IF: 4092000
samplingFreq: 16368000
pseudoranges =
1.0e+007 *
2.0926 2.1584 2.1791 2.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;
% 开始迭代运算
for iter = 1:nmbOfIterations
for i = 1:nmbOfSatellites
if iter == 1
%--- Initialize variables at the first iteration --------------
Rot_X = X(:, i);
trop = 2;
else
%--- Update equations -----------------------------------------
rho2 = (X(1, i) - pos(1))^2 + (X(2, i) - pos(2))^2 + ...
(X(3, i) - pos(3))^2;
traveltime = sqrt(rho2) / settings.c ;%补充code
%--- Correct satellite position (do to earth rotation) --------
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)
%--- Calculate tropospheric correction --------------------
trop = tropo(sin(el(i) * dtr), ...
0.0, 1013.0, 293.0, 50.0, 0.0, 0.0, 0.0);
else
% Do not calculate or apply the tropospheric corrections
trop = 0;
end
end % if iter == 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 % for i = 1:nmbOfSatellites
%判断A的稚是否为4，否则方程无解 补充
if rank(A) ~= 4
pos = zeros(1, 4);
return
end
x = A \ omc;
pos = pos + x;
% 位置迭代变换更新--补充
end
pos = pos';
%计算精度因子--见原理
if nargout == 4
dop = zeros(1, 5);
Q = inv(A'*A);
dop(1) = sqrt(trace(Q)); % GDOP
dop(2) = sqrt(Q(1,1) + Q(2,2) + Q(3,3)); % PDOP
dop(3) = sqrt(Q(1,1) + Q(2,2)); % HDOP
dop(4) = sqrt(Q(3,3)); % VDOP
dop(5) = sqrt(Q(4,4)); % TDOP
end
（2）验证实验结果是否和文件说明中的结果是一致的。
settings =
msToProcess: 1000
numberOfChannels: 4
skipNumberOfBytes: 16000000
fileName: 'usbdata.bin'
dataType: 'int8'
IF: 4092000
settings =
msToProcess: 1000
numberOfChannels: 4
skipNumberOfBytes: 16000000
fileName: 'usbdata.bin'
dataType: 'int8'
IF: 4092000
samplingFreq: 16368000
pos = 1.0e+003 *
6.2664 0.0302 0.0287 -0.0895
el =
-14.6876 -56.7839 -33.1457 -46.8281 25.2618
az =
55.4842 312.1270 296.3907 308.4924 306.6239
dop =
23.7428 21.0671 20.1687 6.0867 10.9499
结论：最小二乘法求接收机位置得到用户的坐标位置和用户时钟偏差，由最终的矩阵A计算误差放大因子DOP。
个人总结：
技术方面:
1) 在解码星历表（数制转换）、求伪距观测量、求卫星位置和卫星钟差实验中能较好的理解和对应实验内容，可以通过对理论知识的研究可以彻底明白其方法与步骤。
2) 对于最小二乘法求接收机位置这个实验中，我想我还是欠缺一定的理解，一方面是对理论知识的掌握不够透彻，另一方面我想对利用最小二乘法求接收机位置的方法和步骤没有很好的研究。
3) 我想通过这次实验仿真让我对解码星历表（数制转换）、求伪距观测量、求卫星位置和卫星钟差、最小二乘法求接收机位置等问题有更深的理解。
4) 当然对于利用matlab仿真导航技术原理的程序步骤有所掌握。
个人收获：
1） 通过这次对解码星历表（数制转换）、求伪距观测量、求卫星位置和卫星钟差、最小二乘法求接收机位置的实验我大致了解了GPS接收的部分工作原理，也让我更加清楚如何实现接收，解算，编译导航电文，叫我不仅仅巩固了理论知识，更让我感性的认识到了关于GPS接收的有关知识。
2） 当然在此次实验中我想我还有以下不足，一、对理论知识的掌握不够娴熟，二、对实验具体细节问题没有很好的认真研究，三、理论联系实际做的不是很好。最后我想通过这次对导航技术综合实验的仿真我还是有很大的收获，尤其是对接收机部分有了更加深入的了解。