王海棟，柴洪洲

（1．信息工程大學(xué) 測(cè)繪學(xué)院，河南 鄭州 450052；2．海軍出版社，天津 300450；3．海軍海洋測(cè)繪研究所，天津 300061）

基于CUBE算法的多波束測(cè)深數(shù)據(jù)自動(dòng)處理研究

王海棟1,2，柴洪洲1,3

（1．信息工程大學(xué) 測(cè)繪學(xué)院，河南 鄭州 450052；2．海軍出版社，天津 300450；3．海軍海洋測(cè)繪研究所，天津 300061）

對(duì)CUBE算法自動(dòng)處理多波束測(cè)深數(shù)據(jù)的模型建立、格網(wǎng)節(jié)點(diǎn)的多重估計(jì)和最優(yōu)估值選取準(zhǔn)則進(jìn)行了詳細(xì)介紹，深入分析了多重估計(jì)的實(shí)用性，并通過(guò)實(shí)測(cè)數(shù)據(jù)對(duì)該算法進(jìn)行實(shí)現(xiàn)。利用了抗差Kalman濾波改進(jìn)CUBE算法。通過(guò)模擬數(shù)據(jù)對(duì)改進(jìn)的CUBE算法進(jìn)行實(shí)驗(yàn)，驗(yàn)證了算法改進(jìn)的必要性。

多波束測(cè)深；CUBE算法；抗差估計(jì)；Kalman濾波；異常值檢測(cè)；規(guī)則格網(wǎng)

近年來(lái)，多波束測(cè)深系統(tǒng)測(cè)得的數(shù)據(jù)密度明顯提高，這對(duì)精確繪制海底地形具有十分重要的意義，但同時(shí)也給數(shù)據(jù)處理帶來(lái)了許多問(wèn)題[1]。針對(duì)多波束數(shù)據(jù)異常值檢測(cè)和海底地形格網(wǎng)估計(jì)的復(fù)雜性，美國(guó)新罕布什爾大學(xué)的B.R.Calder和L.A.Mayer[2]提出了一種自動(dòng)處理多波束數(shù)據(jù)的CUBE（Combined Uncertainty and Bathymetry Estimator）算法。該算法通過(guò)建立貝葉斯動(dòng)態(tài)線性模型BDLM（Bayesian Dynamic Linear Model），利用Kalman濾波和多重估計(jì)，并結(jié)合測(cè)深數(shù)據(jù)的水平和垂向不確定度來(lái)計(jì)算格網(wǎng)節(jié)點(diǎn)，具有良好的數(shù)據(jù)處理效果。

實(shí)質(zhì)上，結(jié)合異常值檢測(cè)的海底地形網(wǎng)格化思想國(guó)外已應(yīng)用于多波束測(cè)深數(shù)據(jù)的處理中，而且用地形格網(wǎng)代替測(cè)深值作為水深的存檔數(shù)據(jù)[3]。CUBE算法正是針對(duì)該工序提出的，算法具有高效性、客觀性、穩(wěn)健性及高精度等特點(diǎn)，已集成于許多商業(yè)軟件[4]中，被廣泛應(yīng)用，但國(guó)內(nèi)研究較少。本文主要針對(duì)該算法的理論實(shí)現(xiàn)與改進(jìn)方法進(jìn)行了研究。

CUBE算法通過(guò)計(jì)算每個(gè)測(cè)深點(diǎn)的水平和垂向不確定度，利用信息傳播模型，將測(cè)深信息進(jìn)行傳遞。在所測(cè)區(qū)域，利用估計(jì)節(jié)點(diǎn)及其聯(lián)合精度來(lái)表示該區(qū)域水深。當(dāng)測(cè)深值可用時(shí)，節(jié)點(diǎn)就接收該點(diǎn)信息并估計(jì)水深，進(jìn)一步通過(guò)收集更多的數(shù)據(jù)進(jìn)行及時(shí)更新。

1 CUBE算法模型的建立

為反映各水深值在節(jié)點(diǎn)的估值，綜合考慮節(jié)點(diǎn)和測(cè)深點(diǎn)間距、水平和垂向測(cè)量精度，建立測(cè)深點(diǎn)i對(duì)節(jié)點(diǎn)j的預(yù)報(bào)信息模型為：

由于只預(yù)報(bào)格網(wǎng)節(jié)點(diǎn)的水深，而且格網(wǎng)間距設(shè)置較小，因此可將節(jié)點(diǎn)估值假定為一個(gè)常量，通過(guò)貝葉斯動(dòng)態(tài)模型理論處理序列化的測(cè)深值，考慮最簡(jiǎn)單的常均值模型[7]：

式中，z[n]為節(jié)點(diǎn)水深估值，d[n]為測(cè)深值，為服從方差為的零均值高斯白噪聲，為系統(tǒng)狀態(tài)噪聲，表示模型的可信程度，并假定其方差W[n]為常值0[2]。該模型也是離散線性系統(tǒng)模型[8]的簡(jiǎn)化，作為CUBE算法的Kalman濾波模型。利用該模型計(jì)算，只保留當(dāng)前的水深估值，待新數(shù)據(jù)進(jìn)入后進(jìn)行更新，減小了計(jì)算過(guò)程的存儲(chǔ)。

2 CUBE算法的格網(wǎng)節(jié)點(diǎn)估計(jì)

為保證輸入測(cè)深值的初始數(shù)據(jù)不為異常值，利用中值濾波對(duì)數(shù)據(jù)序列進(jìn)行排序，按照靠近中位數(shù)的情況將輸入數(shù)據(jù)序列的可信值前移，將異常值延后，再利用Kalman濾波對(duì)排序后的數(shù)據(jù)進(jìn)行迭代，計(jì)算節(jié)點(diǎn)估值。

圖 1 測(cè)深點(diǎn)在節(jié)點(diǎn)處的信息預(yù)報(bào)Fig.1 Prediction information from sounding to node

如圖1為測(cè)深點(diǎn)對(duì)格網(wǎng)節(jié)點(diǎn)的信息預(yù)報(bào)示意圖，運(yùn)用Kalman濾波給出了狀態(tài)參數(shù)的遞推公式作為格網(wǎng)節(jié)點(diǎn)估值，將測(cè)深點(diǎn)的信息傳遞給格網(wǎng)節(jié)點(diǎn)，由于算法主要對(duì)格網(wǎng)節(jié)點(diǎn)進(jìn)行估計(jì)，因此計(jì)算過(guò)程簡(jiǎn)化了濾波的狀態(tài)方程。

3 多重估計(jì)及最優(yōu)估值選取準(zhǔn)則

在格網(wǎng)節(jié)點(diǎn)的估計(jì)過(guò)程中，輸入數(shù)據(jù)常含有一定量的測(cè)深異常。若將其和正常值分別進(jìn)行估計(jì)，得到節(jié)點(diǎn)的多個(gè)估值，利用最優(yōu)估值選取準(zhǔn)則對(duì)其作進(jìn)一步判斷，可提高結(jié)果的可信度。

3.1 算法的多重估計(jì)及其實(shí)現(xiàn)

圖 2 多重估計(jì)及最有效模型選取Fig.2 Multiple estimations and available model selection

利用存在跳變的測(cè)深數(shù)據(jù)序列進(jìn)行多重估計(jì)的實(shí)現(xiàn)，結(jié)果如圖3所示，該圖同時(shí)比較了多重估值與單個(gè)估值的計(jì)算結(jié)果。從圖中可見(jiàn)，當(dāng)測(cè)深數(shù)據(jù)存在明顯跳變時(shí)，節(jié)點(diǎn)的多重估值未被跳開(kāi)的測(cè)深數(shù)據(jù)破壞，當(dāng)380個(gè)測(cè)深點(diǎn)附近的連續(xù)異常數(shù)據(jù)消失時(shí)，估計(jì)又繼續(xù)連接第一個(gè)估值，而單個(gè)估值的計(jì)算未正確估計(jì)節(jié)點(diǎn)水深。

圖 3 多重估值與單個(gè)估值的實(shí)驗(yàn)結(jié)果比較Fig.3 Comparison of multiple estimates with single one

由此可見(jiàn)，當(dāng)存在連續(xù)異常的情況下，對(duì)節(jié)點(diǎn)做多重估計(jì)要明顯好于單個(gè)估值的計(jì)算，多重估計(jì)是將不一致的測(cè)深數(shù)據(jù)得到各自內(nèi)部一致的估值，這在一定程度上使算法具有穩(wěn)健性。

3.2 最優(yōu)估值選取準(zhǔn)則

通過(guò)以上的步驟，每一個(gè)節(jié)點(diǎn)都能夠得到相應(yīng)的節(jié)點(diǎn)估計(jì)信息，其包括水深估值、驗(yàn)后方差以及估值的數(shù)量。當(dāng)節(jié)點(diǎn)具有多重估值時(shí)，利用最優(yōu)估值選取準(zhǔn)則來(lái)確定最優(yōu)估值的序號(hào)k，使同時(shí)將其他估值對(duì)應(yīng)的測(cè)深點(diǎn)檢測(cè)為異常值。

最簡(jiǎn)單的準(zhǔn)則是通過(guò)用于估計(jì)的測(cè)深點(diǎn)數(shù)量來(lái)確定k值，公式如下：

一般情況下，上式在選定估值時(shí)速度較快而且準(zhǔn)確性較高，然而，在有連續(xù)脈沖狀異常時(shí)，有可能使異常值的數(shù)量比正常數(shù)據(jù)高，這一準(zhǔn)則可能不適用，顧及海底曲面的連續(xù)性，節(jié)點(diǎn)j的水深和其鄰近節(jié)點(diǎn)應(yīng)接近，因此，可利用局部區(qū)域水深信息的相關(guān)性來(lái)確定k值。

可依據(jù)處理時(shí)間的要求和數(shù)據(jù)的復(fù)雜程度選擇以上3種方法之一進(jìn)行判斷。

格網(wǎng)節(jié)點(diǎn)多重估計(jì)使算法合理地處理了測(cè)深數(shù)據(jù)異常值的問(wèn)題，但也說(shuō)明具有多重估值的節(jié)點(diǎn)周圍存在測(cè)深異常，因此，實(shí)際處理多波束數(shù)據(jù)時(shí)，需對(duì)這些區(qū)域的處理結(jié)果進(jìn)行審核。

4 基于抗差估計(jì)的CUBE算法改進(jìn)

由于每個(gè)多重估計(jì)模型內(nèi)部，CUBE算法假設(shè)測(cè)深數(shù)據(jù)服從正態(tài)分布。但當(dāng)測(cè)深數(shù)據(jù)所含的異常值偏離程度較小，不足以建立多重估計(jì)時(shí)，模型內(nèi)部的數(shù)據(jù)分布受到污染，Kalman濾波估計(jì)節(jié)點(diǎn)的精度下降，雖然程度較小，但若要建立高精度的海底地形，還是不可忽視的。因此，提出利用抗差Kalman濾波對(duì)該情況下的CUBE算法進(jìn)行改進(jìn)。

抗差Kalman濾波常用于對(duì)動(dòng)態(tài)系統(tǒng)含有異常的觀測(cè)值進(jìn)行估計(jì)[10]，來(lái)彌補(bǔ)Kalman濾波對(duì)污染分布數(shù)據(jù)處理的不足。假設(shè)測(cè)深值相互獨(dú)立，并可能含有異常，且服從污染正態(tài)分布：

一步預(yù)測(cè)方差：

利用濾波獲得測(cè)深值殘差，并通過(guò)權(quán)函數(shù)確定其對(duì)應(yīng)的等價(jià)權(quán)，再次對(duì)數(shù)據(jù)進(jìn)行濾波，各測(cè)深值第m步抗差迭代組成的方差陣與等價(jià)權(quán)陣的關(guān)系為，等價(jià)權(quán)值可利用IGGⅢ方案確定[11]，計(jì)算中可將0等價(jià)權(quán)對(duì)應(yīng)的測(cè)深值標(biāo)定為異常。

5 CUBE算法及其改進(jìn)方法的實(shí)現(xiàn)

實(shí)驗(yàn)的測(cè)深數(shù)據(jù)來(lái)自NOAA測(cè)量船在2001年測(cè)得的美國(guó)新罕布什爾州Portsmouth市入?？跀?shù)據(jù)，使用SeaBat 8101多波束測(cè)深儀，已做聲速剖面、潮汐等改正。

圖 4 測(cè)深數(shù)據(jù)分布圖Fig.4 Bathymetry data distribution

取20個(gè)無(wú)異常值的扇區(qū)，由于通常情況下邊緣波束質(zhì)量較差，所以舍棄條帶兩邊各10個(gè)波束，并對(duì)數(shù)據(jù)加入15個(gè)異常值。波束點(diǎn)和異常值的分布如圖4所示，星形點(diǎn)為測(cè)深異常位置，將數(shù)據(jù)進(jìn)行CUBE算法編程實(shí)現(xiàn)。格網(wǎng)間距取0.7 m，共有節(jié)點(diǎn)949個(gè)，計(jì)算結(jié)果如圖5所示。算法在得到海底地形規(guī)則格網(wǎng)圖之前，會(huì)得到格網(wǎng)多重估值個(gè)數(shù)的情況，如圖5(a)。比較圖4和圖5(a)可以發(fā)現(xiàn)，估值個(gè)數(shù)在2個(gè)以上的區(qū)域同時(shí)也存在測(cè)深異常。

圖 5 CUBE算法計(jì)算結(jié)果Fig.5 Results of CUBE algorithm

實(shí)驗(yàn)利用式(16)準(zhǔn)則進(jìn)行最優(yōu)估值的選取，獲得圖5(b)所示的海底地形規(guī)則格網(wǎng)圖，其準(zhǔn)確反映了海底的地形情況。

對(duì)改進(jìn)方法實(shí)驗(yàn)，模擬95個(gè)水深為20.5 m，精度為0.05 m的測(cè)深數(shù)據(jù)，并加入3個(gè)殘差為0.2 m的異常值，用于估計(jì)一個(gè)節(jié)點(diǎn)（圖6(a)）。為說(shuō)明改進(jìn)算法的優(yōu)點(diǎn)，分別通過(guò)以下兩種方案進(jìn)行計(jì)算：

方案1：利用原CUBE算法進(jìn)行多重估值計(jì)算；

方案2：基于抗差Kalman濾波的CUBE算法進(jìn)行節(jié)點(diǎn)估計(jì)。

由于方案2用到IGGⅢ權(quán)函數(shù)，而濾波過(guò)程中，方差為等價(jià)權(quán)的倒數(shù)，因此在實(shí)際計(jì)算過(guò)程中，將權(quán)函數(shù)等于0的值用一小正數(shù)代替。

兩方案的計(jì)算結(jié)果如圖6(a)所示，從圖中可見(jiàn)，所加入的異常值由于偏離程度較小，并未使CUBE算法建立多重估計(jì)。但計(jì)算過(guò)程中，方案1的濾波過(guò)程受到了異常值的影響，致使節(jié)點(diǎn)估值在異常值處有較大的變化，而方案2未受到影響，估計(jì)比較穩(wěn)定，抗干擾能力較好。濾波過(guò)程中兩種方案水深估值的差值如圖6(b)所示，實(shí)際測(cè)量過(guò)程中，由于儀器的測(cè)量精度不同，該差值可能進(jìn)一步放大。因此，CUBE算法的改進(jìn)是有必要的。另外，CUBE算法具體可通過(guò)CARIS軟件實(shí)現(xiàn)，操作過(guò)程可參見(jiàn)文獻(xiàn)[4]和[12]，本文不再贅述。

6 結(jié) 論

本文研究了CUBE算法的模型建立、格網(wǎng)節(jié)點(diǎn)估計(jì)以及多重估計(jì)的思想。該算法有效解決了自動(dòng)處理的可靠性、穩(wěn)健性及海道測(cè)量的安全性等問(wèn)題，利用其估計(jì)的規(guī)則格網(wǎng)容易生成海底三維地形圖，人工處理時(shí)只需少量工作，提高了多波束處理的效率。

圖 6 兩種方案計(jì)算結(jié)果及比較Fig.6 Results and comparison of two schemes

在多個(gè)可信水深估值存在的區(qū)域，算法根據(jù)自定義的準(zhǔn)則或人工方式選擇最佳估值，是自動(dòng)選擇和人工判斷最優(yōu)水深的結(jié)合。

CUBE算法是個(gè)相對(duì)開(kāi)放的算法，本文提出利用抗差Kalman濾波代替該算法中的Kalman濾波，減弱了較小異常對(duì)估值的影響，改進(jìn)了該算法，進(jìn)一步提高了格網(wǎng)節(jié)點(diǎn)估計(jì)的精度。也可以在其他方面改進(jìn)，包括Kalman濾波模型的狀態(tài)方程和觀測(cè)方程，多重估值的選擇準(zhǔn)則等。

[1]李家彪, 王小波, 華祖根, 等.多波束勘測(cè)原理技術(shù)與方法 [M].北京: 海洋出版社, 1999.

[2]Calder B, Mayer L A.Automatic processing of high-rate, highdensity multibeam echosounder data [J].Geochem.Geophys.Geosyst , 2003, 4(6): 24-48.

[3]Armstrong A, Brennan R, Smith S.Implications of the Navigation Surface Approach for Archiving and Charting Shallow Survey Data [C].Shallow Survey 2003.Sydney Australia: 2003.

[4]CARIS.HIPS and SIPS 6.1.1 User’s Guide [M].Fredericton, NB: Universal System Ltdp, 2007.

[5]IHO Standards for Hydrographic Surveys, Special Publication No.44 [S].5th.International Hydrographic Bureau, 2008.

[6]Calder B, Mayer L A.Robust Automatic Multi-Beam Bathymetric Processing [C].U.S.Hydrographic Conference.Norfolk, VA, USA: 2001.

[7]張孝令, 劉福升.貝葉斯動(dòng)態(tài)模型及其預(yù)測(cè) [M].濟(jì)南: 山東科學(xué)技術(shù)出版社, 1992.

[8]楊元喜.自適應(yīng)動(dòng)態(tài)導(dǎo)航定位 [M].北京: 測(cè)繪出版社, 2006.

[9]Vásquez M E.Tuning the CARIS implementation of CUBE for Patagonian Waters [D].New Brunswick, Canada: University of New Brunswick, 2007.

[10]柴洪洲, 崔岳.動(dòng)態(tài)系統(tǒng)的抗差Kalman濾波及其影響函數(shù) [J].中國(guó)慣性技術(shù)學(xué)報(bào), 2002, 10(3): 26-30.

[11]王海棟, 柴洪洲, 翟天增, 等.多波束測(cè)深異常的兩種趨勢(shì)面檢測(cè)算法比較 [J].海洋通報(bào), 2010, 29(2): 182-186.

[12]Mallace D, Robertson P.Alternative Use of CUBE: How to Fita Square Pegin a Round Hole [C].U.S.Hydrographic Conference.Norfolk, Virginia: 2007.

Research on multibeam bathymetry data automatic processing based on CUBE algorithm

WANG Hai-dong1,2，CHAI Hong-zhou1,3

（1.Institute of Surveying and Mapping, Information Engineering University, Zhengzhou 450052, China; 2.Navy Press, Tianjin 300450, China; 3.Naval Institute of Hydrographic Surveying and Charting, Tianjin 300061, China）

The CUBE algorithm in automatic processing multibeam bathymetry data is introduced in detail by its model construction, multiple estimations of grid node and the rules of selecting best available estimation.The efficiency of multiple estimations is analyzed intensively, and the algorithm is implemented using multibeam real data.The Robust Kalman filter is presented to improve the CUBE algorithm.The improvement of the algorithm is compared by the experiment of synthetic data, and its necessity is showed.

multibeam bathymetry; CUBE algorithm; robust estimation; Kalman filter; outlier detection; regular square grid

P229

A

1001-6932(2011)03-0246-06

2010-09-27；收修改稿日期：2011-02-24

中國(guó)博士后科學(xué)基金項(xiàng)目 (20080431342)。

王海棟 ( 1983－ )，男，碩士，主要從事測(cè)量數(shù)據(jù)處理理論與方法研究。

柴洪洲，教授。電子郵箱：chaihz1969@yahoo.com.cn。