Matlab与测量数据处理实验报告

isCon=false; return; else isCon=true; return; end end end %利用要素名查找对应序号 function [index]=nameToindex(name,table) [row,col]=size(table); for i=1 : row if table(i,1)==name index=i; return; end end index=0;
lineIndex=i; return; end end end lineIndex=0; %由起算点开始计算附和水准路线程序 function [ok,ptTab,lineTab]=caculateHO(ptTab,lineTab,lineIndex) lbname=lineTab(lineIndex,2); lename=lineTab(lineIndex,3); lhh=lineTab(lineIndex,4); lbindex=nameToindex(lbname,ptTab); leindex=nameToindex(lename,ptTab); lbH=ptTab(lbindex,2); leH=ptTab(leindex,2); if lbH==0 && leH~=0 ptTab(lbindex,2)=ptTab(leindex,2)-lhh; nextpt=lbindex; finishpt=leindex; else
%以上为初始化变量 [row,col]=size(lineTab); for i=1 : row
if lineTab(i,5)==0 jsq=jsq+1;ineTab(i,8)=jsq;
else lineTab(i,8)=0;
end end [row,col]=size(angleTab); for i=1 : row
%水准网间接平差 clear; load('ptTab1','ptTab1'); ptTab=ptTab1; lineTab=lineTab1; [isok,ptTab,lineTab]=getHO(ptTab,lineTab); [P,B,l,ptTab,lineTab]=getV(ptTab,lineTab); X=-inv（B'*P*B)*B'*P*l; V=B*X+l; [row,col]=size(B); mo=sqrt(V'*P*V)/(row-col); M=mo^2*inv(B'*P*B); msgbox "计算正确！" %查找起算边函数 function [lineIndex]=getBegin(ptTab,lineTab) [rowNo,colNo]=size(lineTab); for i=1 : rowNo
Al(i,1)lf+wangle+pi()-k; end end if rbname==lename if lwz~=0 && rwzห้องสมุดไป่ตู้=0
AB(i,lwz)=-1;AB(i,rwz)1; Al(i,1)=lfw+angle-pi()-k; end if lwz==0 && rwz~=0 AB(i,rwz)=1; Al(i,1)=lfw+angle-pi()-k; end if lwz~=0 && rwz==0 AB(i,lwz)=-1; Al(i,1)=angle-rfw-pi()-k; end end end AX=inv(AB'*AB)*AB'*Al; AV=AB*AX-Al; [row_1,col_1]=size(AB); m0=sqrt((AV'*AV)/(row_1-col_1)); QQ=inv(AB'*AB); M=m0^2*QQ;AB=0; Al=0; [row1,col1]=size(AV);[row3,col3]=size(AX);[row4,col3]=size(lineTab); maxdet=100;maxindex=0; for mm=1 : row1 getvalue=AV(mm); if getvalue<maxdet && angleTab(mm,5)<2*pi() maxindex=mm;maxdet=getvalue; end end if maxindex~=0 && maxdet<pi()/1800 && angleTab(maxindex,5)<=3*pi()/2 angleTab(maxindex,5)=pi()+angleTab(maxindex,5); end end for mm=1 : row3 getvalue=AX(mm); for nn=1 :row4 if mm==lineTab(nn,8) lineTab(nn,5)=hdtodfm(getvalue); while lineTab(nn,5)>=360 lineTab(nn,5)=lineTab(nn,5)-360; end
leH=ptTab(lbindex,2); if ptTab(lbindex,4)~=0 && ptTab(leindex,4)~=0
B(i,ptTab(lbindex,4))=-1; B(i,ptTab(leindex,4))=1; l(i,1)=(-1)*lhh; end if ptTab(lbindex,4)==0 && ptTab(leindex,4)~=0 B(i,ptTab(leindex,4))=1; l(i,1)=(-1)*(lhh+lbH); end if ptTab(lbindex,4)==0 && ptTab(leindex,4)~=0 B(i,ptTab(lbindex,4))=-1; l(i,1)=lbH-lhh; end if ptTab(lbindex,4)==0 && ptTab(leindex,4)==0 msgbox "建立法方程出错！" return; end P(i,i)=10/lineTab(i,5); end %判断当前点要素是否是已知点 function [isCon]=isContr(ptIndex,ptTab) [rowNo,colNo]size(ptTab); for i=1 : rowNo if ptIndex==ptTab(i,1) if ptTab(i,2)==0
bx=ptTab(pbindex,2);by=ptTab(pbindex,3); ex=ptTab(peindex,2);ey=ptTab(peindex,3); detx=ex-bx;dety=ey-by; if lineTab(i,5)==0
tem=getfw(detx,dety); lineTab(i,5)=hdtodfm(tem); end else %判断数据是否完整 end end %调用子程序，开始计算各边要素近似方位角 [ptTab,angleTab,lineTab,QQ,ok]=getTrueFW03(ptTab,angleTab,lineTab); if not(ok) msgbox "计算近似方位失败！" return; end msgbox "计算成功！" %计算近似方位角子函数 function [ptTab,angleTab,lineTab,QQ,ok]=getTrueFW03(ptTab,angleTab,lineTab) JSQ=0;ab=0;lbzh=0;det=10; maxdet=10;maxindex=0;
beginIndex=lineTab(i,2); endIndex=lineTab(i,3); if isContr(beginIndex,ptTab) && not (isContr(endIndex,ptTab))
lineIndex=i; return; else if not (isContr(beginIndex,ptTab)) && isContr(endIndex,ptTab)
if ptTab(i,4)==0 ptTab(i,5)==jsq;jsq=jsq+1;
else ptTab(i,5)=0;
end end %查找计算起点终点都是已知点的边的方位角，统一方位角已知形式 [row,col]=size(lineTab); for i=1 : row
bname=lineTab(i,2);ename=lineTab(i,3); pbindex=nameToindex(bname,ptTab); peindex=nameToindex(ename,ptTab); if ptTab(pbindex,4)==1 && ptTab(peindex,4)==1
if lbH~=0 && leH==0 ptTab(leindex,2)=ptTab(lbindex,2)+lhh; nextpt=leindex; finishpt=lbindex;
else ok=true; return;
end end [row,col]=size(lineTab); nextbegin=0; for i=1 : row
if lineTab(i,2)==nextpt && finishpt~=lineTab(i,3) nextbegin=i; break;
else if lineTab(i,3)==nextpt && finishpt~=lineTab(i,2) nextbegin=i; break; end
end end ok=true; if nextbegin~=0
边号 1 2 3 4 5 6
P1 P2
点要素表
点号 A P1 P2 B
高程/m 15.016 —— —— 16.016
已知点标志 1 0 0 1
起点点号 A A B B P1 P1
边要素表
终点点号 P1 P2 P1 P2 P2 B
高差/m 1.359 2.009 0.363 1.012 0.657 -0.357
angleTab(i,5)=0; end while abs(maxdet)>pi()/1800
for i=1 :row angle=dfmtohd(angleTab(i,4)); k=angleTab(i,5); llname=angleTab(i,2);rlname=angleTab(i,3); llindex=nameToindex(llname,lineTab); rlindex=nameToindex(rlname,lineTab); lbname=lineTab(llindex,2);lename=lineTab(llindex,3); rbname=lineTab(rlindex,2);rename=lineTab(rlindex,3); lfw=dfmtohd(lineTab(llindex,5)); rfw=dfmtohd(lineTab(rlindex,5)); lwz=lineTab(llindex,8);rwz=lineTab(rlindex,8); if rbname==lbname||rename==lename if lwz~=0 && rwz~=0 AB(i,lwz)=-1;AB(i,rwz)=1; Al(i,1)=angle-k; end if lwz==0 && rwz~=0% AB(i,rwz)=1;Al(i,1)=angle+lfw-k; end if lwz~=0 && rwz==0 AB(i,lwz)=-1;Al(i,1)=angle-k-rfw; end end if rename==lbname if lwz~=0 && rwz~=0 AB(i,lwz)=-1;AB(i,rwz)=1; Al(i,1)=angle+pi()-k; end if lwz==0 && rwz~=0 AB(i,rwz)=1;
[ok,ptTab,lineTab]=caculateHO(ptTab,lineTab,nextbegin); end %计算误差方程系数矩阵函数 function [P,Bl,ptTab,lineTab]=getV(ptTab,lineTab) [row,col]size(ptTab); jsq=1; for i=1 : row
if ptTab(i,3)==0 ptTab(i,4)=jsq; jsq=jsq+1;
else ptTab(i,4)=0;
end end [row,col]=size(lineTab); for i=1 : row
lbname=lineTab(i,2);lename=lineTab(i,3);lhh=lineTab(i,4); lbindex=nameToindex(lbname,ptTab); leindex=nameToindex(lename,ptTab); lbH=ptTab(lbindex,2);
线路长/m 1.1 1.7 2.3 2.7 2.4 4.0
高程平差值 6.3748 7.0279
计算结果
单位权中误差 2.2248
高差平差值中误差 2.2080
高程平差值中误差 1.6208 1.9429
%导线网间接平差 %clear; format long; load('angleTab','angleTab'); load('lineTab','lineTab'); load('ptTab','ptTab'); [row,col]=size(ptTab); jsq=1; for i=1 : row