畢業(yè)設計（論文）基于MATLAB的水準網(wǎng)和測邊網(wǎng)平差程序設計

1、基于 matlab 的水準網(wǎng)和測邊網(wǎng)平差程序設計 摘 要 matlab 是目前在研究機構(gòu)廣泛應用的一種數(shù)值計算及圖形工具軟件，它的 特點是語法結(jié)構(gòu)簡明、數(shù)值計算高效、圖形功能完備，特別適合非專業(yè)編程員完成 數(shù)值計算、科學試驗處理等任務。以往的測量數(shù)據(jù)處理方法需要編制特定的處理矩 陣運算程序，而且程度復雜，難度大。 本文介紹一種基于 matlab 的水準網(wǎng)和測邊網(wǎng)的程序設計方法，與其它算法 語言相比，具有編程簡單，運算速度快的特點。文中分別闡述了水準網(wǎng)和測邊網(wǎng)程 序的理論基礎、實現(xiàn)步驟和運行結(jié)果。通過實例的分析，總結(jié)出利用 matlab 對 測量數(shù)據(jù)處理有很大的應用價值，它縮短了編程的時間，提高

2、工作效率。 關(guān)鍵詞：matlab； 水準網(wǎng)； 測邊網(wǎng)； 程序設計 abstract matlab is one species of numerical-values calculation and graphic tools software which is widely used to apply at research institutions at present. the particularities are: concise grammar-structure、highly efficient in numerical values calculating、complete f

3、unction of graphs、especially it is adapted to evildoing professional programmer to accomplish the tasks that are numerical-values calculating and scientific experiments treating. the ancient methods of measured data-processing need establishing special proceedings of treating matrices operation, mor

4、eover, it is complex and greatly difficult. this article introduces one programming method dealing with leveling and measuring edge network based on matlab. compared with other algorithm language, it has particularities which are simply programming and quickly operating. the article separately expat

5、iate the theories basics、realizing steps and running results at leveling and measuring edge network. with the analysis of examples, it has prodigious application value in measured data-processing by use of matlab. moreover, it shortens programming time and improves working effectiveness. key words：m

6、atlab；leveling network；measuring edge network；programming 目目 錄錄 緒論緒論.4 4 1.1. matlabmatlab 軟件簡介軟件簡介 .5 5 2 2matlabmatlab 在測量平差中在測量平差中的的應應用用.6 6 2.1 測量平差原理的概述 .6 2.2 平差程序總體方案 .7 3 3. . 水準網(wǎng)平差程序水準網(wǎng)平差程序.8.8 3.1 程序的功能 .8 3.2 水準模型網(wǎng)的間接平差 .8 3.2.1 “權(quán)”值的確定 .8 3.2.2 水準路線的平差計算 .9 3.2.3 精度評定 .11 3.3 水準網(wǎng)間接平差程序信息

7、設計 .11 3.4水準網(wǎng)程序與使用說明.12 3.4.1 水準網(wǎng)程序流程圖 .12 3.4.2水準網(wǎng)程序的使用 .12 3.5 案例 .13 4.4. 測邊網(wǎng)平差程測邊網(wǎng)平差程序序設計設計.1616 4.1 數(shù)學模型 .15 4.1.1 誤差方程和法方程的組成 .15 4.1.2邊長觀測的權(quán) .15 4.1.3 解算法方程 .16 4.1.4精度評定 .19 4.2測邊網(wǎng)平差信息設計.20 4.2.1 主要的技術(shù)要求 .21 43 利用 matlab 的繪圖語句繪制網(wǎng)圖 .21 4.44.4 測邊網(wǎng)程序測邊網(wǎng)程序和和使使用說明用說明 .22 4.5程序代碼說明：.23 4.6 程序的使用算例

8、 .25 結(jié)結(jié) 論論 .29 致致 謝謝.3030 參考文獻參考文獻.3131 附錄附錄一一.3232 附錄附錄二二.3636 附錄三附錄三.4646 緒論 作為一名測量技術(shù)人員，如果不掌握一門 pc 機編程語言與便攜計算工具，要 想提高測量工作的效率幾乎寸步難行。測量需求的多樣性與復雜性，造就了測量計 算鮮明的個性化特點，這就是在商業(yè)測量計算軟件高度發(fā)達的今天，掌握一種實用 的程序語言進行編程計算仍有廣泛的市場需求的重要原因。 當今較流行的計算機程序語言基本上都是基于 windows 的，例如 turbo pascal,visual basic,visual c,borland c+等，這些

9、程序語言的優(yōu)勢是基于對象及可利 用 windows 豐富的系統(tǒng)資源，應用它們可以開發(fā)出界面非常豐富和友好的應用程序， 其劣勢主要有以下幾點： 1.windows 程序都非常龐大，學習并熟練掌握它們并非易事。 2.雖然市場上已有的多種專用的測量平差軟件都是采用 c 語言開發(fā)的，但這些軟 件價格都比較貴，而且都帶有加密狗，一次只能供一個用戶使用。出于商業(yè)目的， 開發(fā)商不會公開程序源代碼，這為修改程序功能以適應用戶的特殊需求帶來了不便。 3.在測量生產(chǎn)中，經(jīng)常需要根據(jù)工程的實際情況進行一些個性化的數(shù)值計算工作， 這些數(shù)值計算工作無固定模式，這就需要求測量技術(shù)人員最好能熟練掌握一種適用 于數(shù)值計算的程

10、序語言，以便提高測量計算的效率。 4.c 語言的數(shù)值計算語句不夠豐富，例如，在測量平差計算中，經(jīng)常需要進行的 矩陣運算，尤其是解法方程的矩陣求逆不能直接使用語句實現(xiàn)，而必須應用計算機 算法編程實現(xiàn)。 如果不是基于商業(yè)軟件開發(fā)，只為滿足實際測量工作計算需要，則 c 語言的劣 勢就變成了 matlab 語言的優(yōu)勢。 1. matlab軟件簡介 matlab是從matrix(矩陣)和laboratory(實驗室)各取前3個字母組成的,意思是 矩陣實驗室,是美國mathworks公司于20世紀80年代中期推出的一種交互式、面向?qū)?象的科技應用軟件,是一個為科學和工程計算而專門設計的高級交互式軟件包。

11、matlab 集成了圖示與精確的數(shù)值計算,是一個可以完成各種計算和數(shù)據(jù)可視化的強 有力工具,其優(yōu)秀的數(shù)值計算能力和卓越的數(shù)據(jù)可視化能力使其很快在數(shù)學軟件中脫 穎而出,成為以矩陣運算為主要工作方式的線性代數(shù)、概率論和數(shù)理統(tǒng)計、自動控制、 數(shù)字信號處理、動態(tài)系統(tǒng)仿真等領(lǐng)域教學和科研工作者的有力武器。隨著該軟件自 身的發(fā)展及市場的需求,其功能日趨完善,其最高版本7.0版已經(jīng)推出,隨著版本的不 斷升級,它的數(shù)值計算及符號計算功能得到了進一步完善。 matlab是以矩陣作為數(shù)據(jù)操作的基本單位,矩陣的生成、運算、轉(zhuǎn)置、求逆等非 常簡單。在matlab環(huán)境中,不需要對創(chuàng)建的變量對象給出類型說明和維數(shù),所有的

12、變 量都作為雙精度數(shù)來分配內(nèi)存空間,matlab將自動地為每一個變量分配內(nèi)存。matlab 語言起源于矩陣運算,并已經(jīng)發(fā)展成為一種高度集成的計算機語言,它提供了強大的 科學運算、靈活的程序設計流程、高質(zhì)量的圖形可視化與界面設計、便捷的與其他 程序和語言接口的功能。matlab系統(tǒng)主要包含5 部分的內(nèi)容:matlab 工作環(huán)境、 mablab 數(shù)學函數(shù)庫、matlab語言體系、句柄圖形、matlab應用程序接口(api)。 matlab系統(tǒng)主要功能包括:數(shù)值計算功能、符號計算功能、數(shù)據(jù)分析和可視化、文字 處理功能、simulink動態(tài)仿真功能。同時,matlab又是開放的,除了內(nèi)部函數(shù)之外,所

13、有的matlab 主包文件和各工具包文件都是可讀可改的源文件,用戶可以作為參考掌 握其用法,并可對其修改以適應自己的需要,也可加入自己編寫的文件構(gòu)成新的工具 包。例如,隨著gps 的廣泛應用,orion dynamics and con2t rol corporation、constell inc. gpssoft llc、navsyscorporation等多家公司都相應 開發(fā)出了適于gps數(shù)據(jù)處理的matlab 工具箱。 matlab 是一個集數(shù)值計算、圖形管理、程序開發(fā)于一體的功能十分強大的系統(tǒng)。 將 matlab 應用于測量數(shù)據(jù)的處理是一件非常有意義的工作。mo2hamed 等曾成功地

14、 在 matlab 系統(tǒng)中利用白濾波技術(shù)研究動態(tài)解算 gps 載波相位信號的模糊度問題。因 為測量數(shù)據(jù)的處理特別是測量平差主要應用矩陣運算,而 matlab 又特別易于做矩陣 運算,因此,研究開發(fā)基于 matlab 的測量平差方法具有極好的應用價值。 2matlab 在測量平差中的應用 測量平差數(shù)據(jù)處理主要是基于矩陣的運算,常用的矩陣運算主要是矩陣的生成、 轉(zhuǎn)置、求逆和矩陣求廣義逆等。在matlab環(huán)境中,不需要對創(chuàng)建的變量對象給出類型 說明和維數(shù),所有的變量都作為matlab中的m文件的語法與其他的高級語言類似,是一 種程序化的編程語言,同時也是一種解釋性的編程語言,即逐行解釋運行程序,使程

15、序 容易調(diào)試,計算更為簡捷,而且對于平差原理理解和掌握變得更容易。另外,matlab語 言與數(shù)學語言比較接近,更容易掌握和理解。 2.1測量平差原理的概述 測量平差的函數(shù)模型有條件方程和觀測方程。以條件方程為函數(shù)的模型的最小 二乘平差稱為條件平差；在條件方程中，根據(jù)需要如果還設有一定數(shù)量的未知數(shù)， 則稱為附有參數(shù)的條件平差；以觀測方程為函數(shù)模型的最小二乘平差稱為間接平差； 如果觀測方程中的某些參數(shù)不獨立，則這些不獨立參數(shù)必然存在一些條件，稱這種 平差模型為附有條件的間接平差。本文的兩個程序都采用間接平差模型。 對于一個實際平差問題,根據(jù)所選參數(shù)的個數(shù)、選什么量為參數(shù)以及參數(shù)之間是 否函數(shù)獨立,

16、經(jīng)過仔細推敲可以發(fā)現(xiàn)附有條件的間接平差模型本身就是各種經(jīng)典平差 模型的概括模型,其余的經(jīng)典平差模型,如條件平差模型、間接平差模型、附有未知 數(shù)的條件平差模型和附有限制條件的條件平差模型都是它的特例。 間接平差的公式匯集： 間接平差模型為 min t vbxl v pv （21） 系數(shù)矩陣 b 滿秩，即 rank(b)=t 法方程及解為： （2- 1 0,(,) tt bbebbe nxfnb pb fb pl 2） （2- 1 bbe xnf 3） 參數(shù)的平差值： （2- 0 xxx 4） 觀測量的平差值： （2- llv 5） 單位權(quán)中誤差： （2- 0 t v pv nt 6） 平差參數(shù)的

17、協(xié)方差陣： （2- 21 0bb x x dn 7） 平差函數(shù)的協(xié)方差陣： （2- 21 0 t bb qf nf 8） 2.2 平差程序總體方案 matlab 號稱為全球工程師的共同語言，其語法和 c 語言相似，但它有強大的數(shù) 值計算和繪圖功能，這使之在工程應用方面的計算更出色，本文就基于這種程序設 計語言環(huán)境設計一個控制網(wǎng)平差程序。該程序包含了一個高程控制網(wǎng)平差程序和測 邊網(wǎng)平差程序。 本程序適用于各種等級的高程網(wǎng)和測邊網(wǎng)，程序在設計過程中，始終考慮數(shù)據(jù) 的儲存量。因而本程序不儲存誤差方程的系數(shù)和常數(shù)項，對待定點數(shù)較多的平差網(wǎng)， 組成法方程的系數(shù)矩陣是個稀疏矩陣，如待定點的編號恰當，法方程

18、的系數(shù)會集中 在主元系數(shù)的兩側(cè)形成帶狀。為減少法方程系數(shù)的儲存量，只要按行儲存下三角陣 或按列儲存上三角陣中第一個非零系數(shù)起的系數(shù)，就是通常叫做維變帶寬儲存方法。 3. 水準網(wǎng)平差程序 3.1 程序的功能 本程序適用于二、三、四等水準網(wǎng)平差計算，平差的水準網(wǎng)可以是獨立的、也 可以是附合網(wǎng)，其主要功能是完成水準網(wǎng)的平差計算和精度評定計算。平差計算采 用間接平差法，以歸心的觀測值為高差，以未知點高程為未知參數(shù)。精度評定計算 包括計算單位權(quán)中誤差和每個待定點的高程中誤差。 3.2 水準模型網(wǎng)的間接平差 3.2.1 “權(quán)”值的確定 當在相同的條件下進行水準測量時,其精度是相同的,因而觀測結(jié)果的可靠性也

19、 是同樣的。但如果在不同的條件下進行水準測量時,高程的精度就有所不同,此時稱 為不等精度觀測,所求出的未知量的值、高程的最或是值并對其精度進行評定時,就 需要“權(quán)”了。 由于觀測的不等精度,因而觀測值的可靠程度不同,求未知量的最或是值時,這樣 的一個因素就必須考慮了,這個因素是:可靠性大的某觀測值,其精度高,對測量的最 后結(jié)果的影響也就越大。此時用“權(quán)”值來表示觀測值的可靠程度,那么,“權(quán)”值 愈大,觀測值的可靠程度就愈高。另外,在觀測過程中,觀測值的中誤差愈小,觀測結(jié) 果愈可靠,它的“權(quán)”值就愈大。因而,根據(jù)中誤差來確定“權(quán)”值是非常適當?shù)摹?設以pi 表示觀測值li 的“權(quán)”,m 為中誤差

20、,則“權(quán)”值的定義為: (3-1) 式中:a 為任意的正常數(shù),在一組觀測值中為一個定數(shù)。 在實際測量中,通常是觀測值的中誤差事先并不知道,因而必須先確定觀測值的“權(quán)” ,然后才能求出未知量的最或是值。此時可以利用距離(s)或測站數(shù)(n)來確定觀測值 高程的“權(quán)”。 根據(jù)偶然誤差傳播定律,各觀測點高程hi的中誤差mi由測站數(shù)ni確定時,則有: (3-2) ii mm n 式中:m 為一組觀測值的中誤差,為一個定數(shù). 由(3-2-1)、(3-2-2)兩式可得: (3-3) 同樣可得出: （3- 4） 式中:c為定數(shù),為測距. i s 由(3-3)、（3-4）兩式可以得出這樣一個結(jié)論:當測站觀測高差

21、等精度時,觀測總高 差的“權(quán)”與測站數(shù)或距離成反比。 3.2.2 水準路線的平差計算 1附合路線的平差計算 假定在圖1 示的a 、b 兩水準點之間布設一條水準路線,a 、b 兩水準點的高程 為已知,分別設為、 、c為中間水準點。假定觀測了所有的點的高 a h b h 1 n 2 n 程,現(xiàn)擬求c點的高程的最或是值。 c h 可由水準路線a c、b c分別觀測的高差、計算得出,由此而 c h ac h bc h 得到的觀測高程分別設為hc1、hc2,其值為: hc1=+;hc2=+ a h ac h b h bc h 當hc1、hc2在不等精度條件下觀測得出時,它們的“權(quán)”也不同,分別設為 pc

22、1、pc2,這樣c點的高程的最或是值為: c h (3-5) 1122 12 cccc c cc p hp h h pp 根據(jù)a 點的高程,a c水準路線觀測的高差以及b c水準路線觀測 a h ac h 的高差,可推算出b點的觀測高程為: bc h b h =+- b h a h ac h bc h 水準路線a b 的高程閉合差為: =-= (3-6) h f b h b h 12cc hh 由（3-6)式得到: =- 2c h 1c h h f 由（3-3)式得到: 、 （、分別表示水準路線 1c ac c p n 2c bc c p n ac n bc n ac、bc 的測站數(shù),水準路線

23、ab的測站數(shù)) abacbc nnn 將上述表達式代入(3-2-5)式中,得到: (3-7) 如果以水準路線ac的距離、b c的距離、ab的距離 ( bc s bc s ab s )來確定高程觀測值的“權(quán)”值時,同樣可以得到: abacbc sss (3-8) 圖3-1 水準路線圖 2閉合路線的平差計算 閉合路線的平差計算原理與附合路線相同,因而(3-7)、(3-8)兩式的結(jié)論適用于 閉合路線的平差計算。 (3) 具有一個結(jié)點的水準網(wǎng)的平差計算 如圖2所示為具有一個結(jié)點的水準網(wǎng),b,c,d,為已知高程水準點,ba,ca,d a,為水準路線,則接點a的高程最或是值為: (3-9) 1122331

24、 123 1 n aiai aaaaaai an aaa ai i p h p hp hp h h ppp p 式中分別為水準路線ba,ca,da,計算a的觀測高程,各高 123 , aaa hhh 程相應的“權(quán)”值為 123 , aaa ppp 設的算術(shù)平均值為,各高程觀測值與的差值分別為 123 , aaa hhh 0 a h 0 b h a1,a2,a3,則有: (3-10) 將(3-10)式代入(3-9)式得到: (3-11) 當以測站數(shù)和距離來確定“權(quán)”值時,(3-11)分別可以轉(zhuǎn)化為: (3-12) (3-13) 上述結(jié)論也可應用于小三角水準網(wǎng)平差計算。 3.2.3 精度評定 單位

25、權(quán)中誤差： （3- 0 t v pv nt 12） 平差參數(shù)的協(xié)方差陣： （3- 21 0bb x x dn 13） 平差函數(shù)的協(xié)方差陣： （3- 21 0 t bb qf nf 14） 3.3 水準網(wǎng)間接平差程序信息設計 1數(shù)據(jù)文件的組織 下面給出一個水準網(wǎng)輸入數(shù)據(jù)文件的例子： 3 3 6（已知點個數(shù)、未知點個數(shù)、觀測值個數(shù)） 101 102 103 104 105 106 （點號） 34.788 35.259 37.825 （已知點高程） 104 101 1.625 4.5（起點點號、終點點號、高差觀測值、距離觀測值） 101 102 -0.418 3.1 105 102 0.714 3.

26、4 102 103 1.243 3.8 106 103 -0.577 4.3 103 101 -0.786 2.5 （其中編號數(shù)組未知點在前，已知點在后） 2水準網(wǎng)平差變量約定 表 3-1 變量約定表 變量名說明 ed 已知點個數(shù) dd 未知點個數(shù) sd 總點數(shù) gd 觀測值個數(shù) pn 點號 h0 已知點高程 be 起點點號 en 終點點號 h1 高差觀測值 s 距離觀測值 3.4水準網(wǎng)程序與使用說明 3.4.1 水準網(wǎng)程序流程圖 圖 3-2 水準網(wǎng)流程圖 程序的全部代碼見附錄一。 3.4.2水準網(wǎng)程序的使用 本程序使用 matlab 的矩陣功能計算法方程，在運行程序前首先要有其始數(shù)據(jù)。 其始

27、數(shù)據(jù)是一文件的形式保存在磁盤中，文件的格式在上文已經(jīng)說明過，編好文件 后，以后綴名為.txt 的形式保存。執(zhí)行時在 matlab 命令窗口直接鍵入文件名即 可。 3.5 案例 如下圖水準網(wǎng)，104、105、106 為已知點，101、102、103 為待定點，已知點的 高程分別為 34.788，35.259 ，37.825 。 觀測高差和觀測路線長度分別為： h1=1.652, h2=-0.418 ,h3=0.714 ,h4=1.243, h5=-0.577, h6=-0.786. s1=4.5, s2=3.1, s3=3.4, s4=3.8, s5=4.3, s6=2.5. 圖3-3 水準網(wǎng)圖

28、 首先編數(shù)據(jù)文件，命名為data1.txt. 數(shù)據(jù)的格式如下： 3 3 6 101 102 103 104 105 106 34.788 35.259 37.825 104 101 1.652 4.5 101 102 -0.418 3.1 105 102 0.714 3.4 102 103 1.243 3.8 106 103 -0.577 4.3 103 101 -0.786 2.5 進入matlab界面，在命令窗口直接輸入level3運行程序。彈出如下窗口 圖3-4 數(shù)據(jù)讀入文件 選擇data1.txt即可運行出如下結(jié)果： 圖3-5 計算結(jié)果 在圖 3-5 中，分別輸出了高程的平差值及精度。

29、結(jié)果是一文本的形式保存，用 戶可對它進行編輯。 4. 測邊網(wǎng)平差程序設計 4.1 數(shù)學模型 4.1.1 誤差方程和法方程的組成 控制網(wǎng)中的觀測值為邊長，誤差方程非零項最多為 4 個，所以誤差方程系數(shù)矩 陣采用壓縮格式進行儲存。 可采用以下的方法： a(m,n)=a(m,9) 其中，m 為觀測值個數(shù)，n 為未知點個數(shù)的兩倍。 改進后的 a 陣格式為 =（編號 1，系數(shù) 1，編號 2，系數(shù) 2，編號 4，系數(shù) 4，常數(shù)項） i a 共 9 列。即只存儲誤差方程的 4 個非零參數(shù)系數(shù)。 法方程系數(shù)陣為對稱陣，在存儲時，只需要存其上三角部分就可以了。其占 a n 用的空間為： (1) 2 n n su

30、m 現(xiàn)有 a 陣： a=（編號 1，系數(shù) 1，編號 2，系數(shù) 2，編號 4，系數(shù) 4，常數(shù)項） 其中偶數(shù)項為系數(shù)，加上最后的 a9 為常數(shù)項，在組成法方程時，從 a2 開始分別與 剩下的偶數(shù)項以及常數(shù)項相乘，然后再用 a4 與剩余的項相乘，一直到 a8 為止，這 樣就完成了的過程。需要注意的是：若 a1，a3，a5，a7 小于零，則表 t a na pa 示該點已知點，不參與法方程的組成。 4.1.2 邊長觀測的權(quán) 邊長觀測的精度一般與其長度有關(guān)，定權(quán)公式為 2 0 2 (1,2, ) i i s s pin 式中為所測邊長的方差，為任意選定的單位權(quán)方差。 2 i s i s 2 0 為了定權(quán)

31、必須已知測邊的先驗方差，但精確的已知是十分困難的，一般 i s p 2 i s 采用廠方給定的測距儀精度，即 i si abs 式中，a 為固定誤差（單位 mm），b 為比例誤差（單位：ppm），為邊長（單 i s 位 km）。 4.1.3 解算法方程 由于法方程是對稱正定陣，因此，可采用改進的平方根法進行解算。 平方根法是對稱正定矩陣非常有效的三角分解方法，設 a 為 n 階方陣，如果其 所有順序主子式均不為零，則其存在唯一的分解式： a=ldr 其中 l =，d=，r = 2 1,1 1 1 nn n l ll 1 0 0 n d d 121 1, 1 1 n nn rr r 由于此住 a

32、 對稱性，得，又根據(jù) a 陣正定的性質(zhì)，可證明 d 均為正數(shù)。 t lr 現(xiàn)在設 d= 1 n d d 1 1 d d d 1 1 d dd 即 11 22 dd d 則 1111 2222 ()() t ttt aldlld d lldd lll 為方便，記為 t all 稱為 cholesky 分解，即正定對稱矩陣的平方根分解法。解等階于求解兩個三axb 角方程組： 和lyb t l xy 在用平方根分解法計算時，需要進行 n 次開方運算。為了避免開方，可以直接 采用對稱正定的分解式對平方根法進行改進。從而解方程組可以按 t aldlaxb 如下步驟進行：把 a 分解成，則變成，即等價于

33、t aldlaxb() t ldlxb 1t lyb l xd y 由此可以解出 x 和 y。這稱為改進的平方根法，在計算中避免了開方運算。 平方根法和改進的平方根法的計算量和存儲量比消去法節(jié)約近一半，而且不需 要選主元，能得到比較精確的數(shù)值解。 法方程用改進平方根法解算的過程如下： (1)分解： 1t cs d s 其中 11111 , n mnmn sss sd ss 11111112 11 1 1, 1 111,1 1 1 nn nnnmn nn n nn ccsss s ccs ss ss 11 12 112 1 2222 1111 13 123 213 123 2 3333 1122

34、1122 ,2 ,3 jj jj jjjj jjjj jjjj sc s ss s csscj ss s ss ss ss s csscj ssss 純量計算公式為 1 1 1 1 1,2 ij i jkikj ij k kk ci ss s cij s （2）求逆 1 rs 111n nn rr r r 由 rs=i 得 1112111121 222222 100 010 001 nn nn nnnn rrrsss rrss rs 純量計算公式： 11 12 12 22 11 131223 13 33 11 11(1)(1) 1 22 12 23 33 222423 34 24 44 222

35、2(1)(1) 2 1 () ( () ( ii ii nnnn n nn nnnn n nn r s r s r s r sr s r s r srs r s r s r s r sr s r s r srs r s 通式為 1 1 1 , ii ii j ikkj k ij ij r s r s rji s （3）求積 1111 ()() ttt qs d ss d srdr 1111111 1 n t nnnnnnn rrsr qrdr rsrr 11 1112 121111 222222 1 nn nnnnnn nnnnnnn nn s rs rs rr rrs rs r rrrs r

36、 4.1.4 精度評定 （1） 坐標改正數(shù)以及單位權(quán)中誤差的計算 0 m 使用上三角一維數(shù)組形式存儲坐標改正數(shù)的公式為： 1 11 ,1, in ijijijj jj xq wq win 其中，n=2 dd,的單位是 cm. ix 平差值： 0 xxx 寫成分量的形式，為 0 iii xxx 如果近似坐標的誤差較大，或網(wǎng)形較大，平差的結(jié)果不會精確，這時，就需要進行 迭 代平差，直到兩次平差間互差在允許值內(nèi)。 由測量平差理論： t v pv nt 同樣可得到單位權(quán)中誤差： 0 pvv m mn 其中，m-n 觀測個數(shù)減去未知點個數(shù)； 123, 2, mmmm nddst st pvvpllx 方

37、向觀測的測站數(shù) （2）點位誤差橢圓 誤差橢圓表示了網(wǎng)中點或點與點之間的誤差分布情況如圖。在測量工作中，常 用的誤差橢圓對布網(wǎng)方案作精度分析。繪制誤差橢圓只需要三個數(shù)據(jù)：橢圓長半軸 a，短半軸 b 和主軸方向 ,其求法為 2222222 22 2222222 1 ()4) 2 2 ,tan2 1 ()4) 2 xyxyxy xy xy xyxyxy a b 圖 4-1 誤差橢圓的表達 顧及方差與權(quán)倒數(shù)的關(guān)系，得 2 222 0 2 222 0 ()4) 2 2 ,tan2 ()4) 2 xxyyxxyyxy xy xxyy xxyyxxyyxy aqqqqq q qq bqqqqq 根據(jù)上述的理

38、論，我們實際要求的是、。只要得到了這些元素， xi m yi m xiyi m 就能依照上面的公式來求得橢圓的元素了。 4.2測邊網(wǎng)平差信息設計 外業(yè)測量的數(shù)據(jù)首先應進行預處理，包括測站平差、歸心計算、觀測值歸化到 橢球面的改正、橢球面歸化到高斯平面的改正等，然后將預處理后的數(shù)據(jù)輸入到以 后綴名為 txt 的文本文件中，該數(shù)據(jù)文件的組織格式如下所列： 表 4-1 數(shù)據(jù)組織格式 次序內(nèi)容 1已知點個數(shù) ed,未知點個數(shù) dd,控制參數(shù) 2點號 pn。先輸入已知點編號，各點輸入順序無要求 3已知點坐標,x0,y0,x1,y1,x2,y2 4測量邊的個數(shù) m1 5 邊長的固定誤差 ms,邊長的比例誤

39、差 pp。單位分別是：cm 和 0.000001 6邊長的起始點號 e,終點點號 d,邊長 sid。每一條邊一行，依次列出 7推算近似坐標的路線經(jīng)過的邊數(shù) 8推算近似坐標的起算已知點坐標（按順時針） 9推算近似坐標的路線經(jīng)過的邊的邊號 4.2.1 主要的技術(shù)要求 表 4-2 測邊網(wǎng)的技術(shù)要求 等級平均邊長（km） 測距中誤差 （mm） 測距相對中誤差 二等9 30 130 萬 三等5 30116 萬 四等2 16 112 萬 一級小三角1 1616 萬 二級小三角0.5 16 13 萬 43 利用 matlab 的繪圖語句繪制網(wǎng)圖 1. 網(wǎng)形的繪制 由于網(wǎng)形圖與誤差橢圓繪制在同一個圖形上，因此

40、必須對誤差橢圓進行放大， 在本文的程序中，使用了 inputdlg 對話框輸入語句，其中，確省的放大倍數(shù)為 100。 在程序中，使用了 ed dd pn m1 x y e d sid ai bi fi 等變量，其意義與前面的變量相同。 對繪制的網(wǎng)圖有效放大和縮小功能，即點擊鼠標左鍵放大圖形，點擊右鍵縮小圖形， 利用 matlab 菜單本身的功能，可以將該圖形輸出為各種圖形文件格式。 2 誤差橢圓的繪制 無論多么復雜的圖形，其基本單元還是點和線。換句話說，只需要利用基本元 素的點或線，通過各種組合，也能畫出復雜的圖形。matlab 中沒有提供直接繪制 橢圓的命令，因此可以直接利用連線來畫橢圓。

41、測量中描述誤差橢圓用長半軸 a，短半軸 b 和方位角 fi 三個量。在如圖 4-1 的 xoy中直角坐標系中，橢圓的標準方程為 （4- 22 1 xy ab 1） 如果以角度 i 為變量，則橢圓的標準參數(shù)方程為 （4- cosi sini x y a 2） 設在測量坐標系 xoy 中橢圓長半軸的方位角為，則有 0 （4- 00 00 cossin sincos xx yy 3） 用參數(shù)方程代入，得到 （4- 00 00 coscossinsin sincoscossin xaibi yaibi 4） 測量坐標系與 matlab 或 autocad 繪制的數(shù)學坐標系的 x,y 坐標軸不同， 繪圖

42、時，需要調(diào)換 x，y 坐標。 在上式中，i 取不同的值，就有一組（x，y） ，只需要將這些點連接起來，就可以繪 出一個橢圓。 圖 4-2 誤差橢圓的參數(shù)方程 4.44.4 測邊網(wǎng)程序和使用說明測邊網(wǎng)程序和使用說明 使用本程序的全部數(shù)據(jù)都按規(guī)定的格式編輯成數(shù)據(jù)文件儲存在磁盤上。數(shù)據(jù)文 件的編輯取決于平差網(wǎng)型和觀測值的編號，為此，先繪制平差網(wǎng)的略圖，在圖上標 明各項數(shù)據(jù)的信息。以下就測邊網(wǎng)的點號，邊編號，輸入數(shù)據(jù)，輸出成果，運行程 序等問題作簡明說明。 1 點號和觀測邊的編號 已知點和待定點的編號為 pn（取三為數(shù)），已知點在前，未知點在后，其順序 無要求。但為了減小法方程系數(shù)的帶寬，應使相臨的待

43、定點編號的差數(shù)盡可能小。 平差網(wǎng)的編號見圖 4-3。 2.推算近似坐標的路線 近似坐標的路線是用戶在測邊網(wǎng)略圖上指定出來的。如圖 4-3 的箭頭，就是表 示推算路線。路線的兩個起算點必須為已知點，從兩個已知推算出的第一個未知點 開始，選擇觀測邊，由觀測邊和已求得的近似坐標或已知坐標推算出觀測邊所對的 未知點。本程序是按推算的三個點 a、b、p 順序為順時針。 圖 4-3 測邊網(wǎng)略圖 2 數(shù)據(jù)的輸入 （1）簡單變量 為了在程序運行中數(shù)據(jù)的傳遞，定義了一些全局變量,參見 4-1 表。 （2）數(shù)據(jù)文件 外業(yè)測量的數(shù)據(jù)首先應進行預處理，包括測站平差、歸心計算、觀測值歸化到 橢球面的改正、橢球面歸化到高

44、斯平面的改正等，然后將預處理后的數(shù)據(jù)輸入到以 后綴名為 txt 的文本文件中，該數(shù)據(jù)文件的組織格式如表 4-1。 3 輸出成果 本程序的計算成果是以文件的形式輸出到文本文件中，用戶可以在文本中查看 和編輯。 4.5 程序代碼說明： 程序的總體框架 圖 4-4 測邊網(wǎng)總體流程圖 1數(shù)據(jù)讀入塊 本模塊的功能是打開一個*.txt 的數(shù)據(jù)文件，同時生成一個*out.txt 的文本文件， 記錄用戶數(shù)據(jù)和輸出成果用。程序調(diào)用 fscanf 函數(shù)把文件中的數(shù)據(jù)賦值給相應的變 量，這些變量是后面計算的數(shù)據(jù)依據(jù)。 2近似坐標的計算 測邊網(wǎng)的觀測數(shù)據(jù)是邊長，所以在近似坐標的計算時只能用測邊交會計算。以 及 a、b

45、 兩點坐標及 a、b 到 p 點的距離 b、a，c 為 a、b 兩點的距離，a、b、p 三點按順時針排列。則 p 點的坐標計算公式如下： 222 22 , 2 cossin sincos paabab paabab bca efbe c xxef yxef 計算近似坐標的流程圖入圖 4-5。 圖 4-5 計算近似坐標流程圖 3誤差方程和法方程的形成 程序中，用數(shù)組 a 來存儲誤差方程的編號和系數(shù)，a(i,9)存儲常數(shù)項，w 和 c 分 別存儲法方程系數(shù)和法方程常數(shù)。 4解算法方程。 函數(shù)流程圖如圖 4-6。函數(shù)的代碼見附錄三。 5精度評定 本模塊包括坐標改正數(shù)、單位權(quán)中誤差和誤差橢圓的計算。程

46、序中定義了 dxy（坐標改正數(shù)） 、pvv（即存儲） 、uw0(單位權(quán)中誤差)等。 pvvpllx 同時計算出誤差橢圓的三個參數(shù)長半軸 ai，短半軸 bi 和主軸方向 fi。 圖 4-6 平方根法求逆程序框架圖 6控制網(wǎng)的繪制 程序中調(diào)用 inputdlg 函數(shù)來打開一個對話框來輸入誤差橢圓的參數(shù)；用 text 函 數(shù)來對點號標注；用 ploth 函數(shù)及控制參數(shù)繪制線和誤差橢圓。 程序的完全代碼見附錄二。 4.6 程序的使用算例 有測邊網(wǎng)如圖 4-3 所示。網(wǎng)中 a、b、c、c 點已知，其余為未知點，現(xiàn)用某測距 儀觀測了 13 條邊長，測距精度。起算數(shù)據(jù)和觀測數(shù)據(jù)如下： 6 (31 10) s

47、 mms 表 4-3 已知點坐標 點名xy點名xy a 53743.13661003.826c40049.22953782.79 b 47943.00266225.854d36924.72861027.086 表 4-4 觀測值 編號邊觀測值編號邊觀測值編號邊觀測值 15760.70668720.162115487.073 25187.34275598.57128884.587 37838.8887494.881137228.367 45483.15897493.323 55731.788105438.382 1.編輯數(shù)據(jù)文本文件如下圖 圖 4-7 數(shù)據(jù)文件 2.在 matlab 命令窗口鍵入

48、 nnbb 執(zhí)行程序，運行中會彈出一個對話匡提示用戶 輸入誤差橢圓的放大比例，默認為 100，本例選擇 500。如下圖 圖圖 4-8a4-8a 對話框?qū)υ捒?圖圖 4-8b4-8b 對話框?qū)υ捒?3.計算出的結(jié)果如下： 在圖 4-9a 中，第 1、2 行分別是已知點和未知的個數(shù)；第 4 行是點的編號；第 5 至第 8 行是已知點的坐標；第 10 行是觀測值的個數(shù)；第 12 行是測距的固定誤差 和比例誤差；13 至 25 行是觀測邊的起點號、終點號和觀測邊長；27 至 39 行（接到 圖 4-9b）是點號轉(zhuǎn)換為計算機順序后的觀測邊起點號、終點號和觀測邊長。 接圖 4-9a，在圖 4-9b 中，第

49、 42 至 49 行是推算的近似坐標；第 52 至 64 行是 計算的誤差方程系數(shù)和常數(shù)。 在圖 4-9 中，第 67 至 70 行是法方程的系數(shù)（上三角一維存儲）；第 73 至 76 行是求逆后的法方程系數(shù)；第 79 至 86 是坐標改正數(shù)和坐標平差值；第 89 至 92 行 是誤差橢圓的參數(shù)；第 95 行是單位權(quán)中誤差。 圖 4-9a 計算結(jié)果 圖 4-9b 計算結(jié)果 圖 4-9c 計算結(jié)果 4.輸出的控制網(wǎng)圖和誤差橢圓圖如下: 圖 4-10 控制網(wǎng)圖 結(jié)論 由本論文的兩個程序例子解算過程不難發(fā)現(xiàn)，程序的平原理和其它高級語言編 寫的程序大致相同的。但matlab 語言本身就有獨特的數(shù)值計算

50、功能和圖形繪制功 能，這使它不用編寫專門的函數(shù)而直接處理和計算測量的數(shù)據(jù)，這就大大縮短了編 程的時間。matlab 軟件將人們從乏味的 fortran、c 編程中解放出來，使他們真正的把精力 放在科學研究的核心問題上。但用 matlab 語言編程本身也有它的不足之處： 1.由于 matlab 是一個解釋器，會逐行對程序代碼進行解釋后執(zhí)行，當要處理的 數(shù)據(jù)量很大時，計算機的運行速度明顯變慢了好多。 2.用 matlab 進行界面開發(fā)時沒有像其他面向?qū)ο笳Z言方便。 針對上面的問題，提出一些見解： 1.為了提高整體程序的執(zhí)行效率，應盡量多使用向量化的運算，而避免或少用 for 循環(huán)、while 循環(huán)

51、。事實上，如果程序是純粹的數(shù)值運算、而沒有使用 for 循環(huán) 或 while 循環(huán)，那么其執(zhí)行速度將會接近于純粹用 c 語言來寫的程序代碼。 2.雖然 matlab 是一個完整的科學計算與可視化環(huán)境，但在多數(shù)情況下，如 matlab 編程的代碼執(zhí)行效率不佳、你希望隱藏你的程序代碼以保護產(chǎn)權(quán)、或您想要 進行任何只有在某些高級語言程序才能做到的事，這時就需要用 matlab 所提供的應 用程序接口。 comment 楊楊1: 定 義 讀 取 水 準 網(wǎng) 數(shù) 據(jù) 輸 出 形 參： ed d d s d gd p n h0 k1 k2 h 1 s comment 楊楊2: 定 義 全 局 變 量 co

52、mment 楊楊3: 定 義 計 算 水 準 網(wǎng) 函 數(shù) comment 楊楊4: 將矩陣 a 轉(zhuǎn)化為稀 疏矩陣形式，即矩陣 a 中任何 0 元素 被去除，非零元素及其下標組成 s comment 楊楊5: 在 a 矩陣中豎著數(shù) 第 n 個是-1 或 1 comment 楊楊6: 表示取 a 第 1 到 dd 列的元素 comment 楊楊7: 行列為 0 comment 楊楊8: x=diag(v,k) 以向量 v 的元素作為矩陣 x 的第 k 條對角線元 素 comment 楊楊9: x 的改正數(shù) comment 楊楊10: a(1:2)取 a 的從 1 到 2 的值 附錄一 水準網(wǎng)平差程

53、序代碼： function level3,ed,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata; global pathname net_name s_datafile a1_datafile; global ed dd sd pn gd h0 k1 k2 h1 s dh; dh,h,v,l,uw0,uwh,uw1=calculatelevelnet(ed,dd,sd,pn,gd,h0,k1,k2,h1,s); writelevelnetdata(pn,k1,k2,h1,v,l,h0,dh,h,uwh,uw0); %輸出水準網(wǎng)計算結(jié)果 return fun

54、ctiondh,h,v,l,uw0,uwh,uw1=calculatelevelnet(ed,dd,sd,pn,h0,be,en,hd, distance); % 水準平差網(wǎng) a=sparse(zeros(sd,gd); %求解系數(shù)陣 b=(0:(gd-1)*sd); a(be+b)=-1; a(en+b)=1; a=a; a=a(:,1:dd); l=zeros(gd,1); %求解常數(shù)項 l=h0(be)-h0(en)+hd; p=diag(1./distance); %權(quán)陣 dh=inv(a*p*a)*a*p*1; %高程改正數(shù) h00=h0(dd+1:sd); h0=h0(1:dd);

55、 %待定點高程近似值 h=h0+dh; %待定點高程平差值 v=a*dh-1; %高差觀測值改正數(shù) l=hd+v; %高差觀測值平差值 %精度評定 comment 楊楊11: 單位權(quán)中誤差 comment 楊楊12: comment 楊楊13: 定義空矩陣 comment 楊楊14: 空集 comment 楊楊15: or 或者 comment 楊楊16: %獲取所選文件信 息（文件名、路徑等） comment 楊楊17: b=find(a),a 是一 個矩陣,查詢非零元素的位置 b=find(a,2),找出 a 中最先出現(xiàn)的 2 個不為零的位置 comment 楊楊18: fileid =

56、fopen(filename,permission) fprintf 函數(shù)可以將數(shù)據(jù)按指定格 式寫入到指定的文本文件中基本 格式 fprintf(fid, format, variables) comment 楊楊19: 表示文件打開不成 功 comment 楊楊20: a = fscanf(fileid, format, sizea) uw0=sqrt(v*p*v/(gd-dd); %單位權(quán)中誤差 qxx=inv(a*p*a); uwh=uw0*sqrt(diag(qxx); %待定點高程平差值中誤差 uwh(dd+1:ed+dd)=0.0; qff=a*qxx*a; uwh=uw0*sqr

57、t(diag(qff); %高差平差值中誤差 h=h;h00; %所有點高程 h0=h0;h00; dh=dh;zeros(ed,1); return function ed,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata; global pathname net_name s_datafile b_datafile; global ed dd sd pn gd h0 k1 k2 h1 s k11 k12; k1=;k2=;h=;s=; if(isempty(pathname)|isempty(net_name) filename,pathname=uig

58、etfile(*.txt,input filename); i=find(.=filename); net_name=filename(1:i-1); end fid1=fopen(strcat(pathname,net_name,s_datafile),rt); if(fid1=-1) msgbox(input file or path is not correct,warning,warn); return; end %open afile to read %open afile to read ed=fscanf(fid1,%f,1); %已知點個數(shù) dd=fscanf(fid1,%f,

59、1); %未知點個數(shù) sd=ed+dd; %總點數(shù) gd=fscanf(fid1,%f,1); %觀測點個數(shù) pn=fscanf(fid1,%f,sd); %點號 %known data h0=fscanf(fid1,%f,ed); %已知點高程 h0(dd+1:ed+dd)=h0(1:ed) comment 楊楊21: comment 楊楊22: 重組 heightdiff=fscanf(fid1,%f,4,gd); heightdiff=heightdiff; k1=heightdiff(:,1); %起點 k2=heightdiff(:,2); %終點 k11=heightdiff(:,

60、1); %起點 k12=heightdiff(:,2); %終點 h1=heightdiff(:,3); %高差 s=heightdiff(:,4); %距離 fclose(all); %點號轉(zhuǎn)換 k1,k01=chkdat(sd,pn,k1); k2,k02=chkdat(sd,pn,k2); h0(1:dd)=20000.; ie=0; while(1)%計算近似高程 for k=1:gd i=k1(k); j=k2(k); if(h0(i)le4) h0(j)=h0(i)+h1(k); ie=ie+1; end if(h0(i)le4 a=fscanf(fit1,%f,2); ms=a(