改進(jìn)的布料模擬算法在多波束點(diǎn)云粗差剔除中的應(yīng)用

劉 毅， 王 磊， 王勝利， 胡亮亮， 高文龍

(1.山東科技大學(xué)測(cè)繪與空間信息學(xué)院， 青島 266590； 2.自然資源部，海洋測(cè)繪重點(diǎn)實(shí)驗(yàn)室， 青島 266590；3.山東科技大學(xué)海洋科學(xué)與工程學(xué)院， 青島 266590)

多波束探測(cè)技術(shù)已經(jīng)成為目前水下地形探測(cè)的主要方法[1]。近些年隨著近岸淺水調(diào)查及淺水水下地形項(xiàng)目需求的日益增長(zhǎng)，淺水多波束測(cè)深系統(tǒng)得到了廣泛的應(yīng)用[2-3]。在水下地形測(cè)量過(guò)程中，淺水多波束安裝于船體側(cè)舷容易受到噪聲、復(fù)雜的海洋環(huán)境及設(shè)備參數(shù)設(shè)置不合理等多種因素的影響，測(cè)深數(shù)據(jù)中難免含有粗差[4]。

為了獲得精確的水下地形，需要剔除多波束數(shù)據(jù)粗差，現(xiàn)有的濾波方法有人工剔除粗差和自動(dòng)濾波兩種。人工剔除粗差方法準(zhǔn)確度較高，但是需要消耗大量的勞動(dòng)力，且效率較低，無(wú)法滿(mǎn)足大量數(shù)據(jù)處理[5]。因此許多專(zhuān)家學(xué)者提出自動(dòng)濾波方法。李宏武等[6]利用基于最小二乘支持向量機(jī)(least squares support vector machines，LS-SVM)算法探測(cè)深度異常值；Calder等[7]提出CUBE(combined uncertainty and bathymetry estimator)算法；趙荻能等[8]提出了一種CUBE、曲面濾波參數(shù)聯(lián)合選優(yōu)流程；趙祥鴻等[9]提出利用BP神經(jīng)網(wǎng)絡(luò)剔除多波束測(cè)深數(shù)據(jù)粗差；楊安秀等[10]提出面向多波束測(cè)深數(shù)據(jù)的雙向布料模擬自動(dòng)濾波方法。

Zhang等[11]提出了基于布料模擬濾濃(cloth simulation filtering，CSF)算法對(duì)陸地激光點(diǎn)云濾波，其方法需要調(diào)整參數(shù)較少，精度高，適用不同復(fù)雜地形，給激光點(diǎn)云濾波提供了新的思路。本文中則利用CSF算法剔除多波束點(diǎn)云粗差，如果將CSF算法直接用于多波束點(diǎn)云濾波，會(huì)導(dǎo)致對(duì)凸起礁石和凹陷溝槽等水下地形過(guò)度濾波。因此在CSF算法基礎(chǔ)上，提出一種基于改進(jìn)的CSF算法多波束數(shù)據(jù)濾波來(lái)解決對(duì)于復(fù)雜海底地形過(guò)度濾波問(wèn)題。首先，對(duì)多波束點(diǎn)云數(shù)據(jù)進(jìn)行點(diǎn)云分割；然后，根據(jù)點(diǎn)云數(shù)據(jù)的面積和標(biāo)準(zhǔn)差調(diào)整每塊布料濾波算法的參數(shù)；其次，構(gòu)建布料濾波模型剔除粗差點(diǎn)；最后，將分割后的點(diǎn)云拼接生成水下地形點(diǎn)云數(shù)據(jù)，解決水下地形起伏較大時(shí)過(guò)度濾波的問(wèn)題。多波束測(cè)量系統(tǒng)是在不斷運(yùn)動(dòng)狀態(tài)下進(jìn)行的一種動(dòng)態(tài)測(cè)量，由于測(cè)量平臺(tái)受海浪起伏、風(fēng)、流等海洋環(huán)境效應(yīng)的干擾，多波束測(cè)深數(shù)據(jù)采集過(guò)程難免出現(xiàn)假信號(hào)，從而采集的數(shù)據(jù)與實(shí)際地形存在差異，這就是海洋測(cè)量信息獲取過(guò)程中的粗差問(wèn)題[12]。因此，粗差的判別非常重要，通常水下地形是連續(xù)變化，本文中利用布料模擬連續(xù)變化的水下地形，將離散點(diǎn)和異常線狀點(diǎn)云判定為粗差。

1 改進(jìn)的布料模擬濾波算法

1.1 布料模擬濾波算法原理

Zhang等[10]提出的CSF算法是利用布料模擬技術(shù)進(jìn)行濾波。布料模擬技術(shù)可以想象成一塊布放在一個(gè)水下地形的上方，布?jí)K在重力作用下落下，如果布料足夠柔軟，則可以貼在地形的表面上，最后布?jí)K的形狀就是地形。基于布料模擬技術(shù)，Zhang等[10]開(kāi)發(fā)了CSF算法用于激光雷達(dá)提取地面點(diǎn)，首先將原始點(diǎn)云翻轉(zhuǎn)，模擬一塊布落到翻轉(zhuǎn)后的地面上，最終確定布?jí)K的形狀，并將其形狀作為點(diǎn)云分類(lèi)的依據(jù)。本文中則利用CSF算法剔除多波束點(diǎn)云數(shù)據(jù)粗差點(diǎn)(圖1)。

圖1 布料模擬濾波算法原理

在布料模擬過(guò)程中，布料是由具有質(zhì)量和互連的粒子組成的網(wǎng)格(圖2)。網(wǎng)格節(jié)點(diǎn)上的粒子沒(méi)有大小，質(zhì)量相同且不變，粒子所在空間中的位置決定了布料的形狀。布料格網(wǎng)中連接兩個(gè)粒子的結(jié)構(gòu)稱(chēng)之為“虛擬彈簧”并且遵循胡克定律。

圖2 布料格網(wǎng)模型

布料格網(wǎng)中粒子的位置和速度由作用于粒子上的力決定，由牛頓第二定律可得，粒子位置與力之間的關(guān)系為

(1)

式(1)中：X(t)表示粒子在t時(shí)刻的位置；Fe(X,t)表示粒子所受外力，由粒子自身的重力G和與地形接觸的碰撞力Fc組成，粒子在運(yùn)動(dòng)方向上碰撞到一些物體則會(huì)產(chǎn)生碰撞力；Fi(X,t)表示粒子在X位置和t時(shí)間的內(nèi)力，由粒子間相互作用的“虛擬彈簧”產(chǎn)生(圖3)。

圖3 虛擬彈簧

建立布料格網(wǎng)模型時(shí)，首先，粒子的運(yùn)動(dòng)被限制在垂直方向，因此可以通過(guò)比較粒子和地形的高差來(lái)檢測(cè)碰撞；然后，當(dāng)一個(gè)粒子到達(dá)地形時(shí)，這個(gè)粒子就被設(shè)定為不可移動(dòng)；最后，將力的分成離散的兩個(gè)步驟，先計(jì)算重力對(duì)粒子的影響，后根據(jù)“虛擬彈簧”產(chǎn)生的內(nèi)力修改粒子的位置(圖4)。

圖4 布料格網(wǎng)節(jié)點(diǎn)位移計(jì)算

僅從重力計(jì)算每個(gè)粒子的位置，即內(nèi)力得零時(shí)可得

(2)

式(2)中：m表示粒子質(zhì)量；g表示重力加速度；X(t)表示t時(shí)刻節(jié)點(diǎn)粒子位置；Δt表示時(shí)間長(zhǎng)。當(dāng)已知節(jié)點(diǎn)粒子初始位置與時(shí)間步長(zhǎng)，則可以解算出格網(wǎng)節(jié)點(diǎn)的位置。

為了限制粒子在反轉(zhuǎn)的點(diǎn)云空隙區(qū)域的位移，需要考慮到粒子間的內(nèi)力。因?yàn)閮?nèi)力作用，粒子會(huì)停留在網(wǎng)格中，并且回到初始位置。遍歷所有“虛擬彈簧”，對(duì)于每個(gè)彈簧比較形成這個(gè)彈簧的兩個(gè)粒子之間的高差。如果兩個(gè)粒子都是可以移動(dòng)，在相反方向上移動(dòng)相同的量；如果一個(gè)粒子不可移動(dòng)，另一個(gè)將被移動(dòng)。每個(gè)粒子的位移計(jì)算公式為

(3)

式(3)中：d為粒子的位移量；當(dāng)粒子可以移動(dòng)時(shí)b=1，不可移動(dòng)時(shí)b=0；pi為p0的相鄰粒子的位置，n為點(diǎn)進(jìn)行標(biāo)準(zhǔn)化到垂直方向上的單位向量(0,0,1)T。

1.2 算法改進(jìn)

1.2.1 水下地形點(diǎn)云分割

水下地形是連續(xù)變化的，由于海洋水體的運(yùn)動(dòng)會(huì)不斷地侵蝕水下地形，并且一條測(cè)線覆蓋范圍相對(duì)較廣，水下會(huì)有凹陷溝槽或者凸起礁石等地形。如果對(duì)一整條測(cè)線都使用相同參數(shù)進(jìn)行CSF算法濾波，則容易在地形起伏較大處剔除水下地形點(diǎn)。因此需要對(duì)水下地形點(diǎn)云分割，當(dāng)分割后點(diǎn)云面積大于S1=10 000 m2，則再對(duì)該點(diǎn)云進(jìn)行分割，分割完成后可以針對(duì)不同區(qū)域調(diào)整參數(shù)進(jìn)行CSF算法濾波(圖5)。

圖5 點(diǎn)云分割

1.2.2 自適應(yīng)參數(shù)選取

水下地形數(shù)據(jù)濾波時(shí)，CSF參數(shù)調(diào)整是決定濾波質(zhì)量的關(guān)鍵因素。CSF算法需要調(diào)整參數(shù)為6個(gè)分別為：是否有斜坡St、布料格網(wǎng)分辨率GR、布料剛性RI、時(shí)間步長(zhǎng)dT、距離閾值hcc、迭代次數(shù)n。

(4)

1.3 算法流程

基于改進(jìn)的CSF算法流程(圖6)，可以概括為以下步驟。

圖6 算法流程

步驟1基于PCL對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行分割，使每個(gè)點(diǎn)云面積小于S1。

步驟2基于點(diǎn)云數(shù)據(jù)的標(biāo)準(zhǔn)差設(shè)置濾波參數(shù)：是否有斜坡St；布料剛性RI；布料格網(wǎng)分辨率GR。

步驟3計(jì)算布料在重力作用下格網(wǎng)節(jié)點(diǎn)高程值，計(jì)算節(jié)點(diǎn)之間彈簧力作用下的節(jié)點(diǎn)高程值。

步驟4重復(fù)步驟(3)當(dāng)所有迭代次數(shù)大于設(shè)定n次，或者未達(dá)到n次且當(dāng)所有節(jié)點(diǎn)粒子迭代后均未移動(dòng)則迭代停止。

步驟5計(jì)算水下地形點(diǎn)云點(diǎn)與節(jié)點(diǎn)的高程差，若高程差小于或等于距離閾值hcc，則被分類(lèi)為水下地形點(diǎn)，否則為粗差點(diǎn)。

步驟6將所有水下地形點(diǎn)云拼接。

2 實(shí)驗(yàn)與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)中采用兩處多波束點(diǎn)云數(shù)據(jù)驗(yàn)證該算法有效性，分別為溫州測(cè)區(qū)與青島測(cè)區(qū)中條帶數(shù)據(jù)，其中青島測(cè)區(qū)平均水深41.15 m，數(shù)據(jù)共包含水深點(diǎn)9 204 749，條帶跨度較大，水下地形起伏較大；溫州測(cè)區(qū)數(shù)據(jù)平均水深9.25 m，數(shù)據(jù)共包含水深點(diǎn)3 551 232，條帶跨度較小，水下地形起伏較小。選擇這兩數(shù)數(shù)據(jù)主要由于分別代表了兩種典型的多波束點(diǎn)云數(shù)據(jù)粗差，兩次數(shù)據(jù)粗差點(diǎn)均由噪聲引起且無(wú)法避免，其中溫州數(shù)據(jù)粗差集中于條帶中央，呈線狀，而青島數(shù)據(jù)粗差主要位于條帶兩側(cè)，呈毛刺狀(圖7)。所有數(shù)據(jù)均利用Reson SeaBat T50-P多波束測(cè)深系統(tǒng)采集，表1為Reson SeaBat T50-P多波束技術(shù)參數(shù)指標(biāo)。

圖7 實(shí)驗(yàn)區(qū)域

表1 T50-P多波束測(cè)深系統(tǒng)主要技術(shù)參數(shù)指標(biāo)

2.2 結(jié)果分析

為了驗(yàn)證本文算法有效性，分別利用CSF算法、人工剔除粗差與本文改進(jìn)的CSF算法對(duì)這兩處多波束點(diǎn)云數(shù)據(jù)進(jìn)行濾波。在測(cè)試中，CSF算法參數(shù)設(shè)定為St=true，GR=1，RI=1，n=500，dT=0.5，hcc=0.65。

圖8、圖9分別采用三種濾波方式對(duì)兩處多波束點(diǎn)云得出的結(jié)果，通過(guò)對(duì)比分析可得：

圖8 基于溫州數(shù)據(jù)多種濾波方法比對(duì)

圖9 基于青島數(shù)據(jù)多種濾波方法比對(duì)

(1)三種濾波方式均可對(duì)多波束點(diǎn)云數(shù)據(jù)粗差點(diǎn)剔除，可以看出所有方法對(duì)平坦地形的水下粗差剔除效果較好。

(2)CSF算法由于參數(shù)固定，雖然對(duì)相對(duì)平坦地形濾波效果較好，但是對(duì)于復(fù)雜的水下地形難以達(dá)到理想的效果，產(chǎn)生了一定程度的過(guò)度濾波，導(dǎo)致部分水下地形數(shù)據(jù)的丟失，并且當(dāng)測(cè)線面積較大時(shí)，所需計(jì)算時(shí)間和計(jì)算機(jī)內(nèi)存成倍增加。

(3)人工剔除粗差時(shí)會(huì)遺漏少量粗差點(diǎn)。本文算法則彌補(bǔ)了CSF算法的劣勢(shì)，針對(duì)不同水下地形進(jìn)行自動(dòng)參數(shù)調(diào)整，凸顯了水下地形濾波的優(yōu)勢(shì)。通過(guò)對(duì)整體數(shù)據(jù)對(duì)比分析，驗(yàn)證了本文算法不僅能有效剔除多波束點(diǎn)云數(shù)據(jù)粗差也能更好地保留較大起伏的水下地形特征。

為了更好地量化本文算法優(yōu)勢(shì)，將這兩條測(cè)線分別采用人工剔除粗差、CSF算法與本文算法濾波，濾波前后的各項(xiàng)指標(biāo)如表2所示。與人工方法粗差點(diǎn)剔除比NR是指的利用算法剔除粗差點(diǎn)個(gè)數(shù)與人工剔除粗差點(diǎn)個(gè)數(shù)的比，即NR=Nn/Na，其中Na為人工剔除粗差點(diǎn)個(gè)數(shù)，Nn為算法剔除的粗差點(diǎn)個(gè)數(shù)[10]。

由表2可以看出，采用本文算法和CSF算法得到的NR，對(duì)于相對(duì)平坦區(qū)域和起伏較大區(qū)域分別從1.03下降到1.01和從1.46下降到1.19。從表2中可以看出人工剔除粗差效果雖然很好，但是最小水深均要比兩種算法要小，經(jīng)檢查發(fā)現(xiàn)有遺漏粗差點(diǎn)；CSF算法在平坦地形條件下濾波效果較好，但是在起伏較大地形條件下存在過(guò)度濾波問(wèn)題，導(dǎo)致局部地形數(shù)據(jù)丟失；本文算法表現(xiàn)較好，平坦和起伏較大地形均有良好的濾波效果，且多波束數(shù)據(jù)完整。

表2 多種方法濾波性能比對(duì)

3 結(jié)論

(1)人工剔除粗差效果較好，但是會(huì)遺漏粗差點(diǎn)，表現(xiàn)為最小水深不符合實(shí)際水下地形。

(2)CSF濾波算法對(duì)相對(duì)平坦的水下地形濾波有較好的效果，但是遇到有較大坡度的水下地形時(shí)，會(huì)出現(xiàn)過(guò)度濾波的現(xiàn)象。

(3)基于點(diǎn)云分割和自適應(yīng)參數(shù)調(diào)整提出了改進(jìn)的CSF算法，與CSF算法相比，克服了過(guò)度濾波現(xiàn)象，而且對(duì)于不同地區(qū)與人工方法粗差點(diǎn)剔除比分別從1.03下降到1.01和從1.46下降到1.19，保證了數(shù)據(jù)的完整性。與人工剔除粗差相比，避免了遺漏粗差點(diǎn)，人工干預(yù)很少。但是對(duì)于起伏更大地形，需進(jìn)一步驗(yàn)證。