抗差推估GNSS水準(zhǔn)高程擬合方法研究

郭范春

（遼寧省測(cè)繪產(chǎn)品質(zhì)量監(jiān)督檢驗(yàn)站，遼寧 沈陽(yáng) 110034）

1 概述

GNSS水準(zhǔn)高程擬合所用的數(shù)據(jù)不可避免的包含一定的誤差，誤差產(chǎn)生的原因是多方面的同時(shí)也包含多方面的粗差。清楚認(rèn)識(shí)誤差對(duì)擬合結(jié)果的影響規(guī)律，找到探測(cè)粗差的方法，采取必要措施將其削弱或消除，提高成果的可靠性和精確性十分必要。

目前高程擬合技術(shù)主要有整體擬合、隨機(jī)插值、分區(qū)擬合、逐點(diǎn)內(nèi)插等方法，這些方法對(duì)待定點(diǎn)周圍的參考點(diǎn)質(zhì)量和數(shù)量有嚴(yán)格的要求。由于參考點(diǎn)大地高、待定點(diǎn)大地高及其平面坐標(biāo)這些數(shù)據(jù)本身的誤差大小不盡一致，已知數(shù)據(jù)由于非同一年代實(shí)測(cè)、非同網(wǎng)平差或等諸多原因，致使擬合時(shí)在局部范圍內(nèi)這些數(shù)據(jù)不能滿足要求或兼容，加之有些測(cè)區(qū)范圍較大，或地形起伏較大，這種情況下就對(duì)高程擬合結(jié)果影響更為嚴(yán)重［1－6］。某些數(shù)據(jù)對(duì)于待定點(diǎn)擬合影響較大，變成了粗差，必須剔除或削弱這些點(diǎn)對(duì)周圍待定點(diǎn)的影響。為此，論文提出抗差推估方法結(jié)合其它擬合方法進(jìn)行GNSS水準(zhǔn)擬合計(jì)算，從而提高擬合結(jié)果的質(zhì)量和可靠性。

2 抗差推估GNSS水準(zhǔn)高程擬合

2.1 最小二乘推估擬合模型

設(shè)GNSS點(diǎn)高程異常ζ與平面坐標(biāo)（x，y）有如下式所示的函數(shù)關(guān)系

式中，f（x，y）為ζ中趨勢(shì)值，ε為誤差。

在GNSS水準(zhǔn)高程擬合時(shí)，按常規(guī)擬合方法，則上式函數(shù)關(guān)系的矩陣形式為

因?yàn)閷?duì)于每一個(gè)點(diǎn)都可以列立這樣的方程，在∑ε2＝min條件下，利用最小二乘理論解出式（2）方程組，再按式（1）求出待求點(diǎn)的ζ。

由大地高H（GNSS測(cè)得）與正常高Hr以及高程異常ζ的關(guān)系式H＝Hr＋ζ，從而可以求出正常高Hr。

高程異常是一隨機(jī)變量，若能知道觀測(cè)區(qū)域內(nèi)高程異常分布的隨機(jī)統(tǒng)計(jì)信息，即已知參考點(diǎn)間高程異常的方差－協(xié)方差陣

以及參考點(diǎn)與待定點(diǎn)間高程異常的方差－協(xié)方差陣

就可采用最小二乘推估法確定其趨勢(shì)值。其推估式為

式（3）、式（4）及式（5）中，ζxi（i＝1，…，m）為m個(gè)己知點(diǎn)上的高程異常，ζpj（j＝1，…，n）為n個(gè)待定點(diǎn)上的高程異常，δi2、δij為ζxi的方差、ζxi與ζxj間的協(xié)方差，δpjxi為ζpj與ζxi間的協(xié)方差。Qxx、Qpx的確定借助于觀測(cè)區(qū)域內(nèi)高程異常協(xié)方差函數(shù)而建立。

如果能確定描述高程異常隨機(jī)信息的協(xié)方差函數(shù)，最小二乘推估擬合就是首選方法，比一般擬合法要精確，而且已知點(diǎn)的分布沒(méi)有具體要求。但要建立較為精確的協(xié)方差函數(shù)比較困難。為此采用抗差推估方法進(jìn)行高程擬合。

2.1 抗差推估擬合

抗差推估擬合法是在原最小二乘擬合推估法的基礎(chǔ)上發(fā)展起來(lái)的一種選權(quán)迭代法［9］。它通過(guò)驗(yàn)后方差估計(jì)求出觀測(cè)值的驗(yàn)后方差。然后利用方差檢驗(yàn)找出方差異常大的觀測(cè)值，根據(jù)方差與權(quán)成反比的關(guān)系，給它一個(gè)相應(yīng)小的權(quán)，進(jìn)行下一步的迭代平差計(jì)算，重復(fù)以上過(guò)程，使含有粗差的觀測(cè)值的權(quán)越來(lái)越小甚至于等于0，從而使其對(duì)平差結(jié)果的影響很小，在某種意義上說(shuō)，當(dāng)觀測(cè)值的權(quán)很小或者等于0時(shí)，也就相當(dāng)于從觀測(cè)序列中剔除了該觀測(cè)值，因此在迭代的過(guò)程中，逐步發(fā)現(xiàn)粗差并將其 “剔除”。對(duì)于GPS高程擬合推估，具體的計(jì)算步驟是首先進(jìn)行最小二乘平差，即傳統(tǒng) GNSS高程擬合［10－12］。在本文研究中，式（1）的數(shù)學(xué)模型代換的是多項(xiàng)式曲面擬合模型，如下式

第一次抗差時(shí)，初權(quán)為1，然后根據(jù)式（7）進(jìn)行驗(yàn)后方差估計(jì)：

其中ri為多余觀測(cè)分量，ri＝（Qvvp）ii。

根據(jù)驗(yàn)后方差計(jì)算等價(jià)權(quán)，計(jì)算等價(jià)權(quán)的方法很多，如Huber函數(shù)：

c0一般取1.5～2.0，σ0是單位權(quán)標(biāo)準(zhǔn)差。

或者取指數(shù)函數(shù)為

重復(fù)以上步驟，隨著迭代次數(shù)的增加，如果某一觀測(cè)值的權(quán)變得越來(lái)越小甚至為0，這就說(shuō)明該觀測(cè)值含有粗差，或者與其它觀測(cè)值不相容，不能納入己知點(diǎn)列用于計(jì)算擬合函數(shù)系數(shù)；實(shí)際上，當(dāng)權(quán)降低到很小時(shí)，其對(duì)平差結(jié)果的影響也就很小，甚至沒(méi)有影響，從而計(jì)算結(jié)果抗差化。

迭代次數(shù)的增加，并非無(wú)限次循環(huán)下去，這里給出判定標(biāo)準(zhǔn)。以多項(xiàng)式相關(guān)平面數(shù)學(xué)模型說(shuō)明，如式（11）所示

一般情況下，設(shè)參考點(diǎn)的數(shù)目為n（n大于5），所以要進(jìn)行最小二乘法解算。判定的第一個(gè)標(biāo)準(zhǔn)是，在迭代到k次時(shí)，根據(jù)精度評(píng)定，可以計(jì)算出中誤差σ0，如果

式中μ為工程精度要求，σ0為擬合后的參考點(diǎn)精度。式（12）成立，則表明抗差滿足要求，迭代完畢。其二，當(dāng)在迭代到k次時(shí)，隨著權(quán)不斷變小，某些權(quán)可能為0，如果此時(shí)權(quán)不為0的參考點(diǎn)的個(gè)數(shù)：

即n′小于式（11）所示方程式中的系數(shù)個(gè)數(shù)，則方程不可解。這時(shí)可以降低工程精度的要求，重新解算。

具體算法實(shí)現(xiàn)時(shí)，可以在成果文件中查看每一次迭代的結(jié)果，包括權(quán)、誤差、精度和標(biāo)準(zhǔn)差，這樣就可以根據(jù)具體的要求，考察每一個(gè)參考點(diǎn)的精確性和可靠性，以作為排除點(diǎn)的依據(jù)。

3 工程應(yīng)用及分析

為了考察抗差推估GNSS水準(zhǔn)高程擬合方法的可行性和準(zhǔn)確性，本文采用某GPS高程控制網(wǎng)進(jìn)行計(jì)算分析。該GPS高程控制網(wǎng)所覆蓋的區(qū)域東西長(zhǎng)12km，南北長(zhǎng)18km，總面積約為214km2。圖1為全部GPS點(diǎn)平面分布情況，圖2為27個(gè)GPS試驗(yàn)點(diǎn)分布情況圖。GPS測(cè)點(diǎn)高程異常最大值為6.816m，最小值為6.410m；同時(shí)，水準(zhǔn)高差最大值為6.115m，最小值為3.376m。

3.1 一次抗差推估結(jié)果分析

根據(jù)前面的理論和通過(guò)編程實(shí)現(xiàn)了相關(guān)平面、二次曲面、三次曲面三種方法抗差推估GNSS水準(zhǔn)高程擬合。表1給出了全部控制點(diǎn)一次抗差推估的結(jié)果，其中差值是擬合后的GPS水準(zhǔn)高與幾何水準(zhǔn)高的之差。

圖1 全部GPS點(diǎn)平面分布圖

圖2 27個(gè)GPS試驗(yàn)點(diǎn)分布圖

由表1可以看出，二次曲面抗差結(jié)果符合點(diǎn)較少，且接近0的點(diǎn)較多。但是從超出誤差限的個(gè)數(shù)統(tǒng)計(jì)來(lái)看，相關(guān)平面不符合點(diǎn)個(gè)數(shù)最多，二次曲面卻最少。而且發(fā)現(xiàn)，某些點(diǎn)經(jīng)過(guò)抗差后，三種結(jié)果有不同的差值，例如B972、B842。對(duì)此作為研究，將它們保留。由多項(xiàng)式的特性，再根據(jù)該市地勢(shì)平坦特點(diǎn)，作為檢測(cè)，綜合分析三種抗差推估結(jié)果，將誤差超過(guò)5.0cm或權(quán)接近0的點(diǎn)認(rèn)為含有粗差，剔除掉。由此找出14個(gè)點(diǎn)符合要求，剔除掉后剩下37個(gè)點(diǎn)符合要求。

3.2 二次抗差推估結(jié)果分析

在以上計(jì)算分析的基礎(chǔ)上，進(jìn)行二次抗差推估，其結(jié)果見(jiàn)表2。表2列出了相關(guān)平面和二次曲面抗差的結(jié)果。根據(jù)第一次剔除誤差的方法，綜合分析，再次剔除粗差點(diǎn)10個(gè)，這樣剩下27個(gè)符合要求，它們也是以后進(jìn)行擬合分析的試驗(yàn)點(diǎn)，分布如圖2所示。

3.3 最終結(jié)果分析

本文再將控制點(diǎn)權(quán)接近1或誤差在1cm以下，結(jié)合圖形分布選擇以下首級(jí)控制點(diǎn)，它們?yōu)椋築344、B833、I407、I409、I412、I415、I417、I426、I429、I428共10個(gè)點(diǎn)，其余17個(gè)點(diǎn)則作為待定點(diǎn)（擬合點(diǎn)）。進(jìn)行了計(jì)算。

表1 三種整體抗差推估方法比較

續(xù)表

表2 二次整體抗差推估結(jié)果

續(xù)表

4 結(jié)論

（1）在高程控制網(wǎng)中，已知高程點(diǎn)數(shù)目比較多，且分布情況不理想，在進(jìn)行GPS水準(zhǔn)高程擬合時(shí)，先進(jìn)行抗差計(jì)算，剔除個(gè)別含有粗差的點(diǎn)能有效的提高高程擬合精度。

（2）在一定的條件下，二次曲面抗差的結(jié)果精度要好于相關(guān)平面與三次曲面抗差的結(jié)果精度，經(jīng)過(guò)二次曲面抗差后符合要求的點(diǎn)數(shù)最多，說(shuō)明二次曲面抗差模型的有效性比其它兩種模型要好。

（3）抗差模型的精度與點(diǎn)位分布有關(guān)系，個(gè)別含有粗差的點(diǎn)經(jīng)過(guò)三種模型進(jìn)行抗差計(jì)算后，結(jié)果各不相同，這要求在進(jìn)行抗差計(jì)算時(shí)要結(jié)合多種模型進(jìn)行比較分析，特別點(diǎn)要進(jìn)行多次抗差計(jì)算。

［1］ 李德仁、袁修孝.誤差處理與可靠性理論［M］.武漢：武漢大學(xué)出版社，2002：235－294.

［2］ SHRESTHA R L，NAZIR A，DEWITT B A，etal.Surface Interpolation Techniques to Convert GPS Ellipsoid Heights to Elevations［J］.Surveying and Land Information Systems，1993，53（3）：133－144.

［3］ 李建成.最新中國(guó)陸地?cái)?shù)字高程基準(zhǔn)模型：重力似大地水準(zhǔn)面CNGG2011［J］.測(cè)繪學(xué)報(bào)，2012，41（5）：651－660.

［4］ 蔣 濤.利用航空重力測(cè)量數(shù)據(jù)確定區(qū)域大地水準(zhǔn)面［J］.測(cè)繪學(xué)報(bào)，2013，42（1）：152.

［5］ 李建成.我國(guó)現(xiàn)代高程測(cè)定關(guān)鍵技術(shù)若干問(wèn)題的研究及進(jìn)展［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2007，32（11）：980－987.

［6］ 馬洪濱，董仲宇.基于DEM與重力場(chǎng)模型的GPS水準(zhǔn)高程處理方法研究［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2007，32（9）：767－769.

［7］ 章傳銀，黨亞民，柯寶貴，等.高精度海岸帶重力似大地水準(zhǔn)面的若干問(wèn)題討論［J］.測(cè)繪學(xué)報(bào)，2012，41（5）：709－714.

［8］ 陳俊勇，李建成，寧津生，等.全國(guó)及部分省市地區(qū)高精度、高分辨率大地水準(zhǔn)面的研究及其實(shí)施［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2006，34（4）：283－288.

［9］ 馬洪濱，董仲宇.多面函數(shù)GPS水準(zhǔn)高程擬合中光滑因子求定方法［J］.東北大學(xué)學(xué)報(bào)：自然科學(xué)版，2008，29（8）：1176－1178.

［10］ 陳俊勇，李健成，寧津生，等.我國(guó)大陸高精度、高分辨率大地水準(zhǔn)面的研究和實(shí)施［J］.測(cè)繪學(xué)報(bào)，2001，30（2）：95－100.

［11］ 郭東美，許厚澤.應(yīng)用GPS水準(zhǔn)與重力數(shù)據(jù)聯(lián)合解算大地水準(zhǔn)面［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（5）：621－624.

［12］ 申文斌，韓建成.利用新方法確定的30′×30′全球大地水準(zhǔn)面模型及其精度評(píng)估［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2012，37（10）：1135－1139.