趙祥鴻，暴景陽，黃辰虎，黃賢源，陸秀平

（1.海軍大連艦艇學(xué)院海洋測繪系，遼寧大連116018；2.海軍海洋測繪研究所，天津300061）

基于經(jīng)驗(yàn)?zāi)B(tài)分解削弱多波束殘余誤差的方法

趙祥鴻1，暴景陽1，黃辰虎2，黃賢源2，陸秀平2

（1.海軍大連艦艇學(xué)院海洋測繪系，遼寧大連116018；2.海軍海洋測繪研究所，天津300061）

針對多波束水深數(shù)據(jù)中存在的系統(tǒng)性殘余誤差，提出了基于經(jīng)驗(yàn)?zāi)B(tài)分解方法來削弱殘余誤差的方法：首先利用經(jīng)驗(yàn)?zāi)B(tài)分解方法對多波束測深數(shù)據(jù)作一維分解，將非線性、非平穩(wěn)的多波束測深數(shù)據(jù)分解成準(zhǔn)線性子波，然后構(gòu)建水深數(shù)據(jù)趨勢項(xiàng)與殘余項(xiàng)，利用中央波束趨勢項(xiàng)建立整體數(shù)據(jù)趨勢項(xiàng)，最后加以水深數(shù)據(jù)殘余項(xiàng)還原海底地形，削弱殘余誤差影響。通過實(shí)測多波束測深數(shù)據(jù)驗(yàn)證方法的有效性。

經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）；多波束；殘余誤差；趨勢項(xiàng)；殘余項(xiàng)

多波束測深系統(tǒng)是一個(gè)全覆蓋式聲納測深系統(tǒng)，除自身測量誤差之外，還受聲速變化、定位和測船姿態(tài)、換能器安裝偏差等誤差影響（趙建虎等，2008）。目前國內(nèi)外針對各種誤差影響的改正措施已經(jīng)有很多，趙建虎等（2001）分析了船姿對多波束測量的影響；陽凡林等（2010）分析了多波束測深瞬時(shí)姿態(tài)對水深的影響并提出相應(yīng)的改正措施；孫文川等（2016）提出一種基于最小二乘的多波束換能器橫搖偏差二次校正方法；張志偉（2016）等綜合分析了多波束換能器安裝偏差對海底地形測量的影響。

盡管在進(jìn)行多波束數(shù)據(jù)處理的過程中，針對各項(xiàng)誤差都有相應(yīng)的改正措施，但是處理后的結(jié)果中仍然包含著殘余系統(tǒng)性偏差，表現(xiàn)為相鄰條帶之間重疊區(qū)域內(nèi)測深數(shù)據(jù)出現(xiàn)不匹配的情況。針對此問題，國內(nèi)外專家學(xué)者提出了許多方法：吳自銀等（2005）運(yùn)用GIS面向?qū)ο蠓椒ǚ治隽嗽斐啥嗖ㄊ吘壊ㄊ鴾y量誤差的原因；丁繼勝等（2004）提出了一種誤差強(qiáng)壓制方法來解決相鄰條帶水深不符情況。但是，這些方法致力于解決相鄰條帶測深一致性的問題，沒能從根本上削弱殘余誤差的影響。趙建虎等（2013）提出了一種基于地形變化長波項(xiàng)和短波項(xiàng)相結(jié)合的削弱殘余誤差方法，對于殘余誤差的削弱有很好的效果，但是其殘余誤差削弱原理較為復(fù)雜，實(shí)際操作計(jì)算量較大。

目前小波分析方法廣泛應(yīng)用于處理非平穩(wěn)信號中，但是其在實(shí)際應(yīng)用中存在許多問題，例如小波基函數(shù)的選擇、分解層數(shù)的設(shè)置等，都需要根據(jù)以往的經(jīng)驗(yàn)或者多次實(shí)驗(yàn)來確定，自適應(yīng)性較差。經(jīng)驗(yàn)?zāi)B(tài)分解方法是Huang等（1998）于1998年提出的一種數(shù)據(jù)處理方法，該方法能夠基于數(shù)據(jù)本身的特點(diǎn)，將信號中不同尺度的波動(dòng)或趨勢逐級分解開來，產(chǎn)生一系列本征模態(tài)函數(shù)。其特點(diǎn)是可以處理非線性、非平穩(wěn)過程的數(shù)據(jù)，適合處理多波束測深數(shù)據(jù)，并且該方法自適應(yīng)性較好，任何復(fù)雜的信號都能被分解為有限數(shù)量的簡單分量之和。因此，本文將經(jīng)驗(yàn)?zāi)B(tài)分解方法引入到削弱多波束殘余誤差的工作中，利用經(jīng)驗(yàn)?zāi)７纸夥椒▽Χ嗖ㄊ鴾y深數(shù)據(jù)進(jìn)行一維分解，建立經(jīng)驗(yàn)?zāi)B(tài)函數(shù)與水深數(shù)據(jù)關(guān)系，削弱殘余誤差影響。

1 基于經(jīng)驗(yàn)?zāi)B(tài)分解削弱殘余誤差的方法

1.1 經(jīng)驗(yàn)?zāi)B(tài)分解方法

經(jīng)驗(yàn)?zāi)B(tài)分解方法從本質(zhì)上講是對一個(gè)信號進(jìn)行平穩(wěn)化處理，其結(jié)果是將信號中不同尺度的波動(dòng)或趨勢逐級分解開來，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，每一個(gè)序列稱為一個(gè)本征模態(tài)函數(shù)（IMF）（Echeverria et al，2001；Agarwal et al，2007；Battista et al，2007）。最低頻率的IMF分量通常情況下代表原始信號的趨勢。其大體思路是用波動(dòng)上、下包絡(luò)的平均值去確定“瞬時(shí)平衡位置”，進(jìn)而提取出本征模態(tài)函數(shù)。對數(shù)據(jù)提取若干本征模態(tài)函數(shù)直至滿足終止條件，其處理過程如圖1：

本征模態(tài)函數(shù)需滿足的兩個(gè)條件：（1）在整個(gè)數(shù)據(jù)集中，極值數(shù)和過零數(shù)相等或最多差一個(gè)；（2）在任意點(diǎn)，由局部極大值構(gòu)成的包絡(luò)和由局部極小值構(gòu)成的包絡(luò)均值為零。為了使本征模態(tài)函數(shù)足夠反映實(shí)際的振幅和頻率，設(shè)定篩選停止條件：

循環(huán)終止條件：通常有兩個(gè)判斷標(biāo)準(zhǔn)，一是所得到的本征模態(tài)函數(shù)或者殘余分量小于預(yù)設(shè)的閾值時(shí)；或者當(dāng)殘余分量是單調(diào)函數(shù)時(shí)，不能再篩分出本征模態(tài)函數(shù)。

圖1 經(jīng)驗(yàn)?zāi)B(tài)分解流程

對于原始信號進(jìn)行EMD運(yùn)算，首先就要得到其上下包絡(luò)，即根據(jù)信號極大極小值，用三次樣條曲線進(jìn)行逼近。在進(jìn)行三次樣條插值的時(shí)候，必然使得信號的上下包絡(luò)在信號的兩端附近嚴(yán)重扭曲，即端點(diǎn)飛翼問題。當(dāng)處理低頻或者多分量復(fù)雜信號時(shí)，會(huì)嚴(yán)重淹沒信號的端部特征。為了抑制樣條插值的端點(diǎn)飛翼，又不扭曲原始信號的端部特征，采用端點(diǎn)延拓的方法加以抑制（舒忠平等，2006；孫暉等，2005）。

1.2 殘余誤差削弱方法

多波束水深數(shù)據(jù)X（t）可以表示成由趨勢項(xiàng)Q和殘余項(xiàng)R組成，其中趨勢項(xiàng)包括水深數(shù)據(jù)的主成分，即海底地形的總體趨勢。

對多波束水深數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，將水深數(shù)據(jù)分解為若干經(jīng)驗(yàn)?zāi)B(tài)函數(shù)（IMFi）和殘余量R′。構(gòu)建經(jīng)驗(yàn)?zāi)B(tài)函數(shù)與水深數(shù)據(jù)趨勢項(xiàng)關(guān)系如下：

進(jìn)一步建立經(jīng)驗(yàn)?zāi)B(tài)函數(shù)與水深數(shù)據(jù)的關(guān)系：

表達(dá)式中，n為經(jīng)驗(yàn)?zāi)B(tài)分解的次數(shù)。多波束測深數(shù)據(jù)的殘余誤差主要表現(xiàn)為：兩條相鄰條帶的重疊部分出現(xiàn)水深不符值。多波束測深系統(tǒng)進(jìn)行水深數(shù)據(jù)采集，通常中央波束采集的數(shù)據(jù)質(zhì)量較好。利用兩相鄰條帶中央波束的趨勢項(xiàng)來擬合邊緣波束（重疊部分），加以殘余項(xiàng)還原海底地形，具體步驟如下：

（1）提取多波束水深數(shù)據(jù)中相鄰兩條帶數(shù)據(jù)，以其中一個(gè)條帶的正中波束為基點(diǎn)，采用最小距離法尋找相鄰條帶直線距離最近的點(diǎn)，找出對應(yīng)的另一條帶ping斷面的ping斷面。由于航速和航向等問題，導(dǎo)致相鄰條帶的ping斷面并不嚴(yán)格對應(yīng)，針對這種情況本文選取相鄰的ping斷面與之相對應(yīng)。

（2）對兩條帶之間同一ping斷面上的數(shù)據(jù)X（t）1、X（t）2分別進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，通過上述方法得到兩條帶水深數(shù)據(jù)趨勢項(xiàng)Q1、Q2和殘余項(xiàng)R1、R2。分別截取趨勢項(xiàng)Q1、Q2中的中央波束部分M1、M2和邊緣波束部分L1、L2（中央波束的波束入射角選擇為0°≤φ≤35°，通常認(rèn)為在這個(gè)波束入射角范圍內(nèi)，水深數(shù)據(jù)質(zhì)量最好，可靠性最高），如圖2所示：

圖2 兩條帶水深數(shù)據(jù)趨勢項(xiàng)

（3）對中央波束趨勢項(xiàng)M1、M2擬合出對應(yīng)的邊緣波束的趨勢項(xiàng)L，逐個(gè)循環(huán)各個(gè)ping對，從而獲得兩條帶的總體趨勢項(xiàng)，具體擬合方程如下：

（4）構(gòu)建總體趨勢項(xiàng)Q=Q1+Q2+L，加上殘余項(xiàng)還原海底地形。

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

2.1 實(shí)驗(yàn)準(zhǔn)備

實(shí)驗(yàn)數(shù)據(jù)來源于某海域?qū)崪y數(shù)據(jù)，多波束測深系統(tǒng)每ping包含193個(gè)波束，縱橫波束角均為1.5°，波束開角為140°，采用等角度發(fā)射接收的模式。實(shí)驗(yàn)選取兩條測線各包含500 ping水深數(shù)據(jù)，如圖3所示，可以看出兩條帶重疊部分出現(xiàn)明顯不符情況（為了更好的表現(xiàn)殘余誤差的削弱效果，前期已經(jīng)進(jìn)行了測深異常值的剔除工作）。

圖3 兩條測線原始水深數(shù)據(jù)及對應(yīng)ping斷面

2.2 實(shí)驗(yàn)過程及分析

依據(jù)前文所述方法，分別對兩條帶數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)正交分解，提取本征模態(tài)函數(shù)，構(gòu)建水深數(shù)據(jù)趨勢項(xiàng)和殘余項(xiàng)，具體步驟如下：

（1）先識別出每ping水深數(shù)據(jù)X（t）中的所有極大值點(diǎn)和極小值點(diǎn)，分別利用插值擬合出信號的上包絡(luò)線Xmax（t）和下包絡(luò)線Xmin（t），計(jì)算上下包絡(luò)線的平均值m（t），如圖4所示：

圖4 單ping水深數(shù)據(jù)及擬合的上下包絡(luò)線與平均值

圖5 單ping水深數(shù)據(jù)及所得4個(gè)模態(tài)函數(shù)IMFi和殘余量R′

（2）用原始水深數(shù)據(jù)X（t）減去平均值m（t）得到類距平均值h（t），如果h（t）滿足IMF所規(guī)定的條件，那么h（t）即為分解出的第一個(gè)模式分量，否則，將h（t）作為原始數(shù)據(jù)繼續(xù)重復(fù)上述計(jì)算，直到滿足IMF的條件。

（3）得到第一個(gè)經(jīng)驗(yàn)?zāi)B(tài)函數(shù)，將原始水深數(shù)據(jù)剝離經(jīng)驗(yàn)?zāi)B(tài)函數(shù)作為新的數(shù)據(jù)，重復(fù)步驟（1）和（2），直到滿足循環(huán)結(jié)束條件，得到4個(gè)模態(tài)函數(shù)IMFi以及殘余量R′，如圖5所示：

圖5中，橫坐標(biāo)表示為ping斷面上各個(gè)水深點(diǎn)與最邊緣波束點(diǎn)的距離，縱坐標(biāo)表示水深值。圖5（a）是多波束單ping水深數(shù)據(jù)，圖5（b）-圖2.3（e）是經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解得到的本征模態(tài)函數(shù)，得到的殘余量R′如圖5（f）所示。經(jīng)驗(yàn)?zāi)B(tài)分解將信號中不同尺度的波動(dòng)逐級分解開來，分解產(chǎn)生的本征模態(tài)函數(shù)代表水深數(shù)據(jù)的不同頻率特征，由高頻至低頻依次為IMF1、IMF2、IMF3、IMF4，殘余分量R′中的頻率最低。

分別對選取的兩條測線進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解，將得到的兩組本征模態(tài)函數(shù)（IMFi1、IMFi2）及殘余量R1′、R2′，依據(jù)式（4）構(gòu)建兩條帶水深數(shù)據(jù)趨勢項(xiàng)Q和殘余項(xiàng)R如圖6所示。

圖6 兩條帶水深數(shù)據(jù)趨勢項(xiàng)及殘余項(xiàng)

圖6（a）中藍(lán)色曲線代表兩條帶水深數(shù)據(jù)邊緣波束部分，利用兩條帶中央波束趨勢項(xiàng)擬合邊緣波束趨勢項(xiàng)求得表達(dá)式的系數(shù)，擬合圖如圖7，具體方程表達(dá)式如下：

將重構(gòu)的總體水深數(shù)據(jù)趨勢項(xiàng)加以總體殘余項(xiàng)還原水深數(shù)據(jù)，還原后的單ping水深數(shù)據(jù)效果如圖8，與圖3（b）比較發(fā)現(xiàn)，兩條帶水深數(shù)據(jù)中存在的殘余誤差得到削弱。

圖7 水深數(shù)據(jù)趨勢項(xiàng)

圖8 融合后的單ping對水深數(shù)據(jù)

（4）采用最小距離法尋找條帶之間對應(yīng)的ping斷面。截取相鄰兩條帶測線數(shù)據(jù)，分別提取兩條航跡線正中波束對應(yīng)的位置信息；以其中一條測線的波束位置為基點(diǎn)，在相鄰測線中尋找直線距離最近的點(diǎn)，其對所在的波束即為對應(yīng)的波束。在兩條帶逐個(gè)循環(huán)各個(gè)ping對重復(fù)上述工作，那么整個(gè)水深數(shù)據(jù)的殘余誤差得以削弱。削弱殘余誤差影響后的海底地形，效果如圖9。

比較圖9與圖3，水深數(shù)據(jù)中殘余誤差得到較好的削弱，證明了基于經(jīng)驗(yàn)?zāi)B(tài)分解削弱殘余誤差方法的可行性。量化反映殘余誤差削弱效果如表1所示。

表1中可以看出，經(jīng)過殘余誤差改正后的相鄰條帶不符值的范圍比改正前有所縮??；不符值均值的絕對值和均方誤差值為改正前的一半左右；表明經(jīng)過經(jīng)驗(yàn)?zāi)B(tài)分解方法可以有效構(gòu)建多波束水深數(shù)據(jù)趨勢項(xiàng)及殘余項(xiàng)，削弱殘余誤差的影響。

圖9 削弱殘余誤差影響后的海底地形

表1 相鄰條帶不符值統(tǒng)計(jì)

3 結(jié)論

多波束殘余誤差表現(xiàn)為相鄰條帶重疊水深數(shù)據(jù)的不符，本文利用經(jīng)驗(yàn)?zāi)７纸獾姆椒▽Χ嗖ㄊ鴾y深數(shù)據(jù)進(jìn)行一維分解，將非線性的多波束測深數(shù)據(jù)分解成準(zhǔn)線性的子波。通過提取測深數(shù)據(jù)趨勢項(xiàng)構(gòu)建經(jīng)驗(yàn)?zāi)B(tài)函數(shù)與水深數(shù)據(jù)關(guān)系，并通過實(shí)驗(yàn)驗(yàn)證該方法可有效削弱殘余誤差的影響。

Agarwal,V,Tsoukalas,L.H,2007.Denoising electrical signal via empirical mode decomposition.2007 Irep Symposium－Bulk Power SystemDynamicsandControl－ViiRevitalizingOperational Reliability,1－2:72－77.

Battista.B.M,Knapp,C,McGee,T,et al,2007.Application of the empirical mode decomposition and Hilbert－Huang transform to seismic reflection data.Geophysics,72(2):H29－H37.

Echeverria,J.C.,Crowe,J.A.,Woolfson,M.S,et al,2001.Application of empirical mode decomposition to heart rate variability analysis.Medical&Biological Engineering&Computing,39(4):471－479.

Huang N E,Shen Z,Long S R,et al,1998.The empirical mode decompositionandtheHilbertspectrumfornonlinearand nonstationary time series analysis.Pro.Roy.Soc.Lond,Ser.A454:903－995.

丁繼勝,2004.基于等效聲速剖面法的多波束測深系統(tǒng)聲線折射改正技術(shù).海洋測繪,24(6):27－30.

黃辰虎,2016.聯(lián)合XBT和WOA13模型鹽度信息的深水走航聲速準(zhǔn)確確定.海洋通報(bào),2016(35):554－561.

舒忠平,楊智春,2006.抑制經(jīng)驗(yàn)?zāi)B(tài)分解邊緣效應(yīng)的極值點(diǎn)對稱延拓法.西北工業(yè)大學(xué)學(xué)報(bào),24(5):639－642.

孫暉,2005.經(jīng)驗(yàn)?zāi)B(tài)分解理論與應(yīng)用研究.杭州,浙江大學(xué).

孫文川,暴景陽,金紹華,等,2016.多波束換能器橫搖偏差二次校準(zhǔn)的一種方法.武漢大學(xué)學(xué)報(bào),41(11):1－5.

吳自銀,金祥龍,鄭玉龍,等,2005.多波束測深邊緣波束誤差的綜合校正.海洋學(xué)報(bào),27(4):88－94.

陽凡林,盧秀山,李家彪,等,2010.多波束勘測運(yùn)動(dòng)傳感器偏移的改正方法.武漢大學(xué)學(xué)報(bào),35(7):816－820.

張志偉,暴景陽,肖付民,2016.多波束換能器安裝偏差對海底地形測量的影響.海洋測繪,36(1):51－54.

趙建虎,劉經(jīng)南,2008.多波束測深及圖像數(shù)據(jù)處理.武漢:武漢大學(xué)出版社:101－110.

趙建虎,劉經(jīng)南,張紅梅,2001.船姿分析及其對多波束測量的影響.武漢大學(xué)學(xué)報(bào),26(2):144－149.

趙建虎,張紅梅,嚴(yán)峻,等,2013.削弱殘余誤差對多波束測深綜合影響的方法研究.武漢大學(xué)學(xué)報(bào),38(10):1184－1187.

Weakening effect of residual error based on EMD

ZHAO Xiang-hong1,BAO Jing-yang1,HUANG Chen-hu2,HUANG Xian-yuan2,LU Xiu-ping2
（1.Department of Hydrography and Cartography,Dalian Naval Academy,Dalian 116018,China;2.Naval Institute of Hydrographic Surveying and Charting,Tianjin 300061,China）

This paper puts forward a method for weakening the effect of residual errors of multibeam bathymetric data based on Empirical Mode Decomposition（EMD）.Multibeam bathymetric data is deposed into one-dimensional linear waves by EMD.In that way,the trend and residual of seabed topography is constructed.The trend of whole data is constructed using the center beam trend.Then submarine topography is restored by merging of the data trend and the residual one,which weakens the influence of residual errors.The paper also shows an experiment to prove effectiveness of the method.

Empirical Mode Decomposition（EMD）;multi-beam bathymetric data;residual errors;trend of seabed topography;residual of seabed topography

P733.2；P229

A

1001－6932（2017）06－0662－06

10.11840/j.issn.1001-6392.2017.06.008

2016-09-28；

2017-05-26

國家自然科學(xué)基金（41474012）；國家重大科學(xué)儀器設(shè)備開發(fā)專項(xiàng)（2011YQ12004503）。

趙祥鴻（1989-），博士研究生，主要從事海底地形測量數(shù)據(jù)處理研究。電子郵箱：zhaoxianghong0@163.com。

袁澤軼）