CN103198483B - 基于边缘和光谱反射率曲线的多时相遥感图像配准方法 - Google Patents

Abstract

本发明公开了一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，主要解决现有技术对存在变化区域和整体亮度差异的两时相多光谱图像集进行配准时稳定性不足、精度低的问题。其实现步骤为：对输入两时相多光谱图像集进行预处理后的同时相的小波高频系数矩阵相加并分割，得到两时相强边缘图像；计算两时相强边缘图像的边缘对齐度的配准参数，以对两时相图像集进行配准，得到粗配准图像集；将粗配准图像集转换为相对地物光谱反射率图像集，对该图像集中的时相1的小波高频系数重构得到粗配准边缘图像，用该图像对时相1和时相2的相对地物光谱反射率图像进行局部匹配，得到最终配准图像集。本发明可用于变化检测、融合、及镶嵌。

Description

基于边缘和光谱反射率曲线的多时相遥感图像配准方法

技术领域

本发明属于光学遥感图像处理技术领域，涉及同一多光谱传感器在不同时间获取的相同地理位置的遥感图像间存在较大平移和旋转错位情况下的配准，可用于多时相多光谱遥感图像变化检测、融合、镶嵌等方面的前期配准工作。

背景技术

近年来，随着多光谱遥感图像数据的公开化程度和处理技术水平的不断提高，多光谱遥感图像数据已经被广泛用于城市规划、资源调查、自然灾害预防、农作物长势监测、军事目标探查等多种应用领域。

多光谱遥感图像的广泛应用首先需要进行图像数据的几何校正、辐射校正、配准、去噪、分割、镶嵌、融合、目标检测、变化检测等多项数字图像处理技术。其中，在进行多光谱遥感图像变化检测、融合以及镶嵌等工作前，首先需要对不同传感器获取的同一地区的图像或相同传感器在不同时间获取的同一地区的图像进行配准，以消除因获取图像的时间、角度、环境和传感器成像机理的不同造成图像间的平移、旋转、伸缩及局部形变问题。

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

根据特征空间的不同，现有的多光谱遥感图像配准方法通常可以分为基于像素灰度特征和基于图像几何特征两大类。大部分基于像素灰度特征的方法一般直接利用整幅图像的灰度信息，通过建立某种像素间的相似性测度来衡量两幅图像重叠部分地表反射属性的匹配程度，进而寻找到最优匹配时的平移、旋转和伸缩等变换参数。常见的灰度相似性测度如平均平方差MSD、平均绝对差MAD、归一化互相关NCC等通常要求待配准图像对应像素的灰度值具有线性关系，因而在利用它们进行不同时相的图像配准时容易受到噪声和图像间整体亮度差异的影响，而对于不同传感器或不同波段图像的配准则基本不适用。最初作为多模态医学图像配准相似性测度的最大互信息准则由于不要求待配准图像对应像素的灰度信息具有严格的线性关系，而认为只要两幅图像中主要的几何结构特征对准时即达到最佳匹配，因此已被用于不同时相和不同波段图像间的配准。在对不同时相的多光谱遥感图像的配准中，由于基于像素灰度特征的配准方法在搜索最佳匹配变换参数时通常是由全体像素共同决定图像间的对准程度，计算量大，并且难免会受到图像中不同位置大小的云、阴影及变化地物的干扰，造成配准精度和稳定性的下降。

基于图像几何特征的配准方法则首先需要对待配准图像提取角点、边缘等明显的几何特征，然后进行对应特征对或特征集之间的搜索匹配，进而得到图像变换参数。基于边缘特征的配准方法常见的如利用一些传统空域边缘检测算子进行边缘特征提取，但往往由于提取的边缘过于琐碎，不能很好地体现图像中的结构特征，对后续匹配精度和速度造成影响，因而需要对提取到的边缘进行优化处理。近年来，小波变换wavelet、曲波变换curvelet和轮廓波变换contourlet等多尺度几何分析方法由于在描述和分析图像中的边缘特征时具有优异的性能，也常被用于边缘特征的提取。1997年J.W.Hsieh.H.Y.M.Liao、K.C.Fan等在文献”Imageregistration using a new edge-based approach”（Computer Vision and ImageUnderstanding,1997,vol.67,no.2,pp.112-130）中首次提出利用小波变换得到系数的模局部极大值进行边缘特征点提取，然后以互相关作为相似性测度进行边缘特征点对的搜索匹配，其中还使用了一种非迭代的策略消除误匹配点对，达到了较好的配准效果。2008年孙淑一.吴勇和吴建民在文献”一种基于边缘特征的图像配准方法”（计算机工程与应用,2008,vol.44,no.7,pp.94-96）中、2010年孙雅琳和王菲在文献”基于边缘和互信息的红外与可见光图像配准”（电子科技,2010,vol.23,no.4,pp.80-82）中也都采用了相同的边缘特征提取方法，不同的是，这两种方法在搜索匹配时前者计算了两幅边缘图像重叠部分的交互方差，而后者计算了两幅边缘图像重叠部分的归一化互信息。如果将这两种方法用于不同时相多光谱遥感图像的配准，则容易受到图像中变化部分对应特征点的干扰而造成误配准。2009年陈志刚、尹福昌和孙孚在文献”基于非采样Contourlet变换高分辨率遥感图像配准”（光学学报,2009,vol.29,no.10,pp.2744-2750）中提出首先由图像通过非下采样Contourlet变换得到的各方向子带计算得到梯度矢量模，阈值分割得到边缘特征点后由归一化互相关匹配搜索匹配控制点，并利用概率支持算法消除误匹配点对，降低了变化部分对应特征点错误匹配的概率。总体来说，由于不直接采用像素的灰度值作为相似性测量的对象，因而该类方法受不同时相间的亮度差异、不同传感器以及不同波段图像间的成像特性差异的影响较小，适合多时相多光谱遥感图像的配准，但要求待配准图像中含有比较明显的边缘信息，否则提取不到理想的边缘将很可能导致配准失败。另外，其中边缘特征提取和匹配方法的准确性及稳定性也是影响其配准性能的主要因素。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，以消除不同时相图像间存在的整体亮度差异和变化部分对配准精度及稳定性造成的影响，提高配准的精度和稳定性。

为实现上述目的，本发明的技术方案包括如下步骤：

（1）输入在两个时相获取的同一地区的多光谱图像集I₁和I₂，其中，时相1的波段序列图像集I₁={A₁ ^b}，时相2的波段序列图像集I₂={A₂ ^b}，A_t ^b为两个图像集中的每一幅单波段图像，其中，上标b表示波段序号，b=1,2,…,B，B为总波段数，下标t为时相序号，t={1,2}，每一幅单波段图像A_t ^b均由n行n列像素构成，n的取值为32的倍数；

（2）对两时相多光谱图像集I₁和I₂分别进行窗口大小为3×3像素的中值滤波去噪，并归一化处理，得到归一化图像集和

（3）将归一化图像集和中的每一幅图像均进行多层二维离散Haar小波变换，得到第l层的低频系数矩阵为垂直方向的高频系数矩阵为水平方向的高频系数矩阵为对角线方向的高频系数矩阵为其中分解层数l=1,…,L，L为最大分解层数，其取值范围为1～4；

（4）将时相1和时相2的所有B个波段图像的第L层垂直方向的高频系数矩阵水平方向的高频系数矩阵和对角线方向的高频系数矩阵取绝对值后空间位置对应元素相加，分别得到时相1和时相2的高频系数图像HVD₁和HVD₂；

（5）对高频系数图像HVD₁和HVD₂采用最大类间方差法分别进行阈值分割，得到时相1和时相2的边缘图像和

（6）对边缘图像和统计每一条不与其他边缘有8邻域连接的边缘线段在8邻域下所包含的像素个数，并将其中像素个数小于长度阈值Tc的边缘线段删除，对应得到时相1和时相2的强边缘图像C₁和C₂，其中

（7）将强边缘图像C₁和C₂进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si(a,q,p)，得到边缘对齐度集合SI={Si(a,q,p)}，其中a为旋转角度，q为横移参量，p为纵移参量，a=0,1,2,…,180，q=0,1,2,…,n/2^L+1，p=0,1,2,…,n/2^L+1；

（8）用边缘对齐度集合SI中的最大值所对应的旋转角度a_M、横移参量q_M和纵移参量p_M对时相1的各波段归一化图像均进行相应的旋转和平移，得到时相1的各波段粗配准图像将时相1的所有波段粗配准图像的集合与步骤（2）得到的时相2的归一化图像集合并，得到粗配准图像集合

（9）对粗配准图像集合中两个时相各波段图像分别采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集合R；

（10）对时相1的各波段相对地物光谱反射率图像R₁ ^b，分别进行1层二维离散Haar小波变换，得到对应波段的低频系数矩阵AR^b、垂直方向的高频系数矩阵VR^b、水平方向的高频系数矩阵HR^b、对角线方向的高频系数矩阵DR^b；

（11）将时相1各波段低频系数矩阵AR^b中的元素全部置为零，并与高频系数矩阵VR^b、HR^b和DR^b进行二维离散Haar小波逆变换，得到时相1的各波段重构图像；

（12）将时相1所有B个波段的重构图像取绝对值后对应像素值相加，得到一幅模糊边缘图像；

（13）对时相1的模糊边缘图像采用最大类间方差法进行阈值分割，得到粗配准边缘图像CR；

（14）将粗配准边缘图像CR从第1行第1列像素开始，由左至右、由上至下分为互不重叠且大小均为K列K行的u×u个子图像CRi_j，其中i和j分别为粗配准边缘子图像的行块序号和列块序号，i=1,2,…,u，j=1,2,…,u，u为粗配准边缘图像CR分块后的水平方向和垂直方向的子图像数目，且满足n=K×u，K为偶数，40≤K≤100；

（15）将时相1的各波段相对地物光谱反射率图像R₁ ^b中与粗配准边缘子图像CRi_j相同空间位置范围的图像块记为时相1的各波段相对地物光谱反射率子图像并将时相1的所有波段相对地物光谱反射率子图像与时相2的所有波段相对地物光谱反射率图像{R₂ ^b}采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij、匹配旋角参量θm_ij、匹配横移参量Qm_ij以及匹配纵移参量Pm_ij；

（16）将时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij按行块序号和列块序号的顺序作为对应行列的元素构成矩阵，利用最大类间方差法对该矩阵进行阈值分割，得到差异二值矩阵D_uc；

（17）分别统计差异二值矩阵D_uc中0值元素对应的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij的数目，将各参量的数目最大值记为最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm；

（18）用最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm对粗配准图像集合中时相2的各波段归一化图像进行相应的旋转和平移，得到时相2的各波段最终配准图像

（19）将时相2的所有波段最终配准图像与步骤（8）中得到的时相1的所有波段粗配准图像合并，构成最终配准图像集合

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

a、本发明与现有的多时相多光谱遥感图像配准方法相比，由于采用了多层小波变换并进行多波段合并，能得到分辨率低，且综合所有波段边缘的强边缘图像，避免了单波段高分辨率图像对几何特征的稳定性不足；同时由于本发明采用了先进行两幅低分辨率的强边缘图像的配准，后进行不同时相多波段图像的光谱反射率曲线差异量局部配准的两步配准策略，在有效地减小局部配准搜索次数的同时，不仅避免了不同时相的图像集之间存在的变化部分的干扰，而且提高了配准精度、速度及稳定性。

b、本发明与现有基于边缘特征的多时相多光谱遥感图像配准方法相比，由于在利用边缘特征进行配准的基础上，设计并采用了一种光谱反射率曲线匹配的局部图像块搜索匹配方法，进一步提高了配准精度，避免了不同时相图像间存在的整体亮度差异的干扰。

附图说明

图1是本发明的整体流程图；

图2是本发明的3层二维离散小波变换示意图；

图3是本发明的旋转扩展图像示意图；

图4是本发明的实验图像；

图5是本发明的实验结果图。

具体实施方式

参照图1，本发明的实现步骤如下：

步骤1，输入在两个时相获取的同一地区的待配准多光谱图像集I₁和I₂。

输入的待配准的多光谱图像集为I₁和I₂，其中，时相1的波段序列图像集I₁={A₁ ^b}，时相2的波段序列图像集I₂={A₂ ^b}，A_t ^b为两个图像集中的每一幅单波段图像，其中，上标b表示波段序号，b=1,2,…,B，B为总波段数，下标t为时相序号，t={1,2}，每一幅单波段图像A_t ^b均由n行n列像素构成，n的取值为32的倍数。

步骤2，对两时相多光谱图像集I₁和I₂分别进行中值滤波去噪，并归一化处理，得到归一化图像集和

2.1）将两时相多光谱图像集I₁和I₂中每一幅图像A_t ^b的灰度值区间从0～255转换到0～1之间，然后采用窗口大小为3×3像素的中值滤波进行去噪，得到相应的去噪后图像

2.2）对去噪后图像进行归一化处理，得到归一化图像

${\hat{A}}_t^b=\left\{{\hat{A}}_t^b(x,y)\right|x=\mathrm{1,2},...,n,y=\mathrm{1,2},...,n\},$

其中，为去噪后图像中的像素点(x,y)的灰度值，其计算公式如下：

${\hat{A}}_t^b(x,y)=\frac{{\stackrel{\·\·}{A}}_t^b(x,y)-\min \left(\right\{{\stackrel{\·\·}{A}}_t^b(x,y)\left\}\right)}{\max \left(\right\{{\stackrel{\·\·}{A}}_t^b(x,y)\left\}\right)-\min \left(\right\{{\stackrel{\·\·}{A}}_t^b(x,y)\left\}\right)},$

表示取去噪后图像中所有灰度值中的最大值，表示取去噪后图像中所有灰度值中的最小值，x和y分别为去噪后图像中像素的行序号和列序号；

2.3）用时相1和时相2的归一化图像分别构成归一化图像集和即

${\hat{I}}_1=\{{{\hat{A}}_1}^1,{{\hat{A}}_1}^2,\·\·\·{{\hat{A}}_1}^b,\·\·\·,{{\hat{A}}_1}^B\},$

${\hat{I}}_2=\{{{\hat{A}}_2}^1,{{\hat{A}}_2}^2,\·\·\·,{{\hat{A}}_2}^b,\·\·\·,{{\hat{A}}_2}^B\}.$

步骤3，将归一化图像集和中的每一幅图像均进行二维离散小波变换。

如图2所示，将归一化图像集和中的每一幅图像均进行多层二维离散Haar小波变换，得到第l层的低频系数矩阵为垂直方向的高频系数矩阵为水平方向的高频系数矩阵为对角线方向的高频系数矩阵为其中分解层数l=1,…,L，L为最大分解层数，其取值范围为1～4，本发明实例中取值为3。

由于小波变换对图像中的边缘信息具有良好的描述性，且不会像传统基于梯度的边缘检测算子产生过于琐碎的边缘线，因此本发明采用小波变换得到的高频系数进行处理以获取图像中的主要边缘。然而，当图像尺寸较大时，对整幅图像进行配准的搜索量是非常大的，为了降低搜索耗时，本发明采用了L=3层二维离散Haar小波对图像进行变换，使变换得到的最高层小波系数矩阵只有原图像大小的 $\frac{1}{2^L}\×\frac{1}{2^L}=\frac{1}{64}.$

步骤4，将时相1和时相2的所有B个波段图像的第L层垂直方向的高频系数矩阵为水平方向的高频系数矩阵为和对角线方向的高频系数矩阵为

取绝对值后对应元素相加，分别得到时相1和时相2的高频系数图像HVD₁和HVD₂：

${\mathrm{HVD}}_1=\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\left(\right|H_{1L}^b|+|V_{1L}^b|+|D_{1L}^b\left|\right),$

${\mathrm{HVD}}_2=\underset{b=1}{\overset{B}{\mathrm{\Σ}}}\left(\right|H_{2L}^b|+|V_{2L}^b|+|D_{2L}^b\left|\right),$

其中，|·|表示取绝对值操作。

通过小波变换得到的垂直方向的高频系数矩阵、水平方向的高频系数矩阵和对角线方向的高频系数矩阵对不同方向边缘线段的敏感程度不同，而不同波段的图像反映不同类型地物边界的效果也不同。为了增强不同方向和不同地物边界类型的边缘信息，按照如下公式计算得到高频系数图像；

步骤5，对高频系数图像HVD₁和HVD₂分别采用最大类间方差法进行阈值分割，得到时相1的边缘图像和时相2的边缘图像

步骤6，对边缘图像和统计每一条不与其他边缘有8邻域连接的边缘线段在8邻域下所包含的像素个数，将其中像素个数小于长度阈值Tc的边缘线段删除，得到时相1和时相2的强边缘图像C₁和C₂，其中强边缘图像C₁和C₂的总行数和总列数均为n2^L。

步骤7，将强边缘图像C₁和C₂进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si(a,q,p)，得到边缘对齐度集合SI={Si(a,q,p)}。

7.1）将强边缘图像C₁和C₂相匹配的旋转角度a、横移参量q和纵移参量p的初始值均赋为0；将强边缘图像C₁的行序号和列序号均为ζ的像素点记为Ce，以像素点Ce为中心点O，水平方向为X轴、垂直方向为Y轴建立XOY直角坐标系，其中ζ=n/2^L+1，a=0,1,2,…,180，q=0,1,2,…,n/2^L+1，p=0,1,2,…,n/2^L+1；

7.2）将强边缘图像C₁在XOY坐标系中以O点为中心逆时针旋转a-90度，采用最近邻线性插值方法进行图像插值，得到旋转图像C_1a，则旋转图像C_1a的横坐标和纵坐标范围均为 $[-\mathrm{round}\left(\sqrt{2}(\mathrm{\ζ}-1)\cos (a-135)\right),\mathrm{round}\left(\sqrt{2}(\mathrm{\ζ}-1)\cos (a-135)\right)],$ 这里，round(·)表示对括号内的数值四舍五入；

7.3）在XOY坐标系中，如图3所示，对横坐标和纵坐标范围均在[-2ζ,2ζ]内且旋转图像C_1a以外的像素点的值赋为0，得到大小为(4ζ+1)×(4ζ+1)像素的旋转扩展图像CE_1a；

7.4）将强边缘图像C₂与旋转扩展图像CE_1a的第1+p行至第n/2^L+p行、第1+q列至第n/2^L+q列范围内的像素进行对应空间位置的像素值与运算，将旋转扩展图像CE_1a的行序号在[1+p,n/2^L+p]范围外且列序号在[1+q,n/2^L+q]范围外的像素值赋为0，得到一幅与旋转扩展图像CE_1a相同大小的边缘交集图像CA_aqp；

7.5）在边缘交集图像CA_aqp中，对任意一条不与其他边缘有8邻域连接的边缘线段，计算该边缘线段中包含的像素个数Lc与长度阈值Tc的比值Lc/Tc，若该比值小于1，则将该边缘线段包含的所有像素的值都赋为1，否则，赋值为Lc/Tc；

7.6）对边缘交集图像CA_aqp中的每一条独立的边缘线段都执行步骤7.5），得到增强边缘交集图像CA′_aqp；

7.7）将增强边缘交集图像CA′_aqp中所有像素的值求和，得到边缘对齐度Si(a,q,p)；

7.8）若横移参量q<2ζ+1，则纵移参量p保持不变，对横移参量q的值加1，返回步骤7.4）；若横移参量q=2ζ+1且纵移参量p<2ζ+1，则将横移参量q的值置为0，纵移参量p的值加1，返回步骤7.4）；若横移参量q=2ζ+1，纵移参量p=2ζ+1且旋转角度a<180，则将旋转角度a的值加3，横移参量q和纵移参量p的值都置为0，返回步骤7.2）；若横移参量q=2ζ+1，纵移参量p=2ζ+1，且旋转角度a=180，则强边缘图像C₁和C₂的旋转和平移匹配完成，得到所有的边缘对齐度{Si(a,q,p)}，由此得到边缘对齐度集合SI={Si(a,q,p)}；

7.9）计算边缘对齐度集合SI中的最大值，将该最大值所对应的旋转角度、横移参量和纵移参量分别记为a_M、q_M和p_M。

步骤8，用边缘对齐度集合SI中最大值所对应的旋转角度a_M、横移参量q_M和纵移参量p_M对时相1的各波段归一化图像均进行相应的旋转和平移，得到时相1的各波段粗配准图像将时相1的所有波段粗配准图像的集合与步骤2得到的时相2的归一化图像集合并，得到粗配准图像集合

8.1）将时相1的各波段归一化图像中的第n/2行第n/2列像素点均记为Co，以像素点Co为中心点O，并以水平方向为X轴、垂直方向为Y轴建立XOY直角坐标系；

8.2）将时相1的各波段归一化图像在XOY坐标系中以Co点为中心逆时针旋转a_M-90度，用最近邻线性插值方法进行图像插值，得到时相1的各波段归一化旋转图像

8.3）在XOY坐标系中，对横坐标和纵坐标范围均在[-n,n]内且归一化旋转图像以外的像素点的值赋为0，得到大小为(2n+1)×(2n+1)像素的各波段归一化旋转扩展图像

8.4）将时相1的各波段归一化旋转扩展图像中的第n-2^L×p_M行至第-2^L×p_M行、第-n+2^L×q_M列至第2^L×q_M列范围内的像素形成的n×n大小的图像块作为时相1的各波段粗配准图像

8.5）将时相1的所有波段粗配准图像的集合与步骤（2）得到的时相2的归一化图像集合并，得到粗配准图像集合

步骤9，对粗配准图像集合中两个时相的各波段图像采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集R。

9.1）将时相1的第b波段粗配准图像中像素点(x,y)的灰度值记为再将转换为相对地物光谱反射率值r₁ ^b(x,y)的对数：

$\log r_1^b(x,y)=\log {\tilde{A}}_1^b(x,y)-\log \frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y)-\log \frac{1}{n^2}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y);$

$+\log \frac{1}{n^2\&times;B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y)$

9.2）将时相2的第b波段归一化图像中像素点(x,y)的灰度值记为再将转换为相对地物光谱反射率值r₂ ^b(x,y)的对数：

$\log r_2^b(x,y)=\log {\tilde{A}}_2^b(x,y)-\log \frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\tilde{A}}_2^b(x,y)-\log \frac{1}{n^2}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_2^b(x,y);$

$+\log \frac{1}{n^2\&times;B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_2^b(x,y)$

9.3）分别对步骤（9.1）得到的logr₁ ^b(x,y)和步骤（9.2）得到的logr₂ ^b(x,y)做指数变换，得到该像素(x,y)在时相1和时相2的相对地物光谱反射率值r₁ ^b(x,y)和r₂ ^b(x,y)；

9.4）对第b波段的时相1的粗配准图像和时相2的归一化图像的所有像素点均重复步骤（9.1）至步骤（9.3），得到第b波段时相1的相对地物光谱反射率图像R₁ ^b={r₁ ^b(x,y)}和时相2的相对地物光谱反射率图像R₂ ^b={r₂ ^b(x,y)}；

9.5）对时相1的所有B个波段的粗配准图像重复步骤（9.1）至（9.4），分别得到时相1的相对地物光谱反射率图像R₁ ¹,R₁ ²,…,R₁ ^b,…,R₁ ^B，再将它们合并构成时相1的相对地物光谱反射率图像集R₁；对时相2的所有B个波段的归一化图像重复步骤（9.1）至步骤（9.4），分别得到时相2的相对地物光谱反射率图像R₂ ¹,R₂ ²,…,R₂ ^b,…,R₂ ^B，再将它们合并构成时相2的相对地物光谱反射率图像集R₂。

步骤10，对时相1的各波段相对地物光谱反射率图像R₁ ^b，分别进行1层二维离散Haar小波变换，得到对应波段的低频系数矩阵AR^b、的高频系数矩阵VR^b、水平方向的高频系数矩阵HR^b、对角线方向的高频系数矩阵DR^b。

步骤11，将时相1各波段低频系数矩阵AR^b中的元素全部置为零，并与高频系数矩阵VR^b、HR^b和DR^b进行二维离散Haar小波逆变换，得到时相1的各波段重构图像。

步骤12，将时相1所有B个波段的重构图像取绝对值后对应像素值相加，得到一幅模糊边缘图像。

步骤13，对时相1的模糊边缘图像采用最大类间方差法进行阈值分割，得到粗配准边缘图像CR。

步骤14，将粗配准边缘图像CR从第1行第1列像素开始，由左至右、由上至下分为互不重叠且大小均为K列K行的u×u个子图像CR_ij，其中i和j分别为粗配准边缘子图像的行块序号和列块序号，i=1,2,…,u，j=1,2,…,u，u为粗配准边缘图像CR分块后的水平方向和垂直方向的子图像数目，且满足n=K×u，K为粗配准边缘子图像CR_ij的总行数或总列数，K为偶数，40≤K≤100，本发明实例中K=80。

步骤15，将时相1的各波段相对地物光谱反射率图像R₁ ^b中选取与粗配准边缘子图像CR_ij相同空间位置范围的图像块记为时相1的各波段相对地物光谱反射率子图像并将时相1的所有波段相对地物光谱反射率子图像与时相2的所有波段相对地物光谱反射率图像{R₂ ^b}采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij、匹配旋角参量θm_ij、匹配横移参量Qm_ij以及匹配纵移参量Pm_ij。

相同的地物具有相似的光谱反射率曲线，而不同地物间的光谱反射率曲线则相差很大，故一般只有当不同的物体相邻时，才会在图像中产生明显的边缘。如果两幅图像之间具有较大的空间位置差异时，同一空间位置的像素由于在不同的时相上对应着不同的地物，从而使各像素的光谱反射率曲线差异量较大；反之，则同一位置的像素对应着同一地物，因此各像素的光谱反射率曲线差异量较小，甚至为0。因此，本发明提出了光谱反射率曲线差异量，通过光谱反射率曲线差异量的最小化，达到图像间的最优配准。具体步骤如下：

15.1）将粗配准边缘子图像CR_ij的行块序号i和列块序号j初始值赋为2，CR_ij的中心点为该图像的第v行第v列的像素，其中v=ceil(K/2)，ceil(·)表示取大于或等于括号内数值的最小整数；

15.2）若粗配准边缘子图像CR_ij中像素值为1的边缘像素点的总数小于λ×K，执行步骤（15.3），否则，将时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的旋角参量θ、水平参量Q以及垂直参量P初始值赋为0，执行步骤（15.4），其中边缘约束系数λ的取值范围为2≤λ≤4，本发明实例中λ=2.5；

15.3）若列块序号j<u-1时，则行块序号i保持不变，列块序号j的值加1，返回步骤（15.2）；若列块序号j=u-1且行块序号i<u-1时，则行块序号i的值加1，列块序号j的值赋为2，返回步骤（15.2）；若列块序号j=u-1且行块序号i=u-1时，则将所有行块序号i=1或i=u-1，或列块序号j=1或j=u-1的粗配准边缘子图像CR_ij的光谱反射率曲线总差异量DRC_ij置为空值，执行步骤16；

15.4）将时相2的各波段相对地物光谱反射率图像R₂ ^b中与粗配准边缘子图像CR_ij中心点空间位置对应的像素点记为Ro，其横坐标为(i-1)×K+v，纵坐标为(j-1)×K+v；以像素点Ro为中心点，分别将时相2的各波段相对地物光谱反射率图像R₂ ^b逆时针旋转θ-5度，采用最近邻线性插值方法进行图像插值，得到各波段的相对地物光谱反射率旋转图像在该图像中将以像素点Ro为中心、大小为(3v+1)×(3v+1)像素的图像块作为时相2的各波段窗图像；

15.5）分别将时相1的各波段相对地物光谱反射率图像R₁ ^b中与粗配准边缘子图像CR_ij相同空间位置范围的图像块作为时相1的各波段搜索图像

15.6）将时相1的各波段搜索图像中的第1行第1列像素点与时相2的各波段窗图像中的第1+Q行第1+P列像素点空间位置对应，则时相1的各波段搜索图像中的像素点(x′,y′)与时相2的各波段窗图像中的像素点(x′+Q,y′+P)空间位置对应，计算时相1的各波段搜索图像与时相2的各波段窗图像之间的光谱反射率曲线差异量DRC(θ,Q,P)：

$\mathrm{DRC}(\mathrm{\&theta;},Q,P)=\underset{x^{\&prime;}=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{y^{\&prime;}=1}{\overset{K}{\mathrm{\&Sigma;}}}{V}^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})\&CenterDot;M^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;});$

其中，V^θQP(x′,y′)为时相1的各波段搜索图像与时相2的各波段窗图像之间对应像素的反射率变化矢量的方差：

$V^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\left\{\right[{\stackrel{\&prime;}{R}}_{1\mathrm{ij}}^b(x^{\&prime;},y^{\&prime;})-{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P)]-\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}[{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P)\left]\right\}}^2---4)$

M^θQP(x′,y′)为时相1的各波段搜索图像与时相2的各波段窗图像之间对应像素的反射率变化矢量的模值：

$M^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})=\sqrt{\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}[{\stackrel{\&prime;}{R}}_{1\mathrm{ij}}^b(x^{\&prime;},y^{\&prime;})-{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P){]}^2};$

其中，x′和y′分别为时相1的各波段搜索图像中的像素的行序号和列序号，x′=1,2,…,K，y′=1,2,…,K；

15.7）若水平参量Q<3v+1-K，则垂直参量P保持不变，水平参量Q的值加1，返回步骤（15.6）；若水平参量Q=3v+1-K且垂直参量P<3v+1-K，则将水平参量Q的值置为0，垂直参量P的值加1，返回步骤（15.6）；若水平参量Q=3v+1-K且垂直参量P=3v+1-K且旋角参量θ<10，则将旋角参量θ的值加1，水平参量Q、垂直参量P的值都置为0，返回步骤（15.4）；若水平参量Q=3v+1-K、垂直参量P=3v+1-K且旋角参量θ=10，执行步骤（15.8）；

15.8）将计算得到的光谱反射率曲线差异量DRC(θ,Q,P)中的最小值作为时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像集{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij，将该最小值对应的旋角参量、水平参量和垂直参量分别作为时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij，返回步骤（15.3）循环计算不同行块序号或不同列块序号的时相1的相对地物光谱反射率子图像R_1ij与时相2的相对地物光谱反射率图像R₂的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij。

步骤16，时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij按行块序号和列块序号的顺序作为对应行列的元素构成矩阵，利用最大类间方差法对该矩阵进行阈值分割，得到差异二值矩阵D_uc。

步骤17，分别统计差异二值矩阵D_uc中0值元素对应的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij的数目，将各参量的数目最高者记为最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm。

17.1）统计差异二值矩阵D_uc中0值元素对应的匹配旋角参量θm_ij，将最优旋角参量θmm赋值为max({θm_ij|i=1,2,…,u，j=1,2,…,u})，其中max(·)表示取括号内数值的最大值；

17.2）统计差异二值矩阵D_uc中元素值为0且匹配旋角参量θm_ij=θmm对应的匹配水平参量Qm_ij和匹配垂直参量Pm_ij，将匹配水平参量Qm_ij中的最大值作为最优水平参量Qmm，将匹配垂直参量Pm_ij中的最大值作为最优垂直参量Pmm。

步骤18，用最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm对粗配准图像集合中时相2的各波段归一化图像进行相应的旋转和平移，得到时相2的各波段最终配准图像

18.1）将时相2的各波段归一化图像中第n/2行第n/2列像素点记为Co′，以该像素点为中心点O′，经点O′的水平方向为X′轴，垂直方向为Y′轴建立X′O′Y′直角坐标系；

18.2）在X′O′Y′坐标系中，对时相2的各波段归一化图像以O′点为中心逆时针旋转θmm-5度，用最近邻线性插值方法进行图像插值，得到时相2的各波段归一化旋转图像

18.3）在X′O′Y′坐标系中，对时相2的各波段归一化旋转图像对横坐标和纵坐标范围均在[-3n/4,3n/4]内且时相2的各波段归一化旋转图像以外的像素点的值赋为0，得到大小为(3n/2+1)×(3n/2+1)像素的各波段归一化旋转扩展图像

18.4）将时相2的各波段归一化扩展图像中的第-Pmm+3n/4行至第-Pmm-n/4行、第Qmm-3n/4列至第Qmm+n/4列范围内的像素形成的n×n大小的图像块作为时相2的最终配准图像

步骤19，将时相2的所有波段最终配准图像与步骤（8）中得到的时相1的所有波段粗配准图像合并，构成最终配准图像集合

本发明的效果可通过如下仿真实验结果进一步说明：

1.实验数据

本发明仿真实验所使用的实验数据均是从1999年11月20日和2001年10月8日拍摄的增强专题制图仪ETM图像集中选取的位于哥伦比亚比查达省地区中的一个区域，大小为800×800，灰度级为256级。增强专题制图仪图像集的第6、第8波段与其他波段分辨率不同，因此只选取了除第6、第8波段外的第1至第7共6个波段图像。其中，两个时相的图像集间存在较大的平移、旋转角度误差以及少量的云和阴影。为清晰直观展示两时相间的不同，将实验数据中的两个时相各6波段的图像集内的12幅单波段图像分别显示，如图4所示，其中：

图4（a1）为时相1的第1波段图像；

图4（a2）为时相1的第2波段图像；

图4（a3）为时相1的第3波段图像；

图4（a4）为时相1的第4波段图像；

图4（a5）为时相1的第5波段图像；

图4（a6）为时相1的第7波段图像；

图4（b1）为时相1的第1波段图像；

图4（b2）为时相1的第2波段图像；

图4（b3）为时相1的第3波段图像；

图4（b4）为时相1的第4波段图像；

图4（b5）为时相1的第5波段图像；

图4（b6）为时相1的第7波段图像。

2.实验内容与结果

用本发明的方法对上述实验数据进行配准，结果如图5，其中：

图5（a）为时相1的强边缘图像，该图中白色的像素表示为边缘，黑色的像素表示为背景；

图5（b）为时相2的强边缘图像，该图中白色的像素表示为边缘，黑色的像素表示为背景；

由于实验数据集中同一时相的各波段图像之间是同时采集的，在同一时相的各个波段图像之间是无需配准的，所以，同一波段不同时相之间的两幅图像的配准就完全可以反映出不同时相各个波段的图像集的配准。对图4所示的实验数据进行配准，

得到的时相1的粗配准图像集中的第3波段图像，如图5（c）所示；

得到的时相2的最终配准图像集中的第3波段图像，如图5（d）所示。

结论：通过上述实验结果图发现，本发明方法有效避免了不同时相图像间存在的整体亮度差异和变化部分对配准精度和稳定性的影响；通过对比粗配准图像与最终配准图像发现，基于边缘图像的粗配准过程大大降低了时相1与时相2原图像之间的平移和旋转角度误差，而通过光谱反射率曲线匹配得到的最终配准图像则将粗配准图像中剩余的配准误差有效地进一步降低，使配准精度达到了亚像素级。

Claims (4)

1.一种基于边缘和光谱反射率曲线的多时相遥感图像配准方法，包括步骤如下：

(1)输入在两个时相获的同一地区的多光谱图像集I₁和I₂，其中，时相1的波段序列图像集I₁＝{A₁ ^b}，时相2的波段序列图像集I₂＝{A₂ ^b}，A_t ^b为两个图像集中的每一幅单波段图像，其中，上标b表示波段序号，b＝1,2,…,B，B为总波段数，下标t为时相序号，t＝{1,2}，每一幅单波段图像A_t ^b均由n行n列像素构成，n的取值为32的倍数；

(2)对两时相多光谱图像集I₁和I₂分别进行窗口大小为3×3像素的中值滤波去噪，并归一化处理，得到归一化图像集和

(3)将归一化图像集和中的每一幅图像均进行多层二维离散Haar小波变换，得到第l层的低频系数矩阵为垂直方向的高频系数矩阵为水平方向的高频系数矩阵为对角线方向的高频系数矩阵为其中分解层数l＝1,…,L，L为最大分解层数，其取值范围为1～4；

(4)将时相1和时相2的所有B个波段图像的第L层垂直方向的高频系数矩阵水平方向的高频系数矩阵对角线方向的高频系数矩阵取绝对值后空间位置对应元素相加，分别得到时相1和时相2的高频系数图像HVD₁和HVD₂；

(5)对高频系数图像HVD₁和HVD₂采用最大类间方差法分别进行阈值分割，得到时相1和时相2的边缘图像和

(6)对边缘图像和统计每一条不与其他边缘有8邻域连接的边缘线段在8邻域下所包含的像素个数，并将其中像素个数小于长度阈值Tc的边缘线段删除，对应得到时相1和时相2的强边缘图像C₁和C₂，其中

(7)将强边缘图像C₁和C₂进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si(a,q,p)，得到边缘对齐度集合SI＝{Si(a,q,p)}，其中a为旋转角度，q为横移参量，p为纵移参量，a＝0,1,2,…,180，q＝0,1,2,…,n/2^L+1，p＝0,1,2,…,n/2^L+1；

(8)用边缘对齐度集合SI中的最大值所对应的旋转角度a_M、横移参量q_M和纵移参量p_M对时相1的各波段归一化图像均进行相应的旋转和平移，得到时相1的各波段粗配准图像将时相1的所有波段粗配准图像的集合与步骤(2)得到的时相2的归一化图像集合并，得到粗配准图像集合

(9)对粗配准图像集合中两个时相各波段图像分别采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集合R；

(10)对时相1的各波段相对地物光谱反射率图像R₁ ^b，分别进行1层二维离散Haar小波变换，得到对应波段的低频系数矩阵AR^b、垂直方向的高频系数矩阵VR^b、水平方向的高频系数矩阵HR^b、对角线方向的高频系数矩阵DR^b；

(11)将时相1各波段低频系数矩阵AR^b中的元素全部置为零，并与高频系数矩阵VR^b、HR^b和DR^b进行二维离散Haar小波逆变换，得到时相1的各波段重构图像；

(12)将时相1所有B个波段的重构图像取绝对值后对应像素值相加，得到一幅模糊边缘图像；

(13)对时相1的模糊边缘图像采用最大类间方差法进行阈值分割，得到粗配准边缘图像CR；

(14)将粗配准边缘图像CR从第1行第1列像素开始，由左至右、由上至下分为互不重叠且大小均为K列K行的u×u个子图像CR_ij，其中i和j分别为粗配准边缘子图像的行块序号和列块序号，i＝1,2,…,u，j＝1,2,…,u，u为粗配准边缘图像CR分块后的水平方向和垂直方向的子图像数目，且满足n＝K×u，K为偶数，40≤K≤100；

(15)将时相1的各波段相对地物光谱反射率图像R₁ ^b中与粗配准边缘子图像CR_ij相同空间位置范围的图像块记为时相1的各波段相对地物光谱反射率子图像并将时相1的所有波段相对地物光谱反射率子图像与时相2的所有波段相对地物光谱反射率图像{R₂ ^b}采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij、匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij；

(16)将时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij按行块序号和列块序号的顺序作为对应行列的元素构成矩阵，利用最大类间方差法对该矩阵进行阈值分割，得到差异二值矩阵D_uc；

(17)分别统计差异二值矩阵D_uc中0值元素对应的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij的数目，将各参量的数目最大值记为最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm；

(18)用最优旋角参量θmm、最优水平参量Qmm以及最优垂直参量Pmm对粗配准图像集合中时相2的各波段归一化图像进行相应的旋转和平移，得到时相2的各波段最终配准图像

(19)将时相2的所有波段最终配准图像与步骤(8)中得到的时相1的所有波段粗配准图像合并，构成最终配准图像集合

2.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(7)所述的将强边缘图像C₁和C₂进行旋转和平移匹配，并计算每次匹配的边缘对齐度Si(a,q,p)，通过如下步骤进行：

(7a)将强边缘图像C₁和C₂相匹配的旋转角度a、横移参量q和纵移参量p的初始值均赋为0；将强边缘图像C₁的行序号和列序号均为ζ的像素点记为Ce，以像素点Ce为中心点O，水平方向为X轴、垂直方向为Y轴建立XOY直角坐标系，其中ζ＝n/2^L+1；

(7b)将强边缘图像C₁在XOY坐标系中以O点为中心逆时针旋转a-90度，采用最近邻线性插值方法进行图像插值，得到旋转图像C_1a；

(7c)在XOY坐标系中，对横坐标和纵坐标范围均在[-2ζ,2ζ]内且旋转图像C_1a以外的像素点的值赋为0，得到大小为(4ζ+1)×(4ζ+1)像素的旋转扩展图像CE_1a；

(7d)将强边缘图像C₂与旋转扩展图像CE_1a的第1+p行至第n/2^L+p行、第1+q列至第n/2^L+q列范围内的像素进行对应空间位置的像素值与运算，将旋转扩展图像CE_1a的行序号在[1+p,n/2^L+p]范围外且列序号在[1+q,n/2^L+q]范围外的像素值赋为0，得到一幅与旋转扩展图像CE_1a大小相同的边缘交集图像CA_aqp；

(7e)在边缘交集图像CA_aqp中，对任意一条不与其他边缘有8邻域连接的边缘线段，计算该边缘线段中包含的像素个数Lc与长度阈值Tc的比值Lc/Tc，若该比值小于1，则将该边缘线段包含的所有像素的值都赋为1，否则，赋值为Lc/Tc；

(7f)对边缘交集图像CA_aqp中的每一条独立的边缘线段都执行步骤(7e)，得到增强边缘交集图像CA'_aqp；

(7g)将增强边缘交集图像CA'_aqp中所有像素的值求和，得到边缘对齐度Si(a,q,p)；

(7h)若横移参量q<2ζ+1，则纵移参量p保持不变，对横移参量q的值加1，返回步骤(7d)；若横移参量q＝2ζ+1且纵移参量p<2ζ+1，则将横移参量q的值置为0，纵移参量p的值加1，返回步骤(7d)；若横移参量q＝2ζ+1，纵移参量p＝2ζ+1且旋转角度a<180，则将旋转角度a的值加3，横移参量q和纵移参量p的值都置为0，返回步骤(7b)；若横移参量q＝2ζ+1，纵移参量p＝2ζ+1，且旋转角度a＝180，则强边缘图像C₁和C₂的旋转和平移匹配完成，得到所有的边缘对齐度{Si(a,q,p)}。

3.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(9)所述的对粗配准图像集合中两个时相各波段图像分别采用对数残差修正模型将像素的灰度值转换为相对地物光谱反射率值，得到相对地物光谱反射率图像集合R，通过如下步骤进行：

(9a)将时相1的第b波段粗配准图像中像素点(x,y)的灰度值记为再将转换为相对地物光谱反射率值r₁ ^b(x,y)的对数：

$\begin{array}{c}{{\log r}_1}^b(x,y)=\log {\tilde{A}}_1^b(x,y)-\log \frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y)-\log \frac{1}{n^2}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y)\\ +\log \frac{1}{n^2\&times;B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\tilde{A}}_1^b(x,y)\end{array}---1)$

其中x和y分别为各波段粗配准图像中像素的行序号和列序号，x＝1,2,…,n，y＝1,2,…,n；

(9b)将时相2的第b波段归一化图像中像素点(x,y)的灰度值记为再将转换为相对地物光谱反射率值r₂ ^b(x,y)的对数：

$\begin{array}{c}{{\log r}_2}^b(x,y)=\log {\hat{A}}_2^b(x,y)-\log \frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\hat{A}}_2^b(x,y)-\log \frac{1}{n^2}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\hat{A}}_2^b(x,y)\\ +\log \frac{1}{n^2\&times;B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{n}{\mathrm{\&Sigma;}}}\underset{x=1}{\overset{n}{\mathrm{\&Sigma;}}}{\hat{A}}_2^b(x,y)\end{array}---2)$

(9c)分别对步骤(9a)得到的logr₁ ^b(x,y)和步骤(9b)得到的logr₂ ^b(x,y)做指数变换，得到该像素(x,y)在时相1和时相2的相对地物光谱反射率值r₁ ^b(x,y)和r₂ ^b(x,y)；

(9d)对第b波段的时相1的粗配准图像和时相2的归一化图像的所有像素点均重复步骤(9a)至(9c)，得到第b波段时相1的相对地物光谱反射率图像R₁ ^b＝{r₁ ^b(x,y)}和时相2的相对地物光谱反射率图像R₂ ^b＝{r₂ ^b(x,y)}；

(9e)对时相1的所有B个波段的粗配准图像重复步骤(9a)至(9d)，分别得到时相1的相对地物光谱反射率图像R₁ ¹,R₁ ²,…,R₁ ^b,…,R₁ ^B，再将它们合并构成时相1的相对地物光谱反射率图像集R₁；对时相2的所有B个波段的归一化图像重复步骤(9a)至(9d)，分别得到时相2的相对地物光谱反射率图像R₂ ¹,R₂ ²,…,R₂ ^b,…,R₂ ^B，再将它们合并构成时相2的相对地物光谱反射率图像集R₂。

4.根据权利要求1所述的基于边缘和光谱反射率曲线的多时相遥感图像配准方法，其中步骤(15)所述的将时相1的所有波段相对地物光谱反射率子图像与时相2的所有波段相对地物光谱反射率图像{R₂ ^b}采用光谱反射率曲线差异最小方法进行匹配，得到时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij、匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij，通过如下步骤进行：

(15a)将粗配准边缘子图像CR_ij的行块序号i和列块序号j初始值赋为2，CR_ij的中心点为该子图像的第v行第v列的像素，其中v＝ceil(K/2)，ceil(·)表示取大于或等于括号内数值的最小整数，K为粗配准边缘子图像CR_ij的总行数或总列数；

(15b)若粗配准边缘子图像CR_ij中像素值为1的边缘像素点的总数小于λ×K，执行步骤(15c)，否则，将时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的旋角参量θ、水平参量Q以及垂直参量P初始值赋为0，执行步骤(15d)，其中边缘约束系数λ的取值范围为2≤λ≤4；

(15c)若列块序号j<u-1时，则行块序号i保持不变，列块序号j的值加1，返回步骤(15b)；若列块序号j＝u-1且行块序号i<u-1时，则行块序号i的值加1，列块序号j的值赋为2，返回步骤(15b)；若列块序号j＝u-1且行块序号i＝u-1时，则将所有行块序号i＝1或i＝u-1，或列块序号j＝1或j＝u-1的粗配准边缘子图像CR_ij的光谱反射率曲线总差异量DRC_ij置为空值，执行步骤(16)；

(15d)将时相2的各波段相对地物光谱反射率图像R₂ ^b中与粗配准边缘子图像CR_ij中心点空间位置对应的像素点记为Ro，其横坐标为(i-1)×K+v，纵坐标为(j-1)×K+v；以像素点Ro为中心点，分别将时相2的各波段相对地物光谱反射率图像R₂ ^b逆时针旋转θ-5度，采用最近邻线性插值方法进行图像插值，得到各波段的相对地物光谱反射率旋转图像在该图像中将以像素点Ro为中心、大小为(3v+1)×(3v+1)像素的图像块作为时相2的各波段窗图像

(15e)分别将时相1的各波段相对地物光谱反射率图像R₁ ^b中与粗配准边缘子图像CR_ij相同空间位置范围的图像块作为时相1的各波段搜索图像

(15f)将时相1的各波段搜索图像中的第1行第1列像素点与时相2的各波段窗图像中的第1+Q行第1+P列像素点空间位置对应，则时相1的各波段搜索图像中的像素点(x’,y’)与时相2的各波段窗图像中的像素点(x’+Q,y’+P)空间位置对应，计算时相1的各波段搜索图像与时相2的各波段窗图像之间的光谱反射率曲线差异量DRC(θ,Q,P)：

$\mathrm{DRC}(\mathrm{\&theta;},Q,P)=\underset{x^{\&prime;}=1}{\overset{K}{\mathrm{\&Sigma;}}}\underset{y^{\&prime;}=1}{\overset{K}{\mathrm{\&Sigma;}}}{V}^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})\&CenterDot;M^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})---3)$

其中，V^θQP(x′,y′)为时相1的各波段搜索图像与时相2的各波段窗图像之间对应像素的反射率变化矢量的方差：

$V^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})=\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{\left\{\right[{\stackrel{\&prime;}{R}}_{1\mathrm{ij}}^b(x^{\&prime;},y^{\&prime;})-{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P)]-\frac{1}{B}\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}[{\stackrel{\&prime;}{R}}_{1\mathrm{ij}}^b(x^{\&prime;},y^{\&prime;})-{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P)\left]\right\}}^2---4)$

M^θQP(x′,y′)为时相1的各波段搜索图像与时相2的各波段窗图像之间对应像素的反射率变化矢量的模值：

$M^{\mathrm{\&theta;QP}}(x^{\&prime;},y^{\&prime;})=\sqrt{\underset{b=1}{\overset{B}{\mathrm{\&Sigma;}}}{[{\stackrel{\&prime;}{R}}_{1\mathrm{ij}}^b(x^{\&prime;},y^{\text{\&prime;}})-{\stackrel{\&prime;}{R}}_{2\mathrm{\&theta;ij}}^b(x^{\&prime;}+Q,y^{\&prime;}+P)]}^2}---5)$

其中，x’和y’分别为时相1的各波段搜索图像中的像素的行序号和列序号，x’＝1,2,…,K，y’＝1,2,…,K；

(15g)若水平参量Q<3v+1-K，则垂直参量P保持不变，水平参量Q的值加1，返回步骤(15f)；若水平参量Q＝3v+1-K且垂直参量P<3v+1-K，则将水平参量Q的值置为0，垂直参量P的值加1，返回步骤(15f)；若水平参量Q＝3v+1-K且垂直参量P＝3v+1-K且旋角参量θ<10，则将旋角参量θ的值加1，水平参量Q、垂直参量P的值都置为0，返回步骤(15d)；若水平参量Q＝3v+1-K、垂直参量P＝3v+1-K且旋角参量θ＝10，执行步骤(15h)；

(15h)将计算得到的光谱反射率曲线差异量DRC(θ,Q,P)中的最小值作为时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像集{R_1ij}相匹配的光谱反射率曲线总差异量DRC_ij，将该最小值对应的旋角参量、水平参量和垂直参量分别作为时相2的相对地物光谱反射率图像集{R₂}与时相1的相对地物光谱反射率子图像{R_1ij}相匹配的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij，返回步骤(15c)循环计算不同行块序号或不同列块序号的时相1的相对地物光谱反射率子图像R_1ij与时相2的相对地物光谱反射率图像R₂的匹配旋角参量θm_ij、匹配水平参量Qm_ij以及匹配垂直参量Pm_ij。