基于特征点几何拓扑关系的SAR图像自动配准方法与流程

本发明涉及一种图像自动配准方法，特别是涉及具有仿射变换关系的sar图像自动配准方法。

背景技术：

图像配准是一个将不同时期、不同视点或不同传感器获得的同一地域或地物的两幅或多幅图像进行线性或非线性变换，使两幅图像在几何关系上达到匹配的过程。它的工作原理是利用参考图像与待配准图像中的同名特征点对获取两幅图像之间的变换关系，从而将待配准图像的坐标系变换到参考图像的坐标系中。

变换关系的参数估计为图像自动配准的关键技术。目前图像自动配准主要有两种途径：(1)基于图像灰度的方法；(2)基于特征的方法。方法(1)直接利用灰度信息寻找图像间的变换关系，对辐射特性比较敏感，不太适合辐射失真较大的情况，更不适合斑点噪声较强的sar图像之间的配准。方法(2)首先从两幅图像中分别提取一些几何特征(如区域、曲线、直线以及角点等)作为基元，寻找基元之间的匹配关系，而后将匹配基元对转换成同名点对，最终实现变换参数估计。这类方法克服了基于灰度方法的缺点，已成为目前图像自动配准研究的主要发展方向。但此类方法大多存在几何基元提取易受图像质量影响的不足，并且基元特征提取算法的鲁棒性差、复杂度高。

鉴于上述分析，有必要研究一种针对sar图像的自动配准方法，以解决sar图像配准在非同源数据、低信噪比、场景对比度变化剧烈等实际情况下的适用性问题，提高其准确性和稳健性。

技术实现要素：

本发明针对现有sar图像自动配准技术中存在的特征点抗噪能力低以及同名点计算复杂度高等不足，提出一种基于特征点几何拓扑关系的sar图像自动配准方法。该方法提取图像中固定区域的质心作为特征点，通过特征点间的几何拓扑关系(geomatrictopologicalrelationship，gtr)寻找所有可能的完备匹配点对，利用从完备匹配点对中获得的变换参数对待配准图像进行逆变换，以参考图像与逆变换图像之间的配准精度为依据，筛选出正确的同名点对集合，并获得图像间最佳的仿射变换参数，最终实现sar图像的配准。

本发明所提出的解决方案包括以下步骤：

步骤一：提取参考图像与待配准图像的固定区域及质心

利用区域生长的方法对水域、草坪、农田以及道路等在图像上表现为平滑暗区的固定区域进行提取，提取固定区域后，区域中所有像素点的坐标均值即为区域的质心，将质心作为特征点用于后续的图像配准。

步骤二：确定具有相同仿射变换关系的所有非共线点对之间的几何拓扑关系，包括以下步骤：

2.1)分别从参考图像与待配准图像中提取特征点集并构建基准三角形点对

令p＝{p1，p2，...，pn}与q＝{q1，q2，...，qn}分别为从参考图像与待配准图像中提取的特征点集，将任意三个非共线点对(pi，qi′)、(pj，qj′)与(pk，qk′)组成基准三角形点对，其中，pi、pj与pk为点集p中的特征点，下标i、j与k代表特征点在点集中的序号，qi′、qj′与qk′分别表示点集q中的特征点，下标i′、j′与k′代表特征点在点集中的序号。

2.2)确定基准三角形点对的支持集合

基准三角形点对的支持集合为其中所有任意两个点对的支持集合的交集，即

ssii′jj′kk′＝ssii′jj′∩ssii′kk′∩ssjj′kk′(1)

其中，ssii′jj′kk′是基准三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}的支持集合，其中的任一个点对与基准三角形点对中任两个点对组成的三角形点对，均与基准三角形点对具有相同的面积比。ssii′jj′是点对(pi，qi′)与(pj，qj′)的支持集合，ssii′jj′中每一个点对与(pi，qi′)、(pj，qj′)组成的三角形点对均与基准三角形点对具有相同的面积比，ssii′kk′与ssjj′kk′分别为点对(pi，qi′)、(pk，qk′)与点对(pj，qj′)、(pk，qk′)的支持集合。

2.3)确定符合基准三角形点对确定的仿射变换关系的非共线点对之间的几何拓扑关系

若基准三角形点对的支持集合为空集，则只有该基准三角形中的三个点对符合自身确定的仿射关系；否则，从基准三角形点对与其支持集合中任取一点对，两点对的支持集合与基准三角形点对的支持集合的交集中的点均符合由基准三角形点对确定的仿射变换关系。

步骤三：基于几何拓扑关系寻找完备匹配点集合

定义完备匹配点集合为符合同一个仿射变换关系的所有点对集合。由步骤二得到具有相同仿射变换关系的特征点间的几何拓扑关系可知，若从基准三角形点对与其支持集合中任取一点对，则两点对的支持集合与基准三角形点对的支持集合的交集与基准三角形点对共同构成一个完全匹配点集，即：

其中，cmsii′jj′kk′为包含基准三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}的完备匹配点集，sshh′tt′为基准三角形点对中的任两个元素(ph，qh′)与(pt，qt′)的支持集合。因此，完备匹配点集合寻找通过以下步骤实现：

3.1)寻找匹配三角形点对集合

任意三个非共线的特征点对组成三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}，若各点对所对应的区域面积比值均与三角形面积比值相等，则该三角形点对被划入匹配三角形点对集合c中。

3.2)匹配三角形点对分组

将匹配三角形点对集合中的元素根据面积比值进行分组，比值相近的被分为同一组。分组后得到其中i＝1，2，...，n，ci中所有匹配三角形点对的面积比值都在一定的阈值范围内，n为根据实际情况划分的组数。

3.3)寻找完备匹配点集合

对于步骤3.2)得到的每一个匹配三角形点对子集，令其中每一个匹配三角形点对作为基准三角形点对，根据完备匹配点集中特征点间的几何拓扑关系，寻找包含基准三角形点对的完备匹配点集，通过以下步骤实现：

3.3.1)寻找基准三角形点对中任两个点对的支持集合

对于一组匹配三角形点对cn，其中一个基准三角形点对中的任两个点对(ps，qs′)与(pt，qt′)，其支持集合sstt′ss′可通过寻找cn中包含点对(ps，qs′)、(pt，qt′)的元素的并集获得。

3.3.2)寻找基准三角形点对的支持集合

通过3.3.1)获得基准三角形点对中任两个点对的支持集合后，可根据公式(1)得到基准三角形点对的支持集合。

3.3.3)确定完备匹配点集合

若基准三角形点对的支持集合为空集，则包含该基准三角形点对的完备匹配点集合即为其自身；若支持集合非空，则根据公式(2)计算包含基准三角形点对的完备匹配点集合。

3.4)重复步骤3.3)，直至获得全部完备匹配点集合。

步骤四：筛选确定图像间正确仿射变换关系的同名点集

虽然完备匹配点集合中所有点对均符合一个相同的仿射变换关系，但此仿射变换关系未必与图像间真实的变换关系一致。在所有的完备匹配点集合中，只有一个完备匹配点集合能够确定图像间正确的仿射变换关系，该完备匹配点集合为同名点集。

为了筛选出同名点集，获得图像间的形变关系，通过以下步骤实现：

4.1)对于每一个完备匹配点集{(pi，qi′)，(pj，qj′)，...，(pk，qk′)}，利用最小二乘法寻找一个满足最小平方距离匹配准则的仿射变换t。

4.2)对经区域提取操作后的二值待配准图像进行逆变换。

4.3)计算参考图像与逆变换图像之间的配准精度：

分别从二值参考图像与逆变换图像中提取完备匹配点集中匹配点对应的区域切片，并保证切片像素尺寸完全相同；计算所有切片的图像相似度，并取均值作为参考图像与逆变换图像的配准精度。

4.4)重复步骤4.1)、4.2)、4.3)，对全部完备匹配点集合进行相同操作，最后将具有最高配准精度的完备匹配点集合作为同名点集。

步骤五：图像配准

由同名点集确定的特征点间的映射关系t作为图像间的仿射变换关系，将待配准图像进行逆变换，即可得到与参考图像配准后的结果。

采用本发明可以达到以下技术效果：

1、本发明提取固定区域质心作为特征点，避免了sar图像中边缘检测或直线提取等操作，降低了sar图像噪声的影响；既适用于同源sar图像的配准，也适用于数据特性存在显著差异的非同源图像的配准；

2、本发明根据仿射变换下同名点间的几何拓扑关系，设计了搜索同名点集的方案，具有准确率高、复杂度低的优点。

附图说明

图1是sar图像自动配准方法的流程图；

图2是固定区域及质心提取示例图；

图3是寻找完备匹配点集的流程图；

图4是筛选同名点集的流程图；

图5是配准精度计算过程中所提图像区域切片的示例图；

图6是实测图像的配准结果图。

具体实施方式

下面结合说明书附图对本发明的实施方式进行详细说明。

本发明涉及的一种基于特征点几何拓扑关系的sar图像自动配准方法，首先提取固定区域的质心作为特征点，然后利用仿射变换下特征点的几何拓扑关系寻找完备匹配点集合，利用从完备匹配点对中获得的变换参数对待配准图像进行逆变换，以参考图像与逆变换图像之间的配准精度为依据，筛选出正确的同名点对集合，并获得图像间最佳的仿射变换参数，最终实现sar图像的配准。

图1为该方法的流程图，包含图像固定区域及其质心提取、寻找完备匹配点集合、筛选同名点集以及图像配准四个步骤。具体的操作方法如下：

步骤一：提取参考图像及待配准图像中的固定区域及其质心

由于水域、草坪、农田以及道路等固定区域散射能力弱，在sar图像中表现为平滑的暗区，即区域中各像素点的灰度值较低且差异较小，因此非常适合利用区域生长的方法实现固定区域提取。由于区域生长需要提供种子点，在不具有先验信息的条件下，本发明采用遍历的方法进行种子点搜索。提取固定区域后，区域中所有像素点的坐标均值即为区域的质心。

图2(a)与图2(c)分别为原始的参考sar图像与待配准sar图像，图2(b)与图2(d)两幅二值图像分别为利用区域生长的方法对上述两幅图像的水域提取结果，其中所提区域的像素值为0，其余像素为1。

步骤二：确定具有相同仿射变换关系的所有非共线点对之间的几何拓扑关系，包括以下步骤：

2.1)分别从参考图像与待配准图像中提取特征点集并构建基准三角形点对

令p＝{p1，p2，...，pn}与q＝{q1，q2，...，qn}分别为从参考图像与待配准图像中提取的特征点集，将任意三个非共线点对(pi，qi′)、(pj，qj′)与(pk，qk′)组成基准三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}。由于一个仿射变换关系可通过三个非共线点对确定，因此，一个基准三角形点对即对应一个仿射变换关系。

2.2)确定基准三角形点对的支持集合

基准三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}的支持集合ssii′jj′kk′中，每一个点对(ph，qh′)与基准点对中任两个点对(ps，qs′)与(pt，qt′)组成的三角形点对均为匹配三角形点对，且与基准三角形点对具有相同的面积比。因此，ssii′jj′kk′中的点对必包含于基准三角形点对中任意两点对的支持集合，即ssii′jj′kk′为基准三角形点对中任意两点对支持集合的交集，也即

ssii′jj′kk′＝ssii′jj′∩ssii′kk′∩ssjj′kk′(4)

2.3)确定符合基准三角形点对确定的仿射变换关系的非共线点对之间的几何拓扑关系

若基准三角形点对的支持集合为空集，表明只有该基准三角形中的三个点对符合自身确定的仿射关系。否则，由于点对(pi，qi′)、(pj，qj′)与(pk，qk′)非共线，可分别由(pi，pj，pk)与(qi′，qj′，qk′)组成坐标系其中表示由点pi至点pj的向量。从基准三角形点对的支持集合中任取一点对(ph，qh′)，若两点在各自坐标系下的坐标分别为(αh，βh)与(α′h，β′h)，即

由于(ph，qh′)是ssii′jj′kk′中的一个元素，因此(ph，qh′)与基准三角形点对中任两个点对组成的三角形点对与基准三角形点对相同，即

由公式(5)与公式(6)，可以得到

若点对(pt，qt′)是ssii′jj′kk′与ssii′hh′的交集中的一个元素，且pt与qt′在坐标系中的坐标分别为(αt，βt)与(α′t，β′t)，则根据三角形点对面积比值不变性，有

由公式(5)、(6)、(7)、(8)得到

αt＝α′t，βt＝β′t(9)

公式(9)表明pt与qt′在坐标系中具有完全相同的坐标，因此，点对(pt，qt′)符合由基准三角形点对确定的仿射变换关系。综上可以得到，从基准三角形点对与其支持集合中任取一点对，两点对的支持集合与基准三角形点对的支持集合的交集中的点对均符合由基准三角形点对确定的仿射变换关系。

步骤三：基于几何拓扑关系寻找完备匹配点集合

定义完备匹配点集合为符合同一个仿射变换关系的所有点对集合。由步骤二得到具有相同仿射变换关系的特征点间的几何拓扑关系可知，若从基准三角形点对与其支持集合中任取一点对，则两点对的支持集合与基准三角形点对的支持集合的交集与基准三角形点对共同构成一个完全匹配点集，即：

因此，完备匹配点集合寻找通过以下步骤实现：

令p＝{p1，p2，...，pn}与q＝{q1，q2，...，qm}分别为通过步骤一从参考图像(图2(a))与待配准图像(图2(c))中提取的特征点集。本步骤利用公式(11)，从p与q中寻找所有的完备匹配点集合，其流程图如图(3)所示，且步骤如下：

3.1)寻找匹配三角形点对集合

由于仿射变换具有面积比值不变的性质，匹配点对在各自图中对应的区域面积比值必然为一个固定值，此约束同样适用于由匹配点对作为顶点形成的区域面积比值。因此，为了获得匹配三角形点对，分别从两点集中任取一点，组成特征点对(pi，qi′)，任意三个非共线的特征点对组成三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}，将各点作为三角形顶点并计算两三角形面积比值(代表由点pi，pj与pk作为顶点组成的三角形面积，亦然)，若各点对所对应的区域面积比值均与三角形面积比值相等，则该三角形点对被划入匹配三角形点对集合c中。

3.2)匹配三角形点对分组

将匹配三角形点对集合中的元素根据面积比值进行分组，比值相近的被分为同一组。由于在实际情况下，误差的存在使得计算得到的比值不会完全相同，因此，在分组过程中可设定一个误差范围ξ(一般ξ取所有区域面积均值的3％～6％，此处ξ取为所有区域面积均值的4％)。如果三角形点对t1的面积比值rt1与三角形点对t2的面积比值rt2满足

|rt1-rt2|≤ξ(11)

则将t1与t2划分为同一组。

经上述步骤后，由步骤3.1)得到的匹配三角形点对集合c将被划分为多个子集，即

其中，ci中所有匹配三角形点对具有相同的面积比值，n为根据实际情况划分的组数。

3.3)寻找完备匹配点集合

对于步骤3.2)得到的每一个匹配三角形点对子集cn，其中n∈[1，n]，令其中每一个匹配三角形点对{(pi，qi′)，(pj，qj′)，(pk，qk′)}作为基准三角形点对，根据完备匹配点集中特征点间的几何拓扑关系，寻找包含其中三个点对的完备匹配点集cmsii′jj′kk′，包括以下步骤：

3.3.1)寻找基准三角形点对中任两个点对的支持集合

对于基准三角形点对中的任两个点对(ps，qs′)与(pt，qt′)，计算其支持集合sstt′ss′。由于sstt′ss′中的任一个点对(ph，qh′)与(ps，qs′)、(pt，qt′)组成的三角形点对均为匹配三角形点对，且与基准三角形点对具有相同的面积比，即{(ps，qs′)，(pt，qt′)，(ph，qh′)}∈cn，因此，sstt′ss′可通过寻找cn中包含点对(ps，qs′)、(pt，qt′)的元素的并集获得，即

其中，ck是包含于cn中的匹配三角形点对，‘\’代表不包含。

3.3.2)寻找基准三角形点对的支持集合

对于基准三角形点对中的任意两个点对，根据公式(13)计算其支持集合，而后利用公式(4)计算基准三角形点对的支持集合。

3.3.3)确定完备匹配点集合

若基准三角形点对的支持集合为空集，则包含该基准三角形点对的完备匹配点集合即为其自身；若支持集合非空，则利用公式(11)确定包含基准三角形点对的完备匹配点集合。

3.4)重复步骤3.3)，直至获得全部完备匹配点集合。

步骤四：筛选确定图像间正确仿射变换关系的同名点集

虽然完备匹配点集合中所有点对均符合一个相同的仿射变换关系，但此仿射变换关系未必与图像间真实的变换关系一致。在所有的完备匹配点集合中，只有一个完备匹配点集合能够确定图像间正确的仿射变换关系，该完备匹配点集合为同名点集。

为了筛选出同名点集，获得图像间的形变关系，通过以下步骤实现(图(4)所示)：

对于每一个完备匹配点集，以包含匹配点对(pi，qi′)，(pj，qj′)与(pk，qk′)的完备匹配点集cmsii′jj′kk′为例进行说明。

4.1)利用最小二乘法寻找一个满足最小平方距离匹配准则的仿射变换t；

4.2)利用仿射变换t将经区域提取后获得的二值待配准图像(如图2(d))进行逆变换；

4.3)计算参考图像与逆变换图像之间的配准精度，具体包含以下步骤：

4.3.1)对于完备匹配点集中的每一个匹配点对(pt，qt′)，分别提取其作为质心点对应的二值区域切片图像。qt′对应的区域切片图像需从经逆变换操作后的二值待配准图像(即由步骤4.2)获得的逆变换图像)中提取。提取参考图像中的区域切片时，切片四个顶点在原图中的坐标由公式(15)确定，

其中，lu、ru、ld与rd分别表示切片的左上、右上、左下以及右下四个顶点在原图中的位置。为获得一定的像素冗余，可设定冗余参数a与b，其值一般为1～2个像素，此处均取为2。(xt，yt)是质心点pt对应区域的像素点坐标，min(·)与max(·)分别为取最小值与最大值操作，而‘[]’则表示矩阵。

为了保证两幅切片图像尺寸相同，质心qt′对应的区域切片顶点在逆变换图像中的坐标则为

其中，(xt，yt)是质心点pt对应的区域像素点坐标，(xt′，yt′)是质心点qt′对应的区域在逆变换后的待配准图像中的像素点坐标，参数a与b的取值需与公式(14)相同。

4.3.2)计算完备匹配点集中所有匹配特征点对应的二值区域切片图像相似度，将其作为配准精度acc，

其中，it(xi，yi)为质心点pt对应的区域切片图像，it′(xi，yi)则为质心点qt′对应的区域局部图像，上标t表示该切片图像来源于逆变换后的二值待配准图像，m为cmsii′jj′kk′包含的点对数目。

针对图2所示的实测图像，经步骤三，可得到多组完全匹配点集。对于其中一组完全匹配点集，根据其确定的仿射变换矩阵将待配准图像的二值图像(图2(d))进行逆变换，得到图5(a)。图(b)和(c)分别为其中一组匹配点对对应的区域切片，二者的顶点在原图中的坐标分别由公式(15)与公式(16)确定，因此两幅区域切片图像具有相同的尺寸。由公式(16)得到图5(b)和图5(c)的相似度为0.7318。

4.4)重复步骤4.1)、4.2)、4.3)，对全部完备匹配点集合进行相同操作。最后将具有最高配准精度的完备匹配点集合作为同名点集。

步骤五：图像配准

由同名点集确定的特征点间的仿射变换关系t作为图像间的仿射变换关系，将待配准图像进行逆变换，即可得到与参考图像配准后的结果。

图6是利用本说明的算法将待配准图像(图6(c))与参考图像(图6(a))进行配准，得到的图像配准结果。图6(a)为参考图像，图6(b)为待配准图像逆变换后的结果，根据公式(16)估算得此次实验的配准精度约为0.989。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。