基于D-InSAR技术的同震位移场反演及验证

郭世鹏，康伟，黄继茂，张庭苇，李云，张王菲

(西南林业大学 地理与生态旅游学院，云南 昆明 650224)

0 引 言

SAR干涉测量技术(interferometric synthetic aperture radar，InSAR)可以监测地球陆地表面的微小变化，精度可达毫米级，同时SAR技术具有全天时、全天候和大面积监测等优点。由于使用该技术不需要人进入监测地区，因此逐渐成为获取地震发生区震后位移测量的首选方法[1-2]。D.Massonnet等[3]采用时隔几个月的ERS-1卫星干涉像对揭示了Landers地震的同震位移场，基于此获得的形变图与GPS监测结果非常一致；张红等[1]采用ERS-1/2干涉像对提取了张北地区的地震同震位移场，进一步证实了差分干涉SAR(D-InSAR)技术在地震监测和研究中的潜力；臧妻斌等[4]利用两轨法对意大利拉奎拉地震形变展开研究，成功获取了同震位移场在雷达视线方向(light of sight,LOS)的位移；张瑞等[5]应用D-InSAR技术得到了玉树地震区域的地表形变,并对地震的震源机制和发震机理展开分析。D-InSAR技术由于干涉像对获取的容易性和获取位移场算法的简单性，在未来的地震监测研究和监测中具有极大潜力。这些研究证实了D-InSAR技术在地震形变场监测中的有效性，同时也指出由于SAR成像方式的影响，即部分区域影像噪声较高，相干性低，使得监测结果精度降低，因此需要地表测量数据进行相应的校正或验证[1-2,6-8]。由于地震产生的同震位移场的区域通常较大，地表合适的测量数据(如GPS、水准测量等测量数据)获取困难，难以检验采用D-InSAR获取的地震形变信息的不确定性和精度[2]。近年来针对升降轨数据获取的地表形变特征及采用弹性半空间模型(Okada模型)模拟三维位移场相融合的方法为以上问题提供了解决思路[9-13]。本文基于D-InSAR获取地震位移场的成像几何关系和测量原理，建立Okada模型模拟地震三维位移场与D-InSAR在LOS方向的数学关系，通过升降轨相关性分析、Okada模型模拟地震位移场来间接验证D-InSAR在地震引发的同震位移反演中的有效性。

1 研究区概况及数据获取

1.1 研究区概况

研究区位于中国台湾省高雄市(纬度22°09′43″N～23°21′20″N，经度120°06′25″E～120°53′38″E)，震中位于美浓区(22°55′48″N，120°32′24″E)，如图1所示。地震区地处左镇断层、潮州断层及旗山断层之间的地震密集带，距高雄市区约30 km。此次地震达到里氏6.7级，震源深度为15 km。

图1 研究区地理位置

1.2 研究区D-InSAR数据获取

1.2.1 D-InSAR数据获取

本研究获取了升轨和降轨的Sentinel-1A卫星InSAR数据各一景。Sentinel-1是由位于同一轨道平面内，轨道相位差为180°的1A和1B组成的双星星座，双星星座重访周期为6 d，A，B星各自重访周期为12 d，高纬度地区最快重访时间可达1～3 d。Sentinel-1A卫星搭载一台C波段(频率为5.404 GHz)合成孔径雷达，包含条带(SM)、干涉宽幅(IW)、超宽幅(EW)、波浪(WV)4种成像模式，最大覆盖幅宽可达400 km。其中IW模式为Sentinel-1A的主要模式，该模式成像由9个burst拼接而成，幅宽为375 km。[14]本研究获取了Sentinel-1A卫星IW模式的2对较短时间间隔的干涉SAR影像对，详细参数信息见表1。升降轨覆盖研究区的位置见图1(b)。

表1 SAR影像信息

1.2.2 其他辅助数据获取

本次研究同时获取了辅助的DEM数据(图1(a))，用于模拟地形相位和地理编码，该数据为美国航空航天局(NASA)提供的SRTM DEM(分辨率为30 m)数据。为了提高轨道误差对基线估算精度的影响，本文同时获取了两景影像的精密轨道数据(precise orbit data，POD)，用于轨道误差改正，进而提高基线估算的精度。

2 InSAR数据处理

2.1 InSAR数据处理流程

InSAR数据处理流程如图2所示，包括两个步骤：第一步，采用二轨差分法分别对升降轨数据进行LOS方向地震前后二维位移量的提取；第二步，采用升降轨数据和Okada模型反演形变信息并进行二维位移场提取结果验证。具体处理环节如下。

图2 InSAR数据处理流程图

第一，首先利用ENVI软件处理InSAR单视复数(single look complex，SLC)数据；第二，根据获取的精密轨道数据对用于基线估算的轨道信息进行修正，得到更精确的轨道参数；第三，通过距离向为1和方位向为5的多视处理提高相干性的统计精度；第四，采用Goldstein滤波对多视后的影像进行去噪[15]，最小费用流(minimum cost network flow phase unwrapping algorithm，MCF)进行相位解缠获取形变相位；第五，经过地理编码获得在LOS向的位移信息；第六，分别使用升降轨数据获得的位移信息进行相关性分析；第七，将降轨数据获得的位移信息带入Okada模型模拟地震同震场的信息,以此验证采用D-InSAR反演的位移场位移结果是否具有现实指导意义[8]。

2.2 Okada模型算法

Okada通过分析各种类型的断层错动引起位移的解析表达式，剔除了通用点源及有限面源弹性模型，以走滑方向、倾滑方向和张性方向断层单位位错描述的弹性半无限空间表征地表位移公式，进而模拟地震场的时空演化过程[9-10,12]。该算法描述的位移场ui可表示为

(1)

图3 Okada模型坐标系

(2)

(3)

(4)

(5)

(6)

(7)

式(5)中,

Sloc=aU1+bU2+cU3，

(8)

其中，a,b,c可以通过雷达成像几何关系和U1,U2,U3几何关系计算获得。

3 结果与讨论

3.1 InSAR技术提取结果

图4为覆盖整个地震形变区域的干涉纹图，上盘干涉条纹密集且闭合，呈现4个同心圆环状，一个同心圆代表半个雷达波长的形变量。下盘干涉条纹稀疏且条纹不太明显，存在非闭合条纹，右下部分地区呈现较大失相干现象。从整体看，升降轨获得的干涉图大致相同，反演的同心圆数目也相同。但由于卫星拍摄角度等多种因素和植被山体叠掩现象，升轨干涉纹图和降轨干涉纹图略有区别。虽然上盘和下盘均形成干涉条纹，但两者差异明显。上盘干涉条纹和下盘条纹的差异表明此次地震为逆冲型地震，该结果与全球多家监测机构给出的相关震源机制解释一致。

图4 地震形变区域干涉纹图

笔者将升降轨两组干涉纹图经过相位解缠、地理编码到相同地理坐标系下，并将相位值转为位移量(图5)。图5(a)为升轨影像提取的LOS方向的位移信息，图5(b)为降轨影像提取的位移信息。从图5 D-InSAR反演的同震位移场可以看出，两组D-InSAR结果较为一致，上盘主要表现为抬升，下盘则为下沉。升降轨位移场不仅具有相同的位移趋势，并且上盘和下盘的位移值也基本吻合，其中升降轨监测的位移图中在断层的上盘区域分别抬升14.3 cm和17.2 cm，下盘分别下降10.1 cm和15.7 cm，两者在上盘和下盘的抬升和下降主要是由于卫星获取影像时间差异和升降轨入射角度差异造成的。由于SAR斜距成像，对距离的变化较敏感，因此断层西部上盘在LOS方向远离升降轨，表现为较大抬升；相反，东部下盘LOS方向远离降轨，靠近升轨，使得降轨下沉值大于升轨下沉值。

图5 InSAR同震位移场提取结果

从地质构造看，台湾西南部地处复杂的破碎活动断层，并且西南地区发生的地震基本都与这些断层走滑相关。本次地震震中位于旗山断层和潮州断层的交界附近，结合升降轨LOS位移场的位置关系可以看出，同震位移场的沉降位移区位于地震带的东南盘，而隆升位移区位于地震带的西北盘，这一特征与逆冲断层的性质相符[14-16]。此外，沉降区与隆升区分界明显，并且分界位置处于旗山断层带附近。因此，推测出此次地震很有可能是由多个断层带位移综合引起的。由于研究区所在的台湾岛地处欧亚板块与菲律宾板块的相互挤压带,构造复杂，在北部和东部，菲律宾海板块沿琉球海沟向西北偏北的欧亚大陆俯冲；岛屿以南，南海(在欧亚板块上)在马尼拉海沟向东延伸至菲律宾海板块之下，板块的挤压作用导致台湾地形南北狭长。在地震发生的地方，菲律宾板块以每年8 cm左右的速率向西北-东南方向推挤，因此推测此次地震是菲律宾海板块与东亚大陆在台东纵谷碰撞仰冲而强烈向西推压影响的结果[17-18]。

3.2 升降轨互补验证位移精度

由于缺乏当地实际的水准测量数据，本次研究基于统计方法从沿最大抬升区和最大沉降区验证D-InSAR提取的震后位移信息的精度。消除地形误差后，在升降轨数据的反演结果中共选出54对典型位移点进行相关性分析，分析结果如图6所示。两组数据之间相关性较高，R2=0.978 4。升降轨数据的高相关性表明采用Sentinel-1A数据的D-InSAR技术在反演地震同震位移场中LOS方向位移的有效性。

图6 升降轨位移结果相关性分析

3.3 Okada模型反演位移场验证精度分析

为进一步分析地震机制与D-InSAR反演地震位移场的精确性，研究选用降轨数据反演的位移场提取用于Okada模型的反演震源参数。将D-InSAR同震位移场进行采样后，参考CMT(centroid moment tensor)提供的经度、纬度、走向、倾向、滑动角、深度、断层长度、宽度以及滑动量等求解非线性反演断层。

依据CMT官网推荐的地震参数推荐区间反复迭代，不断缩小区间范围，直至得到合适搜索区间，并根据搜索区间范围确定参数值进行反演，将反演结果的参数设为最优值。经过Okada模型的非线性反演，得到如表2所示的参数值。表2同时比较了反演结果与CMT及USGS(united states geological survey)公布的震源参数，从表2可知，反演结果与CMT接近，与USGS公布的震源深度和倾角有一定差异，但从地震海滩球、震级、地震矩综合分析可知，三者一致性较好。

表2 CMT、USGS和Okada形变场参数

采用Okada模型反演出的地震产生的破裂长度为9 890.4 m，宽度为3 976.2 m，震源深度为14.59 km，震级为MW6.4，累计释放地震矩为5.25×1018N·m，非线性正演断层模型反演结果如图7所示。

由图7断层滑动模型结果可知，地震形成的断层面主要位于震中西北方向，并处在潮州断层和旗山断层附近，断裂发生在向西北-东南方向。向东北方向浅倾，或在向西陡倾的南北走向构造上。此外，地震还伴有走滑现象，断层的滑动角为1.26°，倾角为22°，滑动量为4.47 m。由此可以断定，本次地震的发震机制为逆冲断层，这与D-InSAR获取的逆冲走滑断层发震机制相符，进一步说明了D-InSAR反演同震位移场的可行性。

图7 Okada模型地震参数反演结果

由以上分析可以确定，高雄地震由该区域断层带引起，发震机制为逆冲走滑断层，滑动量为4.47 m，InSAR得到上盘最大抬升量约为17 cm，下盘最大沉降量约为 12 cm。已有研究采用ALOS-2 PALSAR-2的D-InSAR的结果表明，该区域的最大抬升高度约为13 cm，与本研究的研究结果相差在5 cm内[8]。而同样采用2016年2月14日和2月02日的Sentinel-1A的反演结果表明研究区总体位移量为-4～12 cm，与本研究的结果基本一致[7]。

4 结 语

基于Sentinel-A升降轨数据，利用D-InSAR技术得到LOS向同震位移场，基于同震位移场及干涉相位图，确定了研究区地震场的范围，地震引起的地面抬升和沉降程度，结合研究区的地质构造特征，揭示了地震的形成机制。为了验证D-InSAR技术得到同震位移场的有效性，采用54个典型点的升降轨形变量进行了相关性比较，其中升降轨典型特征点形变量线性拟合R2为0.978 4，说明了D-InSAR技术在地震测量中有较大的应用潜力。采用D-InSAR技术降维获得的参数输入Okada模型反演位移场特征，反演结果与已有CMT和USGS官方公布结果一致性较好，Okada模型反演参数显示，地震位于旗山断层和潮州断层的断裂带上，深度为14.59 km，发震机制为走滑逆冲断层，与D-InSAR反演位移场分析结果基本一致。研究表明，利用D-InSAR可以有效进行地震同震位移场的反演，而在没有对应区域测量验证数据时，可以采用升降轨技术验证研究区同震位移场的反演结果，也可采用Okada模型的反演结果定量和定性分析D-InSAR技术获得的位移场有效性。