基于相位一致性计算和RS-GLOH描述子的SAR和光学图像匹配*

嵇宏伟,刘畅,潘志刚,申芳瑜

(中国科学院空天信息创新研究院,北京 100190; 中国科学院大学,北京 101408) (2021年11月11日收稿; 2021年12月30日收稿修改)

异源遥感图像匹配是将不同传感器获取的两幅或多幅图像中的相似特征或区域进行空间位置上的匹配对应过程[1],其中应用最广的是合成孔径雷达(synthetic aperture radar, SAR)和光学图像匹配。光学成像符合人眼视觉,图像信息易于解读,SAR为主动式微波遥感成像传感器,能够进行全天时、全天候的对地观测,可以弥补光学图像易受天气影响的缺点,因此SAR和光学图像匹配结果能够有效发挥各自的优势[2]。但是由于成像机理不同,SAR与光学图像之间存在较大的非线性辐射和几何差异,使得两者间的匹配成为一个难点[3]。所以,实现SAR和光学图像匹配在飞行器匹配导航、SAR与光学的信息融合等应用中具有重大意义。

SAR和光学图像匹配算法主要分为基于区域和基于特征的两种方法[4]。基于区域的方法是从参考图像中提取一定大小的模板子图像并使用相似性度量函数在另一幅图像中搜寻最优匹配区域[5],相似性度量函数主要包括互信息[6]和归一化互相关[7]等。在实际应用中,基于区域的方法计算量大,不适用于两幅图像之间存在旋转和尺度差异的情况,有时还需要提供待匹配区域的地理参考信息[8]。基于特征的方法主要对提取出的点[9]、线[3]和轮廓[10]等显著特征信息建立特征描述向量后进行匹配,与基于区域的方法相比,基于特征的方法计算量小,鲁棒性强,并且具有尺度和旋转不变性,因此本文使用基于特征的方法。

基于特征的方法分为基于空间域和基于频域的两种方法[11]。基于空间域的方法中应用最广的是尺度不变特征变换(scale-invariant feature transform,SIFT)[12]。但是将SIFT运用于SAR和光学图像匹配中很难得到良好的匹配结果,所以一些学者提出了针对SAR和光学图像匹配的改进SIFT算法。例如Ma等[13]提出将每个特征点的位置、尺度和方向联合起来的增强型特征匹配算法,增加了特征点匹配数量。Dellinger等[14]和Xiang等[15]使用指数加权均值比率算子(ratio of exponentially weighted average, ROEWA)建立多尺度Harris空间提取SAR图像的角点,并且后者提出一种特征点位置精调的方法实现SAR和光学图像的精确匹配。但是SAR和光学图像之间的非线性辐射差异导致两者在同一区域表现出不同的灰度特征[16],使得基于空间域的方法很难提取出数量充足且位置对应性好的初始特征点,所以抗辐射差异性更好的频域方法更多地用于多源遥感图像匹配中。基于频域方法中常用的为相位一致性(phase congruency, PC)[17-18],例如Fan等[19]提出相位一致性结构描述子对不同尺度的角点特征进行描述并匹配。Ye等[8]提出相位一致性方向直方图对Harris角点特征进行描述。Li等[20]和Xie等[21]都使用相位一致性模型计算得到的最大矩图作为特征点提取的基础图像,并且前者提出了最大索引图(maximum index map, MIM)对特征点进行描述。但是SAR成像会使得SAR图像中存在透视缩短、叠掩和阴影等几何失真现象,导致两者之间具有非线性几何差异,使得基于相位一致性计算构建的描述子很难得出两者在同一区域的结构相似性信息,所以本文构建了一种新的基于空间域的RS-GLOH特征描述子,可以有效抑制SAR和光学图像间非线性几何差异的对特征描述的影响。

综上所述,本文提出一种结合频域和空间域方法的SAR和光学图像匹配框架,对应流程如图1所示。在特征提取阶段,从基于相位一致性模型计算得到的最小矩图和最大矩图上获得分布均匀且位置对应性好的最小矩点和最大矩点;在特征描述阶段提出RS-GLOH,使用ROEWA和Sobel建立梯度定位和方向直方图(gradient location orientation histogram, GLOH)[22]分别对SAR和光学图像中的特征进行描述,得出两者在同一区域具有尺度和旋转不变性的结构相似性信息;在特征匹配阶段使用最近邻距离比(nearest neighbor distance ratio, NNDR)和快速抽样一致性(fast sample consensus, FSC)分别对最小矩点和最大矩点进行匹配。

1 本文方法

1.1 特征提取

1.1.1 相位一致性理论

Morrone和Owens[23]证明相比于信号的幅度信息,相位信息更有利于特征的提取,并提出初始的相位一致性模型。随后,Kovesi[17]提出改进的相位一致性模型提取图像的结构特征。首先,使用二维Log-Gabor滤波器(2D-LGF)与图像进行卷积运算得出图像的局部频域信息,2D-LGF的频域公式如下

(1)

其中:ω0为滤波器的中心频率;κ/ω0为保持滤波器带宽不变的参数,使得滤波器波形不会产生畸变;θo为滤波器的方向;σθ表示角度带宽。在实际的程序运算中,需要将2D-LGF从频域变换到空间域

ζno(x,y)=ζno-even(x,y)+iζno-odd(x,y),

(2)

其中:ζno-even(x,y)和ζno-odd(x,y)分别为2D-LGF在尺度n和方向o中的偶对称和奇对称滤波器,然后将ζno-even和ζno-odd分别与图像I(x,y)进行卷积运算

(3)

得出I(x,y)在尺度n和方向o的频域幅度和相位分量分别为

(4)

进一步得出相位一致性模型

PC(x,y)=

(5)

ΔΦno(x,y)=

cos(φno-φno-avg)-|cos(φno-φno-avg)|.

(6)

其中:Wo(x)是一维权重函数,ΔΦno(x)为在尺度n和方向o的相位分量;常数ε(ε>0)可以防止分母数值变为0;To是噪声补偿量,满足To=μR+kσR,其中μR和σR为瑞利分布的均值和方差;φno-avg为平均相位角度。

1.1.2 特征点提取

公式(5)所示的相位一致性模型没有考虑方向变化对图像特征的影响,所以Kovesi计算多个方向的PC(θ),并根据经典的矩分析方程,得出具有亮度和对比度不变性的最小矩图m和最大矩图M[18]

(7)

其中3个矩分量为

a=∑θ(PC(θ)cos(θ))2,

b=2∑θ(PC(θ)cos(θ))(PC(θ)sin(θ)),

c=∑θ(PC(θ)sin(θ))2.

(8)

根据上述得出的SAR和光学图像最大矩图和最小矩图如图2所示,可以看出SAR图像相比于光学图像,亮度和对比度较低,具有明显的灰度差异,但是得出的最小矩图和最大矩图能够清晰地表现出两者在同一地区的边缘轮廓信息。因为最小矩图和最大矩图之间具有协同性与互补性[17],所以本文使用极值点检测和非极大值抑制从两种矩图中得到最小矩点和最大矩点,对这两种点特征集合分别进行特征匹配。

本文将最小矩点和最大矩点提取结果与Harris角点进行对比,如图3所示。根据图3(a)所示Harris角点检测结果,在区域A1和A2中没有检测出角点特征,因为区域A1和A2的结构信息较少,并且亮度和对比度较低,很难提取出Harris角点。虽然在结构信息丰富的区域B1和B2中提取出了Harris角点,但是点的位置分布不均匀,而且B1和B2角点的位置对应性并不是很好,不利于后续的匹配。图3(b)所示的本文方法在区域A1和A2、B1和B2中都能够得出分布均匀且位置对应性好的特征点。为区分两种点,本文使用黄点表示最大矩点,红点表示最小矩点,可以看出两种点具有分明的位置区分度。所以,本文将对这两种点分别进行匹配,提高相位一致性

图3 特征点检测对比结果Fig.3 Feature point detection and results comparison

特征点的利用效率。根据上述理论叙述与实验结果,基于相位一致性计算的方法在特征提取阶段能够抑制SAR和光学图像之间辐射差异的影响,得出稳定的初始特征点。

1.2 特征描述

接下来需要为每个特征点建立描述子,考虑到SAR和光学图像之间的非线性几何差异,需要一种特征描述方法能准确描述出SAR和光学图像在同一区域的结构相似性信息。受到文献[11,15]的启发,本文提出一种基于空间域的RS-GLOH描述子,构建流程如图4所示。与原有的单一方法构建特征描述子不同,本文使用ROEWA和Sobel对SAR和光学图像的梯度信息进行提取,然后使用GLOH描述子分别对SAR和光学图像的特征点建立如图4所示的RS-GLOH描述子。

图4 RS-GLOH描述子构建流程Fig.4 RS-GLOH descriptor construction process

ROEWA是一种基于比率的算子,对于相干成像的SAR图像具有很好的梯度检测效果,并且对乘性斑点噪声也具有抑制作用[24]。ROEWA梯度的计算主要依赖于无限对称指数滤波器(infinite symmetric exponential filter,ISEF),其中水平方向的ISEF如图5(a)所示。基于ROEWA的SAR图像梯度求解公式如下

图5 SAR和光学图像的梯度幅值和梯度方向图Fig.5 Gradient amplitude and gradient direction map of SAR and optical images

Gsar-h(x,y)=

(9)

Gsar-v(x,y)=

(10)

其中:K和J分别为ISEF的长和宽,β为适用于SAR 图像的尺度参数。根据式(9)和式(10),得到图像的梯度幅值和方向信息,相应结果如图5(b)和5(c)所示

(11)

Sobel算子是光学图像常用的梯度检测算子[15]。常用的Sobel模板如下所示

(12)

其中:Gh和Gv分别为水平和垂直方向的矩阵模板。考虑到光学图像加性噪声的影响,本文采用高斯加权Sobel模板计算光学图像的梯度,相应的计算过程如下:

(13)

(14)

其中:α为适用于光学图像的尺度参数。图5(d)为水平方向Sobel模板。根据上述结果得出如图5(e)和5(f)所示的光学梯度幅值与方向,对应计算公式为

(15)

在求解对应关键点梯度的过程中需要注意两点:

1)在SAR和光学图像匹配过程中,为了两者之间的梯度幅度一致性,将梯度幅度图进行归一化处理[25]。

2)SAR和光学图像之间的灰度差异会导致在求解梯度的过程中出现梯度反转现象[26],所以本文将原始梯度方向O进行如下处理

(16)

其中Onew为最终使用的梯度方向。

得出图像的梯度信息后,需要对每个SAR和光学图像的特征点分别建立如图4(a)和4(b)所示的ROEWA-GLOH和Sobel-GLOH描述子:首先以某一特征点为圆心,建立一定半径范围的梯度分布直方图并提取特征点对应区域的主方向,让特征点具有旋转不变性。然后基于上述信息构建特征点的GLOH描述子。相比于传统4×4的正方形SIFT特征描述子,对数极坐标描述子GLOH具有更加鲁棒的匹配性能[22]。GLOH的建立标准参照文献[13],假设GLOH的3个圆形区域的半径比率因子为γ,3个圆形区域的半径分别为R={3γ,4.11γ,12γ},根据图5所示的GLOH描述子,本文将特征点对应区域划分为17个独立区域,在每个区域中将梯度方向量化为8个方向,根据梯度直方图统计各个方向的梯度幅值,进而得到对应特征点的17×8=136维特征描述向量。

为说明RS-GLOH描述子构建方法相比于同一种方法构建GLOH描述子的优越性,接下来将RS-GLOH与同一方法得到的GLOH描述子进行对比实验。其中实验需要注意的是,同一方法是指对SAR和光学图像都使用基于Sobel的方法构建GLOH描述子。对应结果如图6所示。

图6 梯度提取和描述子提取对比结果Fig.6 Results comparison of gradient extraction and descriptor extraction

本文主要采用图6(a)和6(d)所示的区域A和区域B作为梯度提取和描述子提取实验的样本。其中区域A为平坦区域,结构信息较少;区域B为建筑物区域,结构信息丰富。根据图6(b)和图6(e)的梯度提取结果,本文方法得到的SAR和光学图像在同一区域的梯度直方图具有更好的重合度,由于SAR图像中存在几何失真现象,使得同一方法得到的SAR图像直方图具有明显的形状畸变。从图6(c)和6(f)得到的描述子结果可以看出,RS-GLOH得到的SAR和光学图像的特征描述子具有更好的相似性。综上所述,本文提出的RS-GLOH能够有效地抑制SAR和光学图像之间的非线性几何差异,从而提取出SAR和光学图像在同一区域的结构相似性信息。

1.3 特征匹配

原有的匹配算法是以特征描述子间的欧式距离为基础来选择匹配点,常用方法为NNDR[12],主要目的是将最近邻和次近邻距离点对应的最近邻距离与次近邻距离的比值与阈值tNNDR进行比较,若比值小于tNNDR,则当前最近邻点为初始匹配点。但是使用NNDR得到的匹配结果仍然存在误匹配,所以本文使用FSC作为外点剔除算法,相比于常用的RANSAC(random sample consensus)算法[27],FSC在同样的迭代次数下可以去除更多的误匹配点[28]。由于本文需要分别匹配最小矩点和最大矩点,得到两种点的匹配结果会存在一部分重复,需要进行去重处理。

2 实验

2.1 实验准备

本文方法在特征提取阶段需要设置4个阈值提取SAR最小矩点、SAR最大矩点、光学最小矩点和光学最大矩点,设置的原则是使提取出来的SAR最小矩点数和光学最小矩点数尽可能接近,最大矩点遵循同样的原则。在特征描述阶段,将ROEWA和Sobel中的尺度参数β和α设定为常用值β=α=2。在特征匹配阶段,将最近邻距离比阈值tNNDR设置为0.9。需要注意的是,本文将分别使用红线和黄线表示最小矩点和最大矩点匹配结果。实验将高分3号(GF-3)、TerraSAR以及中国科学院电子学研究所研制的某机载SAR上获取的SAR图像和Google Earth获取的光学图像作为实验数据,相关信息如表1所示。

表1 本文使用的SAR和光学图像对数据Table 1 The SAR and optical image pair used in our paper

本文使用正确点匹配数量(correct matching point,CMN)、均方根误差(root mean square error,RMSE)和匹配时间(time)作为匹配算法的性能评估指标,接下来主要介绍CMN和RMSE。

1)CMN 在本文中指的是经过NNDR初始匹配、FSC外点筛选和重复点对去除后的匹配点对数量。在保证较高匹配精度的前提下,匹配点数越多越好。

2)RMSE 主要是用来评估仿射变换矩阵H的精度,公式如下

(17)

(18)

RMSE越小,H越精确,匹配的精度越高。从式(17)看出,CMN是RMSE的某个参数,CMN越高,得到的RMSE越精确。

2.2 实验结果

2.2.1 验证尺度和旋转不变性

本节实验将主要对本文算法的尺度和旋转不变性进行验证。使用的是P-A图像对,因为P-A是同一地理区域范围,同一尺度和同一旋转方向的SAR和光学图像对,很容易将其进行后期处理以便于实现SAR和光学图像对之间不同旋转方向和不同尺度比下的匹配验证实验。

本文算法在不同旋转方向图像对的匹配结果如图7所示,对应CMN结果如表2所示。根据实验结果,本文方法在最大角度差异为90°的情况下仍然能够得到101个正确匹配点数,验证了本文方法的旋转不变性。本文方法在不同尺度比的匹配结果如图8所示,在不同尺度比下的CMN结果如表3所示。根据实验结果,本文算法在尺度比范围为0.7～1.2的P-A图像对中最低能够提取出15对正确匹配点对,验证了本文方法的尺度不变性。

表2 本文算法在不同旋转方向图像对的CMNTable 2 CMN of image pairs in different rotation directions of our method

表3 本文算法在不同尺度比图像对的CMNTable 3 CMN of image pairs in different rotation directions of our method

2.2.2 与其他方法的对比

本文所采用的对比算法分别为PSO-SIFT、SAR-SIFT和OS-SIFT。其中PSO-SIFT是一种专用于遥感图像匹配的改进SIFT算法[14]。SAR-SIFT是专用于SAR图像的匹配算法[13]。OS-SIFT是用于SAR和光学图像的匹配算法,是目前针对SAR和光学图像匹配的最优算法[15]。并且实验对比了结合最大矩点和最小矩点匹配(即本文方法1)和分别采用最大矩点和最小矩点匹配(即本文方法2)的实验结果。为公平起见,实验将对特征点提取阈值进行精调,保证上述方法在特征提取阶段得到数量相同的特征点。图9显示了不同方法的匹配结果,表4对比给出不同匹配方法在5对SAR和光学图像对上的具体性能指标。

表4 不同方法的匹配性能对比Table 4 Comparison of matching performance of different methods

接下来对结果进行具体分析,P-A为城市区域图像对,结构信息最为丰富,包含许多道路和建筑物,但是对应SAR图像的一些区域亮度和对比度较低,使得基于空间域的方法在相应的区域中较难提取到特征,所以PSO-SIFT匹配失败。使用Harris角点的SAR-SIFT和OS-SIFT仅得到12和20对匹配点,而本文方法1和2分别得出116 和214对匹配点,远超其他方法,并且RMSE也低于其他方法。因为本文方法2分别对最小矩点和最大矩点进行匹配,对两种点集进行充分利用,在保证匹配精度的前提下比本文方法1得出了更多的匹配点对数。

P-B为郊外机场区域,图像对具有一定的旋转差异,图像中包括机场跑道和零星的建筑物。结构性区域主要集中于图像中央区域,左下方区域为欠纹理区域。其中左下方的平坦区域明显受到SAR斑点噪声的影响,使得灰度分布不均匀。虽然PSO-SIFT、SAR-SIFT和OS-SIFT能够得到正确的匹配点,但是在左下方区域只有极少的正确匹配点,而本文方法在这两个区域中都能得到更好的匹配结果,并且本文方法2相比于方法1得到了更多的匹配点数和更好的匹配精度。

P-C为山区区域,结构性区域主要集中于中央,在此图像对中PSO-SIFT匹配失败,其他方法都在中央的结构性区域得到了正确的匹配结果,但是相比于SAR-SIFT、OS-SIFT和本文方法1,本文方法2得到更多的正确匹配点对数和较好的匹配精度。

P-D为农田区域,根据P-D图像对的匹配结果,PSO-SIFT匹配失败,SAR-SIFT和OS-SIFT仅得到极少的匹配点对,虽然本文方法1得到了较好的匹配精度,但是本文方法2得到了更多的匹配点对数和更好的匹配精度。

P-E为纹理信息较少的郊区区域,在此类图像对中,本文方法得到了更多的匹配点数和更好的匹配精度,并且本文方法2在具有较低RMSE的情况下得到了更多的匹配点数。

本文方法在通过仿射变换矩阵得到的棋盘拼接图如图10所示,可以看出SAR和光学图像拼接整齐,这将为之后的图像融合和变化检测等工作打好基础。

图10 本文方法2的棋盘拼接图Fig.10 Checkerboard mosaic map of our method

3 结论

本文提出一种联合频域方法的相位一致性计算和空间域方法的RS-GLOH的SAR和光学图像匹配算法,能够对具有非线性辐射和几何差异的SAR和光学图像实现高性能的匹配,主要贡献有:

1)在特征提取阶段从基于相位一致性模型得到的最小矩图和最大矩图中得出分布均匀且位置对应性好的最小矩点和最大矩点,有效抑制了SAR和光学图像之间辐射差异对特征提取的影响。由于两种点的位置区分度较为明显,所以本文对这两种点集合分别进行匹配,能够得到更多的正确匹配点数;

2)在特征描述阶段提出了RS-GLOH描述子,能够有效抑制SAR图像几何失真对特征描述的影响,提取出SAR和光学图像在同一区域的特征描述向量,反映出两者之间的结构相似性。

最终,在特征匹配阶段,使用NNDR和FSC得出最终的匹配结果,考虑到最小矩点和最大矩点匹配结果可能存在重复,进行了去重处理。在5组不同区域的SAR和光学图像对匹配实验中,验证了本文方法的旋转和尺度不变性,相比于PSO-SIFT、SAR-SIFT和OS-SIFT具有更好的匹配精度。