引言

页岩气资源是我国重要的战略能源之一，四川盆地有丰富的页岩气资源，然而大规模水力压裂开采使得开采区周围的地震数量急剧增加（Lei et al，2017，2019b；Meng et al，2019），一些较大地震的发生给开采区带来了新的灾害问题。因此，在保证能源生产的同时，也要保证开采区免受人工诱发地震的影响。近年来，诱发地震在全球都是一个重要科学问题（Jackson et al，2014；Norris et al，2016；Keranen，Weingarten，2018）。四川盆地由于特殊的地质构造条件，这个问题更加突出，一些较大的地震事件被怀疑是工业开采所诱发，甚至可能是全球最大的诱发事件（Lei et al，2019a；Schultz et al，2020）。通常的诱发地震主要是依靠地震观测网络监测，即通过地震学的方法进行诱发地震的识别和定位，找出潜在的活化断层，评估其地震风险（Zeng et al，2014；易桂喜等，2019；Liu，Zahradník，2020）。通过干涉合成孔径雷达（interferometric synthetic aperture radar，缩写为InSAR）形变手段能否在诱发地震活动区监测到有意义的地表形变，并进一步解释诱发地震的发生机理和评估诱发地震的危险性，已经成为一个新的研究方向（Bao，Eaton，2016；Barbour et al，2016）。特别是对于页岩气开采的监测，能否得到开采相关的形变信息，提供开采中流体活动和断层活化的证据，为规避诱发地震灾害风险提供新的观测手段，更值得深入探讨。Wang等（2020）以及Li等（2021）都利用ALOS-2和Sentinel-1数据观测到了长宁页岩气区块有明显的地表形变，但其重点给出的是2019年6月17日长宁M_S6.0地震与盐矿注水的关系，并未对长宁页岩气区块的地表形变特征展开分析。

鉴于此，本文基于已有研究，排除长宁地震对研究区形变信息的干扰，利用永久散射体（persistent scatterer，缩写为PS）时间序列方法对长宁地震前所有可用的哨兵数据进行处理，并对地表形变特征进行系统性分析，以解释开采活动造成地面形变的基本特征及其机制。

1.   四川盆地南部页岩气开采与地震活动监测

1.1   四川盆地的构造活动

四川盆地位于青藏高原、秦岭与东扬子丘陵之间，是一个菱形的山间盆地。该盆地的形成受控于东北边界的米仓山—大巴山逆断裂带、西北边界的龙门山逆断裂带、西南边界左旋走滑的鲜水河—小江断裂系和东南边界的华蓥山褶皱带（图1a，1b）。这些构造形成于中生代（Burchfiel et al，1995），新生代以来随着青藏高原向东扩张又被重新活化（Burchfiel et al，1995；Xu，Kamp，2000；Wang et al，2014；Yang et al，2017；Tian et al，2018）。强震记录表明，四川盆地西部和西南边界的现今变形仍然非常强烈，近期曾发生1974年M_S7.1大关地震、2008年M_W7.9汶川地震、2013年M_W6.6芦山地震等强破坏性地震。相对而言，盆地东部边界和东北边界的现今地震活动和构造变形较为微弱（Wang，Shen，2020）。四川盆地南部长宁页岩气区块发育有建武向斜，该向斜开阔舒缓，其两翼地层倾角为5°—10°，核部地层近水平，向北与长宁背斜相连（图1c）。

图  1  四川盆地构造背景及地震活动性时空分布

（a） 四川盆地地震构造和1970年之后的地震分布，地震目录数据来源于国家地震科学数据中心（2021），GPS数据来源于Wang和Shen （2020）

Figure  1.  Tectonic settings of Sichuan basin and spatio-temporal distribution of its seismicity

（a） Seismotectonics of Sichuan basin and earthquake distribution since 1970，where earthquake catalog from National Earthquake Data Center （2021），GPS data from Wang and Shen （2020）

下载: 全尺寸图片 幻灯片

图  1  四川盆地构造背景及地震活动性时空分布

（b） 2009年以来威远和长宁地区地震的空间分布及M4以上地震的震源机制解（GCMT，2021），其中红色震源机制球代表疑似采盐注水引起的地震事件，蓝色代表疑似水力压裂引起的事件，黑色代表天然地震；（c） 长宁页岩气开采区块附近褶皱迹线和地表断层分布；（d） 2009年以来长宁区块M＞0地震的时间演化图

Figure  1.  Tectonic settings of Sichuan basin and spatio-temporal distribution of its seismicity

（b） Spatial distribution of earthquakes in Weiyuan and Changning areas since 2009 and focal mechanisms of earthquakes larger than M4 （GCMT，2021），where red beach balls represent seismic events suspected to be caused by water injection and salt extraction，blue ones represent events suspected to be caused by hydraulic fracturing，and black ones represent natural earthquakes；（c） Fold traces and surface faults around Changning shale gas block. Solid black lines are Changning anticline and Jianwu syncline；（d） Evolution of M＞0 earthquakes in Changning block since 2009

下载: 全尺寸图片 幻灯片

1.2   页岩气开采与地震活动

页岩气开采通常采用水力压裂技术，通过向页岩层中注入数十个MPa的水、砂和其它化学物质组成的压裂液和支撑剂等，使得致密的页岩层在某些方向产生裂隙，把蕴藏在裂缝中的天然气释放出来（Schultz et al，2020）。在脆性岩石被压裂的过程中，随着裂缝打开通常伴有微震的产生，震级通常＜0，这种被称为“压裂微震”的地震活动并不会造成灾害，相反会给压裂操作提供有用的压裂效果信息（Lei et al，2020）。水力压裂作业也可能使得附近不同规模的先存断层活化并在一定应力条件下发生黏滑，从而诱发灾害性地震，这是诱发地震研究的重点。在美国中部的俄亥俄州和俄克拉荷马州以及加拿大西部的页岩气压裂开采作业区域，压裂致震的情况较为普遍，一些地震学的观测证据表明水力压裂过程确实会引发有感地震（Atkinson et al，2016；Bao，Eaton，2016；Schultz et al，2020）。

四川盆地多期构造沉降和造山运动使得该地区富有油、气、煤炭和其它矿产资源，是中国页岩气储量最大的地区（蒲泊伶等，2010；Jiang et al，2016）。2014年之前，整个四川盆地的非天然地震活动主要受长期废水回注（Lei et al，2008，2013）和井盐开采注水的影响（阮祥等，2008；朱航，何畅，2014；Sun et al，2017）。自2014年开始，四川盆地内先后开发了威远、长宁、富顺永川和涪陵焦石坝等页岩气开采区块（董大忠等，2014），为期几个月的短期大量高压注液可能是该地区高水平小地震活动的主要原因。不同构造环境的地震活动性不同，其中长宁区块的诱发地震活动水平较其它区块更高。长宁区块内的上罗页岩气区块位于建武向斜（图1c），目标页岩层为志留系龙马溪组页岩，埋藏深度为2.3—3.0 km （Lei et al，2020）。我们统计了长宁区块内2009—2020年间的地震分布情况，如图1d所示，可以看出：在2014年以前，研究区的地震活动水平较低；但自从规模化水力压裂活动以来，地震活动频率显著上升，其空间分布与压裂井分布有很高的相关性（图1b和1d）。长宁页岩气区块内2018年以前未记录到超过M5的地震，2018年12月16日四川省宜宾市兴文县发生M_L5.7地震，震源矩心深度为3 km左右（Lei et al，2019b），这次地震造成了多人受伤。Lei等（2019b）认为这次地震可能是由长宁上罗页岩气区块水力压裂操作引起，可能是迄今为止全球最大的页岩气水力压裂诱发的地震活动（Schultz et al，2020）。兴文地震源于一个NNW向的单向破裂，长宁背斜的NW向伸展构造可能起到了阻止该破裂扩展的作用，并有可能促进2019年6月17日长宁M_S6.0地震的发生（Lei et al，2019a，2020）。一些研究认为长宁M_S6.0地震及其后的震群活动可能由井盐开采诱发（Lei et al，2019a；Liu，Zahradník，2020；Wang et al，2020；Li et al，2021），而非页岩气开采诱发的结果。与页岩气短期压裂注水不同，采盐时从注水井向地下2.7—3.0 km深处长期注入淡水，溶解盐矿并从对应的抽水井抽回盐水。对于长宁地震而言，注水井底部很可能与InSAR反演推断出来的发震断层面接近相交，而这种对流井技术可能通过注水提高了断层带内的孔隙压力，从而诱发了长宁地震（Wang et al，2020；Li et al，2021）。本研究将主要讨论页岩气开采对诱发地震活动的影响，两种不同作业区域大致以28.3°N为界，空间上基本没有重叠，易于分析不同诱发源及其活动机制。

目前对于水力压裂引起诱发地震活动的机制主要有两种：其一是压裂井与断层之间存在流体压力扩散通道，流体扩散作用增大了断层上的孔隙压，从而导致先存断层发生滑动（Zhang et al，2013；Bao，Eaton，2016；Schultz et al，2017；Goebel，Brodsky，2018）；其二是压裂井与断层之间无流体压力扩散通道，压裂井附近的孔隙压力上升，造成周边区域的岩石变形，进而改变周边区域内与压裂井不连通的先存断层上的剪应力和正应力（Bao，Eaton，2016；Deng et al，2016；Goebel，Brodsky，2018）。当剪应力大于摩擦系数乘以有效正应力时，地震就会随之发生，这一机理一般称为孔隙弹性效应。

1.3   诱发地震InSAR形变监测

目前对四川盆地页岩气诱发地震活动的研究主要使用地震学方法（Lei et al，2019b；Tan et al，2020）。四川盆地的地质情况比北美中西部复杂得多，地震学观测结果受限于台站分布和数据质量，往往不能给出较为确切的反演结果（张捷等，2021）。与天然地震观测一样，监测地表形变也是获取中型或大型诱发地震震源断层破裂信息的有效手段。水力压裂采用多段压裂，每段向地下注入的压裂液为大约1 800 m³的滑溜水和约100 t的支撑剂（任勇等，2015）。然而这些压裂液不是百分之百返排到地表，返排率通常为15%—30% （韩慧芬等，2017），残留在地下的压裂液及其扩散过程可能会引起地层体应变变化或者断层附近孔隙压的增加，造成地表形变并诱发地震。因此，本文首先使用常规InSAR技术分析单时相长波ALOS-2雷达卫星数据，然后利用永久散射体干涉雷达技术（PS-InSAR）分析Sentinel-1多时相雷达卫星数据，对ALOS-2的结果进行验证并进行更高精度的形变分析，进而结合当地井位数据分析该地区地表发生形变的原因。

日本宇航局（Japan Aerospace Exploration Agency，缩写为JAXA）的L波段ALOS-2数据（波长为24.2 cm，重访周期为14 d），相比欧州航天局（European Space Agency，缩写为ESA）的Sentinel-1数据（波长为5.6 cm，最短重访周期为6 d），能够在较长的时间跨度内保持较高的干涉相干性，获得足够信噪比的空间连续地表形变场。由于ALOS-2卫星在本文研究区获得的观测数据很少，所以只能通过常规差分干涉处理，尚无法形成时间序列。四川盆地内植被茂盛，湿度大，雷达电磁波的散射特性随时间的推移变化较快，因此会引起时间去相干现象（Zebker，Villasenor，1992）。而ALOS-2数据能在很大程度上克服该问题，是比较理想的观测数据。常规InSAR方法处理Sentinel-1数据由于相干性的降低而很难得到足够信噪比的信号，因此，我们利用ALOS-2升轨扫描模式合成孔径雷达数据（ScanSAR，地面覆盖范围350 km×350 km）和降轨条带模式合成孔径雷达数据（Stripmap，地面覆盖范围70 km×125 km）对长宁及其周边地区展开形变观测，期望获得研究区在观测时段内完整连续的地表形变场，以便评估页岩气开采过程中可能造成形变分布的总体特征。使用Sentinel-1数据虽然无法通过单幅干涉对有效地获得较长时段内连续的InSAR形变场，但是通过永久散射体技术对Sentinel-1数据进行时间序列处理，仍然有望获得地面一些相位稳定目标（离散点）的形变信息，可以用来验证由ALOS-2数据获得的地表形变信息。此外，关于不同轨道Sentinel-1数据长时间序列的InSAR分析，有助于我们寻找可能存在的隐伏断层，对诱发地震机制的研究及其灾害预防具有重要意义。对于开采部门来说，寻找隐伏断层也有利于避免生产中套管变形等问题带来的经济损失（陈朝伟等，2019）。

大地测量手段越来越多地被应用于长期的形变监测中。Vasco等（2010）使用基于卫星的观测手段监测注入地下大量二氧化碳所引起的地表形变，并将地表形变监测和数值模拟相结合来约束地质力学参数。Comola等（2016）以哈萨克斯坦腾格兹水库为例，建立了油藏的三维地质力学模型，并与InSAR大地测量手段结合起来，推断该地区的地质力学特征。Shirzaei等（2019）开发了基于地表形变和废水回注数据的联合反演方法，反演了目标地层的水力扩散系数，证明了InSAR手段应用于约束地层水文地质参数方面的潜力。Deng等（2020）在美国德克萨斯西部选取了三个地点，利用InSAR时间序列手段测量地表形变，并研究了该地区的地震活动与石油和天然气生产之间的关系。目前也有研究人员对页岩气生产诱发地面形变场进行了初步分析（Kubanek et al，2018；Jordan et al，2019），但对于四川盆地页岩气区块地表形变场的特征缺乏系统性分析和物理机制方面的研究。

2.   实验数据与处理方法

本研究首先基于日本宇航局的L波段ALOS-2 PALSAR数据，使用常规InSAR技术监测长宁页岩气区块两三年间的累积地表形变分布。与C波段雷达数据相比，即使在四川盆地这样地形起伏较大、湿度大、植被茂盛的地区，InSAR数据也能保持良好的相干性。其次，我们使用2015—2019年间获取的升降轨道Sentinel-1数据进行PS-InSAR时间序列分析，以克服常规InSAR技术受时空失相干、地形误差、大气延迟效应等误差的影响，获得较常规InSAR方法更高精度、更长时间序列的形变信息。对于ALOS-2数据来说，由于长波信号会受到电离层中电子密度变化的巨大影响，对于不同观测时间和空间，电离层误差可能达到厘米级甚至是米级（Sandwell et al，2008），这在汶川地震期间ALOS-1卫星的InSAR观测中表现得最为明显（Shen et al，2009）。本研究中使用的ALOS-2数据同样也受到这种误差的严重影响，因此必须考虑这种误差的影响才能获得有效的形变信息。此外，四川盆地常年水汽较重，而且其时空变化非常快，这是影响InSAR观测误差的主要因素，特别是在单时相InSAR观测中对结果精度影响较大，因此我们也尝试使用大气模型辅助对该误差进行校正。

2.1   ALOS-2数据的InSAR处理和误差校正过程

常规InSAR干涉图由同一地区获取的两幅配准后的单视复数（single look complex，缩写为SLC）图像进行复共轭相乘得到，反映的是地面上同一位置在两次成像间隔的相位差，去除地形相位后每一个像素的差分干涉相位由六部分组成：

式中：Δϕ_topo为所使用数字高程模型（digital elevation model，缩写为DEM）误差引起的剩余高程相位；ϕ_tropo为大气对流层延迟相位，是在两次观测时由于对流层差异引起的相位误差；ϕ_ion为两次成像时电离层电子密度变化引入的相位；ϕ_def为两次成像期间内地表移动引起的雷达与地面之间的视线向（light of sight，缩写为LOS）变化；Δϕ_orb为剩余轨道误差相位，对于现代卫星基本可以忽略；随机噪声ϕ_n为失相干噪声。对于地表形变测量来说，只有ϕ_def是我们希望得到的信息，其它相位均为误差项，尤其是大气电离层ϕ_ion和对流层ϕ_tropo的影响，足以把厘米级的形变掩盖。对于常规InSAR处理，我们可以利用误差相位分量的不同空间尺度特点或者外部辅助数据去除它们。对于L波段ALOS-2数据，四川盆地所处的纬度地区受电离层电子密度变化影响较大，需要特别考虑电离层误差的影响。

利用ALOS-2不同成像观测模式，并选取升降轨两种不同成像几何的合成孔径雷达（synthetic aperture radar，缩写为 SAR）数据进行干涉处理并消除地形相位的影响，得到常规InSAR的差分干涉结果，如图2所示。图2a显示整个图像范围均分布着复杂的不规则干涉条纹，由于ALOS-2卫星较ALOS-1卫星在轨道控制方面有显著改善，图像中的轨道误差所产生的影响很小，所以这种大尺度不规则条纹与形变无关，也无法用线性或者非线性轨道模型拟合，主要为电离层相位延迟误差（Liang，Fielding，2017）。图2b显示了布满整个图像的规则干涉条纹，后文基于电离层物理模型的相位分析表明，它们也主要是电离层相位误差。另外，图2a和2b中大尺度的干涉条纹也可能叠加了长波对流层延迟误差。由于Δϕ_topo与垂直基线有关，本文采用1″的航天飞机雷达地形测绘使命（Shuttle Radar Topography Mission，缩写为SRTM）的DEM数据在小基线情况下（表1）引入的误差相对较小，失相干噪声ϕ_n对最后结果的影响较小，在数据相干性较好时可以忽略不计（Hanssen，2001）。对形变结果产生较大影响且较难校正的是大气延迟引起的相位误差，包括对流层延迟和上述电离层延迟两种。对于四川盆地来说，成像期间水汽分布剧烈变化会造成较大的对流层延迟误差，对所有波段的雷达数据均有较大影响。当使用长波雷达数据时，电离层效应在不同纬度、不同时间对合成孔径雷达信号的干扰程度不同，而在四川地区受电离层影响较大（Shen et al，2009，Gomba et al，2016）。对于C波段Sentinel-1数据，电离层在低纬度和高纬度如极地地区对合成孔径雷达信号有严重影响，但本研究所考虑的纬度地区相对小，因此下面重点介绍大气相位延迟的建模和去除方法。

图  2  功率谱滤波后的ALOS-2差分干涉相位结果

（a） 升轨ScanSAR模式差分干涉图，成像时间为2016年8月7日至2019年9月15日；（b） 降轨StripMap模式差分干涉图，成像时间为2017年6月12日至2019年7月8日

Figure  2.  ALOS-2 differential interferometry phases with power spectrum filtering

（a） Differential interferogram in ascending orbit with ScanSAR mode，imaging from August 7，2016 to September 15，2019； （b） Differential interferogram in descending orbit with StripMap mode，imaging from June 12，2017 to July 8，2019

下载: 全尺寸图片 幻灯片

表  1  本研究所用ALOS-2 PALSAR数据的参数

Table  1.  ALOS-2 PALSAR data parameters used in this study

主景日期	从景日期	入射角/°	卫星方向角/°	轨道	垂直基线/m
2016-08-07	2019-09-15	39	−10.2	升轨（146）	75.7
2017-06-12	2019-07-08	40	−169.9	降轨（37）	−27.6

下载: 导出CSV 

| 显示表格

2.1.1   去除电离层相位

电离层是大气中的电离部分，卫星雷达传感器发射的电磁波信号两次穿过电离层，会对SAR信号的相位产生显著影响（Gomba et al，2016；王楠等，2017）。对于ALOS-2卫星数据来说，其L波段信号的频率较Sentinel-1卫星的C波段信号低，电离层对低频雷达系统干涉相位的影响不能忽略。因此，在数据处理中需要估计电离层相位的影响，并将其从干涉图中去除。如果考虑电离层的影响，InSAR干涉相位的组成可表示为

式中：f₀为载波频率；ΔTEC为LOS方向上两次成像的总电子含量（total electron content，缩写为TEC）差值；K为常数，取40.28 m³/s²；c为光速（Liang et al，2018）。式（3）表示干涉相位 ［ 式（2） ］ 中的非频散相位，与频率成正比，主要是剩余地形、形变、对流层和剩余轨道相位，式（4）表示由电离层变化引起的频散相位，与频率成反比。

基于电离层相位的频散特性，可以将SAR距离向信号频谱分割为不同频段来估计电离层相位，即距离向信号频谱分割法（Brcic et al，2010；Rosen et al，2010；Gomba et al，2016；Liang et al，2018）。本文对ALOS-2数据应用距离向频谱分割法得到差分电离层相位（Liang，Fielding，2017），该方法将需要作干涉的两幅SAR图像分割为两个不同距离向中心频率的子带图像：上子带图像（频率为f_u）和下子带图像（频率为f_l）。对两幅上子带图像干涉得到干涉相位Δϕ_u，对两幅下子带图像干涉得到干涉相位Δϕ_l，在不考虑相位展开误差和其它相位误差的情况下，将两子带图像的干涉相位代入

便可以计算得到电离层相位ϕ_ion （Liang et al，2018；Zhang et al，2018），如图3所示。对比原始差分干涉图（图2），可知：形变场中的长波相位与电离层延迟相位非常相似，由此可以确定图2中的误差主要为电离层相位，去除该相位影响的结果如图3c和图3d所示。去除电离层相位后的差分干涉图还主要有对流层延迟误差相位，下面继续消除这部分误差。

图  3  利用距离向信号频谱分割法估计电离层延迟相位

（a） 升轨ScanSAR模式图像的电离层相位估计结果；（b） 降轨StripMap模式图像的电离层相位估计结果；（c） 升轨ScanSAR模式图像去除电离层相位的缠绕相位；（d） 降轨StripMap模式图像去除电离层相位的缠绕相位

Figure  3.  Estimation of the ionospheric phase delay by range split-spectrum method

（a） Ionospheric phase estimated from ascending orbit image with ScanSAR mode；（b） Ionospheric phase estimated from descending orbit image with StripMap mode；（c） Wrapped phase estimated from ascending orbit image with ScanSAR mode after removing ionospheric phase；（d） Wrapped phase estimated from descending orbit image with StripMap mode after removing ionospheric phase

下载: 全尺寸图片 幻灯片

2.1.2   去除对流层相位和长波趋势相位

大气中温度、压力和相对湿度在时间和空间上的变化造成了干涉相位的对流层延迟误差。20%的相对湿度变化即可导致干涉图上的信号产生高达10—14 cm的形变测量误差（Zebker et al，1997；孙广通等，2011），而这个误差可能会掩盖干涉图上的地表微小形变，对于地壳形变场的监测非常不利。根据物理机制的不同可将对流层延迟分为两部分：天顶干延迟（又称垂直分层延迟）和天顶湿延迟（又称水平湍流混合延迟）（Parker et al，2015）。垂直分层延迟相位与高程呈反比关系，即随着高度的增加，水汽含量降低，大气延迟相位减小；水平湍流混合延迟主要是由水汽的含量变化和随机分布引起的水平紊乱延迟（Lohman，Simons，2005；Ebmeier et al，2013）。在平坦地区，水平湍流延迟是主要因素，此时垂直分层延迟受压力主导，在空间上表现得较为平滑；在地形起伏明显的多山区域，由海拔变化引起的垂直分层延迟是主要因素（Doin et al，2009）。

式中：ϕ_tropo0代表从高度h₁到对流层顶部h_top的对流层延迟相位，与大气折射率N有关；N_hydr代表静力学平衡折射率；N_wet代表非静力学折射率；θ为影像获取时的入射角；ϕ_tropo为获取主辅影像时大气折射率的时空变化所导致的对流层延迟相位，其中上标s表示从图像，上标m表示主图像（Bekaert et al，2015）。

目前针对对流层改正的方法有三大类：① 基于外部辅助数据法，外部数据源有GPS数据、大气模型（ERA5，ERA-I，WRF等）、多光谱观测数据（MODIS，MERIS等）；② 基于经验关系的干涉相位法，包括线性相关模型和幂律模型；③ 基于数据本身的校正技术，主要有相位叠加平均法（stacking）和时序InSAR方法等（周洪月，2018）。不同的改正方法对于干、湿延迟分量的敏感度不同，都有各自的优缺点。MERIS和MODIS光谱仪数据只能提供有光照条件下的数据，对云的敏感度很大，并且只能估计湿延迟，而线性相关模型和幂律模型只能估计受地形起伏影响大的干延迟。长宁地区属于多山区域，且水汽含量较大，需要去除对流层的影响。下文将分析如何从干涉相位中去除对流层延迟影响。

由欧洲中尺度天气预报中心（European Centre for Medium-Range Weather Forecasts，缩写为ECMWF）提供的最新一代全球气象再分析产品（ECMWF Re-analysis v5，缩写为 ERA5）既可以估计干延迟分量，又可以估计湿延迟分量，相比于之前的ERA-I模型来说有更好的时空分辨率（Hu，Mallorquí，2019）。本文中对流层延迟改正使用外部辅助数据（ERA5 模型），利用TRAIN软件实现。对流层校正后的图像如图4a所示，可见：在本例中ERA5模型对整个图像大尺度对流层延迟相位的改善明显，但是对小尺度对流层延迟相位作用有限，对我们的研究区仅有几个弧度的改正。

图  4  大气对流层延迟校正及长波趋势相位校正

（a，b） 基于ERA5模型去除升、降轨对流层延迟误差后的图像；（c，d） 长波趋势相位校正后的图像，其中蓝色方框表示本研究的形变区

Figure  4.  Atmospheric troposphere delay correction and long wavelength trend phase correction

（a，b） Images after removing tropospheric delay errors in ascending and descending orbits based on the ERA5 model， respectively；（c，d） Images after long wavelength trend phase correction，where blue boxes denote the deformation area of this study

下载: 全尺寸图片 幻灯片

在四川盆地这种水汽变化较快的地区，ERA5模型可能无法有效地消除长波对流层误差，同时用大气模型进行对流层改正时也可能会引入额外的长波误差，我们将这部分误差称为长波趋势相位误差，如图4b的情况。前文提及的轨道误差在新一代卫星中影响很小，可忽略不计。因此，我们用一阶线性模型对上述校正后的图像进行趋势相位误差校正，去除剩余的长波误差，得到最终校正后的差分干涉图（图4c，4d）。虽然图4c中仍然能看到一些比较明显的误差干涉条纹，其主要来源于大范围（350 km×350 km）对流层延迟误差，但是在我们关注的页岩气开采区块内（图中矩形框区域），这种误差的影响非常有限，不影响对结果的分析。图4d中，通过上述误差校正后已经看不到明显的误差条纹，对比图3d，图像中无论是长宁地震的形变，还是页岩气区块的形变，其信噪比均较高。除了长宁地震，一些较大的地震事件（如M_W＞4.0）的同震形变场也会对InSAR形变场产生影响，但其位置已知，且形变范围有限，比较容易识别，不会影响页岩气区域的形变分析。因此，我们将图4c和图4d的结果作为最终的ALOS-2形变观测结果。

通过上述过程的多步校正，特别是针对ALOS-2数据进行的电离层延迟误差校正，其结果显示利用长波SAR数据进行四川盆地或者龙门山地区的InSAR形变信息提取有重要的应用价值。同时，L波段的ALOS-1或者ALOS-2数据能在长时间内保持较高的相干性，对于复杂环境下，如复杂地面植被和地形条件等，获得高信噪比的InSAR观测结果非常有利。但是L波段的观测数据较少，只有日本宇航局的ALOS-2可以使用，该数据没有针对四川地区的系统性数据获取计划，因此还无法形成比较好的观测时间序列，阻碍了长波InSAR观测的应用。为了获得比较连续的形变时间序列，只能利用欧洲航天局的Sentinel-1数据，具体将在下节中进一步分析。

2.2   Sentinel-1数据的PS-InSAR处理

在ALOS-2卫星观测到的常规干涉形变场中，我们能够清晰地看到长宁页岩气区块存在显著的地表形变（图4），且不同时段的形变存在一定差异，而该区域页岩气开采是唯一可能的形变来源。为了验证ALOS-2数据的观测结果，并获得更高精度的形变信息，我们使用欧洲航天局的Sentinel-1雷达卫星数据，通过时间序列分析技术得到地表形变场。为了避免ALOS-2升降轨数据获取时间的差异，我们选择相同时段的Sentinel-1升降轨数据进行分析。本文选用的ALOS-2数据时间间隔较长，因为使用传统单时相InSAR技术，不但会受到失相干的影响，而且还会受到前述电离层和对流层的干扰。特别是四川盆地浓厚大气所导致的对流层延迟误差的影响，虽然可以通过ERA5天气模型等辅助数据进行校正，但是两者的时空分辨率差异使得校正效果有限（图4）。为了得到更加准确的结果，我们对Sentinel-1数据进行PS-InSAR处理，从经过配准的超过四年累积的Sentinel-1单视复数图像（2015年2月至2019年5月）上提取相位稳定的永久散射体目标进行分析。该时段页岩气开采处于快速发展期，并且形变场中不包括长宁地震的相关形变，能够反映页岩气开采产生的形变特征。

区别于最早的基于幅值和先验形变模型识别稳定点目标的PS-InSAR处理方法（Ferretti et al，2000，2001），我们基于斯坦福大学开发的StaMPS软件（Hooper et al，2004）及其相应的相位分析算法进行点目标的识别和形变分析。该软件基于干涉相位的空间相关性来识别PS点，不需要假设的先验形变模型，适合对非线性的时间变化过程进行监测。该方法中永久散射体的识别不取决于雷达图像幅值，在雷达幅值较低的区域也能发现一定量的稳定相位目标，因此也适用于非城市地区的形变监测。本研究的页岩气开采区块位于四川盆地内部的乡村地区，地面上找不到足够的高幅值散射体，但是用StaMPS方法依然能够识别出足够的PS点进行相位分析。同时，与稳定的构造运动不同，页岩气开采作业的时空变化较大，其造成的地面形变也非稳态变化，无法用简单的形变时间演化模型来表示，但使用基于相位相关性的算法仍然可以得到合适的形变分析结果。

StaMPS方法的主要优势在于PS点的识别，通过定义时间相干性参数γ_x，在一定空间相关性尺度范围内迭代筛选PS点。

式（8）为时间相干性的表达式，其中不含幅值相关信息，所以该方法识别出的PS点目标理论上与幅值无关。另外，式（8）中也未涉及形变的时间模型，因此PS点目标的选取与形变量值高低无关，如果其相位是稳定的，这个点依然可以作为PS点，这也是该算法可以用于非城市地区观测的原因，这有利于时间非稳态信息的提取。

式中，N表示干涉图的数量，ϕ_{int， x， i}表示利用SRTM的DEM地形数据消除地形相位后的干涉相位，$\hat{\phi }_{{\rm{int}}， x， i} $为给定空间相关性范围内对$ {\phi }_{{\rm{int}}， x， i} $进行频率域带通滤波的结果，$\Delta{\hat{\phi }}_{\theta ， x， i}^{{\rm{nc}}}$为空间非相关视角误差相位的最小二乘估计。

筛选出PS点后，从缠绕的干涉相位减去式（8）中估计得到的空间非相关误差，使得相邻PS点之间的相位差小于Π，此时可以对剩余干涉相位进行展开。缠绕的干涉相位为

式中：${\phi }_{{\rm{def}}， x， i}$为两次成像期间的PS点目标沿着LOS方向移动所导致的相位变化；${\phi }_{{\rm{atm}}， x， i}$为两次成像期间的大气分布变化所造成的延迟误差（对于C波段主要是对流层延迟）；$\Delta {\phi }_{{\rm{orb}}， x， i}$为轨道不精确所导致的相位残差；$\Delta {\phi }_{\theta ， x， i}^{{\rm{corr}}}$为空间相关性视角误差；${\phi }_{{\rm{noise}}， x， i}$为时空非相关噪声，由散射机制变化、热噪声和配准误差等引入；x，i代表第x幅干涉图的第i个像素；W{···}代表缠绕相位（Hooper et al，2007）。其中：时空非相关误差项${\phi }_{{\rm{noise}}， x， i}$可以当作噪声处理；空间相关而时间不相关的误差项${\phi }_{{\rm{atm}}， x， i}＋\Delta {\phi }_{{\rm{orb}}， x， i}＋\Delta {\phi }_{\theta ， x， i}^{{\rm{corr}}} $可通过对解缠绕的数据首先进行时间上的高通滤波，再进行空间上的低通滤波而估计；剩余的相位即为空间相关且时间上也相关的形变相位，是我们最终要估计的形变观测结果${\phi }_{{\rm{def}}， x， i} $；大气对流层相位${\phi }_{{\rm{atm}}， x， i} $也可以使用前文给出的ERA5大气模型进行估计，因为C波段数据在四川盆地地区电离层误差相对小，相位模型中未予以考虑。对于最终的平均速度场估计，仍然可能存在一定的长波误差，与前文提到的长波趋势相位类似，主要来自于长波对流层延迟误差。由于本研究涉及的形变区均远小于SAR成像范围，因此这种长波误差可以通过线性坡面进行拟合和消除。

本文搜集了跨越研究区升降两条轨道（A055和D164）的Sentinel-1数据进行PS-InSAR分析，时间跨度为2015年2月至2019年5月底，截止于长宁地震前两周左右。从图像中截选出长宁页岩气开采区及其周围的数据进行单一主景干涉处理（地面覆盖见图5），升轨形成84幅干涉图，降轨形成89幅干涉图，利用StaMPS软件分别对升降轨道作PS-InSAR时序处理。升轨以2017年11月10日获取的数据为主景，最大垂直基线为122 m，最大时间基线为1 116 d；降轨以2017年2月26日获取的数据为主景数据，最大垂直基线为205 m，最大时间基线为972 d。该研究区升降轨道分别识别出20万4 088个和44万4 102个PS点，经过误差去除并进行PS点插值后，得到观测时段内的平均线性速度场（图5），该结果揭示了长宁上罗页岩气开采区在观测时段内的一些典型变形特征。

图  5  基于长宁页岩气开发区Sentinel-1雷达数据得到的升轨（a）和降轨（b） PS-InSAR线性速度场

图中黑点表示该区块中已知页岩气井的空间分布

Figure  5.  PS-InSAR linear velocity field in ascending （a） and descending （b） orbits obtained from Sentinel-1 radar data around Changning shale gas block

The black dots indicate spatial distribution of known shale gas wells in Changning shale gas block

下载: 全尺寸图片 幻灯片

3.   观测结果分析

我们使用两种不同的数据处理方法，分析了两种不同的卫星数据，最终得到了关于页岩气开采区块的一些基本变形特征。这些形变结果虽然受到四川盆地复杂地面观测条件、对流层水汽分布变化等因素的影响，但是通过一些有效的技术手段我们仍然能够得到足够信噪比的观测结果，这一节将结合两种卫星的观测结果分析这种变形的基本特征及其形成机制。

3.1   ALOS-2数据的InSAR结果分析

将图4c和4d中校正完成的干涉图像放大到本研究关注的长宁页岩气开采区块，可以看到长宁页岩气开采区附近不同轨道ALOS-2数据上的地面形变场及其南北向剖面图，如图6所示。图中28°20′N附近的形变场（虚线框内）为2019年6月17日长宁M_W5.8地震的同震形变场（Li et al，2021）；AA′和BB′剖面附近为长宁—上罗页岩气开采区块，整个区块分布有密集的页岩气开采井（图中黑点所示）。在这些开采井附近，ALOS-2形变场显示了显著的地面形变（LOS方向）。由于ALOS-2数据成像时间的限制，升降轨道InSAR观测时间不一致，图6a升轨的观测时间较图6b降轨的观测时间长将近一年（表1），因此两者显示的形变场有一定差异，观测时间较长的升轨数据（图6a）显示了更显著的地面抬升信号（蓝色正值区，假设地面发生垂直形变），这与页岩气水力压裂开采中大量向页岩层注液有关。在ALOS-2的观测期间内，研究区范围内发生过多次M_W＞4.0地震，但是距离本文研究的形变中心区有一定距离，对页岩气形变影响不大，且在观测时间较长的条件下，对于速度场来说，地震形变的影响被削弱，因此本文分析不考虑这些地震带来的影响。另外，开采区的这些井并非同时工作，有些井在水力压裂阶段，而有些井是在采气生产阶段，因此造成了地面形变场的复杂性。在两个观测时段内的一个有趣现象是，跨越形变区相同位置的剖面显示，地面隆升的幅值和峰值位置均不同：对于升轨数据，LOS向形变峰值达到了约120 mm，出现在10 km的位置（图6c）；对于降轨数据，LOS向形变峰值仅约20 mm，出现在13 km的位置（图6c）。假设地面仅为纯垂直形变，这说明在2016年8月至2017年6月这个时段内，压裂生产大规模开展，长宁页岩气区块压裂活动比较剧烈，造成了地面快速大幅度隆升；2017年6月至2019年7月这段时间内，水力压裂的重心已经往南偏移，而且开采的幅度显著下降，对应的地面隆升峰值也大幅降低，仅约20 mm；2017年下半年开始，由于大多数气井压裂工作已完成并转入生产阶段，因此在这段时间内原本上升的区域可能逐渐表现为地面沉降，或者LOS方向上远离降轨卫星的地面水平运动，亦或两种运动兼而有之（图6c中BB′剖面上0—10 km位置的负值区域）。这种形变的变化过程与石油公司公开报道的生产进度以及我们野外实地调查都一致，表明形变剖面很好地反映了这种生产过程的转变。

图  6  长宁页岩气区块近场ALOS-2干涉形变位移场及剖面图

长宁地震序列和水力压裂致震的震源机制解来源于Lei等（2019a），其余2015年和2017年的两次地震的震源机制解来源于GCMT （2021）；白色虚线框表示2019年6月17日长宁地震形变场（a） 升轨ScanSAR模式图像得到的页岩气区块形变位移场，成像间隔为3年；（b） 降轨StripMap模式图像得到的页岩气区块形变位移场，成像间隔为2年；（c） 图（a）中的AA′和图（b）中的BB′形变位移场剖面，两剖面在同一位置

Figure  6.  ALOS-2 interferometry displacement field and profiles in the near field of Changning shale gas block

The focal mechanism solution of the Changning earthquake sequence and hydraulic fracturing are from Lei et al （2019a），the remaining two earthquakes in 2015 and 2017 are from GCMT （2021）. The white dashed box represents the deformation field of the Changning earthquake occurred on June 17，2019 （a） Deformation displacement field of shale gas block obtained from ascending orbit data in ScanSAR mode，and the imaging interval is three years；（b） Deformation displacement field of shale gas block obtained from descending orbit data in StripMap mode，and the imaging interval is two years； （c） Deformation displacements on the profiles AA′ in Fig. （a） and BB′ in Fig. （b），and the two profiles are at the same location

下载: 全尺寸图片 幻灯片

通过对覆盖长宁地区的升降轨ALOS-2数据进行干涉处理，消除电离层误差、大气误差等因素后，得到研究区可靠的差分干涉形变场。结果显示页岩气生产过程中能产生足够大的地面形变，而且也能够通过InSAR技术有效监测到。虽然L波段ALOS-2数据的时间基线较长，但其波长较长，对地面目标的散射机制变化不及C波段敏感，因此即便在四川盆地复杂的地面条件和浓厚多变的水汽观测条件下，两三年内仍能保持较高的干涉相干性，并能从干涉图中清晰分辨研究区的形变及其演化情况。此外，基于ALOS-2的单时相InSAR形变分析也表明，页岩气开采会造成地面的快速隆升或者下降，形变的变化呈非稳态，在压裂阶段由于注入流体而使得地面隆升，在生产阶段可能由于抽取返排液或者天然气而使得地面沉降。受限于ALOS-2数据的可用性，目前尚无法做到长时间序列的InSAR形变观测，这是制约ALOS-2应用于页岩气开采监测的主要问题。

3.2   Sentinel-1数据的PS-InSAR结果分析

Sentinel-1数据的PS-InSAR分析结果（图5）为整个观测时段（2015年2月至2019年5月）的线性速度场，其中建武向斜内的形变特征与Wang等（2020）利用2016年3月至2019年6月的Sentinel-1数据进行时序分析的结果相似。本研究仅使用速度场进行分析，因为基于C波段的InSAR观测会受到四川盆地浓厚水汽的严重干扰以及季节性形变的影响（Li et al，2021），形变时间序列中位移场有很大的不确定性，而速度场是基于所有干涉数据作最小二乘反演所得，能够最大限度地降低时间随机信号，特别是大气对流层信号的干扰。从图7d的形变剖面点的发散程度可以大致估计形变速度场的观测精度，本研究中Sentinel-1的InSAR时间序列能达到的形变观测精度约为1 mm/a，因此它是页岩气开采监测的有力工具，可以为未来开采区地层的水文地质参数反演和人工诱发地震研究提供重要信息。在长宁页岩气开采区附近（图5中黑点分布范围），升降轨数据均显示了较高的信噪比，能客观地反映页岩气开采造成的地面形变。我们同样将图5的速度场结果放大到开采形变区来进一步分析（图7）。

图  7  长宁页岩气区块近场升降轨Sentinel-1数据的LOS向线性速度场

图（a）和（b）中KK′为图8a中二维地震反射剖面位置，OO′为图8b中断层破碎带的北界；图（c）中红色直线为线性速度模型拟合的各点位移。（a） 升轨速度场；（b） 降轨速度场；（c） 图（a）和图（b）中P₁，P₂，P₃和P₄点的位移时间序列；（d） 升轨数据剖面CC′和降轨数据剖面DD′上的形变速率变化

Figure  7.  Linear velocity field of Sentinel-1 data in LOS direction in the near field of Changning shale gas block

In Figs. （a） and （b） KK′ is the location of the two-dimensional seismic reflection profile in Fig. 8a，and OO′ is the north boundary of a fault fracture zone in Fig. 8b；in Fig. （c） the red straight lines are the displacements of each point by linear velocity fitting. （a） Velocity field in ascending orbit；（b） Velocity field in descending orbit；（c） The displacement time series of points P₁，P₂，P₃ and P₄ in Figs. （a） and （b），respectively；(d) Variation of deformation rates on the profile CC′ in ascending orbit and profile DD′ in descending orbit

下载: 全尺寸图片 幻灯片

图7a和7b的范围与图6a和6b相同，图中AA′，BB′，CC′和DD′这几个剖面的位置也相同，这样处理便于我们比较四次不同的观测结果。图7中所使用的SAR数据都在长宁地震之前获取，所以该结果中未包含图6中长宁地震的形变场。另外，选取Sentinel-1升降轨数据的时候尽量保证两者时间接近，以便比较不同观测角度的InSAR形变场。

剖面CC′和DD′ （图7a，7b）的形变速度场的一个主要特征是在升降两条轨道的12—13 km处都存在一个明显的峰值，约为6—8 mm/a （图7d）。这个峰值的位置恰好与图6c降轨ALOS-2数据的峰值位置相同，这反过来说明降轨ALOS-2观测到的峰值是可信的。这个峰值的物理意义在于，在Sentinel-1观测覆盖的长宁页岩气开采时段内，地面形变场以局部隆升为主，隆升主要发生在南北向剖面上10—15 km的位置，图像上显示为一条近东西向的狭长条带（图7a，7b）。在升降轨道不同的观测成像几何条件下，该条带的形变分布和量值都非常接近，因此我们可以确定该形变主要以垂直抬升为主，而两者的微小差别反映了SAR卫星观测入射角的差异（图7d）。形变剖面还反映了一个比较明显的特征是，在0—8 km的范围内存在一个明显的LOS向负值区域，该负值区以5 km的位置为中心并达到峰值约−2 mm/a。如果将图7b中的负值区与图6b中ALOS-2降轨数据的负值区进行对比，可以看到两个负值区域在0—8 km范围内几乎是重叠的，因此可以确定该负值区是真实的地面形变信息，而不是由于Sentinel-1降轨SAR数据中的噪声所致。但是对比升轨Sentinel-1数据可知，相同区域在升轨上并未出现明显的负值，而是大部分接近于零值。在图6c所示的ALOS-2升轨数据AA′剖面上0—5 km的范围内也未见明显的负值分布，这与Sentinel-1数据的观测也是一致的。综合这四种InSAR观测结果推知，图7d中DD′剖面负值出现的原因是存在一个地面水平向的运动分量，其运动方向与升轨卫星飞行方向近乎平行（N12°W向）。这样的运动使我们在升轨图像上很难观测到这个方向的形变，但在降轨条件下观测到了负值形变。出现这样的水平向地面形变的一个可能原因是页岩气开采注入的压裂液产生了NW向的扩散，这种流体的扩散作用也使得地面产生了相同方向的形变。

图7c给出了对应于图7a升轨数据的P₁点、P₂点和对应于图7b降轨数据的P₃点、P₄点的位移时间序列，可以看出：P₁点和P₃点分别为升轨和降轨数据中抬升量较大的位置，在2015—2019年这四年期间，累计抬升量达到约30 mm；P₂点和P₄点分别为升轨和降轨数据中发生沉降较大的位置，在四年期间累计沉降量达到约15 mm。

图7d的形变剖面上还有一个值得讨论的特征是，在剖面15—16 km的位置，形变速率值突然由正变负，该形变转向的位置为图7a和7b的OO′虚线所示，在图像上形变的变化也很明显，而且在升降轨上均呈现相同的特征。此外，从图7a和7b上可以看出，在OO′剖面以南区域没有页岩气开采井的分布。因此，我们推测地面形变的变化与该区域地层结构发生变化有关，需要通过长宁页岩气开采区的地层结构及二维地震反射剖面资料予以进一步分析。

图7中KK′剖面的地层结构如图8a所示，卷入建武向斜变形的沉积地层由下至上依次为震旦系、寒武系、奥陶系、志留系、二叠系、三叠系和侏罗系，缺失泥盆系和石炭系，总厚度可达8000 m，其中震旦系—中三叠系为一套由灰岩、泥岩、泥页岩、白云岩和石膏岩组成的海相沉积，而之上的上三叠系—侏罗系是以砂岩、砾岩为主的陆源碎屑沉积（图8a，8c）。页岩气开采层位为志留系的龙马溪组，为一套富含笔石化石的黑色页岩层，也是页岩气开采的主产气层（何登发等，2019）。图8b给出了长宁页岩气开采区的二维地震反射剖面，在该剖面上可以看到：建武向斜南翼发育有一条宽约2 km的断层破碎带，其在二维石油物探剖面上表现为地震同向轴的杂乱反射；断层带大致倾向南，切穿了震旦系—三叠系的全部沉积盖层。该断层破碎带的北界大概与图7a中的OO′位置对应，它的存在可能导致OO′剖面以南区域页岩气存储层的缺失。从形变场来看，该破碎带以北有密集的页岩气开采井分布，地面隆升明显，而其北界（OO′剖面）可能成为流体扩散的障碍体（图8b），这使得开采产生的大量流体扩散不到OO′剖面以南区域，因此该区域不会发生隆升，其变形行为与弹性体相似，在北侧快速隆升的过程中弹性应力场造成南侧地块的相对下降，其下降幅值明显小于北侧的隆升区。此外，我们在图6c中AA′和BB′剖面上发现，类似的特征对于ALOS-2数据也是存在的，但是由于ALOS-2单时相InSAR观测的噪声干扰以及干涉滤波操作的影响，这个特征并不明显。

图  8  跨越长宁开采区的二维地震反射剖面

（a） KK′剖面（图7）对应的地层结构和近年来长宁附近发生的较大地震事件及其震源机制解（Li et al，2021）；（b） 图7d形变剖面中15—17 km处断层破碎带附近的二维地震反射剖面，其位置见图1c；（c） 长宁地区的沉积地层序列（Li et al，2021）

Figure  8.  Two-dimensional seismic reflection profile across the Changning shale gas block

（a） The stratigraphic structure corresponding to the profile KK′ （Fig. 7） and several large-magnitude events near Changning area in recent years and their focal mechanism solutions （Li et al，2021）；（b） The two-dimensional seismic reflection profile near a fault fracture zone at 15−17 km of the deformation profile in Fig. 7d，and the location of the profile is shown in Fig. 1c；（c） Sedimentary stratigraphic sequence in Changning area （Li et al，2021）

下载: 全尺寸图片 幻灯片

综上所述，基于Sentinel-1数据的PS-InSAR时间序列观测结果表明，多时相InSAR观测技术对于页岩气形变监测是可行的，它能够提供比ALOS-2单时相InSAR观测更高的精度，也能更好地反映页岩气开采的一些形变细节特征。

4.   讨论与结论

本文基于两种SAR卫星数据，使用两种InSAR技术得到了长宁页岩气开采区块的地表形变场。针对L波段ALOS-2卫星数据提供了电离层校正、对流层校正等一系列处理方案，得到了在严重电离层和对流层干扰下长达两三年的InSAR干涉结果。该卫星数据的分析结果表明，页岩气开采能够产生非常显著的地面形变场，而且随着开采的进展，地表形变也随之变化，呈现出非稳态的形变特征。由于页岩气开采需要向页岩层注入大量的流体，而且这些流体大部分会留在地层中而无法返排到地面，因此造成地层体应变变化和地面在开采过程中局部强烈隆升，本研究中隆升量达到了120 mm。在页岩气压裂完成转入生产阶段，生产过程中驻留在地层中的流体逐步向外扩散，会使得地面下沉并产生特定方向的水平运动。

长宁页岩气区块内近年来发生了一些较大地震事件，特别是2019年1月3日珙县M_S5.3地震、2018年12月16日M_S5.7兴文地震以及更早的2017年1月27日M_S5.0筠连地震等（图8），这些地震均发生在建武向斜内，并与长宁页岩气开采区块非常接近。从形变场的角度如何认识这些地震与页岩气开采的关系，对后续诱发地震模型的分析具有重要意义。为了避免InSAR时间序列中同震形变场信号被时间序列反演过程削弱的影响，我们重新处理了ALOS-2降轨数据，并尽可能降低InSAR相位功率谱滤波器的强度，以保证这些小地震的形变信息能够被清晰保留（图9a）。除了筠连地震由于震级较小形变难于识别外，珙县地震和兴文地震的形变信号比较清楚。我们也用Sentinel-1数据处理获得了珙县地震的同震形变场，得到了间隔12天的形变场图像（图9b），但兴文地震的InSAR结果在Sentinel-1数据上失相干较重，仅能对形变的范围进行精确定位。从图9a中可以清楚地识别出三个主要的形变区A，B和C，同时也标记出了可能与珙县地震相关的N201区块的H18井的位置，这口井自2018年11月开始直到珙县地震发生时一直在进行压裂作业（Lei et al，2019b）。形变区A位于页岩气开采产生的主要形变区C的西北边缘，距离C较近。通过查阅地震目录及野外调查我们发现，形变区A与珙县地震的位置相符。为了判断形变区A主要源于水力压裂还是珙县地震弹性变形，我们利用珙县地震前后2018年12月30日和2019年1月11日的Sentinel-1降轨数据作干涉，并转换成LOS向地表形变场，如图9b所示。其形变分布清楚地揭示了一条近南北向断层的弹性形变场，且LOS向形变量值可达61 mm，而水力压裂在12天内很难达到该量值的形变，空间分布也很难达到图9b所示的形变范围，所以尽管H18井位于形变区A内，我们认为其主要是珙县地震所致。同样，形变区B主要是由兴文地震断层滑动引起，附近的开采井距离形变区均较远。需要指出的是，我们在形变图上观测到的形变位置与Lei等（2019b）给出的兴文地震定位结果有较大偏差，而地震定位结果附近无明显的地表形变（图9a）。此外，经精定位的筠连地震位置距离我们的研究区C较远，且在形变图中未观测到明显的地表形变，因此也不会对研究区C的形变有影响。通过以上分析，我们认为研究区C内的形变分布特征主要反映了两三年内水力压裂作业向目标地层注水及返排操作引起的地表形变。图9a也很好地揭示了一些较大地震事件与页岩气开采形变之间的空间分布关系，对于这些地震的诱发模型分析非常重要。

图  9  M_S≥5.0地震引起的地表形变与页岩气开采区形变的关系

（a） 2017年6月12日至2019年7月8日期间的ALOS-2降轨页岩气区块形变场，图中矩形A和B分别为珙县地震和兴文地震引起的地表形变，区域C为页岩气开采引起的地表形变，地震定位结果来自Lei等（2019b），H18井位置源于本研究实地野外调查；（b） 2018年12月30日至2019年1月11日期间的Sentinel-1 LOS向InSAR形变图

Figure  9.  The relationship between the surface deformation caused by the M_S≥5.0 earthquakes and the deformation of shale gas exploitation

（a） ALOS-2 descending deformation field of shale gas block from June 12，2017 to July 8，2019，where the rectangles A and B are the surface deformation caused by the Gongxian and Xingwen earthquakes，C is the surface deformation caused by shale gas exploitation，earthquake location results are from Lei et al （2019b），and the location of well H18 is from our field investigation； （b） Sentinel-1 LOS InSAR deformation map from December 30，2018 to January 11，2019

下载: 全尺寸图片 幻灯片

为了验证ALOS-2数据得到结果的正确性并提高形变的观测精度，我们利用2015年至2019年获取的Sentinel-1 SAR数据进行时序PS-InSAR处理，结果显示两个传感器得到的形变特征有很好的一致性。虽然采用了不同的数据和不同的处理方法，地表观测到的形变特征都能反映页岩气开采造成的地面形变及其变化过程。结合四种不同观测角度的数据，我们也可以分析得到地表形变的主要活动方式，初步揭示出页岩气生产过程中造成地面形变的复杂性。Sentinel-1数据提供了更好的观测精度，也提供了页岩气区块形变的一些细节特征。通过形变与地震反射剖面的比较，我们初步确认在开采区南侧边界存在一条较宽的东西向断层破碎带，其北界如果对应于一条隐伏断层，在持续的开采和注液条件下，极有可能成为具有较高诱发地震风险的地震断层，应在未来的页岩气开采和诱发地震研究中予以重点关注。

ALOS-2数据由JAXA （JAXA 2nd Research Announcement on the Earth Observations，ER2A2N171）提供，中国科学技术大学李俊伦教授提供了页岩气开采井位信息，研究中与日本地质调查局雷兴林教授、四川省地震局易桂喜研究员、中国地震局地球物理研究所房立华研究员、中国地震局地质研究所鲁人齐研究员进行了有益讨论，审稿专家给出了建设性的修改意见，作者在此一并表示感谢。