申仲舒,項(xiàng)謙和,董思學(xué),韋天赦

(1.浙江省測繪科學(xué)技術(shù)研究院,浙江 杭州 311122; 2.中國地質(zhì)大學(xué)(武漢),湖北 武漢 430074;3. 中國煤炭地質(zhì)總局浙江煤炭地質(zhì)局,浙江 杭州 310021)

在大地測量、攝影測量和計(jì)算機(jī)視覺等領(lǐng)域,三維坐標(biāo)變換是將觀測數(shù)據(jù)從原始坐標(biāo)系轉(zhuǎn)換至目標(biāo)系統(tǒng)的常用方法。三維坐標(biāo)轉(zhuǎn)換最重要的步驟是確定7個(gè)轉(zhuǎn)換參數(shù)(1個(gè)尺度比參數(shù)、3個(gè)平移參數(shù)和3個(gè)旋轉(zhuǎn)角參數(shù))。當(dāng)旋轉(zhuǎn)角足夠小且尺度比接近1時(shí),三維坐標(biāo)轉(zhuǎn)換數(shù)學(xué)模型可以簡化為著名的布爾莎模型,通常用線性Gauss-Markov(GM)模型來描述。然而,實(shí)際中旋轉(zhuǎn)參數(shù)和尺度參數(shù)可能是任意的,因此,一些研究采用最小二乘法求解非線性系統(tǒng)[1-2]。

由于實(shí)際中原始坐標(biāo)通常也是觀測所得,矩陣系數(shù)中的某些元素不可避免地會(huì)受隨機(jī)誤差的影響。因此,轉(zhuǎn)換模型不能用GM模型描述,此時(shí)應(yīng)引入變量誤差(EIV)模型[3]。EIV模型的參數(shù)估計(jì)方法首先由文獻(xiàn)[4]在數(shù)值分析領(lǐng)域提出,并命名為TLS,此后,又提出了許多算法,尤其是WTLS算法[5-6]。從應(yīng)用的角度看,TLS的研究已經(jīng)被引用到許多實(shí)際的問題中,如直線擬合、坐標(biāo)轉(zhuǎn)換、GNSS單點(diǎn)定位、遙感影像配準(zhǔn)。近年來,TLS方法被用來解決三維坐標(biāo)轉(zhuǎn)換問題(任意旋轉(zhuǎn)參數(shù)和尺度參數(shù))已經(jīng)成了研究熱點(diǎn)[2,7]。

上述的TLS算法只有在原始坐標(biāo)和目標(biāo)坐標(biāo)中僅包含偶然誤差時(shí)適用。然而,如果觀測值受到粗差污染,用這些方法獲取的轉(zhuǎn)換參數(shù)也會(huì)受到不利影響甚至嚴(yán)重失真。目前,對(duì)于TLS解法的粗差處理問題也有一些研究[1]。然而,由于標(biāo)準(zhǔn)EIV模型形式上的限制,所有的這些方法在大旋轉(zhuǎn)角和任意尺度比的三維坐標(biāo)轉(zhuǎn)換問題中不適用[8-10]。

本文將三維坐標(biāo)轉(zhuǎn)換數(shù)學(xué)模型抽象為非線性EIV模型,并對(duì)Amiri-Simkooei and Jazaeri[11]的思路進(jìn)行推廣,將數(shù)據(jù)探測法拓展至該模型中。首先,利用Euler-Lagrange方法推導(dǎo)出非線性EIV模型的GTLS解。然后,將其等價(jià)轉(zhuǎn)化為經(jīng)典最小二乘問題,分別在方差分量已知和未知方差的條件下,基于經(jīng)典最小二乘理論,構(gòu)造兩類數(shù)據(jù)探測的檢驗(yàn)統(tǒng)計(jì)量。最后,采用試驗(yàn)數(shù)據(jù)對(duì)三維基準(zhǔn)轉(zhuǎn)換新算法的粗差探測性能進(jìn)行驗(yàn)證[7]。

1 非線性EIV模型的WTLS和數(shù)據(jù)探測算法

1.1 非線性EIV模型的GTLS解法

三維坐標(biāo)轉(zhuǎn)換數(shù)學(xué)模型公式為

(1)

轉(zhuǎn)換公式可寫為非線性EIV模型的形式

L-eL=f(a-e1,ξ)=φ(e1,ξ)

(2)

式中,f和φ都是n×1維的抽象向量函數(shù);a和L分別代表在原始和目標(biāo)系統(tǒng)中m×1維和n×1維的觀測向量;e1和eL分別為a和L的隨機(jī)誤差向量;ξ為一個(gè)j×1維的未知參數(shù)向量。

相應(yīng)的隨機(jī)模型為

(3)

其中

(4)

式中,Q為粗差向量e的正定協(xié)因數(shù)陣;QLL和Qaa分別為eL和e1的協(xié)因數(shù)陣;QLa和QaL為eL和e1之間的相關(guān)協(xié)因數(shù)陣;σ02為方差分量。上述模型的估計(jì)準(zhǔn)則為

eTQ-1e=min

(5)

(6)

其中

(7)

A和B需要在迭代過程中不斷更新。以此構(gòu)造拉格朗日目標(biāo)函數(shù)為

(8)

式中,λ為拉格朗日乘子的n×1維向量。利用Euler-Lagrange必要條件可解出

(9)

(10)

(11)

式中,“～”和“＾ ”分別表示預(yù)測值和評(píng)估值。

從式(9)可以得到誤差向量

(12)

(13)

(14)

將式(12)代入式(11)

(15)

其中

QT=GQGT=QLL+QLaBT+BQaL+BQaaBT

(16)

將式(15)代入式(10),可以獲得δξ的表達(dá)式

(17)

1.2 GTLS算法重構(gòu)

式(6)可以寫成如下形式

l=Aδξ+el

(18)

將Ql作為l的協(xié)因數(shù)陣,兩邊對(duì)式(18)取方差

(19)

(20)

(21)

基于上述證明,非線性EIV模型可以利用經(jīng)典最小二乘理論知識(shí)。誤差向量el估計(jì)為

(22)

因此式(13)、式(14)可寫成如下形式

(23)

(24)

此外,估計(jì)參數(shù)的后驗(yàn)方差因子和協(xié)方差矩陣可估計(jì)為

(25)

(26)

由于非線性EIV模型的GTLS估計(jì)是一個(gè)非線性問題,其參數(shù)估計(jì)及其精度評(píng)估是有偏差的。盡管在實(shí)踐中偏差可能很小,為了獲得偏差的影響,可以參考文獻(xiàn)[5]提供的方法,此方法已應(yīng)用于標(biāo)準(zhǔn)EIV模型和約束TLS問題。

1.3 數(shù)據(jù)探測算法應(yīng)用于非線性EIV模型

如果觀測數(shù)據(jù)受粗差污染,參數(shù)估值將不可靠?？稍诜讲罘至恳阎臀粗那闆r下,分別采用正態(tài)檢驗(yàn)和t檢驗(yàn)來檢測和消除粗差[12]。當(dāng)方差分量已知時(shí),第i個(gè)觀測方程的檢測統(tǒng)計(jì)量為

(27)

當(dāng)方差分量未知時(shí),第i個(gè)觀測方程的檢測統(tǒng)計(jì)量為

(28)

基于局部敏感性分析,可將這兩種檢驗(yàn)統(tǒng)計(jì)量推導(dǎo)出來。在相關(guān)觀測數(shù)據(jù)的情況下,它們對(duì)檢測粗差更為敏感[11]。

在出現(xiàn)多個(gè)粗差的情況下,應(yīng)持續(xù)實(shí)施探測過程,這個(gè)過程叫作“數(shù)據(jù)探測”。當(dāng)單位權(quán)方差分量已知時(shí),如果

(29)

或單位權(quán)方差分量未知時(shí),

(30)

即消除觀測方程。如果滿足以上等式,則證明相應(yīng)的觀測方程中存在粗差。粗差可能在觀測向量L或向量a中,甚至在兩者中。應(yīng)將標(biāo)記的觀測方程從方程表中刪除,然后調(diào)整剩余的數(shù)據(jù),重復(fù)執(zhí)行此過程,直至所有統(tǒng)計(jì)檢驗(yàn)都被接受。

2 數(shù)據(jù)探測算法在三維坐標(biāo)轉(zhuǎn)換中的應(yīng)用

在三維坐標(biāo)轉(zhuǎn)換算法的應(yīng)用中,需要確定矩陣A和B的具體表達(dá)式。

假設(shè)共有k個(gè)公共點(diǎn),則相應(yīng)維數(shù)分別為m=n=3k和t=7。

對(duì)于第i個(gè)公共點(diǎn)

(31)

(32)

其中

(33)

(34)

=-μ0·M0

(35)

(36)

(37)

通過將每個(gè)點(diǎn)的相應(yīng)表達(dá)式組合在一起,得到矩陣A和B為

(38)

3 試驗(yàn)分析

假設(shè)有36個(gè)點(diǎn)分布在一個(gè)空間域中。已知源坐標(biāo)系(a系統(tǒng))和目標(biāo)坐標(biāo)系(b系統(tǒng))中每個(gè)點(diǎn)的實(shí)際坐標(biāo)。從系統(tǒng)a轉(zhuǎn)換到系統(tǒng)b的7個(gè)參數(shù)如下:Δx=1000, Δy=1000, Δz=1000,μ=1.5,β1=1.0 rad,β2=1.5 rad,β3=-0.5 rad。

本文將26個(gè)點(diǎn)作為公共點(diǎn),其余10個(gè)點(diǎn)是檢核點(diǎn)。根據(jù)已知的協(xié)方差矩陣,每次模擬均會(huì)產(chǎn)生相關(guān)隨機(jī)誤差。在含有隨機(jī)誤差的坐標(biāo)中引入5個(gè)粗差,其大小介于(-30,-5)和(5,30)倍的先驗(yàn)標(biāo)準(zhǔn)差之間,位置隨機(jī)產(chǎn)生,反復(fù)進(jìn)行500次模擬。

設(shè)計(jì)以下3種計(jì)算方案來求解轉(zhuǎn)換參數(shù):①在引入粗差之前,對(duì)非線性EIV模型進(jìn)行了GTLS估計(jì),并在引入粗差之后進(jìn)行了以下兩種方案;②非線性EIV模型的GTLS解法;③非線性EIV模型的數(shù)據(jù)探測算法,由于方差分量已知,本方案中采用正態(tài)分布檢驗(yàn)(α=0.05)。

(39)

總中誤差為

(40)

上述方案的統(tǒng)計(jì)結(jié)果見表1,包含σx、σy、σz和σp的平均值和最大值,此外,圖1為方案2和方案3的試驗(yàn)對(duì)比結(jié)果。

表1 不同轉(zhuǎn)換方案中誤差統(tǒng)計(jì) m

圖1 轉(zhuǎn)換方案2和3的試驗(yàn)結(jié)果對(duì)比

(1)當(dāng)粗差引入前,GTLS估計(jì)(方案1)得到了3種方案之間的最優(yōu)轉(zhuǎn)換參數(shù)。

(2)引入粗差后,GTLS(方案2)得到的轉(zhuǎn)換參數(shù)嚴(yán)重失真,因此,在模擬過程中,方案2中存在許多大的均方根誤差,此外,σx、σy、σz和σp的對(duì)應(yīng)平均值均增加了約3倍。

(3)與GTLS算法相比,非線性EIV模型的數(shù)據(jù)探測方法(方案3)能有效地提高轉(zhuǎn)換坐標(biāo)的精度。所有方向的最大中誤差明顯小于方案2的最大中誤差。還可看出,方案3的σx、σy、σz和σp的對(duì)應(yīng)平均值分別僅為方案2的39.4%、36.9%、41.6%和38.4%。盡管方案3的結(jié)果不如方案1,但它們之間的差異并不顯著。因此,非線性EIV模型的數(shù)據(jù)探測方法得到的轉(zhuǎn)換參數(shù)仍然是十分可靠的[13-14]。

4 結(jié) 語

為了處理嚴(yán)密三維坐標(biāo)轉(zhuǎn)換的粗差處理問題,本文將數(shù)據(jù)探測算法應(yīng)用于非線性EIV模型中。首先用Euler-Lagrange方法推導(dǎo)出了非線性EIV模型的GTLS迭代解,然后用經(jīng)典最小二乘法對(duì)其進(jìn)行了重構(gòu)。分別在已知和未知方差分量因子的條件下,基于經(jīng)典最小二乘理論,構(gòu)造了數(shù)據(jù)探測的兩個(gè)檢驗(yàn)統(tǒng)計(jì)量。這里需要指出該算法可以看作是Amiri Simkooei和Jazaeri算法的一種推廣形式[8,11,15]。

試驗(yàn)及其結(jié)果表明,該算法能有效削弱粗差的影響,獲得可靠的轉(zhuǎn)換參數(shù)。與WTLS估計(jì)相比,所得參數(shù)的精度和準(zhǔn)確度有了顯著提高。