CN104867126A - 基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法 - Google Patents

Abstract

本发明公开了一种基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，主要解决不采用地面控制点数据的有变化区域的SAR图像配准问题。其实现步骤为：对于输入的浮动图像和基准图像进行滤波后分别用SURF算法提取特征点；对得到的基准图像特征点和浮动图像特征点进行归一化描述符矩阵的相似度量匹配后进行距离和连线方向角约束，得到初始候选匹配点对集；分别在浮动图像和基准图像中利用匹配特征点构建Delaunay三角形网，在两幅图像的三角形中寻找同名三角形对更新候选匹配点对集，将最小仿射误差对应的所有同名三角形对中不重复的顶点对作为最终匹配点对集；最后计算仿射变换矩阵，对浮动图像作仿射变换并双立方插值，得到最终配准结果图。

Description

基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法

技术领域

本发明属于图像处理技术领域，涉及图像配准，特别是一种基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，可用于SAR影像的变化检测、融合、拼接等方面的前期配准工作。 

背景技术

合成孔径雷达（SAR）技术以其对目标的良好识别特性及高分辨、全天候、全天时等特点在军事、环境监测等方面有着重要应用。在对SAR图像进行拼接、融合、变化检测等操作前，需将来自同一地区、在相同时间不同视点或同一传感器不同时间获取的图像在空间上进行配准，消除因获取图像的时间、角度、环境和传感器成像机理的不同造成图像间的平移、旋转、伸缩及局部形变等问题。 

1992年L.G.Brown在综述“A Survey of Image Registration Techniques”（ACM Computing Surveys,1992,vol.24,no.4,pp.325-376）中指出，所有的图像配准技术都应包含特征空间、搜索空间、搜索策略和相似性测度四个方面。特征空间指图像中用于搜索匹配的对象集合；而每一次搜索的匹配程度是以相似性测度作为评价标准的，因而相似性测度直接关系到配准的精确度；搜索空间是搜索匹配得到的所有图像变换参数所组成的空间；搜索策略则是在相似性测度下搜索最优匹配参数时所采取的优化策略。SAR图像配准的过程一般为首先确定特征空间和搜索空间，然后以特征空间为基础设计并采用某种相似性测度，按照一定的搜索策略搜索最佳匹配变换参数，并将其带入图像变换模型从而达到配准的目的。 

根据特征空间的不用，SAR配准算法一般可以分为两类：即基于灰度特征的配准和基于特征的配准。大部分基于灰度特征的方法一般直接利用整幅图像的灰度信息，通过建立某种像素间的相似性测度来衡量两幅图像重叠部分地表反射属性的匹配程度，进而寻找到最优匹配时的平移、旋转和伸缩等变换参数。在利用它们进行SAR图像的配准时，容易受到SAR图像成像原理导致的灰度差异及噪声的影响。 

基于特征的图像配准方法能够解决存在较大几何畸变及灰度信息差别的SAR图像间的配准问题。常用的图像特征有：特征点（包括角点、高曲率点等）、直线段、边缘、轮廓、闭合区域、特征结构等。基于特征的配准方法则首先需要对待配准图像提取特征点、边缘等，然后进行对应特征对或特征集之间的搜索匹配，进而得到图像变换参数。其中，基于边缘特征的配准方法需利用一些边缘检测算 子进行边缘特征提取，再对提取的边缘进行匹配，但往往由于提取的边缘过于琐碎，不能很好地体现图像中的结构特征，对后续匹配精度和速度造成影响。 

点特征相对于线、面等特征有着易于提取、速度快、精度高、稳定性好等许多方面的优势而有着更广泛的应用。基于特征点的方法基本步骤可概括为：特征点的提取、特征空间的建立、特征点的匹配、误匹配去除、变换模型选择及参数求取及图像变换与插值。常用点特征方法如Harris、SUSAN等角点，及局部不变特征算子提取的特征点等。其中局部不变特征具有或尺度不变或仿射不变的特性，对图像的缩放、旋转、仿射扭曲、光照变换等都具有较强的稳定性。 

1999年D.Lowe在文献“Object recognition from local scale-invariant features”（International Conference on Computer Vision,Corfu,Greece,pp.1150-1157）首次提出的SIFT（Scale-Invariant Feature Transform）算法，并于2004在文献“Distinctive image features from scale-invariant keypoints”（Int.J.Comput.Vis.,vol.60,no.2,pp.91–110,Nov.2004.）中进行完善。许多局部特征描述符算法都是在SIFT基础上改进的。SIFT生成图像特征的主要步骤概括为：1）尺度空间极值点提取；2）极值点定位；3）方向分配；4）关键点描述。SIFT在较小旋转角度下具有较强的旋转、尺度不变性，归一化后还对光照变化鲁棒。但由于SIFT采用了一个128维的向量来描述特征点，在特征点较多的情况下增加了运算代价，且在寻找最佳匹配时需要特征间具有较好的区分度。 

2006年Herbert Bay等在文献“SURF:Speeded Up Robust Features”（9^th European Conference on Computer Vision，ECCV2006,PartI,LNCS3951,pp.404–417）中首次提出SURF（Speeded Up Robust Features）描述子，并于2008年在文献“SURF:Speeded Up Robust Features”（International Journal on Computer Vision and Image Understanding，IJCV2008，vol.110,no.3,pp.346–359）中对SURF进行了总结和改进，将SURF的主要步骤概括为：（1）获取积分图；（2）求取近似Hessian矩阵；（3）定位特征点；（4）确定特征点的方向特征；（5）构建特征描述向量。它利用Haar小波来近似SIFT方法中的梯度操作，同时利用积分图技术进行快速计算，得到了具有64维的描述符。SURF的速度是SIFT的3-7倍，大部分情况下它和SIFT的性能相当，因此它在很多应用中得到了应用，尤其是对运行时间要求高的场合。但在旋转角度太大情况下，SURF的特征点匹配效果较差。 

通过相关算子提取特征点后，通常需要构造相应的匹配准则来匹配特征点，并尽可能地在保持正确匹配特征点对的同时去除误匹配。常用的特征点匹配准则 主要有欧式距离最近（NN）准则、最近距离比值（NNDR）准则等。其中NN匹配准则容易产生特征点间“多对一”情况；而NNDR准则对阈值参数设置敏感，有时得不到足够多的匹配点对。常用的误匹配去除方法为随机采样一致性（RANSAC），它通过对特征点进行多次迭代拟合来去除离群的特征点所在的特征点对，需要保证的初始正确匹配特征点对比例足够大，在特征点较多情况下匹配精度提升有限。 

上述多数方法仅对两幅图中的特征点进行匹配，却忽视了同一图像内的特征点间的几何分布特性。2006年XinKang等人在文献“Automatic SAR Image Registration by Using Element Triangle Invariants”（In proceedingof:9th International Conference on Information Fusion,ICIP2006，Aug.2006）中将SAR图像中的强散射点作为特征点，构建Delaunay三角形网，对每个三角形分别计算其质心，将该三角形分为3个小三角形，将每个小三角形的仿射时不变量作为该质心的一种描述，用于匹配，得到最终的匹配特征点对用于求取变换矩阵和配准图。但由于特征点位置一旦确定，对应的三角形网中的三角形也随之确定，这限制了特征点的匹配精度，达不到最佳的配准结果。 

2010年Ruirui Wang等人在文献“The normalized SIFT based on visual matching window and structural information for multi-optical imagery registration”(International Geoscience and Remote Sensing Symposium,IGARSS2010，pp.999-1002)中介绍了一种基于视觉匹配窗和结构信息的归一化SIFT特征点配准方法，对光学图像的配准得到了较好的效果。其主要思想是用三角形的仿射变换后对应边成比例的特性对特征点进行阈值约束，但它需保证有两个正确的匹配点对，然后用其来寻找满足约束的其他点对，这里阈值参数的设定需要相关经验。 

2012年Yong-MeiZhang等人在文献“A new fast multi-source remote sensing image registration algorithm”（Advanced Materials Research,v452-453,pp.950-953,2012）中介绍了一种利用SIFT提取特征点，通过Delaunay三角网去除误匹配点，并引入距离计算来决定相似控制点对的方法用于含有较少平移、旋转和背景噪声的多光谱或全彩色图的配准，提升了配准速度和精度。该方法在较大尺度和旋转情况或噪声的SAR图像配准情况下的性能仍需进一步验证。 

以上配准方法主要是针对不同照明条件、不同视角、不同场景或不同尺度下、不同数据类型的自然图像或光谱遥感图像，要求图像间具有容易提取的几何特征 或区分度较高的特征点用于匹配，对于有着大量相干斑噪声的、且不同时相间图像存在刚性的地物变化的SAR图像配准，则以上方法难以适用。常见的有变化区域的遥感图像的配准大多是基于地面控制点的配准，但是地面控制点数据的获取受各种因素影响，对于已有的多时相遥感图像常常是没有或者缺少地面控制点数据。不采用地面控制点数据的有变化区域的遥感图像配准方法： 

2010年MingtaoDing等人在文献“Registration usingrobust kernel principal component for object-based change detection”（IEEE Geoscience and Remote Sensing Letters,v7,n4,pp.761-765,October2010）中、2012年Zhaoyang Zhang等人在文献“Improved Robust Kernel Subspace for Object-Based Registration and Change Detection”（IEEE Geoscience and Remote Sensing Letters,January21,2013）中都针对多光谱图像中分类后配准算法由于分类不当导致错误控制点，使用鲁棒的核子空间方法提取图像目标的共有模式，然后利用特征空间模式匹配进行配准。 

2010年ChenboShi等人在文献“Topology based affine invariant descriptor for MSERs”（ImageProcessing(ICIP),201017th IEEE International Conference on Publication Year:2010,pp.133-136）中用椭圆拟合于最大稳定极值区域（MSER），将椭圆的仿射不变量作为区域描述，通过基于投票选择的拓扑来作为最佳匹配，用于具有视角和光照变化的多传感图像配准。2012年Xiaoming Li等人在文献“An experimental comparison of localization accuracy of affine region detectors”（Image and Signal Processing(CISP),20125th International Congress on Publication Year:2012,pp.411–414）中比较不同变化条件下(如视角、尺度、光照、散焦和图像压缩），基于检测区域重心（CGR）与仿射不变区域检测子（MSER，Hessian仿射和Harris仿射）的定位精度。 

2012年Boli Xiong等人在文献“A Method Of Acquiring Tie Point Based On Closed Regions In SAR Images”(IGARSS2012)中提出了一种基于提取的闭合区域SAR图像自动寻找控制点方法。这方法主要分三步。第一步是用测地自动轮廓模型（GAC）扣除的图像分割方法提取SAR图像闭合区域。然后，一种多边形逼近处理被用来定位这些区域边界上的特征点。根据得到的特征点，几何哈希理论用于匹配特征点成为控制点。 

2013年Al-Ruzouq等人在文献“Multiple source imagery and linear features for detection of urban expansion in Aqaba City,Jordan”（International Journal of Remote Sensing,v33,n8,pp.2563-2581,April2012）中指出现存基于特征的方法存在的不 足：（1）提出的用于匹配基元的相似度量基于经验和有时主观的；（2）包含的成像必须在自动配准前进行近似配准。线性特征（直线段）具有较高语义，能从图像中可靠提取，故可以从图像中提取的具有显著不同几何和辐射性质的线性特征作图像配准。 

从上面可以看出，不采用地面控制点数据的有变化区域的遥感图像配准方法，一般通过选择图像的几何特征如边缘、线段、轮廓等进行匹配，再将具有一定特性的奇异点作为控制点，往往定位精度不高，在处理具有较大变化区域的图像对时，有时无法找到合适的匹配模板；而通过相关算子直接在图上提取特征点用于匹配，如Harris、SUSAN等角点，SIFT、SURF等仿射不变特征点，一般基于一定的经验假设，有时特征点较多，难以找到合适的匹配点作为控制点。 

发明内容

本发明在于针对上述已有方法的不足，提出了一种基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，通过对特征点匹配准则获取的匹配点对进行约束和利用构建三角网来检测局部结构一致性，提高了特征点的分布质量，利用采样三角形的仿射特性进一步提升了配准精度，能够处理具有尺度变化和旋转的、包含丰富特征的单传感器多时相或多传感器SAR图像配准问题，适用条件为：两幅图像尺寸均在1500×1500像素以下、相对旋转角度小于90度和相对变化区域在50%以下。 

为实现上述目的，本发明的技术方案：基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于，包括下列步骤： 

（1）输入来自同一传感器对同一地区在不同时间获取的两幅SAR图像I₁和I₂，它们的大小均为M×N，每幅图像分别以其左上角作为坐标原点，为了便于描述，称图像I₁为初始浮动图，称图像I₂为初始基准图，则上述两幅图像可分别表示为I_t＝{I_t(x,y)|t＝1,2;1＜x≤M;1＜y≤N}，其中x和y分别为图像的行序号和列序号，M和N分别为图像I₁的最大行序号和最大列序号； 

（2）分别对输入的两幅图像I₁和I₂进行邻域块大小为b×b、搜索窗大小为w×w的Pretest滤波，得到浮动图像F和基准图像R，其中，邻域块大小b的取值范围为3～9，搜索窗大小w的取值范围为5～31； 

（3）对基准图像R用步骤（3a）至（3f）所述的SURF方法，计算得到N_R个特征点，及这N_R个特征点的1×N_R维的横坐标向量X_R和纵坐标向量Y_R、1×N_R 维的特征点所在尺度向量S_R、1×N_R维的主方向向量O_R，以及64×N_R维的归一化描述符矩阵D_R，具体步骤为： 

(3a)计算基准图像的积分图像； 

（3b）计算近似海森矩阵Ha； 

（3c）建立基准图像的尺度空间； 

（3d）定位特征点； 

（3e）计算特征点方向向量； 

（3f）计算特征描述向量； 

（4）对浮动图像F重复步骤（3a）至（3f）所述的SURF方法，计算得到N_F个特征点，及这N_F个特征点的1×N_F维的横坐标向量X_F和纵坐标向量Y_F、1×N_F维的特征点所在尺度向量S_F、1×N_F维的主方向向量O_F，以及64×N_F维的归一化描述符矩阵D_F； 

（5）对步骤（3）得到的基准图像R的N_R个特征点和步骤（4）得到的浮动图F的N_F个特征点进行归一化描述符矩阵的相似度量匹配，得到匹配特征点对集{P₁,P₂}，并对匹配特征点对进行距离和连线方向角约束，得到初始候选匹配点对集{P₁′,P₂′}， 

（5a）计算计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵； 

（5b）根据相似性度量矩阵和匹配准则得到匹配特征点对集{P₁,P₂}； 

（5c）计算按行并列拼接后浮动图像和基准图像中匹配特征点对间的连线距离和与水平方向的夹角； 

（5d）分别计算步骤（5c）中求得的距离和方向夹角的中心和偏离误差； 

（5e）对匹配特征点对集{P₁,P₂}进行距离和方向夹角约束得到初始候选匹配点对集{P₁′,P₂′}； 

（6）对初始候选匹配点对集{P₁′,P₂′}中的特征点矩阵P₁′和P₂′，分别在浮动图像F和基准图像R中构建Delaunay三角形网，在两幅图像的三角形网TIN₁和TIN₂中寻找来自于不同三角形网的、三角形顶点集合相同的三角形构成三角形对：若三角形对的每个三角形的三个顶点不共线，则当前三角形对构成一组同名三角形对，记录当前同名三角形对中的顶点对分别在初始候选点对集{P₁′,P₂′}中的特征点矩阵P₁′和P₂′中的序号，每个记录的序号表示在特征点矩阵P₁′和P₂′中找到的同名三 角形对的顶点对，将在初始候选匹配点对集{P₁′,P₂′}中被记录的所有同名三角形对顶点对去掉重复的点对后构成候选匹配点对集{Q₁,Q₂}； 

（6a）根据初始候选匹配点对集中的初始候选特征点矩阵分别构建浮动图像和基准图像的三角网； 

（6b）寻找同名三角形对； 

（6c）根据同名三角形对得到候选匹配点对集{Q₁,Q₂}； 

（7）从候选匹配点对集{Q₁,Q₂}中，依次选择3对匹配点对构成一个同名三角形对，若此同名三角形对中任意一个三角形的三个顶点不共线，则计算此同名三角形对的仿射误差，将最小仿射误差对应的所有同名三角形对中不重复的顶点对作为最终匹配点对集{Q₁′,Q′₂}， 

（7a）初始化仿射误差矩阵； 

（7b）计算候选匹配点对集{Q₁,Q₂}构成三角形对的仿射误差； 

（7c）重复步骤（7b）直至遍历候选匹配点对集{Q₁,Q₂}中所有的候选匹配点对； 

（7d）计算误差矩阵的最小值对应的候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}； 

（8）用最终匹配点对集{Q₁′,Q′₂}计算齐次坐标矩阵Z_R和Z_F，并计算仿射变换矩阵M，对浮动图像F进行仿射变换并作双立方插值，得到最终配准结果图。 

上述的步骤（5）所描述的初始候选匹配点对集的生成过程，具体如下： 

（5a）计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵AS＝D_F ^T×D_R，其中AS的矩阵元素AS_ij表示浮动图像F中的第i个特征点的归一化描述符向量与基准图像R中的第j个特征点的归一化描述符向量的内积，浮动图像F中第i个特征点的坐标为(x_1i,y_1i)，i=1,2,…,N_F，基准图像R中第j个特征点的坐标为j=(x_2j,y_2j)，j=1,2,…,N_R； 

（5b）根据相似性度量矩阵和匹配准则得到匹配特征点对集{P₁,P₂}，当相似性度量矩阵的元素AS_ij为第i行最小值且是第j列的最小值时，当前浮动图像F中的特征点i与基准图像R中的特征点j才是对应的，得到匹配点对{i,j}；所有N_p对相匹配的特征点的坐标构成匹配特征点对集{P₁,P₂}，其中P₁为浮动图像F中的匹配特征点坐标矩阵，该矩阵中的第i行元素为第i个匹配特征点的横坐标x_1i和纵 坐标y_1i，P₂为基准图像R中的匹配特征点坐标矩阵，该矩阵中的第i行元素为与P₁中第i个匹配特征点相匹配的基准图像R的第i个匹配特征点的横坐标x_2j和纵坐标y_2j，P₁和P₂的维数均为N_p×2； 

（5c）计算按行并列拼接后浮动图像和基准图像中匹配特征点对间的连线距离和与水平方向的夹角，沿浮动图像F的水平方向在其左侧并排拼接基准图像R，得到拼接图像I，则拼接图像的第1行的前一半为浮动图像的第1行，后一半为基准图像的第1行，其它各行均如此，则该拼接图像的总行数不变，总列数为基准图的列数的一倍。因此，在拼接图像I中，P₁中的点坐标(x_1k,y_1k)值不变，而在P₂中的点坐标(x_2k,y_2k)在拼接图像中变为(x_2k,y_2k+N)，N为浮动图像F的总列数。计算在拼接图像I上的匹配特征点间的欧式距离 并计算匹配特征点间的连线与水平方向的方向夹角得到匹配特征点对集{P₁,P₂}的距离向量D＝{d_k}和连线与水平方向的方向角向量A＝{α_k}。其中，tan^-1()表示反正切函数，k＝1,2,...,N_p，α_k的取值范围为0～π； 

（5d）分别计算（5c）中求得的距离和方向夹角的中心和偏离误差，分别统计距离向量D和方向角向量A的直方图，得到距离频数最大值对应的中心距离μ_D和方向角频数最大值对应的中心角度μ_A，计算匹配特征点对集{P₁,P₂}中每组匹配特征点对的距离d_k相对于中心距离μ_D的偏离误差计算每组匹配特征点对的方向夹角αk相对于中心夹角μ_A的偏离误差  ${\mathrm{\σ}}_A=\sqrt{{\mathrm{\Σ}}_{k=1}^{N_p}|{\mathrm{\α}}_k-{\mathrm{\μ}}_A}/N_p;$

（5e）对匹配特征点对集{P₁,P₂}进行距离和方向夹角约束得到初始候选匹配点对集{P₁′,P₂′}，对于匹配特征点对集{P₁,P₂}中的匹配特征点对(x_1k,y_1k)和(x_2k,y_2k)，当其在拼接图I中的距离d_k和方向α_k满足|d_k-μ_D|＜λ_Dσ_D且|α_k-μ_A|＜λ_Aσ_A时，保留该对特征点，否则删除，其中参数λ_D和λ_A是约束范围参数，其取值范围分别为λ_D∈[0.05,0.5]，λ_A∈[0.1,1)，本发明实例中当图像尺寸大于500×500时取λ_D=0.1和λ_A=0.1，当图像尺寸小于500×500时取λ_D＝0.2和 λ_A＝0.2；最后将剩余的特征点对的集合作为初始候选匹配点对集{P₁′,P₂′}，该集合也是按行号一一对应的。 

上述的步骤（6）所描述的候选匹配点对集的生成过程，具体如下： 

（6a）根据初始候选匹配点对集中的初始候选特征点矩阵分别构建浮动图像和基准图像的三角网，分别在浮动图像F中对初始候选匹配特征点矩阵P₁′、基准图像R中对初始候选匹配特征点矩阵P₂′按照Delaunay方法构建三角形网，得到浮动图F的t₁×3维的Delaunay三角形网矩阵TIN₁和基准图R的t₂×3维的Delaunay三角形网TIN₂，其中t₁和t₂分别表示浮动图像F和基准图像R中三角形的总数，TIN₁的每一行中的三个列元素分别代表浮动图像F中一个三角形的三个顶点在初始候选匹配点矩阵P₁′中的行序号，TIN₂的每一行中的三个列元素分别代表基准图像R中一个三角形的三个顶点在初始候选匹配点矩阵P₂′中的行序号； 

（6b）寻找同名三角形对，将TIN₁的第一行与TIN₂的每一行逐行比较，当TIN₁和TIN₂中存在相同的行序号组合时，则找到了一个同名三角形对，保留TIN₁的当前行和TIN₂中相同行序号组合的行，否则删除TIN₁的当前行； 

（6c）根据同名三角形对得到候选匹配点对集{Q₁,Q₂}，对TIN₁的每一行都执行步骤（6b）直到TIN₁的所有行都检测完毕，得到了TIN₁与TIN₂中所有三角形顶点行序号组合相同的三角形组成的同名三角形网络TIN_s；将浮动图像F的初始候选匹配点矩阵P₁′中的TIN_s中所有节点行序号构成浮动图像F中的同名匹配点矩阵Q₁，同样地，将基准图像R的初始候选匹配点矩阵P₂′中与P₁′中的行序号对应的同名匹配点矩阵Q₂，组成候选匹配点对集{Q₁,Q₂}。 

上述的步骤（6b）中所描述的同名三角形对判定准则，具体如下： 

将TIN₁的第一行与TIN₂的每一行逐行比较，当TIN₁和TIN₂中存在相同的行序号组合时，则找到了一个同名三角形对，保留TIN₁的当前行和TIN₂中相同行序号组合的行，否则删除TIN₁的当前行。 

上述的步骤（7）所描述的最终匹配点对集生成过程，具体如下： 

（7a）初始化仿射误差矩阵，构建一个大小为Nm×Nm×Nm的三维矩阵TRIerr，此三维矩阵中的每个矩阵元素初值均为1000，其中Nm为候选匹配点对集{Q₁,Q₂}的点对数目， 

在候选匹配点对集{Q₁,Q₂}中任意选择3对候选匹配点对，第1对的行序号作为三维矩阵TRIerr的第一维行序号，第2对的行序号作为三维矩阵TRIerr的第 二维列序号，第3对的行序号作为三维矩阵TRIerr的第三维页序号，则此矩阵元素的值赋值为仿射误差值，仿射误差值由下述步骤计算得到； 

（7b）计算候选匹配点对集{Q₁,Q₂}构成三角形对的仿射误差，对于候选匹配点对集{Q₁,Q₂}中每次选取三对候选匹配点对的行序号组合为(n₁,n₂,n₃)，判断所对应的浮动图像F中的候选特征点集Q1的点A=(x₁,y₁)、B=(x₂,y₂)和C=(x₃,y₃)是否满足构成三角形的条件(y₁-y₂)(x₁-x₃)≠(y₁-y₃)(x₁-x₂)，如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算对应的基准图像R中的候选特征点集Q₂的点、和是否满足构成三角形的条件如果不满足^w件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算三对候选匹配点对所形成的三角形ΔABC和ΔA′B′C′的仿射误差值  $\mathrm{TRIerr}(n_1,n_2,n_3)=\mathrm{std}\left(\right[\frac{\left|\mathrm{AB}\right|}{|A\′B\′|},\frac{\left|\mathrm{BC}\right|}{|B\′C\′|},\frac{\left|\mathrm{AC}\right|}{|A\′C\′|}\left]\right),$ 其中，|AB|、|BC|、|AC|分别表示三角形ΔABC的三条边长，|A'B'|、|B'C'|、|A'C'|分别表示与三角形ΔABC组成同名三角形对的三角形ΔA′B′C′的对应顶点组成的三条边长，符号“[]”表示元素按行排列成为一个向量，函数std()表示标准差函数，n₁＝1,...,Nm-2，n₂＝n₁+1,...,Nm-1，n₃＝n₂+1,...,Nm；n₃＝n₂+1,...,Nm； 

（7c）重复步骤（7b）直至遍历候选匹配点对集{Q₁,Q₂}中所有的候选匹配点对； 

（7d）计算误差矩阵的最小值对应的候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}。计算误差矩阵TRIerr的最小值，将所有取到最小值的候选匹配点坐标逐行存入三角形顶点矩阵TRI中，将TRI中出现的对应候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}。 

本发明的与其他技术相比具有以下优点： 

1）本发明提出了一种描述符矩阵的内积矩阵用于度量描述符相似性，并以行列交叉位置是否同时为所在行和列最小值来选择初始候选匹配点对，增加了特征的区分度和匹配点对数目，在无需设定阈值参数的同时解决了“多对一”情况。 

2）本发明通过对提取的特征点点对进行距离和连线方向的一致性约束在保留了较多正确匹配点对的同时减少大量输入特征点。 

3）利用具有仿射不变性质的Delaunay三角形网中的三角形对匹配特征点进行局部约束，较好的再次去除了大误差点对；利用三角形的仿射性质构建仿射误差矩阵来自动寻找最佳匹配仿射三角形，可以在特征点匹配精度较低的情况下找到最合适的对应点对集来求仿射变换矩阵和配准图。 

4）本发明使用SURF算法子自动提取特征点，可在硬件配置较低的情况下较快处理大幅的相同传感器拍摄的同一地区不同时相的SAR图像配准问题，对具有旋转、尺度变换的情形或包含变化区域的情形仍具有较高的配准精度。 

附图说明

图1是本发明的实现流程图。 

图2是用于实验的第一组SAR图像对，其中图2（a）是对应地区农田时相1的SAR图像，图2（b)是对应地区农田时相2的SAR图像。 

图3是用于实验的第二组SAR图像对，其中图3（a）是对应地区道路时相1的SAR图像，图3（b)是对应地区道路时相2的SAR图像。 

图4是对第一组SAR图像对采用最近邻准则得到的匹配点对连线示意图。 

图5是对第一组SAR图像对采用比阈值准则得到的匹配点对连线示意图。 

图6是对第一组SAR图像对采用本方法准则得到的初始匹配点对连线示意图。 

图7是对图6的结果进行长度约束参数为0.2和夹角约束参数为0.2进行范围约束后得到的匹配点对连线示意图。 

图8为图7的匹配点对得到的三角网，其中，图8(a)为浮动图像的特征点构造三角网，图8（b）为基准图像特的征点构造的三角网。 

图9对图8的三角网进行同名三角形对检测和顶点删除后得到的三角网，其中，图9（a）为浮动图像同名三角形网，图9（b）为基准图像同名三角形网。 

图10为几种方法对数据1进行处理得到的最终配准结果，其中图10（a）为JH1方法结果，图10（b）为DD方法结果，图10（c）为本发明结果，图10（d）为数据1的基准图。 

图11为几种方法对数据2进行处理得到的最终配准结果，其中图11（a）为JH1方法结果，图11（b）为DD方法结果，图11（c）为本发明结果，图11（d）为数据1的基准图。 

具体实施方式

参照图1，本发明的基于点对约束和三角形网的有变化区域的合成孔径雷达 图像配准方法，具体实施步骤如下： 

（1）输入在两个不同时间对同一地区由同一传感器获取的SAR图像I₁和I₂，大小均为M×N，可表示为I_t＝{I_t(x,y)|t＝1,2;1＜x≤M;1＜y≤N}，其中x和y分别为图像的行序号和列序号，每幅图像均以其左上角作为坐标原点，M和N分别为图像的最大行序号和最大列序号。 

（2）分别对图像I₁和I₂进行邻域块大小为b×b、搜索窗大小为w×w的Pretest滤波，得到浮动图像F和基准图像R，其中，邻域块大小b的取值范围为3～9，本发明实例中取值为3；搜索窗大小w的取值范围为5～31，本发明实例中取值为7。 

（3）对基准图像R用SURF方法计算得到N_R个最终特征点，及这N_R个最终特征点的1×N_R维的横坐标向量X_R和纵坐标向量Y_R、1×N_R维的最终特征点所在尺度向量S_R、1×N_R维的主方向向量O_R，以及64×N_R维的归一化描述符矩阵D_R。具体步骤如下： 

（3a）计算基准图像的积分图像。对基准图像R中像素点G=(x,y)计算其积分图像R_Σ对应位置的像素值R_Σ(G)，计算公式如下 

$R_{\mathrm{\Σ}}\left(G\right)=\underset{p=1}{\overset{p\≤x}{\mathrm{\Σ}}}\underset{q=1}{\overset{q\≤y}{\mathrm{\Σ}}}R(p,q)$

其中，R(p,q)表示基准图像R中对应像素点(p,q)的像素值。 

由于对整幅基准图像进行积分图像遍历后，基准图像中任意一个矩形区域的像素值之和可以通过加减法运算来完成，与矩形的面积无关，图像越大，节省的时间越多，这就为后续步骤中矩形滤波器（boxfilter）与积分图像的卷积运算大大减少了计算量，从而提高了运算效率。 

（3b）对步骤（3a）中得到的基准图像的积分图像R_Σ中的像素点G=(x,y)，按照如下公式计算在s尺度上的近似海森（Hessian）矩阵Ha(G,s) 

$\mathrm{Ha}(G,s)=\left[\begin{array}{cc}{D}_{\mathrm{xx}}(G,s)& D_{\mathrm{xy}}(G,s)\\ D_{\mathrm{xy}}(G,s)& D_{\mathrm{yy}}(G,s)\end{array}\right]$

其中，D_xx(G,s)、D_xy(G,s)、D_yy(G,s)为矩形滤波器与基准图像的积分图像R_Σ卷积后在像素点G=(x,y)处的卷积值。近似海森矩阵Ha(G,s)是用D_xx(G,s)、D_xy(G,s)、D_yy(G,s)替代了海森矩阵 $H(G,\mathrm{\σ})=\left[\begin{array}{cc}{L}_{\mathrm{xx}}(G,\mathrm{\σ})& L_{\mathrm{xy}}(G,\mathrm{\σ})\\ L_{\mathrm{xy}}(G,\mathrm{\σ})& L_{\mathrm{yy}}(G,\mathrm{\σ})\end{array}\right]$ 中的高斯二阶偏导与基 准图像卷积后在像素点G=(x,y)处的卷积值L_xx(G,σ)、L_xy(G,σ)、L_yy(G,σ)得到，其中式中g(σ)表示高斯核函数，σ为高斯尺度。 

（3c）建立基准图像的尺度空间。保持基准图像R大小不变，通过改变矩形滤波器的大小来对基准图像的积分图像R_Σ进行滤波，从而形成图像的尺度空间。 

矩形滤波器的初始尺度为9×9像素大小，对应的高斯尺度σ等于1.2，通过逐步放大滤波矩形的大小来计算其后的滤波矩形大小。本发明实例中，图像金字塔共分5组，每组4层，第一组的增长因子为6，所以4层的滤波矩形大小分别为9×9、15×15、21×21、27×27。第二组的增长因子为第一组增长因子的2倍，且第二组的初始值对应到第一组的第二层，所以第二组4层滤波方框大小分别为15×15、27×27、39×39、51×51。第三组的增长因子是第二组增长因子的2倍，且第三组的初始值是对应第二组的第二层，所以第三组4层滤波方框大小分别为27×27、51×51、75×75、99×99。第四组的增长因子是第三组增长因子的2倍，且第四组的初始值是对应第三组的第二层，所以第四组4层滤波方框大小分别为51×51、99×99、147×147、195×195。第五组的增长因子是第四组增长因子的2倍，且第五组的初始值是对应第四组的第二层，所以第五组4层滤波方框大小分别为99×99、195×195、291×291、387×387。将不同大小的矩形滤波器对基准图像的积分图像进行滤波就得到了尺度图像的金字塔，因而得到了由图像空间位置G=(x,y)与尺度s构成的三维尺度空间(G,s)。 

（3d）定位特征点。依据近似海森矩阵的行列式的局部极值定位特征点的位置。 

（3d1）对步骤（3c）中的三维尺度空间中的每个像素点，按照如下公式计算步骤(3b)得到的近似海森矩阵的行列式的值， 

det(Ha(G,s))=D_xx(G,s)×D_yy(G,s)-(0.9×D_xy(G,s))²； 

如果det(Ha(G,s))大于阈值T_g时，执行步骤（3d2）至步骤（3d4），小于阈值T_g的像素点舍弃不用。其中，det(·)表示计算矩阵的行列式的值，阈值T_g的取值范围为0～1，本发明实例中取T_g=0.0001。 

（3d2）计算近似海森矩阵Ha在三维尺度空间中的极值，从而确定候选兴趣点的位置和尺度。对每个处于中间两层的极值点，判断该极值点的值是否大于（或小于）其3×3×3立体邻域内的26个邻域点的值，如果是则将此极值点作为候选兴趣点，否则删除该极值点。得到各个候选兴趣点的位置G_c和对应的尺度s_c； 

（3d3）计算候选兴趣点的精确位置。由于极值点的搜索是在离散的尺度空间进行的，检测到的极值点与连续空间的极值点之间存在差别。为了得到特征点的精确位置，需要对所在尺度sc上的离散的二维函数Ha(G)进行插值，利用泰勒级数将其展开为二次拟合函数 $\mathrm{Ha}\left(G\right)\≈\mathrm{Ha}\left(0\right)+\frac{\∂\mathrm{Ha}{\left(G\right)}^T}{\∂G}G+\frac{1}{2}{G}^T\frac{{\∂}^2\mathrm{Ha}\left(G\right)}{\∂G^2}G,$ Ha(G)的导数为0时，即得到极值点所在的位置坐标为 得到候选兴趣点的精确位置其中，上标T表示矩阵的转置操作，为求偏导符号；这里通过近邻像素点的差分来近似表示偏导数，如果极值点的位置在x、y、σ的任何方向的偏移大于0.5时，意味着插值中心点已经偏移到其邻近点上，对于这样的点需要进行位置调整和再次插值，直到所有方向的偏移小于0.5，仅保留候选兴趣点集中最稳定和可重复的候选点； 

（3d4）删除低对比度的特征点和不稳定的边缘响应点。当候选兴趣点的精确位置处的极值小于阈值Tg时，该点容易受噪声的干扰而变得不稳定，则删除；对每一个候选兴趣点计算是否成立，若成立，则删除该候选兴趣点。其中，Tr(·)表示计算矩阵的迹，即求对角线元素之和，det(·)表示计算矩阵的行列式的值，r为边缘阈值参数，其取值范围为0～100，本发明实例中取值为10； 

（3d5）将保留下来的所有N_R个候选兴趣点做为最终特征点，得到基准图像R中所有最终特征点的1×N_R维的横坐标向量X_R和1×N_R维的纵坐标向量Y_R及1×N_R维的最终特征点所在尺度向量S_R。其中，为最终特征点所在尺度。 

（3e）计算最终特征点方向向量。 

（3e1）对每一个最终特征点，将其当前尺度下半径为的圆形邻域内像素进行步长为的下采样，对所有采样点计算尺度为的x和y方向上的Harr小波响应值d_x和d_y，得到该邻域内每个点对应x和y方向的响应值组成的响应向量h=(d_x,d_y)，采用以该最终特征点为中心的尺度为的高斯函数对响应向量h进行高斯加权，得到高斯加权响应向量hg； 

（3e2）对每一个最终特征点，以其为中心点，水平向右方向为0度方向，在-30度至30度方向所形成的夹角为60度的扇形区域内，计算所有采样点的高斯加 权响应向量的和hg(0)；然后，将夹角为60度的扇形区域逆时针旋转5度，即在-25度至35度方向所形成的扇形区域内，计算所有采样点的高斯加权响应向量的和hg(5)；依次将夹角为60度的扇形区域逆时针旋转5度直到旋转一周，依次得到所成区域内所有采样点的高斯加权响应向量的和hg(10)、hg(15)、hg(20)、……、hg(350)、hg(355)；将各次得到的扇形区域内所有采样点的高斯加权响应向量和响应向量和组成和向量集合{hg}，将和向量集合中幅度最大时对应的向量的方向记为该最终特征点对应的主方向。 

（3e3）对所有的最终特征点重复步骤（3e1）和步骤（3e2），得到每个最终特征点的主方向，将所有最终特征点的主方向构成1×N_R维的主方向向量O_R。 

（3f）计算特征描述向量。 

（3f1）以最终特征点为中心，将最终特征点的主方向定为垂直的y轴方向，在该最终特征点的边长为的正方形邻域内，将此正方形区域分为4×4=16个子区域。在此正方形邻域内共有400个采样点，每个子区域内有25个采样点； 

（3f2）在每个子区域中，采用尺寸为的Harr模板计算每个子区域中的25个采样像素点中每一个采样像素点的水平方向响应值d_x和垂直方向响应值d_y，采用以该采样像素点为中心的尺度为的高斯函数计算该子区域内的所有Harr小波响应值加权，得到加权水平方向响应值|d_x|和加权垂直方向响应值|d_y|，再对该子区域内的所有采样点的响应值d_x、d_y、|d_x|、|d_y|分别求和，得到该子区域的一个由Σd_x，Σ|d_x|，Σd_y，Σ|d_y|构成的1×4维的子区域响应向量。 

（3f3）对16个子区域均重复步骤（3f2），得到16个1×4维的子区域响应向量，构成一个1×64维的描述符向量，该向量具有旋转不变性和尺度不变性，再对此向量进行归一化，得到了具有光照不变性的归一化描述符向量。将基准图像R中所有N_R个最终特征点的归一化描述符向量合并，得到64×N_R维的归一化描述符矩阵D_R。 

（4）对浮动图像F重复步骤（3a）至（3f）所述的SURF方法，计算得到N_F个最终特征点，及这N_F个最终特征点的1×N_F维的横坐标向量X_F和纵坐标向量Y_F、1×N_F维的最终特征点所在尺度向量S_F、1×N_F维的主方向向量O_F，以及64×N_F维的归一化描述符矩阵D_F。 

（5）对步骤（3）得到的基准图像R的N_R个最终特征点和步骤（4）得到的浮动图F的N_F个最终特征点进行归一化描述符矩阵的相似度量匹配，得到匹配特 征点对集{P₁,P₂}，并对匹配特征点对进行距离和连线方向角约束，得到初始候选匹配点对集{P₁′,P₂′}。 

（5a）计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵AS＝D_F ^T×D_R，其中AS的矩阵元素AS_ij表示浮动图像F中的第i个特征点的归一化描述符向量与基准图像R中的第j个特征点的归一化描述符向量的内积，浮动图像F中第i个特征点的坐标为(x_1i,y_1i)，i=1,2,…,N_F，基准图像R中第j个特征点的坐标为j=(x_2j,y_2j)，j=1,2,…,N_R。 

（5b）根据相似性度量矩阵AS和匹配准则得到匹配特征点对集{P₁,P₂}。当相似性度量矩阵的元素AS_ij为第i行最小值且是第j列的最小值时，当前浮动图像F中的特征点i与基准图像R中的特征点j才是对应的，得到匹配点对{i,j}；所有N_p对相匹配的特征点的坐标构成匹配特征点对集{P₁,P₂}，其中P₁为浮动图像F中的匹配特征点坐标矩阵，该矩阵中的第i行元素为第i个匹配特征点的横坐标x_1i和纵坐标y_1i，P₂为基准图像R中的匹配特征点坐标矩阵，该矩阵中的第i行元素为与P₁中第i个匹配特征点相匹配的基准图像R的第i个匹配特征点的横坐标x_2j和纵坐标y_2j，P₁和P₂的维数均为N_p×2 

（5c）计算按行并列拼接后浮动图像和基准图像中匹配特征点对间的连线距离和与水平方向的夹角。沿浮动图像F的水平方向在其左侧并排拼接基准图像R，得到拼接图像I，则拼接图像的第1行的前一半为浮动图像的第1行，后一半为基准图像的第1行，其它各行均如此，则该拼接图像的总行数不变，总列数为基准图的列数的一倍。因此，在拼接图像I中，P₁中的点坐标值不变，而在P₂中的点坐标(x_2k,y_2k)在拼接图像中变为(x_2k,y_2k+N)，N为浮动图像F的总列数。计算在拼接图像I上的匹配特征点间的欧式距离并计算匹配特征点间的连线与水平方向的方向夹角得到匹配特征点对集{P₁,P₂}的距离向量D＝{d_k}和连线与水平方向的方向角向量A＝{α_k}。其中，tan^-1()表示反正切函数，k＝1,2,...,N_p，α_k的取值范围为0～π； 

（5d）分别计算（5c）中求得的距离和方向夹角的中心和偏离误差。分别统计距离向量D和方向角向量A的直方图，得到距离频数最大值对应的中心距离μ_D 和方向角频数最大值对应的中心角度μ_A。计算匹配特征点对集{P₁,P₂}中每组匹配特征点对的距离d_k相对于中心距离μ_D的偏离误差计算每组匹配特征点对的方向夹角α_k相对于中心夹角μ_A的偏离误差  ${\mathrm{\σ}}_A=\sqrt{{\mathrm{\Σ}}_{k=1}^{N_p}|{\mathrm{\α}}_k-{\mathrm{\μ}}_A}/N_p.$

（5e）对匹配特征点对集{P₁,P₂}进行距离和方向夹角约束得到初始候选匹配点对集{P₁′,P₂′}。对于匹配特征点对集{P₁,P₂}中的匹配特征点对(x_1k,y_1k)和(x_2k,y_2k)，当其在拼接图I中的距离d_k和方向α_k满足|d_k-μ_D|＜λ_Dσ_D且|α_k-μ_A|＜λ_Aσ_A时，保留该对特征点，否则删除，其中参数λ_D和λ_A是约束范围参数，其取值范围分别为λ_D∈[0.05,0.5]，λ_A∈[0.1,1)，本发明实例中当图像尺寸大于500×500时取λ_D=0.1和λ_A=0.1，当图像尺寸小于500×500时取λ_D＝0.2和λ_A＝0.2。最后将剩余的特征点对的集合作为初始候选匹配点对集{P₁′,P₂′}，该集合也是按行号一一对应的。 

（6）在浮动图像F和基准图像R中分别用初始候选匹配点坐标矩阵P₁′和P₂′构建Delaunay三角形网TIN₁和TIN₂，在TIN₁和TIN₂中寻找来自于不同三角形网的、三角形顶点集合相同的三角形，这些对三角形构成若干组三角形对；若一组三角形对的每个三角形的三个顶点不共线，则该组三角形对构成一组同名三角形对；记录当前同名三角形对中的顶点对分别在初始候选点对集{P₁′,P₂′}中的特征点矩阵P₁′和P₂′中的序号，每个记录的序号表示在特征点矩阵P₁′和P₂′中找到的同名三角形对的顶点对，将在初始候选匹配点对集{P₁′,P₂′}中被记录的所有同名三角形对顶点对去掉重复的点对后构成候选匹配点对集{Q₁,Q₂}； 

（6a）根据初始候选匹配点对集中的初始候选特征点矩阵分别构建浮动图像和基准图像的三角网。分别在浮动图像F中对初始候选匹配特征点矩阵P₁′、基准图像R中对初始候选匹配特征点矩阵P₂′按照Delaunay方法构建三角形网，得到浮动图F的t₁×3维的Delaunay三角形网矩阵TIN₁和基准图R的t₂×3维的Delaunay三角形网TIN₂，其中t₁和t₂分别表示浮动图像F和基准图像R中三角形的总数，TIN₁的每一行中的三个列元素分别代表浮动图像F中一个三角形的三个顶点在初始候选匹配点矩阵P₁′中的行序号，TIN₂的每一行中的三个列元素分别代表基准图 像R中一个三角形的三个顶点在初始候选匹配点矩阵P₂′中的行序号； 

（6b）寻找同名三角形对。将TIN₁的第一行与TIN₂的每一行逐行比较，当TIN₁和TIN₂中存在相同的行序号组合时，则找到了一个同名三角形对，保留TIN₁的当前行和TIN₂中相同行序号组合的行，否则删除TIN₁的当前行； 

（6c）根据同名三角形对得到候选匹配点对集。对TIN₁的每一行都执行步骤（6b）直到TIN₁的所有行都检测完毕，得到了TIN₁与TIN₂中所有三角形顶点行序号组合相同的三角形组成的同名三角形网络TIN_s；将浮动图像F的初始候选匹配点矩阵P₁′中的TIN_s中所有节点行序号构成浮动图像F中的同名匹配点矩阵Q₁，同样地，将基准图像R的初始候选匹配点矩阵P₂′中与P₁′中的行序号对应的同名匹配点矩阵Q₂，组成候选匹配点对集{Q₁,Q₂}。 

（7）从候选匹配点对集{Q₁,Q₂}中，依次选择3对匹配点对构成一个同名三角形对，若此同名三角形对中任意一个三角形的三个顶点不共线，则计算此同名三角形对的仿射误差，将最小仿射误差对应的所有同名三角形对中不重复的顶点对作为最终匹配点对集{Q₁′,Q′₂}。 

（7a）初始化仿射误差矩阵。构建一个大小为Nm×Nm×Nm的三维矩阵TRIerr，此三维矩阵中的每个矩阵元素初值均为1000，其中Nm为候选匹配点对集{Q₁,Q₂}的点对数目。 

在候选匹配点对集{Q₁,Q₂}中任意选择3对候选匹配点对，第1对的行序号作为三维矩阵TRIerr的第一维行序号，第2对的行序号作为三维矩阵TRIerr的第二维列序号，第3对的行序号作为三维矩阵TRIerr的第三维页序号，则此矩阵元素的值赋值为仿射误差值，仿射误差值由下述步骤计算得到。 

（7b）计算候选匹配点对集{Q₁,Q₂}构成三角形对的仿射误差。在候选匹配点对集{Q₁,Q₂}中每次选取三对候选匹配点对的行序号组合为(n₁,n₂,n₃)，判断所对应的浮动图像F中的候选特征点矩阵Q₁的点A=(x₁,y₁)、B=(x₂,y₂)和C=(x₃,y₃)是否满足构成三角形的条件(y₁-y₂)(x₁-x₃)≠(y₁-y₃)(x₁-x₂)，如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算对应的基准图像R中的候选特征点矩阵Q₂的点、和是否满足构成三角形的条件如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算三对候选匹配点对所形成的三角形ΔABC和ΔA′B′C′的仿射误差值  $\mathrm{TRIerr}(n_1,n_2,n_3)=\mathrm{std}\left(\right[\frac{\left|\mathrm{AB}\right|}{|A\′B\′|},\frac{\left|\mathrm{BC}\right|}{|B\′C\′|},\frac{\left|\mathrm{AC}\right|}{|A\′C\′|}\left]\right),$ 其中，|AB|、|BC|、|AC|分别表示三角形ΔABC的三条边长，|A'B'|、|B'C'|、|A'C'|分别表示与三角形ΔABC组成同名三角形对的三角形ΔA′B′C′的对应顶点组成的三条边长，符号“[]”表示元素按行排列成为一个向量，函数std()表示标准差函数，n₁＝1,...,Nm-2，n₂＝n₁+1,...,Nm-1，n₃＝n₂+1,...,Nm； 

（7c）重复步骤（7b）直至遍历候选匹配点对集{Q₁,Q₂}中所有的候选匹配点对。 

（7d）计算误差矩阵的最小值对应的候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集。计算误差矩阵TRIerr的最小值，将所有取到最小值的候选匹配点坐标逐行存入三角形顶点矩阵TRI中，将TRI中出现的对应候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}。 

（8）用最终匹配点对集{Q₁′,Q′₂}计算齐次坐标矩阵Z_R和Z_F，并计算仿射变换矩阵M，对浮动图像F作仿射变换并双立方插值，得到最终配准结果图。设基准图像R中的坐标(x,y)经仿射变换M后在浮动图像F中的坐标为故步骤（7b）中ΔABC与ΔA'B'C'的顶点的仿射变换模型可以用如下的矩阵形式来描述  $M=\left[\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ 0& 0& 1\end{array}\right],$ 其中，参数m₁₁,m₁₂,m₂₁,m₂₂表示尺度缩放和旋转；m₁₃为水平位移，m₂₃为垂直位移，，故可利用最终匹配点对集来求仿射变换矩阵M。 

（8a）计算齐次坐标矩阵Z_R和Z_F。分别将最终匹配点对集{Q₁′,Q′₂}中的最终匹配特征点坐标矩阵Q₁′和Q′₂均由原来的两列扩展为三列，其中新增的第三列的值均置为1，得到齐次坐标矩阵 $Z_R=\left[\begin{array}{ccc}{x}_A& x_B& x_C\\ y_A& y_B& y_C\\ 1& 1& 1\end{array}\right]$ 和 $Z_F=\left[\begin{array}{ccc}{\hat{x}}_A& {\hat{x}}_B& {\hat{x}}_C\\ {\hat{y}}_A& {\hat{y}}_B& {\hat{y}}_C\\ 1& 1& 1\end{array}\right].$

（8b）计算M＝Z_R ^T/Z_F ^T得到仿射变换矩阵M，其中上标T表示对矩阵进行转置，符号“/”表示矩阵的除运算。 

（8c）对浮动图像F按照仿射变换矩阵M进行图像仿射变换和双立方插值得到浮动图像I₁相对于基准图像I₂的最终配准结果图。 

本算法的效果可以通过如下仿真试验进一步说明。 

1）仿真条件： 

硬件平台为：IntelCore2DuoCPUE65502.3GHz、2GB RAM Win7OS 

软件平台为：Matlab2012b 

2）所用实验数据 

数据1：分别从ALOS PALSAR卫星于2000年6月18日和2009年6月19对黄河入海口地区拍摄的两幅一景大图中对应地区农田图像，如图2所示，其中图2（a）为时相1图像，图2（b）为时相2图像，图像尺寸均为400*400，显著特征为农田的用地变化，具有明显的旋转和对比度差异。 

数据2：分别从ALOS PALSAR卫星于2000年6月18日和2009年6月19对黄河入海口地区拍摄的两幅一景大图中对应地区道路图像，如图3（a）为时相1图像，图3(b)为时相2图像，图像尺寸均为1000*1000，显著特点为道路和耕地的纹理丰富，具有与数据1一样的旋转参数。 

在本发明中，除非特别说明，上述数据均以时相1图像作为浮动图像，时相2图像作为基准图像。 

3）对比试验 

为了说明本发明的有效性，本发明与如下两个对比方法进行了对比： 

对比方法1：采用本章提出的匹配准则，仅作范围约束得到控制点的方法，对保留的控制点构造三角网后直接寻找误差最小的三角形的方法记为JH1； 

对比方法2：采用比阈值准则（0.7）得到初始匹配后，进行顶点三角数目统计，去除数目不一致的顶点后保留的点作为控制点的方法，记为DD； 

4）评价标准 

应用本算法对给定的两组数据中的两幅多时相图进行试验，应用最终的特征点数目、最终的配准图像与基准图的比较及变换矩阵参数等作评价。 

（1）匹配点数与同名点数的差距越小，说明特征点对的正确对应越多。 

（2）设点集P₁＝{(x_l,y_l)|l＝1,2,...,LN}与是一对匹配点对集，它们是一一对应的，其中点对数目为LN。根据下面的变换关系 

$M=\left[\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ 0& 0& 1\end{array}\right],\left[\begin{array}{c}{\hat{x}}_l\\ {\hat{y}}_1\\ 1\end{array}\right]=M\×\left[\begin{array}{c}{x}_l\\ y_l\\ 1\end{array}\right],l=\mathrm{1,2},...,\mathrm{LN}$

可由下式计算那么P₁到P₂变换误差 

$\mathrm{RMSE}=\sqrt{\frac{\underset{l=1}{\overset{\mathrm{LN}}{\mathrm{\Σ}}}[{({\hat{x}}_l-m_{11}{x}_l-m_{12}{y}_l-m_{13})}^2+{({\hat{y}}_l-m_{21}{x}_l-m_{22}{y}_l-m_{23})}^2}{\mathrm{LN}}}$

RMSE的值越小说明点对集的仿射变换一致性越好。 

（3）由于数据2和数据3是来自同一对数据源，故可相互作为参考进行比较。若同一算法对数据2和数据3分别处理得到的最终变换矩阵的旋转参数分别为{m₁₁,m₁₂,m₂₁,m₂₂}和那它们之间旋转参数方差计算公式  $\mathrm{\σ}=\sqrt{{({\hat{m}}_{11}-m_{11})}^2+{({\hat{m}}_{12}-m_{12})}^2+{({\hat{m}}_{21}-m_{21})}^2+{({\hat{m}}_{22}-m_{22})}^2}.$ 其值越小说明算法的性能越好。 

仿真内容 

仿真1匹配准则的对比 

对数据1采用SURF算法提取特征点后，分别按照最近邻准则（NN）、比阈值准则（阈值取0.7）和本发明的匹配准则对比，实验结果如图4到图7所示。 

从图4到图7可知：NN准则得到了具有歧义性的匹配点对，即“多对一”情况；而NSR准则在获得匹配点对较少的同时，正确匹配率也较低，具体表现在图5中得到的11组匹配点对，其中只有3对为正确匹配；而本发明提出的匹配准则，在初始匹配阶段就获得了大量的候选匹配对，经过范围约束后在保证足够的正确匹配对的同时去除了大量误匹配特征点对。 

仿真2同名三角网示例 

对数据1采用本发明方法处理，得到的初始三角形网和同名三角形网如图8到图9所示。从中可以看出本发明提出的方法再次剔除了部分误匹配点对，使浮动图和基准图中的特征点的虚拟几何结构趋于一致。 

仿真3算法性能对比 

对数据1和数据2分别利用JH1方法、DD方法和本发明的方法进行处理，其实验结果如图10到图11所示。具体试验数据如表1及表2所示。 

表1、匹配点对数目、变换矩阵参数及其RMSE、运行时间等实验数据对比 

表2旋转参数误差对比 

	JH1	DD	ours
				ΔM11	6.96	-0.02	0.11
ΔM12	3.41	0.2	-0.05
				ΔM21	-14.97	0.33	0.07
ΔM22	7.16	0.04	-0.13
				σ	18.31	0.39	0.19

由图10、图11可知 

（1）相对于其他算法，本发明中的通过寻找同名三角网络构建候选点对集和仿射误差最小三角形顶点组合，得到了较高精度的配准结果。 

（2）通过寻找同名网络，能够有效去除一些分布仿射性质较差的节点对应的点对，可以更具有较好的区分相似三角和仿射三角形，提升了最终的配准精度。 

由表1、表2可知： 

（1）由表1可知，相对于最近邻与次近邻比阈值准则，本发明的准则能获得更多的初始匹配点对，这就会加大寻找数量较多的同名点对和候选匹配点对的可能，从而提升最终匹配点对的匹配精度。 

（2）数据2中JH1方法出现了较大的误差的原因是未经过同名顶点确认，不能有效去除位置偏移较大的匹配点对，仅使用原始的不规则三角网中的三角形而未进行采样三角形不能克服三角形顶点定位精度的限制来进一步提升匹配精度，从而影响配准精度；而DD方法仅比较三角网中每个顶点的三角形数目，同样未能很好地解决匹配精度不高的问题；本发明的方法较好地解决了这一问题，可以突破三角网的限制来进一步去除匹配精度较差的匹配点对。 

（3）根据表2，通过对相同方法对数据1和数据2进行处理，比较其仿射变换矩阵M,的旋转参数变化可知，本发明对来自同一数据源的两组数据分别进行处理得到的旋转参数误差是所列三种方法中最小的，这进一步验证了本发明在处理具有旋转和变化区域的图像对时，仍具有良好的性能。 

综上所述，本发明中的候选点对集获取方法得到了较多的匹配点对；通过Delaunay三角形网可以很好地刻画一幅图像中的特征点分布情况；而通过本发明可以至少找到一组三角形对仿射误差最小的一对三角形来求仿射变换，具有比对比方法更高的配准精度。 

Claims (7)

1.基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于，包括下列步骤：

（1）输入来自同一传感器对同一地区在不同时间获取的两幅SAR图像I₁和I₂，它们的大小均为M×N，每幅图像分别以其左上角作为坐标原点，为了便于描述，称图像I₁为初始浮动图，称图像I₂为初始基准图，则上述两幅图像可分别表示为I_t＝{I_t(x,y)|t＝1,2;1＜x≤M;1＜y≤N}，其中x和y分别为图像的行序号和列序号，M和N分别为图像I₁的最大行序号和最大列序号；

（2）分别对输入的两幅图像I₁和I₂进行邻域块大小为b×b、搜索窗大小为w×w的Pretest滤波，得到浮动图像F和基准图像R，其中，邻域块大小b的取值范围为3～9，搜索窗大小w的取值范围为5～31；

（3）对基准图像R用步骤（3a）至（3f）所述的SURF方法，计算得到N_R个特征点，及这N_R个特征点的1×N_R维的横坐标向量X_R和纵坐标向量Y_R、1×N_R维的特征点所在尺度向量S_R、1×N_R维的主方向向量O_R，以及64×N_R维的归一化描述符矩阵D_R，具体步骤为：

（3a）计算基准图像的积分图像；

（3b）计算近似海森矩阵Ha；

（3c）建立基准图像的尺度空间；

（3d）定位特征点；

（3e）计算特征点方向向量；

（3f）计算特征描述向量；

（4）对浮动图像F重复步骤（3a）至（3f）所述的SURF方法，计算得到N_F个特征点，及这N_F个特征点的1×N_F维的横坐标向量X_F和纵坐标向量Y_F、1×N_F维的特征点所在尺度向量S_F、1×N_F维的主方向向量O_F，以及64×N_F维的归一化描述符矩阵D_F；

（5）对步骤（3）得到的基准图像R的N_R个特征点和步骤（4）得到的浮动图F的N_F个特征点进行归一化描述符矩阵的相似度量匹配，得到匹配特征点对集{P₁,P₂}，并对匹配特征点对进行距离和连线方向角约束，得到初始候选匹配点对集{P₁′,P₂′}，

（5a）计算计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵；

（5b）根据相似性度量矩阵和匹配准则得到匹配特征点对集{P₁,P₂}；

（5c）计算按行并列拼接后浮动图像和基准图像中匹配特征点对间的连线距离和与水平方向的夹角；

（5d）分别计算步骤（5c）中求得的距离和方向夹角的中心和偏离误差；

（5e）对匹配特征点对集{P₁,P₂}进行距离和方向夹角约束得到初始候选匹配点对集{P₁′,P₂′}；

（6）对初始候选匹配点对集{P₁′,P₂′}中的特征点矩阵P₁′和P₂′，分别在浮动图像F和基准图像R中构建Delaunay三角形网，在两幅图像的三角形网TIN₁和TIN₂中寻找来自于不同三角形网的、三角形顶点集合相同的三角形构成三角形对：若三角形对的每个三角形的三个顶点不共线，则当前三角形对构成一组同名三角形对，记录当前同名三角形对中的顶点对分别在初始候选点对集{P₁′,P₂′}中的特征点矩阵P₁′和P₂′中的序号，每个记录的序号表示在特征点矩阵P₁′和P₂′中找到的同名三角形对的顶点对，将在初始候选匹配点对集{P₁′,P₂′}中被记录的所有同名三角形对顶点对去掉重复的点对后构成候选匹配点对集{Q₁,Q₂}；

（6a）根据初始候选匹配点对集中的初始候选特征点矩阵分别构建浮动图像和基准图像的三角网；

（6b）寻找同名三角形对；

（6c）根据同名三角形对得到候选匹配点对集{Q₁,Q₂}；

（7）从候选匹配点对集{Q₁,Q₂}中，依次选择3对匹配点对构成一个同名三角形对，若此同名三角形对中任意一个三角形的三个顶点不共线，则计算此同名三角形对的仿射误差，将最小仿射误差对应的所有同名三角形对中不重复的顶点对作为最终匹配点对集{Q₁′,Q′₂}，

（7a）初始化仿射误差矩阵；

（7b）计算候选匹配点对集{Q₁,Q₂}构成三角形对的仿射误差；

（7c）重复步骤（7b）直至遍历候选匹配点对集{Q₁,Q₂}中所有的候选匹配点对；

（7d）计算误差矩阵的最小值对应的候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}；

（8）用最终匹配点对集{Q₁′,Q′₂}计算齐次坐标矩阵Z_R和Z_F，并计算仿射变换矩阵M，对浮动图像F进行仿射变换并作双立方插值，得到最终配准结果图。

2.根据权利要求1所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于步骤（5）所描述的初始候选匹配点对集的生成过程，具体如下：

（5a）计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵AS＝D_F ^T×D_R，其中AS的矩阵元素AS_ij表示浮动图像F中的第i个特征点的归一化描述符向量与基准图像R中的第j个特征点的归一化描述符向量的内积，浮动图像F中第i个特征点的坐标为(x_1i,y_1i)，i=1,2,…,N_F，基准图像R中第j个特征点的坐标为j=(x_2j,y_2j)，j=1,2,…,N_R；

（5b）根据相似性度量矩阵和匹配准则得到匹配特征点对集{P₁,P₂}，当相似性度量矩阵的元素AS_ij为第i行最小值且是第j列的最小值时，当前浮动图像F中的特征点i与基准图像R中的特征点j才是对应的，得到匹配点对{i,j}；所有N_p对相匹配的特征点的坐标构成匹配特征点对集{P₁,P₂}，其中P₁为浮动图像F中的匹配特征点坐标矩阵，该矩阵中的第i行元素为第i个匹配特征点的横坐标x_1i和纵坐标y_1i，P₂为基准图像R中的匹配特征点坐标矩阵，该矩阵中的第i行元素为与P₁中第i个匹配特征点相匹配的基准图像R的第i个匹配特征点的横坐标x_2j和纵坐标y_2j，P₁和P₂的维数均为N_p×2；

（5c）计算按行并列拼接后浮动图像和基准图像中匹配特征点对间的连线距离和与水平方向的夹角，沿浮动图像F的水平方向在其左侧并排拼接基准图像R，得到拼接图像I，则拼接图像的第1行的前一半为浮动图像的第1行，后一半为基准图像的第1行，其它各行均如此，则该拼接图像的总行数不变，总列数为基准图的列数的一倍。因此，在拼接图像I中，P1中的点坐标(x_1k,y_1k)值不变，而在P₂中的点坐标(x_2k,y_2k)在拼接图像中变为(x_2k,y_2k+N)，N为浮动图像F的总列数。计算在拼接图像I上的匹配特征点间的欧式距离并计算匹配特征点间的连线与水平方向的方向夹角得到匹配特征点对集{P₁,P₂}的距离向量D＝{d_k}和连线与水平方向的方向角向量A＝{α_k}。其中，tan^-1()表示反正切函数，k＝1,2,...,N_p，α_k的取值范围为0～π；

（5d）分别计算（5c）中求得的距离和方向夹角的中心和偏离误差，分别统计距离向量D和方向角向量A的直方图，得到距离频数最大值对应的中心距离μ_D和方向角频数最大值对应的中心角度μ_A，计算匹配特征点对集{P₁,P₂}中每组匹配特征点对的距离dk相对于中心距离μ_D的偏离误差计算每组匹配特征点对的方向夹角αk相对于中心夹角μ_A的偏离误差 ${\mathrm{\σ}}_A=\sqrt{{\mathrm{\Σ}}_{k=1}^{N_p}|{\mathrm{\α}}_k-{\mathrm{\μ}}_A}/N_p;$

（5e）对匹配特征点对集{P₁,P₂}进行距离和方向夹角约束得到初始候选匹配点对集{P₁′,P₂′}，对于匹配特征点对集{P₁,P₂}中的匹配特征点对(x_1k,y_1k)和(x_2k,y_2k)，当其在拼接图I中的距离d_k和方向α_k满足|d_k-μ_D|＜λ_Dσ_D且|α_k-μ_A|＜λ_Aσ_A时，保留该对特征点，否则删除，其中参数λ_D和λ_A是约束范围参数，其取值范围分别为λ_D∈[0.05,0.5]，λ_A∈[0.1,1)，本发明实例中当图像尺寸大于500×500时取λD=0.1和λA=0.1，当图像尺寸小于500×500时取λ_D＝0.2和λ_A＝0.2；最后将剩余的特征点对的集合作为初始候选匹配点对集{P₁′,P₂′}，该集合也是按行号一一对应的。

3.根据权利要求2所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于：步骤（5a）所描述的描述符相似度量矩阵和步骤（5b）所示的特征点匹配准则，具体如下：

（5a）计算浮动图像F的归一化描述符矩阵D_F与基准图像R的归一化描述符矩阵D_R之间的相似性度量矩阵AS＝D_F ^T×D_R，其中AS的矩阵元素AS_ij表示浮动图像F中的第i个特征点的归一化描述符向量与基准图像R中的第j个特征点的归一化描述符向量的内积，浮动图像F中第i个特征点的坐标为(x_1i,y_1i)，i=1,2,…,N_F，基准图像R中第j个特征点的坐标为j=(x_2j,y_2j)，j=1,2,…,N_R；

（5b）根据相似性度量矩阵和特征点匹配准则得到匹配特征点对集{P₁,P₂}，当相似性度量矩阵的元素AS_ij为第i行最小值且是第j列的最小值时，当前浮动图像F中的特征点i与基准图像R中的特征点j才是对应的，得到匹配点对{i,j}；所有N_p对相匹配的特征点的坐标构成匹配特征点对集{P₁,P₂}，其中P₁为浮动图像F中的匹配特征点坐标矩阵，该矩阵中的第i行元素为第i个匹配特征点的横坐标x_1i和纵坐标y_1i，P₂为基准图像R中的匹配特征点坐标矩阵，该矩阵中的第i行元素为与P₁中第i个匹配特征点相匹配的基准图像R的第i个匹配特征点的横坐标x_2j和纵坐标y_2j，P₁和P₂的维数均为N_p×2。

4.根据权利要求1所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于步骤（6）所描述的候选匹配点对集的生成过程，具体如下：

（6a）根据初始候选匹配点对集中的初始候选特征点矩阵分别构建浮动图像和基准图像的三角网，分别在浮动图像F中对初始候选匹配特征点矩阵P₁′、基准图像R中对初始候选匹配特征点矩阵P₂′按照Delaunay方法构建三角形网，得到浮动图F的t₁×3维的Delaunay三角形网矩阵TIN₁和基准图R的t₂×3维的Delaunay三角形网TIN₂，其中t₁和t₂分别表示浮动图像F和基准图像R中三角形的总数，TIN₁的每一行中的三个列元素分别代表浮动图像F中一个三角形的三个顶点在初始候选匹配点矩阵P₁′中的行序号，TIN₂的每一行中的三个列元素分别代表基准图像R中一个三角形的三个顶点在初始候选匹配点矩阵P₂′中的行序号；

（6b）寻找同名三角形对，将TIN₁的第一行与TIN₂的每一行逐行比较，当TIN₁和TIN₂中存在相同的行序号组合时，则找到了一个同名三角形对，保留TIN₁的当前行和TIN₂中相同行序号组合的行，否则删除TIN₁的当前行；

（6c）根据同名三角形对得到候选匹配点对集{Q₁,Q₂}，对TIN₁的每一行都执行步骤（6b）直到TIN₁的所有行都检测完毕，得到了TIN₁与TIN₂中所有三角形顶点行序号组合相同的三角形组成的同名三角形网络TIN_s；将浮动图像F的初始候选匹配点矩阵P₁′中的TIN_s中所有节点行序号构成浮动图像F中的同名匹配点矩阵Q₁，同样地，将基准图像R的初始候选匹配点矩阵P₂′中与P₁′中的行序号对应的同名匹配点矩阵Q₂，组成候选匹配点对集{Q₁,Q₂}。

5.根据权利要求4所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于步骤（6b）中所描述的同名三角形对判定准则，具体如下：

将TIN₁的第一行与TIN₂的每一行逐行比较，当TIN₁和TIN₂中存在相同的行序号组合时，则找到了一个同名三角形对，保留TIN₁的当前行和TIN₂中相同行序号组合的行，否则删除TIN₁的当前行。

6.根据权利要求1所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于步骤（7）所描述的最终匹配点对集生成过程，具体如下：

（7a）初始化仿射误差矩阵，构建一个大小为Nm×Nm×Nm的三维矩阵TRIerr，此三维矩阵中的每个矩阵元素初值均为1000，其中Nm为候选匹配点对集{Q₁,Q₂}的点对数目，

在候选匹配点对集{Q₁,Q₂}中任意选择3对候选匹配点对，第1对的行序号作为三维矩阵TRIerr的第一维行序号，第2对的行序号作为三维矩阵TRIerr的第二维列序号，第3对的行序号作为三维矩阵TRIerr的第三维页序号，则此矩阵元素的值赋值为仿射误差值，仿射误差值由下述步骤计算得到；

（7b）计算候选匹配点对集{Q₁,Q₂}构成三角形对的仿射误差，对于候选匹配点对集{Q₁,Q₂}中每次选取三对候选匹配点对的行序号组合为(n₁,n₂,n₃)，判断所对应的浮动图像F中的候选特征点集Q₁的点A=(x₁,y₁)、B=(x₂,y₂)和C=(x₃,y₃)是否满足构成三角形的条件(y₁-y₂)(x₁-x₃)≠(y₁-y₃)(x₁-x₂)，如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算对应的基准图像R中的候选特征点集Q₂的点、和是否满足构成三角形的条件如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算三对候选匹配点对所形成的三角形ΔABC和ΔA′B′C′的仿射误差值 $\mathrm{TRIerr}(n_1,n_2,n_3)=\mathrm{std}\left(\right[\frac{\left|\mathrm{AB}\right|}{|A\′B\′|},\frac{\left|\mathrm{BC}\right|}{|B\′C\′|},\frac{\left|\mathrm{AC}\right|}{|A\′C\′|}\left]\right),$ 其中，|AB|、|BC|、|AC|分别表示三角形ΔABC的三条边长，|A'B'|、|B'C'|、|A'C'|分别表示与三角形ΔABC组成同名三角形对的三角形ΔA′B′C′的对应顶点组成的三条边长，符号“[]”表示元素按行排列成为一个向量，函数std()表示标准差函数，n₁＝1,...,Nm-2，n₂＝n₁+1,...,Nm-1，n₃＝n₂+1,...,Nm；n₃＝n₂+1,...,Nm；

（7c）重复步骤（7b）直至遍历候选匹配点对集{Q₁,Q₂}中所有的候选匹配点对；

（7d）计算误差矩阵的最小值对应的候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}。计算误差矩阵TRIerr的最小值，将所有取到最小值的候选匹配点坐标逐行存入三角形顶点矩阵TRI中，将TRI中出现的对应候选匹配点对集{Q₁,Q₂}中的点对作为最终匹配点对集{Q₁′,Q′₂}。

7.根据权利要求6所述的基于点对约束和三角形网的有变化区域的合成孔径雷达图像配准方法，其特征在于步骤（7b）所描述的仿射误差计算过程及仿射三角形对仿射误差计算过程，具体如下：

对于候选匹配点对集{Q₁,Q₂}中每次选取三对候选匹配点对的行序号组合为(n₁,n₂,n₃)，判断所对应的浮动图像F中的候选特征点集Q₁的点A=(x₁,y₁)、B=(x₂,y₂)和C=(x₃,y₃)是否满足构成三角形的条件(y₁-y₂)(x₁-x₃)≠(y₁-y₃)(x₁-x₂)，如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算对应的基准图像R中的候选特征点集Q2的点 $\hat{B}=({\hat{x}}_2,{\hat{y}}_2)$ 和 $\hat{C}=({\hat{x}}_3,{\hat{y}}_3)$ 是否满足构成三角形的条件 $({\hat{y}}_1-\widehat{y_2})({\hat{x}}_1-{\hat{x}}_3)\≠({\hat{y}}_1-{\hat{y}}_3)({\hat{x}}_1-{\hat{x}}_2),$ 如果不满足条件则选取下一个的三对候选匹配点对重复判断构成三角形的条件，如果满足条件则计算三对候选匹配点对所形成的三角形ΔABC和ΔA′B′C′的仿射误差值 $\mathrm{TRIerr}(n_1,n_2,n_3)=\mathrm{std}\left(\right[\frac{\left|\mathrm{AB}\right|}{|A\′B\′|},\frac{\left|\mathrm{BC}\right|}{|B\′C\′|},\frac{\left|\mathrm{AC}\right|}{|A\′C\′|}\left]\right),$ 其中，|AB|、|BC|、|AC|分别表示三角形ΔABC的三条边长，|A'B'|、|B'C'|、|A'C'|分别表示与三角形ΔABC组成同名三角形对的三角形ΔA′B′C′的对应顶点组成的三条边长，符号“[]”表示元素按行排列成为一个向量，函数std()表示标准差函数，n₁＝1,...,Nm-2，n₂＝n₁+1,...,Nm-1，n₃＝n₂+1,...,Nm；n₃＝n₂+1,...,Nm。