利用球諧分析方法研究2017年九寨溝MS7.0地震前重力場變化

李忠亞 胡敏章 郝洪濤 李 輝

1 中國地震局地震研究所，武漢市洪山側(cè)路40號，430071 2 中國地震局地震大地測量重點實驗室，武漢市洪山側(cè)路40號，430071 3 湖北省地震局，武漢市洪山側(cè)路48號，430071

根據(jù)中國地震臺網(wǎng)測定，2017-08-08四川省西北部九寨溝縣發(fā)生MS7.0地震，震中位于33.2°N、103.82°E，震源深度約20 km。中國地震局地球物理研究所、美國地質(zhì)調(diào)查局和相關(guān)文獻(xiàn)計算結(jié)果均表明，本次九寨溝地震為走滑型地震[1-2]。九寨溝地震周邊地形和構(gòu)造復(fù)雜(圖1)，近10 a來，其南部和北部先后發(fā)生2008年汶川MS8.0地震、2013年蘆山MS7.0地震、2013年岷縣-漳縣MS6.6地震和2016年門源MS6.4地震等破壞性強(qiáng)震。本次地震的發(fā)震斷層為樹正斷裂[1,3]，其位于塔藏斷裂、岷江斷裂和虎牙斷裂西北延長線的交會區(qū)域，也是東昆侖斷裂帶由走滑到虎牙斷裂逆沖構(gòu)造的過渡區(qū)域[4]。

地震孕震過程中的質(zhì)量遷移必然伴隨著重力場的改變，震前重力場變化可以由GRACE重力衛(wèi)星和地表重力觀測得到[5-6]。通過絕對重力儀和相對重力儀組合的方式可以進(jìn)行流動重力觀測，不同期觀測得到的重力變化特征可以作為強(qiáng)震發(fā)生前的可靠前兆[7]。目前利用流動重力資料判斷地震發(fā)生的可能地點和時間主要依據(jù)重力場時空動態(tài)演化特征和相關(guān)經(jīng)驗信息。識別重力場動態(tài)演化規(guī)律時，目標(biāo)區(qū)域經(jīng)常會出現(xiàn)多個“正負(fù)重力變化過渡帶”等震前重力場變化特征，但受限于目前對地震孕育機(jī)制的認(rèn)識和觀測點位分布等因素，多數(shù)情況下難以對其進(jìn)行取舍。本文嘗試通過將重力變化數(shù)據(jù)進(jìn)行球諧展開并濾波，壓制局部重力場變化特征，突出主要演化規(guī)律，以便于目標(biāo)特征的識別。

F1：東昆侖斷裂帶;F2：阿萬倉斷裂；F3：塔藏斷裂;F4：光蓋山-迭山斷裂;F5：西秦嶺斷裂帶;F6：岷江斷裂;F7：虎牙斷裂;F8：龍門山斷裂帶;F9：鮮水河斷裂帶;F10：龍日壩斷裂圖1 青藏高原東緣地形、主要斷裂帶和城市分布Fig.1 The distribution of topography, main fault zones and cities in eastern Tibetan plateau

1 重力變化球諧展開方法

空間任意一點(θ,λ)的引力位球諧展開表達(dá)式[8]為:

(1)

(2)

重力儀在地球表面觀測到的重力變化δg為:

δg=Δg+βur

(3)

式中，Δg為空間固定點的重力變化，β為布格梯度，ur為徑向位移。采用球模型近似，空間固定點重力變化Δg和引力位變化ΔV有如下關(guān)系：

(4)

當(dāng)Δg已知時，可以直接進(jìn)行球諧展開：

(5)

式(4)和式(5)的球諧系數(shù)滿足：

(6)

2 重力變化計算結(jié)果

流動重力觀測結(jié)果可以很好地反映地震孕育過程中產(chǎn)生的重力變化[9-11]。本文重力變化數(shù)據(jù)來自中國地震局重力學(xué)科管理部，時間變化范圍為2015-04～2017-04。每期重力數(shù)據(jù)處理過程是采用絕對重力結(jié)果作為起算基準(zhǔn)，然后將絕對重力和相對重力進(jìn)行聯(lián)合平差，計算過程中采用絕對重力基準(zhǔn)對相對重力儀格值進(jìn)行解算以削弱格值誤差。不同期數(shù)據(jù)之間作差即可得到相應(yīng)時間段內(nèi)的原始重力變化。

為了獲取研究區(qū)域重力場動態(tài)變化圖像，需對重力變化數(shù)據(jù)進(jìn)行后處理，流程如下：1)剔除極少數(shù)誤差較大的數(shù)據(jù)，最終使用的數(shù)據(jù)點位平均精度優(yōu)于15 μGal。2)由于提供的重力變化數(shù)據(jù)是地表值，而根據(jù)上節(jié)方法計算處理時需要空間固定點重力變化值，為此進(jìn)行布格梯度歸算。布格梯度改正包含垂直形變效應(yīng)和布格層效應(yīng)。垂直梯度可取 -3.086 μGal/cm，地殼密度采用平均值2.67 g/cm3，由布格改正公式得到布格改正為1.1 μGal/cm。布格梯度是垂直改正和布格改正之和，大小為 -2.0 μGal/cm。研究區(qū)域每年隆升速率變化范圍在1～6 mm/a[12]，以最大隆升速率6 mm/a計算，每年因為垂直形變引起的重力變化絕對值不超過1.2 μGal，該數(shù)值小于相對重力儀觀測到的重力變化，因此本文不區(qū)分重力儀觀測的地表重力變化和空間固定點重力變化。3)實際數(shù)據(jù)處理時，選擇整個南北地震帶重力觀測數(shù)據(jù)(重力觀測點位分布見文獻(xiàn)[13])。將區(qū)域內(nèi)重力變化數(shù)據(jù)進(jìn)行格網(wǎng)化，區(qū)域外重力變化值設(shè)為0，根據(jù)式(5)進(jìn)行球諧分析，得到60階球諧系數(shù)。利用球諧系數(shù)進(jìn)行求和并作3°范圍平滑處理，得到的重力變化見圖2。

圖2顯示，九寨溝地震前2 a青藏高原東緣重力場動態(tài)變化劇烈。下面具體分析重力場變化特征：

圖2 青藏高原東緣重力場動態(tài)演化Fig.2 Gravity field changes in eastern Tibetan plateau

1)2015-04～2015-10期間，全區(qū)域絕大部分地區(qū)為正重力變化。在九寨溝、黑水、汶川、丹巴和雅安一帶是重力增加峰值區(qū)域，最大正重力變化達(dá)到20 μGal。除莊浪東北地區(qū)重力變化接近0外，其余地區(qū)重力變化幅度沿著峰值區(qū)域向西北和東南兩側(cè)遞減。峰值區(qū)域走向與龍門山斷裂帶走向基本一致，說明重力場變化與構(gòu)造運動密切相關(guān)。

2)2015-10～2016-04期間，大部分地區(qū)重力變化趨勢與2015-04～2015-10相反，表現(xiàn)為負(fù)重力變化。巴顏喀拉塊體東部和柴達(dá)木塊體的貴南至澤庫一帶是主要的重力減小地區(qū)，周邊的川滇塊體東北部和華南地塊西北部重力變化平緩，但在漳縣和西和東北地區(qū)出現(xiàn)正重力變化。區(qū)域整體重力場變化特征表現(xiàn)為，巴顏喀拉塊體東緣沿東北方向至莊浪附近，負(fù)重力變化逐漸過渡為正重力變化，形成負(fù)-正重力變化過渡帶。

3) 2016-04～2016-10期間，重力場變化最明顯的特點是巴顏喀拉塊體負(fù)重力變化在2015-10～2016-04期間變化基礎(chǔ)上繼續(xù)減小，并且負(fù)重力變化分布向西南延伸至川滇塊體東北部。阿壩、馬爾康至新龍一帶是負(fù)重力變化的峰值區(qū)域，最大負(fù)重力變化達(dá)到-35 μGal。在柴達(dá)木塊體和華南塊體內(nèi)部，重力變化分布特征分別是向北側(cè)和東側(cè)逐漸由負(fù)重力變化向正重力變化過渡，在貴南附近、資陽和遂寧東南區(qū)域出現(xiàn)重力正變化特征。區(qū)域整體重力場變化特征為川滇塊體東北部負(fù)重力變化向東部、東南部和北部逐漸過渡為正重力變化，形成重力變化過渡帶。

4)2016-10～2017-04期間，重力場變化在前2期基礎(chǔ)上出現(xiàn)反轉(zhuǎn)回調(diào)，具體表現(xiàn)為2個特征：一是負(fù)重力變化，在阿壩、若爾蓋、岷縣和漳縣西北地區(qū)，重力變化由0逐漸減小，在貴南附近數(shù)值達(dá)到 -20 μGal；二是正重力變化，在汶川、綿陽、平武和廣元一帶出現(xiàn)正重力變化峰值區(qū)域，最大重力變化為18 μGal，并且正重力峰值區(qū)域走向與龍門山斷裂帶走向一致。在重力變化峰值區(qū)域兩側(cè)，即西北方向至阿壩、若爾蓋、岷縣和漳縣一帶，東南方向至遂寧和資陽一帶，正重力變化幅度逐漸減小。區(qū)域重力場變化特征是貴南附近的負(fù)重力變化沿東南方向逐漸過渡至正重力變化，形成重力變化過渡帶。

3 重力變化與九寨溝地震

地震孕育過程中伴隨著質(zhì)量遷移和能量交換，可通過分析震前重力場動態(tài)演化特征來捕捉地震前兆信息。分析圖2可知，九寨溝地震前，震源周邊重力場演化可以分為3個階段。第1階段為重力增大階段，2015-04～2015-10巴顏喀拉塊體東部邊界重力場均表現(xiàn)為正變化，沿著邊界向四周正變化幅度逐漸減小。第2階段為重力減小階段，2015-10～2016-04重力場開始出現(xiàn)負(fù)變化，2016-04～2016-10重力場繼續(xù)減小，并且變化幅度增大。該階段內(nèi)，九寨溝地震震源均位于正負(fù)重力變化過渡帶上。第3階段為重力場恢復(fù)調(diào)整階段，巴顏喀拉塊體東部表現(xiàn)出與第2階段截然不同的重力變化，即重力增加。并且在該階段，震源仍然位于正負(fù)重力變化過渡帶上。仔細(xì)對比第2階段和第3階段的正負(fù)過渡帶可以發(fā)現(xiàn)，前者由負(fù)向正過渡的方向為東北向，后者為東南向，兩者方向發(fā)生約90°旋轉(zhuǎn)。正負(fù)重力變化過渡帶旋轉(zhuǎn)，表明巖石圈內(nèi)部物質(zhì)遷移劇烈。祝意青等[13]指出，九寨溝地震之前震源周邊重力場變化原因有2個：1)2013年四川蘆山7.0級地震和甘肅岷縣-漳縣6.6級地震使青藏高原東部巖石圈物質(zhì)重新分布；2)川西高原及鄰區(qū)深部殼幔邊界與上地幔物質(zhì)和能量強(qiáng)烈交換。九寨溝地震前重力場變化體征顯示，區(qū)域正負(fù)重力變化過渡帶空間走向發(fā)生旋轉(zhuǎn)后，緊接著發(fā)生了九寨溝地震。因此，捕捉正負(fù)重力變化過渡帶空間走向旋轉(zhuǎn)變化信息對于尋找地震可能發(fā)生的時間具有重要前兆意義。

本文分析重力場動態(tài)演化規(guī)律時，將經(jīng)重力網(wǎng)平差后的數(shù)據(jù)進(jìn)行球諧展開至60階，并進(jìn)行3°范圍平滑濾波，這樣處理后可以在圖2(b)、2(c)、和2(d)所示的3個時間段重力場變化圖像中清晰顯示重力正負(fù)變化過渡帶。3期正負(fù)重力變化過渡帶的交會區(qū)域基本集中在四川西北部，九寨溝地震震源即位于該區(qū)域內(nèi)部。破壞性強(qiáng)震孕育會引起重力場在較大范圍內(nèi)出現(xiàn)改變，本文中所用的球諧分析方法可以較好地突出大尺度重力場變化特征，抑制局部細(xì)節(jié)特征，這種特點有利于確定地震可能發(fā)生的大概位置。已有研究表明，強(qiáng)震易發(fā)生在重力變化正負(fù)梯度過渡帶的零值線附近[9,11,14-15]，本文給出的結(jié)果顯示，九寨溝MS7.0地震震源位于正負(fù)過渡帶內(nèi)，但與零值線存在一定距離，這與本文中采用的球諧分析技術(shù)突出主要重力場變化特征、觀測誤差、點位分布等因素有關(guān)。實際采用流動重力觀測資料進(jìn)行地震前兆分析時，建議采用重力網(wǎng)平差后直接得到的結(jié)果和球諧分析方法得到的結(jié)果同時進(jìn)行分析。

4 結(jié) 語

本文采用球諧分析方法研究九寨溝地震前重力場動態(tài)變化特征，取得的主要認(rèn)識有：

1)2015-04～2017-04期間，震源周邊區(qū)域重力場變化經(jīng)歷增大、減小和恢復(fù)調(diào)整3個階段，后2個階段在震源附近出現(xiàn)正負(fù)重力變化過渡帶，且重力變化過渡帶空間走向出現(xiàn)旋轉(zhuǎn)現(xiàn)象。重力場演變特征較好地反映了九寨溝地震的前兆現(xiàn)象。

2)本文使用的球諧分析方法，可以有效突出區(qū)域重力場演化的主要特征，便于地震前兆特征識別。但同時也會抑制局部細(xì)節(jié)變化，實際分析重力場變化數(shù)據(jù)時，建議與傳統(tǒng)方法同時使用。