多轨道InSAR震间形变速率场拼接方法

华 俊 龚文瑜* 单新建 王振杰 季灵运 刘传金 李永生

1)中国地震局地质研究所，地震动力学国家重点实验室，北京 100029 2)中国石油大学(华东)，青岛 266580 3)中国地震局第二监测中心，西安 710054 4)国家自然灾害防治研究院，北京 100085

0 引言

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)是一种高精度、 高效率的地表形变测量方法，常用于地表形变的长期观测，如地震构造观测(徐小波，2016； Li Betal.，2020； Li Yetal.，2020)、 火山活动观测、 地表沉降监测(张永红等，2016)、 冰川漂移和山体滑坡(王桂杰等，2010)等。但单一观测轨道SAR影像的标准幅宽范围有限，常常无法覆盖整个形变区域。以青藏高原地区为例，其内部分布很多EW走向的大型走滑断裂(如甘孜-玉树断裂带、 嘉黎断裂带等)，走向长度往往达上百千米(Armijoetal.，1986； 闻学泽等，2003； 闫欢欢，2017)，需要多个同一轨道(升轨或降轨)的SAR数据条带(以下均简称为多轨道InSAR)才能完全覆盖。近年来，关于多轨道InSAR形变拼接的研究主要集中于以垂直运动为主要观测目标的地表沉降监测应用领域(Ketelaaretal.，2007； Geetal.，2010； 葛大庆，2013； 罗三明等，2014； 熊思婷等，2014； 孙赫，2015)，而运动较为复杂的构造形变观测领域的研究相对较少。在构造形变观测方面，以Chaussard等(2015)的研究成果为例，一般采用时间序列方法分别计算2个相邻降轨数据的视线向(Line of Sight，LOS)速率后，直接对2个InSAR数据进行拼接。Wang等(2007)首先将多个条带模式的InSAR LOS向形变场转换到水平方向，然后再恢复为统一入射角的LOS干涉图，重建了玛尼地震的同震形变场。Wang等(2019)在藏南地区的研究中未直接对InSAR数据进行拼接，而是联合GPS速度场并结合物理模型重建了大尺度应变率分布和速度场。Weiss等(2020)通过SAR卫星观测几何参数，将全球导航卫星系统(Global Navigation Satellite System，GNSS)的三维形变转换到卫星视线向，并计算其与对应区域的InSAR观测值的差值，通过二次多项式拟合统一InSAR和GNSS资料的参考基准，然后对多个轨道进行直接拼接。这些研究对获取大范围InSAR形变场的应用研究提供了很好的思路，但并没有针对现代宽幅卫星数据定量考虑入射角等其他导致差异的信号源对多轨道 InSAR 数据拼接的影响。

聚焦于震间形变研究中广域形变场重建的应用需求，本文以哨兵1号卫星数据为例，首先分析了多轨道InSAR数据重叠区域产生形变差异的原因，包括入射角、 参考点、 非形变信号等； 然后基于模拟数据集，定量分析入射角对多轨道InSAR震间形变速率场拼接的影响，提出基于多项式估计的比值法修正入射角； 讨论了针对青藏高原走滑断裂的震间形变特点，将数据转换到地距向(Azimuth look direction，ALD)后进行拼接的可能性。最后，本文以青藏高原东南部为实验区域，对多轨道单一方向InSAR震间形变速率场进行了拼接并开展定量分析。

1 多轨拼接的误差分析

以相邻轨道重叠区域形变速率差为例，同一位置形变信号差φdiff的相位贡献可以用数学模型表示为

φdiff=φl-φr=φunw-ref+φinc-diff+φatm+φt-diff+φres

(1)

式中，φl和φr分别表示为图1 中对应的左侧轨道(以下称为主轨道)和右侧轨道(以下称为辅轨道)的形变速率相位。

图1 哨兵1号卫星数据相邻轨道重叠区域入射角变化示意图

φunw-ref表示不同轨道空间参考带来的差异，这是由于各独立轨道内的干涉解缠相位是相对于影像覆盖范围内某一空间的起算点，因此，不同轨道的空间起算点的差异导致各个轨道解算的形变速率存在明显的空间基准偏差。这种参考点不一致带来的偏差可以通过联合GPS水平速度场成果进行统一。目前，已有一些基于GPS数据校正InSAR形变速率场的方法。例如，Meng等(2010)提出了一种去除、 滤波、 恢复方法，结合GPS和InSAR数据，获得了高信噪比的震间地壳形变场。另一种方法是以GPS和InSAR之间的差值作为观测值，使用二阶多项式消除InSAR形变速率中的长波长误差，从而统一参考基准(Hussainetal.，2018； Neelyetal.，2020； Weissetal.，2020)。在下文的研究中，我们应用二阶多项式拟合GPS和InSAR形变速率之间的差值，以统一参考基准。

φinc-diff表示侧视条件下入射角的差异。SAR卫星的侧视成像模式导致不同轨道保持独立的侧视空间几何参考系统。尤其是随着现代宽幅SAR传感器的发展，SAR卫星入射角从近距离到远距离的变化增大，进一步加大了相邻轨道InSAR形变速率场的差异。表1 给出了不同SAR卫星传感器的多种模式数据近距离到远距离的入射角范围(Torresetal.，2012； 朱建军等，2017； Wangetal.，2020)。其中，ERS等历史卫星的入射角仅为20°～26°； 欧洲航空局 Sentinel-1卫星(1)https： ∥earth.esa.int/eogateway/。TOPS模式数据的入射角为29°～46°； 而日本ALOS-2卫星(2)https： ∥www.eorc.jaxa.jp/ALOS/en/。ScanSAR模式数据的入射角可达32°～54°，其入射角有22°的差异。图1 以哨兵1号卫星数据为例，直观地显示了SAR卫星侧视成像导致的卫星视线向的空间参考差异所带来的影响。

表1 SAR卫星传感器的主要参数

φatm表示残余的大气误差引起的信号差异。这是由于不同轨道的SAR影像获取时间存在差异，大气条件发生变化，在提取形变信息时大气误差无法完全消除(刘斌，2013)，导致残余的大气误差与形变相位耦合在一起，大气校正的好坏直接影响多轨道InSAR形变速率场的拼接效果。φt-diff表示SAR数据由于获取的时间跨度不完全一致引起的形变差异，尤其是在时间上非线性的地表形变会加剧这一差异。φres表示失相干噪声等其他残余误差引起的形变差异。

因此，为了获取大范围InSAR形变场，考虑到构造形变可能在时间上存在非线性变化(如同震、 震后位移等)，应该尽量使用具有相似时间跨度和时间采样的多轨道数据集，以减弱如残余大气误差、 时间不一致等引起的形变差异影响。下文将首先针对入射角差异讨论并分析多轨道拼接的策略。

2 入射角效应研究

卫星视线向InSAR形变速率场实际上是三维(three dimensional，3D)形变速率到LOS向的投影，两者之间的关系可由式(2)表示(Hanssen，2001)：

(2)

式中，dLOS表示LOS向形变速率，du、dn和de分别表示垂直地表、 SN和EW向的形变速率，θinc和αh分别表示雷达卫星相对目标位置的入射角和方位角。

显然，如果研究区域存在至少3个或以上不同视角的InSAR形变速率场，可通过式(2)将其分解为三维形变速率，然后分别对三维形变速率进行多轨道拼接。但在实际应用中，往往很难获取研究区域多个不同视角且大面积的星载SAR观测数据。更常见的情况是，仅有单轨道(升轨或降轨)或2个轨道(升、 降轨道组合)观测数据，在这种条件下，完成三维形变的分解往往需要假设垂直或水平位移可以忽略(Mohretal.，1998)，或者通过外部数据加以约束(例如GNSS形变观测产品)(Gudmundssonetal.，2002； Mengetal.，2010； 宋小刚等，2015； 曹海坤，2017)。对于多轨道InSAR速率场拼接，本文首先通过震间形变场模拟数据定量分析多轨道拼接中的误差特点并讨论有关的解决方案。

2.1 震间形变场模拟数据和统计分析

模拟数据集的构建以甘孜-玉树断裂带(左旋走滑断裂)为例。这里使用二维反正切弹性位错模型(Savageetal.，1973)模拟走滑断裂的震间形变速率场(图2a，b)：

(3)

式中，x表示与断层间的垂直距离，s表示断层深部的滑动速度，d表示断层的闭锁深度，V(x)表示距离断层x处的滑动速度，这里假设断层长度为75km。孟国杰等(2016)和吕丽星等(2017)根据构造地貌调查得到的长期滑动速率约为7mm/a，整体闭锁深度约为20km。

我们根据式(3)计算获取了沿着断层走向的滑动速率，并基于断层走向分解得到了SN和EW向的滑动速率(图2a，b)。将二维震间形变速率场分割为带有重叠区域的3部分，基于哨兵1号卫星的空间几何参数(入射角范围为28.9°～43.0°)，使用式(2)分别投影到相邻2个轨道的LOS向(图2c，d)。可见，地表形变速率大的区域，2个轨道间信号的目视差异也较大。

图2 模拟得到的震间形变速率场

基于以上模拟数据构建了走滑断裂的震间形变场，忽略垂直形变速率，则式(2)可改写为

dLOS=-sin(θinc)DN(rangeline)-sin(θinc)DE(rangeline)

(4)

(5)

图3 邻轨重叠区域形变速率差的比值和直方图分布

2.2 入射角影响的修正策略

2.2.1 基于多项式改正的方法

图4 比值与入射角的关系

(6)

式中，(x，y)表示重叠区域像素点EW向和SN向距离，a、b、c为待估参数。

针对图2 中的模拟数据，采用比值法进行入射角差异改正，改正前、 后的结果如图5a、 5b所示。其中，图5a 为2个相邻轨道直接叠加后原形变速率场的模拟结果，图5b 为入射角改正后的结果，图5c 和图5d 分别为加入大气湍流后的原形变速率场和经入射角校正后的形变速率场。对比图5a 和图5b 可知，在对入射角进行修正后，相邻轨道的InSAR形变速率场在重叠区域具有较好的一致性，且均值和标准差都为0，说明比值法可以有效地减弱入射角的影响。

图5 基于模拟数据的邻轨InSAR形变速率场拼接

此外，考虑到实际应用时InSAR形变速率场中的大气误差无法完全消除，残余大气误差仍然存在于InSAR数据中(Mengetal.，2010； 刘斌，2013； 占文俊等，2015； 姜宇，2017； Neelyetal.，2020)。我们基于功率谱法(Hanssen，2001)向InSAR震间形变速率场中增加了模拟的大气湍流信号(图5c)，其均值为0.38mm，标准差为0.30mm，以探讨残余大气误差对邻轨InSAR形变速率拼接的影响。可见，重叠区域InSAR形变速率的不一致性进一步增大。采用同样的方法进行入射角改正后，发现虽然比值法能够对入射角的差异进行一定程度上的修正，但由于大气信号的存在，拼接质量也受到了一定程度的影响。因此，在进行多轨道InSAR形变速率场拼接之前，尽可能减弱大气等误差具有重要意义。此外，比值法是将相邻轨道纳入到统一的侧视观测坐标系内的手段，在该坐标系下远场区域的入射角相对于原侧视坐标系增大，将导致信号发生变化(如图5d 中的红框所示)，对目视效果带来影响。

2.2.2 ALD形变转换与多轨道形变场拼接

针对研究目标——构造走滑断裂震间形变的应用场景，这里讨论在地距向(Azimuth look direction，ALD)进行多轨道InSAR震间形变速率场拼接的策略。

为了减少入射角带来的影响，将LOS向形变转换至ALD方向，即将式(2)改写为式(7)。如果研究区域地表形变以水平形变为主，假设垂直贡献可忽略不计，进一步可得到式(8)。可见，ALD向形变速率仅与卫星的方位角有关，因此可以在ALD向对InSAR形变产品进行拼接。

(7)

(8)

该策略与Wang等(2007)介绍的方式类似，其在对玛尼地震的同震形变场进行拼接时，首先将LOS形变场投影到ALD方向，然后再将其全部投影到23°入射角方向，恢复为统一入射角的干涉形变场。值得注意的是，他们在转换过程中并未忽略垂直位移，而是通过统一入射角重投影削减了垂直贡献导致的差异。由于其处理历史卫星条带模式数据的入射角近-远距离变化较小，因此相位差异本身影响也较小。考虑到哨兵1号卫星相邻轨道间入射角差异较大，采用类似的方法进行处理时，垂直位移导致的差异系数最大可约达0.2(假设统一入射角为30°)； 而当地表垂直形变不大时，2种方法均能产生较好的拼接结果。考虑本研究的目标场景(以走滑运动为主的断裂)，在下文研究中除将形变场转换至ALD方向外，同时假设垂直位移可以忽略。

该处理策略方式简单，实现过程中除需获得各个轨道的形变场外，还需逐点获取对应的入射角。一般在计算InSAR震间形变速率场的同时，可基于卫星的几何参数获取对应点的入射角以及方位角。本研究中，由于是直接针对形变速率场结果进行分析和处理，且考虑到较小的入射角变化并不会引起较大的LOS向形变误差(Shirzaeietal.，2015)，因此我们计算哨兵卫星数据的最小入射角和最大入射角，再利用线性多项式拟合的方式获取对应各像素点的入射角。

3 基于真实数据的多轨道InSAR形变速率场的拼接

为了检验上述策略的可行性，本文以青藏高原东南部地区为实验区，进行大范围InSAR形变速率场拼接实验和分析。实验区域范围为(27°～35°N，91°～99°E)(图6)，总体以水平运动为主(Liangetal.，2013； Zhengetal.，2017)。该地区是目前中国地震灾害较为频繁的地区，且于2010年4月14日发生了玉树MW7.1 地震(图6)。

图6 研究区域示意图

实验区域包括 3 个二级块体，从北向南分别为巴颜喀拉、 羌塘和拉萨块体(Lovelessetal.，2011)。其中，巴颜喀拉块体是青藏高原隆起发育的前缘带，是中国南北地震带的中段，构造变形极其复杂，断层活动强烈，控制了一系列历史强震的发生(Armijoetal.，1986； Liangetal.，2013)。嘉黎断裂带是一条近EW走向的右旋走滑断裂，全长约500km，呈向N突出的弧形，为青藏高原主体的南部构造边界，以1.5～6mm/a的速度向E挤压(Armijoetal.，1986； Liangetal.，2013； 陈长云，2014)。甘孜-玉树断裂带自晚第四纪以来以左旋走滑运动为主，滑动速率为(6.6±0.1)～(7.4±1.2)mm/a(王阎昭等，2011)。澜沧江断裂带以右旋走滑运动为主，滑动速率约为2mm/a(闫欢欢，2017)。此外，东喜马拉雅造山带是青藏高原构造应力较大、 冰川活动强烈和形变最强的区域，也是地球上变化剧烈、 构造类型保存较为完整的地区，运动呈现强烈隆升和向S挤压的特征(唐方头等，2010； Guptaetal.，2015)。

3.1 数据集

本文采用Liang等(2013)处理得到的GPS形变速率场，并基于此统一基准。他们采用了750个GPS观测数据得到了相对于欧亚参考框架下的青藏高原三维形变速度场。如图6 所示，黑色三角形和红色三角形是InSAR覆盖范围内GPS点的空间分布。本文使用的3个降轨 InSAR 形变成果基于哨兵1号卫星数据处理得到，其时间跨度大致为2015年6月—2019年1月。图9a 显示了原始InSAR形变速率场的空间分布，其中，每个轨道由5、 6个标准景组成。干涉图的处理采用哨兵1号卫星精密轨道数据进行了基线精化，并利用30m分辨率的SRTM DEM数据进行模拟并去除平地和地形相位； 最后利用相位叠加平均法(stacking)(Sandwelletal.，1998)求解该区域高精度的地表形变速率场(图9a)。由于大气信号的影响，各InSAR形变速率场的精度存在差异。以T77和T4为例，两者都捕捉到了2017年11月17日西藏米林6.9级地震的同震形变信号(图9a 中的红圈)，且邻轨重叠区域形变速率存在不一致的现象，其主要原因是各轨道的参考点不一致，且入射角存在差异。

如图6 中的黑色方框所示，其内包括T77、 T4和T1063个轨道。其中自西向东，3个轨道内分别有11、 10和20个GPS点，可用作地面控制点(Ground Control Point，GCP)以统一 InSAR 数据的参考基准。此外，我们在1984年WGS84坐标系下对多轨道InSAR形变速率场进行拼接。

为了获取GPS观测站所在位置的InSAR形变速率，本文以GPS点为中心，半径为5km，采用反距离加权方法获取各GPS点对应位置的InSAR形变速率值。

3.2 大范围形变场的拼接成果与分析

我们首先基于研究区域的GPS水平速度统一了InSAR形变速率场的参考基准(图9b)。在处理过程中，基于T77、 T4和T106的GPS和InSAR形变速率的差值拟合的二阶多项式的R2分别为0.7、 0.7和0.9。为了评估使用GPS速度校正InSAR数据参考基准的精度，抽取了5个未参与校正的GPS点(图6 中的红色三角形)对改正前、 后的结果进行比较(图7)，改正前、 后GPS和InSAR之间形变差的平均值约为14.22mm/a和0.01mm/a，改正效果明显。

图7 GPS与InSAR形变速率差值基准统一改正前、 后的比较

总体而言，统一基准后，相邻轨道重叠区域InSAR形变速率差的差异明显减小，但仍然存在InSAR形变速率不一致的现象，原因可能是相邻轨道重叠区域存在入射角差异或其他数据处理误差(如相位解缠误差等)。

如图8 所示，我们对所有GPS点与InSAR统一基准前、 后形变速率之间的差异进行了比较。在改正前，GPS与InSAR形变速率存在明显的信号差异，且不同轨道的差异量不同(图8a)，这也说明各个轨道之间的参考点存在差异。统一基准后GPS和InSAR形变速率之间的差异显著减小，校正后GPS和InSAR差值的数学平均值为0.003mm/a，标准差为1.5mm/a(图8b)，表明基于GPS形变速率统一InSAR参考基准具有较好的效果。

图8 基准统一前、 后GPS点形变速率与InSAR点形变速率的比较

在基准统一的基础上，分别采用比值多项式改正法(图9c)和ALD转换的方式(图9d)进行改正，各轨道重叠区域的形变速率总体具有较好的目视效果，比值法和将InSAR数据转换到ALD向拼接对于该区域都是可行的。其中使用入射角改正法时，将T106条带作为主轨道，并将其他轨道统一纳入该轨道的侧视参考坐标系。因此，随着加入条带数量的增加，远场入射角不断增大，将导致远场形变垂直向的贡献被压缩、 水平向贡献被拉伸，在解译的过程中需要加以注意。此外，T4和T106重叠区域南部(图9a 中的红框)的形变速率场在拼接后呈现明显不一致的现象，一方面可能是由于此处相干性差导致形变场解算精度低，另一方面可能是附近区域垂直形变较大(如冰川活动等)带来的影响所致。

图9 多轨道InSAR形变速率场拼接

为了更好地进行定量讨论，在2个重叠区域作了2条剖面(图9)，并计算相邻轨道间沿剖面的InSAR形变速率差(图10a，b)及对应的直方图(图10c—j)。结果显示，经过基于多项式改正的比值法和ALD向形变转换策略处理后，重叠区域形变差的平均值和标准差均较小，2种策略都能够比较有效地减少不同轨道间的偏差。其中ALD向形变转换策略略好于多项式改正法，可见ALD转换的策略更适用于该地区。

图10 剖面AA′和BB′覆盖下的InSAR形变速率差及统计分析

3.3 讨论

前文以覆盖青藏高原南缘的真实InSAR形变速率条带产品为例，利用本文所提出的比值法和ALD方向转换策略进行了拼接测试。总体而言，两者均能对多条带数据进行拼接。根据统计精度可知，ALD转换策略的效果略好于比值法，说明目标区域的实际运动情况以水平运动为主，ALD转换策略的前提假设基本成立。此外，ALD转换策略不存在未知参数估计，其转换方式更为简单直接，而比值法则需要对重叠区域进行多项式拟合，该过程将受到形变速率场中残余误差的影响。

总体而言，ALD方向转换策略对于青藏高原大部分走滑断裂是适用的，但对于具有不可忽略的垂直位移的区域则并不适用，在该情景下可采用比值法进行拼接，或可参考Wang等(2007)提出的将其恢复为统一入射角的LOS干涉图的方式，以减少垂直位移的影响。比值法在处理中将多个轨道统一至单一侧视参考坐标系，随着纳入轨道的增多，重建得到的远场入射角将增大，和近场位置形变相比，远场位置的形变存在较大拉伸，在结果分析中需要加以注意。此外，InSAR形变场靠近山区的位置可能本身就存在解算结果不可靠的问题，需要对低质量像素进行掩膜。例如，本研究中采用基于相干性的掩膜处理，未来还可加入阴影-叠掩图进行掩膜，这对于使用比值法减弱入射角的影响具有重要意义。

4 结论

本研究针对震间形变研究中广域形变场重建的应用场景，首先对邻轨重叠区域InSAR形变速率的差异来源进行了分析，包括空间参考点、 侧视入射角、 其他非形变信号(如残余大气误差、 卫星数据获取时间差异等)。总体而言，在保证各个轨道SAR数据的时间覆盖基本一致的条件下，应首先解决参考基准和入射角的差异问题，从而实现多轨道InSAR形变速率场的融合拼接。本文基于模拟数据的研究，讨论了削弱入射角的影响，提出了更好的实现多轨道数据拼接的办法，包括基于多项式估计的比值法和转换至ALD方向后拼接的策略。通过真实数据的研究，发现基于多项式估计的比值法能够较好地修正入射角的影响，并获取高精度、 大范围的InSAR震间形变速率场，但最终结果均位于单侧的侧视参考坐标系下，容易导致远场形变产生目视效果失真的问题。将数据转换至ALD方向后，也可进行大范围的InSAR形变速率场拼接，该方法适用于以水平位移为主的大型走滑断裂区域。在实际应用中，还应该根据区域目标的运动特点和轨道数量等因素选择以上介绍的2种策略，以实现多轨道InSAR数据的拼接。本文以青藏高原东南部大范围震间形变场重建为主要目标，基于GPS形变速率获取了高精度、 统一参考基准的InSAR形变速率场，并验证了2种多轨道拼接解决策略的可用性，最终采用青藏高原东南部地区的真实数据实现了高精度、 大范围InSAR形变速率场重建，这对研究大范围、 高精度构造形变的空间变化具有重要意义。

致谢欧洲航空局为本研究提供了哨兵1号卫星数据； 国家自然灾害防治研究院李永生副研究员提供了T77轨道数据； 中国地震局第二监测中心季灵运研究员和刘传金助理研究员提供了T4和T106 InSAR数据； GMT中文社区为本研究提供了断层数据支持； 文中部分图件采用GMT绘制。在此一并表示感谢！