EXCEL内置函数在测量平差中的应用.doc

EXCEL内置函数在测量平差中的应用
　　　　　　　　　　　　　　　　　　　　　　　　王海英 陈世梅
　　　　　　　　　　　　　　　　　　　（南宁市勘测院 广西南宁 530021）
【摘　要】　测量平差是测量数据处理的重要环节，然而平差计算过程中，线性方程解算的问题却比较麻烦。本文利用Excel的内置函数简化了线性方程解算的问题，并结合具体的实例，讲述利用Excel内置函数进行测量平差的方法。
【关键词】　矩阵 Excel 测量平差

　　0 引言
　　在现代测量中，对通过各种测量方法所采集得到的原始数据，往往需要根据误差理论的方法对存在的各类误差进行平差处理从而取得最或然结果。而对于偶然误差的处理，利用最小二乘法准则进行平差计算的过程，通常都要对误差方程式或条件方程式进行整合处理求得法方程组，然后解算法方程、计算改正数、精度评定等等，一系列计算步骤的进行，如果没有现成的专业软件，而用手工的办法去完成平差计算工作，将是一件很繁琐的事情。
　　其实，测量平差的过程简单的来说就是解算线性方程组的过程，而最令人烦恼的部分是线性方程的求解。在处理一些小的工程项目平差计算的时候，如果手头没有像MATLAB等具有矩阵运算功能的应用软件，利用Microsoft Office Excel的内置函数，同样的可以帮助计算者比较轻松的完成计算任务。
　　以下要阐述的就是利用Excel转置粘贴功能以及矩阵计算的函数MMULT（矩阵乘）、MINVERSE（矩阵求逆），实现测量平差之线性方程组解算的过程。
　　1 平差的数学模型
　　为了加快平差解算的作业效率，应该根据实际情况选择适当的数学模型。实际工作中，有两种数学模型得到了较为广泛的应用，即间接平差模型和条件平差模型；它们的应用公式[1]分别如下：
　　1.1 间接平差
　　误差方程：V = B δx+ J ；J = B X0+ B0—L ；
　　法方程：BT P B δx+ BT P J = 0 ；
　　改正数计算：δx = —（ BT P B） -1 BT P J ；
　　PVV计算：V-1 P V = JT P V + δxT BT P V ；
　　权函数式：δΦ = F δx ；
　　平差值的权倒数计算： 1 / PΦ= FT （BT P B）-1 F ；
　　公式中各种符号的含义：观测向量 L；相应的权阵 P；误差方程的系数矩阵B；未知数的改正数δx；常数阵 J；常数向量B0；未知数的近似值X0；权函数的系数阵F。
　　1.2 条件平差
　　条件方程式：A V + W= 0 ；W = A L + A0；
　　法方程：A P-1ATK + W = 0 ；
　　改正数计算：V = P-1 AT K ；
　　PVV计算：V-1 P V = —WT K ；
　　权函数式：δΦ = F δi；
　　平差值的权倒数计算：1 / PΦ= FTP-1 F — （A P-1 F）T （A P-1AT ）-1 A P-1 F ；
　　公式中各种符号的含义：观测向量 L；相应的权阵 P；条件方程的系数矩阵A；条件改正数 V；条件闭合差 W；法方程联系数矩阵 K；权函数的系数阵F。
　　2 计算方法
　　不难看出，以上平差计算的数学模型，全都是以矩阵的形式给定的。如果掌握了矩阵计算的方法，所有的计算将迎刃而解。
　　由于测量计算作业过程中采用不同的数学模型，以及图形条件的多样化，对于误差方程或者条件方程的产生，在本文所述的解算方案当中，还不能自动完成；当然，可以利用Excel中的VBA功能编程处理，但将可能给计算者增加了难度。在此仅就手工生成误差方程组或者条件方程组之后，如何使用Excel内置函数进行矩阵计算的方案进行阐述。
　　2.1 数据的输入
　　将系数矩阵的数据填入电子表格的相应区域中。
　　2.2 矩阵的转置
　　（1）选中电子表格中系数矩阵的数据区域，点击右键，选击‘复制’菜单选项。
　　（2）在电子表格中点击系数矩阵数据区域以外的任意单元格，点击右键，选击“选择性粘帖”菜单项。
　　（3）点击“选择性粘帖”菜单项，出现“选择性粘帖”对话框；该对话框有三部分选项，第一部分“粘帖”，选中“全部”选项，第二部分“运算”，选中“无”选项，第三部分选中“转置”选项，单击确定，即可完成系数矩阵数据的转置操作。
　　2.3 矩阵乘运算
　　（1）选中电子表格中的空白区域，该区域的行数、列数等于矩阵相乘所得矩阵的行数、列数，点击工具栏中的fx（函数）工具按钮。
　　（2）弹出插入函数的对话框，在该对话框的左边的函数分类中选择“数学与三角函数”，在右边的函数名中选中“MMULT”，单击确定。
　　（3）弹出函数参数输入的对话框，该对话框提示输入两组参数，第一个参数Array1栏中输入第4步骤中转置矩阵区域的行列号，在Array2栏中输入系数矩阵区域的行列号——可以不用键盘输入，而用Array*内右端的按钮，回到表格视图中用鼠标涂选表格区域，按一下回车键即可；在表格视图的公式栏应该有“=MMULT（‘转置矩阵区域’：‘系数矩阵区域’）”的描述。
　　（4）同时按住Shift和Ctrl两键，按下回车键，这时，第1步骤选定的区域内的单元格所显示的结果即为联系方程的系数阵。
　　2.4 矩阵求逆运算
　　（1）选中电子表格中的空白区域，其行列数与系数矩阵相同，点击工具栏中的fx工具按钮。
　　（2）在弹出的对话框中左边的函数分类选取“数学与三角函数”，右边选择函数名“MINVERSE”，单击确定。
　　（3）在弹出的参数选择对话框中Array1中输入系数矩阵，在表格视图的公式栏中也应该有“=MINVERSE （系数矩阵区域）”的描述。
　　（4）同时按住Shift和Ctrl两键，按下回车键，第1步骤选定区域内的单元格所显示的结果即为系数矩阵的逆矩阵。
熟练掌握上述矩阵计算的步骤、方法，解决测量平差之线性方程组解算的问题就会变得简单、容易了。
　　3 注意事项
　　在计算过程中，受Excel软件本身的约束和限制，可能会带来不可预想的结果，因此需要注意一下几个事项：
　　（1）受电子表格列数的限制，Excel最多可以计算有256列的矩阵。若要求解有更多列的矩阵，可以利用分块矩阵的办法，或者VBA 语言编程进行解决。
　　（2）改变计算结果的精度，可作以下操作：‘格式’菜单→单元格→数字→数值→小数位数，选定相应的数值即可。若需用双精度计算，设定小数位为15位，即可得到双精度的求解结果。
　　（3）有关矩阵的计算，还有其他的运算功能，如相加、相减等，均可以在Excel编辑功能的‘复制’与‘选择性粘贴’中实现。
　　4 实例应用
　　根据以上提供的计算步骤，本文以“后方交会、附加测边条件”的一个实际应用为例，进一步说明Excel内置函数在测量平差中的应用。其中涉及矩阵运算的环节均利用Excel的内置函数求解（计算过程中边长观测以及角度观测值的权系数均取1；单位为米、秒）。
　　 4.1 已知数据及观测成果（待定点：P）
测站
观测方向
观测点纵坐标
观测点横坐标
观测方向角
观测边长

P
仙鹤顶
2431773.087
491172.970
　00　00　00

龙珠塔
2427738.078
494836.670
　30　35　02

南剑
2426330.488
493695.187
　41　37　50

YN256
2430554.354
485826.957
　175　31　47
147.405
　　4.2 近似值计算
　　取待定点P的近似坐标为：X=2430597.925；Y=485967.775；并解算得以下近似的边长、方位角：
待定点
已知点
边长
方位角
P
仙鹤顶
5336.203
77　16　40.03
龙珠塔
9318.585
107　52　20.41
南剑
8827.452
118　54　34.29
YN256
147.405
252　48　26.00
　　4.3 方向、边长误差方程系数及常数项
　　由公式aij = Δy*ρ″/Sij2 ；bij = -Δx*ρ″/Sij2可计算出Δx项和Δy项的系数。
　　对S = SQR（Δx 2+Δy 2）求偏导并推导出边长条件的线性方程，求得边长条件的系数。
　　由观测值和近似方位角，解算得误差方程的常数项J = B X0+ B0—L.
　　根据误差方程的数学模型 V = B δx + J ;其系数矩阵和常数阵如下：

待定点
已知点
Δx 项系数
Δy 项系数
常数项
V01
P
仙鹤顶
37.7
－8.5
0.0
V02
龙珠塔
21.1
6.8
38.4
V03
南剑
20.5
11.3
4.3
V04
YN256
－1336.8
413.6
－1.0
Vs
（边长条件）
YN256
0.3
1.0
0.0
　　4.4 改化后的角度、边长误差方程系数及常数项

待定点
（自）方向
（到）方向
δx项系数
δy项系数
常数项
改正数
V1
P
仙鹤顶
龙珠塔
－16.6
15.3
38.4
0
V2
仙鹤顶
南剑
－17.2
19.8
4.3
0
V3
仙鹤顶
YN256
－1374.5
422.1
－1.0
0
Vs
（边长）
P
YN256
0.3
1.0
0.0
0.0
（*从此开始，把方程组的系数阵、常数阵的系数键入EXCEL的相应单元格中，并进行相应计算。）
　　4.5 矩阵转置运算
　　（1）将系数阵的数据填入电子表格的（A1：B3）区域中；将常数项填入（D1：D3）区域中；如下图：

　　（2）选取（A1：B3）区域，做‘编辑’→‘复制’的菜单操作；
　　（3）鼠标点击一下A6单元格；
　　（4）做‘编辑’菜单下的‘选择性粘贴’操作
　　（5）在‘选择性粘贴’的对话框中，选择：粘贴的‘全部’，运算的‘无’，以及‘转置’共三项；然后确定，便完成了矩阵的转置。如下图：

　　4.6 法方程组生成（矩阵乘运算）
　　（1）选中（A9：B10）区域；点击插入函数的工具按钮fx；
　　（2）弹出插入函数的对话框，在对话框的左边的函数分类中选择“数学与三角函数”，在右边的函数名中选中“MMULT”，单击确定。
　　（3）根据BT P B，输入函数的参数：在Array1栏中输入A6：C7；在Array2栏中输入A1：B3；此时表格视图的公式栏中应有MMULT（A6:C7，A1:B3）的描述。同时按住SHIFT和CTRL键不放，按下回车键，这时选中的（A9：B10）区域便有了结果。如下图

　　（4）、计算BT P J时使用类似的方法，即：选中（A17：A18）区域，在Array1栏中输入A6：C7；在Array2栏中输入D1：D3；此时表格视图的公式栏中应有MMULT（A6:C7，D1:D3）的描述。同时按住SHIFT和CTRL键不放，按下回车键，这时选中的（A17：A18）区域便有了结果。如下图

　　4.7 法方程组解算（矩阵求逆）
　　顾及全站仪的观测精度以及现场的情况，拟将边长观测值作为附加约束条件，进行有附加条件的间接平差。则有附加条件间接平差的联系方程的系数为：
δx
δy
K
L
1889821.65
－580770.99
0.3
663.10
－580770.99
178794.54
1.0
250.56
0.3
1.0
0.0
0.0
　　那么，须作：
　　（1）在单元格A19中填入0.0
　　（2）在（C9：C11）；（A11：C11）两区域键入K的系数。如下图

　　（3）选中（A13：C15）区域，点击插入函数的工具按钮fx；在弹出的对话框中左边函数分类选择“数学与三角函数”，右边选择函数名“MINVERSE”，单击确定。
　　输入函数的参数：在Array1栏中输入A9：C11；此时表格视图栏中应有MMULT（A9:C11）的描述。同时按住SHIFT和CTRL键不放，按下回车键，这时选中的（A9：A10）区域便得到计算的结果。如下图

　　这就是联系（法）方程组系数阵的逆阵（即权逆阵）。
　　4.8 方程的解
　　根据公式：δx = —（ BT P B） -1（BT P J）进行计算。
　　选中（A21：A23）区域；在Array1栏中输入A13：C15；在Array2栏中输入A17：A19；此时表格视图的公式栏中应有MMULT（A13：C15，A17:A19）的描述。同时按住SHIFT和CTRL键不放，按下回车键，这时选中的（A21：A23）区域即得到结果。如下图

　　那么方程的解就是
　　δx = 0.0003；δy = － 0.0001；K = －416.0114；
　　4.9 待定点P的坐标
　　由方程的解计算待定点P改正后的坐标值为：
　　XP = 2430597.925 ＋ 0.0003 = 2430597.925；
　　YP =485967.775 － 0.0001 =485967.775；
　　4.10 精度评定
　　Mo = SQR（（38.42＋4.32）/2） = ±27.3
　　待定点P的点位中误差为：
　　Ms = ±27.3 * SQR（0.000000444＋0.000000040） = ±0.019米。