顧及2 套坐標(biāo)誤差的三維坐標(biāo)變換方法

李微曉 ,沈云中 ,李博峰

(1.同濟(jì)大學(xué) 測(cè)量與國(guó)土信息工程系, 上海200092;2.淮海工學(xué)院 測(cè)繪工程學(xué)院, 江蘇連云港222001;3.現(xiàn)代工程測(cè)量國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室, 上海200092)

運(yùn)用全球定位系統(tǒng)(GPS)技術(shù)建立控制網(wǎng)涉及到三維坐標(biāo)變換模型.三維坐標(biāo)轉(zhuǎn)換通常采用7 個(gè)轉(zhuǎn)換參數(shù), 即3 個(gè)平移參數(shù)、3 個(gè)旋轉(zhuǎn)參數(shù)和1 個(gè)尺度參數(shù).根據(jù)旋轉(zhuǎn)與尺度參數(shù)參考點(diǎn)的不同定義,三維坐標(biāo)轉(zhuǎn)換模型可分為:Bursa-Wolf 模型(簡(jiǎn)稱Bursa 模型)、Molodensky 模型及武測(cè)模型等[1].然而,利用最小二乘方法求解這些模型時(shí)通常只采用1組公共點(diǎn)坐標(biāo)改正數(shù)向量, 而不是對(duì)2 套坐標(biāo)的公共點(diǎn)誤差分別引入改正數(shù)進(jìn)行解算,因此, 求得的改正數(shù)向量不能用于對(duì)任何一套坐標(biāo)進(jìn)行誤差改正,以免引起公共點(diǎn)與待轉(zhuǎn)換點(diǎn)坐標(biāo)之間的縫隙, 除非對(duì)待轉(zhuǎn)換點(diǎn)坐標(biāo)的誤差也作相應(yīng)的改正.事實(shí)上,若能同時(shí)考慮2 套坐標(biāo)的誤差, 在求解轉(zhuǎn)換參數(shù)時(shí),求得2 套坐標(biāo)的改正數(shù)并進(jìn)行改正, 并用最小二乘配置法(least-squares collocation,LSC)推估待轉(zhuǎn)換點(diǎn)的坐標(biāo)改正數(shù)并進(jìn)行改正, 有望提高轉(zhuǎn)換坐標(biāo)的精度.LSC 由K rarup 于1969 年在研究重力異常時(shí)提出[2] .1973 年奧地利的Moritz 進(jìn)行了深入研究, 提出帶系統(tǒng)參數(shù)的LSC[3].從20 世紀(jì)80 年代以來(lái),我國(guó)學(xué)者對(duì)該理論進(jìn)行了廣泛深入的研究.周江文提出了LSC 模型的兩步解法[4];楊元喜等討論了基于方差分量估計(jì)的自適應(yīng)LSC 模型,并將其應(yīng)用于高程系統(tǒng)轉(zhuǎn)換和地理信息系統(tǒng)(GIS)誤差糾正[5];李博峰指出LSC 模型的理論基礎(chǔ)是極大驗(yàn)后估計(jì),實(shí)質(zhì)是采用觀測(cè)值對(duì)隨機(jī)參數(shù)先驗(yàn)信息的改造,其改造程度由觀測(cè)值的方差和觀測(cè)值與參數(shù)的相關(guān)程度共同決定[6].

本文以Bursa 模型為例研究顧及2 套坐標(biāo)誤差的三維坐標(biāo)轉(zhuǎn)換, 用2 套坐標(biāo)改正數(shù)的加權(quán)平方和最小為準(zhǔn)則導(dǎo)出了轉(zhuǎn)換參數(shù)的計(jì)算公式;并根據(jù)公共點(diǎn)的坐標(biāo)改正數(shù)以及公共點(diǎn)與轉(zhuǎn)換點(diǎn)的協(xié)方差陣,推估轉(zhuǎn)換點(diǎn)的坐標(biāo)改正數(shù)并對(duì)其改正, 從而求得近似無(wú)縫的坐標(biāo)轉(zhuǎn)換結(jié)果.采用模擬數(shù)據(jù)分析了新方法的有效性, 與直接將非公共點(diǎn)轉(zhuǎn)換到新坐標(biāo)系的傳統(tǒng)方法相比, 新方法能有效地提高坐標(biāo)轉(zhuǎn)換精度.

1 七參數(shù)坐標(biāo)轉(zhuǎn)換模型

Bursa-Wolf 模型以坐標(biāo)原點(diǎn)為旋轉(zhuǎn)和尺度的參考點(diǎn),對(duì)應(yīng)的7 個(gè)參數(shù)轉(zhuǎn)換模型為

式中:dx,dy,dz為3 個(gè)平移參數(shù);εx,εy,εz為3 個(gè)旋轉(zhuǎn)參數(shù);Δu為尺度參數(shù);下標(biāo)Ⅰ與Ⅱ表示2 套坐標(biāo)系;i表示第i個(gè)公共點(diǎn).假設(shè)n個(gè)公共點(diǎn)在2 套坐標(biāo)系下對(duì)應(yīng)的協(xié)方差陣分別為Σ11,Ⅰ和Σ11,Ⅱ, 傳統(tǒng)的求解方法只考慮公共點(diǎn)在第1 套坐標(biāo)系下的誤差,對(duì)應(yīng)的轉(zhuǎn)換參數(shù)的最小二乘解為

2 顧及2 套坐標(biāo)誤差的坐標(biāo)轉(zhuǎn)換模型

觀測(cè)向量L 等于第2 套坐標(biāo)系公共點(diǎn)的坐標(biāo)向量LⅡ與對(duì)應(yīng)第1 套坐標(biāo)系坐標(biāo)向量LⅠ之差,即

采用最小二乘求解,同時(shí)考慮觀測(cè)向量L 中公共點(diǎn)在2 套坐標(biāo)系下的誤差.令vⅡ, vⅠ分別為L(zhǎng)Ⅱ,LⅠ的改正數(shù)向量,有

誤差方程式可表示為

顯然, 觀測(cè)向量LⅡ與LⅠ相互獨(dú)立, 以2 套公共點(diǎn)坐標(biāo)改正數(shù)的加權(quán)平方和最小為準(zhǔn)則

根據(jù)拉格朗日條件極值原理,構(gòu)造極值函數(shù)

式中:K為聯(lián)系數(shù)向量, ξ為7 參數(shù)向量, 分別對(duì)殘差向量vⅠ, vⅡ和7 參數(shù)向量ξ求一階導(dǎo)數(shù), 并令其為零, 得

由式(8)中前2 式得

兩式相減, 得

顧及式(4)和式(5), 則

將式(11)代入式(8)中第3 式, 得

因此, 考慮二套坐標(biāo)誤差,對(duì)應(yīng)的轉(zhuǎn)換參數(shù)的最小二乘解ξ^N為

求得轉(zhuǎn)換參數(shù)后, 由式(10)計(jì)算K,再代入式(9)計(jì)算觀測(cè)值改正數(shù).

3 轉(zhuǎn)換點(diǎn)的坐標(biāo)改正數(shù)算法

4 算例與分析

4.1 算例設(shè)計(jì)和協(xié)方差的構(gòu)造

算例旨在驗(yàn)證本文坐標(biāo)轉(zhuǎn)換方法能否提高坐標(biāo)轉(zhuǎn)換精度.在經(jīng)度110°～130°、緯度20°～40°范圍內(nèi)設(shè)計(jì)2 個(gè)方案.方案1 :每隔4°取一個(gè)點(diǎn), 共6 ×6 個(gè)點(diǎn),見(jiàn)圖1a ;方案2 :每隔2°取一個(gè)點(diǎn), 共11 ×11 個(gè)點(diǎn),見(jiàn)圖1b .

圖1 選點(diǎn)方案Fig .1 Point distribution in different control networks

根據(jù)圖1 中各點(diǎn)的經(jīng)緯度,設(shè)其大地高為0 m ,采用WGS84 橢球參數(shù)計(jì)算每個(gè)點(diǎn)的空間三維坐標(biāo)作為近似坐標(biāo).

在2 種方案中, 第1 套坐標(biāo)系中的所有點(diǎn)的協(xié)方差陣根據(jù)三維空間距離觀測(cè)值的間接平差計(jì)算得到,測(cè)距精度模擬為10-2+10-8D,D是2 點(diǎn)之間的三維空間距離, 以米為單位.列出點(diǎn)位坐標(biāo)的估值與觀測(cè)值L 的改正數(shù)v 之間的函數(shù)關(guān)系為

其中,B 為誤差方程的系數(shù)陣.則X 的估值協(xié)因數(shù)陣為

其中,Nbb=BTPB ,P 為權(quán)陣.協(xié)方差陣為

第2 套坐標(biāo)系中全部點(diǎn)的協(xié)方差陣用GPS 網(wǎng)平差的方式計(jì)算.首先根據(jù)模擬點(diǎn)坐標(biāo)計(jì)算出任意2點(diǎn)的基線, 給基線模擬標(biāo)準(zhǔn)差為10-2+10-8D的正態(tài)分布隨機(jī)誤差, 然后列誤差方程式并附加全部點(diǎn)坐標(biāo)改正數(shù)平方和最小為約束條件, 以位于選點(diǎn)范圍中部的公共點(diǎn)為基準(zhǔn)點(diǎn)進(jìn)行平差計(jì)算, 其中, 3 個(gè)基線分量的權(quán)取P x=(Δx)2/D2,P y=(Δy)2/D2,P z=(Δz)2/D2.

4 .2 解算與分析

方案1 :根據(jù)4 .1 所述, 對(duì)圖1a 中的所有36 個(gè)點(diǎn)分別按照距離觀測(cè)和按GPS 網(wǎng)平差觀測(cè)進(jìn)行,對(duì)應(yīng)得到全部點(diǎn)在第1 套和第2 套坐標(biāo)系下的協(xié)方差陣.平差后考慮網(wǎng)中與外圍點(diǎn)相聯(lián)的觀測(cè)數(shù)較少,點(diǎn)位誤差遠(yuǎn)大于內(nèi)部點(diǎn)位誤差, 如果對(duì)全部點(diǎn)位誤差統(tǒng)計(jì), 則統(tǒng)計(jì)結(jié)果不能客觀反映內(nèi)部點(diǎn)位精度,故僅取16 個(gè)中間點(diǎn)比較,其中5 個(gè)是公共點(diǎn)、11 個(gè)是轉(zhuǎn)換點(diǎn).將5 個(gè)公共點(diǎn)按式(13)計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù)后,按式(9)計(jì)算第1 套坐標(biāo)系和第2 套坐標(biāo)系中公共點(diǎn)的坐標(biāo)改正數(shù)vⅠ, vⅡ,再分別用式(14)和(15)計(jì)算第1 套坐標(biāo)系和第2 套坐標(biāo)系中11 個(gè)轉(zhuǎn)換點(diǎn)的坐標(biāo)改正數(shù)v2Ⅰ, v2Ⅱ, 并改正這些轉(zhuǎn)換點(diǎn)坐標(biāo)(記改正后的第1 套和第2 套轉(zhuǎn)換點(diǎn)坐標(biāo)分別為～X2,Ⅰ和～X2,Ⅱ);然后用轉(zhuǎn)換參數(shù)ξ轉(zhuǎn)換改正后的第1 套轉(zhuǎn)換點(diǎn)坐標(biāo)～X2,Ⅰ(記為～X2,Ⅱ),并與改正后的第2 套轉(zhuǎn)換點(diǎn)坐標(biāo)～X2,Ⅱ比較,按式(19)計(jì)算11 個(gè)轉(zhuǎn)換點(diǎn)的坐標(biāo)精度

以及按(20)計(jì)算轉(zhuǎn)換的點(diǎn)位精度

該過(guò)程稱為L(zhǎng)SC 法.傳統(tǒng)方法直接用式(1)計(jì)算轉(zhuǎn)換參數(shù)ξ^L并轉(zhuǎn)換11 個(gè)轉(zhuǎn)換點(diǎn), 并與X2,Ⅱ比較按式(19)和式(20)計(jì)算轉(zhuǎn)換精度, 該過(guò)程稱為L(zhǎng)S 法.為了在統(tǒng)計(jì)意義上評(píng)價(jià)LSC 方法的優(yōu)越性,按4 .1 中的協(xié)方差模擬過(guò)程計(jì)算1 000 次,并按式(20)分別計(jì)算LSC 和LS 坐標(biāo)轉(zhuǎn)換點(diǎn)位精度,如圖2 所示.顯然,LSC 的精度比LS 高.類似于方案1 的實(shí)驗(yàn)過(guò)程, 對(duì)方案2 作1 000 次計(jì)算,結(jié)果如圖3 所示.

圖2 6×6 網(wǎng)中LSC 法與直接轉(zhuǎn)換的誤差比較Fig .2 Comparison of errors between LSC method and the direct conversion method for a 6×6 grid

圖3 11×11 網(wǎng)LSC 法與直接轉(zhuǎn)換的誤差比較Fig .3 Comparison of errors between the LSC method and thedirect conversion method for a 11×11 grid

對(duì)2 種方案的1 000 次試驗(yàn)坐標(biāo)轉(zhuǎn)換點(diǎn)位精度和3 個(gè)坐標(biāo)分量精度統(tǒng)計(jì)如表1 所示.表2 給出了精度提高百分比.結(jié)果表明:同一范圍方案1 選36個(gè)點(diǎn),方案2 選121 個(gè)點(diǎn),當(dāng)控制點(diǎn)不變, 其余點(diǎn)的點(diǎn)位間距越小, 點(diǎn)與點(diǎn)之間坐標(biāo)相關(guān)性越大, 考慮坐標(biāo)改正的LSC 方法較LS 方法能有效提高坐標(biāo)轉(zhuǎn)換精度.

表1 2 種方案的轉(zhuǎn)換精度比較Tab.1 Comparison of conversion accuracy in the two schemes cm

表2 LSC 相對(duì)于LS 的精度提高率Tab.2 Rate of increased relativeaccuracy of LSC to LS

5 結(jié)論

研究了顧及2 套坐標(biāo)誤差的三維坐標(biāo)轉(zhuǎn)換模型,以2 套坐標(biāo)改正數(shù)的加權(quán)平方和最小為準(zhǔn)則導(dǎo)出了轉(zhuǎn)換參數(shù)的計(jì)算公式;用推估方法根據(jù)公共點(diǎn)與轉(zhuǎn)換點(diǎn)的協(xié)方差陣計(jì)算轉(zhuǎn)換點(diǎn)的坐標(biāo)改正數(shù)并對(duì)其改正,求得了無(wú)縫的坐標(biāo)轉(zhuǎn)換結(jié)果.得到結(jié)論:

(1)坐標(biāo)轉(zhuǎn)換的2 期坐標(biāo)通常來(lái)自不同測(cè)量手段,若已知其協(xié)方差陣, 用LSC 法先對(duì)轉(zhuǎn)換點(diǎn)坐標(biāo)修正再進(jìn)行坐標(biāo)轉(zhuǎn)換能有效提高轉(zhuǎn)換精度.

(2)當(dāng)點(diǎn)的坐標(biāo)之間相關(guān)性較大, 直接將轉(zhuǎn)換點(diǎn)的坐標(biāo)換算成新坐標(biāo)會(huì)產(chǎn)生較大的誤差, 采用LSC 法改正,能很好地解決殘差問(wèn)題.

(3)采用本文模型求解轉(zhuǎn)換參數(shù)并用LSC 對(duì)轉(zhuǎn)換點(diǎn)改正, 能實(shí)現(xiàn)整體、無(wú)縫的坐標(biāo)轉(zhuǎn)換結(jié)果.

[1] 劉大杰, 施一民, 過(guò)靜珺.全球定位系統(tǒng)(GPS)的原理與數(shù)據(jù)處理[M] .上海:同濟(jì)大學(xué)出版社, 1996.LI U Dajie,SH I Yimin, GUO Jingjun.The theory and data processing of global position system (GPS)[M] .Shanghai :Tongji University Press, 1996.

[2] Krarup T .Acontribution to the mathematical foundation of phy sical geodesy [R] . C openhagen : Danish Geodetic Institude,1969.

[3] Mo ritz H .Least squares collocation[M] .Munich:Deutsche Geodaetische Kommission(DGK), 1973.

[4] 周江文.擬合推估新解——兩步解法[J] .測(cè)繪學(xué)報(bào), 2002, 31(2):189.ZHO U Jiangw en.Atwo-step solution of collocation[J] .Acta Geodaetica et Cartographica Sinica,2002, 31(2):189.

[5] 楊元喜, 張菊清, 張亮.基于方差分量估計(jì)的擬合推估及其在GIS 誤差糾正的應(yīng)用[J] .測(cè)繪學(xué)報(bào), 2008, 37(1):152.YANG Yuanxi, ZH ANG Juqing,ZHANG Liang .Estimation based collocation and its in GIS error variance fitting component application[J] .Acta Geodaetica et Cartographica Sinica,2008, 37(1):152.

[6] 李博峰.混合整數(shù)GNSS 隨機(jī)模型及函數(shù)模型參數(shù)估計(jì)理論與方法[D] .上海:同濟(jì)大學(xué), 2010.LI Bofeng .Theory and method of parameter estimation for GNSS stochastic and function models[D] .Shanghai:Tongji University,2010.

[7] Schaffrin B .T LS-collocation:The total least-squares approach to EIV-models with prior information [C/ CD] ∥18th International Workshop on Matrices and Statistics.Smolenice Castle :Institute of Measurement Science, 2009.