基于matlab的测量导线描绘

1测绘学院学号:2009301610208夏少波手机:基于MATLAB的测量导线描绘摘要:本文从导线的描绘出发，提出了该问题的MATLAB解决方法。给出了数据的处理结果和图表。并且，设计了二维和三维的代码描述。其中对待处理的数据有一些优化，但计算方法和计算的逻辑是没有问题的。最后，MATLAB对这个问题的解决也较为令人满意。关键字:控制测量导线代码运行结果MATLAB任何测量工作都会产生误差，而误差在绝对意义上是不可避免的，所以在测量工作中必须采用一定的方法来控制误差的范围，以期达到工作需求。测量工作中遵循“先整体后局部，先控制后碎部”的原则。简而言之就是在工作区域内建立一定精度的由一定数目的已知坐标点构成的几何图形，然后在几何图形内进行测量。在测量过程中，我们可以通过各个已知点的坐标来验算测得的未知点的坐标。这样通过不断地检核、测量、检核，成果的误差就可以控制在一定范围内。其中，已知点被称作“控制点”，几何图形被称为“控制网”，这种测量方法被称为“控制测量”。一般意义上的控制测量可以分为“平面控制测量”和“高程控制测量”。平面控制测量是基平面坐标系X-Y，而高程控制测量是基于高程H的。目前，随着测绘技术的发展，也出现了把两种控制测量结合起来的“三维控制测量”。控制网根据其精度从高到低大致可以分为一、二、三、四等，精度越高的点数目越少。点数不够自然无法形成有效地控制网。因此，测量工作中需要通过已知的几个点，“一生二，二生三”，绘出一个可用于实际工程的控制网。“一生二，二生三”的过程需要已知点前后通视（即互相看得见），然而，在城市中，由于建筑物较多，视线不好，往往会出现已知点之间无法通视的现象。所以就引入了导线测量（traversesurvey）。导线测量又有很多分类，在此就不在赘述了。下面来介绍导线测量中的附和导线。所谓附和导线就是：“导线起始于一个已知控制点而终至于另一个已知控制点”。通过测量导线上未知点与已知控制点坐标的对比、计算，进行平差（控制误差的一种手段），得出一条导线上的各个待测点的坐标。这个过程就是导线测量的主体步骤了。导线示意图如下2根据导线示意图，这里有五点需要说明：一、测绘中，坐标系很多，一般导线测量使用的是“平面直角坐标系”，之所以打上了双引号，是因为，测量工作中的坐标系是以正北方向（竖轴）为X方向，横轴为Y方向，极坐标中顺时针旋转方向为正方向。二、图中2、3等点要安放仪器进行测量，称这些点为测站。三、A、B、C、D为已知点。P2、P3等为待测点。导线测量中需要计算出各个点的坐标。四、关于P1左、P2左等角度，为测绘中所谓的左角。所谓的左角可以通过导线的前进方向来确定。比如，在这张图中，测量的顺序为A-B-…-C-D,那么前进方向即为A-B-…-C-D,而前进方向的左手边的角度就是左角了（右角亦然，并不规定一定要测哪个角，计算方法略有不同而已）。这些角一般称作转折角。五、平面的导线控制测量需要在一测站测相邻的左右测站之间的夹角（如图中的P1），该测站到前进方向上的下一测站的距离（如图中的D(1-2)，也称作边长）。如果把高程算进去，还需要获得两测站之间的高差。而获得高差数据的方法很多，一般可以通过水准测量或三角高程测量获得。在这里介绍一下三角高程的基本原理。三角高程测量的基本原理就是利用三角形各边角关系来获得所需数据。如图二;已知A、B两点的水平距离D，和一个夹角α，通过计算D×tan(α)，再通过加减测量仪器高就可以获得A、B两点的高程差。当然实际运作中还要考虑地球曲率、大气温度等误差。在本文中，我对模型作了简化，对仪器高和地球曲率等因素做了默认设置。图二以上为本文要讨论的问题背景做了一个比较完整的描述下面进入问题的主体。一、提出问题有一测量小组通过查询获得了三个已知点的坐标，又通过野外测量，做了导线，并获得了导线中的转折角和边长，同时他们也获得了各个测站的α角。数据如下：3A列是边长，B列是左角值A列是α角值现在要求获得各个待测点的坐标，并描绘出导线的形状，对比已知点得出测量值与实际获得坐标值的偏差。二、对问题的分析和方案的产生现在获得的数据已经可以画出导线了。问题是坐标系是不同的，那么可以进行转换。MATLAB数学坐标系的Y轴是正北方向，那么我们规定：Y轴北方向，X轴为东方向。化为极坐标时，顺时针为角的方向。确定坐标系后，我们可以根据已知坐标把做边长和转折角转化为坐标值。计算方法：α1=α2+α（转折角）－180°;X1=X2＋D(边长)×sin(α1);Y1=Y2＋D(边长)×cos(α1);注：此处的坐标转化考虑到坐标系的转化，故sin、cos互换了。通过以上计算获得了这些数据还是不够的。还需要用户输入已知点的坐标来对比测量结果。同时为了增加代码的通用性，可以由用户输入所需处理的文件名称，以及初始坐标值。同时，还可以通过用户选择，来实现计算和描绘平面的导线或者是三维的导线。通过用户选择那个函数，调用了二维描述和三维描述这两个函数。从而实现整个功能。如上所示，整个框架已经构成了。三、描述解决方案首先给出可供用户选择的函数：functionbegin()%这个函数给了用户选择。2Dor3D?input('welcometousemyprogram£¡pressEntertostar£¡')N=input('choose1for2D,choose2for3D')switch(N)case1daoxian();case2gaocheng();end4再给出描画二维导线的函数描述：functiondaoxian()%thisfunctionisfor2DFILENAME=input('whatisthenameofdatefile?')%获取用户所要处理的文件名X0=input('inputthex：')Y0=input('inputthey：')%用户输入已知坐标。Angle=atan(Y0/X0);%用户输入已知的坐标点后，进行角度的转化，使之能在坐标系中正常实现。A=xlsread(FILENAME);%获取用户所要处理的文件内容L=A(:,1);L=L';CT=A(:,2);CT=CT';%以上是对获取数据的处理，使之符合矩阵运算。GT=[];CT(1)=Angle;CT=(CT./(360))*2*pi;%把一般测得的360制角转化为2Pi制X=[0,X0,0,0,0,0,0];Y=[0,Y0,0,0,0,0,0];%初始化X,YGT(3)=CT(1)+CT(2)-pi;I=4;fork=1:4%转化各个转角，使结果符合坐标系的运算GT(I)=CT(I-1)+GT(I-1)-pi;I=I+1;endforI=1:7ifGT(I)0%对小于0度的角进行转换GT(I)=GT(I)+2*pi;endifGT(I)2*pi%对大于360度的角进行转换GT(I)=GT(I)-2*pi;endendGT(1)=0;GT(2)=0;I=3;forJ=1:5%核心计算，获得坐标。X(I)=X(I-1)+L(I-1).*sin(GT(I));Y(I)=Y(I-1)+L(I-1).*cos(GT(I));I=I+1;endOUTPUT=[X',Y']%输出数据saveOUTPUT%保存数据。plot(X,Y,'b-O');%画图。axis('equal');5xlabel('X');ylabel('Y,北方向');title('2D-----------map');%以上为坐标的优化gridon;holdon;%对比已知点，求出坐标差，并描点。CONx=input('pleaseenterthexofcontrolpoint：');CONy=input('pleaseentertheyofcontrolpoint：');plot(CONx,CONy,'ro');%描出用于检核的点。delx=X(7)-CONx%输出相差X,Y坐标大小。dely=Y(7)-CONy函数的描述和解释代码的后面已经阐述，现在直接给出三维的函数描述：functiongaocheng()%Thisfunctionisfor3DFILENAME1=input('whatisthenameoflengt-hdatefile?')%获得用户所需要处理的文件名。FILENAME2=input('whatisthenameofhigh-datefile?')X0=input('inputthex：')%获得初始的起算坐标Y0=input('inputthey：')Z0=input('inputthez：')Angle=atan(Y0/X0);%计算出起算的两点间在坐标系中的夹角。A=xlsread(FILENAME1);%获得用户提供的数据L=A(:,1);L=L';CT=A(:,2);CT=CT';%以上是对获取数据的处理，使之符合矩阵运算。GT=[];CT(1)=Angle;%用户输入已知的坐标点后，进行角度的转化，使之能在坐标系中正常实现。CT=(CT./(360))*2*pi;%转化角度X=[0,X0,0,0,0,0,0];%初始化坐标值Y=[0,Y0,0,0,0,0,0];CT(1)=Angle;GT(3)=CT(1)+CT(2)-pi;I=4;fork=1:4%计算获得可计算的坐标方位角GT(I)=CT(I-1)+GT(I-1)-pi;I=I+1;endforI=1:7ifGT(I)0%对小于0度的角进行转换GT(I)=GT(I)+2*pi;endifGT(I)2*pi%对大于360度的角进行转换GT(I)=GT(I)-2*pi;6endendGT(1)=0;GT(2)=0;I=3;forJ=1:5%核心计算，获得坐标值。X(I)=X(I-1)+L(I-1).*sin(GT(I));Y(I)=Y(I-1)+L(I-1).*cos(GT(I));I=I+1;endHT=xlsread(FILENAME2);%获得用户提供的高程α角。HT=HT';HT=(HT./(360))*2*pi;%α角转化；LT=[];fork=1:5LT(k)=L(k+1);endT=LT.*tan(HT);%通过tan函数获得H的初始值。Z=[];Z(1)=0;Z(2)=Z0;fork=3:7%获得H的值。Z(k)=Z(k-1)+T(k-2);endOUTPUT=[X',Y',Z']%输出计算结果saveOUTPUT;plot3(X,Y,Z,'k-O');%画出图像。axis('equal');xlabel('X');ylabel('Y,北方向');zlabel('H,高程');gridon;title('3D--------map');%以上为坐标的优化holdon;%获得已知点，进行比较，描绘，输出结果CONx=input('pleaseenterthexofcontrolpoint：');CONy=input('pleaseentertheyofcontrolpoint：');CONz=input('pleaseenterthezofcontrolpoint：');plot3(CONx,CONy,CONz,'ro')%描出用于检核的点。delx=X(7)-CONx%输出相差的坐标值dely=Y(7)-CONydelyz=Z(7)-CONz现在整个函数的实现已经给出来了。其中，functionbegin()放在一个M文件里，functiondaoxian()放在一个M文件里，functiongaocheng()放在一个M文件里。7其实，第2D和3D就差了一个坐标元素。主体函数还是差不多的。不过，为了表达清楚，我还是把它们都写出来了。我在写这些函数的时候，遇到了一些矩阵数组计算的小问题，最后虽然得到了解决，但我一直怀疑是不是有更好的描述方法。四、运行MATLAB,获得结果首先，我要说明的是，本来我是想在word里直接运行代码的，但是，由于不同的M文件，这样的表达，我觉得不是很清晰（可能我方法不对），所以就选择了截图运行过程的方法。运行时首