單波束水深約束的側掃聲吶圖像三維微地形反演

王愛學，趙建虎，尚曉東，張紅梅

(1.武漢大學 測繪學院，湖北 武漢430079； 2.東華理工大學 江西省數字國土重點實驗室，江西 南昌 330013； 3.武漢大學 動機學院，湖北 武漢 430072)

?

單波束水深約束的側掃聲吶圖像三維微地形反演

王愛學1，2，趙建虎1，尚曉東1，張紅梅3

(1.武漢大學 測繪學院，湖北 武漢430079； 2.東華理工大學 江西省數字國土重點實驗室，江西 南昌 330013； 3.武漢大學 動機學院，湖北 武漢 430072)

高分辨率海底微地形可呈現更豐富的海床細節特征，對水下目標正確判讀具有重要意義。除去底質因素影響，側掃聲吶圖像明暗變化與海底地形起伏變化有關，挖掘側掃聲吶圖像中的微地形信息，可以彌補現有測深手段分辨率上的不足。針對SFS方法從側掃聲吶圖像反演出的相對地形缺少起伏尺度約束和起算基準等問題，本文通過小波分析方法從單波束數據中分別提取與側掃聲吶圖像反演地形相關的高頻地形及表征地形趨勢變化的低頻地形，以此構建側掃聲吶反演地形的尺度約束及基準改正模型。試驗獲得的反演地形從整體趨勢到細節特征均與真實地形具有較好的一致性，統計結果表明，反演地形內、外符合均方根誤差均小于15 cm，反演地形整體精度與單波束地形相當。

側掃聲吶圖像；微地形反演；明暗形狀恢復；約束方法；地形頻譜分析

隨著我國海洋開發向遠海深海推進，各種海洋活動對海床地形的質量和分辨率也提出了更高的要求[1-2]。精細的海床地形可表征海床表面“麻坑”或“丘頂”等微觀特征[3-4]，對發現天然氣水合物疑似埋藏區域、分析海底礦藏成因機制、研究海床沉積物搬運機理均具有重要意義。目前，高分辨率海底地形主要通過大比例尺單波束測量和多波束測量來獲取[5]。單波束獲取高分辨率海底地形，需要密集布線，測量效率低、成本高；多波束系統測量效率高，但隨著測量深度增加，其測深點間距增大，邊緣地形分辨率快速下降。與上述測深系統不同，側掃聲吶采用拖曳模式貼近海底，可獲得到厘米級分辨率海底地貌圖像，盡管其不具有直觀的地形信息，但去除底質因素的影響，聲吶圖像的陰影及明暗強度均與地形變化存在明顯依賴關系。若能挖掘這些微地形信息，并將其與常規測深數據相結合，即可彌補現有測深系統在分辨率上的缺陷。計算機視覺領域采用明暗修復形狀(shape from shading，SFS)方法由圖像明暗強度反演目標形狀[6]，該方法基于Lambert法則和一系列的約束條件求解目標在特定尺度下的最優相對高度。為了求解SFS反演地形的絕對尺度并實現其與實測地形的融合，本文采用小波分析方法來尋找側掃聲吶圖像反演地形與實際地形的相關性，并以此為基礎提出了一種基于單波束水深約束的側掃聲吶圖像微地形反演算法。

1 SFS地形反演基本原理

SFS地形反演方法以Lambert漫反射法則為基礎，認為入射強度一定的情況下，能量在粗糙表面的漫反射強度僅與能量入射方向、物體表面法線形成夾角θ的余弦有關[7-9]。物體或海底表面高度可用空間曲面函數z=f(x，y)描述(如圖1)，物體表面法向量N(nx，ny，nz)可表示為高度的梯度形式N(p，q， -1)，根據Lambert法則可列出聲能在海床表面的輻照度方程：

I(x，y) =R(p(x，y)，q(x，y))=I0cosθ=

(1)

式中：I0為入射能量強度，S(ps，qs， -1)為能量入射方向，式(1)表征了聲能反射強度與地形梯度之間的關系，SFS問題的基本任務就是對給定的圖像I，根據每個像素灰度值I(x，y)求定被測物體表面上對應位置的表面相對地形參數，進而重構物體的表明形狀z(x，y)。該問題求解還需一系列約束條件，式(2)為Zheng等[10]綜合了強度約束、強度梯度約束及可積性約束等條件形成的一種經典最小化SFS泛函方程：

J=?F(p，q，z)=?Ω[(I-R)2+(Rx-Ix)2+(Ry-Iy)2+μ((zx-p)2+(zy-q)2))]dxdy=min

(2)

上述泛函極值問題可通過變分法轉化為歐拉方程，并采用有限差分方法建立迭代過程，近而求得圖像各處的最優地形。大量試驗表明，SFS方法在計算機合成圖像的三維恢復中可反演出滿意的目標形狀[11]，但在真實的光學圖像或聲學影像三維恢復中，反演結果往往與真實形狀有一定的差異。其原因除了漫反射模型的理想化外，聲吶發射聲波的頻率和脈沖寬度一定，其采集的聲吶圖像明暗強度難以反映海底地形的整體變化趨勢，若無真實地形做初始約束時，SFS最小化法僅能解算出特定尺度(起伏范圍)和基準(起算基值)下的相對微地形。

圖1 海床表面形狀模型Fig.1 Sea floor surface model

2 反演地形與真實地形相關頻段提取

若將真實海床地形看成一個全頻段非平穩隨機信號，那么該信號可經小波變換在空間域分解為多個頻段地形，在底質相同區域，側掃聲吶圖像的反演地形應與其中某一頻段地形具有相關性。確定小波變換的最優分解層數，是尋找與反演地形最相關的真實地形頻段，并利用外部水深數據構建反演地形約束模型的關鍵。

(3)

在真實地形逐層分解過程中，當提取的高頻部分中心頻率與聲吶圖像反演地形的主要頻率相當時，二者地形起伏變化達到最大相關性。因此，可以先對聲吶圖像的反演地形進行傅里葉分析，獲得其地形頻譜分布特征，并選其半功率譜對應頻率作為小波變換預提取高頻信息的中心頻率f′，進而利用式(4)，可求出需要分解的尺度a:

a=f*·fs/f′

(4)

N=lba=lb(f*fs/f′)

(5)

根據分解層數與分解尺度的關系按式(5)可求出最佳的分解層數N。

3 單波束地形與反演地形的約束融合

單波束地形分辨率主要受測線布設寬度限制，但其測線剖面具有較高的采樣率，可較好描述海底地形的整體趨勢和細節特征。因此，可選取少量單波束測線，采用小波分解方法分別提取與側掃聲吶圖像反演地形相關的高頻地形和低頻地形，利用高頻地形為同區域反演地形提供尺度約束，而低頻地形作為起算基準，將上述尺度和基準用于無單波束水深的區域，從而得到整個區域的恢復地形。該地形反演約束方法的主要流程如圖2所示。

圖2 單波束地形約束的SFS地形反演方法流程圖Fig.2 Flow chart about SFS terrain recovering method constrained by single-beam soundings

3.1 尺度約束

忽略海底底質、噪聲等因素的影響，側掃聲吶圖像的反演地形與單波束高頻地形的關系可用線性方程描述，若令單波束采樣序列為X，側掃聲吶圖像反演地形為f(X)，單波束高頻地形為hhig(X)，g(X)為尺度約束后的反演地形，三者的關系可用下式表示：

g(X)≈hhig(X)=kf(X)+d

(6)

尺度約束的目的是使反演地形與單波束高頻地形的起伏程度一致，而地形的標準差是衡量地形起伏程度的重要指標[13-14]。因此，可利用地形的標準差來求解尺度縮放系數k。兩組地形序列在統計規律上應該滿足如下關系：

(7)

3.2 基準改正

結合式(6)、(7)，可以求出反映真實尺度的反演地形。該地形還需添加正確的起算基準，才能得到最終的反演地形。根據小波分解原理，真實地形H(X)經小波多層分解后，得到目標頻段的高頻地形hhig(X)，同時也獲得反映真實地形趨勢變化的低頻地形hlow(X)，該低頻地形即為反演地形的起算基準：

(8)

整個約束反演的過程可用式(8)表示，H′(X)為經過尺度和基準約束后的最終反演地形。上述尺度與基準求解方法適用于與單波束測線同區域的側掃聲吶圖像地形反演，在沒有單波束測線的區域，則需要采用插值方法，如移動曲面擬合、反距離加權、三角剖分等，對前文求得的尺度參數和起算基準進行插值或延拓。

4 微地形反演試驗

本文利用實測的單波束水深數據和同區域側掃聲吶圖像設計了外部水深約束的SFS地形反演方法試驗。選用深圳以南某水域為試驗區，區域內海底底質相近，地形特征變化明顯；采用單波束測量對該區域海底地形進行精細測量，單波束采樣間隔為0.5 m，測線間隔為10 m，圖3為試驗區單波束測量地形圖；采用Edgetech 4200側掃聲吶采集同區域海底表面的聲吶圖像，圖4為一塊150 m×85 m的側掃聲吶圖像，圖像分辨率為300×170，像素尺寸為0.5 m，該聲吶圖像已經過必要的畸變改正、灰度均衡等處理，聲吶圖像在試驗區中范圍如圖3所示。

圖3 所選區域海底地形及測線分布情況Fig.3 Seabed topography of test site and layout of survey lines

圖4 實測側掃聲吶圖像Fig.4 Side-scan sonar image measured in test area

在與聲吶圖像重合的區域內共有9條單波束測線，如圖3所示，選取其中5條，即圖中的虛線測線，作為內符合測線，用于構建側掃聲吶圖像反演地形的約束模型，剩下的4條線，即圖中的實線測線，作為外符合檢測線。

圖5為側掃聲吶圖像經最小化SFS數值算法計算得到的反演地形圖，從中可以看出，SFS算法根據圖像的明暗變化得到圖像區域的相對地形，反演地形的高度在0上下浮動，尺度浮動范圍為-1×10-4～1×10-4m，顯然該反演地形并不具有正確的起伏尺度和絕對的起算基準。接下來用5條內符合測線對其進行尺度和起算基準修正。

圖5 側掃聲吶圖像的SFS反演地形圖Fig.5 SFS recovering topography from side-scan sonar image

根據5條內符合單波束測線及每條測線上各點的坐標，在反演地形中插值同位置處的反演值，從而得到5對位置對應的實測單波束地形剖面和反演地形剖面序列。分別對每條單波束地形剖面進行小波分解，本文選用bior3.3小波函數，按前所述方法，對各單波束地形分解至第7層。圖6分別列出了5條內符合測線的實際地形剖面、7層小波分解后的低頻地形、7層小波分解后的高頻地形以及對應位置上的反演地形剖面，圖中橫坐標為剖面采樣序號。

從圖6中可看出，經7層小波分解，從5條內符合單波束測線中均能較好的分離出高、低頻地形。從各單波束高頻地形及其對應的反演地形起伏變化來看，在同一位置處，反演地形與高頻地形具有一致的變化趨勢，這說明側掃聲吶圖像近似反映了海底地形的高頻細節特征，基于聲吶圖像反演的海底地形與實際地形的高頻部分存在較強的線性相似關系。

利用式(6)、(7)求出尺度縮放參數k及尺度平移參數d，如表1所示。將上述5條測線求解得到的尺度參數，通過合理的插值方法延拓至整個反演區域。本試驗采用線狀單波束測深作為外部約束數據，測線與測線之間近似平行分布，測線間隔約為20 m，鑒于外源數據的分布特征，可以根據插值點至相鄰測線之間的垂直距離，采用反距離加權法，求解相鄰測線之間區域的尺度參數。

表1 各測線的尺度參數

根據各位置上求得的尺度參數，對整個反演地形進行尺度修正，尺度約束后的反演地形圖如圖7所示。從該圖可以看出，經過尺度約束以后，反演地形的浮動范圍從-1×10-4～1×10-4m修正至-0.4～0.4 m的正常尺度范圍。

接下來，利用從5條單波束水深中提取的低頻地形數據，采用帶線性插值的三角剖分法構建整個反演區域的起算基準面，如圖8所示，該起算基準面反映了海底地形的整體變化趨勢。

將圖7中經過尺度修正后的反演地形與起算基準面相加，得到最終的海底反演地形，如圖9所示。圖10為經過密集單波束精細測量獲得的試驗區域海底地形圖，其分辨率與側掃聲吶圖像相同，可以認為該地形是真實的海底的地形，與之相比，圖9中的反演地形不僅具備正確的整體變化趨勢，而且保留了側掃聲吶圖像豐富的細節特征，反演地形中的細節特征與真實地形的細節特征具有相同的起伏程度和較好的相似性。

為了定量的分析該算法的地形恢復效果，對上述9條測線序列，計算約束反演地形與同位置處單波束實測地形之間的差值，并統計每條測線和所有測線總體的均方差、平均誤差、最大誤差、最小誤差等相關精度參數，內、外符合精度統計情況分別見表2和表3。

圖6 單波束測線小波分解及與反演地形剖面的比較Fig.6 Wavelet decomposition of single-beam soundings and comparison with recovering terrain

圖7 尺度約束后的反演地形圖Fig.7 Recovering topography constrained by scale correction

圖8 低頻地形生成的起算基準面Fig.8 Reference surface generated by low frequency terrain

圖9 尺度和基準約束后的最終恢復地形圖Fig.9 Final recovering topography with scale and reference correction

圖10 同區域單波束精細測量地形圖Fig.10 Elaborate terrain surveyed by single-beam sonar

表3 外符合精度統計表

從內符合精度統計表可以看出(見表2)，最終的反演地形在5條內符合測線處的均方根誤差皆小于15 cm，平均誤差均為零，這主要是因為本文采用地形的標準差來計算尺度縮放參數，采用線性變換進行尺度修正，這保證了尺度約束前后反演地形的統計規律不會發生變化。從外符合精度統計表可以看出(見表3)，在沒有單波束水深區域，通過尺度和基準的延拓，最終的反演地形精度與內符合區域相當，均方根誤差同樣小于15 cm，每條外符合測線的平均誤差也較小，達到了厘米級甚至毫米級。圖11和12分別為內、外符合概率分布曲線，從中可以看出，整個反演地形的誤差近似服從標準正態分布，基本上不存在系統誤差。上述統計指標均滿足IHO海道測量規范的要求。

圖11 內符合概率密度曲線Fig.11 Internal probability density curve

圖12 外符合概率密度曲線Fig.12 External probability density curve

5 結論

1)試驗及統計結果表明，本文所提出的采用少量單波束測深數據對聲吶圖像的反演地形進行約束改正的思路是可行的、正確的，經過尺度和基準約束以后的反演地形從整體趨勢到細節特征均與真實地形保持了較好的一致性和準確性，達到了由側掃聲吶圖像恢復高分辨率海底微地形的目的。

2)該方法僅需要少量外部水深數據，便可利用側掃聲吶圖像反演出精度可靠的、高分辨率的海底地形，大大較少了大比例尺單波束測量的成本，也可彌補多波束邊緣測深分辨率隨水深增加而下降的缺陷。

由于側掃聲吶圖像的明暗變化除受地形因素影響外，還與底質、入射角、環境噪聲等有關，海底并非是理想的漫反射體，SFS方法用于聲吶圖像海底地形反演有其局限性，尤其是側掃聲吶近場區域圖像并不滿足朗伯體法則。因此，基于側掃聲吶圖像開展微地形恢復，應避開天底大掠射角區域。采用文中所提方法開展大區域微地形反演時，側掃聲吶數據測線間隔應保證50%以上覆蓋率，而單波束測線間隔應根據實際海底地形變化程度合理選擇，一般一條側掃聲吶圖像應在左右兩側各配合一條單波束測線進行建模。

[1]金翔龍. 海洋地球物理技術的發展[J]. 東華理工學院學報， 2004， 27(1)： 6-13.

JIN Xianglong. The development of technique of marine geophysics[J]. Journal of East China Institute of Technology， 2004， 27(1)： 6-13.

[2]陸秀平， 黃謨濤， 翟國君， 等.多波束測深數據處理關鍵技術研究進展與展望[J]. 海洋測繪， 2016， 36(4)： 1-6.

LU Xiuping，HUANG Motao， ZHAI Guo jun， et al. Development and prospect of key technologies formultibeam echosounding data processing[J]. Hydrographic surveying and charting， 2016， 36(4)： 1-6.

[3]王健， 邱文弦， 趙俐紅. 天然氣水合物發育的構造背景分析[J]. 地質科技情報， 2010， 29(2)： 100-106.

WANG Jian， QIU Wenxian， ZHAO Lihong. Tectonic settings analysis of gas hydrate deposits development[J]. Geological science and technology information，2010， 29(2)： 100-106.

[4]鄧希光， 吳廬山， 付少英，等. 南海北部天然氣水合物研究進展[J]. 海洋學研究， 2008， 26(2)： 67-74.

DENG Xiguang， WU Lushan， FU Shaoying， et al. The research advances of natural gas hydrates in northern south China Sea[J]. Journal of marine science， 2008， 26(2)： 67-74.[5]董慶亮， 歐陽永忠， 陳岳英，等. 側掃聲納和多波束測深系統組合探測海底目標[J]. 海洋測繪， 2009， 29(5)： 51-53.DONG Qinliang， OUYANG Yongzhong， CHEN Yueying， et al. Measuringbottom of sea target with side-scansonarandmultibeam soundingsystem[J]. Hydrographic surveying and charting，2009， 29(5)： 51-53.

[6]HORNB K. Obtaining shape from shading information[M]. Massachusetts: MIT press， 1989.

[7]HORN B K. Shape from shading: a method for obtaining the shape of a smooth opaque object from one view[R]. US-MIT， 1970， 232.

[8]HORNBK.Height and gradient from shading[J]. International journal of computer vision， 1990， 5(1)： 37-75.

[9]HORNBK，SZELISKIRS， YUILLEAL. Impossible shaded images[J]. IEEE transactions on pattern analysis and machine intelligence， 1993， 15(2)： 166-170.

[10]ZHENG Qinfen， CHELLAPPA R.Estimation of illuminant direction， albedo， and shape from shading[J]. IEEE transactions on pattern analysis & machine intelligence， 1991: 13(7)： 680-702.

[11]ZHANG Ruo， TSAI Pingsing， CRYER J E， et al. Shapefromshading: a survey[J]. Pattern analysis and machine intelligence， 1999， 21(8)： 690-706.

[12]樊計昌， 劉明軍， 海燕，等. 計算尺度函數和小波函數中心頻率的 GUI 及其應用[J]. 科技導報， 2008， 25(24)： 36-39.

FAN Jichang， LIU Mingjun， HAI Yan， et al. GUI for computing center-frequency of the scale and waveletfunctionsandits applications[J].Science & technology review，2008， 25(24)： 36-39.

[13]隋剛， 郝兵元， 彭林. 利用高程標準差表達地形起伏程度的數據分析[J]. 太原理工大學學報， 2010， 41(4)： 381-384.

SUI Gang， HAO Bingyuan， PENG Lin. Data analysis of elevation standard deviation classifying relief degree of land surface [J]. Journal of Taiyuan University of Technology， 2010， 41(4)： 381-384.

[14]劉曉萌， 常占強. 數字高程模型數據分布規律的研究[J]. 河北師范大學學報: 自然科學版， 2006， 30(4)： 473-477.

LIU Xiaomeng， CHANG Zhanqiang. Research on the distribution law of DEM data [J]. Journal of Hebei Normal University: natural science edition，2006， 30(4)： 473-477.

本文引用格式：

王愛學， 趙建虎， 尚曉東， 等. 單波束水深約束的側掃聲吶圖像微地形反演[J]. 哈爾濱工程大學學報， 2017， 38(5)： 739-745.

WANG Aixue，ZHAO Jianhu，SHANG Xiaodong，et al. Recovery of seabed 3D micro-topography from side-scan sonar image constrained by single-beam soundings [J]. Journal of Harbin Engineering University， 2017， 38(5)： 739-745.

Recovery of seabed 3D micro-topography from side-scan sonar image constrained by single-beam soundings

WANG Aixue1，2，ZHAO Jianhu1，SHANG Xiaodong1，ZHANG Hongmei3

(1.School of geodesy and Geomatics， Wuhan University， Wuhan 430079， China; 2.Jiangxi Province Key Lab for Digital Land， East China Institute of Technology， Nanchang 330013， China; 3.School of Power and Mechanical Engineering， Wuhan 430072， China)

Without considering the impact of marine sediment， the gray intensity of a side-scan sonar image mainly depends on the changing seabed micro-topography， which can be used to compensate for the low-resolution deficiency of traditional bathymetry methods. In order to recover the micro-topography and transform it into a valid scale and absolute reference， we used the wavelet analysis method to extract high- and low-frequency terrains from single-beam soundings， in which the high-frequency terrains correspond with topography recovered from side-scan sonar images and the low-frequency terrains represent topography change. On this basis， we constructed a size restraint and reference correction model for side-scan sonar topography recovery. The test results show that the final recovered topography is consistent with the actual terrain in overall trend and local detail. Statistical results show that the internal and external root mean square error (RMSE) are both less than 15 cm， and the final recovery accuracy is the same as that for single-beam soundings.

side-scan sonar image; micro-topography recovery; shape from shading; constrain method; spectrum analysis of topography

2015-12-17.

日期：2017-04-26.

國家自然科學基金項目(41376109，41576107，41606114)；江西數字國土重點實驗室開放基金項目(DLLJ201502)；精密工程與工業測量國家測繪局重點實驗室開放基金項目(PF2015-7).

王愛學(1984-)， 男， 講師， 博士.

王愛學，E-mail: axwang@sgg.whu.edu.cn.

10.11990/jheu.201512057

P229.1

A

1006-7043(2017)05-0739-07

網絡出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170426.1031.010.html