引言

海底沉积物的特性，对于军用反潜作战、智能化武器的发展、探测水雷，以及民用上的海洋工程选址、海底资源开发都起到重要作用。其中海底微地形粗糙度是海底沉积物性质中的重要参数。因此，研究海底粗糙度的表达和探测有着重要意义。由于海底沉积物往往具有松散颗粒的结构特点^[1]，打捞后原始结构会遭到影响，所以多采用近距离、原位观测海底微地形粗糙度的技术方法。

海底微地形粗糙度的测量随着上世纪七十年代海底声散射研究的开始而兴起。作为传统的测量方法，手工测量和电导率探针方法属于接触式测量，这种方式存在分辨率低、准确度差的问题。1978年，Akal和Hovem作为早期的研究人员，曾使用立体摄影测量法统计和量化海底微地形^[2]，基于多个不同角度的水下照片对摄像机进行校准，以获得海底三维数字高程模型^[3]。但是这种方法成像时需要校准，实现复杂，并由于海洋复杂的环境中难以实施单一位置多镜头复现的问题，导致该方法在实际应用中难以开展。十九世纪以来，测量海底底质粗糙度参数的方法得到极大发展，2016年，Isakson提出了一种移动式激光测量方法，将激光测量装置安装在ROV(remote operated vehicle)上，可以在无需海底设施的条件下实现大面积测量，并于佛罗里达部署实验^[4]。

作为三维重建方法——SFS算法，成像时不需要校准，且可以通过单张照片中的亮度信息得到海底微地形的数据。SFS算法目前也已经发展了多种理论和应用，并被应用于地形观测、缺陷检测、人体建模等^[5-6]。然而，其水下光学应用相对较少。本文将SFS算法拓展到水下，并通过该方法表示三维海底微地形，提取粗糙度参数。

1   理论模型

1.1   微地形测量

工件表面粗糙度是一个重要的参数^[7-8]，而海底微地形粗糙度对于声纳的性能也具有不可忽视的影响。获取海底微地形的数字表示是计算粗糙度参数的第一步。用函数z(x, y)来描述海底微地形成像结果，x、y对应于二维成像中的横纵坐标，z是像素点(x, y)对应的真实海底点的相对高度，也就是需要求解的变量。具有亮度I(x, y)的图像中的每个点p(x, y, z)对应于真实海底中的点P(X, Y, Z)。考虑到海水对光线具有强吸收、散射作用，依据光学规律建立光线传播模型，如图 1所示。

图  1  水下光传播路径

Figure  1.  Underwater light propagation

下载: 全尺寸图片 幻灯片

根据Beers法则，水下光吸收的规律符合^[9]：

其中:B₀代表发射光强度; d是通过介质传播的距离; B是光线在水中经过d距离后的光强; 参数α是介质中的光吸收系数。通常，吸收系数以cm^-1为单位。考虑到上述海水对光的吸收、散射等物理现象，在假设透视投影的情况下，根据Shaomin Zhang的理论，水下散射规律可写为^[10]

其中：I_W(x, y)表示归一化后每个像素的灰度值；P_S是从成像点P到水下光源的距离；P_C是从成像点P到相机的距离；设$p = \frac{{\partial Z}}{{\partial X}}, q = \frac{{\partial Z}}{{\partial Y}}$, k是坐标和梯度p和q的函数；φ是单位法向量n与光源方向之间的夹角。假设与到成像区域的距离相比，光源和相机方向上的分离可以忽略不计。根据Shaomin Zhang ^[11]，(2)式简化为

其中${P_C} = \sqrt {{{\left( {X - {X_c}} \right)}^2} + {{\left( {Y - {Y_c}} \right)}^2} + Z{{\left( {Z - {Z_c}} \right)}^2}} $。即在Lambert模型I(x, y)=σcosφ基础上，水下光散射模型引入了因子c_f进行修正:

c_f是一个依赖于P_C和α的修正因子。点P(X, Y, Z)处的单位法线可以表示为$\mathit{\boldsymbol{n = }}{\left[ { - p, - q, 1} \right]^{\rm{T}}}/\sqrt {\left( {{p^2} + {q^2} + 1} \right)} $。光源方向表示为$\mathit{\boldsymbol{s = }}{\left[ { - {s_X}, - {s_Y}, 1} \right]^{\rm{T}}}/\sqrt {\left( {s_X^2 + s_X^2 + 1} \right)} $。对于n和s，有^[11]：

代入方程(3)，有：

I_W(x, y)为水下照片里各像素点的归一化亮度，方程(6)是一个含有p、q两个未知数的方程。在不添加约束的情况下，这是一个病态方程，因此需要对方程添加约束才能进行求解。Horn加入了平滑、可积和强度约束，将方程转换为功能极值的解；Zheng和Chellappa用图像梯度代替了平滑的约束条件，用函数极值的方法求解该方程^[12-13]。Tsai引入了牛顿迭代法^[14]，离散化之后，梯度p和q为

用Gr(p, q)表示$\left( {\frac{1}{{P_C^2}}} \right){e^{ - 2\alpha {P_C}}})\cos \gamma $，待求解方程(6)可以表示为

将k设置为迭代次数，对上式在z(x, y)=z^k-1(x, y)处进行泰勒展开，并进行一阶近似，我们可以得到：

一阶展开后，(9)式可由(11)式近似：

其中f_z是z(i, j)的导数：

经过迭代后，进行标定，即可得到微地形的三维数字模型。

1.2   粗糙度求解

经过上一节求解，可以得到水下微地形。测量时水下装置配置如图 2所示，将摄像机放置于待测区域上方的ROV或水下框架上操作。待测参数包括从相机到该区域中心的距离以及光吸收系数。

图  2  水下装置配置示意图

Figure  2.  Underwater installation schematic

下载: 全尺寸图片 幻灯片

成像结果如图 3(a)所示，成像范围为0.4 m×0.4 m，经过三次迭代后，获得的数字高程模型(digital elevation model, DEM)如图 3(b)所示。

图  3  图像及其DEM

Figure  3.  Underwater micro-topography image and its DEM

下载: 全尺寸图片 幻灯片

海底散射的研究与三维地形的具体表征无关，它依赖于频域的统计信息。在图像处理之后，获得界面高度函数z=f(R)，R =(x, y)。海底散射需要空间频率f(周期/单位长度)的统计信息来表示，$f = \frac{1}{\lambda }$。实际上，波数k是更为常用的表征粗糙度谱的自变量。k和f的关系是

海底粗糙度的功率谱强度近似符合幂律形式，因此可以表示为

其中：W(K)是粗糙度谱；K表示二维波长向量；参数ω₂称为谱强度；参数γ₂称为谱指数。用B(R)表示高程数字信息的协方差函数。根据Wiener-Khinchin定理，粗糙度谱W(K)可以由B(R)的傅里叶变换求解得到：

设H(K_x, K_y)为微地形的频谱，z=f(R)在采样数为m/n的x/y方向采样，采样间隔为l_x/l_y，有：

在估算其粗糙度谱之前，为减少频谱泄漏，将二维Hamming窗应用于高度场数据，高度频谱强度等级以dB表示。

图  4  计算得到的功率谱密度

Figure  4.  Cabculated 2-D power spectra density

下载: 全尺寸图片 幻灯片

对海底粗糙度功率谱进行幂律拟合，结果如图 5所示，其中W(K)和K以对数形式进行表示。拟合结果为γ₂=2.581，ω₂=10^-4.768m^4-γ₂，置信度为95%。

图  5  拟合后的功率谱强度

Figure  5.  Roughness power spectra intensity after fitting

下载: 全尺寸图片 幻灯片

2   实验验证

验证实验于2017年在青岛市即墨海域附近部署，其海图位置与海底拍摄情况如图 6所示。

图  6  即墨附近水域海图及典型海底沙坡地形

Figure  6.  Nautical chart near Jimo and seabed ripples diagram

下载: 全尺寸图片 幻灯片

利用三角架架设水下相机，对海底沙坡进行成像，得到图 7(a)中拍摄结果和图 7(b)处理结果。将原始波纹图像与其功率谱强度进行比较，二维功率谱强度在空间频率域内具有整体倾斜的特点，如图 7(b)所示，对应着图 7(a)中沙坡的分布方向。根据图 8中拟合的结果，观察得到粗糙度参数主要是由较高空间频率部分，即微地形的小尺度起伏决定的。对海底微地形的粗糙度功率谱数据进行幂律拟合，γ₂=1.404, ω₂=10^-8.295m^4-γ₂，具有95%的置信度。所得的结果和Briggs在美国东南海岸测量海底粗糙度的结果相当^[15]。如图 9所示，水下SFS算法是一种有效的海底粗糙度测量算法。与Ispkson提出的方法相比，本方法更加便捷、价格低廉。

图  7  海底沙坡及其对应的功率谱强度

Figure  7.  Simulated seabed ripples and its power spectrum intensity

下载: 全尺寸图片 幻灯片

图  8  拟合后的功率谱强度

Figure  8.  Roughness power spectra intensity after fitting

下载: 全尺寸图片 幻灯片

图  9  手工测量和立体相机测量方法

Figure  9.  Manual and stereo camera measurement methods

下载: 全尺寸图片 幻灯片

3   结论

本文提出一种基于SFS的光学海底微地形粗糙度测量方法，结合水下光传播的特点，利用水下SFS得到海底微地形的数字高程模型进行仿真，并在即墨海域部署实验进行验证，证明了该算法有效性。该方法可以通过单张水下照片获得海底微地形特征。对于动态海底而言，这种可安装在ROV或框架上的光学系统非常有效，并且有利于长期和大面积观测。