秦 鋒 張振虎
(1. 中國核工業(yè)集團三門核電有限公司 浙江 臺州 317112;2. 上??辈煸O(shè)計研究院(集團)有限公司 上海 200093)
空間直角坐標(biāo)轉(zhuǎn)換需要計算七參數(shù),包括尺度系數(shù)、旋轉(zhuǎn)角度(旋轉(zhuǎn)矩陣)及平移參數(shù)。目前坐標(biāo)轉(zhuǎn)換的算法主要分為兩類,一類是基于經(jīng)典最小二乘(least square,LS)的坐標(biāo)轉(zhuǎn)換算法,另一類是基于變量含誤差(errors in variables,EIV)模型的整體最小二乘(total least squares,TLS)或加權(quán)整體最小二乘(weighted total least squares,WTLS)坐標(biāo)轉(zhuǎn)換算法。
尺度系數(shù)通常認為是兩坐標(biāo)系距離之比。在經(jīng)典最小二乘準(zhǔn)則下,小角度坐標(biāo)轉(zhuǎn)換一般采用線性最小二乘一次性計算出包括尺度系數(shù)在內(nèi)的坐標(biāo)轉(zhuǎn)換七參數(shù)。對于基于經(jīng)典最小二乘準(zhǔn)則的大角度坐標(biāo)轉(zhuǎn)換,文獻[1]論證了尺度系數(shù)和旋轉(zhuǎn)矩陣(或旋轉(zhuǎn)角度)之間不相關(guān),和平移參數(shù)之間呈強相關(guān),故大部分基于經(jīng)典最小二乘準(zhǔn)則的大角度坐標(biāo)轉(zhuǎn)換算法在計算坐標(biāo)轉(zhuǎn)換參數(shù)時先計算尺度系數(shù),再計算旋轉(zhuǎn)矩陣和平移參數(shù)。如文獻[2]通過計算坐標(biāo)轉(zhuǎn)換前后所有距離比值的平均值來計算尺度系數(shù),文獻[3]對坐標(biāo)數(shù)據(jù)重心化處理后采用最大似然估計計算尺度系數(shù),文獻[4]對坐標(biāo)數(shù)據(jù)重心化處理后通過計算坐標(biāo)轉(zhuǎn)換前后重心化坐標(biāo)平方和比值的開方值來計算尺度系數(shù)。另有部分基于經(jīng)典最小二乘準(zhǔn)則的大角度坐標(biāo)轉(zhuǎn)換算法一次性計算出坐標(biāo)轉(zhuǎn)換七參數(shù),如文獻[5]對誤差方程線性化處理后通過迭代一次性求解坐標(biāo)轉(zhuǎn)換七參數(shù),不需要分別計算尺度系數(shù)、旋轉(zhuǎn)矩陣及平移參數(shù)。但在實際應(yīng)用中發(fā)現(xiàn),以上計算坐標(biāo)轉(zhuǎn)換尺度系數(shù)的方法計算結(jié)果互有差異,難以確定何種方法尺度系數(shù)計算結(jié)果最優(yōu)。
在TLS或WTLS準(zhǔn)則下,坐標(biāo)轉(zhuǎn)換參數(shù)中尺度系數(shù)的求解較LS準(zhǔn)則下復(fù)雜。目前已有的TLS或WTLS坐標(biāo)轉(zhuǎn)換方法主要分為兩類,一類是分別計算尺度系數(shù)、旋轉(zhuǎn)矩陣及平移參數(shù),另一類通過迭代一次性計算出全部參數(shù)。第一類方法中,如文獻[6]采用多元總體最小二乘方法分別求解坐標(biāo)轉(zhuǎn)換參數(shù),文獻[7]推導(dǎo)了特定加權(quán)矩陣下采用多元總體最小二乘求解坐標(biāo)轉(zhuǎn)換參數(shù)的方法。第二類方法中,如文獻[8]采用提出了一種基于TLS準(zhǔn)則的小角度坐標(biāo)轉(zhuǎn)換算法,文獻[9]采用改進的加權(quán)整體最小二乘算法用于小角度坐標(biāo)轉(zhuǎn)換,文獻[10]提出了一種通用的加權(quán)整體最小二乘坐標(biāo)轉(zhuǎn)換方法。
本文通過奇異值分解算法(singular value decomposition,SVD)推導(dǎo)了LS及TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換尺度系數(shù)的計算公式,并結(jié)合工程實例,與其他尺度系數(shù)計算方法進行了對比。同時,為驗證本文算法,采用了布羅伊登-弗萊徹-戈德法布-香農(nóng)(Broyden-Fletcher-Goldfarb-Shanno,BFGS)優(yōu)化算法進行結(jié)果驗證。對于WTLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換尺度系數(shù)的計算及其與旋轉(zhuǎn)角度、平移參數(shù)之間的相關(guān)性,本文基于文獻[10]中的方法進行了數(shù)據(jù)驗證,并得出了相應(yīng)結(jié)論。
坐標(biāo)轉(zhuǎn)換模型可表示為
為方便參數(shù)計算,根據(jù)文獻[11]對坐標(biāo)進行重心化處理。數(shù)據(jù)重心化處理后,旋轉(zhuǎn)矩陣和平移參數(shù)可分開進行計算,其計算公式為
分別定義矩陣A和B為
重心化后的坐標(biāo)轉(zhuǎn)換模型就可表示為:A=λRB
求取正交矩陣R,使A-λRB的弗羅貝紐斯(Frobenius)范數(shù)最小時所得的R即為最佳轉(zhuǎn)換矩陣,即求解以下問題,即
(5)
(6)
式中,tr是trace的簡寫,表示矩陣的跡。
因λ>0,故在tr(BATR)取最大值時,式(6)得到最小值。
對BAT進行奇異值分解,得
(7)
式中,U是左奇異矩陣;V是右奇異矩陣;Σ=diag(σ1,σ2,σ3),σ1,σ2,σ3為BAT的奇異值。根據(jù)文獻[12],R=VUT時,tr(BATR)取最大值。將求得的R代入式(6),此時式(6)為關(guān)于λ的一元二次方程求極值問題,當(dāng)λ為式(6)時,式(6)得到最小值。
(8)
平移參數(shù)可根據(jù)式(4)進行計算。
考慮到矩陣A、B中的隨機誤差EA、EB,構(gòu)建EIV模型為
(9)
假設(shè)隨機誤差EA、EB獨立同分布,基于拉格朗日乘子法,構(gòu)建整體最小二乘優(yōu)化函數(shù)為
(10)
函數(shù)f分別對EA、EB、μ求偏導(dǎo)并置0,解出EA、EB、μ后代入式(10),得
(11)
同式(7),對BAT進行奇異值分解;同理,當(dāng)且僅當(dāng)R=VUT時,f取極小值。將R=VUT代入式(11)中,可得
(12)
式(2)對λ求偏導(dǎo)并置0,因λ>0,tr(Σ)>0,故可解出
(13)
平移參數(shù)可根據(jù)式(4)進行計算。
分別定義以下矩陣:
其中,Ln表示元素全是1的列向量。
坐標(biāo)轉(zhuǎn)換模型又可表示為:
(14)
構(gòu)建基于LS準(zhǔn)則的坐標(biāo)轉(zhuǎn)換優(yōu)化函數(shù):
(15)
參照式(9)至式(11),構(gòu)建基于TLS準(zhǔn)則的坐標(biāo)轉(zhuǎn)換優(yōu)化函數(shù)為
(16)
對上述兩個函數(shù)分別采用BFGS算法進行優(yōu)化,迭代求取坐標(biāo)轉(zhuǎn)換七參數(shù)。具體BFGS算法可參照文獻[13]。
坐標(biāo)轉(zhuǎn)換數(shù)據(jù)選用文獻[10]中的數(shù)據(jù),詳細數(shù)據(jù)見表1。
表1 坐標(biāo)轉(zhuǎn)換原始數(shù)據(jù) 單位:m
分別采用文獻[2-5]中算法計算尺度系數(shù),同時采用本文算法及驗證算法分別計算LS準(zhǔn)則下及TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換七參數(shù),解算結(jié)果見表2。
表2 LS及TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換參數(shù)解算結(jié)果
從以上計算結(jié)果可以看出:
(1)本文LS坐標(biāo)轉(zhuǎn)換算法與基于BFGS優(yōu)化算法的LS坐標(biāo)轉(zhuǎn)換驗證算法解算的七參數(shù)完全一致,而文獻[2-5]尺度系數(shù)計算結(jié)果均與本文LS算法結(jié)果有差異,說明LS準(zhǔn)則下本文提出的尺度系數(shù)計算方法更準(zhǔn)確。
(2)在TLS準(zhǔn)則下,本文算法解算的尺度系數(shù)也與驗證算法解算結(jié)果一致,但與本文LS算法計算的尺度系數(shù)不一致,說明LS準(zhǔn)則下尺度系數(shù)的計算方法并不能用于TLS準(zhǔn)則下尺度系數(shù)的計算。
(3)比較上述案例LS準(zhǔn)則下和TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換參數(shù)計算結(jié)果,可以發(fā)現(xiàn)兩種準(zhǔn)則下旋轉(zhuǎn)角度計算結(jié)果完全一致,這與1.1節(jié)和1.2節(jié)推導(dǎo)結(jié)果一致。
從LS準(zhǔn)則下構(gòu)建的優(yōu)化函數(shù)式(6)及TLS準(zhǔn)則下構(gòu)建的優(yōu)化函數(shù)式(11)可以看出,無論尺度系數(shù)取何值,函數(shù)取最小值時所求取的旋轉(zhuǎn)矩陣均不變,說明在LS準(zhǔn)則和TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換尺度系數(shù)和旋轉(zhuǎn)矩陣(旋轉(zhuǎn)角度)不相關(guān);另根據(jù)式(4)可以看出,計算平移參數(shù)需要尺度系數(shù)及旋轉(zhuǎn)矩陣,故平移參數(shù)與尺度系數(shù)是相關(guān)的,這與文獻[1]中結(jié)論一致。
“坐標(biāo)轉(zhuǎn)換尺度系數(shù)與旋轉(zhuǎn)矩陣不相關(guān),與平移參數(shù)相關(guān)”這一結(jié)論適用LS準(zhǔn)則及TLS準(zhǔn)則,但是否同樣適用于WTLS準(zhǔn)則,目前的研究較少,故對于WTLS準(zhǔn)則下該結(jié)論是否成立需要進一步的驗證。
為驗證WTLS準(zhǔn)則下尺度系數(shù)與旋轉(zhuǎn)矩陣(旋轉(zhuǎn)角度)、平移參數(shù)的相關(guān)性,設(shè)計了以下六組權(quán)陣,每組權(quán)陣下再分別計算三種坐標(biāo)轉(zhuǎn)換參數(shù)(第一種計算全部七參數(shù);第二種和第三種分別指定尺度系數(shù)為1和2,計算坐標(biāo)轉(zhuǎn)換六個參數(shù)),具體權(quán)陣數(shù)據(jù)見表3。
表3 WTLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換設(shè)計權(quán)陣
坐標(biāo)轉(zhuǎn)換數(shù)據(jù)采用第3節(jié)案例分析中數(shù)據(jù),并采用文獻[10]中算法進行坐標(biāo)轉(zhuǎn)換參數(shù)解算,解算結(jié)果見表4。
表4 WTLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換參數(shù)解算結(jié)果
續(xù)表4
續(xù)表4
備注:所求得的角度均為弧度。
從表4計算結(jié)果可以看出,第一組權(quán)陣下,即對1.3節(jié)中矩陣C、D進行列加權(quán)、行不加權(quán)且列加權(quán)矩陣相同,當(dāng)尺度系數(shù)變化時,旋轉(zhuǎn)角度不變,平移參數(shù)改變,說明在第一組權(quán)陣下尺度系數(shù)與旋轉(zhuǎn)角度不相關(guān),與平移參數(shù)相關(guān)。第二組至第六組權(quán)陣下,尺度系數(shù)變化時,旋轉(zhuǎn)角度及平移參數(shù)均改變,說明在這些權(quán)陣下尺度系數(shù)與旋轉(zhuǎn)角度及平移參數(shù)均相關(guān)。
本文基于奇異值分解(SVD)算法推導(dǎo)了LS準(zhǔn)則下和TLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換尺度系數(shù)計算方法,并結(jié)合案例采用BFGS優(yōu)化算法進行了驗證,同時分析了LS、TLS及WTLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換尺度系數(shù)與旋轉(zhuǎn)角度、平移參數(shù)之間的相關(guān)性,并得出了以下結(jié)論:
(1)本文提出的LS準(zhǔn)則下及TLS準(zhǔn)則下尺度系數(shù)公式與最優(yōu)化驗證算法計算結(jié)果完全一致,說明本文方法可用于尺度系數(shù)的精確計算。
(2)對比LS及TLS準(zhǔn)則下的坐標(biāo)轉(zhuǎn)換七參數(shù)計算結(jié)果,其旋轉(zhuǎn)矩陣(旋轉(zhuǎn)角度)相同,尺度系數(shù)不同,因尺度系數(shù)不同導(dǎo)致平移參數(shù)不同。
(3)在LS準(zhǔn)則及TLS準(zhǔn)則下,坐標(biāo)轉(zhuǎn)換尺度系數(shù)和旋轉(zhuǎn)矩陣(旋轉(zhuǎn)角度)不相關(guān),與平移參數(shù)相關(guān);但在WTLS準(zhǔn)則下,除個別特殊權(quán)陣外,尺度系數(shù)與旋轉(zhuǎn)矩陣和平移參數(shù)均相關(guān)。故在LS準(zhǔn)則下及TLS準(zhǔn)則下,坐標(biāo)轉(zhuǎn)換尺度系數(shù)、旋轉(zhuǎn)矩陣及平移參數(shù)可分別計算,而在WTLS準(zhǔn)則下坐標(biāo)轉(zhuǎn)換七參數(shù)不能分別計算(個別特殊權(quán)陣除外),只能通過優(yōu)化算法迭代一次性求取全部七個參數(shù)。