王延存，俞家勇，田茂義，周茂倫，李國玉

(1.山東科技大學(xué) 測繪科學(xué)與工程學(xué)院，山東 青島 266590；2.青島秀山移動測量有限公司，山東 青島 266590)

0 引言

車載移動測量系統(tǒng)[1]由三維激光掃描儀、慣性測量單元(inertial measurement unit，IMU)、全球定位系統(tǒng)(global positioning system，GPS)等傳感器組成，該系統(tǒng)能夠快速獲取高精度點云數(shù)據(jù)，可廣泛應(yīng)用于三維建模、大比例尺地形圖、城市部件普查等方面,可以有效降低測量人員的勞動強度，提高工作效率。車載移動測量系統(tǒng)憑借自身的優(yōu)勢，成為當(dāng)前研究的熱點。然而點云數(shù)據(jù)的高程為大地高H，常規(guī)工程測量的高程是正常高h(yuǎn)，為了能將車載激光點云數(shù)據(jù)運用于一般工程測量中，必須把大地高H轉(zhuǎn)化為正常高h(yuǎn)；大地高與正常高之間的差值為高程異常ζ，即ζ=H-h。高程異常如圖1所示。因此如何獲取高精度的高程異常改正值，成為點云數(shù)據(jù)處理的一個研究重點。目前確定高程異常的方法有地球重力場模型法[2-4]和GPS水準(zhǔn)法。重力場模型法是通過EGM2008地球重力場模型計算待定點的高程異常值，該方法可以獲取高精度高程異常值，但是重力數(shù)據(jù)獲取比較困難且成本較高，因此在工程中應(yīng)用比較少。GPS水準(zhǔn)法是在測區(qū)內(nèi)布設(shè)一定數(shù)量且分布均勻的控制點，進行GPS和水準(zhǔn)聯(lián)測，得到控制點的高程異常，進而通過曲面擬合獲得該區(qū)域的高程異常數(shù)學(xué)表達(dá)式，再內(nèi)插得到該區(qū)域其他點的高程異常值。常用的曲面擬合方法包括多項式曲面擬合法、多面函數(shù)法、移動曲面法、最小二乘配置法、薄板樣條函數(shù)法(TPS)以及BP神經(jīng)網(wǎng)絡(luò)法等[5-15]。地形變化較大區(qū)域，文獻(xiàn)[16-17]通過結(jié)合不同的擬合方法來適應(yīng)不同的地區(qū)；針對小區(qū)域范圍，文獻(xiàn)[18-19]進行了研究，利用云模型對高程異常進行定性定量分析，形成云形圖，直觀獲得該小區(qū)域的高程異常值。以上研究都取得了較好的效果。

圖1 高程異常示意圖

在實際的工程應(yīng)用中，由于外界條件的影響，組合導(dǎo)航定位精度較差區(qū)域會導(dǎo)致點云高程誤差大，此時利用文獻(xiàn)[5-19]所提擬合方法求得的模型參數(shù)的可靠性降低，不能滿足高精度測圖要求。對于高程異常值含有粗差的情況，文獻(xiàn)[20-24]進行了研究，但是局限于通過降低粗差的權(quán)值來提高模型參數(shù)的穩(wěn)定性，而沒有剔除粗差。且文獻(xiàn)[20-24]都是針對大區(qū)域含粗差情況下的高程異常擬合，小區(qū)域范圍內(nèi)的穩(wěn)健高程異常擬合相關(guān)研究較少。針對這一情況，本文提出了一種穩(wěn)健高程異常改正方法，該方法以薄板樣條函數(shù)為基礎(chǔ)，結(jié)合RANSAC算法，并以狄克松判別法為準(zhǔn)則，獲取高程異常值。通過實驗證明該方法可以有效的剔除粗差，獲得高精度的高程異常改正值。

1 薄板樣條函數(shù)原理

薄板樣條函數(shù)插值[25-26]是一種二維的插值方法，在具有3個以上非共線點的情況下，薄板樣條函數(shù)具有唯一解。假設(shè)在空間區(qū)域R2內(nèi)有n個非共線點(t1，t2，…，tn)，每個點的坐標(biāo)為(xi，yi，zi)(i=1，2，…，n)，其中zi=z(xi，yi)。

(1)

式中：δi和bj為待定系數(shù)；‖·‖為歐氏范數(shù)，在二維空間即為兩點間的距離。ω是TPS的核函數(shù)，ω(r)=r2ln r/16π。當(dāng)滿足條件ωT=0時，式(1)就是空間區(qū)域R2上的自然薄板樣條函數(shù)。

2 穩(wěn)健自然薄板樣條函數(shù)法

在小區(qū)域范圍內(nèi)，高程異?？梢钥闯晒饣B續(xù)曲面，非常適合通過自然薄板樣條函數(shù)進行擬合。但當(dāng)擬合數(shù)據(jù)中偶然誤差較大或含有粗差時，直接利用自然薄板樣條函數(shù)擬合得到的參數(shù)δi和bj不是最優(yōu)解，進行插值時，距離粗差點越近的區(qū)域受到的影響越大。在這種情況下，需要引進一種穩(wěn)健的擬合方法。本文在自然薄板樣條函數(shù)的基礎(chǔ)上，結(jié)合RANSAC算法，通過狄克松判別法尋找并剔除粗差點，進而得到該區(qū)域的高程異常值。

2.1 RANSAC算法

RANSAC算法是最有效的穩(wěn)健估計算法之一，當(dāng)數(shù)據(jù)中的粗差數(shù)據(jù)超過50%時，使用該算法仍然可以得到理想的估計參數(shù)。RANSAC算法首先設(shè)定一個準(zhǔn)則，根據(jù)這個準(zhǔn)則將數(shù)據(jù)分為局內(nèi)點(可以用數(shù)學(xué)模型表達(dá)的數(shù)據(jù))和局外點(不可以用數(shù)學(xué)模型表達(dá)的數(shù)據(jù))，最后利用保留下來的局內(nèi)點進行參數(shù)計算，從而得到最優(yōu)參數(shù)估值。最小迭代數(shù)k通過公式(2)計算得到。

(2)

式中：m為計算模型需要的最小數(shù)據(jù)量；ε為數(shù)據(jù)錯誤率。通常情況下，ε、P和k是根據(jù)實際情況動態(tài)確定的。

2.2 狄克松判別法

狄克松判別法采用極差比的思想，可以快速有效檢測粗差的存在。RANSAC算法中，判定條件的選取直接決定最終結(jié)果的精度，因此在本文中，利用狄克松判別法作為局內(nèi)點和局外點的判定條件。

通過RANSAC算法隨機選取待擬合點，利用自然薄板樣條函數(shù)進行擬合，得到每個待擬合點的高程異常殘差Δξi，將殘差Δξi按照從小到大的順序排列。

Δξ1≤Δξ2≤…≤Δξn

(3)

(4)

2.3 穩(wěn)健擬合方法

利用自然薄板樣條函數(shù)進行曲面擬合至少需要4個非共線點，根據(jù)RANCAC算法的思想，即需要在待擬合點中隨機選取4個點作為種子點，根據(jù)式(2)計算最小迭代次數(shù)k，并計算自然薄板樣條函數(shù)模型參數(shù)δi和bj的初始值，利用該初始值計算其余待擬合點的高程異常值，并與已知高程異常值做差得到高程異常值差值Δξi，利用狄克松判別法判定尋找該組擬合點對應(yīng)的粗差值。將粗差值判定為局外點，其余的作為局內(nèi)點，并統(tǒng)計局內(nèi)點的數(shù)量。迭代統(tǒng)計最大局內(nèi)點數(shù)量，將該組數(shù)據(jù)作為最終的待擬合點。該穩(wěn)健擬合方法的流程如圖2所示。

①利用ε、P、m根據(jù)式(2)計算出最小迭代次數(shù)k，本文中m=4。

②隨機選取m個點，計算出模型參數(shù)δi和bj的初始值，并利用自然薄板樣條函數(shù)的特性檢核參數(shù)是否正確。

③根據(jù)②得到的參數(shù)計算Δξi，根據(jù)公式(3)、公式(4)，進行自殘序列兩端逐點判別，如果小于臨界值，將其判定為局內(nèi)點，反之為局外點。

④重復(fù)②至③k次，統(tǒng)計每次迭代的局內(nèi)點數(shù)量及局內(nèi)點點號。

⑤選取最大局內(nèi)點，根據(jù)公式(1)計算出參數(shù)δi和bj的最優(yōu)解。

⑥通過最優(yōu)解計算檢核點殘差。

圖2 穩(wěn)健擬合方法流程圖

3 實驗及結(jié)果分析

本文采用的數(shù)據(jù)為南方某地的車載移動測量系統(tǒng)獲取的點云數(shù)據(jù)，該地區(qū)為丘陵地區(qū)。高程變化比較平緩。在測區(qū)內(nèi)通過隨機抽取，共選擇了26個均勻分布的地面點，根據(jù)某勘測院1∶500水準(zhǔn)精化模型獲得對應(yīng)的正常高，在該區(qū)域選擇7個高程不同的地面點作為檢核點，檢核點高程分布在32～43 m范圍內(nèi)。高程異常擬合點點位分布如圖3所示。

圖3 待擬合點點位分布

針對移動測量作業(yè)過程中，組合導(dǎo)航系統(tǒng)定位精度受環(huán)境影響較大，采集前需要對作業(yè)區(qū)域進行規(guī)劃，確保系統(tǒng)精度能夠滿足地籍測量要求。但在實際作業(yè)過程中，依然會存在局部小范圍區(qū)域GPS定位精度較差情況，一般此區(qū)域高程誤差在5～10 cm。一般誤差較大區(qū)域占總的采集范圍大約在5%～10%。在選擇地面觀測數(shù)據(jù)過程中，為保證擬合點均勻分布，采樣擬合點難免會有少數(shù)部分在誤差較大區(qū)域，需要在擬合方法中考慮粗差剔除，因此粗差取值范圍在5～10 cm。本文考慮粗差點在擬合點數(shù)據(jù)中所占比例較小(考慮極端條件10%，擬合點數(shù)為26)，故分別取1個、3個粗差點進行分析，闡述本文方法的穩(wěn)健性。

目前最常用的高程異常擬合方法是多項式曲面法，本文采用三次曲面法進行擬合作為對比實驗。三次曲面法是把高程異常值構(gòu)成的曲面看作標(biāo)準(zhǔn)的三次曲面，而在地形變化較大的地區(qū)，高程異常值構(gòu)成的曲面比較復(fù)雜，三次曲面法不能真實地描述該曲面，所得參數(shù)誤差較大。因此該方法適用于高程異常變化比較平緩，且變化趨勢近似于三次曲面的平原、低矮丘陵地區(qū)。使用該方法時，擬合點位要分布均勻且不含有粗差，否則局部地區(qū)可能會有較大誤差。利用三次曲面法進行高程異常擬合，只需求解10個待定系數(shù)，計算效率較高，且方程中最高次項為三階，有效避免了“龍格”現(xiàn)象。本文實驗數(shù)據(jù)所在區(qū)域為丘陵地區(qū)，高程異常值變化比較平緩，適合使用三次曲面法進行擬合。

在實驗中，分別利用三次曲面法、自然薄板樣條函數(shù)法、穩(wěn)健自然薄板樣條函數(shù)法進行高程異常擬合，統(tǒng)計3種方法的檢核點高程異常值的殘差Δξi，結(jié)果如圖4所示。

圖4 3種方法各檢核點殘差比較

根據(jù)圖4有以下分析結(jié)果：

①從圖4(a)可以看出，在數(shù)據(jù)不包含粗差的情況下，其中三次曲面擬合的殘差最大值為3.4 cm，而通過薄板樣條函數(shù)進行擬合時，殘差在2.5 cm以下，說明在該區(qū)域薄板樣條函數(shù)法的擬合精度優(yōu)于三次曲面法。實驗表明，2種方法在擬合點不含有粗差的情況下，其精度都可以滿足1∶500地形圖要求。

②從圖4(b)、圖4(c)可以看出，在數(shù)據(jù)中通過數(shù)值模擬加入一個粗差點時，用薄板樣條函數(shù)法擬合得到參數(shù)，求得的第3檢核點的殘差較大。加入3個粗差點時，薄板樣條函數(shù)法對應(yīng)的第3、4、6檢核點的殘差較大，最大值達(dá)到了5.2 cm，且只有一個檢核點超過了1：500地形圖的高程精度要求；利用三次曲面進行擬合得到的檢核點的殘差較大，最大值達(dá)到了6.8 cm，且在粗差較多的情況下，50%以上的檢核點的高程精度都不能滿足1∶500地形圖的要求。因此，三次曲面法在該區(qū)域進行高程異常擬合的抗差性要低于薄板樣條函數(shù)法。

③利用本文提出的方法，雖然隨著待擬合點中粗差的增多，各檢核點的殘差增大，但是各檢核點的殘差都控制在3 cm以下，即在模擬的極端條件下(含有10%粗差)，本文提出的方法在該區(qū)域依然滿足1∶500地形圖的精度要求。導(dǎo)致這一現(xiàn)象的原因是，傳統(tǒng)的自然薄板樣條函數(shù)法是將所有的數(shù)據(jù)不加區(qū)分地全部用于擬合中，這導(dǎo)致求解得到的模型參數(shù)不能反映真實的高程異常值，根據(jù)自然薄板樣條函數(shù)的特性可知，距離粗差點越近的檢核點，其殘差值越大；而穩(wěn)健薄板樣條函數(shù)法對應(yīng)的殘差較小且穩(wěn)定，是因為穩(wěn)健法并使用RANSAC思想，并根據(jù)狄克松準(zhǔn)則將粗差點剔除，尋找到最大局內(nèi)點，最后利用保留下來的待擬合點(局內(nèi)點)進行擬合，體現(xiàn)出該方法具有較強的穩(wěn)健性。

4 結(jié)束語

在點云數(shù)據(jù)含有粗差的情況下，利用傳統(tǒng)的擬合方法得到的模型參數(shù)是不可靠的。針對這一情況，本文在自然薄板樣條函數(shù)的基礎(chǔ)上，提出了一種穩(wěn)健高程異常改正方法。該方法的核心思想是結(jié)合RANSAC算法，并根據(jù)狄克松判別法剔除粗差點，從而達(dá)到穩(wěn)健的效果。實驗結(jié)果表明，本文提出的方法可以在保證效率的前提下，有效地消除粗差點帶來的影響，提高高程異常改正精度，解決了傳統(tǒng)方法的不足。同時排除的粗差點還可以作為測圖誤差較大區(qū)域的標(biāo)志，有助于車載激光點云數(shù)據(jù)質(zhì)量把控，提高整體作圖精度；該方法不僅可以應(yīng)用于點云數(shù)據(jù)的高程異常擬合，還可以應(yīng)用于其他的包含粗差的曲面擬合和制作DEM。需要指出的是，由于自然薄板樣條函數(shù)擬合得到的曲面具有連續(xù)、光滑的特性，因此本方法在分區(qū)域擬合中也具有一定的優(yōu)勢。