-
遥感影像配准是遥感影像变化检测技术中至关重要的数据预处理环节，遥感载荷技术指标差异、观测角度、时相、地形起伏等内外部因素造成的几何畸变，给遥感影像变化检测带来极大挑战。
Dai等[1]研究影像配准误差对变化检测精度的影响，实验结果表明影像配准误差小于0.2个像元，变化检测错误率小于10%；史文中等[2]采用缓冲区分析的方法研究了影像配准误差对变化检测误差分布的模式，认为配准误差在0~1个像元, 影像变化检测误差像元基本上位于影像中目标边缘位置；邬建伟等[3]对遥感影像多项式配准精度影响因素进行分析，平坦地区TM影像的多项式配准精度可达到0.2个像元，但对于地形起伏区域和多源遥感影像配准很难达到0.2个像元的精度。针对影像配准中的局部几何畸变问题，Goshtasby等[4]、Ehlers[5]、Wiemker等[6]利用具有一定空间分布、数量足够的控制点，在局部变形区域通过局部距离加权、二次曲面函数、薄板样条函数等进行局部配准；不少学者采用基于不规则三角网的多源遥感影像配准方法，通过自动提取多源影像同名点并构建三角网(triangulated irregular network，TIN)，对每个三角形建立配准几何变换模型，较全局遥感影像配准方法大幅提高了影像的整体和局部配准精度，但该类算法对多源影像同名点提取的精度以及空间分布有较高的要求[7-10]。在基于像元的局部配准方面，张泽旭等[11]用全局光流场配准红外图像的背景区域，再以层次匹配算法精确配准目标，该方法对场景的运动有很好的鲁棒性；李峰等[12-13]将弹性配准算法应用于遥感影像处理领域，成功用于影像超分辨率重建以及航空大视场多光谱线阵扫描仪影像多光谱波段间数据的高精度匹配。
本文在文献[12-13]的基础上提出一种基于加速抗差特征(speed up robust feature，SURF)算法[14]的全局匹配和局部配准模型相结合的面向影像变化检测的弹性配准方法，在全局初步配准基础上，以像元为单位建立局部配准参数解算和平滑模型，考虑到遥感影像变化检测中不同时相影像之间存在的局部辐射亮度差异, 本文采用不同时相遥感影像差值特征影像各像元正态分布密度函数构建像元局部参数解算权重进行局部参数解算，采用影像金字塔由粗到细的数据处理策略，完成不同时相影像弹性配准。
HTML
-
面向影像变化的弹性配准方法主要包括影像初步全局配准、局部平移参数解算、局部平移参数平滑3个主要步骤。
-
文献[12-13]中，采用基于像元灰度的最小二乘仿射变换法进行全局自动配准，对影像的重叠度和几何变换的初始值有一定的要求(初始值误差过大或影像重叠度太低均会造成匹配失败)，且计算量较大。本文采用SURF算法实现影像全局仿射变换配准，主要包括兴趣点检测、描述子生成、兴趣点匹配、粗差点剔除、影像全局仿射变换等过程，效率和适用性明显提高。在粗差匹配点剔除方面，先采用随机抽样一致性方法(random sample consensus，RANSAC)[15]初步剔除错误匹配点，采用最小二乘法求解仿射变换参数，计算匹配点均方根误差(root mean squared error，RMSE)，当RMSE小于给定的阈值时，参数为最终全局仿射变换参数；当RMSE大于给定的阈值时，剔除匹配点中误差最大的同名点，重新计算仿射变换参数，如此循环运算直至RMSE小于给定阈值，确定最终全局仿射变换参数[16]。
-
1) 局部平移模型
经过基于特征点的影像全局匹配，参考影像与浮动影像之间实现粗略匹配。考虑到遥感影像“刚性”程度较强，配准影像间像元通过水平和垂直方向的小范围平移可实现更高精度配准，基于这一假设以像元为单位在一定窗口范围W(以像元为中心，长宽分别为n个像元的区域像元集合，一般取5≤n≤50)内建立局部平移模型：
式中，WT(x, y)为参考影像窗口内坐标为(x, y)处像元的灰度值(范围为0~255)；WS(u, v)为浮动影像相同窗口大小内对应像元(u, v)的灰度值；u=x+m1，v=y+m2，m1和m2为参考影像窗口中心像元像元匹配至浮动影像像元的水平和垂直方向的位移，构成平移参数向量m =[m1, m2]T，m为2维向量，式(1)可变换为：
构建窗口像元局部平移代价函数为：
利用式(3)最小化为准则，求m。
式(3)一阶泰勒级数展开近似表达为：
式中, WTu和WTv为参考影像窗口内像元WT(x, y)在水平和垂直方向的一阶偏导数; WT、WS为WT(x, y)、WS(u, v)简写。
式(4)可简写为：
式中, k为窗口集合W内对应像元的灰度差，c为2维梯度向量。
式(5)对m求导数得：
由式(8)得：
式(9)实际上是由窗口W内n2个像元通过最小二乘法估计中心像元平移参数m的过程。其中C为由参考影像窗口W内n2个像元的梯度向量组成的2×n2维矩阵，k为n2个对应像元的灰度差构成的n2×1维向量。m向量为参考影像像元配准至浮动影像像元的平移参数，而浮动影像像元配准至参考影像的平移参数向量为Δ m =[-m1, -m2]T，采用这种逆向匹配策略[12-13]可使矩阵C在迭代计算过程中保持不变，提高参数解算效率。
2) 局部加权平移模型
由式(9)可知，平移参数与窗口内各像元的灰度差值相关，不同时相遥感影像局部存在辐射亮度差异时必然会对平移参数的准确性产生影响。通常参考影像与配准浮动影像的像元灰度差近似服从正态分布，如式(10)、式(11)所示：
Z为像元灰度差，其中μ和σ分别为灰度差分布的均值和标准差，f(Z)为灰度差正态分布密度函数，一般情况下参考影像和浮动影像中的辐射亮度差异较大的像元较少，即变化像元出现的概率较低。根据这一假设，可以各像元灰度差构建各像元参与局部匹配运算的权w(Z)：
对式(9)局部窗口参与解算各像元赋予不同的权重：
式中, Wz为n2×n2对角矩阵，对角线元素为相应位置像元的权值记为wz(x, y)。
-
浮动影像各像元经过局部平移参数解算后，每个像元会求得一组平移参数，通过平滑约束处理可使各像元平移匹配参数与相邻像元参数之间具备较好的连续性和一致性。根据这一准则建立平滑模型为：
式中，Epix(m)为每个像元整体匹配的代价函数; Exy(m)为平移代价函数; Esm(m)为平滑约束代价函数; 式(17)中ζe为约束权重(0≤ζe≤1，其值越大平滑性越强)；以Epix(m)取最小值为准则进行平滑约束平移匹配参数解算，详细推导过程参见文献[12-13]。
-
本文弹性配准算法，首先，采用SURF算法提取同名特征点实现影像全局匹配；其次，建立影像金字塔从分辨率最低图层开始进行局部弹性配准。对当前图层参考影像和浮动影像计算差值特征图像，求出各像元的权值。进入局部加权参数计算环节，在一定大小的窗口内计算每个像元的平移参数，对全部像元的平移参数进行平滑，判断平移参数解算是否达到精度要求，未达到则利用当前参数进行影像变换和重采样，进行新一轮局部平移参数解算和平滑处理，直至满足精度要求。对当前各像元的解算参数采用双线性插值，获得下一分辨率图层匹配的初始值。重复上述过程，最后利用各像元解算参数对浮动影像进行最终配准，算法流程如图 1所示。
Figure 1. Flowchart of Elastic Registration Data Preprocessing
1.1. 基于SURF算法的全局初步配准
1.2. 局部平移参数解算
1.3. 局部平移参数平滑
1.4. 弹性配准变化检测数据预处理流程
-
本文实验数据采用中国镇江市2003年QuickBird多光谱数据作为参考影像，2004年QuickBird多光谱数据、2005年SPOT-5全色与多光谱融合数据、2006年ALOS全色与多光谱融合数据分别作为浮动影像，影像空间分辨率均采样为2.5 m。实验区域影像中含有山地、建筑物、道路等，有较多的变化信息，数据来自三种不同的遥感载荷，对本文算法验证具有一定的代表性。实验数据分为三组，第一组参考影像为2003年QuickBird影像，浮动影像为2004年QuickBird影像，影像大小为948×734像元，该组数据由于属于同一传感器且影像差异较小，采用10×10窗口大小进行弹性配准运算；第二组数据参考影像为2003年QuickBird影像，浮动影像为2005年SPOT-5影像，影像大小为1 002×1 002像元；第三组数据参考影像为2003年QuickBird影像，浮动影像为2006年ALOS影像，影像大小为1 002×1 002像元；第二、三组影像由于差异较大采用20×20窗口进行弹性配准运算。实验分别对三组数据进行全局仿射变换配准、文献[12-13]方法、本文方法弹性配准，采用在典型位置(地形起伏误差较大的区域)人工选取15对同名点计算x、y方向的均方根误差RMSE(式(18))的方法进行配准精度评价，并将浮动影像配准结果分别与参考影像、文献[12-13]与本文方法配准结果之间进行求差取绝对值运算，再进行密度分割，目视分析两种方法的配准效果。此外对于空间分辨率相差较大的影像，将第二组数据中SPOT-5数据采样成10 m空间分辨率数据进行弹性配准实验，通过密度分割效果图可直观反映弹性配准方法对于尺度差异较大影像的适用性。
式中, zi为配准影像在同名点i处(x, y)坐标; zi*为参考影像在在同名点i处(x, y)坐标。
-
在全局SURF仿射变换全局配准环节特征点提取阈值为0.000 1；粗差剔除环节，RANSAC粗差点剔除阈值为0.01，RMSE阈值为3像元；弹性配准环节，平移参数迭代收敛条件为：方向平移参数前后两次迭代平均误差小于0.02个像元为局部运算结束条件；参数平滑收敛条件为：方向平移参数平滑前后两次迭代平均误差小于0.001为平滑迭代结束条件。本文实验由于文献方法在全局初步配准环节几何形变较大时无法实现全局仿射变换配准，故全局配准环节均采用SURF全局配准方法，算法时间从局部弹性配准开始计时，直至算法结束。实验的软件和硬件条件为：MATLAB 2010b / 32-bit Windows xp / Intel Core(TM) i3-M 380/ CPU2.53GHz / 2 GB RAM。
-
如表 1所示，三组数据实验数据表明文献[12-13]和本文弹性配准方法较全局配准方法配准精度有较大的提升，影像之间配准误差可由5~10个像元降低至1个像元以内(包括地形起伏区域)；本文加权弹性配准方法，通过在局部参数解算环节，采用加权方法，降低影像辐射亮度差异较大像元对局部参数解算的影响，局部差异较大区域配准精度和配准效果显著提高(如图 2(c)和图 2(d)可看出，在图 2(c)中新增道路有显著扭曲现象，而图 2(d)中道路扭曲程度不大)，且较文献[12-13]方法速度有一定的提高。本文以与参考影像差异最大的第二组数据为例进行讨论。
Figure 2. Experiment Result Images of the Second Set of Images
图 3(a)、图 3(b)、图 3(c)分别为第二组2.5 m分辨率影像全局配准影像与参考影像、本文方法配准影像与参考影像、本文方法与文献[12-13]方法配准影像的差值密度分割特征图，图 3(d)、图 3(e)、图 3(f)分别为第二组采样为10 m分辨率浮动影像全局配准影像与参考影像(2.5 m分辨率)、本文方法配准影像与参考影像、与文献[12-13]方法配准影像的差值密度分割特征图(差值密度分割特征图为影像之间差值绝对值图像进行32个灰度等级的图像分割，由蓝到深红表示灰度值由低至高，选用上述配准影像和参考影像红波段数据构建图 3(a)至图 3(f)特征图)。由图 3(a)和图 3(d)可以看出，全局配准影像特征图由于存在较大配准误差，高亮噪声显著；图 3(b)和图 3(e)图显示出本文配准方法可有效降低配准误差产生的噪声；图 3(c)和图 3(f)显示出本文方法和文献[12-13]方法除在辐射亮度显著差异区域存在高亮噪声，其它区域配准效果保持着较好的一致性，图 3(d)至图 3(f)也说明本文方法对于空间分辨率相差4倍的影像仍然适用。
Figure 3. Difference Indesity Sliced Feature Images of Registration Images
本文弹性配准方法中两个重要的参数即局部窗口大小和影像金字塔的层数需针对不同的情况进行设置。当影像局部存在较大的配准误差或影像局部辐射亮度差异较大，应选取大窗口，一般为20以上，一方面同一窗口尽可能多地容纳影像局部同名点，另一方面抑制少数噪声像元的影响，配准影像“刚性”较强；大窗口具有上述优势的同时，也存在着计算量大，无法实现小目标的精确配准的问题，因此窗口大小的选择仍是一个难题。在影像金字塔层数方面，Bruzzone等[17]和Celik等[18]曾采用多尺度变化检测的方法降低配准误差产生的伪变化，本文局部加权过程本质上属于检测变化的过程，当像元之间存在较大的配准误差时，此时加权方法会将配准误差较大像元判别为变化像元而导致该像元无法精确配准，降低影像分辨率可减少这种情况的发生；本文参考影像与其它时相的影像之间配准误差大概在10个像元左右，因此配准时建立3层影像金字塔，空间分辨率分别为2.5、5、10 m。本文弹性配准方法局部参数平滑环节一方面使配准影像看起来更加平滑连续，另一方面也存在像元局部真实位移与邻近像元发生突变时，由于平滑的原因无法获得像元真实平移参数的问题。
DownLoad:
