探析基于Excel实现工程测量中的坐标转换

罗孝楠 周广勇 谢全远

摘要：因各个国家卫星导航使用的坐标系不一致，GPS工程测量中坐标转换是实际应用的重要步骤，针对目前坐标转换软件数据录入不便、过程不可控等问题，文中探析四参数法坐标转换参数求解方法，并在Excel下利用少量内部函数，设计实现了四参数求解和坐标转换电子表格。通过工程实例坐标进行验证，求解结果准确，计算过程数据可视，编辑录入数据方便，方便在测量工程、教学实验中使用。

关键词：卫星导航；四参数；坐标转换；Excel

中图分类号：TP311.1    文献标识码：A

文章编号：1009-3044（2023）35-0118-03

开放科学（资源服务）标识码（OSID）

0 引言

伴随着卫星技术深入发展，在实践中得到了广泛应用，特别是美国全球定位系统（GPS） 的成功建立和应用，我国北斗卫星导航测量技术也取得了前所未有的进步和发展。卫星定位测量技术具有精准、快速、高效和全天候等显著的特点，使得测绘行业经历了一场深刻的变革。丰富的应用场景已经被应用到了现代社会经济中的很多领域，尤其是在军事、交通、测量等方面和领域取得了显著的成效，凸显了重要的应用价值。在因此全球导航定位卫星系统，在全球范围内的众多行业中，尤其是在测绘行业有着越来越广泛的应用。

全球定位系统（GPS） 所采用的坐标系是WGS-84坐标系，因此直接采集的点位坐标属于WGS-84坐标系，而目前我们国家工程测量成果普遍使用的是以CGCS2000坐标系或是地方（任意）独立坐标系为基础的坐标数据。因此必须将WGS-84坐标转换到CSCS2000坐标系或地方（任意）独立坐标系[1]。坐标转换是从一种坐标系统映射到另一种坐标系统的过程。通过建立两个坐标系统之间一一对应关系来实现。一般四参数应用于平面转换和七参数应用于三维转换。

在许多测量数据处理软件中，如徕卡LGO、天宝TBC、南方SGO、Pinncle等数据处理软件，都有坐标系转换模块，有些功能比较齐全，如在TBC软件中包含了七参数法、格网法、多元回归法；LGO软件中有格网法、七参数法、三参数法、格网与参法结合法，有三维转换也有二维转换。在实际应用中，可以结合测区大小、范围内已知点的数量与分布情况决定采用哪一种方法。这些软件在学习、教学过程中还是存在一些问题，软件比较复杂庞大、学习成本比较高，计算过程不可见，不利于对原理的理解。

Excel是微软公司推出的一个广泛应用于商业和科学领域的电子表格程序。它具有一些强大的功能和特点，使其成为各种应用程序中不可缺少的工具。简洁的界面、优秀的计算功能和图表工具，再加上出色的市场营销，使Excel成为流行的个人计算机电子表格软件。它的主要功能有制作表格、数字格式的转换、制作图表、函数功能、数据处理功能等[2]。强大的数据编辑功能，使得利用Excel编制坐标转换程序表格方便可行。

1 四参数坐标转换原理

四参数是用于两个平面直角坐标系之间的互相转换，而七参数是用于两个三维空间直角坐标系之间的转换。四参数可以由两个及以上具有平面二维坐标的已知等级控制点求出，高程可以用拟合方法求出，求解较为简单，也较容易理解；七参数转换是空间坐标系的转换需要在测区布设一定密度的控制网点，需要三个以上已知三维坐标的控制点，利用整个网的WGS-84坐标系下的三维约束平差结果和当地坐标系统的二维约束平差结果及各点的高程解算，参数求解比较烦琐，理解起来相对困难。

1.1 平面坐标转换

点位的平面坐标可以用高斯投影取得。当地面两点的距离小于十千米时，地球曲率的影响可以忽略，此时可以采用四参数来完成两种坐标系的转换，精度也能达到预定的需求。四参数转换计算如公式（1） 。

[x2y2=∆x∆y+mcosα-sinαsinαcosαx1y1]    （1）

其中，x1、y1为原始坐标，x2、y2为目标坐标，单位米，Δx、Δy为两个坐标平移参数，即两个平面坐标系的坐标原点之间的坐标差值。α为两个坐标系的坐标轴间的旋转角度，通过旋转α，可以使两个坐标系的X和Y轴重合在一起。M为尺度缩放因子，即两个坐标系内的同一段直线的长度比值，实现尺度的比例转换。通常m值非常接近于1。因为有四个参数所以称之为四参数坐标转换。这样通过坐标的平移、坐标轴的旋转和尺度的缩放，可以实现平面坐标的转换。而这四个参数如何求解是坐标转换的关键步骤。

1.2 四参数的求解

四参数求解是坐标转换的一个逆过程，求解需要至少两对已知点，分别知道其在目标平面的坐标和原始平面的坐标，根据间接平差方法，令[a=m*cosα-1] ， [b=m*sinα]， 按最小二乘原理推导[3]，限于篇幅直接给出公式如下：

[L=x21-x11y21-y11…x2i-x1iy2i-y1i…x2n-x1ny2n-y1n，B=1001…10…10…01…01x11-y11y11x11…x1iy1i…x1ny1n…-y1ix1i…-y1nx1n，]

[X=?x∆yab，P=1001000000001001]   （2）

因为一个点只能建立两个误差方程，要解四个未知数，至少需要两个点建立四个误差方程，如果有多的点，就使用最小二乘法。

求得结果X=（BTPB）-1（BTPL），可以得到Δx、Δy、a、b的值，进而可以解得α=acrtan[b/（a+1）]、m=sqrt[（a+b）2+b2]，得到Δx、Δy、α、m四个参数值[4-5]。

1.3 四参数转换的可靠精度

转换参数的求取时，所采用的已知点的位置分布需要注意如下几个方面：已知点的选取位置应分布合理均匀，分布于测区四周和中心。为达到既定转换精度，尽量采用多个已知点，让这些点位能完全并均匀覆盖整个转换区域。留取几个检查点，作为检核。为了保证参数求解精度，需要注意以下几点：

1） 在参数求解前，对已知坐标进行预处理，保证数据准确。

2） 控制坐标分布、密度。

3） 残差分析。参数求解后对残差进行综合分析，剔除误差较大的点，提高可靠性。

2 Excel的电子表格编辑设计

在 Excel中，函数实际上是一个预先编写的特定计算公式。按照这个特定的计算公式对一个或多个值执行计算，并得出一个或多个运算结果，叫作函数值。使用这些函数不仅可以完成许多复杂的计算，而且还可以简化公式的繁杂程度。使用函数可以简化或缩短工作表中的公式，使数据处理简单方便。

Excel每一个函数都由一系列的参数来控制其输入输出结果。函数的参数可以是数值、文本、单元格引用等数据类型，将函数的参数传入函数中进行处理得到结果。其中，一些函数会有循环嵌套，即通过一个函数调用另一个函数的方式进行计算，从而实现更为复杂的计算需求。

2.1 所用函数

实现坐标转换需要用到矩阵运算，Excel中有矩阵运算相关函数可以方便得到所需结果，常见的矩阵函数有：1） 矩阵转置：TRANSPOSE（原矩阵）；2） 矩阵乘法：A、B矩阵相乘MMULT（矩阵A，矩阵B）；3） 求逆矩阵：MINVERSE（原矩阵）；4） 求矩阵行列式的值 MDETERM（矩阵）；5） INDIRECT函数。在使用矩阵函数时候在公式编辑栏中输入公式：如“=MINVERSE（矩阵A）”，按下“Shift+Ctrl+Enter”组合键，即可得到对应的逆矩阵。

2.2 具体编辑设计

为了防止数据编辑出错、界面简洁，在Excel中新建两个工作表并命名为“输入”和“过程”。在“输入”工作表中编辑输入已知点和坐标转换如图1：

界面数据录入分为目标坐标和源坐标，各自对应X、Y坐标。按行从上到下排列分为P1、P2、P3…Pn点。右侧F列编辑残差列，对应每个点转换后的残差，来监测已知点的数据质量。残差的计算方式为：1） 由已知点求得转换参数；2） 源坐标按第一步的转换参数求出转换后的坐标；3） 转换后的坐标和已知目标坐标X、Y进行求差，再求平方和、开方即为残差。

在“过程”工作表中编辑生成过程矩阵，其中L为2n行1列的矩阵，B为2n行4列的矩阵（n为已知点的数量）。实际应用中已知点数量为2个及以上，数量并不是固定的，因此L、B两个矩阵宜竖向排列，以便在点数增加时进行矩阵扩展如图2所示。

以L矩阵第1、2行为例，第一行输入“=IF（F1>=1，INDIRECT（"输入！D3"）-INDIRECT（"输入！B3"），0）”，第二行输入“=IF（F1>=1，INDIRECT（"输入！E3"）-INDIRECT（"输入！C3"），0）”。F1单元格为已知点数量可以用COUNT函数取得，即已知点数量大于等于1时第一行值为x21-x11否则为0。因为Excel的公式会随着复制剪切同步变化，这会导致已知点编辑时引起后面公式变化，所以加入INDIRECT函数防止剪切单元格时公式跟随变化。以此类推根据公式2中L矩阵编辑其余单元格。

矩阵B第一行第一列输入“=IF（F1>=1，1，0）”即已知点数量大于等于1，第一列值为1，否则为0;第二列固定输入0，第三列输入“=IF（F1>=1，INDIRECT（"输入！B3"），0）”，即已知点大于等于1时，第三列单元格值为“输入”B3单元格的值，否则为零；第四列输入“=IF（F1>=1，INDIRECT（"输入！C3"）*-1，0）”，即已知点F1数量大于1时，第四列单元格值为“输入”C3单元格值乘以-1，否则为零。据此根据需要扩展矩阵的行数，此例为5个点，当输入点数量少于5时，空白的地方利用IF函数根据F1单元格的点数自动填充为0。即已知点数量小于等于5个时，矩阵会自动把不足5的地方填充为零，不影响最终结果的计算。可以达到已知点数量动态使用的效果。如果需要更多的已知点数量，只需要手动把矩阵范围扩展到相应点数即可。

编辑求得B矩阵的转置矩阵BT，因Excel的TRANSPOSE函数使用时需要固定单元格数量，所以矩阵BT采用手动引用编辑生成，以便矩阵的跟随矩阵B动态扩展。

矩阵BT第一行从左到右等于矩阵B的第一列从上到下，矩阵BT第二行从左到右等于矩阵B的第二列从上到下，依次编辑好矩阵BT，结果如图3所示。

编辑构造矩阵BTPB，权系数矩阵P默认为1的单位阵，在此忽略，BTPB结果为4行4列的矩阵，它的大小是不变的，和已知点数量无关。在空白处选中4行4列单元格，在编辑栏输入“=MMULT（矩阵BT，矩阵B）”，然后同按Ctrl+Shift+回车键返回得到结果图4所示：

求逆矩阵（BTPB）-1，逆矩阵大小和源矩阵大小相同，在空白处选中4行4列单元格，在编辑栏输入“=MINVERSE（矩阵BTPB）” 然后同按Ctrl+Shift+回车键返回得到结果图5所示：

编辑构造矩阵BTPL， 权系数矩阵P默认为1的单位阵，在此忽略。矩阵BTPL结果为4行1列的矩阵，它的大小也是不变的，和已知点数量无关。空白处选中4行1列单元格，在编辑栏输入“=MMULT（矩阵BT，矩阵L）”，然后同按Ctrl+Shift+回车键返回得到结果图6左边BTPL：

求方程的最小二乘解X=（BTPB）-1*BTPL，其结果为4行1列固定大小矩阵。空白处选中4行1列单元格，在编辑栏输入“=MMULT（矩阵（BTPB）-1， 矩阵BTPL）”，引用好各个参数矩阵的位置，然后同时按Ctrl+Shift+Enter键返回得到矩阵X结果图6右边所示，分别为Δx、Δy、a、b的值。按照公式2最终解得α=acrtan[b/（a+1）]、m=sqrt[（a+b）2+b2]，得到Δx、Δy、α、m四个参数值。

3 实例与结果分析

现有某工程已知点坐标如表1。

如图7所示，输入P1、P2、P3三个点得到四个转换参数平移Δx、Δx、α、m分别为53.193m、-50.134m、6.578E-6弧度、1.0000003。残差均小于1mm。并有残差列可供使用者评估各个已知点的精度。转换后的参数可以用于实际测量中。

4 结束语

本文针对工程测量中重要的坐標转换实现问题进行研究，探讨二维平面坐标转换的实现原理、精度控制和注意事项。经过矩阵运算、实现使用Excel 内部函数进行坐标转换。经页面设计和调试测试达到精度要求。该方法相对于专用软件和Excel VBA编程方法，实现简单，无需复杂的编程过程、学习成本小；过程数据可见，过程中生成的矩阵可以分析数据，更能方便理解其转换的意义、录入保存数据方便，适用于多种测量、教学场景。

参考文献：

[1] 李征航，黄劲松.GPS测量与数据处理[M].3版.武汉：武汉大学出版社，2016.

[2] 冯注龙，丘荣茂.Excel之光：高效工作的Excel完全手册[M].北京：电子工业出版社，2019.

[3] 武汉大学测绘学院测量平差学科组.误差理论与测量平差基础[M].2版.武汉：武汉大学出版社，2009.

[4] 王萌，高远，范咪娜.区域坐标转换模型对比分析[J].矿山测量，2021，49（2）：71-76.

[5] 谢伟杰.平面四参数法坐标转换的实例应用探讨[J].工程建设与设计，2020（13）：46-47，50.

【通联编辑：梁书】