基于MATLAB的測(cè)量平差數(shù)據(jù)處理講解

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

2、率。關(guān)鍵詞：MATLAB；平差；程序設(shè)計(jì)ABSTRACTMATLAB is one species of numeiical-values calculation and grapluc tools software wluch is widely used to apply at reseaich institutions at present. The paiticulaiities are: concise gnunmai-structure、lughly efficient in numeiical values calculating、complete finiction of gi

3、aplis especially it is adapted to evildoing professional piogranunei- to accomplish the tasks that are nunieiical-values calculating and scientific expenments treating. The ancient methods of measured data-processmg need establislung special pioceedmgs of treating matiices operation, moieovei; it is

4、 complex and gieatlv difficultThis article mtroduces one piogianuning method dealing with leveling and measunng edge network based on MATLAB. Compared with other algoiitlmi language, it has paiticulaiities which are simply progianuiinig and quickly operating. The article separately expatiate the the

5、ories basics fealizuig steps and mimmg results at leveling and measuring edge netwoik. With the analysis of examples, it has prodigious application value in measured data-processing by use of MATLAB Moreover, it shortens piogianuning time and unpioves working effectivenessKey words： MATLAB； progianm

6、img緒論作為一名測(cè)量技術(shù)人員，如果不掌握一門PC機(jī)編程語言與便攜計(jì)算工具，要想 提高測(cè)量工作的效率幾乎寸步難行。測(cè)量需求的多樣性與復(fù)雜性，造就了測(cè)量計(jì)算鮮 明的個(gè)性化特點(diǎn)，這就是在商業(yè)測(cè)量計(jì)算軟件高度發(fā)達(dá)的今天，掌握一種實(shí)用的程序 語言進(jìn)行編程計(jì)算仍有廣泛的市場(chǎng)需求的重要原因。當(dāng)今較流行的計(jì)算機(jī)程序語言基本上都是基于Wmdows的，例如Turbo Pascal,Visual Basic,Visual C,Borland C十十等，這些程序語言的優(yōu)勢(shì)是基于對(duì)象及可利用 Wmdows豐富的系統(tǒng)資源，應(yīng)用它們可以開發(fā)出界面非常豐富和友好的應(yīng)用程序，其 劣勢(shì)主要有以下幾點(diǎn)：1 .Windows程序都

7、非常龐大，學(xué)習(xí)并熟練掌握它們并非易事。2. 雖然市場(chǎng)上己有的多種專用的測(cè)量平差軟件都是采用C語言開發(fā)的，但這些軟 件價(jià)格都比較貴，而且都帶有加密狗，一次只能供一個(gè)用戶使用。出于商業(yè)目的，開 發(fā)商不會(huì)公開程序源代碼，這為修改程序功能以適應(yīng)用戶的特殊需求帶來了不便。3. 在測(cè)量生產(chǎn)中，經(jīng)常需要根據(jù)工程的實(shí)際情況進(jìn)行一些個(gè)性化的數(shù)值計(jì)算工作, 這些數(shù)值計(jì)算工作無固定模式，這就需要求測(cè)量技術(shù)人員最好能熟練掌握一種適用于 數(shù)值計(jì)算的程序語言，以便提高測(cè)量計(jì)算的效率。4. C語言的數(shù)值計(jì)算語句不夠豐富，例如，在測(cè)量平差計(jì)算中，經(jīng)常需要進(jìn)行的 矩陣運(yùn)算，尤其是解法方程的矩陣求逆不能直接使用語句實(shí)現(xiàn)，而必須應(yīng)

8、用計(jì)算機(jī)算 法編程實(shí)現(xiàn)。如果不是基于商業(yè)軟件開發(fā)，只為滿足實(shí)際測(cè)量工作計(jì)算需要，則C語言的劣勢(shì) 就變成了 MATLAB語言的優(yōu)勢(shì)。內(nèi)蒙古科技大學(xué)畢業(yè)設(shè)計(jì)MATLAB軟件簡(jiǎn)介1.MATLAB軟件簡(jiǎn)介MATLAB是從Matrix （矩陣）和Laboratory （實(shí)驗(yàn)室）各取前3個(gè)字母組成的，意思是 矩陣實(shí)驗(yàn)室,是美國(guó)MathWorks公司于20世紀(jì)80年代中期推出的一種交互式、面向?qū)ο?的科技應(yīng)用軟件，是一個(gè)為科學(xué)和工程計(jì)算而專門設(shè)計(jì)的高級(jí)交互式軟件包。MATLAB 集成了圖示與精確的數(shù)值計(jì)算,是一個(gè)可以完成各種計(jì)算和數(shù)據(jù)可視化的強(qiáng)有力工具, 其優(yōu)秀的數(shù)值計(jì)算能力和卓越的數(shù)據(jù)可視化能力使其很快在

9、數(shù)學(xué)軟件中脫穎而出，成 為以矩陣運(yùn)算為主要工作方式的線性代數(shù)、概率論和數(shù)理統(tǒng)計(jì)、自動(dòng)控制、數(shù)字信號(hào) 處理、動(dòng)態(tài)系統(tǒng)仿真等領(lǐng)域教學(xué)和科研工作者的有力武器。隨著該軟件自身的發(fā)展及 市場(chǎng)的需求，其功能日趨完善，其最高版本7. 0版己經(jīng)推出，隨著版本的不斷升級(jí)，它的 數(shù)值計(jì)算及符號(hào)計(jì)算功能得到了進(jìn)一步完善。MATLAB是以矩陣作為數(shù)據(jù)操作的基本單位，矩陣的生成、運(yùn)算、轉(zhuǎn)置、求逆等非 常簡(jiǎn)單。在MATLAB環(huán)境中，不需要對(duì)創(chuàng)建的變量對(duì)象給出類型說明和維數(shù)，所有的變量 都作為雙精度數(shù)來分配內(nèi)存空間,MATLAB將自動(dòng)地為每一個(gè)變量分配內(nèi)存。MATLAB語 言起源于矩陣運(yùn)算，并己經(jīng)發(fā)展成為一種高度集成的計(jì)算

10、機(jī)語言，它提供了強(qiáng)大的科 學(xué)運(yùn)算、靈活的程序設(shè)計(jì)流程、高質(zhì)量的圖形可視化與界面設(shè)計(jì)、便捷的與其他程序 和語言接口的功能。MATLAB系統(tǒng)主要包含5部分的內(nèi)容:MATLAB工作環(huán)境31ablab數(shù) 學(xué)函數(shù)庫、MATLAB語言體系、句柄圖形、MATLAB應(yīng)用程序接口 （API） MATLAB系統(tǒng)主 要功能包括：數(shù)值計(jì)算功能、符號(hào)計(jì)算功能、數(shù)據(jù)分析和可視化、文字處理功能、 SIMULINK動(dòng)態(tài)仿真功能。同時(shí),MATLAB乂是開放的，除了內(nèi)部函數(shù)之外，所有的MATLAB 主包文件和各工具包文件都是可讀可改的源文件，用戶可以作為參考掌握其用法，并 可對(duì)其修改以適應(yīng)自己的需要,也可加入自己編寫的文件構(gòu)成新

11、的工具包。例如，隨著 GPS 的廣泛應(yīng)用，Orion Dynamics and Con21 rol Corporation、Constell Inc. GPSSoft LLC、NavsysCorporation等多家公司都相應(yīng)開發(fā)出了適于GPS數(shù)據(jù)處理的 MATLAB工具箱。MATLAB是一個(gè)集數(shù)值計(jì)算、圖形管理、程序開發(fā)于一體的功能十分強(qiáng)大的系統(tǒng)。 將MATLAB應(yīng)用于測(cè)量數(shù)據(jù)的處理是一件非常有意義的工作。Mo2hamed等曾成功地在 MATLAB系統(tǒng)中利用白濾波技術(shù)研究動(dòng)態(tài)解算GPS載波相位信號(hào)的模糊度問題。因?yàn)?測(cè)量數(shù)據(jù)的處理特別是測(cè)量平差主要應(yīng)用矩陣運(yùn)算，而MATLAB 乂特別易于做矩

12、陣運(yùn) 算，因此，研究開發(fā)基于MATLAB的測(cè)量平差方法具有極好的應(yīng)用價(jià)值。畢業(yè)設(shè)計(jì)水準(zhǔn)網(wǎng)平差程序2. MATLAB在測(cè)量平差中的應(yīng)用測(cè)量平差數(shù)據(jù)處理主要是基于矩陣的運(yùn)算，常用的矩陣運(yùn)算主要是矩陣的生成、 轉(zhuǎn)置、求逆和矩陣求廣義逆等。在MATLAB環(huán)境中，不需要對(duì)創(chuàng)建的變量對(duì)象給出類型 說明和維數(shù)，所有的變量都作為MATLAB中的M文件的語法與其他的高級(jí)語言類似，是一 種程序化的編程語言，同時(shí)也是一種解釋性的編程語言，即逐行解釋運(yùn)行程序,使程序 容易調(diào)試,計(jì)算更為簡(jiǎn)捷，而且對(duì)于平差原理理解和掌握變得更容易。另外,MATLAB語 言與數(shù)學(xué)語言比較接近,更容易掌握和理解。2. 1測(cè)量平差原理的概述測(cè)

13、量平差的函數(shù)模型有條件方程和觀測(cè)方程。以條件方程為函數(shù)的模型的最小二 乘平差稱為條件平差；在條件方程中，根據(jù)需要如果還設(shè)有一定數(shù)量的未知數(shù)，則稱 為附有參數(shù)的條件平差；以觀測(cè)方程為函數(shù)模型的最小二乘平差稱為間接平差；如果 觀測(cè)方程中的某些參數(shù)不獨(dú)立，則這些不獨(dú)立參數(shù)必然存在一些條件，稱這種平差模 型為附有條件的間接平差。本文的兩個(gè)程序都采用間接平差模型。對(duì)于一個(gè)實(shí)際平差問題，根據(jù)所選參數(shù)的個(gè)數(shù)、選什么量為參數(shù)以及參數(shù)之間是 否函數(shù)獨(dú)立，經(jīng)過仔細(xì)推敲可以發(fā)現(xiàn)附有條件的間接平差模型本身就是各種經(jīng)典平差 模型的概括模型,其余的經(jīng)典平差模型，如條件平差模型、間接平差模型、附有未知數(shù) 的條件平差模型和附

14、有限制條件的條件平差模型都是它的特例。間接平差的公式匯集：間接平差模型為(2-1)(2-2)(2-3)V = BX-l 1 VTPV f min系數(shù)矩陣B滿秩，即rank(B)= t法方程及解為：gx-fc = Wb = BTPBJe = BrPl)參數(shù)的平差值：X = X + x(2-4)觀測(cè)量的平差值：L=L+V(2-5)單位權(quán)的中誤差：b。=(2-6)(2-7)(2-8)平差參數(shù)的協(xié)方差陣：Dxx=crQN平差函數(shù)的協(xié)方差陣：0廠元嚴(yán)滬2. 2平差程序總體方案MATLAB號(hào)稱為全球工程師的共同語言，其語法和C語言相似，但它有強(qiáng)大的數(shù)值 計(jì)算和繪圖功能，這使之在工程應(yīng)用方面的計(jì)算更出色，本

15、文就基于這種程序設(shè)計(jì)語 言環(huán)境設(shè)計(jì)一個(gè)控制網(wǎng)平差程序。該程序包含了一個(gè)高程控制網(wǎng)平差程序和測(cè)邊網(wǎng)平 差程序。本程序適用于各種等級(jí)的高程網(wǎng)和測(cè)邊網(wǎng)，程序在設(shè)計(jì)過程中，始終考慮數(shù)據(jù)的 儲(chǔ)存量。因而本程序不儲(chǔ)存誤差方程的系數(shù)和常數(shù)項(xiàng)，對(duì)待定點(diǎn)數(shù)較多的平差網(wǎng)，組 成法方程的系數(shù)矩陣是個(gè)稀疏矩陣，如待定點(diǎn)的編號(hào)恰當(dāng)，法方程的系數(shù)會(huì)集中在主 元系數(shù)的兩側(cè)形成帶狀。為減少法方程系數(shù)的儲(chǔ)存量，只要按行儲(chǔ)存下三角陣或按列 儲(chǔ)存上三角陣中第一個(gè)非零系數(shù)起的系數(shù)，就是通常叫做維變帶寬儲(chǔ)存方法。3. 水準(zhǔn)網(wǎng)平差程序 3. 1程序的功能本程序適用于二、三、四等水準(zhǔn)網(wǎng)平差計(jì)算，平差的水準(zhǔn)網(wǎng)可以是獨(dú)立的、也可 以是附合網(wǎng)，

16、其主要功能是完成水準(zhǔn)網(wǎng)的平差計(jì)算和精度評(píng)定計(jì)算。平差計(jì)算采用間 接平差法，以歸心的觀測(cè)值為高差，以未知點(diǎn)高程為未知參數(shù)。精度評(píng)定計(jì)算包括計(jì) 算單位權(quán)中誤差和每個(gè)待定點(diǎn)的高程中誤差。3. 2水準(zhǔn)模型網(wǎng)的間接平差3. 2. 1“權(quán)”值的確定當(dāng)在相同的條件下進(jìn)行水準(zhǔn)測(cè)量時(shí)，其精度是相同的，因而觀測(cè)結(jié)果的可靠性也 是同樣的。但如果在不同的條件下進(jìn)行水準(zhǔn)測(cè)量時(shí)，高程的精度就有所不同，此時(shí)稱為 不等精度觀測(cè)，所求出的未知量的值、高程的最或是值并對(duì)其精度進(jìn)行評(píng)定時(shí)，就需要 “權(quán)” 了。由于觀測(cè)的不等精度，因而觀測(cè)值的可靠程度不同，求未知量的最或是值時(shí)，這樣 的一個(gè)因素就必須考慮了，這個(gè)因素是:可靠性大的某觀

17、測(cè)值，其精度高，對(duì)測(cè)量的最 后結(jié)果的影響也就越大。此時(shí)用“權(quán)”值來表示觀測(cè)值的可靠程度，那么，“權(quán)”值愈 大，觀測(cè)值的可靠程度就愈高。另外，在觀測(cè)過程中，觀測(cè)值的中誤差愈小，觀測(cè)結(jié)果愈 可靠，它的“權(quán)”值就愈大。因而，根據(jù)中誤差來確定“權(quán)”值是非常適當(dāng)?shù)?。設(shè)以刃 表示觀測(cè)值的“權(quán)” 5為中誤差，則“權(quán)”值的定義為：叫(3-1)式中必為任意的正常數(shù)，在一組觀測(cè)值中為一個(gè)定數(shù)。在實(shí)際測(cè)量中，通常是觀測(cè)值的中誤差事先并不知道，因而必須先確定觀測(cè)值的“權(quán)”, 然后才能求出未知量的最或是值。此時(shí)可以利用距離(5)或測(cè)站數(shù)(A)來確定觀測(cè)值一 一一高程的“權(quán)”。根據(jù)偶然誤差傳播定律，各觀測(cè)點(diǎn)高程旳的中誤差

18、加由測(cè)站數(shù)M確定時(shí),則有：叫=n?(3-2)式中:加為一組觀測(cè)值的中誤差,為一個(gè)定數(shù).由(3-2-1) (3-2-2)兩式可得：P,=亠=丄1礦 MM(3-3)式中：C=氣r由于A、加都為定數(shù)、故C也為定數(shù). m同樣可得出:=f(3-4)式中：勸定數(shù)，耳為測(cè)距.由(3-3) (3-4)兩式可以得出這樣一個(gè)結(jié)論：當(dāng)測(cè)站觀測(cè)高差等精度時(shí)，觀測(cè)總高差 的“權(quán)”與測(cè)站數(shù)或距離成反比。3. 2. 2水準(zhǔn)路線的平差計(jì)算1. 附合路線的平差計(jì)算假定在圖1示的力、萬兩水準(zhǔn)點(diǎn)之間布設(shè)一條水準(zhǔn)路線，/、方兩水準(zhǔn)點(diǎn)的高程 為已知，分別設(shè)為日人、hb、“、C為中間水準(zhǔn)點(diǎn)。假定觀測(cè)了所有的點(diǎn)的高 程,現(xiàn)擬求Q點(diǎn)的高程H

19、e的最或是值?？捎伤疁?zhǔn)路線A -a B _C分別觀測(cè)的高差hAC 4九C計(jì)算得出，由此而 得到的觀測(cè)高程分別設(shè)為弘1、c2,其值為：Hc-Ha + A hAC ； Hc2=孔 + 4 hBC當(dāng)尿1、*2在不等精度條件下觀測(cè)得出時(shí),它們的“權(quán)”也不同，分別設(shè)為斤1、 P&,這樣C點(diǎn)的高程日。的最或是值為：(3-5)根據(jù)力點(diǎn)的高程Ha,A 一 CK準(zhǔn)路線觀測(cè)的高差力g以及方一 CK準(zhǔn)路線觀測(cè)的高差/ /張，可推算出乃點(diǎn)的觀測(cè)高程/為：H，B = HA + AhAC-hBC水準(zhǔn)路線力-B的高程閉合差為：(3-6)fh = HfB-HB=Hcl-Hc2由(3-6)式得到：Hc2 = Hcl-fhcc由

20、3)式得到心喬石“、心分別表示水準(zhǔn)路線I B-C的測(cè)站數(shù)，水準(zhǔn)路線力f別勺測(cè)站數(shù)Nab = Nac + NQ將上述表達(dá)式代入(3-2-5)式中，得到：耳1 PqP + Pd辦(3-7)如果以水準(zhǔn)路線力一來J距離片。、方一期距離Spe/f尿勺距離S朋(SAB = SAC + SBC) 來確定高程觀測(cè)值的“權(quán)”值時(shí)，同樣可以得到：(3-8)B圖3T水準(zhǔn)路線圖2. 閉合路線的平差計(jì)算閉合路線的平差計(jì)算原理與附合路線相同，因而(3-7). (3-8)兩式的結(jié)論適用于 閉合路線的平差計(jì)算。(3)具有一個(gè)結(jié)點(diǎn)的水準(zhǔn)網(wǎng)的平差計(jì)算如圖2所示為具有一個(gè)結(jié)點(diǎn)的水準(zhǔn)網(wǎng),B, C, D,為己知高程水準(zhǔn)點(diǎn),B-A, C

21、-Ay D 一力,為水準(zhǔn)路線，則接點(diǎn)/的高程最或是值為：(3-9)_ EuH知 + Pa：Ha2 + PaHA3 + _f=l式中粘皿皿分別為水準(zhǔn)路線萬一4 C-AyD-Ay 計(jì)算力的觀測(cè)高程，各高程相應(yīng)的“權(quán)”值為PAl,PA2,PAi-設(shè)比肆心的算術(shù)平均值為日各高程觀測(cè)值與日鴛的差值分別為 劉1,丹2, M3,則有：1 2 3 + + + &4垃 一一=1 2 3 或- - 12 3=將(3-10)式代入(3-9)式得到:(3-11)當(dāng)以測(cè)站數(shù)和距離來確定“權(quán)”值時(shí)，（3-11）分別可以轉(zhuǎn)化為:r=l(3-12)f 召（成）憶=氏+(3-13)1=1上述結(jié)論也可應(yīng)用于小三角水準(zhǔn)網(wǎng)平差計(jì)算。

22、3. 2. 3精度評(píng)定單位權(quán)中誤差:(3-12)平差參數(shù)的協(xié)方差陣：D, 小一寫（3-13）A A平差函數(shù)的協(xié)方差陣：0廠元嚴(yán)AT1（3-14）3. 3水準(zhǔn)網(wǎng)間接平差程序信息設(shè)計(jì)1. 數(shù)據(jù)文件的組織下面給出一個(gè)水準(zhǔn)網(wǎng)輸入數(shù)據(jù)文件的例子：3 3 6 （己知點(diǎn)個(gè)數(shù)、未知點(diǎn)個(gè)數(shù)、觀測(cè)值個(gè)數(shù)）101 102 103 104 105 106 （點(diǎn)號(hào)）34.788 35.259 37.825 （己知點(diǎn)高程）104 101 1.625 4.5 （起點(diǎn)點(diǎn)號(hào)、終點(diǎn)點(diǎn)號(hào)、高差觀測(cè)值、距離觀測(cè)值）101 102 -0.418 3.1105 102 0.714 3.4102 103 1.243 3.8106 103

23、-0.577 4.3103 101 -0.786 2.5（其中編號(hào)數(shù)組未知點(diǎn)在前，己知點(diǎn)在后）2.水準(zhǔn)網(wǎng)平差變量約定表3-1變量約定表變量名說明ed已知點(diǎn)個(gè)數(shù)dd未知點(diǎn)個(gè)數(shù)sd總點(diǎn)數(shù)gd觀測(cè)值個(gè)數(shù)pn點(diǎn)號(hào)hO已知點(diǎn)高程be起點(diǎn)點(diǎn)號(hào)en絞占占號(hào)hl高差觀測(cè)值s距離觀測(cè)值3.4水準(zhǔn)網(wǎng)程序與使用說明3. 4. 1水準(zhǔn)網(wǎng)程序流程圖圖3-2水準(zhǔn)網(wǎng)流程圖程序的全部代碼見附錄一。3. 4. 2水準(zhǔn)網(wǎng)程序的使用本程序使用MATLAB的矩陣功能計(jì)算法方程，在運(yùn)行程序前首先要有其始數(shù)據(jù)。 其始數(shù)據(jù)是一文件的形式保存在磁盤中，文件的格式在上文己經(jīng)說明過,編好文件后, 以后綴-名為.TXT的形式保存。執(zhí)行時(shí)在MATL

24、AB命令窗口直接鍵入文件名即可。3.5案例如下圖水準(zhǔn)網(wǎng)，104、105、106為己知點(diǎn)，101、102 103為待定點(diǎn)，己知點(diǎn)的 高程分別為 34. 788, 35. 259 , 37. 825 。觀測(cè)高差和觀測(cè)路線長(zhǎng)度分別為：hl二 1.652, h2=-0. 418 ,h3=0.714 , h4=l. 243, h5=-0. 577, h6=-0. 786.sl=4 5, s2=3. 1, s3=3. 4, s4=3. & s5=4. 3, s6二2. 5.圖3-3水準(zhǔn)網(wǎng)圖首先編數(shù)據(jù)文件，命夕I為datal. txt.數(shù)據(jù)的格式如下：3 3 6101 102 103 104 105 106

25、34. 788 35. 259 37. 8251041011.6524. 5101102-0.4183. 11051020. 7143.41021031.2433.8106103-0. 5774. 3103 101 -0.786 2. 5進(jìn)入MATLAB界面，在命令窗口直接輸入level3運(yùn)行程序。彈出如下窗口Input filenameIZI兇查找范圍(I):dadadahawo航 旨自色am.txttxt.txt閏新淫文本文檔(3) out新建文本文擋文本文檔文件類型d):I*, txt打開迫)I 取消 I圖3-4數(shù)據(jù)讀入文件 選擇datal. txt即可運(yùn)行出如下結(jié)果：圖3-5計(jì)算結(jié)果在

26、圖3-5中，分別輸出了高程的平差值及精度。結(jié)果是一文本的形式保存，用戶可 對(duì)它進(jìn)行編輯。水準(zhǔn)網(wǎng)平差程序代碼function level3ed, dd, sd, gd, pn, hO, kl, k2, hl, s =readlevelnetdata;global pathname net_name s_datafile al_datafile;global ed dd sd pn gd hO kl k2 hl s dh;dh, h, V, L, uwO, uwh, uwl=calculatelevelnet (ed, dd, sd, pn, gd, hO, kl, k2, hl, s);writ

27、elevelnetdata(pn, kl, k2, hl, V, L, hO, dh, h, uwh, uwO) ;%輸岀水準(zhǔn)網(wǎng)計(jì)算結(jié)果returnfunctionLdh, h, V, L, uwO, uwh, uwl=calculatelevelnet (ed, dd, sd, pn, hO, be, en, hd, distance);%水準(zhǔn)平差網(wǎng)A二sparse (zeros (sd, gd) ;%求解系數(shù)陣b=(0: (gd-l)*sd);A (be +b)=l;A(en，+b)=l;A二A;A=A(:, 1:dd);l=zeros (gd, 1) ;%求解常數(shù)項(xiàng)l=h0(be)h0(

28、en)+hd; p=diag(l. /distance);%權(quán)陣dh=inv(A *p*A)*A，*p*l; hOO=hO(dd+l:sd); hO=hO(l:dd);h=hO+dh;V=A*dh-l;L=hd+V;%高程改正數(shù)%待定點(diǎn)高程近似值%待定點(diǎn)高程平差值%高差觀測(cè)值改正數(shù)%高差觀測(cè)值平差值%精度評(píng)定uwO=sqrt (V *p*V/ (gddd) ;%單位權(quán)中誤差Qxx二 inv(A *p*A);uwh=uwO*sqrt (diag(Qxx) ;%待定點(diǎn)高程平差值中誤差uwh(dd+1:ed+dd)二0. 0;Qff=A*Qxx*A，;%高差平差值中誤差%所有點(diǎn)高程uwh=uw0*s

29、qrt(diag(Qff); h二h;h00;ho二h0;h00;dh=dh;zeros(ed, 1);returnfunction ed, dd, sd, gd, pn, hO, kl, k2, hl, s =readlevelnetdata;global pathname net_name s_datafile b_datafile;global ed dd sd pn gd hO kl k2 hl s kll kl2: kl= ;k2= ;h= ;s=;if(isempty(pathname) isempty(net_name)filename, pathname =uigetfile

30、(J * txt , Input filename); i=find(，. =f ilename);net_name二filename (1: il);endf idl=fopen (strcat (pathname, net_name, s_dataf ile), rt)； if(fidl=-l)msgbox C Input File or Path is not correct ,，Warning, warn，);return;end%open afile to read %open afile to readed二fscanf (fidl, %f , 1) ;%已知點(diǎn)個(gè)數(shù)dd二fscan

31、f (fidl, %f , 1) ;%未知點(diǎn)個(gè)數(shù)sd 二 ed+dd;%總點(diǎn)數(shù)gd=fscanf (fidl, %F , 1) ;%觀測(cè)點(diǎn)個(gè)數(shù)pn=fscanf (f idl, %f, sd) ;%點(diǎn)號(hào)%known datahO二fscanf (fidl, %f, ed) ;%己知點(diǎn)高程hO(dd+1:ed+dd)=h0(1:ed) heightdiff=fscanf (fidl, %f , 4, gd); heightdiff二heightd辻f;kl=heightdiff (:, 1) ;%起點(diǎn)k2=heightdiff(:,2) ;%終點(diǎn)kllheightdiff (:, 1) ;%起點(diǎn)

32、%終點(diǎn)%高差%距離kl2=heightdiff(:, 2); hl=heightdiff (:,3); s=heightdiff (:,4); fclose ( all);%點(diǎn)號(hào)轉(zhuǎn)換kl, k01=chkdat (sd, pn, kl);k2, k02=chkdat(sd, pn, k2); hO(l:dd)=20000.;ie二0;while(l)%計(jì)算近似高程for k=l:gdi=kl (k);j=k2 (k);if(h0(i)le4) hO(j)=hO(i)+hl(k);ie二ie+1;if(h0(i)le4&h0(j)A（m,9）其中，m為觀測(cè)值個(gè)數(shù)，n為未知點(diǎn)個(gè)數(shù)的兩倍。改進(jìn)后的A

33、陣格式為4=（編號(hào)1,系數(shù)1,編號(hào)2,系數(shù)2,編號(hào)4,系數(shù)4,常數(shù)項(xiàng)）共9列。即只存儲(chǔ)誤差方程的4個(gè)非零參數(shù)系數(shù)。法方程系數(shù)陣為對(duì)稱陣，在存儲(chǔ)時(shí)，只需要存其上三角部分就可以了。其占用的空間為：n(n +1)sum =現(xiàn)有A陣A=（編號(hào)1,系數(shù)1,編號(hào)2,系數(shù)2,編號(hào)4,系數(shù)4,常數(shù)項(xiàng)）其中偶數(shù)項(xiàng)為系數(shù)，加上最后的A9為常數(shù)項(xiàng)，在組成法方程時(shí)，從A2開始分別與 剩下的偶數(shù)項(xiàng)以及常數(shù)項(xiàng)相乘，然后再用A4與剩余的項(xiàng)相乘，一直到A8為止，這 樣就完成了心=屮只4的過程。需要注意的是：若Al, A3, A5, A7小于零，則表 示該點(diǎn)己知點(diǎn)，不參與法方程的組成。4.1.2邊長(zhǎng)觀測(cè)的權(quán)邊長(zhǎng)觀測(cè)的精度一般與

34、其長(zhǎng)度有關(guān)，定權(quán)公式為Ps = -（j = l，2,J2）*叮式中為所測(cè)邊長(zhǎng)q的方差，1小5其中L =12 ,D= ,R = 0 d/心1)n 7由于此住A對(duì)稱性，得L = R, 乂根據(jù)A陣正定的性質(zhì)，可證明D均為正數(shù)?，F(xiàn)在設(shè)陣1申.= .d何丿d風(fēng)丿即D = D1 丄11T貝 ijA = LDlI = LD2D2lI = (LD2)(D2lI)=L L為方便，記為a=ll!稱為Cholesky分解，即正定對(duì)稱矩陣的平方根分解法。解AX=b等階于求解兩個(gè)三角方程組：LY = bL! X =Y在用平方根分解法計(jì)算時(shí)，需要進(jìn)行n次開方運(yùn)算。為了避免開方，可以直接釆用對(duì)稱正定的A = LDC分解式對(duì)

35、平方根法進(jìn)行改進(jìn)。從而解方程組AX=b可以按如下步驟進(jìn)行：把A分解成A = LDCf則AX=b變成(LDU)X=b,即等價(jià)于LY = bHX = DlY由此可以解出X和Y。這稱為改進(jìn)的平方根法，在計(jì)算中避免了開方運(yùn)算。平方根法和改進(jìn)的平方根法的計(jì)算量和存儲(chǔ)量比消去法節(jié)約近一半，而且不需要選主元，能得到比較精確的數(shù)值解。法方程用改進(jìn)平方根法解算的過程如下:(1)分解：C = STDlS(sSau3 In其中s =,D =Snm nn /S12Sljj2SllS22511S22i=l純量計(jì)算公式為5二 5, 5X -E il，ag skk(2)求逆iiInnn /由RS=I得純量計(jì)算公式:nn i

36、tn_ r22S12r23 =S33*_ r22S24 + 23534)r24 通式為1rii =；-1防-，ji（3）求積Usii_ rilrin - 44_ r22S2nh 空(一 1)$(-1”?r2n =12ri2 =；S22” _一(佔(zhàn)3 + 佔(zhàn)3)r!3 _ 一(人21“ + +Q = (S7 S 尸=SlD(Sl)T = RDRT1111Q = RDR=nn /4.1.4精度評(píng)定(1) 坐標(biāo)改正數(shù)以及單位權(quán)中誤差的計(jì)算叫使用上三角一維數(shù)組形式存儲(chǔ)坐標(biāo)改正數(shù)的公式為：Hln心=工你 -工的旳，心匕/j=其中,11=2X dd, Xi的單位是dll.平差值：X = X + x寫成分量

37、的形式，為X： = X： + 5兀如果近似坐標(biāo)的誤差較大，或網(wǎng)形較大，平差的結(jié)果不會(huì)精確，這時(shí)，就需要進(jìn)行迭代平差，直到兩次平差間互差在允許值內(nèi)。由測(cè)量平差理論：同樣可得到單位權(quán)中誤差：其中，mm觀測(cè)個(gè)數(shù)減去未知點(diǎn)個(gè)數(shù)；m = “ + m2 +加3, = 2 x dd + ST, ST -方向觀測(cè)的測(cè)站數(shù)PW = pll + Sxco(2) 點(diǎn)位誤差橢圓誤差橢圓表示了網(wǎng)中點(diǎn)或點(diǎn)與點(diǎn)之間的誤差分布情況如圖。在測(cè)量工作中，常用的誤差橢圓對(duì)布網(wǎng)方案作精度分析。繪制誤差橢圓只需要三個(gè)數(shù)據(jù)：橢圓長(zhǎng)半軸a,短半軸b和主軸方向4),其求法為宀 + （時(shí) + b J + g2-b：），+4吋）2b 卜 tan

38、 20 二jGj + _ J（Qxx-QQ+4Q）根據(jù)上述的理論，我們實(shí)際要求的是5、加、八 加如。只要得到了這些元素， 就能依照上面的公式來求得橢圓的元素了。4. 2 測(cè)邊網(wǎng)平差信息設(shè)計(jì)外業(yè)測(cè)量的數(shù)據(jù)首先應(yīng)進(jìn)行預(yù)處理，包括測(cè)站平差、歸心計(jì)算、觀測(cè)值歸化到橢 球面的改正、橢球面歸化到高斯平面的改正等，然后將預(yù)處理后的數(shù)據(jù)輸入到以后綴 名為TXT的文本文件中，該數(shù)據(jù)文件的組織格式如下所列：表4-1數(shù)據(jù)組織格式次序內(nèi)容1已知點(diǎn)個(gè)數(shù)ed,未知點(diǎn)個(gè)數(shù)dd,控制參數(shù)2點(diǎn)號(hào)pn先輸入已知點(diǎn)編號(hào)，各點(diǎn)輸入順序無要求3己知點(diǎn)坐標(biāo) 5xO,yO,xl,yl,x2,y2 4測(cè)量邊的個(gè)數(shù)ml5邊長(zhǎng)的固定誤差ms,邊

39、長(zhǎng)的比例誤差pp。單位分別是:cm和 X 0.0000016邊長(zhǎng)的起始點(diǎn)號(hào)6終點(diǎn)點(diǎn)號(hào)d、邊長(zhǎng)skL每一條邊一行，依次列岀7推算近似坐標(biāo)的路線經(jīng)過的邊數(shù)8推算近似坐標(biāo)的起算已知點(diǎn)坐標(biāo)（按順時(shí)針）9推算近似坐標(biāo)的路線經(jīng)過的邊的邊號(hào)4.2.1主要的技術(shù)要求表4-2測(cè)邊網(wǎng)的技術(shù)要求等級(jí)平均邊長(zhǎng)（km）測(cè)距中誤差（mm）測(cè)距相對(duì)中誤差二等9301/30 萬三等5301 / 16 萬四等2161/12 萬一級(jí)小三角1161/6萬二級(jí)小三角0.5161/3萬4. 3利用MATLAB的繪圖語句繪制網(wǎng)圖1. 網(wǎng)形的繪制由于網(wǎng)形圖與誤差橢圓繪制在同一個(gè)圖形上，因此必須對(duì)誤差橢圓進(jìn)行放大，在 本文的程序中，使用了

40、inputdig對(duì)話框輸入語句，其中，確省的放大倍數(shù)為100。在 程序中，使用了 ed dd pnml xy e d sid ai bi fi等變量，其意義與前面的變量相同。對(duì) 繪制的網(wǎng)圖有效放大和縮小功能，即點(diǎn)擊鼠標(biāo)左鍵放大圖形，點(diǎn)擊右鍵縮小圖形，利 用MATLAB菜單本身的功能，可以將該圖形輸出為各種圖形文件格式。2誤差橢圓的繪制無論多么復(fù)雜的圖形，其基本單元還是點(diǎn)和線。換句話說，只需要利用基本元素 的點(diǎn)或線，通過各種組合，也能畫出復(fù)雜的圖形。MATLAB中沒有提供直接繪制橢 圓的命令，因此可以直接利用連線來畫橢圓。測(cè)量中描述誤差橢圓用長(zhǎng)半軸A,短半軸B和方位角FI三個(gè)量。在如圖4-1的

41、x Oyf中直角坐標(biāo)系中，橢圓的標(biāo)準(zhǔn)方程為(4-1)(4-2)如果以角度1為變量，則橢圓的標(biāo)準(zhǔn)參數(shù)方程為 = Acosfy = Bsnu設(shè)在測(cè)量坐標(biāo)系xOy中橢圓長(zhǎng)半軸的方位角為，則有COS %,Sill %-S1110OCOS 00(4-3)2用參數(shù)方程代入，得到(4-4)x= A cos(pQ cos z B siii(p0 sin iy = A sin(pQ cosz + B cos(pQ sin i測(cè)量坐標(biāo)系與MATLAB或AUTOCAD繪制的數(shù)學(xué)坐標(biāo)系的x,y坐標(biāo)軸不同，繪 圖時(shí)，需要調(diào)換x, y坐標(biāo)。在上式中，1取不同的值，就有一組（x, y）,只需要將這些點(diǎn)連接起來，就可以繪出

42、一個(gè)橢圓。Y圖4-2誤差橢圓的參數(shù)方程4.4測(cè)邊網(wǎng)程序和使用說明使用本程序的全部數(shù)據(jù)都按規(guī)定的格式編輯成數(shù)據(jù)文件儲(chǔ)存在磁盤上。數(shù)據(jù)文件 的編輯取決于平差網(wǎng)型和觀測(cè)值的編號(hào)，為此，先繪制平差網(wǎng)的略圖，在圖上標(biāo)明各 項(xiàng)數(shù)據(jù)的信息。以下就測(cè)邊網(wǎng)的點(diǎn)號(hào)，邊編號(hào)，輸入數(shù)據(jù)，輸出成果，運(yùn)行程序等問 題作簡(jiǎn)明說明。1點(diǎn)號(hào)和觀測(cè)邊的編號(hào)己知點(diǎn)和待定點(diǎn)的編號(hào)為pn （取三為數(shù)），己知點(diǎn)在前，未知點(diǎn)在后，其順序 無要求。但為了減小法方程系數(shù)的帶寬，應(yīng)使相臨的待定點(diǎn)編號(hào)的差數(shù)盡可能小。平 差網(wǎng)的編號(hào)見圖4-3。2. 推算近似坐標(biāo)的路線近似坐標(biāo)的路線是用戶在測(cè)邊網(wǎng)略圖上指定出來的。如圖4-3的箭頭，就是表示 推算路線

43、。路線的兩個(gè)起算點(diǎn)必須為己知點(diǎn)，從兩個(gè)己知推算出的第一個(gè)未知點(diǎn)開始, 選擇觀測(cè)邊，由觀測(cè)邊和己求得的近似坐標(biāo)或己知坐標(biāo)推算出觀測(cè)邊所對(duì)的未知點(diǎn)。 本程序是按推算的三個(gè)點(diǎn)A、B、P順序?yàn)轫槙r(shí)針。圖4-3測(cè)邊網(wǎng)略圖2數(shù)據(jù)的輸入(1) 簡(jiǎn)單變量為了在程序運(yùn)行中數(shù)據(jù)的傳遞，定義了一些全局變量，參見4-1表。(2) 數(shù)據(jù)文件外業(yè)測(cè)量的數(shù)據(jù)首先應(yīng)進(jìn)行預(yù)處理，包括測(cè)站平差、歸心計(jì)算、觀測(cè)值歸化到橢 球面的改正、橢球面歸化到高斯平面的改正等，然后將預(yù)處理后的數(shù)據(jù)輸入到以后綴 名為TXT的文本文件中，該數(shù)據(jù)文件的組織格式如表4-1。3輸出成果本程序的計(jì)算成果是以文件的形式輸出到文本文件中，用戶可以在文本中查看和

44、 編輯。4.5程序代碼說明：程序的總體框架數(shù)據(jù)的諛入近似坐標(biāo)計(jì)算誤差方程和法方程的組成解算法方程精度評(píng)定輸出計(jì)算結(jié)杲繪制控制網(wǎng)圖和誤差橢圓圖4-4測(cè)邊網(wǎng)總體流程圖1. 數(shù)據(jù)讀入塊本模塊的功能是打開一個(gè)*.txt的數(shù)據(jù)文件，同時(shí)生成一個(gè)Sllt.txt的文本文件， 記錄用戶數(shù)據(jù)和輸出成果用。程序調(diào)用fscanf函數(shù)把文件中的數(shù)據(jù)賦值給相應(yīng)的變 量，這些變量是后面計(jì)算的數(shù)據(jù)依據(jù)。2. 近似坐標(biāo)的計(jì)算測(cè)邊網(wǎng)的觀測(cè)數(shù)據(jù)是邊長(zhǎng)，所以在近似坐標(biāo)的計(jì)算時(shí)只能用測(cè)邊交會(huì)計(jì)算。以及 A、B兩點(diǎn)坐標(biāo)及A、B到P點(diǎn)的距離b、a, c為A、B兩點(diǎn)的距離，A、B、P三點(diǎn) 按順時(shí)針排列。則P點(diǎn)的坐標(biāo)計(jì)算公式如下： rJ

45、= y/b2-e2Zr + L-cr2cXp = XA+ecos a * - f sin aAB Yp = Xa +esinaAB + f cosaAB計(jì)算近似坐標(biāo)的流程圖入圖4-5o輸入A、E點(diǎn)號(hào)for i = l: length (m)一 一一一蘭 1-嚴(yán)巡之|巳（m （匕書亠彳f*dm (il) )=A| dF二d(m(i)P二d(m(i)B=PA=PJ一仝）二AP=d(w(i)P=e (m(i)計(jì)算近似坐標(biāo)h0（F）,70（P）圖4-5計(jì)算近似坐標(biāo)流程圖3. 誤差方程和法方程的形成程序中，用數(shù)組a來存儲(chǔ)誤差方程的編號(hào)和系數(shù)，a（i,9）存儲(chǔ)常數(shù)項(xiàng)，w和c分別 存儲(chǔ)法方程系數(shù)和法方程常數(shù)

46、。4. 解算法方程。函數(shù)流程圖如圖4-6。函數(shù)的代碼見附錄三。5. 精度評(píng)定本模塊包括坐標(biāo)改正數(shù)、單位權(quán)中誤差和誤差橢圓的計(jì)算。程序中定義了 dxy（坐 標(biāo)改正數(shù)）、pw （即存儲(chǔ)PW = pll + Sx.co ）、uwO（單位權(quán)中誤差）等。同時(shí)計(jì)算 出誤差橢圓的三個(gè)參數(shù)長(zhǎng)半軸ai,短半軸bi和主軸方向fiofor i=l:ndi=Crri/2. 0)Yfor j=l:nss=C(di+j)for k=l:i-l dk=(lri)*Gn-V2. 0) ss=S-C(dk+i)*C(dij)/C(dki-k)L=J?c (di4j)=l/ssc (d計(jì) j)二 ee*CsCCdij)for k

47、=j-*-l:nss=C(di+j)*C(dj+j)dkz(lrl)*(ii-k/2.0)s 5=55吃di恢桁時(shí)林)札(血松)c (di4-j)=ssfor i=l:nl di=(i-1) *(ni/2. 0)fcr i=l:irldi=Crri/2. 0)for i=j+l:n ss=CCdi+j)foi k=i+l: j-1 dk=(Jrl)*(n-V2. 0) sszss-Cc (di+k) *C (dk+j)圖4-6平方根法求逆程序框架圖6.控制網(wǎng)的繪制程序中調(diào)用inputdig函數(shù)來打開一個(gè)對(duì)話框來輸入誤差橢圓的參數(shù)；用text函數(shù) 來對(duì)點(diǎn)號(hào)標(biāo)注；用ploth函數(shù)及控制參數(shù)繪制線和

48、誤差橢圓。程序的完全代碼見附錄二。4.6程序的使用算例有測(cè)邊網(wǎng)如圖4-3所示。網(wǎng)中A、B、C、C點(diǎn)己知，其余為未知點(diǎn)，現(xiàn)用某測(cè)距 儀觀測(cè)了 13條邊長(zhǎng)，測(cè)距精度s = (3m + 1x10Ys)。起算數(shù)據(jù)和觀測(cè)數(shù)據(jù)如下：表4-3已知點(diǎn)坐標(biāo)點(diǎn)需XY點(diǎn)名XYA53743.13661003.826C40049.22953782.79B47943.00266225.854D36924.72861027.086表4-4觀測(cè)值編號(hào)邊觀測(cè)值編號(hào)邊觀測(cè)值編號(hào)邊觀測(cè)值15760.70668720.162115487.07325187.34275598.57128884.58737838.8887494.881137228.36745483.15897493.32355731.788105438.3821 編輯數(shù)據(jù)文本文件如下圖圖4-7數(shù)據(jù)文件2.在MATLAB命令窗口鍵入nnbb執(zhí)行程序，運(yùn)行中會(huì)彈出一個(gè)對(duì)話匡提示用戶輸入誤差橢圓的放大比例，默認(rèn)為100,本例選擇500。如下圖竺i叩山the cab10D-OK Cancel圖4-8b對(duì)話框圖4-8a對(duì)話框3. 計(jì)算出的結(jié)果如下：在圖4-9a中，第1、2行分別是己知點(diǎn)