专利名称：基于边缘和光谱反射率曲线的多时相遥感图像配准方法
技术领域：
本发明属于光学遥感图像处理技术领域，涉及同一多光谱传感器在不同时间获取的相同地理位置的遥感图像间存在较大平移和旋转错位情况下的配准，可用于多时相多光谱遥感图像变化检测、融合、镶嵌等方面的前期配准工作。
背景技术：
近年来，随着多光谱遥感图像数据的公开化程度和处理技术水平的不断提高，多光谱遥感图像数据已经被广泛用于城市规划、资源调查、自然灾害预防、农作物长势监测、军事目标探查等多种应用领域。多光谱遥感图像的广泛应用首先需要进行图像数据的几何校正、辐射校正、配准、去噪、分割、镶嵌、融合、目标检测、变化检测等多项数字图像处理技术。其中，在进行多光谱遥感图像变化检测、融合以及镶嵌等工作前，首先需要对不同传感器获取的同一地区的图像或相同传感器在不同时间获取的同一地区的图像进行配准，以消除因获取图像的时间、角度、环境和传感器成像机理的不同造成图像间的平移、旋转、伸缩及局部形变问题。1992 年 L.G.Brown 在文献”A Survey of Image Registration Techniques'ACMComputing Surveys, 1992, vol.24, n0.4, pp.325-376)中指出，所有的图像配准技术都是由特征空间、搜索空间、搜索策略和相似性测度四个方面构成。特征空间指图像中用于进行搜索匹配的信息集合，是搜索匹配的对象；而每一次搜索的匹配程度是以相似性测度为评价标准的，因而相似性测度直接关系到配准的精确度；搜索空间是搜索匹配得到的所有图像变换参数所组成的空间；搜索策略则是搜索最优匹配参数时所采取的优化策略。多光谱遥感图像配准的过程一般为首先确定特征空间和搜索空间，然后以特征空间为基础设计并采用某种相似性测度，进而按一定的搜索策略搜索最佳匹配变换参数，并将其带入图像变换模型从而达到配准的目的。根据特征空间的不同，现有的多光谱遥感图像配准方法通常可以分为基于像素灰度特征和基于图像几何特征两大类。大部分基于像素灰度特征的方法一般直接利用整幅图像的灰度信息，通过建立某种像素间的相似性测度来衡量两幅图像重叠部分地表反射属性的匹配程度，进而寻找到最优匹配时的平移、旋转和伸缩等变换参数。常见的灰度相似性测度如平均平方差MSD、平均绝对差MAD、归一化互相关NCC等通常要求待配准图像对应像素的灰度值具有线性关系，因而在利用它们进行不同时相的图像配准时容易受到噪声和图像间整体亮度差异的影响，而对于不同传感器或不同波段图像的配准则基本不适用。最初作为多模态医学图像配准相似性测度的最大互信息准则由于不要求待配准图像对应像素的灰度信息具有严格的线性关系，而认为只要两幅图像中主要的几何结构特征对准时即达到最佳匹配，因此已被用于不同时相和不同波段图像间的配准。在对不同时相的多光谱遥感图像的配准中，由于基于像素灰度特征的配准方法在搜索最佳匹配变换参数时通常是由全体像素共同决定图像间的对准程度，计算量大，并且难免会受到图像中不同位置大小的云、阴影及变化地物的干扰，造成配准精度和稳定性的下降。
基于图像几何特征的配准方法则首先需要对待配准图像提取角点、边缘等明显的几何特征，然后进行对应特征对或特征集之间的搜索匹配，进而得到图像变换参数。基于边缘特征的配准方法常见的如利用一些传统空域边缘检测算子进行边缘特征提取，但往往由于提取的边缘过于琐碎，不能很好地体现图像中的结构特征，对后续匹配精度和速度造成影响，因而需要对提取到的边缘进行优化处理。近年来，小波变换wavelet、曲波变换curvelet和轮廓波变换contourlet等多尺度几何分析方法由于在描述和分析图像中的边缘特征时具有优异的性能，也常被用于边缘特征的提取。1997年J.W.Hsieh.H.Y.M.Liao、K.C.Fan 等在文献，，Image registration using a new edge-based approach，，(ComputerVision and Image Understanding, 1997，vol.67，n0.2，pp.112-130)中首次提出利用小波变换得到系数的模局部极大值进行边缘特征点提取，然后以互相关作为相似性测度进行边缘特征点对的搜索匹配，其中还使用了一种非迭代的策略消除误匹配点对，达到了较好的配准效果。2008年孙淑一.吴勇和吴建民在文献”一种基于边缘特征的图像配准方法”(计算机工程与应用，2008，vol.44，n0.7，pp.94-96)中、2010年孙雅琳和王菲在文献”基于边缘和互信息的红外与可见光图像配准”(电子科技，2010, vol.23，n0.4，pp.80-82)中也都采用了相同的边缘特征提取方法，不同的是，这两种方法在搜索匹配时前者计算了两幅边缘图像重叠部分的交互方差，而后者计算了两幅边缘图像重叠部分的归一化互信息。如果将这两种方法用于不同时相多光谱遥感图像的配准，则容易受到图像中变化部分对应特征点的干扰而造成误配准。2009年陈志刚、尹福昌和孙孚在文献”基于非采样Contourlet变换高分辨率遥感图像配准”(光学学报，2009, vol.29，n0.10，pp.2744-2750)中提出首先由图像通过非下采样Contourlet变换得到的各方向子带计算得到梯度矢量模，阈值分割得到边缘特征点后由归一化互相关匹配搜索匹配控制点，并利用概率支持算法消除误匹配点对，降低了变化部分对应特征点错误匹配的概率。总体来说，由于不直接采用像素的灰度值作为相似性测量的对象，因而该类方法受不同时相间的亮度差异、不同传感器以及不同波段图像间的成像特性差异的影响较小，适合多时相多光谱遥感图像的配准，但要求待配准图像中含有比较明显的边缘信息，否则提取不到理想的边缘将很可能导致配准失败。另外，其中边缘特征提取和匹配方法的准确性及稳定性也是影响其配准性能的主要因素
发明内容
本发明的目的在于针对上述已有技术的不足，提出一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，以消除不同时相图像间存在的整体亮度差异和变化部分对配准精度及稳定性造成的影响，提高配准的精度和稳定性。为实现上述目的，本发明的技术方案包括如下步骤:(I)输入在两个时相获取的同一地区的多光谱图像集IdP I2，其中，时相I的波段序列图像集I1=IA11I，时相2的波段序列图像集I2= {A2b}，Atb为两个图像集中的每一幅单波段图像，其中，上标b表示波段序号，b=l, 2，…，B, B为总波段数，下标t为时相序号，t={l, 2}，每一幅单波段图像Atb均由η行η列像素构成，η的取值为32的倍数；(2)对两时相多光谱图像集I1和I2分别进行窗口大小为3X3像素的中值滤波去噪,并归一化处理,得到归一化图像集IjPf2;
(3)将归一化图像集Ijpf2中的每一幅图像均进行多层二维离散Haar小波变换，
得到第I层的低频系数矩阵为AP、垂直方向的高频系数矩阵为V、水平方向的高频系数
矩阵为对角线方向的高频系数矩阵为D:,其中分解层数1=1，...，L，L为最大分解层数，其取值范围为I 4 ;(4)将时相I和时相2的所有B个波段图像的第L层垂直方向的高频系数矩阵V:、
水平方向的高频系数矩阵和对角线方向的高频系数矩阵Dfi取绝对值后空间位置对应元素相加，分别得到时相I和时相2的高频系数图像HVD1和HVD2 ；(5)对高频系数图像HVD1和HVD2采用最大类间方差法分别进行阈值分割，得到时相I和时相2的边缘图像0和 2;(6)对边缘图像G和02，统计每一条不与其他边缘有8邻域连接的边缘线段在8邻域下所包含的像素个数，并将其中像素个数小于长度阈值Tc的边缘线段删除，对应得到
时相I和时相2的强边缘图像C1和C2，其中:
权利要求
1.一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，包括步骤如下: (1)输入在两个时相获取的同一地区的多光谱图像集I1和I2，其中，时相I的波段序列图像集I1=IA11I，时相2的波段序列图像集I2= {A2b}，Atb为两个图像集中的每一幅单波段图像，其中，上标b表示波段序号，b=l, 2，…，B, B为总波段数,下标t为时相序号，t={l, 2},每一幅单波段图像Atb均由η行η列像素构成，η的取值为32的倍数； (2)对两时相多光谱图像集I1和I2分别进行窗口大小为3X3像素的中值滤波去噪，并归一化处理，得到归一化图像集IjPf2; (3)将归一化图像集^和12中的每一幅图像均进行多层二维离散Haar小波变换，得到第I层的低频系数矩阵为Ag、垂直方向的高频系数矩阵为Vfj、水平方向的高频系数矩阵为对角线方向的高频系数矩阵为1^，其中分解层数1=1，...，L，L为最大分解层数，其取值范围为I 4 ; (4)将时相I和时相2的所有B个波段图像的第L层垂直方向的高频系数矩阵%、水平方向的高频系数矩阵Hf,、对角线方向的高频系数矩阵私取绝对值后空间位置对应元素相加，分别得到时相I和时相2的高频系数图像HVD1和HVD2 ； (5)对高频系数图像HVD1和HVD2采用最大类间方差法分别进行阈值分割，得到时相I和时相2的边缘图像 和0 (6)对边缘图像C,和62,统计每一条不与其他边缘有8邻域连接的边缘线段在8邻域下所包含的像素个数，并将其中像素个数小于长度阈值Tc的边缘线段删除，对应得到时相I和时相2的强边缘图像C1和C2，其中%:1.IUu (7)将强边缘图像C1和C2进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si (a, q，P)，得到边缘对齐度集合SI= {Si (a，q，P) }，其中a为旋转角度，q为横移参量，P为纵移参量，a=0, 1,2, -,180, q=0, 1,2,…，n/2L+l，p=0，1,2,...，n/2L+l ； (8)用边缘对齐度集合SI中的最大值所对应的旋转角度aM、横移参量qM和纵移参量pM对时相I的各波段归一化图像I/均进行相应的旋转和平移，得到时相I的各波段粗配准图像ΑΛ将时相I的所有波段粗配准图像的集合与步骤(2)得到的时相2的归一化图像集 2合并，得到粗配准图像集合I; (9)对粗配准图像集合I中两个时相各波段图像分别采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集合R ; (10)对时相I的各波段相对地物光谱反射率图像R113,分别进行I层二维离散Haar小波变换，得到对应波段的低频系数矩阵ARb、垂直方向的高频系数矩阵VRb、水平方向的高频系数矩阵HRb、对角线方向的高频系数矩阵DRb ； (11)将时相I各波段低频系 数矩阵ARb中的元素全部置为零，并与高频系数矩阵VRb、HRb和DRb进行二维离散Haar小波逆变换，得到时相I的各波段重构图像； (12)将时相I所有B个波段的重构图像取绝对值后对应像素值相加，得到一幅模糊边缘图像； (13)对时相I的模糊边缘图像采用最大类间方差法进行阈值分割，得到粗配准边缘图像CR; (14)将粗配准边缘图像CR从第I行第I列像素开始，由左至右、由上至下分为互不重叠且大小均为K列K行的uXu个子图像CRij，其中i和j分别为粗配准边缘子图像的行块序号和列块序号，i=l，2，...，u，j=l，2，…，u，u为粗配准边缘图像CR分块后的水平方向和垂直方向的子图像数目，且满足n=KX U，K为偶数，40彡K彡100 ； (15)将时相I的各波段相对地物光谱反射率图像R113中与粗配准边缘子图像CRij相同空间位置范围的图像块记为时相I的各波段相对地物光谱反射率子图像.并将时相I的所有波段相对地物光谱反射率子图像丨丨与时相2的所有波段相对地物光谱反射率图像IR2bI采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R2}与时相I的相对地物光谱反射率子图像{RUj}相匹配的光谱反射率曲线总差异量DRCiP匹配旋角参量Θ mi」、匹配横移参量Qmij以及匹配纵移参量Pmij ； (16)将时相2的相对地物光谱反射率图像集{RJ与时相I的相对地物光谱反射率子图像{R j}相匹配的光谱反射率曲线总差异量DRCij按行块序号和列块序号的顺序作为对应行列的元素构成矩阵，利用最大类间方差法对该矩阵进行阈值分割，得到差异二值矩阵Duc ; (17)分别统计差异二值矩阵du。中O值元素对应的匹配旋角参量emu、匹配水平参量Qmij以及匹配垂直参量Pmij的数目，将各参量的数目最大值记为最优旋角参量Θ mm、最优水平参量Q_以及最优垂直参量P_ ； (18)用最优旋角参量Θ_、最优水平参量Q_以及最优垂直参量P_对粗配准图像集合！中时相2的各波段归一化图像进行相应的旋转和平移，得到时相2的各波段最终配准图像 (19)将时相2的所有波段最终配准图像{A/j与步骤(8)中得到的时相I的所有波段粗配准图像丨A。合并，构成最终配准图像集合L
2.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(7)所述的将强边缘图像C1和C2进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si (a, q, p),通过如下步骤进行: (7a)将强边缘图像C1和C2相匹配的旋转角度a、横移参量q和纵移参量p的初始值均赋为O ;将强边缘图像C1的行序号和列序号均为ζ的像素点记为Ce，以像素点Ce为中心点O，水平方向为X轴、垂直方向为Y轴建立XOY直角坐标系，其中ζ=η2 +1 ； (7b)将强边缘图像C1在XOY坐标系中以O点为中心逆时针旋转a-90度，采用最近邻线性插值方法进行图像插值，得到旋转图像Cla ； (7c)在XOY坐标系中，对横坐标和纵坐标范围均在[-2 ζ，2 ζ ]内且旋转图像Cla以外的像素点的值赋为0，得到大小为(4 ζ +1) X (4 ζ +1)像素的旋转扩展图像CEla ； (7d)将强边缘图像C2与旋转扩展图像CEla的第1+p行至第η/2ι+ρ行、第Ι+q列至第n/2L+q列范围内的像素进行对应空间位置的像素值与运算，将旋转扩展图像CEla的行序号在[Ι+ρ,η/^+ρ]范围外且列序号在[l+q’n/^+q]范围外的像素值赋为O，得到一幅与旋转扩展图像CEla大小相同的边缘交集图像CAaqp ； (7e)在边缘交集图像CAaqp中，对任意一条不与其他边缘有8邻域连接的边缘线段，计算该边缘线段中包含的像素个数Lc与长度阈值Tc的比值Lc/Tc，若该比值小于1，则将该边缘线段包含的所有像素的值都赋为1，否则，赋值为Lc/Tc ； (7f)对边缘交集图像CAaqp中的每一条独立的边缘线段都执行步骤(7e)，得到增强边缘交集图像CA' aqp ； (7g)将增强边缘交集图像CA' _中所有像素的值求和，得到边缘对齐度Si (a, q, p)； (7h)若横移参量q〈2 ζ +1，则纵移参量P保持不变，对横移参量q的值加1，返回步骤(7d);若横移参量q=2(+l且纵移参量ρ〈2ζ+1，则将横移参量q的值置为O，纵移参量p的值加1，返回步骤(7d);若横移参量q=2ζ+l，纵移参量p=2ζ+l且旋转角度a〈180，则将旋转角度a的值加3，横移参量q和纵移参量P的值都置为O，返回步骤(7b);若横移参量q=2 ζ +1，纵移参量ρ=2 ζ +1，且旋转角度a=180，则强边缘图像C1和C2的旋转和平移匹配完成，得到所有的边缘对齐度{Si(a，q，p)}。
3.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(9)所述的对粗配准图像 集合I中两个时相各波段图像分别采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集合R，通过如下步骤进行: (9a)将时相I的第b波段粗配准图像A,中像素点(x，y)的灰度值记为/ΤΛυ),再将转换为相对地物光谱反射率值Ab(Xj)的对数:
4.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(15)所述的将时相I的所有波段相对地物光谱反射率子图像与时相2的所有波段相对地物光谱反射率图像{R2b}采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R2}与时相I的相对地物光谱反射率子图像{Rnj}相匹配的光谱反射率曲线总差异量DRCij、匹配旋角参量Qmij、匹配横移参量Qmij以及匹配纵移参量Pmu，通过如下步骤进行: (15a)将粗配准边缘子图像CRij的行块序号i和列块序号j初始值赋为2，CRu的中心点为该子图像的第V行第V 列的像素，其中V=Ceil (K/2)，ceil (.)表示取大于或等于括号内数值的最小整数，K为粗配准边缘子图像CRij的总行数或总列数； (15b)若粗配准边缘子图像CRij中像素值为I的边缘像素点的总数小于λ ΧΚ，执行步骤(15c)，否则，将时相2的相对地物光谱反射率图像集{R2}与时相I的相对地物光谱反射率子图像{RUj}相匹配的旋角参量Θ、水平参量Q以及垂直参量P初始值赋为O,执行步骤(15d)，其中边缘约束系数λ的取值范围为2彡λ < 4 ; (15c)若列块序号j〈u-l时，则行块序号i保持不变，列块序号j的值加1，返回步骤(15b);若列块序号j=u-l且行块序号i〈u_l时，则行块序号i的值加1，列块序号j的值赋为2，返回步骤(15b);若列块序号j=u-l且行块序号i=u-l时，则将所有行块序号i=l或i=u-l,或列块序号j=l或j=u-l的粗配准边缘子图像CRij的光谱反射率曲线总差异量DRCij置为空值，执行步骤(16); (15d)将时相2的各波段相对地物光谱反射率图像R2b中与粗配准边缘子图像CRij中心点空间位置对应的像素点记为Ro，其横坐标为(1-l)XK+v，纵坐标为(j-l)XK+v;以像素点Ro为中心点，分别将时相2的各波段相对地物光谱反射率图像R2b逆时针旋转Θ -5度，采用最近邻线性插值方法进行图像插值，得到各波段的相对地物光谱反射率旋转图像在该图像中将以像素点Ro为中心、大小为(3v+l) X (3v+l)像素的图像块作为时相2的各波段窗图像 (15e)分别将时相I的各波段相对地物光谱反射率图像R113中与粗配准边缘子图像CRij相同空间位置范围的图像块作为时相I的各波段搜索图像64; (15f)将时相I的各波段搜索图像中的第I行第I列像素点与时相2的各波段窗图像中的第1+Q行第1+P列像素点空间位置对应，则时相I的各波段搜索图像^^中的像素点(X' ,1')与时相2的各波段窗图像中的像素点(X' +Q，y' +P)空间位置对应，计算时相I的各波段搜索图像Ι 〗,与时相2的各波段窗图像食 ,之间的光谱反射率曲线差异量 DRC(0，Q，p):
全文摘要
本发明公开了一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，主要解决现有技术对存在变化区域和整体亮度差异的两时相多光谱图像集进行配准时稳定性不足、精度低的问题。其实现步骤为对输入两时相多光谱图像集进行预处理后的同时相的小波高频系数矩阵相加并分割，得到两时相强边缘图像；计算两时相强边缘图像的边缘对齐度的配准参数，以对两时相图像集进行配准，得到粗配准图像集；将粗配准图像集转换为相对地物光谱反射率图像集，对该图像集中的时相1的小波高频系数重构得到粗配准边缘图像，用该图像对时相1和时相2的相对地物光谱反射率图像进行局部匹配，得到最终配准图像集。本发明可用于变化检测、融合、及镶嵌。
文档编号G06T5/00GK103198483SQ20131011788
公开日2013年7月10日 申请日期2013年4月7日 优先权日2013年4月7日
发明者王桂婷, 焦李成, 孙一博, 公茂果, 钟桦, 王爽, 张小华 申请人:西安电子科技大学