一种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软件和程序语言专业语法命令的情况下即可使用该程序进行平差，操作简单并且易用，实例证明了方法的可行性。