机载激光雷达测深系统的研发开始于20世纪60年代末期，其中以澳大利亚研制成功的WRELADS-1为代表，主要进行了激光测深的机理研究，到了20世纪90年代瑞典研制的HAWE EYE系统标志着机载激光雷达测绘近海海底地貌进入实用化^[1]。机载激光雷达集成全球卫星定位系统、惯性导航系统、扫描系统和激光测距系统，通过各个子系统得到位置、姿态和距离信息来解算地物的空间位置^[2]。机载海洋激光雷达在海洋地形测绘中相对于传统的人工测量加船载测量方式，可以克服对于浅水地区船只无法通行的困难，具有机动性高、测绘面积大、成本低等优势^[3]。因此机载海洋激光雷达测深在海洋浅海深度提取、海岸线调查、军事侦查等调查具有重要地位^[4]。

机载海洋激光雷达测深原理是从机载激光全回波波形中提取出海面回波和海底回波时刻，通过光速乘以时间计算出海水深度，因此如何正确提取海底回波时刻是机载激光雷达的关键技术。为了提升机载海洋激光雷达的最大可探测深度，国产机载海洋激光雷达采用敏感度强的PMT装置，用于放大强度低的深水回波信号，但同时也放大了噪声。因此对于深水的回波波形，海底的微弱信号在回波波形上幅度小，而且容易被噪声淹没。常见的海底回波提取算法包括波峰检测、解卷积、波形分解^[5]，且算法主要包括两大步骤，初始参数确定和参数优化^[6]。由于这些方法都对于单条波形进行处理，参数的确定易受到噪声的影响，进而影响到下一步参数优化，因此容易提取出错误的海底回波时刻，甚至提取不到海底回波。为此需要一种能够有效从噪声中提取出海底回波的处理方法。

综上，由于噪声在空间上不具有相关性，而海底回波在空间上具有相关性，文中将波形按采集顺序组合成二维的回波强度图像，通过图像处理方法，可以有效区分海底回波信号和噪声部分，进而提取出海底回波廓线，提升海底回波的提取能力。

实验采用的机载海洋激光雷达设备是中国科学院上海光学精密机械研究所研制的新型机载双频雷达测深系统（Mapper 5000），如图1所示，实验区域为中国海南省三亚市南海某海岛附近海域。

图  1  机载双频激光雷达Mapper 5000系统

Figure 1.  Airborne dual-frequency lidar Mapper 5000 system

该系统具有海洋测绘和陆地测绘两种模块，海洋测绘模块使用1064 nm波长激光和532 nm波长激光，其中1064 nm波长激光探测器为APD，该回波信息用于记录海表回波信号，532 nm波长激光探测器为PMT，该回波信息用于提取海底回波信息。海洋测绘模块采用卵形扫描器和高速波形采集卡，可实现最大30°角，采样间隔1 ns的探测采集^[7]。为了能够采集到不同海水深度的海底回波信号，海洋测绘模块采用了分视场的接收光路设计，如图2所示，总的视场角为40 mrad，其中中间的6 mrad视场作为浅水通道，6~40 mrad的环形视场作为深水通道。对于深水通道削去了中间视场，可以避免海面和浅水区域的强回波信号，而随着激光光束随海水深度逐渐展宽，深水的海底回波信号将被环形视场接收。文中选取实验数据为新型机载双频雷达测深系统的深水通道数据。

图  2  分视场接收光路示意图

Figure 2.  Schematic diagram of receiving light path for sub-field of view

新型机载双频雷达测深系统的全回波波形会因为激光器能量的抖动或是探测器响应饱和使得回波波形能量基线值不同，如图3(a)所示，因此为了增大信号和噪声的对比度，避免在图像上因为回波基线值不同导致后续图像分割错误的情况，需要对单一回波进行基值校正。文中选取每条波形最后100 ns范围内的回波均值作为基线值，将回波减去基线值，完成基值校正，如图3(b)所示。

图  3  未校正基值波形之间对比和校正基值后波形之间对比。（a）未校正基值波形之间对比；（b）校正基值波形之间对比

Figure 3.  Comparisons of uncorrected base waveform and corrected base waveform. (a) Uncorrected base value waveform comparison; (b) Corrected base value waveform comparison

为了能够拼接出连贯的海底回波波形廓线，选取新型机载双频雷达测深系统扫描电机编码数相同的波形进行拼接，保证图像拼接列数方向为飞机飞行方向。将拼接图像利用公式(1)进行灰度归一化：

式中：$ A $为输入的原始图像；$ B $为归一化后图像。归一化后拼接的海水回波能量剖面图，如图4所示。

图  4  海水回波能量剖面图

Figure 4.  Seawater echo energy profile

图中每一列代表一条回波波形，灰度值代表该点的回波能量，从上至下每一行代表一个采样时刻，采样时刻单位为ns。由图4可见一条明显的海底回波廓线，但在500行往下廓线逐渐消失，不再明显，相应波形如图3中波形1所示。为了减小高回波能量对图像灰度值显示范围的拉伸，截取图像401~600行并通过公式(1)进行归一化，用于下一步的图像处理，如图5所示。图像中可以较明显的看到海底回波廓线，相对于图3中波形1无明显海底回波，初步验证了图像处理提取海水深度的可行性。

图  5  401~600 ns范围的海水回波能量剖面图

Figure 5.  Seawater echo energy profile in the range of 401-600 ns

图像预处理包括灰度化、高通滤波、对比度增强等^[8]。为了凸显海底回波廓线区域，消除噪声对后续海底回波廓线的影响需要对图像进行滤波处理。对于任意深度的海底回波信号能量应当与其相邻的海底回波能量相似，因此选择采样双边滤波对图像进行平滑操作并突出海底回波廓线区域，双边滤波在高斯滤波考虑像素空间距离的基础上还考虑了灰度的相似性，在消除噪声的同时保持图像的边缘细节^[9]。双边滤波由公式(2)定义：

其中,$ {W}_{p} $为：

式中：$ p $为输入像素的像素位置；$ q $为滑动窗口$ S $中像素的像素位置；$ {G}_{\delta s} $为空间域核；$ {G}_{\delta \gamma } $为图像像素域核；$ I\left(p\right) $为输入像素灰度值；$ \bar {I}\left(p\right) $为输出像素灰度值；$ I\left(q\right) $为滑动窗口中像素灰度值。进一步，$ {G}_{\delta s} $空间域核如下式：

$ {G}_{\delta \gamma } $图像像素域核如下式：

空间域核和图像像素域核可以视为二维高斯函数，$ \delta s $为空间域核标准差，$ \delta \gamma $为像素域核标准差。在图像的平坦区域，图像灰度值变换小，像素域权重约为1，此时空间域权重起主要作用，相当于进行高斯滤波。在图像边缘区域，空间小范围内灰度值剧烈变化，空间域权重约为1，此时像素域权重变大起主要作用，因而保持了边缘信息。

实验表明，当窗口大小取22，$ \delta s $取5，$ \delta \gamma $取0.7时滤波效果好，如图6所示。

图  6  经过双边滤波后的海水回波能量剖面图

Figure 6.  Seawater echo energy profile after bilateral filtering

对经过双边滤波后的海水回波能量剖面图进行海底回波廓线提取之前，需要对图像进行分割，将海底回波廓线区域和海水散射层回波信号区域分割开来。再单独对感兴趣的区域做进一步的处理。常见的图像分割方法可以分为基于阈值的分割方法、基于区域的分割方法和基于边缘的分割方法等。文中选用基于阈值的分割方法中的Niblack方法。

Niblack算法是一种常见的局部阈值二值化算法，局部阈值通过计算以像素点为中心，窗口大小范围$ r $内的灰度平均值和标准差得到，如果该像素值大于计算的局部阈值，视该像素点为前景，赋值为1，反之视为背景赋值为0，具体阈值计算公式见公式(6)：

式中：$ T $是阈值；$ m $是邻域$ S $中的均值；$ s $是邻域$ S $中的方差；$ k $是预先设定的修正值。

由于局部阈值二值化算法能够根据像素点所在位置不同，确定不同的二值化阈值，实施自适应的二值化处理，因而对背景灰度变换较大，存在突发噪声的图像处理也有较好的效果。但Niblack对于边界区域$ (r-1)/2 $内的像素无法求取阈值，同时如果整个窗口内的像素都为背景，通过Niblack计算后必有一部分像素被误认为前景^[10]。文中选取的Niblack局部阈值窗口大小为13，修正值为0.2，同时为防止海底廓线区域出现断裂，用半径为3的圆结构对图像进行闭操作，处理后图像见图7（a）。

图中除了海底廓线区域外还存在较多错误提取的前景目标，需要对分割的前景目标做进一步的判别以从分割区域中选择出海底回波廓线区域。

通过观察分割目标所在位置对应的双边滤波后海水回波能量剖面图，如图7(b)所示，发现海底回波廓线区域灰度值普遍比背景区域大，因此统计经过双边滤波后的图像上相同位置区域的均值，找到均值最大的区域视为海底廓线区域，至此提取出海底回波廓线区域。对寻找到的海底回波廓线区域按列进行处理，找到每一列第一个为目标的索引值$ t_1 $和最后一个为目标的索引值$ t_2 $，通过公式(7)计算该列的海底回波对应时刻$ t $为：

绘制出海底回波廓线，如图8所示。

图  7  (a) Niblack局部阈值化后的海水回波能量剖面图； (b) 分割目标区域对应海水回波能量

Figure 7.  (a) Seawater echo energy profile after partial thresholding of Niblack; (b) Segmented area corresponds to the seawater echo energy

图  8  海底回波廓线

Figure 8.  Seafloor echo profile

将图5对应数据通过传统波形处理方法进行处理，不进行激光入射角度修正，比对文中方法与传统波形处理方法所求401 ns后的激光斜距，如图9(a)所示。文中方法与传统波形处理在海水斜距50 m之前有很好的一致性，且图像处理方法得到斜距与波形处理方法得到斜距相差如图9(b)所示，两者斜距的平均偏差为0.5148 m，充分说明了文中方法的优越性。

图  9  （a）图像处理方法与波形处理方法比对结果；（b）图像处理斜距减去波形处理斜距后的差值

Figure 9.  (a) Comparison result of image processing method and waveform processing method; (b) Difference between the image processing slope distance and the waveform processing slope distance

绘制另外四组扫描电机编码数不同的海水回波能量剖面图，并处理得到相应的海底回波廓线，如图10所示，对应四组图像处理和波形处理比对结果如图11所示，斜距差值如图12所示。

图  10  不同电机编码数下的海底回波廓线

Figure 10.  Seafloor echo profile under different motor code numbers

图  11  对应图10区域图像处理和波形处理比对结果

Figure 11.  The comparison result of image processing and waveform processing corresponds to the area in Figure 10

图  12  对应图10区域图像处理和波形处理斜距差值

Figure 12.  Image processing and waveform processing slope distance difference corresponds to Figure 10

统计这四组数据通过图像处理方法和波形处理方法探测到海底回波个数、最大探测斜距以及平均偏差见表1。

由表可见，文中方法与传统波形处理方法的斜距平均偏差为0.4547 m，文中图像处理方法在探测海底回波个数和最大探测斜距上均优于传统波形处理方法，在图10（d）区域的图像处理海底回波提取率约是波形处理方法的三倍，在图10（c）区域图像处理方法相对传统波形处理的最大探测斜距提升了约12 m。

为了提升机载海洋激光雷达深水区域的海底回波提取能力，文中提出了基于图像处理提取海水深度的算法，利用双边滤波消除噪声，突出海底回波廓线区域，利用Niblack局部阈值二值化对图像进行分割，提取出海底回波廓线区域，有效提升了海底回波提取能力。文中方法与传统波形处理方法的斜距平均偏差为0.4547 m，且在深水区域，最大探测距离相较波形处理提升了12 m。

下一步，对图像处理方法的参数选择做进一步研究，找到一种自适应的滤波参数设置办法，以避免在海底回波能量变化大的情况下将背景与前景混合的情况。在海水深度提取精度方面，将通过边缘检测的算法，拟合海底廓线，使其得到亚像素的海水深度精度。