非下采样剪切波域抗错性尺度不变的遥感图像匹配算法

刘志强，朱立谷

(中国传媒大学 计算机学院，北京 100024)

1 引 言

图像匹配是指将两幅或多幅图像进行比较[1]，找到它们的共有景物。遥感图像匹配在军事和民用方面都具有重要的价值，但由于高分辨率遥感图像信息量巨大，一般都有几何畸变、地貌特征变化、地面杂波混入等影响，所以遥感图像匹配更加困难。

随着计算机技术的不断发展，图像匹配速度有很大程度上的提高，自动的归纳图像匹配技术弥补了人工方法的不足，已经被广泛地应用于遥感图像匹配处理领域。

过去几十年中，各种图像匹配算法相继出现[2-4]，且人们在结合许多数学理论和方法后不断提出新的匹配算法，一般可分为基于灰度的匹配和基于特征的匹配。基于灰度的匹配方法直接根据图像的灰度定义参考图像和待匹配图像之间的相似性度量。这些基于灰度的匹配算法计算简单、易于实现，但是对图像畸变的适应能力较弱，对遥感图像匹配不能产生很好的作用。基于特征的匹配方法是先从待匹配的图像和参考图像中分别提取特征，然后在两幅图像的特征之间建立对应关系。匹配中常用的特征主要有点、边缘、区域特征3种。点特征易于标示和操作，同时也反映了图像的本质特征，所选取的兴趣点是指相对于邻域表现出某种奇异性的像素点。点特征容易提取，但所含信息量较少，所以建立两幅图像中特征点对应关系和保证适当的特征点数目是其难点所在。吴一全等[2]提出了一种利用Contourlet变换、Tsallis熵和改进粒子群优化的多源遥感图像匹配算法。叶沅鑫等[3]针对多源遥感图像间几何变形和灰度差异造成的匹配难题，提出了一种结合SIFT和边缘信息的匹配方法，相比直接利用SIFT自带的特征描述向量，该算法有效地提高了匹配正确率。陈飒等[4]提出了一种基于Contourlet域Hausdorff距离和粒子群的多源遥感图像匹配算法，实现了全分辨率情况下遥感图像的匹配。为了进一步提高遥感图像匹配的精度和运算效率。Zhang等[5]提出了一种基于特征点匹配的遥感图像配准方法，利用特征点的三角形局部特征区域表示遥感图像的局部特征，再利用KNN算法进行相似区域匹配，遥感图像的特征表示鲁棒性得到了提升，匹配精度也随之得到改善。Shen等[6]利用以SIFT和归一化互相关技术为基础提出了一种新颖的遥感图像匹配技术，提高了遥感图像匹配精度基础上，降低了时间复杂度。

基于此分析，本文提出一种非下采样剪切波域抗错性尺度不变的遥感图像匹配算法。首先，利用非下采样剪切波变换分解待匹配图像获得低频子带，采用抗错性尺度不变检测器提取分布均匀且稳定的特征点；然后，利用四元数指数矩提取出抗错性尺度不变特征点局部特征区域的特征，并通过相似度计算进行预匹配。最后，通过随机抽样一致去除掉存在的错误匹配对。通过仿真实验，证明了提出算法的有效性和实用性，显著地提高了遥感图像匹配的准确率。

2 提出的算法

2.1 特征点与特征区域

近些年，很多特征点检测器被相继提出，并广泛应用于目标识别及图像处理领域，例如SIFT(Scale-Invariant Feature Transform)特征点[7]、SURF(Speeded Up Robust Feature)特征点[8]等。上述特征点在目标识别和图像配准领域得到了广泛应用，但其鲁棒性仍然偏低，尤其在图像匹配中会严重降低精度。为解决此问题，本文将使用NSST域的SIFER特征点[9]作为遥感图像感兴趣点，提取出稳定的特征点，具体步骤如下：

1)对遥感图像进行二级NSST分解，并提取低频子带；

2)在低频子带上创建尺度空间。在x，y方向上，对图像进行n尺度滤波，先需要计算n+2个波长位置，公式如下所示：

(1)

若所求的连续波长之间差值小于2，则λj=λj-1+2，这样做是为了避免在较小的波长位置为亚像素级采样特征尺度而失去意义。计算与λ0，λ1，...，λn+1相应的滤波尺度σ0，σ1，...，σn+1。

(2)

其中，λ是基本正弦波的波长，σ是高斯窗口宽度，b表示半幅带宽。

3)尺度空间滤波。调制高斯函数与余弦函数构造一维余弦调制的高斯滤波(Cosine Modulated Gaussian，CMG)表示为：

(3)

其中，φ表示滤波方向。这里将一维CMG拓展到二维CMG滤波，可以对图像分别在x，y方向上进行滤波处理，得到的二维CMG滤波器：

CMGX(x，y；σ，λ)=e-(x2+y2)/2σ2cos(2πx/λ)

(4)

CMGY(x，y；σ，λ)=e-(x2+y2)/2σ2cos(2πy/λ)

(5)

4)尺度空间响应计算。利用公式(4)与公式(5)求出尺度σ0，σ1，...，σn+1上分别在x，y方向上的CMG响应，并分别表示为CMGX及CMGY，计算出在不同尺度上CMGX和CMGY的加和响应值Sj：

Sj=CMGXj+CMGYj

(6)

5)尺度空间极值检测。在CMG加和响应的尺度空间内利用3×3×3的窗口提取极大值和极小值点：

Sj(×)>∀Sj(·)，Sj(×)>∀Sj-1(·)
Sj(×)>∀Sj+1(·) ，Sj(×)<∀Sj(·)
Sj(×)<∀Sj-1(·)，Sj(×)<∀Sj+1(·)

(7)

其中，“×”表示当前位置的响应之和，“·”表示相邻位置的响应之和，通过公式(7)将判断为极值点的位置作为候选的特征点。

6)去除边缘响应。图像边缘区域通常会出现响应和的极值点，但这些位置并不是很好的特征点位置.步骤5就是用来去除边缘响应的影响，实现去除边缘区域的目的：

(8)

其中，x=(x，y)，gx=CMGX，gy=CMGY，W=λ，计算公式(8)里矩阵R的特征值的比：

eratio=γmax/γmin

(9)

(10)

(11)

通过eratio值来判断出CMG加和响应中的极值点是否为出现在边缘区域，如果候选特征点的eratio大于某个阈值时(本文算法中边缘阈值选择为6)，该点即被判定为边缘点，则该点需要被剔除掉，而剩余的响应极值点就是最终剩余的SIFER特征点。

2.2 特征点提取与局部特征区域构造

通常局部区域可以构造成很多形状，例如长方形、正方形、圆形和椭圆形等，但都需要能够使构造的区域对局部去同步攻击有很好的抵抗能力。本文自适应地构造出椭圆特征区域，相比较其他算法具有更好的鲁棒性。椭圆区域的构造可以分成以下环节：

1)利用SIFER检测器在宿主遥感图像中提取出特征点，X(x，y)表示椭圆区域中心点。

2)通过公式(8)可以计算出二阶自相关矩阵R(x)，然后根据公式(12)和(13)自相关矩阵R(x)的特征值(λ1，λ2)：

(12)

(13)

3)分别利用特征值(λ1，λ2)和特征向量(v1，v2)计算出局部椭圆区域的长轴ra、短轴rb和方向角θ：

(14)

2.3 四元数指数矩及稳定性能分析

本文中期望借助指数矩的优秀描述能力对遥感图像局部特征区域进行描述，所以这里将结合四元数理论对指数矩阵进行改进，从而推到出适合于彩色图像的四元数指数矩。

四元数作为复数的一种推广[10]，由爱尔兰数学家哈密尔顿正式提出。四元数表现为一种超复数形式，由一个实部和三个虚部组成，可以写成如下表达式：

q=a+bi+cj+dk

(15)

其中，a、b、c和d都表示为实数，i、j和k表示三个虚数单位，并需要满足以下规则：

i2=j2=k2=ijk=-1
ij=k，jk=i，ki=j
ji=-k，kj=-i，ik=-j

(16)

由上面的公式，可以把四元数的共轭和幅值的定义成如下形式：

(17)

假设一个大小为M×N的彩色图像f(x，y)具有红色、绿色、蓝色三个颜色通道，则可以用四元数形式表示为纯四元数(没有实数部分)，具体表示为：

f(x，y)=fR(x，y)i+fG(x，y)j+fB(x，y)k

(18)

其中，fR(x，y)、fG(x，y)和fB(x，y)分别代表了彩色图像f(x，y)的红色、绿色及蓝色通道。这里面的i，j和k为虚数单位。

(19)

根据正交完整函数系理论，彩色图像f(r，θ)可以通过有限阶数的指数矩(n≤nmax，l≤lmax)近似重构回来，如果重构中矩的阶数越多，重构回的图像就更加接近原图像，具体可以表示为如下形式：

(20)

结合指数矩和四元数理论，可以推导出四元数指数矩，同理，也可以得到四元数Zernike矩、四元数伪Zernike矩、四元数OFMMs，为了充分对比每种四元数矩的重构能力，进行了下面的重构实验。

本文中利用四元数指数矩作为SIFER特征点的局部特征，除了良好的特征描述能力，四元数指数矩对常规信号攻击和几何攻击均具有良好的抵抗能力。这部分将对四元数指数矩的几何抵抗原理进行分析，具体如下所示。

2.3.1 平移不变性

首先计算一幅N*N图像f(i，j)的质心(xc，yc)，并将图像中心移到质心坐标位置，再计算其四元数指数矩，就得到了指数矩的平移不变性。

xc=m1，0/m0，0，yc=m0，1/m0，0

2.3.2 缩放不变性

设图像函数为g(r′，θ)，找到图像的半径k，则r′的变化范围为0≤r′≤k，其归一化灰度图像函数可表示为：f(r，θ)=g(kr，θ)=g(r′，θ)。其中，r=r′/k的取值范围是0≤r≤1，f(r，θ)是经过归一化的极坐标下图像表示。此时，在极坐标系下通过f(r，θ)获得的矩是缩放不变的，对于相同的f(r，θ)，在0≤r≤1范围内，经过缩小或放大，对于任意f(r′/k，θ)都可以最终归一化为同一个图像函数f(r，θ)，所以归一化后的彩色图像四元数指数矩就具有了缩放不变性。

2.3.3 旋转不变性

图1给出了四元数指数矩(幅值)抵抗几何攻击能力的测试结果(128×128×8bit灰度图像Barbara)。实验结果表明，四元数指数矩的幅值具有很好的几何不变特性。

2.4 提出的算法步骤

本文的核心步骤包括遥感图像NSST低频的SIFER特征点提取，局部区域四元数指数矩特征提取环节以及图像配准环节。利用SIFER算子从NSST低频提取出稳定性更高的遥感图像特征点。考虑颜色通道间的相关性，充分利用了彩色图像中重要的颜色信息，利用四元数指数矩幅值表示每个特征点的局部区域的特征，并将其构成特征向量进行初步相似度匹配。对于匹配出现的错误匹配对，通过随机抽样一致进行去除，有效地提高了遥感图像匹配的精度。该算法引入了稳定性及特征表示能力的SIFER算子提取特征点，针对彩色的遥感图像，利用颜色不变量理论充分考虑了重要的颜色内容，使提出的SIFER特征点更加稳定，利用指数矩作为每个特征点局部区域的特征提高了区域特的刻画能力。该方案的基本流程图如图2所示。

(a)无攻击 (b)放大1.2

(c)缩小0.5 (d)旋转15度

(e)旋转30度并缩小0.7倍 (f)放大1.2、旋转15度图1 四元数指数矩(幅值)抵抗几何攻击能力的测试结果

3 实验结果分析

分别利用SIFT遥感图像匹配算法[11]、SURF遥感图像匹配算法[12]和本文提出的算法，针对实际的遥感图像进行实验对比分析，所有算法均采用相同的特征表示方法以及RANSAC算法。实验环境为8.0GB的RAM，3.80GHz的CPU，实验平台为Matlab 2015a。

评价标准使用检测率用来衡量匹配性能的好坏，通过检测率柱状图来客观形象的展示检测率的高低。检测率(R)定义为正确匹配点对数和待匹配图的点个数的比值：

(21)

实验的参考图像为谷歌地球下载的4幅可见光卫星遥感图像，涵盖北京地区4个不同类型的场景，分别如图3所示，大小均为512×512。评价性能的指标包括匹配算法的抗击何攻击能力和抗常规攻击能力。抗击何攻击能力中主要验证算法的旋转不变性和尺度不变性，在抗常规攻击能力中主要验证在椒盐噪声、高斯噪声以及亮度变化下算法的匹配性能。

3.1 抗几何攻击能力测试

以图3中居民区2图像为例，将待匹配图像分别旋转20度，45度，70度，如图4所示。将图3中的参考图像分别与攻击后遥感图像进行匹配，分析3种算法的抗旋转能力，详细对比结果见表1。

图2 遥感图像匹配流程图

(a)居民区1 (b)体育场 (c)居民区2 (d)河流图3 遥感图像参考图像

(a)旋转20度 (b)旋转45度 (c)旋转70度图4 旋转后图像

旋转角度算 法待匹配图像特征点个数预匹配点对数RANSAC后匹配点对数匹配率20度SIFT算子2751258430.18%SURF算子29714410936.70%本文算子31016913242.58%45度SIFT算子2891466422.15%SURF算子3121759129.17%本文算子33819110531.07%70度SIFT算子2681329133.96%SURF算子28915411439.45%本文算子30318314146.53%

以图3中居民区2图像为例，将待匹配图像分别缩放0.5，0.7，1.5，如图5所示。将参考图像分别与攻击后遥感图像进行匹配，分析3种算法的抗缩放能力，详细对比结果见表2。

由表1和表2综合看出，本文算法抗旋转能力最优，明显高于SIFT算法和SURF算法。

3.2 抗常规攻击能力测试

对三种算法的抗常规攻击能力进行测试，如图6中分别给出了添加椒盐噪声0.02，0.05，0.1以及高斯噪声标准方差0.02，0.05，0.1的6幅测试图像。

将参考图像分别与图7中各幅图像进行匹配，分析3种算法的抗噪声能力，结果如图7所示。图7(a)是3种算法对含高斯噪声遥感图像的匹配性能对比，横坐标是噪声方差，纵坐标是匹配率；图7(b)是3种算法对含椒盐噪声遥感图像的匹配性能对比，横坐标是噪声密度，纵坐标是匹配率。

(a)缩放0.5 (b)缩放0.7 (c)缩放1.5图5 缩放后图像

缩放尺度算法待匹配图像特征点个数预匹配点对数RANSAC后匹配点对数匹配率0.5SIFT算子2751258430.18%SURF算子29714410936.70%本文算子31016913242.58%0.7SIFT算子2891466422.15%SURF算子3121759129.17%本文算子33819110531.07%1.5SIFT算子2681329133.96%SURF算子28915411439.45%本文算子30318314146.53%

(a)椒盐噪声0.02 (b)椒盐噪声0.05 (c)椒盐噪声0.1

(d)高斯噪声0.02 (e)高斯噪声0.05 (f)高斯噪声0.1图6 添加椒盐噪声和高斯噪声后图像

(a) (b)图7 不同噪声下匹配结果对比

由图7可知，虽然在不含噪声的情况下，SURF算法的匹配率低于SIFT算法，但随着高斯噪声方差/椒盐噪声密度的增加，SIFT算法的匹配率在大幅下降，幅度远高于SURF算法，甚至最后匹配率会低于SURF算法，抗噪能力明显不如SURF算法；而本文算法不仅在无噪声的情况下匹配率最高，优于SIFT算法和SURF算法，而且抗噪性能也最好，匹配率始终保持对其他两种算法的绝对优势。这是因为NSST把图像分解成1个低频分量和多个不同方向的高频分量，低频分量基本保持了原图的整体特征，而高频反映了图像的细节，含有大部分的噪声，如果匹配的时候把高频也考虑进去，不但常常受高频噪声的影响而产生误匹配，而且匹配分析的时间消耗也很大。而本文只对图像低频分量进行匹配，可以大大降低噪声等细节的影响，使得算法既提高了匹配的速度也提高了匹配的精度。

对待匹配图像处理，构造中、低、高3种不同亮度，如图8所示。其中，图8(b)由(a)亮度降低4倍所得，(c)为(a)亮度提高2倍所得。将参考图像分别与图8中各幅图像进行匹配，分析3种算法的抗亮度变化能力，详细对比结果见表3。

由表3可知，在不同的亮度情况下，3种算法提取到的特征点数目基本保持不变，只是随着亮度增加而略有增加；当亮度较低时，只有SIFT算法的匹配率略有下降，SURF算法和本文算法基本不变；当亮度较高时，3种算法的匹配率都有所降低，但SIFT算法下降幅度最大，SURF算法次之。由此得出，本文算法抗亮度变化性能最好，SURF算法次之，SIFT算法最差。

图8 不同亮度下遥感图像

亮度变化算 法待匹配图像特征点个数预匹配点对数RANSAC后匹配点对数匹配率亮度分量×0.5SIFT算子2751258430.18%SURF算子29714410936.70%本文算子31016913242.58%亮度分量×1.2SIFT算子2891466422.15%SURF算子3121759129.17%本文算子33819110531.07%亮度分量×1.5SIFT算子2681329133.96%SURF算子28915411439.45%本文算子30318314146.53%

4 总结

本文提出一种非下采样剪切波域抗错性尺度不变的遥感图像匹配算法。首先，利用非下采样剪切波变换分解待匹配图像获得低频子带，采用抗错性尺度不变检测器提取分布均匀且稳定的特征点；然后，利用四元数指数矩提取出抗错性尺度不变特征点局部特征区域的特征，并通过相似度计算进行预匹配。最后，通过随机抽样一致去除掉存在的错误匹配对。