单波束水深约束的侧扫声呐图像三维微地形反演

王爱学，赵建虎，尚晓东，张红梅

(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