GPS计算例题

例1 图为一简单GPS网，用两台GPS接收机观测，测得5条基线向量，n?15，每一个基线向量中三个坐标差观测值相关，由于两台GPS接收机接收机观测，所以各观测基线向量相互独立，网中点LOC1其三维坐标已知，其余三个为待定点，参数个数t?9。

图1

已知点信息（单位：m） LOC1

1

X -1974638.7340 Y 4590014.8190 Z 3953144.9235

2

编号 起点 终点 ?X -1218.561 ?Y -1039.227 ?Z 1737.720 基线方差阵 1 LOC2 LOC1 对?2.320999E-7??-5.097008E-71.339931E-6? 称?-4.371401E-71.109356E-61.008592E-6???对?1.044894E-6??-2.396533E-66.341291E-6? 称?-2.319683E-65.902876E-66.035577E-6???对?5.850064E-7??-1.329620E-63.362548E-6? 称?-1.252374E-63.069820E-63.019233E-6???2 LOC4 LOC1 270.457 -5.3.208 1879.923 3 LOC4 LOC2 1489.013 536.030 142.218 4 LOC3 LOC2 1405.531 -178.157 1171.380 对?1.205319E-6??-2.636702E-66.858585E-6? 称?-2.174106E-65.480745E-64.820125E-6???对?9.662657E-6??-2.175476E-55.194777E-5? 称?-1.971468E-54.633565E-54.324110E-5???5 LOC4 LOC3 83.497 714.153 -1029.199

3

编程代码： % 输入基线阵

[Filename pathName]=uigetfile({'*.txt','Txt Files|*.txt'},'输入基线差'); L=length(Filename); if L<5

msgbox('wrong File','File open Error'); end

fid = fopen([pathName Filename]); %打开数据文件 tline1 = fgetl(fid); str = '[-]?\\d[\\.]?\\d*';

f1=str2double(regexp(tline1,str,'match')); m_Pnumber=f1(1,1);%总点数 m_knPnumber=f1(1,2);%已知点数 m_Lnumber=f1(1,3);%总的观测值 tline2= fgetl(fid);

st = '[\\w][\\w][\\w][-]?\\d[\\.]?\\d*';

pname=regexp(tline2,st,'match');%提取点名

[P1,P2,P3,P4,P5]=textread([pathName Filename],'%s%s%s%s%s','headerlines',2); %输入基线方差阵

[Filename1 pathName1]=uigetfile({'*.txt','Txt Files|*.txt'},'输入基线方差阵'); L=length(Filename); if L<5

msgbox('wrong File','File open Error'); end

a=textread([pathName1 Filename1]);

for i=1:m_knPnumber %提取与已知点相关数据 pname1{i}=P1{i};

point2(i)=str2num(cell2mat(P2(i))); point3(i)=str2num(cell2mat(P3(i))); point4(i)=str2num(cell2mat(P4(i))); end

XYZ=zeros(1,3*m_Pnumber);%命名点坐标大小 for i=1:m_knPnumber

c=strmatch(pname1(i),pname);%将点名转换成点号 XYZ(1,3*c-2)=point2(i); XYZ(1,3*c-1)=point3(i); XYZ(1,3*c)=point4(i); end

%将未知点的x坐标赋值为10^30 for i=1:(m_Pnumber-m_knPnumber)

XYZ(1,3*(m_knPnumber+i)-2)=10^30; end

for i=1:m_Lnumber

dir1(i)=strmatch(P1(i+m_knPnumber),pname);%起点点号

4

dir2(i)=strmatch(P2(i+m_knPnumber),pname);%终点点号 XX(i)=str2num(cell2mat(P3(i+m_knPnumber)));%x坐标 YY(i)=str2num(cell2mat(P4(i+m_knPnumber)));%y坐标 ZZ(i)=str2num(cell2mat(P5(i+m_knPnumber)));%z坐标 end

%计算未知点近似坐标 for i=1:m_Lnumber

k1=dir1(i);%起点点号 k2=dir2(i);%终点点号

if XYZ(1,3*k1-2)<=10^29 & XYZ(1,3*k2-2)>=10^29 %起点已知，终点未知 XYZ(1,3*k2-2)=XYZ(1,3*k1-2)+XX(i); XYZ(1,3*k2-1)=XYZ(1,3*k1-1)+YY(i); XYZ(1,3*k2)=XYZ(1,3*k1)+ZZ(i);

elseif XYZ(1,3*k1-2)>=10^29 & XYZ(1,3*k2-2)<=10^29%起点未知，终点已知 XYZ(1,3*k1-2)=XYZ(1,3*k2-2)-XX(i); XYZ(1,3*k1-1)=XYZ(1,3*k2-1)-YY(i); XYZ(1,3*k1)=XYZ(1,3*k2)-ZZ(i); end end

%计算法方程系数

for i=1: m_Lnumber k1=dir1(i); k2=dir2(i);

B(3*i-2,3*k1-2)=-1;%X相关系数 B(3*i-2,3*k2-2)=1;

if k1<=m_knPnumber% k1点系数已知，该系数为0 B(3*i-2,3*k1-2)=0; end

if k2<=m_knPnumber% k2点系数已知，该系数为0 B(3*i-2,3*k2-2)=0; end

B(3*i-1,3*k1-1)=-1;%Y相关系数 B(3*i-1,3*k2-1)=1;

if k1<=m_knPnumber% k1点系数已知，该系数为0 B(3*i-1,3*k1-1)=0; end

if k2<=m_knPnumber% k2点系数已知，该系数为0 B(3*i-1,3*k2-1)=0; end

B(3*i,3*k1)=-1;%Z相关系数 B(3*i,3*k2)=1;

if k1<=m_knPnumber% k1点系数已知，该系数为0 B(3*i,3*k1)=0; end

5