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

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

2、計ABSTRACTMATLABisonespeciesofnumerical-valuescalculationandgraphictoolssoftwarewhichiswidelyusedtoapplyatresearchinstitutionsatpresent.Theparticularitiesare:concisegrammar-structure、highlyefficientinnumericalvaluescalculating、completefunctionofgraphs、especiallyitisadaptedtoevildoingprofessionalprogr

3、ammertoaccomplishthetasksthatarenumerical-valuescalculatingandscientificexperimentstreating.Theancientmethodsofmeasureddata-processingneedestablishingspecialproceedingsoftreatingmatricesoperation,moreover,itiscomplexandgreatlydifficult.Thisarticleintroducesoneprogrammingmethoddealingwithlevelingandm

4、easuringedgenetworkbasedonMATLAB.Comparedwithotheralgorithmlanguage,ithasparticularitieswhicharesimplyprogrammingandquicklyoperating.Thearticleseparatelyexpatiatethetheoriesbasics、realizingstepsandrunningresultsatlevelingandmeasuringedgenetwork.Withtheanalysisofexamples,ithasprodigiousapplicationval

5、ueinmeasureddata-processingbyuseofMATLAB.Moreover,itshortensprogrammingtimeandimprovesworkingeffectiveness.Keywords：MATLAB；programming緒論作為一名測量技術(shù)人員，如果不掌握一門PC機編程語言與便攜計算工具，要想提高測量工作的效率幾乎寸步難行。測量需求的多樣性與復(fù)雜性，造就了測量計算鮮明的個性化特點，這就是在商業(yè)測量計算軟件高度發(fā)達的今天，掌握一種實用的程序語言進行編程計算仍有廣泛的市場需求的重要原因。當今較流行的計算機程序語言基本上都是基于Windows的，例如T

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

7、測量技術(shù)人員最好能熟練掌握一種適用于數(shù)值計算的程序語言，以便提高測量計算的效率。C語言的數(shù)值計算語句不夠豐富，例如，在測量平差計算中，經(jīng)常需要進行的矩陣運算，尤其是解法方程的矩陣求逆不能直接使用語句實現(xiàn)，而必須應(yīng)用計算機算法編程實現(xiàn)。如果不是基于商業(yè)軟件開發(fā)，只為滿足實際測量工作計算需要，則C語言的劣勢就變成了MATLAB語言的優(yōu)勢。內(nèi)蒙古科技大學(xué)畢業(yè)設(shè)計MATLAB軟件簡介1.MATLAB軟件簡介MATLAB是從Matrix（矩陣）和Laboratory（實驗室）各取前3個字母組成的,意思是矩陣實驗室,是美國MathWorks公司于20世紀80年代中期推出的一種交互式、面向?qū)ο蟮目萍紤?yīng)用軟件

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

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

10、內(nèi)部函數(shù)之外，所有的MATLAB主包文件和各工具包文件都是可讀可改的源文件，用戶可以作為參考掌握其用法，并可對其修改以適應(yīng)自己的需要，也可加入自己編寫的文件構(gòu)成新的工具包。例如，隨著GPS的廣泛應(yīng)用,OrionDynamicsandCon21rolCorporation、ConstellInc.GPSSoftLLC、NavsysCorporation等多家公司都相應(yīng)開發(fā)出了適于GPS數(shù)據(jù)處理的MATLAB工具箱。MATLAB是一個集數(shù)值計算、圖形管理、程序開發(fā)于一體的功能十分強大的系統(tǒng)。將MATLAB應(yīng)用于測量數(shù)據(jù)的處理是一件非常有意義的工作Mo2hamed等曾成功地在MATLAB系統(tǒng)中利用白

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

12、原理理解和掌握變得更容易。另外,MATLAB語言與數(shù)學(xué)語言比較接近，更容易掌握和理解。2.1測量平差原理的概述測量平差的函數(shù)模型有條件方程和觀測方程。以條件方程為函數(shù)的模型的最小二乘平差稱為條件平差；在條件方程中，根據(jù)需要如果還設(shè)有一定數(shù)量的未知數(shù)，則稱為附有參數(shù)的條件平差；以觀測方程為函數(shù)模型的最小二乘平差稱為間接平差；如果觀測方程中的某些參數(shù)不獨立，則這些不獨立參數(shù)必然存在一些條件，稱這種平差模型為附有條件的間接平差。本文的兩個程序都采用間接平差模型。對于一個實際平差問題，根據(jù)所選參數(shù)的個數(shù)、選什么量為參數(shù)以及參數(shù)之間是否函數(shù)獨立，經(jīng)過仔細推敲可以發(fā)現(xiàn)附有條件的間接平差模型本身就是各種經(jīng)典

13、平差模型的概括模型，其余的經(jīng)典平差模型，如條件平差模型、間接平差模型、附有未知數(shù)的條件平差模型和附有限制條件的條件平差模型都是它的特例。間接平差的公式匯集：間接平差模型為V二BX-l21）（2-2）（2-3）VTPVTmin系數(shù)矩陣B滿秩，即rank（B）=t法方程及解為：Nx-f二0,（N-i二BTPB,f二BtPI）bbebbex=N-ifbbe參數(shù)的平差值：X二X+x0（2-4）觀測量的平差值：L二L+V（2-5）單位權(quán)的中誤差：花占(2-6)平差參數(shù)的協(xié)方差陣：D=b2N-1(2-7)XX0bb平差函數(shù)的協(xié)方差陣：Q2FtN-iF(2-8)榊0bb2.2平差程序總體方案MATLAB號稱

14、為全球工程師的共同語言，其語法和C語言相似，但它有強大的數(shù)值計算和繪圖功能，這使之在工程應(yīng)用方面的計算更出色，本文就基于這種程序設(shè)計語言環(huán)境設(shè)計一個控制網(wǎng)平差程序。該程序包含了一個高程控制網(wǎng)平差程序和測邊網(wǎng)平差程序。本程序適用于各種等級的高程網(wǎng)和測邊網(wǎng)，程序在設(shè)計過程中，始終考慮數(shù)據(jù)的儲存量。因而本程序不儲存誤差方程的系數(shù)和常數(shù)項，對待定點數(shù)較多的平差網(wǎng)，組成法方程的系數(shù)矩陣是個稀疏矩陣，如待定點的編號恰當，法方程的系數(shù)會集中在主元系數(shù)的兩側(cè)形成帶狀。為減少法方程系數(shù)的儲存量，只要按行儲存下三角陣或按列儲存上三角陣中第一個非零系數(shù)起的系數(shù)，就是通常叫做維變帶寬儲存方法。3.水準網(wǎng)平差程序3.1

15、程序的功能本程序適用于二、三、四等水準網(wǎng)平差計算，平差的水準網(wǎng)可以是獨立的、也可以是附合網(wǎng)，其主要功能是完成水準網(wǎng)的平差計算和精度評定計算。平差計算采用間接平差法，以歸心的觀測值為高差，以未知點高程為未知參數(shù)。精度評定計算包括計算單位權(quán)中誤差和每個待定點的高程中誤差。3.2水準模型網(wǎng)的間接平差3.2.1“權(quán)”值的確定當在相同的條件下進行水準測量時，其精度是相同的，因而觀測結(jié)果的可靠性也是同樣的。但如果在不同的條件下進行水準測量時，高程的精度就有所不同，此時稱為不等精度觀測，所求出的未知量的值、高程的最或是值并對其精度進行評定時，就需要“權(quán)”了。由于觀測的不等精度，因而觀測值的可靠程度不同，求未

16、知量的最或是值時，這樣的一個因素就必須考慮了，這個因素是:可靠性大的某觀測值，其精度高，對測量的最后結(jié)果的影響也就越大。此時用“權(quán)”值來表示觀測值的可靠程度那么，“權(quán)”值愈尢觀測值的可靠程度就愈高。另外，在觀測過程中,觀測值的中誤差愈小，觀測結(jié)果愈可靠,它的“權(quán)”值就愈大。因而根據(jù)中誤差來確定“權(quán)”值是非常適當?shù)摹ＴO(shè)以Pi表示觀測值Li的“權(quán)”,m為中誤差，則“權(quán)”值的定義為：P.=呂(3-1)式中:A為任意的正常數(shù)，在一組觀測值中為一個定數(shù)。在實際測量中，通常是觀測值的中誤差事先并不知道，因而必須先確定觀測值的“權(quán)”,然后才能求出未知量的最或是值。此時可以利用距離(S)或測站數(shù)(N)來確定觀

17、測值一高程的“權(quán)”。根據(jù)偶然誤差傳播定律,各觀測點高程Hi的中誤差mi由測站數(shù)Ni確定時，則有：(3-2)式中:m為一組觀測值的中誤差,為一個定數(shù).由(3-2T)、(3-2-2)兩式可得：P_C(3-3)式中上二亠.由于衛(wèi)5邯為定數(shù)故C也為定數(shù).7U同樣可得出:(3-4)式中:C為定數(shù)，s為測距.i由(3-3)、(3-4)兩式可以得出這樣一個結(jié)論:當測站觀測高差等精度時，觀測總高差的“權(quán)”與測站數(shù)或距離成反比。3.2.2水準路線的平差計算1.附合路線的平差計算假定在圖1示的A、B兩水準點之間布設(shè)一條水準路線,A、B兩水準點的高程為已知，分別設(shè)為H、H、n、n、C為中間水準點。假定觀測了所有的點

18、的高TOC o 1-5 h zAB12程,現(xiàn)擬求C點的高程H的最或是值。CH可由水準路線AfC、BfC分別觀測的高差A(yù)h、Ah計算得出，由此而CACBC得到的觀測高程分別設(shè)為Hcl、Hc2，其值為：Hc1=H+Ah;Hc2=H+AhAACBBC當Hcl、Hc2在不等精度條件下觀測得出時,它們的“權(quán)”也不同，分別設(shè)為Pel、Pc2,這樣C點的高程H的最或是值為：C(3-5)”PHPHTOC o 1-5 h zH=clclc21cP+Pclc2根據(jù)A點的高程H,AfC水準路線觀測的高差A(yù)h以及BfC水準路線觀測AAC的高差A(yù)h，可推算出B點的觀測高程H為：BCBH=H+Ah-AhBAACBC水準路

19、線AfB的高程閉合差為：f=H-H=HH（3-6）hBBclc2由（3-6）式得到：H=H-fc2clh HYPERLINK l bookmark82 o Current Document CC由（3-3）式得到：P=亠P=亠（N、N分別表示水準路線AfC、BclN、c2NACBCfC的測站數(shù)，水準路線AfB的測站數(shù)N二N+N）ABACBC將上述表達式代入（3-2-5）式中,得到：(3-7)如果以水準路線AfC的距離S、BfC的距離S、AfB的距離SBCBCAB(S二S+S)ABACBC來確定高程觀測值的“權(quán)”值時，同樣可以得到:H.圖3-1水準路線圖2閉合路線的平差計算閉合路線的平差計算原理

20、與附合路線相同，因而(3-7)、(3-8)兩式的結(jié)論適用于閉合路線的平差計算。(3)具有一個結(jié)點的水準網(wǎng)的平差計算工PHAiAiXpAii=1(3-9)如圖2所示為具有一個結(jié)點的水準網(wǎng),B,C,D,為已知高程水準點,B-A,C-A,D-A,為水準路線，則接點A的高程最或是值為：口PH+PH+PH+H二一A2A2A3A3AP+P+P+A1A2A3式中H,H,H分別為水準路線B-A,C-A,D-A,計算A的觀測高程，各高A1A2A3程相應(yīng)的“權(quán)”值為P,P,PA1A2A3設(shè)H,H,H的算術(shù)平均值為H0，各高程觀測值與H0的差值分別為SA1,A1A2A3ASA2,5A3,則有：41=耳-找Jtl=戌

21、+越1狂=驗-戌或毘2=&2&2=屁_出H石=止+(3-10)將(3-10)式代入(3-9)式得到:(3-11)當以測站數(shù)和距離來確定“權(quán)”值時，（3-11）分別可以轉(zhuǎn)化為：耳=疋+于當，（3-12）耳=此十占；（3-13）上述結(jié)論也可應(yīng)用于小三角水準網(wǎng)平差計算。3.2.3精度評定單位權(quán)中誤差：bVTPV0Vnt(3-12)平差參數(shù)的協(xié)方差陣：D-b2N-1XX0bb(3-13)平差函數(shù)的協(xié)方差陣：Q-b2FtN-1F林0bb(3-14)3.3水準網(wǎng)間接平差程序信息設(shè)計1.數(shù)據(jù)文件的組織下面給出一個水準網(wǎng)輸入數(shù)據(jù)文件的例子：336（已知點個數(shù)、未知點個數(shù)、觀測值個數(shù)）1011021031041

22、05106（點號）34.78835.25937.825（已知點高程）1011.6254.5（起點點號、終點點號、高差觀測值、距離觀測值）101102-0.4183.11020.7143.41031.2433.8106103-0.5774.3101-0.7862.5（其中編號數(shù)組未知點在前，已知點在后）2.水準網(wǎng)平差變量約定表3-1變量約定表變量名說明ed已知點個數(shù)dd未知點個數(shù)sd總點數(shù)gd觀測值個數(shù)pn點號h0已知點高程be起點點號en終占占號k、八、八、Jh1咼差觀測值s距離觀測值3.4水準網(wǎng)程序與使用說明3.4.1水準網(wǎng)程序流程圖圖3-2水準網(wǎng)流程圖程序的全部代碼見附錄一。3.4.2水準

23、網(wǎng)程序的使用本程序使用MATLAB的矩陣功能計算法方程，在運行程序前首先要有其始數(shù)據(jù)。其始數(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.&s5=4.3,

24、s6=2.5.圖3-3水準網(wǎng)圖首先編數(shù)據(jù)文件，命名為data1.txt.數(shù)據(jù)的格式如下：33610210310410510634.78835.25937.8251011.6524.5102-0.4183.11020.7143.41031.2433.8103-0.5774.3101-0.7862.5進入MATLAB界面，在命令窗口直接輸入level3運行程序。彈出如下窗口Inputfilename查我范團Q：口work_T十匡呼datal.如七data3.txt0新建文本文襠1.3Jout0新建丈本文擋。毗U1打開|*.txtJ取消文件名逼：文件類型(I):圖3-4數(shù)據(jù)讀入文件選擇datal.t

25、xt即可運行出如下結(jié)果：FileEdi+TeKtCellToolsDebugDesktopWindowHelp站a#4A4*VSStack：Base田待定點高程值及中誤差:點號近似高程咼差高程平差値高程平差値中誤差pnh0(m)dh(m)h(m)vwh10136.4400-0.008036.4320C.016210236.0220-0.029535.9925C.015910337.2650-0.034137.2309C.016310434.78800.000034.7880C.000010535.25900.000035.25900.000010637.82500.000037.82500.0

26、000高差觀測値平差値:plainfileLn1Col1DV1pnpnh0(m)V(m)H(m)1041011.6520-0.00801.6440101102-0.4180-0.0215-0.43951051020.71400.01950.73351021031.2430-0.00461.238106103-0.5770-0.0171-0.5941103101-0.7860-0.0129-0.7989D12345678910111213141516171819單位權(quán)中誤差：0.0118m圖3-5計算結(jié)果在圖3-5中，分別輸出了高程的平差值及精度。結(jié)果是一文本的形式保存，用戶可對它進行編輯。水準

27、網(wǎng)平差程序代碼functionlevel3ed,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata;globalpathnamenet_names_datafilea1_datafile;globaledddsdpngdh0k1k2h1sdh;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é)果returnfunctiondh,h,V,L,uw0,uwh,

28、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;%求解常數(shù)項A=A(:,1:dd);l=zeros(gd,1);%權(quán)陣%高程改正數(shù)%待定點高程近似值%待定點高程平差值%高差觀測值改正數(shù)%高差觀測值平差值l=h0(be)-h0(en)+hd;p=diag(1./distance);dh=inv(A*p*A)*A*p*1;h00=h0(dd+1:sd);h0=h0(1:dd);h=h0+

29、dh;V=A*dh-1;L=hd+V;%精度評定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*sqrt(diag(Qff);%高差平差值中誤差h=h;h00;%所有點高程h0=h0;h00;dh=dh;zeros(ed,1);returnfunctioned,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata;globalpathnamenet_names_datafi

30、leb_datafile;globaledddsdpngdh0k1k2h1sk11k12;k1=;k2=;h=;s=;if(isempty(pathname)|isempty(net_name)filename,pathname=uigetfile(*.txt,Inputfilename);i=find(.=filename);net_name=filename(1:i-1);endfid1=fopen(strcat(pathname,net_name,s_datafile),rt);if(fid1=-1)msgbox(InputFileorPathisnotcorrect,Warning,w

31、arn);return;end%openafiletoread%openafiletoreaded=fscanf(fid1,%f,1);%已知點個數(shù)dd=fscanf(fid1,%f,1);%未知點個數(shù)sd=ed+dd;%總點數(shù)gd=fscanf(fid1,%f,1);%觀測點個數(shù)pn=fscanf(fid1,%f,sd);%點號%knowndatah0=fscanf(fid1,%f,ed);%已知點高程h0(dd+1:ed+dd)=h0(1:ed)heightdiff=fscanf(fid1,%f,4,gd);heightdiff=heightdiff;k1=heightdiff(:,1);

32、%起點k2=heightdiff(:,2);%終點k11=heightdiff(:,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(l)%計算近似高程fork=1:gdi=k1(k);j=k2(k);if(h0(i)le4)h0(j)=h0(i)+h1(k);ie=ie+1;endif(h0(i)le4&h0(j)A（

33、m,9）其中，m為觀測值個數(shù)，n為未知點個數(shù)的兩倍。改進后的A陣格式為A=（編號1,系數(shù)1,編號2,系數(shù)2,編號4,系數(shù)4,常數(shù)項）i共9列。即只存儲誤差方程的4個非零參數(shù)系數(shù)。法方程系數(shù)陣N為對稱陣,在存儲時,只需要存其上三角部分就可以了。其占A用的空間為：n（n+1）sum=2現(xiàn)有A陣：A=（編號1,系數(shù)1,編號2,系數(shù)2,編號4,系數(shù)4,常數(shù)項）其中偶數(shù)項為系數(shù)，加上最后的A9為常數(shù)項，在組成法方程時，從A2開始分別與剩下的偶數(shù)項以及常數(shù)項相乘,然后再用A4與剩余的項相乘,一直到A8為止,這樣就完成了N=AtPA的過程。需要注意的是：若A1,A3,A5,A7小于零，則表A示該點已知點,不

34、參與法方程的組成。4.1.2邊長觀測的權(quán)邊長觀測的精度一般與其長度有關(guān),定權(quán)公式為C2P二0-siQ2si（i二12，n）式中O2為所測邊長s的方差，a2為任意選定的單位權(quán)方差。si0i為了定權(quán)P必須已知測邊的先驗方差a2sisi但精確的已知是十分困難的,一般采用廠方給定的測距儀精度,即a=a+bSsii式中，a為固定誤差（單位mm）,b為比例誤差（單位：ppm）,S為邊長（單i位km)。4.1.3解算法方程由于法方程是對稱正定陣，因此，可采用改進的平方根法進行解算。平方根法是對稱正定矩陣非常有效的三角分解方法，設(shè)A為n階方陣，如果其所有順序主子式均不為零，則其存在唯一的分解式：A=LDR其中

35、L=ri1.2，D=rdi0、，R=rirr12.in1ni11n,n-1丿10d丿nrn-i,n1丿由于此住A對稱性，得L=RT，又根據(jù)A陣正定的性質(zhì)，可證明D均為正數(shù)?，F(xiàn)在設(shè)(ddiD=d即111丄，A=LDLt=LD2D2Lt=(LD2)(D2Lt)=LLT則為方便A=LLt稱為Cholesky分解，即正定對稱矩陣的平方根分解法。解AX=b等階于求解兩個三角方程組：LY=b和LtX=Y在用平方根分解法計算時，需要進行n次開方運算。為了避免開方，可以直接采用對稱正定的A=LDL分解式對平方根法進行改進。從而解方程組AX=b可以按如下步驟進行：把A分解成A=LDLt，則AX=b變成(LDLT

36、)X二b，即等價于LY=bLtX=D-iY由此可以解出X和Y。這稱為改進的平方根法，在計算中避免了開方運算。平方根法和改進的平方根法的計算量和存儲量比消去法節(jié)約近一半，而且不需要選主元，能得到比較精確的數(shù)值解。法方程用改進平方根法解算的過程如下(1)分解：C=StD-is其中1S11S1n,D=1S11kS丿mnkS丿mn1c11c)1ns-2-s111S11S)1ncn1C丿nnS丿mns-ns11n1,72sn-1,n-1s1j=c1jssc2j12_Ljsc3js11ss131j2js2jc2js11s11s+232j+ss223j,s3j=c3js131js1122cijs1jcij-

37、1k=1ssskki1,j2純量計算公式為i=12)求逆R=S-11r11r、1nR=r丿nn由RS=I得1r11r12r22r1nr2nsiis12s22rnn227純量計算公式：r=一iisiirsr=ii1212s22(rs+rs)r二1113122313s33(rs+.+rsr=111n1(n1)(n1)n1nsnnrsr=221223s33(rs+rs)r=22_24233424s44(rs+rsr=222n2(n1)(n1)n2nsnn通式為1r=iisii玄rsikkjr二一i,jiijsij3）求積Q=(StD-1S)-1=S-1D(S-1)t=RDRtQ=RDRt廠sr111

38、1s12sr11r1n(s11rnn八snnr11r1nr12r22sr1n1nsr2n2nr11rnnr丿nn、rnnsrnnnnr1nrnnrnn4.1.4精度評定（1）坐標改正數(shù)以及單位權(quán)中誤差的計算m0使用上三角一維數(shù)組形式存儲坐標改正數(shù)的公式為：5x=一可qw-工qw,i=1,nijijjjj=1j=1其中，n=2xdd,x的單位是cm.i平差值：X=Xo+x寫成分量的形式，為X=X0+5xiii如果近似坐標的誤差較大，或網(wǎng)形較大，平差的結(jié)果不會精確，這時，就需要進行迭代平差，直到兩次平差間互差在允許值內(nèi)。由測量平差理論：=Vtpv.n-1同樣可得到單位權(quán)中誤差：fiPVVTm=v-

39、0m一n其中，m-n觀測個數(shù)減去未知點個數(shù)；m=m+m+m,n=2xdd+ST,ST一方向觀測的測站數(shù)123PVV=pll+5x（2）點位誤差橢圓誤差橢圓表示了網(wǎng)中點或點與點之間的誤差分布情況如圖。在測量工作中，常用的誤差橢圓對布網(wǎng)方案作精度分析。繪制誤差橢圓只需要三個數(shù)據(jù)：橢圓長半軸a，短半軸b和主軸方向，其求法為2otan2p=xyo2o2xya2=(O2+O2+2o2)2+402)2xyxyxyb2=丄(02+o2一(o2-o2)2+4o2)2xy耳xyxy1793.00266225.85440049.22953782-79036924.72861027.S8613760.7061055

40、187.342065483.1587838.S80061051017493.3237494.8810304108058884.5：871QS081021087228.36721012圖4-7數(shù)據(jù)文件2.在MATLAB命令窗口鍵入nnbb執(zhí)行程序，運行中會彈出一個對話匡提示用戶輸入誤差橢圓的放大比例，默認為100,本例選擇500。如下圖圖4-8a對話框圖4-8b對話框3計算出的結(jié)果如下：在圖4-9a中，第1、2行分別是已知點和未知的個數(shù)；第4行是點的編號；第5至第8行是已知點的坐標；第10行是觀測值的個數(shù)；第12行是測距的固定誤差和比例誤差；13至25行是觀測邊的起點號、終點號和觀測邊長；27至

41、39行（接到圖4-9b）是點號轉(zhuǎn)換為計算機順序后的觀測邊起點號、終點號和觀測邊長。接圖4-9a，在圖4-9b中，第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é)果-lai如IjjslsM|豪詆耳53me|jH1出畫1程）51netp田m曰日百32百3斷述IGZ3331?wfit.&ra3447TW.EE35S7T4SL32336T5開遲3B33?35翻1&5&2E

42、wsF.oraS$34,5&r72S,33FtoUii鼻的腔貝止標125370.1360JDDXB2(3a丁齟乳D02S6225.S2.WO15；淞如.理Q閃亞.DM16isseo.2?0和5汕.5師4：.390SSO30.2Mta43762225T?6S.592戲E師.BM3D51訝農(nóng)方程弔加和糟功52-&.-D.11L-a.DDDD.SS4I.DOQ乩ILL2-HH-CL&Mdcen53-r.DDDD.QQE-C.ODDD.D9T1.000-O.fliK2-QXi-CLOS?CLCOD54-r.ODDD.NE-fi.ODDD.Tfi4n.oao-O.SJC帚cm-

43、clr&i-CLCOD551.-D.ma2.DDD1.DDdn.oao-DLS帚cm-L-iXO-CLCOD5G3.DDDD.B574.DDD-D.ELEE.aao-D.3E?fi.SBX53EQ.CGD573.ODD.99Da.DDDD.42-.oaa-D.9CO-2.C833-OlL程-H.知50-3.ODD-D.GB4-2.non-n.6-0000=6&16-0000=MS3LS3F59-l.DOO-D.0130.000L刖85-0110flB9L35-0000=40890,000-D.39flb.noaD.92L乩QilD乩臨6-000-CLSeinmoSJ5.000f除6.DOO-n

44、.Jflfl1.000a.ass2rijjOL-lif-5O.mn02-LD00-d.rwn.noa-n.7(iaTb0iOa.TL-iaoooa-LS.ldS3LDOO0.OFL2.0Q0-G.XTMH-fcrfH：gx，allipc*稠jbi.23T2.521.2.5902.20?Ilfl.QLL90106人6323.025i.1853.2022.劉J2F.33T9L1QT人觸*舊岡I:：.袱:.0TT2.2LQ62104.57&3.355l.呦3.53d2.SSL11-4.QPi93IrSftOlL1J1114Jdith=lEearn94thjivtA2statifiisityxflir

45、ru=3.93Z旳LM匸3.523796-電1Jn|U1?OS3圖4-9c計算結(jié)4.輸出的控制網(wǎng)圖和誤差橢圓圖如下:卜Figure1FileEditViewInsertToolsDesktopWindowHelp勺q包背可口eT測邊網(wǎng)平差程序代碼：globaledddsdddlpnx0y0mlm2m3msppedsidmdgfaabbccrtrrttqlwglobalpathnamenet_names_datafileb_datafilefilename;xO=;yO=;e=;d=;sid=;g=;f=;pn=;%數(shù)據(jù)的讀入if(isempty(pathname)lisempty(net_na

46、me)filename,pathname=uigetfile(*.txt,Tnputfilename);i=find(.=filename);net_name=filename(1:i-1);endfit1=fopen(strcat(pathname,filename),rt);if(fit1=-1)msgbox(InputFileorPathisnotcorrect,Waming,warn);%openafilereturn;endtoreadb_datafile=out.txtfit2=fopen(strcat(pathname,net_name,b_datafile),wt);if(fi

47、t1=-1)msgbox(InputFileorPathisnotcorrect,Warning,warn);return;enda=fscanf(fit1,%d,2);ed=a(1);dd=a(2);sd=ed+dd;fprintf(fit2,%5dn,ed,dd,dd1);fprintf(fit2,n);pn=fscanf(fit1,%d,sd);pn1=pn;pn1,i2=chkdat(sd,pn,pn1);fprintf(fit2,%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5dn,pn);fprintf(fit2,n);a=fscanf(fit1,%f,2*ed);

48、fori=1:edx0(i)=a(2*i-1);y0(i)=a(2*i);fprintf(fit2,%15.3f%15.3fn,x0(i),y0(i);endfprintf(fit2,n);m1=fscanf(fit1,%d,1);fprintf(fit2,%5dn,m1);isid=0;if(m10)fprintf(fit2,siden);a=fscanf(fit1,%f,2);ms=a(1);pp=a(2);fprintf(fit2,%6.2f%6.2fn,ms,pp);a=fscanf(fit1,%d%d%f,3*m1);fori=1:m1e(i)=a(3*i-2);d(i)=a(3*i

49、-1);sid(i)=a(3*i);fprintf(fit2,%5d%5d%15.3fn,e(i),d(i),sid(i);ende,i1=chkdat(sd,pn,e);d,i2=chkdat(sd,pn,d);i3=0;fori=1:m1;if(e(i)=d(i)i3=1;fprintf(fit2,*%5d%5d%5d*n,i,e(i),d(i);endfprintf(fit2,%5d%5d%15.3fn,e(i),d(i),sid(i);endisid=i1+i2+i3;endidir=0;kk=isidif(kk0)msgbox(Errorbyfunctionrddat1,Warnin

50、g,warn);return;end%近似坐標計算xyknow(1:ed)=pn1(1:ed);xyunknow(1:dd)=pn1(ed+1:sd);aa=fscanf(fit1,%d,2);A=aa(1);B=aa(2);mn=fscanf(fit1,%d,1);bb=fscanf(fit1,%d,mn);m=bb;fori=1:length(m)ifi=1ife(m(i)=A|e(m(i)=BP=d(m(i);elseP=e(m(i);endelseife(m(i-1)=A|d(m(i-1)=AB=P;elseA=P;endife(m(i)=A|e(m(i)=BP=d(m(i);else

51、P=e(m(i);endendforj=1:m1if(e(j)=P|d(j)=P)&(d(j)=A|e(j)=A)s1=sid(j);endif(e(j)=P|d(j)=P)&(d(j)=B|e(j)=B)s2=sid(j);endenddeltx=x0(B)-x0(A);delty=y0(B)-y0(A);alfaAB=alfa(deltx,delty);ss=sqrt(deltx*deltx+delty*delty);ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee);x0(P)=x0(A)+ee*cos(alfaAB)-ff*sin(al

52、faAB);y0(P)=y0(A)+ee*sin(alfaAB)+ff*cos(alfaAB);endfprintf(fit2,n);fprintf(fit2,計算的近似坐標n);fori=1:sdfprintf(fit2,%15.3f%15.3fn,x0(i),y0(i);end%k=1;while(k)lo=2062.648062470964;n=2*dd;sum=n*(n+1)/2.0;sd=ed+dd;a=;a(1:m1,1:9)=0.0;c(1:sum)=0.0;w(1:n)=0.0;fprintf(fit2,誤差方程系數(shù)和常數(shù)項n);fori=1:m1%邊長觀測誤差方程dx=x0(

53、d(i)-x0(e(i);dy=y0(d(i)-y0(e(i);ss=sqrt(dx*dx+dy*dy);cosa=dx/ss;sina=dy/ss;a(i,1)=2*e(i)-1-2*ed+1.0e-9;a(i,2)=-cosa;a(i,3)=a(i,1)+1;a(i,4)=-sina;a(i,5)=2*d(i)-1-2*ed+1.0e-9;a(i,6)=cosa;a(i,7)=a(i,5)+1;a(i,8)=sina;a(i,9)=10.0*(ss-sid(i);%ql(i)=(ms2+(ss*pp*0.0001)人2);ql(i)=(ms+ss*pp*0.001);ql(i)=100/q

54、l(i)A2;fprintf(fit2,%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3fn,a(i,:);end%fprintf(fit2,%15.3fn,ql);fori=1:m1%形成法方程forj=1:4jj=fix(a(i,2*j-1);if(jj=0)continue;endw(jj)=w(jj)+a(i,2*j)*a(i,9)*ql(i);di=(jj-1)*(n-jj/2.0);fork=1:4kk=fix(a(i,2*k-1);if(kkkk)continue;endc(di+kk)=c(di+kk)+a(i,2*k)

55、*a(i,2*j)*ql(i);endendendfprintf(fit2,n);fprintf(fit2,%15.3fn,w);fprintf(fit2，法方程系數(shù)n);fprintf(fit2,%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3fn，c)%坐標改正數(shù)和單位權(quán)中誤差的計算sd=ed+dd;n=2*dd;k=0;fori=1:ndxy(i)=0.0;di=(i-1)*(n-i/2.0);forj=1:ndj=(j-1)*(n-j/2.0);if(j1.0)%k=1;%enddxy(i)=-dxy(i)/100.0;end%

56、end%while(k11)while(k11)x(1:ed)=x0(1:ed);y(1:ed)=y0(1:ed);fori=1:ddx(ed+i)=x0(ed+i)+dxy(2*i-1);y(ed+i)=y0(ed+i)+dxy(2*i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);fori=1:sdif(i=ed)vx(i)=0.0;vy(i)=0.0;elsevx(i)=dxy(2*(i-ed)-1);vy(i)=dxy(2*(i-ed);endendvy%l4.4fnfprintf(fit2,坐標改正數(shù)和坐標平差值n);fprintf(fit2,pnvxxyn);fori=l:sdfprintf(fit2,%6d%8.4f%14.4f%8.4f,pn(i),vx(i),x(i),vy(i),y(i);endpvv=0.0;fori=l:np