吳繼忠 王安怡

1 南京工業(yè)大學(xué)測繪學(xué)院，南京市浦珠南路30號，211816

2 大連海洋大學(xué)海洋與土木工程學(xué)院，大連市黑石礁街52號，116023

根據(jù)空間維數(shù)的不同，空間直角坐標(biāo)轉(zhuǎn)換一般分為二維平面和三維空間坐標(biāo)轉(zhuǎn)換。二維平面坐標(biāo)轉(zhuǎn)換可直接變換成線性模型，其解算過程容易實現(xiàn)。三維空間直角坐標(biāo)轉(zhuǎn)換中，在旋轉(zhuǎn)角度很小的情況下可以對旋轉(zhuǎn)矩陣作一定的近似處理，從而轉(zhuǎn)化為線性模型進(jìn)行求解。當(dāng)旋轉(zhuǎn)角為大角度時，在小范圍測區(qū)內(nèi)7參數(shù)之間存在強(qiáng)相關(guān)性［1］，更難以處理的是旋轉(zhuǎn)矩陣的各個參數(shù)關(guān)系復(fù)雜，通常進(jìn)行線性化或迭代運(yùn)算［2－3］進(jìn)行求解，這些方法在實際應(yīng)用中均有一定的條件限制或者使用復(fù)雜。本文提出了空間直角坐標(biāo)轉(zhuǎn)換的統(tǒng)一模型，結(jié)合空間直角坐標(biāo)轉(zhuǎn)換中旋轉(zhuǎn)矩陣的正交特性，將模型求解轉(zhuǎn)換為正交Procrustes問題，實現(xiàn)了二維和三維空間直角坐標(biāo)轉(zhuǎn)換的直接線性解算。

1 空間直角坐標(biāo)轉(zhuǎn)換的數(shù)學(xué)模型

二維平面坐標(biāo)轉(zhuǎn)換模型可以表示為：

式中，λ為尺度參數(shù)；為平移參數(shù)；旋轉(zhuǎn)矩陣，φ是旋轉(zhuǎn)角。顯然，R2是正交矩陣，滿足

三維空間直角坐標(biāo)的轉(zhuǎn)換模型表示為：

式中，λ為尺度參數(shù)；為平移參數(shù)；，其中為3個軸向旋轉(zhuǎn)對應(yīng)的旋轉(zhuǎn)矩陣，其構(gòu)成可參見相關(guān)文獻(xiàn)，均為正交矩陣，因此R3也是正交矩陣，滿足

對比式（1）和式（2）可知，二維平面坐標(biāo)轉(zhuǎn)換可以看成是三維空間直角坐標(biāo)轉(zhuǎn)換的一種特殊形式，二者性質(zhì)完全相同，區(qū)別僅在于矩陣的維數(shù)不同。以三維空間直角坐標(biāo)轉(zhuǎn)換為例，設(shè)兩套坐標(biāo)系的公共點數(shù)目為m，各點在不同坐標(biāo)系下的坐標(biāo)記為：

B為公共點在目標(biāo)坐標(biāo)系統(tǒng)下的坐標(biāo)，A為公共點在原坐標(biāo)系統(tǒng)下的坐標(biāo)。顯然，每個公共點均滿足式（2）的條件，將第i（i＝1，2，…，m）個點代入式（2）并將左右兩邊矩陣進(jìn)行轉(zhuǎn)置可得：

其中J＝［1，1，…，1］T為m×1的常數(shù)向量。由上述分析，空間直角坐標(biāo)轉(zhuǎn)換的統(tǒng)一模型可寫為如下形式：

式中，k為常量，當(dāng)k＝2或3時分別對應(yīng)二維平面坐標(biāo)轉(zhuǎn)換和三維空間直角坐標(biāo)轉(zhuǎn)換，m為轉(zhuǎn)換中公共點的數(shù)目，R∈Rk×k是轉(zhuǎn)置后的旋轉(zhuǎn)矩陣，是滿足正交特性的方陣。式（6）即為空間直角坐標(biāo)轉(zhuǎn)換的統(tǒng)一模型。

2 正交Procrustes分析

利用公共點求解轉(zhuǎn)換參數(shù)實際是要求出參數(shù)λ、R、T，使得λAR＋JT盡可能接近B，即求解如下問題：

其中‖·‖F(xiàn)表示Frobenius范數(shù)，R滿足RTR＝I，式（7）本質(zhì)上屬于正交Procrustes問題。正交Procrustes問題可表述為：給定矩陣A∈Rn×p，B∈Rn×q，其中n＞q，求正交矩陣R∈Rp×q滿足‖AR－B‖F(xiàn)＝min。

根據(jù)正交Procrustes分析的解算方法［4］，針對式（6）的求解過程如下：記ε＝λAR＋JT－B，根據(jù)Lagrange原理構(gòu)造極值函數(shù)：

式中，L為乘數(shù)因子矩陣。將式（8）分別對λ、R、T求偏導(dǎo)數(shù)：

由式（11）可得：

由式（10）左乘RT得：

由于RTATAR和（L＋LT）／2均是對稱矩陣，RTATB－RTATJTT也必須是對稱矩陣，即

將式（12）代入式（14）得：

將上式展開，根據(jù)對稱矩陣的特點同理可得：

記D＝AT（I－JJT／m）B，要使式（16）成立，必須

滿足：

將式（17）分別左乘R、右乘RT，得：

將式（18）左右相乘得：

將D進(jìn)行奇異值分解，svd（D）＝UΣVT，其中U和V均為正交矩陣，代入式（19）得：

故旋轉(zhuǎn)矩陣R＝UVT，將式（12）代入式（9）得：

將R、λ代入式（12）即可得到平移參數(shù)T。根據(jù)旋轉(zhuǎn)矩陣R的展開式，可以反算出不同軸向上的旋轉(zhuǎn)角，其過程不再贅述。最后用單位權(quán)中誤差來評定精度，單位權(quán)中誤差計算公式為：

式中，ε＝λAR＋JT－B；f為自由度，在二維轉(zhuǎn)換和三維轉(zhuǎn)換的情況下分別等于2m－4和3m－7，m為公共點數(shù)量。

在計算過程中，首先將公共點坐標(biāo)按式（6）構(gòu)建矩陣A、B和常量矩陣J，再根據(jù)A、B、J構(gòu)建矩陣D＝AT（I－JJT／m）B。對矩陣D進(jìn)行奇異值分解，依次計算得到旋轉(zhuǎn)矩陣R、尺度參數(shù)λ和平移參數(shù)T，最后進(jìn)行精度評定。

3 算例分析

3.1 算例1

以文獻(xiàn)［5］中的數(shù)據(jù)為例，5 個公共點在WGS－84和獨(dú)立坐標(biāo)系下的坐標(biāo)如表1所示。為便于比較，取1、2、3號點求轉(zhuǎn)換參數(shù)，再利用轉(zhuǎn)換參數(shù)將4、5號點轉(zhuǎn)換到地方坐標(biāo)系，并與其已知值進(jìn)行比較。

表1 平面坐標(biāo)轉(zhuǎn)換公共點坐標(biāo)Tab.1 Cartesian coordinates of common points

利用本文模型和方法，首先計算出旋轉(zhuǎn)矩陣，得到的旋轉(zhuǎn)角為359°01′20.575 4″，進(jìn)一步得到尺度因子為3.064 875×10－6，最后由旋轉(zhuǎn)矩陣和尺度參數(shù)得到平移參數(shù)Δx＝－2 465 635.256 m，Δy＝－433 223.055m。利用上述轉(zhuǎn)換參數(shù)將4、5號點轉(zhuǎn)換到地方坐標(biāo)系，結(jié)果見表2。

表2 平面坐標(biāo)轉(zhuǎn)換公共點坐標(biāo)Tab.1 Cartesian coordinates of common points

由表2可知，本文計算結(jié)果與已知結(jié)果間僅存在微小的差異，最大坐標(biāo)較差為3 mm。造成這一差異的原因是點位坐標(biāo)是實測結(jié)果，均包含一定的觀測誤差。

3.2 算例2

以4個點的三維空間直角坐標(biāo)為原始坐標(biāo)，模擬10 套不同轉(zhuǎn)換參數(shù)（表3），利用式（2）將4個控制點的原始坐標(biāo)轉(zhuǎn)換到目標(biāo)坐標(biāo)系下，最后利用本文算法計算出轉(zhuǎn)換參數(shù)并與其模擬的真值進(jìn)行比較，得出各個轉(zhuǎn)換參數(shù)的計算誤差及轉(zhuǎn)換參數(shù)求解的中誤差，結(jié)果列于表4。表4 最后一列為各個方案計算轉(zhuǎn)換參數(shù)的中誤差。

表3 模擬參數(shù)Tab.3 Simulation parameters and solutions

表4 模擬參數(shù)的求解誤差Tab.4 Errors of simulation parameters

從表4可以看出，本文計算方法對旋轉(zhuǎn)角、尺度參數(shù)和平移參數(shù)的數(shù)值范圍沒有限制，因而適合于任意條件下的三維空間直角坐標(biāo)轉(zhuǎn)換。其計算過程為矩陣運(yùn)算，是一種純粹的線性計算，突破了以往三維坐標(biāo)轉(zhuǎn)換方法中對旋轉(zhuǎn)角的限制及計算方法的瓶頸問題，實現(xiàn)了三維坐標(biāo)轉(zhuǎn)換的直接解算。由于模擬數(shù)據(jù)的點位坐標(biāo)不存在誤差，在不同方案中，旋轉(zhuǎn)角的計算誤差小于10－13，尺度參數(shù)計算誤差小于10－14，平移參數(shù)計算誤差小于10－8，計算中誤差的數(shù)量級也在10－8以下，這些計算誤差是計算機(jī)運(yùn)算過程中數(shù)據(jù)截斷誤差的累計所造成的，其大小對實際應(yīng)用而言完全可以忽略不計。

4 結(jié) 語

本文提出了空間直角坐標(biāo)轉(zhuǎn)換的統(tǒng)一方法。首先建立二維平面坐標(biāo)轉(zhuǎn)換和三維空間直角坐標(biāo)轉(zhuǎn)換的統(tǒng)一模型，采用正交Procrustes分析對模型進(jìn)行解算。解算過程不需要參數(shù)的先驗近似值，也不需要線性化處理和迭代計算，對轉(zhuǎn)換參數(shù)的數(shù)值范圍沒有限制，實現(xiàn)了空間直角坐標(biāo)的轉(zhuǎn)換在表示模型和解算方法上的統(tǒng)一。

［1］王解先.七參數(shù)轉(zhuǎn)換中參數(shù)之間的相關(guān)性［J］.大地測量與地球動力學(xué)，2007，27（2）：43－46（Wang Jiexian.Correlations among Parameters in Seven－Parameter Transformation Model［J］.Journal of Geodesy and Geodynamics，2007，27（2）：43－46）

［2］曾文憲，陶本藻.三維坐標(biāo)轉(zhuǎn)換的非線性模型［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2003，28（5）：566－568（Zeng Wenxian，Tao Benzao.Non－Linear Adjustment Model of Three－Dimensional Coordinate Transformation［J］.Geomatics and Information Science of Wuhan University，2003，28（5）：566－568）

［3］姚宜斌，黃承猛，李程春，等.一種適用于大角度的三維坐標(biāo)轉(zhuǎn)換參數(shù)求解算法［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2012，37（3）：253－256（Yao Yibin，Huang Chengmeng，Li Chengchun et al.A New Algorithm for Solution of Transformation Parameters of Big Rotation Angle’s 3D Coordinate［J］.Geomatics and Information Science of Wuhan University，2012，37（3）：253－256）

［4］Crosilla F，Beinat A.Use of Generalised Procrustes Analysis for the Photogrammetric Block Adjustment by Independent Models［J］.ISPRS Journal of Photogrammetry ＆ Remote Sensing，2002，56（2）：195－209

［5］姚宜斌.平面坐標(biāo)系統(tǒng)相互轉(zhuǎn)換的一種簡便算法［J］.測繪信息與工程，2001（1）：1－3（Yao Yibin.A Simple and Convenient Arithmetic for the Transformation of Plane Coordinate System［J］.Journal of Geomatics，2001（1）：1－3