levelnet.m

水准网平差MATLAB程序

function [h,uwl,uwh]=levelnet(pathname,net_name);
clear
global ed dd sd pn gd h0 k1 k2 h1 s dh;% k11 k12
global pathname net_name; %s_datafile a1_datafile b_datafile
[ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata();
[dh,h,V,L,uw0,uwh,uwl]=calculatelevelnet(ed,dd,sd,gd,pn,h0,k1,k2,h1,s);
% writelevenetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0);
return
% 平差计算
function [dh,h,V,L,uw0,uwh,uwl]=calculatelevelnet(ed,dd,sd,gd,pn,h0,be,en,hd,distance);
A=sparse(zeros(sd,gd));
b=(0:(gd-1))*sd;
A(be'+b)=-1;
A(en'+b)=1;
A=A';
A=A(:,1:dd);                   %系数矩阵
l=zeros(gd,1);
l=h0(be)-h0(en)+hd;            %L矩阵
dia=distance;
p=inv(diag(dia));              %权矩阵
dh=inv(A'*p*A)*A'*p*l;         %X矩阵
h00=h0(dd+1:sd);               %已知点向量
h0=h0(1:dd);                   %待定点高程近似值
h=h0+dh;                       %待定点高程平差值
V=A*dh-l                       %高差观测值改正数
L=hd+V                         %观测高差平差值
uw0=sqrt(V'*p*V/(gd-dd));      %单位权中误差
Qxx=inv(A'*p*A);               %协因数矩阵
uwh=uw0*sqrt(diag(Qxx));       %待定点高程平差值中误差
uwh(dd+1:ed+dd)=0.0;
Qff=A*Qxx*A';                  
uwl=uw0*sqrt(diag(Qff));       %高差平差值中误差
h=[h;h00]                      %所有点高程
h0=[h0;h00]
ccc=dot(h,h0);
dh=[dh;zeros(ed,1)]
return

function [ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata();
% global ed dd sd pn gd h0 k1 k2 h1 s dh k11 k12;
% global pathname net_name s_datafile a1_datafile b_datafile;
k1=[];k2=[];h=[];s=[];
% if(isempty(pathname)|isempty(net_name))

    [filename,pathname]=uigetfile('*.txt','input filename');
%     i=find('.'==filename);
    net_name=filename;
% end
% pname='E:\program\C#\DataProcess\matlabDLL\';
% netname='aaa.txt';

fid1=fopen(strcat(pathname,net_name),'rt');
% fid1=fopen(strcat(pname,netname),'rt');
if(fid1==-1)
    msgbox('input file or path is not correct','warning','warn');
    return
end
ed=fscanf(fid1,'%f',1);
dd=fscanf(fid1,'%f',1);
sd=ed+dd;
gd=fscanf(fid1,'%f',1);
pn=fscanf(fid1,'%f',sd);
h0=fscanf(fid1,'%f',ed);
h0(dd+1:ed+dd)=h0(1:ed);
heightdiff=fscanf(fid1,'%f',[4,gd]);
heightdiff=heightdiff';
k1=heightdiff(:,1);
k2=heightdiff(:,2);
k11=heightdiff(:,1);
k12=heightdiff(:,2);
h1=heightdiff(:,3);
s=heightdiff(:,4);
fclose('all');
% 近似高程计算
[k1,k01]=chkdat(sd,pn,k1);
[k2,k02]=chkdat(sd,pn,k2);
h0(1:dd)=20000;
ie=0;
while(1)
    for k=1:gd
        i=k1(k);
        j=k2(k);
        if(h0(i)<1e4&h0(j)>1e4)
            h0(j)=h0(i)+h1(k);
            ie=ie+1;
        end
        if(h0(i)>1e4&h0(j)<1e4)
            h0(j)=h0(i)-h1(k);
            ie=ie+1;
        end
    end
    if(ie==dd)
        break;
    end
end
h0=reshape(h0,length(h0),1);
return

function [h1,k]=chkdat(sd,pn,h1);
n=length(h1);
k=0;
for i=1:n
    i1=0;
    for j=1:sd
        if(h1(i)==pn(j))
            i1=1;
            h1(i)=j;
            break;
        end
    end
    if(i1==0)
%         fprinf(fit2,'%5d %5d\n',i,n1(i));
        k=1;
    end
end
return


% function writelevenetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0)
% fid2