基于MATLAB的GPS网平差

测绘数据繁杂，数据处理很麻烦，运用简单编程设计程序实现普通测绘数据的平差，成图，自动结算，为测绘行业提供简便。本代码只针对简易GPS网进行平差，主要是秩亏gps网，无约束gps网，含有一个，两个或者无已知点的情况，其他情况不做讨论。

clc;
clear;
disp（‘*********************************************测绘工程16-2杨帆*****************************************************‘）
%读取txt文本
[filename1filepath]=uigetfile（‘*.txt‘‘请选择GPS数据文件‘）;%uigetfile显示检索文件的对话框，filename1为文件名
yf=fopen（strcat（filepathfilename1）‘rt‘）;%strcat:横向连接字符串
if （yf==-1）
   msgbox（‘出错！出错！出错！‘）; %magbox为弹出警告的指令
   return;
end

yzds=fscanf（yf‘%f‘1）;%已知点数
if yzds==0
    x01=fscanf（yf‘%f‘1）;%已知点x坐标
    y01=fscanf（yf‘%f‘1）;%已知点y坐标
    z01=fscanf（yf‘%f‘1）;%已知点z坐标
    wzds=fscanf（yf‘%f‘1）;%未知点数
    jxs=fscanf（yf‘%f‘1）;%基线数
    sigma0=fscanf（yf‘%f‘1）;%中误差
    i=1;
  while （i    qd（i）=fscanf（yf‘%f‘1）;%起点
    zd（i）=fscanf（yf‘%f‘1）;%终点
    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无约束网平差结果‘）
  %权阵P的形成
  disp（‘**********************************************************P阵**********************************************************‘）
  T=zeros（jxs*3jxs*3）;%建立行数列数为基线向量条数3倍的零矩阵，用来存储方差协方差阵
  for i=1:jxs %存协方差
        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）;
    %系数矩阵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]）;%对角阵
  end
  for i=1:jxs
        B（3*i-2:3*i3*zd（i）-2:3*zd（i））=diag（[111]）;
  end    
    disp（B）;
  %常数项矩阵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（‘-------------------------------------------单位权中误差（单位：m）--------------------------------------------‘）
    disp（strcat（‘σ0=‘num2str（delta0））） %num2str为把数值转换成字符串，转换后可以使用disp函数输出
    Q0=inv（NBB）;%σBbb-1
    Q=Q0*B‘*P*B*Q0;
    %由假定已知点计算其他未知点平差后坐标
    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