關于坐標轉換尺度系數計算的分析

秦 鋒 張振虎

(1. 中國核工業集團三門核電有限公司 浙江 臺州 317112;2. 上海勘察設計研究院(集團)有限公司 上海 200093)

0 引言

空間直角坐標轉換需要計算七參數,包括尺度系數、旋轉角度(旋轉矩陣)及平移參數。目前坐標轉換的算法主要分為兩類,一類是基于經典最小二乘(least square,LS)的坐標轉換算法,另一類是基于變量含誤差(errors in variables,EIV)模型的整體最小二乘(total least squares,TLS)或加權整體最小二乘(weighted total least squares,WTLS)坐標轉換算法。

尺度系數通常認為是兩坐標系距離之比。在經典最小二乘準則下,小角度坐標轉換一般采用線性最小二乘一次性計算出包括尺度系數在內的坐標轉換七參數。對于基于經典最小二乘準則的大角度坐標轉換,文獻[1]論證了尺度系數和旋轉矩陣(或旋轉角度)之間不相關,和平移參數之間呈強相關,故大部分基于經典最小二乘準則的大角度坐標轉換算法在計算坐標轉換參數時先計算尺度系數,再計算旋轉矩陣和平移參數。如文獻[2]通過計算坐標轉換前后所有距離比值的平均值來計算尺度系數,文獻[3]對坐標數據重心化處理后采用最大似然估計計算尺度系數,文獻[4]對坐標數據重心化處理后通過計算坐標轉換前后重心化坐標平方和比值的開方值來計算尺度系數。另有部分基于經典最小二乘準則的大角度坐標轉換算法一次性計算出坐標轉換七參數,如文獻[5]對誤差方程線性化處理后通過迭代一次性求解坐標轉換七參數,不需要分別計算尺度系數、旋轉矩陣及平移參數。但在實際應用中發現,以上計算坐標轉換尺度系數的方法計算結果互有差異,難以確定何種方法尺度系數計算結果最優。

在TLS或WTLS準則下,坐標轉換參數中尺度系數的求解較LS準則下復雜。目前已有的TLS或WTLS坐標轉換方法主要分為兩類,一類是分別計算尺度系數、旋轉矩陣及平移參數,另一類通過迭代一次性計算出全部參數。第一類方法中,如文獻[6]采用多元總體最小二乘方法分別求解坐標轉換參數,文獻[7]推導了特定加權矩陣下采用多元總體最小二乘求解坐標轉換參數的方法。第二類方法中,如文獻[8]采用提出了一種基于TLS準則的小角度坐標轉換算法,文獻[9]采用改進的加權整體最小二乘算法用于小角度坐標轉換,文獻[10]提出了一種通用的加權整體最小二乘坐標轉換方法。

本文通過奇異值分解算法(singular value decomposition,SVD)推導了LS及TLS準則下坐標轉換尺度系數的計算公式,并結合工程實例,與其他尺度系數計算方法進行了對比。同時,為驗證本文算法,采用了布羅伊登-弗萊徹-戈德法布-香農(Broyden-Fletcher-Goldfarb-Shanno,BFGS)優化算法進行結果驗證。對于WTLS準則下坐標轉換尺度系數的計算及其與旋轉角度、平移參數之間的相關性,本文基于文獻[10]中的方法進行了數據驗證,并得出了相應結論。

1 數學模型及算法

坐標轉換模型可表示為

為方便參數計算,根據文獻[11]對坐標進行重心化處理。數據重心化處理后,旋轉矩陣和平移參數可分開進行計算,其計算公式為

分別定義矩陣A和B為

重心化后的坐標轉換模型就可表示為:A=λRB

1.1 基于LS準則的坐標轉換尺度系數計算

求取正交矩陣R,使A-λRB的弗羅貝紐斯(Frobenius)范數最小時所得的R即為最佳轉換矩陣,即求解以下問題，即

(5)

(6)

式中,tr是trace的簡寫,表示矩陣的跡。

因λ>0,故在tr(BATR)取最大值時,式(6)得到最小值。

對BAT進行奇異值分解,得

(7)

式中,U是左奇異矩陣;V是右奇異矩陣;Σ=diag(σ1,σ2,σ3),σ1,σ2,σ3為BAT的奇異值。根據文獻[12],R=VUT時,tr(BATR)取最大值。將求得的R代入式(6),此時式(6)為關于λ的一元二次方程求極值問題,當λ為式(6)時,式(6)得到最小值。

(8)

平移參數可根據式(4)進行計算。

1.2 基于TLS準則的坐標轉換尺度系數計算

考慮到矩陣A、B中的隨機誤差EA、EB,構建EIV模型為

(9)

假設隨機誤差EA、EB獨立同分布,基于拉格朗日乘子法,構建整體最小二乘優化函數為

(10)

函數f分別對EA、EB、μ求偏導并置0,解出EA、EB、μ后代入式(10),得

(11)

同式(7),對BAT進行奇異值分解;同理,當且僅當R=VUT時,f取極小值。將R=VUT代入式(11)中,可得

(12)

式(2)對λ求偏導并置0,因λ>0,tr(Σ)>0,故可解出

(13)

平移參數可根據式(4)進行計算。

1.3 驗證算法

分別定義以下矩陣:

其中,Ln表示元素全是1的列向量。

坐標轉換模型又可表示為:

(14)

構建基于LS準則的坐標轉換優化函數:

(15)

參照式(9)至式(11),構建基于TLS準則的坐標轉換優化函數為

(16)

對上述兩個函數分別采用BFGS算法進行優化,迭代求取坐標轉換七參數。具體BFGS算法可參照文獻[13]。

2 案例分析

坐標轉換數據選用文獻[10]中的數據,詳細數據見表1。

表1 坐標轉換原始數據 單位：m

分別采用文獻[2-5]中算法計算尺度系數,同時采用本文算法及驗證算法分別計算LS準則下及TLS準則下坐標轉換七參數,解算結果見表2。

表2 LS及TLS準則下坐標轉換參數解算結果

從以上計算結果可以看出:

(1)本文LS坐標轉換算法與基于BFGS優化算法的LS坐標轉換驗證算法解算的七參數完全一致,而文獻[2-5]尺度系數計算結果均與本文LS算法結果有差異,說明LS準則下本文提出的尺度系數計算方法更準確。

(2)在TLS準則下,本文算法解算的尺度系數也與驗證算法解算結果一致,但與本文LS算法計算的尺度系數不一致,說明LS準則下尺度系數的計算方法并不能用于TLS準則下尺度系數的計算。

(3)比較上述案例LS準則下和TLS準則下坐標轉換參數計算結果,可以發現兩種準則下旋轉角度計算結果完全一致,這與1.1節和1.2節推導結果一致。

3 相關性分析

3.1 LS及TLS準則下尺度系數分析

從LS準則下構建的優化函數式(6)及TLS準則下構建的優化函數式(11)可以看出,無論尺度系數取何值,函數取最小值時所求取的旋轉矩陣均不變,說明在LS準則和TLS準則下坐標轉換尺度系數和旋轉矩陣(旋轉角度)不相關;另根據式(4)可以看出,計算平移參數需要尺度系數及旋轉矩陣,故平移參數與尺度系數是相關的,這與文獻[1]中結論一致。

3.2 WTLS準則下尺度系數分析

“坐標轉換尺度系數與旋轉矩陣不相關,與平移參數相關”這一結論適用LS準則及TLS準則,但是否同樣適用于WTLS準則,目前的研究較少,故對于WTLS準則下該結論是否成立需要進一步的驗證。

為驗證WTLS準則下尺度系數與旋轉矩陣(旋轉角度)、平移參數的相關性,設計了以下六組權陣,每組權陣下再分別計算三種坐標轉換參數(第一種計算全部七參數;第二種和第三種分別指定尺度系數為1和2,計算坐標轉換六個參數),具體權陣數據見表3。

表3 WTLS準則下坐標轉換設計權陣

坐標轉換數據采用第3節案例分析中數據,并采用文獻[10]中算法進行坐標轉換參數解算,解算結果見表4。

表4 WTLS準則下坐標轉換參數解算結果

續表4

續表4

備注：所求得的角度均為弧度。

從表4計算結果可以看出,第一組權陣下,即對1.3節中矩陣C、D進行列加權、行不加權且列加權矩陣相同,當尺度系數變化時,旋轉角度不變,平移參數改變,說明在第一組權陣下尺度系數與旋轉角度不相關,與平移參數相關。第二組至第六組權陣下,尺度系數變化時,旋轉角度及平移參數均改變,說明在這些權陣下尺度系數與旋轉角度及平移參數均相關。

4 結束語

本文基于奇異值分解(SVD)算法推導了LS準則下和TLS準則下坐標轉換尺度系數計算方法,并結合案例采用BFGS優化算法進行了驗證,同時分析了LS、TLS及WTLS準則下坐標轉換尺度系數與旋轉角度、平移參數之間的相關性,并得出了以下結論:

(1)本文提出的LS準則下及TLS準則下尺度系數公式與最優化驗證算法計算結果完全一致,說明本文方法可用于尺度系數的精確計算。

(2)對比LS及TLS準則下的坐標轉換七參數計算結果,其旋轉矩陣(旋轉角度)相同,尺度系數不同,因尺度系數不同導致平移參數不同。

(3)在LS準則及TLS準則下,坐標轉換尺度系數和旋轉矩陣(旋轉角度)不相關,與平移參數相關;但在WTLS準則下,除個別特殊權陣外,尺度系數與旋轉矩陣和平移參數均相關。故在LS準則下及TLS準則下,坐標轉換尺度系數、旋轉矩陣及平移參數可分別計算,而在WTLS準則下坐標轉換七參數不能分別計算(個別特殊權陣除外),只能通過優化算法迭代一次性求取全部七個參數。