基于選權(quán)迭代的總體最小二乘算法在三維坐標轉(zhuǎn)換中的應用

徐 波 高井祥 李增科 劉 洋

1 中國礦業(yè)大學國土環(huán)境與災害監(jiān)測國家測繪地理信息局重點實驗室，徐州市大學路1號，221116

在進行三維坐標轉(zhuǎn)換時，利用LS 算法求解轉(zhuǎn)換參數(shù)只考慮觀測向量中的誤差，而沒有考慮系數(shù)矩陣中的誤差。總體最小二乘（TLS）算法將觀測向量和系數(shù)矩陣中的誤差都考慮在一起。文獻［1－2］介紹了總體最小二乘算法在坐標轉(zhuǎn)換中的應用，并與最小二乘算法進行比較，認為TLS算法較LS算法更加合理。但是當控制點中存在粗差時，LS 和TLS 都不能得到有效的參數(shù)值。將抗差理論與LS和TLS相結(jié)合，可有效定位存在粗差的控制點，并對粗差控制點進行定權(quán)，求得未知參數(shù)值［3］。本文根據(jù)實例得出，基于抗差的LS算法和基于抗差的TLS算法能夠得到有效的未知參數(shù)，但是抗差TLS 較抗差LS 收斂更快，并且隨著存在粗差控制點數(shù)目的增加，TLS算法優(yōu)勢更明顯。

1 坐標轉(zhuǎn)換模型

空間三維坐標轉(zhuǎn)換通常采用七參數(shù)模型：

式中，［XT YT ZT］T和［XS YS ZS］T分別是控制點在目標坐標系和原坐標系的空間三維坐標，［ΔX ΔY ΔZ］T為平移參數(shù)，［εX εY εZ］T為旋轉(zhuǎn)參數(shù)，μ為尺度參數(shù)［1，4］。

當控制點數(shù)目大于2時上式的線性模型為：

2 模型參數(shù)解算方法

2.1 最小二乘與總體最小二乘算法

當觀測數(shù)大于必要觀測數(shù)時，可以根據(jù)最小二乘算法求解未知參數(shù)的最優(yōu)解［2］。但最小二乘算法只考慮了觀測量的誤差，沒有考慮系數(shù)矩陣的誤差，在解決高精度坐標轉(zhuǎn)換時嚴密性較差。

總體最小二乘算法在最小二乘算法基礎(chǔ)上，將系數(shù)矩陣的誤差加入到平差過程中：

vec（ΔB）為矩陣的列向量，ΔB、ΔL為隨機誤差矩陣。式（3）可以化為：

對式（5）采用奇異值分解SVD 來求解參數(shù)的TLS解。求出未知參數(shù)X的估值其 中為矩陣［B，L］T·個特征向量組成的正交矩陣。單位權(quán)方差及參數(shù)的協(xié)方差矩陣為：

利用總體最小二乘算法求解空間三維坐標轉(zhuǎn)換參數(shù)時，可以有效地控制觀測量和系數(shù)矩陣的誤差，求得未知參數(shù)的最優(yōu)解。但是，當控制點坐標含有粗差時，僅利用總體最小二乘算法并不能有效地定位并剔除含有粗差的控制點，需要將抗差和總體最小二乘算法聯(lián)合起來［5－9］。

2.2 基于選權(quán)迭代的總體最小二乘算法

基于選權(quán)迭代的總體最小二乘算法的關(guān)鍵在于定權(quán)，定權(quán)的方法有Hample權(quán)函數(shù)、周江文（IGG）權(quán)函數(shù)、Hubert權(quán)函數(shù)等。試驗證明，在參數(shù)選取合理且粗差大于最小可探測偏差時，利用上述幾種選權(quán)函數(shù)并根據(jù)本文的抗差TLS方法探測和定位粗差均能達到要求［10］。

借助IGG I權(quán)函數(shù)，并取k0＝max（abs（v（i））／σ，為防止由于k0取值原因而產(chǎn)生秩虧現(xiàn)象，若k0＝max（abs（v（i））／σ就不會剔除任何觀測值，取權(quán)函數(shù)為，σ為方差因子，式中取k0＝3.0。

將IGG I權(quán)函數(shù)和總體最小二乘算法結(jié)合起來，IGG I權(quán)函數(shù)需要分別對觀測值向量和系數(shù)矩陣求權(quán)。設(shè)系數(shù)矩陣B的權(quán)為，觀測向量L的權(quán)為，表示為：

設(shè)原坐標系和目標坐標系下觀測值的協(xié)因數(shù)陣分別為Q1和Q2，對應的權(quán)陣為P1和P2。假設(shè)控制點i在原坐標系和目標坐標系下的坐標分別為和，對應的隨機誤差分布為和，則坐標變化關(guān)系式可以表示如下：

因此可以將誤差向量v寫成如下形式：

則選權(quán)迭代的總體最小二乘方程式可以寫成：

坐標轉(zhuǎn)換方程式可以寫成：

設(shè)誤差向量和未知參數(shù)的初值分別為v0和ξ0，對上式線性化得：

可以求得未知參數(shù)和拉格朗日乘數(shù)的改正數(shù)為：

通過迭代計算直到參數(shù)向量滿足‖ξv－ξv－1‖≤η即可，通常取η＝10－10。

3 算例分析

為驗證基于選權(quán)迭代的總體最小二乘算法在求解三維坐標轉(zhuǎn)換參數(shù)中的優(yōu)越性，隨機生成10個點。設(shè)七參數(shù)為（ΔX，ΔY，ΔZ）＝（100，50，20），（ωx，ωy，ωz）＝（10，15，20），旋轉(zhuǎn)參數(shù)單位為（″），μ＝1.001，在原坐標系統(tǒng)中的位置和轉(zhuǎn)換后在目標坐標系中的位置可以在圖1、圖2看出。

圖1 控制點在原坐標系中的分布Fig.1 The distribution of the control points in the original coordinate system

給2號點在原坐標和目標坐標中分別加入粗差（0.5，0.5，0.5）和（－0.5，－0.5，－0.5），分別用上述方法算出兩套坐標的轉(zhuǎn)換參數(shù)（表1）。

由表1看出，當控制點存在粗差時，利用LS和TLS算法求解三維坐標轉(zhuǎn)換參數(shù)會出現(xiàn)數(shù)值失真的情況，而抗差LS和抗差TLS算法求解的未知參數(shù)與真值相近。在收斂速度方面，抗差TLS算法的收斂速度較抗差LS算法快［10］。

圖2 控制點在目標坐標系中的分布Fig.2 The distribution of the control points in the target coordinate system

表1 1個點存在粗差時4種方法計算出的轉(zhuǎn)換參數(shù)Tab.1 Four methods to calculate the conversion parameters when one－point with gross errors

假設(shè)有2個控制點存在粗差，將2號點和8號點分別在原坐標和目標坐標中加入粗差（0.5，0.5，0.5）和（－0.5，－0.5，－0.5），并進行計算（表2）。

表2 2個點存在粗差時4種方法計算出的轉(zhuǎn)換參數(shù)Tab.2 Four methods to calculate the conversion parameters when two－point with gross errors

表2同樣反映出，當控制點存在粗差時，LS和TLS算法求解三維坐標轉(zhuǎn)換參數(shù)會出現(xiàn)失真情況，抗差LS和抗差TLS算法求解的轉(zhuǎn)換參數(shù)與真值差值較小。

假設(shè)有5個控制點存在粗差，將2、5、6、8、9號點在原坐標和目標坐標中分別加入粗差（0.5，0.5，0.5）和（－0.5，－0.5，－0.5），并進行計算（表3）。

表3 5個點存在粗差時4種方法計算出的轉(zhuǎn)換參數(shù)Tab.3 Four methods to calculate the conversion parameters when five－point with gross errors

表3同樣反映了表1、2中LS和TLS存在的問題，但是隨著存在粗差的控制點數(shù)目的增加，抗差LS算法求解的三維坐標轉(zhuǎn)換參數(shù)與真值之差也變大，而抗差TLS算法求解的未知參數(shù)并沒有明顯變大，表明隨著存在粗差的控制點數(shù)目的增加，抗差TLS算法較其他3種方法更加有效。

4 結(jié) 語

1）當控制點存在粗差時，利用LS和TLS求解轉(zhuǎn)換參數(shù)時會發(fā)生嚴重的失真，而利用抗差LS和抗差TLS求解參數(shù)能夠有效消除粗差對轉(zhuǎn)換參數(shù)的影響。通過模擬數(shù)據(jù)分析，基于抗差TLS求解轉(zhuǎn)換參數(shù)比基于抗差LS求解的效果更好一些。

2）隨著存在粗差控制點數(shù)目的增多，LS、抗差LS和TLS所計算出的轉(zhuǎn)換參數(shù)與真值的偏差越來越大，根據(jù)IGG I的定權(quán)結(jié)合TLS算法計算出的轉(zhuǎn)換參數(shù)無明顯變化。

3）對于空間轉(zhuǎn)換模型而言，需要控制點在參心坐標系下的三維空間直角坐標求解模型參數(shù)。由于參心坐標系是通過光學觀測建立的，控制點缺乏精確的大地高信息，使得控制點在參心坐標系下的空間三維直角坐標無法精確得到。這方面的問題還需繼續(xù)研究。

［1］施一民.現(xiàn)代大地控制測量［M］.北京：測繪出版社，2003（Shi Yimin.Modern Geodetic Control Survey［M］.Beijing：Surveying and Mapping Press，2003）

［2］Golub H G，Vanl F C.An Analysis of the Total Least Squares Problem［J］.SIAM Journal on Numerical Analysis，1980，17（6）：883－893

［3］周江文.經(jīng)典誤差理論與抗差估計［J］.測繪學報，1989，18（2）：116－120（Zhou Jiangwen.Classical Theory of Errors and Robust Estimation［J］.Acta Geodaetica et Cartographica Sinica，1989，18（2）：116－120）

［4］陶本藻，邱衛(wèi)寧.誤差理論與測量平差基礎(chǔ)［M］.武漢：武漢大學出版社，2010（Tao Benzao，Qiu Weining.Error Theory and Fundation of Surveying Adjustment［M］.Wuhan，Wuhan University Press，2010）

［5］黃令勇，呂志平，任雅奇，等.多元總體最小二乘在三維坐標轉(zhuǎn)換中的應用［J］.武漢大學學報：信息科學版，2014（7）：793－798（Huang Lingyong，LüZhiping，Ren Yaqi，et al.Application of Multivariate Total Least Square In Three－Dimensional Coordinate Transformation［J］.Geomatics and Information Science of Wuhan Unversity，2014（7）：793－798）

［6］Tao Y Q，Gao J X，Yao Y F.TLS Algorithm for GPS Height Fitting Based on Robust Estimation［J］.Survey Review，2014，46（336）：184－188

［7］王樂洋，許才軍.總體最小二乘研究進展［J］.武漢大學學報：信息科學版，2013，38（7）：850－856（Wang Leyang，Xu Caijun.Progress in Total Least Squares［J］.Geomatics and Information Science of Wuhan Unversity，2013，38（7）：850－856）

［8］楊仕平，范東明，龍玉春.基于整體最小二乘法的任意旋轉(zhuǎn)角度三維坐標轉(zhuǎn)換［J］.大地測量與地球動力學，2013，33（2）：114－119（Yang Shiping，F(xiàn)an Dongming，Long Yuchun.Three－Dimensional Coordination Transformation Adapted to Arbitrary Rotation Angle Based on Total Least Squares Method［J］.Journal of Geodesy and Geodynamics，2013，33（2）：114－119）

［9］陸玨，陳義，鄭波.總體最小二乘算法在三維坐標轉(zhuǎn)換中的應用［J］.大地測 量與地球動 力學，2008，28（5）：77－81（Lu Jue，Chen Yi，Zheng Bo.Application Total Least Squares To Three Dimensional Datum Transformation［J］.Geomatics and Information Science of Wuhan Unversity，2008，28（5）：77－81）

［10］龔循強，李志林.一種利用IGG II方案的穩(wěn)健混合總體最小二乘方法［J］.武漢大學學報：信息科學版，2014，39（4）：462－466（Gong Xunqiang，Li Zhilin.A Robust Mixed LSTLS Based on IGG II Scheme［J］.Geomatics and Information Science of Wuhan Unversity，2014，39（4）：462－466）