本发明属于图像处理
技术领域：
，更进一步涉及遥感图像处理
技术领域：
中的一种基于结构信息和空间约束的遥感图像配准方法。本发明可用于不同时间、不同角度、不同传感器得到的灰度非线性变化的遥感图像配准。
背景技术：
：图像配准是将不同时间、不同角度或不同传感器拍摄的同一景物的两幅或多幅图像进行对准叠加的过程。给定参考图像和待配准图像，由于图像成像条件不同，图像在光照、对比度、尺度、旋转、灰度变化等方面可能都有很大差异，图像配准的主要目的就是要克服这些困难，找到一个从待配准图像到参考图像的最佳空间几何变换，将两幅图像变换到同一坐标系下，使得变换后的待配准图像相对于某一相似度量函数来说与参考图像具有最大程度的相似性。西安电子科技大学在其申请的专利“基于互信息选块和sift特征的遥感图像配准方法”(专利申请号：cn201410379927，公开号：cn104200461a)中提出了一种基于尺度不变性特征转换sift(scaleinvariantfeaturetransform)和互信息的遥感图像配准方法。该方法首先对图像进行分块，随机从参考遥感图像和待配准遥感图像中选取图像对，计算每一对图像的互信息；然后对互信息进行降序排列；选取前n个互信息较大的图像块提取尺度不变性特征转换sift特征，并进行粗匹配，去除错误匹配，计算两幅图像之间的变换矩阵参数和变换后的互信息值；最终选取互信息值最大的变换矩阵参数作为图像配准结果。该方法能够加快图像配准的速度，但是该方法在特征点检测、特征点主方向分配和特征点描述子生成阶段采用尺度不变性特征转换sift方法，在遥感图像对的灰度非线性变化较大的情况下正确匹配率会迅速下降，当遥感图像中存在较多相似特征时，该方法仅仅采用最近邻方法，得到较多错误匹配点对。floradellinger在其发表的论文“sar-sift:asift-likealgorithmforsarimages”(《ieeetransactionsongeoscienceandremotesensing》，2015，53(1)：453-466)中提出了一种构建基于合成孔径雷达的尺度不变性特征转换sar-sift(syntheticapertureradar-scaleinvariantfeaturetransform)特征描述子的方法。该方法利用指数加权均值比算子roewa(ratioofexponentiallyweightedaverages)定义了一种新的梯度计算方法比例梯度gr(gradientbyratio)，通过计算图像在各尺度上的基于合成孔径雷达的尺度不变性特征转换sar-sift矩阵得到基于合成孔径雷达的尺度不变性特征转换sar-sift尺度空间，然后在尺度空间上寻找极值点，最后提取特征点周围邻域局部特性，生成基于合成孔径雷达的尺度不变性特征转换sar-sift特征描述符。由于提取的基于合成孔径雷达的尺度不变性特征转换sar-sift特征具有尺度和旋转不变性，且对光照变化具有不变性，而且考虑到合成孔径雷达sar(syntheticapertureradar)图像中斑点噪声的统计特征，利用指数加权均值比算子而非差分法求图像梯度，对噪声较为鲁棒，因而被成功应用于合成孔径雷达sar图像配准领域。但是，该方法仍然存在的不足之处是，采用基于合成孔径雷达的尺度不变性特征转换sar-sift对遥感图像进行特征提取时，在提取图像特征过程中，当图像灰度变化不规律时，利用梯度信息生成特征点的特征描述符会丢失特征点的部分邻域信息，而且未考虑特征点的空间信息，具有相似特征向量的特征点容易导致正确匹配点对被遗漏，因此针对同一区域灰度变化很大的遥感图像配准问题，基于合成孔径雷达的尺度不变性特征转换sar-sift算法不能提供大量且准确的特征匹配点对，从而导致配准精度较低或失配。技术实现要素：本发明的目的在于针对上述现有技术的不足，提出一种基于结构信息和空间约束的遥感图像配准方法，可以利用遥感图像结构信息提升特征点特征描述符的有效性和稳定性，利用空间约束减少误匹配，增加正确匹配点对数，提高遥感图像配准的精度，解决了现有技术对遥感图像配准时效果不好的问题。实现本发明目的的思路是，将参考遥感图像和待配准遥感图像中所有像素点的梯度幅值和相位一致性值相加得到遥感图像像素点的结构信息值，采用结构信息生成特征点的特征描述符，通过最近邻方法筛选出初次匹配点对，利用空间约束筛选匹配点对，增加正确匹配点对数量。本发明的具体步骤包括如下：(1)读取数据：输入两幅遥感图像，将其中一幅遥感图像作为参考遥感图像，将另一幅遥感图像作为待配准遥感图像；(2)检测特征点：(2a)利用sar-harris尺度空间公式，分别生成参考遥感图像和待配准遥感图像sar-harris尺度空间中的幅值图像；(2b)利用阈值比较法，分别检测参考遥感图像和待配准遥感图像的sar-harris尺度空间幅值图像上的特征点；(3)采用结构信息生成特征描述符：(3a)利用相位一致性公式，分别计算参考遥感图像和待配准遥感图像的所有像素点的相位一致性值；(3b)利用梯度幅值公式，分别计算参考遥感图像和待配准遥感图像所有像素点的梯度幅值；(3c)利用梯度幅值公式，分别计算参考遥感图像和待配准遥感图像所有像素点的梯度幅值；(3d)通过将具有相同坐标的像素点的相位一致性值和梯度幅值相加，分别计算出参考遥感图像和待配准遥感图像所有像素点得结构信息值；(3e)利用直方图统计法，分别生成参考遥感图像和待配准遥感图像中初始特征点的描述符；(4)利用空间约束筛选匹配点对；(4a)利用最近邻方法，从参考遥感图像和待配准遥感图像的特征点中筛选出初次匹配点对；(4b)利用空间约束法，从参考遥感图像和待配准遥感图像的特征点中二次筛选匹配点对；(5)获得几何形变参数：采用随机抽样一致性算法，计算参考遥感图像和待配准遥感图像之间的变换矩阵，得到待配准遥感图像的几何形变参数；(6)将待配准遥感图像进行几何变换：利用得到的几何形变参数，将待配准遥感图像进行几何变换，得到配准结果。本发明与现有技术相比具有以下优点：第一，由于本发明采用结构信息生成特征点的特征描述符，充分反映特征点的邻域特征信息，克服了现有技术在在图像灰度变化不规律时，利用梯度信息生成的特征点的特征描述符会丢失特征点的部分邻域特征信息的问题，使得本发明提高了获取特征点邻域特征信息的能力。第二，由于本发明利用空间约束筛选匹配点对，增加正确匹配点对，克服了现有技术中采用最近邻方法，得到较多错误匹配点对的问题，使得本发明提高了滤除错误匹配点对的能力。附图说明图1是本发明的流程图；图2是本发明仿真实验1的仿真图；图3是本发明仿真实验2的仿真图。具体实施方式下面结合附图对本发明做进一步描述。参照图1，本发明的实施步骤做进一步描述。步骤1，读取数据。输入两幅遥感图像，将其中一幅遥感图像作为参考遥感图像，将另一幅遥感图像作为待配准遥感图像。步骤2，检测特征点。利用sar-harris尺度空间公式，分别生成参考遥感图像和待配准遥感图像sar-harris尺度空间中的幅值图像。所述的sar-harris尺度空间公式如下：r＝det(c)-0.04·tr(c)其中，r表示一幅遥感图像的sar-harris尺度空间，det表示矩阵的行列式操作，c表示一幅遥感图像的sar-harris矩阵，tr表示对矩阵主对角线元素的求和操作。利用阈值比较法，分别检测参考遥感图像和待配准遥感图像的sar-harris尺度空间幅值图像上的特征点。所述的阈值比较法的步骤如下：第1步，在sar-harris尺度空间幅值图像中选取一个像素点，以所选取的像素点为中心，取3x3邻域，得到8个邻域点。第2步，重复第一步，依次从sar-harris尺度空间幅值图像选择像素点，直到得到所有像素点的8个邻域点。第3步，比较sar-harris尺度空间幅值图像中每个像素点与其相邻8个像素点的灰度值。第4步，将中心点的灰度值大于其周围8个邻域点的灰度值的点作为极值点。第5步，将所得到的极值点中灰度值大于0.8的像素点作为特征点。步骤3，采用结构信息生成特征描述符。利用相位一致性公式，分别计算参考遥感图像和待配准遥感图像所有像素点的相位一致性值。相位一致性公式如下：ω(x,y)＝cos(φn(x,y)-ψ(x,y))-|sin(φn(x,y)-ψ(x,y))|其中，pc表示遥感图像像素点的相位一致性值，x和y分别表示遥感图像中像素点的行坐标和列坐标，∑表示求和操作，n表示加伯gabor滤波器的尺度标号，取值范围为1～8，w表示频率扩展的权重系数，取值范围为0～1，表示向上取整操作，an表示加伯gabor滤波器尺度标号为n时的幅度响应值，ω表示加伯gabor滤波器的相位偏移量，t表示加伯gabor滤波器的噪声阈值，取值为0.1，cos表示余弦操作，φn表示加伯gabor滤波器尺度标号为n的相位响应值，|·|表示绝对值操作，sin表示正弦操作，ψ表示加伯gabor滤波器在各个尺度标号的平均相位值。利用梯度幅值公式，分别计算参考遥感图像和待配准遥感图像所有像素点的梯度幅值。所述的梯度幅值公式如下：其中，g表示遥感图像像素点的梯度幅值，x和y分别表示遥感图像中像素点的行坐标和列坐标，i表示遥感图像，表示开平方操作，表示偏导操作。利用梯度方向公式，分别计算参考遥感图像和待配准遥感图像所有像素点的梯度幅值。所述的梯度方向公式如下：其中，θ表示遥感图像像素点的梯度方向角度值，x和y分别表示遥感图像中像素点的行坐标和列坐标，tan-1表示反正切操作，表示偏导操作，i表示遥感图像，/表示相除操作。通过将具有相同坐标的像素点的相位一致性值和梯度幅值相加，分别得到参考遥感图像和待配准遥感图像所有像素点的结构信息值。利用直方图统计法，分别生成参考遥感图像和待配准遥感图像中所有特征点的描述符。直方图统计法的步骤如下：第1步，在参考遥感图像中，选取一个特征点。第2步，将所选取的特征点作为中心点，在其周围取一个16x16的邻域。第3步，用所选特征点16×16邻域内所有像素点的梯度方向角度值，构造梯度方向直方图。第4步，选取梯度方向直方图的最大值作为所选特征点的主方向角度值。第5步，用主方向角度值更新所选特征点梯度方向的角度。；第6步，将0°～360°梯度方向角度区间，划分为0°～44°、45°～89°、90°～134°、135°～179°、180°～224°、225°～269°、270°～314°、315°～359°八个角度子区间。第7步，将所选特征点的16×16邻域窗口平均划分成16个4×4的子窗口。第8步，依据每个子窗口中的每个像素点的梯度方向角度值大小，将其划分到八个角度子区间内。第9步，将每个子窗口内每个子区间中所有像素点的结构信息值相加，得到八个窗口结构信息值，将其存放在一个128维的列矢量中。第10步，将这个128维的列矢量作为所选特征点的特征描述符。第11步，重复执行第1步至第10步，依次选取特征点，直到得到参考遥感图像中所有特征点的特征描述符。第12步，在待配准图像中选取一个特征点。第13步，执行第1步至第12步操作一次，得到待配准遥感图像中所选特征点的特征描述符。第14步，重复12步和13步，依次选取特征点，直到得到待配准遥感图像中所有特征点的特征描述符。步骤4，利用空间约束筛选匹配点对。利用最近邻方法，从参考遥感图像和待配准遥感图像的特征点中筛选出初次匹配点对。所述的最近邻方法的步骤如下：第1步，分别将参考遥感图像和待配准遥感图像中的特征点存放在参考图像特征点集合u和待配准遥感图像特征点集合v。第2步，在参考遥感图像特征点集合u中选取一个特征点p。第3步，计算参考遥感图像所选取的特征点p的特征描述符和待配准遥感图像特征点集合v中所有特征点的特征描述符的欧氏距离。第4步，将参考遥感图像中所选取的特征点p与待配准遥感图像中所有特征点的特征描述符的欧式距离进行排序，得到最小值d1和次最小值d2，以及取得最小值时待配准遥感图像的特征点o，若d1/d2＜0.6，将p和o作为初次匹配点对。第5步，重复执行第2步至第4步，直到参考遥感图像特征点集合u中所有特征点都被选取。第6步，将所有匹配点对合并成一个匹配点对集合k。利用空间约束法，从参考遥感图像和待配准遥感图像的特征点中二次筛选匹配点对。所述的空间约束法的步骤如下：第1步，依次计算匹配点对集合k中所有匹配点对两个像素点的主方向角度值之差，得到匹配点对集合k的主方向角度差均值ω。第2步，在匹配点对集合k中任意选取两对匹配点，计算两对特征点中位于参考遥感图像中的两个特征点之间距离与两对特征点中位于浮动遥感图像中两个特征点之间距离的比值r。第3步，从参考遥感图像和浮动遥感图像中各自任取一个特征点，组成一个特征点对。第4步，重复第3步，直到参考遥感图像任意一个特征点都和浮动遥感图像中所有特征点都组成特征点对，所有特征点对合并成特征点对集合l。第5步，从所得到的特征点对集合l中选取一个特征点对，计算所选特征点对两个像素点的主方向角度值之差，判断是否满足点约束公式，若满足该公式，执行第六步，否则执行第七步，点约束公式如下。|ω-e|＜τ其中，|·|表示绝对值操作，e表示所选取的特征点对的主方向差，τ表示点约束阈值，取值为3。第6步，从匹配点对集合k任意选取一个匹配点对，判断所选取的特征点对两个像素点的坐标是否满足线约束公式，若满足该公式，将所选取特征点对作为最终匹配点对，否则执行第七步，线约束公式如下：其中，|·|表示求绝对值操作，r表示匹配点对集合k中任意选取两对匹配点距离比值，表示开平方操作，(a,b)和(c,d)分别表示所选匹配点对两个像素点的坐标，(x,y)和(u,v)分别是两幅遥感图像中所选取的特征点对两个像素点的坐标，ρ表示线约束阈值，取值为0.1。第7步，重复第五步和第六步，依次从特征点对集合l中选取特征点对，直到所有特征点对都被选取。第8步，将所有最终匹配点对合并成最终匹配点对集合。步骤5，获得几何形变参数。采用随机抽样一致性算法，计算参考遥感图像和待配准遥感图像之间的变换矩阵，得到待配准遥感图像的几何形变参数。步骤6，将待配准遥感图像进行几何变换。利用得到的几何形变参数，将待配准遥感图像进行几何变换，得到配准结果。下面结合实验仿真对本发明效果做进一步说明。1.仿真条件：本发明的仿真实验平台采用intel(r)core(tm)cpui5-65003.20ghz，内存为8gb，运行windows7的pc机，编程语言为matlab2016b。2.仿真内容与结果分析：图2为本发明仿真实验1的仿真图。其中图2(a)和图2(b)分别是usgs项目中的两幅多光谱遥感图像，图像大小均为614x611个像素。图2(a)为本发明用于遥感图像配准的参考图，由landsat-7etm+传感器于2000年7月拍摄，图2(b)为本发明用于遥感图像配准的待配准图，由landsat4-5tm传感器于1999年6月拍摄，可以看出图2(a)和图2(b)对应区域存在着非线性灰度变化。图2(c)为利用本发明的方法对图2(a)和2(b)进行配准得到的配准结果融合图。从图2(c)可以看出，本发明方法得到的待配准图与参考图中对应区域的图像纹理信息准确重合在一起，说明本发明对遥感图像进行了精确配准。本发明的仿真实验中利用空间约束法，筛选得到107对匹配点对。为了评价遥感图像配准的精确度，通过下述的均方根误差公式，对107对匹配点对进行加权平均，计算得到均方根误差值rmse(rootmeansquareerror)。所述的均方根误差计算公式如下：其中，rmse表示均方根误差值，表示开平方操作，n表示筛选得到的匹配点对个数，∑表示求和操作，i表示匹配点对的标号，xi和yi表示标号为i的两个特征点其中一个特征点的行坐标值和列坐标值，ui和vi表示标号为i的两个特征点中另一个特征点的行坐标值和列坐标值。通过上述均方根误差计算公式，计算得到107对匹配点对之间的均方根误差值为0.5997。现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法得到0对匹配点对，匹配点对之间的均方根误差值为0。将本发明与现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法分别得到的匹配点对数和均方根误差rmse值填入表1中，得到本发明与现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法的对比实验结果一览表。从表1可以看出，现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法不能得到匹配点对，表明该方法对图2(a)和图2(b)配准失效，而本发明能够成功配准该组图像，这些数据说明本发明方法能够成功对具有不规律灰度变化的遥感图像进行配准。表1本发明与sar-sift方法实验结果比较一览表sar-sift方法本发明匹配点对数0107均方根误差值00.5997图3为本发明仿真实验2的仿真图。其中图3(a)和3(b)分别为同一传感器alos-palsar分别于2006年6月和2009年3月拍摄的武汉地区的两幅合成孔径雷达sar图像，图像大小均为400x400个像素。图3(a)为本发明用于遥感图像配准的参考图，图3(b)为本发明用于遥感图像配准的待配准图，图3(c)为利用本发明方法对图3(a)和3(b)所示的图像进行配准，得到的配准结果融合图。从图3(c)可以看出，本发明方法得到的待配准图与参考遥感图中对应区域的图像纹理信息准确重合在一起，说明本发明对遥感图像进行了精确配准。本发明的仿真实验中利用空间约束法，筛选得到67对匹配点对。通过下述的均方根误差公式，对67对匹配点对进行加权平均，计算得到均方根误差值rmse(rootmeansquareerror)。所述的均方根误差计算公式如下：其中，rmse表示均方根误差值，表示开平方操作，n表示筛选得到的匹配点对个数，∑表示求和操作，i表示匹配点对的标号，xi和yi表示标号为i的两个特征点其中一个特征点的行坐标值和列坐标值，ui和vi表示标号为i的两个特征点中另一个特征点的行坐标值和列坐标值。通过上述均方根误差计算公式，计算得到67对匹配点对之间的均方根误差值为0.5090。现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法得到8对匹配点对，匹配点对之间的均方根误差值为0.5239。将本发明与现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法分别得到的匹配点对数和均方根误差rmse值填入表2，得到本发明与现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法的对比实验结果一览表。从表2可以看出，相比于现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法，本发明方法对图3(a)和图3(b)配准得到的匹配点对更多，均方根误差值更小，说明本发明增加了匹配点对，提高了遥感图像配准的精度。表2本发明与传统sar-sift方法实验结果比较一览表sar-sift方法本发明匹配点对数867均方根误差值0.52390.5090综上所述，本发明分别对两组具有非线性灰度变化的遥感图像配准进行了仿真实验，相比于现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法，本发明方法均能得到较多的匹配点对和较小的均方根误差值，配准精度较高。两个仿真实验均表明，本发明采用结构信息生成特征点的特征描述符，利用空间约束筛选正确匹配点对，有效改善了现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法中由于非线性灰度变化导致筛选出错误匹配点对，增加了正确匹配点对，解决了现有技术基于合成孔径雷达的尺度不变性特征转换sar-sift方法应用于遥感图像配准失效或精度较低的问题。当前第1页12