基于各向异性梯度尺度空间的遥感图像配准方法技术

本发明专利技术公开了一种基于各向异性梯度尺度空间的遥感图像配准方法，主要解决遥感图像之间亮度非线性变化较大情况时正确匹配率较低的问题。其实现步骤为：(1)输入遥感图像对；(2)构造各向异性扩散的尺度空间；(3)计算梯度幅度图像；(4)检测特征点；(5)生成特征点主方向；(6)生成特征点描述子；(7)特征点匹配；(8)删除错误匹配的特征点对；(9)配准参考图像和待配准图像。本发明专利技术在各向异性尺度空间的梯度幅度图像上进行特征点检测、特征点主方向生成和特征点描述子生成，能有效应对图像之间亮度有较大分线性变化的情况，可用于复杂多源和多光谱遥感图像配准。

【技术实现步骤摘要】
基于各向异性梯度尺度空间的遥感图像配准方法
本专利技术属于图像处理
，更进一步涉及对遥感图像进行配准处理
中的一种基于各向异性梯度尺度空间的遥感图像配准方法。本专利技术可应用于对多光谱、多源不同时刻获取的遥感图像的配准。
技术介绍
图像配准是指在不同时段，对同一场景从不同视角使用相同或者不同传感器拍摄的有重叠区域图像进行几何矫正的过程。图像配准技术是近年来发展迅速的图像处理技术之一。图像配准技术已经广泛用于各个领域，如航空航天技术、图像镶嵌、地理信息系统、图像融合、三维重建、目标识别和变化检测等等。目前随着遥感技术的发展，由于传感器的不同物理特性所产生的遥感图像不断增多，因此综合利用各种图像进行数据提取和分析已成为遥感领域的一个重要手段。同时由于各种传感器之间物理特性和成像方式的不同，在数据应用和数据融合时不同几何特性和不同分辨率的图像间必须进行严格的配准。目前遥感图像配准主要分为两类：基于区域灰度的配准方法和基于特征的配准方法。其中，基于区域灰度的常用图像配准方法有：互相关(CC)、不变矩、基于FFT的相位相关法和互信息(MI)等。但是基于灰度的图像配准方法存在着以下缺点：①对图像灰度变化比较敏感，尤其是非线性的光照变化，这大大降低了算法的性能；②计算复杂度高；③对目标的旋转、形变和遮挡比较敏感。为了克服其存在的缺点，人们提出了基于特征的图像配准方法。基于图像特征的方法首先从图像中提取边缘、角点、轮廓和区域中心等特征，然后使用特征之间的相关决定图像的最优对齐。这些显著特征可以大大压缩图像的信息量，使得计算量减小，速度更快，而且对图像的灰度变化具有鲁棒性。目前来说点特征是最为常用且效率较高的一种方法。二维图像中的角点、拐点、交叉点等是图像的明显特征。这些点具有平移、旋转和缩放不变性，几乎不受光照条件的影响。它只需用图像中大约0.05％的像素点就可以表示整幅图像的数据信息，其信息含量高、计算速度快，使实时处理成为可能。因此，基于特征的图像配准方法在遥感图像配准领域也得到了广泛的应用。经典的SIFT算法具有对可见光图像旋转、缩放、部分的光照和视角变换保持不变性的优点。当SIFT方法检测特征点数目增加时，生成特征点主方向和生成特征点描述子步骤花费时间迅速增加，限制了SIFT算法在实际中的应用。由于SIFT算法的优点，国内外学者提出了许多基于SIFT的提升算法。然而当使用基于SIFT的方法匹配遥感图像时出现了许多不正确的匹配点，正确匹配率(CMT)迅速下降。原因在于拍摄时间、光谱和获取传感器的不同导致遥感图像对同一区域的像素亮度显著不同，同时像素对之间的亮度映射可能是线性、非线性，甚至会出现亮度反转情况。Sedaghat在其发表的论文“UniformRobustScale-InvariantFeatureMatchingforOpticalRemoteSensingImages”(《IEEETransactionsonGeoscienceandRemoteSensing》,2011,pages4516-4527)中提出了一种UR-SIFT算法。该算法较好的控制了尺度空间提取的特征点的数量、质量和分布，提高了具有局部变换的部分遥感图像的配准精度，但是该算法仍然存在的不足之处是，不能精确的配准图像亮度有较大非线性变化的遥感图像对。Gong在其发表的论文“ANovelCoarse-to-FineSchemeforAutomaticImageRegistrationBasedonSIFTandMutualInformation”(《IEEETransactionsonGeoscienceandRemoteSensing》,2014,pages4328-4338)中公开了一种由粗到精的配准方法。该方法先使用SIFT算法获得图像对之间的初始变换关系，然后结合初始变换关系和互信息方法获得更加精确的图像配准。该方法仍然存在的不足之处是，由于该方法是基于互信息的方法，所以计算复杂度高，而且配准精度依懒于sift算法获得的初始解，当图像亮度存在较大非线性变化的情况下，sift算法并不能获得很好的初始解，依赖sift算法获得的初始解的互信息方法也不能实现精确的配准。西安电子科技大学在其申请的专利“基于互信息选块和sift特征的遥感图像配准方法”(专利申请号：CN201410379927，公开号：CN104200461A)中提出了一种基于sift和互信息的遥感图像配准方法。该方法的实现过程为：随机从参考遥感图像和待配准遥感图像中选取图像对，计算每一对图像的互信息；对互信息进行降序排列；选取前n个互信息较大的子图像对；对选取的前n个互信息较大的图像对提取sift特征，并进行粗匹配；删除错误匹配点，进行细匹配；计算配准参数和互信息值；选取互信息最大的配准参数作为最终的配准参数。该方法虽然能够加快图像配准的速度，但是，该方法仍然存在的不足之处是，在特征点生成、特征点主方向生成和特征点的描述符生成阶段仍采用sift方法，在遥感图像对的亮度非线性变化较大的情况下正确匹配率迅速下降，同时该方法采用传统的最近邻和次近邻距离比匹配准则，当遥感图像中存在较多重复特征情况下正确配对的特征点数目也迅速下降。
技术实现思路
本专利技术的目的是克服上述已有技术的不足，提出一种基于各向异性梯度尺度空间的遥感图像配准方法，以提高特征匹配的准确度，实现对图像亮度有较大非线性变化的遥感图像对的配准。本专利技术实现上述目的的思路是：首先根据非线性扩散原理构建参考遥感图像和待配准遥感图像的各向异性尺度空间，再对参考遥感图像和待配准遥感图像的各向异性尺度空间的梯度幅度图像上使用Harris算子进行角点检测，接着生成点特征的对图像亮度的非线性变化有较强鲁棒性的描述子，最后使用改进的特征点匹配准则进行特征点匹配，得到最后的遥感图像对配准结果图。本专利技术的步骤包括如下：(1)输入参考遥感图像和待配准遥感图像；(2)构造各向异性扩散的尺度空间：(2a)计算参考遥感图像和待配准遥感图像的各向异性尺度空间各层的尺度值；(2b)将尺度值转换到时间度量值；(2c)对输入的参考遥感图像和待配准遥感图像，采用标准差为σ0的高斯滤波，得到参考遥感图像和待配准遥感图像各向异性尺度空间第0层图像；(2d)将各向异性尺度空间层的序号i初始化为零；(2e)按照下式，分别计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵：其中，ci表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii应用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％百分位的梯度幅度值；(2f)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做行扩散：其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii的行扩散后的图像，I1表示行数和列数均与参考遥感图像或者待配准遥感图像的列数相同的单位方阵，ti和ti+1分别表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值本文档来自技高网...

【技术保护点】
一种基于各向异性梯度尺度空间的遥感图像配准方法，包括如下步骤：(1)输入参考遥感图像和待配准遥感图像；(2)构造各向异性扩散的尺度空间：(2a)计算参考遥感图像和待配准遥感图像的各向异性尺度空间各层的尺度值；(2b)将尺度空间值转换到时间度量值；(2c)对输入的参考遥感图像和待配准遥感图像，采用标准差为σ0的高斯滤波，得到参考遥感图像和待配准遥感图像各向异性尺度空间第0层图像；(2d)将各向异性尺度空间层的序号i初始化为零；(2e)按照下式，分别计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵：ci=11+(|▿Igi|K)2]]>其中，ci表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii应用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％百分位的梯度幅度值；(2f)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做行扩散：I1i+1=(I-2(ti+1-ti)·A1(Ii))-1Ii]]>其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii的行扩散后的图像，I表示行数和列数均与参考遥感图像或者待配准遥感图像的列数相同的单位方阵，ti和ti+1分别表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值，A1(Ii)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像扩散系数ci的矩阵，Ii表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像，(·)‑1表示逆矩阵操作；(2g)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做列扩散：I2i+1=(I-2(ti+1-ti)·A2(Ii))-1Ii]]>其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii的列扩散后的图像，I表示行数和列数均与参考遥感图像或者待配准遥感图像的行数相同的单位方阵，ti和ti+1分别表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值，A2(Ii)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像扩散系数ci的矩阵，Ii表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像，(·)‑1表示逆矩阵操作；(2h)按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i+1层图像：Ii+1=12(I1i+1+I2i+1)]]>其中，Ii+1表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i+1层图像，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像经过行扩散后的图像，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像经过列扩散后的图像；(2i)判断i≥N‑2是否成立，若成立，得到参考遥感图像和待配准遥感图像的各向异性尺度空间，否则，令i＝i+1，执行步骤(2e)；其中，i表示各向异性尺度空间层的序号，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；(3)计算梯度幅度图像：使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像；(4)计算梯度幅度图像的差分：(4a)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的x轴正方向的差分：GIxn=-101-202-101⊕GIn]]>其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIn沿水平向右的x轴正方向的差分，GIn表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作，n＝0,1,...,N‑1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；(4b)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的y轴正方向的差分：GIyn=-1-21000121⊕GIn]]>其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIn沿竖直向下的y轴正方向的差分，GIn表示输入参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作，n＝0,1,...,N‑1，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；(5)按照下式，计算梯度幅度图像的梯度幅度：GGIn=(GIxn)2+(GIyn)2]]>其中，GGIn表示参考遥感图像或...

【技术特征摘要】
1.一种基于各向异性梯度尺度空间的遥感图像配准方法，包括如下步骤：(1)输入参考遥感图像和待配准遥感图像；(2)构造各向异性扩散的尺度空间：(2a)计算参考遥感图像和待配准遥感图像的各向异性尺度空间各层的尺度值；(2b)将尺度值转换到时间度量值；(2c)对输入的参考遥感图像和待配准遥感图像，采用标准差为σ0的高斯滤波，得到参考遥感图像和待配准遥感图像各向异性尺度空间第0层图像；(2d)将各向异性尺度空间层的序号i初始化为零；(2e)按照下式，分别计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵：其中，ci表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像的扩散系数矩阵，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii应用标准差为1的高斯滤波后的图像，表示高斯滤波后的图像的梯度幅度，|·|表示取模操作，K表示对比度因子，K的取值为梯度幅度的统计直方图70％百分位的梯度幅度值；(2f)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做行扩散：其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii的行扩散后的图像，I1表示行数和列数均与参考遥感图像或者待配准遥感图像的列数相同的单位方阵，ti和ti+1分别表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层和第i+1层的时间度量值，A(Ii)表示编码参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像扩散系数ci的矩阵，Ii表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像，(·)-1表示逆矩阵操作；(2g)按照下式，对参考遥感图像和待配准遥感图像各向异性尺度空间的第i层图像做列扩散：其中，表示沿参考遥感图像或者待配准遥感图像各向异性尺度空间的第i层图像Ii的列扩散后的图像，I2表示行数和列数均与参考遥感图像或者待配准遥感图像的行数相同的单位方阵；(2h)按照下式，计算参考遥感图像和待配准遥感图像各向异性尺度空间的第i+1层图像：其中，Ii+1表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第i+1层图像；(2i)判断i≥N-2是否成立，若成立，得到参考遥感图像和待配准遥感图像的各向异性尺度空间，否则，令i＝i+1，执行步骤(2e)；其中，N表示参考遥感图像和待配准遥感图像各向异性尺度空间层的总数；(3)计算梯度幅度图像：使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像；(4)计算梯度幅度图像的差分：(4a)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的x轴正方向的差分：其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIn沿水平向右的x轴正方向的差分，GIn表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像，表示相关操作；(4b)按照下式，使用索贝尔算子Sobel，计算参考遥感图像和待配准遥感图像各向异性尺度空间的梯度幅度图像的y轴正方向的差分：其中，表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像GIn沿竖直向下的y轴正方向的差分；(5)按照下式，计算梯度幅度图像的梯度幅度：其中，GGIn表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIn的梯度幅度，表示开根号操作；(6)按照下式，计算梯度幅度图像的梯度角度：其中，AGGIn表示参考遥感图像或者待配准遥感图像各向异性尺度空间的梯度幅度图像的第n层图像GIn的梯度角度，arctan(·)表示四象限反正切操作；(7)检测特征点：使用Harris角点检测算子，在参考遥感图像和待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测特征点，得到特征点集合：其中，CIR表示在参考遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点集合，CIS表示在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点集合，R表示在参考遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点的总数，S表示在待配准遥感图像的各向异性尺度空间的梯度幅度图像上检测的特征点的总数；(8)生成特征点主方向：(8a)将水平方向角度在[0,2π]内36等分；(8b)由特征点集合的半径为6σn的圆形邻域区域内像素的梯度方向AGGIn(X)，确定圆形区域内像素所在的等分角度范围，其中，σn表示参考遥感图像或者待配准遥感图像各向异性尺度空间的第n层的尺度值，AGGIn(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度角度；(8c)分别累加位于每个等分角度范围内的所有像素的梯度幅度值GGIn(X)，形成梯度方向直方图，其中，GGIn(X)表示特征点周围圆形邻域内位置坐标为X的像素的梯度幅度；(8d)将梯度方向直方图中大于最大值0.8倍的数值对应的梯度方向，作为特征点的主方向；(9)生成特征点描述子：(9a)将半径为ρ的特征点的圆形邻域逆时针旋转θm度，圆形邻域沿半径方向将ρ分为3个区间，内圆邻域半径为0.25ρ，中间圆邻域半径为0.73ρ，外圆邻域半径为ρ，圆形邻域沿角度方向将[0,2π]等分为8个区间，内圆作为一个整体，特征点周围圆形邻域被划分为了17个不同位置的面积相等的子区域，其中，ρ的取值为12σn，θm表示该特征点的主方向；(9b)将特征点圆形邻域内像素的笛卡尔坐标转换为对数极坐标，对数极坐标角度方向水平向右，在[0,2π]的范围内等分为8个区间，对数极坐标对数长度方向竖直向下，在范围内非等分为3个区间，其中，ρ表示特征点周围圆形邻域的半径；(9c)对数极坐标网格中每个子区域内的所有像素根据其梯度幅度GGIn(X)和梯度方向AGGIn(X)计算梯度方向直方图，每一个子区域形成了一个8维的梯度方向向量，依次拼接17子区域的梯度方向向量就形成了一个136维的特征点的描述子，其中，梯度方向直方图的角度在[0,2π]的范围内等分为8个区间；(9d)采用步骤(9a)、步骤(9b)、步骤(9c)的相同方法，生成特征点集合CIR的描述子集合DIR和特征点集合CI...

【专利技术属性】
技术研发人员：马文萍，闻泽联，焦李成，武越，马进，任琛，
申请(专利权)人：西安电子科技大学，
类型：发明
国别省市：陕西;61