孫瑞舉

(山東正元地理信息工程有限責(zé)任公司，山東濟(jì)南250101)

一、引 言

坐標(biāo)變換在測(cè)繪學(xué)科中應(yīng)用廣泛，自2008年我國(guó)啟用2000國(guó)家大地坐標(biāo)系后，要求把以往的1954北京坐標(biāo)系、1980西安坐標(biāo)系下的各類(lèi)測(cè)繪成果轉(zhuǎn)換成2000國(guó)家大地坐標(biāo)系，這就涉及到坐標(biāo)系的變換問(wèn)題。

二、二維坐標(biāo)的轉(zhuǎn)換

點(diǎn)的二維坐標(biāo)可以用向量［x y］T表示，假設(shè)某點(diǎn)坐標(biāo)轉(zhuǎn)換前后的坐標(biāo)分別為［x1y1］T和［x2y2］T，則二維坐標(biāo)的坐標(biāo)轉(zhuǎn)換公式為

式中，m為坐標(biāo)尺度參數(shù)，表示坐標(biāo)轉(zhuǎn)換前后坐標(biāo)系的尺度變化情況;θ為坐標(biāo)旋轉(zhuǎn)角，表示轉(zhuǎn)換后的坐標(biāo)系繞前一個(gè)坐標(biāo)系逆時(shí)針旋轉(zhuǎn)的角度;Δx，Δy為平移參數(shù)，表示坐標(biāo)原點(diǎn)的平移量。

在實(shí)際工作中坐標(biāo)轉(zhuǎn)換參數(shù)都是未知的，需要根據(jù)一定數(shù)量的具有新舊坐標(biāo)的控制點(diǎn)來(lái)求得。為了計(jì)算方便，將式(1)改寫(xiě)為如下形式

式中，x0和y0為平移量;a=k·cosθ;b=k·sinθ;k=1+m。

假設(shè)有n個(gè)控制點(diǎn)，其轉(zhuǎn)換前的坐標(biāo)分別為(x1，y1)，(x2，y2)，…，(xn，yn)，轉(zhuǎn)換后的坐標(biāo)為(x1'，y1')，(x2'，y2')…(xn'，yn')，根據(jù)式(2)可以列出如下方程

寫(xiě)成矩陣形式為

式中，

式中有4個(gè)未知量，當(dāng)n=2時(shí)，有唯一解;當(dāng)n＞2時(shí)，可以求得最小二乘解

三、三維坐標(biāo)的轉(zhuǎn)換

三維坐標(biāo)轉(zhuǎn)換的情況與二維坐標(biāo)轉(zhuǎn)換類(lèi)似，也是通過(guò)尺度參數(shù)、旋轉(zhuǎn)角參數(shù)和平移參數(shù)來(lái)確定，只是旋轉(zhuǎn)角有3個(gè)分別繞x軸旋轉(zhuǎn)、繞y軸旋轉(zhuǎn)和繞z軸旋轉(zhuǎn)的旋轉(zhuǎn)角，平移參數(shù)也有3個(gè)，三維坐標(biāo)的坐標(biāo)轉(zhuǎn)換公式為

式中，m為坐標(biāo)尺度參數(shù)，表示坐標(biāo)轉(zhuǎn)換前后坐標(biāo)系的尺度變化情況;Δx，Δy，Δz為平移參數(shù)，表示坐標(biāo)原點(diǎn)在3個(gè)方向的平移量。

角度很小時(shí)取

于是式(7)可以簡(jiǎn)化為

假設(shè)有n個(gè)控制點(diǎn)，其轉(zhuǎn)換前的坐標(biāo)分別為(x1，y1，z1)，(x2，y2，z2)，…，(xn，yn，zn)，轉(zhuǎn)換后的坐標(biāo)為(x'1，y'1，z'1)，(x'2，y'2，z'2)，…，(x'n，y'n，z'n)，根據(jù)式(8)可以列出如下方程

寫(xiě)成矩陣形式為

式中

式中，k=1+m，a=kεx，b=kεy，c=kεz

式中有7個(gè)未知量，當(dāng)n=3時(shí)，有唯一解;當(dāng)n＞3時(shí)，可以求得最小二乘解

四、求解轉(zhuǎn)換參數(shù)

通過(guò)上面的推導(dǎo)，已得到通過(guò)n個(gè)控制點(diǎn)來(lái)求解二維和三維坐標(biāo)的轉(zhuǎn)換參數(shù)公式，即式(5)和式(11)。下面我們討論如何利用VB來(lái)求解參數(shù)，即線性方程組求解。求解線性方程組是測(cè)量程序設(shè)計(jì)中的經(jīng)常遇到的問(wèn)題，本次討論的求解方法是直接法中的列選主元Guass約化法。求解轉(zhuǎn)換參數(shù)的編程思路如下:

1)通過(guò)讀取控制點(diǎn)坐標(biāo)得到系數(shù)項(xiàng)矩陣A及常數(shù)項(xiàng)矩陣L，用數(shù)組A()和L()表示。

2)求得系數(shù)矩陣A的轉(zhuǎn)置矩陣AT，用數(shù)組At()表示。

3)編寫(xiě)矩陣相乘的通用程序，得到AT·A的矩陣(用數(shù)組Aaa()表示)，及得到AT·L的矩陣(用數(shù)組Atl()表示)。

4)編寫(xiě)列選主元Guass約化法求解線性方程組的通用過(guò)程，來(lái)求得未知量

5)根據(jù)未知量得到二維或三維坐標(biāo)的轉(zhuǎn)換參數(shù)。

五、列選主元Guass約化法求解的數(shù)學(xué)過(guò)程

本程序的難點(diǎn)是列選主元Guass約化法的通用程序的編寫(xiě)，下面我們討論其求解思路。

1.Gauss約化法

對(duì)于一般的線性方程組

式中，

使用Gauss約化法求解，就是將上述方程組通過(guò)n-1步約化，轉(zhuǎn)化為上三角方程組

再回代，求此方程組的解。

則方程組得Gauss約化的具體步驟是

用―lik乘第1行再加到第k行，可以消去，具體公式是

4)將其直接回代得

即可以求出線性方程組的解。

2.列選主元Gauss消去法

六、列選主元Guass約化法求解的VB序通用過(guò)程

根據(jù)上述數(shù)學(xué)求解線性方程組的過(guò)程原理，用VB編寫(xiě)列選主元Guass約化法求解通用過(guò)程如下:

七、結(jié)束語(yǔ)

實(shí)際工作中經(jīng)常會(huì)遇到坐標(biāo)系的變換問(wèn)題，尤其是2000國(guó)家大地坐標(biāo)系的使用，如何把以往坐標(biāo)系的成果轉(zhuǎn)換到2000國(guó)家大地坐標(biāo)系，涉及的工作量大，格式不一，通過(guò)本程序可以靈活解決。

［1］孔祥元，梅是義.控制測(cè)量學(xué)［M］.武漢:武漢大學(xué)出版社，2003.

［2］佟彪.VB語(yǔ)言與測(cè)量程序設(shè)計(jì)［M］.北京:中國(guó)電力出版社，2007.