一種Matlab與C#混合編程實現(xiàn)測量平差程序的方法

林敬娜,程 鋼,路曉明，張啟華,劉玲玲

(1.河南理工大學(xué) 測繪與國土信息工程學(xué)院，河南 焦作 454003;2.河南省自然資源調(diào)查規(guī)劃院，鄭州 450052;3.江蘇省地質(zhì)測繪院，南京 211102)

隨著經(jīng)濟社會發(fā)展和新一輪科技革命，測繪科學(xué)與技術(shù)得到了長足的發(fā)展，測繪的內(nèi)涵及形式也發(fā)生巨大的變化。測量平差作為測繪學(xué)科的基礎(chǔ)理論與方法，是測繪科學(xué)研究及實際測繪工作的基礎(chǔ)和必備工具。測量平差的理論性和綜合性比較強，計算過程中涉及矩陣計算，公式較多，計算量大，僅僅依靠人工計算過程繁雜，簡單易用的計算機軟件平差符合當(dāng)前的用戶需求。

Matlab具有強大的數(shù)值計算性能，在涉及大量數(shù)據(jù)處理分析特別是矩陣計算方面具有其他一些程序語言難以超越的便捷性，在測量平差中進行矩陣計算和線性方程組解算時具有獨特優(yōu)勢。C#具有良好的圖形交互界面和窗體應(yīng)用程序開發(fā)功能，利用其與Matlab混合編程，能夠快速實現(xiàn)測量數(shù)據(jù)平差的相關(guān)功能。

潘雄、白征東等人闡述了Matlab在測量平差中的優(yōu)勢，并以間接平差為例介紹了實現(xiàn)過程[1-2]。杜澤明等人給出了C#與Matlab混合進行平差程序設(shè)計的基本框架，但未詳述具體平差方法的實現(xiàn)過程[3]。朱杰等人基于C#提出矩陣類的設(shè)計并探討其在平差程序設(shè)計中的應(yīng)用[4]。劉淼等人以間接平差為例探討基于Matlab Web Server實現(xiàn)平差計算的方法，該程序只能在網(wǎng)絡(luò)支持條件下使用[5]。文中在綜合以上研究成果優(yōu)點的基礎(chǔ)上，基于Matlab與C#設(shè)計實現(xiàn)一套計算簡單方便的測量平差系統(tǒng)，以條件平差和間接平差為例，對實現(xiàn)過程的關(guān)鍵方法進行探討，為相關(guān)工作提供參考。

1 原理介紹

1.1 Matlab與C#簡介

1.1.1 Matlab簡介

Matlab語言是一種功能十分強大的計算機語言，數(shù)值計算性能非常高，是一種可以應(yīng)用于數(shù)值分析、程序開發(fā)以及數(shù)據(jù)可視化的高級技術(shù)計算語言和環(huán)境。Matlab以矩陣計算作為基礎(chǔ)，擁有完備的函數(shù)集，將數(shù)值計算、程序設(shè)計與數(shù)據(jù)可視化集于一體，可以進行數(shù)學(xué)計算、符號運算、可視化建模、圖形處理等等。無論是在科學(xué)研究還是工程領(lǐng)域中都具有功能豐富的應(yīng)用工具箱[6]，應(yīng)用范圍廣泛，能夠幫助用戶解決復(fù)雜的數(shù)學(xué)計算與分析問題，大大縮減用戶編程計算的工作量。

1.1.2 C#簡介

C#是一種C和C++衍生出來的基于.NET框架的面向?qū)ο蟮木幊陶Z言，由微軟公司發(fā)布。它擁有非常多的類庫，使用方便，利于學(xué)習(xí)，并且具有非常友好的用戶界面，可用于Windows桌面應(yīng)用程序、Web應(yīng)用程序、Web服務(wù)程序、Windows服務(wù)程序等程序開發(fā)。C#推出后憑借其強大的操作能力、面向組件、面向?qū)ο蟆⒗陂_發(fā)的特點受到廣大程序員的喜愛與支持。

1.1.3 Matlab與C#混合編程

Matlab與C#混合編程，充分利用Matlab強大的計算能力和C#中Winform直觀的顯示效果，進行優(yōu)勢互補?；旌暇幊虒崿F(xiàn)方法多種多樣，文中采用C#調(diào)用Matlab編寫的動態(tài)鏈接dll文件的方法，首先將平差處理函數(shù)在Matlab中實現(xiàn)并封裝為dll文件，然后在VS中建立C#窗體應(yīng)用程序進行dll文件的引用，實現(xiàn)Matlab與C#的混合編程。

1.2 測量平差原理

在測量過程中，由于測量儀器精度、人為因素和外界條件的一些影響，存在不可避免的測量誤差，測量平差的目的即在于通過平差方法消除觀測誤差產(chǎn)生的矛盾，計算出觀測值最可靠的結(jié)果，并且進行測量成果的精度評定。

測量學(xué)中涉及多種類型控制網(wǎng)的平差，有導(dǎo)線網(wǎng)平差、水準網(wǎng)平差、GPS網(wǎng)平差等，測量平差方法也有多種：間接平差，條件平差，附有限制條件的間接平差，附有未知參數(shù)的條件平差等。每一種平差方法都有其特定的平差模型和精度評定的方法。在進行測量平差計算時，要根據(jù)不同計算模型的特點建立函數(shù)組，并選擇相應(yīng)的計算方法，可以生成函數(shù)庫，通過函數(shù)調(diào)用完成測量數(shù)據(jù)分析和平差的處理工作。下文以間接平差和條件平差在水準網(wǎng)平差中的應(yīng)用為例，介紹開發(fā)過程。

1.2.1 間接平差

間接平差法又叫參數(shù)平差法，該方法選定合適數(shù)量并且與觀測值存在一定數(shù)學(xué)關(guān)系的獨立未知量作為參數(shù)，函數(shù)模型是將每個觀測值都分別表達成選定參數(shù)的函數(shù)[7]。間接平差數(shù)學(xué)模型比較簡單，便于評定平差值及其函數(shù)的精度。模型求解方法是按照最小二乘原理的方法，求自由極值，解出參數(shù)的最或然值，從而求得各觀測值的平差值。

(1)

式中：B表示誤差方程的系數(shù)陣；l表示常數(shù)項；V表示觀測值的改正數(shù)。

Nbb=BTPB，W=BTPl.

(2)

(3)

解求出平差值之后，還要分析平差值的中誤差，即為精度評定：

(4)

1.2.2 條件平差

在測量工作中，為了能及時發(fā)現(xiàn)錯誤，盡可能地減少測量誤差，提高測量精度和效率，測量人員往往會進行必要的多余觀測。一個幾何模型中條件方程個數(shù)由多余觀測數(shù)決定，有幾個多余觀測，就會產(chǎn)生幾個條件方程，條件平差法就是以條件方程為函數(shù)模型的平差方法。因此，條件平差以觀測量間的函數(shù)關(guān)系為條件構(gòu)造條件方程，不含獨立參數(shù)。這些條件方程間需要獨立，保證列出的方程組線性無關(guān)[8]。

設(shè)平差問題中有n個觀測值L，已知觀測值協(xié)因數(shù)陣Q=P-1，必要觀測數(shù)為t。條件平差的方程模型為：

AV-W=0.

(5)

式中：A表示條件方程的系數(shù)陣；W表示常數(shù)項；V表示觀測值的改正數(shù)。

實際上，條件方程求解的關(guān)鍵是改正數(shù)V的求解。根據(jù)最小二乘原理，V必須滿足VTPV=min的要求，用拉格朗日求極值的方法可得到觀測值改正數(shù)：

V=P-1ATK.

(6)

其中

則觀測值的平差值

(7)

就可以得到各個觀測量的平差值。然后進行精度評定，分析平差值函數(shù)的中誤差：

(8)

2 測量平差系統(tǒng)實現(xiàn)

2.1 Matlab函數(shù)設(shè)計

每種測量平差方法都有其自身的平差模型與解算步驟，文中將每種測量平差方法分別編寫成為Matlab中的function函數(shù)，保證方法的獨立性，便于后續(xù)調(diào)用。下文以水準平差為例，編寫相關(guān)函數(shù)。

2.1.1 間接平差

間接平差的函數(shù)形式為function [X1,L1,D1]=jianjieAdjustment(a1,b1,c1)。其中function是函數(shù)聲明，jianjieAdjustment是函數(shù)名稱。傳入?yún)?shù)a1，b1，c1為3個矩陣：a1矩陣行數(shù)為水準網(wǎng)測段數(shù)(觀測值總數(shù))，列數(shù)為4，分別表示每測段路線的起點編號、終點編號、高差和水準路線長度；b1矩陣行數(shù)為已知點個數(shù)，列數(shù)為2，分別表示已知高程點編號和高程數(shù)值；c1為1*4矩陣，分別表示水準網(wǎng)測段數(shù)、水準點個數(shù)、未知高程點個數(shù)和已知高程點個數(shù)。輸出參數(shù)X1，L1，D1為3個矩陣，分別表示高程平差值、高差平差值和高程平差值的協(xié)方差陣。注意水準網(wǎng)中的水準點點號按照先已知水準點，后未知水準點的規(guī)則依次進行編號。

函數(shù)設(shè)計主要分為三部分：

第一，誤差方程系數(shù)矩陣和常數(shù)項矩陣的獲取。系數(shù)矩陣B為n×t矩陣，常數(shù)項矩陣l為n×1矩陣，n表示測段個數(shù)，t表示未知點個數(shù)，將每測段的起點和終點高程是否已知作為判定規(guī)則，逐測段來計算B和l中的每個元素B(i,j)與l(i,j)。過程代碼如下：①當(dāng)?shù)趇測段起點和終點高程均未知時，B(i,qi-y)=-1;B(i,zh-y)=1；②當(dāng)?shù)趇測段起點高程已知，終點高程未知時，B(i,zh-y)=1，且可計算終點近似高程X0(zh,1)=X0(qi,1)+L(i,1)；③當(dāng)?shù)趇測站起點高程未知，終點高程已知時B(i,qi-y)=-1，且可計算起點近似高程X0(qi,1)=X0(zh,1)-L(i,1)；B矩陣中其余元素均為0。常數(shù)項矩陣元素l(i,1)=L(i,1)-X0(point2,1)+X0(point1,1)。其中qi表示起點編號，zh表示終點編號，y表示已知高程點個數(shù)，L表示觀測值。

第二，誤差方程的求解。借助Matlab的矩陣計算方法，設(shè)置輔助量Nbb=B′*P*B，W=B′*P*l；計算未知參數(shù)x=inv(Nbb)*W，高差改正數(shù)V=B*x-l，高差平差值L1=L+V，高程平差值X1=X0((y+1):end,1)+x。

第三，精度評定，即參數(shù)平差值協(xié)方差的計算。首先計算單位權(quán)中誤差：d0=sqrt ((V′*P*V)/(n-t))，n-t表示多余觀測數(shù)。其次，求取參數(shù)平差值協(xié)方差陣D1=d0*d0 *Q1，其中Q1=inv(Nbb)表示參數(shù)平差值的協(xié)因數(shù)陣。最后求取參數(shù)及觀測值的中誤差。

2.1.2 條件平差

條件平差的函數(shù)形式為function [L1,H0,D1]=tiaojianAdjustment(a1,b1,c1)。其中傳入?yún)?shù)與間接平差一致，輸出參數(shù)L1，H0，D1為3個矩陣，分別表示高差平差值、高程平差值和高差平差值協(xié)方差陣。水準點點號編號規(guī)則不變。

條件平差的函數(shù)設(shè)計分為誤差方程系數(shù)矩陣和常數(shù)項矩陣的獲取、條件方程的求解、未知點高程計算、精度評定4部分。

條件方程的系數(shù)矩陣A和常數(shù)項矩陣W運用文獻[9]的方法，通過引入一個傳遞矩陣t來獲取，t為N*n的矩陣，n表示測段個數(shù)，N表示水準點個數(shù)。根據(jù)參與計算的測段路線的起點和終點高程是否已知，未知水準點的近似高程是否已經(jīng)計算得到為判定規(guī)則，計算傳遞矩陣t中的元素t(i,j)[9]。傳遞矩陣t計算完成后，再根據(jù)傳遞矩陣t計算系數(shù)矩陣A和常數(shù)項矩陣W，A為r*n矩陣，W為r*1矩陣，n表示測段個數(shù)，r表示多余觀測數(shù)，若某測段路線高差未用于計算傳遞矩陣t，則第m個條件方程的系數(shù)A(m,1:n)=t(qi,1:n)-t(zh,1:n);A(m,k)=1;W(m,1)=0(j,1)-H0(i,1)-L(k,1)，H0和L分別表示高程和高差觀測值。

條件方程的求解根據(jù)條件平差公式計算，改正數(shù)V=inv(P)*A′*K，其中聯(lián)系數(shù)K=inv (Naa)*W，Naa=A*inv(P)*A′，高差平差值L1=L+V。

未知高程點的高程計算根據(jù)已知高程點的高程與條件方程求解得到的高差平差值計算，其值由測段起點高程與測段高差之和，或測段終點高程與測段高差之差依次求取，直至所有水準點高程計算完成。

精度評定實現(xiàn)了觀測值平差值的精度計算。單位權(quán)中誤差：d0=sqrt((V′*P*V)/r)，r表示多余觀測數(shù)。觀測值平差值協(xié)方差陣Ql=Q-Q*A′*inv(Naa)*A*Q，其中Q表示觀測值的協(xié)因數(shù)陣。觀測值平差值協(xié)方差陣D1=d0*d0*Q1，D1矩陣各對角線即為高差平差值的中誤差。

2.2 C#系統(tǒng)實現(xiàn)

Matlab與C#混合編程利用了Matlab強大的計算能力和C#友好的用戶界面的優(yōu)勢，制作一個測量平差系統(tǒng)平臺，能夠更加方便地進行平差計算，進行平差網(wǎng)的直觀圖形顯示，并且可以脫離Matlab與VS的系統(tǒng)軟件，也不需要掌握Matlab與C#的軟件使用和程序語法規(guī)則，而且此平臺操作簡單，可完成不同測量平差方法進行平差的功能。

2.2.1 系統(tǒng)界面搭建

C#的窗體應(yīng)用程序可以通過工具箱中的控件搭建用戶界面，直接拖拽到窗體中就可完成測量平差平臺所需界面的設(shè)計，并為控件按鈕添加事件相應(yīng)程序。測量平差系統(tǒng)平臺搭建界面如圖1所示，其中文件菜單提供了數(shù)據(jù)輸入、保存和輸出等功能。測量數(shù)據(jù)輸入文件存儲測段起點、終點及其高差和水準路線長度等觀測數(shù)據(jù)；高程數(shù)據(jù)輸入文件存儲已知高程點號和高程等數(shù)據(jù)；點號數(shù)據(jù)輸入文件存儲測段個數(shù)、水準點個數(shù)、未知點個數(shù)和已知點個數(shù)等數(shù)據(jù)。平差網(wǎng)圖形輸入文件可以導(dǎo)入測量平差網(wǎng)的拓撲結(jié)構(gòu)，以供實現(xiàn)圖形可視化，方便計算者直觀地進行分析。平差方法菜單提供5種平差方法可供選擇，分別實現(xiàn)不同平差方法的平差計算。

圖1 系統(tǒng)界面

2.2.2 動態(tài)鏈接庫(DLL)生成

動態(tài)鏈接庫(DLL)是一種可實現(xiàn)代碼共享和重用的基于Windows的程序模塊，包含可執(zhí)行代碼、數(shù)據(jù)、資源等[10]。Matlab有自帶編譯器，可創(chuàng)建function函數(shù)的動態(tài)鏈接庫，即.dll文件。在Matlab Compiler窗口，選擇Library Compiler，在設(shè)置窗口中選擇.NET Assembly，加載編寫好的測量平差的M文件，設(shè)置類名和路徑，即可利用Package進行打包，完成動態(tài)鏈接庫創(chuàng)建。

2.2.3 C#對DLL文件的調(diào)用

1)項目準備。實現(xiàn)C#與Matlab的混合編程，需要在C#項目中添加上述生成的DLL文件以及MWArray.dll環(huán)境文件，然后才可以實現(xiàn)Matlab函數(shù)方法調(diào)用及數(shù)據(jù)參數(shù)轉(zhuǎn)換。引用添加成功后，需要為項目代碼添加必要的命名空間，例如數(shù)據(jù)轉(zhuǎn)換需要：using MathWorks.Matlab.NET.Arrays。

2)控件代碼編寫。由于測量平差需要輸入的數(shù)據(jù)繁多，規(guī)定以txt文本文檔的形式輸入。控制程序通過輸入數(shù)據(jù)文件獲取相關(guān)的參數(shù)，并存入數(shù)組之中以供計算使用。

C#中調(diào)用Matlab函數(shù)的關(guān)鍵是掌握數(shù)據(jù)類型一致的矩陣與數(shù)組之間的轉(zhuǎn)換，將C#中的參數(shù)輸入到Matlab函數(shù)時，要將參數(shù)轉(zhuǎn)化為Matlab的參數(shù)形式，通常是MWArray類型。例如條件平差參數(shù)a1數(shù)據(jù)類型轉(zhuǎn)換：MWNumericArraya1=new MWNumericArray(p1)，其中MWNumericArray是Matlab與C#的中間數(shù)據(jù)類型。Matlab函數(shù)返回的參數(shù)，也要轉(zhuǎn)化為C#用的類型，比如數(shù)組或者數(shù)值類型。C#中調(diào)用條件平差函數(shù)封裝成的類方法的調(diào)用代碼如下，以agrsIn表示輸入的矩陣數(shù)組，agrsOut表示返回的矩陣數(shù)組：

MWArray[] agrsIn = new MWArray[] { a1, b1,c1 };

MWArray[] agrsOut = new MWArray[3];

TiaojianAdjustment.Class1 level = new TiaojianAdjustment.Class1();

level.tiaojianAdjustment(3, ref agrsOut, agrsIn);

平差計算完成后，可借助文本框完成結(jié)果輸出。

2.2.4 可執(zhí)行文件生成

將C#中的測量平差數(shù)據(jù)處理窗體程序生成可執(zhí)行的exe文件，便于程序分發(fā)和共享。將C#窗體程序中引用的DLL文件放入項目文件夾下，利用VS自帶的打包項目功能或者第三方打包軟件將此項目生成exe文件即可，程序可不依賴于Matlab和VS軟件環(huán)境運行。

2.3 實例驗證

文中以水準平差為例說明平差程序的應(yīng)用，在某校園內(nèi)布置水準網(wǎng)如圖2所示，利用自動安平水準儀(KL80)進行四等水準觀測，水準網(wǎng)中水準路線閉合差、前后視距差、紅黑面讀數(shù)差、紅黑面所測高差之差等均符合國家四等水準測量規(guī)范。A和D兩點表示已知水準點，HA=94.957 m，HD=94.919 m，B，C，E，F(xiàn)4點表示待定水準點，其余的觀測高差和相應(yīng)的水準路線長度見表1，求解待定點的高程平差值及中誤差。

圖2 水準網(wǎng)路線圖

表1 實驗相關(guān)數(shù)據(jù) m

將此水準網(wǎng)中的水準點從已知到未知依次編號，A,D分別為1、2號點，B,C,E,F分別為3、4、5、6號點，整理觀測數(shù)據(jù)，在系統(tǒng)中輸入整理好的已知數(shù)據(jù)文件，分別進行條件平差(見圖3)和間接平差(見圖4)。

將其計算結(jié)果保留三位小數(shù)如表2所示。

由表2結(jié)果可見，在此水準網(wǎng)中，待定點B的高程中誤差最小，精度最高，E的高程中誤差最大，精度最低。并且此例采用兩種平差方法計算的結(jié)果相同，印證了同一個平差問題無論是采用間接平差還是條件平差，其結(jié)果是一致的，因此可以驗證文中平差算法具有可行性。利用Matlab與C#混合編程實現(xiàn)的測量平差數(shù)據(jù)處理系統(tǒng)軟件使用方便，計算準確，解決了人工計算所帶來的麻煩。

圖3 條件平差計算結(jié)果

圖4 間接平差計算結(jié)果

表2 計算結(jié)果

3 結(jié)束語

文中以水準網(wǎng)條件平差和間接平差為例，說明Matlab與C#混合編程實現(xiàn)測量平差程序的方法，體現(xiàn)了Matlab具有強大數(shù)據(jù)計算以及C#具有良好的用戶界面開發(fā)功能的優(yōu)勢，探究和驗證了Matlab與C#混合編程在測量平差數(shù)據(jù)處理中的實用性。文中設(shè)計了一套測量平差簡單計算的應(yīng)用系統(tǒng)，使得用戶在不需要掌握Matlab與VS軟件和程序語言專業(yè)語法命令的情況下即可使用該程序進行平差，操作簡單并且易用，實例證明了方法的可行性。