王 洒,宫辉力,杜钊峰,顾兆芹
1.首都师范大学 三维信息获取与应用教育部重点实验室,北京100048;2.首都师范大学 资源环境与旅游学院,北京100048;3.郑州轻工业学院 计算机与通信工程院,河南 郑州450002;4.河南城建学院 测绘与城市空间信息系,河南平顶山467044
随着对地观测技术的发展,干涉测量技术得到不断的改进和优化,在InSAR基础上发展起来的永久散射体干涉测量技术(persistent scatterers for SAR interferometry,PSInSAR)的应用范围逐步扩大[1-3]。它在同一地区的多景SAR图像的基础上,利用数学算法把在长时间内散射特征较为稳定并且回波信号较强的点(有角反射效应的人工建筑物、裸露的岩石等)统计出来,这些稳定的点就是永久散射体(permanent scatterer,PS),它们的相干性较高,并且不受大气变化影响。相对于水准测量等常规监测方法,PSInSAR技术不仅可以大范围获取高时空分辨率的地表形变细节信息,并且在理论上有着毫米级监测精度,因此被广泛应用于地表形变场监测(如地震、地面沉降以及山体滑坡引起的地表位移)[4-7]。
使用PSInSAR监测地表形变场变化时,为了提高监测的精度,在处理过程中往往会使用十几幅甚至几十幅图像[5]。选取其中一幅为主图像,其他所有图像为从图像,这些从图像分别与主图像配准并重采样,与主图像形成干涉对。如果主图像选取不恰当,就会影响到图像干涉质量,进而导致PS点选取不准确,测量结果不可靠。为了保证数据的正确性,就要重新选取主图像,重复进行SAR数据处理,这样既费时又费力。关于主图像的选取方法,国内外学者有一系列的研究:文献[8]提出SAR图像的相干性与时空基线、多普勒质心频率以及热噪声等因素有关;文献[9]在时间基线、空间有效基线以及多普勒质心频率差3个因素的基础上建立综合相关函数求取最优公共主图像;文献[10]把时间基线和空间基线作为选取主图像的主要影响因子;文献[11]把GPS数据融合在主图像选取模型中以便于去除大气延迟影响。上述研究可概括为通过建立模型获取理论上最优主图像,但研究中并未关注主从图像间的干涉效果,且无后续试验来验证所选取的主图像是否为最优。文献[12]提出把聚类分析用于主图像的选取,该方法虽然从考虑干涉对集合入手,但带有一定的随意性,需要人为设定阈值来判断影像是否相干。
基于此,本文研究内容主要包括3个方面:首先在前人研究的基础上,除了综合考虑时间基线、空间基线、多普勒质心频率以及热噪声等因素外,把季节性变化也融合到最佳主图像选取函数中去;接着利用此函数计算试验数据中主从图像之间的相干性以及每幅影像为主图像时可相干从影像个数,在此基础上预测最佳主图像;最后用软件实际处理结果来验证最佳主图像选取方法的有效性。
根据PSInSAR的基本原理,地表同一分辨单元在不同的成像时刻所形成的回波信号存在着很多去相干因素。主从图像之间的相干性受雷达成像(图像热噪声、基线去相关、时间去相干、多普勒质心去相干等)和一系列处理过程的去相干因素(数据处理去相干)的影响[8-13,15-16]。这些因素降低了干涉图的质量,增加了信息反演的难度。
最佳主图像选取为永久散射体干涉测量的第1步,不涉及数据处理去相干的影响[13],只需考虑初始状态下主从图像之间的相干性。因此,图像之间的相干性在各因子的综合作用下可以表示为[14-15]
式中,rtot为总相干性;rSNR、rgeom、rDC、rtemporal分别表示系统热噪声、影像间几何、多普勒质心以及时间相干程度。
在主图像的选取中,必须综合考虑到所有SAR图像之间的系统热噪声、时间基线、空间基线以及多普勒质心频率等因素,从中选取与其他所有图像相干程度最高的SAR图像作为最佳主图像,以达到最好的干涉效果。上述4个影响因子都可以用函数表达。
2.1.1 系统热噪声函数
信噪比(SNR)与星载SAR系统特征有关,当两幅SAR图像的SNR值不同时,系统热噪声函数表示为[9,13,15]
SNR的值可以根据卫星的系统参数确定,对于固定模式的卫星可以视为一个常数。比如ERS-1/2卫星的SNR值为11.7。因此f(rSNR)也视为一个常数。
2.1.2 SAR图像间几何相干性函数
SAR图像间的几何关系如图1所示,其几何相干性通常与垂直基线有关。在一定的垂直基线距内,SAR图像之间可以产生干涉。垂直基线越长,相干性越低,当垂直基线超过临界范围,SAR图像 完 全 不 相 干[9,11,13,15]。因此,图像间的几何相干性表示为一个分段函数
式中,B⊥为SAR图像之间的垂直基线;B⊥C为临界垂直基线,表示为
式中,W为雷达带宽;R是斜距;θ是入射角;Ωslope为斜面的坡度;f0是频率。
图1 SAR图像间几何关系Fig.1 The geometrical relationship between SAR images
2.1.3 SAR图像间多普勒质心频率相干性函数
对于SAR图像来说,多普勒质心频率差越大,相干性越低,当多普勒质心频率差超过临界值时,图像完全不相干[9-11]。和SAR图像间几何相干函数一样,图像间多普勒质心频率相干性可表示为分段函数[9-11]
式中,f为图像之间的多普勒质心频率差;fC为多普勒质心频率差临界值。
2.1.4 考虑季节性因素的时间相干性函数
时间去相干是一个比较复杂的因素,主要取决于环境因子,例如降雨和植被变化等[13]。SAR图像成像时间间隔越长,时间失相关越严重[9]。在以往的研究中,通常会把时间相干表示为[9-11,13,16]
式中,T为图像之间的时间基线差;TC为临界相干时间基线。
仅靠SAR图像成像时间间隔的长短来评判时间失相关,往往会忽略季节性影响。事实上,季节性变化更易导致图像失相干[17-18]。地球上大部分区域都存在着季节性变化,由此导致的环境差异会对干涉效果产生很大的影响。因此对于地表建筑群无明显变化的地区,在时间相干函数中也要考虑到季节性因素的影响。本文在图像成像时间差别和季节性差别的基础上,建立融合季节性因素的时间相干函数。具体推导过程如下:
步骤1 对时间基线取整,获得两幅干涉图像之间的相隔年份整数值
步骤3 时间基线实际相隔年份数减去年份的中点值,获取实际值和中点值的距离,并求其绝对值
步骤4 离中点值越远,季节性变化影响越小,中点值左右相同间隔内季节性图像可以看做近似相等,表示为
步骤5 季节性变化较小的图像之间有较大的相干性,表示为
因此融合季节性因素的时间相干函数表示为
简化为
根据SAR图像之间相干性函数(公式(1)),对于2幅图像来说,图像间的相干可以表示为
在前文中已知f(rSNR)为一个常数,它的变化对主图像的选取并没有决定性的影响,所以在表达式中将此项省略。本文研究区存在季节性变化,用融合季节性变化的时间相干函数rseason来表达时间相干性。优化后的图像间相干函数为
同理,在最佳主图像的选取过程中,对于N幅图像,每幅SAR图像都要依次作为主图像按照公式(9)分别计算与其他N-1幅从图像间的相关性,最终生成一个N×N的相干系数矩阵。
图像的干涉质量要通过可干涉图像数目的多少以及干涉图像相干性的高低来衡量。因此在最佳主图像的选取过程中,要根据相干系数矩阵提供的信息,考虑两个因素,完全不相干图像数目以及相干图像的平均相干系数。衡量是否为最佳主图像,要遵守以下两个规则:
(1)SAR图像的完全不相干图像(相干系数为零)数目最少,为了在后文中表述方便,完全不相干图像数目用字母D来表示;
(2)在满足规则(1)的情况下,图像平均相干系数最高,相干系数均值R可以通过以下公式求得
式中,D越小表明该图像为主图像时与其从图像总体相干性越好;R越大表明相干系数越高。因此,D越小,R越大,被选为最佳主图像的概率越大。
本文选取2002—2007年拉斯维加斯地区13幅Envisat ASAR升轨数据来验证最佳主图像选取方法的有效性。研究区范围为10km×10km,具体如图2所示。
图2 试验数据范围Fig.2 Area of SAR images
通常情况下,对于依靠实用时间模型来选择PS像素,最少需要25幅干涉图像[19]。文中使用的软件为 StaMPS/MTI(Stanford Method for Persistent Scatterers/Multi Temporal InSAR),该软件仅用12幅SAR图像就可以确定PS点的网络分布[20]。表1为2002-12-12的SAR 影像为主图像时试验数据参数。
表1 试验数据参数Tab.1 Parameters of experiment data
前述的最佳主图像影响因子所涉及的相干性函数表达式中存在分段性,需要依靠其临界值确定函数最终表达状态。根据试验数据参数特征,影响因子临界值计算方法如下。
3.2.1 垂直基线距
只有图像间的垂直基线在临界范围内,图像才可相干。ASAR数据在Image模式下波长为0.056m,入射角在19.2°~26.7°之间,带宽为1.6×107Hz,频率为5.33×109Hz,星下点距离为242~347km。由临界垂直基线计算公式(公式(4))得知,在坡度为0°的时候,B⊥C有最大值,约为586m。因此,文中SAR图像间垂直基线临界值取为586m。
3.2.2 多普勒质心频率差
在本文中,试验数据多普勒质心频率差的最大值为56.3Hz。因此项差值较小,在文中取质心频率差最大值临界值,fC为56.3Hz。
3.2.3 时间基线距
永久散射体干涉测量可以监测长时间序列下的地表形变,对于图像间时间基线几乎没有上限限制,这就意味着无需设定时间基线临界值。时间基线对主图像选取的影响只需使用式(7)进行计算即可。
确定影响因子临界值后,分别把13幅Envisat ASAR图像作为主图像,依据其数据参数经函数计算后,生成一个相干系数的对称矩阵(表2),图像编号所代表相应日期的SAR图像(表1),矩阵中行列值为主从图像间的相干系数,零值表示图像完全不相干。根据前述最佳主图像选取规则,依次计算矩阵中13幅主图像对应的零值图像的个数(D)及相干系数均值,并按照D值进行排序。如果出现完全不相干图像数目相同的情况,则依据相干系数均值进行排序。
表2 各图像依次为主图像时的相干系数矩阵Tab.2 Coherence matrix of different master images
续表2
图3为依据最佳主图像选取方法所获取的完全不相干图像个数和相干系数均值的预测结果。图中横轴为图像编号,左侧纵轴为以对应编号为主图像时,主从图像完全不相干个数D,用柱状图表示;右侧纵轴为主从图像的相干系数均值R,用曲线表示。以编号为7(2005-12-01)的影像为主图像时,完全不相干数目为最大值9。以编号为1(2002-12-12)、11(2007-01-25)、12(2007-03-01)的影像为主图像时,它们的完全不相干图像数目为最小值3,对应的相干系数R分别为0.173 362、0.267 507、0.322 989。
图3 相干系数均值和不相干图像数目预测结果Fig.3 Predictions of average coherence and the number of incoherent images
根据前述最佳主图像选取规则,最佳主图像的选取要满足完全不相干图像个数最小且相关系数最大的条件,因此通过该方法选取的最佳主图像是编号为11(2007-01-25)的SAR图像。
为验证最佳主图像选取方法的有效性,分别以试验数据中13幅SAR图像作为主图像,借助StaMPS/MTI软件生成13组干涉图组合。图4~图7分 别 是 编 号 为1(2002-12-12)、11(2007-01-25)、12(2007-03-01)、编号为7(2005-12-01)的影像为主图像时,经软件处理后的主从图像干涉结果。图4~图6干涉效果在这13组干涉图组合中效果较好,都只有3幅干涉图像无明显干涉条纹;图7干涉效果最差,有9幅干涉图像无明显干涉条纹。除了图4~图7外,还有另外9幅SAR影像做主图像时生成的干涉图组合,因篇幅所限,就不在此逐一列出。这些图像的干涉结果和和图3中使用最佳主图像选取方法的预测结果高度吻合。
图4 2002-12-12图像为主图像时干涉效果Fig.4 Interferograms with respect to the master image Dec 12,2002
图5 2007-01-25图像为主图像时干涉效果Fig.5 Interferograms with respect to the master image Jan 25,2007
图6 2007-03-01图像为主图像时干涉效果Fig.6 Interferograms with respect to the master image Mar 1,2007
图7 2005-12-01图像为主图像时干涉效果Fig.7 Interferograms with respect to the master image Dec 1,2005
由软件中干涉结果输出文件计算可得,图4~图6所对应的相干图像的相干系数均值分别为0.48、0.58和0.52,这也和图3中预测的相干系数之间的关系一样:编号为11时相干系数最大,编号为12时次之,编号为1时最小。所以软件处理结果与预测结果相同,均以2007-01-25的SAR影像作为主图像时干涉效果最佳。
从软件实际处理结果和预测结果对比可知,本文提出的最佳主图像选取方法不仅能预测到主从图像完全不相干图像的数目,并且还能够在完全不相干图像数目相同的情况下根据相干系数进一步区分图像干涉效果的好坏。因此该方法在主图像选取中有很好的效果,可以快速、方便地挑选出最佳主图像。
最佳主图像选取是永久散射体干涉测量处理过程中首要考虑的问题。根据SAR图像干涉原理和永久散射体干涉测量中SAR图像的处理流程,本文以系统热噪声、多普勒质心频率以及时空基线作为影响主图像选取因子,在此基础上融合季节性变化的影响,建立PSInSAR最佳主图像选取函数,生成相干系数矩阵,并按照主图像选取规则,预测最佳主图像。通过对比验证,发现使用该方法得到的预测结果与软件处理结果存在着很好的一致性。
同前人未考虑主从图像干涉图最优性或人为设定阈值的方法相比,本文所提出的最佳主图像选取方法不仅可以有效地预测干涉结果中的完全不相干图像个数,并且能够计算出可干涉图像之间的相关性。因此该主图像选取方法可有效避免在永久散射体干涉处理过程中由于主图像选取不当而导致的重复处理SAR数据的现象,有助于提高长时间序列下永久散射体干涉测量的图像处理效率。
[1] FERRETTI A,PRATI C,ROCCA F.Permanent Scatterers in SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2001,39(1):8-20.
[2] FERRETTI A.Submillimeter Accuracy of InSAR Time Series:Experimental Validation[J].IEEE Transactions Geoscience and Remote Sensing,2007,45(5):1142-1153.
[3] ZHU Yefei,YU Jun,WU Jianqiang,et al.Application Comparison of D-InSAR and PS-InSAR Techniques in Suzhou Land Subsidence Monitoring[J].Journal of Geology,2010,34(3):289-294.(朱叶飞,于 军,武健强,等.D-InSAR与PS-InSAR技术应用于苏州地面沉降监测之比较[J].地质学刊,2010,34(3):289-294.)
[4] ZHANG Jingfa,GONG Lixia,JIANG Wenliang.Application of PS InSAR Technique to Measurement of Long Time Crustal Deformation[J].Recent Developments in World Seismology,2006(6):1-6.(张景发,龚利霞,姜文亮.PS InSAR技术在地壳长期缓慢形变监测中的应用[J].国际地震动态,2006(6):1-6.)
[5] GE Daqing,WANG Yan,GUO Xiaofang,et al.Surface Deformation Monitoring with Multi-Baseline D-InSAR Based on Coherent Point Target[J].Journal of Reomte Sensing,2007,11(4):574-580.(葛大庆,王艳,郭小方,等.基于相干点目标的多基线D-InSAR技术与地表形变监测[J].遥感学报,2007,11(4):574-580.)
[6] GONG Wenyu,ZHANG Jixian,ZHANG Yonghong.Detecting Ground Deformation with Permanent Scatterer of Suzhou Region[C]∥International Workshop on Earth Observation and Remote Sensing Applications.Beijing:IEEE,2008:1-5.
[7] JIANG Liming,LIN Hui.Integrated Analysis of SAR Interferometric and Geological Data for Investigating Longterm Reclamation Settlement of Chek Lap Kok Airport,Hong Kong[J].Engineering Geology,2010,110(3):77-92.
[8] ZEBKER H A,VILLASENOR J.Decorrelation in Interferometric Radar Echoes[J].IEEE Transactions on Geoscience and Remote Sensing,1992,30(5):950-959.
[9] CHEN Qiang,DING Xiaoli,LIU Guoxiang.Method for Optimum Selection of Common Master Acquisiton for PSDInSAR[J].Acta Geodaetica et Cartographica Sinica,2007,36(4):395-399.(陈强,丁晓利,刘国祥.PSDInSAR公共主影像的优化选取[J].测绘学报,2007,36(4):395-399.)
[10] ZHU Zhengwei,ZHOU Jianjiang.Optimum Selection of Common Master Image for Ground Deformation Monitoring Based on PS-DInSAR Technique[J].Journal of Systems Engineering and Electronics,2009,20(6):1213-1220.
[11] LONG Sichun,LIU Jingnan,LI Tao,et al.Method for Optimum Selection of Common Master Acquisition for PSDInSAR Fusing GPS Data[J].Journal of Tongji University:Natural Science,2010,38(3):453-458.(龙四春,刘经南,李陶,等.融合GPS数据的PS-DInSAR公用主影像的优化选取[J].同济大学学报:自然科学版,2010,38(3):453-458.)
[12] PAN Bin,SHU Ning.Cluster Analysis for Selection of Time Series Interferometric SAR Imagery[J].Journal of Applied Sciences,2010,28(5):501-506.(潘斌,舒宁.聚类分析用于序列SAR干涉像对选取[J].应用科学学报,2010,28(5):501-506.)
[13] WANG Chao,ZHANG Hong,LIU Zhi.Satellite Synthetic Aperture Radar Inteffefometry[M].Beijing:Science Press,2002.(王超,张红,刘智.星载合成孔径雷达干涉测量[M].北京:科学出版社,2002.)
[14] SHI Shiping.DEM Generation Using ERS-1/2Interferometric SAR Data[J].Acta Geodaetica et Cartographica Sinica,2000,11(4):317-322.(史世平.使用 ERS-1/2干涉测量SAR数据生成DEM[J].测绘学报,2000,11(4):317-322.)
[15] TAO Qiuxiang ,LIU Guolin,HAO Huadong et al.A Study on the Method to Optimize and Select Common Master Images in PS InSAR[J].Journal of Shandong University of Science and Technology,2009,28(5):10-15.(陶秋香,刘国林,郝华东,等.PS InSAR公共主图像优化选取方法研究[J].山东科技大学学报,2009,28(5):10-15.)
[16] ZHANG Hong,WANG Chao,WU Tao,et al.Research on DInSAR Methods Based on Coherent Target[M].Beijing:Science Press,2002.(张红,王超,吴涛,等.基于相干目标的DInSAR方法研究[M].北京:科学出版社,2009.)
[17] BELL J W,AMELUNG F,FERRETIA A,et al.Permanent Scatterer InSAR Reveals Seasonal and Long-term Aquifersystem Response to Groundwater Pumping and Artificial Recharge[J].Water Resources Research,2008,44:1-18.
[18] HERRERA G,TOMáS R,LOPEZ-SANCHEZ J M,et al.Advanced DInSAR Analysis on Mining Areas:La Union Case Study(Murcia,SE Spain)[J].Engineering Geology,2007,90(3-4):148-159.
[19] COLESANTI C,FERRETTI A,PRATI C,et al.Monitoring Landslides and Tectonic Motions with the Permanent Scatterers Technique[J].Engineering Geology,2003,68(1-2):3-14.
[20] HOOPER A,SEGALL P,ZEBKER H.Persistent Scatterer Interferometric Synthetic Aperture Radar for Crustal Deformation Analysis,with Application to Volcan Alcedo,Galapagos[J].Journal of Geophysical Research:B:Solid earth,2007,112(b7):1-21.