基于K-means聚類與RBFNN的點云DEM構建方法

趙慶展 李沛婷 馬永建 田文忠

(1.石河子大學信息科學與技術學院， 石河子 832003； 2.國家遙感中心新疆兵團分部， 石河子 832003；3.兵團空間信息工程技術研究中心， 石河子 832003； 4.石河子大學機械電氣工程學院， 石河子 832003)

0 引言

數字高程模型(Digital elevation model，DEM)是地形高程信息的數字化表示，是重要的基本地形產品之一[1]，不僅可以非常直觀地顯示地形和地貌，而且為各種地形特征的定量分析和不同類型專題圖的繪制提供了基本數據[2]。無人機機載激光雷達(Light detection and ranging，LiDAR)是集激光測距技術、計算機技術、高精度動態慣性導航系統(Inertial navigation system，INS)和高精度動態差分全球定位系統(Differential global positioning system，DGPS)等于一體的攝影測量新技術[3]，可提供地物三維空間坐標、地物激光回波強度、獲取時間等信息[4]。因無人機機載LiDAR具有起降靈活、機動性好等優點而被廣泛應用，可直接獲取高精度、高密度點云數據以生成DEM。

同時，點云濾波是生成DEM的基礎，即把點云分成地面點與非地面點。林祥國等[5]采用多基元三角網漸進加密方法實現點云濾波，使濾波誤差降低。李鵬程等[6]基于表面實現點云濾波，即采用點云的波形信息進行加權曲面擬合以獲取點云表面。張繼賢等[7]利用K-means聚類方法快速分割點云中的電線。周曉明[8]結合點云全波形屬性信息和聚類方法實現濾波。綜上，具有代表性的點云濾波方法可以分為基于逐漸加密、基于表面和基于聚類等方法[9]。本文選擇原理簡單、且便于處理海量數據的K-means聚類方法實現點云濾波。

神經網絡不需建立精確的數學模型即可實現非線性映射，以預測相關數據。KU?AK等[10]利用自組織映射神經網絡和K-means聚類兩種方法，結合點云法向量、強度和曲率3種屬性信息，實現點云分割，以生成DEM。ZUO等[11]采用神經網絡實現城市點云分類，以驗證神經網絡可以成功分離地面點、建筑物、樹木和裸露土等地物。徑向基函數神經網絡(Radical basis function neural network，RBFNN)以其學習速度快、不易陷入局部極小值等優點，常被應用在數據插值預測中。陳昌華等[12]利用主成分分析方法消除徑向基函數神經網絡(RBFNN)輸入層數據的相關性，即以主成分分析模型的結果作為RBFNN的輸入，建立土壤水分的預測模型，結果表明，RBFNN預測精度比誤差反向傳播神經網絡預測精度高。周仲禮等[13]采用模擬退火蟻群算法改進徑向基函數神經網絡模型，分別對地層高程進行面插值和對礦體品位進行空間體插值，且與普通克里金(Kriging)插值方法進行交叉驗證，結果表明，改進徑向基函數神經網絡插值效果明顯優于克里金插值方法。目前使用徑向基函數神經網絡實現點云預測的研究較少[14]，本文使用徑向基函數神經網絡預測地面點云高程值，且采用Delaunay三角網生成DEM。

1 數據獲取

以AeroScout B1-100型單旋翼油動無人機作為飛行平臺，搭載Riegl VUX-1型激光掃描儀、OxTS Survey+2型慣性導航系統和Sony可見光相機獲取研究區LiDAR數據和可見光數據。無人機的長、寬、高為3.3 m×1.0 m×1.3 m，尾旋翼直徑0.65 m，空機質量50 kg，有效載荷18 kg，發動機功率13.23 kW，標準油箱容積10 L。Riegl VUX-1型激光掃描儀是專業測繪級激光掃描儀，掃描方式為線性掃描，近紅外波長，最大轉速200 r/s，激光脈沖頻率高達550 kHz，記錄16 bit的回波強度。Riegl VUX-1型激光掃描儀詳細參數見表1。

表1 Riegl VUX-1型激光掃描儀產品參數Tab.1 Product parameters of Riegl VUX-1 laser scanner

研究區域位于新疆兵團第八師一五零團五連周邊荒漠植被區，地理位置為44°57′29″～44°58′0″N，85°58′35″～85°59′4″E，激光雷達數據基本包含了荒漠植被區的部分植物類型，主要以梭梭、駱駝蓬、麻黃、堿蓬草、駝絨藜為主。獲取數據的飛行任務參數設置如下：航高60 m，巡航速度6 m/s，航線3條。在采用Riegl LMS配套軟件和OxTS NAVgraph配套軟件對激光雷達數據進行配準、校正、平差和刪除噪聲點等預處理后，截取部分區域數據作為試驗區，且以WGS-1984-UTM-Zone-44 N為投影坐標系導出.las格式的點云數據。為方便后期描述，令X、Y、Z、P分別表示由點云三維坐標和回波強度組成的N×1維矩陣，其中X和Y分別表示東向和北向的位置，單位為(°)，Z表示點云高程，單位為m。試驗區點云總數N為69 544，點云X中元素的范圍在-46.065°～-15.000°之間，Y中元素的范圍在83.958°～112.145°之間，Z中元素的范圍在275.863～280.535 m之間，P中元素的范圍在12 036～49 850之間。

2 研究方法

2.1 零-均值標準化

因為K-means聚類方法需使用不同距離來度量數據相似性，所以為了消除數據量綱和取值范圍差異的影響，需要將數據縮放到相同區間，即標準化[15]。同時，標準化后的數據也利于加快RBFNN的訓練速度，提高預測的準確性[16]。常用的標準化方法有3種：最大-最小標準化、零-均值標準化和小數定位標準化。其中，最大-最小標準化容易受到數據最大值和最小值的影響，小數定位標準化中移動的小數位數取決于絕對值最大的樣本點，故本文選擇零-均值標準化，與其它兩種標準化方法對比，零-均值標準化保持了異常值所包含的有用信息，標準化公式為

(1)

(2)

(3)

式中Snormal——點云屬性t標準值

Sorigin——點云屬性t原始值

μt——點云屬性t的均值

σt——點云屬性t的標準差

ti——點云屬性t的第i個點

2.2 K-means聚類方法

2.2.1K-means聚類原理

K-means聚類方法是把數據分成k簇，使簇內數據具有較高的相似度，而不同簇之間數據相似度較低[17]，采用K-means聚類方法實現點云聚類獲得地面點云的基本步驟如下：

(1)根據采集得到的研究區數據，確定需要的聚類數目k。

(2)確定k個聚類對象的初始聚類中心{Ui|i=1，2，…，k}。

(3)分別計算每個樣本數據到k個聚類中心的距離。

(4)篩選樣本數據到聚類中心的最小距離，且將此樣本劃分到聚類最小的類別中。

(5)當所有樣本點都劃分完畢，重新計算k個對象的聚類中心。

(6)與前一次獲取的k個聚類中心比較，判斷聚類中心是否發生變化或達到預定的最大迭代次數。如果聚類中心發生變化，轉到步驟(3)；如果聚類中心發生變化但已經達到最大迭代次數或者聚類中心沒有發生變化，則算法終止[18]。

通常，采用不同的距離標準可以得到不同的聚類結果。本文選擇適用范圍較廣的平方歐幾里距離作為距離標準，其計算公式為

(4)

式中d(x,y)——屬性t中樣本點云x與樣本點云y的歐幾里距離

xt——屬性t的樣本點云x

yt——屬性t的樣本點云y

2.2.2肘方法確定最佳聚類數目

K-means聚類方法的難點之一是必須事先指定聚類數目k，k值的選擇會嚴重影響聚類結果，故本文使用肘方法確定最佳聚類數目?；舅枷胧沁x擇簇內誤差平方和(Within-cluster sum of squared errors，SSE)(即聚類偏差)驟變時的k值[19]，此時對應的k值為最佳聚類數目。即先對不同聚類數目進行K-means聚類，然后根據聚類偏差得到不同k值的曲線，且繪制不同k值對應的聚類偏差圖，最后找到聚類偏差最顯著拐點處對應的k值，其中聚類偏差的計算公式為

(5)

其中

式中DSSE——簇j的聚類偏差

Bi——第i個點

U′j——簇j中心點

l——指定的最大聚類數目

K-means聚類方法的另一個難點是初始聚類中心的選擇。為了得到最優的聚類結果，本文在隨機選擇初始聚類中心和采用肘方法最佳聚類數目的基礎上，針對不同的隨機初始中心獨立運行K-means聚類方法多次，然后從中選擇DSSE最小的模型作為最終模型。

2.3 點云抽稀

海量地面點云給后期處理帶來較大影響，尤其是在數據處理速度方面[20]，故在保證研究對象的必要信息下，對高密度的點云進行抽稀有一定的實際意義。為了簡化抽稀過程，本文按照比例對點云進行抽稀，且僅對比抽稀率為20%和80%的抽稀結果。

Kriging插值是一種統計插值方法，可以準確對表面預測，不僅考慮了待估點與已知點位置的相互關系，而且考慮了變量的空間相關性[21]，所以采用Kriging方法插值抽稀結果，獲取較好的抽稀率。隨機抽取100個地面點云作為檢驗數據集，然后對剩余地面點云進行抽稀率為20%和80%的抽稀，對抽稀后的點云進行點云高程插值。最后從插值結果中提取檢驗數據集的高程預測值，且通過點云的預測值與實測值進行對比分析。

均方根誤差(RMSE)能夠評判插值結果，均方根誤差越小，表明插值精度越高，誤差越小[22]，其計算公式為

(6)

式中RRMSE——點云高程的均方根誤差

Zi——檢驗數據集中點云高程實測值

N′——檢驗數據集中點云數量

2.4 徑向基函數神經網絡模型

2.4.1RBFNN模型建立

RBFNN是應用多變量的徑向基函數設計而成，由輸入層、隱含層和輸出層組成。RBFNN的輸入層為n維向量I=(I1，I2，…，In)；隱含層為m維向量D=(D1，D2，…，Dm)；輸出層為f(I)，即輸出層是RBFNN的輸入層對應的實測值，以此實現RBFNN模型的建立[23]，其輸出公式為

(7)

其中

(8)

式中f(I)——徑向基函數神經網絡的輸出層函數

Wi——第i個隱含層神經元到輸出層神經元的權值

Di(I)——隱含層徑向基函數，采用高斯函數

Ci——第i個隱含層神經元中基函數的中心

ri——第i個隱含層神經單元的寬度，調節網絡的靈敏度

RBF神經網絡模型設計結構如圖1所示。

圖1 RBF神經網絡模型結構圖Fig.1 RBFNN model structure

采用軟件Matlab神經網絡工具箱的newrb函數設計RBFNN，在進行徑向基函數逼近時，newrb函數可以自動添加隱含層的神經元數量直到滿足訓練誤差要求，其實現公式為

Nnet=newrb(I，T，Ggo，Ssp，Mmn，Ddf)

(9)

式中Nnet——訓練得到的RBFNN

T——預測的目標數據

Ggo——RBFNN的目標誤差

Ssp——徑向基函數散布常數，默認值1

Mmn——最大神經元數量

Ddf——在網絡訓練過程中顯示的頻率，newrb中默認值為25

2.4.2RBFNN模型驗證

為了驗證模型的準確性，用測試點云對訓練完成后的RBFNN模型進行測試，且與實測高程值進行對比分析，來檢驗RBFNN對點云高程預測的準確性。采用軟件Matlab中sim函數，來實現測試點云高程的預測，計算公式為

A=sim(Nnet,Q)

(10)

式中A——點云高程預測值矩陣

Q——由測試數據集中點云坐標X、Y和回波強度P組成的矩陣

通過線性回歸法對預測值和實測值進行分析，線性回歸分析法是將實測值作為自變量，預測值作為因變量，建立一元線性回歸方程，檢驗回歸方程的決定系數R2。決定系數R2越接近1，則說明實測值和預測值吻合度越高，即預測效果越好。

3 結果分析

3.1 K-means聚類結果

采用軟件Matlab提取研究區點云三維坐標和回波強度值，且對點云進行零-均值標準化處理。采用Python語言編程確定最佳聚類數目和實現K-means點云聚類方法，圖2為聚類偏差結果，橫坐標表示聚類數目，縱坐標表示聚類偏差。分析圖2可知，當k為4時，聚類偏差呈現肘型，即對于研究區點云，最佳聚類數目為4，此時的聚類偏差DSSE為4 304.32。

圖2 聚類偏差可視化Fig.2 Within-cluster sum of squared errors visualization

采用K-means聚類方法對回波強度進行k為4的聚類，設置聚類最大迭代次數為3 000、聚類誤差為0.000 1。圖3a為采用K-means聚類得到不同類別可視化結果，數據1、數據2、數據3和數據4分別對應簇1、簇2、簇3和簇4。圖3b是簇1灰度圖，通過目視檢查分析可得，簇1包含研究區地面點云和離群點云，其點云總數為50 728，點云高度范圍在275.863～279.636 m之間。K-means聚類得到的地面點云完整地保留了試驗區真實地表起伏情況，但仍然存在少量離群點云，因此利用反復建立三角網濾出簇1中的少量離群點云，從而獲取地面點云數為48 722，此時點云高度范圍為275.870～277.600 m。

圖3 K-means聚類三維可視化點云圖Fig.3 Visualization point clouds by using K-means clustering method

3.2 克里金插值結果

通過計算均方根誤差評價抽稀精度，可知采用20%和80%的抽稀率對應的均方根誤差分別為0.021 m和0.023 m，即采用20%抽稀率抽稀可大大減少點云數據量。本文對地面點云進行抽稀處理后得到點云數為9 724。圖4為20%抽稀率抽稀處理后的克里金插值可視化結果，其中綠色點表示經過抽稀后的點云，紅色點表示100個檢驗點云。

圖4 20%抽稀率抽稀后的克里金插值的可視化結果Fig.4 Visual Kriging interpolation result at 20% point cloud thinning rate

3.3 徑向基函數神經網絡訓練結果

圖5 RBFNN訓練結果和線性回歸分析Fig.5 Training results of RBFNN and regression results of RBFNN

隨機采取抽稀后點云的70%作為RBFNN的訓練數據集，30%作為評估網絡的測試數據集。選擇訓練數據集中點云的X、Y坐標矩陣和回波強度矩陣P作為訓練RBFNN的輸入，即n=3，輸入向量I=(I1,I2,I3)。把訓練數據集中點云的高程Z作為網絡的輸出向量，即Z=f(I)。根據給定的目標誤差，多次訓練網絡確定隱含層數量。建立輸入層神經元數量為3，輸出層神經元數量為1，訓練誤差為0.02，散布常數為1的RBFNN。圖5a為RBFNN的訓練結果，可以看出，訓練次數達到100即隱含層神經元數量為100時，訓練誤差0.020 2很逼近目標誤差0.02，對應的RBF訓練時間為56 s。

利用訓練完成的RBFNN對測試數據集進行插值預測高程，記錄預測結果。同時采用線性回歸法進行分析，圖5b為RBFNN模型預測高程與實測高程的線性回歸分析，對應的決定系數R2為0.887，說明預測值與實測值的擬合度較高。均方根誤差RMSE為0.168 m，說明插值預測精度較高。

3.4 內插結果

根據訓練完成的徑向基函數神經網絡，對抽稀后的地面點云高程值進行內插，且根據高程值構建Delaunay三角網以生成試驗區DEM。圖6為采用Delaunay三角網生成的DEM，共構建97 393個三角網。

圖6 基于Delaunay三角網生成的DEM圖Fig.6 Generated DEM diagram based on Delaunay triangulation network

4 討論

本文目的是驗證利用RBF神經網絡內插值實現點云高程的預測，所以對K-means聚類結果僅采用了目視檢查來定性評價聚類結果。在后續研究中，可以計算漏分誤差、錯分誤差和總體誤差來定量評價聚類結果。另外，如果需要提高K-means聚類速度，可以先采用K-means聚類對標準化處理的點云三維坐標進行聚類分簇，然后針對不同的簇，再次使用K-means聚類對回波強度標準值進行聚類，且合并每類地面點云。

本文隨機選取100個點云測試Kriging插值結果，可知測試點分布較好，但是存在部分區域測試點分布較少，故后期可先計算點云曲率再抽稀。在建立RBF神經網絡模型時，選擇輸入變量為點云坐標矩陣X、Y和點云回波強度矩陣P，并沒有充分考慮到點云之間的相關性，因此RBF神經網絡輸入層數據的相關性還需進一步研究。

5 結論

(1)采用K-means聚類方法實現點云濾波獲得試驗區48 722個地面點云，采用肘方法確定最佳聚類數目為4，采用Python語言編程設置聚類數目、最大迭代次數3 000和聚類誤差0.000 1等參數，不僅可以快速得到聚類結果，而且可以將K-means聚類獨立運行多次，選擇最小聚類偏差DSSE為4 304.32作為最終結果。

(2)在保證用最少的數據量表示地面真實狀況下，選擇簡單、運算速度快的抽稀方法。在ArcGIS軟件中按比例對點云進行抽稀，且采用Kriging插值方法對不同抽稀后的點云進行插值。通過比較插值結果，得到較優抽稀率為20%。將RBF神經網絡引入到點云處理過程中，嘗試建立一種快速高效的空間預測方法，從誤差分析可以看出，RBF神經網絡對空間數據具有較好的預測能力，預測時間為56 s，其預測的決定系數R2為0.887。