基于多波束測深聲吶的水中氣體目標(biāo)檢測方法

張萬遠(yuǎn)， 王雪斌， 周天， 張宏偉， 李東洋

(1.哈爾濱工程大學(xué) 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001；2.海洋信息獲取與安全工信部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué))，黑龍江 哈爾濱 150001；3.哈爾濱工程大學(xué) 水聲工程學(xué)院，黑龍江 哈爾濱 150001；4.中海石油深海開發(fā)有限公司 深水工程建設(shè)中心，廣東 深圳 518000；5.天津大學(xué) 機(jī)械工程學(xué)院，天津 300072)

近年來隨著海洋資源開發(fā)和海洋環(huán)境保護(hù)的迫切需要，水中氣體目標(biāo)探測已經(jīng)成為國內(nèi)外學(xué)者重要的研究方向與熱點(diǎn)問題[1-3]。多波束測深聲吶以其搭載平臺多樣化、低可視環(huán)境下適應(yīng)性強(qiáng)等優(yōu)勢，已成為海底目標(biāo)探測、海洋資源勘探和海洋環(huán)境監(jiān)測的重要設(shè)備之一[4]。與光學(xué)探測設(shè)備相比，多波束測深聲吶可以實(shí)現(xiàn)高精度、高分辨率海底大范圍的地形地貌、底質(zhì)分類以及水體目標(biāo)的探測[5-6]，同時(shí)兼具水體目標(biāo)的二維成像能力[7-9]。多波束測深聲吶利用底檢測算法，能夠有效檢測海底地形并估計(jì)出高精度深度值，但對于同時(shí)存在水體目標(biāo)的多回波環(huán)境，常規(guī)的底檢測算法不能滿足復(fù)雜海洋環(huán)境下目標(biāo)探測的需求。針對這一問題，Christoffersen[10]提出波束域多次檢測算法，雖然該算法能夠較好地檢測多目標(biāo)有效回波，但是容易受不均勻背景噪聲和目標(biāo)回波時(shí)寬展比等因素的影響，算法穩(wěn)定性有待提升?；诙嗖ㄊ鴾y深聲吶水體圖像可以提取紋理特征、幾何特征、矩不變特征、延伸率特征、平均回波強(qiáng)度、方差和對比度等特征，進(jìn)行水下靜止目標(biāo)的探測和識別[11-12]。Leighton等[4]提出一種基于水聽器陣列被動探測理論并總結(jié)主被動聲學(xué)監(jiān)測方法的適用性。低頻被動聲學(xué)監(jiān)測系統(tǒng)具有低功耗、定位準(zhǔn)等特點(diǎn)，主動成像聲吶探測效率高、范圍廣。Schneider等[13]將粒子圖像測速法(particle imaging velocimetry, PIV)應(yīng)用在多波束測深聲吶圖像序列，在檢測水中氣泡群的同時(shí)估計(jì)其上升速度，但該方法對于高壓管道泄漏產(chǎn)生的大量密集氣泡群存在圖像序列相關(guān)性低等問題，具有一定的應(yīng)用局限性。

本文提出了一種基于多波束測深聲吶的水中氣體目標(biāo)檢測方法。該方法基于目標(biāo)動態(tài)和靜態(tài)特征，對水中氣體目標(biāo)進(jìn)行綜合檢測。在一維波束域上對氣體目標(biāo)進(jìn)行回波檢測；在二維聲吶單幀圖像上檢測氣體目標(biāo)輪廓和尺度不變特征(scale invariant feature transform flow, SIFT)；在三維聲吶圖像序列中，應(yīng)用SIFT特征流估計(jì)氣體目標(biāo)運(yùn)動特征。在水池條件下和外場水中氣體目標(biāo)檢測實(shí)驗(yàn)中驗(yàn)證該方法的可行性與有效性。

1 水中氣體目標(biāo)綜合檢測

多波束測深聲吶通過接收水下目標(biāo)的回聲信號實(shí)現(xiàn)目標(biāo)檢測與參數(shù)估計(jì)?；芈曅盘柕膹?qiáng)弱與目標(biāo)的聲反射特性密切相關(guān)?？紤]到海洋環(huán)境復(fù)雜，常有魚群、水草、浮游動植物等目標(biāo)干擾，僅從回波強(qiáng)度檢測氣體目標(biāo)是不全面也是不可靠的，本文提出如圖1所示的水中氣體目標(biāo)檢測流程。

圖1 水中氣體目標(biāo)檢測流程Fig.1 Detection flow of underwater gas targets

首先獲取第一組聲吶通道數(shù)據(jù)，對一維波束數(shù)據(jù)進(jìn)行自適應(yīng)閾值多次回波檢測，判斷是否存在水體目標(biāo)。不存在則重新獲取下一組通道數(shù)據(jù)，反之則獲取二維聲吶圖像?？紤]到多波束測深聲吶圖像普遍分辨率低、噪聲嚴(yán)重等特點(diǎn)，采用中值濾波去除聲吶圖像的斑點(diǎn)噪聲，然后應(yīng)用數(shù)學(xué)形態(tài)學(xué)開閉運(yùn)算提取目標(biāo)輪廓特征。在二維聲吶圖像中，水草和繩纜等目標(biāo)同樣具有氣體的形態(tài)等靜態(tài)特征，因此基于多波束測深聲吶圖像序列的氣體目標(biāo)動態(tài)特征檢測是實(shí)現(xiàn)氣體目標(biāo)檢測的一種新穎方法。文中聲吶圖像序列的動特征檢測通過圖像序列特征匹配來實(shí)現(xiàn)。以聲吶圖像全局像素點(diǎn)的SIFT特征向量實(shí)現(xiàn)對應(yīng)點(diǎn)匹配，在匹配過程中采用由粗到精的層次匹配策略和嚴(yán)格的匹配約束，從整體細(xì)化到局部，最后通過求解能量泛函最小化問題來獲取全局最優(yōu)解，實(shí)現(xiàn)全局像素點(diǎn)最優(yōu)匹配。最終估計(jì)目標(biāo)的運(yùn)動特征來判別是否為水中氣體目標(biāo)。

1.1 自適應(yīng)閾值多次檢測

在一維波束域中，本文結(jié)合恒虛警檢測器對多次檢測算法進(jìn)行改進(jìn)，提出一種自適應(yīng)閾值多測檢測算法，提高了算法對復(fù)雜背景環(huán)境的適應(yīng)性與穩(wěn)定性。

1.1.1 多次檢測基本原理

多次檢測算法包括預(yù)檢測和底檢測2部分[10]，預(yù)檢測是在波束形成的基礎(chǔ)上自動計(jì)算檢測閾值，確定有效的目標(biāo)回波區(qū)間，結(jié)合底檢測算法對目標(biāo)回波區(qū)間進(jìn)行處理分析，完成對水體和海底目標(biāo)的檢測。底檢測和預(yù)檢測處理相對獨(dú)立，且國內(nèi)外已經(jīng)進(jìn)行了大量研究，本文不再贅述。圖2是預(yù)檢測處理示意圖。

圖2 預(yù)檢測示意Fig.2 Pre-detection diagram

多次檢測算法基于波束包絡(luò)計(jì)算檢測閾值，從而確定預(yù)檢測區(qū)間：

(1)

其中：

(2)

式中：A為一個(gè)波束的波束幅度；M為波束總采樣點(diǎn)數(shù)；參數(shù)K為比例因子，用于調(diào)整閾值的計(jì)算。該方法基于全局波束數(shù)據(jù)計(jì)算閾值，當(dāng)目標(biāo)回波時(shí)寬占總回波時(shí)寬比例發(fā)生變化或者背景分布不均勻時(shí)，使用該固定閾值不能穩(wěn)定檢測目標(biāo)回波區(qū)間。考慮到多波束測深聲吶不均勻背景噪聲分布和邊緣波束目標(biāo)回波展寬的特點(diǎn)，檢測算法需要一種恒虛警檢測器自適應(yīng)計(jì)算檢測閾值。

1.1.2 可變性指示恒虛警檢測器

恒虛警檢測(constant false-alarm rate, CFAR)是根據(jù)檢測單元的局部背景自適應(yīng)地計(jì)算檢測門限，保持虛警概率恒定，被廣泛應(yīng)用于聲吶信號檢測。單元平均(cell averaging, CA)CFAR檢測器在背景分布均勻的情況下，擁有較強(qiáng)的檢測性能，但其在非均勻背景中檢測性能下降；有序統(tǒng)計(jì)(ordered statistic, OS)CFAR檢測器在已知干擾目標(biāo)數(shù)目時(shí)，其檢測性能優(yōu)于CA-CFAR，但在干擾目標(biāo)數(shù)目未知時(shí)，檢測性能下降；自動刪除平均(automatic censoring cell averaging, ACCA)CFAR自適應(yīng)檢測器通過自動刪除背景估計(jì)單元中的干擾目標(biāo)，提高了檢測器在干擾目標(biāo)未知時(shí)的適應(yīng)性，但其他情況下，有較大的檢測損失；Verma[14]將基于可變性指示(variability index, VI)CFAR檢測器用于聲吶信號檢測，根據(jù)前沿滑窗和后沿滑窗估計(jì)單元的變化性，自適應(yīng)選擇CA-CFAR、最大選擇(greatest of, GO)CFAR和最小選擇(smallest of, SO)CFAR檢測器，體現(xiàn)了檢測器的自適應(yīng)能力。

VI-CFAR檢測器[15]檢測原理如圖3所示，圖中檢測器以平方檢波結(jié)果作為輸入，D為檢測單元，檢測單元兩邊為用于背景噪聲功率估計(jì)的參考滑窗，參考滑窗與檢測單元之間存在長度為lng的保護(hù)單元防止目標(biāo)能量進(jìn)入?yún)⒖蓟?。VI-CFAR通過前沿和后沿滑窗的可變性指數(shù)VI和統(tǒng)計(jì)均值比(mean ratio, MR)分別與閾值KVI和KMR比較，判斷檢測器前沿和后沿滑窗是否均勻和相同，并選擇相應(yīng)背景噪聲估計(jì)單元計(jì)算VI并求和和門限系數(shù)CN計(jì)算自適應(yīng)閾值[14]。如果檢測單元的值超過自適應(yīng)閾值，判斷為目標(biāo)存在，用H1表示；反之，判斷為目標(biāo)不存在，用H2表示。

圖3 VI-CFAR檢測器原理Fig.3 VI-CFAR detector block diagram

1.1.3 水體目標(biāo)檢測

由于水體環(huán)境復(fù)雜，有鰾魚類和浮游動植物與氣泡有相似的回波強(qiáng)度，對水中氣體目標(biāo)檢測造成嚴(yán)重的干擾。本文根據(jù)多次檢測結(jié)果，預(yù)先判斷并提取水體目標(biāo)。如圖4所示，根據(jù)多次檢測結(jié)果在水平方向上設(shè)置寬度為W的水平滑動窗口，然后將滑動窗口沿水平方向從左到右進(jìn)行移動，每次移動距離為W/2。在每次移動中，統(tǒng)計(jì)滑動窗口內(nèi)的測深點(diǎn)數(shù)目，若測深點(diǎn)數(shù)大于設(shè)定閾值，則認(rèn)為該滑動窗口內(nèi)存在水體目標(biāo)，并記錄該滑動窗口位置以及該窗口內(nèi)的最大、最小測深值，

圖4 水體目標(biāo)檢測Fig.4 Water column target detection

1.2 水中氣體目標(biāo)圖像序列檢測

對二維聲吶圖像運(yùn)用數(shù)學(xué)形態(tài)學(xué)處理，去除圖像中目標(biāo)內(nèi)部較小的空洞、平滑目標(biāo)邊緣，最終提取目標(biāo)輪廓特征，進(jìn)一步判斷水中氣體目標(biāo)。為了剔除水草等目標(biāo)的干擾，提取圖像SIFT特征并基于聲吶圖像序列進(jìn)行稠密匹配，估計(jì)目標(biāo)運(yùn)動特征，為氣體目標(biāo)的分布規(guī)模定量評估奠定基礎(chǔ)。

1.2.1 數(shù)學(xué)形態(tài)學(xué)處理

在圖像處理領(lǐng)域，利用數(shù)學(xué)形態(tài)學(xué)方法進(jìn)行圖像處理具有簡化圖像數(shù)據(jù)、保持目標(biāo)的基本形態(tài)特征等優(yōu)點(diǎn)?？紤]到圖像中水中氣體邊緣模糊，內(nèi)部有細(xì)小空洞等特點(diǎn)，運(yùn)用形態(tài)學(xué)開運(yùn)算去除圖像中較小的突出，使目標(biāo)輪廓平滑，斷開圖像中的細(xì)小連接；形態(tài)學(xué)閉運(yùn)算填充圖像中較小的空洞，目標(biāo)輪廓存在的細(xì)小斷裂得以彌合。

形態(tài)學(xué)運(yùn)算定義如下：假設(shè)圖像用f(x,y)表示，結(jié)構(gòu)元素用k(x,y)表示，則膨脹的表達(dá)式為：

(f⊕k)(x,y)=max{f(x-i,y-j)-k(i,j)|

k(i,j)∈K,f(x-i,y-j)∈F}

(3)

腐蝕的表達(dá)式為：

(fΘk)(x,y)=max{f(x-i,y-j)+k(i,j)|

k(i,j)∈K,f(x-i,y-j)∈F}

(4)

對圖像先進(jìn)行腐蝕后膨脹為開運(yùn)算：

F°K=(AΘK)⊕K

(5)

對圖像先進(jìn)行膨脹后腐蝕為閉運(yùn)算：

F·K=(A⊕K)ΘK

(6)

式中：符號⊕表示膨脹運(yùn)算，符號Θ表示腐蝕運(yùn)算，符號°表示開運(yùn)算，符號·表示閉運(yùn)算，F(xiàn)為圖像集合；K為結(jié)構(gòu)元素集合。

聲吶原始圖像如圖5所示，從圖中可以清晰得到在聲吶圖像中直立向上的部分為水體氣體目標(biāo)、一條水平的亮線為海底，而圖中明顯的圓弧由旁瓣泄漏引起。選取合適的結(jié)構(gòu)元素對目標(biāo)圖像進(jìn)行數(shù)學(xué)形態(tài)學(xué)運(yùn)算，再將得到的結(jié)果與原圖像相減得到目標(biāo)的輪廓特征。結(jié)果如圖6所示，圖6(a)為水體目標(biāo)圖像，圖6(b)為經(jīng)過形態(tài)學(xué)處理后得到的水體目標(biāo)輪廓。

圖5 聲吶原始圖像Fig.5 Sonar original image

圖6 水體目標(biāo)、輪廓提取Fig.6 The water column target and its contour

1.2.2 尺度不變特征(SIFT)向量集

Lowe[16]提出SIFT局部特征，在不同尺度空間上查找特征點(diǎn)，計(jì)算出關(guān)鍵點(diǎn)的方向。查找到的關(guān)鍵點(diǎn)是一些十分突出，不會因光照，放射變換和噪聲等因素而變化的點(diǎn)，如角點(diǎn)、邊緣點(diǎn)、亮區(qū)的暗點(diǎn)及暗區(qū)的亮點(diǎn)等。Liu等[17]基于稠密匹配理論提出的尺度不變特征變化流，在光學(xué)圖像配準(zhǔn)和人臉識別領(lǐng)域已有成功的先例。SIFT描述子具有旋轉(zhuǎn)、縮放和平移不變性，對入射角度、噪聲和復(fù)雜環(huán)境多目標(biāo)的干擾保持一定程度的魯棒性。

為了更好地解決密集目標(biāo)的匹配問題，本文僅采用其特征提取部分。對聲吶原始圖像采用中值濾波和動態(tài)亮度分配進(jìn)行預(yù)處理，提高圖像質(zhì)量，將圖像中某一個(gè)像素點(diǎn)鄰域16×16窗口分成4×4單元陣列，每個(gè)單元以45°為間隔量化為8個(gè)方向，從而生成一個(gè)128維特征向量作為該像素的SIFT描述子。遍歷圖像中所有像素點(diǎn)，全部像素點(diǎn)的SIFT描述子構(gòu)成一幀圖像的SIFT特征向量集。

1.2.3 稠密匹配

在連續(xù)多幀聲吶圖像SIFT特征向量集中，本文借鑒光流法中能量泛函的思想[18]，基于SIFT流場光滑和目標(biāo)小位移運(yùn)動2條基本假設(shè)進(jìn)行SIFT特征的匹配。能量函數(shù)為：

min(α|v(p)-v(q)|,d)

(7)

式中：s1和s2為連續(xù)2幀聲吶圖像的SIFT特征向量集。w(p)=(u(p),v(p))為p點(diǎn)流場速度向量，p(x,y)為圖像網(wǎng)格坐標(biāo)；u(p)為水平方向速度矢量；v(p)為垂直方向速度矢量，由連續(xù)2幀聲吶圖像像素點(diǎn)的SIFT描述子匹配計(jì)算得出；ε為p點(diǎn)四鄰域空間；t和d分別為截?cái)郘1范數(shù)的閾值。式中第1項(xiàng)為數(shù)據(jù)項(xiàng)，限定SIFT描述子沿著速度矢量方向進(jìn)行匹配；第2項(xiàng)為小位移項(xiàng)，對于特征比較弱的區(qū)域點(diǎn)，保證其偏移量盡可能??；其他項(xiàng)為平滑項(xiàng)，表示相鄰特征點(diǎn)的速度矢量相近。對于該能量泛函最小化問題，通過循環(huán)信念傳播優(yōu)化算法[19](loopy belief propagation, LBP)進(jìn)行求解，最終實(shí)現(xiàn)聲吶圖像逐點(diǎn)稠密匹配。

1.2.4 由粗到精金字塔匹配策略

由聲吶圖像稠密匹配可知，假設(shè)聲吶圖像像素點(diǎn)數(shù)為L2個(gè)，則LBP算法的時(shí)間復(fù)雜度為O(L4)。為了減少計(jì)算時(shí)間，使用由粗到精的匹配策略進(jìn)行優(yōu)化。如圖7所示在金字塔匹配策略中，將待匹配的連續(xù)2幀聲吶圖像的SIFT特征向量集s1和s2通過高斯函數(shù)進(jìn)行模糊以及降采樣處理得到聲吶圖像金字塔。在金字塔頂層上粗略地估計(jì)SIFT流，即每個(gè)像素點(diǎn)的速度矢量，最后逐步地從粗到細(xì)向下進(jìn)行傳播和細(xì)化。經(jīng)過實(shí)驗(yàn)驗(yàn)證，處理1對256×256圖像，直接的LBP算法處理時(shí)間需要超過2 h，僅式(7)中的數(shù)據(jù)項(xiàng)所需運(yùn)行內(nèi)存已經(jīng)超過16 GB，而由粗到精的金字塔匹配算法的復(fù)雜度為O(L2lbL)，處理時(shí)間為30 s左右，相比于直接的匹配算法運(yùn)算速度取得了顯著的提升。

圖7 由粗到精金字塔匹配示意Fig.7 Coarse-to-fine pyramid matching diagram

2 水池與外場實(shí)驗(yàn)

2.1 水池實(shí)驗(yàn)

為了驗(yàn)證本文所提出方法的有效性，選用實(shí)驗(yàn)設(shè)備為哈爾濱工程大學(xué)研制的多波束測深聲吶。聲吶系統(tǒng)工作頻率為300 kHz，脈沖寬度為35 μs，波束寬度為1°×1°。

在哈爾濱工程大學(xué)水聲技術(shù)國家重點(diǎn)實(shí)驗(yàn)室信道水池對水中氣體目標(biāo)綜合檢測算法進(jìn)行驗(yàn)證。在水池實(shí)驗(yàn)中將多波束測深聲吶傾斜來更好地獲取水中氣體目標(biāo)的整體形態(tài)，實(shí)驗(yàn)裝置布放如圖8所示。

圖8 信道水池實(shí)驗(yàn)示意Fig.8 Channel pool experiment

空氣壓縮機(jī)、氣管和噴嘴等裝置模擬水中氣體目標(biāo)實(shí)驗(yàn)環(huán)境，進(jìn)行多波束測深聲吶水中氣體目標(biāo)綜合檢測實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果如圖9所示，為連續(xù)3幀聲吶圖像的檢測結(jié)果以及SIFT特征流算法估計(jì)出水體目標(biāo)向上的上升速度。在靜水條件下，氣體的上升速度方向垂直向上，伴有左右微弱振蕩。

圖9 信道水池氣體目標(biāo)速度Fig.9 Gas targets velocity in channel pool

2.2 松花湖實(shí)驗(yàn)

為了進(jìn)一步驗(yàn)證算法的穩(wěn)定性與有效性，在外場條件下模擬水中氣體目標(biāo)綜合檢測場景，以實(shí)驗(yàn)船為平臺對水中氣體目標(biāo)進(jìn)行綜合檢測實(shí)驗(yàn)，實(shí)驗(yàn)船和實(shí)驗(yàn)裝置如圖10所示，首先利用三角鐵架將水管及噴嘴固定，然后將其放入湖底，在岸上水管的另一端與空氣壓縮機(jī)連接作為氣體發(fā)生裝置，通過閥門開關(guān)控制壓力大小。實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果如圖11所示，圖11為連續(xù)3幀聲吶圖像的檢測結(jié)果以及SIFT特征流算法估計(jì)出目標(biāo)運(yùn)動特征，在松花湖外場條件下，水中氣體目標(biāo)的在上升的同時(shí)受水流等因素影響產(chǎn)生水平方向的速度矢量。速度方向可分為左上、右上和垂直向上3種情況。

圖10 實(shí)驗(yàn)船和實(shí)驗(yàn)裝置Fig.10 Experiment vessel and device

圖11 松花湖氣體目標(biāo)速度Fig.11 Gas targets velocity in Songhua River

2.3 氣體目標(biāo)速度統(tǒng)計(jì)分析

在外場實(shí)驗(yàn)中，改變氣體的泄漏壓力和噴嘴孔徑來探究泄漏壓力和噴嘴孔徑對泄漏氣體上升速度的影響。選取直徑為0.50 mm的噴嘴，泄漏壓強(qiáng)分別采用0.15、0.30、0.50和0.70 MPa 4個(gè)檔位。圖12(a)給出了不同泄漏壓強(qiáng)下氣體上升速度的盒形圖。從統(tǒng)計(jì)分析結(jié)果可得知，在噴嘴孔徑不變條件下，泄漏氣體的上升速度隨著泄漏壓力的增加而增大。在0.30 MPa泄漏壓強(qiáng)條件下，噴嘴孔徑分別采用0.50，1.00和2.00 mm進(jìn)行3組實(shí)驗(yàn)，統(tǒng)計(jì)分析結(jié)果如圖12(b)所示。壓力恒定時(shí)，泄漏氣體上升速度隨著噴嘴孔徑的增大而增大。

圖12 壓力與噴嘴孔徑對泄漏氣體上升速度的影響Fig.12 The influence of pressure and nozzle aperture on the rising speed of leakage gas

3 結(jié)論

1)在復(fù)雜背景條件下，自適應(yīng)閾值多次檢測算法有效穩(wěn)定地檢測多目標(biāo)回波以及水中氣體目標(biāo)，實(shí)現(xiàn)一維波束域水中氣體目標(biāo)檢測。

2)基于多波束測深聲吶二維圖像，數(shù)學(xué)形態(tài)學(xué)處理可以平滑氣體目標(biāo)邊緣較小的突出，填充目標(biāo)內(nèi)部的空洞。

3)基于圖像金字塔和稠密匹配理論的SIFT流檢測算法可以有效的估計(jì)出水中氣體的運(yùn)動特征。

該方法具有較強(qiáng)的實(shí)用性和重要的工程應(yīng)用性，下一步工作考慮搭載水下移動平臺實(shí)現(xiàn)對近海海域油氣管道泄漏的自動檢測，進(jìn)工程應(yīng)用嘗試。