基于小波分析的多波束声呐图像角度响应改正方法

朱正任，冯成凯，徐琪尧，韩 冰，付桂合，阳凡林,2

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590； 2.海岛(礁)测绘技术国家测绘地理信息局重点实验室，山东 青岛 266590)

多波束测深系统(multibeam echo sounder, MBES)不仅可以得到高精度的水深地形数据，还能获取丰富的反向散射强度数据。反向散射强度在较大程度上可确定沉积物类型，是识别海底底质类型的重要参数[1-3]。但是，在实际测量过程中，影响多波束反向散射强度的因素较多，除海底底质的物理特性外，声呐系统参数、海水特性(吸收和折射)、测量模式(发射距离和角度)等因素也会对反向散射强度产生影响[4-5]。因此，国内外学者对反向散射强度数据处理进行了广泛的研究，并形成了基于主动声呐方程的较为完善的改正模型。反向散射强度数据经过模型改正后，其中大部分影响可以消除，但角度响应由于受到声波散射机理的影响，表现为反向散射强度随入射角而变化，传统的角度响应改正模型在复杂底质环境下适应性较差，改正效果不理想。

针对多波束声呐图像角度响应改正问题，国内外学者从两种角度出发，针对AR改正后存在的残差进行改正。一种是从数理统计的角度出发，如唐秋华等[6]通过条带两侧强度数据加权内插得到中央波束区域强度数据进行中央区域残差改正；Yang等[7]采用等均值方差拟合的方法，通过对中央区域数据的压缩和移动进行改正；王煜[8]通过构建中央波束区与非中央波束区模型进行残差改正，上述方法均未考虑真实的声波散射机理，改正效果往往与真实底质有较大差异，造成一定的精度损失。另一种是从声学机理角度出发，完善AR改正模型，如Hemmerstand等[9]建立了Lambertian模型进行AR改正；Hellequin等[10]建立了综合声学模型进行AR改正；严俊等[11]通过提取角度响应曲线的AR参数，进行AR改正模型的重构。上述方法均有效减弱了角度响应对声呐图像的影响，但未考虑不同底质对AR改正模型的影响。Zhao等[12]建立了基于预分类的AR改正模型，但未给出高效的分类数的确定方法且计算量较大，存在一定局限性。AR曲线本身亦可以用于分类，但由于其以ping为单位，分类结果的精度低于声呐图像，故本研究有效结合两者的优势，提出一种基于小波分析的多波束声呐图像角度响应改正方法，通过小波分解技术，将AR曲线分解为长波项与短波项，利用AR曲线长波项三维概率密度图先进行同底质区域划分，然后在同底质区域下完成AR改正。该方法可以在不需要确定分类数的前提下，有效消除不同底质对AR改正的影响，获取高质量的声呐图像，为后续高精度的海底底质分类提供基础。

1 角度响应曲线的小波分解及其三维概率密度表示

角度响应曲线是反向散射强度随入射角变化的曲线，由于海底底质的不同，声呐系统接收到的回波强度信号严格地说是非平稳信号[13]，因此，需要用非平稳随机信号的方法进行处理。小波变换在空间域和频率域同时具有良好的分析特性[14]，可以对非平稳信号进行有效分析。故将角度响应曲线进行小波分解，分解后得到其长波项(低频信号)和短波项(高频信号)。然而单ping角度响应曲线的长波项和短波项很难看出其规律性变化，为直观地显示长波项和短波项的趋势变化和分布情况，绘制长波项和短波项的三维概率密度图。

1.1 角度响应曲线的小波分解

(1)

(2)

1.2 短波项三维概率密度图及其特点

短波项代表角度响应曲线的声强特征变化，为直观地显示其反向散射强度随入射角变化的整体趋势，绘制短波项的三维概率密度图。由于反射散射强度随入射角的变化曲线呈对称分布，选取[0°, 60°]进行分析。

图1(a)和1(b)分别是同底质区域和不同底质区域下相邻50 ping角度响应曲线短波项的三维概率密度图(这里以小波参数：小波基coif5、分解层数1为例，说明短波项三维概率密度图的特点)。其中，横坐标代表入射角，纵坐标表示反向散射强度，在指定的入射角处，反向散射强度出现的频率以图右侧的颜色条表示。可以看出，无论在同底质区域下还是不同底质区域下短波项的反向散射强度均随入射角的变化在0附近波动，这是由于短波项代表角度响应曲线变化的细节信息，不包含其整体变化趋势。

图1 短波项的三维概率密度图

1.3 长波项三维概率密度图及其特点

长波项代表角度响应曲线的轮廓线，包含角度响应曲线的主要变化趋势。由于不同的底质类型对应不同的角度响应曲线[15]，因此，长波项三维概率密度图可以直观地显示某一区域是否为同一底质类型。

如图2(a)和图2(b)所示，分别是同底质区域下和不同底质区域下角度响应曲线长波项的三维概率密度图。可以看出图2(a)只包含1种类型的角度响应曲线，认为该区域是同一底质类型；而图2(b)中包含2种类型的角度响应曲线，故该区域存在不同底质类型。

2 多波束声呐图像角度响应改正方法

多波束角度响应曲线经小波分解为长波项和短波项。其中，长波项包含角度响应曲线的主要变化趋势。因此，声呐图像的角度响应改正即角度响应曲线长波项的改正，将改正后的长波项与短波项进行信号重构，得到改正后的多波束声呐图像。

2.1 最优小波分解参数的获取

小波基和分解层数是小波分解的重要参数，不同分解参数对信号分解的影响不同[16]。因此，最优参数的选择是进行角度响应曲线小波分解的重要前提。本研究基于小波分解后的短波项不受角度响应影响的特点，提出了利用短波项三维概率密度曲线的特征参数来确定最优小波分解参数的方法(图3)。

案例教学法还能实现教学相长。在案例教学中，教师不仅是引导者而且也是学习者。一方面，教师掌握着教学进程，引导学生思考、组织讨论研究，并进行总结、归纳。另一方面，收集案例的过程就是一个学习的过程，教师在教学中与学生共同讨论同一话题，也可以从中获得大量感性材料。教师在课堂上不再“独唱”，既能更合理地利用体力和脑力，也能从中获得一种成就感。

图3 最优小波参数获取的流程图

图4(a)是短波项三维概率密度图的剖面图，表示在不同入射角处反向散射强度出现的频率(横坐标为反向散射强度，纵坐标为出现的频率，这里以15°和60°为例)，由图可以看出，在指定入射角处，角度响应曲线的短波项近似服从N(0,σ2)的正态分布，其中σ2为方差。依据中心极限定理，角度响应曲线序列数N>30时，子样BSθ(N)、BSN(θ)序列均服从正态分布[17]，且均值在u=0处波动，因此，其概率密度函数为

(3)

(4)

(5)

图4 短波项三维概率密度的剖面图及其概率密度曲线

图5 统计参数对小波基和分解层数的评价结果

以短波项概率密度曲线的特征参数符合均值与方差近似为0的原则下(确保角度响应曲线的特征趋势不被分解到短波项)，尽可能地保证较多的分解层数(确保角度响应曲线分解得彻底)。由图5(a)和图5(b)可以看出，在不同小波基的情况下，均值绝对值与方差在分解层数为5时产生较大跳变，即分解层数大于5时，角度响应曲线的特征趋势会被分解到短波项中，因此，最优分解层数为5层，结合表1，在相同分解层数为5时，小波基coif5的均值绝对值和方差分别为0.007 8和0.171 9，均小于小波基sym4(0.024 7,0.300 8)和小波基db6(0.011 5,0.317 5)的均值绝对值和方差。为此，将小波基coif5、分解层数为5层作为角度响应曲线小波分解的最优参数。图6(b)和图6(c)，为1 ping角度响应曲线在小波基coif5和分解层数5时，小波分解得到的长波项和短波项，其中，短波项不受角度响应的影响，在强度值0附近波动；长波项变化趋势明显，包含角度响应曲线的主要信息，验证了该参数选择的有效性。

图6 1 ping角度响应曲线的小波分解

图7 多波束条带同底质区域划分流程图

2.2 多波束条带同底质区域划分

在获取角度响应曲线小波分解的最优参数后，还需对多波束条带进行同底质区域划分，来消除不同底质类型对AR改正的影响，进行自适应分区改正。由1.3节可知，长波项三维概率密度图可以判别某一区域是否为同一底质类型，因此，采用二分法对多波束条带进行同底质区域划分，整体流程如图7所示。

由于长波项三维概率密度图无法进行定量分析，绘制同底质区域下和不同底质区域下长波项三维概率密度图(图2)的剖面图。在常见底质下，入射角θ较小时(一般不超过15°)，换能器接收的多为镜面反射，称为D1区；当θ变大时(一般在15°～60°之间)换能器接收的多为后向散射部分，称为D2区，即漫反射区；θ进一步增大时(一般超过60°)，称为D3区或高入射区[12]。因此，选取绘制D1、D2、D3区的边界角度15°和60°处的反向散射强度出现的频率，如图8(a)和图8(b)所示。

图8 长波项三维概率密度图的剖面图

分析图8可得，同底质区域下，长波项三维概率密度图的剖面图只存在1个波峰，当某一区域存在不同底质时，长波项三维概率密度图的剖面图存在2个或多个波峰。故以长波项三维概率密度在入射角15°和60°处的剖面图波峰的个数为判定条件，采用二分法对多波束条带进行同底质区域划分(图9)。

由于条带左右两部分可能存在不同的角度响应曲线，将条带划分为左右两部分进行同底质区域划分(海底底质复杂，以角度响应曲线按块划分无法保证某一区域完全是同一底质，这里对条带进行粗划分，为后续AR改正提供同底质区域范围，提高改正效率)。

2.3 长波项的AR改正及强度信号重构

对于任意一条测线，采用2.1节中的方法可以获得角度响应曲线小波分解的最优参数：小波基coif5和分解层数5；采用2.2节中的方法可以得到多波束条带左右舷的同底质区域。由上述分析可得对声呐图像角度响应改正即对角度响应曲线长波项的改正。因此，算法整体流程如图10所示。

图9 二分法同底质区域划分流程图

为避免单ping信号改正的偶然性，在同底质区域下，将角度响应曲线的长波项叠加取平均，获得该区域长波项的均值曲线BSMean；然后将该区域每1 ping的长波项减去该均值曲线，并统一归化到长波项漫反射区的平均值强度BSM，从而实现长波项的AR改正。即：

(6)

3 实验与分析

为验证本研究算法的有效性，采用高分辨率浅水多波束系统R2SONIC2024在我国某附近海域获取的相邻4条带数据进行实验，该次实验多波束的采样模式为平均强度采样，波束个数为256个，波束开角120°，声呐频率400 kHz，实验区水深20～30 m。选取4条带中的1条测线采用本算法进行AR改正。

首先，将原始多波束反向散射强度数据进行传播损失改正、声照区面积改正等一系列预处理改正；其次，将条带左舷和右舷的角度响应曲线进行小波分解，分解参数采用2.1节中得到的最优分解参数：小波基coif5、分解层数5层，分解后得到角度响应曲线的长波项和短波项；然后，绘制长波项的三维概率密度图，并通过二分法对条带进行同底质区域的划分，划分结果如图11所示，在划分的条带小区域内基本属于同一底质类型。最后，在同底质区域下，进行角度响应曲线长波项的AR改正，将改正后的长波项与短波项进行信号重构，得到改正后的声呐图像(图12(c))。

为了对比现有方法，结果中添加了未进行AR改正的原始声呐图像(图12(a))和传统Lambert法则改正后的声呐图像(图12(b))。通过对比分析图12(a)和图12(c)可以看出，经过AR改正后的多波束声呐图像底质分布均匀，条带中央区域“高亮”的灰度异常现象被消除，同时边缘区域的数据也被归化到平均回波强度上，表明了本算法的有效性。对比图12(b)，由于Lambert法则改正算法只是固定的改正模型，没有考虑不同底质类型AR曲线的不同，因此，在复杂海底底质情况下，无法进行自适应改正，会出现过度改正或改正不彻底的现象，同时也验证了本方法的优越性。

图11 二分法划分区域后的结果

以测线为单位依次消除AR影响，并在地理编码下进行拼接，形成整个区域的海底声呐图像，拼接结果如图13(c)所示。与原始声呐图像拼接后的结果(图13(a))相比，声呐图像中条带中央和边缘区域“高亮”现象得到了有效改善，整个测区灰度变化均匀，底质分界清晰。图13(b)为Lambert法则改正后地理编码的结果，可以看出，该方法在一定程度上减弱了AR的影响，但仍存有一定残差，影响了声呐图像的质量。为进一步量化AR改正的效果，选取同底质区域连续50 ping数据绘制反向散射强度散点图并拟合趋势线，获取其MIC(maximal information coefficient)指数[18]。

图12 AR改正前后声呐图像比较

图13 整个测区的AR改正比较

MIC指数是衡量变量之间相关性的参数，MIC指数越小，说明反向散射强度与入射角的相关性越小，即AR改正效果越好。如图14所示，本算法的MIC指数为0.092，相较于AR改正前的0.433降低了0.341，小于Lambert法则改正的0.274；且改正后的强度点较为集中，趋势线基本呈直线，证明了本算法的有效性。

图14 AR改正效果比较

4 结论

针对多波束声呐图像角度响应改正，本研究综合考虑了不同底质类型对AR改正的影响，提出一种基于小波分析的多波束声呐图像改正方法，该方法利用信号分解与重构的思想，将多波束角度响应曲线分解为长波项和短波项，利用长波项和短波项各自的特点，绘制三维概率密度图，得到小波分解的最优参数和多波束条带同底质区域的划分范围，然后在同底质区域下进行AR改正得到改正后的声呐图像。实验结果表明，对多波束条带同底质区域划分，进行自适应分区改正，改善了传统改正模型适应性差、改正不精确等问题，减弱了角度响应对多波束声呐图像的影响，提高了声呐图像的质量，有效适用于浅水海域地形平坦和底质较复杂的区域，为多波束声呐图像处理提供了参考。