劉浩林，覃珊珊，李鵬鴿，張曉暉，劉青

(西安理工大學(xué) 自動(dòng)化與信息工程學(xué)院，陜西 西安 710048)

0 引言

側(cè)掃聲吶效率高及價(jià)格低，尤其具有高分辨率的特點(diǎn)，是海洋開(kāi)發(fā)，海底地形勘察，考古調(diào)查等行業(yè)領(lǐng)域的有力工具[1]。然而，因其工作原理制約，其獲取的聲吶圖像無(wú)法直觀得到被測(cè)區(qū)域高度信息。為解決這一問(wèn)題，國(guó)內(nèi)外學(xué)者在聲吶圖像的三維重構(gòu)方面展開(kāi)了大量研究，并取得了一些成就[2–3]。

明暗恢復(fù)形狀（Shape from shading,SFS）方法為海底地形三維重構(gòu)提供了思路，它可利用圖像的明暗變化來(lái)獲取物體表面的相對(duì)高度[4]。該方法有最小化法、演化法、局部分析法及線性化法4 種。其中最小化法為多數(shù)文獻(xiàn)公認(rèn)魯棒性及抗噪性較好的一類[5–6]。最小化法最開(kāi)始由Horn 等[7]提出，將SFS 問(wèn)題描述為一個(gè)能量函數(shù)，并加入光滑性約束。Zheng[8]并沒(méi)有用光滑性約束條件，而采用圖像梯度約束，從而消除光滑約束的局限性。在此基礎(chǔ)上，Johnson 等[9]為了保證最終的物體表面的相對(duì)高度連續(xù)可積，提出了圖像可積性約束。之后，Leclerc 和Bobick[10]提出了離散化的直接重構(gòu)物體表面相對(duì)高度的迭代算法。然而，現(xiàn)有的SFS 最小化方法是在所研究的對(duì)象為光滑表面的假設(shè)上提出的，對(duì)于實(shí)際的側(cè)掃聲吶圖像的恢復(fù)精度不高。

為此，本文從側(cè)掃聲吶工作原理出發(fā)，首先對(duì)側(cè)掃聲吶原始數(shù)據(jù)進(jìn)行解碼及可視化研究，然后利用小波變換對(duì)SFS 最小化方法進(jìn)行改進(jìn)，對(duì)側(cè)掃聲吶二維瀑布圖進(jìn)行三維重構(gòu)，實(shí)現(xiàn)海底地形反演，并對(duì)反演結(jié)果進(jìn)行分析。

1 側(cè)掃聲吶工作原理

側(cè)掃聲吶工作時(shí)聲吶探測(cè)裝置安裝在拖魚(yú)上，聲吶發(fā)射陣列會(huì)沿垂直于拖魚(yú)前進(jìn)方向的一側(cè)或兩側(cè)發(fā)射一束垂直開(kāi)角很大但水平開(kāi)角很小的聲波脈沖，脈沖遇到海底表面或者水下物體被不斷反射，可按反射信號(hào)的強(qiáng)弱程度畫(huà)出灰度變化不均的聲吶圖像[11]。圖1為側(cè)掃聲吶的工作原理示意圖。

圖1 側(cè)掃聲吶工作原理示意圖Fig.1 Schematic diagram of the working principle of side scan sonar

從聲吶圖像中可以觀察出海底地貌變化，是否有礙航物和海底底質(zhì)類型等信息。對(duì)于凸起、堅(jiān)硬或者粗糙的海底區(qū)域，一般反射回波信號(hào)強(qiáng)；反之，凹陷、細(xì)軟或平緩的海底表面反射回波信號(hào)弱，而被遮擋的海底表面不會(huì)產(chǎn)生反射回波信號(hào)，距離越近的海底反射回波信號(hào)越強(qiáng)[12]。

2 原始數(shù)據(jù)解析

現(xiàn)有的大部分側(cè)掃聲吶原始數(shù)據(jù)處理軟件與聲吶設(shè)備捆綁，也有少數(shù)是國(guó)外的商業(yè)軟件，如Sonar WaveLite 軟件等，但由于軟件側(cè)重的功能不同，缺乏側(cè)掃聲吶全部參數(shù)提取及解析功能，不能滿足海洋地質(zhì)調(diào)查的需求[13–14]。

本文選用C++作為數(shù)據(jù)解析系統(tǒng)的開(kāi)發(fā)語(yǔ)言，利用opencv 作為圖像處理工具，在Visual Studio 2015 中實(shí)現(xiàn)了側(cè)掃聲吶XTF 格式原始數(shù)據(jù)文件的信息提取及可視化，其程序流程圖如圖2 所示。

圖2 XTF 文件數(shù)據(jù)解析流程圖Fig.2 XTF file data analysis flow chart

可視化程序部分利用opencv 視覺(jué)庫(kù)的函數(shù)分別對(duì)通道數(shù)據(jù)中的聲吶回波強(qiáng)度信息進(jìn)行提取，并將其賦值給灰度圖像矩陣，進(jìn)行聲吶瀑布圖顯示，如圖3（a）所示。同時(shí)為了方便對(duì)結(jié)果進(jìn)行驗(yàn)證，在生成圖像時(shí)進(jìn)行反色運(yùn)算，如圖3（b）所示。為驗(yàn)證程序的準(zhǔn)確性，將本文解析軟件結(jié)果與國(guó)外軟件SonarWaveLite 得到的側(cè)掃聲吶瀑布圖進(jìn)行比照，對(duì)比圖如圖4 所示。其中圖左為國(guó)外軟件生成的圖，圖右為程序運(yùn)行結(jié)果圖。

圖3 側(cè)掃聲吶瀑布圖Fig.3 Side scan sonar waterfall chart

圖4 側(cè)掃聲吶瀑布圖對(duì)比圖Fig.4 Side scan sonar waterfall chart comparison

由于國(guó)外軟件會(huì)對(duì)原始瀑布圖進(jìn)行簡(jiǎn)單圖像處理后再顯示，如圖像增強(qiáng)等操作，所以相較于程序運(yùn)行結(jié)果圖對(duì)比度會(huì)較強(qiáng)，但是2 張圖片的基本灰度變化趨勢(shì)是相同的，由此可驗(yàn)證程序的有效性。

3 側(cè)掃聲吶圖像三維重構(gòu)

在側(cè)掃聲吶的三維重構(gòu)部分，本文提出了基于小波變換的SFS 最小化算法，二維圖像在經(jīng)小波分解后的低頻部分采用了基于Priwitt 算子的三維重構(gòu)算法，而在圖像的高頻部分很好地利用了SFS 最小化方法的優(yōu)勢(shì)。實(shí)驗(yàn)表明該方法有效可行，且恢復(fù)結(jié)果的精度得到提高。

3.1 圖像的小波變換

利用小波的的多分辨率特點(diǎn)，可以對(duì)圖片進(jìn)行分解及重構(gòu)，如圖5 所示。其中圖像分解時(shí)，首先對(duì)圖像進(jìn)行行分解，獲得圖像在水平方向上的低頻分量L 和高頻分量H，然后在其基礎(chǔ)上再進(jìn)行列分解，最終獲得一個(gè)低頻分量LL 及LH，HL 和HH 三個(gè)高頻分量圖像。圖像重構(gòu)時(shí)，首先對(duì)分解結(jié)果的每一列進(jìn)行一維離散小波逆變換，再對(duì)所得數(shù)據(jù)的每一行進(jìn)行一維離散小波逆變換，即可獲得重構(gòu)圖像。

圖5 圖像的二維離散小波分解和重構(gòu)過(guò)程Fig.5 Image decomposition and reconstruction process of two-dimensional discrete wavelet

小波變換用于圖像分解時(shí)，最重要的就是小波基的選擇。同一幅圖像，用不同的小波基進(jìn)行分解得到的效果都是不一樣的。本文選擇文獻(xiàn)[14]中提到的具有良好平滑性的小波基db2，來(lái)實(shí)現(xiàn)高質(zhì)量的圖像分解及重構(gòu)。

3.2 基于Priwitt 算子的三維重構(gòu)算法

對(duì)于小波變換后低頻部分圖像的處理，本文提出基于Priwitt 算子的重構(gòu)算法。首先利用Priwitt 算子計(jì)算海底表面每一點(diǎn)的法矢向量，然后利用數(shù)值積分計(jì)算海底表面點(diǎn)的三維高度值，并利用雙線性插值得到海底表面點(diǎn)更精確的高度信息。

3.2.1 表面法式的計(jì)算

海底表面法矢可以通過(guò)表面的傾角和偏角計(jì)算得到。

聲波脈沖入射方向可以用單位向量表示如下：

N(x,y)為像素點(diǎn)l(x,y)處的法向單位向量，其可表示如下：

SFS 方法普遍假設(shè)入射光為無(wú)限遠(yuǎn)處點(diǎn)光源，即S=(0,0,?1)，且物體表面反射模型為朗伯模型[4]，由朗伯模型可知海底某一點(diǎn)l的反射強(qiáng)度值E(x,y)由海底表面反射率 ρ、聲源強(qiáng)度值E0及傾角 θ的余弦值決定，即

海底反射模型幾何關(guān)系如圖6 所示，傾角 θ是法向量N和入射向量S的夾角。

圖6 海底反射模型幾何關(guān)系Fig.6 Geometric relationship of seabed reflection model

式（3）中，ρ及E0一般是常數(shù)，因此可以把E(x,y)歸一化為：

由式（4）可得出傾角的計(jì)算公式為：

從圖6 可看出，偏角 φ為海底表面法線矢量N在xoy平面的投影與x軸的夾角。而瀑布圖像中灰度梯度方向?yàn)楹５妆砻娣ň€投影的方向，所以可通過(guò)灰度梯度來(lái)計(jì)算海底表面法線矢量。

通常，在點(diǎn)(x,y)處的灰度E(x,y)的梯度值為一個(gè)二維列向量，定義如下：

式中，Gx和Gy分別為E(x,y)在x方向和y方向的1 階偏導(dǎo)數(shù)。其中，偏導(dǎo)數(shù)的計(jì)算可以用可以用差分運(yùn)算來(lái)表示。Priwitt 算子就是基于差分運(yùn)算的一個(gè)梯度算子，其水平和垂直梯度算子的模板如圖7 所示。

圖7 Priwitt 算子（8 鄰域3×3）Fig.7 Priwitt operator (8 neighborhood 3×3)

對(duì)于圖中的任一點(diǎn) (i,j)，利用Priwitt 算子對(duì)每一像素點(diǎn)相鄰8 個(gè)像素進(jìn)行加權(quán)運(yùn)算，分別得到該像素點(diǎn)在x方向和y方向的梯度，表達(dá)式為：

在利用式（7）計(jì)算出各個(gè)點(diǎn)的梯度后，根據(jù)式（6）即可以計(jì)算得到偏角 φ的值。

由圖6 中幾何關(guān)系，N(x,y)可表示為：

至此，可以計(jì)算得到所有點(diǎn)的表面法矢向量N(x,y)的值。

3.2.2 數(shù)值積分

在獲得了瀑布圖像的各個(gè)表面點(diǎn)的表面法矢向量N(x,y)的值后，還需根據(jù)N(x,y)計(jì)算海底表面的高度值z(mì)(x,y)。本文利用數(shù)值積分運(yùn)算來(lái)計(jì)算z(x,y)的值。

本文采用牛頓-柯特斯公式的梯度變形公式來(lái)計(jì)算積分。設(shè)[a,b]為一個(gè)有限區(qū)間，有

由于海底表面能夠用含有表面高度信息的方程z=f(x,y)表示，且表面梯度p和q又可以由下式表示：

根據(jù)式（10）可得到z的積分表達(dá)式：

式中，z0表示參考點(diǎn) (x0,y0)處的深度值，結(jié)合式（10）所示的數(shù)值積分，最終計(jì)算出z值。

在一般的離散數(shù)據(jù)計(jì)算中，獲得的數(shù)據(jù)不具有連續(xù)性，因此本文采用雙線性插值的方法，來(lái)獲得除了離散點(diǎn)處函數(shù)值以外的更多數(shù)據(jù)。

3.2.3 基于小波變換的SFS 最小化方法

本文提出的基于小波變換的SFS 最小化方法步驟如下：

1）將二維瀑布圖像經(jīng)小波二次分解后得到1 個(gè)低頻子帶LL 和3 個(gè)高頻子帶。

2）在圖像的低頻部分采用基于Priwitt 算子梯度算法，計(jì)算出這部分物體表面三維高度值。

3）在圖像的3 個(gè)高頻部分采用文獻(xiàn)[10]中經(jīng)典最小化方法迭代求解，計(jì)算出對(duì)應(yīng)點(diǎn)的高度值。

4）將第2 和第3 步得到的實(shí)驗(yàn)結(jié)果結(jié)合起來(lái)，恢復(fù)出物體表面三維形狀。

3.3 實(shí)驗(yàn)結(jié)果比較與分析

利用對(duì)比實(shí)驗(yàn)對(duì)本文提出的基于小波變換的SFS 最小化算法進(jìn)行驗(yàn)證分析。

3.3.1 實(shí)驗(yàn)說(shuō)明

本文關(guān)于三維重構(gòu)算法的評(píng)價(jià)標(biāo)準(zhǔn)采用2 種方式：

1）視覺(jué)效果

以三維重構(gòu)的圖形結(jié)果作為SFS 算法的評(píng)判標(biāo)準(zhǔn)，主要是水下物體的恢復(fù)形狀、海底表面紋路是否清晰及高度起伏是否明顯。

2）量化表示

本文為評(píng)價(jià)SFS 算法在三維重構(gòu)中的效果，將重構(gòu)結(jié)果利用朗伯反射模型投影得到對(duì)應(yīng)的平面灰度圖。再將其與原始輸入圖像利用平均絕對(duì)誤差（MAE）、相關(guān)系數(shù)（CC）及平均結(jié)構(gòu)相似性（MSSIM）3 個(gè)評(píng)價(jià)指標(biāo)進(jìn)行圖像相似性分析。其中，MAE的值越小，表明圖像之間的誤差越小，圖像相似度也越高。相關(guān)系數(shù)的絕對(duì)值|CC|越接近于1 時(shí)，圖像相關(guān)程度越好。MSSIM 越接近1，則2 張圖像越相似。

3.3.2 驗(yàn)證實(shí)驗(yàn)結(jié)果

以Imagenex SportScan 側(cè)掃聲吶實(shí)測(cè)圖像中截取的部分左舷圖像數(shù)據(jù)為例，大小為446×376像素，利用文獻(xiàn)[7]和文獻(xiàn)[10]提出的2 種經(jīng)典SFS 最小化算法及本文提出的算法進(jìn)行對(duì)比實(shí)驗(yàn)，得到三維重構(gòu)結(jié)果如圖8 所示。

分別將三種重構(gòu)結(jié)果進(jìn)行投影，在朗伯反射模型下得到對(duì)應(yīng)的平面灰度圖，如圖9 所示。

將圖9 中三維結(jié)構(gòu)對(duì)應(yīng)的平面灰度圖像與原始輸入圖像進(jìn)行相似度對(duì)比分析，得到的分析結(jié)果如表1所示。

圖9 側(cè)掃聲吶實(shí)地測(cè)試圖像的三維重構(gòu)Fig.9 Three-dimensional reconstruction of field test images of side scan sonar

從圖8 可以觀察到，從視覺(jué)效果來(lái)看，對(duì)于側(cè)掃聲吶實(shí)測(cè)圖像，本文提出的算法恢復(fù)效果最好，三維重構(gòu)圖可清楚顯示海底表面高度波動(dòng)以及水下地形的細(xì)微特征。而文獻(xiàn)[10]的算法效果次之，文獻(xiàn)[7]的算法最差。從表1的評(píng)價(jià)指標(biāo)可看出，本文提出的算法與輸入圖像間的MAE 值最小，表明兩圖像之間的絕對(duì)誤差相對(duì)較小。本文提出算法得到的對(duì)比圖像間CC 及MSSIM 值相比之下是最接近于1的，說(shuō)明圖像間相似度之高，即本文提出算法的三維重構(gòu)結(jié)果精度高。

表1 重構(gòu)結(jié)果對(duì)應(yīng)的平面灰度圖評(píng)價(jià)指標(biāo)對(duì)照表Tab.1 Corresponding table of evaluation indexes of plane gray image corresponding to the reconstruction results

圖8 側(cè)掃聲吶圖像的三維重構(gòu)結(jié)果圖Fig.8 The 3D reconstruction result of the side scan sonar image

4 結(jié)語(yǔ)

本文首先深入了解側(cè)掃聲吶工作原理及作業(yè)模式，在此基礎(chǔ)上，利用Visual Studio 平臺(tái)實(shí)現(xiàn)對(duì)XTF 格式原始數(shù)據(jù)進(jìn)行解析，獲取各種有效信息及回波強(qiáng)度值，利用opencv 工具實(shí)現(xiàn)側(cè)掃聲吶二維瀑布圖的可視化，獲得了理想的結(jié)果。利用小波變換對(duì)SFS 最小化法進(jìn)行改進(jìn)，以二維瀑布圖作為輸入，實(shí)現(xiàn)了海底地形的三維重構(gòu)。最后通過(guò)對(duì)比實(shí)驗(yàn)對(duì)本文提出的算法進(jìn)行分析，結(jié)果表明，與2 種經(jīng)典的SFS 算法相比，基于小波變換的SFS 最小化法性能明顯提高，可清晰顯示海底表面起伏變化和水下地形的細(xì)微特征，從而驗(yàn)證了算法的有效性及準(zhǔn)確性。