總體最小二乘平差方法在GPS高程擬合中的應(yīng)用研究

龔循強(qiáng)，陳 磬，周秀芳

(1. 西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756； 2. 江西省數(shù)字國(guó)土重點(diǎn)實(shí)驗(yàn)室(東華理工大學(xué))，江西 南昌 330013； 3. 四川省應(yīng)急測(cè)繪保障與地質(zhì)災(zāi)害監(jiān)測(cè)工程技術(shù)研究中心，四川 成都 610041； 4. 江西核工業(yè)測(cè)繪院，江西 南昌 330038； 5. 四川農(nóng)業(yè)大學(xué) 經(jīng)濟(jì)管理學(xué)院，四川 成都 611130)

一、引 言

GPS高程擬合是將高程異常值與平面坐標(biāo)近似描述為多項(xiàng)式關(guān)系，利用同時(shí)已知大地高和正常高的公共點(diǎn)組成誤差方程，進(jìn)而求解多項(xiàng)式系數(shù)。目前，普遍采用最小二乘(LS)求解GPS高程擬合的多項(xiàng)式系數(shù)，該方法具有簡(jiǎn)單、方便等優(yōu)勢(shì)，但也同樣存在未顧及系數(shù)矩陣誤差這一客觀現(xiàn)象?；诖?，已有相關(guān)學(xué)者采用能夠同時(shí)顧及系數(shù)矩陣和觀測(cè)向量誤差的總體最小二乘[1-2](TLS)平差方法進(jìn)行參數(shù)估計(jì)，從而使得估計(jì)結(jié)果更加精確。文獻(xiàn)[3]采用基于迭代法的TLS求解參數(shù)估值，其實(shí)質(zhì)等同于基于奇異值分解的總體最小二乘(SVD-TLS)，該方法存在未考慮系數(shù)矩陣部分有誤差的情況，且現(xiàn)今較多文獻(xiàn)[4-6]對(duì)TLS的求解都是考慮系數(shù)矩陣中的全部誤差，而在測(cè)繪實(shí)踐中，經(jīng)常會(huì)遇到方程式系數(shù)矩陣部分有誤差的情況，這種情況需采用混合總體最小二乘(LS-TLS)求解。

近年來，關(guān)于GPS高程擬合方法的研究較多[7-8]，而對(duì)擬合模型的最佳解算方法卻較少涉及。本文在對(duì)TLS平差方法進(jìn)行理論研究的基礎(chǔ)上，針對(duì)LS和SVD-TLS存在對(duì)模型系數(shù)矩陣和觀測(cè)向量誤差估計(jì)不合理，造成估計(jì)結(jié)果不嚴(yán)密的現(xiàn)象，嘗試采用LS-TLS進(jìn)行GPS高程擬合，在合理考慮系數(shù)矩陣和觀測(cè)向量誤差的情況下，利用本文提出的基于TLS平差的粗差探測(cè)方法進(jìn)行粗差剔除后，對(duì)GPS高程擬合參數(shù)進(jìn)行求解。最后結(jié)合真實(shí)數(shù)據(jù)，分別利用LS、SVD-TLS和LS-TLS進(jìn)行GPS高程平面擬合，并對(duì)擬合結(jié)果進(jìn)行了比較分析。

二、總體最小二乘平差方法

目前，LS是對(duì)擬合函數(shù)參數(shù)求解的最常用方法，但這種方法不能顧及系數(shù)矩陣中的誤差。而TLS平差方法能夠同時(shí)顧及觀測(cè)向量和系數(shù)矩陣誤差，使得函數(shù)的解算結(jié)果更合理。

1. 基于奇異值分解的總體最小二乘

奇異值分解是線性代數(shù)中一種重要的矩陣分解，同時(shí)也是對(duì)TLS參數(shù)估計(jì)的常用方法，能較好地處理病態(tài)矩陣問題。在函數(shù)模型中同時(shí)顧及系數(shù)矩陣和觀測(cè)向量誤差，設(shè)矩陣方程為

(A+ΔA)X=L+ΔL

(1)

假設(shè)觀測(cè)次數(shù)為n，參數(shù)個(gè)數(shù)為m，通常情況下n>m，將增廣矩陣B=[AL]奇異值分解為B=UΣVT，SVD-TLS的解可由增廣矩陣右奇異向量的最后一列求得，進(jìn)一步可以計(jì)算得到殘差平方和及單位權(quán)中誤差。

2. 混合總體最小二乘

由式(1)可以看出，SVD-TLS求解參數(shù)的同時(shí)考慮了觀測(cè)向量和系數(shù)矩陣中的全部誤差。然而在測(cè)繪實(shí)踐中，經(jīng)常會(huì)遇到系數(shù)矩陣部分有誤差的情況，若采用SVD-TLS求解顯然是不嚴(yán)密的，目前解決這類問題的適宜方法是LS-TLS。假設(shè)A1沒有誤差影響，則顧及系數(shù)矩陣和觀測(cè)向量誤差，LS-TLS函數(shù)模型為

A1X1+(A2+ΔA2)X2=L+ΔL

(2)

3. 基于總體最小二乘平差的粗差探測(cè)方法

TLS平差是研究數(shù)據(jù)擬合的有效方法，其不僅考慮到觀測(cè)向量含有誤差，同時(shí)也考慮到系數(shù)矩陣含有誤差的情況，反映了觀測(cè)向量和系數(shù)矩陣的總體數(shù)據(jù)特征。因此，不管對(duì)于二維數(shù)據(jù)擬合還是對(duì)于多維數(shù)據(jù)擬合，總體數(shù)據(jù)的態(tài)勢(shì)都將在直線或平面一定范圍內(nèi)出現(xiàn)。對(duì)于總體觀測(cè)數(shù)據(jù)來說，如果有某一個(gè)觀測(cè)量含有粗差或異常值，則該觀測(cè)數(shù)據(jù)必定相對(duì)于其他數(shù)據(jù)離擬合的直線或平面的距離要大。

由于TLS平差所形成的擬合直線或擬合平面反映的是總體觀測(cè)數(shù)據(jù)的變化趨勢(shì)，因此，如果觀測(cè)數(shù)據(jù)中沒有粗差或異常值，基于TLS平差所形成的擬合殘差必定分布于某個(gè)誤差橢圓內(nèi)。若某個(gè)觀測(cè)數(shù)據(jù)擬合殘差超出了這個(gè)誤差橢圓，則可以認(rèn)為此觀測(cè)數(shù)據(jù)有可能含有粗差或異常值?；谠撍枷肟梢岳肨LS擬合殘差及極限誤差的相關(guān)特性探尋觀測(cè)數(shù)據(jù)中的粗差或異常值。下面以平面擬合為例進(jìn)行具體詳解，步驟如下：

1) 設(shè)有平面方程z=a+bx+cy，利用TLS平差對(duì)其進(jìn)行擬合，求出參數(shù)的估計(jì)值。

2) 求出中誤差σ，并取置信區(qū)間為(-2σ,2σ)。

3) 求出擬合殘差Vi，判定每一個(gè)擬合殘差是否在置信區(qū)間內(nèi)。

4) 若擬合殘差超出置信區(qū)間，則可認(rèn)為形成該擬合殘差的觀測(cè)數(shù)據(jù)中含有粗差，可以將其剔除；反之，則認(rèn)為數(shù)據(jù)中不存在粗差。

三、實(shí)例分析

以上簡(jiǎn)要闡述了采用LS、SVD-TLS和LS-TLS進(jìn)行參數(shù)估計(jì)及精度評(píng)定的方法，并給出了基于TLS平差的粗差探測(cè)方法。下面采用真實(shí)數(shù)據(jù)進(jìn)行GPS高程平面擬合，分別利用LS、SVD-TLS和LS-TLS進(jìn)行計(jì)算，并從模型、參數(shù)估值、殘差平方和及單位權(quán)中誤差等幾個(gè)方面對(duì)各自的GPS高程擬合結(jié)果進(jìn)行比較分析。

以文獻(xiàn)[11]例5.8中的GPS水準(zhǔn)擬合觀測(cè)數(shù)據(jù)作為樣本，并在數(shù)據(jù)中加入一定的粗差，觀測(cè)數(shù)據(jù)見表1，選取測(cè)站點(diǎn)的坐標(biāo)為自變量(x,y)，各點(diǎn)的高程異常值作為因變量δh，則有δh=a+bx+cy。

表1 加入粗差后的GPS水準(zhǔn)擬合觀測(cè)數(shù)據(jù)

由于該平面方程系數(shù)矩陣部分含有誤差，因此采用LS-TLS求解擬合參數(shù)，根據(jù)本文提出的基于TLS平差的粗差探測(cè)方法進(jìn)行粗差的探測(cè)與剔除。由參數(shù)估計(jì)值求得中誤差σ=1.695 8，即得置信區(qū)間為(-3.391 6,3.391 6)，并求出擬合殘差Vi，各測(cè)站點(diǎn)的擬合殘差如圖1所示。

判定每個(gè)測(cè)站點(diǎn)的擬合殘差是否在置信區(qū)間內(nèi)，對(duì)存在粗差的測(cè)站數(shù)據(jù)進(jìn)行剔除，根據(jù)圖1及求出的置信區(qū)間可將5、8、17、20測(cè)站數(shù)據(jù)剔除。同時(shí)采用LS、SVD-TLS和LS-TLS 3種方法對(duì)剔除粗差后的GPS水準(zhǔn)擬合數(shù)據(jù)進(jìn)行平面擬合，求出擬合參數(shù)、殘差平方和及單位權(quán)中誤差，結(jié)果見表2。

圖1 測(cè)站點(diǎn)的LS-TLS擬合殘差

求解方法abcS殘/m^σ0/mLS-144.95040.12000.19052.79500.4317SVD-TLS-147.93200.12250.19392.71000.4250LS-TLS-145.08030.12010.19072.66000.4211

表2中分別采用LS、SVD-TLS和LS-TLS求解了擬合平面參數(shù)，由于TLS平差方法同時(shí)顧及了系數(shù)矩陣x、y和觀測(cè)向量δh的誤差，因此建立的模型更加符合實(shí)際情況，其模型的殘差平方和及單位權(quán)中誤差優(yōu)于LS。

從平面方程δh=a+bx+cy可以看出，部分系數(shù)矩陣和觀測(cè)向量都是由GPS數(shù)據(jù)組成的矩陣，采用LS進(jìn)行平面擬合，其只顧及了觀測(cè)向量δh的誤差，因此擬合的平面并不精確；采用SVD-TLS進(jìn)行平面擬合，所得擬合平面與LS和LS-TLS的擬合平面相差較大，通過分析構(gòu)建的函數(shù)模型和解算方法發(fā)現(xiàn)，SVD-TLS將系數(shù)矩陣中不存在誤差的常數(shù)項(xiàng)部分也考慮了進(jìn)去，因而造成擬合結(jié)果變異；LS-TLS合理地顧及了系數(shù)矩陣x、y和觀測(cè)向量δh的誤差，因此建立的模型更加合理，所得擬合精度最優(yōu)，較LS和SVD-TLS擬合平面結(jié)果更加接近實(shí)際，故可認(rèn)為合理考慮系數(shù)矩陣和觀測(cè)向量誤差的LS-TLS擬合平面結(jié)果最為可靠。

四、結(jié) 論

本文采用TLS平差方法進(jìn)行GPS高程擬合，在利用提出的基于TLS平差的粗差探測(cè)方法進(jìn)行粗差剔除的基礎(chǔ)上，結(jié)合真實(shí)數(shù)據(jù)，進(jìn)行擬合平面的參數(shù)計(jì)算及精度評(píng)定。通過LS、 SVD-TLS和LS-TLS的比較分析研究，得出如下相關(guān)結(jié)論：

1) 采用本文提出的基于TLS平差的粗差探測(cè)方法能較好地對(duì)數(shù)據(jù)中存在的粗差進(jìn)行探測(cè)與剔除。

2) 在GPS高程擬合參數(shù)估計(jì)中，由于LS沒有考慮系數(shù)矩陣存在誤差的情況，因此擬合結(jié)果并不嚴(yán)密。SVD-TLS 將系數(shù)矩陣中不含誤差的常數(shù)矩陣進(jìn)行考慮，造成擬合結(jié)果變異。

3) 采用LS-TLS進(jìn)行GPS高程擬合參數(shù)估計(jì)，它同時(shí)合理地顧及系數(shù)矩陣和觀測(cè)向量誤差，所得擬合精度最優(yōu)，因此可認(rèn)為L(zhǎng)S-TLS擬合平面的結(jié)果最為可靠。

參考文獻(xiàn)：

[1] GOLUB G H, VAN LOAN C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal Numerical Analysis, 1980, 17(6): 883-893.

[2] VAN HUFFEL S，LEMMERLING P. Total Least Squares and Error-in-variables Modeling: Analysis, Algorithms and Applications [M]. Netherlands:Springer,2002.

[3] 趙輝, 張書畢, 張秋昭. 基于加權(quán)總體最小二乘法的GPS高程擬合 [J]. 大地測(cè)量與地球動(dòng)力學(xué), 2011, 31(5): 88-90.

[4] 龔循強(qiáng), 李通, 陳西江. 總體最小二乘法在曲線擬合中的應(yīng)用[J]. 地礦測(cè)繪, 2012, 28(3): 4-6.

[5] 官云蘭, 周世健, 張立亭，等. 穩(wěn)健整體最小二乘直線擬合[J]. 工程勘察, 2012,40(2): 60-62.

[6] 官云蘭, 劉紹堂, 周世健，等. 基于整體最小二乘的穩(wěn)健點(diǎn)云數(shù)據(jù)平面擬合 [J]. 大地測(cè)量與地球動(dòng)力學(xué), 2011, 31(5): 80-83.

[7] 叢康林, 岳建平. 基于SVR的GPS高程擬合模型研究[J]. 測(cè)繪通報(bào), 2011(2): 8-11.

[8] 朱東銘, 秦偉, 范百興. GPS橋梁控制網(wǎng)高程擬合及其精度分析[J]. 測(cè)繪通報(bào), 2009(9): 42-44.

[9] VAN HUFFEL S, VAN DEWALLE J. The Total Least Squares Problem：Computational Aspects and Analysis [M].Philadelphia: Society for Industrial and Applied Mathematics,1991.

[10] GOLUB G H, HOFFMAN A, STEWART G W. A Generalization of the Eckart-Young-Mirsky Matrix Approximation Theorem [J]. Linear Algebra and Its Applications, 1987, 88/89: 317-327.

[11] 邱衛(wèi)寧, 陶本藻, 姚宜斌，等. 測(cè)量數(shù)據(jù)處理理論與方法[M]. 武漢: 武漢大學(xué)出版社, 2008.