基于構(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，二者之間的差值稱(chēng)為高程異常 ζi［1～2］。

目前實(shí)現(xiàn)GPS大地高與正常高之間的擬合方法有很多，包括:加權(quán)均值法、多項(xiàng)式曲線(xiàn)擬合、多項(xiàng)式曲面擬合、多面函數(shù)曲面擬合、線(xiàn)性移動(dòng)擬合法、神經(jīng)網(wǎng)絡(luò)法、支持向量基回歸(SVM)模型等［3～7］。其中以多項(xiàng)式曲面擬合中的平面擬合法最為簡(jiǎn)單，容易編程實(shí)現(xiàn)，但是單一的平面擬合法卻難以滿(mǎ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)不能擬合?；诖朔N原因，提出構(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ù)的求?。?］，進(jìn)而得出待擬合點(diǎn)的高程異常。本文利用MATLAB編程建立GPS高程擬合系統(tǒng)，通過(guò)實(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)的算法主要分為三類(lèi):分治算法;逐點(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ù)滿(mǎn)足所給的閥值。在每個(gè)子區(qū)域里構(gòu)建凸包生成Delaunay三角網(wǎng)，尋找相鄰?fù)拱g兩凸包的頂支撐線(xiàn)和底支撐線(xiàn)，對(duì)兩線(xiàn)與凸包邊界圍成的多邊形區(qū)域三角網(wǎng)化，并用局部?jī)?yōu)化準(zhǔn)則對(duì)其進(jìn)行優(yōu)化。

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

若通過(guò)重構(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)。以此類(lèi)推，直到剩下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))，無(wú)法體現(xiàn)出地形變化，具體數(shù)據(jù)如表1所示。

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

通過(guò)表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)通過(guò)實(shí)例計(jì)算表明，該方法容易編程實(shí)現(xiàn)，可進(jìn)行大批量GPS水準(zhǔn)擬合，并且通過(guò)對(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］張小紅，程世來(lái)，許曉東．基于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］高偉，徐紹銓?zhuān)瓽PS高程分區(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］劉萬(wàn)選，劉加生．兩種常用的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］陳鵬，陳正陽(yáng)．廣義延拓法在 GPS高程轉(zhuǎn)換中的應(yīng)用［J］．大地測(cè)量與地球動(dòng)力學(xué)，2010，30(2):125～128．