加權總體最小二乘方法在 ITRF轉換中的應用＊

陸 玨 陳 義 鄭 波

(1)同濟大學測量與國土信息工程系,上海 200092 2)現代工程測量國家測繪局重點實驗室,上海 200092 3)上海市地籍事務中心上海市土地登記事務中心,上海 200003)

加權總體最小二乘方法在 ITRF轉換中的應用＊

陸 玨1)陳 義1,2)鄭 波3)

(1)同濟大學測量與國土信息工程系,上海 200092 2)現代工程測量國家測繪局重點實驗室,上海 200092 3)上海市地籍事務中心上海市土地登記事務中心,上海 200003)

針對 ITRF轉換中,兩套坐標系下的點坐標值均存在誤差,且各點之間精度不等、甚至相關的情況,提出利用加權總體最小二乘方法對轉換參數進行解算。通過模擬數據和真實數據的解算證明了加權總體最小二乘方法在 ITRF轉換中的適用性,與其他方法相比,利用加權總體最小二乘方法能夠得到準確的、更為合理的轉換參數。

加權總體最小二乘(WTLS);誤差(EI V)模型;協因數陣;坐標轉換;布爾莎模型

1 引言

總體最小二乘方法在求解觀測方程 Ax=b時,除了按經典最小二乘方法求解時要考慮觀測向量 b的誤差外,還需考慮系數矩陣 A的誤差[1,2]。系數矩陣與各種模型改正、先驗信息的誤差有關,因此,總體最小二乘法相對于最小二乘法而言更合理和完整。

近年來,Schaffrin對約束總體最小二乘法的計算方法[3]、線性回歸加權總體最小二乘法的計算方法[4]、以及附有線性和二次約束條件的總體最小二乘法的計算方法[5]進行了探討。同時,在總體最小二乘方法的應用方面,Felus[6]和 Schaffrin[7,8]分別討論了基準轉換中系數矩陣存在誤差時,如何合理地確定基準轉換的參數,不過,他們討論的僅僅是 2維的情況,并且認為兩個坐標系統的坐標服從相同的誤差分布,且各個誤差相互獨立。陸玨[9]和陳義[10]討論了三維線性基準轉換、空間后方交會的總體最小二乘估計方法,但解算時也認為誤差服從相同的分布,且相互獨立。

在 ITRF轉換的觀測方程中,觀測向量 b和系數矩陣A分別由目標坐標系和原坐標系下的點坐標值構成。這些點不僅都存在誤差,且由于觀測手段、觀測距離、網型結構等原因使得它們的精度并不相等,且存在相關。

為了同時考慮觀測向量和系數矩陣的誤差,以及它們之間的不等精度和相關性,本文采用加權總體最小二乘 (WTLS)方法來解算 ITRF轉換問題。在對該方法原理和算法研究的基礎上,通過對模擬數據和 ITRF2005至 ITRF2000真實坐標轉換數據進行計算,并與加權最小二乘(WLS)方法和總體最小二乘(TLS)方法的計算結果進行比較,證明了利用WTLS方法計算 ITRF轉換參數能夠得到準確的、更合理的結果。

2 加權總體最小二乘方法

在變量誤差 (EI V)模型中,觀測向量 b和系數矩陣A均包含隨機誤差,因此觀測方程可定義為[]：

觀測向量 b∈Rn×1,它包含著 e。參數向量 x∈Rm×1。系數矩陣 A∈Rn×m,其中 n＞m=rank(A)。EA代表存在于 A陣中的隨機誤差矩陣,在最小二乘中 EA≡0[11]。

可見,與僅考慮觀測值向量 b中的誤差、忽略系數矩陣A中誤差的LS方法相比,TLS方法關心的則是當觀測向量 b和系數矩陣A中都含有誤差,需要同時平差的參數估計問題。

進一步地,隨機誤差 e和 EA的數字特征可表示為：

“vec”是指將矩陣從左至右按列拉直所得到的列向量,這里定義 vecEA=eA。為單位權方差,Qb和QA分別為 e和 eA的權逆陣,可表示為：

Q0∈Rm×m,QE∈Rn×n,?代表的是矩陣之間的直積,即“Kronecker-Zehfuss積”,具體表現為：M?N：=[mijN],M=[mij]。通過矩陣間的直積,可以得到系數矩陣A中誤差的權逆陣。此時,WTLS方法的目標函數為：

根據拉格朗日極值條件,建立方程：

其中λ為拉格朗日乘算子,EAx=(xT?In)eA。在拉格朗日條件方程的基礎上,對方程中的每個參數進行求導解算,最終可以得到WTLS方法的解算過程為[4]：

第一步,解算初值：

第二步,迭代解算：

有偏的單位權方差為：

3 加權總體最小二乘方法在 ITRF轉換中的應用

兩個不同時間下的 ITRF轉換是一個小角度的三維坐標轉換,若對原坐標系和目標坐標系下的點坐標都進行重心化,則可以將平移參數略去,只解算尺度參數和旋轉參數[12]。因此重心化后,新舊坐標系重心化后的布爾莎模型可表示為：其中 (x2i,y2i,z2i)為第 i個控制點在目標坐標系下的坐標,它們組成了觀測向量 b。(x1i,y1i,z1i)為對應點在原坐標系下的坐標,它們構成了系數矩陣A中的變量。這里共有 k個對應點。μ、(εx,εy,εz)分別為坐標轉換參數中的 1個尺度參數以及 3個分別對x軸、y軸和 z軸的旋轉角參數。

對于具有 k個已知點的三維坐標轉換觀測方程而言,n=3k,m=4,Q0=I4×4。令原坐標系下各點坐標之間的協因數陣為：

若已知原坐標系和目標坐標系下各點之間的方差協方差陣,則Qorg(Qorg∈R3k×3k)和Qb已知。下面介紹如何解算QE：

由于

G∈Rn×4n,由協因數傳播定律[13]可知：

通過系數矩陣 A中逐個元素對原坐標系下每個點的三維坐標進行求導得到：

γ∈R4n×n,可見 eA=γdh,根據協因數傳播定律[13]：

因此,將式 (12)和 (16)代入式 (18)可得 QeA,再將QeA代入式 (15)即可計算出矩陣 QE。而對于式(15)的矩陣 G中的坐標轉換參數可先由WLS方法求出,之后再由WTLS方法不斷迭代更新代入解算。

4 計算步驟與算例

4.1 計算步驟

應用WTLS方法進行 ITRF的轉換可以按以下步驟實現：

1)對兩套坐標系下的點坐標進行重心化;

2)根據目標坐標系下點坐標的方差協方差陣Qb計算觀測向量 b的權矩陣 Pb;

3)利用WLS方法根據式 (11)解算出坐標轉換參數初值,從而計算式(15)中的矩陣 G;

4)根據原坐標系下各點坐標的協因數陣列出Qh,如式 (12)、(13)所示;

5)由式 (16)計算γ,并將γ和 Qh代入式 (18)計算QeA;

6)將步驟 3)和 5)計算出的 G和 QeA代入式(15)得到 QE;

7)由式(6)和(7)計算 ITRF轉換參數的WTLS解,再次帶入式(15)更新 G和QE的值;

8)最后,由式 (8)和式 (9)分別計算觀測向量和系數矩陣的改正數,并由式(10)計算WTLS的單位權方差。

4.2 算例

為了驗證上述算法的正確性,分別利用模擬數據和實際數據,采用WLS、TLS和WTLS方法進行計算,并從參數估值、單位權中誤差等幾個方面對結果進行對比和分析。

例 1 先采用模擬三維小角度坐標數據進行解算,兩套坐標系下對應點的坐標如表 1所示。

表1 原坐標系和目標坐標系下對應點的坐標(單位:m)Tab.1 Coordi nates of corresponding poi nts in two coordi nate system s(un it:m)

表1中,模擬的坐標轉換參數為 (x0,y0,z0)= (10,10,10),μ=1.01,(εx,εy,εz)=(0.05°,0.02°, 0.08°)。為了驗證WTLS方法能夠準確地反應出兩套坐標系下點的誤差情況,并進行改正,采用了1 000組模擬數據。

第一次,對原坐標系及目標坐標系下點的 X、Y、Z坐標值分別加入單位權中誤差為 1 mm、均值為零,和單位權中誤差為 1 cm、均值為零的隨機誤差。WTLS方法計算出的觀測向量單位權中誤差和系數矩陣的單位權中誤差見圖 1和圖 2。

圖1和圖 2所表示的觀測向量單位權中誤差和系數矩陣的單位權中誤差均值分別為 0.988 2 cm和 0.983 mm,與加入的誤差相匹配。

利用WLS、TLS和WTLS方法計算出的坐標轉換參數平均值和模擬的坐標轉換參數之差如表 2所示。

由表 2可看出,相比于WLS和 TLS方法,利用WTLS計算的結果更加接近模擬值。

第二次,對原坐標系及目標坐標系下點的 X、Y、Z坐標值分別加入單位權中誤差為 5 mm、均值為零,和單位權中誤差為 2 cm、均值為零的隨機誤差。

圖1 WTLS方法估計的觀測向量的單位權中誤差(方案1)Fig.1 Unit weight mean square error of observation vector estimated byWTLS solution(caseⅠ)

此時,WTLS方法計算出的觀測向量單位權中誤差和系數矩陣的單位權中誤差分別如圖 3和圖 4所示。

圖3和圖 4計算的觀測向量單位權中誤差和系數矩陣的單位權中誤差均值分別為 1.960 1 cm和4.907 nn,同樣與加入的誤差相匹配。

再次利用WLS、TLS和WTLS方法計算坐標轉換參數平均值和模擬值之差如表 3所示。

表3的結果再次證明了利用WTLS方法不僅可以準確地計算出坐標的轉換參數,且結果更加接近模擬值。

圖2 WTLS方法估計的系數矩陣的單位權中誤差(方案1)Fig.2 Unit weight mean square error of coefficient matrix estimated byWTLS solution(caseⅠ)

表2 由WLS、TLS和W TLS方法計算出的坐標轉換參數與模擬值之差Tab.2 D ifferences of transformation parameters between calculation results and si mulation values

圖3 WTLS方法估計的觀測向量的單位權中誤差(方案2)Fig.3 Unit weight mean square error of observation vector estimated byWTLS solution(caseⅡ)

圖4 WTLS方法估計的系數矩陣的單位權中誤差(方案2)Fig.4 Unit weight mean square error of coefficient matrixestimated byWTLS solution(caseⅡ)

表3 由WLS、TLS和W TLS方法計算出的坐標轉換參數與模擬值之差Tab.3 D ifferences of transformation parameters betweencalculation results and si mulation values

而從圖 1至圖 4的結果可以看出,當給兩套坐標系下的點坐標值分別加入不同量級的誤差時, WTLS計算出的 e和 EA都能夠準確地予以反應。

例 2 采用國際地球自轉服務(IERS)網站上公布的全球較穩定的 70個站點在 ITRF2005以及ITRF2000坐標框架下的坐標、速度以及它們之間的方差協方差陣,分別采用WLS、TLS和WTLS方法進行計算。3種方法的計算結果以及與網站公布的坐標轉換參數見表 4,其中,D=μ-1。

表4 WLS、TLS以及W TLS計算坐標轉換參數結果Tab.4 Comparison among transformation parameters calculated byWLS,TLS andW TLS solutions

可以看出,相比WLS和 TLS方法,利用WTLS方法計算出的坐標轉換參數與 IERS網站公布的結果更接近。而WLS和WTLS方法計算的單位權中誤差分別為 1.498 51和 1.342 82,這也再次證明了WTLS方法可以對 EI V模型進行合理的解算,從而計算出正確的、精度更高的坐標轉換參數。

5 結論

1)總體最小二乘方法不僅顧及觀測向量中的誤差,而且顧及包含變量的系數矩陣中的誤差,因此理論上更加嚴密,估計效果更佳,可用于 ITRF轉換的計算中;

2)當觀測向量以及系數矩陣中的各元素不等精度、甚至相關時,應用WTLS方法解算 ITRF轉換參數更為合理;

3)通過模擬數據的計算可以看到,WTLS方法能夠對原坐標系和目標坐標系下點的誤差分別改正,且改正量與加入的誤差相符合。

1 Golub H G and Van Loan F C.An analysis of the total least squares problem[J].SI AM Journal on Numerical Analysis, 1980,17(6)：883-893.

2 魏木生.廣義最小二乘問題的理論和計算[M].北京：科學出版社,2006.(WeiMusheng.Generalized least squares theory and calculation[M].Beijing：Science Press,2006)

3 Burkhard Schaffrin.A note on constrained total leastsquares estimation[J].LinearAlgebra and itsApplications, 2006,417：245-258.

4 Schaffrin B andW ieserA.On weighted total least-squares adjust ment for linear regression[J].Journal of Geodesy, 2008,82,415-421.

5 Schaffrin B and Felus Y A.An algorithmic approach to the tatal least-squares problem with linear and quadratic constraints[J].Stud Geophys Geod.,2009,53：1-16.

6 Yaron A and FelusBurkhard Schaffrin.Perfor ming similarity transformations using the error-in-variables model[R]. ASPRS 2005 Annual Conference Balti more,Maryland, March 7-11,2005.

7 Schaffrin B and Felus YA.On the multivariate total leastsquares approach to empirical coordinatetransformation：Three algorithms[J].J Geodesy,2008,82：373-383.

8 SchaffrinB,et al.Total least-squares for geodetic straightline and plane adjustment[J].Boll Geod Sci Aff.,2006, 65：141-168.

9 陸玨,陳義,鄭波.總體最小二乘方法在三維坐標轉換中的應用[J].大地測量與地球動力學,2008,(5)：77-81.(Lu Jue,Chen Yi and ZhengBo.Research study on three -dimensional datum transformation using total least squares [J].Journal of Geodesy and Geodynamics,2008,(5)：77 -81)

10 陳義,陸玨,鄭波.總體最小二乘方法在空間后方交會中的應用[J].武漢大學學報 (信息科學版),2008,33 (12)：1 271-1 274. (Chen Yi,Lu Jue and Zheng Bo. Application of total least squares to space resection[J]. Geomatics and Information Science of Wuhan University 2008,33(12)：1 271-1 274)

11 Van Huffel S and Vandeualle J.The total least squares problem[J].Computational Aspects and Analysis,Frontiers in Applied,Mathematics,1991,(9)：1-87.

12 Shen Y Z,Chen Y and Zheng D H.A quaternion-based geodetic datum transformation algorithm[J]. Springer-Verlag J Geod,2006,80：233-239.

13 樊功瑜.誤差理論與測量平差[M].上海：同濟大學出版社,1998.(Fan Gongyu.Error theory and surveying adjustment[M].Shanghai：TongjiUniversity Press,1998)

APPL ICATI ON OFW EIGHTED TOTAL LEAST SQUARES IN ITRF TRANSFORMATI ON

Lu Jue1),Chen Yi1,2)and ZhengBo3)

(1)The Departm ent of Surveying and Geo-infor m atics,Tongji University,Shanghai 200092 2)Key Laboratory of Advanced Surveying Engineering of SBSM,Shanghai 200092 3)Shanghai CadastreM aragenent Center,Shanghai 200003)

According to the fact that the points in two transformational coordinate systems are all affected by random errors,which make the observation vector and coefficientmatrix in error equations both include errors,and even,these observations and elements in coefficient matrix may be heteroscedastic and correlated.We employ the weighted total least squares solution to calculate the parametersof ITRF transformation.To demonstrate the performance and efficiency of the application ofWTLS solution in coordinate transformation,the simulation and real experiments are i mplemented.The results show that theWTLS solution is correct and more reasonable.

weighted total least squares;errormodel;cofactormatrix;coordinate transformation;Bursa model

1671-5942(2011)04-0084-06

2011-01-24

國家自然科學基金(41074017)

陸玨,女,1985年生,博士生,研究方向：測量數據處理、攝影測量.E-mail：lujue1985@126.com

P207

A