利用VB解決不同坐標系的變換問題

孫瑞舉

(山東正元地理信息工程有限責任公司，山東濟南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.