二次CUBE濾波算法及在邊坡亂石區(qū)測深數(shù)據(jù)處理中的應(yīng)用

鄭永新，張紅梅，趙建虎

1. 武漢大學(xué)測繪學(xué)院，湖北 武漢 430079； 2. 武漢大學(xué)海洋研究院，湖北 武漢 430079； 3. 武漢大電氣與自動化學(xué)院，湖北 武漢 430072

多波束測深系統(tǒng)(multibeam bathymetric system,MBS)因具有高密度、全覆蓋測深特點[1-2]，在港口、航道等工程中有著廣泛的應(yīng)用。實際應(yīng)用中，特別在海岸防波堤拋石工程中需要邊施工邊測量，對拋石間較大縫隙進行填補，這就需要比較準(zhǔn)確和清晰的水下地形信息。然而，施工區(qū)域水下地形測量受懸浮物、石塊遮擋、系統(tǒng)噪聲和測量船噪聲等影響，MBS獲取的測深點云存在大量粗差，嚴(yán)重影響了測深數(shù)據(jù)對海底測量對象的準(zhǔn)確描述，不能給予拋石施工準(zhǔn)確的指導(dǎo)信息，必須給予有效濾除。據(jù)此，國內(nèi)外學(xué)者開展了大量研究。目前，常見濾波方法主要有COP法[1]、Ware法[1]、Knight&Wells法[1]、Eag(RDANH)法[3]、趨勢面法[1,4-5]、抗差估計法[6-7]、Bayes估計[8]、中值濾波[9]、均值濾波、局部方差檢測和小波分析相結(jié)合的濾波方法[10]以及選權(quán)迭代加權(quán)平均[11]等方法。上述方法對于海量MBS點云數(shù)據(jù)處理存在速度慢，無法滿足邊施工邊測量，近實時獲取水下地形的需求，且普適性不高，亂石區(qū)等復(fù)雜海床測深數(shù)據(jù)濾波性能欠佳等不足，導(dǎo)致無法準(zhǔn)確獲取水下地形信息，無法正確指導(dǎo)拋石施工。在文獻[12]提出MBS測深點的水平和垂直不確定度理論計算方法的基礎(chǔ)上，文獻[13—15]提出了CUBE濾波方法。由于該方法具有濾波高效、可靠、抗差、穩(wěn)定等特點[16-19]，目前在CARIS和HYPACK等國際商用MBS測深數(shù)據(jù)處理軟件中廣泛使用[20]。當(dāng)測深點在格網(wǎng)節(jié)點周圍均勻分布，存在少量粗差時，CUBE算法可準(zhǔn)確估計格網(wǎng)節(jié)點水深值。然而，在亂石區(qū)，現(xiàn)有CUBE算法會遇到如下問題：其一，MBS測量會因亂石遮擋產(chǎn)生測量盲區(qū)，導(dǎo)致盲區(qū)測深數(shù)據(jù)空白，波束可達區(qū)密集，測深數(shù)據(jù)分布(特別是在邊坡亂石區(qū))不均勻性問題突出，導(dǎo)致現(xiàn)有CUBE算法對邊坡亂石區(qū)測深數(shù)據(jù)處理時估計節(jié)點的深度不準(zhǔn)確；其二，現(xiàn)有CUBE算法僅開展深度估計，未顧及平面位置估計，亂石區(qū)的上述數(shù)據(jù)特征會引起估計的水深點位置不準(zhǔn)確，從而引起地形特征混亂；其三，邊坡亂石區(qū)地形變化復(fù)雜，波束在亂石間會產(chǎn)生多次回波后到達換能器，或因環(huán)境復(fù)雜，水中存在大量遮擋物時導(dǎo)致測深數(shù)據(jù)中存在大量粗差，基于現(xiàn)有CUBE算法的最優(yōu)水深估值選擇原則[13]將很難得到準(zhǔn)確的水深，還需要借助手工編輯剔除粗差，也因此顯著降低了數(shù)據(jù)濾波的效率，嚴(yán)重影響施工的進度。

為此，針對CUBE算法不足，在深入研究的基礎(chǔ)上，提出了顧及平面位置估計和開展二次CUBE估計的濾波算法，以期解決上述問題，實現(xiàn)MBS測深數(shù)據(jù)對亂石區(qū)地形的真實反映。

1 CUBE算法原理及其不足

1.1 CUBE算法

CUBE是基于規(guī)則格網(wǎng)估計各格網(wǎng)節(jié)點最優(yōu)水深的濾波方法，主要包括如下3個部分。

1.1.1 測深點深度不確定度傳遞

在根據(jù)水深確定格網(wǎng)節(jié)點搜索范圍內(nèi)，將測點水平和垂直不確定度σH,j、σV,j傳遞到待估格網(wǎng)點上[13]

(1)

式中，δij為測深點到待估格網(wǎng)點的平面距離；SH為水平不確定度的影響因子，默認(rèn)為1.96；α為地形起伏影響因子，默認(rèn)為2；Δmin為格網(wǎng)大小。

1.1.2 多重估計

以Kalman濾波方法建立估計方程[13]，式(2)和式(3)分別給出了狀態(tài)方程與觀測方程

z[n+1]=z[n]+w[n]w[n]～N(0,W[n])

(2)

(3)

式中，w[n]、v[n]為系統(tǒng)噪聲和觀測噪聲，均服從正態(tài)分布。由于模型是為了獲取格網(wǎng)節(jié)點的唯一水深，因此有W[n]=0 m2(文獻[13])。

式(4)—式(9)給出了濾波估計過程[13]

(4)

(5)

(6)

(7)

(8)

(9)

定義圖1(b)中σf為預(yù)測誤差，圖1(c)中d1和d2分別為輸入水深值與估值1、估值2差值，則CUBE算法多重估計更新過程為：

(1) 初始化一個估值后或者只存在一個估值時，進行如圖1(a)或圖1(b)所示的過程。

(2) 存在多個估值時(如圖1(c))選擇與輸入水深值差異最小估值，重復(fù)圖1(a)或1(b)。

(3) 重復(fù)步驟(1)和(2)，對所有數(shù)據(jù)格網(wǎng)節(jié)點進行水深估計。

圖1 CUBE多重估計更新Fig.1 Multiple estimation and update in CUBE filtering

1.1.3 最優(yōu)水深估值選擇原則[13]

原則1：選取估值納入點數(shù)最多的一個估值作為最優(yōu)估值。

原則2：在待估格網(wǎng)點周圍一定范圍內(nèi)，尋找單個估計值的格網(wǎng)點水深為參考水深，尋找待估格網(wǎng)節(jié)點估值中與該參考水深最接近的估值作為最優(yōu)估值。

原則3：綜合原則1和2。

1.2 邊坡亂石區(qū)CUBE算法的不足

(1) 忽略了邊坡亂石區(qū)測深點的分布特點。

邊坡亂石區(qū)MBS測量時，因大坡度地形、石塊(圖2(c))遮擋，導(dǎo)致大量測量盲區(qū)(圖2(a))，測深點在垂直方向上間斷分布(圖2(a)、(c))?，F(xiàn)有CUBE算法基于測點均勻分布(圖2(b))和多重估計原理開展估值圖2(c)中的測點分布不均勻問題會導(dǎo)致對節(jié)點Node的估計產(chǎn)生兩個水深估值A(chǔ)和B。無論是選擇A估值還是B估值，最終會賦值給Node作為該節(jié)點的估值，進而造成節(jié)點平面位置和水深估值不準(zhǔn)，引起地形特征錯位甚至錯誤。

(2) 無法消除大量粗差影響。

傳統(tǒng)CUBE濾波給出了3個水深估值優(yōu)選原則，其優(yōu)缺點如下：

原則1：簡單，易于實施，當(dāng)?shù)匦纹骄徍蜏y深數(shù)據(jù)中存在少量粗差時，可有效剔除粗差，但當(dāng)粗差較多且集中分布時，易獲得錯誤水深估值。

原則2：當(dāng)?shù)匦纹骄?、粗差較少或者離散分布時，可有效消除粗差，具有抗差性，但當(dāng)?shù)匦螐?fù)雜、數(shù)據(jù)質(zhì)量較差時，多重估計后將有大量的格網(wǎng)節(jié)點出現(xiàn)多個估值情況，此時很難找到只有單個估值情況。如取舍不當(dāng)，會導(dǎo)致待估格網(wǎng)點水深錯誤。

原則3：情況同原則1和2。

實際數(shù)據(jù)處理時，如圖3(a)中紅色區(qū)域為邊坡亂石區(qū)，采用CUBE算法對該區(qū)水深數(shù)據(jù)進行濾波，每個格網(wǎng)節(jié)點將獲取多個估值，因此原則2、3無法實施，只能采用原則1選擇最終估值。然而，針對圖3(b)矩形區(qū)測深數(shù)據(jù)中存在大量粗差，借助原則1獲取的估值如圖3(c)，可看出傳統(tǒng)CUBE并未完全消除粗差的影響。

2 顧及數(shù)據(jù)特征的二次CUBE濾波算法

在邊坡亂石區(qū)，為解決因原始波束數(shù)據(jù)中存在大量粗差產(chǎn)生的由傳統(tǒng)CUBE濾波帶來的節(jié)點估值不準(zhǔn)問題(如圖3(c)所示)，提出一種快速二次CUBE濾波方法。

2.1 顧及位置估計的一次CUBE濾波

針對1.2節(jié)中問題(1)，本文提出在傳統(tǒng)CUBE算法的基礎(chǔ)上，在進行水深估計的同時對水深估值的平面位置估計。平面位置估計與水深估計類似，將水深值z用平面坐標(biāo)(x、y)代替，用每個波束的水平不確定度代替水深不確定度。采用類似式(5)、式(7)和式(8)的形式，則有

(10)

(11)

(12)

在傳統(tǒng)CUBE濾波方法的基礎(chǔ)上，增加位置估計，對原始的測深數(shù)據(jù)開展一次濾波。顧及平面估計的一次CUBE濾波所得散點(圖4中所示的三角網(wǎng)節(jié)點)不再是每個格網(wǎng)節(jié)點。

2.2 二次CUBE濾波

提取一次CUBE濾波所得散點水深值、不確定度及對應(yīng)的每個點上多重估計水深。如圖4所示，圖中中心點(xm,yn)為待檢測點，以該點為中心，直徑為a的圓范圍內(nèi)測深點對該點估值進行檢測，a=2x×Δmin(x≥1,x∈N)。節(jié)點上估值的垂直不確定度為

(13)

式中，k表示節(jié)點(xi,yj)的第k個估計值；nk表示該估值的吸收的測深點數(shù)；σv,j表示該估值納入的第j個測深點的垂直不確定度。

將周圍點(xi,yj)的不確定度σ(i,j)傳遞到待估點(xm,yn)上。本文顧及平面位置影響和亂石區(qū)地形梯度變化顯著特點，給出傳遞模型如下

(14)

當(dāng)水深小于50 m，波束大小為0.5°×1°，儀器量程分辨率為1.25 cm，條帶覆蓋角120°，假定節(jié)點不確定度σ(i,j)為0.15 m，則根據(jù)式(14)，表1列出了當(dāng)周圍點距待估點距離為1Δmin時不同μ將對應(yīng)不同大小的σi,j值。由表1可知，假設(shè)相鄰節(jié)點水深值差異為0.5 m，為容許該水深差，σi,j應(yīng)該大于0.5 m，此時μ應(yīng)大于3，取值4。μ不能過大，否則可能會導(dǎo)致粗差濾除失敗。μ值參數(shù)需根據(jù)當(dāng)前測量水深的不確定度和地形起伏來綜合確定。

σ(i,j)/mμσi,j/m0.1530.4740.1540.6180.1550.765

然后，采用多重估計對水深點值與周圍一定范圍內(nèi)的水深值進行一致性檢測(如圖1)。若發(fā)現(xiàn)該水深值帶有粗差，利用周圍水深值采用多重估計和最優(yōu)參考水深選擇方法獲取最優(yōu)水深參考值dr。dr的選取過程如下：

(1) 將一次CUBE濾波獲得的待估點水深作為初始估值d0，利用直徑為a區(qū)域內(nèi)的所有點的平面和水深一次估值結(jié)果，對待估點進行多重估計，并統(tǒng)計每個估值納入的點數(shù)ξi。

(3) 根據(jù)各估值的百分比εi，尋找最大的εmax。若εmax≥ε′，則εmax對應(yīng)的估值即最優(yōu)水深參考值dr；否則，表明未能找到任何可靠的水深參考值，刪除或標(biāo)記待檢測節(jié)點為不可靠點。

最后，給出如下方法，開展二次水深估計:

(1) 從一次CUBE濾波多重估計結(jié)果中獲取與dr差異最小的水深估值d(m,n,index)

(15)

式中，index為對應(yīng)索引數(shù);count為待估點上的一次CUBE濾波獲取的多重估值個數(shù)。

(2) 然后計算地形梯度變化量ddiff

(16)

(3) 顧及地形梯度ddiff，基于式(17)的原則確定最終水深估值df

(17)

式中，df為最終水深優(yōu)選估值。

獲取可靠水深值的同時，更新該點平面位置。對所有一次CUBE濾波后點進行上述二次濾波處理，最終完成亂石區(qū)測深數(shù)據(jù)的濾波。

3 閾值參數(shù)的確定

3.1 搜索范圍直徑a的確定

如圖4所示，極端情況下，假設(shè)黑色點及中心點全為粗差點，隨著直徑a增大，區(qū)域粗差率降低。理論上，粗差率越低，越易獲得無粗差影響的參考值，但對于復(fù)雜地形，a過大，最終獲得的參考值與待估點的地形相關(guān)性將大大降低，使得待估點水深估值優(yōu)選結(jié)果不可靠。所以參數(shù)a取值時既要考慮粗差率，又要顧及地形起伏。

當(dāng)直徑a取2、4、6倍Δmin，圖4中圓形區(qū)包含的粗差格網(wǎng)點數(shù)占該區(qū)域總點數(shù)(包含區(qū)域邊緣點)的百分比分別為100%、42.9%和20.5%。粗差率一定的情況下，地形變化較大，a應(yīng)取較小值，這樣獲取參考值時納入的節(jié)點與待估點越靠近，參考值與待估點水深值在地形上相關(guān)性更大，最終優(yōu)選得到的水深值越可靠。粗差率通常是未知的，但與地形坡度、水深和測量設(shè)備等相關(guān)。試驗區(qū)大量試驗表明，水深小于50 m，地形坡度大于30°，Δmin=0.5 m，粗差大量聚集時，a取6Δmin較為合適；地形坡度小于30°，a取6～8倍Δmin較為合適。

3.2 有效點數(shù)閾值η給定

η是確保以待估點為中心，直徑為a范圍內(nèi)有足夠多的有效水深值點數(shù)。如果有效點數(shù)太少，獲取的參考水深可靠性將大大降低。地形邊緣或數(shù)據(jù)稀疏地方，有效點數(shù)較少，這種點的可靠性較差(往往表現(xiàn)為破碎地形)，可以將其刪除或者標(biāo)記。有效點數(shù)取值可以根據(jù)直徑為a的搜索區(qū)域的外接正方形包含的點數(shù)得到，即η=ceil[(a/Δmin+1)2/2]，其中ceil函數(shù)為向上取整。

3.3 閾值ε′的確定

閾值ε′用于判斷待檢測節(jié)點水深值和選擇可靠的參考水深估值。如圖4所示，若圖中較大的黑色格網(wǎng)圓節(jié)點及中心點皆為異常點，由3.1節(jié)可得，當(dāng)a=6Δmin時，粗差率為20.5%。如果ε′小于粗差率，可能將待估點水深值判斷為可靠值，從而導(dǎo)致錯誤檢測，因此，ε′必須大于粗差率，才能保證獲得正確的檢測結(jié)果。同時，為了獲得可靠的水深參考值，閾值ε′不宜過大。因此，ε′應(yīng)滿足20.5%<ε′<50%。

4 試驗及分析

為驗證本文方法對亂石區(qū)MBS測深點濾波的有效性，試驗采用以色列海岸防波堤拋石工程施工區(qū)測量數(shù)據(jù)。該數(shù)據(jù)采用Sonic2024多波束測量儀器[21]采集，表2為Sonic2024各項誤差分量[13]。測區(qū)水深1～25 m，最大浪高小于0.25 m，海洋潮汐影響0.1 m左右。邊坡拋石區(qū)石塊大小不一，呈現(xiàn)不規(guī)則分布。根據(jù)施工要求，每次拋石后都要對水下拋石區(qū)域進行測量，準(zhǔn)確掌握每次拋石情況，為下一次拋石提供指導(dǎo)。受亂石區(qū)多次反射波等影響，測量數(shù)據(jù)中存在大量粗差，如圖5中橢圓形虛線區(qū)域所示。

以下采用3種數(shù)據(jù)處理方案開展試驗，分析本文方法的性能。

方案1：采用Caris中傳統(tǒng)CUBE算法濾波。基于CARIS軟件(內(nèi)置CUBE濾波)對測深數(shù)據(jù)處理，結(jié)果如圖6(a)和圖7(a)所示。方案1基于格網(wǎng)數(shù)據(jù)估值，估值結(jié)果仍存在大量粗差，錯誤地反映了邊坡亂石區(qū)地形。

表2 Sonic2024各項誤差分量

方案2：二次CUBE算法濾波。試驗中Δmin=0.5 m、a=3 m、η=25、ε′=25%、μ=4。試驗結(jié)果如圖6(b)、(d)和圖7(b)所示。因同時開展了水深和位置估計，因此估值結(jié)果不再是格網(wǎng)節(jié)點，而是一系列散點，采用三角網(wǎng)構(gòu)建地形模型更為恰當(dāng)(圖6(d))。利用該方案檢測出了大量的異常測深點如圖7(d)所示(圖7(a)、(b)、(c)均為(d)中紅色矩形區(qū)域側(cè)視圖)，由圖7(a)、(d)可知，檢測出的異常點分布與實際吻合，最終濾波結(jié)果粗差剔除徹底。抽取圖5中邊坡亂石區(qū)矩形區(qū)域內(nèi)(圖5(b))的原始測深數(shù)據(jù)，采用方案1和方案2處理，濾波結(jié)果分別如圖6(a)、(c)和(b)、(d)所示。可以看出：

(1) 方案2塊石地形特征更加明顯，石頭間隙更加分明，粗差被徹底消除。

(2) 比較圖6(c)、(d)，方案2得到的估計點在平坦地區(qū)幾乎與方案1得到的格網(wǎng)節(jié)點分布一致，在復(fù)雜地形區(qū)則存在明顯差異，但對亂石區(qū)地形描述更真實。

(3) 方案2顧及了亂石區(qū)測深數(shù)據(jù)特點，同時開展了位置和水深估計，并基于二次濾波消除了粗差影響，因此為測深點賦予了準(zhǔn)確的三維估值，從而也實現(xiàn)了對地形的真實表達。

圖2 邊坡亂石區(qū)測深點分布及對估值的影響Fig.2 Effects of sounding point distribution on sounding estimationin slope riprap area

圖3 亂石區(qū)傳統(tǒng)CUBE濾波結(jié)果Fig.3 Traditional CUBE filtering results in underwater riprap area

圖4 二次CUBE濾波Fig.4 Second CUBE filtering

方案3：CARIS手工編輯。方案3基于CARIS濾波后出現(xiàn)粗差時常采用的一種補救措施，因邊坡亂石區(qū)地形復(fù)雜，有時難以判斷準(zhǔn)確的地形信息，需要有豐富的水深數(shù)據(jù)處理經(jīng)驗。如圖7(c)所示，異常測深點被剔除，基本反映了實際地形變化，但過程非常耗時。以方案3結(jié)果為參考，方案1和方案2濾波獲得的水深值與之比較，差值分布如圖8所示，統(tǒng)計結(jié)果如表3所示。

圖5 邊坡拋石區(qū)地形及大量粗差點云Fig.5 Slope riprap point cloud with a large number of gross errors

圖6 方案1和方案2濾波地形及格網(wǎng)Fig.6 Filter terrain and grid obtained by solution 1 and solution 2

圖7 大量粗差剔除Fig.7 Removing a lots of gross error

Tab.3 Statistics on the estimation error of water depth in slope riprap of different methods m

從圖8和表3可知：

(1) 方案2偏差最小，且偏差分布符合正態(tài)分布。

(2) 方案1偏差均值較大，主要因其機理不足所致，也進一步表明了本文工作的必要性。

(3) 方案1中的問題未在方案2結(jié)果中出現(xiàn)(如圖7(b)和圖8(b)、(c))；方案2實現(xiàn)了自動處理，且將方案1的成果精度提高了近5倍，從而也驗證了本文方法的有效性。

5 結(jié) 論

本文基于數(shù)據(jù)特征提出的二次CUBE濾波算法，給出了位置和水深的同步估計方法、顧及地形梯度的平面和水深不確定度傳遞模型、參考水深的多重估計和優(yōu)選算法并顧及地形梯度的二次水深估計模型，解決了傳統(tǒng)CUBE濾波算法在亂石區(qū)出現(xiàn)的水深估計不準(zhǔn)確、地形特征模糊和粗差剔除能力不足等問題，實現(xiàn)了邊坡亂石區(qū)水深的準(zhǔn)確估計。試驗表明，二次CUBE濾波算法實現(xiàn)了粗差的自動檢測和測深點的準(zhǔn)確估值；在亂石區(qū)，大大提高了傳統(tǒng)CUBE算法水深估值的精度，形成的地形特征更加準(zhǔn)確和明晰；在平坦海床，取得了與傳統(tǒng)CUBE濾波算法近似相同的地形結(jié)果。相對傳統(tǒng)CUBE濾波算法，二次CUBE算法的抗差性更強，數(shù)據(jù)處理效率更高。