廣義整體最小二乘粗差探測的三維坐標轉換方法

申仲舒,項謙和,董思學,韋天赦

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

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

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

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

本文將三維坐標轉換數學模型抽象為非線性EIV模型,并對Amiri-Simkooei and Jazaeri[11]的思路進行推廣,將數據探測法拓展至該模型中。首先,利用Euler-Lagrange方法推導出非線性EIV模型的GTLS解。然后,將其等價轉化為經典最小二乘問題,分別在方差分量已知和未知方差的條件下,基于經典最小二乘理論,構造兩類數據探測的檢驗統計量。最后,采用試驗數據對三維基準轉換新算法的粗差探測性能進行驗證[7]。

1 非線性EIV模型的WTLS和數據探測算法

1.1 非線性EIV模型的GTLS解法

三維坐標轉換數學模型公式為

(1)

轉換公式可寫為非線性EIV模型的形式

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

(2)

式中,f和φ都是n×1維的抽象向量函數;a和L分別代表在原始和目標系統中m×1維和n×1維的觀測向量;e1和eL分別為a和L的隨機誤差向量;ξ為一個j×1維的未知參數向量。

相應的隨機模型為

(3)

其中

(4)

式中,Q為粗差向量e的正定協因數陣;QLL和Qaa分別為eL和e1的協因數陣;QLa和QaL為eL和e1之間的相關協因數陣;σ02為方差分量。上述模型的估計準則為

eTQ-1e=min

(5)

(6)

其中

(7)

A和B需要在迭代過程中不斷更新。以此構造拉格朗日目標函數為

(8)

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

(9)

(10)

(11)

式中,“～”和“＾ ”分別表示預測值和評估值。

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

(12)

(13)

(14)

將式(12)代入式(11)

(15)

其中

QT=GQGT=QLL+QLaBT+BQaL+BQaaBT

(16)

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

(17)

1.2 GTLS算法重構

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

l=Aδξ+el

(18)

將Ql作為l的協因數陣,兩邊對式(18)取方差

(19)

(20)

(21)

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

(22)

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

(23)

(24)

此外,估計參數的后驗方差因子和協方差矩陣可估計為

(25)

(26)

由于非線性EIV模型的GTLS估計是一個非線性問題,其參數估計及其精度評估是有偏差的。盡管在實踐中偏差可能很小,為了獲得偏差的影響,可以參考文獻[5]提供的方法,此方法已應用于標準EIV模型和約束TLS問題。

1.3 數據探測算法應用于非線性EIV模型

如果觀測數據受粗差污染,參數估值將不可靠。可在方差分量已知和未知的情況下,分別采用正態檢驗和t檢驗來檢測和消除粗差[12]。當方差分量已知時,第i個觀測方程的檢測統計量為

(27)

當方差分量未知時,第i個觀測方程的檢測統計量為

(28)

基于局部敏感性分析,可將這兩種檢驗統計量推導出來。在相關觀測數據的情況下,它們對檢測粗差更為敏感[11]。

在出現多個粗差的情況下,應持續實施探測過程,這個過程叫作“數據探測”。當單位權方差分量已知時,如果

(29)

或單位權方差分量未知時,

(30)

即消除觀測方程。如果滿足以上等式,則證明相應的觀測方程中存在粗差。粗差可能在觀測向量L或向量a中,甚至在兩者中。應將標記的觀測方程從方程表中刪除,然后調整剩余的數據,重復執行此過程,直至所有統計檢驗都被接受。

2 數據探測算法在三維坐標轉換中的應用

在三維坐標轉換算法的應用中,需要確定矩陣A和B的具體表達式。

假設共有k個公共點,則相應維數分別為m=n=3k和t=7。

對于第i個公共點

(31)

(32)

其中

(33)

(34)

=-μ0·M0

(35)

(36)

(37)

通過將每個點的相應表達式組合在一起,得到矩陣A和B為

(38)

3 試驗分析

假設有36個點分布在一個空間域中。已知源坐標系(a系統)和目標坐標系(b系統)中每個點的實際坐標。從系統a轉換到系統b的7個參數如下:Δx=1000, Δy=1000, Δz=1000,μ=1.5,β1=1.0 rad,β2=1.5 rad,β3=-0.5 rad。

本文將26個點作為公共點,其余10個點是檢核點。根據已知的協方差矩陣,每次模擬均會產生相關隨機誤差。在含有隨機誤差的坐標中引入5個粗差,其大小介于(-30,-5)和(5,30)倍的先驗標準差之間,位置隨機產生,反復進行500次模擬。

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

(39)

總中誤差為

(40)

上述方案的統計結果見表1,包含σx、σy、σz和σp的平均值和最大值,此外,圖1為方案2和方案3的試驗對比結果。

表1 不同轉換方案中誤差統計 m

圖1 轉換方案2和3的試驗結果對比

(1)當粗差引入前,GTLS估計(方案1)得到了3種方案之間的最優轉換參數。

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

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

4 結 語

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

試驗及其結果表明,該算法能有效削弱粗差的影響,獲得可靠的轉換參數。與WTLS估計相比,所得參數的精度和準確度有了顯著提高。