水准网平差程序Matlab实现 全部代码,详细教程
生活随笔
收集整理的這篇文章主要介紹了
水准网平差程序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 dataline = fgetl(fid);observs(:, iii) = str2num(line); end2 初始化高程向量,若其為未知點,則高程設為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_totalif(sum(kkk == idx_known) == 1)h0(kkk) = h_known(ttt);ttt = ttt+1;elseh0(kkk) = 1e6;end end3 利用while函數和for函數求其初始高程:
idx_unknown = find(h0 == 1e6); %calculate the origin height of each points while(sum(h0 == 1e6)>=1)for ttt = 1:num_observstarts= 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);endend end4 求解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 -1starts_index = find(idx_unknown == starts_) ; B(hhh, starts_index) = -1;endif(sum(find( idx_known == ends_ )) == 0)%if the end point is unknown set it to 1ends_index = find(idx_unknown == ends_) ;B(hhh, ends_index) = 1;endl(hhh) = -h0(starts_) + h0(ends_) - h(hhh); end5 利用得到的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_totalif(sum(find( idx_known == sss )) == 0 )sigma(sss) = sigma0*sqrt(Qxx(ttt,ttt));dh(sss) = x_(ttt);ttt = ttt+1;elsesigma(sss) = 0;dh(sss) = 0;end end6 輸出
%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_totalfid3 = 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總結
以上是生活随笔為你收集整理的水准网平差程序Matlab实现 全部代码,详细教程的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Origin-画局部放大图,并和原图放在
- 下一篇: redis win10开机自启