一种基于高斯分布聚类的水下声呐图像匹配方法，通过图像配准和优化声纳三维运动参数，对二维声纳图像进行精准的三维重建，包括以下步骤：步骤A：提取特征，建立特征间的匹配关系，获得地图的高程信息；步骤B：进行位姿估计，将特征地图更新，从而生成三维空间地图。本发明实现了对声纳图像从特征提取到涵盖高程信息的环境地图的重建，同时提供一定的运动姿态估计信息，可用于水下机器人声纳图像处理与建图领域。
1.一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：所述方法包括以下步骤：
步骤A，提取特征，建立特征间的匹配关系，获得地图的高程信息；具体包括以下3个步骤：
步骤A3中计算场景高程图包括以下步骤：步骤A3‑1，根据平面性假设的方法，估计场景相对平坦部分和三维物体上的点的仰角；
基于投射的阴影估算场景中相对平坦部分和3‑D物体上的点的仰角；
步骤A3‑2，使用最小和最大距离值的图像点来固定平面，即前端点和后缘点的三维坐标是固定的，在不同方位扫描出类似的点对，通过投射阴影确定遮挡轮廓点建立目标物体仰角与阴影的对应关系，每个三维对象从前端到后端的高程值由线性插值填充；
步骤A3‑3，通过低通滤波对平滑后的数据进行处理，以减小图像噪声对局部峰值的影响；通过k‑means分割，建立背景、对象和阴影区域的阈值，定位物体到阴影和阴影到地面的转换；
步骤A3‑4，通过后缘点的仰角和遮挡轮廓的角度估算出三维物体的大小，利用下式计算出场景高程图：
式中，Hs为声纳探测到海底的距离，Rs为斜距，Ls为阴影长度，Ht为目标物高度，Lt为目标物长度；
步骤B，进行位姿估计，将特征地图更新，从而生成三维空间地图；具体包括以下步骤：步骤B1，位姿估计；
步骤B3，生成三维地图，包括以下步骤：根据步骤B2优化参数后，计算速度加快得到最优解，计算相对于初始位置的声纳轨迹，完成精确的三维重建。
2.根据权利要求1所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤A1中采集声纳数据，进行特征提取包括以下2个步骤：步骤A1‑1，采集声纳数据，从声纳数据图像中选取像素明显的区域为目标处理区域；
步骤A1‑2，利用目标特征区域灰度分布信息并通过场景分析，将目标区域分为三个类别，分辨出目标像素区域和其他非目标区域，根据灰度分布信息，消除由于噪声产生的像素区域，实现特征点的描述。
3.根据权利要求2所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤A1‑2包括以下2个步骤：步骤A1‑2a，利用目标特征区域灰度分布信息并通过场景分析，将声纳图像中相对明显的灰度值分为三种类别：组成3D目标或结构的明亮小斑点、与3D目标相邻的表面投射的阴影区、介于前两个灰度值之间的平坦表面；
步骤A1‑2b，通过场景分析和像素强度，将目标像素区域和其他非目标区域分辨出，去除由于噪声产生的像素区域，减少图像误差。
4.根据权利要求1所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤A2中图像配准，生成高斯地图包括以下步骤：步骤A2‑1，通过K‑means聚类方法将像素强度较强的目标像素聚类；
步骤A2‑2，采用低通滤波去除噪声，减少对特征提取的误差。
5.根据权利要求4所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤A2‑1包括以下步骤：步骤A2‑1a，对得到的目标特征点网格化，基于网格内的点，通过K‑means聚类法，将每个网格化目标像素进行网格内聚类；
步骤A2‑1b，用高斯分布表示k‑means聚类过的目标特征点，计算协方差矩阵的特征值和特征向量反应表面的方向和光滑性信息，用聚类均值和协方差表示任一聚类，生成高斯地图；均值向量和协方差矩阵对应于每个二维图像区域的位置、大小和方向。
6.根据权利要求5所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤A2‑1a包括以下步骤：选取1％‑2％具有高亮度的像素，剔除小于8的像素，使任何子区域不大于32像素，将大于32像素的区域划分为较小的区域，得到合适的特征区域。
7.根据权利要求1所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤B1位姿估计包括以下步骤：步骤B1‑1，在得到的场景高程图中，选取第一帧声纳图像作为参考图像，第二帧图像作为待匹配图像，每个图像点都通过声纳刚性运动后到新的位置，通过计算帧到帧的运动参数来对两张声纳图像进行配准；
步骤B1‑2，计算声纳运动参数，寻求最适合在转换特性从参考图像到待匹配图像，用空间转换函数进行位姿估计，对两幅图像的变换即旋转平移对比，得到最好的变换参数；然后每次将前一次得到的配准图像作为参考，对下一帧声纳图像进行配准，直至整个图像序列所有待匹配声纳图像配准完成，将所有相邻帧的配准统一到相邻的参考系中，减少成对配准造成的积累误差。
8.根据权利要求1所述的一种基于高斯分布聚类的水下声呐图像匹配方法，其特征在于：步骤B2参数优化包括以下步骤：步骤B2‑1，由步骤A3‑4得到的高程图每个图像点都投影到相应的空间点，经过声纳刚性运动后计算出新的位置三维场景点，得到两个声纳视图；
步骤B2‑2，对两个声纳视图参数优化，将第二个视图的三维点变换到第一个视图的坐标系中，最佳配准评估所有转换的特征高斯分布。
一种基于高斯分布聚类的水下声呐图像匹配方法
技术领域
[0001] 本发明涉及船舶与海洋技术和声纳图像处理领域，特别涉及水下探测与感知以及目标三维重建的方法，具体提出一种基于高斯分布聚类的水下声呐图像匹配方法。
背景技术
[0002] 在过去的十年中，各种基于特征、基于模板、基于区域以及基于傅立叶的方法的二维FS声纳图像配准技术已经被探索，基于各种技术的高分辨率二维(2D)多波束前向扫描
(FS)声纳视频系统已经商业化，可用于浑浊水域的操作和成像。但是以往的研究一般基于
单个的特征点，然后图像配准是制定作为一个优化问题,采用正态分布转换。但是网格大小
的选择涉及分辨率和计算时间之间的权衡，并不一定是自动的。点的分布不服从单变量高
斯分布，以往的高斯分布计算会导致真实分布的错误表示，复杂网格结构计算优化时间过
长。生成高程图时3D投射到2D声纳图像上仰角会丢失，使重建的目标误差大，以及优化方法
的不理想影响图像配准和精确的三维重建。
发明内容
[0003] 本发明的目的是开发一种基于高斯分布聚类的水下声呐图像匹配方法，通过图像配准和优化声纳三维运动参数，通过二维声纳图像进行水下环境三维重建。
[0004] 一种基于高斯分布聚类的水下声呐图像匹配方法，包括以下步骤：
[0005] 步骤A，提取特征，建立特征间的匹配关系，获得地图的高度信息；具体包括以下3个步骤：
[0006] 步骤A1，采集声纳数据，进行特征检测；
[0007] 步骤A2，图像配准，生成高斯地图；
[0008] 步骤A3，计算出场景高程地图；
[0009] 步骤B，进行位姿估计，将特征地图更新，从而生成三维空间地图；具体包括以下步骤：
[0010] 步骤B1，位姿估计；
[0011] 步骤B2，参数优化；
[0012] 步骤B3，生成三维地图。
[0013] 进一步地，步骤A1中采集声纳数据，进行特征提取包括以下4 个步骤：
[0014] 步骤A1‑1，采集声纳数据，从声纳数据图像中选取像素明显的区域为目标处理区域；
[0015] 步骤A1‑2，利用目标特征区域灰度分布信息并通过场景分析，将目标区域分为三个类别，分辨出目标像素区域和其他非目标区域，根据灰度分布信息，消除由于噪声产生的
像素区域，实现特征点的描述。
[0016] 进一步地，其特征在于，步骤A‑1‑2包括以下2个步骤：
[0017] 步骤A1‑2a，利用目标特征区域灰度分布信息并通过场景分析，将声纳图像中相对明显的灰度值分为三种类别：组成3D目标或结构的明亮小斑点、与3D目标相邻的表面投射
的阴影区、介于前两个灰度值之间的平坦表面；
[0018] 步骤A1‑2b，通过场景分析，将组成3D目标或结构的明亮小斑点作为目标区域，将目标像素和其他非目标区域分辨出，去除由于噪声产生的像素区域，减少图像误差。
[0019] 进一步地，步骤A2中图像配准，生成高斯地图包括以下步骤：
[0020] 步骤A2‑1，通过K‑means聚类方法将像素强度较强的目标像素聚类；
[0021] 步骤A2‑2，采用低通滤波去除噪声，减少对特征提取的误差。
[0022] 进一步地，步骤A2‑1包括以下步骤：
[0023] 步骤A2‑1a，对步骤A1‑2b得到的目标特征点网格化，基于网格内的点，通过K‑means聚类法，将每个网格化目标像素进行网格内聚类；
[0024] 步骤A2‑1b，用高斯分布表示k‑means聚类过的目标特征点，计算协方差矩阵的特征值和特征向量反应表面的方向和光滑性等信息，用聚类均值和协方差表示任一聚类，生
成高斯地图。均值向量和协方差矩阵对应于每个二维图像区域的位置、大小和方向。
[0025] 进一步地，步骤A2‑1b包括以下2个步骤：
[0026] 步骤A2‑1b‑1，选取1％‑2％具有高亮度的像素，剔除小于8的像素，选取k值使任何子区域不大于32像素，将大于32像素的区域划分为较小的区域，得到合适的特征区域；
[0027] 步骤A2‑1b‑2，计算均值向量和协方差矩阵对应于每个二维图像区域的位置、大小和方向，生成高斯地图。
[0028] 进一步地，步骤A3中计算场景高程图包括以下步骤：
[0029] 步骤A3‑1，根据平面性假设的方法，计算扫描方位的投射阴影点的仰角的测量值；
[0030] 步骤A3‑2，使用最小和最大距离值的图像点来固定平面，即前端点和后缘点的三维坐标是固定的，在不同方位扫描出类似的点对，通过投射阴影确定遮挡轮廓点建立目标
物体仰角与阴影的对应关系，每个三维对象从前端到后端的高程值由线性插值填充；
[0031] 步骤A3‑3，通过低通滤波对平滑后的数据进行处理，以减小图像噪声对局部峰值的影响；通过k‑means分割，建立背景、对象和阴影区域的阈值，即聚类后根据灰度值分为三
种类别，划分图像类别定位物体到阴影和阴影到地面的转换；
[0032] 步骤A3‑4，通过后缘点的仰角和遮挡轮廓的角度估算出三维物体的大小，利用下式计算出场景高程图：
[0033]
[0034] 式中，Hs为声纳探测到海底的距离，Rs为斜距，Ls为阴影长度， Ht为目标物高度，Lt为目标物长度。
[0035] 进一步地，步骤B1位姿估计包括以下步骤：
[0036] 步骤B1‑1，在步骤A3‑4得到的场景高程图中，选取第一帧声纳图像作为参考图像，第二帧图像作为待匹配图像，每个图像点都通过声纳刚性运动后到新的位置，通过计算帧
到帧的运动参数来对两张声纳图像进行配准；
[0037] 步骤B1‑2，计算声纳运动参数，寻求最适合在转换特性从参考图像到待匹配图像，用空间转换函数进行位姿估计，对两幅图像的变换即旋转平移等对比，得到最好的变换参
数；然后每次将前一次得到的配准图像作为参考，对下一帧声纳图像进行配准，直至整个图
像序列所有待匹配声纳图像配准完成，将所有相邻帧的配准统一到相邻的参考系中，减少
成对配准造成的积累误差。
[0038] 进一步地，步骤B2参数优化包括以下步骤：
[0039] 步骤B2‑1，由步骤A3‑4得到的高程图每个图像点都投影到相应的空间点，经过声纳刚性运动后计算出新的位置三维场景点，得到两个声纳视图；
[0040] 步骤B2‑2，对两个声纳视图参数优化，将第二个视图的三维点变换到第一个视图的坐标系中，最佳配准评估所有转换的特征高斯分布。
[0041] 进一步地，步骤B3生成三维地图包括以下步骤：
[0042] 根据步骤B2优化参数后，计算速度加快得到最优解，计算相对于初始位置的声纳轨迹，完成精确的三维重建。
[0043] 本发明达到的有益效果为：实现了对声纳图像从特征提取到涵盖高程信息的环境地图的重建，同时提供一定的运动姿态估计信息，可用于水下机器人声纳图像处理与建图
领域。
附图说明
[0044] 图1为本发明实施例中整个声纳图像处理的流程图。
[0045] 图2为本发明实施例中声纳数据图像中像素明显的目标处理区域。
[0046] 图3为本发明实施例中通过场景分析得到的三个级别像素区。
[0047] 图4为本发明实施例中海底声纳图和k‑means聚类处理过的图像。
[0048] 图5为本发明实施例中用高斯分布表示k‑means聚类处理目标像素。
[0049] 图6为本发明实施例中利用仰角信息计算高程图。
[0050] 图7为本发明实施例中声纳视图1和通过刚性运动得到的视图2。
具体实施方式
[0051] 下面结合说明书附图对本发明的技术方案做进一步的详细说明。
[0052] 本发明的目的是开发一种基于高斯分布聚类的水下声呐图像匹配方法，通过图像配准和优化声纳三维运动参数，通过二维声纳图像进行水下环境三维重建，整体步骤如图1
所示，包括以下步骤：
[0053] 步骤A：提取特征，建立特征间的匹配关系，获得地图的高度信息；
[0054] 步骤B：进行位姿估计，将特征地图更新，从而生成三维空间地图。
[0055] 步骤A特征检测包括以下步骤：步骤A1采集声纳数据，进行特征检测；步骤A2图像配准，生成高斯地图；步骤A3计算出场景高程地图。
[0056] 步骤A1中采集声纳数据，进行特征提取包括以下2个步骤：
[0057] 步骤A1‑1，采集声纳数据，从声纳数据图像中选取像素明显的区域为目标处理区域，如图2所示。
[0058] 步骤A1‑2，利用目标特征区域灰度分布信息并通过场景分析，将目标区域分为三个类别，分辨出目标像素区域和其他非目标区域，根据灰度分布信息，消除由于噪声产生的
像素区域，实现特征点的描述。
[0059] 步骤A1‑2包括以下2个步骤：
[0060] 步骤A1‑2a，利用目标特征区域灰度分布信息并通过场景分析，将声纳图像中相对明显的灰度值分为三种类别，如图3所示：1.组成 3D目标或结构的明亮小斑点；2.与3D目标
相邻的表面投射的阴影区； 3.介于前两个灰度值之间的平坦表面。
[0061] 步骤A1‑2b，通过场景分析和像素强度，将目标像素和其他非目标区域分辨出，去除由于噪声产生的像素区域，减少图像误差。
[0062] 步骤A2中图像配准，生成高斯地图包括以下步骤：
[0063] 步骤A2‑1，通过K‑MEANS聚类方法将像素强度较强的目标像素聚类，如图4所示。
[0064] 步骤A2‑2，采用低通滤波去除噪声，减少对特征提取的误差。
[0065] 作为本发明的进一步限定，步骤A2‑1包括以下步骤：
[0066] 步骤A2‑1a，对步骤A1‑2b得到的目标特征点网格化，基于网格内的点，通过K‑means聚类法，将每个网格化目标像素进行网格内聚类，如图5，选取1％‑2％具有高亮度的
像素，剔除小于8的像素，选取k值使任何子区域不大于32像素，将大于32像素的区域划分为
较小的区域，得到合适的特征区域。
[0067] 步骤A2‑1b，用高斯分布表示k‑means聚类过的目标特征点，由公式1和公式2计算协方差矩阵的特征值和特征向量反应表面的方向和光滑性等信息，用聚类均值和协方差表
示任一聚类，生成高斯地图。均值向量和协方差矩阵对应于每个二维图像区域的位置、大小
和方向。
[0068]
[0069]
[0070] 其中， 示一个网格内所有扫描点，通过上下标以及取值j j j
范围表示区域j内N 点的集合(共R个区域)中的第 i个点。特征区域j用均值μ和方差Σ高
斯表示如图3所示不同于正态分布，基于集群的高斯映射(统计均值，协方差)表示，消除了
保持一个繁琐的网格结构的必要性，而不失去精度，减少优化过程中的计算时间。
[0071] 步骤A3中计算场景高程图包括以下步骤：
[0072] 步骤A3‑1，根据平面性假设的方法进行计算，估计场景相对平坦部分和三维物体上的点的仰角。平面性假设主要是估算场景中相对平坦部分和3‑D物体上的点的仰角(基于
它们投射的阴影)。
[0073] 步骤A3‑2，使用最小和最大距离值的图像点来固定平面，即前端点和后缘点的三维坐标是固定的，如图6所示，在不同方位扫描出类似的点对，通过投射阴影确定遮挡轮廓
点建立目标物体仰角与阴影的对应关系，每个三维对象从前端到后端的高程值由线性插值
填充。
[0074] 步骤A3‑3，通过低通滤波对平滑后的数据进行处理，以减小图像噪声对局部峰值的影响；通过k‑means分割，建立背景、对象和阴影区域的阈值(在步骤A1‑2中建立)定位物
体到阴影和阴影到地面的转换。
[0075] 步骤A3‑4，利用公式3通过后缘点的仰角和遮挡轮廓的角度估算出三维物体的大小，计算出场景高程图。小物体可能没有明显的阴影，对于这些物体，将物体高度为设为零。
[0076]
[0077] 式中，Hs：声纳探测到海底的距离，Rs：斜距，Ls：阴影长度， Rh：平距，Ht：目标物高度，Lt：目标物长度。
[0078] 步骤B包括以下步骤：步骤B1位姿估计；步骤B2参数优化；步骤B3生成三维地图。
[0079] 步骤B1位姿估计包括以下步骤：
[0080] 步骤B1‑1，在步骤A3‑4得到的场景高程图中，选取第一帧声纳图像作为参考图像，第二帧图像作为待匹配图像，每个图像点都可以通过声纳刚性运动后到新的位置，如图7，
通过计算帧到帧的运动参数来对两张声纳图像进行配准。
[0081] 三维世界投影到二维声纳图像上丢失仰角，可表示为三维场景点 Ps沿三维点的距离和方位角的零高程(Xs，Ys)平面上的映射h(Ps)，计算公式为：
[0082]
[0083] 式中，R为声纳束到反射目标的距离，Φ为方位角。
[0084] 其中声纳坐标系中点P的三维坐标系Ps＝(Xs Ys Zs)T可以用 (R,θ,Φ)T表示，笛卡尔和球面声纳坐标换算公式：
[0085]
[0086]
[0087] 用6个分量描述声纳的三维运动，T＝[tx,ty,tz]T和 W＝[wx,wy,wz]T，二维图像中不包括转动分量[wx,wy,]，通过计算帧到帧的运动参数[tx,ty,tz,wz]对两幅图进行配准。
[0088]
[0089] 步骤B1‑2，计算声纳运动参数，寻求最适合在转换特性从参考图像到待匹配图像，用空间转换函数进行位姿估计，利用公式7将第一幅图投影到第二幅图S′，对两幅图像的变
换即旋转平移等对比。声纳位置运动产生两个场景图(在同一范围R内，不同高度角上)图像
配准通过场景特征的仰角改变通过公式8估计两帧声纳的运动参数[tx,ty,tz,wz]，求出最
小S和S′(公式9)，得到最好的变换参数；然后每次将前一次得到的配准图像作为参考，对下
一帧声纳图像进行配准，直至整个图像序列所有待匹配声纳图像配准完成，将所有相邻帧
的配准统一到相邻的参考系中，减少成对配准造成的积累误差。
[0090]
[0091] S′＝HS (9)
[0092] 式中，H是一个包含平移和平面内旋转的变换矩阵。
[0093] 步骤B2参数优化包括以下步骤：
[0094] 步骤B2‑1，由步骤A3‑4得到的高程图每个图像点都可以用公式 5投影到相应的空间点Ps′在分量m＝(T,W)声纳刚性运动后，新的位置P′三维场景点由公式10和公式11确定。
[0095]
[0096]
[0097] 式中，R为旋转矩阵，W为旋转速度向量。
[0098] 声纳图像对应的S＝h(Ps)和S′＝h(Ps′)由公式9和公式12确定
[0099]
[0100] 经过声纳刚性运动后计算出新的位置三维场景点，得到两个声纳视图。
[0101] 步骤B2‑2包括以下步骤：
[0102] 对步骤B2‑1得到的两个声纳视图由公式13进行参数优化，将第二个视图的三维点变换到第一个视图的坐标系中，最佳配准在评估所有转换的特征高斯分布时产生最大函数
值。
[0103]
[0104] 式中，Gj(s)表示第一个视图中的高斯分布(对应第j个特征)，G′j(s)表示第二个视‑1
图中高斯分布(对应第j个特征)，Psi和P′si表示两个视图特征点集合对应的三维场景，M 将
第二个视图的三维点变换到第一个视图的坐标系中。
[0105] 步骤B3生成三维地图包括以下步骤：
[0106] 根据步骤B2优化参数后，第一幅图象直接映射到S′的第二幅视图中，利用公式8中‑1
的变换矩阵H＝(M,R,Φ)，cosΦ≈1，S＝h(MPs)和 S′＝h(M Ps)由公式14计算。
[0107]
[0108] 齐次变换M，根据公式15递归计算计算相对于初始位置的声纳轨迹，完成精确的三维重建。
[0109] kM0＝kMk‑1k‑1M0 (公式15)
[0110] 式中，k表示帧数(时间指数)。
[0111] 以上所述仅为本发明的较佳实施方式，本发明的保护范围并不以上述实施方式为限，但凡本领域普通技术人员根据本发明所揭示内容所作的等效修饰或变化，皆应纳入权
利要求书中记载的保护范围内。
