水准网平差程序Matlab实现 全部代码，详细教程

一、程序结构：

输入文件格式如下所示：

 

1 读取文件：主要用到fopen, str2num函数

fid = fopen('input_leveling.txt','r');
line1 =  fgetl(fid); 
vec1 = str2num(line1);
%read the first line and change it to number vector (the number of known points, unknown points, and obervations)
num_known = vec1(1);
num_unknown = vec1(2);
num_observ = vec1(3);
line2 = fgetl(fid); %read the second line(the index of known point)
idx_known = str2num(line2);
line3 = fgetl(fid);% read the third line(the height of known point)
h_known = str2num(line3);
line4 = fgetl(fid)
flag = str2num(line4);
for iii = 1:num_observ %  read the obersvation data
        line = fgetl(fid);
        observs(:, iii) = str2num(line);
end

2 初始化高程向量，若其为未知点，则高程设为1e8，若为已知点，则为其高程。

start_pt = observs(1,:); % start point
end_pt = observs(2,:);% end point
h = observs(3,:);% oberved height
dist = observs(4,:);% distances of each baseline
num_total = num_known+num_unknown;
ttt = 1;
%create the origin height vector 
%for unknown points -> 1e6 
%for known points -> height
for kkk = 1:num_total
        if(sum(kkk == idx_known) == 1)
                h0(kkk) = h_known(ttt);
                ttt = ttt+1;
        else
                h0(kkk) = 1e6;
        end
end

3 利用while函数和for函数求其初始高程：

idx_unknown =  find(h0 == 1e6);
%calculate the origin height of each points
while(sum(h0 == 1e6)>=1)
        for ttt = 1:num_observ
                starts= start_pt(ttt);
                ends = end_pt(ttt);
                if((h0(starts)>= 1e6)&&(h0(ends)<1e6))
                        h0(starts) = h0(ends) - h(ttt);
                elseif((h0(starts)< 1e6)&&(h0(ends)>=1e6))
                        h0(ends) = h0(starts) + h(ttt);
                end
        end
end

4 求解B矩阵, l 矩阵以及权重矩阵：

p = diag(1./dist);% p is the weight matrix 
B = zeros(num_observ, num_unknown);% B 
for hhh = 1:num_observ;
        starts_ = start_pt(hhh);
        ends_ = end_pt(hhh);
        if(sum(find( idx_known == starts_ )) == 0)%if the start point is unknown set it to -1
                starts_index = find(idx_unknown == starts_) ;      
                B(hhh, starts_index) = -1;
        end
        if(sum(find( idx_known == ends_ )) == 0)%if the end point is unknown set it to 1
                ends_index = find(idx_unknown == ends_) ;
                B(hhh, ends_index) = 1;
        end
        l(hhh) = -h0(starts_) + h0(ends_) - h(hhh);
end

5 利用得到的B, p, l进行平差：

%adjustment
l = -l;

Nbb = B'*p*B;
w = B'*p*l';
x_ = inv(Nbb)*w;
v = B*x_ - l';
num_r = num_observ-(num_unknown);
sigma0 = sqrt(v'*p*v/num_r);
Qxx = inv(B'*p*B);
ttt = 1;
for sss = 1: num_total
        if(sum(find( idx_known == sss )) == 0 )
        sigma(sss) = sigma0*sqrt(Qxx(ttt,ttt));
        dh(sss) = x_(ttt);
        ttt  = ttt+1;
        
        else
                sigma(sss) = 0;
                dh(sss) = 0;
        end
end

6 输出

%outputs the results
fid2 =  fopen('adjustment_results.txt','wt')
fprintf(fid2, '%4s   %4s        %4s       %4s         %4s\n', 'idx','h0','dh','h','theta_h')
for ppp = 1: num_total
       fid3 =  fopen('adjustment_results.txt','a') 
        fprintf(fid3, '%2d   %4.4f   %4.4f   %4.4f   %4.4f\n', ppp, h0(ppp),dh(ppp),h0(ppp)+dh(ppp),sigma(ppp))
end