基于構(gòu)建Delaunay三角網(wǎng)和平面擬合法組合的GPS高程轉(zhuǎn)換

王新 ，張紹良

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

1 引言

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

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

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

2.1 算法模型

(1)平面擬合模型

多項(xiàng)式曲面擬合的一般公式為［10］:

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

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

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

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

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

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

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

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

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

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

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

3 實(shí)例計(jì)算與分析

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

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

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

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

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

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

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

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

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

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

圖5 殘差比較圖

4 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

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

［10］劉萬選，劉加生．兩種常用的GPS高程擬合模型擬合精度研究［J］．測(cè)繪與空間地理信息，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］．工程勘測(cè)，2006(2):64～67．

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

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