基于丹麥法計算權(quán)因子的GPS高程擬合方法研究

張艷兵,田林亞,畢繼鑫,王永敏

(河海大學(xué) 地球科學(xué)與工程學(xué)院,江蘇 南京 211100)

基于丹麥法計算權(quán)因子的GPS高程擬合方法研究

張艷兵,田林亞,畢繼鑫,王永敏

(河海大學(xué) 地球科學(xué)與工程學(xué)院,江蘇 南京 211100)

針對GPS高程擬合這一問題,介紹了等值線法、解析擬合法等傳統(tǒng)GPS高程擬合方法,并選用丹麥法計算權(quán)因子,構(gòu)造等價權(quán)矩陣,結(jié)合總體最小二乘法抵抗粗差點對GPS高程擬合結(jié)果的影響。結(jié)果表明:等值線法與其它兩種方法對不含粗差的等精度觀測數(shù)據(jù)擬合的效果相當;基于丹麥法計算權(quán)因子的抗差總體最小二乘法對不等精度觀測數(shù)據(jù)的偶然誤差處理效果更佳,且通過權(quán)因子可以有效處理粗差點對高程擬合的影響,擬合效果更佳。

GPS高程擬合;等值線法;解析擬合法;丹麥法;抗差總體最小二乘法

0 引 言

傳統(tǒng)高程測量是指利用水準儀進行測量,獲取所需高程,即正常高[1]。隨著現(xiàn)代化進程的推進,社會建設(shè)的規(guī)模越來越大,對高程系統(tǒng)快速獲取要求更高,使得傳統(tǒng)測量手段在大范圍、地形起伏較大、通視性不良的區(qū)域難以實施。而GPS的高精度定位和快速獲取使其得到了廣泛關(guān)注,結(jié)合水準測量體系,形成了GPS高程測量理論[2],降低了水準測量在地形起伏較大區(qū)域的作業(yè)難度。

GPS直接獲取的高程是大地高,而在具體工程應(yīng)用中所使用的高程系統(tǒng)是正常高,兩者的參考面并不是同一基準面[3],所以GPS高程要在獲取的大地高基礎(chǔ)上,減去兩個參考面的差異值。由于地球密度分布不均勻,使得似大地水準面成為一個連續(xù)變化的不規(guī)則面,進而導(dǎo)致高程異常值的不規(guī)則變化,難于準確獲取[4-5]。如圖1所示,由于θ很小,所以控制點的高程異常值ξi可以根據(jù)ξi=Hi-hi直接獲取,然后采用某種方法建立控制點高程異常值與待求點高程異常值的相互關(guān)系,最后根據(jù)這種相互關(guān)系確定待求點的高程異常值估值。因此,如何精確測定控制點與待求點的相互關(guān)系是提高GPS高程測量精度的一個關(guān)鍵問題。

圖1 觀測點大地高與正常高的示意圖

1 等值線法

由于高程異常值隨著平面坐標不斷變化,故可以將其轉(zhuǎn)化為三維模型,根據(jù)繪制等高線原理,繪制該區(qū)域高程異常值的等值線畫圖,然后根據(jù)待求點平面坐標在等值線畫圖中拾取高程異常值估值。

具體過程為首先在測區(qū)選取適當?shù)腉PS控制網(wǎng),但需要注意:1) GPS控制點需均勻分布;2) GPS控制點要布設(shè)在測區(qū)最高點、最低點、地性線上,使其可以充分反映測區(qū)范圍內(nèi)高程異常面的起伏變化;3) 在地勢起伏較大的區(qū)域,需加密GPS控制點的數(shù)量。其次選取其中m個控制點作為已知點,并對其進行水準聯(lián)測,獲取已知點的高程異常值,其它控制點作為檢核點,然后利用matlab畫出該區(qū)域高程異常值的等值線畫圖,并在圖上標出檢核點的位置。最后利用matlab中點拾取工具,拾取檢核點的高程異常值信息,將其作為檢核點的高程異常值估值。

2 解析擬合法

解析擬合法,就是將已知點的高程異常值與檢核點的高程異常值的相關(guān)關(guān)系用數(shù)學(xué)知識精確表達為函數(shù)模型,然后根據(jù)函數(shù)模型計算檢核點的高程異常值估值。常見的解析擬合法有多項式擬合法和抗差總體最小二乘法。

2.1多項式擬合法

多項式擬合法的基本原理是首先利用GPS測量數(shù)據(jù)與水準聯(lián)測高程計算出已知點的高程異常值ξi,然后建立高程異常值ξi與平面坐標x,y的函數(shù)模型f(x,y),求得函數(shù)模型參數(shù),最后根據(jù)函數(shù)模型與待求點的平面坐標,計算待求點的高程異常值估值。由于二次多項式擬合模型[6-9]精度較高,且所需的已知點數(shù)量較少,在實際應(yīng)用中得到了廣泛使用,故本文將采用二次多項式模型擬合GPS高程異常值。

設(shè)已知點的高程異常ξ與其平面坐標x,y之間關(guān)系為

ξ=f(x,y)+ε，

(1)

式中:ξ為高程異常值;f(x,y)為ξ的近似表達式;ε為誤差值。

二次多項式模型可表示為

f(x,y)=a0+a1x+a2y+a3xy+

a4x2+a5y2，

(2)

寫成矩陣形式:

ξ=BX+ε，

(3)

改寫為誤差方程:

V=BX-l.

(4)

要使擬合曲面最佳,只需滿足∑vivi=min,即:

VTPV=min .

(5)

在聯(lián)測的n個控制點中,選取m(m

2.2抗差總體最小二乘法

系數(shù)矩陣A和觀測向量ξ都是根據(jù)觀測數(shù)據(jù)獲得,當觀測數(shù)據(jù)誤差較大甚至含有粗差時,A和ξ也會受到很大影響,使平差結(jié)果出現(xiàn)較大偏差,擬合效果不佳。而抗差總體最小二乘法針對這一特性,引入權(quán)因子,優(yōu)化了平差模型,降低了觀測數(shù)據(jù)中粗差對擬合效果的影響。為便于比較抗差總體最小二乘法與多項式擬合法的區(qū)別,繼續(xù)采用二次多項式構(gòu)建總體最小二乘模型。

采用總體最小二乘模型[10],其誤差方差為

ξ+eξ=(A+EA)X,

(6)

式中:A∈Rn×m,ξ∈Rn,X∈Rm,EA為系數(shù)陣A的誤差,eξ為觀測向量ξ的觀測誤差。

若令eA=vec(EA),QA和Qξ分別為eA和eξ的協(xié)因數(shù)陣,則有:

(7)

要求式(6)的最佳估計值,需要滿足:

(8)

由于觀測值中含有粗差,所以要對總體最小二乘模型進行改正,本文通過丹麥法計算權(quán)因子[11],抵制觀測值中粗差的影響。

(9)

計算步驟:

1) 計算初始值;

v0=0,x0=(ATPξ0A)-1ATPξ0ξ.

(10)

2) 迭代過程,計算中間值;

Qξi=QξI-1ωi,QAi=QAi-1ωi,

(11)

(12)

(13)

3) 判斷迭代終止條件為

‖Xi+1-Xi‖<ε.

(14)

根據(jù)式(14)判斷迭代過程是否結(jié)束,若不成立,重復(fù)迭代;若成立,迭代結(jié)束,計算單位權(quán)中誤差,其估值為

(15)

式中:r為自由度;λ=(Qξ+(XTQ0X)QA)-1(ξ-AX)。

3 算例分析

為了驗證抗差總體最小二乘法在GPS高程擬合中的有效性,選取云南某地采集數(shù)據(jù)與某城市控制網(wǎng)數(shù)據(jù)進行GPS高程異常值擬合計算。如圖2所示,數(shù)據(jù)1:云南某地區(qū)共布設(shè)了27個GPS控制點,GPS網(wǎng)為C級網(wǎng),控制點全部為三等水準點,剔除3個帶有明顯粗差的控制點,選取其中的2、6、12、14號點作為檢核點,其它點作為已知點。數(shù)據(jù)2:某中型城市的城市控制網(wǎng),共已知20個GPS-E級控制點,為研究GPS擬合原理,并對以上控制點進行了3等水準測量。選取其中的2、15、20號點作為檢核點,其它點作為已知點。

圖2 兩組實驗數(shù)據(jù)點位空間分布圖

3.1等值線法

由兩組實驗數(shù)據(jù)分別生成間距為10 cm和20 cm的等值線畫圖,如圖3所示,然后根據(jù)檢核點平面坐標在圖3上拾取其高程異常值估值,結(jié)果如表1、表2所示。

圖3 兩組實驗數(shù)據(jù)生成的等值線畫圖

點號X坐標/mY坐標/m高程異常值/m高程異常值估值/m偏差/m22763781.594573220.44032.61632.5970.01962768522.200593641.08432.03232.064-0.032102766456.717613160.71432.04832.135-0.103122776204.448553188.40332.47132.2430.228

表2 數(shù)據(jù)2中檢核點高程異常值計算結(jié)果

從表1與表2中可以看出,由等值線法內(nèi)插出的檢核點高程異常值與真值(這里把聯(lián)測的水準高程與大地高高差作為檢核點高程異常值的真值)之差最小可以達到cm級精度,最大可以達到厘米級,所以用等值線法來內(nèi)插高程異常值的結(jié)果并不穩(wěn)定,且要求內(nèi)插點盡可能分布在中心位置,同時,該方法也缺少精度評價指標,不易判斷生成的等值線正確與否,也不能判斷內(nèi)插點的高程異常值正確與否。

3.2解析擬合法

由于等值線法缺少精度衡量指標,很難判斷檢核點內(nèi)插出的高程異常值是否正確,因此,引入多項式擬合法與抗差總體最小二乘法來對兩組數(shù)據(jù)進行擬合,并在數(shù)據(jù)1中添加粗差點(2 762 850.290 8,566 913.499 9,32.324 8)與(2 777 857.3 74 6,594 527.091 8,32.276 4),分析多項式擬合法與抗差總體最小二乘法的擬合效果。

利用多項式擬合法與抗差總體最小二乘法分別對數(shù)據(jù)1與數(shù)據(jù)2中的已知點進行擬合處理,得到各點高程異常值ξi的殘差值vi,如圖4所示。其中,圖4(a)是利用多項式擬合法對數(shù)據(jù)1的處理結(jié)果,從中可以看出粗差點14,26不能得到很好的處理,對其它已知點的擬合結(jié)果造成了一定程度上的影響;圖4(b)是利用抗差總體最小二乘法對數(shù)據(jù)1的處理結(jié)果,從中可以粗差點對其它已知點的擬合結(jié)果影響較小;圖4(c)、4(d)示出了利用兩種方法對數(shù)據(jù)2的處理結(jié)果,可以看出兩種方法擬合效果相差不大。

圖4 兩組實驗數(shù)據(jù)高程異常值擬合殘差 (a)多項式擬合法擬合殘差(數(shù)據(jù)1);(b) RWTLS法擬合殘差(數(shù)據(jù)1);(c)多項式擬合法擬合殘差(數(shù)據(jù)2);(d) RWTLS法擬合殘差(數(shù)據(jù)2)

表3 兩組實驗數(shù)據(jù)中檢核點在不同擬合方法中高程異常值差值計算結(jié)果

根據(jù)已知點的殘差值與檢核點的高程異常值差值進行精度分析,其中,內(nèi)符合精度計算公式為

(16)

式中,v為已知點殘差值向量,r為自由度。

外符合精度計算公式為

(17)

式中： Δ為檢核點的高程異常值真值與擬合值之差；n為檢核點的個數(shù)。

表4示出了兩組實驗數(shù)據(jù)在不同擬合方法下的精度。

表4 兩組實驗數(shù)據(jù)在不同擬合方法下的計算精度

4 結(jié)束語

本文采用等值線法、多項式擬合法、抗差總體最小二乘法分別對數(shù)據(jù)1與數(shù)據(jù)2進行GPS高程擬合。當觀測數(shù)據(jù)不含粗差,且觀測數(shù)據(jù)質(zhì)量較高時,可以看出多項式擬合法與總體最小二乘擬合法都能取得很好的擬合效果,同時,通過檢核點在三種方法下擬合結(jié)果的比較,可以得到等值線法內(nèi)插出的未知點高程異常值精度也很高,且該方法計算方便,便于提取區(qū)域中的任意點的高程異常值,此時等值線法更加適用;當觀測數(shù)據(jù)中含有少量粗差時,前兩種方法的擬合精度有所下降,且越靠近粗差點,擬合效果越差,而抗差總體最小二乘法擬合效果依然可以保持穩(wěn)定,綜上,抗差總體最小二乘法適用范圍更廣,擬合效果也比較穩(wěn)定。

擬合計算中,本文只采用丹麥法來計算權(quán)因子,且取c=σ0,并未就c的取值對抗差效果進行分析;也未引入其他模型與丹麥法計算權(quán)因子進行比較,分析其對抗差總體最小二乘法擬合效果的影響,有待進一步的研究。

[1] 孔祥元. 大地測量學(xué)基礎(chǔ)[M]. 武漢:武漢大學(xué)出版社, 2010.

[2] 陳為民,張旭東,符華年,等. GPS高程測量代替等級水準測量的應(yīng)用研究[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版),2013(7):828-831.

[3] 王珍,程壘. GPS高程擬合方法的研究和應(yīng)用[J]. 測繪通報, 2015(s1):102-105.

[4] 田曉,鄭洪艷,許明元,等. 一種改進的適用于不同地形的GPS高程擬合模型[J]. 測繪通報, 2017(1):35-38,64.

[5] 吳迪軍,熊偉. 跨海橋梁GPS高程擬合方法[J]. 測繪科學(xué),2017(6):1-6.

[6] 張潘,余代俊,張玉剛,等. GPS高程擬合方法研究及精度對比試驗[J]. 測繪通報, 2015(9):54-56.

[7] 李宏博,史先琦,萬奇靈,等. GPS高程擬合在丘陵地區(qū)地形測量中的應(yīng)用研究[J]. 測繪地理信息, 2015, 40(4):44-47.

[8] 劉斌,郭際明,史俊波,等. 利用EGM2008模型與地形改正進行GPS高程擬合[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2016, 41(4):554-558.

[9] 王昶,王旭. GPS高程擬合分區(qū)方法研究[J]. 大地測量與地球動力學(xué), 2016, 36(1):26-29.

[10] 劉亞彬,鄭南山,張旭,等. GPS高程擬合的加權(quán)總體最小二乘抗差估計[J]. 大地測量與地球動力學(xué), 2016, 36(1):30-34.

[11] 林國標,劉立龍,蔡成輝,等. 一種基于丹麥法的改進型雙步M估計[J]. 大地測量與地球動力學(xué),2015(2):235-238，247.

ResearchonGPSHeightFittingMethodBasedonDanishMethod

ZHANGYanbing,TIANLinya,BIJixin,WANGYongmin

(CollegeofEarthScienceandEngineering,HohaiUniversity,Nanjing211100，hina)

In this paper, the traditional GPS height fitting methods, such as contour method and analytic fitting method, are introduced for GPS elevation fitting. The Danish method is used to calculate the weight factor for the first time to construct the equivalent weight matrix, and combined with the overall least squares method to resist the influence of the rough point on the GPS elevation fitting results. The results show that the effect of the contour method is equivalent to the other two methods for fitting the equal precision observation data without gross error. Based on the weight factor calculated by Danish method, The deal result of the robust overall least squares method is better in dealing with the accidental error of the unequal accuracy observation data,and the effect of the rough points on the elevation fitting can be effectively dealt by using the weighting factor, and the fitting effect is better.

GPS elevation fitting; contour method; analytic fitting method; Danish method; robust overall least squares method

10.13442/j.gnss.1008-9268.2017.04.017

P228.4

A

1008-9268(2017)04-0096-06

2017-04-06

聯(lián)系人: 張艷兵E-mail:1821708477@qq.com

張艷兵(1992-),男,碩士研究生，研究方向為測量數(shù)據(jù)處理理論與方法。

田林亞(1963-),男,教授,博士,研究方向為精密工程測量理論與技術(shù)。

畢繼鑫(1994-),男,碩士研究生,研究方向為測量數(shù)據(jù)處理理論與方法。