整體最小二乘和最小二乘擬合空間直線的比較

胡 川 陳 義，2 朱衛(wèi)東 錢承軍

1 同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海市四平路1239號，200092

2 現(xiàn)代工程測量國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室，上海市四平路1239號，200092

3 上海海洋大學(xué)海洋科學(xué)學(xué)院，上海市浦東新區(qū)南匯新城鎮(zhèn)滬城環(huán)路999號，200120

三維空間直線擬合是工業(yè)測量中一個(gè)常見問題［1］。空間直線不能直接應(yīng)用最小二乘法（least squares，LS）［2］，而是采用點(diǎn)到直線距離平方和最小的擬合準(zhǔn)則［3］進(jìn)行擬合。該準(zhǔn)則是一個(gè)非線性估計(jì)問題，文獻(xiàn)［4］將其線性化后采用帶約束條件的間接平差模型進(jìn)行迭代求解。這種基于高斯迭代解的算法對初值非常敏感，可能出現(xiàn)不收斂情況。文獻(xiàn)［3］用特征值分解法解決上述不收斂問題，并提出基于選權(quán)迭代的穩(wěn)健空間直線擬合算法。因點(diǎn)到直線的距離不是直接觀測值，該準(zhǔn)則的最大困難是無法直接采用觀測精度定權(quán)。

空間直線可用與坐標(biāo)平面垂直的兩個(gè)相交平面來表示。假設(shè)其中一個(gè)平面上的坐標(biāo)觀測值無誤差，用LS法估計(jì)兩個(gè)平面的參數(shù)，最后用兩個(gè)擬合的平面恢復(fù)空間直線［5］。如果僅擬合其中一個(gè)平面，然后擬合空間直線在該平面上的投影線，即是文獻(xiàn)［6］提出的無迭代算法。文獻(xiàn)［7］將空間直線垂直投影到平面上，進(jìn)行兩兩組合，選擇其中點(diǎn)到直線距離和最小的組合作為空間直線參數(shù)的擬合結(jié)果。上述3 種方法都假設(shè)某些坐標(biāo)無誤差，這與實(shí)際不符。文獻(xiàn)［2］提出采用整體最小二乘（total least squares，TLS）法擬合兩個(gè)垂直投影平面，考慮了所有坐標(biāo)誤差，更符合實(shí)際，但沒有考慮加權(quán)情況。

本文模擬各坐標(biāo)點(diǎn)等精度、非等精度和各坐標(biāo)分量非等精度的觀測數(shù)據(jù)。將空間直線投影到坐標(biāo)平面上，采用整體最小二乘法和最小二乘法對投影直線進(jìn)行擬合。比較3 種模擬場景下TLS和LS估計(jì)的參數(shù)和驗(yàn)后方差，并比較三維激光掃描數(shù)據(jù)的擬合結(jié)果。

1 TLS擬合空間直線

1.1 擬合公式推導(dǎo)

假設(shè)一條空間直線通過點(diǎn)P0（x0，y0，z0），方向向量為（F，G，H），則直線的對稱式方程為：

將上述空間直線投影到坐標(biāo)平面上，有：

以XOY投影面上直線為例。設(shè)ey和ex分別是y和x的誤差矢量，大小為m×1。它們的隨機(jī)屬性可以表達(dá)為［8，9］：

式中，L是由m個(gè)y分量組成的m×1矢量；A是m×2系數(shù)矩陣，其第一列由與y分量相對應(yīng)的x分量組成，第二列全為常數(shù)1；ξ是待估計(jì)參數(shù)矢量；轉(zhuǎn)換矩陣K計(jì)算公式為：

式（7）為矢量導(dǎo)數(shù)公式。引入TLS平差準(zhǔn)則［8］：

將式（9）分別對ey、ex、λ和ξ求偏導(dǎo)，并令其等于零，有［8］：

求解式（10）～（13），可以得到：

驗(yàn)后單位權(quán)方差為［8］：

1.2 空間直線重建

根據(jù)擬合的平面直線，可以找到一個(gè)過該直線且垂直該坐標(biāo)平面的平面，3個(gè)坐標(biāo)平面得到3個(gè)垂直平面，平面兩兩組合可以重建出3條空間直線［7］。平面y－a1x－b1＝0和z－a2y－b2＝0的交線為：

平面z－a2y－b2＝0和x－a3z－b3＝0的交線為：

平面x－a3z－b3＝0和y－a1x－b1＝0的交線為：

根據(jù)交叉出的3條空間直線的對稱式，可以得到直線上的一個(gè)已知點(diǎn)和方向。可以采用下式計(jì)算點(diǎn)到直線的距離：

Δxi、Δyi和Δzi是測量點(diǎn)與對稱方程中已知點(diǎn)的差值。

2 比較分析

2.1 模擬對比

設(shè)有已知空間直線：

在空間直線上任取10個(gè)精確坐標(biāo)點(diǎn)。設(shè)置如下3種場景：

1）設(shè)所有測量點(diǎn)等精度觀測，給各坐標(biāo)點(diǎn)附加期望為零、方差分別為0.000 001、0.000 1、0.001、0.01、1.0m2的隨機(jī)誤差；

2）設(shè)所有測量點(diǎn)非等精度觀測，按點(diǎn)的先后順序附加上期望為零、方差從0.1m2增加到1m2的隨機(jī)誤差；

3）設(shè)所有測量點(diǎn)的坐標(biāo)分量非等精度觀測，給x、y、z坐標(biāo)分量附加上期望為零、方差分別為0.8、0.4和0.2m2的隨機(jī)誤差。

在場景1中，協(xié)因數(shù)矩陣Qx、Qy和Qz都是單位陣，先驗(yàn)單位權(quán)方差分別為0.000 001、0.000 1、0.001、0.01、1.0 m2。在不同方差條件下，用LS和TLS各計(jì)算1 000次。參數(shù)估計(jì)值和驗(yàn)后單位權(quán)方差估計(jì)值的平均值列于表1。可以發(fā)現(xiàn)，在誤差比較小時(shí)，LS和TLS估計(jì)的參數(shù)值幾乎相同。當(dāng)誤差增大時(shí)，沒有證據(jù)表明TLS獲得的參數(shù)估計(jì)值比LS 更接近真實(shí)值。但是，多數(shù)情況下TLS的結(jié)果更接近真實(shí)值。不管在哪種方差條件下，TLS估計(jì)的驗(yàn)后單位權(quán)方差都比LS更接近先驗(yàn)值。

在場景2中，協(xié)因數(shù)矩陣Qx、Qy和Qz的對角線元素對應(yīng)于各點(diǎn)模擬方差值，非對角線元素全為零，先驗(yàn)單位權(quán)方差為1。用LS和TLS分別計(jì)算10 000次。表2描述了場景2中LS和TLS計(jì)算的參數(shù)和驗(yàn)后單位權(quán)方差估計(jì)值。可以看出，不管是參數(shù)估計(jì)還是單位權(quán)方差估計(jì)結(jié)果，TLS明顯要比LS更接近真實(shí)值，特別是方差估計(jì)值。

在場景3中，協(xié)因數(shù)矩陣Qx、Qy和Qz分別是0.8、0.4和0.2與單位矩陣的乘積，先驗(yàn)單位權(quán)方差為1m2，同樣用LS和TLS分別模擬計(jì)算1 000次，將計(jì)算的平均值和平均值與真值的差值列于表3。可以發(fā)現(xiàn)，TLS 獲得的參數(shù)估計(jì)值比LS更接近真實(shí)值。選擇不同的坐標(biāo)平面作為投影面，TLS 估計(jì)的參數(shù)與真值的差異非常小，但是LS對應(yīng)不同的投影面，其估計(jì)結(jié)果與真值的接近程度不相同。換言之，TLS估計(jì)參數(shù)的準(zhǔn)確性對投影面的依賴性比較小，而LS 依賴性較大。TLS獲得的驗(yàn)后單位權(quán)方差比LS更接近真值。

表1 場景1中LS和TLS計(jì)算的參數(shù)和驗(yàn)后單位權(quán)方差平均值Tab.1 Comparisons of the means of estimated parameters and variance of unit weight with LS and TLS in scenario 1

表2 場景2中LS和TLS計(jì)算的參數(shù)和驗(yàn)后單位權(quán)方差平均值以及與真值之差Tab.2 Comparisons of the mean of estimated parameters and variance of unit weight with LS and TLS，and the differences between the estimated value and real value in scenario 2

表3 場景3中LS和TLS計(jì)算的參數(shù)和驗(yàn)后單位權(quán)方差平均值以及與真值之差Tab.3 Comparisons of the mean of estimated parameters and variance of unit weight with LS and TLS，and the differences between the estimated value and real value in scenario 3

2.2 工程實(shí)踐

用Faro三維激光掃描儀對新建的南京青奧步行橋進(jìn)行三維激光掃描，提取六面體鋼結(jié)構(gòu)橋梁上兩條邊相交處的部分掃描點(diǎn)（圖1）。菱形點(diǎn)是空間點(diǎn)在XOY面上的投影，五角星是空間點(diǎn)在XOZ面上的投影，小圓圈是空間點(diǎn)在YOZ面上的投影。

圖1 掃描的三維坐標(biāo)點(diǎn)Fig.1 Laser scanned data points

假設(shè)各點(diǎn)是等精度觀測，權(quán)陣為單位矩陣，用LS法和TLS法分別擬合3條平面投影直線，將擬合的參數(shù)和驗(yàn)后單位權(quán)方差列于表4。可以看出，在XOY投影面的直線，兩種方法擬合結(jié)果和驗(yàn)后方差完全相同。在其他兩個(gè)投影面上，LS和TLS估計(jì)的參數(shù)和單位權(quán)方差出現(xiàn)較大差異，這與前面的模擬計(jì)算相符。但是，TLS估計(jì)的驗(yàn)后單位權(quán)方差也出現(xiàn)較大值，即明顯較大。

表4 LS和TLS擬合掃描點(diǎn)的參數(shù)估計(jì)值Tab.4 Estimated parameters and variance of unit weight with LS and TLS using the laser scanned data points

將擬合的平面直線參數(shù)按照式（20）、（21）和（22）重建空間直線，分別稱為空間直線（20）、（21）和（22）。采用式（23）分別計(jì)算各點(diǎn)到擬合空間直線的距離。圖2描述了各點(diǎn)到LS和TLS擬合的空間直線（20）、（21）和（22）的距離。可以發(fā)現(xiàn)，掃描點(diǎn)到TLS擬合的空間直線的距離都非常接近，相反到LS擬合的空間直線的距離存在較大的差異。這說明，不管選擇怎樣的投影面，TLS法都可以得到近似相同的一條空間直線，而LS法不行。

圖2 掃描點(diǎn)到擬合空間直線的距離Fig.2 The distance between the measured points and the fitted spatial lines

將各點(diǎn)到擬合空間直線距離之和的最大距離、最小距離和最大最小距離差列于表5。可以看出，TLS 重建的3 條直線的距離總和非常接近，而LS法重建的3 條直線距離總和存在較大差異，而且LS 的距離總和大于TLS 的距離總和。可以發(fā)現(xiàn)，TLS結(jié)果比LS 更加穩(wěn)定。這與前面的結(jié)論相同，LS擬合結(jié)果對投影面的選擇依賴性較強(qiáng)，而TLS的依賴性較弱。

表5 各點(diǎn)到空間直線的距離總和、最大值、最小值和最大小值之差／mTab.5 The sum of d，max d，min dand the differences between max and min，where dis the distance between the measured points and the fitted spatial lines／m

3 結(jié) 語

1）當(dāng)已知各點(diǎn)坐標(biāo)等精度觀測的情況下，LS和TLS估計(jì)的空間直線參數(shù)幾乎完全相同，但是TLS的驗(yàn)后方差分量估計(jì)值比LS的更接近先驗(yàn)值。此時(shí)，如果僅僅是為了估計(jì)參數(shù)，為了簡便可以直接采用LS擬合空間直線；如果需要對參數(shù)估計(jì)結(jié)果進(jìn)行統(tǒng)計(jì)評價(jià)，建議采用TLS擬合法。

2）當(dāng)各點(diǎn)非等精度觀測時(shí)，TLS估計(jì)的參數(shù)和驗(yàn)后單位權(quán)方差都比LS 的要更接近真實(shí)值。因此，此種情況建議使用TLS擬合法。

3）當(dāng)各坐標(biāo)分量非等精度觀測時(shí)，LS擬合結(jié)果對投影面的選擇具有較強(qiáng)的依賴性，而TLS的依賴性較弱。此時(shí)建議采用TLS進(jìn)行空間直線擬合，選擇任何兩個(gè)投影面皆可。

4）對于三維激光掃描數(shù)據(jù)，TLS的擬合結(jié)果比LS更優(yōu)。此時(shí)建議采用TLS法擬合，如果需要進(jìn)一步提高擬合精度，需要合理確定各坐標(biāo)分量的權(quán)值。

［1］陳基偉.工業(yè)測量數(shù)據(jù)擬合研究［D］.上海：同濟(jì)大學(xué)，2005（Chen Jiwei.Research on Industrial Measurement Data Fitting［D］.Shanghai：Tongji University，2005）

［2］姚宜斌，黃書華，孔建，等.空間直線擬合的整體最小二乘算法［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2014，39（5）：571－574（Yao Yibin，Huang Shuhua，Kong Jian，et al.Total Least Squares Algorithm for Fitting Spatial Straight Line［J］.Geomatics and Information Science of Wuhan University，2014，39（5）：571－574）

［3］潘國榮，唐杭.特征分解與選權(quán)迭代在空間直線擬合中的應(yīng)用［J］.東南大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，43（A2）：250－255（Pan Guorong，Tang Hang.Application of Eigen Decomposition and Selecting Weight Iteration in Spatial Line Fitting［J］.Journal of Southeast University：Natural Science Edition，2013，43（A2）：250－255）

［4］王解先，季凱敏.工業(yè)測量擬合［M］.北京：測繪出版社，2008（Wang Jiexian，Ji Kaimin.Industrial Measurement Fitting［M］.Beijing：Surveying and Mapping Press，2008）

［5］襲楊.空間直線擬合的一種方法［J］.齊齊哈爾大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，25（2）：64－68（Xi Yang.A Method for Fitting of a Space Straight Line［J］.Journal of Qiqihar University：Natural Science Edition，2009，25（2）：64－68）

［6］郭際明，向巍，尹洪斌.空間直線擬合的無迭代算法［J］.測繪通報(bào)，2011（2）：24－25（Guo Jiming，Xiang Wei，Yin Hongbin.Three－Dimensional Line Fitting without Iteration［J］.Bulletin of Surveying and Mapping，2011（2）：24－26）

［7］王偉鋒，溫耐.空間直線擬合研究［J］.許昌學(xué)院學(xué)報(bào)，2010（5）：37－39（Wang Weifeng，Wen Nai.Research on Fitting Method of Space Straight Line［J］.Journal of Xuchang University，2010，9（5）：37－39）

［8］Schaffrin B，Wieser A.On Weighted Total Least－Squares Adjustment for Linear Regression［J］.Journal of Geodesy，2008.82（7）：415－421

［9］胡川，陳義.解加權(quán)總體最小二乘平差問題的一種新方法［J］.大地測量與地球動力學(xué)，2012，32（6）：106－110（Hu Chuan，Chen Yi.An Innovation Algorithm for Solution of Weighted Total Least Squares Adjustment［J］.Journal of Geodesy and Geodynamics，2012，32（6）：106－110）