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

(西安理工大學 自動化與信息工程學院，陜西 西安 710048)

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 圖像的小波變換

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

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

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

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

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

3.2.1 表面法式的計算

海底表面法矢可以通過表面的傾角和偏角計算得到。

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

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

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

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

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

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

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

從圖6 可看出，偏角 φ為海底表面法線矢量N在xoy平面的投影與x軸的夾角。而瀑布圖像中灰度梯度方向為海底表面法線投影的方向，所以可通過灰度梯度來計算海底表面法線矢量。

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

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

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

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

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

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

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

3.2.2 數(shù)值積分

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

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

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

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

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

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

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

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

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

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

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

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

3.3 實驗結(jié)果比較與分析

利用對比實驗對本文提出的基于小波變換的SFS 最小化算法進行驗證分析。

3.3.1 實驗說明

本文關(guān)于三維重構(gòu)算法的評價標準采用2 種方式：

1）視覺效果

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

2）量化表示

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

3.3.2 驗證實驗結(jié)果

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

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

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

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

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

表1 重構(gòu)結(jié)果對應(yīng)的平面灰度圖評價指標對照表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é)語

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