随着我国海洋开发向远海深海推进，各种海洋活动对海床地形的质量和分辨率也提出了更高的要求^[1-2]。精细的海床地形可表征海床表面“麻坑”或“丘顶”等微观特征^[3-4]，对发现天然气水合物疑似埋藏区域、分析海底矿藏成因机制、研究海床沉积物搬运机理均具有重要意义。目前，高分辨率海底地形主要通过大比例尺单波束测量和多波束测量来获取^[5]。单波束获取高分辨率海底地形，需要密集布线，测量效率低、成本高；多波束系统测量效率高，但随着测量深度增加，其测深点间距增大，边缘地形分辨率快速下降。与上述测深系统不同，侧扫声呐采用拖曳模式贴近海底，可获得到厘米级分辨率海底地貌图像，尽管其不具有直观的地形信息，但去除底质因素的影响，声呐图像的阴影及明暗强度均与地形变化存在明显依赖关系。若能挖掘这些微地形信息，并将其与常规测深数据相结合，即可弥补现有测深系统在分辨率上的缺陷。计算机视觉领域采用明暗修复形状(shape from shading，SFS)方法由图像明暗强度反演目标形状^[6]，该方法基于Lambert法则和一系列的约束条件求解目标在特定尺度下的最优相对高度。为了求解SFS反演地形的绝对尺度并实现其与实测地形的融合，本文采用小波分析方法来寻找侧扫声呐图像反演地形与实际地形的相关性，并以此为基础提出了一种基于单波束水深约束的侧扫声呐图像微地形反演算法。

1 SFS地形反演基本原理

SFS地形反演方法以Lambert漫反射法则为基础，认为入射强度一定的情况下，能量在粗糙表面的漫反射强度仅与能量入射方向、物体表面法线形成夹角θ的余弦有关^[7-9]。物体或海底表面高度可用空间曲面函数z=f(x，y)描述(如图 1)，物体表面法向量N(n_x，n_y，n_z)可表示为高度的梯度形式N(p，q，-1)，根据Lambert法则可列出声能在海床表面的辐照度方程：

$ \begin{array}{l} I\left( {x, y} \right) = R\left( {p\left( {x, y} \right), q\left( {x, y} \right)} \right) = {I_0}\cos \theta = \\ {I_0}\frac{{pp_s + qq_s + 1}}{{\sqrt {{p^2} + {q^2} + 1} \sqrt {p_s^2 + q_s^2 + 1} }} \end{array} $

(1)

式中：I₀为入射能量强度，S(p_s，q_s，-1) 为能量入射方向，式(1) 表征了声能反射强度与地形梯度之间的关系，SFS问题的基本任务就是对给定的图像I，根据每个像素灰度值I(x，y)求定被测物体表面上对应位置的表面相对地形参数，进而重构物体的表明形状z (x，y)。该问题求解还需一系列约束条件，式(2) 为Zheng等^[10]综合了强度约束、强度梯度约束及可积性约束等条件形成的一种经典最小化SFS泛函方程：

$ \begin{gathered} J = \iint {F\left( {p, q, z} \right) = \iint_\mathit{\Omega } {\left[{{{\left( {I-R} \right)}^2}} \right. + \left( {{R_x}-{I_x}} \right)}}2 + \hfill \\ \left. {\left. {{{\left( {{R_y}-{I_y}} \right)}^2} + \mu \left( {{{\left( {{z_x} - p} \right)}^2} + {{\left( {{z_y} - q} \right)}^2}} \right)} \right)} \right]{\text{d}}x{\text{d}}y = \min \hfill \\ \end{gathered} $

(2)

上述泛函极值问题可通过变分法转化为欧拉方程，并采用有限差分方法建立迭代过程，近而求得图像各处的最优地形。大量试验表明，SFS方法在计算机合成图像的三维恢复中可反演出满意的目标形状^[11]，但在真实的光学图像或声学影像三维恢复中，反演结果往往与真实形状有一定的差异。其原因除了漫反射模型的理想化外，声呐发射声波的频率和脉冲宽度一定，其采集的声呐图像明暗强度难以反映海底地形的整体变化趋势，若无真实地形做初始约束时，SFS最小化法仅能解算出特定尺度(起伏范围)和基准(起算基值)下的相对微地形。

2 反演地形与真实地形相关频段提取

若将真实海床地形看成一个全频段非平稳随机信号，那么该信号可经小波变换在空间域分解为多个频段地形，在底质相同区域，侧扫声呐图像的反演地形应与其中某一频段地形具有相关性。确定小波变换的最优分解层数，是寻找与反演地形最相关的真实地形频段，并利用外部水深数据构建反演地形约束模型的关键。

小波变换采用尺度函数$ {\hat \varphi _{a, b}}$和小波函数${{\hat \psi }_{a, b}} $将原始信号分别在尺度空间和小波空间中进行逐层分解，以此实现信号不同频带的分类。在不同频带分离的过程中，尺度函数和小波函数的中心频率随层数不断变化，小波函数和尺度函数的中心频率 f_a^* 与层数的关系如下^[12]

$ f_a^* = \frac{{2{\rm{ \mathsf{ π} }}{f_s}}}{{\left\| {{{\mathit{\boldsymbol{\hat \xi }}}_{a, b}}} \right\|_2^2}}\int_{{f_l}}^{{f_h}} {f{{\left| {{{\mathit{\boldsymbol{\hat \xi }}}_{a, b}}\left( f \right)} \right|}^2}} {\rm{d}}f = \frac{{{f^*}{f_s}}}{a} = \frac{{{f^*}}}{{a\Delta }} $

(3)

式中： f^*为小波母函数的中心频率，f_s为信号采样频率，Δ 为时间采样间隔，a为分解尺度。当$ {{\hat \xi }_{a, b}}\left( f \right) = {{\hat \varphi }_{a, b}}$时，f_a^*为缩放后尺度函数的中心频率；当$ {{\hat \xi }_{a, b}}\left( f \right) = {{\hat \psi }_{a, b}}$时，f_a^* 为缩放后小波函数中心频率。

在真实地形逐层分解过程中，当提取的高频部分中心频率与声呐图像反演地形的主要频率相当时，二者地形起伏变化达到最大相关性。因此，可以先对声呐图像的反演地形进行傅里叶分析，获得其地形频谱分布特征，并选其半功率谱对应频率作为小波变换预提取高频信息的中心频率 f′，进而利用式(4)，可求出需要分解的尺度a:

$ a = {f^*} \cdot {f_s}/f' $

(4)

$ N = {\text{lb}}a = {\text{lb}}\left( {{f^*}{f_s}/f'} \right) $

(5)

根据分解层数与分解尺度的关系按式(5) 可求出最佳的分解层数N。

3 单波束地形与反演地形的约束融合

单波束地形分辨率主要受测线布设宽度限制，但其测线剖面具有较高的采样率，可较好描述海底地形的整体趋势和细节特征。因此，可选取少量单波束测线，采用小波分解方法分别提取与侧扫声呐图像反演地形相关的高频地形和低频地形，利用高频地形为同区域反演地形提供尺度约束，而低频地形作为起算基准，将上述尺度和基准用于无单波束水深的区域，从而得到整个区域的恢复地形。该地形反演约束方法的主要流程如图 2所示。

3.1 尺度约束

忽略海底底质、噪声等因素的影响，侧扫声呐图像的反演地形与单波束高频地形的关系可用线性方程描述，若令单波束采样序列为X，侧扫声呐图像反演地形为f(X)，单波束高频地形为h_hig(X)，g(X)为尺度约束后的反演地形，三者的关系可用下式表示：

$ \begin{gathered} {h_{{\text{hig}}}}\left( X \right)-{{\bar h}_{{\text{hig}}}} = k\left( {f\left( X \right)-\bar f} \right) \hfill \\ g\left( X \right) \approx {h_{{\text{hig}}}}\left( X \right) = kf\left( X \right) + d \hfill \\ \;\;\;\;\;d = {{\bar h}_{{\text{hig}}}}-k\bar f \hfill \\ \end{gathered} $

(6)

式中： k为尺度缩放参数，d为尺度平移参数，认为在一定区域内k和d都为一定值。$ {\bar f}$和$ {{\bar h}_{{\rm{hig}}}}$ 分别为同区域内反演地形和同频段单波束高频地形的期望，式(6) 保证了反演地形经尺度约束后，与从单波束高频地形具有相同的高度期望。

尺度约束的目的是使反演地形与单波束高频地形的起伏程度一致，而地形的标准差是衡量地形起伏程度的重要指标^[13-14]。因此，可利用地形的标准差来求解尺度缩放系数k。两组地形序列在统计规律上应该满足如下关系：

$ \sqrt {\frac{{\sum\limits_{i = 1}^N {{{\left( {h_{{\text{hig}}}^i-{{\bar h}_{{\text{hig}}}}} \right)}^2}} }}{{N-1}}} = k\sqrt {\frac{{\sum\limits_{i = 1}^N {{{\left( {{f^i}-\bar f} \right)}^2}} }}{{N - 1}}} $

(7)

3.2 基准改正

结合式(6)、(7)，可以求出反映真实尺度的反演地形。该地形还需添加正确的起算基准，才能得到最终的反演地形。根据小波分解原理，真实地形H(X)经小波多层分解后，得到目标频段的高频地形h_hig(X)，同时也获得反映真实地形趋势变化的低频地形h_low(X)，该低频地形即为反演地形的起算基准：

$ \left\{ {\begin{array}{*{20}{c}} {H\left( X \right) = {h_{{\text{low}}}}\left( X \right) + {h_{{\text{hig}}}}\left( X \right)} \\ {H'\left( X \right) = {h_{{\text{low}}}}\left( X \right) + g\left( X \right)} \end{array}} \right. $

(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所示。

在与声呐图像重合的区域内共有9条单波束测线，如图 3所示，选取其中5条，即图中的虚线测线，作为内符合测线，用于构建侧扫声呐图像反演地形的约束模型，剩下的4条线，即图中的实线测线，作为外符合检测线。

图 5为侧扫声呐图像经最小化SFS数值算法计算得到的反演地形图，从中可以看出，SFS算法根据图像的明暗变化得到图像区域的相对地形，反演地形的高度在0上下浮动，尺度浮动范围为-1×10^-4~1×10^-4 m，显然该反演地形并不具有正确的起伏尺度和绝对的起算基准。接下来用5条内符合测线对其进行尺度和起算基准修正。

根据5条内符合单波束测线及每条测线上各点的坐标，在反演地形中插值同位置处的反演值，从而得到5对位置对应的实测单波束地形剖面和反演地形剖面序列。分别对每条单波束地形剖面进行小波分解，本文选用bior3.3小波函数，按前所述方法，对各单波束地形分解至第7层。图 6分别列出了5条内符合测线的实际地形剖面、7层小波分解后的低频地形、7层小波分解后的高频地形以及对应位置上的反演地形剖面，图中横坐标为剖面采样序号。

从图 6中可看出，经7层小波分解，从5条内符合单波束测线中均能较好的分离出高、低频地形。从各单波束高频地形及其对应的反演地形起伏变化来看，在同一位置处，反演地形与高频地形具有一致的变化趋势，这说明侧扫声呐图像近似反映了海底地形的高频细节特征，基于声呐图像反演的海底地形与实际地形的高频部分存在较强的线性相似关系。

利用式(6)、(7) 求出尺度缩放参数k及尺度平移参数d，如表 1所示。将上述5条测线求解得到的尺度参数，通过合理的插值方法延拓至整个反演区域。本试验采用线状单波束测深作为外部约束数据，测线与测线之间近似平行分布，测线间隔约为20 m，鉴于外源数据的分布特征，可以根据插值点至相邻测线之间的垂直距离，采用反距离加权法，求解相邻测线之间区域的尺度参数。

根据各位置上求得的尺度参数，对整个反演地形进行尺度修正，尺度约束后的反演地形图如图 7所示。从该图可以看出，经过尺度约束以后，反演地形的浮动范围从-1×10^-4~1×10^-4 m修正至-0.4~0.4 m的正常尺度范围。

接下来，利用从5条单波束水深中提取的低频地形数据，采用带线性插值的三角剖分法构建整个反演区域的起算基准面，如图 8所示，该起算基准面反映了海底地形的整体变化趋势。

将图 7中经过尺度修正后的反演地形与起算基准面相加，得到最终的海底反演地形，如图 9所示。图 10为经过密集单波束精细测量获得的试验区域海底地形图，其分辨率与侧扫声呐图像相同，可以认为该地形是真实的海底的地形，与之相比，图 9中的反演地形不仅具备正确的整体变化趋势，而且保留了侧扫声呐图像丰富的细节特征，反演地形中的细节特征与真实地形的细节特征具有相同的起伏程度和较好的相似性。

为了定量的分析该算法的地形恢复效果，对上述9条测线序列，计算约束反演地形与同位置处单波束实测地形之间的差值，并统计每条测线和所有测线总体的均方差、平均误差、最大误差、最小误差等相关精度参数，内、外符合精度统计情况分别见表 2和表 3。

从内符合精度统计表可以看出(见表 2)，最终的反演地形在5条内符合测线处的均方根误差皆小于15 cm，平均误差均为零，这主要是因为本文采用地形的标准差来计算尺度缩放参数，采用线性变换进行尺度修正，这保证了尺度约束前后反演地形的统计规律不会发生变化。从外符合精度统计表可以看出(见表 3)，在没有单波束水深区域，通过尺度和基准的延拓，最终的反演地形精度与内符合区域相当，均方根误差同样小于15 cm，每条外符合测线的平均误差也较小，达到了厘米级甚至毫米级。图 11和12分别为内、外符合概率分布曲线，从中可以看出，整个反演地形的误差近似服从标准正态分布，基本上不存在系统误差。上述统计指标均满足IHO海道测量规范的要求。

5 结论

1) 试验及统计结果表明，本文所提出的采用少量单波束测深数据对声呐图像的反演地形进行约束改正的思路是可行的、正确的，经过尺度和基准约束以后的反演地形从整体趋势到细节特征均与真实地形保持了较好的一致性和准确性，达到了由侧扫声呐图像恢复高分辨率海底微地形的目的。

2) 该方法仅需要少量外部水深数据，便可利用侧扫声呐图像反演出精度可靠的、高分辨率的海底地形，大大较少了大比例尺单波束测量的成本，也可弥补多波束边缘测深分辨率随水深增加而下降的缺陷。

由于侧扫声呐图像的明暗变化除受地形因素影响外，还与底质、入射角、环境噪声等有关，海底并非是理想的漫反射体，SFS方法用于声呐图像海底地形反演有其局限性，尤其是侧扫声呐近场区域图像并不满足朗伯体法则。因此，基于侧扫声呐图像开展微地形恢复，应避开天底大掠射角区域。采用文中所提方法开展大区域微地形反演时，侧扫声呐数据测线间隔应保证50%以上覆盖率，而单波束测线间隔应根据实际海底地形变化程度合理选择，一般一条侧扫声呐图像应在左右两侧各配合一条单波束测线进行建模。