王新 ，張紹良

(1.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇徐州 221008;2.江蘇省資源環(huán)境信息工程重點實驗室，江蘇徐州 221008)

1 引言

與常規(guī)水準(zhǔn)測量相比，GPS水準(zhǔn)具有費用低、效率高、觀測時間短、操作方便、全天候等特點，因此在為研究地表形變的變化機(jī)理提供動態(tài)監(jiān)測數(shù)據(jù)有著獨特的優(yōu)勢。但GPS測量是在WGS－84坐標(biāo)系中進(jìn)行的，所提供的高程為相對于WGS－84橢球的大地高H，而我國在工程實際應(yīng)用中使用的是正常高Hr。為此，必須將GPS大地高H轉(zhuǎn)換成正常高Hr，二者之間的差值稱為高程異常 ζi［1～2］。

目前實現(xiàn)GPS大地高與正常高之間的擬合方法有很多，包括:加權(quán)均值法、多項式曲線擬合、多項式曲面擬合、多面函數(shù)曲面擬合、線性移動擬合法、神經(jīng)網(wǎng)絡(luò)法、支持向量基回歸(SVM)模型等［3～7］。其中以多項式曲面擬合中的平面擬合法最為簡單，容易編程實現(xiàn)，但是單一的平面擬合法卻難以滿足較大范圍和地形變化比較大的地區(qū)，文獻(xiàn)［8］指出高程異常由中長波段的趨勢項(地球重力場和重力異常)和短波項(地形起伏引起的高程異常)組成，而在 20 km～200 km內(nèi)中長波段的趨勢項變化趨勢和大地水準(zhǔn)面變化趨勢基本相同且變化趨勢平緩，即只與坐標(biāo)有關(guān)，可用一般擬合函數(shù)進(jìn)行擬合。而短波項與地形信息相關(guān)，變化比較復(fù)雜，大區(qū)域內(nèi)不能擬合?；诖朔N原因，提出構(gòu)建Delaunay三角網(wǎng)和平面擬合法的組合模型，將聯(lián)測點構(gòu)建成Delaunay三角網(wǎng)，使整個區(qū)域分成若干三角形小區(qū)域，然后插入待擬合點進(jìn)行Delaunay重構(gòu)，選擇與待擬合點相關(guān)聯(lián)的聯(lián)測點，并且使用逐點剔除法優(yōu)化這些關(guān)聯(lián)點，選擇最佳擬合組合進(jìn)行平面擬合函數(shù)的求取［9］，進(jìn)而得出待擬合點的高程異常。本文利用MATLAB編程建立GPS高程擬合系統(tǒng)，通過實例驗證得出:該方法可以提高擬合精度，并且能為GPS水準(zhǔn)聯(lián)測點的分布密度和位置選取提供建議。

2 算法模型和程序?qū)崿F(xiàn)流程

2.1 算法模型

(1)平面擬合模型

多項式曲面擬合的一般公式為［10］:

其中 ζi(xi，yi)為點(xi，yi)的高程異常，ai待求參數(shù)，當(dāng)參數(shù)只有3個時，即為平面擬合。

平面擬合法數(shù)學(xué)表達(dá)式為:

式中xi和yi為聯(lián)測點的平面坐標(biāo)，ζi為相應(yīng)的高程異常，a0、a1、a2為待定系數(shù)，因此為求出待定系數(shù)的值必須要有3個或3個以上的聯(lián)測點。

(2)Delaunay三角網(wǎng)算法

生成Delaunay三角網(wǎng)的算法主要分為三類:分治算法;逐點插入法，三角網(wǎng)生長法［10］。本文采用分治算法，該算法首先由 Shamos和 Hoey提出［11］，Lewis和Robinson將該算法首先應(yīng)用于Delaunay三角網(wǎng)。

分治算法的具體步驟為:首先把所有的離散點進(jìn)行以X坐標(biāo)為主Y坐標(biāo)為輔升序排序，將點集平均分成數(shù)據(jù)點個數(shù)相等的兩個子集，在每個子集里再次分成數(shù)量相等的子集，直至每個局部區(qū)域里的數(shù)據(jù)點數(shù)滿足所給的閥值。在每個子區(qū)域里構(gòu)建凸包生成Delaunay三角網(wǎng)，尋找相鄰?fù)拱g兩凸包的頂支撐線和底支撐線，對兩線與凸包邊界圍成的多邊形區(qū)域三角網(wǎng)化，并用局部優(yōu)化準(zhǔn)則對其進(jìn)行優(yōu)化。

(3)區(qū)域GPS水準(zhǔn)優(yōu)化模型

若通過重構(gòu)Delaunay三角形得到的關(guān)聯(lián)點為n個，本文使用平面擬合模型，即選取3個或3個以上的點擬合即可。使用逐點剔除法從n個關(guān)聯(lián)點中逐點剔除對擬合精度貢獻(xiàn)最小的點，最后剩下的3個關(guān)聯(lián)點就是GPS水準(zhǔn)擬合計算的最優(yōu)擬合方案。

從n個點中剔除一點，共有n種剔除方案，選擇其中使得擬合精度最高的一個方案作為n－1個點時的最佳擬合方案，這樣就剔除了一個點。以此類推，直到剩下3個點，即為最優(yōu)的GPS水準(zhǔn)點。這樣該算法總共需要計算n+(n+1)+…+(3+1)中方案。然后比較3～n點擬合方案中擬合精度最高的作為最終擬合方案，所選中的點為最佳擬合點，相應(yīng)分布為最佳分布［12］。

2.2 程序?qū)崿F(xiàn)流程

圖1 程序?qū)崿F(xiàn)流程圖

3 實例計算與分析

某GPS控制網(wǎng)有37個GPS控制點(圖2)，聯(lián)測四等水準(zhǔn)［13］。將 2、20、31、53、87、96、116、143、157、187、190、225、242、243、255、265、291 號點作為初始構(gòu)網(wǎng)擬合點，其余點為檢核點。當(dāng)插入檢核點(以50號點為例)后構(gòu)成新的三角網(wǎng)。從圖3可看到與檢核點構(gòu)成三角網(wǎng)的擬合點，對這些擬合點進(jìn)行優(yōu)化，選出最優(yōu)的擬合組合進(jìn)行擬合參數(shù)的求取。

圖2 擬合點分布圖及初始Delaunay三角網(wǎng)

圖3 重構(gòu)后的Delaunay三角網(wǎng)

當(dāng)插入高程異常構(gòu)成三維Delaunay三角網(wǎng)后(圖4)，可以清楚地看到高程異常面變化非常不規(guī)則，常規(guī)的曲面函數(shù)擬合勢必造成較大的誤差，因此使用Delaunay三角剖分，將大區(qū)域按照地形變化分成多個三角形小區(qū)域可以更為精確地描述出高程異常的變化趨勢，提高擬合精度。因此，地形變化較大地區(qū)，為提高擬合精度必須增加聯(lián)測點的數(shù)量和制定合理的聯(lián)合點分布位置。

圖4 三維Delaunay三角網(wǎng)

為進(jìn)一步提高擬合精度，使用逐點剔除法選擇最佳擬合方案，使擬合結(jié)果最優(yōu)。計算得出擬合中誤差為±1.9 cm，擬合誤差最大的兩個點分別為246、93，擬合誤差分別為 4.8 cm、4.6 cm。從圖5和原數(shù)據(jù)分析具體原因為:246、93兩點所處位置地形變化較大(相關(guān)聯(lián)擬合點高程異常最大差距分別為21 cm、27 cm)，且所測聯(lián)測點比較稀疏(只有3個關(guān)聯(lián)點)，無法體現(xiàn)出地形變化，具體數(shù)據(jù)如表1所示。

幾種擬合方法所得結(jié)果如表2，具體殘差比較如圖5所示。從計算結(jié)果可以看到，本文方法能提高擬合精度。

通過表1和圖5可以看出，各種擬合算法的誤差走向相似，最大擬合誤差均出現(xiàn)在93、246兩點，且本文方法在各點的擬合精度均高于其他方法。

誤差最大點及相關(guān)聯(lián)擬合點數(shù)據(jù)表 表1

四種方法計算結(jié)果比較表 表2

圖5 殘差比較圖

4 結(jié)論

(1)利用構(gòu)建Delaunay三角網(wǎng)和平面擬合組合法進(jìn)行GPS水準(zhǔn)擬合時，是將大區(qū)域以地形變化劃分為多個三角形小區(qū)域，并以待擬合點內(nèi)插重構(gòu)三角形為基準(zhǔn)進(jìn)行擬合，避免了大區(qū)域單一擬合函數(shù)造成的誤差。

(2)利用構(gòu)建Delaunay三角網(wǎng)和平面擬合組合法進(jìn)行GPS高程轉(zhuǎn)換時，其精度與測區(qū)的似大地水準(zhǔn)面的復(fù)雜程度、水準(zhǔn)點的密度和分布有關(guān)，足夠的密度和合理的分布，是保證轉(zhuǎn)換精度的必要條件［14］。

(3)使用GPS水準(zhǔn)優(yōu)化法選擇關(guān)聯(lián)擬合點的最優(yōu)組合能剔除可能造成擬合值發(fā)生突變的擬合點，提高擬合精度。

(4)通過實例計算表明，該方法容易編程實現(xiàn)，可進(jìn)行大批量GPS水準(zhǔn)擬合，并且通過對比構(gòu)建Delaunay三角網(wǎng)和平面擬合組合法有較高精度。

［1］ 李征航，黃勁松．GPS測量與數(shù)據(jù)處理［M］．武漢:武漢大學(xué)出版社，2005．

［2］寧津生，劉經(jīng)南，陳俊勇等．現(xiàn)代大地測量理論與技術(shù)［M］．武漢:武漢大學(xué)出版社，2006．

［3］張小紅，程世來，許曉東．基于Kriging統(tǒng)計的GPS高程擬合方法研究［J］．大地測量與地球動力學(xué)，2007(2):47～51．

［4］楊明清，朱達(dá)成，陳現(xiàn)春．用神經(jīng)網(wǎng)絡(luò)方法轉(zhuǎn)換GPS高程［J］．測繪學(xué)報，1999(4):301～307．

［5］王繼剛，胡永輝，孔令杰．基于最小二乘支持向量機(jī)的區(qū)域GPS高程轉(zhuǎn)換組合［J］．大地測量與地球動力學(xué)，2009(5):99～102．

［6］高偉，徐紹銓．GPS高程分區(qū)擬合轉(zhuǎn)換正常高的研究［J］．武漢大學(xué)學(xué)報(信息科學(xué)版)，2004(10):908～911．

［7］ 叢康林，岳建平．基于SVR的GPS高程擬合模型研究［J］．測繪通報，2011(2):8～11．

［8］王旭，劉文生．GPS高程擬合方法的研究［J］．測繪科學(xué)，2010(35):28～30．

［9］沈云中，高達(dá)凱．GPS水準(zhǔn)點優(yōu)化選擇法［M］．兩岸重力及水準(zhǔn)面研討會，臺北，2003:131～135．

［10］劉萬選，劉加生．兩種常用的GPS高程擬合模型擬合精度研究［J］．測繪與空間地理信息，2009，32(3):147～150，156．

［11］Shamos M I．a(chǎn)nd Hoey D．Closest－point Problems，In:Proceedings of the 16th Annual Symposium on the Foundations of Computer Science［C］．1975，151 ～162．

［12］趙超英，劉雷．GPS水準(zhǔn)擬合優(yōu)化法探討［J］．工程勘測，2006(2):64～67．

［13］金雪漢，尹長林．GPS高程轉(zhuǎn)換中的函數(shù)模型逼近研究［J］．長沙電力學(xué)院學(xué)報(自然科學(xué)版)，2005(2):75～77．

［14］陳鵬，陳正陽．廣義延拓法在 GPS高程轉(zhuǎn)換中的應(yīng)用［J］．大地測量與地球動力學(xué)，2010，30(2):125～128．