基于Matlab 的七參數(shù)空間坐標(biāo)轉(zhuǎn)換研究與實(shí)現(xiàn)

劉 平，段志強(qiáng)，謝 超

（1.湖北省測繪成果檔案館，湖北 武漢 430074）

基于Matlab 的七參數(shù)空間坐標(biāo)轉(zhuǎn)換研究與實(shí)現(xiàn)

劉 平1，段志強(qiáng)1，謝 超1

（1.湖北省測繪成果檔案館，湖北 武漢 430074）

分析了七參數(shù)空間坐標(biāo)轉(zhuǎn)換模型和Matlab語言在矩陣運(yùn)算方面的優(yōu)勢。在此基礎(chǔ)上開發(fā)實(shí)現(xiàn)了七參數(shù)空間坐標(biāo)轉(zhuǎn)換程序。簡要介紹了利用Matlab語言開發(fā)此程序的一點(diǎn)體會。

Matlab；七參數(shù)；空間坐標(biāo)轉(zhuǎn)換

1 七參數(shù)空間坐標(biāo)轉(zhuǎn)換基本原理

坐標(biāo)系統(tǒng)間的差異主要來自于坐標(biāo)系統(tǒng)的定義差，即原點(diǎn)位置、坐標(biāo)軸向的定向和尺度的定義差。坐標(biāo)轉(zhuǎn)換時應(yīng)優(yōu)先考慮坐標(biāo)系統(tǒng)定義差異的轉(zhuǎn)換，即首先應(yīng)完成相似變換。在相似變換的基礎(chǔ)上，再考慮對剩余誤差進(jìn)行擬合，使精度較低的坐標(biāo)框架點(diǎn)附合到精度較高的坐標(biāo)系統(tǒng)的框架點(diǎn)坐標(biāo)，使統(tǒng)一后的坐標(biāo)系框架點(diǎn)坐標(biāo)具有較好的一致性。

七參數(shù)空間坐標(biāo)轉(zhuǎn)換，實(shí)際上是2個不同基準(zhǔn)面之間的空間直角坐標(biāo)的相似變換。其基本原理是：假設(shè)2個坐標(biāo)基準(zhǔn)都滿足“剛體”（基準(zhǔn)無變形）條件，則2個不同基準(zhǔn)間存在3個平移量和3個旋轉(zhuǎn)量及1個尺度變換量的差別。

對于既有旋轉(zhuǎn)、縮放又有平移的2個空間直角坐標(biāo)系的坐標(biāo)換算，存在3個平移參數(shù)和3個旋轉(zhuǎn)參數(shù)以及1個尺度變化參數(shù)，共7個參數(shù)。其坐標(biāo)轉(zhuǎn)換數(shù)學(xué)模型為：

一般εX、εY、εZ為微小轉(zhuǎn)角，可取：

cosεX=cosεY=cosεZ=1

sinεX= εX，sinεY= εY，sinεZ= εZsinεXsinεY=sinεXsinεZ=sinεYsinεZ=0于是有：

式（1）即為布爾沙-沃爾夫七參數(shù)空間坐標(biāo)轉(zhuǎn)換模型。此模型的轉(zhuǎn)換步驟是首先作三軸旋轉(zhuǎn)，再統(tǒng)一尺度，最后再平移。如果能求解出這7個參數(shù)值，則所有待轉(zhuǎn)換點(diǎn)的坐標(biāo)轉(zhuǎn)換也就可以求得了。公式中含有7個轉(zhuǎn)換參數(shù)，為了求得7個轉(zhuǎn)換參數(shù)，至少需要3個公共點(diǎn)（也稱為重合點(diǎn)），當(dāng)多于3個公共點(diǎn)時，可按最小二乘法求得7個參數(shù)的最或是值。

則可寫出如下形式的誤差方程：改寫成矩陣形式為：

式中，δX =（?X0，?Y0，?Z0，a1，a2，a3，a4）T，為待求點(diǎn)轉(zhuǎn)換參數(shù)向量；V =（VX2i，VY2i，VZ2i）T為改正數(shù)向量；L=（X2i，Y2i，Z2i）T已知值；B為系數(shù)矩陣。

根據(jù)最小二乘法原理，可列出法方程為：

其解為：

當(dāng)已知公共點(diǎn)為等精度時，上述權(quán)矩陣為單位權(quán)矩陣，于是有：

單位權(quán)方差：

式中，n為參與求解轉(zhuǎn)換參數(shù)的點(diǎn)數(shù)。

2 七參數(shù)空間坐標(biāo)轉(zhuǎn)換技術(shù)流程

1）收集整理（或施測）坐標(biāo)轉(zhuǎn)換區(qū)域范圍的已知公共點(diǎn)空間坐標(biāo)成果；

2）分析選取用于計算坐標(biāo)轉(zhuǎn)換參數(shù)的公共點(diǎn)，原則是選取能覆蓋轉(zhuǎn)換區(qū)域且精度較高、分布均勻的已知公共點(diǎn)成果；

3）根據(jù)上述七參數(shù)模型計算坐標(biāo)轉(zhuǎn)換參數(shù)；

4）分析用于計算坐標(biāo)轉(zhuǎn)換參數(shù)的公共點(diǎn)的坐標(biāo)轉(zhuǎn)換殘差和殘差中誤差，以及未參與求解參數(shù)的公共檢查點(diǎn)坐標(biāo)轉(zhuǎn)換殘差和殘差中誤差，查找并剔除用于計算坐標(biāo)轉(zhuǎn)換參數(shù)的粗差公共點(diǎn)；

5）用可靠公共點(diǎn)重新計算坐標(biāo)轉(zhuǎn)換參數(shù)，直到精度滿足規(guī)定限差為止。

6）根據(jù)最終合格的轉(zhuǎn)換參數(shù)，代入式（2），求解出所有待轉(zhuǎn)換點(diǎn)的轉(zhuǎn)換坐標(biāo)。

如果已知公共點(diǎn)是高斯平面坐標(biāo)，則在作七參數(shù)坐標(biāo)轉(zhuǎn)換前需要將其經(jīng)高斯投影反算得到大地經(jīng)緯度B、L，而大地高H=正常高+高程異常。然后再根據(jù)同基準(zhǔn)下的大地坐標(biāo)與空間坐標(biāo)轉(zhuǎn)換關(guān)系，將大地坐標(biāo)轉(zhuǎn)換成空間坐標(biāo)。

3 七參數(shù)空間坐標(biāo)轉(zhuǎn)換研究與實(shí)現(xiàn)

3.1 問題的分析處理

根據(jù)矩陣?yán)碚撝匠虃€數(shù)大于未知量個數(shù)的方程組稱為“超定方程組”。求解坐標(biāo)系轉(zhuǎn)換參數(shù)的計算中，一般都要求有多余的公共點(diǎn)參與，即在解求7個轉(zhuǎn)換參數(shù)（7個未知數(shù)）時，方程個數(shù)往往遠(yuǎn)大于7，顯然坐標(biāo)轉(zhuǎn)換計算過程所涉及的都是“超定方程組”。求解“超定方程組”比較常用的方法是最小二乘法。形象地說，就是在無法完全滿足所給定條件的情況下，求解“超定方程組”的一個最或是解。

從上述坐標(biāo)轉(zhuǎn)換數(shù)學(xué)模型可知，在合理選擇了參與求解轉(zhuǎn)換參數(shù)的公共點(diǎn)和外部檢查點(diǎn)后，坐標(biāo)轉(zhuǎn)換過程其實(shí)質(zhì)是根據(jù)參與求解轉(zhuǎn)換參數(shù)的公共點(diǎn)數(shù)量，動態(tài)組成B矩陣和L矩陣（隨著參與求解坐標(biāo)轉(zhuǎn)換參數(shù)的公共點(diǎn)增加，B矩陣和L矩陣也不斷增大）。然后對數(shù)值矩陣進(jìn)行一系列的運(yùn)算，而其中關(guān)鍵運(yùn)算是矩陣求逆。在Matlab中，對形如Ax=b的矩陣方程作求逆運(yùn)算，一般可用“左除”、內(nèi)部數(shù)值求逆函數(shù)求逆及符號運(yùn)算求逆函數(shù)求逆等。

上述式（7）或式（8）中矩陣B為3n×7階的矩陣（n為參與求解參數(shù)的公共點(diǎn)數(shù)），當(dāng)參與求解參數(shù)的公共點(diǎn)為5個時，其B矩陣是一個15×7階的矩陣。顯然，對B矩陣求逆運(yùn)算將是很復(fù)雜且運(yùn)算量很大的工作。在矩陣求逆過程中，如果矩陣接近奇異，Matlab將給出警告信息。

由于參與求解坐標(biāo)轉(zhuǎn)換參數(shù)點(diǎn)的坐標(biāo)值X1i、Y1i、Z1i都很大（一般都有十進(jìn)制的7位整數(shù)），而系數(shù)矩陣B中其他元素僅為1或0，即系數(shù)矩陣B的向量之間長度過于懸殊，致使矩陣B的條件數(shù)過大，且BTB的條件數(shù)更大，即此時矩陣呈病態(tài)趨勢（或稱矩陣接近奇異），求逆時容易產(chǎn)生擾動，造成求解的轉(zhuǎn)換參數(shù)可能不正確。

對上述警告性提示，筆者采用了2種方法處理。一是首先對坐標(biāo)進(jìn)行重心化處理，將所有點(diǎn)的坐標(biāo)平移到重心化坐標(biāo)系中，減小坐標(biāo)的絕對值。求逆完成后，在適當(dāng)處再反算還原到原坐標(biāo)系中。經(jīng)過這樣處理后，警告提示就再沒有出現(xiàn)了。另一種方法是采用Matlab中符號矩陣求逆運(yùn)算，因矩陣中變量元素是符號變量，即變量沒有賦值，避開了數(shù)值矩陣求逆而進(jìn)行符號矩陣求逆，從而避免因數(shù)值計算誤差導(dǎo)致的矩陣接近奇異的情況發(fā)生，完成求逆運(yùn)算后再用相關(guān)函數(shù)將符號矩陣轉(zhuǎn)換成數(shù)值矩陣，也可避免出現(xiàn)“矩陣接近奇異…”的警告提示。

3.2 七參數(shù)空間坐標(biāo)轉(zhuǎn)換的實(shí)現(xiàn)

先用自編程序?qū)^(qū)域散點(diǎn)展繪到圖上，根據(jù)圖上點(diǎn)的分布，均勻選擇參與求解轉(zhuǎn)換參數(shù)的重合點(diǎn)，再將待轉(zhuǎn)換坐標(biāo)系下所有點(diǎn)的空間坐標(biāo)和目標(biāo)坐標(biāo)系下的重合點(diǎn)空間坐標(biāo)建立一個文本文件保存，然后啟動程序。程序設(shè)計上首先通過界面提示輸入?yún)⑴c求解轉(zhuǎn)換參數(shù)的重合點(diǎn)個數(shù)，并默認(rèn)重合定向點(diǎn)（即參與求解參數(shù)的重合點(diǎn)）為5，且最少要大于或等于3個，如輸入重合點(diǎn)數(shù)小于3，程序會提示輸入錯誤。再通過界面分別選擇待轉(zhuǎn)換坐標(biāo)系所有點(diǎn)坐標(biāo)文件和目標(biāo)坐標(biāo)系公共點(diǎn)坐標(biāo)文件，程序?qū)?個文件中的空間坐標(biāo)數(shù)據(jù)分別保存到2個矩陣變量中。通過對坐標(biāo)系重心化處理及一系列運(yùn)算，計算出轉(zhuǎn)換坐標(biāo)、各點(diǎn)殘差及殘差中誤差以及單位權(quán)中誤差和7個轉(zhuǎn)換參數(shù)等，并將所有結(jié)果保存到Excel表格中。

3.3 程序正確性驗(yàn)證

1）經(jīng)多個區(qū)域?qū)嵗龜?shù)據(jù)計算（網(wǎng)上很多實(shí)例數(shù)據(jù)及結(jié)果）對比，程序計算結(jié)果完全正確。

2）將七參數(shù)空間坐標(biāo)轉(zhuǎn)換出的結(jié)果再轉(zhuǎn)換到大地坐標(biāo)，然后繼續(xù)經(jīng)投影正算到高斯坐標(biāo)，對比其原始高斯坐標(biāo)成果，結(jié)果所有同名點(diǎn)差值也都小于1 mm。

3）還開發(fā)了一個平面四參數(shù)坐標(biāo)轉(zhuǎn)換程序，對同一區(qū)域作平面四參數(shù)坐標(biāo)轉(zhuǎn)換，并與七參數(shù)坐標(biāo)轉(zhuǎn)換程序結(jié)果對比，轉(zhuǎn)換計算出的坐標(biāo)值幾乎完全一致，證明了2個獨(dú)立開發(fā)的坐標(biāo)轉(zhuǎn)換程序完全正確。

[1] 孔祥元．大地測量學(xué)基礎(chǔ)[M]．武漢:武漢大學(xué)出版社,2006

[2] 張志涌．精通Matlab[M]．北京:北京航空航天大學(xué)出版社,2003

[3] Moler C．Matlab 數(shù)值計算[M]．北京:機(jī)械工業(yè)出版社,2006

[4] 肖筱南．現(xiàn)代數(shù)值計算方法[M]．北京:北京大學(xué)出版社,2003

[5] 同濟(jì)大學(xué)數(shù)學(xué)教研室．工程數(shù)學(xué)線性代數(shù)[M]．北京:高等教育出版社,1982

[6] 成英燕,李夕銀．適用于不同橢球的高斯平面坐標(biāo)正反算的實(shí)用算法[J]．測繪科學(xué),2004,29(4):26-27

[7] 劉陶勝,黃聲享,羅力,等．基于重心基準(zhǔn)的平面坐標(biāo)轉(zhuǎn)換研究[J]．大地測量與地球動力學(xué),2011,31(2):102-106

[8] 陳宇,白征東,羅騰,等．基于改進(jìn)的布爾沙模型的坐標(biāo)轉(zhuǎn)換方法[J]．大地測量與地球動力學(xué),2010,30(3):71-75

[9] 李瀟,尹暉．基于最小二乘配置的三維空間坐標(biāo)轉(zhuǎn)換[J]．測繪工程,2008,17(2):16-18,29

[10] 劉平,呂海兵,張峰,等．?dāng)?shù)值計算誤差初步理解及簡單分析[J]．地理空間信息,2011,9(3): 156-158

P226

B

1672-4623（2014）06-0062-03

10.3969/j.issn.1672-4623.2014.06.022

劉平，正高職高級工程師，主要從事測繪地理信息內(nèi)業(yè)技術(shù)和質(zhì)量管理工作。

2014-01-26。