基于matlab的测量导线描绘

测绘学院学号:2009301610208 夏少波手机:

基于MATLAB的测量导线描绘

摘要: 本文从导线的描绘出发，提出了该问题的MATLAB解决方法。给出了数据的处理结果和图表。并且，设计了二维和三维的代码描述。其中对待处理的数据有一些优化，但计算方法和计算的逻辑是没有问题的。最后，MATLAB对这个问题的解决也较为令人满意。

关键字: 控制测量导线代码运行结果MATLAB

任何测量工作都会产生误差，而误差在绝对意义上是不可避免的，所以在测量工作中必须采用一定的方法来控制误差的范围，以期达到工作需求。测量工作中遵循“先整体后局部，先控制后碎部”的原则。简而言之就是在工作区域内建立一定精度的由一定数目的已知坐标点构成的几何图形，然后在几何图形内进行测量。在测量过程中，我们可以通过各个已知点的坐标来验算测得的未知点的坐标。这样通过不断地检核、测量、检核，成果的误差就可以控制在一定范围内。其中，已知点被称作“控制点”，几何图形被称为“控制网”，这种测量方法被称为“控制测量”。

一般意义上的控制测量可以分为“平面控制测量”和“高程控制测量”。平面控制测量是基平面坐标系X-Y，而高程控制测量是基于高程H的。目前，随着测绘技术的发展，也出现了把两种控制测量结合起来的“三维控制测量”。

控制网根据其精度从高到低大致可以分为一、二、三、四等，精度越高的点数目越少。点数不够自然无法形成有效地控制网。因此，测量工作中需要通过已知的几个点，“一生二，二生三”，绘出一个可用于实际工程的控制网。“一生二，二生三”的过程需要已知点前后通视（即互相看得见），然而，在城市中，由于建筑物较多，视线不好，往往会出现已知点之间无法通视的现象。所以就引入了导线测量（traverse survey）。

导线测量又有很多分类，在此就不在赘述了。下面来介绍导线测量中的附和导线。所谓附和导线就是：“导线起始于一个已知控制点而终至于另一个已知控制点”。通过测量导线上未知点与已知控制点坐标的对比、计算，进行平差（控制误差的一种手段），得出一条导线上的各个待测点的坐标。这个过程就是导线测量的主体步骤了。导线示意图如下

根据导线示意图，这里有五点需要说明：

一、测绘中，坐标系很多，一般导线测量使用的是“平面直角坐标系”，之所以打上了双引

号，是因为，测量工作中的坐标系是以正北方向（竖轴）为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两点的高程差。当然实际运作中还要考虑地球曲率、大气温度等误差。在本文中，我对模型作了简化，对仪器高和地球曲率等因素做了默认设置。

图二

以上为本文要讨论的问题背景做了一个比较完整的描述下面进入问题的主体。

一、提出问题

有一测量小组通过查询获得了三个已知点的坐标，又通过野外测量，做了导线，并获得了导线中的转折角和边长，同时他们也获得了各个测站的α角。数据如下：

A列是边长，B列是左角值A列是α角值

现在要求获得各个待测点的坐标，并描绘出导线的形状，对比已知点得出测量值与实际获得坐标值的偏差。

二、对问题的分析和方案的产生

现在获得的数据已经可以画出导线了。问题是坐标系是不同的，那么可以进行转换。MATLAB数学坐标系的Y轴是正北方向，那么我们规定：Y轴北方向，X轴为东方向。化为极坐标时，顺时针为角的方向。

确定坐标系后，我们可以根据已知坐标把做边长和转折角转化为坐标值。计算方法：

α1=α2+α（转折角）－180°;

X1=X2＋D(边长) ×sin(α1);

Y1=Y2＋D(边长) ×cos(α1);

注：此处的坐标转化考虑到坐标系的转化，故sin、cos互换了。

通过以上计算获得了这些数据还是不够的。还需要用户输入已知点的坐标来对比测量结果。同时为了增加代码的通用性，可以由用户输入所需处理的文件名称，以及初始坐标值。同时，还可以通过用户选择，来实现计算和描绘平面的导线或者是三维的导线。通过用户选择那个函数，调用了二维描述和三维描述这两个函数。从而实现整个功能。

如上所示，整个框架已经构成了。

三、描述解决方案

首先给出可供用户选择的函数：

f unction begin ()

%这个函数给了用户选择。2D or 3D?

input('welcome to use my program￡?press Enter to star￡?')

N = input('choose 1 for 2D,choose 2 for 3D')

switch(N)

case 1

daoxian();

case 2

gaocheng();

end

再给出描画二维导线的函数描述：

function daoxian()

% this function is for 2D

FILENAME = input('what is the name of date file?')%获取用户所要处理的文件名

X0 =input('input the x：')

Y0 =input('input the y：')%用户输入已知坐标。

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,Y

GT(3)=CT(1)+CT(2)-pi;

I=4;

for k=1:4 %转化各个转角，使结果符合坐标系的运算

GT(I)=CT(I-1)+GT(I-1)-pi;

I=I+1;

end

for I=1:7

if GT(I)<0 %对小于0度的角进行转换

GT(I)=GT(I)+2*pi;

end

if GT(I)>2*pi %对大于360度的角进行转换

GT(I)=GT(I)-2*pi;

end

end

GT(1)=0;

GT(2)=0;

I=3;

for J=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;

end

OUTPUT = [X' ,Y'] % 输出数据

save OUTPUT %保存数据。

plot(X,Y,'b - O');%画图。

axis('equal');

xlabel('X');

ylabel('Y,北方向');

title('2D-----------map'); %以上为坐标的优化

grid on;

hold on;

% 对比已知点，求出坐标差，并描点。

CONx = input('please enter the x of control point：');

CONy = input('please enter the y of control point：');

plot(CONx,CONy,'r o'); %描出用于检核的点。

delx =X(7)- CONx %输出相差X,Y坐标大小。

dely =Y(7)- CONy

函数的描述和解释代码的后面已经阐述，现在直接给出三维的函数描述：

function gaocheng()

%This function is for 3D

FILENAME1 = input('what is the name of lengt-hdate file?')%获得用户所需要处理的文件名。FILENAME2 = input('what is the name of high-date file?')

X0 =input('input the x：')%获得初始的起算坐标

Y0 =input('input the y：')

Z0 =input('input the z：')

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;

for k=1:4 %计算获得可计算的坐标方位角

GT(I)=CT(I-1)+GT(I-1)-pi;

I=I+1;

end

for I=1:7

if GT(I)<0 %对小于0度的角进行转换

GT(I)=GT(I)+2*pi;

end

if GT(I)>2*pi %对大于360度的角进行转换

GT(I)=GT(I)-2*pi;

end

end

GT(1)=0;

GT(2)=0;

I=3;

for J=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;

end

HT=xlsread(FILENAME2);% 获得用户提供的高程α角。

HT=HT';

HT=(HT./(360))*2*pi;% α角转化；

LT=[];

for k=1:5

LT(k)=L(k+1);

end

T=LT.*tan(HT);%通过tan函数获得H的初始值。

Z=[];

Z(1)=0;

Z(2)=Z0;

for k=3:7 %获得H的值。

Z(k)=Z(k-1)+T(k-2);

end

OUTPUT = [X' ,Y',Z'] %输出计算结果

save OUTPUT;

plot3(X,Y,Z,'k - O'); %画出图像。

axis('equal');

xlabel('X');

ylabel('Y,北方向');

zlabel('H,高程');

grid on;

title('3D--------map');%以上为坐标的优化

hold on;

%获得已知点，进行比较，描绘，输出结果

CONx = input('please enter the x of control point：');

CONy = input('please enter the y of control point：');

CONz = input('please enter the z of control point：');

plot3(CONx,CONy,CONz,'r o') %描出用于检核的点。

delx=X(7)- CONx % 输出相差的坐标值

dely=Y(7)- CONy

delyz=Z(7)- CONz

现在整个函数的实现已经给出来了。其中，function begin ()放在一个M文件里，function daoxian()放在一个M文件里，function gaocheng()放在一个M文件里。

其实，第2D和3D就差了一个坐标元素。主体函数还是差不多的。不过，为了表达清楚，我还是把它们都写出来了。我在写这些函数的时候，遇到了一些矩阵数组计算的小问题，最后虽然得到了解决，但我一直怀疑是不是有更好的描述方法。

四、运行MATLAB,获得结果

首先，我要说明的是，本来我是想在word里直接运行代码的，但是，由于不同的M文件，这样的表达，我觉得不是很清晰（可能我方法不对），所以就选择了截图运行过程的方法。

运行时首先运行begin函数，选1，再按要求输入所要求调用的文件名：导线数据.xlsl，再输入起算点的X,Y坐标，此处是100,200。运行后，再根据提示输入控制点的坐标进行检核，

这里是（-36,350）。整个输入过程就结束了。

运行过程如下：

启动界面

下图显示的是输出的X,Y的坐标值

下图显示的delx，dely是与控制点相差的数值

下图显示的是导线的图像，其中红点为检核所用的控制点（此处几乎重合）

现在来看看三维导线的运行。首先调用begin函数，选2，按要求输入文件名，再输入起算点的坐标，这里是（100,200,300）；运行如下

即输出了计算结果，然后按要求说如已知控制点的坐标用于检核。这里输入（170,390,640），则其输出了测算点与控制点差值。结果如下：

最后我们可以看一下三维图像的运行结果（如下图），注意，红点是用于检核的已知坐标的点，不得不提一下的是，实际操作中，精度往往比较高，这里为了突出误差，便于观察和理解图像，把误差放大了。

运行到这里可以说所有的结果都获得了，问题得到了解决。

五、总结

这个问题的解决其实很不完美。首先，这个模型不是实际的模型，我对其中的数据进行了优化（如角度和长度的位数限制），并且祛除了一些比较复杂的能够减少误差的计算，只保留基本的几种思想方法，对导线的测量进行了描绘。但是问题的全貌，和数据的主要处理方法都到位了。这个问题当然可以拓展，比如，直接画出初步平差后的图像。这个问题就复杂多了，我这些天的算法设计总是有问题，所以就不好拿出来了。

另外，这个程序可重用性不高，虽然能实现对不同数据源文件的读取，用户自主输入起算和检核坐标，但是，如果测站数不一样，还需要修改函数内部的部分数据。这个方面的问题，我目前还不知道如何解决更好。我以后会边学习，边思考，找出好方法。

还要说明的是，对于导线的描绘，测绘有专门的软件如CASS等，其对于专业内的绘图功能自然比MATLAB好些。此处的学习意义大于运用意义。

最后不得不说一下我对这门课的感受，因为MATLAB这门选修课是我修过的“最有技术含量”的公选课。只可惜我目前只学了半个学期的专业课，所以很多专业的应用我都不能“立即”去开发。但是通过这门课的学习和这次结课作业的练习，我认识了MATLAB这种工具,知道了他的强悍功能。同时，我还发现了和自己专业相关的工具箱，如“Mapping Toolbox”等。通过互联网，我还了解到MATLAB在“GPS”的计算还有“平差”的计算中有较多的的应用。我想以后有了机会，具备了能力，我是会尝试一下各种应用的。

参考文献

《MATLAB7及工程问题解决方案》（美）Delores M.Etter 邱李华译2006.3 《MATLAB基础及其应用》北京大学周开利邓春晖2009.11

MATLAB7.6.0帮助

《数字测图原理与方法》武汉大学2010 .7

参考网站

MATLAB中文论坛：https://www.360docs.net/doc/7116252685.html,

MATLAB百度贴吧：https://www.360docs.net/doc/7116252685.html,/f?kw=matlab

MATLAB教程网：https://www.360docs.net/doc/7116252685.html,

编程爱好者MATLAB讨论区：https://www.360docs.net/doc/7116252685.html,/club/bbs.asp?bbsid=57

2010年11月12日星期五