谭宏伟 周芳
(合肥工业大学计算机与信息学院 安徽省合肥市 230009)
合成孔径雷达干涉测量技术[1-2](interferometric synthetic aperture radar,InSAR)利用SAR 传感器获取地面同一区域的两幅或多幅干涉图像,并对图像中对应地面同一散射单元的像素进行相位差分反演来获得地面的高程或形变信息。然而,由于获取干涉图像时SAR 传感器空间位置以及方位向时间的差异,图像间存在平移、缩放和旋转等变形[3],因此,需要对干涉图像进行配准。
SAR 干涉图像配准时通常在参考图像(主图像)上布置一系列控制点,通过求待配准图像(辅图像)在这些控制点上的亚像素级配准偏移量,并对配准偏移量与控制点坐标间的关系进行多项式拟合从而得到主辅图像间的坐标映射函数。目前,根据求取控制点配准偏移量时的测度函数的不同,SAR 干涉图像配准算法主要可以分为互相关函数法[4-6]、频谱比值法[7-8]以及相位平均波动函数法[9]三类,其中互相关函数法是最常采用的配准算法,它具有操作简单、稳健性强的特点。目前全球测绘任务SRTM 及TanDEM-X 的干涉图像配准均是采用互相关函数算法,一些知名的InSAR 数据处理软件也都采用这种算法[10]。
图1:配准算法流程图
互相关函数配准算法根据数据处理的对象是SAR 幅度数据还是包含相位的复数数据可以相应的分为实互相关函数配准和复互相关函数配准。研究表明,实互相关函数配准算法更加稳健,在图像低相干区域也能表现出一定的配准精度,然而其在干涉图像相干性较高且散射特性较为一致的区域配准精度不如复互相关函数配准;复互相关函数配准算法易受图像相干性以及干涉相位变化的影响,对于干涉图上相干性较低或地形变化剧烈的区域其配准性能不如实互相关函数配准[10]。对于SAR 干涉图像,其配准精度一般应优于1/8 像元[11],若配准误差大于这一数值,则干涉图像对的相干性将会受到影响,从而使生成的干涉图相位中含有较多噪声。
为进一步提高SAR 干涉图像的配准精度,本文结合实、复互相关函数配准算法的优势提出一种两级配准算法。提出算法首先采用实互相关函数配准算法对干涉图像进行初步配准,从而得到粗配准辅图像以及低精度的干涉图和相干系数图;然后在相干系数图上提取一系列均匀分布的高相干点作为控制点,并定义相位可靠性指标,以便根据初配准干涉图上的控制点相位可靠性为每个控制点选择相应的精配准方法;最后在初配准的基础上采用选择的精配准方法对控制点配准,得到控制点精确的配准偏移量。通过对实测干涉图像进行配准实验,结果表明提出算法能有效提高SAR 干涉图像的配准精度。
SAR 干涉图像配准中实互相关测度函数定义如下:
式中,s1和s2分别为主辅复数SAR 图像上截取的图像块,(m,n)为坐标位置,||表示求复数模,M、N 为互相关运算的图像大小。
复互相关测度函数定义为:
图2:图像网格区间划分
其中,*表示复数的共轭操作。
主辅SAR 图像块的实互相关函数或复互相关函数越大,两图像块越匹配,相应的图像块中心位置越对齐。
2.2.1 控制点配准偏移量计算
对SAR 干涉图像对配准时首先在主图像上布置一系列均匀分布的控制点,并以这些控制点为中心设置一系列一定大小的目标窗口。在辅图像相应位置处截取同样大小的匹配窗口,通过移动匹配窗口的位置并计算其与目标窗口的互相关函数,当互相关函数达到最大时匹配窗口的中心位置即为控制点的配准位置。由于干涉图像的配准精度需要达到亚像素精度,因此,在控制点亚像素精度配准时需要对匹配窗口进行插值,进而以亚像素网格间距移动匹配窗口,获得亚像素精度的配准偏移量。
2.2.2 偏移模型参数拟合
在得到一系列控制点的坐标偏移量后,就可以对辅图像相对于主图像的坐标偏移模型进行参数估计。偏移模型一般采用如下所示的多项式表示:
2.2.3 辅图像插值重采样
利用坐标偏移模型就可以计算辅图像相对于主图像上每一像素的坐标偏移量,进而得到主图像上每一点在辅图像上的对应坐标位置。由于坐标值往往不是整数,因此,需要对辅图像进行插值重采样。
受SAR 干涉图像成像条件的限制,主辅SAR 图像中往往存在一些低相干区域,当互相关函数配准中控制点位于低相干区域时,配准偏移量的计算精度会降低。然而由于缺乏干涉图像相干性的先验信息,互相关函数配准一般将控制点均匀分布在图像上,无法自主选择高相干的控制点进行配准。另一方面,对于实、复互相关函数,实互相关函数配准的稳健性更强,在低相干区域也能表现出一定的配准精度,然而由于没有利用复数图像的相位信息,其在干涉图像相干性较高且散射特性较为一致的区域配准精度不如复互相关函数配准。复互相关函数配准算法的配准精度容易受图像相干性以及干涉相位变化的影响,对于干涉图上相干性较低或地形变化剧烈的区域其配准性能不如实互相关函数配准[10]。
基于以上分析,本文结合实、复互相关函数配准算法的优势提出一种SAR 干涉图像两级配准算法,其配准过程可以分为以下三个步骤:
(1)利用实互相关函数配准算法计算主辅SAR 图像上一系列控制点的低精度坐标偏移量,并利用这些坐标偏移量进行偏移模型估计和辅图像插值重采样,得到初步配准的辅图像;
(2)由主图像和上一步中生成的初配准辅图像计算低精度的干涉图和相干系数图,然后利用得到的相干系数图在主图像上提取一系列均匀分布的高相干点作为下一步精配准的控制点,并计算干涉图上控制点处的相位可靠性,以便根据相位可靠性在每个控制点处为下一步精确配准选择相应的配准方法;
图3:SIR-C/X-SAR 干涉图像幅度图
图4:SIR-C/X-SAR 干涉图像配准结果
(3)在每个控制点处利用选择的配准方法计算控制点精确的配准偏移量,进而重新进行偏移模型拟合和辅图像插值重采样,得到配准辅图像。
算法流程如图1所示。
图5:ALOS PALSAR 干涉图像幅度图
图6:ALOS PALSAR 干涉图像配准结果
分别用符号S1和S2表示主辅复数SAR 图像。对S1、S2取模,获得其对应的幅度图像。传统基于窗口互相关函数的配准算法在计算辅图像上控制点的配准位置时是通过在辅图像上移动寻找与主图像窗口取得最大互相关函数值的位置,这种方法计算量较大。由于两幅图像互相关函数的傅里叶变换等于一幅图像的傅里叶变换与另一幅图像的傅里叶变换的共轭相乘,因此,可以通过傅里叶变换的方法将窗口最大互相关函数移动查找的过程转变成两图像窗口互相关函数最大值定位的过程,这使得计算效率大大提高。两图像窗口互相关函数的快速傅里叶变换计算方法如下式:
式中,FFT2()和FFT2-1()分别为二维傅里叶正变换和逆变换,FFTSHIFT()表示将图像频谱的零频成分移到频谱中间,y1和y2为两图像窗口数据。通过搜索ρr取得最大值的位置即可得到窗口中心点的配准偏移量:
其中,、 分别为方位向和距离向配准偏移量,M、N 分别为数据窗口的方位向和距离向尺寸,max_x、max_y 分别为ρr取得最大值时的方位向和距离向位置坐标。另外,可以通过傅里叶变换补零的方法来实现亚像素精度的配准偏移量计算。
利用上述方法就可以计算一系列控制点处的辅图像配准偏移量,然后进行偏移模型拟合和辅图像插值重采样,得到初步配准的辅图像,用符号S'2表示。
在配准时,若控制点落入图像的阴影区域,则配准结果可能出现错误,因此还需要进行粗差剔除。本文采用以下两个准则来进行粗差剔除:
(1)粗配准后以控制点为中心的主辅SAR 幅度图像块的互相关函数小于某一阈值(如0.8),则将此控制点配准结果剔除;
(2)此控制点与周围控制点配准偏移矢量长度差的平均值大于周围控制点配准偏移矢量长度的2 倍标准差,则将此控制点配准结果剔除。
将主图像S1与第一步中初配准的辅图像S'2共轭相乘即可得到粗精度的干涉相位图φ:
式中,arg()为求复数相位值。并计算其对应的相干系数图ρ:
式中,(i,j)为当前计算相干系数的像素位置,Ω 为以(i,j)为中心的矩形窗口,(k,l)为窗口内的像素坐标。
本文设置粗配准相干系数图的均值作为进一步配准的控制点相干系数阈值,提取相干系数图上高于这一阈值的像素作为预选控制点。InSAR 图像配准中,需要控制点在图像上趋于均匀分布,这样有利于偏移模型的正确估计。因此还需对预选控制点进行离散化,避免控制点集中分布在某一局部高相干区域。本文通过以下步骤来使预选控制点离散化:
(1)选择1/10 主图像图幅大小作为控制点离散化的窗口大小(当图像尺寸较大时,可以缩小离散化窗口相对于图像的比例,以增加控制点选择的个数),并以此离散化窗口将原图像依次划分为一系列网格区间,如图2所示。将每一网格区间内相干系数最大的预选控制点作为控制点,并记录这些控制点的坐标位置。
(2)设置相邻网格区间内的控制点行列方向的间距最小值分别为离散化窗口行列尺寸的1/2,若当前网格区间内的控制点与其相邻网格内的控制点在行方向或列方向的间距小于这一设定值,则向使行或列间距增大的方向选取下一个预选控制点作为控制点。重复执行相邻网格区间内所选控制点的间距检测和重新设置工作,直至间距满足设定值。若当前网格区间中没有点满足要求,则此网格区间不选择控制点。
表1:SIR-C/X-SAR 干涉图像配准性能分析
表2:ALOS PALSAR 干涉图像配准性能分析
(3)在每个网格区间内重复步骤(2),直至完成所有网格内控制点的选取。
选出一系列高相干控制点后,由于控制点的相干性较高,且粗配准得到了较为准确的干涉相位图,因此可以在这些控制点所在的局部邻域内估计粗略的干涉相位项,并将估计的干涉相位项从主SAR 图像中去除,进而利用复互相关函数配准算法就可以得到控制点更加精确的配准偏移量。然而在实际数据处理时,当主辅SAR图像相干质量很差时,选择的控制点相干性可能依然不佳,复互相关函数配准不如实互相关函数;另一方面,当控制点处于非匀质区域时,粗配准得到的干涉相位图上控制点邻域范围内含有较多的相位噪声,此时难以由含噪干涉相位估计真实相位,前述方法的配准误差较大。
在理想干涉图中,局部干涉图频谱会出现较明显的峰值,并且频谱峰值会随着噪声等级的增加逐渐变小。据此,本文定义一个相位可靠性指标来评价粗配准干涉相位图在控制点邻域内的相位可靠性。相位可靠性计算方法如下:
式中,max()为求矩阵最大值,sum()为对矩阵求和,P(u,v)为控制点邻域范围内的干涉图频谱。设置相位可靠性阈值为rth,当rp≥rth时,粗配准控制点邻域范围内的干涉相位比较可靠,采用前述方法计算控制点处精确的配准偏移量,当rp 利用为每一控制点选择的配准方法就可以在初配准的基础上计算精确的配准偏移量。当采用实互相关函数进行配准时,配准过程与第一步类似,不再赘述,当采用相位补偿的复互相关函数进行配准时,数据处理过程如下: (1)控制点邻域范围内的干涉相位估计。由含噪干涉相位估计实际相位的操作即为干涉图滤波,然而本步只是估计大致的干涉相位项,从而将其从主图像中去除,以减弱干涉相位对复互相关函数配准精度的影响,因此为了减小计算量,本文采用多视滤波的方法进行相位估计,计算方法如下: (2)主图像干涉相位去除。在主图像相应控制点位置处的同等大小窗口内将估计的干涉相位去除,可以得到去除干涉相位后的主图像窗口数据: 式中,c1为主图像控制点处的数据窗口,φp为估计的干涉相位。 (3)精确配准偏移量计算。同样,在粗配准的辅图像相同位置处提取同等大小的数据窗口c2,然后采用复互相关函数配准算法计算c'1和c2的中心像素配准偏移量,并将此偏移量的值与第一步中估计的粗偏移模型计算的偏移量相加,得到控制点精确的配准偏移量。 得到一系列控制点精确的配准偏移量后,采用最小二乘方法对主辅图像间的多项式偏移模型进行参数估计,在此基础上计算主图像每一点在辅图像上的相对坐标位置,进而对辅图像插值重采样,得到配准辅图像。 为验证提出算法的有效性,本文对两组不同数据源的干涉图像对进行配准实验,并将本文算法与实互相关函数配准算法、复互相关函数配准算法进行对比。 第一组实验数据是由SIR-C/X-SAR 系统获取的意大利Etna 火山区域图像,主辅图像成像日期分别为1996年10月8日和1996年10月9日,图像大小均为1024×1024 像素,地面分辨率约为30 米,主辅SAR 幅度图像如图3所示。由于成像区域为山区,地物多为植被且高程变化剧烈,因此理论上图像对相干性较低且干涉相位变化剧烈。 分别采用实互相关函数、复互相关函数以及本文算法进行干涉图像配准,参数的设置情况为:配准窗口大小设为51 像素×51 像素,其中实、复互相关配准以及本文算法的初步配准阶段均在主图像上均匀布置400 个控制点,而在本文算法的高相干控制点查找阶段设置离散化窗口的大小为图像大小的1/20。生成的干涉相位图以及相干系数图如图4所示。由图可以看出,复互相关函数配准得到的干涉图的相位条纹比较模糊,在图像右下部分几乎没有相位条纹出现,对应区域的相干系数图也较暗,配准效果较差。从目视效果来看,实互相关函数配准算法以及本文算法的配准结果较为相似,都得到了较为清晰的干涉相位图以及亮度更高的相干系数图。 分别计算不同算法的配准偏移量多项式拟合时的距离向、方位向均方根误差以及配准结果的残差点数和平均相干系数,结果如表1所示。可以看出,提出算法的均方根误差最小、平均相干系数最大,残差点数目也最少,具有最优的配准性能。 第二组实验数据是ALOS 卫星上的PALSAR 传感器获取的西安地区图像,主辅图像成像日期分别为2009年8月7日和2009年9月22日,图像大小均为30927×4903 像素,方位向、距离向分辨率分别约为3m 和15m。为了便于图像显示,分别在主辅图像中心位置处截取1024×1024 像素大小的图像块进行配准实验,其幅度图像如图5。 采用本文算法以及两种对比算法对干涉图像进行配准,配准参数设置情况与实验1 相同,得到的干涉相位图以及相干系数图如图6所示。可以看出在本实验中,实互相关函数配准算法对应的相干系数图在图像右上部分较暗,相应区域相位条纹模糊,配准效果不佳,而复互相关函数配准以及本文算法配准得到的结果总体一致。 计算各种算法配准结果的性能,结果如表2所示。可以看出,提出算法在三个性能指标上同样表现最优。 本文提出一种联合实、复互相关函数的干涉图像两级配准算法,算法首先利用较为稳健的实互相关函数配准算法对图像初配准,然后利用初配准的结果选择高相干的点作为控制点并定义控制点相位可靠性,从而利用相位可靠性为精配准选择相应的配准方法。由于初配准提供了配准图像的基本信息,因此在此基础上进行进一步的配准能得到更高的配准精度。通过对两组干涉图像进行配准实验,并从目视效果以及定量性能指标上评价配准的效果,结果证明了提出算法的有效性。3.3 控制点精确配准偏移量计算和辅图像插值重采样
4 实验结果与分析
4.1 SIR-C/X-SAR干涉图像配准
4.2 ALOS PALSAR干涉图像配准
5 结论