加權(quán)總體最小二乘方法在 ITRF轉(zhuǎn)換中的應(yīng)用＊

陸 玨 陳 義 鄭 波

(1)同濟(jì)大學(xué)測量與國土信息工程系,上海 200092 2)現(xiàn)代工程測量國家測繪局重點(diǎn)實(shí)驗(yàn)室,上海 200092 3)上海市地籍事務(wù)中心上海市土地登記事務(wù)中心,上海 200003)

加權(quán)總體最小二乘方法在 ITRF轉(zhuǎn)換中的應(yīng)用＊

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

(1)同濟(jì)大學(xué)測量與國土信息工程系,上海 200092 2)現(xiàn)代工程測量國家測繪局重點(diǎn)實(shí)驗(yàn)室,上海 200092 3)上海市地籍事務(wù)中心上海市土地登記事務(wù)中心,上海 200003)

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

加權(quán)總體最小二乘(WTLS);誤差(EI V)模型;協(xié)因數(shù)陣;坐標(biāo)轉(zhuǎn)換;布爾莎模型

1 引言

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

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

在 ITRF轉(zhuǎn)換的觀測方程中,觀測向量 b和系數(shù)矩陣A分別由目標(biāo)坐標(biāo)系和原坐標(biāo)系下的點(diǎn)坐標(biāo)值構(gòu)成。這些點(diǎn)不僅都存在誤差,且由于觀測手段、觀測距離、網(wǎng)型結(jié)構(gòu)等原因使得它們的精度并不相等,且存在相關(guān)。

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

2 加權(quán)總體最小二乘方法

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

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

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

進(jìn)一步地,隨機(jī)誤差 e和 EA的數(shù)字特征可表示為：

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

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

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

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

第一步,解算初值：

第二步,迭代解算：

有偏的單位權(quán)方差為：

3 加權(quán)總體最小二乘方法在 ITRF轉(zhuǎn)換中的應(yīng)用

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

對(duì)于具有 k個(gè)已知點(diǎn)的三維坐標(biāo)轉(zhuǎn)換觀測方程而言,n=3k,m=4,Q0=I4×4。令原坐標(biāo)系下各點(diǎn)坐標(biāo)之間的協(xié)因數(shù)陣為：

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

由于

G∈Rn×4n,由協(xié)因數(shù)傳播定律[13]可知：

通過系數(shù)矩陣 A中逐個(gè)元素對(duì)原坐標(biāo)系下每個(gè)點(diǎn)的三維坐標(biāo)進(jìn)行求導(dǎo)得到：

γ∈R4n×n,可見 eA=γdh,根據(jù)協(xié)因數(shù)傳播定律[13]：

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

4 計(jì)算步驟與算例

4.1 計(jì)算步驟

應(yīng)用WTLS方法進(jìn)行 ITRF的轉(zhuǎn)換可以按以下步驟實(shí)現(xiàn)：

1)對(duì)兩套坐標(biāo)系下的點(diǎn)坐標(biāo)進(jìn)行重心化;

2)根據(jù)目標(biāo)坐標(biāo)系下點(diǎn)坐標(biāo)的方差協(xié)方差陣Qb計(jì)算觀測向量 b的權(quán)矩陣 Pb;

3)利用WLS方法根據(jù)式 (11)解算出坐標(biāo)轉(zhuǎn)換參數(shù)初值,從而計(jì)算式(15)中的矩陣 G;

4)根據(jù)原坐標(biāo)系下各點(diǎn)坐標(biāo)的協(xié)因數(shù)陣列出Qh,如式 (12)、(13)所示;

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

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

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

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

4.2 算例

為了驗(yàn)證上述算法的正確性,分別利用模擬數(shù)據(jù)和實(shí)際數(shù)據(jù),采用WLS、TLS和WTLS方法進(jìn)行計(jì)算,并從參數(shù)估值、單位權(quán)中誤差等幾個(gè)方面對(duì)結(jié)果進(jìn)行對(duì)比和分析。

例 1 先采用模擬三維小角度坐標(biāo)數(shù)據(jù)進(jìn)行解算,兩套坐標(biāo)系下對(duì)應(yīng)點(diǎn)的坐標(biāo)如表 1所示。

表1 原坐標(biāo)系和目標(biāo)坐標(biāo)系下對(duì)應(yīng)點(diǎn)的坐標(biāo)(單位:m)Tab.1 Coordi nates of corresponding poi nts in two coordi nate system s(un it:m)

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

第一次,對(duì)原坐標(biāo)系及目標(biāo)坐標(biāo)系下點(diǎn)的 X、Y、Z坐標(biāo)值分別加入單位權(quán)中誤差為 1 mm、均值為零,和單位權(quán)中誤差為 1 cm、均值為零的隨機(jī)誤差。WTLS方法計(jì)算出的觀測向量單位權(quán)中誤差和系數(shù)矩陣的單位權(quán)中誤差見圖 1和圖 2。

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

利用WLS、TLS和WTLS方法計(jì)算出的坐標(biāo)轉(zhuǎn)換參數(shù)平均值和模擬的坐標(biāo)轉(zhuǎn)換參數(shù)之差如表 2所示。

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

第二次,對(duì)原坐標(biāo)系及目標(biāo)坐標(biāo)系下點(diǎn)的 X、Y、Z坐標(biāo)值分別加入單位權(quán)中誤差為 5 mm、均值為零,和單位權(quán)中誤差為 2 cm、均值為零的隨機(jī)誤差。

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

此時(shí),WTLS方法計(jì)算出的觀測向量單位權(quán)中誤差和系數(shù)矩陣的單位權(quán)中誤差分別如圖 3和圖 4所示。

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

再次利用WLS、TLS和WTLS方法計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù)平均值和模擬值之差如表 3所示。

表3的結(jié)果再次證明了利用WTLS方法不僅可以準(zhǔn)確地計(jì)算出坐標(biāo)的轉(zhuǎn)換參數(shù),且結(jié)果更加接近模擬值。

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

表2 由WLS、TLS和W TLS方法計(jì)算出的坐標(biāo)轉(zhuǎn)換參數(shù)與模擬值之差Tab.2 D ifferences of transformation parameters between calculation results and si mulation values

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

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

表3 由WLS、TLS和W TLS方法計(jì)算出的坐標(biāo)轉(zhuǎn)換參數(shù)與模擬值之差Tab.3 D ifferences of transformation parameters betweencalculation results and si mulation values

而從圖 1至圖 4的結(jié)果可以看出,當(dāng)給兩套坐標(biāo)系下的點(diǎn)坐標(biāo)值分別加入不同量級(jí)的誤差時(shí), WTLS計(jì)算出的 e和 EA都能夠準(zhǔn)確地予以反應(yīng)。

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

表4 WLS、TLS以及W TLS計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù)結(jié)果Tab.4 Comparison among transformation parameters calculated byWLS,TLS andW TLS solutions

可以看出,相比WLS和 TLS方法,利用WTLS方法計(jì)算出的坐標(biāo)轉(zhuǎn)換參數(shù)與 IERS網(wǎng)站公布的結(jié)果更接近。而WLS和WTLS方法計(jì)算的單位權(quán)中誤差分別為 1.498 51和 1.342 82,這也再次證明了WTLS方法可以對(duì) EI V模型進(jìn)行合理的解算,從而計(jì)算出正確的、精度更高的坐標(biāo)轉(zhuǎn)換參數(shù)。

5 結(jié)論

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

2)當(dāng)觀測向量以及系數(shù)矩陣中的各元素不等精度、甚至相關(guān)時(shí),應(yīng)用WTLS方法解算 ITRF轉(zhuǎn)換參數(shù)更為合理;

3)通過模擬數(shù)據(jù)的計(jì)算可以看到,WTLS方法能夠?qū)υ鴺?biāo)系和目標(biāo)坐標(biāo)系下點(diǎn)的誤差分別改正,且改正量與加入的誤差相符合。

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 魏木生.廣義最小二乘問題的理論和計(jì)算[M].北京：科學(xué)出版社,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 陸玨,陳義,鄭波.總體最小二乘方法在三維坐標(biāo)轉(zhuǎn)換中的應(yīng)用[J].大地測量與地球動(dòng)力學(xué),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 陳義,陸玨,鄭波.總體最小二乘方法在空間后方交會(huì)中的應(yīng)用[J].武漢大學(xué)學(xué)報(bào) (信息科學(xué)版),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].上海：同濟(jì)大學(xué)出版社,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

國家自然科學(xué)基金(41074017)

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

P207

A