地球物理学报  2016, Vol. 59 Issue (4): 1371-1382   PDF    
基于子带干涉测量技术的巴基斯坦地震形变获取研究
庾露1,2, 单新建1, 宋小刚1, 屈春燕1    
1. 中国地震局地质研究所, 地震动力学国家重点实验室, 北京 100029;
2. 广西壮族自治区遥感中心, 南宁530023
摘要: 2013年9月24日发生在巴基斯坦俾路支省(Balochistan)境内阿瓦兰县(Awaran)的MW7.7级地震,在地表产生了最大达10 m的滑动量.利用TerraSAR-X短波雷达数据获取的InSAR同震形变场产生了密集且大范围的干涉条纹,给后续的相位解缠带来困难.而子带干涉法是一种无需或只需进行少量相位解缠,即可获得绝对相位差的新方法.其主要思路是通过缩减带宽以增长波长,从而减少干涉条纹数,降低解缠难度或不需解缠直接得到绝对相位差.但由于带宽的缩减,导致噪声的增大和旁瓣带来的额外干扰,使干涉图质量下降,因此在子带干涉参数选取、噪声滤波以及处理流程等方面需要特殊处理,特别是子带的中心频率和带宽的选取会很大程度地影响测量精度.首先选取典型DEM实验区,以干涉图相干性和误差为评价指标,利用逐步参数选取法,研究相关参数对子带干涉测量的影响,制定最优的参数方案,认识参数选取的原则和方法.在此基础上,将子带干涉应用于巴基斯坦地震的同震形变场获取.最后,将子带干涉、Landsat 8光学影像的交叉频谱相关法、offset-tracking、常规DInSAR获取的同震形变场进行比较,并与模型拟合的形变场进行对比分析.结果表明,子带干涉虽然会受失相干的影响,其提取的形变场范围相较于Landsat 8和offset-tracking有所缺失,但在共同覆盖的区域其精度和噪声水平更优,相比较于常规DInSAR,更适用于条纹密集和形变量大的地区.
关键词: 子带干涉     巴基斯坦地震     同震形变    
Deformation of the 2013 Pakistan MW7.7 earthquake derived from sub-band InSAR
YU Lu1,2, SHAN Xin-Jian1, SONG Xiao-Gang1, QU Chun-Yan1    
1. State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China;
2. Guangxi Remote Sensing Center, Nanning 530023, China
Abstract: On September 24, 2013, an MW7.7 earthquake strok Awaran, Balochistan Province, Pakistan, causing large strike slip up to 10 m. The coseismic deformation field, produced by TerraSAR-X short wave radar data, shows dense and extensive interference fringes which bring difficulties to thephase unwrapping. The sub-band InSAR is a new method for obtaining absolute phase, with a little or without phase unwrapping. This method shortens the bandwidth to increase the wavelength, and could reduce the number of fringes, which decreases the difficulty of phase unwrapping. However the noise and the extra interference of sidelobes will increase and the interferogram quality will decline, while the bandwidth reduces. More consideration should be taken into the parameter selection, noise filtering and processing flow, especially the selection of center frequency and bandwidth. In this paper, we first select the typical DEM experiment zone as the research area, and study the effect of parameter selection on sub-band InSAR with coherences and elevation errors as the evaluation index. Then, we apply the optimal parameter scheme to derive the coseismic deformation field of Pakistan. The coseismic deformation, extracted by sub-band interference, Landsat 8 optical images cross-spectrum correlation, offset-tracking, and conventional DInSAR separately, are analyzed and compared with the deformation field from modeling. It indicates that sub-band InSAR performs batter in precision and noise level, in spite of the slight influence of decorrelation and loss of the extracted deformation field, compared with Landsat-8 and offset tracking. Therefore this method is more suitable to be applied in areas with large deformation and dense fringes.
Key words: Sub-band InSAR     Pakistan earthquake     Coseismic deformation    
1 引言

常规合成孔径雷达干涉测量技术(Interferometry Synthetic Aperture Radar,InSAR)在地形测绘(杜亚男等,2015占文俊等,2015)和形变监测方面(刘云华等,2014;李永生等,2015单新建等,2015)的发展已经非常成熟,其成功应用的关键取决于数据处理中的一些关键步骤,相位解缠便是其中之一,解缠结果的质量直接影响测量的精度(廖明生和林珲,2003).针对解缠问题,众多学者提出了不同的解缠算法,例如:Costantini(1998)提出的最小费用流(minimum cost flow,MCF)法;Goldstein等(1998)提出的枝切法.这些解缠算法的提出,极大地推动了InSAR测量技术的发展.但在一些地表环境复杂地区,大部分的解缠方法存在一定的局限性.例如:一方面由于解缠算法的不同,会造成InSAR的结果各异,从而带来了不确定性;另一方面,对于地形复杂或相位梯度过大的区域,干涉条纹会过于密集,再加上去相干因素的影响,滤波会出现条纹混叠现象,导致后续错误的解缠结果.如何在应用中解决这些问题,已成为当今研究的热点问题之一,在这一需求下,子带干涉成为全球关注的焦点.子带干涉可以增加有效波长,大大减少条纹数量,从而降低解缠难度或直接得到绝对相位.该方法最早由Madsen等(1993)针对当时NASA DC-8空基SAR数据提取地形提出,其基本原理是利用带通滤波器,在距离向的频率域上将主副两景单视复数影像(single look complex,SLC)分别分割成上下两个子带,再通过上-上子带、下-下子带干涉处理后,再次进行上下子带二次干涉,从而模拟获得较长的波长.在其后的一段时间内,受制于天基雷达波带宽的限制,这一技术一直未能得到广泛的应用.随着新一代高分辨率,宽带宽的雷达卫星如TerraSAR、 CosmoSkyMed的相继发射并投入使用,雷达波带宽已达到300 MHz,已可以提取出信噪比很高的子带.在这一背景下,又有一批学者重新对该方法进行了研究.如Derauw等(2010)分析了子带干涉法下地面散射体的相干性;Brcic等(2008)通过子带干涉法利用TerraSAR-X数据,提取了南美Salar de Uyuni盐湖区域的DEM,并与常规InSAR的结果进行了对比.目前对于子带干涉技术的应用,主要集中在提取DEM数据上,而对于地表形变的提取则少有涉及.

而对于本文所研究的2013年发生在巴基斯坦境内的MW7.7级地震的同震地表形变场提取.目前使用遥感手段,主要是利用Landsat 8光学卫星震前和震后影像,基于交叉频谱相关(cross-spectrum correlation)算法(Leprince et al.,2007a)获得(冯光财等,2015Avouac et al.,2014).雷达影像方面,主要使用了TerraSAR-X StripMap模式和ScanSAR模式像对提取同震形变场.其中StripMap模式像对覆盖了形变中心,但是因为条纹分布密集和较大范围失相干等原因,无法通过常规DInSAR提取出与真实情况相符的形变场,一些学者通常采用offset-tracking等方法获得(Jolivet et al.,2014).另一组TerraSAR-X ScanSAR模式的像对其覆盖范围远离形变中心区域,形变量的显著下降,Yague-Martinez等(2014)采取常规DInSAR法成功提取了远离形变中心东北部的同震形变场,并与Avouac等(2014)的同震模型进行了对比,二者结果吻合较高.

本文在前期研究的成果上,首先通过选取典型DEM示范区,研究和讨论子带干涉处理过程中参数的选取问题,并尝试总结出最优化的参数方案.然后在此基础上,将子带干涉法应用于2013年巴基斯坦地震形变中心附近(StripMap模式像对所覆盖的范围)同震形变场的获取,力求通过该方法,更准确地提取干涉条纹密集区的形变量.

2 子带干涉测量法 2.1 基本原理

根据雷达干涉测量原理(李德仁等,2000),干涉相位φ

其中fc是雷达传感器的中心频率,ΔR是斜距差,c是光速.

为了计算高程,首先需要去除平地相位,去平后的干涉相位φ′与高程差Δh之间的关系(王超等,2002)为

其中R为主影像斜距,θ为入射角,B为垂直基线距,λ为波长(λ=c/fc).

当式(2)中φ′=2π时的高程差为Hamb=λRsinθ/(2B)称为模糊高.模糊高是指一个去平后的干涉相位周期所能代表的最大高程差.对于相同的高程差,越小的模糊高会产生越多的干涉条纹.在地形起伏剧烈的区域,较小的模糊高会产生过于密集的干涉条纹,进而影响滤波和相位解缠的结果.由式(2)可知,模糊高与频率成反比,与垂直基线距成反比.由于一组干涉像对其基线距是已定的,如需减少条纹数增加模糊高,只能通过增长波长实现.

而对于地表形变的探测,λ/2称为模糊形变周期,是一个去平并去地形相位后的干涉相位周期所能代表的最大形变差.与模糊高类似,模糊形变周期与波长成正比,对于大形变的区域,过小的模糊形变周期也会产生过于密集的干涉条纹,如需减少条纹数,也只有通过增长波长实现.

子带干涉法可通过模拟的手段获得较长的波长.其原理是:

1)首先利用带通滤波器,将主影像和辅影像在距离向频率域上,分割出上下两个互不重叠的子带像对,如图 1所示,其中上子带和下子带的距离向中心频率分别为fup=fc+f0flo=fc-f0fc为原中心频率,f0为偏移量,B为全频带宽,b为子带带宽.

图 1 距离向频谱分割Fig. 1 Range spectrum divided into upper sub-band and lower sub-band

2)对上子带像对和下子带像对分别进行干涉处理,根据式(1)相应地可得到上子带干涉相位φup和下子带干涉相位φlo,公式为

3)再将上子带干涉相位与下子带干涉相位进行差分,得到子带干涉相位Δφ,公式为

图 1可知,f0的取值范围在B以内,按照一般SAR传感器的设计B<<fc,因此f0<<fc.根据式(5)以及波长和频率的关系(λ=c/f)可知,经过子带干涉处理后,极大地增长了波长,从而在高程计算或地表形变提取中极大地增加模糊高或模糊形变周期.相应地就能减少甚至消除干涉图中的条纹数量,使得相位解缠变得简单甚至无需解缠而得到绝对相位差.

2.2 最优参数选取

选取用于提取DEM数据的影像对为TerraSAR-X Spotlight模式,覆盖范围为澳大利亚中部的北领地南部,影像覆盖澳大利亚乌鲁汝-卡塔楚塔国家公园(Uluru-Kata Tjuta National Park)中心和周边地区,影像中心是一巨型的独块砂岩岩体,又被称为艾尔斯岩(Young et al.,2002).ASTER GDEM 30 m分辨率的高程数据显示该地区最大高程差约为357 m,该地区除艾尔斯岩体外,周边是较为平坦的戈壁,像对获取时间间隔仅11天,可认为期间不存在地表形变,干涉相位主要由地形相位和平地相位组成,影像相干性非常好适合做处理流程和参数选取方面的讨论.因此,采用逐步参数筛选方法,在DEM提取过程中,针对处理流程和参数的选取进行讨论,得出最优化的方案.

2.2.1 处理流程

图 2给出了子带干涉法的基本处理流程:

图 2 子带干涉法处理流程Fig. 2 Processing flow of sub-band InSAR

1)距离向频谱分割:首先将主辅影像分别在距离向频率域上分割出互不重叠的上下子带,如果子带之间频率出现重叠,则重叠部分的能量会对结果产生严重的干扰.在带通滤波器的选取上,为了集中主瓣能量,并能够控制旁瓣能量快速衰减,本文采用了基于Kaiser窗函数的FIR滤波器(Kaiser,1974).

2)子带配准:重采样后的辅影像其频谱会产生一定的改变,如果首先进行配准再进行距离向频谱分割,则有可能在某个频段上引入不必要的误差,表现在干涉图背景上出现了一些周期性的线性条纹.因此,子带配准应在距离向频率分割之后进行.

3)子带干涉:分别对上子带和下子带进行干涉和去除地形相位等处理,得到上下子带去平干涉相位后,再进行一次干涉得到子带干涉相位.

4)相位平滑降噪:由于子带带宽只是全带宽的一部分,而带宽与分辨率成正比(Cumming and Wong,2005),因此子带干涉后的结果会降低分辨率,引入额外的噪声并放大原有的噪声,需要对其进行一定程度的平滑降噪处理.本文采用了基于Goldstein的自适应滤波器(Goldstein and Werner,1998),采用两次迭代以进一步减小噪声.

5)相位解缠:对于一些高程差或形变量很大的区域,子带干涉所产生的模糊高或模糊形变周期有可能仍无法完全涵盖最大高程差或最大形变量并形成相位缠绕,因此仍需进行解缠以求出绝对相位差.反之,如果模糊高或模糊形变周期已涵盖最大高程差或最大形变量,子带干涉相位差就是绝对值,这一步则可以省略.

2.2.2 参数选取

子带的中心频率会影响波长,子带带宽的划分则会影响精度,由于已有ASTER GDEM等数字高程模型作为先验参考,因此可设计一系列配置方案,研究这两个参数对高程测量结果的影响程度,选取线性误差和统计误差作为指标来对不同参数产生的结果进行评价,从而实现最优化参数的选取.如表 1图 3所示,本文以0.05B为间隔,从0.1B开始逐步增大子带带宽,fc为原SLC的中心频率.同时,为了降低谱间干扰的影响,所有方案如图 3所示,均将上下子带带宽的间隔设置最大,子带中心频率均处于子带带宽中心.

表 1 子带中心频率和带宽参数配置方案 Table 1 The parameter setting scheme of central frequency and bandwidth of sub-band

图 3 子带中心频率和带宽配置示意图Fig. 3 Sketch of central frequency and bandwidth of sub-band

线性误差对比分析:首先将上述10组方案的干涉图进行平滑降噪,然后计算得到高程,并与ASTER GDEM高程数据做差分,得到高程残差图(图 4).可看出从方案1至方案10,在高程残差图上出现线性趋势条纹的程度越为严重.分析原因,是因为从方案1至方案10的子带带宽是不断增大的,上下子带带宽的间隔则逐步缩小.带通滤波器自身存在不可避免的旁瓣能量泄漏,会随着间隔的逐步缩小对子带干涉相位产生越为严重的谱间干扰,即出现距离向的线性趋势条纹.而将这些相位条纹转换为高程后,也成为高程上的误差,反映在残差图上就是呈线性趋势性的条纹误差.

图 4 方案1至方案10的高程残差对比Fig. 4 Comparison of elevation residuals among scheme 1 to scheme 10

统计误差对比分析:我们采用了均方根误差(root mean square error,RMSE)来检测方案1至方案10计算得到的高程值相对于ASTER GDEM高程值的偏离程度.由方案1至方案10计算得到的RMSE分别为{47.5844,29.6204,23.0572,24.3553,23.3479,25.7188,23.8768,28.9921,32.4175,37.7373}(m),与各自对应的子带带宽可绘制出如图 5所示的关系曲线.如该曲线所示,RMSE随子带带宽的增大呈现先下降后上升的趋势.其中最小值出现在0.2B的位置,对应方案3.而在此之前方案1和方案2的RMSE均较大,是因为其对应的子带带宽0.1B和0.15B过小,分辨率损失过大所带来的相位模糊造成的,而旁瓣能量泄漏此时造成的影响在此处非常轻微,基本可以忽略.当子带带宽大于0.2B且继续增大时,对应的RMSE在0.25B~0.35B处出现缓慢抬升,而当>0.35B时则出现较快抬升.分析原因,虽然子带带宽的不断增大会对分辨率以及干涉质量(相干性)有所提升,但是也会使谱间干扰问题不断加重,给测量结果带来的负面影响甚至超过了分辨率提升带来的增益.

图 5 方案1至方案10的均方根误差对比Fig. 5 Comparison of RMSE among scheme 1 to scheme 10

综合上述分析,由于子带带宽和分辨率成正比,当设置较小带宽以避免谱间干扰影响的同时,会出现因分辨率降低过大而造成的模糊和噪声放大等误差,反之如果设置了较大的带宽虽然能提高子带影像的分辨率,却会带来谱间干扰的严重影响.因此,可看出谱间干扰的减少和分辨率的提高是“此消彼长”的相互矛盾的关系,无法同时兼顾.通过在真实数据基础上的一系列实验,可得到子带带宽为0.2~0.3B,且上下子带带宽间隔最大时,分辨率可获得较可靠的精度保证,同时受谱间干扰的影响也较为轻微,是理想的子带带宽选择区间.

3 巴基斯坦地震子带干涉形变场获取 3.1 震区概况与数据选取

2013年9月24日,在巴基斯坦俾路支省(Balochistan)境内的阿瓦兰县(Awaran)北北东66 km处发生了一次震级为MW7.7的地震,时隔4天后,又一个震级为MW6.6级的地震在该地区发生.冯光财等(2015)Avouac等(2014),均采用交叉频谱相关(cross-spectrum correlation)算法(Leprince et al.,2007b),利用Landsat 8光学卫星震前和震后影像,计算像对的偏移量获取了此次地震大范围的地表形变场.这些研究显示,此次地震形变沿着地表破裂迹线呈弧形分布,地震类型为左旋走滑兼少许逆冲,破裂迹线北部造成的形变总体上大于南部.本文选取的像对为TerraSAR-X StripMap模式,影像覆盖此次地震的形变中心,参数如表 2所示.ASTER GDEM资料显示,该区域在破裂迹线北部山岭分布密集,多数山峰在800~1000 m以上,地形起伏情况相较于南部更为复杂.图 6b是二轨法差分干涉相位与强度图的叠加,由于X波段的模糊形变周期只有1.6 cm,对形变和地形都十分敏感,易产生密集的干涉条纹或失相干.干涉条纹在分布形态上,破裂迹线南北两侧各有不同.南侧以大量密集且不连续的干涉条纹为主(图 7-2);在影像中部迹线南侧位置,存在一处大面积的失相干带,由于此处是居民地,可能受地震影响建筑物大面积倒塌,或是人类活动改变了地表地面状况而造成失相干.迹线北侧则因过大的形变和复杂的地形影响造成了大范围的失相干,只有影像西北角残存了部分干涉条纹(图 7-1).

图 6 (a) Landsat 8提取的2013年巴基斯坦地震地表形变(冯光财等(2015)),已转换为雷达地距方向;正值为接近卫星方向; (b) TerraSAR-X StripMap常规DInSAR干涉相位与强度图叠加;其覆盖区域为图6a中绿色矩 形框所示.图6a和图6b中的黑色实线为地表破裂迹线Fig. 6 (a) The surface deformation of 2013 Pakistan earthquake from Landsat 8 satellite images (Feng, et al, 2015), which is transformed into the ground-range direction. Positive values mean close to the satellite; (b) The stacking chart of conventional DinSAR interferometric phase and magnitude of TerraSAR-X StripMap images. The green rectangular of Fig.6a indicates its coverage region. The black solid lines in Fig.6a and Fig.6b are surface rupture zone

图 7 图6b的局部放大,其覆盖区域为.图 6b中红色矩形框所示Fig. 7 The partial enlargement of Fig.6b. The red rectangular of Fig.6b indicates its coverage region

表 2 阿瓦兰地区TerraSAR-X StripMap影像参数 Table 2 The parameters of TerraSAR-X Stripmap images
3.2 子带干涉形变场提取与验证

由Landsat 8提取的形变场显示(图 6a),本文选用的TerraSAR-X影像覆盖范围内的地距向形变量为-5.64~2.88 m,对应的斜距向形变约为-4.69~2.4 m.使用子带干涉法需要模拟产生至少18.76 m的波长,即子带中心频率约为0.16B=16 MHz,才不会产生相位缠绕.由3.3.1节分析可知,这样的参数配置会使上下子带的中心频率过于接近,带来谱间干扰的影响,从而使结果存在较为严重的误差.另一方面,常规二轨法差分干涉相位显示(图 6b),紧邻破裂迹线区域由于形变量过大和地形复杂等因素,导致失相干.这些在常规InSAR中失相干区域的区域,其干涉条纹不可识别,干涉相位被严重的噪声淹没,即使经过子带干涉处理,其结果也是不正确的,而实际有干涉条纹覆盖的区域,其斜距向形变量约为-4.0~2.4 m.综合上述分析,由于此次地震产生的地表形变量远大于一般子带干涉所能模拟的模糊形变周期,因此在子带中心频率和带宽的配置上,应以减少干涉条纹数为主要目标.

通过上文的一系列实验,已经获得了子带参数配置的最优方案.该方案虽然基于地形提取实验,但实验中讨论的分辨率缩减和谱间干扰这两方面问题实际是由图像自身性质变化造成的,与试验区的选取以及子带干涉的用途无关.换句话说,子带干涉无论用途为何,其最终目的都是为了减少干涉条纹数以降低解缠难度,因此上文的最优参数配置方案,也同样适用于子带干涉形变场的提取.综合上述分析,本文将采用上文中子带中心频率为0.8B,子带带宽为0.2B的最优方案,模拟产生的波长为3.75 m,模糊形变周期为1.87 m,提取该震例的形变场.

图 8a是经过子带干涉处理、5(距离向)×6(方位向)的多视处理,以及3次迭代自适应滤波处理(滤波窗口分别为128、96、64)后得到的干涉相位,其干涉条纹分布在破裂迹线南侧十分稀疏,相位在大部分地区处在一个形变周期范围内,只在靠近破裂迹线附近出现了少量的缠绕现象;北侧的干涉条纹出现在影像西北角,并出现了明显的相位缠绕.参考Landsat 8提取的形变场(图 6a),破裂迹线南北两侧的形变量由内向外逐步递减,且两侧的形变运动方向相反,因此需要对破裂迹线南北两侧的相位分别解缠.图 8b是分别对断层两侧子带干涉相位采用MCF方法解缠,计算地距向形变,并拼接得到的结果,红色圆点所在位置为解缠起始点.

图 8 (a) 子带干涉相位; (b) 分别对破裂迹线南北两侧的干涉相位解缠,并计算得到的地距向形变场;红色圆点为解缠起始点.图8a和图8b中的黑色实线为地表破裂迹线Fig. 8 (a) Sub-band interferometric phase; (b) Phase unwrapping from north and south side of surface rupture zone respectively, and obtaining deformation field in slant direction. Red point is the beginning of phase unwrapping. Note: The black solid lines in Fig.8a and Fig.8b are surface rupture zone
4 形变场比较与验证

由于阿瓦兰地震发生在偏远的山区,GPS台站分布稀疏,地面调查资料缺乏,因此本文使用的基准参考形变场来自冯光财等(2015)的研究成果,其主要基于Landsat 8、USGS(美国地质调查局)的观测结果,经Okada弹性半空间形变模型(Okada,1985)反演出此次地震的断层滑动分布,最后正演得到(图 9a).本文将这一结果分别与子带干涉法、Landsat 8光学影像的交叉频谱相关法、offset-tracking法和常规DInSAR各自提取的地表形变场(图 9b-e)进行差分,得到各自的残差图(图 9f-i),并分别在图 9a-e的相同位置从南至北绘制了3条横跨断层的剖面p1,p2,p3(图 10).

图 9 不同手段获得的同震形变场及残差分布图
(a) 模型拟合得到的地表形变场; (b) 子带干涉法提取的地表形变场,(低相干区已被掩膜去除); (c) Landsat 8提取的地表形变场; (d) TerraSAR-X基于offset-tracking算法提取的地表形变场; (e) 常规DInSAR提取的地表形变场,(低相干区已被掩膜去除); (f-i) 子带 干涉法、Landsat 8、offset-tracking法和常规DInSAR各自的的残差分布.
Fig. 9 The coseismic deformation field and residuals distribution map obtained by different approaches.
(a) Model prediction, (b) sub-band InSAR (low coherence area has been removed by mask), (c) Landsat 8, (d) TerraSAR-X based on offset-tracking, (e) conventional DInSAR (low coherence area has been removed by mask); (f-i) The residuals distribution of sub-band interferometric, Landsat 8, offset-tracking and conventional DInSAR.

图 10 剖面地距向形变量Fig. 10 The deformation value in the profile ground-range direction

对比分析图 9图 10,可看出以上4种方法均不同程度地获得了此次地震的地表形变场.其中:

1)offset-tracking法和交叉频谱相关法,都是一种基于强度信息的子像素精密匹配算法,适用于大形变量的探测(Gray et al.,1998),offset-tracking法还可避免InSAR测量中低相干或失相干而造成无法解缠等问题(刘云华等,2012李佳等,2013).如图 9c图 9d所示,上述两种方法相比子带干涉法(图 9b)和常规DInSAR(图 9e),均获得了更完整的地表形变场和更清晰的地表破裂迹线,且迹线南北两侧的运动方向也与模型拟合的结果(图 9a)相吻合.但这两种方法会受到地形、阴影和数据质量等因素的影响,其精度一般要比InSAR低,offset-tracking法的理论精度在分米级(Strozzi et al.,2002),而交叉频谱相干法的理论最高精度为像元大小的1/50(Leprince et al.,2007a),另外二者均存在细节分辨力偏低和噪声偏大等问题.图 9c图 9d显示,在破裂迹线南侧,存在大量的噪声异常值;与之对应的残差分布(图 9g图 9h)上,可看到北侧在紧邻迹线的区域形变值偏大;在远离破裂迹线的北侧区域,Landsat 8表现出总体吻合局部偏大的特点,offset-tracking则呈现由近到远先偏大后偏小的趋势.从剖面中可看出,上述2种方法均存在形变值在各处均呈现波动较大的现象,说明受噪声干扰较为严重.

2)子带干涉法和常规DInSAR法,受到影像中2条失相干带的影响,无法计算出正确的相位,进而在最终结果中(图 9b图 9e)该处的形变值为幅值抖动剧烈的噪声,故本文已将其掩膜去除.如图 9f图 9i图 10所示,2条失相干带将影像分为南部、中部和西北部3块相对独立的区域(分别在图 9b图 9e中以I、II和III标出),其中:

I号区域,位于影像南部是远离形变中心的区域,上述2种方法在此处获得的形变场均与模型拟合结果拟合较好,其平均误差分别在5 cm和8 cm的范围内;剖面线(图 10)显示与模型拟合结果相吻合,且噪声相较于offset-tracking法和交叉频谱相关法少,形变趋势变化平稳.

II号区域,在影像中部,破裂迹线的南侧,是形变剧烈的区域.在图 9b-d中黑色实线矩形框标出了一处形变值急剧增大区域,对应图 10b剖面线可看出相较于模型拟合结果,形变值发生了急剧的抬升,在残差图上,此处出现了较大的负值误差,而常规DInSAR的结果(图 9e)在此处与模型拟合结果吻合较好.本文认为,上述出现较大误差的3种方法的结果,具有一定的一致性,且是由不同数据源和不同处理手段得到的,因此这一形变值剧烈增大的区域可能客观存在.如果这一假设成立,则模型拟合结果此处形变值变化平缓,可能与四叉树降采样的密度设置有关,而DInSAR的结果也未出现形变值剧增的趋势,则可能是因为对密集条纹做解缠时,因噪声和低相干性的影响,造成条纹混叠并丢失导致.

III号区域,位于影像西北部,由模型拟合得到的形变场(图 9a)及对应的剖面显示(图 10),该地区形变量约为-1.4~-3.7 m,由南至北基本呈线性递减,递减速率(斜率)约为0.15 m·km-1.Landsat 8和offset-tracking提取的形变场在剖面上,与该结果基本吻合,但存在振幅波动较大的噪声.而常规DInSAR计算得到的形变,与其他方法得到的形变差异很大,这些差异同时反映在形变速率以及形变值分布这2个方面:剖面显示,其形变递减速率十分平缓只有2.6 cm·km-1;形变分布(图 9e)显示,其计算得到的地距向形变值约为-0.15~-0.36 cm,与模型拟合结果的误差很大(图 9i).造成上述2方面误差的原因,也与该处干涉条纹密集,解缠时条纹丢失有关.对于子带干涉,在该处只发生了一个周期的相位缠绕,可十分容易地.解缠得到更接近真实情况的绝对相位差.图 9b所示,由子带干涉获取的西北部的地距向形变值约为-0.87~-3.77 m与模型拟合结果相接近.残差分布(图 9f)显示,残差由北至南,由形变远场至近场呈现递减的趋势;剖面(图 10)显示,最远端与模型拟合结果的误差最大达到0.7 m,随着逐步靠近形变近场,误差逐步控制在0.5~0.2 m或更低的水平,相较于offset-tracking和Landsat 8提取的形变值的平均误差(分别为1.5 m和1.8 m)更小,且趋势平稳增大,噪声水平更低.我们注意到,子带干涉法在此处出现了相较于影像南部更大的误差,其原因是因为采用了常规的多项式全局配准策略,在处理影像局部存在的大偏移、复杂地形情况时出现的局部偏差,这一配准处理方法上的局限,在常规InSAR中也是存在且难以避免的.

5 结论

子带干涉法作为一种传统InSAR技术的衍生,其减少干涉条纹数的特性,可有效避免对条纹密集区域进行解缠出现的误差;也可避免因绝对相位差过大而待解缠区域过小,造成条纹丢失,产生错误的解缠结果.因此可直接用于测量,也可作为一种手段用于检验传统InSAR解缠结果的准确性.本文通过子带干涉法以提取澳大利亚艾尔斯岩区域的DEM为基础,分析了子带中心频率和带宽参数的选取,并将最优化方案应用于2013年巴基斯坦地震的同震地表形变场提取中,通过对这实例数据进行处理和对比分析,得出如下结论:

1)子带中心频率和带宽参数配置的不同会对最终结果产生较大的影响.在保证波长足够长的前提下,应选择高相干点分布比例大的配置方案,可保证干涉图的总体质量;应将上下子带的中心频率和带宽的间隔设置尽量大,可减少谱间干扰产生的伪信号混淆测量结果以提高测量精度.

2)为避免子带干涉图上产生斜向条纹影响测量结果,应首先对主辅影像做距离向频率分割,生成上下子带影像,再分别对其进行配准和重采样.在对子带干涉图进行降噪平滑时,可将搜索窗口适当增大,并进行二次或更多次的迭代滤波,以达到更好的降噪效果.

3)同常规InSAR一样,子带干涉法也依赖于相干性.常规InSAR中失相干的区域,在子带干涉中该区域的干涉相位也可能是不准确.因此,可利用常规InSAR的相干性分布图,结合其他观测数据,预先估计出影像中有效干涉相位所覆盖的区域内的值域范围,以选取合适的子带中心频率和带宽等参数.在对常规InSAR干涉条纹密集的区域做解缠时,易发生条纹丢失使参与计算的条纹数低于实际情况,从而造成较大的误差.而子带干涉法能有效减少或消除干涉条纹数,降低解缠难度甚至无需解缠,可有效地避免条纹丢失问题,获得更为准确的结果.

4)通过对2013年巴基斯坦震例进行的研究,可看出子带干涉法在提取大尺度地表形变上具有很大的优势.首先,虽然受失相干的影响其提取的形变场范围相较于交叉频谱相关法和offset-tracking法有所缺失,但在共同覆盖的区域,子带干涉法所测量的结果均显示出了更高的精度和更低的噪声水平.其次,在DInSAR存在干涉条纹且十分密集的区域(如影像中部区域),子带干涉法不但能更为准确地获得绝对相位差,而且与交叉频谱相关法或offset-tracking法相比,其反映的形变场细节更丰富.

致谢 感谢夏耶教授生前对本文完成给予的十分宝贵的帮助.

参考文献
[1] Avouac J P, Ayoub F, Wei S J, et al. 2014. The 2013, MW7.7 Balochistan earthquake, energetic strike-slip reactivation of a thrust fault. Earth and Planetary Science Letters, 391: 128-134.
[2] Brcic R, Eineder M, Bamler R. 2008. Absolute phase estimation from TerraSAR-X acquisitions using wideband interferometry.//Proceedings of IEEE Radar Conference. Pasadena, CA: IEEE.
[3] Costantini M. 1998. A novel phase unwrapping method based on network programming. IEEE Transactions on Geoscience and Remote Sensing, 36(3): 813-821.
[4] Cumming I G, Wong F H. 2005. Digital Processing of Synthetic Aperture Radar Data: Algorithms and Implementation. Boston: Artech House.
[5] Derauw D, Orban A, Barbier C. 2010. Wide band SAR sub-band splitting and inter-band coherence measurements. Remote Sensing Letters, 1(3): 133-140.
[6] Du Y N, Feng G C, Li Z W et al .2015.Generation of high precision DEM from TerraSAR-X/TanDEM-X. Chinese Journal of Geophysics (in Chinese),58(9): 3089-3102,doi: 10.6038/cjg20150907.
[7] Feng G C, Xu B, Shan X J, et al. 2015. Coseismic deformation and source parameters of the 24 September 2013 Awaran, Pakistan MW7.7 Earthquake derived from optical Landsat 8 satellite images. Chinese Journal of Geophysics (in Chinese), 58(5): 1634-1644, doi: 10.6038/cjg20150515.
[8] Goldstein R M, Zebker H A, Werner C L. 1988. Satellite radar interferometry: Two-dimensional phase unwrapping. Radio Science, 23(4): 713-720.
[9] Goldstein R M, Werner C L. 1998. Radar interferogram filtering for geophysical applications. Geophysical Research Letters, 25(21): 4035-4038.
[10] Gray A L, Mattar K E, Vachon P W, et al. 1998. InSAR results from the RADARSAT Antarctic Mapping Mission data: estimation of glacier motion using a simple registration procedure. 1998 IEEE International Geoscience and Remote Sensing Symposium, 1998, 3: 1638-1640.
[11] Jolivet R, Duputel Z, Riel B, et al. 2014. The 2013 MW7.7 Balochistan Earthquake: Seismic potential of an accretionary wedge. Bulletin of the Seismological Society of America, 104(2): 1020-1030.
[12] Kaiser J F. 1974. Nonrecursive digital filter design using the I0-Sinh window function.//Proceedings of the 1974 IEEE International Symposium on Circuits and Systems. New York: IEEE, 20-23.
[13] Leprince S, Barbot S, Ayoub F, et al. 2007a. Automatic and precise orthorectification, coregistration, and subpixel correlation of satellite images, application to ground deformation measurements. IEEE Transactions on Geoscience and Remote Sensing, 45(6): 1529-1558.
[14] Leprince S, Ayoub F, Klingert Y, et al. 2007b. Co-registration of optically sensed images and correlation (COSI-Corr): An operational methodology for ground deformation measurements.//Proceedings of IEEE International Geoscience and Remote Sensing Symposium. Barcelona: IEEE, 1943-1946.
[15] Li D R, Zhou Y Q, Ma H C. 2000. Principles and applications of interferometry SAR. Science of Surveying and Mapping (in Chinese), 25(1): 9-12.
[16] Li J, Li Z W, Wang C C et al .2013.Using SAR offset-tracking approach to estimate surface motion of the South Inylchek Glacier in Tianshan. Chinese Journal of Geophysics,56(4): 1226-1236,doi: 10.6038/cjg20130417.
[17] Li Y S, Feng W P, Zhang J F,et al. 2015.Coseismic slip of the 2014 MW6.1 Napa, California Earthquake revealed by Sentinel-1A InSAR.Chinese Journal of Geophysics (in Chinese),58(7): 2339-2349,doi: 10.6038/cjg20150712.
[18] Liao M S, Lin H. 2003. Synthetic Aperture Radar Interferometry: Principle and Signal Processing (in Chinese). Beijing: Surveying and Mapping Press.
[19] Liu Y H, Qu C Y, Shan X J. 2012.Two-dimensional displacement field of the Wenchuan earthquake inferred from SAR intensity offset-tacking. Chinese Journal of Geophysics (in Chinese),55(10): 3296-3306,doi: 10.6038/j.issn.0001-5733.2012.10.012.
[20] Liu Y H, Wang C S, Shan X J, et al. 2014.Result of SAR differential interferometry for the co-seismic deformation and source parameter of the MS7.0 Lushan Earthquake.Chinese Journal of Geophysics (in Chinese), 57(8): 2495-2506,doi: 10.6038/cjg20140811.
[21] Madsen S N, Zebker H A, Martin J. 1993. Topographic mapping using radar interferometry: Processing techniques. IEEE Transactions on Geoscience and Remote Sensing, 31(1): 246-256.
[22] Okada Y. 1985. Surface deformation due to shear and tensile faults in a half-space. Bulletin of the Seismological Society of America, 75(4): 1135-1154.
[23] Shan X J, Zhang G H, Wang C S, et al. 2015.Joint inversion for the spatial fault slip distribution of the 2015 Nepal MW7.9 earthquake based on InSAR and GPS observations.Chinese Journal of Geophysics (in Chinese),58(11): 4266-4276,doi: 10.6038/cjg20151131.
[24] Strozzi T, Luckman A, Murray T, et al. 2002. Glacier motion estimation using SAR offset-tracking procedures. IEEE Transactions on Geoscience and Remote Sensing, 40(11): 2384-2391.
[25] Wan Y G. 2012. Digital Signal Processing-Using MATLAB (in Chinese). 2nd ed. Beijing: Science Press.
[26] Wang C, Zhang H, Liu Z. 2002. Spaceborne Synthetic Aperture Radar Interferometry (in Chinese). Beijing: Science Press.
[27] Wang X P. 2010. Minimum cost flow and its ameliorative algorithm for phase unwrapping of InSAR image. Science of Surveying and Mapping (in Chinese), 35(4): 129-131.
[28] Yague-Martinez N, Fielding E, Haghshenas-Haghighi M, et al. 2014. Ground displacement measurement of the 2013 M7.7 and M6.8 Balochistan earthquake with TerraSAR-X ScanSAR data.//Proceedings of IEEE International Geoscience and Remote Sensing Symposium. Quebec City, QC: IEEE, 950-953.
[29] Young D N, Duncan N, Camacho A, et al. 2002. Ayers Rock, Northern Territory, Map Sheet GS52-8 (2nd edition) (Map). 1:250 000. Northern Territory Geological Survey. Geological Map Series Explanatory Notes.
[30] Zhan W J, Li Z W, Wei J C, et al. 2015.A strategy for modeling and estimating atmospheric phase of SAR interferogram. Chinese Journal of Geophysics (in Chinese), 58(7): 2320-2329,doi: 10.6038/cjg20150710.
[31] 杜亚男, 冯光财, 李志伟等. 2015.TerraSAR-X/TanDEM-X获取高精度数字高程模型技术研究. 地球物理学报,58(9): 3089-3102,doi: 10.6038/cjg20150907.
[32] 冯光财, 许兵, 单新建等. 2015. 基于Landsat 8光学影像的巴基斯坦Awaran MW7.7地震形变监测及参数反演研究. 地球物理学报, 58(5): 1634-1644, doi: 10.6038/cjg20150515.
[33] 李德仁, 周月琴, 马洪超. 2000. 卫星雷达干涉测量原理与应用. 测绘科学, 25(1): 9-12.
[34] 李佳, 李志伟, 汪长城等. 2013.SAR偏移量跟踪技术估计天山南依内里切克冰川运动. 地球物理学报,56(4): 1226-1236,doi: 10.6038/cjg20130417.
[35] 李永生, 冯万鹏, 张景发等. 2015. 2014年美国加州纳帕MW6.1地震断层参数的 Sentinel-1A InSAR反演. 地球物理学报, 58(7): 2339-2349, doi: 10.6038/cjg20150712.
[36] 廖明生, 林珲. 2003. 雷达干涉测量: 原理与信号处理基础. 北京: 测绘出版社.
[37] 刘云华, 屈春燕, 单新建 .2012.基于SAR影像偏移量获取汶川地震二维形变场. 地球物理学报,55(10): 3296-3306,doi: 10.6038/j.issn.0001-5733.2012.10.012.
[38] 刘云华, 汪驰升, 单新建等. 2014.芦山MS7.0级地震InSAR形变观测及震源参数反演. 地球物理学报,57(8): 2495-2506,doi: 10.6038/cjg20140811.
[39] 单新建,张国宏,汪驰升等. 2015. 基于InSAR和GPS观测数据的尼泊尔地震发震断层特征参数联合反演研究. 地球物理学报, 58(11): 4266-4276, doi: 10.6038/cjg20151131.
[40] 占文俊, 李志伟, 韦建超等. 2015.一种InSAR大气相位建模与估计方法. 地球物理学报,58(7): 2320-2329,doi: 10.6038/cjg20150710.