基于SFS 方法的側掃聲吶圖像三維重構

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

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

0 引言

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

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

為此，本文從側掃聲吶工作原理出發，首先對側掃聲吶原始數據進行解碼及可視化研究，然后利用小波變換對SFS 最小化方法進行改進，對側掃聲吶二維瀑布圖進行三維重構，實現海底地形反演，并對反演結果進行分析。

1 側掃聲吶工作原理

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

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

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

2 原始數據解析

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

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

圖2 XTF 文件數據解析流程圖Fig.2 XTF file data analysis flow chart

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

圖3 側掃聲吶瀑布圖Fig.3 Side scan sonar waterfall chart

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

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

3 側掃聲吶圖像三維重構

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

3.1 圖像的小波變換

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

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

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

3.2 基于Priwitt 算子的三維重構算法

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

3.2.1 表面法式的計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2.2 數值積分

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

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

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

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

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

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

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

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

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

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

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

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

3.3 實驗結果比較與分析

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

3.3.1 實驗說明

本文關于三維重構算法的評價標準采用2 種方式：

1）視覺效果

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

2）量化表示

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

3.3.2 驗證實驗結果

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

分別將三種重構結果進行投影，在朗伯反射模型下得到對應的平面灰度圖，如圖9 所示。

將圖9 中三維結構對應的平面灰度圖像與原始輸入圖像進行相似度對比分析，得到的分析結果如表1所示。

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

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

表1 重構結果對應的平面灰度圖評價指標對照表Tab.1 Corresponding table of evaluation indexes of plane gray image corresponding to the reconstruction results

圖8 側掃聲吶圖像的三維重構結果圖Fig.8 The 3D reconstruction result of the side scan sonar image

4 結語

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