作者：易衡　郑茂仪

来源：《科技视界》2015年第23期

        【摘 要】以相似变换四参数模型、仿射变换六参数模型和正形变换六参数模型，建立最小二乘法为准则进行间接平差的数学模型，利用Excel软件的宏程序VBA进行编程，求局部的不同坐标系间的转换参数，最后利用这些参数实现不同坐标系之间的相互转换。利用高斯投影正反算公式和换带计算公式，利用Excel的公式计算功能实现高斯投影的正反算和换带计算。

        【关键词】Excel宏程序；最小二乘算法；高斯投影正反算；相似变换；仿射变换；正形变换；换带计算

        0 引言

        随着GPS，测量机器人等测量仪器的不断更新，工程建设中也越来越多的出现不同复杂情况，这就使坐标转换的问题凸显出来。例如，北京54和西安80坐标换算是工程建设中经常遇见的问题，可使不同坐标系统下的测绘成果得到充分利用，因此快速准确地进行坐标转换是非常重要的。再如，在GPS测量中，很多实测的点可能是大地坐标的形式输出，那么要将大地坐标转换为空间直角坐标也是要非常快速和重要的。Excel作为常用的办公软件，几乎在每台微机上都有配备，所以使用Excel做相应的数据转换是非常方便且容易获取的。

        在数字化测图时，应当是先进行控制测量，后进行地形图测量。但是有许多时候，例如当要求在短期内要完成测图，而当地国家坐标又不能查到，为了赶时间，只得先采取假定坐标系进行控制和测图，最后采用坐标转换的方法转换成国家标准坐标系地形图。在种情况下，可以采用“坐标换算”、“图上转换”等方法解决此类问题。借助Excel中的VBA编程，就能方便的批量转换坐标，使得工作变得更加方便和快捷。

        1 不同坐标系之间的相互转换

        1.1 数学模型的建立

        1.1.1 相似变换四参数模型

        x =x Scosa+y Ssina+△x y =-x Ssina+y Scosa+△y （1）

        其中，（x2，y2）原坐标系下的坐标，（x1，y1）为目的坐标系下的坐标，（△x0， △y0）表示平移参数，a表示旋转参数，S表示尺度参数。

        1.1.2 仿射变换六参数模型

        x =ax +by +cy =dx +ey +f（2）

        其中a和e分别确定点（x2，y2）在输出坐标中x1方面和y1方向上的缩放尺度。而b和d确定旋转角度，c和f分别确定在x1方向和y1方向上的水平平移尺寸。

        1.1.3 正形变换六参数模型

        x =p +x p +（x 2-y 2）p -2x y q y =q +y p +x q +2x y p +（x 2-y 2）q （3）

        1.2 转换模型的平差矩阵

        对于n个已知点对，公式（1）的相似变换模型可以建立的误差方程为

        V V = 0 1 -x y2 △x △y δ μ-x y （4）

        对于n个已知点对，公式（2）的仿射变换模型可以建立的误差方程为

        V V = 0 0 0 x y2 1abcdef-x y （5）

        对于n个已知点对，公式（2）的正形变换模型可以建立的误差方程为

        V V = 0 y2 2x y2 1 x （x -y ）p p p q q q -x y （6）

        2 VBA程序设计

        2.1 建立矩阵运算模块

        因为在运算中要进行矩阵的计算，所以要建立一个新的模块，这个模块有一系列的函数可以调用，从而实现矩阵求逆、转置、相乘等相关矩阵运算。

        这是相关函数的自定义：

        Public Sub jzzz（a（） As Double， at（） As Double） '矩阵转置

        Public Sub jzxc（a（） As Double， b（） As Double， c（） As Double） '矩阵相乘

        Private Function h（a（） As Double） As Integer '计算矩阵行数

        Private Function l（a（） As Double） As Integer '计算矩阵列数

        Public Function jzqn（mtxAR（） As Double） As Boolean '矩阵求逆

        2.2 界面设计

        在Excel表格中设置已知的参考坐标系下的已知坐标和两个坐标系的同点坐标对，并固定至少三个同点坐标对如图1所示：

        图1 Excel表格输入设置

        运行宏命令，在宏编辑器中设计操作界面如图2所示Z

        在计算按钮双击事件中添加运算代码，

        Private Sub CommandButton1_Click（）

        。。 。。

        jzzz a（）， at（） ‘矩阵转置

        jzxc at（）， a（）， c（） ‘矩阵相乘

        jzqn c（） ‘矩阵求逆

        jzxc c（）， at（）， d（）‘矩阵求逆

        jzxc d（）， l（）， cs（）‘矩阵求逆

        。。 。。

        图2 图形界面

        以上代码能够自动读取Excel表格中已知的同点坐标点和已知点数据，并进行矩阵间接平差计算，利用仿射变换六参数解算，最终得到转换后的计算结果如图3：

        图3 计算结果

        设计好的Excel宏，在使用的时候，只要点击Excel菜单栏上的工具>>宏（M）>>运行宏或者直接使用快捷操作Alt+F8就能启动宏。此宏命令在平面坐标转换中如北京54平面坐标转西安80平面坐标可以相当快捷方便的使用，为工程计算带来便利。

        3 高斯投影转换

        在工程中经常会遇见这样的问题，很多GPS测量设备或遥感数据只能提供当地的大地坐标，而真正工程使用却是平面坐标，那么就要通过确定的经线来进行高斯平面投影或换带计算，那么这些计算是相当复杂和繁重的。这时，我们可以利用Excel公式计算功能，将既定的高斯投影公式输入到Excel表格中，从而实现准确迅速的计算。

        3.1 坐标正算

        x=X+Ntcos Bl

        + cos B（5-t +9η +4η ）l + cos B（61-58t +t ）l

        y=Ncos Bl

        1+ cos B（1-t +η ）l + cos B（5-18t +t +14η -58t η ）l （7）

        为了计算方便，改写高斯正算公式如公式（7）

        将公式输入Excel表中X=B9+D10*F9*（0.5*POWER（D9，2）+（1/24）*（5-POWER（F9，2）+9*B10+4*POWER（B10，2））*POWER（D9，4）+（1/720）*（61-58*POWER（F9，2）+POWER（F9，4））*POWER（D9，6））

        Y=D10*（D9+（1/6）*（1-POWER（F9，2）+B10）*POWER（D9，3）+（1/120）*（5-18*POWER（F9，2）+POWER（F9，4）+14*B10-58*B10*POWER（F9，2））*POWER（D9，5））

        Excel计算流程， 1>输入已知点的大地经纬度坐标B、L、L0及椭球参数2>计算参数X、m、t、l、η2、N、L3>计算待求数据项X，Y。

        计算结果如图4。

        3.2 高斯投影坐标反算

        本实例设计目的在于将投影平面坐标反算到大地坐标。

        为了计算方便，改写高斯反算公式如下：

        l=

        1- （1+2t +η ） + （5+28t +24t +6η -8η t ）

        B=B y

        1- （5+3t +η -9η t ） + （61+90t +45t ） （8）

        将公式（8）输入Excel表格中

        L=F6+F10

        B=B9-（（1+B10）/PI（））*F9*（90*D9*D9-7.5*（5+3*F9*F9+B10-9*B10*F9*F9）*POWER（D9，4）+0.25*（61+90*F9*F9+45*POWER（F9，4））*POWER（D9，6））

        Excel计算流程，1>输入已知点的空间直角坐标X，Y，L0及椭球参数 2>计算参数Bf、n、t 、n 、X、l3>计算待求数据项B，L。

        计算结果如图4。

        3.3 高斯换带坐标换算

        高斯投影虽然保证了角度没有变形这一优点，但其长度变形较为严重。为了高斯投影的长度变形，必须依子午线进行分带，把投影范围在子午线东、西两侧一定的狭长带内分别进行。但这又使得统一的坐标系分割成各带的坐标系。于是，因分带的结果产生了新的矛盾，即在生产建设中提出各相邻带的相互联系问题。这个问题是通过由一个带的平面坐标换算到相邻带的平面坐标。

        利用高斯投影的正反算公式，亦可进行不同投影带坐标的换带计算。其计算步骤如下：

        （1）根据高斯投影坐标x，y，反算得纬度B和经度差l；

        （2）由子午线的经度L0，求得经度 L=L0 +l；

        （3）根据换带后新的子午线经度L0'，计算相应的经差；

        （4）由高斯投影正算，求得新的高斯投影坐标x'，y'。

        Excel计算流程，1>输入已知点的原坐标、原坐标带L0、拟转坐标带L0’及椭球参数2>计算参数项B η 、B、l、X、m、η 、n、L、t 、n 、l’、X’、t、N3>计算待求数据项X，Y。计算结果如图4。

        4 结论与建议

        将Excel方便的宏VBA和公式计算功能应用于测量快速计算中，不仅有良好的界面，而且加快了数据处理的进程，为工程提供可视化计算和表格处理，让非专业的测量人员也可以方便使用。

        本文仅在不同参考平面坐标系下的相互转换和高斯投影的相关坐标转换中做了简要的方法和理论讲述，在实际运用中，可以根据不同的工程项目特点设计Excel计算命令程序，使得工程计算中大大提高坐标转换的计算效率。

        【参考文献】

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

        [2]孔祥元，郭际明，刘宗泉.大地测量学基础[M].武汉大学出出版社.

        [3]李庆扬，等，编.数值分析[M].施普林格出版社.

        [4]徐振明，等，编.Visual Basic.Net程序设计与应用[M].中国水利水电出版社.

        [5]赵丹.基于vc++的常用大地坐标转换程序实现[J].铁道勘测与设计，2009.

        [6]严蔚敏，等，编.数据结构[M].清华大学出版社.

        [7]李宝玉，等，编.大比例尺数字化测图技术[M].西南交通大学出版社.

        [8]同济大学数学系编著.线性代数[M].高等教育出版社.

        [9]同济大学应用数学系主编.高等数学[M].高等教育出版社.

        [10]潘正风，等，编.数字测图原理与方法[M].武汉大学出版社.

        [11]杨晓明，等，编.数字测绘基础[M].测绘出版社.

        [责任编辑：刘展]下载本文