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

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

0 引言

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

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

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

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

1.1 研究區(qū)域與實驗震例

青藏高原位于亞洲大陸中部(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ì)呈較強的黏塑性，且高原持續(xù)受到印度板塊向N的推擠作用力，使得下地殼不斷驅(qū)動著上地殼脆性層的運動與變形，由此進(jìn)一步形成了青藏高原地區(qū)復(fù)雜的孕震環(huán)境(李吉均等，1998)。

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

1.2 數(shù)據(jù)

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

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

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

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

2 方法

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

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

本方法具體包括3個步驟：

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

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

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

2.1 去除水文信息

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

(1)

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

(2)

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

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

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

(3)

(4)

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

MSH=λmax-λmin

(5)

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

2.3 重力異常的刻畫

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

(6)

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

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

3 實驗結(jié)果與分析

本研究以汶川地震和尼泊爾地震為例，根據(jù)GRACE擾動引力位增量的2階梯度計算得到了重力應(yīng)變張量，進(jìn)而獲得各像素在震前2a的重力最大切應(yīng)變時間序列，并進(jìn)一步通過偏移指數(shù)K探測時間序列中存在的異常。基于此，本研究在斷裂帶觀測尺度上完成了震前重力異常時空變化特征的分析，并以此探討震前斷裂帶可能存在的構(gòu)造活動。此外，為驗證重力異常為震前特有的現(xiàn)象，本研究還開展了非震年份的同期對比實驗。其中，非震年份的偏移指數(shù)值是利用GRACE衛(wèi)星月數(shù)據(jù)無缺失值的2005—2010年的同期數(shù)據(jù)計算得到的偏移指數(shù)K值的平均值，詳細(xì)計算過程為： 基于重力最大切應(yīng)變分別計算斷裂帶區(qū)域中每個0.5°×0.5°網(wǎng)格上在2005—2006年、 2007—2008年和2009—2010年的K值，然后可得到這3個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ā)震時間相隔在3個月以內(nèi)，因此將以上4個地震視作一個地震序列，詳細(xì)信息見表1。圖 3 顯示了震前2a龍門山斷裂帶區(qū)域的重力異常值(>3的K值)總和的變化曲線。從中發(fā)現(xiàn)，在汶川地震前3個月即地震序列震前半年，斷裂帶上出現(xiàn)了重力異常值迅速上升的現(xiàn)象，并達(dá)到了2年內(nèi)的最大值，由此推測在2008年2—4月期間，斷裂帶區(qū)域的構(gòu)造活動尤為活躍。為了更加清晰直觀地獲取震前重力異常演變規(guī)律，本研究還顯示了重力異常震前半年的時空分布情況，如圖 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ā)震時間

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

圖 4b 表明，在汶川MW7.9 地震及其2個余震(汶川MW6.1 地震和青川MW6.1 地震)發(fā)震前的1～3個月，也即為青川MW6.0 地震震前的4～6個月間，龍門山斷裂帶的中段和北段出現(xiàn)了大面積異常區(qū)，且異常區(qū)的面積為該地區(qū)上個時間段(即圖4a)異常面積的3倍。由此可知，此時龍門山斷裂帶的中段和北段出現(xiàn)了顯著的重力異?，F(xiàn)象，結(jié)合本研究計算得到的該時段的GRACE重力變化值為正值這一事實，可推測出該時間段里龍門山斷裂帶已進(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個月內(nèi)的重力異常值的空間分布，從圖中可知龍門山斷裂帶上再次出現(xiàn)了明顯的重力異?，F(xiàn)象，主要表現(xiàn)為高異常值沿斷裂帶走向分布，并多出現(xiàn)在斷裂帶北部。由圖4d 再結(jié)合圖4e 中的重力異常值的空間分布可推測，此時龍門山斷裂帶的主要構(gòu)造活動向N遷移，而1個月后在龍門山斷裂帶的東北端發(fā)生的青川MW6.0 地震驗證了我們的推測。這也表明相比于單個地震而言，針對地震序列的時空分析可以更加全面地認(rèn)識斷裂帶的構(gòu)造活動特征。圖 5 顯示，在非震年份的同期中，龍門山斷裂帶上出現(xiàn)的最大K值為2.8，并未出現(xiàn)重力異常。

圖 5 龍門山斷裂帶在非震年份的同期偏移指數(shù)K的時空變化結(jié)果Fig. 5 Temporal and spatial variation results of the simultaneous offset index K of the Longmenshan fault zone in non-earthquake years.圖中圓點的含義與圖 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ā)震時間

3.2 尼泊爾震例

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

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

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

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

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

4 討論

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

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

圖 9 汶川地震前2a龍門山斷裂帶上最大切應(yīng)變和擾動位的K值時間序列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的時間序列計算的龍門山斷裂帶區(qū)域里所有像素的K值(K>1、 2、 3)總和。b—d 基于重力值本身(相鄰月份擾動位增量)在震前2a的時間序列計算的龍門山斷裂帶區(qū)域里所有像素的K值(K>1、 2、 3)總和，其各自對應(yīng)的高斯濾波平滑半徑分別為250km、 300km和500km。圖中紅色箭頭表示汶川地震的發(fā)震時間。由于K值是由GRACE衛(wèi)星的3個相鄰月重力場數(shù)據(jù)計算得到，為了便于表示，在圖中的時間坐標(biāo)軸上僅標(biāo)注這3個月中最后1個月的時間

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

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

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

5 結(jié)論

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

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

(2)在震前，重力異常值呈現(xiàn)出沿斷裂帶分布的現(xiàn)象，而且大范圍地存在于震中附近，最大異常值出現(xiàn)在距離震中約50km的范圍內(nèi)，且上述顯著異常發(fā)生的時間為震前半年內(nèi)。

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