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

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

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

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

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

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

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

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

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

坐標(biāo)系統(tǒng)間的差異主要來(lái)自于坐標(biāo)系統(tǒng)的定義差，即原點(diǎn)位置、坐標(biāo)軸向的定向和尺度的定義差。坐標(biāo)轉(zhuǎn)換時(shí)應(yīng)優(yōu)先考慮坐標(biāo)系統(tǒng)定義差異的轉(zhuǎn)換，即首先應(yīng)完成相似變換。在相似變換的基礎(chǔ)上，再考慮對(duì)剩余誤差進(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個(gè)不同基準(zhǔn)面之間的空間直角坐標(biāo)的相似變換。其基本原理是：假設(shè)2個(gè)坐標(biāo)基準(zhǔn)都滿足“剛體”（基準(zhǔn)無(wú)變形）條件，則2個(gè)不同基準(zhǔn)間存在3個(gè)平移量和3個(gè)旋轉(zhuǎn)量及1個(gè)尺度變換量的差別。

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

一般εX、εY、εZ為微小轉(zhuǎn)角，可?。?/p>

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個(gè)參數(shù)值，則所有待轉(zhuǎn)換點(diǎn)的坐標(biāo)轉(zhuǎn)換也就可以求得了。公式中含有7個(gè)轉(zhuǎn)換參數(shù)，為了求得7個(gè)轉(zhuǎn)換參數(shù)，至少需要3個(gè)公共點(diǎn)（也稱為重合點(diǎn)），當(dāng)多于3個(gè)公共點(diǎn)時(shí)，可按最小二乘法求得7個(gè)參數(shù)的最或是值。

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

式中，δ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)為等精度時(shí)，上述權(quán)矩陣為單位權(quán)矩陣，于是有：

單位權(quán)方差：

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

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

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

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

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

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

5）用可靠公共點(diǎn)重新計(jì)算坐標(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 問(wèn)題的分析處理

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

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

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

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

對(duì)上述警告性提示，筆者采用了2種方法處理。一是首先對(duì)坐標(biāo)進(jìn)行重心化處理，將所有點(diǎn)的坐標(biāo)平移到重心化坐標(biāo)系中，減小坐標(biāo)的絕對(duì)值。求逆完成后，在適當(dāng)處再反算還原到原坐標(biāo)系中。經(jīng)過(guò)這樣處理后，警告提示就再?zèng)]有出現(xiàn)了。另一種方法是采用Matlab中符號(hào)矩陣求逆運(yùn)算，因矩陣中變量元素是符號(hào)變量，即變量沒(méi)有賦值，避開(kāi)了數(shù)值矩陣求逆而進(jìn)行符號(hào)矩陣求逆，從而避免因數(shù)值計(jì)算誤差導(dǎo)致的矩陣接近奇異的情況發(fā)生，完成求逆運(yùn)算后再用相關(guān)函數(shù)將符號(hào)矩陣轉(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)建立一個(gè)文本文件保存，然后啟動(dòng)程序。程序設(shè)計(jì)上首先通過(guò)界面提示輸入?yún)⑴c求解轉(zhuǎn)換參數(shù)的重合點(diǎn)個(gè)數(shù)，并默認(rèn)重合定向點(diǎn)（即參與求解參數(shù)的重合點(diǎn)）為5，且最少要大于或等于3個(gè)，如輸入重合點(diǎn)數(shù)小于3，程序會(huì)提示輸入錯(cuò)誤。再通過(guò)界面分別選擇待轉(zhuǎn)換坐標(biāo)系所有點(diǎn)坐標(biāo)文件和目標(biāo)坐標(biāo)系公共點(diǎn)坐標(biāo)文件，程序?qū)?個(gè)文件中的空間坐標(biāo)數(shù)據(jù)分別保存到2個(gè)矩陣變量中。通過(guò)對(duì)坐標(biāo)系重心化處理及一系列運(yùn)算，計(jì)算出轉(zhuǎn)換坐標(biāo)、各點(diǎn)殘差及殘差中誤差以及單位權(quán)中誤差和7個(gè)轉(zhuǎn)換參數(shù)等，并將所有結(jié)果保存到Excel表格中。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

P226

B

1672-4623（2014）06-0062-03

10.3969/j.issn.1672-4623.2014.06.022

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

2014-01-26。