单点定位程序说明(优选.)
卫星导航定位算法与程序设计_单点定位程序流程
GPS 单点定位程序流程、计算流程读取RINEX O 文件,读取一个历元观测值epoch数据预处理根据epoch 中的卫星号和历元时刻T R 在ephlst 查找相应的卫星星历, 准则 T R TOE 3600.0S 。
对论改正;d ）对卫星位置X Si T Si进行地球自转改正，得到e ）根据X W i T Si 和测站概略位置X r 计算卫星和测站的几何距离 R Sf ）根据几何距离R Si求信号传播时间1SR S/c 。
1、读取RINEX N 文件，将所有星历放到一个列表（数组） ep hist 中。
2、 3、 4、 程序初始化, 置测站概略位置为 X r ，接收机钟差初值dt r 。
X OX r cdt r OXO YO第一次迭代，取X O ZOcdt5、选择epoch 中一颗卫星S i 观测值，设其伪距为S& 计算卫星S i 的信号发射的概略时刻T s方法如下: a ）卫星S i 的信号传播时间:SS iO/c dt r dt Si ;dt Si 为卫星钟差，需要进行相b ）卫星S i 的信号发射时刻: T ST RS iOc ）卫星S i 在T Si 时刻的位置 X ST SiX SYZ SiTTSiT S ；17、 8、 10、 11、 12、g ）如果107，则退出迭代。
T ,Sih ）否则S 1Si，回带到b ）进行迭代。
求卫星s 方向余弦b Si^R ^,b 1SY Y S求卫星s 在观测方程式中的余数项:其中:T RlSiSiR Si c dt Sid trop Si——卫星S 的伪距观测值； 卫星S i 到测站的几何距离； —以米表示的卫星S i 的钟差； 对流层延迟改正量，单位米，用简化的R Si — c dt Sdtro pd iono - D RTCM1S即为卫星信号发射时刻。
Z Z S,b 3S 1d iono D R TCMhop field 模型计算;电离层延迟改正量，单位米，采用无电离层伪距组合观测值时，此项为 —对伪距的差分改正值，此处为 0;选择epoch 中下一颗卫星S j 观测值，设其伪距为 重复第6—9步，计算每颗卫星的系数和余数项 将所有卫星的系数组成误差方程， x,y,z,cdt rS j为未知参数进行求解，形式应该是：AX Li 0,1,L 0；b20b： MbjTx y z cdtl S 0 l S 1 LlS i,svnum求解法方程刃 A TPA 1 A TPL，求出定位结果13、14、15、16、XYX i ZcdtrX0Y0Z0Cdt r0xxcdt与X o进行比较，判断位置差值,a）如果各分量差值＞0.001m，则令X ob）若小于则退出迭代。
绝对单点定位的计算流程
绝对单点定位的计算流程下载温馨提示:该文档是我店铺精心编制而成，希望大家下载以后，能够帮助大家解决实际的问题。
文档下载后可定制随意修改，请根据实际需要进行相应的调整和使用，谢谢!并且，本店铺为大家提供各种各样类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，如想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by theeditor. I hope that after you download them,they can help yousolve practical problems. The document can be customized andmodified after downloading,please adjust and use it according toactual needs, thank you!In addition, our shop provides you with various types ofpractical materials,such as educational essays, diaryappreciation,sentence excerpts,ancient poems,classic articles,topic composition,work summary,word parsing,copy excerpts,other materials and so on,want to know different data formats andwriting methods,please pay attention!绝对单点定位的计算流程一、准备工作阶段。
在进行绝对单点定位计算之前，需做好充分准备。
GPS单点定位实验报告
GPS原理与应用实验题目：GPS单点定位专业：测绘工程班级：12-01学号：2012212600姓名：王威指导教师：陶庭叶时间：2014.11目录一、实验目的 (3)二、实验原理 (3)三、实验内容 (3)四、实验效果图 (9)五、实验总结 (9)一．实验目的1.深入了解单点定位的计算过程；2.加强单点定位基本公式和误差方程式，法线方程式的记忆；3.通过上机调试程序加强动手能力的培养。
二．实验原理一个接收机接受三个火三个以上卫星信号，得出卫星坐标和伪距，利用间接平差计算接收机的坐标。
三．实验内容1.程序流程图2、实验数据3、实验程序代码Private Sub Command1_Click()CommonDialog1.Filter = "TXT files|*.txt|" CommonDialog1.FilterIndex = 1CommonDialog1.ShowOpenOpen monDialog1.FileName For Input As #1 Do While Not EOF(1)Line Input #1, Texttextbuff = textbuff + Text + vbCrLfLoopClose #1kk = MSFlexGrid1.Rows - 1Dim aReDim a(kk - 1)a = Split(textbuff, vbCrLf)For j = 1 To kkMSFlexGrid1.TextMatrix(j, i) = a(j - 1 + 5 * (i - 1)) Next iNext jFor k = 1 To kkMSFlexGrid1.TextMatrix(k, 0) = "第" & k & "个点" Next kMSFlexGrid1.TextMatrix(0, 1) = "X"MSFlexGrid1.TextMatrix(0, 2) = "Y"MSFlexGrid1.TextMatrix(0, 3) = "Z"MSFlexGrid1.TextMatrix(0, 4) = "伪距" MSFlexGrid1.TextMatrix(0, 5) = "钟差"End SubPrivate Sub Command2_Click()kk = MSFlexGrid1.Rows - 1X0 = 0: Y0 = 0: Z0 = 0c = 299792458Dim a()ReDim a(kk - 1, 3)Dim ll()ReDim ll(kk - 1, 0)For i = 1 To kkl = (MSFlexGrid1.TextMatrix(i, 1) - X0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)m = (MSFlexGrid1.TextMatrix(i, 2) - Y0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)n = (MSFlexGrid1.TextMatrix(i, 3) - Z0) / Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2)a(i - 1, 0) = la(i - 1, 1) = ma(i - 1, 2) = na(i - 1, 3) = -1lk = MSFlexGrid1.TextMatrix(i, 4) - Sqr((MSFlexGrid1.TextMatrix(i, 1) - X0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 2) - Y0) ^ 2 + (MSFlexGrid1.TextMatrix(i, 3) - Z0) ^ 2) + c * MSFlexGrid1.TextMatrix(i, 5)ll(i - 1, 0) = lkNext igzs = xc(qiuni(xc(zz(a), a)), xc(zz(a), ll))X0 = X0 - gzs(0, 0)Y0 = Y0 - gzs(1, 0)Z0 = Z0 - gzs(2, 0)j = j + 1Next iiText2.Text = "X=" & X0 & vbCrLf & vbCrLf & "Y=" & Y0 & vbCrLf & vbCrLf & "Z=" & Z0V = jian(ll, xc(a, gzs))zjl = xc(zz(V), V)σ0 = Sqr(zjl(0, 0)) / (kk - 3)Qx = qiuni(xc(zz(a), a))Text3.Text = "σX=" & σ0 * Sqr(Qx(0, 0)) & vbCrLf & vbCrLf & "σY=" & σ0 * Sqr(Qx(1, 1)) & vbCrLf & vbCrLf & "σZ=" & σ0 * Sqr(Qx(2, 2))End SubPrivate Sub Form_Load()MSFlexGrid1.ColWidth(1) = 1300MSFlexGrid1.ColWidth(2) = 1300MSFlexGrid1.ColWidth(3) = 1300MSFlexGrid1.ColWidth(4) = 1300Text2.Text = ""Text3.Text = ""End SubPublic Function jian(m, n)Dim i, j As IntegerIf UBound(m, 1) <> UBound(n, 1) Or UBound(m, 2) <> UBound(n, 2) Then MsgBox ("请确认输入数组是否可以相减！")ElseDim c()ReDim c(UBound(m, 1), UBound(n, 2))For i = 0 To UBound(c, 1)For j = 0 To UBound(c, 2)c(i, j) = m(i, j) - n(i, j)Next jNext ijian = cEnd IfEnd Function'矩阵的转置Public Function zz(a)Dim i As Integer, j As Integer, t As Integer, b()If UBound(a, 1) = UBound(a, 2) ThenFor i = 0 To UBound(a, 1)For j = 0 To UBound(a, 2)t = a(i, j)a(i, j) = a(j, i)a(j, i) = tEnd IfNext jNext izz = aElseReDim b(UBound(a, 2), UBound(a, 1))For i = 0 To UBound(a, 2)For j = 0 To UBound(a, 1)b(i, j) = a(j, i)Next jNext izz = bEnd IfEnd Function'两矩阵相乘Public Function xc(a, b)Dim i As Integer, j As Integer, k As Integer If UBound(a, 2) <> UBound(b, 1) ThenMsgBox ("这两个矩阵不能够相乘") Exit FunctionEnd IfReDim sd(UBound(a, 1), UBound(b, 2)) For i = 0 To UBound(a, 1)For j = 0 To UBound(b, 2)For k = 0 To UBound(b, 1)sd(i, j) = sd(i, j) + a(i, k) * b(k, j)Next kNext jNext ixc = sdEnd FunctionPublic Function qiuni(a)Dim c, m%, n%, p#, l%, i%, j%, ab#m = UBound(a, 1)n = UBound(a, 2)If m <> n ThenMsgBox ("该矩阵不可逆！！！")Exit FunctionEnd IfReDim c(m, 2 * n + 1)For i = 0 To mFor j = 0 To nc(i, j) = a(i, j)Next jNext iFor i = 0 To mFor j = m + 1 To 2 * m + 1 c(i, j) = 0Next jNext ii = 0For j = m + 1 To 2 * m + 1 c(i, j) = 1i = i + 1Next jFor k = 0 To nIf c(k, k) = 0 ThenFor i = k + 1 To nIf c(i, k) <> 0 ThenGoTo thisEnd IfNext iIf i = n + 1 ThenMsgBox ("该矩阵不可逆!!!") Exit FunctionEnd Ifthis:For j = 0 To 2 * m + 1p = c(k, j)c(k, j) = c(i, j)c(i, j) = pNext jEnd Ifab = 1# / c(k, k)For j = 0 To 2 * m + 1c(k, j) = c(k, j) * abNext jFor i = 0 To nIf i <> k ThenFor j = 0 To 2 * m + 1If j <> k Thenc(i, j) = c(i, j) - c(i, k) * c(k, j)End IfNext jc(i, k) = 0End IfNext iNext kFor i = 0 To mFor j = 0 To ma(i, j) = c(i, j + n + 1) a(i, j) = Round(a(i, j), 4) Next jNext iqiuni = aEnd Function四．实验结果图五．实验总结此次实验让我深入了解单点定位的计算过程，加强了对单点定位基本公式和误差方程式，法线方程式的记忆。
RTKLIB源码解析（一）——单点定位（pntpos.c）
RTKLIB源码解析（⼀）——单点定位（pntpos.c）RTKLIB源码解析（⼀）——单点定位(pntpos.c)标签： GNSS RTKLIB 单点定位前段时间⼀直忙着写毕业论⽂，所以也没有太多时间来阅读 RTKLIB源码，最近好⽍把 pntpos中的相关代码看了⼀遍，知道了 RTKLIB是如何实现单点伪距定位的。
这⾥把每⼀个函数都做成了⼩卡⽚的形式，每个函数⼤都包含函数签名、所在⽂件、功能说明、参数说明、处理过程、注意事项和我的疑惑这⼏个部分，介绍了阅读代码时我⾃⼰的看法和疑惑。
所以希望诸位看官能帮忙解答我的疑惑，与我交流，也希望能帮助后来也有需要阅读 RTKLIB源码的⼈，给他们多提供⼀份思路。
总⽽⾔之，既为⼈，也为⼰。
这份⽂档是使⽤ Cmd Markdown完成的，在作业部落上其格式显式的⾮常完整，但是在博客园中⽬录、代码块和流程图似乎都没有显⽰出来，所以这⾥也贴上本⽂在作业部落上的链接，对格式“零容忍”的同学请移步那⾥。
⽬录pntposint pntpos (const obsd_t *obs, int n, const nav_t *nav, const prcopt_t *opt, sol_t *sol,double *azel, ssat_t *ssat, char *msg)所在⽂件：pntpos.c功能说明：依靠多普勒频移测量值和伪距来进⾏单点定位，给出接收机的位置、速度和钟差参数说明：函数参数，8个：obsd_t *obs I observation dataint n I number of observation datanav_t *nav I navigation dataprcopt_t *opt I processing optionssol_t *sol IO solutiondouble *azel IO azimuth/elevation angle (rad) (NULL: no output)ssat_t *ssat IO satellite status (NULL: no output)char *msg O error message for error exit返回类型:int O (1:ok,0:error)调⽤关系：如⽆特别说明，本⽂所出现的流程图中，纵向代表时间上的先后调⽤顺序，横向代表层次上的逐级调⽤顺序。
利用MATLAB编制的GPS单点定位程序
function t=TimetoJD(Y,M,D,h,f,s)if(Y>=80)Y=Y+1900;elseY=Y+2000;endif M<=2Y=Y-1;M=M+12;endJD=fix(365.25*Y)+fix(30.6001*(M+1))+D+h/24+f/1440+s/24/3600+1720981.5;t=JD-2444244.5;function [head,obs]=ReadObsData%读接收机观测数据文件%HeadODat :a structure stores header information if o-file% .ApproXYZ[3]; //approximate coordinate% .interval; //intervals of two adjacent epochs% .SiteName; //point name% .Ant_H; //antenna height% .Ant_E; //east offset of the antenna center% .Ant_N; //north offset of then antenna center% .TimeOB; //first epoch time in modifuied Julian time% .TimeOE; //last epoch time in modifuied Julian time% .SumOType; //number of types of observables% .SumOO[SumOType]; //type of observables 0-L1,1-L2,2-C1,3-P1,4-P2,5-D1,6-D2,%ObsODat :a structure stores observables by one and one epoch% .TimeOEpp[2]; //reciever time of epoch% .SatSum; //number of satellites% .SatCode[SatSum]; //satellites' PRN% .Obs_FreL1[SatSum];% .Obs_FreL2[SatSum];% .Obs_RangeC1[SatSum];% .Obs_RangeP1[SatSum];% .Obs_RangeP2[SatSum];global HeadODat;global ObsODat;[fname,fpath]=uigetfile('*.*','选择一个O文件');O_filename=strcat(fpath,fname);fid1=fopen(O_filename,'rt');if (fid1==-1)msgbox('file invalide','warning','warn');return;end%将文件头数据存入结构体HeadODat中t=0;while(t<100)s=fgets(fid1);t=t+1;L=size(s,2);if L<81s(L+1:81)=' ';endinstrS=s(61:81);%测站点近似坐标if strncmp(instrS,'APPROX POSITION XYZ',19)HeadODat.ApproXYZ=zeros(1,3);HeadODat.ApproXYZ(1,1)=str2num(s(1:14));HeadODat.ApproXYZ(1,2)=str2num(s(15:28));HeadODat.ApproXYZ(1,3)=str2num(s(29:42));%历元间隔;elseif strncmp(instrS,'INTERVAL',8)HeadODat.interval=str2num(s(5:11));%测站点号elseif strncmp(instrS,'MARKER NAME',11)HeadODat.SiteName=s(1:4)%天线中心改正elseif strncmp(instrS,'ANTENNA: DELTA H/E/N',20)HeadODat.Ant_H=str2num(s(1:14));HeadODat.Ant_E=str2num(s(15:28));HeadODat.Ant_N=str2num(s(29:42));%第一个历元时间elseif strncmp(instrS,'TIME OF FIRST OBS',17)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOB=TimetoJD(year,month,day,hour,minute,second);%最后一个历元时间elseif strncmp(instrS,'TIME OF LAST OBS',16)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOE=TimetoJD(year,month,day,hour,minute,second);%观测值类型elseif strncmp(instrS,'# / TYPES OF OBSERV',19)HeadODat.SumOType=str2num(s(1:6));HeadODat.SumOO=ones(1,HeadODat.SumOType)*-1;for k=1:HeadODat.SumOTypef=s(k*6+5:k*6+6);if strcmp(f,'L1')HeadODat.SumOO(1,k)=0;elseif strcmp(f,'L2')HeadODat.SumOO(1,k)=1;elseif strcmp(f,'C1')HeadODat.SumOO(1,k)=2;elseif strcmp(f,'P1')HeadODat.SumOO(1,k)=3;elseif strcmp(f,'P2')HeadODat.SumOO(1,k)=4;elseif strcmp(f,'D1')HeadODat.SumOO(1,k)=5;elseif strcmp(f,'D2')HeadODat.SumOO(1,k)=6;endend%头文件结束elseif strncmp(instrS,'END OF HEADER',13)break;elsecontinue;endend%观测数据结构体%观测数据结构t=0;while ~feof(fid1)%每个历元的第一行数据，时间和观测到的卫星号s=fgets(fid1);t=t+1;year=str2num(s(1:3));month=str2num(s(4:6));day=str2num(s(7:9));hour=str2num(s(10:12));minute=str2num(s(13:15));second=str2num(s(16:26));%历元时间ObsODat(t).TimeOEp=[year,month,day,hour,minute,second];ObsODat(t).TimeOEpp=TimetoJD(year,month,day,hour,minute,second);%该历元观测到的卫星数ObsODat(t).SatSum=str2num(s(30:32));%该历元观测到的卫星号ObsODat(t).SatCode=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL2=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeC1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP2=zeros(1,ObsODat(t).SatSum);for k=1:ObsODat(t).SatSumf=s(31+k*3:32+k*3);ObsODat(t).SatCode(1,k)=str2num(f);end%每个历元的观测数据，按卫星号先后顺序分行存for k=1:ObsODat(t).SatSums=fgets(fid1);%判断一个卫星的观测数据是否分两行存储，如果为两行，则再读入一行if HeadODat.SumOType>5sg=fgets(fid1);s=strcat(s,sg);endL=size(s,2);%补充数据长度if L<HeadODat.SumOType*16s(L+1:HeadODat.SumOType*16)=' ';end%对观测数据判断其类型，并存储到相应的数组中for j=1:HeadODat.SumOTypestemp=s(j*16-15:j*16);stemp=deblank(stemp);if isempty(stemp)continue;endstempNum=str2num(stemp);stempLength=size(stempNum,2);if stempLength>1stempNum=stempNum(1,1);endif HeadODat.SumOO(1,j)==0ObsODat(t).Obs_FreL1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==1ObsODat(t).Obs_FreL2(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==2ObsODat(t).Obs_RangeC1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==3ObsODat(t).Obs_RangeP1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==4ObsODat(t).Obs_RangeP2(1,k)=stempNum;elsecontinue;endend%完成一个卫星的所有观测数据存储end%完成一个历元的数据存储end%完成所有历元的数据存储head=HeadODat;obs=ObsODat;returnfunction EphDat=ReadGpsDataglobal EphDatEphDat=struct;[filename,pathname]=uigetfile('*.**N','读取GPS广播星历文件'); fid1=fopen(strcat(pathname,filename),'rt');if(fid1==-1)msgbox('Input File or Path is not correct','warning ','warn');return;endwhile(1)temp=fgetl(fid1);header=findstr(temp,'END OF HEADER');if(~isempty(header))break;endendi=1;while(1)temp=fgetl(fid1);break;endEphDat(i).PRN=str2double(temp(1:2));year=str2double(temp(3:5));EphDat(i).toc=TimetoJD(year,str2double(temp(6:8)),str2double(temp(9:11)),str2double(temp(12: 14)),str2double(temp(15:17)),str2double(temp(18:22)));EphDat(i).a0=str2num(temp(23:41));EphDat(i).a1=str2num(temp(42:60));EphDat(i).a2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).idoe=str2num(temp(4:22));EphDat(i).Crs=str2num(temp(23:41));EphDat(i).dn=str2num(temp(42:60));EphDat(i).M0=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).Cuc=str2num(temp(4:22));EphDat(i).e=str2num(temp(23:41));EphDat(i).Cus=str2num(temp(42:60));EphDat(i).sqa=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).toe=str2num(temp(4:22));EphDat(i).Cic=str2num(temp(23:41));EphDat(i).omg0=str2num(temp(42:60));EphDat(i).Cis=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).i0=str2num(temp(4:22));EphDat(i).Crc=str2num(temp(23:41));EphDat(i).w=str2num(temp(42:60));EphDat(i).omg=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).iDot=str2num(temp(4:22));EphDat(i).cflgl2=str2num(temp(23:41));EphDat(i).weekno=str2num(temp(42:60));EphDat(i).pflgl2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).svacc=str2num(temp(4:22));EphDat(i).svhlth=str2num(temp(23:41));EphDat(i).tgd=str2num(temp(42:60));EphDat(i).aodc=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).ttm=str2num(temp(4:22));if(i~=1) %删除重复数据if (EphDat(i-1).PRN~=EphDat(i).PRN)break;elseif abs(EphDat(i-1).toc-EphDat(i).toc)<0.001i=i - 1;endendendi=i + 1;endfunction DDDWformat longclear allticglobal HeadODat;global ObsODat;global EphDat;%先读N文件，再读O文件EphDat=ReadGpsData;[HeadODat,ObsODat]=ReadObsData;%多个历元加权平均求测站点坐标N = size(EphDat,2);c=299792458;epochNUM=size(ObsODat,2);%观测数据的个数Xk0=HeadODat.ApproXYZ(1,1); %测站点的近似坐标Yk0=HeadODat.ApproXYZ(1,2);Zk0=HeadODat.ApproXYZ(1,3);for k=1:epochNUMtr=ObsODat(k).TimeOEp;%历元接收数据时间Snum=ObsODat(1,k).SatSum; %该历元观测到的卫星数if Snum<4break; %去除无解的历元endCode=ObsODat(k).SatCode;%该历元观测到的卫星号组R=ObsODat(k).Obs_RangeC1 ; %取出C1观测值,列向量%卫星发射时间的迭代计算for j=1:Snumif R(j) == 0continue;endt=TimetoJD(tr(1),tr(2),tr(3),tr(4),tr(5),tr(6));t2 = mod(t,7)*24*3600;%gps second%卫星钟差dd=-1;for i=1:Nif(EphDat(i).PRN==Code(j)&abs(t-EphDat(i).toc)<0.0417*2)%读入时间相近的星历文件a0= EphDat(i).a0;a1=EphDat(i).a1;a2= EphDat(i).a2;toe=EphDat(i).toe;dt=a0+(t2-toe)*a1+(t2-toe)^2*a2;%卫星钟差dd=1;break;endendif dd==-1msgbox('没有相关星历文件');return;endtt = tr;ts = tr(6) - R(j)/c;%用秒进行迭代sign = 1;while(sign>1E-8)tt(6) = ts;[Xs1,Ys1,Zs1]= CalPos(Code(j),tt);ss1 = sqrt((Xk0-Xs1)^2+(Yk0-Ys1)^2+(Zk0-Zs1)*(Zk0-Zs1));%卫星位置加地球自转改正qe=0.000072921151467; %地球自转角速度delt=ss1/c;Rotation=[cos(qe*delt),sin(qe*delt),0;-sin(qe*delt),cos(qe*delt),0;0,0,1];h=Rotation*[Xs1;Ys1;Zs1] ;Xs=h(1);Ys=h(2);Zs=h(3);ss = sqrt((Xk0-Xs)^2+(Yk0-Ys)^2+(Zk0-Zs)*(Zk0-Zs));ts = tr(6)-ss/c;sign = norm(ts-tt(6));endaxk = (Xk0-Xs)/ss;ayk = (Yk0-Ys)/ss;azk = (Zk0-Zs)/ss;EAB=CAL2POL(Xk0,Yk0,Zk0,Xs,Ys,Zs);if j==1A = [axk,ayk,azk,1];L = R(j)-ss++c*dt-2.47/(sin(EAB)+0.0121);elseA = [A;axk,ayk,azk,1]; %构造误差方程L = [L;(R(j)-ss++c*dt)-2.47/(sin(EAB)+0.0121)];%常数项endendif Snum==4dx=inv(A)*L;aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);elseif Snum > 4dx = inv(A'*A)*A'*L; %构造法方程并求解aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);endendXp=sum(Px.*(Xk0+xchange))/sum(Px)Yp=sum(Py.*(Yk0+ychange))/sum(Py)Zp=sum(Pz.*(Zk0+zchange))/sum(Pz)figure(1);subplot(3,1,1);plot(xchange,'black--');subplot(3,1,2);plot(ychange,'black--');subplot(3,1,3);plot(zchange,'black--');tocfunction [x,y,z]= CalPos(PRN,time)global EphDatt1=TimetoJD(time(1),time(2),time(3),time(4),time(5),time(6));t2=TimetoJD(time(1),time(2),time(3),0,0,0);if isempty(EphDat)EphDat=ReadGpsData;endsz=size(EphDat,2);gg=0;%判断最近的时间for i=1:szif(EphDat(i).PRN==PRN & abs(EphDat(i).toc-t1)<0.0417*2)G=EphDat(i);gg=1;break;endendt3=t2-G.weekno*7;% 减整周数t=t3*24*3600+time(4)*3600+time(5)*60+time(6);%化为GPS秒u=3.986004418E+14; %地球引力常数we=7.2921151467E-5; %地球自转速度a=G.sqa^2; %地球长半轴n0=sqrt(u/a^3); %计算的平运动tk=t-G.toe; %从参考历元开始的时间if tk>=302400tk=tk-604800;elseif tk<-302400tk=tk + 604800;endn=n0+G.dn; %改正后的运动Mk=G.M0+n*tk;Ek=Mk; %偏近点角弧度for i=1:4Ek=Mk+G.e*sin(Ek);endfk=2*atan((sqrt((1+G.e)/(1-G.e))*tan(Ek/2))); %真近点角if(fk<0)fk=fk+2*pi;endcosf=(cos(Ek)-G.e)/(1-G.e*cos(Ek));sinf=(sqrt(1-G.e^2)*sin(Ek))/(1-G.e*cos(Ek));uk=fk+G.w; %升交角矩及轨道摄动改正项iuk=G.Cuc*cos(2*uk)+G.Cus*sin(2*uk);irk=G.Crc*cos(2*uk)+G.Crs*sin(2*uk);iik=G.Cic*cos(2*uk)+G.Cis*sin(2*uk);uk=uk+iuk ; %改正后的升角距r=a*(1-G.e*cos(Ek))+irk; %改正后的地心相径i=G.i0+iik+G.iDot*tk; %改正后的倾角xx=r*cos(uk); %在轨道面中的位置yy=r*sin(uk);omg=G.omg0+(G.omg-we)*tk-we*G.toe; %改正后的升交点经度x=xx*cos(omg)-yy*cos(i)*sin(omg);y=xx*sin(omg)+yy*cos(i)*cos(omg);z=yy*sin(i);CalPos=[x,y,z];%打印结果% x=num2str(x,'%12.6f');% y=num2str(y,'%12.6f');% z=num2str(z,'%12.6f');% fprintf('卫星号PRN:%2.0f',PRN); fprintf('\n');% fprintf('时间time:%4.0f%4.0f%4.0f%4.0f%4.0f%4.1f',time); fprintf('\n');% fprintf('所求卫星在地心坐标系中的空间坐标为：\n');% fprintf('X = ');fprintf(x,'\n');fprintf('m\n');% fprintf('Y = ');fprintf(y,'\n');fprintf('m\n');% fprintf('Z = ');fprintf(z,'\n');fprintf('m\n');%returnfunction EAB=CAL2POL(X,Y,Z,XB,YB,ZB)format longL=atan(Y/X);R=sqrt(X*X+Y*Y+Z*Z);a=6378137;e=sqrt(0.0066943799013);seta=atan(Z/sqrt(X*X+Y*Y));B=asin(sqrt((a*a-R*R*cos(seta)*cos(seta))/(a*a-R*R*cos(seta)*cos(seta)*e*e)))B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)*sin(B))));while abs(B1-B)>0.01*pi/180/3600B=B1;B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)^2)));endH=R*cos(seta)/cos(B)-(a/sqrt(1-e*e*sin(B)*sin(B)));L=L*180/pi;B=B*180/pi;DX=XB-X;DY=YB-Y;DZ=ZB-Z;NB=-sin(B)*cos(L)*DX-sin(B)*sin(L)*DY+cos(B)*DZ; EB=-sin(L)*DX+cos(L)*DY;UB=cos(B)*cos(L)*DX+cos(B)*sin(L)*DY+sin(B)*DZ; SAB=sqrt(NB*NB+UB*UB+EB*EB);EAB=asin(UB/SAB);。
第6章 单点定位
• 单位权中误差，其受伪距测量精度、星历精度及大气延迟 影响；
• 对应的协因数矩阵，它由卫星的空间几何分布决定
6.1 伪距单点定位
协因数矩阵中各个元素反映了在特定的卫星空间几何分布 下，不同参数的定位精度及其相关性信息。因此，利用这些信 息即可描述卫星空间几何分布对定位精度影响的精度因子: 常用的精度因子有： （1）几何精度因子（Geometric Dilution of Precision, GPOP）
式中 ， Qx
qYX
qYY
qYZ
R
sin L0
cos L0
0
为坐标转换矩阵，
qZX qZY qZZ
cos B0 cos L0 cos B0 sin L0 sin B0
B0和L0分别为对应的大地纬度和大地经度。
6.1 伪距单点定位
由此可得到另两个常用的精度因子 （4）水平精度因子（Horizontal Dilution of Precision, GPOP）
23934824.154
23978631.5766+43
822.577
24181945.803 24298916.7595+116 967.755
22957572.280 22965399.9529+7 834.145
22385541.968 22355209.7858 30 330.506
ni
1
VX VY
VZ
=i
i0
+cV t
i
S
Ii Ti i
cVtR
若接收机同时接收n（ n ≥4）颗卫星，则上式可写为：
6.1 伪距单点定位
l1
l2
l3
单点定位设备小博士使用流程
单点定位设备小博士使用流程
一、翻页键：
1按动此键将循环显示各个主页面2从某种操作中退出到主页面电源键：
1持续按住此键将开机或关机。
2短时间按下此键将打开或关闭背景光。
二、上下键：
1在各页面或菜单中，上下移动光标。
2在卫星状态页面中，调节屏幕显示对比度
3在航迹导航页面中，放大或缩小比例尺。
4在罗盘导航页面中，查看各种数据。
输入键：
1激活光标所在选项。
2确认菜单选项。
3在可以进行输入操作的地方输入数据。
ATTENTION：etrex（小博士）被设计为左手握机并操作，但用右手也可以很好地握机和操作，这取决于您的个人习惯。
将接收机拿到室外开阔的地点，显示屏向使其内置天线朝向开阔的天空。
按住电源您将看到欢迎画面，随后进入“卫星状态ATTENTION：当您第一次使用该设备时，置，以后将只需要15-45秒时间来定位。
当足够的卫星被锁定时，页面顶部的窗口航”，如果您的“小博士”无法接收到足够如您正在室内或者当前的GPS卫星信号很好。
GPS单点定位算法及实现
GPS单点定位算法及实现GPS单点定位算法是通过接收来自卫星的信号，通过计算接收信号到达时间差以及接收信号强度等信息，确定自身的位置坐标。
常见的GPS单点定位算法包括最小二乘法定位算法、加权最小二乘法定位算法、无拓扑算法等。
最小二乘法定位算法是一种基本的GPS定位算法，通过最小化测量误差的平方和，求得位置坐标最优解。
该算法假设接收器没有任何误差，并且卫星几何结构是已知的。
具体实现步骤如下：1.收集卫星信息：获取可见卫星的位置和信号强度信息。
2.数据预处理：对接收信号进行滤波和数据处理，例如去除离群点、噪声滤除等。
3.卫星定位计算：根据接收器和可见卫星之间的距离和相对几何关系，计算每颗卫星与接收器之间的距离。
4.平面定位计算：根据卫星位置和距离信息，使用最小二乘法求取接收器的经度和纬度。
5.高度定位计算：根据卫星位置和距离信息，使用最小二乘法或其他方法求取接收器的高度。
加权最小二乘法定位算法在最小二乘法定位算法的基础上加入对测量数据的加权处理，以提高定位精度。
加权最小二乘法定位算法的实现步骤与最小二乘法定位算法类似，只是在卫星定位计算和平面定位计算中，对每个测量值进行加权处理。
无拓扑算法是一种基于统计的定位算法，不需要事先知道接收器和卫星的几何关系，而是通过分析多个卫星的信息来确定接收器的位置。
其实现步骤如下：1.收集卫星信息：获取可见卫星的位置和信号强度信息。
2.数据预处理：对接收信号进行滤波和数据处理，例如去除离群点、噪声滤除等。
3.卫星选择：选择可见卫星中信号强度最强的几颗卫星。
4.定位计算：根据已选择的卫星信息，使用统计模型或其他算法计算接收器的位置。
1.数据采集与处理：获取和处理接收信号、卫星信息和测量数据，对数据进行有效的滤波和预处理。
2.算法选择与优化：根据定位精度和计算效率的要求，选择合适的算法，并进行算法优化和参数调整。
3.数据处理与结果可视化：对定位结果进行处理和分析，可通过地图等方式可视化结果，以便用户更直观地了解定位情况。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
GPS单点定位程序文档说明
一程序说明
本程序的基本功能是利用测站接收机获取的观测文件（O文件），结合相应历元下的卫星导航文件（N文件）来计算测站点在WGS-84坐标系下的坐标。
所用编程语言为C++，编程环境为VC++6.0。
二单历元卫星坐标计算流程图
三程序设计流程图
四程序模块说明
（一）对类的说明
1.PointPosition类
double GetDelta_t(CommonTime Toe,CommonTime Toc);
用于计算两个历元时刻的时间间隔
CoorCartesian ComputeSatCoor(double Tk,OneNavData oneNaveData);
利用导航数据计算卫星坐标
Bool GetPreciseSatCoor(CommonTimeObsTime,CoorCartesian&Sitecoor , CoorCartesian &Satcoor,OneNavData oneNavData,double SatClkBais);
通过迭代得到新号发送时刻的卫星坐标
double ComputeSatClkBias(double SendTime_Tk,OneNavData oneNavData);
计算卫星的钟差改正
Factor ComputeFactors(CoorCartesian &Sitecoor, CoorCartesian &Satcoor, double &satClkBias,double &Tr, double TropDelay, double ionDelay, double p1);
计算组建法方程需要的各个元素
double ComputeTropDelay(CoorCartesian SatCoor,CoorCartesian SitCoor);
计算对流层误差
double ComputeIonDelay(const double L1,const double L2);
计算电离层误差
2．ReadObsData类
ReadObsFile(const string& FileName);
ObsFileHeader ReadObsHeader(const string &FileName);
读取观测文件的头文件部分
EntirObsData ReadObsData(const string &FileName);
读取观测文件的数据部分
3．ReadNavData类
NavFileHeader ReadNavHeader(const string &FileName);
读取导航文件的头文件部分
AllNavData ReadNavData(const string &FileName);
读取导航文件的数据部分
OneNavData SelectEpochNavData(AllNavData allNavData,string SatPrn,CommonTime ObsTime)
查找并获取要计算的观测历元下的导航数据
（二）程序模块的连接关系
1.分别用文件流打开相应的观测数据文件和星历文件；
2.调用ReadObsFile.ReadObsHeader()读取观测文件头文件；
3.调用ReadObsFile.ReadObsData()按照观测历元读取观测数据；
4调用ReadNavFile.ReadNavHeader()读取导航文件的头文件部分；
5调用ReadNavFile.ReadNavData()读取导航文件的数据部分；
6 利用循环并调用ReadNavFile.SelectEpochNavData()获取观测历元的导航数据
7 调用puteSatCoor()计算卫星坐标
8利用puteFactor()计算法方程的各个元素
9迭代得到结果。
五程序的不足之处
1.程序较多的出现类调其他类中的方法，从而类的独立性较差。
2.在坐标系统中，只进行了空间直角坐标系和大地坐标系的转换，空间直角坐
标系和测站坐标系的转换，没有完全整理坐标转换。
3.计算的精度较低。
4.只利用GPS卫星的观测数据并没有利用Glonass卫星的观测数据，从而使数
据为得到充分利用。
最新文件---------------- 仅供参考--------------------已改成word文本--------------------- 方便更改。