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

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

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

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

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

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

1 原理介紹

1.1 Matlab與C#簡介

1.1.1 Matlab簡介

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

1.1.2 C#簡介

C#是一種C和C++衍生出來的基于.NET框架的面向對象的編程語言，由微軟公司發布。它擁有非常多的類庫，使用方便，利于學習，并且具有非常友好的用戶界面，可用于Windows桌面應用程序、Web應用程序、Web服務程序、Windows服務程序等程序開發。C#推出后憑借其強大的操作能力、面向組件、面向對象、利于開發的特點受到廣大程序員的喜愛與支持。

1.1.3 Matlab與C#混合編程

Matlab與C#混合編程，充分利用Matlab強大的計算能力和C#中Winform直觀的顯示效果，進行優勢互補。混合編程實現方法多種多樣，文中采用C#調用Matlab編寫的動態鏈接dll文件的方法，首先將平差處理函數在Matlab中實現并封裝為dll文件，然后在VS中建立C#窗體應用程序進行dll文件的引用，實現Matlab與C#的混合編程。

1.2 測量平差原理

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

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

1.2.1 間接平差

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

(1)

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

Nbb=BTPB，W=BTPl.

(2)

(3)

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

(4)

1.2.2 條件平差

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

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

AV-W=0.

(5)

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

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

V=P-1ATK.

(6)

其中

則觀測值的平差值

(7)

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

(8)

2 測量平差系統實現

2.1 Matlab函數設計

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

2.1.1 間接平差

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

函數設計主要分為三部分：

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

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

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

2.1.2 條件平差

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

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

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

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

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

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

2.2 C#系統實現

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

2.2.1 系統界面搭建

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

圖1 系統界面

2.2.2 動態鏈接庫(DLL)生成

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

2.2.3 C#對DLL文件的調用

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

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

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

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

MWArray[] agrsOut = new MWArray[3];

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

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

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

2.2.4 可執行文件生成

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

2.3 實例驗證

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

圖2 水準網路線圖

表1 實驗相關數據 m

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

將其計算結果保留三位小數如表2所示。

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

圖3 條件平差計算結果

圖4 間接平差計算結果

表2 計算結果

3 結束語

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