張秋明

（福建省港航管理局勘測(cè)中心（福建省港航勘察科技有限公司），福建 福州 350009）

0.引言

多波束測(cè)深系統(tǒng)可以有效地獲得密集的三維空間數(shù)據(jù)，被廣泛應(yīng)用于海洋工程、海底形變監(jiān)測(cè)和海洋資源開發(fā)等領(lǐng)域[1，2]?？紤]到海洋環(huán)境的復(fù)雜性，如魚群、儀器噪聲和可變觀測(cè)條件，MBS不能完全準(zhǔn)確地獲取實(shí)際的海底空間數(shù)據(jù)，在多波束測(cè)量數(shù)據(jù)中存在著異常值，導(dǎo)致其不能充分反映實(shí)際海底地形，給航行安全帶來潛在威脅。為了滿足國(guó)際海道測(cè)量規(guī)范[3]的要求，需要后處理方法以減少這些異常值的影響。

異常值的檢測(cè)是多波束測(cè)深數(shù)據(jù)自動(dòng)處理的關(guān)鍵因素。由于目前的軟件不能有效地解決這一問題，常用的方法是人工編輯法。隨著MBS的不斷發(fā)展，采集數(shù)據(jù)集變得越來越密集，需要花費(fèi)大量時(shí)間的人工編輯法變得更加笨拙[4]。研究人員已經(jīng)提出了許多方法來解決這個(gè)問題。Shaw等[5]提出了一種自回歸算法來檢測(cè)輕微異常值。根據(jù)相鄰點(diǎn)的情況，Du等[6]提出應(yīng)用聚類技術(shù)對(duì)點(diǎn)進(jìn)行分類。Mann等[7]應(yīng)用中值濾波探測(cè)和去除淺水環(huán)境中的異常值?；谪惾~斯分析，Calder等[8]建立了CUBE算法，將穩(wěn)健估計(jì)與實(shí)測(cè)數(shù)據(jù)進(jìn)行比較，根據(jù)兩者存在的差異大小確定異常值。Debese等[9]利用局部魯棒表面估計(jì)和分層自適應(yīng)方法對(duì)不同尺寸的海底網(wǎng)格進(jìn)行分類。Rezvani等[2]應(yīng)用M估計(jì)來衰減異常值的影響。本文結(jié)合將小波變換的定位能力與卡爾曼濾波的魯棒動(dòng)態(tài)預(yù)測(cè)相結(jié)合，提出了一種新的異常值自動(dòng)檢測(cè)算法，并利用數(shù)據(jù)對(duì)其性能進(jìn)行檢測(cè)。

1.基本原理

1.1 離散小波變換

小波變換是一種能同時(shí)進(jìn)行局部分析和多分辨率分析的小波變換。母小波是各種小波變換的基礎(chǔ)如公式（1）所示。通過尺度擴(kuò)展和平移，可以得到一族小波函數(shù)如公式（2）所示：

其中，參數(shù)a是尺度因子，參數(shù)b是平移因子。離散小波變換（DWT）是對(duì)基本小波的尺度和平移的離散化。

DWT顯示出非常適合于信號(hào)，因?yàn)樗臅r(shí)間尺度是由信號(hào)的大小決定的。當(dāng)數(shù)據(jù)在不同頻率下分解成一定的水平時(shí)，就會(huì)識(shí)別出隱藏的特征。

1.2 Kalman濾波

離散卡爾曼濾波的數(shù)學(xué)模型由狀態(tài)方程和觀測(cè)方程兩部分組成。其離散形式可通過公式（3）表示：

其中，φk,k－1表示時(shí)間k-1到k的傳遞矩陣；Γk,k－1表示系統(tǒng)噪聲矩陣；Hk表示k時(shí)刻的系統(tǒng)觀測(cè)矩陣；Vk表示觀測(cè)噪聲；Xk表示狀態(tài)變量的向量；sdk表示系統(tǒng)的觀測(cè)向量矩陣。

2.人工編輯法

現(xiàn)今海洋測(cè)繪工程領(lǐng)域，絕大多數(shù)的工程技術(shù)人員采用的是caris軟件進(jìn)行內(nèi)業(yè)處理。數(shù)據(jù)處理流程主要包括編輯船型配置文件、聲速剖面改正、潮位改正、剔除異常值、CUBE生成、合并數(shù)據(jù)、按比例提取數(shù)據(jù)、圖上展點(diǎn)和繪圖等過程。編輯船型配置文件主要是使用各模塊的相對(duì)位置信息，以及計(jì)算出來的校準(zhǔn)參數(shù)等數(shù)據(jù)構(gòu)造測(cè)船的模型；聲速剖面改正主要是將以波束角和聲波來回傳播時(shí)間格式記錄的多波束測(cè)量原始數(shù)據(jù)轉(zhuǎn)換成：沿航跡、垂直航跡及深度格式的數(shù)據(jù)；剔除“飛點(diǎn)”（也就是異常值）是人工去除錯(cuò)誤“跳點(diǎn)”的過程。每道工序也都需要一定的計(jì)算時(shí)間。特別是剔除異常值需要占據(jù)三分之二的時(shí)間，它需要一個(gè)條帶一個(gè)條帶的人工剔除，非常耗時(shí)。例如2km2面積的多波束掃測(cè)圖，外業(yè)約一天，內(nèi)業(yè)也占據(jù)了一天的時(shí)間，其中剔除異常值需要三分之二的時(shí)間，而單波束測(cè)圖僅需要幾分鐘可出數(shù)據(jù)。

3.組合算法

3.1 構(gòu)建條帶數(shù)據(jù)序列

在本研究中，我們提出了一種空間模式下的異常值檢測(cè)方案。在聲速剖面校正和潮汐校正后，對(duì)MBS觀測(cè)數(shù)據(jù)進(jìn)行變換，得到在局部地理坐標(biāo)系中定義的三維坐標(biāo)。由于所獲得的點(diǎn)云數(shù)據(jù)分布不均勻，不能直接作為后續(xù)處理的輸入數(shù)據(jù)，為了滿足算法的處理要求，將點(diǎn)云數(shù)據(jù)在相同距離中分成幾條，這類似于將表面數(shù)據(jù)“切割”成多條數(shù)據(jù)序列（如圖1所示）。條帶與航行方向垂直。天然海底的形態(tài)幾乎是連續(xù)的，當(dāng)條帶寬度較小時(shí)，單個(gè)條帶中的地形點(diǎn)可以看成是一個(gè)窄的連續(xù)數(shù)據(jù)序列。然而，異常值會(huì)破壞系統(tǒng)的穩(wěn)定性，這是利用條帶數(shù)據(jù)序列檢測(cè)異常值的理論基礎(chǔ)。

圖1 多波束測(cè)深數(shù)據(jù)劃分示意圖

3.2 組合算法處理

當(dāng)穩(wěn)定的數(shù)據(jù)序列只包含孤立的異常值時(shí)，小波分析就能準(zhǔn)確地檢測(cè)出異常值的位置。然而，當(dāng)存在多個(gè)連續(xù)離群點(diǎn)時(shí)，小波檢測(cè)結(jié)果與異常值的位置并不完全一致。換句話說，檢測(cè)到的異常區(qū)域可能同時(shí)包含異常值和真實(shí)探測(cè)數(shù)據(jù)。此時(shí)，我們引入一個(gè)閾值來判斷高頻信號(hào)是否溢出。如果高頻信號(hào)連續(xù)溢出，則可確定為異常區(qū)域。同時(shí)，對(duì)這一異常區(qū)域的范圍進(jìn)行了標(biāo)記，濾波后得到了所有條帶的新的高頻信號(hào)。對(duì)新的高頻信號(hào)和低頻信號(hào)進(jìn)行重構(gòu)，得到新的條帶數(shù)據(jù).它既能消除隨機(jī)噪聲，又能減弱孤立點(diǎn)的影響，大大提高了卡爾曼濾波模型的精度。根據(jù)安裝在船上的不同傳感器的信息，確定卡爾曼濾波的參數(shù)。聯(lián)合算法檢測(cè)和刪除異常值的基本過程。信號(hào)對(duì)新的高頻信號(hào)和低頻信號(hào)進(jìn)行重構(gòu)，得到新的條帶數(shù)據(jù)。它既能消除隨機(jī)噪聲，又能減弱異常值的影響，大大提高了卡爾曼濾波模型的精度。根據(jù)安裝在船上的不同傳感器的信息，確定卡爾曼濾波的參數(shù)。組合算法檢測(cè)和刪除異常值的基本過程（如圖2所示）：

圖2 組合算法探測(cè)和移除異常值流程圖

4.實(shí)測(cè)案例

利用仿真數(shù)據(jù)和實(shí)際數(shù)據(jù)對(duì)組合算法的性能進(jìn)行了測(cè)試。利用仿真數(shù)據(jù)進(jìn)行測(cè)試的目的是驗(yàn)證該算法的有效性。實(shí)際數(shù)據(jù)用于在實(shí)際環(huán)境中測(cè)試應(yīng)用程序。測(cè)量誤差通常有三種：系統(tǒng)誤差、不規(guī)則噪聲和異常值。由于這些數(shù)據(jù)是在早期處理的，因此將消除MBS的系統(tǒng)誤差。不規(guī)則噪聲雖然不可避免且不可預(yù)測(cè)，但可以用正態(tài)分布來表示。因此，這些數(shù)據(jù)的測(cè)試結(jié)果是可靠的。

4.1 模擬數(shù)據(jù)測(cè)試

本文利用SeaBat T20-P寬帶MBS采集的點(diǎn)云數(shù)據(jù)對(duì)算法進(jìn)行了測(cè)試。由于邊緣數(shù)據(jù)質(zhì)量差，只能用中間數(shù)據(jù)進(jìn)行測(cè)試。該地區(qū)地形復(fù)雜，但地形起伏不大。水深在66m-71m之間，數(shù)據(jù)集包含10952個(gè)三維點(diǎn)。

測(cè)量環(huán)境中的異常值可能以兩種形式出現(xiàn)：孤立點(diǎn)或連續(xù)區(qū)域。為了驗(yàn)證該方法的可靠性和穩(wěn)定性，我們還將在這兩種方法中添加異常值。模擬數(shù)據(jù)（如圖3（b）所示）是通過添加異常值（總數(shù)的8%）生成的，范圍從1至3m。這些條帶沿航行方向以一定距離（測(cè)區(qū)的測(cè)深平均距離）構(gòu)造?？偣伯a(chǎn)生120條帶。每個(gè)條帶中的水深被重新排序，以產(chǎn)生新的數(shù)據(jù)序列。然后，創(chuàng)建新的索引，并與原始數(shù)據(jù)建立關(guān)聯(lián)。該區(qū)域原始數(shù)據(jù)經(jīng)小波分解得到的高頻信號(hào)主要是測(cè)深過程中的隨機(jī)誤差，可用于計(jì)算原始均方誤差σ0。截止誤差為σ1=2σ0，由于條帶數(shù)太大，我們選取了六條有代表性的條帶的處理結(jié)果為例。這些條帶共有43個(gè)異常值。利用DB3小波分析了每條帶的高頻信號(hào)和低頻信號(hào)。當(dāng)高頻信號(hào)大于截?cái)嗾`差σ1時(shí)，數(shù)據(jù)被定義為超限。（如圖4所示）除了112條帶中1m的異常值外，所有異常值都已被檢測(cè)到。

圖4為各條帶卡爾曼濾波結(jié)果。濾波值與正常深度基本一致，只是與峰值有明顯偏差。當(dāng)比較偏差和截?cái)嗾`差時(shí)，可以確定超限測(cè)深點(diǎn)為孤立點(diǎn)，并對(duì)其位置和指標(biāo)進(jìn)行標(biāo)注。當(dāng)所有條帶的處理完成后，標(biāo)記的離群點(diǎn)將從計(jì)算機(jī)內(nèi)存中刪除。然后對(duì)所有的條帶進(jìn)行重構(gòu)，得到新的MBS數(shù)據(jù)集。

圖3 第（a）20、（b）55、（c）57、（d）68、（e）93和（f）112條帶數(shù)據(jù)的小波分解和異常區(qū)域的判斷結(jié)果（紅線代表異常區(qū)域）

圖4 第（a）20、（b）55、（c）57、（d）68、（e）93和（f）112條帶組合算法的濾波結(jié)果

4.2 實(shí)測(cè)數(shù)據(jù)案例

測(cè)試案例位于福建省湄洲灣大生島附近海域。該海域水深約為5m-35m，水流流速較大，獨(dú)立漁網(wǎng)較多，水流流速大容易產(chǎn)生噪聲點(diǎn)，漁網(wǎng)多也會(huì)被多波束掃測(cè)到，這些并不是真實(shí)的海底水深數(shù)據(jù)，實(shí)際數(shù)據(jù)集包含約63000個(gè)水深，是在最大開啟角為160°的SeaBat T50P上采集的。由于真實(shí)海底未知，本文從檢測(cè)時(shí)間和探測(cè)位置的符合度兩個(gè)方面對(duì)組合算法與人工編輯過程進(jìn)行了比較。結(jié)果（如表1和圖5所示）：

圖5 實(shí)測(cè)數(shù)據(jù)異常值探測(cè)和移除的結(jié)果比較：（a）組合算法和（b）手工編輯

表1 組合算法和手工編輯在實(shí)際應(yīng)用中比較

5.結(jié)論

該方法不需要迭代運(yùn)算，適用于多核處理。根據(jù)模擬數(shù)據(jù)和實(shí)際數(shù)據(jù)的測(cè)試結(jié)果，該方法具有較強(qiáng)的魯棒性和較高的去除效率。但是，由于測(cè)區(qū)邊界點(diǎn)分布不均勻，難以構(gòu)造連續(xù)的數(shù)據(jù)條帶，無法對(duì)異常點(diǎn)進(jìn)行完全檢測(cè)和剔除。未發(fā)現(xiàn)的異常值主要集中在MBS海底的邊界地區(qū)。此外，但該方法將一些少量的自然尖峰標(biāo)記為異常值。該方法適用于天然水域，地勢(shì)平緩，未施工區(qū)域，若有礁石等區(qū)域會(huì)被誤刪，因此，需要進(jìn)一步優(yōu)化各種參數(shù)，如小波閾值、卡爾曼遞推方程的參數(shù)和條帶寬度。