clc;
clear;
disp（‘*********************************************測(cè)繪工程16-2楊帆*****************************************************‘）
%讀取txt文本
[filename1filepath]=uigetfile（‘*.txt‘‘請(qǐng)選擇GPS數(shù)據(jù)文件‘）;%uigetfile顯示檢索文件的對(duì)話框，filename1為文件名
yf=fopen（strcat（filepathfilename1）‘rt‘）;%strcat:橫向連接字符串
if?（yf==-1）
???msgbox（‘出錯(cuò)！出錯(cuò)！出錯(cuò)！‘）;?%magbox為彈出警告的指令
???return;
end

yzds=fscanf（yf‘%f‘1）;%已知點(diǎn)數(shù)
if?yzds==0
????x01=fscanf（yf‘%f‘1）;%已知點(diǎn)x坐標(biāo)
????y01=fscanf（yf‘%f‘1）;%已知點(diǎn)y坐標(biāo)
????z01=fscanf（yf‘%f‘1）;%已知點(diǎn)z坐標(biāo)
????wzds=fscanf（yf‘%f‘1）;%未知點(diǎn)數(shù)
????jxs=fscanf（yf‘%f‘1）;%基線數(shù)
????sigma0=fscanf（yf‘%f‘1）;%中誤差
????i=1;
??while?（i????qd（i）=fscanf（yf‘%f‘1）;%起點(diǎn)
????zd（i）=fscanf（yf‘%f‘1）;%終點(diǎn)
????dx（i）=fscanf（yf‘%f‘1）;%△x
????dy（i）=fscanf（yf‘%f‘1）;%△y
????dz（i）=fscanf（yf‘%f‘1）;%△z
????sigmaxx（i）=fscanf（yf‘%f‘1）;%σxx
????sigmayy（i）=fscanf（yf‘%f‘1）;%σyy
????sigmazz（i）=fscanf（yf‘%f‘1）;%σzz
????sigmayx（i）=fscanf（yf‘%f‘1）;%σyx
????sigmazy（i）=fscanf（yf‘%f‘1）;%σzy
????sigmazx（i）=fscanf（yf‘%f‘1）;%σzx
????i=i+1;
??end
disp（‘**************************************************************************************************************‘）
disp（‘????????????????????????????????????????????GPS無(wú)約束網(wǎng)平差結(jié)果‘）
??%權(quán)陣P的形成
??disp（‘**********************************************************P陣**********************************************************‘）
??T=zeros（jxs*3jxs*3）;%建立行數(shù)列數(shù)為基線向量條數(shù)3倍的零矩陣，用來(lái)存儲(chǔ)方差協(xié)方差陣
??for?i=1:jxs?%存協(xié)方差
????????T（3*i-13*i-2）=sigmayx（i）;%σyx
????????T（3*i3*i-1）=sigmazy（i）;%σzy
????????T（3*i3*i-2）=sigmazx（i）;%σzx
??end
????T;
????D=T+T‘;%方差陣
????disp（D）;
??for?i=1:jxs?%存方差
????????D（3*i-23*i-2）=sigmaxx（i）;%σxx
????????D（3*i-13*i-1）=sigmayy（i）;%σyy
????????D（3*i3*i）=sigmazz（i）;%σzz
??end
????D;
????P=inv（（D/sigma0^2））;%p=inv（Q）q=D/σ0^2
????disp（P）;
????%系數(shù)矩陣B的形成
????disp（‘**********************************************************B陣**********************************************************‘）
????B=zeros（jxs*3wzds*3）;%nxt
??for?i=1:jxs
????????B（3*i-2:3*i3*qd（i）-2:3*qd（i））=diag（[-1-1-1]）;%對(duì)角陣
??end
??for?i=1:jxs
????????B（3*i-2:3*i3*zd（i）-2:3*zd（i））=diag（[111]）;
??end????
????disp（B）;
??%常數(shù)項(xiàng)矩陣l的形成
??disp（‘**********************************************************l陣**********************************************************‘）
??l=zeros（jxs*31）;
??for?i=1:jxs
??????l（3*i-2:3*i）=[dx（i）;dy（i）;dz（i）];
??end
???disp（l）;
?%間接平差
????NBB=B‘*P*B+G0*G;
????w=B‘*P*l;
????xp=inv（NBB）*w;
????v=B*xp-l;%誤差方程
????delta0=sqrt（（v‘*P*v）/（wzds-1））;%σ0
????disp（‘-------------------------------------------單位權(quán)中誤差（單位：m）--------------------------------------------‘）
????disp（strcat（‘σ0=‘num2str（delta0）））?%num2str為把數(shù)值轉(zhuǎn)換成字符串，轉(zhuǎn)換后可以使用disp函數(shù)輸出
????Q0=inv（NBB）;%σBbb-1
????Q=Q0*B‘*P*B*Q0;
????%由假定已知點(diǎn)計(jì)算其他未知點(diǎn)平差后坐標(biāo)
????jxxl=v+l;
??for?i=1:wzds-1
????????X（1）=x01;Y（1）=y01;Z（1）=z01;
????????X（i+1）=X（i）-jxxl（3*i-2）;
????????Y（i+1）=Y（i）-jxxl（3*i-1）;
????????Z（i+1）=Z（i）-jxxl（3*i）;
??end
????dh=ze