應(yīng)用設(shè)計(jì)矩陣求解不同坐標(biāo)系統(tǒng)轉(zhuǎn)換下的參數(shù)問題

韓 笑，蒲英霞

(1. 南京大學(xué) 地理與海洋科學(xué)學(xué)院，江蘇 南京 210023；2. 江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210023；3. 江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，江蘇 南京 210023)

0 引 言

空間基準(zhǔn)和坐標(biāo)系統(tǒng)是大地測量、地圖制圖和“3S”技術(shù)中非常重要的組成部分，為了獲取不同用途的地理空間數(shù)據(jù)，我們通常需要在不同空間基準(zhǔn)（如北京1954、西安1980、WGS-84、CGCS2000等）和不同坐標(biāo)系統(tǒng)（地理坐標(biāo)、投影坐標(biāo)、空間直角坐標(biāo)等）之間相互轉(zhuǎn)換[1]。如當(dāng)采用GPS數(shù)據(jù)和地面控制點(diǎn)數(shù)據(jù)建立具有空間參考的衛(wèi)星影像時(shí)，需要將GPS數(shù)據(jù)和地面控制點(diǎn)數(shù)據(jù)統(tǒng)一到相同的坐標(biāo)系中；當(dāng)合并不同數(shù)據(jù)來源的地圖數(shù)據(jù)時(shí)，也需要先統(tǒng)一他們的空間參考系然后進(jìn)行合并處理。

坐標(biāo)變換是將空間數(shù)據(jù)從一個(gè)坐標(biāo)系統(tǒng)轉(zhuǎn)換到另一個(gè)坐標(biāo)系統(tǒng)的過程，其實(shí)質(zhì)是根據(jù)兩個(gè)坐標(biāo)系下的公共已知點(diǎn)求解坐標(biāo)變換參數(shù)[2]。程新輝等人采用四參數(shù)法將北京1954轉(zhuǎn)換至西安1980坐標(biāo)系，并進(jìn)行了轉(zhuǎn)換誤差的分析[3]；王解先等人討論了將WGS-84地心坐標(biāo)轉(zhuǎn)換到北京1954直角坐標(biāo)的兩個(gè)模型——平面轉(zhuǎn)換模型和空間轉(zhuǎn)換模型，并且比較了空間轉(zhuǎn)換模型中三參數(shù)和七參數(shù)對(duì)模型轉(zhuǎn)換結(jié)果的影響[4]。黎舒等研究了從西安1980到CGCS2000坐標(biāo)系的變換方法，并根據(jù)具體情況提出了二維四參數(shù)、七參數(shù)和三維七參數(shù)等模型[5]。在坐標(biāo)變換的誤差處理方面，較為成熟且應(yīng)用最廣的是在求解參數(shù)的過程中采用最小二乘法[6]。孔建基于整體最小二乘法推導(dǎo)了四參數(shù)仿射變換的參數(shù)估計(jì)公式，避免了常規(guī)矩陣分解方法計(jì)算的復(fù)雜性[7]；張鵬杰等提出利用加權(quán)整體最小二乘法進(jìn)行參數(shù)求解，同時(shí)考慮原始坐標(biāo)和目標(biāo)坐標(biāo)的誤差，并根據(jù)誤差的影響程度賦予不同的權(quán)值[8]?，F(xiàn)有的坐標(biāo)變換方法研究中，大多是基于給定的兩個(gè)坐標(biāo)系統(tǒng)來研究某一種坐標(biāo)轉(zhuǎn)換方法，往往是獨(dú)立闡述，沒有試圖挖掘不同轉(zhuǎn)換方法之間的聯(lián)系。在利用最小二乘法推導(dǎo)參數(shù)估計(jì)表達(dá)式時(shí)，我們發(fā)現(xiàn)不同坐標(biāo)變換方法的推導(dǎo)過程存在很大的相似性。本文在Allan研究的基礎(chǔ)上[9]，引入設(shè)計(jì)矩陣的概念，以此作為紐帶將不同的坐標(biāo)變換方法聯(lián)系起來，推導(dǎo)出了基于最小二乘法的統(tǒng)一參數(shù)估計(jì)表達(dá)式。

1 不同坐標(biāo)轉(zhuǎn)換方法

在進(jìn)行坐標(biāo)變換的過程中，若已知兩個(gè)不同坐標(biāo)系統(tǒng)下足夠數(shù)量的坐標(biāo)點(diǎn)，即觀測值個(gè)數(shù)大于必要觀測數(shù)時(shí)，則可以采用最小二乘法求解坐標(biāo)轉(zhuǎn)換參數(shù)的最優(yōu)估計(jì)值。一般地，根據(jù)一組已知數(shù)據(jù)集，利用最小二乘法解算未知參數(shù)的過程如下：

首先，建立已知坐標(biāo)數(shù)據(jù)與未知坐標(biāo)轉(zhuǎn)換參數(shù)之間的一般函數(shù)關(guān)系式

在公式（1）中，當(dāng)達(dá)到最優(yōu)估計(jì)時(shí)，可得如下的微分關(guān)系式

將（3）式代入（4）式中，可以簡化為

A，C是雅克比矩陣，在這里稱之為設(shè)計(jì)矩陣[9]。

實(shí)際上，已知點(diǎn)的觀測值總會(huì)因?yàn)榱烤V或者觀測條件的不同而產(chǎn)生系統(tǒng)誤差，為減小誤差，我們?cè)谇蠼膺^程中計(jì)算殘差平方的加權(quán)和，因而引入權(quán)矩陣W，可以推導(dǎo)出以下表達(dá)式：

式中，W是一個(gè)對(duì)角矩陣，且對(duì)角線元素對(duì)應(yīng)著各自方差的倒數(shù)[9]。

基于公式（6）待求參數(shù)的一般表達(dá)式，問題的關(guān)鍵就變成求解不同坐標(biāo)變換模型對(duì)應(yīng)的設(shè)計(jì)矩陣。下面針對(duì)幾個(gè)常用的坐標(biāo)轉(zhuǎn)換模型，給出設(shè)計(jì)矩陣的推導(dǎo)過程。

1.1 相似變換

對(duì)于二維坐標(biāo)系下的相似變換，我們通常采用四參數(shù)模型，即經(jīng)過坐標(biāo)系的平移、旋轉(zhuǎn)和比例變化實(shí)現(xiàn)[10-11]。該方法適用于坐標(biāo)軸正交且各個(gè)方向的尺度因子相同的平面坐標(biāo)系之間的轉(zhuǎn)換，需要至少兩個(gè)控制點(diǎn)來估計(jì)變換參數(shù)，具體表達(dá)如下：

設(shè)（Xs,Ys）表示源坐標(biāo)系下的坐標(biāo)值，（XT,YT）表示目標(biāo)坐標(biāo)系下的坐標(biāo)值，則有函數(shù)關(guān)系式：

式中，a0、b0分別為X、Y方向的平移參數(shù)，a=μcosα，b=μsinα，α為旋轉(zhuǎn)角度，μ為兩坐標(biāo)軸的尺度參數(shù)。

設(shè)計(jì)矩陣A，C可以表示為：

進(jìn)一步解算出：

以上的演算過程是針對(duì)單一坐標(biāo)點(diǎn)的。通常情況下，會(huì)提供多個(gè)已知點(diǎn)的坐標(biāo)值，這時(shí)我們可以用Ai，Ci來表示第i個(gè)點(diǎn)的設(shè)計(jì)矩陣，將所有已知點(diǎn)的設(shè)計(jì)矩陣組合在一起，可得表達(dá)式：

1.2 仿射變換

仿射變換是基于仿射坐標(biāo)系而建立的一種坐標(biāo)變換數(shù)學(xué)模型，是經(jīng)過坐標(biāo)系的平移、比例、旋轉(zhuǎn)、對(duì)稱和錯(cuò)切等復(fù)合變換得到的[12]。仿射變換適用于坐標(biāo)軸不正交或者各個(gè)方向尺度因子不同的參考系，要求至少3個(gè)控制點(diǎn)來求解變換參數(shù)，其基本關(guān)系式為：

式中，a0、b0分別為X、Y方向的平移參數(shù)，a1=μxcosα，b1=μxsinα；a2=μysinα ，b2=μycosα，α為旋轉(zhuǎn)角度，μx、μy分別為X、Y軸的尺度參數(shù)。特別地，當(dāng) μx=μy時(shí)，有a1=b2，a2=b1，此時(shí)的六參數(shù)仿射變換關(guān)系式與四參數(shù)相似變換相同，即四參數(shù)相似變換是六參數(shù)仿射變換的特殊形式。

第i個(gè)點(diǎn)的設(shè)計(jì)矩陣Ai，Ci：

將所有已知點(diǎn)的設(shè)計(jì)矩陣組合在一起：

1.3 七參數(shù)相似變換法

三維地心參考系下的坐標(biāo)變換，通常采用七參數(shù)變換模型[13]，即3個(gè)平移量、3個(gè)旋轉(zhuǎn)角和1個(gè)尺度因子，該變換需要至少3個(gè)以上的控制點(diǎn)，其模型表達(dá)式為：

式中，ΔX、ΔY、ΔZ分別為XS、YS、ZS方向的偏移量，R為旋轉(zhuǎn)矩陣，αx、αy、αz分別為XS、YS、ZS坐標(biāo)軸的旋轉(zhuǎn)角。

一般地，比例因子接近于1，所以我們可以用（1+κ）表示比例因子，κ為極小量。第i個(gè)點(diǎn)的設(shè)計(jì)矩陣為：

這種轉(zhuǎn)換方法最普遍的應(yīng)用是一個(gè)大地基準(zhǔn)向另一個(gè)大地基準(zhǔn)的轉(zhuǎn)化，這些情況下旋轉(zhuǎn)角都非常小（幾乎只是幾秒），因此，XS＞＞ YS，則有：

同理，

將以上代入（16）、（18）式中，有：

閉合差向量：

將所有已知點(diǎn)的設(shè)計(jì)矩陣組合在一起，有：

2 算例驗(yàn)證

本文選取了某地5個(gè)控制點(diǎn)分別在北京1954和西安1980坐標(biāo)系統(tǒng)下的坐標(biāo)值，采用二維坐標(biāo)系統(tǒng)下的四參數(shù)相似變換和六參數(shù)仿射變換模型，推導(dǎo)設(shè)計(jì)矩陣，求解變換參數(shù)。表1給出了上述5個(gè)控制點(diǎn)分別在兩個(gè)不同坐標(biāo)系統(tǒng)下的高斯平面直角坐標(biāo)值。

表1 控制點(diǎn)在北京1954和西安1980坐標(biāo)系下的坐標(biāo)值Tab.1 Coordinates of control points in Beijing 1954 and Xi'an 1980

應(yīng)用公式(8)，可以計(jì)算出相似變換下的設(shè)計(jì)矩陣A、C，代入公式（6）求得相似變換參數(shù)：a0=-229.386 7，b0=-366.942 8，a=1.000 008 3，b=0.000003 84。

同理，應(yīng)用公式(12)，求出仿射變換下的設(shè)計(jì)矩陣A和C，代入公式（6）求得變換參數(shù)：a0=-200.398 66，a1=1.000 006 24，a2=0.000 003 306，b0=-303.194 67，b1=-0.000 008 182，b2=1.000 007 11。

選取CGCS2000和西安1980坐標(biāo)系下的F、G、H、I 4個(gè)控制點(diǎn)坐標(biāo)，需要將CGCS000的大地坐標(biāo)和西安1980的高斯平面坐標(biāo)轉(zhuǎn)化到空間直角坐標(biāo)系中（見表2），然后采用三維坐標(biāo)系下的七參數(shù)變換模型，推導(dǎo)出設(shè)計(jì)矩陣，求解變換參數(shù)。

表2 控制點(diǎn)在CGCS2000和西安1980下的坐標(biāo)值Tab.2 Coordinates of control points in CGCS2000 and Xi'an 1980

根據(jù)公式(19)、（20）,列出七參數(shù)變換的設(shè)計(jì)矩陣A和C，由公式（6）計(jì)算得到相應(yīng)的變換參數(shù)：k=-2.219E-06，αx=-9.506 7E-07，αy=9.218E-06，αz=-1.295E-05，ΔX=190.248，ΔY=103.828，ΔZ=30.551 6。

接下來，將3種變換方法計(jì)算得到的變換參數(shù)應(yīng)用到該區(qū)域內(nèi)的11個(gè)控制點(diǎn)，分別計(jì)算這些點(diǎn)在北京1954轉(zhuǎn)西安1980的相似變換、北京1954轉(zhuǎn)西安1980的仿射變換、CGCS2000轉(zhuǎn)西安1980的七參數(shù)相似變換下的坐標(biāo)值，比較計(jì)算值與真實(shí)值的差異(表3，表4，表5)。相似變化的X、Y均方根誤差分別為0.1353和0.1755，仿射變化的X、Y均方根誤差分別為0.1407和0.0127，七參數(shù)相似變化的X、Y、h上的均方根誤差分別為0.0033、0.0022和0.5178，均在可接受范圍內(nèi)。比較相似變換和仿射變換的均方根誤差可以看出，相似變換與仿射變換在X方向的均方根誤差比較接近，但在Y方向的均方根誤差相差較大。如我們所知，假設(shè)坐標(biāo)軸正交，且各個(gè)方向的尺度因子相同，則仿射變換的表達(dá)式與相似變換相同，即a1=b2=a，a2=-b1=b。回顧上面參數(shù)估計(jì)的結(jié)果發(fā)現(xiàn)，a1、b2的估計(jì)值與a的估計(jì)值相近，a2的估計(jì)值與b的估計(jì)值也較為接近，但-b1的值與b的值相差較大，因此，相似變換與仿射變換在Y坐標(biāo)的計(jì)算上差異較大。不同于前面兩者，七參數(shù)相似變換是適用于三維直角坐標(biāo)系的，其考慮的參數(shù)更加充分，因此，在X、Y方向的誤差值相對(duì)較小，計(jì)算結(jié)果更接近真實(shí)值。

表3 相似變換的計(jì)算值和真值比較Tab.3 Comparison between calculated values and true values in similarity transformation

表4 仿射變換計(jì)算值與真值比較Tab.4 Comparison between calculated values and true values in aきne transformation

表5 七參數(shù)相似變換計(jì)算值與真值比較Tab.5 Comparison between calculated values and true values in sevenparameter similarity transformation

3 結(jié)束語

本文在常用坐標(biāo)轉(zhuǎn)換模型研究的基礎(chǔ)上，引用設(shè)計(jì)矩陣的概念，根據(jù)最小二乘法原理推導(dǎo)了坐標(biāo)變換統(tǒng)一的參數(shù)估計(jì)表達(dá)式。以四參數(shù)相似變換、六參數(shù)仿射變換和三維七參數(shù)相似變換為例，分別闡述了應(yīng)用設(shè)計(jì)矩陣求解變換參數(shù)的過程，并將求得的參數(shù)應(yīng)用于其他控制點(diǎn)，比較計(jì)算值與真實(shí)值的誤差，有效驗(yàn)證了引用設(shè)計(jì)矩陣的合理性。與前人的工作相比，本文的特點(diǎn)在于以設(shè)計(jì)矩陣為紐帶，將不同的坐標(biāo)變換方法聯(lián)系起來，對(duì)于給定的坐標(biāo)變換方法，只要求出其設(shè)計(jì)矩陣，代入?yún)?shù)估計(jì)表達(dá)式進(jìn)行求解，簡化了計(jì)算過程。需要注意的是，在進(jìn)行坐標(biāo)變換時(shí)應(yīng)該根據(jù)點(diǎn)的分布情況，合理地選擇控制點(diǎn)。在滿足控制點(diǎn)合理分布的條件下，控制點(diǎn)個(gè)數(shù)越多，計(jì)算結(jié)果往往更加精確。特別地，自2008年啟用CGCS2000坐標(biāo)系以來，隨著國家經(jīng)濟(jì)建設(shè)的發(fā)展，地方坐標(biāo)系測繪成果向CGCS2000轉(zhuǎn)換的需求越來越大，本文的研究有助于提高地理坐標(biāo)系轉(zhuǎn)換的效率，推廣CGCS2000坐標(biāo)系的使用。