2. 国防科技大学 电子科学与工程学院,湖南 长沙 410073
2. School of Electronic Science and Technology, National University of Defense Technology, Changsha 410073, China
1 引 言
干涉合成孔径雷达技术是以不同观测几何下获取的两幅或多幅SAR复图像数据的干涉相位为信息源进而反演得到地表的三维地形[1, 2, 3]。干涉相位的质量将决定InSAR系统的最终产品精度。SAR复图像的高精度配准是获得高质量干涉相位的前提。因此,复图像配准是干涉SAR数据处理中最为关键的步骤之一。
为了从SAR复图像对中获取干涉相位信息,就要确保用于干涉计算的图像像素对应于地面同一散射单元。由于两幅SAR图像在获取时观测几何存在差异,使得图像之间存在着一定的偏移和扭曲,需要对图像进行配准处理。在干涉SAR数据处理中,一般要求配准误差小于0.1像素。对于全球干涉测量,SAR复图像配准可以分为两步:① 几何配准[4],根据成像几何关系,利用平台轨道(天线相位中心)数据和成像参数计算得到两幅图像的配准偏移量,这一步骤的配准精度与轨道数据质量以及有无辅助DEM数据紧密关联,早期的星载数据几何配准精度一般为20~30像素,目前TerraSAR-X/TanDEM-X的几何配准精度优于10像素;② 图像配准,这一步骤是针对SAR图像进行处理,一般情况下又可分为像素级配准和亚像素级配准两个环节。像素级配准是基于某一种配准测度准则,计算两幅SAR图像在不同偏移位置上的配准度量值,由此得到1像素精度的配准结果。然而,随着雷达观测模式的多样化和图像分辨率的提高,已经很难保证两幅图像的像素级配准结果在整个测绘带内满足1像素的精度,只能通过分块进行像素级配准处理才能达到要求。亚像素级配准和像素级配准较为类似,不同的是亚像素级配准先对SAR图像进行插值处理后再作配准测度计算,或者直接对配准测度计算结果进行插值,从而达到亚像素级的配准精度。
目前,针对干涉SAR图像配准问题的研究主要是针对图像配准这一步骤进行的。根据配准所采用的测度函数,已有的干涉SAR图像配准方法主要可分为3大类:相关函数法[5, 6, 7]、平均波动函数法[8]和最大频谱法[9, 10]。相关函数法是以图像对之间的相关函数作为配准度量,当该值最大时认为图像对配准。平均波动函数法是以干涉条纹两个方向的相位梯度绝对值之和作为配准度量,当该值最小时认为图像对配准。最大频谱法是以干涉图频谱模值最大值与其余频率模值之和的比值作为配准度量,当该比值最大时认为图像对配准。除了这3大类算法以外,还有一些其他的混合改进以及基于图像特征的配准算法[11, 12, 13]。已有的这些配准算法主要面临着一个问题,算法在提出时是针对某一类型图像取得了很好的配准效果,但对其他数据的普适性和稳健性有待进一步考究。
星载InSAR系统的一个重要任务是进行全球测绘,因此各环节处理算法的稳健性至关重要。随着星载雷达技术的发展,各种观测模式的出现,尤其是分辨率的提高,研究快速、高精度、稳健的SAR图像配准算法具有十分重要的意义。基于此,在分析各种配准测度函数特性的基础上,本文提出了一种稳健、高效的星载干涉SAR复图像配准方法,该方法能适用于各种类型的SAR图像配准任务。
2 配准测度函数的选择配准方法的核心在于配准测度函数的选择,目前的配准测度函数主要有:相关函数、平均波动函数和频谱比值。文献[9]针对干涉图在频域的特性提出了基于频谱比值大小的最大频谱配准方法,认为当图像配准时,在干涉图的频谱图上会存在一个明显的峰值。然而,当频谱中存在两个以上的主要频率或干涉条纹质量较低导致没有明显的峰值时,最大频谱法将出现不稳定。从干涉条纹的清晰程度出发,文献[8]提出一种新的配准测度函数,称为平均波动函数(average fluctuation)
式中,φi,j为干涉相位;W配准计算窗口。当两幅图像配准时,干涉条纹最清晰,平均波动函数值最小。但是,在低相干高噪声区域或亚像素级配准时,平均波动函数值变化不够灵敏,配准误差的方差较大。
相关测度是一种最基本的统计方法,广泛应用于各种类型的图像配准,是许多配准算法的基础,它具有操作简单且稳健性强的特点。目前仅有的两次全球干涉测绘任务SRTM[14]及TanDEM-X[15, 16]的干涉数据处理模块的配准算法均采用相关测度函数。另外,一些知名成熟的InSAR数据处理软件的配准算法也是采用相关函数法。因此,相关函数是一种较为稳健的配准测度函数,适合全球测绘任务下的干涉SAR图像配准。在干涉SAR复图像配准时,相关函数的计算又可分为实相关计算和复相关计算两种,两幅SAR复图像的实相关函数定义为
式中,s1和s2表示两幅SAR复图像;M、N表示用于相关计算的窗口大小;(u,v)表示相关计算的滑动位置;| |表示取模操作。复相关函数定义为
式中,*表示复共轭。需要指出的是,当干涉条纹较为密集时,考虑两幅图像之间的相位差异,可以利用辅助DEM或其他技术手段对式(3)进行相位补偿。研究发现,实、复相关测度函数各自的特点如下:
(1) 对于高相干且散射特性较为一致的区域,复相关函数的配准精度优于实相关函数,实相关的配准误差约为复相关配准误差的倍[17]。
(2) 对于存在明显特征地貌的区域,实相关函数的配准精度要高于复相关函数[18]。
(3) 在实际干涉数据的低相干区域,实相关函数的性能要优于复相关函数。
下面利用实测数据来验证实、复相关函数的上述特点。从SIR-C/X-SAR系统获取的实测数据中截取了3块数据(64像素×64像素),分别对应于含有特征地貌区域、不含特征地貌区域以及低相干区域,如图 1(a)、(d)和(g)所示。图 1(b)~(c)给出了含有特征地貌的SAR图像对的滑动实、复相关测度函数计算结果,图 1(e)~(f)给出了不含特征地貌、散射特性较为一致的SAR图像对的滑动实、复相关测度函数计算结果,图 1(h)~(i)给出了低相干区域SAR图像对的滑动实、复相关测度函数计算结果。
在图 1中,实、复相关函数的峰值位置、峰的大小及分布均略有不同,都有各自适应的情况,很难确切地说两种相关函数孰优孰劣。一种合理的思路是将两种相关函数进行联合,不同情况下采取不同的相关函数作为配准测度函数。
因此,设计一种准则,使得在这种准则下算法能自适应地选取合适的相关函数,是配准算法设计的关键。为此,提出了配准灵敏度准则,配准灵敏度包括方位向灵敏度和距离向灵敏度,是指配准测度函数归一化后在其峰值位置沿方位向和距离向剖面的3dB宽度[19],类似于SAR图像分辨率的定义。图 2给出配准灵敏度示意图,图 2(a)是归一化后测度函数,图 2(b)、(c)分别给出了在其峰值位置沿方位向和距离向的剖面图,对应的3dB宽度(半功率点处的宽度),即为配准测度函数的方位向灵敏度和距离向灵敏度。
由此可知,配准灵敏度越小表明该配准测度函数的配准精度越高。因此,在实际SAR图像对配准过程中,可以依据配准灵敏度来选取相应的(实、复)相关函数作为配准测度函数。
3 实、复相关联合配准方法流程在上一节中,分析了各种测度函数的优缺点,并确定了稳健的配准方法应该以相关函数作为其配准测度函数。针对存在实、复两种相关函数,提出配准灵敏度准则,用以选择合适的相关函数。本节给出基于联合实、复相关函数的干涉SAR图像配准方法的详细操作流程。
当SAR图像对精确配准时,其相关函数达到最大值。所以,只需在配准窗口内所有偏移位置上计算相应的相关函数值,等价于计算两配准窗口图像的互相关函数,互相关测度函数的峰值位置决定了配准偏移量。根据相关定理[20],两幅图像互相关函数的傅里叶变换等于一幅图像的傅里叶变换与另一幅图像的傅里叶变换的共轭相乘。因此,实、复相关函数可以通过快速傅里叶变换实现
式中,norm( )表示归一化操作,这样就把逐点的滑动窗相关运算转化为图像块之间的互相关计算,得到数据块间的滑动相关函数,极大地提高了计算效率。综上所述,联合实、复相关函数配准算法的详细步骤如下。
步骤1:SAR图像像素级粗配准。由于实相关函数较为稳健,因此两幅SAR图像之间的像素级配准基于实相关完成。像素级配准的目的是使得后续亚像素级配准的输入窗口能尽可能多的重叠(相关)。根据输入数据的图幅大小、分辨率高低以及成像工作模式决定是否需要分块进行像素级配准。在条带模式下,实测数据处理中一般以5000像素×5000像素(方位向×距离向)分块进行粗配准。
步骤2:SAR图像分块。传统的配准方法一般是通过在SAR图像上布置控制点,对控制点所在子图像窗口进行配准处理,得到所有控制点处的配准偏移量,然后结合多项式模型即可求解得到原始图像任意像素处的配准偏移量。然而,随着星载雷达技术的发展,各种观测模式的出现,尤其是分辨率的提高,多项式模型在某些情况下已经不再满足实际配准要求[16]。此时,可以将SAR图像按照某一大小(如128像素×128像素)进行分块,对每一子块进行亚像素级配准,得到相应的配准偏移量,然后通过局部曲面拟合或内插的方式得到整幅图像任意像素处的配准偏移量。
步骤3:SAR图像子块亚像素级配准。对每一子块图像对,分别计算其实、复相关函数,亚像素级配准是通过傅里叶变换的补零实现的,得到亚像素级实、复相关函数后,计算其相应的配准灵敏度,并基于此选择相应的方位向和距离向配准测度相关函数,从而得到相应的配准偏移量(峰值位置)。
步骤4:子块亚像素级配准偏移量的粗差剔除。得到每一子块图像的亚像素级配准偏移量后,可以通过以下准则剔除粗差:① 方位向或距离向的配准灵敏度大于某一门限值;② 配准偏移量的矢量长度超过周围偏移量矢量长度的2倍标准差范围;③ 配准偏移量的矢量长度超过周围偏移量矢量长度的2倍中值范围。在实际粗差剔除时,任选上述准则的一种或两种进行组合均可。被剔除子块的配准偏移量可以通过对其周围偏移量的局部曲面拟合(或插值)得到。
步骤5:辅图像重采样处理。对子块配准偏移量进行插值处理就可求得整幅图像任意像素处的亚像素级配准偏移量。然后基于某一插值核函数对SAR辅图像进行重采样,完成两幅SAR图像间的精确配准。
4 实测数据处理与分析为验证算法的有效性,利用该方法对实测数据进行处理。所用实测数据是ALOS-PALSAR于2009-07-01和2009-08-16获取的重复轨道干涉数据,基线长度约为2700 m,相干系数约为0.65。图 4(a)显示了该区域的SAR幅度图像,整个观测区域呈山地形貌,且地物类别较为丰富,包括机场(SAR图像左下角暗黑区域)、大面积城市建筑区(SAR图像整个右侧较亮区域),水域(SAR图像中部暗黑区域),比较适合校验配准算法的性能。
首先对两幅SAR图像进行像素级配准,从主、辅图像中间各取一块大小为1024像素×1024像素的区域,基于实相关测度函数进行像素级配准。然后,将SAR图像进行分块,分别计算其方位向和距离向的实、复相关函数配准灵敏度(傅里叶变换时作16倍插值),结果如图 3所示。为了更好地显示配准方案的中间结果,本文在处理这一景数据时,没有按照第3节中步骤2的方法(按128像素×128像素)进行分块,而是在方位向取18个子块,距离向取10个子块,总共180个子块,每个子块的大小均为128像素×128像素。
在图 3中,横坐标表示用于配准的数据子块索引号,纵坐标表示配准灵敏度(像素数目)。可以看到,复相关函数的配准灵敏度均值要略小于实相关函数,这也验证了一般情况下复相关配准测度函数的配准精度略高于实相关函数这一规律。但是,在部分区域,实相关函数的配准灵敏度小于复相关函数,这也说明了联合实、复相关函数的必要性和有效性。
计算出各数据块的实、复相关配准灵敏度后,利用第3节的方法即可确定各块的实、复相关配准测度函数分布,结果如图 4(b)所示。在总共180个子块中,方位向、距离向均采用复相关作为配准测度函数的有111块,方位向、距离向均采用实相关作为配准测度函数的有29块,方位向、距离向分别采用实、复相关函数作为配准测度函数的有6块,方位向、距离向分别采用复、实相关函数作为配准测度函数的有20块,依据粗差剔除准则①、②剔除了14个子块的配准偏移量(方位向距离向配准灵敏度门限值均取46),通过局部曲面拟合的方式得到剔除位置的偏移量,结果如图 4(c)所示。对图 4(c)所示的配准偏移量进行双线性插值可以得到所有像素处精确偏移量,然后对SAR辅图像进行重采样处理。图 4(d)显示的是配准完成后去平地效应后的干涉相位图。
图 4(b)中,■ 表示方位向、距离向均采用复相关函数作为配准测度函数;● 表示方位向、距离向均采用实相关函数作为配准测度函数;▼ 表示方位向、距离向分别采用实、复相关函数作为配准测度函数;▲ 表示方位向、距离向分别采用复、实相关函数作为配准测度函数;★ 表示剔除的粗差。值得注意的是,由图 4(d)可以看出,图像右侧的干涉条纹清晰程度相比于其他区域要差。这是因为该区域是城区,在地形干涉条纹上叠加了城市建筑高度所对应的高程相位,而城区建筑高程变化是非连续的,所以图像右侧的干涉条纹清晰度(连续程度)在视觉效果上差于其他区域。在1∶1比例的干涉相位图上可以清晰地看到城市建筑所对应的干涉相位变化。为了更进一步说明联合实、复相关函数的优点,还分别利用实相关和复相关作为配准测度函数对该图像对进行了配准处理,不同方法配准结果的残差点数目及相干系数统计结果如表 1所示。计算相干系数所采用的方法是Guarnieri相干估计法[21],估计窗口大小为9像素×9像素。
由表 1可知,实、复相关函数联合的配准性能明显优于单一的实相关或复相关配准。而相关函数是一种稳健的配准测度函数,因此,本文联合实、复相关的配准方法能够稳健、高效地适应于全球干涉数据配准处理。
5 结 论本文提出了一种基于联合实、复相关函数的干涉SAR图像配准方法。在分析各种配准测度函数特性的基础上,确定相关函数作为匹配度量,并提出了配准灵敏度概念,从而能够自适应地选取配准测度函数。该方法由于联合了实、复相关函数各自的优点,使得在匹配的准确度和稳定性上较传统方法有所提高。针对ALOS-PALSAR复杂场景实测数据的处理结果表明,该方法能稳健、高效地适应于各种地貌类型的实测数据处理。需要指出的是,在应用过程中,配准偏移量粗差剔除可以选择文中给出准则中的一种或两种进行组合均可,3种准则的门限值可以根据实际情况进行合理的取值,不同的取值会对结果略有影响,本文给出了门限取值的参考值,如何获得最优取值也是下一步需要研究的问题。
致 谢:感谢德国宇航中心(DLR)和日本宇宙航空开发研究机构(JAXA)提供的SIR-C/X-SAR、ALOS-PALSAR数据。
[1] | KRIEGER G,MOREIRA A,FIEDLER H,et al.Interferometric Synthetic Aperture Radar (SAR) Missions Employing Formation Flying[J].Proceedings of IEEE,2010,98(5):816-843. |
[2] | ROSEN P A,HENSELEY S,JOUGHIN I R,et al.Synthetic Aperture radar Interferometry[J].Proceedings of IEEE,2000,88(3):333-382. |
[3] | WERNER M,Shuttle Radar Topography Mission (SRTM) Mission Overview[J].Proceedings of IEEE,2001,55(2):75-79. |
[4] | SANSOSTI E,BERARDINO P,MANUNTA M,et al.Geometrical SAR Image Registration[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(10):2861-2870. |
[5] | LI F K,GOLDSTEIN R M.Studies of Multi-baseline Space-borne Interferometric Synthetic Aperture Radars[J].IEEE Transactions on Geoscience and Remote Sensing,1990,28(1):88-97. |
[6] | STONE H S,ORCHARD M T,CHANG E C,et al.A Fast Direct Fourier-based Algorithm for Subpixel Registration of Images[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(10):2235-2243. |
[7] | ZENG Qiming,XIE Xuetong.A FFT-based Complex Correlation Function Method Applied to Interferometric Complex Image Corregistration[J].Acta Gendaetica et Cartegraphica Sinica,2004,33(2):127-131.(曾琪明,解学通.基于谱运算的复相关函数法在干涉复图像配准中的应用[J].测绘学报,2004,33(2):127-131.) |
[8] | LIN Q,VESECKY J F,ZEBKER H A,et al.New Approaches in Interferometric SAR Data Processing[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(3):560-567. |
[9] | GABRIE1 A K,GOLDSTEIN R M.Crossed Orbit Interferometry:Theory and Experimental Results from SIR-B[J].International Journal of Remote Sensing,1988,9(5):857-872. |
[10] | ZHAO Zhiwei,YANG Ruliang,QI Haiming.An Improved Maximum Spectrum Peak Co-registration Algorithm for Space-borne InSAR Complex Data[J].Acta Gendaetica et Cartographica Sinica,2008,37(1):64-69.(赵志伟,杨汝良,祁海明.一种改进的星载干涉SAR复图像最大频谱配准算法[J].测绘学报,2008,37(1):64-69.) |
[11] | SHI Xiaojin,ZHANG Yunhua.A New Image Registration Method for Repeat-Pass InSAR Based on Fourier-Mellin Transformation and Correlation-Coefficient Algorithm[J].Journal of Electronics and Information Technology,2009,31(4):803-807.(石晓进,张云华.基于Fourier-Mellin变换和相干系数法的重复轨道干涉SAR图像配准新方法[J].电子与信息学报,2009,31(4):803-807.) |
[12] | CHEN Lifu,WEI Lideng,XIANG Maosheng,et al.Auto-registration Imaging Algorithm of Non-linear Approximation for Airborne Dual-antenna InSAR[J].Journal of Electronics and Information Technology,2010,32(9):2208-2214.(陈立福,韦立登,向茂生,等.机载双天线干涉SAR非线性近似自配准成像算法[J].电子与信息学报,2010,32(9):2208-2214.) |
[13] | ZOU Weibao,LI Yan,LI Zhilin,et al.Improvement of the Accuracy of InSAR Image Co-registration Based on Tie Points-A Review[J].Sensors,2009,9:1259-1281. |
[14] | RABUS B,EINEDER M,ROTH A,et al.The Shuttle Radar Topography Mission-a New Class of Digital Elevation Models Acquired by Space-borne Radar[J].ISPRS Journal of Photogrammetry and Remote Sensing,2003,57:241-262. |
[15] | FRITZ T,BREIT H,BALSS U,et al.Processing of Interferometric TanDEM-X Data[C]//Proceedings of the European Conference on Synthetic Aperture Radar.Aachen:[s.n.],2010:412-415. |
[16] | MARTINEZ N Y,EINEDER M,BRCIC R,et al.TanDEM-X Mission:SAR Image Coregistration Aspects[C]//Proceedings of the European Conference on Synthetic Aperture Radar.Aachen:[s.n.],2010:576-579. |
[17] | BAMLER R.Interferometric Stereo Radar-grammetry:Absolute Height Determination from ERS-ENVISAT Interferograms[C]//Proceedings of the European Conference on Synthetic Aperture Radar.Munich:[s.n.],2000:742-745. |
[18] | LEONG K K,EE C C,WANG C A H,et al.DTM Generation from 35-day Repeat Pass ERS-1 Interferometry[C]//Proceedings of IEEE International Geoscience and Remote Sensing Symposium.[S.l]:IEEE,1994:2288-2290. |
[19] | SKOLNIK M I.Handbook of Radar[M].WANG Jun Translate.Beijing:Publishing House of Electronics Industry,2003.(SKOLNIK M I.雷达手册[M].王军译.北京:电子工业出版社,2003.) |
[20] | OPPENHEIM A V,WILLSKY A S.Signals and Systems[M].New Jersey:Prentice-Hall,1997. |
[21] | GUARNIERI A M,PRATI C.SAR Interferometry:A "Quick and Dirty" Coherence Estimator for Data Browsing[J].IEEE Transactions on Geoscience and Remote Sensing,1997,35:660-669. |
[22] | WANG Qingsong.Research on High Efficiency and High Precision Processing Techniques of Spaceborne Interferometric Synthetic Aperture Radar[D].Changsha:National University of Defcnse Technology,2011.(王青松.星载干涉合成孔径雷达高效高精度处理技术研究[D].长沙:国防科学技术大学,2011.) |