CN104732529B - 一种遥感图像形状特征配准方法 - Google Patents

Abstract

一种遥感图像形状特征配准方法，采用对立检测背景模型，从图像中提取稳定的最大显著形状特征，完成仿射归一化，结合以配准精度为指导的误匹配去除方法控制误差。该方法基于反数学问题描述，提供了一种高精度遥感图像配准机制，对不同图像的灰度差异和波段间的场景变化不敏感，克服了背景或纹理变化带来的难题，同时具有仿射不变特性且并行执行效率高，特别适用于尺度变化和旋转条件下的多源遥感图像配准。

Description

一种遥感图像形状特征配准方法

技术领域

本发明属于图像处理领域，涉及一种遥感图像的配准方法。

背景技术

遥感图像的配准是对取自不同时间、不同传感器、或不同视角的同一场景的两幅以上图像进行信息对齐的过程。多时相或不同波段的遥感图像提供了冗余的信息描述，可综合其信息进行整合或专项产品生成，以弥补视场或时间分辨率的不足。

常用的图像配准方法大致分为两类：一类是基于区域的方法，即按区域相似性或相关性对图像进行配准，如Pearson等人最先提出的相关系数法、Barnea等人提出的序贯相似法、Shannon互信息法等。此类方法适用性窄、精度较低，但其简单易实现，在特定的少数场合使用；另一类为基于特征(点、线、面等)的方法，即利用图像中稳定可靠的特征信息，提高特征的区分能力，同时确保对光照、噪声、变换关系等干扰项的鲁棒性。其配准精度和效率较好，是当前主流的图像配准手段，包括Harris等人提出的角点及后续算法(Harris-Laplace、Hessian-Laplace、SUSAN)、Lowe提出的SIFT(Scale invariant featuretransform)及后续算法(ASIFT、PCA-SIFT)、Fourier-Mellin 功率谱函数法、Corvi等人采用的小波分析法等。此类算法采用特征附近的区域信息，形成高维的特征描述符，提高了配准的稳定性，但不具有严格的仿射不变性且计算复杂度较高。

通常，多时相的遥感图像具备以下特性：

1)多时相遥感图像背景通常发生较大的变化，灾害前后的图像更为明显。但整体的形状拓扑结构保存较好；

2)成像条件不一致，包括噪声、侧摆、光照等；

3)像素尺寸巨大，数据处理时效性矛盾突出。

相应地，已公开的配准方法对多时相，尤其是多波段遥感图像配准存在以下问题：

1)区域特征描述符依赖纹理信息，对形状明显的同态区域，如海岸和道路等地物表现不佳，适用性差；

2)理论上对仿射变换缺少不变性支持，对侧摆造成的干扰较为敏感；

3)无法进行全局配准，其计算复杂度高。

发明内容

本发明解决的技术问题是：克服现有基于区域特征描述符进行图像配准的不足，针对大幅遥感图像，提供了一种基于形状特征的并行图像配准方法，通过对子块图像采用反数学模型提取极大显著水平线，抵抗仿射变化、噪声等的影响，归一化后可以实现快速搜索，并且在全局图像变换时能够兼顾配准精度和全局几何关系，实现快速高精度的仿射不变配准。

本发明的技术解决方案是：一种遥感图像形状特征配准方法，包括如下步骤：

(1)对基准图像I_R和待配准图像I_S分别进行相同的分块，得到基准图像分块集合Sub_R和待配准图像分块集合Sub_S，Sub_R和Sub_S中的元素分别记为和i为图像分块的索引号，各图像分块在宽度方向和高度方向的长度分别记为subw和subh；

(2)遍历Sub_R和Sub_S，采用水平线描述各图像分块和中的形状信息，得到Sub_R的水平线全集以及Sub_S的水平线全集

(3)从中挑选出由最大显著水平线构成的最大显著子集从中挑选出由最大显著水平线构成的最大显著子集

(4)提取中各元素的双切线或平直片段，构成形状描述集合提取中各元素的双切线或平直片段，构成形状描述集合

(5)利用步骤(4)得到的和对和进行仿射不变归一化，得到处理后的水平线集合和

(6)遍历和寻找匹配的水平线，得到匹配形状对，上的匹配形状形成集合上的匹配形状形成集合

(7)将与中的水平线按长度进行等分，等分时以匹配形状对的仿射不变归一化后的零点为起始点，并将等分时的起始点与各分点对应的像素点作为与中的匹配点对；中的匹配点集合为中的匹配点集合为

(8)根据匹配精度要求，从和中剔除不满足匹配精度的匹配点对，得到同名点对集合P_R和P_S；

(9)利用匹配点对集合P_R和P_S拟合全局坐标变换映射关系，根据得到的全局坐标变换映射关系，对待配准图像I_S进行整体坐标变换，并通过像素插值得到配准结果。

所述步骤(2)中的水平线为区域χ_λ(u)的包络线像素的集合，

其中λ为离散化给出的数值阈值，u(x，y)表示图像。

所述步骤(3)中获取最大显著水平线与的计算方法相同，具体为：当满足关系式N_llH(min_x，y∈C|Du(x，y)|)^l/2＜ξ时，将长度为l的水平线C为作为ξ水平下的最大显著水平线C_max，其中Du(x，y)为图像u(x，y)的梯度图，ξ为控制检测虚警数的阈值，符号#表示集合的基数，N_ll为全部水平线的个数，C_max的集合即为

所述步骤(4)中获取与的方法相同，具体为：如果C_max上存在拐点，则取拐点对应的双切线作为形状描述的一部分；如果C_max为凸形或拐点较少时，则采用平直片段进行形状描述，具体操作如下：

假定C_max是n个折点组成的曲线，记其始末点弦c，主方向为像素长度为d，为C_max上相邻两点C_i、C_i+1相对c的最大方向偏转，

a)对最大显著水平线上所有的平直片段，寻找且α＜α^*的候选片段，p^*＝10^-3，α^*＝1；

b)寻找并保留候选片段中p(α，d)最小值的片段，去除与最小值片段相连的任意候选片段；

c)重复直到遍历所有候选片段，最终得到C_max的平直片段描述F_max；

F_max与双切线共同构成集合

所述步骤(5)中仿射不变归一化获取和的方法相同，具体为：

对任意最大显著水平线C_any，其上某平直片段或双切线F_any，完成：

a)C_any包围区域的质心记为G，从G作直线平行于F_any，得到的两分块质心分别记为G₁和G₃；

b)连接G₁和G₃，对G₁和G₃连线所分两分块计算其质心，分别记为G₂和G₄；

c)点G、G₁和G₂构成仿射变换基底，对C_any的仿射归一化就是寻找仿射变换实现将{G，G₂，G₁}映射到{(0，0)，(1，0)，(0，1)}；

和归一化后得到的水平线集合记为和

所述步骤(6)中寻找匹配的水平线的方法如下：当且仅当满足关系式时，确定形状与是ε水平显著匹配对，其中sup为集合的上确界，d为归一化曲线间的欧式像素距离，m为形状的集合基数，ε＝1，δ^*∈[0，1]。

所述步骤(8)中剔除不满足匹配精度的匹配点对的方法为：

a)按RANSAC从和找到仿射变换关系T；

b)将按T投影得到计算与每点像素偏差e，并剔除e＞ω的点对，取ω＝0.5～2，其中

c)重复步骤a)直到剩余每个点对都满足精度要求。

本发明与现有技术相比的优点在于：

本发明针对大幅遥感图像，基于反数学问题描述，提供了一种快速高精度遥感图像配准机制。

1.为避免传统方法对局部区域的严重依赖，本发明方法采用形状作为配准依据，以适应多时相和不同波段的精确配准，克服了背景或纹理变化带来的配准难题；

2.本发明方法根据对立检测背景模型与反数学问题描述，采用极大显著水平线作为形状描述，控制形状特征数量并对不同图像的灰度差异和波段间的场景变化不敏感，同时具有仿射不变特性；

3.本发明方法根据对立检测背景模型与反数学问题描述，指导形状描述的匹配，从源头控制误匹配的可能性；

4.本发明结合以配准精度为指导的误匹配去除方法控制配准误差，可以实现亚像素的高精度配准；

5.采用本发明方法对图像进行分块处理的并行执行效率高，特别适用于尺度变化和旋转条件下的多源大幅遥感图像配准。

附图说明

图1为本发明方法的流程框图；

图2为本发明大幅遥感图像分块并行处理示意图；

图3为本发明最大显著水平线上平直片段提取示意图；

图4为本发明平直片段仿射归一化示意图。

具体实施方式

根据Wertheimer等人提出的Gestalt不变定律，形状识别应该满足以下条件：

(1)像素间反差变化不能影响其效果；

(2)能较好地容忍视点变化带来的影响；

(3)对图像中普遍存在的噪声不敏感；

(4)对局部遮挡较为鲁棒，因此其局部表现能力较为重要；

(5)对降采样不敏感，可以理解为尺度变化。

从实验分析上说，有两种意义上的形状：任何非形变物体在不同的视角和光照条件下可以得到不同的图像，如果图像间包含可以识别的形状信息，那么说物体自身是可识别的，其形状信息构成了它的识别码(目标背景变化较小)，通常实验室环境下的近景图像序列就是这种类型；由于标准化生产的普及，通常在不同的生产线上生产同一产品，其形状非常易于识别，特别是商品的商标(目标背景变化很大)人眼非常擅长从背景变化的图像中抓取不变的分段形状信息，这在认知心理学上称为“前景背景问题”，从图像整体出发，把握其中的稳定前景或形状，那么背景中的变化就不会影响到识别的结果。

一幅图像或商标，存在它的独特形状；一个区域的遥感图像，其地物组成的拓扑结构图也成为了这个地区的“商标”，只要背景变化的幅度使得这些形状能够较为稳定的提取出来，就可以采用它们完成配准。本发明就是基于这个认识展开的。

需要说明的是，本发明中所述的图像均为[0，1]空间浮点数据。

以下结合说明书附图详细说明本发明的技术方案。图1为本发明方法的流程图，主要步骤如下：

1)为避免全局图像配准对现有通用平台的数据和内存压力，且便于形状特征均匀分布及并行计算，按图2所示对基准图像I_R和待配准图像I_S分别进行相同的分块，分块大小视数据处理能力和内存要求而定，通常取subw＝subh＝1024，这里subw和subh分别为图像块在宽度方向和高度方向的长度(用像素表示)。最终得到的分块集合记为Sub_R和Sub_S，下标R和S用于区分基准图像和待配准图像，下同；

2)首先必须确定形状基元提取方法，在图像丰富的信息中找到不变的或者稳定的形状特征。

考虑到不变性及鲁棒性要求，本发明采用水平线作为图像的形状描述。

记任意图像为u(x，y)，则λ水平下u(x，y)的水平线定义为区域χ_λ(u)的包络线：

其中λ为使用者离散化给出的数值，λ步长分的越细水平线越多，复杂度和存储压力则越大。根据实际计算，推荐λ步长为0.02。

χ_λ(u)是图像区域，包络线是水平线。水平线的集合为图像的拓扑图，是图像的完备描述且满足以下属性：水平线为Jordan曲线，不会相交且按包含关系成树状结构；水平线表示了图像的全部信息，能够用于原始图像的重构；对图像中的整体反差变化不敏感。

具体地，根据对立检测背景模型，设原始图像的噪声为加性噪声。遍历Sub_R和Sub_S，采用水平线描述子块和中的形状信息，即水平线全集和i为分块索引号(最大值为步骤1中的图像分块数量)。

3)根据对立检测理论和多重假设检验原理，从和挑选出最大显著子集和

以计算为例(计算原理相同)，定义长度为l的水平线C为ξ水平下的最大显著水平线C_max，当：

N_llH(min_x，y∈C|Du(x，y)|)^l/2＜ξ

其中，Du(x，y)为图像u(x，y)的梯度图，通常采用差分模板卷积获得。阈值ξ控制着检测的虚警数，通常不希望出现虚警，取值固定为1， 符号#表示集合的基数(集合中元素的个数)，N_ll为步骤2)获取的全部水平线的个数。

上式描述的就是多重检验中的样本与假设之间的显著性，将显著水平线检测问题归结为与背景假设偏差的小样本检测。

事实上，水平线是对比不变的。另外，可以证明，显著水平线的提取对图像间的仿射变换保持不变性，且对图像的噪声也不敏感，这对图像配准精度极其重要。

即为C_max构成的集合，同理可得

4)提取出最大显著水平线后，需要对其进行解释和描述。

双切线(同时是形状上的不同位置两个点的切线)是一种非常有效的形状描述手段，如果形状上存在有限的拐点，那么也就存在同样多的双切线。但也并不是所有的形状都存在双切线，凸形曲线不存在双切线而简单的曲线只有很少的几个，因此需要引入其他的鲁棒曲线描述。曲线上的各平直部分，记录着曲线的弯折信息，可以用作曲线的描述。

从和中提取最大显著水平线C_max的双切线或平直片段，完成形状描述，记为和

以计算为例：

如果C_max上存在拐点，那么取拐点对应的双切线作为形状描述的一部分；

当C_max为凸形或拐点较少时(程序设计时会优选使用拐点检测。凸形没有拐点；拐点数量少时，仅使用双切线会匹配失败)，采用其平直片段进行形状描述。

假定C_max是n个折点组成的曲线，记其始末点弦c，主方向(图3中的标示方向)为像素长度(起始点的像素长度)为d。定义，

为C_max上相邻两点C_i、C_i+1相对c的最大方向偏转，如图3所示。

a)对最大显著水平线上所有的平直片段，寻找且α＜α^*的候选片段。根据复杂度及数据量分析，取p^*＝10^-3，α^*＝1(其中α^*并不是非常重要，它的改变只会影响到p(α，d)的大小但不影响排序。取α^*＝1弧度简化计算过程。对p^*的选取按容许一个虚警，另外兼顾信息量与复杂度间的平衡取p^*＝10^-3)；

b)寻找并保留候选片段中p(α，d)最小值的片段，去除与最小值片段相连的任意候选片段；

c)重复直到遍历所有候选片段，最终得到C_max的平直片段描述F_max。

上述三个过程提取了反映水平线的重要曲折变化片段，作为图像中景物的不变性描述。

F_max与双切线的整体即为集合。同理可得

5)由于图像间存在未知的仿射变换关系，这些步骤中产生的形状基元不能直接进行比较。处理这个问题的一般思想是归一化，仿射归一化就是寻找形状基元使得它对平面内任意的仿射变化保持稳定。换言之，任意仿射变化下的形状在仿射不变归一化后将得到同样的结果。

根据步骤4)得到的和对和进行仿射不变归一化，如图4所示。

不失一般性，对任意最大显著水平线C_any，其上某平直片段或双切线F_any，完成：

a)C_any包围区域的质心记为G，从G作直线平行于F_any，得到的两分块质心分别记为G₁和G₃；

b)连接G₁和G₃，对G₁和G₃连线所分两分块计算其质心，分别记为G₂和G₄；

c)点G、G₁和G₂构成仿射变换基底，对C_any的仿射归一化就是寻找仿射变换实现将{G，G₂，G₁}映射到{(0，0)，(1，0)，(0，1)}。本步骤的映射就是将已经发生仿射变化的区域按质心进行归一化，映射至标准基底{(0，0)，(1，0)，(0，1)}，作为后续匹配的基础)

上述三个过程的一个基本假设是，经过仿射变换前后的区域质心形成的基底是仿射不变的。

和归一化后得到的水平线集合记为和

6)考察给定的形状S与基数为m的形状集合假定存在一种距离或相似(不相似)度量d，从中可以找出使得d(S，S′)取值很小(以下的讨论另外一层意义就是解释“小”究竟意味着什么)。

S与S′的匹配有两种情况：

1)H₀：S与S′的匹配是偶然的。例如，形状集合可能基数很大，每一条候选形状都可能在其中找到相似的形状；

2)H₁：S与S′的匹配是由于这两者确实相似性很高，例如来自于同一目标的不同图像。

对假设H₁进行充分分析，必须明确要考察的所有目标上的任一细节，这往往需要大量的先验知识，而对一些具体的形状如单个字母等则无法应用，此外这样做还需要额外的学习算法。相反，从假设H₀(背景模型)出发要实际些。

具体地，从中遍历寻找与匹配的水平线。

形状与是ε水平显著匹配对，当且仅当其中sup为集合的上确界，d为归一化曲线间的欧式像素距离，m为形状的集合基数，ε＝1(其控制的是匹配的虚警数，一般不希望出现虚警)。sup为集合的上确界，sup里面的求和符号是δ是变量，用于寻找通常复杂度较高，因此下面给出了快速性要求的设置。

如对快速性要求明显，可设置δ^*为一固定值(形状经过归一化，δ^*∈[0，1])以达到匹配次序化，采用哈希算法实现快速搜索匹配。

上述操作与最大显著水平线的寻找思路是一致的，即在多重假设检验的基础上，寻找满足匹配距离最小的匹配水平线。

与的匹配形状对记为与

7)将与的曲线按长度进行等分，按复杂度需求本发明选五等分(可以任意等分，复杂度随等分数上升)，起始点为归一化后的零点。

起始点与各分点对应的像素点即为与的匹配点对。

8)收集子块和的匹配点对，记为和

子块间变换关系近似为仿射，按照随机抽样估计(RANSAC)方法进行误匹配去除，通过配准精度要求ω迭代去除后向投影误差较大的点对。

a)按RANSAC从和找到仿射变换关系T；

b)将按T投影得到计算与每点像素偏差e，并剔除e＞ω的点对，其中

c)重复步骤a)直到剩余每个点对都满足精度要求。取ω＝0.5～2，其值越小，配准的难度越大；

9)为避免分块进行变换产生裂缝，需利用全部同名点对产生全局变换关系。具体地，按步骤2)到步骤8)遍历Sub_R和Sub_S包含的所有子块，收集全部同名点对，记为P_R和P_S。按最小二乘法得到二次多项式变换模型，使得：

P_R(x)＝a₀P_S(x)²+a₁P_S(x)+a₂P_S(x)P_S(y)+a₃

P_R(y)＝b₀P_S(y)²+b₁P_S(y)+b₂P_S(x)P_S(y)+b₃；

多项式求解可以采用开源的工具实现。

10)按二次多项式变换模型系数a₀～a₃和b₀～b₃，对待配准图像I_S进行整体坐标变换，并进行相应的像素插值，得到配准结果。

为不限定对本发明的理解，步骤9)与10)采用的变换模型可以是任意多项式模型，从工程实践中看，高阶的多项式带来的精度改善很小但计算代价太高。本发明推荐采用二阶多项式。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims (7)

1.一种遥感图像形状特征配准方法，其特征在于包括如下步骤：

(1)对基准图像I_R和待配准图像I_S分别进行相同的分块，得到基准图像分块集合Sub_R和待配准图像分块集合Sub_S，Sub_R和Sub_S中的元素分别记为和i为图像分块的索引号，各图像分块在宽度方向和高度方向的长度分别记为subw和subh；

(2)遍历Sub_R和Sub_S，采用水平线描述各图像分块和中的形状信息，得到Sub_R的水平线全集以及Sub_S的水平线全集

(3)从中挑选出由最大显著水平线构成的最大显著子集从中挑选出由最大显著水平线构成的最大显著子集

(4)提取中各元素的双切线或平直片段，构成形状描述集合提取中各元素的双切线或平直片段，构成形状描述集合

(5)利用步骤(4)得到的和对和进行仿射不变归一化，得到处理后的水平线集合和

(6)遍历和寻找匹配的水平线，得到匹配形状对，上的匹配形状形成集合上的匹配形状形成集合

(7)将与中的水平线按长度进行等分，等分时以匹配形状对的仿射不变归一化后的零点为起始点，并将等分时的起始点与各分点对应的像素点作为与中的匹配点对；中的匹配点集合为中的匹配点集合为

(8)根据匹配精度要求，从和中剔除不满足匹配精度的匹配点对，得到同名点对集合P_R和P_S；

(9)利用匹配点对集合P_R和P_S拟合全局坐标变换映射关系，根据得到的全局坐标变换映射关系，对待配准图像I_S进行整体坐标变换，并通过像素插值得到配准结果。

2.根据权利要求1所述的一种遥感图像形状特征配准方法，其特征在于： 所述步骤(2)中的水平线为区域的包络线像素的集合，

其中λ为离散化给出的数值阈值，u(x，y)表示图像。

3.根据权利要求1所述的一种遥感图像形状特征配准方法，其特征在于：所述步骤(3)中获取最大显著水平线与的计算方法相同，具体为：当满足关系式N_llH(min_x，y∈C|Du(x，y)|)^l/2＜ξ时，将长度为l的水平线C作为ξ水平下的最大显著水平线C_max，其中Du(x，y)为图像u(x，y)的梯度图，ξ为控制检测虚警数的阈值，符号#表示集合的基数，N_ll为全部水平线的个数，C_max的集合即为

4.根据权利要求3所述的一种遥感图像形状特征配准方法，其特征在于：所述步骤(4)中获取与的方法相同，具体为：如果C_max上存在拐点，则取拐点对应的双切线作为形状描述的一部分；如果C_max为凸形或拐点较少时，则采用平直片段进行形状描述，具体操作如下：

假定C_max是n个折点组成的曲线，记其始末点弦主方向为像素长度为d，为C_max上相邻两点C_i、C_i+1相对的最大方向偏转，

a)对最大显著水平线上所有的平直片段，寻找且α＜α^*的候选片段，p^*＝10^-3，α^*＝1；

b)寻找并保留候选片段中p(α，d)最小值的片段，去除与最小值片段相连的任意候选片段；

c)重复直到遍历所有候选片段，最终得到C_max的平直片段描述F_max；

F_max与双切线共同构成集合

5.根据权利要求1所述的一种遥感图像形状特征配准方法，其特征在于：所述步骤(5)中仿射不变归一化获取和的方法相同，具体为：

对任意最大显著水平线C_any，其上某平直片段或双切线F_any，完成：

a)C_any包围区域的质心记为G，从G作直线平行于F_any，得到的两分块质心 分别记为G₁和G₃；

b)连接G₁和G₃，对G₁和G₃连线所分两分块计算其质心，分别记为G₂和G₄；

c)点G、G₁和G₂构成仿射变换基底，对C_any的仿射归一化就是寻找仿射变换实现将{G，G₂，G₁}映射到{(0，0)，(1，0)，(0，1)}；

和归一化后得到的水平线集合记为和

6.根据权利要求1所述的一种遥感图像形状特征配准方法，其特征在于：所述步骤(6)中寻找匹配的水平线的方法如下：当且仅当满足关系式时，确定形状与是ε水平显著匹配对，其中sup为集合的上确界，d为归一化曲线间的欧式像素距离，m为形状的集合基数，ε＝1，δ^*∈[0，1]，是基数为m的形状集合。

7.根据权利要求1所述的一种遥感图像形状特征配准方法，其特征在于：所述步骤(8)中剔除不满足匹配精度的匹配点对的方法为：

a)按RANSAC从和找到仿射变换关系T；

b)将按T投影得到计算与每点像素偏差e，并剔除e＞ω的点对，取ω＝0.5～2，其中

c)重复步骤a)直到剩余每个点对都满足精度要求。