基于InSAR 數(shù)據(jù)分析四川盆地南部長寧頁巖氣區(qū)塊地表形變場基本特征＊

包雨鑫 孫建寶， 李 濤 梁存任詹 艷 韓 靜 李永生 張景發(fā)

1） 中國北京 100029 中國地震局地質(zhì)研究所

2） 中國北京 100029 地震動力學(xué)國家重點實驗室

3） 中國北京 100871 北京大學(xué)遙感與地理信息系統(tǒng)研究所

4） 中國北京 100084 應(yīng)急管理部國家自然災(zāi)害防治研究院

引言

頁巖氣資源是我國重要的戰(zhàn)略能源之一，四川盆地有豐富的頁巖氣資源，然而大規(guī)模水力壓裂開采使得開采區(qū)周圍的地震數(shù)量急劇增加（Leiet al，2017，2019b；Menget al，2019），一些較大地震的發(fā)生給開采區(qū)帶來了新的災(zāi)害問題。因此，在保證能源生產(chǎn)的同時，也要保證開采區(qū)免受人工誘發(fā)地震的影響。近年來，誘發(fā)地震在全球都是一個重要科學(xué)問題（Jacksonet al，2014；Norriset al，2016；Keranen，Weingarten，2018）。四川盆地由于特殊的地質(zhì)構(gòu)造條件，這個問題更加突出，一些較大的地震事件被懷疑是工業(yè)開采所誘發(fā)，甚至可能是全球最大的誘發(fā)事件（Leiet al，2019a；Schultzet al，2020）。通常的誘發(fā)地震主要是依靠地震觀測網(wǎng)絡(luò)監(jiān)測，即通過地震學(xué)的方法進行誘發(fā)地震的識別和定位，找出潛在的活化斷層，評估其地震風(fēng)險（Zenget al，2014；易桂喜等，2019；Liu，Zahradník，2020）。通過干涉合成孔徑雷達（interferometric synthetic aperture radar，縮寫為InSAR）形變手段能否在誘發(fā)地震活動區(qū)監(jiān)測到有意義的地表形變，并進一步解釋誘發(fā)地震的發(fā)生機理和評估誘發(fā)地震的危險性，已經(jīng)成為一個新的研究方向（Bao，Eaton，2016；Barbouret al，2016）。特別是對于頁巖氣開采的監(jiān)測，能否得到開采相關(guān)的形變信息，提供開采中流體活動和斷層活化的證據(jù)，為規(guī)避誘發(fā)地震災(zāi)害風(fēng)險提供新的觀測手段，更值得深入探討。Wang 等（2020）以及Li 等（2021）都利用ALOS-2 和Sentinel-1 數(shù)據(jù)觀測到了長寧頁巖氣區(qū)塊有明顯的地表形變，但其重點給出的是2019 年6 月17 日長寧MS6.0 地震與鹽礦注水的關(guān)系，并未對長寧頁巖氣區(qū)塊的地表形變特征展開分析。

鑒于此，本文基于已有研究，排除長寧地震對研究區(qū)形變信息的干擾，利用永久散射體（persistent scatterer，縮寫為PS）時間序列方法對長寧地震前所有可用的哨兵數(shù)據(jù)進行處理，并對地表形變特征進行系統(tǒng)性分析，以解釋開采活動造成地面形變的基本特征及其機制。

1 四川盆地南部頁巖氣開采與地震活動監(jiān)測

1.1 四川盆地的構(gòu)造活動

四川盆地位于青藏高原、秦嶺與東揚子丘陵之間，是一個菱形的山間盆地。該盆地的形成受控于東北邊界的米倉山—大巴山逆斷裂帶、西北邊界的龍門山逆斷裂帶、西南邊界左旋走滑的鮮水河—小江斷裂系和東南邊界的華鎣山褶皺帶（圖1a，1b）。這些構(gòu)造形成于中生代（Burchfielet al，1995），新生代以來隨著青藏高原向東擴張又被重新活化（Burchfielet al，1995；Xu，Kamp，2000；Wanget al，2014；Yanget al，2017；Tianet al，2018）。強震記錄表明，四川盆地西部和西南邊界的現(xiàn)今變形仍然非常強烈，近期曾發(fā)生1974 年MS7.1 大關(guān)地震、2008 年MW7.9 汶川地震、2013 年MW6.6 蘆山地震等強破壞性地震。相對而言，盆地東部邊界和東北邊界的現(xiàn)今地震活動和構(gòu)造變形較為微弱（Wang，Shen，2020）。四川盆地南部長寧頁巖氣區(qū)塊發(fā)育有建武向斜，該向斜開闊舒緩，其兩翼地層傾角為5°—10°，核部地層近水平，向北與長寧背斜相連（圖1c）。

圖1 四川盆地構(gòu)造背景及地震活動性時空分布（a） 四川盆地地震構(gòu)造和1970 年之后的地震分布，地震目錄數(shù)據(jù)來源于國家地震科學(xué)數(shù)據(jù)中心（2021），GPS 數(shù)據(jù)來源于Wang 和Shen （2020）Fig. 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 四川盆地構(gòu)造背景及地震活動性時空分布（b） 2009 年以來威遠和長寧地區(qū)地震的空間分布及M4 以上地震的震源機制解（GCMT，2021），其中紅色震源機制球代表疑似采鹽注水引起的地震事件，藍色代表疑似水力壓裂引起的事件，黑色代表天然地震；（c） 長寧頁巖氣開采區(qū)塊附近褶皺跡線和地表斷層分布；（d） 2009 年以來長寧區(qū)塊M＞0 地震的時間演化圖Fig. 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 頁巖氣開采與地震活動

頁巖氣開采通常采用水力壓裂技術(shù)，通過向頁巖層中注入數(shù)十個MPa 的水、砂和其它化學(xué)物質(zhì)組成的壓裂液和支撐劑等，使得致密的頁巖層在某些方向產(chǎn)生裂隙，把蘊藏在裂縫中的天然氣釋放出來（Schultzet al，2020）。在脆性巖石被壓裂的過程中，隨著裂縫打開通常伴有微震的產(chǎn)生，震級通常＜0，這種被稱為“壓裂微震”的地震活動并不會造成災(zāi)害，相反會給壓裂操作提供有用的壓裂效果信息（Leiet al，2020）。水力壓裂作業(yè)也可能使得附近不同規(guī)模的先存斷層活化并在一定應(yīng)力條件下發(fā)生黏滑，從而誘發(fā)災(zāi)害性地震，這是誘發(fā)地震研究的重點。在美國中部的俄亥俄州和俄克拉荷馬州以及加拿大西部的頁巖氣壓裂開采作業(yè)區(qū)域，壓裂致震的情況較為普遍，一些地震學(xué)的觀測證據(jù)表明水力壓裂過程確實會引發(fā)有感地震（Atkinsonet al，2016；Bao，Eaton，2016；Schultzet al，2020）。

四川盆地多期構(gòu)造沉降和造山運動使得該地區(qū)富有油、氣、煤炭和其它礦產(chǎn)資源，是中國頁巖氣儲量最大的地區(qū)（蒲泊伶等，2010；Jianget al，2016）。2014 年之前，整個四川盆地的非天然地震活動主要受長期廢水回注（Leiet al，2008，2013）和井鹽開采注水的影響（阮祥等，2008；朱航，何暢，2014；Sunet al，2017）。自2014 年開始，四川盆地內(nèi)先后開發(fā)了威遠、長寧、富順永川和涪陵焦石壩等頁巖氣開采區(qū)塊（董大忠等，2014），為期幾個月的短期大量高壓注液可能是該地區(qū)高水平小地震活動的主要原因。不同構(gòu)造環(huán)境的地震活動性不同，其中長寧區(qū)塊的誘發(fā)地震活動水平較其它區(qū)塊更高。長寧區(qū)塊內(nèi)的上羅頁巖氣區(qū)塊位于建武向斜（圖1c），目標(biāo)頁巖層為志留系龍馬溪組頁巖，埋藏深度為2.3—3.0 km （Leiet al，2020）。我們統(tǒng)計了長寧區(qū)塊內(nèi)2009—2020 年間的地震分布情況，如圖1d 所示，可以看出：在2014 年以前，研究區(qū)的地震活動水平較低；但自從規(guī)?；毫鸦顒右詠?，地震活動頻率顯著上升，其空間分布與壓裂井分布有很高的相關(guān)性（圖1b 和1d）。長寧頁巖氣區(qū)塊內(nèi)2018 年以前未記錄到超過M5 的地震，2018 年12 月16 日四川省宜賓市興文縣發(fā)生ML5.7 地震，震源矩心深度為3 km 左右（Leiet al，2019b），這次地震造成了多人受傷。Lei 等（2019b）認為這次地震可能是由長寧上羅頁巖氣區(qū)塊水力壓裂操作引起，可能是迄今為止全球最大的頁巖氣水力壓裂誘發(fā)的地震活動（Schultzet al，2020）。興文地震源于一個NNW 向的單向破裂，長寧背斜的NW 向伸展構(gòu)造可能起到了阻止該破裂擴展的作用，并有可能促進2019 年6 月17 日長寧MS6.0 地震的發(fā)生（Leiet al，2019a，2020）。一些研究認為長寧MS6.0 地震及其后的震群活動可能由井鹽開采誘發(fā)（Leiet al，2019a；Liu，Zahradník，2020；Wanget al，2020；Liet al，2021），而非頁巖氣開采誘發(fā)的結(jié)果。與頁巖氣短期壓裂注水不同，采鹽時從注水井向地下2.7—3.0 km 深處長期注入淡水，溶解鹽礦并從對應(yīng)的抽水井抽回鹽水。對于長寧地震而言，注水井底部很可能與InSAR 反演推斷出來的發(fā)震斷層面接近相交，而這種對流井技術(shù)可能通過注水提高了斷層帶內(nèi)的孔隙壓力，從而誘發(fā)了長寧地震（Wanget al，2020；Liet al，2021）。本研究將主要討論頁巖氣開采對誘發(fā)地震活動的影響，兩種不同作業(yè)區(qū)域大致以28.3°N 為界，空間上基本沒有重疊，易于分析不同誘發(fā)源及其活動機制。

目前對于水力壓裂引起誘發(fā)地震活動的機制主要有兩種：其一是壓裂井與斷層之間存在流體壓力擴散通道，流體擴散作用增大了斷層上的孔隙壓，從而導(dǎo)致先存斷層發(fā)生滑動（Zhanget al，2013；Bao，Eaton，2016；Schultzet al，2017；Goebel，Brodsky，2018）；其二是壓裂井與斷層之間無流體壓力擴散通道，壓裂井附近的孔隙壓力上升，造成周邊區(qū)域的巖石變形，進而改變周邊區(qū)域內(nèi)與壓裂井不連通的先存斷層上的剪應(yīng)力和正應(yīng)力（Bao，Eaton，2016；Denget al，2016；Goebel，Brodsky，2018）。當(dāng)剪應(yīng)力大于摩擦系數(shù)乘以有效正應(yīng)力時，地震就會隨之發(fā)生，這一機理一般稱為孔隙彈性效應(yīng)。

1.3 誘發(fā)地震InSAR 形變監(jiān)測

目前對四川盆地頁巖氣誘發(fā)地震活動的研究主要使用地震學(xué)方法（Leiet al，2019b；Tanet al，2020）。四川盆地的地質(zhì)情況比北美中西部復(fù)雜得多，地震學(xué)觀測結(jié)果受限于臺站分布和數(shù)據(jù)質(zhì)量，往往不能給出較為確切的反演結(jié)果（張捷等，2021）。與天然地震觀測一樣，監(jiān)測地表形變也是獲取中型或大型誘發(fā)地震震源斷層破裂信息的有效手段。水力壓裂采用多段壓裂，每段向地下注入的壓裂液為大約1 800 m3的滑溜水和約100 t 的支撐劑（任勇等，2015）。然而這些壓裂液不是百分之百返排到地表，返排率通常為15%—30% （韓慧芬等，2017），殘留在地下的壓裂液及其擴散過程可能會引起地層體應(yīng)變變化或者斷層附近孔隙壓的增加，造成地表形變并誘發(fā)地震。因此，本文首先使用常規(guī)InSAR 技術(shù)分析單時相長波ALOS-2 雷達衛(wèi)星數(shù)據(jù)，然后利用永久散射體干涉雷達技術(shù)（PS-InSAR）分析Sentinel-1 多時相雷達衛(wèi)星數(shù)據(jù)，對ALOS-2 的結(jié)果進行驗證并進行更高精度的形變分析，進而結(jié)合當(dāng)?shù)鼐粩?shù)據(jù)分析該地區(qū)地表發(fā)生形變的原因。

日本宇航局（Japan Aerospace Exploration Agency，縮寫為JAXA）的L 波段ALOS-2 數(shù)據(jù)（波長為24.2 cm，重訪周期為14 d），相比歐州航天局（European Space Agency，縮寫為ESA）的Sentinel-1 數(shù)據(jù)（波長為5.6 cm，最短重訪周期為6 d），能夠在較長的時間跨度內(nèi)保持較高的干涉相干性，獲得足夠信噪比的空間連續(xù)地表形變場。由于ALOS-2 衛(wèi)星在本文研究區(qū)獲得的觀測數(shù)據(jù)很少，所以只能通過常規(guī)差分干涉處理，尚無法形成時間序列。四川盆地內(nèi)植被茂盛，濕度大，雷達電磁波的散射特性隨時間的推移變化較快，因此會引起時間去相干現(xiàn)象（Zebker，Villasenor，1992）。而ALOS-2 數(shù)據(jù)能在很大程度上克服該問題，是比較理想的觀測數(shù)據(jù)。常規(guī)InSAR 方法處理Sentinel-1 數(shù)據(jù)由于相干性的降低而很難得到足夠信噪比的信號，因此，我們利用ALOS-2 升軌掃描模式合成孔徑雷達數(shù)據(jù)（ScanSAR，地面覆蓋范圍350 km×350 km）和降軌條帶模式合成孔徑雷達數(shù)據(jù)（Stripmap，地面覆蓋范圍70 km×125 km）對長寧及其周邊地區(qū)展開形變觀測，期望獲得研究區(qū)在觀測時段內(nèi)完整連續(xù)的地表形變場，以便評估頁巖氣開采過程中可能造成形變分布的總體特征。使用Sentinel-1 數(shù)據(jù)雖然無法通過單幅干涉對有效地獲得較長時段內(nèi)連續(xù)的InSAR 形變場，但是通過永久散射體技術(shù)對Sentinel-1 數(shù)據(jù)進行時間序列處理，仍然有望獲得地面一些相位穩(wěn)定目標(biāo)（離散點）的形變信息，可以用來驗證由ALOS-2 數(shù)據(jù)獲得的地表形變信息。此外，關(guān)于不同軌道Sentinel-1 數(shù)據(jù)長時間序列的InSAR 分析，有助于我們尋找可能存在的隱伏斷層，對誘發(fā)地震機制的研究及其災(zāi)害預(yù)防具有重要意義。對于開采部門來說，尋找隱伏斷層也有利于避免生產(chǎn)中套管變形等問題帶來的經(jīng)濟損失（陳朝偉等，2019）。

大地測量手段越來越多地被應(yīng)用于長期的形變監(jiān)測中。Vasco 等（2010）使用基于衛(wèi)星的觀測手段監(jiān)測注入地下大量二氧化碳所引起的地表形變，并將地表形變監(jiān)測和數(shù)值模擬相結(jié)合來約束地質(zhì)力學(xué)參數(shù)。Comola 等（2016）以哈薩克斯坦騰格茲水庫為例，建立了油藏的三維地質(zhì)力學(xué)模型，并與InSAR 大地測量手段結(jié)合起來，推斷該地區(qū)的地質(zhì)力學(xué)特征。Shirzaei等（2019）開發(fā)了基于地表形變和廢水回注數(shù)據(jù)的聯(lián)合反演方法，反演了目標(biāo)地層的水力擴散系數(shù)，證明了InSAR 手段應(yīng)用于約束地層水文地質(zhì)參數(shù)方面的潛力。Deng 等（2020）在美國德克薩斯西部選取了三個地點，利用InSAR 時間序列手段測量地表形變，并研究了該地區(qū)的地震活動與石油和天然氣生產(chǎn)之間的關(guān)系。目前也有研究人員對頁巖氣生產(chǎn)誘發(fā)地面形變場進行了初步分析（Kubaneket al，2018；Jordanet al，2019），但對于四川盆地頁巖氣區(qū)塊地表形變場的特征缺乏系統(tǒng)性分析和物理機制方面的研究。

2 實驗數(shù)據(jù)與處理方法

本研究首先基于日本宇航局的L 波段ALOS-2 PALSAR 數(shù)據(jù)，使用常規(guī)InSAR 技術(shù)監(jiān)測長寧頁巖氣區(qū)塊兩三年間的累積地表形變分布。與C 波段雷達數(shù)據(jù)相比，即使在四川盆地這樣地形起伏較大、濕度大、植被茂盛的地區(qū)，InSAR 數(shù)據(jù)也能保持良好的相干性。其次，我們使用2015—2019 年間獲取的升降軌道Sentinel-1 數(shù)據(jù)進行PS-InSAR 時間序列分析，以克服常規(guī)InSAR 技術(shù)受時空失相干、地形誤差、大氣延遲效應(yīng)等誤差的影響，獲得較常規(guī)InSAR 方法更高精度、更長時間序列的形變信息。對于ALOS-2 數(shù)據(jù)來說，由于長波信號會受到電離層中電子密度變化的巨大影響，對于不同觀測時間和空間，電離層誤差可能達到厘米級甚至是米級（Sandwellet al，2008），這在汶川地震期間ALOS-1 衛(wèi)星的InSAR 觀測中表現(xiàn)得最為明顯（Shenet al，2009）。本研究中使用的ALOS-2 數(shù)據(jù)同樣也受到這種誤差的嚴重影響，因此必須考慮這種誤差的影響才能獲得有效的形變信息。此外，四川盆地常年水汽較重，而且其時空變化非常快，這是影響InSAR 觀測誤差的主要因素，特別是在單時相InSAR 觀測中對結(jié)果精度影響較大，因此我們也嘗試使用大氣模型輔助對該誤差進行校正。

2.1 ALOS-2 數(shù)據(jù)的InSAR 處理和誤差校正過程

常規(guī)InSAR 干涉圖由同一地區(qū)獲取的兩幅配準(zhǔn)后的單視復(fù)數(shù)（single look complex，縮寫為SLC）圖像進行復(fù)共軛相乘得到，反映的是地面上同一位置在兩次成像間隔的相位差，去除地形相位后每一個像素的差分干涉相位由六部分組成：

式中：Δφtopo為所使用數(shù)字高程模型（digital elevation model，縮寫為DEM） 誤差引起的剩余高程相位；φtropo為大氣對流層延遲相位，是在兩次觀測時由于對流層差異引起的相位誤差；φion為兩次成像時電離層電子密度變化引入的相位；φdef為兩次成像期間內(nèi)地表移動引起的雷達與地面之間的視線向（light of sight，縮寫為LOS）變化；Δφorb為剩余軌道誤差相位，對于現(xiàn)代衛(wèi)星基本可以忽略；隨機噪聲φn為失相干噪聲。對于地表形變測量來說，只有φdef是我們希望得到的信息，其它相位均為誤差項，尤其是大氣電離層φion和對流層φtropo的影響，足以把厘米級的形變掩蓋。對于常規(guī)InSAR 處理，我們可以利用誤差相位分量的不同空間尺度特點或者外部輔助數(shù)據(jù)去除它們。對于L 波段ALOS-2 數(shù)據(jù)，四川盆地所處的緯度地區(qū)受電離層電子密度變化影響較大，需要特別考慮電離層誤差的影響。

利用ALOS-2 不同成像觀測模式，并選取升降軌兩種不同成像幾何的合成孔徑雷達（synthetic aperture radar，縮寫為 SAR）數(shù)據(jù)進行干涉處理并消除地形相位的影響，得到常規(guī)InSAR 的差分干涉結(jié)果，如圖2 所示。圖2a 顯示整個圖像范圍均分布著復(fù)雜的不規(guī)則干涉條紋，由于ALOS-2 衛(wèi)星較ALOS-1 衛(wèi)星在軌道控制方面有顯著改善，圖像中的軌道誤差所產(chǎn)生的影響很小，所以這種大尺度不規(guī)則條紋與形變無關(guān)，也無法用線性或者非線性軌道模型擬合，主要為電離層相位延遲誤差（Liang，F(xiàn)ielding，2017）。圖2b 顯示了布滿整個圖像的規(guī)則干涉條紋，后文基于電離層物理模型的相位分析表明，它們也主要是電離層相位誤差。另外，圖2a 和2b 中大尺度的干涉條紋也可能疊加了長波對流層延遲誤差。由于Δφtopo與垂直基線有關(guān)，本文采用1″的航天飛機雷達地形測繪使命（Shuttle Radar Topography Mission，縮寫為SRTM）的DEM 數(shù)據(jù)在小基線情況下（表1）引入的誤差相對較小，失相干噪聲φn對最后結(jié)果的影響較小，在數(shù)據(jù)相干性較好時可以忽略不計（Hanssen，2001）。對形變結(jié)果產(chǎn)生較大影響且較難校正的是大氣延遲引起的相位誤差，包括對流層延遲和上述電離層延遲兩種。對于四川盆地來說，成像期間水汽分布劇烈變化會造成較大的對流層延遲誤差，對所有波段的雷達數(shù)據(jù)均有較大影響。當(dāng)使用長波雷達數(shù)據(jù)時，電離層效應(yīng)在不同緯度、不同時間對合成孔徑雷達信號的干擾程度不同，而在四川地區(qū)受電離層影響較大（Shenet al，2009，Gombaet al，2016）。對于C 波段Sentinel-1 數(shù)據(jù)，電離層在低緯度和高緯度如極地地區(qū)對合成孔徑雷達信號有嚴重影響，但本研究所考慮的緯度地區(qū)相對小，因此下面重點介紹大氣相位延遲的建模和去除方法。

圖2 功率譜濾波后的ALOS-2 差分干涉相位結(jié)果（a） 升軌ScanSAR 模式差分干涉圖，成像時間為2016 年8 月7 日至2019 年9 月15 日；（b） 降軌StripMap 模式差分干涉圖，成像時間為2017 年6 月12 日至2019 年7 月8 日Fig. 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 數(shù)據(jù)的參數(shù)Table 1 ALOS-2 PALSAR data parameters used in this study

2.1.1 去除電離層相位

電離層是大氣中的電離部分，衛(wèi)星雷達傳感器發(fā)射的電磁波信號兩次穿過電離層，會對SAR 信號的相位產(chǎn)生顯著影響（Gombaet al，2016；王楠等，2017）。對于ALOS-2 衛(wèi)星數(shù)據(jù)來說，其L 波段信號的頻率較Sentinel-1 衛(wèi)星的C 波段信號低，電離層對低頻雷達系統(tǒng)干涉相位的影響不能忽略。因此，在數(shù)據(jù)處理中需要估計電離層相位的影響，并將其從干涉圖中去除。如果考慮電離層的影響，InSAR 干涉相位的組成可表示為

式中：f0為載波頻率；ΔTEC 為LOS 方向上兩次成像的總電子含量（total electron content，縮寫為TEC）差值；K為常數(shù)，取40.28 m3/s2；c為光速（Lianget al，2018）。式（3）表示干涉相位 ［ 式（2） ］ 中的非頻散相位，與頻率成正比，主要是剩余地形、形變、對流層和剩余軌道相位，式（4）表示由電離層變化引起的頻散相位，與頻率成反比。

基于電離層相位的頻散特性，可以將SAR 距離向信號頻譜分割為不同頻段來估計電離層相位，即距離向信號頻譜分割法（Brcicet al，2010；Rosenet al，2010；Gombaet al，2016；Lianget al，2018）。本文對ALOS-2 數(shù)據(jù)應(yīng)用距離向頻譜分割法得到差分電離層相位（Liang，F(xiàn)ielding，2017），該方法將需要作干涉的兩幅SAR 圖像分割為兩個不同距離向中心頻率的子帶圖像：上子帶圖像（頻率為fu）和下子帶圖像（頻率為fl）。對兩幅上子帶圖像干涉得到干涉相位Δφu，對兩幅下子帶圖像干涉得到干涉相位Δφl，在不考慮相位展開誤差和其它相位誤差的情況下，將兩子帶圖像的干涉相位代入

便可以計算得到電離層相位φion（Lianget al，2018；Zhanget al，2018），如圖3 所示。對比原始差分干涉圖（圖2），可知：形變場中的長波相位與電離層延遲相位非常相似，由此可以確定圖2 中的誤差主要為電離層相位，去除該相位影響的結(jié)果如圖3c 和圖3d 所示。去除電離層相位后的差分干涉圖還主要有對流層延遲誤差相位，下面繼續(xù)消除這部分誤差。

圖3 利用距離向信號頻譜分割法估計電離層延遲相位（a） 升軌ScanSAR 模式圖像的電離層相位估計結(jié)果；（b） 降軌StripMap 模式圖像的電離層相位估計結(jié)果；（c） 升軌ScanSAR 模式圖像去除電離層相位的纏繞相位；（d） 降軌StripMap 模式圖像去除電離層相位的纏繞相位Fig. 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% 的相對濕度變化即可導(dǎo)致干涉圖上的信號產(chǎn)生高達10—14 cm 的形變測量誤差（Zebkeret al，1997；孫廣通等，2011），而這個誤差可能會掩蓋干涉圖上的地表微小形變，對于地殼形變場的監(jiān)測非常不利。根據(jù)物理機制的不同可將對流層延遲分為兩部分：天頂干延遲（又稱垂直分層延遲）和天頂濕延遲（又稱水平湍流混合延遲）（Parkeret al，2015）。垂直分層延遲相位與高程呈反比關(guān)系，即隨著高度的增加，水汽含量降低，大氣延遲相位減??；水平湍流混合延遲主要是由水汽的含量變化和隨機分布引起的水平紊亂延遲（Lohman，Simons，2005；Ebmeieret al，2013）。在平坦地區(qū)，水平湍流延遲是主要因素，此時垂直分層延遲受壓力主導(dǎo)，在空間上表現(xiàn)得較為平滑；在地形起伏明顯的多山區(qū)域，由海拔變化引起的垂直分層延遲是主要因素（Doinet al，2009）。

式中：φtropo0代表從高度h1到對流層頂部htop的對流層延遲相位，與大氣折射率N有關(guān)；Nhydr代表靜力學(xué)平衡折射率；Nwet代表非靜力學(xué)折射率；θ為影像獲取時的入射角；φtropo為獲取主輔影像時大氣折射率的時空變化所導(dǎo)致的對流層延遲相位，其中上標(biāo)s 表示從圖像，上標(biāo)m 表示主圖像（Bekaertet al，2015）。

目前針對對流層改正的方法有三大類：① 基于外部輔助數(shù)據(jù)法，外部數(shù)據(jù)源有GPS 數(shù)據(jù)、大氣模型（ERA5，ERA-I，WRF 等）、多光譜觀測數(shù)據(jù)（MODIS，MERIS 等）；② 基于經(jīng)驗關(guān)系的干涉相位法，包括線性相關(guān)模型和冪律模型；③ 基于數(shù)據(jù)本身的校正技術(shù)，主要有相位疊加平均法（stacking）和時序InSAR 方法等（周洪月，2018）。不同的改正方法對于干、濕延遲分量的敏感度不同，都有各自的優(yōu)缺點。MERIS 和MODIS 光譜儀數(shù)據(jù)只能提供有光照條件下的數(shù)據(jù)，對云的敏感度很大，并且只能估計濕延遲，而線性相關(guān)模型和冪律模型只能估計受地形起伏影響大的干延遲。長寧地區(qū)屬于多山區(qū)域，且水汽含量較大，需要去除對流層的影響。下文將分析如何從干涉相位中去除對流層延遲影響。

由歐洲中尺度天氣預(yù)報中心（European Centre for Medium-Range Weather Forecasts，縮寫為ECMWF）提供的最新一代全球氣象再分析產(chǎn)品（ECMWF Re-analysis v5，縮寫為 ERA5）既可以估計干延遲分量，又可以估計濕延遲分量，相比于之前的ERA-I 模型來說有更好的時空分辨率（Hu，Mallorquí，2019）。本文中對流層延遲改正使用外部輔助數(shù)據(jù)（ERA5 模型），利用TRAIN 軟件實現(xiàn)。對流層校正后的圖像如圖4a 所示，可見：在本例中ERA5 模型對整個圖像大尺度對流層延遲相位的改善明顯，但是對小尺度對流層延遲相位作用有限，對我們的研究區(qū)僅有幾個弧度的改正。

在四川盆地這種水汽變化較快的地區(qū)，ERA5 模型可能無法有效地消除長波對流層誤差，同時用大氣模型進行對流層改正時也可能會引入額外的長波誤差，我們將這部分誤差稱為長波趨勢相位誤差，如圖4b 的情況。前文提及的軌道誤差在新一代衛(wèi)星中影響很小，可忽略不計。因此，我們用一階線性模型對上述校正后的圖像進行趨勢相位誤差校正，去除剩余的長波誤差，得到最終校正后的差分干涉圖（圖4c，4d）。雖然圖4c 中仍然能看到一些比較明顯的誤差干涉條紋，其主要來源于大范圍（350 km×350 km）對流層延遲誤差，但是在我們關(guān)注的頁巖氣開采區(qū)塊內(nèi)（圖中矩形框區(qū)域），這種誤差的影響非常有限，不影響對結(jié)果的分析。圖4d 中，通過上述誤差校正后已經(jīng)看不到明顯的誤差條紋，對比圖3d，圖像中無論是長寧地震的形變，還是頁巖氣區(qū)塊的形變，其信噪比均較高。除了長寧地震，一些較大的地震事件（如MW＞4.0）的同震形變場也會對InSAR 形變場產(chǎn)生影響，但其位置已知，且形變范圍有限，比較容易識別，不會影響頁巖氣區(qū)域的形變分析。因此，我們將圖4c 和圖4d 的結(jié)果作為最終的ALOS-2 形變觀測結(jié)果。

通過上述過程的多步校正，特別是針對ALOS-2 數(shù)據(jù)進行的電離層延遲誤差校正，其結(jié)果顯示利用長波SAR 數(shù)據(jù)進行四川盆地或者龍門山地區(qū)的InSAR 形變信息提取有重要的應(yīng)用價值。同時，L 波段的ALOS-1 或者ALOS-2 數(shù)據(jù)能在長時間內(nèi)保持較高的相干性，對于復(fù)雜環(huán)境下，如復(fù)雜地面植被和地形條件等，獲得高信噪比的InSAR 觀測結(jié)果非常有利。但是L 波段的觀測數(shù)據(jù)較少，只有日本宇航局的ALOS-2 可以使用，該數(shù)據(jù)沒有針對四川地區(qū)的系統(tǒng)性數(shù)據(jù)獲取計劃，因此還無法形成比較好的觀測時間序列，阻礙了長波InSAR 觀測的應(yīng)用。為了獲得比較連續(xù)的形變時間序列，只能利用歐洲航天局的Sentinel-1 數(shù)據(jù)，具體將在下節(jié)中進一步分析。

2.2 Sentinel-1 數(shù)據(jù)的PS-InSAR 處理

在ALOS-2 衛(wèi)星觀測到的常規(guī)干涉形變場中，我們能夠清晰地看到長寧頁巖氣區(qū)塊存在顯著的地表形變（圖4），且不同時段的形變存在一定差異，而該區(qū)域頁巖氣開采是唯一可能的形變來源。為了驗證ALOS-2 數(shù)據(jù)的觀測結(jié)果，并獲得更高精度的形變信息，我們使用歐洲航天局的Sentinel-1 雷達衛(wèi)星數(shù)據(jù)，通過時間序列分析技術(shù)得到地表形變場。為了避免ALOS-2 升降軌數(shù)據(jù)獲取時間的差異，我們選擇相同時段的Sentinel-1 升降軌數(shù)據(jù)進行分析。本文選用的ALOS-2 數(shù)據(jù)時間間隔較長，因為使用傳統(tǒng)單時相InSAR 技術(shù)，不但會受到失相干的影響，而且還會受到前述電離層和對流層的干擾。特別是四川盆地濃厚大氣所導(dǎo)致的對流層延遲誤差的影響，雖然可以通過ERA5 天氣模型等輔助數(shù)據(jù)進行校正，但是兩者的時空分辨率差異使得校正效果有限（圖4）。為了得到更加準(zhǔn)確的結(jié)果，我們對Sentinel-1 數(shù)據(jù)進行PS-InSAR 處理，從經(jīng)過配準(zhǔn)的超過四年累積的Sentinel-1 單視復(fù)數(shù)圖像（2015 年2 月至2019 年5 月）上提取相位穩(wěn)定的永久散射體目標(biāo)進行分析。該時段頁巖氣開采處于快速發(fā)展期，并且形變場中不包括長寧地震的相關(guān)形變，能夠反映頁巖氣開采產(chǎn)生的形變特征。

圖4 大氣對流層延遲校正及長波趨勢相位校正（a，b） 基于ERA5 模型去除升、降軌對流層延遲誤差后的圖像；（c，d） 長波趨勢相位校正后的圖像，其中藍色方框表示本研究的形變區(qū)Fig. 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

區(qū)別于最早的基于幅值和先驗形變模型識別穩(wěn)定點目標(biāo)的PS-InSAR 處理方法（Ferrettiet al，2000，2001），我們基于斯坦福大學(xué)開發(fā)的StaMPS 軟件（Hooperet al，2004）及其相應(yīng)的相位分析算法進行點目標(biāo)的識別和形變分析。該軟件基于干涉相位的空間相關(guān)性來識別PS 點，不需要假設(shè)的先驗形變模型，適合對非線性的時間變化過程進行監(jiān)測。該方法中永久散射體的識別不取決于雷達圖像幅值，在雷達幅值較低的區(qū)域也能發(fā)現(xiàn)一定量的穩(wěn)定相位目標(biāo)，因此也適用于非城市地區(qū)的形變監(jiān)測。本研究的頁巖氣開采區(qū)塊位于四川盆地內(nèi)部的鄉(xiāng)村地區(qū)，地面上找不到足夠的高幅值散射體，但是用StaMPS 方法依然能夠識別出足夠的PS 點進行相位分析。同時，與穩(wěn)定的構(gòu)造運動不同，頁巖氣開采作業(yè)的時空變化較大，其造成的地面形變也非穩(wěn)態(tài)變化，無法用簡單的形變時間演化模型來表示，但使用基于相位相關(guān)性的算法仍然可以得到合適的形變分析結(jié)果。

StaMPS 方法的主要優(yōu)勢在于PS 點的識別，通過定義時間相干性參數(shù)γx，在一定空間相關(guān)性尺度范圍內(nèi)迭代篩選PS 點。

式（8）為時間相干性的表達式，其中不含幅值相關(guān)信息，所以該方法識別出的PS 點目標(biāo)理論上與幅值無關(guān)。另外，式（8）中也未涉及形變的時間模型，因此PS 點目標(biāo)的選取與形變量值高低無關(guān)，如果其相位是穩(wěn)定的，這個點依然可以作為PS 點，這也是該算法可以用于非城市地區(qū)觀測的原因，這有利于時間非穩(wěn)態(tài)信息的提取。

篩選出PS 點后，從纏繞的干涉相位減去式（8）中估計得到的空間非相關(guān)誤差，使得相鄰PS 點之間的相位差小于Π，此時可以對剩余干涉相位進行展開。纏繞的干涉相位為

本文搜集了跨越研究區(qū)升降兩條軌道（A055 和D164）的Sentinel-1 數(shù)據(jù)進行PS-InSAR 分析，時間跨度為2015 年2 月至2019 年5 月底，截止于長寧地震前兩周左右。從圖像中截選出長寧頁巖氣開采區(qū)及其周圍的數(shù)據(jù)進行單一主景干涉處理（地面覆蓋見圖5），升軌形成84 幅干涉圖，降軌形成89 幅干涉圖，利用StaMPS 軟件分別對升降軌道作PS-InSAR 時序處理。升軌以2017 年11 月10 日獲取的數(shù)據(jù)為主景，最大垂直基線為122 m，最大時間基線為1 116 d；降軌以2017 年2 月26 日獲取的數(shù)據(jù)為主景數(shù)據(jù)，最大垂直基線為205 m，最大時間基線為972 d。該研究區(qū)升降軌道分別識別出20 萬4 088 個和44 萬4 102 個PS 點，經(jīng)過誤差去除并進行PS 點插值后，得到觀測時段內(nèi)的平均線性速度場（圖5），該結(jié)果揭示了長寧上羅頁巖氣開采區(qū)在觀測時段內(nèi)的一些典型變形特征。

圖5 基于長寧頁巖氣開發(fā)區(qū)Sentinel-1 雷達數(shù)據(jù)得到的升軌（a）和降軌（b） PS-InSAR 線性速度場圖中黑點表示該區(qū)塊中已知頁巖氣井的空間分布Fig. 5 PS-InSAR linear velocity field in ascending （a） and descending （b） orbits obtained from Sentinel-1 radar data around Changning shale gas blockThe black dots indicate spatial distribution of known shale gas wells in Changning shale gas block

3 觀測結(jié)果分析

我們使用兩種不同的數(shù)據(jù)處理方法，分析了兩種不同的衛(wèi)星數(shù)據(jù)，最終得到了關(guān)于頁巖氣開采區(qū)塊的一些基本變形特征。這些形變結(jié)果雖然受到四川盆地復(fù)雜地面觀測條件、對流層水汽分布變化等因素的影響，但是通過一些有效的技術(shù)手段我們?nèi)匀荒軌虻玫阶銐蛐旁氡鹊挠^測結(jié)果，這一節(jié)將結(jié)合兩種衛(wèi)星的觀測結(jié)果分析這種變形的基本特征及其形成機制。

3.1 ALOS-2 數(shù)據(jù)的InSAR 結(jié)果分析

將圖4c 和4d 中校正完成的干涉圖像放大到本研究關(guān)注的長寧頁巖氣開采區(qū)塊，可以看到長寧頁巖氣開采區(qū)附近不同軌道ALOS-2 數(shù)據(jù)上的地面形變場及其南北向剖面圖，如圖6 所示。圖中28°20′N 附近的形變場（虛線框內(nèi)）為2019 年6 月17 日長寧MW5.8 地震的同震形變場（Liet al，2021）；AA′和BB′剖面附近為長寧—上羅頁巖氣開采區(qū)塊，整個區(qū)塊分布有密集的頁巖氣開采井（圖中黑點所示）。在這些開采井附近，ALOS-2 形變場顯示了顯著的地面形變（LOS 方向）。由于ALOS-2 數(shù)據(jù)成像時間的限制，升降軌道InSAR 觀測時間不一致，圖6a 升軌的觀測時間較圖6b 降軌的觀測時間長將近一年（表1），因此兩者顯示的形變場有一定差異，觀測時間較長的升軌數(shù)據(jù)（圖6a）顯示了更顯著的地面抬升信號（藍色正值區(qū)，假設(shè)地面發(fā)生垂直形變），這與頁巖氣水力壓裂開采中大量向頁巖層注液有關(guān)。在ALOS-2的觀測期間內(nèi)，研究區(qū)范圍內(nèi)發(fā)生過多次MW＞4.0 地震，但是距離本文研究的形變中心區(qū)有一定距離，對頁巖氣形變影響不大，且在觀測時間較長的條件下，對于速度場來說，地震形變的影響被削弱，因此本文分析不考慮這些地震帶來的影響。另外，開采區(qū)的這些井并非同時工作，有些井在水力壓裂階段，而有些井是在采氣生產(chǎn)階段，因此造成了地面形變場的復(fù)雜性。在兩個觀測時段內(nèi)的一個有趣現(xiàn)象是，跨越形變區(qū)相同位置的剖面顯示，地面隆升的幅值和峰值位置均不同：對于升軌數(shù)據(jù)，LOS 向形變峰值達到了約120 mm，出現(xiàn)在10 km 的位置（圖6c）；對于降軌數(shù)據(jù)，LOS 向形變峰值僅約20 mm，出現(xiàn)在13 km 的位置（圖6c）。假設(shè)地面僅為純垂直形變，這說明在2016 年8 月至2017 年6 月這個時段內(nèi)，壓裂生產(chǎn)大規(guī)模開展，長寧頁巖氣區(qū)塊壓裂活動比較劇烈，造成了地面快速大幅度隆升；2017 年6 月至2019 年7 月這段時間內(nèi)，水力壓裂的重心已經(jīng)往南偏移，而且開采的幅度顯著下降，對應(yīng)的地面隆升峰值也大幅降低，僅約20 mm；2017 年下半年開始，由于大多數(shù)氣井壓裂工作已完成并轉(zhuǎn)入生產(chǎn)階段，因此在這段時間內(nèi)原本上升的區(qū)域可能逐漸表現(xiàn)為地面沉降，或者LOS 方向上遠離降軌衛(wèi)星的地面水平運動，亦或兩種運動兼而有之（圖6c 中BB′剖面上0—10 km 位置的負值區(qū)域）。這種形變的變化過程與石油公司公開報道的生產(chǎn)進度以及我們野外實地調(diào)查都一致，表明形變剖面很好地反映了這種生產(chǎn)過程的轉(zhuǎn)變。

圖6 長寧頁巖氣區(qū)塊近場ALOS-2 干涉形變位移場及剖面圖長寧地震序列和水力壓裂致震的震源機制解來源于Lei 等（2019a），其余2015 年和2017 年的兩次地震的震源機制解來源于GCMT （2021）；白色虛線框表示2019 年6 月17 日長寧地震形變場（a） 升軌ScanSAR 模式圖像得到的頁巖氣區(qū)塊形變位移場，成像間隔為3 年；（b） 降軌StripMap 模式圖像得到的頁巖氣區(qū)塊形變位移場，成像間隔為2 年；（c） 圖（a）中的AA′和圖（b）中的BB′形變位移場剖面，兩剖面在同一位置Fig. 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

通過對覆蓋長寧地區(qū)的升降軌ALOS-2 數(shù)據(jù)進行干涉處理，消除電離層誤差、大氣誤差等因素后，得到研究區(qū)可靠的差分干涉形變場。結(jié)果顯示頁巖氣生產(chǎn)過程中能產(chǎn)生足夠大的地面形變，而且也能夠通過InSAR 技術(shù)有效監(jiān)測到。雖然L 波段ALOS-2 數(shù)據(jù)的時間基線較長，但其波長較長，對地面目標(biāo)的散射機制變化不及C 波段敏感，因此即便在四川盆地復(fù)雜的地面條件和濃厚多變的水汽觀測條件下，兩三年內(nèi)仍能保持較高的干涉相干性，并能從干涉圖中清晰分辨研究區(qū)的形變及其演化情況。此外，基于ALOS-2 的單時相InSAR 形變分析也表明，頁巖氣開采會造成地面的快速隆升或者下降，形變的變化呈非穩(wěn)態(tài)，在壓裂階段由于注入流體而使得地面隆升，在生產(chǎn)階段可能由于抽取返排液或者天然氣而使得地面沉降。受限于ALOS-2 數(shù)據(jù)的可用性，目前尚無法做到長時間序列的InSAR 形變觀測，這是制約ALOS-2 應(yīng)用于頁巖氣開采監(jiān)測的主要問題。

3.2 Sentinel-1 數(shù)據(jù)的PS-InSAR 結(jié)果分析

Sentinel-1 數(shù)據(jù)的PS-InSAR 分析結(jié)果（圖5）為整個觀測時段（2015 年2 月至2019 年5 月）的線性速度場，其中建武向斜內(nèi)的形變特征與Wang 等（2020）利用2016 年3 月至2019 年6 月的Sentinel-1 數(shù)據(jù)進行時序分析的結(jié)果相似。本研究僅使用速度場進行分析，因為基于C 波段的InSAR 觀測會受到四川盆地濃厚水汽的嚴重干擾以及季節(jié)性形變的影響（Liet al，2021），形變時間序列中位移場有很大的不確定性，而速度場是基于所有干涉數(shù)據(jù)作最小二乘反演所得，能夠最大限度地降低時間隨機信號，特別是大氣對流層信號的干擾。從圖7d 的形變剖面點的發(fā)散程度可以大致估計形變速度場的觀測精度，本研究中Sentinel-1 的InSAR時間序列能達到的形變觀測精度約為1 mm/a，因此它是頁巖氣開采監(jiān)測的有力工具，可以為未來開采區(qū)地層的水文地質(zhì)參數(shù)反演和人工誘發(fā)地震研究提供重要信息。在長寧頁巖氣開采區(qū)附近（圖5 中黑點分布范圍），升降軌數(shù)據(jù)均顯示了較高的信噪比，能客觀地反映頁巖氣開采造成的地面形變。我們同樣將圖5 的速度場結(jié)果放大到開采形變區(qū)來進一步分析（圖7）。

圖7a 和7b 的范圍與圖6a 和6b 相同，圖中AA′，BB′，CC′和DD′這幾個剖面的位置也相同，這樣處理便于我們比較四次不同的觀測結(jié)果。圖7 中所使用的SAR 數(shù)據(jù)都在長寧地震之前獲取，所以該結(jié)果中未包含圖6 中長寧地震的形變場。另外，選取Sentinel-1 升降軌數(shù)據(jù)的時候盡量保證兩者時間接近，以便比較不同觀測角度的InSAR 形變場。

圖7 長寧頁巖氣區(qū)塊近場升降軌Sentinel-1 數(shù)據(jù)的LOS 向線性速度場圖（a）和（b）中KK′為圖8a 中二維地震反射剖面位置，OO′為圖8b 中斷層破碎帶的北界；圖（c）中紅色直線為線性速度模型擬合的各點位移。（a） 升軌速度場；（b） 降軌速度場；（c） 圖（a）和圖（b）中P1，P2，P3和P4 點的位移時間序列；（d） 升軌數(shù)據(jù)剖面CC′和降軌數(shù)據(jù)剖面DD′上的形變速率變化Fig. 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 P1，P2，P3 and P4 in Figs. （a） and （b），respectively；(d) Variation of deformation rates on the profile CC′ in ascending orbit and profile DD′ in descending orbit

剖面CC′和DD′ （圖7a，7b）的形變速度場的一個主要特征是在升降兩條軌道的12—13 km 處都存在一個明顯的峰值，約為6—8 mm/a （圖7d）。這個峰值的位置恰好與圖6c 降軌ALOS-2 數(shù)據(jù)的峰值位置相同，這反過來說明降軌ALOS-2 觀測到的峰值是可信的。這個峰值的物理意義在于，在Sentinel-1 觀測覆蓋的長寧頁巖氣開采時段內(nèi)，地面形變場以局部隆升為主，隆升主要發(fā)生在南北向剖面上10—15 km 的位置，圖像上顯示為一條近東西向的狹長條帶（圖7a，7b）。在升降軌道不同的觀測成像幾何條件下，該條帶的形變分布和量值都非常接近，因此我們可以確定該形變主要以垂直抬升為主，而兩者的微小差別反映了SAR 衛(wèi)星觀測入射角的差異（圖7d）。形變剖面還反映了一個比較明顯的特征是，在0—8 km 的范圍內(nèi)存在一個明顯的LOS 向負值區(qū)域，該負值區(qū)以5 km 的位置為中心并達到峰值約-2 mm/a。如果將圖7b 中的負值區(qū)與圖6b 中ALOS-2 降軌數(shù)據(jù)的負值區(qū)進行對比，可以看到兩個負值區(qū)域在0—8 km 范圍內(nèi)幾乎是重疊的，因此可以確定該負值區(qū)是真實的地面形變信息，而不是由于Sentinel-1 降軌SAR 數(shù)據(jù)中的噪聲所致。但是對比升軌Sentinel-1 數(shù)據(jù)可知，相同區(qū)域在升軌上并未出現(xiàn)明顯的負值，而是大部分接近于零值。在圖6c 所示的ALOS-2 升軌數(shù)據(jù)AA′剖面上0—5 km 的范圍內(nèi)也未見明顯的負值分布，這與Sentinel-1 數(shù)據(jù)的觀測也是一致的。綜合這四種InSAR 觀測結(jié)果推知，圖7d 中DD′剖面負值出現(xiàn)的原因是存在一個地面水平向的運動分量，其運動方向與升軌衛(wèi)星飛行方向近乎平行（N12°W 向）。這樣的運動使我們在升軌圖像上很難觀測到這個方向的形變，但在降軌條件下觀測到了負值形變。出現(xiàn)這樣的水平向地面形變的一個可能原因是頁巖氣開采注入的壓裂液產(chǎn)生了NW 向的擴散，這種流體的擴散作用也使得地面產(chǎn)生了相同方向的形變。

圖7c 給出了對應(yīng)于圖7a 升軌數(shù)據(jù)的P1點、P2點和對應(yīng)于圖7b 降軌數(shù)據(jù)的P3點、P4點的位移時間序列，可以看出：P1點和P3點分別為升軌和降軌數(shù)據(jù)中抬升量較大的位置，在2015—2019 年這四年期間，累計抬升量達到約30 mm；P2點和P4點分別為升軌和降軌數(shù)據(jù)中發(fā)生沉降較大的位置，在四年期間累計沉降量達到約15 mm。

圖7d 的形變剖面上還有一個值得討論的特征是，在剖面15—16 km 的位置，形變速率值突然由正變負，該形變轉(zhuǎn)向的位置為圖7a 和7b 的OO′虛線所示，在圖像上形變的變化也很明顯，而且在升降軌上均呈現(xiàn)相同的特征。此外，從圖7a 和7b 上可以看出，在OO′剖面以南區(qū)域沒有頁巖氣開采井的分布。因此，我們推測地面形變的變化與該區(qū)域地層結(jié)構(gòu)發(fā)生變化有關(guān)，需要通過長寧頁巖氣開采區(qū)的地層結(jié)構(gòu)及二維地震反射剖面資料予以進一步分析。

圖7 中KK′剖面的地層結(jié)構(gòu)如圖8a 所示，卷入建武向斜變形的沉積地層由下至上依次為震旦系、寒武系、奧陶系、志留系、二疊系、三疊系和侏羅系，缺失泥盆系和石炭系，總厚度可達8000 m，其中震旦系—中三疊系為一套由灰?guī)r、泥巖、泥頁巖、白云巖和石膏巖組成的海相沉積，而之上的上三疊系—侏羅系是以砂巖、礫巖為主的陸源碎屑沉積（圖8a，8c）。頁巖氣開采層位為志留系的龍馬溪組，為一套富含筆石化石的黑色頁巖層，也是頁巖氣開采的主產(chǎn)氣層（何登發(fā)等，2019）。圖8b 給出了長寧頁巖氣開采區(qū)的二維地震反射剖面，在該剖面上可以看到：建武向斜南翼發(fā)育有一條寬約2 km 的斷層破碎帶，其在二維石油物探剖面上表現(xiàn)為地震同向軸的雜亂反射；斷層帶大致傾向南，切穿了震旦系—三疊系的全部沉積蓋層。該斷層破碎帶的北界大概與圖7a 中的OO′位置對應(yīng)，它的存在可能導(dǎo)致OO′剖面以南區(qū)域頁巖氣存儲層的缺失。從形變場來看，該破碎帶以北有密集的頁巖氣開采井分布，地面隆升明顯，而其北界（OO′剖面）可能成為流體擴散的障礙體（圖8b），這使得開采產(chǎn)生的大量流體擴散不到OO′剖面以南區(qū)域，因此該區(qū)域不會發(fā)生隆升，其變形行為與彈性體相似，在北側(cè)快速隆升的過程中彈性應(yīng)力場造成南側(cè)地塊的相對下降，其下降幅值明顯小于北側(cè)的隆升區(qū)。此外，我們在圖6c 中AA′和BB′剖面上發(fā)現(xiàn)，類似的特征對于ALOS-2 數(shù)據(jù)也是存在的，但是由于ALOS-2 單時相InSAR 觀測的噪聲干擾以及干涉濾波操作的影響，這個特征并不明顯。

圖8 跨越長寧開采區(qū)的二維地震反射剖面（a） KK′剖面（圖7）對應(yīng)的地層結(jié)構(gòu)和近年來長寧附近發(fā)生的較大地震事件及其震源機制解（Li et al，2021）；（b） 圖7d 形變剖面中15—17 km 處斷層破碎帶附近的二維地震反射剖面，其位置見圖1c；（c） 長寧地區(qū)的沉積地層序列（Li et al，2021）Fig. 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 數(shù)據(jù)的PS-InSAR 時間序列觀測結(jié)果表明，多時相InSAR 觀測技術(shù)對于頁巖氣形變監(jiān)測是可行的，它能夠提供比ALOS-2 單時相InSAR 觀測更高的精度，也能更好地反映頁巖氣開采的一些形變細節(jié)特征。

4 討論與結(jié)論

本文基于兩種SAR 衛(wèi)星數(shù)據(jù)，使用兩種InSAR 技術(shù)得到了長寧頁巖氣開采區(qū)塊的地表形變場。針對L 波段ALOS-2 衛(wèi)星數(shù)據(jù)提供了電離層校正、對流層校正等一系列處理方案，得到了在嚴重電離層和對流層干擾下長達兩三年的InSAR 干涉結(jié)果。該衛(wèi)星數(shù)據(jù)的分析結(jié)果表明，頁巖氣開采能夠產(chǎn)生非常顯著的地面形變場，而且隨著開采的進展，地表形變也隨之變化，呈現(xiàn)出非穩(wěn)態(tài)的形變特征。由于頁巖氣開采需要向頁巖層注入大量的流體，而且這些流體大部分會留在地層中而無法返排到地面，因此造成地層體應(yīng)變變化和地面在開采過程中局部強烈隆升，本研究中隆升量達到了120 mm。在頁巖氣壓裂完成轉(zhuǎn)入生產(chǎn)階段，生產(chǎn)過程中駐留在地層中的流體逐步向外擴散，會使得地面下沉并產(chǎn)生特定方向的水平運動。

長寧頁巖氣區(qū)塊內(nèi)近年來發(fā)生了一些較大地震事件，特別是2019 年1 月3 日珙縣MS5.3 地震、2018 年12 月16 日MS5.7 興文地震以及更早的2017 年1 月27 日MS5.0 筠連地震等（圖8），這些地震均發(fā)生在建武向斜內(nèi)，并與長寧頁巖氣開采區(qū)塊非常接近。從形變場的角度如何認識這些地震與頁巖氣開采的關(guān)系，對后續(xù)誘發(fā)地震模型的分析具有重要意義。為了避免InSAR 時間序列中同震形變場信號被時間序列反演過程削弱的影響，我們重新處理了ALOS-2 降軌數(shù)據(jù)，并盡可能降低InSAR 相位功率譜濾波器的強度，以保證這些小地震的形變信息能夠被清晰保留（圖9a）。除了筠連地震由于震級較小形變難于識別外，珙縣地震和興文地震的形變信號比較清楚。我們也用Sentinel-1 數(shù)據(jù)處理獲得了珙縣地震的同震形變場，得到了間隔12 天的形變場圖像（圖9b），但興文地震的InSAR 結(jié)果在Sentinel-1 數(shù)據(jù)上失相干較重，僅能對形變的范圍進行精確定位。從圖9a 中可以清楚地識別出三個主要的形變區(qū)A，B和C，同時也標(biāo)記出了可能與珙縣地震相關(guān)的N201 區(qū)塊的H18 井的位置，這口井自2018 年11 月開始直到珙縣地震發(fā)生時一直在進行壓裂作業(yè)（Leiet al，2019b）。形變區(qū)A位于頁巖氣開采產(chǎn)生的主要形變區(qū)C的西北邊緣，距離C較近。通過查閱地震目錄及野外調(diào)查我們發(fā)現(xiàn)，形變區(qū)A與珙縣地震的位置相符。為了判斷形變區(qū)A主要源于水力壓裂還是珙縣地震彈性變形，我們利用珙縣地震前后2018 年12 月30 日和2019 年1 月11 日的Sentinel-1 降軌數(shù)據(jù)作干涉，并轉(zhuǎn)換成LOS 向地表形變場，如圖9b 所示。其形變分布清楚地揭示了一條近南北向斷層的彈性形變場，且LOS 向形變量值可達61 mm，而水力壓裂在12 天內(nèi)很難達到該量值的形變，空間分布也很難達到圖9b 所示的形變范圍，所以盡管H18 井位于形變區(qū)A內(nèi)，我們認為其主要是珙縣地震所致。同樣，形變區(qū)B主要是由興文地震斷層滑動引起，附近的開采井距離形變區(qū)均較遠。需要指出的是，我們在形變圖上觀測到的形變位置與Lei 等（2019b）給出的興文地震定位結(jié)果有較大偏差，而地震定位結(jié)果附近無明顯的地表形變（圖9a）。此外，經(jīng)精定位的筠連地震位置距離我們的研究區(qū)C較遠，且在形變圖中未觀測到明顯的地表形變，因此也不會對研究區(qū)C的形變有影響。通過以上分析，我們認為研究區(qū)C內(nèi)的形變分布特征主要反映了兩三年內(nèi)水力壓裂作業(yè)向目標(biāo)地層注水及返排操作引起的地表形變。圖9a 也很好地揭示了一些較大地震事件與頁巖氣開采形變之間的空間分布關(guān)系，對于這些地震的誘發(fā)模型分析非常重要。

圖9 MS≥5.0 地震引起的地表形變與頁巖氣開采區(qū)形變的關(guān)系（a） 2017 年6 月12 日至2019 年7 月8 日期間的ALOS-2 降軌頁巖氣區(qū)塊形變場，圖中矩形A 和B 分別為珙縣地震和興文地震引起的地表形變，區(qū)域C 為頁巖氣開采引起的地表形變，地震定位結(jié)果來自Lei 等（2019b），H18 井位置源于本研究實地野外調(diào)查；（b） 2018 年12 月30 日至2019 年1 月11 日期間的Sentinel-1 LOS 向InSAR 形變圖Fig. 9 The relationship between the surface deformation caused by the MS≥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 數(shù)據(jù)得到結(jié)果的正確性并提高形變的觀測精度，我們利用2015 年至2019 年獲取的Sentinel-1 SAR 數(shù)據(jù)進行時序PS-InSAR 處理，結(jié)果顯示兩個傳感器得到的形變特征有很好的一致性。雖然采用了不同的數(shù)據(jù)和不同的處理方法，地表觀測到的形變特征都能反映頁巖氣開采造成的地面形變及其變化過程。結(jié)合四種不同觀測角度的數(shù)據(jù)，我們也可以分析得到地表形變的主要活動方式，初步揭示出頁巖氣生產(chǎn)過程中造成地面形變的復(fù)雜性。Sentinel-1 數(shù)據(jù)提供了更好的觀測精度，也提供了頁巖氣區(qū)塊形變的一些細節(jié)特征。通過形變與地震反射剖面的比較，我們初步確認在開采區(qū)南側(cè)邊界存在一條較寬的東西向斷層破碎帶，其北界如果對應(yīng)于一條隱伏斷層，在持續(xù)的開采和注液條件下，極有可能成為具有較高誘發(fā)地震風(fēng)險的地震斷層，應(yīng)在未來的頁巖氣開采和誘發(fā)地震研究中予以重點關(guān)注。

ALOS-2 數(shù)據(jù)由JAXA （JAXA 2nd Research Announcement on the Earth Observations，ER2A2N171）提供，中國科學(xué)技術(shù)大學(xué)李俊倫教授提供了頁巖氣開采井位信息，研究中與日本地質(zhì)調(diào)查局雷興林教授、四川省地震局易桂喜研究員、中國地震局地球物理研究所房立華研究員、中國地震局地質(zhì)研究所魯人齊研究員進行了有益討論，審稿專家給出了建設(shè)性的修改意見，作者在此一并表示感謝。