多特征多级别的可见光全色与多光谱高精度配准方法

1.一种可见光全色图像与多光谱图像的配准方法，其特征在于，包括如下步骤：
步骤S1、对全色图像和多光谱图像进行多尺度分解，生成低分辨率的全色图像和多光谱图像；
步骤S2、在所述低分辨率的全色图像和多光谱图像上提取、匹配SIFT特征并去除外点，利用匹配的SIFT特征对和迭代重加权最小二乘法得到初始变换模型；
步骤S3、在原始图像对上利用所述初始变换模型提供的几何约束，进行基于图像块对的SIFT特征提取、匹配及外点去除，并在所有的SIFT特征对集合上选择最优的变换类型并利用迭代重加权最小二乘法求得变换参数；
步骤S4：根据所述变换模型对多光谱图像进行变换，得到变换后的多光谱图像，其中，所述步骤S3包括：
步骤S31、根据全色图像和多光谱图像的降采样倍数s1和s2以及最优前向变换模型计算对应原始分辨率全色图像和多光谱图像的变换模型
步骤S32、将全色图像分成均匀子块，对其中的每一子块P，根据初始变换模型得到对应的多光谱图像块M；在图像块P和M上分别提取SIFT特征，然后得到每对子块间的SIFT特征对，收集所有的子块间的SIFT特征对，得到对应整幅全色图像和多光谱图像的SIFT特征对集合；
pq
步骤S33、利用初始变换模型θ 提供的几何约束去除SIFT特征对集合中的外点，然后，在该集合上利用RANSAC算法去除外点，得到全局SIFT特征对集合；
步骤S34、分别以仿射变换和投影变换为最终变换模型，利用步骤S33得到的SIFT特征对，通过叠代重加权最小二乘法进行变换类型选择、求解变换参数。
2.如权利要求1所述的可见光全色图像与多光谱图像的配准方法，其特征在于，所述步骤S1包括：
步骤S11：若已知图像经纬度信息，根据经纬度信息取出全色图像和多光谱图像的重叠区域；若没有图像经纬度信息，利用两幅图像的整幅图像信息进行配准；
步骤S12：根据多光谱图像的尺寸确定多光谱图像的多尺度分解层数n1，根据全色图像尺寸及其与多光谱图像的空间分辨率差异确定全色图像的多尺度分解层数n2；
步骤S13：对全色图像和多光谱图像进行降采样，全色图像降采样倍数为 多光谱图像降采样倍数为
3.如权利要求1所述的可见光全色图像与多光谱图像的配准方法，其特征在于，所述步骤S2包括：
S21、在低分辨率全色图像和低分辨率多光谱图像上分别提取SIFT特征，得到两个SIFT特征集合；
S22、基于SIFT特征描述向量的欧式距离和最近邻比值法对两个SIFT特征集合进行匹配：
S23、将低分辨率下的全色图像和多光谱图像之间的变换类型视为仿射变换，利用RANSAC算法去除外点，得到可信度高的SIFT特征对；
S24、利用所述可信度高的SIFT特征对通过迭代重加权最小二乘法求取仿射变换参数，所述仿射变换的模型作为精配准步骤的初始变换模型。
4.如权利要求3所述的可见光全色图像与多光谱图像的配准方法，其特征在于，步骤S2的提取SIFT特征的步骤包括：
步骤S21’、构建图像的DOG金字塔；
步骤S22’、在DOG金字塔的每层上提取极值点，所谓极值点是指在局部邻域内D(x,y,kσ)取值最大的点；
步骤S23’、对于所提取的极值点，去除局部曲率不对称的极值点；
步骤S24’、计算SIFT特征的亚像素级别的空间位置、尺度，SIFT特征指的是保留下来的极值点；
步骤S25’、确定所述SIFT特征主方向，所谓主方向是指与以SIFT特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向；
步骤S26’、为SIFT特征构造描述向量，所谓描述向量是指用于刻画该SIFT特征周围图像块统计特征的、由梯度方向直方图构成的向量。
5.如权利要求3所述的可见光全色图像与多光谱图像的配准方法，其特征在于，所述步骤S2的迭代重加权最小二乘法求取初始变换参数的步骤包括：
步骤S21”、根据步骤S23的可信度高的SIFT特征对利用最小二乘法求得仿射变换模pq
型，作为所述初始变换模型θ ，采用ICP(Iterative Closest Point)方法通过前向匹配和后向匹配得到新的SIFT特征对；
步骤S22”、根据当前匹配的SIFT特征对利用迭代重加权最小二乘法重新估计前向和后向变换参数；通过ICP和参数估计的交替迭代，确定最优变换参数和最优的匹配。
6.如权利要求1所述的可见光全色图像与多光谱图像的配准方法，其特征在于，所述步骤S4包括：
对多光谱图像的进行双三次样条插值得到变换后的多光谱图像。

多特征多级别的可见光全色与多光谱高精度配准方法\n技术领域\n[0001] 本发明涉及一种用于可见光全色与多光谱图像的配准方法，尤其是航天、航空传感器平台获取的遥感影像的配准方法。\n背景技术\n[0002] 单独依靠高空间分辨率全色图像或低空间分辨率多光谱图像都无法得到目标全面、真实的信息，必须同时提供细节信息和颜色特征，人眼才能准确判断遥感图像中包含的目标。然而，由于硬件的制约，不可能同时获取高空间分辨率的多光谱图像。为此，一般的卫星都同时拍摄高空间分辨率全色图像和低空间分辨率多光谱图像。利用配准、融合技术，可以将空间分辨率全色图像和低空间分辨率多光谱图像生成高空间分辨率的多光谱图像，生成的图像既保持了全色图像的空间分辨率特性，又融入了多光谱图像的颜色信息，从而为判图人员提供了便利。然而，与手持数码细节拍摄的高清图像相比，遥感图像的质量较差，配准难度增加。对于大尺寸遥感图像，重复性结构(如建筑物)的存在给传统的特征匹配方法带来了挑战。然而，传统的图像融合方法对图像配准精度有很高的要求(一般要求配准精度优于0.5像素)。在这种情况下，研发出自动、实用的遥感图像高精度配准系统，变得尤为迫切和必要。\n[0003] 在遥感图像配准方面，国内外研究者和技术人员都进行了大量研究，并推出了一些相关产品如Erdas和Envi。这些软件在QuickBird、WorldView等国外卫星数据上可以得到很好的配准精度。然而，由于国产卫星图像质量的特殊性，这些商业软件对国产卫星的配准性能并不好。因此导致了这样的局面：一方面，国家加大投资力度，大力发展国产卫星，并将民用卫星数据免费提供给国内用户；另一方面，由于国产卫星图像配准难度很大，很多用户宁可花钱购买国外卫星数据，大量国产卫星数据被闲置。这种局面造成了国家资金、资源的巨额浪费；如果这种趋势继续持续下去，还会造成我国对国外的恶性依赖，对国防、经济、科技发展都极为不利。实际上，即使是Erdas和Envi软件，采用的也是传统的配准方法，没有考虑到大尺寸遥感图像配准的特殊要求。因此，对于大尺寸遥感图像，现有的软件处理速度都很慢；甚至因内存问题无法完成大尺寸图像的配准。因此，必须针对国产卫星数据研究具有自主知识产权的高精度配准方法，在精度、速度、自动化程度等方面满足实际应用的需要。\n发明内容\n[0004] (一)要解决的技术问题\n[0005] 本发明所要解决的技术问题是针对国产卫星图像的特殊性和现有商业软件对大尺寸图像配准的不足，提供一种自动、高效的遥感图像高精度配准方法\n[0006] (二)技术方案\n[0007] 本发明提出一种可见光全色图像与多光谱图像的配准方法，包括如下步骤：\n[0008] 步骤S1、对全色图像和多光谱图像进行多尺度分解，生成低分辨率的全色图像和多光谱图像；\n[0009] 步骤S2、在所述低分辨率图像对上提取、匹配SIFT特征并去除外点，利用匹配的SIFT特征对和迭代重加权最小二乘法得到初始变换模型；\n[0010] 步骤S3、在原始图像对上利用所述初始变换模型据提供的几何约束，进行基于图像块对的SIFT特征提取、匹配及外点去除，并在所有的SIFT特征对集合上选择最优的变换类型并利用迭代重加权最小二乘法求得变换参数；\n[0011] 步骤S4：根据所述变换模型对多光谱图像进行变换，得到变换后的多光谱图像。\n[0012] (三)有益效果\n[0013] 本发明通过基于多特征多级别的可见光全色与多光谱高精度自动配准方法，可以有效解决国产卫星图像配准问题，弥补国外商业软件对大尺寸遥感图像处理方面的缺陷。\n本发明的配准方法对国产卫星图像具有很好的通用性和实用性，将大大推动国产卫星的广泛应用，具有很好的经济效益。\n附图说明\n[0014] 图1是本发明的系统框图；\n[0015] 图2是本发明的方法的一个具体实施例的提取SIFT特征的例图；\n[0016] 图3是本发明的方法的一个具体实施例的SIFT特征提取中DOG构建图例；\n[0017] 图4是本发明的方法的一个具体实施例的SIFT特征提取中极值点提取图例；\n[0018] 图5是本发明的方法的一个具体实施例的SIFT特征提取中主方向检测图例；\n[0019] 图6是本发明的一个具体实施例的为特征点构造描述向量的一个示意图。\n具体实施方式\n[0020] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明作进一步的详细说明。\n[0021] 本发明的多特征多级别的全色与多光谱图像的高精度配准方法既可以硬件方式实现，也可以软件方式实现。例如在个人计算机、工控机及服务器上以软件的形式安装并执行，也可本发明的方法做成嵌入式芯片以硬件的形式来体现。下面参照附图对本发明的具体实施方式进行说明。\n[0022] 在下面的描述中，所述的“图像”特指通过遥感设备获得的遥感图像，并且是已进行了数字化的数字图像。然而，本发明并不限于遥感图像，对于其它技术领域中获得的需要配准的全色图像和多光谱图像，本发明也可适用。\n[0023] 图1给出了一个实施例系统框图，如图1所示，总的来说，本发明的方法包括如下步骤：\n[0024] 步骤S1：多尺度分解。对全色图像和多光谱图像进行多尺度分解，生成低分辨率的全色图像和多光谱图像；\n[0025] 根据本发明的一种具体实施方式，该实施例中的多尺度分解通过降采样实现。虽然多尺度分解可以由小波金字塔、高斯金字塔等多种方法实现，但计算量比降采样大的多；\n而在精配准阶段，SIFT(Scale Invariant Feature Transform)特征还要在原始图像上提取。因此，基于降采样的多尺度分解方式在保持配准精度的同时，还大大减少了计算量。\n[0026] 步骤S2：粗配准。在所述低分辨率图像对上提取、匹配SIFT特征并去除外点，利用匹配的SIFT特征对和迭代重加权最小二乘法得到初始变换模型；\n[0027] 根据本发明的一种具体实施方式，采用最近邻比值法来匹配SIFT特征，并去除外点，外点指不满足上述变换模型的SIFT特征对；\n[0028] 变换模型包含变换类型和变换参数两层含义。同样，变换模型的确定也包括变换类型选择和变换参数求解。根据本发明的一种具体实施方式，粗配准阶段的变换类型为仿射变换，利用迭代重加权最小二乘法得到变换参数。所述变换模型是后续配准变换模型的初始值即初始变换模型，在精配准阶段变换模型的类型及其参数都将不断调整。\n[0029] 步骤S3：精配准。在原始图像对上利用所述初始变换模型据提供的几何约束，进行基于图像块对的SIFT特征提取、匹配及外点去除，并在所有的SIFT特征对集合上选择最优的变换类型并利用迭代重加权最小二乘法求得变换参数；\n[0030] 根据本发明的一种具体实施方式，在原始图像对上采用基于图像块对的SIFT特征提取、匹配及外点去除。\n[0031] 根据本发明的一种具体实施方式，精配准的变换类型是仿射变换和投影变换中的一种，利用迭代重加权最小二乘法选择变换类型并求取变换参数。\n[0032] 步骤S4：图像变换。根据所述变换模型对多光谱图像进行变换，得到变换后的多光谱图像。\n[0033] 根据本发明的一种具体实施方式，根据变换模型对多光谱图像进行双三次样条插值，然后根据变换模型及全色图像和多光谱图像的尺寸计算参考图像和变换后的多光谱图像的最大重叠区域的坐标，取出、保存所述最大重叠区域的坐标内的全色图像和多光谱图像。\n[0034] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，对本发明作进一步的详细说明。\n[0035] 在该实施例中，多尺度分解是按照如下流程实现的：\n[0036] 步骤S11：若已知图像经纬度信息，根据经纬度信息取出全色图像和多光谱图像的重叠区域；若没有图像经纬度信息，利用两幅图像的整幅图像信息进行配准；\n[0037] 步骤S12：根据多光谱图像的尺寸确定多光谱图像的多尺度分解层数n1，根据全色图像尺寸及其与多光谱图像的空间分辨率差异确定全色图像的多尺度分解层数n2；\n[0038] 步骤S13：对全色图像和多光谱图像进行降采样，全色图像降采样倍数为多光谱图像降采样倍数为\n[0039] 本发明的步骤S2是在低分辨率全色图像和低分辨率多光谱图像上的配准。在该实施例中，粗配准是按照如下流程实现的：\n[0040] S21：特征提取。在低分辨率全色图像和低分辨率多光谱图像上分别提取SIFT特征，得到两个SIFT特征集合A和B。每个SIFT特征是132维的特征向量，前4维特征分别是特征点的x坐标、y坐标、最优尺度、主方向，后128维特征是该特征点的归一化的SIFT特征表示向量；\n[0041] S22：特征匹配。基于SIFT特征描述向量的欧式距离和最近邻比值法对两个SIFT特征集合进行匹配：对集合A中的某个SIFT特征a，它在集合B中的最近邻SIFT特征b1和次近邻SIFT特征b2的SIFT特征表示向量距离之间的欧式距离分别为t1和t2，若t1/t2＜0.8，则a和b1为SIFT特征对。\n[0042] S23：外点去除。将低分辨率下的全色图像和多光谱图像之间的变换类型视为仿射变换，利用RANSAC(RANdom SAmple Consensus，随机抽样一致性)算法去除外点，得到可信度高的SIFT特征对。\n[0043] S24：仿射变换参数求解。利用所述可信度高的SIFT特征对通过迭代重加权最小二乘法求取仿射变换参数，所述仿射变换模型作为精配准步骤的初始变换模型。\n[0044] 图2是本发明的方法的具体实施例中步骤S2的提取SIFT特征的流程图。在该实施例中，SIFT特征的提取是按照如下流程实现的：\n[0045] 步骤S21’、构建图像的DOG(Difference of Gaussian，高斯差分)金字塔。\n[0046] 设图像为I(x，y)，则DOG金字塔中第k层高斯差分图像为D(x，y，kσ)＝L(x，y，(k+1)σ)-L(x，y，kσ)，其中，L(x，y，σ)＝G(x，y，σ)*I(x，y)， *表示卷积运算。图3是DOG构建图例，如图3所示，图像在5个尺度上进行高斯卷积，得到的DOG金字塔中有4个图像。\n[0047] 步骤S22’、在DOG金字塔的每层上提取极值点，所谓极值点是指在局部邻域内D(x，y，kσ)取值最大的点。\n[0048] 提取极值点的过程为，选择DOG金字塔上的任一点，若该点在该层及上、下相邻两层的26个邻域内不是极值点，则将该点去除，否则将该点作为极值点。图4是提取极值点的示意图，如图4所示，DOG金字塔中的第k层图像的标“×”的点与第k-1层、第k+1层的\n26个标“●”的点进行比较，若标“×”的点对应的D(x，y，kσ)是这26个邻域点中的最大值，则标“×”的点为极值点。\n[0049] 步骤S23’、对于所提取的极值点，去除局部曲率非常不对称的极值点。在该实施例中，计算差分图像D的局部Hessian矩阵H，去除满足下面条件的极值点：tr(H)2/det(H)＞\n10，其中det(H)表示矩阵H的行列式值，tr(H)表示矩阵H的迹；\n[0050] 步骤S24’、计算SIFT特征的亚像素级别的空间位置、尺度，SIFT特征指的是保留下来的极值点。\n[0051] 假设SIFT特征X＝(x，y，σ)，其中x，y，σ分别为极值点X的x、y方向坐标及尺度参数，此时三个方向的坐标都为正数。根据高斯差分图像D(x，y，σ)泰勒展开公式计算出以SIFT特征X为原点相对于X在行、列及层数上的\n偏移量，即 其中D表示泰勒展开公式的一次项，D、 和 可由SIFT特征\nX及其邻域的像素根据差分法计算。然后按照如下规则对特征点进行亚像素级插值：如果三个方向上的偏移量都小于0.5像素，则该点就是特征点 为所求的亚像素级极值点坐标；如果某一方向上的偏移量大于等于0.5像素，如假设x方向的偏移量大于0.5像素，将偏移量四舍五入到一整数值a，将a与x相加得到新的SIFT特征X2＝(x+a，y，s)，然后将极值点X2按上述步骤进行操作，直到三个方向上的偏移量都小于0.5。\n[0052] 步骤S25’、确定所述SIFT特征主方向，所谓主方向是指与以SIFT特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向。\n[0053] 在该实施例中，在以SIFT特征X＝(x，y，σ)为中心、以1.5σ为半径的邻域窗口内采样，计算高斯平滑图像L(x，y，σ)在上述邻域窗口内的每一像素的梯度方向并用直方图统计邻域像素的梯度方向，得到一个梯度方向\n直方图。梯度方向直方图是一种关于梯度方向θ的统计图，其范围是0～360度，其中每\n10度一个柱，总共36个柱。该梯度方向直方图的峰值代表了该特征点处邻域梯度的主方向，即作为该SIFT特征的方向。图5是本发明的该实施例的梯度直方图的一个示例图。在该图中显示了采用7个柱为特征点确定主方向的示例。在该梯度方向直方图中，当存在另一个相当于主峰值80％能量的峰值时，则将这个方向认为是该特征点的辅方向。一个特征点可能会被指定具有多个方向(一个主方向，一个以上辅方向)。\n[0054] 步骤S26’、为SIFT特征构造描述向量，所谓描述向量是指用于刻画该SIFT特征周围图像块统计特征的、由梯度方向直方图构成的向量。\n[0055] 在该实施例中，首先将SIFT特征周围图像块的坐标轴旋转为所述SIFT特征的主方向，以确保旋转不变性；然后将特征点周围16×16像素的窗口内分成4×4像素的小块，计算每个小块的8个方向的梯度方向直方图，将每个小块的梯度方向直方图连接起来形成\n128维的特征向量；最后，将其归一化到单位长度。图6是本发明的为特征点构造描述向量的一个示意图。\n[0056] 在该实施例中，迭代重加权最小二乘法求取初始变换参数是按照如下流程实现的：\n[0057] S21”、根据步骤S23的可信度高的SIFT特征对利用最小二乘法求得仿射变换模pq\n型θ ，采用ICP(Iterative Closest Point)方法通过前向匹配和后向后向匹配得到新的SIFT特征对；\n[0058] 前向匹配以全色图像为参考图像，ICP方法为多光谱图像上的每一个SIFT特征寻找其在全色图像上对应的SIFT特征。对于多光谱图像的任一SIFT特征ci＝(xi，yi，si，pq\nθi)，它在仿射变换模型θ 下对应全色图像的坐标记为Ci，在全色图像上距离Ci最近的SIFT特征di＝(xi，yi，si，θi)为ci的匹配的SIFT特征。\n[0059] 后向匹配以多光谱图像为参考图像，ICP方法为全色图像上的每一个SIFT特征寻找其在多光谱图像上对应的SIFT特征，操作过程与前向匹配类似。\n[0060] 步骤S22”、根据当前匹配的SIFT特征对利用迭代重加权最小二乘法重新估计前向和后向变换参数；通过ICP和参数估计的交替迭代，确定最优变换参数和最优的匹配；\n[0061] 当前变换模型与当前SIFT特征对紧密相关。为此，在给定变换类型的前提下，采用迭代重加权最小二乘法不断调整变换参数及匹配的SIFT特征。假设当前变换模型为 (k表示变换类型，k＝1表示仿射变换、k＝2表示投影变换)，为方便叙述，在步骤pq\nS22”中，当变换类型给定后将变换模型 仍记为θ 。对于多光谱图像的任一SIFT特征及其在全色图像上的匹配SIFT特征 迭代重加权最小二\n乘法的目标函数为：\n[0062] \n[0063] 上式中，\n[0064] df(p，q；θpq) ＝ |(T(p；θpq)-q)Tηq|/sq， wf；i 为 每\n一SIFT特征对的权重。为简化计算，实际求解时将上述目标函数拆分为两个步骤：求解每一SIFT特征对匹配权重wd，i和更新变换参数，这两个步骤循环交叉迭代，直到算法收敛。匹配权重wd，i根据当前变换模型θpq下SIFT特征对(pi，qi)的匹配误差进行计算，其中， 在当前匹配权重wd，i的基础\n上，利用加权最小二乘法估计新通过优化 求解新的\n变换参数。σf为SIFT特征的尺度归一因子，计算公式如下：\n[0065] \n[0066] 上述算法收敛时对应的θpq和θqp为当前匹配区域的最优前向变换模型和后向变换模型。\n[0067] 本发明的步骤S3是在原始全色图像和原始多光谱图像上的配准。在该实施例中，精配准是按照如下流程实现的：\n[0068] S31：初始变换模型转换。根据全色图像和多光谱图像的降采样倍数s1和s2以及pq\n最优前向变换模型θ 计算对应原始分辨率全色图像和多光谱图像的变换变换模型 为pq\n方便叙述，在步骤S3中仍将 记为θ 。\n[0069] S32：分块特征提取、匹配。将全色图像分成1024*1024像素的均匀子块，对其中的pq\n每一子块P，根据初始变换模型θ 可以得到对应的多光谱图像块M；在图像块P和M上分别按照S21所述方法提取SIFT特征，然后按照S22所述方法得到每对子块间的SIFT特征对；收集所有的子块间的SIFT特征对，得到对应整幅全色图像和多光谱图像的SIFT特征对集合S0。\npq\n[0070] S33：外点去除。先利用初始变换模型θ 提供的几何约束去除SIFT特征对集合S0中的外点，具体步骤如下：对于S0中的某SIFT特征对(a，b)，a＝(x1，y1，s1，θ1)，b＝pq\n(x2，y2，s2，θ2)，若(x1，y1)与T(b，θ )之间的欧式距离小于τ，则保留该SIFT特征对；否则，从S0中去除SIFT特征对得到匹配点对集合S1。然后，在集合S1上利用RANSAC(RANdom SAmple Consensus，随机抽样一致性)算法去除外点，得到全局SIFT特征对集合S。T(b，pq pq\nθ )表示SIFT特征坐标(x2，y2)在变换模型θ 下的坐标。τ是一阈值，因为在步骤S34中还要通过迭代的方式重新计算每个SIFT特征对的权重，τ的取值对配准精度影响不大。\n本发明实例中，τ＝50.\n[0071] S34：变换类型选择及精确变换参数求解。分别以仿射变换和投影变换为最终变换模型，利用步骤S33得到的SIFT对通过叠代重加权最小二乘法进行变换类型选择、求解变换参数。\n[0072] 在该实施例中，变换类型选择及精确变换参数求解是按照如下流程实现的：\n[0073] 令变换类型k＝1，2，并分别进行步骤S22操作，得到对应的前向变换模型和后向变换分别为 和 当前最优匹配的前向和后向SIFT特征对集合为 和 在\nAkaike Information准则下的不适度为：\n[0074] \n[0075] 上式中， 1表示变换类型k的自由度，k＝1时，1＝6；k＝2时，\n1＝9；|A|表示集合A中元素个数。与最小的Ii对应的 和 为当前匹配区域的最优前向、后向变换模型。\n[0076] 本发明的该实施例的步骤S4是按如下流程实现的：根据步骤S3求得的变换模型，对多光谱图像的进行双三次样条插值得到变换后的多光谱图像，变换后的多光谱图像与全色图像已经配准。\n[0077] 以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。