引言

合成孔径雷达干涉(interferometry synthetic aperture radar，简写为InSAR)技术作为一种获取大面积地形、地表位移的有效手段被广泛应用于地形图测绘、土地利用类型分类、海洋风浪场反演等领域，特别是在地震(Lindsey et al，2015)、火山活动(Gong et al，2015)、冰川漂移(Li et al，2015)、地表沉降(Zhang et al，2016)与山体滑坡(Wang et al，2016)等引起的地表形变监测方面，InSAR显示出其独特的技术优势.

随着InSAR技术研究的不断深入，尽管其观测精度已有较大提高，在实际应用中却由于大气效应等因素的影响，而使得其广泛应用受到了限制. Massonnet等(1994)首次给出了由于大气误差所产生的干涉条纹；而Zebker等(1997)的研究进一步表明，相对湿度在时间和空间上20%的变化，可能会导致形变测量中10—14 cm的测量误差.

当使用InSAR技术研究地震同震或者火山喷发等剧烈的地球物理过程时，大气误差相比地表形变量较小，大气误差的存在对研究其物理机制不会产生很大影响；但是使用InSAR获取震间形变会变得非常具有挑战性，这是由于震间形变微弱，信噪比小，形变信号很容易被噪声淹没.因此，需要对大气误差进行精细的处理，才能保证震间微小形变获取的精度.目前已有多种技术来消除InSAR大气效应，可以划分为两类：其一为基于影像自身的校准技术，例如永久散射体(permanent scatterers)，小基线集(small baseline subset)等时序InSAR技术；其二为基于外部数据的校准技术，例如GPS，MERIS (medium resolution imaging spectrometer)，MODIS (moderate-resolution imaging spectrometer)，WRF (weather research and forecasting)和ERA-Ⅰ (ERA-Ⅰnterim)等数字大气模型.基于影像自身的校准方法一般是基于时间序列影像，利用先验模型加后续滤波的方法进行大气效应去除.但是对震间形变的获取而言，由于形变波长较长，形变信号往往与轨道和大气误差耦合在一起，纯粹的模型改正和滤波方法常常不起作用，因此需要借助外部的大气数据来加以改正. GPS数据可以给出单点上精确的大气延迟量，但由于一般的GPS测站空间密度过低，无法对干涉图所覆盖的大范围区域内的大气延迟信号给出合理的表达，因而无法满足InSAR的大气改正.而一些星载的成像光谱数据(例如MODIS和MERIS)和数字气象模型数据(例如ERA-Ⅰ和WRF)，由于其稳定的数据来源和空间分辨率或对大气条件的良好估计而得到了InSAR研究人员的青睐(Li et al，2006；Doin et al，2009；Jolivet et al，2014；Jung et al，2014).但是由于大气在时空上的多变性，这些数据的通用性比较差，目前还没有一种外部数据对大气延迟的估计精度可以满足所有区域的InSAR大气校正，最佳的方法是基于研究区数据本身来选取最好的大气校正数据.

本文拟选取青藏高原东北缘的老虎山断裂为研究区，基于研究区本身短时间基线干涉数据，对不同外部大气数据(MERIS, ERA-Ⅰ, WRF)改正进行试验研究，来评价不同大气数据的改正效果，从而为研究区选择出最优的大气改正方法，进而用于后续长时间基线、短空间基线干涉图的大气改正，在此基础上利用层叠法(stacking)获取老虎山断裂精确的震间形变速率场，并通过与已有成果的比较来验证最优方法的可靠性.

1.   3种外部数据大气改正方法比较

1.1   3种外部大气数据简介

本文选取星载成像光谱数据(例如MERIS)和数字气象模型数据(例如ERA-Ⅰ和WRF)作为外部大气数据来源参与改正InSAR的大气效应.

MERIS是被动多光谱扫描仪，与ASAR (advanced synthetic aperture radar)同时搭载在欧洲空间局发射的EnviSAT卫星上. MERIS产品的空间分辨率有300 m和1200 m两种，本文选用后者. MERIS数据由于空间分辨率高，且与ASAR数据获取时间同步，因此通常被认为是改正ASAR影像大气效应的最优数据源(Li et al，2006).然而，MERIS是光学传感器，只能在白天工作，并且MERIS对云的敏感性使其应用范围受限.本文所采用的MERIS数据，其云覆盖率均小于30%，并将数据中云覆盖区域的像元设为空值.

ERA-Ⅰ是欧洲中尺度气象预报中心(European Centre for Medium Range Weather Forecasting)发布的全球再分析资料.该数据包括16个气象参数(气温，相对湿度，水汽分压等)的估值，时间为1989年以来每天UTC 00:00，06:00，12:00，18:00，水平分辨率为0.75°，垂直方向0—50 km被划分成37层.本文采用的SAR影像的获取时刻为UTC 03:15，该时刻的大气延迟是利用ERA-Ⅰ数据计算两个邻近时刻(00:00和06:00) 的大气延迟，然后进行线性插值而得(Bekaert et al，2015).

WRF是由美国多个科研机构合作建立的新一代中尺度模式，集数值天气模式、大气模拟与数据同化于一体，能够更好地改善对中尺度天气的模拟和预报.本文使用了WRF 3.4预测模型，由美国国家环境预报中心提供的FNL数据作为初始边界条件，模拟运行时长为12小时，包含影像获取前后各约6小时.本次WRF模拟试验中采用三重嵌套网格，从外到内网格分辨率分别为27 km，9 km和3 km，最终采用分辨率为3 km的模拟结果参与计算；垂直方向分为31层，最高层气压设为10 hPa.考虑到海原地区的地理条件(例如气候特点、植被覆盖等)(王立稳等，2015)，模拟使用的物理方案为：微物理过程采用Lin方案，长波辐射采用ntm方案，短波辐射采用Dudhia方案，陆面过程采用Noah陆面模式，边界层采用MYJ TKE方案，积云对流参数化方案采用Grell-Devenyi方案.

1.2   大气延迟计算方法

SAR传感器发射的电磁波在卫星与地表之间传播时，其相位会受空气折射率的影响，这是大气误差的主要来源.对于C波段的ASAR数据来说，其大气延迟主要来自对流层，按物理来源可将对流层相位延迟分为两部分，即流体静力学延迟L^hydro和湿延迟L^wet.尽管绝对的流体静力学延迟量能达到约2.3 m，较湿延迟(约0.3 m)大得多(Bevis et al，1992)，但是由于InSAR干涉相位是两次影像的差分相位，相比于相对稳定的静力学延迟，在时间和空间上具有更大多变性的湿延迟对InSAR干涉相位的影响更大.

对于ERA-Ⅰ和WRF大气产品来说，我们可以从中直接提取出一些大气参数，然后基于下式计算大气延迟量：

式中，P表示总气压，单位为Pa；T为气温，单位为K；e表示水汽分压；k₁, k′₂, k₃为通用气体常数，其中k₁=0.776 K/Pa, k′₂=0.716 K/Pa, k₃=3.75×10³ K²/Pa.

而对于MERIS产品，其中除了包括来自ERA-Ⅰ中用于计算流体静力学延迟量的气象参数外，还提供了观测的水汽数据(即大气可降水汽含量PWV).假设压力随着海拔的增加成指数减小，标高为8.34 km (Ramon et al，2003)，则大气延迟量为

式中，R_d为气体常数，R_d=287.05 J/(kg·K)；P₀为海平面气压；h为海拔；Π为无量纲转化参数，此处Π=6.2 (Li et al，2006).

通过上述计算方法，我们可以利用3种大气数据计算出每张SAR影像获取时刻对应的大气延迟分布. 图 1给出了其中某一幅影像(10016) 对应的3种外部大气数据计算得到的流体静力学延迟和湿延迟.卫星相位中记录的大气延迟是沿着卫星视线向的大气相位，可由垂直向的天顶相位延迟投影至卫星视线向得到，即L_inc=L/cosθ，其中θ表示卫星入射角.对于一张干涉图的大气相位来说，可由主影像和辅影像获取对应时刻(t₁和t₂)的大气相位相减得到，即φ^atm=－4π(L₁－L₂)/(λcosθ)，其中L₁和L₂分别为t₁和t₂时刻的垂直向天顶相位延迟，λ为电磁波的波长.

图  1  利用3种不同的大气改正方法获取的同一时刻同一区域的流体静力学延迟(a)和湿延迟(b)

Figure  1.  Hydrostatic delay (a) and wet delay (b) for the same area at the same time based on three different atmospheric correction methods (MERIS, ERA-Ⅰ and WRF)

下载: 全尺寸图片 幻灯片

1.3   3种大气改正方法对比研究

为了评价MERIS, ERA-Ⅰ和WRF这3种大气数据的大气改正效果，本文利用研究区中的短时间基线干涉图进行相关的试验研究.短时间基线干涉图中，由于所覆盖的时间间隔较短，一般为1—2个月，而该地区的震间形变信号大约处于mm/a数量级，所以干涉相位中的形变信号几乎可以忽略不计，大气和轨道信号占据主导地位，而轨道信息一般符合线性趋势，通过模型拟合对其加以去除，剩余的信号则可认为完全是由大气延迟引起.所以首先利用3种大气数据计算的延迟量进行大气去除，然后去除轨道信号，最后对残余的相位(即相位残差)进行统计分析，进而评价3种大气改正方法的有效性. 图 2为整个评价的流程图.

图  2  对比3种大气改正方法的流程图

Figure  2.  Process flowchart for comparing three different atmospheric correction methods

下载: 全尺寸图片 幻灯片

本文在研究区雷达数据集中，挑选了7对短时间基线(1—2个月)干涉对，形成差分干涉图，然后利用相应的MERIS, ERA-Ⅰ和WRF数据计算出大气延迟量，分别用于对干涉图进行大气改正，然后进行轨道误差拟合与去除，最后得到残余相位图.

为了定量化地评价MERIS, ERA-Ⅰ和WRF这3种大气改正方法，本文计算了残余相位的频率分布图，如图 3所示.每个直方图表示7幅短时间基线干涉图的残差统计结果，在理想情况下，所有误差在改正后均应为0，即残余相位为0的频率是100%；但是，实际中大气误差不可能被完全消除.使用正态曲线拟合残差，基于MERIS，ERA-Ⅰ和WRF改正结果的残差标准差分别为0.7135 mm, 0.8083 mm和3.8159 mm.这样看来，基于MERIS和ERA-Ⅰ改正结果的残差标准差都较小，说明这两种改正方法适用于研究区的InSAR大气改正；而基于WRF的大气改正结果的残差标准差远大于其它两种，其残差分布范围更广，这可能是由于青藏高原的地形使得天气预测复杂化所致.青藏高原的海拔超过4000 m，是世界上最高大和最宽广的高原，以复杂的地形著称，分布有大范围的冰雪、高山、冰川、多年冻土和高原湖泊. WRF虽然分辨率较高，但是并不能把复杂的地形和土地覆盖引起的局部动态变化完全考虑在内(Su et al，2013).海原断裂带位于青藏高原东北缘，使用WRF模拟大气环境并非易事(代如锋，2015).另外，如图 3c所示，基于WRF的大气改正结果残差频率分布存在几个极大值，这与误差分布只有一个极大值的假设不相符.考虑到如果所选取的大气改正方法不合适，则可能会引入额外的误差.因此，本文选择改正效果好的MERIS和ERA-Ⅰ大气改正方法，并将其用于老虎山断裂带震间形变速率场的获取中.

图  3  基于3种大气改正方法的残差相位直方图

黑色叉号表示每幅图残差的频率分布(共7幅)，红色曲线表示基于7幅图残差分布拟合的正态曲线

Figure  3.  Histograms of residual phase from the results corrected by three different atmospheric correction methods

Black crosses are the frequency for each interferogram (seven in total), and red curves stand for the fitted normal curve of the residuals of the seven interferogams

下载: 全尺寸图片 幻灯片

2.   老虎山断裂形变速率场获取

2.1   老虎山断裂构造背景

海原断裂带位于青藏高原东北缘，西起祁连山中部，东至六盘山西麓.该断裂带长约100 km，走向约110°，整体呈左旋走滑，部分吸收了印度洋板块向欧亚板块汇聚的位移.老虎山断裂位于海原断裂带的中段、天祝空区的东段，其东侧发生了1920年M_W8.0海原地震，西侧发生了1927年M_W7.8古浪地震(图 4)，近几十年来在该断裂上发生了多次微震(Liu-Zeng et al，2007)，而两边的断裂呈闭锁状态(Jolivet et al，2012)，潜在地震危险性较大(Gaudemer et al，1995)，近年来一直被关注.

图  4  海原断裂带的地质构造图和试验数据分布

背景为来自SRTM (Shuttle Radar Topography Mission)的高程，黑色虚线表示Envisat ASAR数据和MERIS数据的覆盖范围，黑色实线矩形框表示WRF第三层嵌套的范围，红色方块表示ERA-Ⅰ的数据分布

Figure  4.  Tectonic map of the Haiyuan fault zone and distribution of test data

Fault traces are superimposed on SRTM DEM. The black dashed rectangle shows the coverage of Envisat ASAR data and MERIS data. The black solid rectangle stands for the coverage of WRF simulation. Red squares represent data from ERA-Ⅰ. And the epicenter of two large historical earthquakes are marked as magenta points

下载: 全尺寸图片 幻灯片

2.2   老虎山断裂的数据条件及处理流程

本文搜集了覆盖老虎山断裂的雷达数据和大气数据.雷达数据的时间跨度为2003—2010年，共7年的数据资源，为选择最佳数据资源获取可靠的形变数据提供了丰富的数据基础. 图 5给出了干涉数据的基线分布情况，可以看到较多的长时间基线干涉对，这些干涉图是后续进行形变场获取的主要数据源.由于这些干涉对具有较短的空间基线，能够最大程度地减小几何去相关和地形的影响，从而使干涉图保持好的相干性；较长的时间基线，保证了干涉图内积累的形变信号相对较强，这样在后续利用层叠法(stacking)进行累积平均时，更容易捕捉到震间形变信号.

图  5  降轨数据061条带覆盖老虎山断裂的ASAR数据时空分布

红色圆圈表示条带内的雷达数据，其中红色实心圆表示无云情况.红线表示可由MERIS数据改正的长时间基线干涉图，绿线表示可由ERA-Ⅰ数据改正的长时间基线干涉图

Figure  5.  Spatio-temporal baseline plot of Envisat ASAR acquisitions for descending track 061 covering the Laohushan fault

Acquisitions are numbered time sequentially, and 11 of them marked with red solid circle are largely cloud-free. The two long-temporal-baseline inter-ferograms denoted by red line are corrected for atmospheric effects using MERIS data, and the four long-temporal-baseline interferograms marked by green lines are corrected by using ERA-Ⅰ data

下载: 全尺寸图片 幻灯片

基于研究区雷达数据的基线分布，本文选取长时间基线、短空间基线的干涉图进行差分干涉处理，从处理结果中选取相干性较好的作为备选之用. SAR原始数据的聚焦、成像以及影像配准、基线估计、去除平地相位、滤波、解缠、地理编码等处理是利用ROI_PAC软件(Rosen et al，2004)处理完成.地形数据来自SRTM数字高程模型(Farr et al，2007)，精密轨道数据来自于DORIS轨道数据.为提高相干性，在相位解缠前进行了多视处理.

从MERIS和ERA-Ⅰ产品里提取了与每幅SAR影像对应的大气延迟数据，但是由于受MERIS数据可用性和干涉图相干性质量的限制，只有两幅长时间基线干涉图可采用MERIS数据进行大气改正；ERA-Ⅰ不受云的影响，有4幅长时间基线干涉图可参与计算.

上述大气改正方法的比较表明，MERIS和ERA-Ⅰ的改正效果较好，因此本文将采用MERIS和ERA-Ⅰ数据对老虎山断裂的数据进行大气改正.首先对两对无云条件下的干涉对数据(图 5中红色基线所示)利用MERIS数据进行大气改正，然后进行线性轨道误差拟合和去除，最后得到两幅残余相位图，对其进行累积平均得到MERIS数据改正下的形变速率场；对于其余4对长时间基线干涉图(图 5中绿色基线所示)也进行相似的大气和轨道误差改正.不同的是，我们利用ERA-Ⅰ数据进行大气改正，并对得到的4幅残差相位图进行累积平均处理，得到了ERA-Ⅰ改正后的形变速率场.整个数据处理流程如图 6所示. 图 7针对MERIS和ERA-Ⅰ改正，分别以一对干涉图作为实例对改正过程加以说明.

图  6  利用长时间基线干涉图获取震间形变场的流程

Figure  6.  Flowchart of deriving interseismic deformation field based on long-temporal-baseline interferograms

下载: 全尺寸图片 幻灯片

图  7  以一个基于MERIS改正的干涉对(040522-100116，上)和基于ERA-Ⅰ改正的干涉对(041218-100116，下)为例展示图 6所示的形变速率场获取过程

从左至右分别为原始干涉图、湿延迟、流体静力学延迟、经过大气改正的干涉图、拟合得到的轨道误差和轨道改正后的残余误差

Figure  7.  Deriving process of the deformation field from one interferogram pair (040522-100116, upper panels) corrected by MERIS data and the other one (041218-100116, lower panels) corrected by ERA-Ⅰ data as shown in Fig.6

From left to right, shown are original interferograms, the wet delay from MERIS and ERA-Ⅰ data, the hydrostatic delay from MERIS and ERA-Ⅰ data, the atmosphere-corrected interferograms, linear orbital ramp, and residual interferograms after orbital correction

下载: 全尺寸图片 幻灯片

2.3   老虎山断裂形变速率场结果及分析

基于图 6所示的处理流程，本文获取了分别经过MERIS和ERA-Ⅰ数据大气改正的InSAR震间形变速度场，如图 8a和图 8b所示.可以看出，两个结果给出了相似的形变速率场.值得注意的是，对两个相互独立的数据集所采用的大气改正方法并不相同，但得到了相似的形变速率场，验证了本文结果的可靠性.从整体上来看，断层南盘向着卫星视线方向运动，而北盘背离卫星方向运动，显示出了左旋走滑的运动趋势，这与地质研究结果(国家地震局地质研究所，宁夏回族自治区地震局，1990)和GPS观测结果(Gan et al，2007)一致.跨断层剖面显示，上下盘的相对运动速率为视线向2.5 mm/a, 假设其完全由平行于断层的走滑运动引起，则可以转换为平行于断层的走滑速率约6.5 mm/a.相对于南边的老虎山断裂，北边的香山—天景山断裂运动缓慢.另外，一个显著的特征是在老虎山断裂附近存在一个明显的高速形变梯度区(如图 8c中的AA′所示).为定量地分析断层的物理性质，对震间形变场提取两条剖面，一条横跨老虎山断裂(AA′剖面)，另一条(BB′剖面)位于老虎山断裂东边即景泰断裂西段.从AA′剖面(图 8c)可以看出，在老虎山断裂附近2—3 km范围内，形变速率出现了突变，这说明有浅层的断层蠕滑存在，与Jolivet等(2012)的观测结果一致.为了研究其闭锁深度和深部滑动速率，基于arctan模型(Savage，Burford，1973)，反演了4个剖面的数据.对于老虎山断裂：基于MERIS改正和轨道改正的剖面结果反演得到的滑动速率为6.9 mm/a，闭锁深度为4 km；基于ERA-Ⅰ改正和轨道改正的剖面结果反演得到的滑动速率为7.0 mm/a，闭锁深度为3 km.对于景泰断裂：基于MERIS改正和轨道改正的剖面结果反演得到的滑动速率为7.6 mm/a，闭锁深度为5 km；基于ERA-Ⅰ改正和轨道改正的剖面结果反演得到的滑动速率为8.4 mm/a，闭锁深度为7 km.这些结果与已有的InSAR (Cavalié et al，2008)和GPS (Gan et al，2007)研究成果相符合.

图  8  InSAR震间形变场和形变场剖面线

(a)和(b)分别表示基于MERIS改正和ERA-Ⅰ改正的震间形变场，其中红色实线矩形框内的断层为老虎山断裂，红色虚线矩形框内的断层为景泰断裂；(c)和(d)分别表示速度场剖面线AA′和BB′平行于断层的水平运动速率，其中黄点和绿点分别表示基于MERIS和ERA-Ⅰ大气改正后的形变场在剖面线上的投影，而红线和黑线分别表示对应其平均值

Figure  8.  Interseismic velocity field of study area and its two profiles across the fault

(a) The MERIS-corrected rate map; (b) The ERA-Ⅰ-corrected rate map. The red solid rectangle covers the Laohushan fault, and the red dashed rectangle covers the Jingtai fault; (c) The fault-parallel velocity of the AA′ profile from two rate maps; (d) The fault-parallel velocity of the BB′ profile from two rate maps. Black and green points are the projection of MERIS-corrected and ERA-Ⅰ-corrected rate map on the profile, and red line and black line are their mean values

下载: 全尺寸图片 幻灯片

3.   讨论与结论

本文从实际案例出发，通过试验对3种外部大气数据进行了评价，选择最适合的大气数据进行了试验区干涉图的大气改正，并成功地将其应用于老虎山断裂InSAR震间形变速率场的获取.通过试验得出以下结论:

1) 老虎山地区干涉图的大气改正效果显示，通过引入独立的外部大气数据来估计干涉图中的大气延迟，能够很好地解决长波震间形变信号与长波大气信号之间的耦合问题；

2) 试验区最优大气改正方法的选取及其成功应用，说明基于研究区本身短时间基线干涉图可以很好地评价外部大气数据对该区域大气延迟估计的准确性，为准确判断何种数据适用于研究区的大气改正提供了可靠的参考依据；

3) 从老虎山区域的试验结果可以看出，在进行合理的大气改正和轨道改正的基础上，即使利用有限的长时间基线干涉图(几幅)，InSAR层叠法也能给出可靠的震间形变速率场，这给SAR数据资源不是很丰富的区域的InSAR形变场获取带来了希望.

另外，无论是对于EnviSAT和PALSAR等老卫星，还是ALOS-2和Sentinel-1A等新卫星，利用InSAR处理卫星影像获取更精确的震间微小形变，大气改正均必不可少.在充分利用各种类型的水汽数据、地形数据的基础上，探索具有更强适应性的大气改正方法，是InSAR大气改正发展的方向.随着新卫星的发射、影像分辨率的提高和数值大气模型的发展，InSAR大气误差将会有更合理的解决办法(崔喜爱等，2014).