CN112686935B - 基于特征融合的机载测深雷达与多光谱卫星影像配准方法 - Google Patents

Abstract

本发明公开了一种基于特征融合的机载测深雷达与多光谱卫星影像自动配准方法，需要计算机载测深雷达(ALB)点云数据的地形梯度和强度梯度，将三维点云数据转化成二维，同时计算影像数据红、绿、蓝三种波段的梯度图，使用基于卷积神经网络构建图像融合模型分别对两种数据进行特征图融合，利用聚类分析剔除融合图像的干扰特征，然后基于SIFT‑RANSAC方法进行粗配准，然后利用一种改进的Demons算法进行精配准；利用粗配准和精配准的结果对卫星影像进行纠正，最终纠正后的卫星影像即为与ALB数据相配准的卫星影像。

Description

基于特征融合的机载测深雷达与多光谱卫星影像配准方法

技术领域

本发明属于遥感数据处理领域，具体涉及一种基于特征融合的机载测深雷达与多光谱卫星影像自动配准方法。

背景技术

机载激光测深系统(ALB)是近三十年来发展起来的一种海洋测量技术，可用于对海岸线和水下进行无缝测绘。高空间分辨率、高光谱分辨率和高时间分辨率的多光谱图像采集和处理技术的发展，可以为我们提供更多的光谱信息和纹理细节，以便进一步的海岸带监测和分析。ALB数据可以得到精确的地物坐标，而卫星图像具有丰富的光谱信息和纹理特征。两种数据源都有明显的优点和局限性。为了弥补彼此的局限并结合它们各自的优势，需要将这两个数据源融合起来，以提供更丰富和有价值的信息。

三维激光点云数据(3D)与二维光谱图像(2D)的融合与应用，首先要解决的问题是几何配准。如何通过建立一个功能转换模型来对齐两个数据源的元素一直是一个长期的研究问题。尽管有许多经过地理校准和正射校正的图像产品可用，但互补数据集之间经常存在配准错误。因此，在数据融合和后续产品生成之前，需要精确的配准。另外，激光点云与光学图像的巨大差异给配准带来了困难。与传统的图像配准不同，激光点云和光学图像属于跨模态异构数据。为了满足实际应用，目前进行配准的方法可分为基于物体特征的方法和基于物体强度的方法。前者着重于点、线、斑块等显著特征，以建立不同数据集之间的对应关系。基于显著特征的配准方法不仅提高了配准精度，而且降低了计算复杂度。基于强度的配准方法可以通过计算数据集之间的相似度来确定相应的关系。这些方法具有较强的鲁棒性，广泛应用于医学影像领域中由不同类型传感器获取的医学图像配准，不需要复杂的特征提取和特征匹配。用从图像和点云强度数据中提取的尺度不变特征变换(SIFT)特征来确定两者之间的对应关系。此外，互信息已经被用来确定点云强度图和图像之间的统计和功能关系。然而，没有高效、鲁棒的点云预处理，无法获得激光雷达数据的准确特征，这也是当前研究的热点。此外，高质量的空间几何需要足够的数据密度来支持。并不是所有的点云数据都可以使用基于特征提取的配准算法对激光雷达数据和图像进行配准。

对于陆地区域卫星影像与陆地点云之间已经有比较多的模型。但是对于海岛、海岸带区域ALB数据与影像数据配准，存在以下几个问题：(1)卫星影像几何校正不充分，使得影像与点云对应关系并不是一个刚性的，增加了配准难度；(2)二维数据与三维数据是完全不同的空间表达，如何进行不同维度两种数据是需要解决的问题；(3)海岛礁区域显著特征(点、线、面)有限给配准带来困难，如何充分利用有限的特征是值得思考的一个问题；(4)ALB数据水上、水下特征差异太大，如何避免水下特征损失，助力影像配准是需要解决的问题。

发明内容

为解决上述问题，得到准确的配准结果，本研究充分考虑了两种数据的特点和优势，利用深度卷积神经网络的方法得到不同数据的特征图像，提出了一种基于特征融合的机载测深雷达与多光谱卫星影像自动配准方法。

本发明采用的技术方案是：第一步，对ALB点云数据和卫星影像数据进行预处理，包括ALB点云数据栅格化，点云地形和强度的梯度计算，多光谱影像不同波段梯度计算、ALB水下梯度特征调整，得到ALB和卫星影像的梯度图像；第二步，利用深度卷积神经网络计算待融合梯度图像的权值图像，然后利用拉普拉斯金字塔分解、高斯金字塔分解得到源图像分别分解到不同的空间频带上，然后对特征图进行融合，再利用拉普拉斯金字塔重构得特征融合后的图像，同样的方法依次得到ALB空间特征、强度特征的融合图像以及卫星影像不同波段梯度特征的融合图像；第三步是利用聚类分析对融合图像进行滤波，进而用SIFT-RANSAC算法进行粗配准；第四步，利用改进的Demons算法进行融合图像的精配准；第五步，利用粗配准和精配准的结果对卫星影像进行纠正，最终纠正后的卫星影像即为与ALB数据相配准的卫星影像。其具体步骤如下：

步骤一，ALB点云数据和卫星影像数据进行预处理，包括ALB点云地形数据和强度数据栅格化，点云地形和强度的梯度计算，多光谱影像不同波段梯度计算，以及ALB水下梯度特征调整，最终得到ALB梯度图像和卫星影像的梯度图像；

步骤二，利用深度卷积神经网络计算待融合梯度图像的权值图像，然后利用拉普拉斯金字塔分解、高斯金字塔分解得到源图像分别分解到不同的空间频带上，然后对特征图进行融合，再利用拉普拉斯金字塔重构得到特征融合后的图像，利用上述方法分别得到ALB点云地形梯度图像和强度梯度图像的融合图像1，以及卫星影像不同波段梯度图像的融合图像2；

步骤三，利用聚类分析对融合图像1和融合图像2分别进行滤波，得到融合图像1′和融合图像2′，进而执行SIFT-RANSAC算法计算融合图像1′和融合图像2′的旋转矩阵ΔR和偏移矩阵ΔT，实现粗配准；

步骤四，利用Demons算法对粗配准后的融合图像1′和融合图像2′进行精配准，得到最终得到的偏移矩阵为T_final；

步骤五，结合两幅影像的旋转矩阵ΔR和偏移矩阵ΔT和步骤四中计算的T_final对原始的卫星影像I_RS进行纠正得到I′_RS，最终I′_RS即为和ALB数据相配准的卫星影像。

进一步的，步骤一中使用反距离加权插值方法对ALB点云地形数据和强度数据进行插值，得到栅格数据；

以ALB点云地形数据为例，公式(1)计算插值点(x,y)到离散点(x_i,y_i)的距离，利用公式(2)计算每个点的权重，权重是距离的倒数的函数，然后用公式(3)得到插值点的坐标得到栅格数据；

n为离散点个数，z_i为离散点z坐标，用同样的方法插值得到强度的栅格数据。

进一步的，步骤一中利用公式(4)-(6)分别计算ALB点云地形和强度的梯度图像，以点云的地形数据为例，点云空间结构特征，即地形梯度图计算公式如下；

公式(4)-(6)中S_x,S_y分别为计算x、y方向导数值的两个滤波核，f为待计算的栅格数据，利用同样的方法计算强度梯度图和卫星影像不同波段的梯度图。

进一步的，步骤一中利用公式(7)对ALB水下梯度特征进行调整，使其在特征图中不被忽略，以强度特征I_new为例子，I_old为待放大的特征，I_new缩放后的范围为[max_new，min_new]，I_old对应原始数据的范围[max_old,min_old]；

进一步的，步骤二中得到ALB点云地形梯度图像和强度梯度图像的融合图像1的具体实现方式如下；

步骤2.1，通过CNN生成权重图W，将强度梯度图像I和地形梯度图像V分别输入到CNN的两个分支中，获得权值图W，每个分支的权值相同，每个分支由3个卷积层和1个最大池化层组成，拼接后生成的512个特征图直接连接到二维矢量上，并送入双向softmax层，然后生成两类的概率分布，代表两种标准化的权重分布；

步骤2.2，将每个强度梯度图像I和地形梯度图像V分解为拉普拉斯金字塔L{I}和L{V}，L{I}和L{V}分别表示I和V的金字塔，将权重图W分解为高斯金字塔G{W}^l，l表示分解的层数，每个金字塔的分解总水平设置为最高的可能值H×W是强度梯度图像I和地形梯度图像V的空间大小；

步骤2.3，系数融合，利用公式(8)计算L{I}和L{V}的l级局部能量图局部能量图和

m,n滤波器尺度，用公式(9)计算融合模式的相似度：

该测度的范围为[-1，1]，数值越大，相似度越高，设置阈值t以确定所使用的融合模式，融合策略为公式(10)：

L{F}^l(x,y)为融合后的l层拉普拉斯金字塔，G{W}^l(x,y)为l层高斯金字塔；

步骤2.4，拉普拉斯金字塔重构；对融合后的拉普拉斯金字塔，从其顶层开始逐层从上至下按下式进行递推，恢复其对应的高斯金字塔，并得到最终的融合图像G₀；

进一步的，步骤三的具体实现方式如下；

步骤3.1，聚类分析，剔除干扰特征；在增强方法EnFCM框架下，对融合图像的灰度直方图进行聚类，通过迭代运算得到模糊隶属度矩阵，最后，利用滤波器对模糊隶属度矩阵进行修正，聚类灰度直方图的目标函数设计为：

u_kl表示图像灰度l相对于第k个聚类中心v_k的模糊隶属度，ξ是模态学重构后的图像，ξ_l是灰度级，q表示包含在ξ中的灰度级数量，1≤l≤q，通常远小于N，为形态闭合重构，f为输入的融合图像；然后，将上述优化问题转化为目标函数最小的无约束优化问题；

其中q和c分别表示图像中聚类中心的个数和灰度级的个数，γ_l表示灰度值为l的像素个数，表示图像灰度l相对于第k个聚类中心v_k的模糊隶属度，m为模糊参数，λ是拉格朗日乘子；因此，目标函数的最小化问题转化为寻找上面拉格朗日函数的鞍点和衍生品的拉格朗日参数，其中分割结果由模糊隶属度矩阵u_kl和聚类中心v_k共同确定，其数学表达式为：

为了得到一个稳定的模糊隶属度矩阵U＝(u_kl)^c×q，重复执行公式(15-17)直到矩阵{U^(t)-U^(t+1)}<η，t为迭代次数，其中η为最小误差阈值，通过隶属度滤波对u_kl进行修正，得到更好的隶属度划分矩阵U′，此外，还引入了中值滤波算法med以提高隶属度滤波性能和算法速度；

U″＝med{U′} (19)

根据聚类分析设置类别数，得到最终的隶属度划分矩阵U″，对干扰类别进行自动剔除，其中干扰类别的设置需要根据具体情况来确定；融合图像1和融合图像2经过上述聚类分析后得到的结果分别为融合图像1′和融合图像2′；

步骤3.2，执行SIFT-RANSAC计算两幅影像的旋转矩阵ΔR和偏移矩阵ΔT，对融合图像1′和融合图像2′进行粗配准。

进一步的，步骤四的具体实现方式如下；

步骤4.1，在SIFT-RANSAC配准结果的基础上对待配准图像进行格网划分，设粗配准后的融合图像1′进行格网划分后得到I_coarse1，格网框选区域内的I_coarse1为静止图像用f表示，粗配准后的融合图像2′进行格网划分后得到格网框选区域I_coarse2，格网框选区域内的I_coarse2称为运动图像用m表示；

步骤4.2，依次遍历4.1中格网划分的区域，对应格网内的图像利用改进的Demons算法进行非刚性配准，执行如下操作：

第一步：模态变换，具体实现方式如下；

运动图像m的直方图图像为T_m(m)，静止图像直方图图像为T_f(f)，联合灰度直方图是直方图图像为T_mf(m,f)可以通过2-D模态变换矩阵来表达。

G_m和G_f分别为m和f的最大灰度值，T_mf(x,y)是m中灰度值为x，f中对应点灰度值为y的像素对的总个数，即：

T_mf(m,f)＝T_m(m)·T_f(f) (21)

而T_m(m)和T_f(f)可以通过分别对图像m和f进行灰度统计得到。通过联合直方图得到联合概率密度的计算式为：

第二步：计算互信息，具体实现方式如下；

计算运动图像m与静止图像f的互信息为：

I(m,f)＝H(m)+H(f)-H(m,f) (23)

其中，

H(m,f)＝∑_m,fp(m,f)log(p(m,f)) (24)

m的熵H(m)＝∑_mp(m)log(p(m))，同理f的熵为H(f)，H(m,f)是联合熵。p(m)和p(f)是m和f中灰度值的概率，p(m,f)是图像灰度值的联合概率，将(24)带入(23)得到：

将公式(22)带入到公式(25)，计算模态变换矩阵下两个图像的互信息：

NM＝∑_m,fT_mf(m,f)是联合直方图所有元素值的总和；

第三步，计算配准误差函数，具体实现方式如下；

Demons算法的工作框架借鉴了光流算法，光流算法用于估计视频图像中相邻图像目标的位移，称为目标移动速度，为了防止图像梯度为0时计算出现问题，在传统的位移计算公式的分母中加了一项图像灰度差(m-f)，静态图像f中p点与运动图像m中对应点相匹配所需的估计位移u定义如下：

其中u＝(u_x,u_y)，表示m变换成f的位移，包含x和y不同方向的位移；是静态图像的梯度；是移动图像的梯度；配准的目的就是为了得到运动图像的位移，通过引入一个包含运动图像边缘力的方程，以提高配准收敛速度和稳定性，此外应用归一化因子α调整力的强度，则改正后的估计位移u为：

初始u设置为0，此时u的上限为1/α，因此通过设置α的值来调节形变的幅度；

定义的配准误差函数为：

S为形变位移，描述每个像素在x和y方向的位移，代表对图像的变换，U为每次迭代计算后S的偏移量，σ_i和σ_x代表强度不确定性和变换不确定性；根据当前形变参数S变换后的运动图像重新表示为在求新的u时的配准误差函数：

那么配准的过程就是一个寻找E最小的过程，计算误差梯度：

求E最小转换为则：

每一次迭代得到的变化矩阵都要进行高斯滤波，使得变化正则化，抑制噪声，防止出现错误形变，在变换过程中，新图像的边缘会出现一些缺陷，因此，模态变换后的图像在一定程度上失去了作为边缘力的能力，因此加入对误差函数E进行了修改：

为了避免局部极小值，加快配准速度，采用尺度空间方法；

步骤4.3，将遍历所有格网得到的m各像素沿x轴和y轴方向的累积偏移量矩阵∑U按照格网位置拼接得到I_coarse2所有像素点的位移T_all＝∑_格网∑U；

步骤4.4，对I_coarse2进行聚类分析得到M′；

步骤4.5，对图像M′进行二值化处理，将累积偏移量矩阵T_all映射到图像M′上，对应M′上值为1的区域对应的偏移矩阵被保留，对应M′上值为0的区域对应的偏移矩阵被删除，其中被删除的偏移信息由插值得到，最终得到的偏移矩阵为T_final，此外，为了避免局部错误配准点的干扰，采用高斯滤波算法对偏移矩阵T_final进行滤波操作。

本发明具有如下优点：创新性提出一种基于特征融合的点云与影像配准的方法，解决了三维数据与二维数据的配准难题，充分利用了ALB数据空间结构特征和强度特征，充分利用了水下特征，提高了其在配准中的贡献，聚类分析可有效剔除影响配准的特征信息，提高特征点正确匹配的概率，精配准模型的提出可以有效提高点云和影像之间的配准精度。

附图说明

图1为本发明的流程图；

图2为本发明提取的梯度特征图，左上为地形梯度图，右上为地形梯度水上部分和水下部分统计直方图；左下为强度梯度图，右下为强度梯度水上部分统计直方图；

图3为本发明提取的特征融合图像，左边：a1-a4为地形梯度图中的四个区域；b1-b4为强度梯度图中的四个区域，c1-c4为融合后图像四个区域，右边：从左至右每列依次为a，b，c中四个区域的局部放大图；

图4为本发明提取的聚类分析结果，(a1)ALB特征融合图像,(a2)ALB特征融合图像聚类分析结果，(a3)ALB特征融合图像剔除干扰特征后结果，(b1)卫星影像特征融合图像,(b2)卫星影像特征融合图像聚类分析结果，(b3)卫星影像特征融合图像剔除干扰特征后结果；

图5为本发明实验结果与常规方法结果对比。

具体实施方式

如图1所示，本发明提供的一种基于特征融合的机载测深雷达与多光谱卫星影像自动配准方法，具体包括如下步骤：

第一步，对ALB点云数据和卫星影像数据进行预处理，包括ALB点云坐标数据和强度数据栅格化，点云地形和强度的梯度计算，多光谱影像不同波段梯度计算，ALB水下梯度特征调整，得到ALB和卫星影像的梯度图像；

步骤1.1：使用反距离加权插值方法对ALB点云坐标数据和强度数据进行插值，以坐标数据为例，公式(1)计算插值点(x,y)到离散点(x_i,y_i)的距离，利用公式(2)计算每个点的权重：权重是距离的倒数的函数，然后用公式(3)得到插值点的坐标得到栅格数据；

n为离散点个数，z_i为离散点z坐标，用同样的方法插值得到强度的栅格数据。

步骤1.2：梯度图计算，然后用公式(4)-(6)分别计算ALB地形和强度的梯度图像。同样以点云的坐标数据为例，点云空间结构特征，即梯度图计算公式如下；

公式(4)-(6)中S_x,S_y分别为计算x、y方向导数值的两个滤波核，f为待计算的栅格数据，用同样的方法计算强度梯度图和卫星影像(I_RS)不同波段的梯度图；

步骤1.3：水下梯度特征调整，对海水覆盖区域的ALB地形和强度特征进行放大，一般情况下水下特征不太明显，海岛海岸带区域既有水上又有水下区域，如果只关注水上区域显然是不合适的，原则上其在配准中是占据重要位置的，因此对水下特征利用公式(7)进行放大使其在特征图中不被忽略。以强度特征I_new为例子，I_old为待放大的特征，I_new缩放后的范围为[max_new，min_new],I_old对应原始数据的范围[max_old,min_old]。

本发明中设置水下显著特征的前50％与陆地显著特征的前30％范围一致。

第二步，利用基于深度卷积神经网络的特征融合模型对特征进行融合，分别对ALB数据得到的两幅梯度图像(地形梯度图和强度梯度图)进行融合，以及对卫星影像得到的三个不同波段梯度图像进行特征融合。以前者为例，首先计算待融合梯度图像的权值图像，然后利用拉普拉斯金字塔分解、高斯金字塔分解得到源图像分别分解到不同的空间频带上，然后对特征图进行融合，再利用拉普拉斯金字塔重构得特征融合后的图像；

步骤2.1：通过CNN生成权值图W。将强度梯度图像I和地形梯度图像V分别输入到CNN的两个分支中。然后使用下面描述的方法获得权值图W。在CNN图像融合算法中，约束两个分支的权值相同，每个分支由3个卷积层和1个最大池化层组成。拼接后生成的512个特征图直接连接到二分类输出上，并送入双向softmax层。然后生成两类的概率分布，代表两种标准化的权重分布。也就是说，每类的概率值代表了每项权重分配的概率。每个类别的概率对应于输入小块的权重，和为1。用高质量的图像和大小为16×16的模糊图像(不同焦点)训练网络。然后采用基于随机梯度下降算法的softmax最大损失函数作为优化目标进行最小化。由于完全连通层转化为包含两个核(8×8×512)的等效卷积层，因此网络可以处理任意大小的源图像。然后这些值被用作分支位置中所有像素的权重。我们对重叠的像素进行平均，以获得与输入的强度梯度图像I和地形梯度图像V相同大小的权值图W。

步骤2.2：对高斯金字塔的每一层图像进行高斯卷积，然后将前一层图像上采样的预测图像减去所得到的一系列差分图像，即为拉普拉斯金字塔分解图像。其目的是将源图像分解成不同的空间频带，等效于带通滤波的构造过程。将每个强度梯度图像I和地形梯度图像V分解为拉普拉斯金字塔L{I}和L{V}。L{I}和L{V}分别表示I和V的金字塔。将权重图W分解为高斯金字塔G{W}^l。l表示分解的层数，每个金字塔的分解总水平设置为最高的可能值H×W是源图像(强度梯度图像I和地形梯度图像V)的空间大小。

步骤2.3：系数融合，用公式(8)计算L{I}和L{V}的l级局部能量图和

m,n滤波器尺度，用公式(9)计算融合模式的相似度：

该测度的范围为[-1，1]，数值越大，相似度越高，设置阈值t以确定所使用的融合模式。融合策略为公式(10)：

L{F}^l(x,y)为l层拉普拉斯金字塔，G{W}^l(x,y)为l层高斯金字塔。

步骤2.4：拉普拉斯金字塔重构，对融合后的拉普拉斯金字塔，从其顶层开始逐层从上至下按下式进行递推，可以恢复其对应的高斯金字塔，并最终可得到最终的融合图像G₀，该图像作为融合图像1。

步骤2.5：用步骤2.1-2.4中的方法得到ALB坐标梯度图与强度梯度图的融合图像1，同理，利用步骤2.1-2.4中的方法在对卫星影像红、绿、蓝波段的梯度图像进行融合得到一副融合图像2。

第三步是利用聚类分析对融合图像进行滤波，进而用SIFT-RANSAC算法对融合图像1和融合图像2进行粗配准；

步骤3.1：聚类分析，剔除干扰特征。在增强方法(EnFCM)框架下，对融合图像的灰度直方图进行聚类，通过迭代运算得到模糊隶属度矩阵。最后，利用滤波器对模糊隶属度矩阵进行修正。因为聚类灰度直方图与EnFCM相似，目标函数设计为：

u_kl表示图像灰度l相对于第k个聚类中心v_k的模糊隶属度。ξ是模态学重构后的图像，ξ_l是像素点的灰度值，q表示包含在ξ中的灰度级数量,1≤l≤q，通常远小于N。为形态闭合重构，f为输入的融合图像。然后，将上述优化问题转化为目标函数最小的无约束优化问题。

其中q和c分别表示图像中聚类中心的个数和灰度级的个数。γ_l表示灰度值为l的像素个数，m为模糊参数。λ是拉格朗日乘子。因此，目标函数的最小化问题转化为寻找上面拉格朗日函数的鞍点和衍生品的拉格朗日参数。其中分割结果由划分矩阵u_kl和聚类中心v_k共同确定。其数学表达式为:

为了得到一个稳定的模糊隶属度矩阵U＝(u_kl)^c×q，重复执行公式(15-17)直到矩阵{U^(t)-U^(t+1)}<η，t为迭代次数，其中η为最小误差阈值。通过隶属度滤波对u_kl进行修正，得到更好的隶属度划分矩阵U′，加快算法收敛速度。此外，为了提高隶属度滤波性能和算法速度，还引入了中值滤波算法med。

U″＝med{U′}, (19)

根据聚类分析需要设置类别数，得到最终的隶属度划分矩阵U″，对干扰类别进行自动剔除，其中干扰类别的设置需要根据具体情况来确定。融合图像1和融合图像2经过上述聚类分析后得到的结果分别为融合图像1′和融合图像2′。

步骤3.2：执行SIFT-RANSAC计算两幅影像的旋转矩阵ΔR和偏移矩阵ΔT，对融合图像1′和融合图像2′进行粗配准，粗配准在一定程度上减少了两种数据的空间几何距离，但是不能满足精度要求，因此需要进行精配准。

第四步，利用改进的Demons算法对粗配准后的融合图像1′和融合图像2′进行精配准，得到最终的偏移矩阵为T_final。

步骤4.1：在SIFT-RANSAC配准结果的基础上对待配准图像(即粗配准后的融合图像1′和融合图像2′)进行格网划分，划分尺度视情况而定，粗配准后的融合图像1′进行格网划分后得到I_coarse1，粗配准后的融合图像2′进行格网划分后得到I_coarse2。

步骤4.2：依次遍历4.1中格网划分的区域，对应格网内的图像利用改进的Demons算法进行进一步的非刚性精配准，其中执行如下操作，此时的格网框选区域内的I_coarse2称为运动图像用m表示，格网框选区域内的I_coarse1称为静止图像用f表示：

第一：模态变换。

尽管ALB融合图像与光学图像整体的梯度图整体趋势一致，但在细节上存在很大差异，这为配准带来了一定困难。该发明将模态变换应用到Demons配准中，允许对多模态进行配准。为了实现基于强度的局部图像快速配准，利用联合直方图峰值将两种图像表示转换为另一种图像表示，即两幅图像的相互直方图的极大值，可以用于将一幅图像中的区域灰度表示转换为另一个图像的表示。运动图像m的直方图图像为T_m(m)，静止图像直方图图像为T_f(f)，联合灰度直方图是直方图图像为T_mf(m,f)可以通过2-D模态变换矩阵来表达。

G_m和G_f分别为m和f的最大灰度值，T_mf(x,y)是m中灰度值为x，f中对应点灰度值为y的像素对的总个数，即：

T_mf(m,f)＝T_m(m)·T_f(f) (21)

而T_m(m)和T_f(f)可以通过分别对图像m和f进行灰度统计得到。通过联合直方图得到联合概率密度的计算式为：

第二：计算互信息。计算运动图像m与静止图像f的互信息为：

I(m,f)＝H(m)+H(f)-H(m,f) (23)

其中，

H(m,f)＝∑_m,fp(m,f)log(p(m,f)) (24)

m的熵H(m)＝∑_mp(m)log(p(m))，同理f的熵为H(f)，H(m,f)是联合熵。p(m)和p(f)是m和f中灰度值的概率，p(m,f)是图像灰度值的联合概率，将(24)带入(23)得到：

互信息配准算法默认待配准图像具有一定的均匀强度区域，这些区域由彼此对应的边分隔开来，但强度和纹理不相同。

将公式(22)带入到公式(25)，计算模态变换矩阵下两个图像的互信息：

NM＝∑_m,fT_mf(m,f)是联合直方图所有元素值的总和。

第三：计算配准误差函数。

Demons算法的工作框架借鉴了光流算法，光流算法用于估计视频图像中相邻图像目标的位移，称为目标移动速度。为了防止图像梯度为0时计算出现问题，在传统的位移计算公式的分母中加了一项图像灰度差(m-f)，静态图像F中p点与运动图像M中对应点相匹配所需的估计位移u定义如下：

其中u＝(u_x,u_y)，表示m变换成f的位移，包含x和y不同方向的位移；是静态图像的梯度；是移动图像的梯度。配准的目的就是为了得到运动图像的位移。本发明引入了一个包含运动图像边缘力的方程，以提高配准收敛速度和稳定性，此外应用归一化因子α调整力的强度，则改正后的估计位移u为：

初始u设置为0，此时u的上限为1/α,因此可以通过设置α的值来调节形变的幅度,可以在最初的几次迭代中使用较小的α，随着迭代次数的增加，使用大一点的α，α用I(m,f)代替。

定义的配准误差函数为：

S为形变位移，描述每个像素在x和y方向的位移。代表对图像的变换，U为每次迭代计算后S的偏移量。σ_i和σ_x代表强度不确定性(图像噪声)和变换不确定性。根据当前形变参数S变换后的运动图像重新表示为在求新的u时的配准误差函数：

那么配准的过程就是一个寻找E最小的过程，计算误差梯度：

求E最小转换为则：

每一次迭代得到的变化矩阵都要进行高斯滤波，使得变化正则化，抑制噪声，防止出现错误形变。在变换过程中，新图像的边缘会出现一些缺陷。因此，模态变换后的图像在一定程度上失去了作为边缘力的能力。因此我们加入对误差函数E进行了修改：

步骤4.3：将遍历所有格网得到的m各像素沿x轴和y轴方向的累积偏移量矩阵∑U按照格网位置拼接得到I_coarse2所有像素点的位移T_all＝∑_格网∑U。

步骤4.4：对I_coarse2利用步骤3.1进行聚类分析得到M′。

步骤4.5：对图像M′进行二值化处理，将累积偏移量矩阵T_all映射到图像M′上。对应M′上值为1的区域对应的偏移矩阵被保留，对应M′上值为0的区域对应的偏移矩阵被删除。其中被删除的偏移信息由插值得到，最终得到的偏移矩阵为T_final。此外，为了避免局部错误配准点的干扰，采用高斯滤波算法对偏移矩阵T_final进行滤波操作。

第五步：结合步骤3.2计算两幅影像的旋转矩阵ΔR和偏移矩阵ΔT和步骤4.5计算的T_final对卫星影像I_RS进行纠正得到I′_RS，最终I′_RS即为和ALB数据相配准的卫星影像。

本实例选取2012年12月9日在海南岛蜈支洲附近获取的ALB数据和2012年9月28日获取的WorldView2多光谱影像数据，给出本发明的一个实施例；

(一)ALB和影像数据梯度计算；

对ALB点云数据利用公式(1)-(3)进行反距离加权插值后得到栅格数据(3D)，然后利用公式(4)-(6)进行梯度计算得到梯度图(依次计算点云坐标梯度图和点云强度梯度图)(2D)，如图2所示，左上为地形梯度图，右上为地形梯度水上部分和水下部分统计直方图；左下为强度梯度图，右下为强度梯度水上部分统计直方图。然后分别计算WorldView2红、绿、蓝三个波段的梯度特征图(2D)；然后用公式(7)对水下区域的特征信息进行放大。

(二)使用卷积神经网络进行特征融合；

按照步骤二利用基于卷积神经网络构建的融合模型对ALB的两幅特征图像进行融合，使其同时具有空间结构特征和属性特征。同样的利用该模型对多光谱卫星影像的三个波段的特征进行融合，得到最终的融合图像，如图3所示，a为地形梯度图；b为强度梯度图，c为融合后图像，其中四个区域的局部放大图被显示。

(三)粗配准；

按照步骤3.1对融合后图像进行聚类分析，结果如图4，(a1)ALB特征融合图像,(a2)ALB特征融合图像聚类分析结果，(a3)ALB特征融合图像剔除干扰特征后结果，(b1)卫星影像特征融合图像,(b2)卫星影像特征融合图像聚类分析结果，(b3)卫星影像特征融合图像剔除干扰特征后结果，然后进行步骤3.2的粗配准。

图3聚类分析结果

统计发现，加入聚类分析可以有效提高配准点对的正确率。

(四)精配准，利用步骤四中的操作对点云与影像数据进行精配准，得到配准结果，选取了九处进行精度比较，得到的结果如图5；精配准比粗配准在配准精度上提高了3.5m左右。

测试分析：将提取出基于特征融合的点云与影像自动配准方法可以自动的实现两种数据的配准，比传统的配准方法具有更高的精度。

(五)结合粗配准和精配准结果对原始卫星影像I_RS进行纠正得到I′_RS，最终I′_RS即为和ALB数据相配准的卫星影像。

本文中所描述的具体实施例仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于，该方法具体包括以下步骤：

步骤一，ALB点云数据和卫星影像数据进行预处理，包括ALB点云地形数据和强度数据栅格化，点云地形和强度的梯度计算，多光谱影像不同波段梯度计算，以及ALB水下梯度特征调整，最终得到ALB梯度图像和卫星影像的梯度图像；

步骤二，利用深度卷积神经网络计算待融合梯度图像的权值图像，然后利用拉普拉斯金字塔分解、高斯金字塔分解得到源图像分别分解到不同的空间频带上，然后对特征图进行融合，再利用拉普拉斯金字塔重构得到特征融合后的图像，利用上述方法分别得到ALB点云地形梯度图像和强度梯度图像的融合图像1，以及卫星影像不同波段梯度图像的融合图像2；

步骤二中得到ALB点云地形梯度图像和强度梯度图像的融合图像1的具体实现方式如下；

步骤2.1，通过CNN生成权重图W，将强度梯度图像I和地形梯度图像V分别输入到CNN的两个分支中，获得权值图W，每个分支的权值相同，每个分支由3个卷积层和1个最大池化层组成，拼接后生成的512个特征图直接连接到二维矢量上，并送入双向softmax层，然后生成两类的概率分布，代表两种标准化的权重分布；

步骤2.2，将每个强度梯度图像I和地形梯度图像V分解为拉普拉斯金字塔L{I}和L{V}，L{I}和L{V}分别表示I和V的金字塔，将权重图W分解为高斯金字塔G{W}^l，l表示分解的层数，每个金字塔的分解总水平设置为最高的可能值H×W是强度梯度图像I和地形梯度图像V的空间大小；

步骤2.3，系数融合，利用公式(8)计算L{I}和L{V}的l级局部能量图局部能量图和

m,n滤波器尺度，用公式(9)计算融合模式的相似度：

相似度的范围为[-1，1]，数值越大，相似度越高，设置阈值t以确定所使用的融合模式，融合策略为公式(10)：

L{F}^l(x,y)为融合后的l层拉普拉斯金字塔，G{W}^l(x,y)为l层高斯金字塔；

步骤2.4，拉普拉斯金字塔重构；对融合后的拉普拉斯金字塔，从其顶层开始逐层从上至下按下式进行递推，恢复其对应的高斯金字塔，并得到最终的融合图像G₀；

步骤三，利用聚类分析对融合图像1和融合图像2分别进行滤波，得到融合图像1′和融合图像2′，进而执行SIFT-RANSAC算法计算融合图像1′和融合图像2′的旋转矩阵ΔR和偏移矩阵ΔT，实现粗配准；

步骤四，利用Demons算法对粗配准后的融合图像1′和融合图像2′进行精配准，得到最终得到的偏移矩阵为T_final；

步骤五，结合两幅影像的旋转矩阵ΔR和偏移矩阵ΔT和步骤四中计算的T_final对原始的卫星影像I_RS进行纠正得到I′_RS，最终I′_RS即为和ALB数据相配准的卫星影像。

2.根据权利要求1所述的基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于：步骤一中使用反距离加权插值方法对ALB点云地形数据和强度数据进行插值，得到栅格数据；

以ALB点云地形数据为例，公式(1)计算插值点(x,y)到离散点(x_i,y_i)的距离，利用公式(2)计算每个点的权重，权重是距离的倒数的函数，然后用公式(3)得到插值点的坐标得到栅格数据；

n为离散点个数，z_i为离散点z坐标，用同样的方法插值得到强度的栅格数据。

3.根据权利要求1所述的基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于：步骤一中利用公式(4)-(6)分别计算ALB点云地形和强度的梯度图像，以点云的地形数据为例，点云空间结构特征，即地形梯度图计算公式如下；

公式(4)-(6)中S_x,S_y分别为计算x、y方向导数值的两个滤波核，f为待计算的栅格数据，利用同样的方法计算强度梯度图和卫星影像不同波段的梯度图。

4.根据权利要求1所述的基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于：步骤一中利用公式(7)对ALB水下梯度特征进行调整，使其在特征图中不被忽略，以强度特征I_new为例子，I_old为待放大的特征，I_new缩放后的范围为[max_new，min_new]，_old对应原始数据的范围[_old,in_old]；

5.根据权利要求1所述的基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于：步骤三的具体实现方式如下；

步骤3.1，聚类分析，剔除干扰特征；在增强方法EnFCM框架下，对融合图像的灰度直方图进行聚类，通过迭代运算得到模糊隶属度矩阵，最后，利用滤波器对模糊隶属度矩阵进行修正，聚类灰度直方图的目标函数设计为：

u_kl表示图像灰度l相对于第k个聚类中心v_k的模糊隶属度，ξ是模态学重构后的图像，ξ_l是灰度级，q表示包含在ξ中的灰度级数量，1≤l≤q，通常远小于N，为形态闭合重构，f为输入的融合图像；然后，将优化问题即公式(12)转化为目标函数最小的无约束优化问题；

其中q和c分别表示图像中聚类中心的个数和灰度级的个数，γ_l表示灰度值为l的像素个数，表示图像灰度l相对于第k个聚类中心v_k的模糊隶属度，m为模糊参数，λ是拉格朗日乘子；因此，目标函数的最小化问题转化为寻找上面拉格朗日函数的鞍点和衍生品的拉格朗日参数，其中分割结果由模糊隶属度矩阵u_kl和聚类中心v_k共同确定，其数学表达式为：

为了得到一个稳定的模糊隶属度矩阵U＝(u_kl)^c×q，重复执行公式(15-17)直到矩阵{U^(t)-U^(t+1)}<η，t为迭代次数，其中η为最小误差阈值，通过隶属度滤波对u_kl进行修正，得到更好的隶属度划分矩阵U′，此外，还引入了中值滤波算法med以提高隶属度滤波性能和算法速度；

U″＝med{U′} (19)

根据聚类分析设置类别数，得到最终的隶属度划分矩阵U″，对干扰类别进行自动剔除，其中干扰类别的设置需要根据具体情况来确定；融合图像1和融合图像2经过上述聚类分析后得到的结果分别为融合图像1′和融合图像2′；

步骤3.2，执行SIFT-RANSAC计算两幅影像的旋转矩阵ΔR和偏移矩阵ΔT，对融合图像1′和融合图像2′进行粗配准。

6.根据权利要求1所述的基于特征融合的机载测深雷达与多光谱卫星影像配准方法，其特征在于：步骤四的具体实现方式如下；

步骤4.1，在SIFT-RANSAC配准结果的基础上对待配准图像进行格网划分，设粗配准后的融合图像1′进行格网划分后得到I_coarse1，格网框选区域内的I_coarse1为静止图像用f表示，粗配准后的融合图像2′进行格网划分后得到格网框选区域_coarse2，格网框选区域内的I_coarse2称为运动图像用m表示；

步骤4.2，依次遍历4.1中格网划分的区域，对应格网内的图像利用改进的Demons算法进行非刚性配准，执行如下操作：

第一步：模态变换，具体实现方式如下；

运动图像m的直方图图像为T_m()，静止图像直方图图像为T_f()，联合灰度直方图是直方图图像为T_mf(,f)可以通过2-D模态变换矩阵来表达；

G_m和G_f分别为m和f的最大灰度值，T_mf(,y)是m中灰度值为x，f中对应点灰度值为y的像素对的总个数，即：

T_mf(m,f)＝T_m(m)·T_f(f) (21)

而T_m()和T_f()可以通过分别对图像m和f进行灰度统计得到；通过联合直方图得到联合概率密度的计算式为：

第二步：计算互信息，具体实现方式如下；

计算运动图像m与静止图像f的互信息为：

I(m,f)＝H(m)+H(f)-H(m,f) (23)其中，

H(m,f)＝∑_m,fp(m,f)log(p(m,f)) (24)

m的熵H(m)＝∑_mp(m)log(p(m))，同理f的熵为H()，H(,)是联合熵，p(m)和p(f)是m和f中灰度值的概率，p(m,f)是图像灰度值的联合概率，将(24)带入(23)得到：

将公式(22)带入到公式(25)，计算模态变换矩阵下两个图像的互信息：

NM＝∑_m,fT_mf(m,f)是联合直方图所有元素值的总和；

第三步，计算配准误差函数，具体实现方式如下；

Demons算法的工作框架借鉴了光流算法，光流算法用于估计视频图像中相邻图像目标的位移，称为目标移动速度，为了防止图像梯度为0时计算出现问题，在传统的位移计算公式的分母中加了一项图像灰度差(m-f)，静态图像f中p点与运动图像m中对应点相匹配所需的估计位移u定义如下：

其中u＝(u_x,u_y)，表示m变换成f的位移，包含x和y不同方向的位移；是静态图像的梯度；是移动图像的梯度；配准的目的就是为了得到运动图像的位移，通过引入一个包含运动图像边缘力的方程，以提高配准收敛速度和稳定性，此外应用归一化因子α调整力的强度，则改正后的估计位移u为：

初始u设置为0，此时u的上限为1/α，因此通过设置α的值来调节形变的幅度；

定义的配准误差函数为：

S为形变位移，描述每个像素在x和y方向的位移，代表对图像的变换，U为每次迭代计算后S的偏移量，σ_i和σ_x代表强度不确定性和变换不确定性；根据当前形变参数S变换后的运动图像重新表示为在求新的u时的配准误差函数：

那么配准的过程就是一个寻找E最小的过程，计算误差梯度：

求E最小转换为则：

每一次迭代得到的变化矩阵都要进行高斯滤波，使得变化正则化，抑制噪声，防止出现错误形变，在变换过程中，新图像的边缘会出现一些缺陷，因此，模态变换后的图像在一定程度上失去了作为边缘力的能力，因此加入对误差函数E进行了修改：

为了避免局部极小值，加快配准速度，采用尺度空间方法；

步骤4.3，将遍历所有格网得到的m各像素沿x轴和y轴方向的累积偏移量矩阵∑U按照格网位置拼接得到I_coarse2所有像素点的位移T_all＝∑_格网∑U；

步骤4.4，对I_coarse2进行聚类分析得到M′；

步骤4.5，对图像M′进行二值化处理，将累积偏移量矩阵T_all映射到图像M′上，对应M′上值为1的区域对应的偏移矩阵被保留，对应M′上值为0的区域对应的偏移矩阵被删除，其中被删除的偏移信息由插值得到，最终得到的偏移矩阵为T_final，此外，为了避免局部错误配准点的干扰，采用高斯滤波算法对偏移矩阵T_final进行滤波操作。