利用VB解决不同坐标系的变换问题
2013-12-12孙瑞举
孙瑞举
(山东正元地理信息工程有限责任公司,山东济南250101)
一、引 言
坐标变换在测绘学科中应用广泛,自2008年我国启用2000国家大地坐标系后,要求把以往的1954北京坐标系、1980西安坐标系下的各类测绘成果转换成2000国家大地坐标系,这就涉及到坐标系的变换问题。
二、二维坐标的转换
点的二维坐标可以用向量[x y]T表示,假设某点坐标转换前后的坐标分别为[x1y1]T和[x2y2]T,则二维坐标的坐标转换公式为
式中,m为坐标尺度参数,表示坐标转换前后坐标系的尺度变化情况;θ为坐标旋转角,表示转换后的坐标系绕前一个坐标系逆时针旋转的角度;Δx,Δy为平移参数,表示坐标原点的平移量。
在实际工作中坐标转换参数都是未知的,需要根据一定数量的具有新旧坐标的控制点来求得。为了计算方便,将式(1)改写为如下形式
式中,x0和y0为平移量;a=k·cosθ;b=k·sinθ;k=1+m。
假设有n个控制点,其转换前的坐标分别为(x1,y1),(x2,y2),…,(xn,yn),转换后的坐标为(x1',y1'),(x2',y2')…(xn',yn'),根据式(2)可以列出如下方程
写成矩阵形式为
式中,
式中有4个未知量,当n=2时,有唯一解;当n>2时,可以求得最小二乘解
三、三维坐标的转换
三维坐标转换的情况与二维坐标转换类似,也是通过尺度参数、旋转角参数和平移参数来确定,只是旋转角有3个分别绕x轴旋转、绕y轴旋转和绕z轴旋转的旋转角,平移参数也有3个,三维坐标的坐标转换公式为
式中,m为坐标尺度参数,表示坐标转换前后坐标系的尺度变化情况;Δx,Δy,Δz为平移参数,表示坐标原点在3个方向的平移量。
角度很小时取
于是式(7)可以简化为
假设有n个控制点,其转换前的坐标分别为(x1,y1,z1),(x2,y2,z2),…,(xn,yn,zn),转换后的坐标为(x'1,y'1,z'1),(x'2,y'2,z'2),…,(x'n,y'n,z'n),根据式(8)可以列出如下方程
写成矩阵形式为
式中
式中,k=1+m,a=kεx,b=kεy,c=kεz
式中有7个未知量,当n=3时,有唯一解;当n>3时,可以求得最小二乘解
四、求解转换参数
通过上面的推导,已得到通过n个控制点来求解二维和三维坐标的转换参数公式,即式(5)和式(11)。下面我们讨论如何利用VB来求解参数,即线性方程组求解。求解线性方程组是测量程序设计中的经常遇到的问题,本次讨论的求解方法是直接法中的列选主元Guass约化法。求解转换参数的编程思路如下:
1)通过读取控制点坐标得到系数项矩阵A及常数项矩阵L,用数组A()和L()表示。
2)求得系数矩阵A的转置矩阵AT,用数组At()表示。
3)编写矩阵相乘的通用程序,得到AT·A的矩阵(用数组Aaa()表示),及得到AT·L的矩阵(用数组Atl()表示)。
4)编写列选主元Guass约化法求解线性方程组的通用过程,来求得未知量
5)根据未知量得到二维或三维坐标的转换参数。
五、列选主元Guass约化法求解的数学过程
本程序的难点是列选主元Guass约化法的通用程序的编写,下面我们讨论其求解思路。
1.Gauss约化法
对于一般的线性方程组
式中,
使用Gauss约化法求解,就是将上述方程组通过n-1步约化,转化为上三角方程组
再回代,求此方程组的解。
则方程组得Gauss约化的具体步骤是
用―lik乘第1行再加到第k行,可以消去,具体公式是
4)将其直接回代得
即可以求出线性方程组的解。
2.列选主元Gauss消去法
六、列选主元Guass约化法求解的VB序通用过程
根据上述数学求解线性方程组的过程原理,用VB编写列选主元Guass约化法求解通用过程如下:
七、结束语
实际工作中经常会遇到坐标系的变换问题,尤其是2000国家大地坐标系的使用,如何把以往坐标系的成果转换到2000国家大地坐标系,涉及的工作量大,格式不一,通过本程序可以灵活解决。
[1]孔祥元,梅是义.控制测量学[M].武汉:武汉大学出版社,2003.
[2]佟彪.VB语言与测量程序设计[M].北京:中国电力出版社,2007.