基于雙調(diào)和樣條內(nèi)插和高斯曲率極值的多面函數(shù)擬合高程異常方法

孫佳龍，焦明連，郭淑艷

(1. 淮海工學(xué)院 測繪工程學(xué)院，江蘇 連云港 222001；2. 江蘇省海洋資源開發(fā)研究院，江蘇 連云港 222005)

一、引 言

隨著GPS測量技術(shù)的不斷完善和發(fā)展，求取高精度的高程異常在實(shí)際工程應(yīng)用中越來越重要[1-3]。在擬合區(qū)域較小的情況下，一般的擬合方法可以達(dá)到較高的精度，但由于這些方法未顧及似大地水準(zhǔn)面的物理性質(zhì)，擬合出的曲面只能是高程異常的趨勢面，與高程異常的實(shí)際數(shù)值存在一定差異[4-7]。而多面函數(shù)在理論上可以以任意精度逼近任意復(fù)雜曲面，因此可以更加逼近實(shí)際的高程異常曲面[8]。而當(dāng)GPS/水準(zhǔn)點(diǎn)數(shù)量較多時(shí)，如果隨機(jī)選取擬合點(diǎn)，則選擇的方案數(shù)過多，可能超出普通計(jì)算機(jī)的計(jì)算能力。鑒于高程異常曲面具有均衡性，且在數(shù)學(xué)上可用調(diào)和函數(shù)來描述，本文提出了基于雙調(diào)和樣條及高斯曲率極值的多面函數(shù)擬合方法，通過計(jì)算高程異常曲面網(wǎng)格數(shù)據(jù)的高斯曲率及極值，選取距離高斯曲率極值點(diǎn)最近的GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)作為多面函數(shù)的結(jié)點(diǎn)，參與高程異常曲面的擬合，而剩余點(diǎn)作為檢核點(diǎn)，對擬合的曲面進(jìn)行外部檢核。

二、雙調(diào)和樣條函數(shù)

對于一組不規(guī)則的空間數(shù)據(jù)點(diǎn)，通?？梢岳萌螛訔l函數(shù)找出最平滑的曲線或曲面。從力學(xué)的角度看，該過程就是強(qiáng)迫一個(gè)彈性桿或彈性板去匹配各個(gè)數(shù)據(jù)點(diǎn)[9]。

樣條的格林函數(shù)的點(diǎn)壓力滿足雙調(diào)和方程[10]

(1)

式中，x為空間位置；δ(x)是Delta函數(shù)；格林函數(shù)Φ(x)就是方程式(1)的一個(gè)解，表示內(nèi)插曲線或曲面的位移，是空間位置x的函數(shù)，可以寫成

(2)

當(dāng)該格林函數(shù)對N個(gè)在xi的數(shù)據(jù)點(diǎn)ωi進(jìn)行插值時(shí)，式(1)就成為

(3)

ω(xi)=ωi

(4)

式中，ωi是內(nèi)插曲線或曲面在空間位置xi處的位移；αj為加權(quán)系數(shù)。

三、基于高斯曲率極值的結(jié)點(diǎn)判定方法

在曲面論中，Weingarten公式刻畫了曲面單位法向的運(yùn)動(dòng)狀況，可以進(jìn)一步了解曲面的彎曲狀況。Weingarten公式的矩陣形式為[11]

(5)

式中，矩陣ω為曲面在參數(shù)(u1,u2)下的Weingarten矩陣。ω可以表示為

曲面上一點(diǎn)的法曲率關(guān)于切線方向的兩個(gè)最值分別為曲面在該點(diǎn)處的主曲率，而使法曲率達(dá)到最值的兩個(gè)切線方向分別為曲面在該點(diǎn)處的主方向。曲面的主曲率和主方向分別為Weingarten變換的特征值和特征方向。

曲面在一點(diǎn)處的兩個(gè)主曲率的乘積K稱為在該點(diǎn)處的高斯曲率。高斯曲率用Weingarten矩陣表示為[12]

對于曲面上某分片區(qū)域的一點(diǎn)，如果該點(diǎn)的2個(gè)主曲率(K1和K2)中任何一個(gè)沿著對應(yīng)的主方向上為極值，則該點(diǎn)稱為曲率極值點(diǎn)。由于高斯曲率極值可以很好地刻畫曲面的彎曲情況，因此，本文選取高程異常曲面上距離高斯曲率極值點(diǎn)最近的GSP/水準(zhǔn)數(shù)據(jù)點(diǎn)作為最后參與擬合的數(shù)據(jù)點(diǎn)，即多面函數(shù)中核函數(shù)的結(jié)點(diǎn)，參與多面函數(shù)的擬合。

四、算例分析

某市E級GPS控制網(wǎng)，進(jìn)行了三等水準(zhǔn)測量，共包含71個(gè)GPS/水準(zhǔn)點(diǎn)，本文以此數(shù)據(jù)為例，利用雙調(diào)和樣條函數(shù)插值出了該市的高程異常曲面，并計(jì)算了高程異常曲面網(wǎng)格數(shù)據(jù)的高斯曲率及極值點(diǎn)。采用距離高斯曲率極值點(diǎn)最近的原則對所有GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)進(jìn)行篩選，利用剩余數(shù)據(jù)點(diǎn)對擬合的曲面進(jìn)行檢核。將篩選出的數(shù)據(jù)點(diǎn)作為結(jié)點(diǎn)進(jìn)行多面函數(shù)擬合，而其余點(diǎn)作為檢核點(diǎn)對擬合的高程異常曲面進(jìn)行外部精度檢核，各種數(shù)據(jù)點(diǎn)如圖1所示。

在圖1中，星號表示高斯曲率極值點(diǎn)，三角符號為所有的GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)，而三角符號加四邊形符號為從所有GPS/水準(zhǔn)數(shù)據(jù)中選擇的參與擬合的數(shù)據(jù)點(diǎn)，即GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)中距離高斯曲率極值點(diǎn)最近的點(diǎn)。圖中顏色條表示高程異常值，單位為m。從圖1中可以看出，整個(gè)觀測區(qū)域的高程異常曲面呈現(xiàn)北高南低的趨勢，但在中部卻又有一個(gè)低凹的異常區(qū)域出現(xiàn)，因此，觀測區(qū)域的高程異常變化趨勢比較復(fù)雜。所有的GPS/水準(zhǔn)數(shù)據(jù)又呈現(xiàn)出南多北少的分布特點(diǎn)，特別是在東北區(qū)域，沒有一個(gè)數(shù)據(jù)點(diǎn)。通過計(jì)算高斯曲率極值，即該點(diǎn)的高斯曲率均大于(極大值)或小于(極小值)周圍4個(gè)網(wǎng)格點(diǎn)的高斯曲率，如圖1中各種顏色的最深處和兩種顏色之間的過渡處，表示曲面彎曲變化較大較快的情況，因此，利用高斯曲率極值可以很好地反映曲面的變化情況。由于計(jì)算出的高斯曲率均位于插值出的網(wǎng)格點(diǎn)上，而GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)一般并不與這些網(wǎng)格點(diǎn)重合，從圖1中可以看出，高斯曲率極值點(diǎn)有時(shí)與GPS/水準(zhǔn)點(diǎn)相距較近，而有時(shí)卻很遠(yuǎn)，但是，距離這些高斯曲率極值網(wǎng)格點(diǎn)最近的GPS/數(shù)據(jù)點(diǎn)都可以在一定程度上反映曲面的彎曲情況，在圖1中也可以看出，被選擇的數(shù)據(jù)點(diǎn)也比較均勻地分布在觀測區(qū)域里，因此，選擇這些數(shù)據(jù)點(diǎn)參與曲面的擬合可以在一定程度上綜合反映高程異常曲面的變化趨勢。

圖1 高斯曲率極值點(diǎn)、GPS/水準(zhǔn)數(shù)據(jù)點(diǎn)及選擇的擬合點(diǎn)

在71個(gè)GPS/水準(zhǔn)數(shù)據(jù)中，共有44個(gè)高斯曲率極值點(diǎn)，利用高斯曲率極值點(diǎn)作為多面函數(shù)的結(jié)點(diǎn)對高程異常曲面進(jìn)行擬合。當(dāng)高斯曲率極值K>0時(shí)，多面函數(shù)采用橢圓拋物面函數(shù)進(jìn)行擬合；而當(dāng)高斯曲率極值K<0時(shí)，采用雙曲拋物面函數(shù)進(jìn)行擬合。以剩余的27個(gè)數(shù)據(jù)點(diǎn)作為檢核點(diǎn)，二次多項(xiàng)式、移動(dòng)二次曲面、基于雙曲拋物面的多面函數(shù)(hyperbolic paraboloid of polyhedral function，HPPF)、橢圓拋物面的多面函數(shù)(elliptic paraboloid of polyhedral function，EPPF)及本文提出的基于高斯曲率極值的多面函數(shù)(Gaussian curvature extremum of polyhedral function，GCEPF)方法的擬合殘差如圖2所示。

從圖2中可以看出，二次多項(xiàng)式擬合的殘差序列波動(dòng)性較大，在多個(gè)檢核點(diǎn)的殘差都較大，最大絕對值約為10 cm；而移動(dòng)二次曲面和EPPF的殘差序列的趨勢雖較為平緩，但最大殘差絕對值也接近6 cm；HPPF和本文提出的GCEPF方法得到的殘差序列的整體趨勢趨于平緩，最大絕對值約為4 cm。為了更細(xì)致地比較幾種方法的優(yōu)劣性，各種擬合方法擬合得到殘差的最大值、最小值、平均值和標(biāo)準(zhǔn)差見表1。

圖2 5種方法的擬合殘差

表1 殘差序列的統(tǒng)計(jì)特征 cm

從表1中可以看出，GCEPF方法殘差最大值雖然比其他幾種方法(EPPF除外)大1～2 cm，但最小值卻比其他幾種方法增加4～9 cm，使極值的誤差控制在5 cm以內(nèi)(其他4種方法均超過5 cm)，滿足了一般測量工作對高程精度的要求。

HPPF和EPPF殘差序列的標(biāo)準(zhǔn)差比二次多項(xiàng)式分別減小了27.2 mm和25.7 mm，但卻比移動(dòng)二次曲面增加了2.5 mm和4 mm，說明采用單一某種核函數(shù)的多面函數(shù)擬合的精度優(yōu)于二次多項(xiàng)式，但卻不及移動(dòng)二次曲面。可見，由于移動(dòng)二次曲面擬合方法顧及了擬合點(diǎn)周邊的高程異常起伏的影響，因此擬合的曲面更加逼近真值，擬合精度較高。

而以高斯曲率極值為參考的GCEPF方法能更細(xì)致地識(shí)別高程異常曲面的突變情況，其殘差序列的標(biāo)準(zhǔn)差比二次多項(xiàng)式減小了33.5 mm，比其他采用單一某種曲面的多面函數(shù)HPPF和EPPF減少了6.3 mm和7.8 mm，比移動(dòng)二次曲面也減小了3.8 mm，因此，GCEPF殘差序列更趨于平緩，更逼近真值，優(yōu)于以上4種方法。

五、結(jié)束語

對于區(qū)域高程異常曲面，基于某種單一核函數(shù)的多面函數(shù)很難較好地進(jìn)行擬合，而從某種意義上說，基于單一核函數(shù)的多面函數(shù)已經(jīng)失去了多面函數(shù)的優(yōu)勢，而在實(shí)際情況中，想要找到適合不同區(qū)域的核函數(shù)，實(shí)現(xiàn)多個(gè)核函數(shù)的組合擬合，從而真正體現(xiàn)其“多面性”是比較困難的。

本文利用雙調(diào)和樣條函數(shù)對高程異常曲面進(jìn)行插值，利用高斯曲率極值的條件選擇核函數(shù)及其結(jié)點(diǎn)，既顧及了高程異常曲面的均衡性，同時(shí)又能有效地對核函數(shù)及其結(jié)點(diǎn)加以選定，從而實(shí)現(xiàn)了利用“多面”函數(shù)進(jìn)行擬合的目的。最終的擬合結(jié)果也說明，本文提出的基于雙調(diào)和樣條函數(shù)及高斯曲率極值的高程異常擬合方法的擬合精度比以上4種方法得到的精度更高，更逼近真值。

參考文獻(xiàn)：

[1] 陳傳法, 岳天祥. 基于預(yù)處理共軛梯度法的高精度曲面模型解算[J]. 中國礦業(yè)學(xué)報(bào), 2010, 39(2): 290-294.

[2] 程蘆穎, 許厚澤. 顧及測站點(diǎn)上重力場信息的大地水準(zhǔn)面高的擬合方法[J]. 武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2008, 33(7): 701-705.

[3] 車延國. GPS高程水準(zhǔn)擬合模型與精度分析[J]. 遼寧工程技術(shù)大學(xué)學(xué)報(bào), 2008, 27(1): 35-38.

[4] 劉大杰，姚連璧，周全基. GPS水準(zhǔn)的擬合基準(zhǔn)面高程系統(tǒng)[J]. 測繪學(xué)報(bào), 2000, 29(S1):11-15.

[5] 胡伍生, 華錫生, 張志偉. 平坦地區(qū)轉(zhuǎn)換GPS高程的混合轉(zhuǎn)換方法[J]. 測繪學(xué)報(bào), 2002,31(2): 128-133.

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

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

[8] 張菊清, 劉平芝. 抗差趨勢面與正交多面函數(shù)結(jié)合擬合DEM數(shù)據(jù)[J]. 測繪學(xué)報(bào), 2008, 37(4): 526-530.

[9] SANDWELL D. Bihamonic Spline Interpolation of GEOS-3 and SEASAT Altimeter Data[J]. Geophysical Research Letter, 1987, 14(2): 139-142.

[10] JUN Yongtae, CHOI Kuiwon. Automated Feature-based Registration for Reverse Engineering of Human Model[J]. Journal of Mechanical Science and Technology, 2005, 19(12): 2213-2223.

[11] 馬驪溟, 徐毅, 李澤湘. 基于高斯曲率極值點(diǎn)的散亂點(diǎn)云數(shù)據(jù)特征點(diǎn)提取[J]. 系統(tǒng)仿真學(xué)報(bào), 2008, 20(9): 2341-2344.

[12] RABBARI T, DIJKMAN S. An Intergrated Approach for Modeling and Global Registration of Point Clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2007,61(2): 355-370.

[13] YANG M, LEE E. Segmentation of Measured Data Using a Parametric Quadric Surface Approximation [J]. Computer Aided Design, 1999, 31(7): 449-457.