宋冬梅 王 慧 單新建 王 斌 崔建勇

1)中國(guó)石油大學(xué)(華東)，海洋與空間信息學(xué)院，青島 266580 2)中國(guó)石油大學(xué)(華東)，研究生院，青島 266580 3)中國(guó)地震局地質(zhì)研究所，北京 100029

0 引言

破壞性地震會(huì)造成國(guó)家經(jīng)濟(jì)建設(shè)和人民生命財(cái)產(chǎn)安全的巨大損失(Allenetal.，2019； Fanetal.，2019)。因此，地震監(jiān)測(cè)對(duì)于抗震設(shè)防具有重要的現(xiàn)實(shí)意義(Tanetal.，2019)。為此，世界上各個(gè)地震頻發(fā)的國(guó)家或地區(qū)均建立了多個(gè)地震監(jiān)測(cè)臺(tái)站(申俊等，2018)。然而，利用地面臺(tái)網(wǎng)傳統(tǒng)的地震監(jiān)測(cè)手段只能獲取有限區(qū)域范圍內(nèi)的地震前兆信息，難以從宏觀上反映孕震過程(Paveletal.，2017； Ziegeretal.，2018)。20世紀(jì)后期，遙感技術(shù)快速發(fā)展，使得大面積觀測(cè)和獲取實(shí)時(shí)信息得以實(shí)現(xiàn)，這有效地彌補(bǔ)了定點(diǎn)觀測(cè)等傳統(tǒng)地震監(jiān)測(cè)手段的不足(范一大等，2016； Zhaoetal.，2021)，遙感技術(shù)已經(jīng)成為地震監(jiān)測(cè)的重要手段(彭令等，2017； Piscinietal.，2017； Jiaoetal.，2018； Baiketal.，2019)。

Massonnet等(1994)采用衛(wèi)星雷達(dá)干涉測(cè)量技術(shù)獲得了1992年美國(guó)加利福尼亞州蘭德斯地震震后的地表形變信息。單新建等(2014)利用InSAR和GPS測(cè)量技術(shù)獲取了汶川地震垂直連續(xù)形變場(chǎng)，結(jié)果顯示發(fā)震斷層兩側(cè)的垂直形變衰減較快，橫跨斷裂帶的形變量>30cm，寬度≤50km，且沿?cái)鄬哟怪毙巫兏咧祬^(qū)主要集中分布在發(fā)震斷裂的南段、 中段和北端。以上提到的SAR技術(shù)和GPS測(cè)量雖然能夠很好地監(jiān)測(cè)地表的地殼運(yùn)動(dòng)，但是卻難以探測(cè)地殼深部的物質(zhì)運(yùn)移和質(zhì)量變化(王武星等，2014)。而地球重力場(chǎng)可反映地球的內(nèi)部結(jié)構(gòu)和組成，被視作是固體地球動(dòng)力過程的歷史再現(xiàn)(Sunetal.，2017)。GRACE衛(wèi)星作為重力衛(wèi)星技術(shù)的代表，其獲得的重力場(chǎng)數(shù)據(jù)包含了地球的密度結(jié)構(gòu)信息，可用于震前板塊間和板塊內(nèi)部作用機(jī)制等方面的研究(Peidouetal.，2019； Velicognaetal.，2020)。Wahr等(1998)在GRACE衛(wèi)星發(fā)射之前便已給出了時(shí)變重力場(chǎng)的基本理論，并進(jìn)行了數(shù)值模擬計(jì)算。Sun等(2004)基于位錯(cuò)理論就GRACE衛(wèi)星能否探測(cè)到同震重力變化開展了研究，得出結(jié)論： GRACE重力衛(wèi)星能夠探測(cè)出7.5級(jí)以上地震的同震形變信號(hào)。Han等(2006)基于GRACE RL04 Level-1數(shù)據(jù)對(duì)2004年蘇門答臘地震進(jìn)行了研究，獲得了世界上第一個(gè)由重力衛(wèi)星數(shù)據(jù)探測(cè)得到的同震重力變化。

雖然GRACE衛(wèi)星能夠探測(cè)到大地震的重力變化，但由于重力場(chǎng)模型中的位系數(shù)含有較大噪聲，使得反演結(jié)果會(huì)出現(xiàn)嚴(yán)重的南北條帶狀波紋現(xiàn)象，因此許多學(xué)者提出如高斯濾波等多種濾波方法對(duì)GRACE衛(wèi)星數(shù)據(jù)進(jìn)行去噪(Guetal.，2017； 丁一航等，2018； 郭飛霄等，2018； Chaoetal.，2019； 王陳燕等，2019； 劉瀟，2021)。然而，進(jìn)行濾波處理后發(fā)現(xiàn)，盡管噪聲可被去除，但同時(shí)也造成有用信息的丟失，使得地震重力場(chǎng)異常信息的反演精度下降(Peidouetal.，2020)。針對(duì)目前基于GRACE衛(wèi)星數(shù)據(jù)的重力異常信息提取技術(shù)的不足，本研究提出了一種基于最大切應(yīng)變的震前重力異常信息提取方法，該方法在抑制噪聲的同時(shí)，減少了有用信息的損失。本文以2008年汶川地震和2015年尼泊爾地震為例，基于GRACE衛(wèi)星數(shù)據(jù)反演的重力最大切應(yīng)變進(jìn)行了斷裂帶的震前構(gòu)造活動(dòng)分析。

1 研究區(qū)域及相關(guān)數(shù)據(jù)介紹

1.1 研究區(qū)域與實(shí)驗(yàn)震例

青藏高原位于亞洲大陸中部(26°00′～39°47′N，73°19′～104°47′E)，是世界上海拔最高、 面積最大且地形最復(fù)雜的高原(張培震等，2003)。新生代以來，印度板塊與歐亞板塊通過碰撞使得青藏高原大規(guī)模地殼隆起，并由此產(chǎn)生了厚約80km的地殼(Cuietal.，2001)。由于地殼內(nèi)的升溫導(dǎo)致其內(nèi)部介質(zhì)呈較強(qiáng)的黏塑性，且高原持續(xù)受到印度板塊向N的推擠作用力，使得下地殼不斷驅(qū)動(dòng)著上地殼脆性層的運(yùn)動(dòng)與變形，由此進(jìn)一步形成了青藏高原地區(qū)復(fù)雜的孕震環(huán)境(李吉均等，1998)。

青藏高原內(nèi)部存在較多斷裂帶(圖 1)，是全球陸地強(qiáng)震爆發(fā)的主體區(qū)。21世紀(jì)以來，該區(qū)域的地質(zhì)構(gòu)造運(yùn)動(dòng)十分活躍，且強(qiáng)震頻發(fā)(張培震等，2003； 徐錫偉等，2017)。據(jù)統(tǒng)計(jì)，在2000—2021年期間，該區(qū)域共發(fā)生了43個(gè)6級(jí)及以上強(qiáng)震，其中包括2次7.5級(jí)以上的大地震，即MW7.9 汶川地震和MW7.8 尼泊爾地震，各自對(duì)應(yīng)的發(fā)震斷裂帶為龍門山斷裂帶和南喜馬拉雅造山帶，如圖 1 所示(1)http： ∥earthquake.usgs.gov/earthquakes/search/。。其中，龍門山斷裂帶處于青藏高原東部，是一條呈NE-SW向延伸的構(gòu)造帶，長(zhǎng)度>500km，寬30～50km(Dirksetal.，1994； Xuetal.，2009； Zhangetal.，2009； 孫闖，2017)。喜馬拉雅主碰撞造山帶是一條向S突出的EW向弧形構(gòu)造帶，長(zhǎng)約2500km，寬300～500km，屬于青藏高原的南邊界，并位于印度板塊和歐亞板塊的交界處(Burgetal.，1997)。需要說明的是，汶川地震和尼泊爾地震爆發(fā)后的3個(gè)月內(nèi)，其發(fā)震斷裂帶上還發(fā)生了多個(gè)6級(jí)以上地震，各地震的詳細(xì)信息見表1。

圖 1 龍門山斷裂帶(黃色矩形框①)和南喜馬拉雅造山帶(黃色矩形框②)在青藏高原上的位置Fig. 1 The location of the Longmenshan fault zone(yellow rectangle ①) and the southern Himalayan orogenic belt(yellow rectangle ②)in Tibetan plateau.F1喜馬拉雅主沖帶； F2喀喇昆侖-嘉里斷裂帶； F3瑪尼-玉樹-鮮水河斷裂帶； F4昆侖-瑪沁斷裂帶； F5阿爾金-海原斷裂帶； F6金沙江-紅河斷裂帶。上述斷裂帶將青藏高原劃分為6個(gè)不同的活動(dòng)塊體(Taylor et al.，2009)，分別為拉薩塊體 (B1)、 羌塘塊體(B2)、 巴顏喀拉塊體(B3)、 柴達(dá)木塊體(B4)、 祁連塊體(B5)和川滇塊體(B6)

1.2 數(shù)據(jù)

1.2.1 GRACE數(shù)據(jù)

GRACE衛(wèi)星由德國(guó)空間飛行中心DLR和美國(guó)宇航局NASA聯(lián)合設(shè)計(jì)研制，并于2002年3月成功發(fā)射(Inceetal.，2019)。GRACE由2顆低軌衛(wèi)星組成，工作模式采用衛(wèi)星跟蹤衛(wèi)星模式，包括低低模式(SST-LL)和高低模式(SST-HL)。本研究所用的重力數(shù)據(jù)是由美國(guó)得克薩斯大學(xué)空間研究中心(CSR)于2018年4月公開發(fā)布的GRACE RL06版本Level-2產(chǎn)品(2)http： ∥www2.csr.utexas.edu/grace/RL06.html。，該產(chǎn)品包含了2002年4月—2017年7月共計(jì)163個(gè)月的重力場(chǎng)數(shù)據(jù)，重力場(chǎng)模型的最大階次為96，其時(shí)間分辨率為1個(gè)月，空間分辨率約為300km。需要說明的是，該數(shù)據(jù)已經(jīng)消除了由潮汐及非潮汐引起的大氣和海洋質(zhì)量變化方面的影響(Swensonetal.，2008)。由于GRACE為極軌衛(wèi)星，衛(wèi)星軌道傾角較大，使得GRACE重力場(chǎng)模型的2階項(xiàng)精度較低，故本文采用CSR提供的衛(wèi)星激光測(cè)距(SLR)的測(cè)量值對(duì)2階項(xiàng)進(jìn)行了替換(Chengetal.，2004)。

1.2.2 GLDAS數(shù)據(jù)

全球陸地?cái)?shù)據(jù)同化系統(tǒng)(GLDAS)由美國(guó)宇航局(NASA)、 美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(NCEP)和國(guó)家海洋大氣局(NOAA)聯(lián)合研制(Spennemannetal.，2015)。GLDAS采用數(shù)據(jù)同化技術(shù)將衛(wèi)星觀測(cè)數(shù)據(jù)和地表觀測(cè)數(shù)據(jù)整合到統(tǒng)一的模型中。模型以網(wǎng)格數(shù)據(jù)的形式集成了多種陸地表面場(chǎng)信息，如全球的降水量(降雨和降雪)、 蒸散(土壤水分蒸發(fā)和植被蒸騰)、 地表徑流、 地下徑流、 土壤濕度、 地表溫度和地表熱流等。本研究采用的GLDAS-2.1模型數(shù)據(jù)(3)https： ∥ldas.gsfc.nasa.gov/gldas/。包括地下0～2m土壤含水量和雪水，其時(shí)間分辨率為1個(gè)月，空間分辨率為1°。

2 方法

為盡可能地凸顯出構(gòu)造活動(dòng)信息，本文提出了一種基于GRACE衛(wèi)星數(shù)據(jù)的震前重力異常信息提取方法，通過最大切應(yīng)變對(duì)重力異常信息進(jìn)行表征。該方法在抑制GRACE重力衛(wèi)星數(shù)據(jù)中存在的南北條帶誤差的同時(shí)，還避免了與構(gòu)造活動(dòng)相關(guān)的重力信號(hào)的丟失。該方法的詳細(xì)流程如圖 2 所示。

圖 2 本文方法的流程圖Fig. 2 Flow chart of the proposed method.

本方法具體包括3個(gè)步驟：

(1)去除水文信息季節(jié)性變化的影響。GRACE衛(wèi)星觀測(cè)到的重力場(chǎng)變化主要由陸地水儲(chǔ)量變化、 構(gòu)造變形和地下物質(zhì)流動(dòng)等因素引起，其中陸地水儲(chǔ)量變化引起的重力變化相對(duì)較大(Deggimetal.，2021)。因此，首先需要去除由水文的季節(jié)性變化引起的重力場(chǎng)變化，得到去除水文干擾的擾動(dòng)引力位增量。具體過程見2.1節(jié)。

(2)計(jì)算最大切應(yīng)變。本研究通過引入重力最大切應(yīng)變來描述與構(gòu)造活動(dòng)相關(guān)的信息。在步驟(1)的基礎(chǔ)上，計(jì)算了擾動(dòng)引力位增量的2階梯度，并進(jìn)一步根據(jù)重力應(yīng)變張量獲得了重力最大切應(yīng)變，以凸顯出與構(gòu)造活動(dòng)相關(guān)的重力變化。其中，擾動(dòng)引力位的2階梯度對(duì)南北條帶誤差有明顯的壓制作用。具體過程見2.2節(jié)。

(3)基于偏移指數(shù)K進(jìn)行震前異常信息時(shí)空分析。為進(jìn)一步探究震前重力異常的時(shí)空演變規(guī)律，需要將最大切應(yīng)變?cè)跁r(shí)間域上進(jìn)行標(biāo)準(zhǔn)化，得到偏移指數(shù)K，并在空間域上對(duì)K值的分布特征進(jìn)行分析。具體過程見2.3節(jié)。

2.1 去除水文信息

GRACE時(shí)變重力場(chǎng)可以通過球諧系數(shù)恢復(fù)。球諧系數(shù)實(shí)際上是擾動(dòng)位常微分方程的解，滿足以下拉普拉斯方程(Eshaghetal.，2012)：

(1)

然而，GRACE重力衛(wèi)星觀測(cè)到的地殼表層質(zhì)量變化信息是由地質(zhì)構(gòu)造活動(dòng)以及水文的季節(jié)性變化共同作用所引起的。同時(shí)，陸地水儲(chǔ)量的季節(jié)性變化(如四川盆地在夏季降雨量偏多)引起的重力變化比較顯著(王武星等，2014)。因此，本研究利用GLDAS水文模型提供的土壤水和雪水?dāng)?shù)據(jù)可實(shí)現(xiàn)部分水文信號(hào)的分離。為此，首先將GLDAS數(shù)據(jù)進(jìn)行球諧系數(shù)展開并截止到96階(與GRACE的階數(shù)相同)，獲得水文球諧系數(shù)Cgldas和Sgldas，并與GRACE數(shù)據(jù)的球諧系數(shù)Cgrace、Sgrace分別作差，得到去除土壤水和雪水影響后的球諧系數(shù)Ccorrect和Scorrect。然后，對(duì)上述得到的相鄰月份的球諧系數(shù)進(jìn)行差分處理，以進(jìn)一步濾除水文信號(hào)的季節(jié)性變化的影響：

(2)

設(shè)ΔCti、 ΔSti和ΔCti+1、 ΔSti+1為經(jīng)過差分處理后的相鄰月差分球諧系數(shù)，將它們分別代入式(1)中即可得到各自對(duì)應(yīng)的擾動(dòng)引力位增量(即重力變化量)ΔTi和ΔTi+1。

2.2 重力最大切應(yīng)變

去除水文干擾以后的重力變化包含地殼運(yùn)動(dòng)、 變形和物質(zhì)運(yùn)移活動(dòng)信息。為最大限度地凸顯出與構(gòu)造活動(dòng)相關(guān)的重力異常信息，本研究引入重力最大切應(yīng)變。該應(yīng)變通過重力應(yīng)變張量導(dǎo)出，而重力應(yīng)變張量是基于拉格朗日方法由擾動(dòng)引力位增量的2階梯度計(jì)算得到的，其對(duì)南北條帶誤差有明顯的壓制作用(Lietal.，2015； Daietal.，2016)。

(3)

(4)

上述重力應(yīng)變張量E的最大特征值與最小特征值之差即為最大切應(yīng)變MSH：

MSH=λmax-λmin

(5)

MSH在地質(zhì)學(xué)中的物理意義為地殼某點(diǎn)處在應(yīng)力作用下產(chǎn)生的切應(yīng)變的最大值(許才軍等，2006)，因此本研究將重力最大切應(yīng)變作為描述構(gòu)造活動(dòng)強(qiáng)弱的表征量。

2.3 重力異常的刻畫

為進(jìn)一步提取出震前異常變化，更好地認(rèn)識(shí)斷裂帶的孕震過程，本文借助偏移指數(shù)K對(duì)最大切應(yīng)變時(shí)間序列的數(shù)據(jù)進(jìn)行刻畫。通過設(shè)定偏移指數(shù)K的閾值，實(shí)現(xiàn)各像素點(diǎn)位置處不同時(shí)間序列的異常探測(cè)。偏移指數(shù)K的計(jì)算公式為

(6)

式中，(x，y)表示像素位置，t為序列數(shù)據(jù)的時(shí)間維，MSH為重力最大切應(yīng)變值，μMSH為MSH數(shù)據(jù)序列的均值，σMSH為MSH數(shù)據(jù)序列的標(biāo)準(zhǔn)差。

按照統(tǒng)計(jì)學(xué)的3倍標(biāo)準(zhǔn)差原則，當(dāng)一組數(shù)據(jù)中的某個(gè)數(shù)值與均值的差值大于3倍標(biāo)準(zhǔn)差時(shí)，則這個(gè)數(shù)據(jù)被判定為異常數(shù)據(jù)(李光強(qiáng)，2009)。因此，根據(jù)式(6)可知，>3的K值即為最大切應(yīng)變的異常數(shù)據(jù)。需要注意的是，偏移指數(shù)值的計(jì)算針對(duì)圖像中各像素點(diǎn)在其時(shí)間序列上的變化，與空間域上的其他像素點(diǎn)無關(guān)。此外，在對(duì)各像素點(diǎn)的時(shí)間序列進(jìn)行異常探測(cè)的基礎(chǔ)上，還對(duì)每個(gè)像素點(diǎn)的K值在空間域上的分布特點(diǎn)進(jìn)行觀察，從而得到對(duì)異常時(shí)空演變規(guī)律的認(rèn)識(shí)。

3 實(shí)驗(yàn)結(jié)果與分析

本研究以汶川地震和尼泊爾地震為例，根據(jù)GRACE擾動(dòng)引力位增量的2階梯度計(jì)算得到了重力應(yīng)變張量，進(jìn)而獲得各像素在震前2a的重力最大切應(yīng)變時(shí)間序列，并進(jìn)一步通過偏移指數(shù)K探測(cè)時(shí)間序列中存在的異常?；诖?，本研究在斷裂帶觀測(cè)尺度上完成了震前重力異常時(shí)空變化特征的分析，并以此探討震前斷裂帶可能存在的構(gòu)造活動(dòng)。此外，為驗(yàn)證重力異常為震前特有的現(xiàn)象，本研究還開展了非震年份的同期對(duì)比實(shí)驗(yàn)。其中，非震年份的偏移指數(shù)值是利用GRACE衛(wèi)星月數(shù)據(jù)無缺失值的2005—2010年的同期數(shù)據(jù)計(jì)算得到的偏移指數(shù)K值的平均值，詳細(xì)計(jì)算過程為： 基于重力最大切應(yīng)變分別計(jì)算斷裂帶區(qū)域中每個(gè)0.5°×0.5°網(wǎng)格上在2005—2006年、 2007—2008年和2009—2010年的K值，然后可得到這3個(gè)K值的平均值K0，最后選取與震期相同月份的K0表征非震年份的同期狀態(tài)。

3.1 汶川震例

考慮到2008年5月12日、 5月25日及8月5日的MW6.1 、MW6.1 和MW6.0 3次余震和汶川MW7.9 主震皆發(fā)生在龍門山斷裂帶上，且發(fā)震時(shí)間相隔在3個(gè)月以內(nèi)，因此將以上4個(gè)地震視作一個(gè)地震序列，詳細(xì)信息見表1。圖 3 顯示了震前2a龍門山斷裂帶區(qū)域的重力異常值(>3的K值)總和的變化曲線。從中發(fā)現(xiàn)，在汶川地震前3個(gè)月即地震序列震前半年，斷裂帶上出現(xiàn)了重力異常值迅速上升的現(xiàn)象，并達(dá)到了2年內(nèi)的最大值，由此推測(cè)在2008年2—4月期間，斷裂帶區(qū)域的構(gòu)造活動(dòng)尤為活躍。為了更加清晰直觀地獲取震前重力異常演變規(guī)律，本研究還顯示了重力異常震前半年的時(shí)空分布情況，如圖 4 所示。圖中最大的K值高達(dá)4.6，如圖4b 中的紫色框所示。

表 1 本研究中各震例的詳細(xì)信息(數(shù)據(jù)來源： USGS官網(wǎng))Table1 Details of each earthquake case involved in this study(data source： USGS)

圖 3 龍門山斷裂帶上震前2a的重力異常值總和Fig. 3 The sum of gravity anomalies on the Longmenshan fault zone two years before the Wenchuan earthquake.紅色虛線表示 MW7.9 汶川地震的發(fā)震時(shí)間

圖 4 龍門山斷裂帶的偏移指數(shù)K在震前的時(shí)空變化結(jié)果Fig. 4 Pre-earthquake spatio-temporal variation of the offset index K on the Longmenshan fault zone.該圖像的像素是對(duì)去除水文影響后的擾動(dòng)引力位進(jìn)行插值所得的0.5°×0.5°的結(jié)果。圖幅中的最大值位于紫色矩形框所在位置。紅色標(biāo)注日期為發(fā)震期，紅色圓點(diǎn)表示發(fā)震位置

圖 4b 表明，在汶川MW7.9 地震及其2個(gè)余震(汶川MW6.1 地震和青川MW6.1 地震)發(fā)震前的1～3個(gè)月，也即為青川MW6.0 地震震前的4～6個(gè)月間，龍門山斷裂帶的中段和北段出現(xiàn)了大面積異常區(qū)，且異常區(qū)的面積為該地區(qū)上個(gè)時(shí)間段(即圖4a)異常面積的3倍。由此可知，此時(shí)龍門山斷裂帶的中段和北段出現(xiàn)了顯著的重力異?，F(xiàn)象，結(jié)合本研究計(jì)算得到的該時(shí)段的GRACE重力變化值為正值這一事實(shí)，可推測(cè)出該時(shí)間段里龍門山斷裂帶已進(jìn)入應(yīng)力累積過程，積累了大量的應(yīng)變能，失穩(wěn)性已達(dá)到較高水平。而在此之前，僅在龍門山斷裂帶的南部區(qū)域出現(xiàn)少量異常值(圖4a)。值得注意的是，地震序列的震中位置皆在上述顯著異常區(qū)內(nèi)，其中MW6.1 和MW6.0 青川地震的震中位置距研究區(qū)內(nèi)最大切應(yīng)變的最大K值所在的位置均在50km范圍內(nèi)。

圖 4d 為汶川MW7.9 地震震后1個(gè)月內(nèi)的重力異常值的空間分布，從圖中可知龍門山斷裂帶上再次出現(xiàn)了明顯的重力異?，F(xiàn)象，主要表現(xiàn)為高異常值沿?cái)嗔褞ё呦蚍植?，并多出現(xiàn)在斷裂帶北部。由圖4d 再結(jié)合圖4e 中的重力異常值的空間分布可推測(cè)，此時(shí)龍門山斷裂帶的主要構(gòu)造活動(dòng)向N遷移，而1個(gè)月后在龍門山斷裂帶的東北端發(fā)生的青川MW6.0 地震驗(yàn)證了我們的推測(cè)。這也表明相比于單個(gè)地震而言，針對(duì)地震序列的時(shí)空分析可以更加全面地認(rèn)識(shí)斷裂帶的構(gòu)造活動(dòng)特征。圖 5 顯示，在非震年份的同期中，龍門山斷裂帶上出現(xiàn)的最大K值為2.8，并未出現(xiàn)重力異常。

圖 5 龍門山斷裂帶在非震年份的同期偏移指數(shù)K的時(shí)空變化結(jié)果Fig. 5 Temporal and spatial variation results of the simultaneous offset index K of the Longmenshan fault zone in non-earthquake years.圖中圓點(diǎn)的含義與圖 4 一致

圖 6 南喜馬拉雅造山帶上震前2a的重力異常值總和Fig. 6 The sum of gravity anomalies in the southern Himalayan orogenic belt two years before the MW7.8 Nepal earthquake.紅色虛線表示 MW7.8 尼泊爾地震的發(fā)震時(shí)間

3.2 尼泊爾震例

尼泊爾地震序列包含了6個(gè)強(qiáng)震，這6個(gè)地震的發(fā)震時(shí)間間隔不到1個(gè)月，詳細(xì)信息見表1。本研究計(jì)算了震前2年南喜馬拉雅造山帶區(qū)域的所有重力異常值總和，計(jì)算結(jié)果見圖 6。從圖中可發(fā)現(xiàn)，在震前6個(gè)月，區(qū)域上出現(xiàn)了重力異常值顯著升高的現(xiàn)象，并達(dá)到了2013—2015年間的最大值。圖 7 顯示了震前6個(gè)月南喜馬拉雅造山帶上重力異常值的空間分布。

圖 7 南喜馬拉雅造山帶上的偏移指數(shù)K的震前時(shí)空變化結(jié)果Fig. 7 Pre-earthquake spatio-temporal variation of the offset index Kin the southern Himalayan orogenic belt.標(biāo)注*的時(shí)間段內(nèi)存在GRACE原始數(shù)據(jù)缺失的現(xiàn)象，即無2014年12月的GRACE數(shù)據(jù)

如圖7b 所示，在MW7.8 尼泊爾地震的震前6個(gè)月，南喜馬拉雅造山帶區(qū)域出現(xiàn)了大面積異常區(qū)，占所在圖幅面積的三分之一。該異常區(qū)全部圍繞著6個(gè)震中位置展布，并且出現(xiàn)最大K值的位置距離MW7.8 尼泊爾地震的震中約50km。上述異常現(xiàn)象表明，這6個(gè)震中附近區(qū)域皆存在較高的應(yīng)力作用。由此推測(cè)，在2014年10月—2015年1月期間，印度板塊對(duì)青藏高原的擠壓作用明顯增強(qiáng)，斷裂帶兩側(cè)的構(gòu)造活動(dòng)異?；钴S。如圖7c 所示，沿著南喜馬拉雅造山帶出現(xiàn)了多個(gè)異常值，且尼泊爾地震序列里的6個(gè)震中皆位于高異常值出現(xiàn)的位置，結(jié)合該時(shí)段GRACE數(shù)據(jù)中的重力變化值為正這一事實(shí)，可以推斷出此時(shí)的斷裂帶上已完成了應(yīng)力積累。而在震時(shí)斷裂帶及其以南區(qū)域的GRACE重力變化為負(fù)值，由此推測(cè)斷裂帶所受的應(yīng)力達(dá)到臨界值后產(chǎn)生了破裂，前期積累的應(yīng)力得到釋放，并最終導(dǎo)致發(fā)生了尼泊爾地區(qū)的6個(gè)地震。與地震期不同的是，在非震年份的同期中，南喜馬拉雅造山帶上存在的最大K值為2.6，這表明該區(qū)域并未出現(xiàn)重力異?，F(xiàn)象(圖 8)。

圖 8 南喜馬拉雅造山帶在非震年份的同期偏移指數(shù)K的時(shí)空變化結(jié)果Fig. 8 Temporal and spatial variation results of the simultaneous offset index K of the southern Himalayan orogenic belt in non-earthquake years.圖中圓點(diǎn)的含義與圖 7 一致

綜上所述，基于地震序列能夠更好地揭示斷裂帶在震前的構(gòu)造活動(dòng)過程。研究發(fā)現(xiàn)，震前往往伴隨著強(qiáng)烈的斷裂帶構(gòu)造活動(dòng)，應(yīng)力集中現(xiàn)象明顯，具體表現(xiàn)為研究區(qū)中異常的最大值出現(xiàn)在距離震中約50km范圍內(nèi)，且該異常發(fā)生的時(shí)間往往在震前半年，這可能是一種地震前兆信號(hào)。特別值得注意的是，不同于以往研究中將GRACE RL06數(shù)據(jù)60階以后的球諧系數(shù)直接舍棄或者對(duì)其進(jìn)行簡(jiǎn)單濾波處理，本研究所采用的重力最大切應(yīng)變的計(jì)算方法較好地利用了該衛(wèi)星數(shù)據(jù)的中高階信息，與單純借助重力值本身相比更能靈敏地探測(cè)到地球內(nèi)部的密度變化信息，從而更好地凸顯出了由構(gòu)造活動(dòng)所引起的震前重力異常變化。

4 討論

以往對(duì)震前重力異常的研究方法通常直接基于重力值本身(即擾動(dòng)位增量)的變化進(jìn)行分析，但這種方法需要首先利用不同半徑的高斯濾波對(duì)南北條帶噪聲進(jìn)行消除，容易造成數(shù)據(jù)的失真與信息的丟失(Hasegawaetal.，2011； Zhaoetal.，2015)。而本研究所使用的重力最大切應(yīng)變是在擾動(dòng)位增量的基礎(chǔ)上進(jìn)一步通過計(jì)算其2階梯度所得到的，該梯度可以在不損失重力信息的同時(shí)又能夠達(dá)到抑制噪聲的效果(Lietal.，2015； Daietal.，2016)。為比較以上2種方法對(duì)構(gòu)造活動(dòng)探測(cè)能力的差別，選取龍門山斷裂帶為研究區(qū)，對(duì)基于最大切應(yīng)變的異常提取方法和傳統(tǒng)的基于擾動(dòng)位增量的異常提取方法所得的汶川地震震前2a時(shí)間序列里的總K值進(jìn)行了比對(duì)分析。該處的總K值為研究區(qū)內(nèi)每個(gè)像素的K值的總和。實(shí)驗(yàn)結(jié)果如圖 9 所示。

圖 9a 顯示了基于最大切應(yīng)變時(shí)間序列的研究區(qū)總K值的變化情況，可以發(fā)現(xiàn)>3的總K值在震前3個(gè)月內(nèi)迅速增大，結(jié)合此時(shí)K值的分布與斷裂帶的空間展布基本一致這一事實(shí)(圖4b)，上述現(xiàn)象極有可能是斷裂帶強(qiáng)烈的構(gòu)造活動(dòng)所引起的。圖9b、 9c和9d顯示了基于擾動(dòng)引力位增量時(shí)間序列的總K值變化，其各自對(duì)應(yīng)的高斯濾波平滑半徑分別為250km、 300km和500km。在這3張圖中可以發(fā)現(xiàn)，除了以250km為平滑半徑的高斯濾波所得的異常提取結(jié)果在震前2a有>3的K值出現(xiàn)外，其余半徑的濾波結(jié)果均未出現(xiàn)K>3的情況。對(duì)>1和>2的K值變化曲線進(jìn)行分析可以發(fā)現(xiàn)，經(jīng)過平滑半徑為250km和300km的高斯濾波處理后，擾動(dòng)位增量的總K值在震前半年出現(xiàn)了較為明顯的局部增長(zhǎng)現(xiàn)象(圖9b，c)，但隨著濾波半徑的進(jìn)一步增大，該現(xiàn)象隨之消失，如圖9d 所示，這說明隨著濾波半徑的增大，在去噪強(qiáng)度增加的同時(shí)也造成了信號(hào)的丟失。

圖 9 汶川地震前2a龍門山斷裂帶上最大切應(yīng)變和擾動(dòng)位的K值時(shí)間序列Fig. 9 Total K value of maximum shear strain and disturbance potential in the Longmenshan fault zone two years before the Wenchuan earthquake.a 基于最大切應(yīng)變震前2a的時(shí)間序列計(jì)算的龍門山斷裂帶區(qū)域里所有像素的K值(K>1、 2、 3)總和。b—d 基于重力值本身(相鄰月份擾動(dòng)位增量)在震前2a的時(shí)間序列計(jì)算的龍門山斷裂帶區(qū)域里所有像素的K值(K>1、 2、 3)總和，其各自對(duì)應(yīng)的高斯濾波平滑半徑分別為250km、 300km和500km。圖中紅色箭頭表示汶川地震的發(fā)震時(shí)間。由于K值是由GRACE衛(wèi)星的3個(gè)相鄰月重力場(chǎng)數(shù)據(jù)計(jì)算得到，為了便于表示，在圖中的時(shí)間坐標(biāo)軸上僅標(biāo)注這3個(gè)月中最后1個(gè)月的時(shí)間

表 2 震前2a最大切應(yīng)變和擾動(dòng)位時(shí)間序列中各自的最大K值Table2 Maximum K value in the time series of maximum shear strain and disturbance potential two years before the Wenchuan earthquake

此外，單個(gè)像素的重力最大切應(yīng)變時(shí)間序列里的最大K值為4.6，出現(xiàn)在震前3個(gè)月內(nèi)。而擾動(dòng)位增量時(shí)間序列的最大K值出現(xiàn)的時(shí)間在發(fā)震期的1a以前，且由濾波半徑為300km和500km的高斯濾波處理后的擾動(dòng)位增量時(shí)間序列里的最大K值皆<3，詳細(xì)信息見表2。此外，由本文方法計(jì)算得到的最大K值位于龍門山斷裂帶的北端，相比于經(jīng)典方法得到的最大K值，該值所在的位置最接近青川地震的震中。需要說明的是，表2 中最大K值所在像素的顏色取決于其值的大小，顏色標(biāo)尺參照?qǐng)D 5。

綜上所述，相比于常用的擾動(dòng)位方法，本研究基于最大切應(yīng)變的方法對(duì)重力場(chǎng)的異常信息更為敏感，因此對(duì)與構(gòu)造活動(dòng)相關(guān)的重力異常信息具有較強(qiáng)的探測(cè)能力。

5 結(jié)論

GRACE衛(wèi)星能夠探測(cè)到地震重力場(chǎng)的變化，因此被廣泛應(yīng)用于地震監(jiān)測(cè)研究中。然而，以往學(xué)者們?cè)谔崛〉卣鸬闹亓ψ兓畔r(shí)，均采用不同方法對(duì)GRACE衛(wèi)星數(shù)據(jù)進(jìn)行濾波處理，造成了重力場(chǎng)中高階信息的損失，因此很難提取到8級(jí)及以下地震的重力異常。對(duì)此，本研究提出了一種基于最大切應(yīng)變的震前重力異常信息提取方法，該方法既能夠避免重力信息的損失，又能夠達(dá)到抑制噪聲的效果，以凸顯與構(gòu)造活動(dòng)相關(guān)的重力異常信息。本文以汶川地震和尼泊爾地震為例，分析了龍門山發(fā)震斷裂帶的震前重力異常的時(shí)空分布特征。主要研究結(jié)論如下：

(1)本文所提出的基于最大切應(yīng)變的方法比傳統(tǒng)方法對(duì)重力場(chǎng)的異常信息提取能力更強(qiáng)，因此能夠提取到汶川地震和尼泊爾地震這2個(gè)8級(jí)以下地震的震前異常信息，這為重力數(shù)據(jù)應(yīng)用于地震監(jiān)測(cè)研究提供了一個(gè)新方法。

(2)在震前，重力異常值呈現(xiàn)出沿?cái)嗔褞Х植嫉默F(xiàn)象，而且大范圍地存在于震中附近，最大異常值出現(xiàn)在距離震中約50km的范圍內(nèi)，且上述顯著異常發(fā)生的時(shí)間為震前半年內(nèi)。

(3)通過地震期與非震期的對(duì)比分析可以發(fā)現(xiàn)，斷裂帶在震前出現(xiàn)了大面積的重力異常，而非震期則表現(xiàn)得較為平靜，未出現(xiàn)明顯的重力異常，表明斷裂帶在震前的構(gòu)造活動(dòng)比非震期更為活躍。