基于matlab的測量導(dǎo)線描繪

1、測繪學(xué)院 學(xué)號(hào):2009301610208 夏少波 手機(jī): 基于MATLAB的測量導(dǎo)線描繪摘要: 本文從導(dǎo)線的描繪出發(fā)，提出了該問題的MATLAB解決方法。給出了數(shù)據(jù)的處理結(jié)果和圖表。并且，設(shè)計(jì)了二維和三維的代碼描述。其中對(duì)待處理的數(shù)據(jù)有一些優(yōu)化，但計(jì)算方法和計(jì)算的邏輯是沒有問題的。最后，MATLAB對(duì)這個(gè)問題的解決也較為令人滿意。關(guān)鍵字: 控制測量 導(dǎo)線 代碼 運(yùn)行結(jié)果 MATLAB任何測量工作都會(huì)產(chǎn)生誤差，而誤差在絕對(duì)意義上是不可避免的，所以在測量工作中必須采用一定的方法來控制誤差的范圍，以期達(dá)到工作需求。測量工作中遵循“先整體后局部，先控制后碎部”的原則。簡而言之就是在工作區(qū)域內(nèi)建立一定

2、精度的由一定數(shù)目的已知坐標(biāo)點(diǎn)構(gòu)成的幾何圖形，然后在幾何圖形內(nèi)進(jìn)行測量。在測量過程中，我們可以通過各個(gè)已知點(diǎn)的坐標(biāo)來驗(yàn)算測得的未知點(diǎn)的坐標(biāo)。這樣通過不斷地檢核、測量、檢核，成果的誤差就可以控制在一定范圍內(nèi)。其中，已知點(diǎn)被稱作“控制點(diǎn)”，幾何圖形被稱為“控制網(wǎng)”，這種測量方法被稱為“控制測量”。一般意義上的控制測量可以分為“平面控制測量”和“高程控制測量”。平面控制測量是基平面坐標(biāo)系X-Y，而高程控制測量是基于高程H的。目前，隨著測繪技術(shù)的發(fā)展，也出現(xiàn)了把兩種控制測量結(jié)合起來的“三維控制測量”?？刂凭W(wǎng)根據(jù)其精度從高到低大致可以分為一、二、三、四等，精度越高的點(diǎn)數(shù)目越少。點(diǎn)數(shù)不夠自然無法形成有效地控

3、制網(wǎng)。因此，測量工作中需要通過已知的幾個(gè)點(diǎn)，“一生二，二生三”，繪出一個(gè)可用于實(shí)際工程的控制網(wǎng)?！耙簧钡倪^程需要已知點(diǎn)前后通視（即互相看得見），然而，在城市中，由于建筑物較多，視線不好，往往會(huì)出現(xiàn)已知點(diǎn)之間無法通視的現(xiàn)象。所以就引入了導(dǎo)線測量（traverse survey）。導(dǎo)線測量又有很多分類，在此就不在贅述了。下面來介紹導(dǎo)線測量中的附和導(dǎo)線。所謂附和導(dǎo)線就是：“導(dǎo)線起始于一個(gè)已知控制點(diǎn)而終至于另一個(gè)已知控制點(diǎn)”。通過測量導(dǎo)線上未知點(diǎn)與已知控制點(diǎn)坐標(biāo)的對(duì)比、計(jì)算，進(jìn)行平差（控制誤差的一種手段），得出一條導(dǎo)線上的各個(gè)待測點(diǎn)的坐標(biāo)。這個(gè)過程就是導(dǎo)線測量的主體步驟了。導(dǎo)線示意圖如下根

4、據(jù)導(dǎo)線示意圖，這里有五點(diǎn)需要說明：一、 測繪中，坐標(biāo)系很多，一般導(dǎo)線測量使用的是“平面直角坐標(biāo)系”，之所以打上了雙引號(hào)，是因?yàn)?，測量工作中的坐標(biāo)系是以正北方向（豎軸）為X方向，橫軸為Y方向，極坐標(biāo)中順時(shí)針旋轉(zhuǎn)方向?yàn)檎较?。二?圖中2、3等點(diǎn)要安放儀器進(jìn)行測量，稱這些點(diǎn)為測站。三、 A、B、C、D為已知點(diǎn)。P2、P3等為待測點(diǎn)。導(dǎo)線測量中需要計(jì)算出各個(gè)點(diǎn)的坐標(biāo)。四、 關(guān)于P1左、P2左等角度，為測繪中所謂的左角。所謂的左角可以通過導(dǎo)線的前進(jìn)方向來確定。比如，在這張圖中，測量的順序?yàn)锳-B-C-D,那么前進(jìn)方向即為A-B-C-D,而前進(jìn)方向的左手邊的角度就是左角了（右角亦然，并不規(guī)定一定要測哪個(gè)

5、角，計(jì)算方法略有不同而已）。這些角一般稱作轉(zhuǎn)折角。五、 平面的導(dǎo)線控制測量需要在一測站測相鄰的左右測站之間的夾角（如圖中的P1），該測站到前進(jìn)方向上的下一測站的距離（如圖中的D(1-2)，也稱作邊長）。如果把高程算進(jìn)去，還需要獲得兩測站之間的高差。而獲得高差數(shù)據(jù)的方法很多，一般可以通過水準(zhǔn)測量或三角高程測量獲得。在這里介紹一下三角高程的基本原理。三角高程測量的基本原理就是利用三角形各邊角關(guān)系來獲得所需數(shù)據(jù)。如圖二;已知A、B兩點(diǎn)的水平距離D，和一個(gè)夾角，通過計(jì)算D×tan()，再通過加減測量儀器高就可以獲得A、B兩點(diǎn)的高程差。當(dāng)然實(shí)際運(yùn)作中還要考慮地球曲率、大氣溫度等誤差。在本文中，

6、我對(duì)模型作了簡化，對(duì)儀器高和地球曲率等因素做了默認(rèn)設(shè)置。圖二以上為本文要討論的問題背景做了一個(gè)比較完整的描述下面進(jìn)入問題的主體。一、 提出問題有一測量小組通過查詢獲得了三個(gè)已知點(diǎn)的坐標(biāo)，又通過野外測量，做了導(dǎo)線，并獲得了導(dǎo)線中的轉(zhuǎn)折角和邊長，同時(shí)他們也獲得了各個(gè)測站的角。數(shù)據(jù)如下：A列是邊長， B列是左角值 A列是角值現(xiàn)在要求獲得各個(gè)待測點(diǎn)的坐標(biāo)，并描繪出導(dǎo)線的形狀，對(duì)比已知點(diǎn)得出測量值與實(shí)際獲得坐標(biāo)值的偏差。二、 對(duì)問題的分析和方案的產(chǎn)生現(xiàn)在獲得的數(shù)據(jù)已經(jīng)可以畫出導(dǎo)線了。問題是坐標(biāo)系是不同的，那么可以進(jìn)行轉(zhuǎn)換。MATLAB數(shù)學(xué)坐標(biāo)系的Y軸是正北方向，那么我們規(guī)定：Y軸北方向，X軸為東方向。化

7、為極坐標(biāo)時(shí)，順時(shí)針為角的方向。 確定坐標(biāo)系后，我們可以根據(jù)已知坐標(biāo)把做邊長和轉(zhuǎn)折角轉(zhuǎn)化為坐標(biāo)值。計(jì)算方法： 1=2+（轉(zhuǎn)折角）180°X1=X2D(邊長) ×sin(1);Y1=Y2D(邊長) ×cos(1);注：此處的坐標(biāo)轉(zhuǎn)化考慮到坐標(biāo)系的轉(zhuǎn)化，故sin、cos互換了。 通過以上計(jì)算獲得了這些數(shù)據(jù)還是不夠的。還需要用戶輸入已知點(diǎn)的坐標(biāo)來對(duì)比測量結(jié)果。同時(shí)為了增加代碼的通用性，可以由用戶輸入所需處理的文件名稱，以及初始坐標(biāo)值。同時(shí)，還可以通過用戶選擇，來實(shí)現(xiàn)計(jì)算和描繪平面的導(dǎo)線或者是三維的導(dǎo)線。通過用戶選擇那個(gè)函數(shù)，調(diào)用了二維描述和三維描述這兩個(gè)函數(shù)。從而實(shí)現(xiàn)整個(gè)

8、功能。 如上所示，整個(gè)框架已經(jīng)構(gòu)成了。三、 描述解決方案首先給出可供用戶選擇的函數(shù)：function begin ()%這個(gè)函數(shù)給了用戶選擇。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再給出描畫二維導(dǎo)線的函數(shù)描述：function daoxian()

9、% this function is for 2D FILENAME = input('what is the name of date file?')%獲取用戶所要處理的文件名X0 =input('input the x：')Y0 =input('input the y：')%用戶 輸入 已知坐標(biāo)。Angle = atan(Y0/X0); %用戶輸入已知的坐標(biāo)點(diǎn)后，進(jìn)行角度的轉(zhuǎn)化，使之能在坐標(biāo)系中正常實(shí)現(xiàn)。A=xlsread(FILENAME); %獲取用戶所要處理的文件內(nèi)容L=A(:,1); L=L' CT=A(:,2); CT=CT

10、'%以上是對(duì)獲取數(shù)據(jù)的處理，使之符合矩陣運(yùn)算。GT=;CT(1)=Angle;CT=(CT./(360)*2*pi;%把一般測得的360制角轉(zhuǎn)化為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 %轉(zhuǎn)化各個(gè)轉(zhuǎn)角，使結(jié)果符合坐標(biāo)系的運(yùn)算GT(I)=CT(I-1)+GT(I-1)-pi; I=I+1;endfor I=1:7if GT(I)<0 %對(duì)小于0度的角進(jìn)行轉(zhuǎn)換GT(I)=GT(I)+2*pi;endif GT(I)>2*pi %對(duì)大于360度的角進(jìn)行轉(zhuǎn)換

11、 GT(I)=GT(I)-2*pi;endendGT(1)=0;GT(2)=0;I=3;for J=1:5 %核心計(jì)算，獲得坐標(biāo)。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' % 輸出數(shù)據(jù)save OUTPUT %保存數(shù)據(jù)。plot(X,Y,'b - O');%畫圖。axis('equal');xlabel('X');ylabel('Y,北方向');title('2D-map'

12、); %以上為坐標(biāo)的優(yōu)化grid on;hold on;% 對(duì)比已知點(diǎn)，求出坐標(biāo)差，并描點(diǎn)。 CONx = input('please enter the x of control point：');CONy = input('please enter the y of control point：'); plot(CONx,CONy,'r o'); %描出用于檢核的點(diǎn)。 delx =X(7)- CONx %輸出相差X,Y坐標(biāo)大小。dely =Y(7)- CONy函數(shù)的描述和解釋代碼的后面已經(jīng)闡述，現(xiàn)在直接給出三維的函數(shù)描述：function g

13、aocheng()%This function is for 3DFILENAME1 = input('what is the name of lengt-hdate file?')%獲得用戶所需要處理的文件名。FILENAME2 = input('what is the name of high-date file?')X0 =input('input the x：')%獲得初始的起算坐標(biāo)Y0 =input('input the y：')Z0 =input('input the z：')Angle = atan(

14、Y0/X0); %計(jì)算出起算的兩點(diǎn)間在坐標(biāo)系中的夾角。A=xlsread(FILENAME1);%獲得用戶提供的數(shù)據(jù)L=A(:,1); L=L' CT=A(:,2); CT=CT' %以上是對(duì)獲取數(shù)據(jù)的處理，使之符合矩陣運(yùn)算。GT=;CT(1)=Angle; %用戶輸入已知的坐標(biāo)點(diǎn)后，進(jìn)行角度的轉(zhuǎn)化，使之能在坐標(biāo)系中正常實(shí)現(xiàn)。CT=(CT./(360)*2*pi;%轉(zhuǎn)化角度 X=0,X0,0,0,0,0,0;%初始化坐標(biāo)值 Y=0,Y0,0,0,0,0,0;CT(1)=Angle; GT(3)=CT(1)+CT(2)-pi;I=4;for k=1:4 %計(jì)算獲得可計(jì)算的坐標(biāo)方位

15、角GT(I)=CT(I-1)+GT(I-1)-pi; I=I+1;endfor I=1:7if GT(I)<0 %對(duì)小于0度的角進(jìn)行轉(zhuǎn)換GT(I)=GT(I)+2*pi;endif GT(I)>2*pi %對(duì)大于360度的角進(jìn)行轉(zhuǎn)換 GT(I)=GT(I)-2*pi;endend GT(1)=0;GT(2)=0;I=3;for J=1:5 %核心計(jì)算，獲得坐標(biāo)值。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&#

16、39;HT=(HT./(360)*2*pi;% 角轉(zhuǎn)化；LT=;for k=1:5 LT(k)=L(k+1);endT=LT.*tan(HT);%通過tan函數(shù)獲得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' %輸出計(jì)算結(jié)果save OUTPUT;plot3(X,Y,Z,'k - O'); %畫出圖像。axis('equal');xlabel('X');ylabel('Y,北方向');z

17、label('H,高程');grid on;title('3D-map');%以上為坐標(biāo)的優(yōu)化hold on;%獲得已知點(diǎn)，進(jìn)行比較，描繪，輸出結(jié)果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

18、 o') %描出用于檢核的點(diǎn)。delx=X(7)- CONx % 輸出相差的坐標(biāo)值dely=Y(7)- CONydelyz=Z(7)- CONz現(xiàn)在整個(gè)函數(shù)的實(shí)現(xiàn)已經(jīng)給出來了。其中，function begin ()放在一個(gè)M文件里，function daoxian()放在一個(gè)M文件里，function gaocheng()放在一個(gè)M文件里。其實(shí)，第2D和3D就差了一個(gè)坐標(biāo)元素。主體函數(shù)還是差不多的。不過，為了表達(dá)清楚，我還是把它們都寫出來了。我在寫這些函數(shù)的時(shí)候，遇到了一些矩陣數(shù)組計(jì)算的小問題，最后雖然得到了解決，但我一直懷疑是不是有更好的描述方法。四、 運(yùn)行MATLAB,獲得結(jié)果首

19、先，我要說明的是，本來我是想在word里直接運(yùn)行代碼的，但是，由于不同的M文件，這樣的表達(dá)，我覺得不是很清晰（可能我方法不對(duì)），所以就選擇了截圖運(yùn)行過程的方法。運(yùn)行時(shí)首先運(yùn)行begin函數(shù)，選1，再按要求輸入所要求調(diào)用的文件名：導(dǎo)線數(shù)據(jù).xlsl，再輸入起算點(diǎn)的X,Y坐標(biāo)，此處是100,200。運(yùn)行后，再根據(jù)提示輸入控制點(diǎn)的坐標(biāo)進(jìn)行檢核，這里是（-36,350）。整個(gè)輸入過程就結(jié)束了。運(yùn)行過程如下：啟動(dòng)界面下圖顯示的是輸出的X,Y的坐標(biāo)值下圖顯示的delx，dely是與控制點(diǎn)相差的數(shù)值下圖顯示的是導(dǎo)線的圖像，其中紅點(diǎn)為檢核所用的控制點(diǎn)（此處幾乎重合）現(xiàn)在來看看三維導(dǎo)線的運(yùn)行。首先調(diào)用begin

20、函數(shù)，選2，按要求輸入文件名，再輸入起算點(diǎn)的坐標(biāo)，這里是（100,200,300）；運(yùn)行如下即輸出了計(jì)算結(jié)果，然后按要求說如已知控制點(diǎn)的坐標(biāo)用于檢核。這里輸入（170,390,640），則其輸出了測算點(diǎn)與控制點(diǎn)差值。結(jié)果如下：最后我們可以看一下三維圖像的運(yùn)行結(jié)果（如下圖），注意，紅點(diǎn)是用于檢核的已知坐標(biāo)的點(diǎn)，不得不提一下的是，實(shí)際操作中，精度往往比較高，這里為了突出誤差，便于觀察和理解圖像，把誤差放大了。 運(yùn)行到這里可以說所有的結(jié)果都獲得了，問題得到了解決。五、 總結(jié)這個(gè)問題的解決其實(shí)很不完美。首先，這個(gè)模型不是實(shí)際的模型，我對(duì)其中的數(shù)據(jù)進(jìn)行了優(yōu)化（如角度和長度的位數(shù)限制），并且祛除了一些比較復(fù)雜的能夠減少誤差的計(jì)算，只保留基本的幾種思想方法，對(duì)導(dǎo)線的測量進(jìn)行了描繪。但是問題的全貌，和數(shù)