趙 輝 張書(shū)畢 張秋昭
(1)中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院,徐州 221008 2)國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室,徐州221008)
基于加權(quán)總體最小二乘法的GPS高程擬合*
趙 輝1,2)張書(shū)畢1)張秋昭1,2)
(1)中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院,徐州 221008 2)國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)國(guó)家測(cè)繪局重點(diǎn)實(shí)驗(yàn)室,徐州221008)
在GPS高程擬合中,針對(duì)傳統(tǒng)最小二乘方法不能解決系數(shù)矩陣存在誤差的問(wèn)題,提出了一種基于加權(quán)總體最小二乘的擬合方法。對(duì)平面和二次曲面多項(xiàng)式建立更加合理的擬合模型,并給出了相應(yīng)的迭代算法。實(shí)例計(jì)算表明,加權(quán)最小二乘方法能夠得到更好的估計(jì)參數(shù),高程異常值擬合精度也相應(yīng)提高。
GPS高程;多項(xiàng)式擬合;EIV模型;最小二乘;加權(quán)總體最小二乘
在實(shí)際工程應(yīng)用中,我國(guó)高程系統(tǒng)普遍采用正常高。傳統(tǒng)水準(zhǔn)測(cè)量獲取基于似大地水準(zhǔn)面的正常高,而GPS定位技術(shù),也可以得到基于WGS84參考橢球的大地高,且操作簡(jiǎn)單方便[1]。用GPS定位技術(shù)獲取的正常高去滿足工程應(yīng)用需要,需進(jìn)行高程系統(tǒng)的轉(zhuǎn)換,也就必須要知道測(cè)點(diǎn)的高程異常值。
對(duì)于大多數(shù)工程應(yīng)用,由于測(cè)區(qū)范圍不大,高程異常變化平緩,多項(xiàng)式擬合便可滿足精度要求[2]。該方法的基本原理是將高程異常值與平面坐標(biāo)近似描述為多項(xiàng)式關(guān)系,利用同時(shí)已知大地高和正常高的公共點(diǎn)組成誤差方程,根據(jù)最小二乘原理求解多項(xiàng)式系數(shù)。然而采用經(jīng)典的Gauss-Markov模型,并不能解決在系數(shù)矩陣中含有誤差的問(wèn)題。如果觀測(cè)向量和系數(shù)陣都存在誤差,那么最小二乘解將不再最優(yōu),而是有偏的[3]。為此,本文引入基于加權(quán)總體最小二乘(WTLS)的GPS高程多項(xiàng)式擬合,并利用實(shí)例數(shù)據(jù)對(duì)算法可行性進(jìn)行驗(yàn)證。
觀測(cè)變量含有誤差的EIV(Errors-In-Variable)線性函數(shù)關(guān)系式為[4,5]
其中,y為含有隨機(jī)誤差ey的n維觀測(cè)向量,A為含有隨機(jī)誤差EA的n×m維系數(shù)矩陣,x為m維待估參數(shù)向量。隨機(jī)誤差具有如下統(tǒng)計(jì)性質(zhì):
式中,vec表示矩陣列向量化算子,協(xié)因數(shù)陣Qy=,QA=。在總體最小二乘準(zhǔn)則下,可解得顧及系數(shù)陣誤差的參數(shù)估計(jì)
對(duì)于多元同方差的EIV模型,觀測(cè)值權(quán)陣Py和系數(shù)矩陣權(quán)陣PA均為單位陣,可采用迭代法解算參數(shù)估值[6,7]:
對(duì)于廣義總體最小二乘(GTLS)中的系數(shù)陣權(quán)陣,按
通過(guò)迭代可求解參數(shù)的WTLS估計(jì),并進(jìn)行精度評(píng)定。迭代步驟如下[8,9]:
GPS高程擬合多用平面和二次曲面多項(xiàng)式,考慮到由于系數(shù)陣存在誤差,本文引入基于加權(quán)總體最小二乘的解算方法。
二次曲面多項(xiàng)式擬合函數(shù)模型為:
式中,ζi為平面坐標(biāo)(xi,yi)點(diǎn)的高程異常,ai(i=0,…,5)為多項(xiàng)式待估參數(shù)。可寫(xiě)成形如式(1)的矩陣形式,式中
設(shè)平面坐標(biāo)和高程異常值具有同等精度,σx= σy=σζ。對(duì)于平面多項(xiàng)式,只需取系數(shù)陣A的前3列。由于第一列為常數(shù),從而可令P0第一列為零,則
采用某市20個(gè)GPS聯(lián)測(cè)水準(zhǔn)的高程數(shù)據(jù),分析總體最小二乘對(duì)多項(xiàng)式GPS高程擬合的可行性與有效性,取均勻分布的前10個(gè)點(diǎn)為擬合公共點(diǎn),其余10個(gè)點(diǎn)為外部檢核點(diǎn)。點(diǎn)位分布如圖1所示。
為保證計(jì)算精度,先對(duì)坐標(biāo)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理,組成觀測(cè)方程。分別采用傳統(tǒng)最小二乘、總體最小二乘和加權(quán)最小二乘求得平面和曲面多項(xiàng)式的擬合參數(shù),表1和表2中列出了3種方法的參數(shù)估值與單位權(quán)中誤差。
由表1和表2可以看出,3種方法對(duì)平面擬合和二次曲面擬合參數(shù)求解,采用傳統(tǒng)最小二乘和總體最小二乘兩者計(jì)算結(jié)果并無(wú)太大差別。然而加權(quán)總體最小二乘單位權(quán)中誤差都有所提高,參數(shù)估值的精度更好。
圖2為加權(quán)總體最小二乘平面和二次多項(xiàng)式擬合殘差,從圖2可以看出最大殘差為6 cm,大多數(shù)在2 cm附近。圖3為傳統(tǒng)最小二乘與總體最小二乘、加權(quán)總體最小二乘擬合殘差的較差,可以看出總體最小二乘法對(duì)殘差提高不大,而加權(quán)總體最小二乘法可提高2 mm左右。
表1 平面擬合參數(shù)結(jié)果及精度(單位:m)Tab.1 Plane fitting results and their accuracy with some method(unit:m)
圖1 GPS水準(zhǔn)點(diǎn)位分布Fig.1 Distribution of GPS leveling points
圖2 WTLS擬合殘差Fig.2 Residuals of fitting with WTLS
圖3 LS與TLS、WTLS擬合結(jié)果較差Fig.3 Comparison among the fitting results with LS、TLS and WTLS
表2 二次曲面擬合參數(shù)結(jié)果及其精度(單位:m)Tab.2 Quadric surface fitting results and accuracy with three methods(unit:m)
傳統(tǒng)最小二乘方法在系數(shù)陣存在誤差的情況下,并不能很好地估計(jì)多項(xiàng)式系數(shù)參數(shù)??紤]到估計(jì)系數(shù)陣含有誤差和觀測(cè)變量不等精度,采用了更加合理的加權(quán)總體最小二乘方法。根據(jù)實(shí)例計(jì)算結(jié)果,加權(quán)總體最小二乘方法的單位權(quán)中誤差更小,檢核點(diǎn)殘差也相應(yīng)減小。在GPS高程擬合中,加權(quán)總體最小二乘法可以提高精度。
1 李征航,黃勁松.GPS測(cè)量與數(shù)據(jù)處理[M].武漢:武漢大學(xué)出版社,2010.(Li Zhenghang and Huang Jinsong.GPS surveying and data processing[M].Wuhan:Wuhan University Press,2010)
2 徐紹銓?zhuān)?GPS測(cè)量原理及應(yīng)用[M].武漢:武漢大學(xué)出版社,2008.(Xu shaoshuai,et al.The principle and application of GPS surveying[M].Wuhan:Wuhan University Press,2008)
3 邱衛(wèi)寧,等.測(cè)量數(shù)據(jù)處理理論與方法[M].武漢:武漢大學(xué)出版社,2008. (Qiu Weining,et al.The theory and method of surveying data processing[M].Wuhan:Wuhan University Press, 2008)
4 Chi Lun Cheng,Mastrondrdl N and Palge C.Total least squares and errors-in-variables modeling[J].Computational Statistics and Data Analysis,2007,52:1 076-1 079.
5 王樂(lè)洋,許才軍,魯鐵定.病態(tài)總體最小二乘模型的正則化算法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2010,35(11):1 346-1 350.(Wang Leyang,Xu Caijun and Lu Tieding.Ridge estimation method in Ill-posed weighted total least squares adjustment[J],Geomatices and Information Science of Wuhan University,2010,35(11):1 346-1 350)
6 Burkhard Schaffrin.A note on constrained total least-squares estimation[J].Linear Algebra and its Applications,2006,417:245-258.
7 Schaffrin B and Felus Y A.On the multivariate total leastsquares approach to empirical coordinate transformations.three algorithms[J].Journal of Geodesy,2008,82(6):373-383.
8 Schaffrin B and Wieser A.On weighted total least-square adjustment for linear regression[J].Journal of Geodesy,2008,82:415-421.
9 陳瑋嫻,等.加權(quán)總體最小二乘在三維激光標(biāo)靶擬合中的應(yīng)用[J].大地測(cè)量與地球動(dòng)力學(xué),2010,(5):90-96.(Chen Weixian,et al.Application of weighted total leastsquares to target fitting of three-dimensional laser scanning[J].Journal of Geodesy and Geodynamics,2010,(5):90 -96)
GPS HEIGHT FITTING OF WEIGHTED TOTAL LEAST-SQUARES ADJUSTMENT
Zhao Hui1,2),Zhang Shubi1)and Zhang Qiuzhao1,2)
1)School of Environment and Spatial Informatics,China University of Mining and Technology,Xuzhou 221008 2)Key Laboratory for Land Environment and Disaster Monitoring of SBSM,Xuzhou 221008
In GPS height fitting,a new method of Weighted Total Least-Squares adjustment(WTLS)is presented for solving the error of coefficient matrix.A more reasonable fitting model of plane and quadric polynomial is established,and the corresponding iterative algorithm is given.The examples of calculations show that the polynomial parameter is more reasonable and the solved height anomaly is more accurate.
GPS height;polynomial fitting;errors-in-variables model;least square;weighted total least-square
1671-5942(2011)05-0088-04
2011-03-14
江蘇省普通高校研究生科研創(chuàng)新計(jì)劃項(xiàng)目(CXLX11-0323)
趙輝,男,1987年生,碩士,研究方向:GNSS理論及應(yīng)用研究.E-mail:zhaohui@cumt.edu.cn
P207
A