一種基于CSF-WOA-LSSVM的匹配點云土石方量計算方法

何廣煥， 唐詩華*， 王文貫， 張炎， 劉銀濤， 蒙金龍

(1.桂林理工大學測繪地理信息學院， 桂林 541006; 2.廣西空間信息與測繪重點實驗室， 桂林 541006;3.廣西建設職業技術學院土木工程學院， 南寧 530007)

土石方量計算在工程施工項目、國土資源監測中占據著不可忽視地位，其計算結果的準確性是工程項目的施工進度、造價預算以及土地資源儲備量估算的重要影響因素[1]，如何準確高效地計算土石方量已成為當今眾多學者的研究方向。無人機傾斜攝影測量技術可快速獲取測區影像匹配密集點云，進而可利用密集點云建立數字高程模型或提取地面點進行土石方量計算[2-3]。當前，諸多學者將無人機傾斜攝影測量技術應用在土石方量計算上已取得一些階段性的研究成果。李博等[4]提出利用Photoscan軟件對無人機影像進行匹配密集點云生成，然后對通過點云分類和非地面點高程修正法生成的數字高程模型進行土石方量計算，高效、準確地獲取了土石方量計算成果。相詩堯等[5]利用五鏡頭影像對公路建立實景三維模型，通過重新構建模型表面三角網來剔除植被、房屋、樹木等影響土方計算的地物，并利用清表后的實景三維模型生成地面點云，進而計算公路土石方量。高利敏等[6]通過匹配點云生成兩期DEM模型，結合Virtual Surveyor軟件計算土方變化量，為土方監測提供了一種新方法。王春香等[7]利用改進的神經網絡對三維點云模型空洞進行修復，很好地還原了三維點云模型的真實性。李雅盟等[8]運用布料模擬濾波(cloth simulation filtering，CSF)法對激光點云進行地面點提取，通過設置相應的地形閾值取得了很好的效果。刑鵬威等[9]通過不規則三角網法對匹配密集點云進行地面點云提取，并利用最小二乘支持向量機(least squares support vector machine，LSSVM)對地面點云空洞進行修復，獲得了穩定、可靠的地面點云。

現有的研究已很好地促進了無人機傾斜攝影測量技術在土石方量計算中的應用與發展，但直接通過匹配密集點云計算土石方量的研究比較少，目前數字正射影像圖、數字高程模型、實景三維模型等低空攝影測量產品的成果質量均依賴于匹配密集點云的三維點位精度，由此可看出，匹配密集點云在精度上具有更大的優勢[10-11]。為了進一步促進無人機匹配點云在土方計算工作中的作業效率和成果精度，現提出通過CSF濾波法對匹配密集點云進行濾波、降噪、抽稀等處理提取地面點，然后利用經鯨魚優化算法(whale optimization algorithm，WOA)優化后LSSVM對所提取的地面點空洞進行插值修復，并將修復完成的地面點數據導入南方CASS軟件中進行土石方量計算，進而為土方工程提供一種更為快速、精度更為可靠的土方計算方法。

1 土石方量計算流程和算法原理

1.1 土石方量計算流程

為了高效、準確地采集土石方測區的地面點數據，利用無人機傾斜攝影測量技術、PPK技術、RTK技術對測區進行影像、控制點、基站衛星觀測文件等數據采集。利用Photoscan、PEN-PPK軟件對采集完成的數據進行處理生成匹配密集點云，對匹配密集點云進行CSF點云濾波、WOA-LSSVM地面點插值修復，進而獲取測區真實地面點數據，最后將處理完成的地面點數據導入南方CASS軟件進行土石方量計算。具體流程如圖1所示。

圖1 土石方量計算流程Fig.1 Calculation process of earthwork volume

1.2 CSF點云濾波原理

傳統分離地面點和非地面的濾波算法，大多數是依賴坡度、高程等變化因素來提取地面點，但需要用戶對點云數據具有豐富的先驗知識才能很好地設置濾波器中的各種參數。Zhang等[12]提出了一種CSF濾波法，該方法先把測區點云進行180°轉置，然后將一塊網格狀構成的布料貼附在轉置的地形點云上，布料在重力、鄰近節點作用力的相互拉伸下產生位移，當布料靜止時，所覆蓋的區域就代表當前測區的地形，如圖2所示。

圖2 布料模擬濾波基本原理Fig.2 Basic principles of cloth simulation filtering

具體表達公式為

(1)

式(1)中：m為粒子的重量，一般定義為1；X為粒子在時間t時的位置；Fext(X,t)為粒子受到的重力大小；Fint(X,t)為粒子受到鄰近節點作用力大小。重力和節點相互作用力隨時間 變化而變化，所以可用數值積分來求解。

1.3 LSSVM點云空洞插值修復模型

支持向量機屬于深度學習中的線性回歸法，通過訓練學習利用所有的支持向量得出最優分割平面，但在對偶問題上容易出現二次規劃現象，為克服其缺點，LSSVM利用最小二乘法將原不等式約束改為等式約束，化簡了計算流程，提高了原算法的泛化能力和收斂速度[13-14]，能讓支持向量機更適用于樣本多、面積大的點云空洞修復。其原理及修復過程如下。

給定點云訓練樣本集{(x1,y1),(x2,y2),…, (xn,yn)}, 則xn∈Rn表示高程點的平面二維坐標，yn∈Rn為高程一維輸出樣本。利用映射函數φ(x)把輸入的三維地面點數據映射到高維特征空間，并在該空間中構造優化的線性回歸函數，其數學模型為

y(x)=wTφ(x)+b

(2)

式(2)中：y(x)為高程預測值；wT式為訓練權重向量；b為偏差值。根據結構風險最小化原則，將其轉化為以下優化問題：

(3)

約束條件為

s.t.yi=wTφ(xi)+b+ξi,i=1,2，…，n

(4)

式中：c為正則化參數，其作用是為了控制地面點云模型的擬合度；ξ為懲罰變量，表示地面點云模型的擬合誤差允許值。通過拉格朗日函數和KKT條件優化，可得地面點空洞修復模型為

(5)

式(5)中：ai為拉格朗日乘子，其作用是為找到模型最優解；K(x,xi)為核函數，主要作用是減少計算量。根據Mercer條件，利用徑向基函數(RBF)作為地面點空洞修復模型的核函數，具體表達式為

(6)

式(6)中：σ為核函數參數；exp(·)為高斯核函數。

1.4 鯨魚優化算法

鐘明輝等[15]受鯨魚群捕食行為的啟發，提出了一種鯨魚優化算法。WOA通過模仿鯨魚對獵物的3種捕獵方式，以達到搜索獵物位置最優解的目的。該算法結構簡易，并且在收斂速度、全局尋優能力方面具有一定的優勢，其主要作用是為LSSVM點云插值模型提供可靠的正則化參數c和核函數參數σ，進而更好地修復地面空洞點。

1.4.1 環繞獵物

鯨魚環繞獵物過程中，座頭鯨的位置看作待優化問題的解，當座頭鯨確定獵物位置后，靠近獵物，并將其位置當作獵物位置，其他鯨魚不斷更新自己的位置并向座頭鯨位置靠近，形成環繞獵物的狀態，具體公式為

D=|C·X*(d)-X(d)|

(7)

X(d+1)=X*(d)-A·D

(8)

式中：·表示逐元素乘法運算；D為座頭鯨與獵物的距離參數向量；d為當前迭代數；X*(d)為獵物位置向量；X(d)為座頭鯨的位置向量；X(d+1)為其他鯨魚群的位置向量；C和A為兩個控制參數。

(9)

式(9)中：·表示逐元素乘法運算；r為[0,1]內的隨機數；TMaxlter為種群最大迭代次數；a隨鯨魚種群迭代次數的增加從2～0線性遞減。

1.4.2 螺旋氣泡攻擊

鯨魚進行螺旋氣泡攻擊時，先計算自身到獵物位置的距離并慢慢逼近，當到達獵物位置對其進行螺旋氣泡攻擊，具體公式為

D′=|X*(d)-X(d)|

(10)

X(d+1)=D′·ebl·cos(2πl)+X*(d)

(11)

式中：D′為當前鯨魚與獵物的距離參數；b為1的螺旋狀常數；l為[-1,-2]內的隨機數。

1.4.3 全局搜索獵物

為了避免鯨魚在捕食獵物的過程中出現局部最優解的情況，鯨魚群可通過全局隨機搜索的方式實現捕食，根據獵物位置的變化而不斷更新自身位置，具體公式為

D=|C·Xrand-X(d)|

(12)

X(d+1)=Xrand-A·D

(13)

式中：Xrand為當前隨機獵物的位置。

2 工程實例

2.1 測區概況

本工程為廣西陸川縣某花崗巖石礦場的土石方量測量項目，測區南北長約480 m，東西寬約550 m，地形相對高差約為110 m，占地面積約為200畝(1畝=667 m2)，大部分區域地表呈裸露狀態。測區內部含有部分建筑物、棚房、灌木等地物以及石礦生產線中的大型設備。為驗證文章土石方計算方法的可行性，均勻選取16個被地物遮擋的隱蔽高程點作為檢查點，并通過RTK技術獲取其坐標信息。測區及檢查點點位分布的具體情況如圖3所示。

2.2 數據采集與預處理

2.2.1 數據采集

利用大疆精靈4RTK單鏡頭無人機對測區進行影像數據采集，根據對測區現場的踏勘情況同時結合飛行設備的硬件參數，采用傾斜攝影3D技術、PPK定位技術、RTK定位技術對測區進行數據采集，共獲取測區743張相片、5個點位分布均勻的像控點、1個基站架設點的靜態數據及其厘米級的三維坐標。其中航線設計的具體參數為：相對航高為80 m、地面分辨率為2.1 cm、飛行速度為7.1 m/s、航向和旁向重疊度為80%、云臺角度為-60°。像控點及基站點坐標信息通過RTK技術獲取，并對基站點進行靜態數據采集，為后期影像POS數據解算提供衛星觀測數據。

2.2.2 數據預處理

數據預處理的主要目的是獲取測區可靠的匹配密集點云，為后期點云濾波、地面空洞插值修復、土石方量計算提供準確的原始點云數據。首先，利用PEN-PPK軟件對導入的原始影像POS數據、靜態觀測數據、厘米級基站點坐標進行聯合解算，進而獲取高精度的影像POS數據。其次,將準確的POS數據、傾斜相片、像控點導入Photoscan軟件進行相對定向、絕對定向、空三加密等處理，進而生成精度可靠的匹配密集點云，具體情況如圖4所示。

圖4 匹配密集點云Fig.4 Match dense point cloud

2.3 CSF點云濾波降噪

匹配密集點云往往含有許多噪聲點，除了地面點云之外，還含有建筑物、植被、湖泊、路燈等非地面點云。本文研究利用CSF濾波法對測區提取地面點云，該算法只需要設置4個參數：時間步長為0.65 s，地面點閾值為0.2 m、格網分辨率為0.2 m、迭代次數為300次。其中，時間步長指布料粒子位置更新所需要的時間；地面點閾值指將具離地面有一定高度的離散點判斷為地面點；網格分辨率指布料網格大小，一般設置為點云分辨率的2～3倍。經CSF濾波提取的地面點云如圖5所示。

圖5 具有空洞的地面點云Fig.5 Ground point cloud with cavity

2.4 地面點空洞修復

2.4.1 WOA-LSSVM插值

通過CSF濾波提取的地面點云，基本剔除了建筑物、植被、大型設備等非地面點云，同時也造成了部分區域產生地面空洞，直接利用該點云數據進行方量計算會使計算結果遠偏離土石方量真實值，所以需對地面點云空洞進行修復，盡可能還原地表真實地形。采用WOA-LSSVM對地面點空洞進行插值修復,并且與LSSVM、Kriging插值法以及RTK實測數據進行比較分析。由于點云數據量巨大，且為了提高軟件的計算速度，則在修復之前先對提取的地面點云進行抽稀，獲取平均間距為3 m的地面點。在WOA算法優化參數的過程中，將LSSVM插值修復模型的點位中誤差作為獵物位置更新的適應度函數，鯨魚種群數量設置為60頭，最大迭代次數為100，螺旋狀常數設置為1，為LSSVM插值修復模型提供該研究區域最合適的正則化參數和核函數參數。最后把待插修復點的平面坐標導入訓練完成的三維點云擬合模型中進行高程預測。圖6為經WOA-LSSVM修復后地面點三維模型。

圖6 修復后的地面三維模型Fig.6 Repaired 3D ground model

2.4.2 插值修復模型精度分析

由于測區高差過大，并且存在多處高陡崖，所以繪制了全測區的3 m等值線圖和局部測區1 m等值線圖來驗證經WOA-LSSVM插值修復后的地面點可靠性，并且與LSSVM插值法、Kriging插值法以及通過RTK技術按5 m采樣間隔獲取的實測數據做對比分析，如圖7和圖8所示。

由圖7可知，相比RTK實測數據的3 m等值線圖，LSSVM插值法、Kriging插值法存在少量等值線局部突變的情況，并且Kriging插值法在測區中部出現等值線分層，表明兩種方法所插值修復的地面點數據存在的離散誤差點較多。由圖8可知，相比RTK實測數據的1 m等值線圖，Kriging插值法、LSSVM插值法在局部測區的左下方出現大范圍等值線突起，并且在許多地方出現大小不一的等值線突變，表明存在許多面狀區域性的地面點粗差。而WOA-LSSVM修復法在兩種等值線圖則幾乎沒有出現等值線分層、突起、變異的情況，并且等值線的走向與RTK實測數據的等值線圖更為貼近，表明該方法具有一定的可靠性和穩定性。為了進一步驗證WOA-LSSVM插值修復法的精確性，選取16個被地物遮擋的高程點作為檢查點，進行精度統計分析，結果如表1所示。

由表1可知，WOA-LSSVM插值法的檢查點殘差變化區間為-0.317～0.214 m，LSSVM插值法為-0.344～0.224 m，Kriging插值法為-0.519～0.449 m。在整體精度上，WOA-LSSVM插值法的平均誤差和中誤差分別為：0.105 m和0.137 m，LSSVM插值法為：0.148 m和0.172 m，Kriging插值法為：0.135 m和0.205 m。相較于其他兩種方法，經WOA優化后的LSSVM插值法的檢查點殘差變化區間最小、效果更好，在精度方面也更具優越性。為了更為清晰直觀地展現三種插值修復法的變化趨勢和穩定性，繪制了檢查點殘差的變化折線圖，結果如圖9所示。

圖7 4種方法的全測區3 m等值線圖Fig.7 3 m contour map of the whole survey area by four methods

圖8 4種方法的局部區域1 m等值線圖Fig.8 Local 1 m contour map of four methods

表1 3種修復法的檢查點精度統計表Table 1 Statistical table of checkpoint accuracy of three repair methods

圖9 檢查點殘差對比圖Fig.9 Comparison diagram of checkpoint residuals

由圖9中的高程殘差波動趨勢可知，相對于LSSVM插值法、WOA-LSSVM插值法，利用Kriging插值法修復的檢查點高程殘差波動范圍明顯更大，而經過WOA優化的LSSVM插值法所修復的隱蔽點高程殘差相比原插值方法更為穩定、波動性更小、更趨近于真值。

2.5 土石方量計算

為驗證本文方法在土石方量計算中的準確性和可靠性，利用RTK技術按5 m間隔對測區進行地面點數據采集，并在研究區植被較茂密、地勢起伏大的區域進行適當的地面點加密采集，然后將CSF-WOA-LSSVM法、CSF-LSSVM法、CSF-Kriging法、RTK技術所處理和實測的地面點數據導入南方CASS軟件中，設置相同的標高及相關參數，以及確定相同的區域邊界線，通過方格網法計算測區的土石方量，結果如表2所示。

通過對4種方法所計算的土石方量結果進行分析比較可知：

(1)通過對無人機匹配點云進行濾波、空洞插值修復所完成的土方計算工作，在工作效率方面要比傳統RTK土石方測量更具優越性，將工作時長縮短了近3倍，并且大大減少了外業工作時間。

(2)與通過傳統RTK技術實測計算所獲取的土石方量結果相比，經過濾波、插值的3種方法均能滿足土石方量計算的規范要求，其中CSF-Kriging法的土石方量和差值比分別為：98 437.4 m3和2.03%；CSF-LSSVM法為90 150.2 m3和1.86%；CSF-WOA-LSSVM法為：67 073.6 m3和1.38%，并且在3種方法中準確性最高，更適用于在土方工程。

表2 土方計算結果對比Table 2 Comparison of earthwork calculation results

(3)利用WOA-LSSVM插值法修復的地面點數據進行在土石方量計算，所計算的結果相比LSSVM插值法的準確性提高了26%，這主要是因為WOA算法為LSSVM插值法提供了更為準確的正則化參數和核函數參數，從而提高了土石方量計算的準確性。

3 結論

在土石方工程中，如何快速、準確地計算土石方量，對整個項目的施工進度和成本控制具有重要影響。利用無人機攝影測量技術生成研究區的匹配點云，通過CSF濾波法對其進行地面點提取，利用WOA-LSSVM插值法、LSSVM插值法、Kriging插值法對所提取的地面點空洞進行修復，最后將3組插值修復完成的地面點導入南方CASS軟件進行土石方量計算，并結合RTK實測土方數據對3種方法進行了精度驗證和適用性分析，得出以下結論：①利用無人機影像匹配點云進行土石方量計算，得益于其不受研究區域環境障礙的影響，使得外業數據采集的時間大大縮短，進而有效地提高了工作效率；②CSF濾波算法整體參數設置簡單，處理過程中無需過多的先驗性知識，并且可較好地提取匹配點云中的地面點；③對具有孔洞的匹配地面點云，利用WOA-LSSVM插值法對其進行修復可獲得到穩定、可靠的完整地面匹配點云，具體精度達到了0.137 m；④將WOA-LSSVM插值法所修復完成的地面匹配點云用于土石方量計算，與RTK技術獲取的土石方量計算結果相比更為接近，差值比為1.38%，完全滿足土方計算的規范要求，可為土石方工程提供一定的參考依據。