鄒正波 李 輝 吳云龍 吳桂桔 康開軒

1 中國地震局地震研究所（地震大地測量重點實驗室），武漢市洪山側路40號，430071

2 中國地震局地殼應力研究所武漢科技創(chuàng)新基地，武漢市洪山側路40號，430071

3 武漢大學測繪學院，武漢市珞瑜路129號，430079

2015－04－25尼泊爾M8.1級地震發(fā)生在喜馬拉雅斷裂帶中段，是地震研究的熱點區(qū)域。重力觀測對研究巨大地震震源區(qū)內部物質變化具有重要意義。自20世紀60年代邢臺地震之后，我國就開始地面定點和流動重力測量，并在多次大地震測報中發(fā)揮重要作用［1－2］。然而地面重力測量受觀測條件限制存在空間分布不均勻和時間分辨率低等不足，難以追蹤大地震造成的大范圍重力變化及空間遷移過程，而大范圍、長時間、近連續(xù)的衛(wèi)星重力觀測成為地面觀測的有力補充。研究表 明［1，3－6］，GRACE衛(wèi)星有能力為大地震孕震階段的重力場變化提供證據。因此，本文利用GRACE重力衛(wèi)星觀測資料研究此次地震前震源區(qū)及周邊的重力變化特征，并結合球形位錯模型模擬該地震同震重力場變化，為研究該地震孕震過程和發(fā)震機理提供地球動力學依據。

1 衛(wèi)星重力數據與方法

1.1 衛(wèi)星重力數據

本文利用美國德州大學奧斯汀分?？臻g研究中心公布的2002－04～2015－01 共141個最大階次達60階的RL05月重力場模型數據，該模型扣除了大氣、海洋、固體潮、海潮及極潮等影響，其中部分月份由于缺少觀測數據而無重力場模型數據?？紤]到GRACE在C20項的不足，利用由5顆SLR衛(wèi)星按同樣時間長度解算得到的C20項替換。

為研究地球動力學相關信息，需分離水文因素的影響。選取全球陸面數據同化系統(tǒng)GLDAS（global land data assimilation system）模型1°×1°數據，將陸地水變化信息（土壤水、雪水及總植冠水）轉化為與GRACE數據相當的60階月水文模型（球諧系數）序列［7］，并將其從衛(wèi)星信號中扣除（GRACE－GLDAS）。

1.2 重力變化信號提取

為消除南北異常條帶的影響，對GRACE 數據選取400km 扇形濾波，并選取2003～2010的月重力場模型作為背景，計算該區(qū)域重力場的時空變化；而水文信號利用GLDAS月數據（對其進行球諧轉化后采用GRACE 數據相同的處理方法）加以扣除，以分析地球動力學對區(qū)域重力場的影響?？紤]到GRACE的空間分辨率及該地震孕震與中國大陸存在的聯(lián)系，故而選取20°～45°E、65°～110°N作為研究區(qū)域。

1.3 長期變化率、周期變化及時間序列重力變化

由于GRACE提取的重力變化時間序列反映了全球質量重新分布的過程，其中海洋、大氣及水文等重力場信號均具有周期性特征，且模型可能扣除不完全，可利用周期函數進行分析。將信號分為周期信號及線性信號兩部分，采用最小二乘法對未知參數進行求解，分離點位重力長期變化及周期性信號影響，為后續(xù)討論提供依據。

2 衛(wèi)星重力變化結果

2.1 年度重力變化

利用GRACE扣除GLDAS陸地水質量變化影響后的長期重力資料，獲得平均年度重力變化。以此為背景，得到2003～2014年累計重力場變化。圖1中紅點為相應年度發(fā)生的M7.0以上地震震中分布（表1）。

表1 2002－04～2015－04研究區(qū)大于M7.0地震信息Tab.1 Earthquakes with magnitude larger than M7.0 in study area during 2002April to 2015April

分析圖1、表1可以發(fā)現(xiàn)：

1）重力變化存在明顯的分時段特征。以2008年為界，前后區(qū)域重力變化特征不同。2003～2008年重力變化以增加為主，且7級以上地震很少，僅在2008－03、05發(fā)生兩個地震；2009～2014以重力減小占主導，7級以上地震發(fā)生5個。這一特征可能顯示地震對區(qū)域重力場的影響。

2）重力變化存在顯著的正負分區(qū)特征。①印度北部重力負變化顯著，重力值逐年減少，從5 μGal變化到－1μGal，再到－5μGal；②青藏高原東角重力變化也出現(xiàn)由正變負的特征，特別是緬甸弧地區(qū)，重力變化由－1μGal到－5μGal，2011年后開始有所恢復；③青藏高原中部（新疆西藏青海交界）自2009年開始出現(xiàn)持續(xù)增加現(xiàn)象；④板塊邊界成為這些區(qū)域分布的重要依據。

3）衛(wèi)星重力變化為地震預測的地點判定提供了大尺度的依據。已有的時變重力資料的震例研究表明，地震多發(fā)生在重力變化與主斷裂構造走向一致的重力梯度帶上，發(fā)生在重力正負變化的轉化帶上即重力變化的零值線附近［2］，如圖2。

4）重力變化為尼泊爾地震的孕震過程及孕震環(huán)境提供了有力證據。衛(wèi)星重力變化表征的是空間固定點的重力變化，不受地面高程影響，其變化來源于質量遷移，且與質量變化成正比。2014年重力變化顯示，尼泊爾地震震中位于重力正負變化的交界區(qū)域，其孕震背景為印度板塊沿北東方向以45mm／a的速率擠壓歐亞板塊。由該地震的震源機制可知，發(fā)震斷層的南部區(qū)域處于斷層的壓縮區(qū)，質量增加、重力為正；北部位于斷層的拉張區(qū)，重力減少。

圖1 2003～2014年尼泊爾及周緣重力變化Fig.1 Annual gravity variation in Nepal and its periphery

2.2 重力長期變化趨勢

對研究區(qū)GRACE－GLDAS結果進行時間序列分析，扣除年、半年周期重力變化后得到的線性變化率即為該點2002－04～2015－01 的重力長期變化率，如圖2。圖中紅點表示相應年度發(fā)生的7級地震震中位置（表1），黑點代表2015年M8.1尼泊爾地震震中，藍色線條為一級塊體及二級塊體邊界。

圖2 尼泊爾及周緣重力長期變化率Fig.2 Gravity long term change rate in Nepal and its periphery

由圖2（a）不難發(fā)現(xiàn)：1）在尼泊爾震中區(qū)兩側近東西向存在顯著的重力變化，印度北部地區(qū)重力變化最大，其最大速率為－0.7μGal／a，青藏高原東角負變化特征顯著；2）其他地區(qū)重力變化相對平穩(wěn)，青藏高原中部以0.2μGal／a的速度增加，新疆天山以－0.1μGal／a的速率減小，四川、貴州等地以正變化特征為主；3）結合斷層分布發(fā)現(xiàn)，多數斷層邊界區(qū)域重力變化均極小或者為0，例如南北地震帶地區(qū)。

年重力變化（圖1）具有分段特征，同時與地震有很好的對應關系，因此，將長期變化以2008－05 M8.0 汶川地震為界進行時間分段研究：1）2002－04～2008－04（圖2（b））青藏高原西南部、印度北部以－1μGal／a的速率減小，而川青甘交界及周邊重力增速也達到0.3μGal／a；2）2008－05～2015－01（圖2（c））發(fā)生M7.0以上地震7次，青藏高原東角及周邊的重力長期較小，尼泊爾地震震中位于重力變化高梯度帶的零值線附近；3）7 級以上地震震中空間分布及塊體邊界顯示，玉樹地震、蘆山地震、汶川地震、緬甸地震等均發(fā)生在塊體邊界且重力變化緩慢區(qū)域或零值線附近，而非極值區(qū)。

2.3 震中時間序列變化

以震中點的GRACE－GLDAS重力長期變化為基礎，扣除年、半年周期重力信號，給出點位時間序列。由于GRACE 空間分辨率小，因此該點時間序列反映的是所在區(qū)域重力變化特征。結果顯示，震源區(qū)在13a中基本保持重力減?。▓D3，藍色點線為時間變化序列，紅色曲線是藍色曲線的多項式擬合），峰值達5～－4μGal，且從2012年開始經歷震前重力增加－減少－增加的動力變化過程。鑒于唐山地震［8］、日本地震［9］在地震發(fā)生前數年經歷了重力減?。黾樱l(fā)震－減小的過程，雖然GRACE數據未能更新到4月，但其震前的基本特征已經具備，因此可將此規(guī)律應用到日后的地震預報實踐中。

圖3 震中點重力時間序列變化Fig.3 Time series gravity change of the epicenter

3 理論模擬同震重力場變化

選取USGS公布的有限斷層模型，該模型的發(fā)震斷層參數為：傾角為10°，斷層走向為295°，沿斷層走向與傾向分成11×11個子斷層，子斷層尺度為20km×15km。利用球形位錯理論模擬計算尼泊爾M8.1地震震源區(qū)及周邊的空間固定點同震重力場變化信息［10］。

賈民育等［11］認為，監(jiān)測8級地震孕震重力變化，要求地面重力測網的空間分辨率達到1 160 km；而褚慶忠等［1］總結多個地震發(fā)現(xiàn)，7～8級地震震前地表重力變化的范圍可達500km，時間長達0.5～3a?？紤]到地震的影響范圍及GRACE重力衛(wèi)星的分辨率，以震中點為中心區(qū)域（23°～33°E，80°～90°N），網格間距為0.1°×0.1°。模擬結果（圖4）證實，斷層重力場變化基本沿斷層兩側對稱分布，且在壓縮區(qū)增加（斷層下盤），拉張區(qū)（上盤）減少。

為研究GRACE月重力場模型能否探測到這一變化，將同震重力變化進行分析處理（圖5）。第一步，球諧展開后截斷成60×60的球諧系數；第二步，對球諧系數進行400km 扇形濾波處理。結果發(fā)現(xiàn)，通過球諧展開截斷，重力變化由原來的－50～70μGal減少至－0.7～0.8μGal（圖5）；而經第二步處理后，截斷重力變化降低了近一個數量級（圖5（b）），空間范圍也發(fā)生很大變化。兩步處理均使模型分辨率由原來的0.1°×0.1°變?yōu)?°×3°乃至更大，數值更小。因此，以現(xiàn)有的GRACE衛(wèi)星的月尺度全球重力場模型提取該同震變化是難以實現(xiàn)的。采用局部重力場方法改善現(xiàn)有模型的空間分辨率及觀測精度，或直接利用地震前后經過震中點附近的星間距離及其變率數據提取同震變化可能是更有效的途徑。此外，2017年發(fā)射的裝載有星間激光測距儀的GRACE Follow On衛(wèi)星也很值得期待。

圖4 位錯模型模擬的尼泊爾大地震空間固定點同震大地水準面及重力變化Fig.4 The coseismic geoid and gravity changes of the points fitting on the space simulated by dislocation model

圖5 處理后的同震重力變化Fig.5 Processed coseismic gravity changes

4 結 語

利用衛(wèi)星重力資料對2015尼泊爾地震震前區(qū)域重力場變化進行計算與分析，得出如下結論：

1）區(qū)域重力變化表現(xiàn)出顯著的分區(qū)分段特點，這些特點與地震發(fā)生及構造分布存在很好的對應關系，體現(xiàn)了地震對區(qū)域重力場的影響。

2）震中區(qū)附近M7.0 級以上歷史地震及尼泊爾M8.1地震震前衛(wèi)星重力結果顯示，地震基本發(fā)生在具有構造背景的斷層邊界重力變化梯度帶的過渡帶上，且位于重力變化的零值線附近。

3）震前衛(wèi)星重力變化及理論模擬的同震重力變化顯示出重力變化與動力學的關系，壓縮區(qū)重力增加，拉張區(qū)重力減少。

4）多個地震震例顯示，重力變化在震前1～2 a出現(xiàn)減?。黾樱l(fā)震－增加－減小的現(xiàn)象。

［1］褚慶忠，武澤，邵先杰，等.地震重力異常形成機理及其在地震預報中的作用［J］.地震工程學報，2014，36（1）：201－206（Chu Qingzhong，Wu Ze，Shao Xianjie，et al.Mechanism of Gravity Anomaly before Earthquake and Its Role in Earthquake Prediction［J］.China Earthquake Engineering Journal，2014，36（1）：201－206）

［2］祝意青.青藏高原東北緣強震前兆特征研究——流動重力方法［J］.國際地震動態(tài)，2007（5）：1－6（Zhu Yiqing.Precursory Characteristics of Stronger Earthquakes in Northeastern Edge of Qinghai－Tibet Plateau by Using Mobile Gravity Technique［J］.Recent Developments and World Seismology，2007（5）：1－6）

［3］鄒正波，羅志才，吳海波，等.日本Mw9.0地震前GRACE衛(wèi)星重力變化［J］.測繪學報，2012，41（2）：171－176（Zou Zhengbo，Luo Zhicai，Wu Haibo，et al.Gravity Changes Observed by GRACE before the Japan Mw9.0Earthquake［J］.Acta Geodaetica et Cartographica Sinica，2012，41（2）：171－176）

［4］鄒正波，羅志才，李輝，等.GRACE 探測強地震重力變化［J］.大地測量與地球動力學，2010，30（2）：6－9（Zou Zhengbo，Luo Zhicai，Li Hui，et al.Detection of Gravity Changes of Strong Earthquake Using GRACE［J］.Journal of Geodesy and Geodesy，2012，41（2）：171－176）

［5］王武星，顧國華，陳石.利用GRACE 觀測資料分析日本MW9.0地震前區(qū)域重力變化特征［J］.地震地質，2014，36（2）：523－535（Wang Wuxing，Gu Guohua，Chen Shi.Study on Regional Gravity Changes before the MW9.0Japan Eartqhuake Detected by GRACE［J］.Seismology and Geology，2014，36（2）：523－535）

［6］姜磊，李德慶，徐志萍，等.蘆山Ms7.0地震震前GRACE衛(wèi)星時變重力場特征研究［J］.地震學報，2014，36（1）：84－94（Jiang Lei，Li Deqing，Xu Zhiping，et al.Characteristics of Gravity Field before the Lushan Ms7.0 Earthquake Based on Time Variable Gravity Data from GRACE Satellite［J］.Acta Seismologica Sinica，2014，36（1）：84－94）

［7］Wahr J，Molenaar M，Bryan F.Time Variability of the Earth’s Gravity Field：Hydrological and Oceanic Effects and Their Possible Detection Using GRACE［J］.Journal of Geophysical Research：Solid Earth（1978－2012），1998，103（B12）：30 205－30 229

［8］李瑞浩，黃建梁，李輝，等.唐山地震前后區(qū)域重力場變化機制［J］.地震學報，1997，19（4）：399－407（Li Ruihao，Huang Jianliang，Li Hui，et al.The Mechanism of Regional Gravity Variation before and after the Tangshan Earthquake［J］.Acta Seismologica Sinica，1997，19（4）：399－407）

［9］鄒正波，羅志才，吳海波，等.日本Mw9.0級地震前GRACE衛(wèi)星重力變化［J］.測繪學報，2012，41（2）：171－176（Zou Zhengbo，Luo Zhicai，Wu Haibo，et al.Gravity Changes Observed by GRACE before the Japan Mw9.0 Earthquake［J］.Acta Geodaetica et Cartographica Sinica，2012，41（2）：171－176）

［10］付廣裕，孫文科.球體位錯理論計算程序的總體設計與具體實現(xiàn)［J］.地震，2012，32（2）：73－87（Fu Guangyu，Sun Wenke.Overall Design and Specific Structures of the Computing Codes for Coseismic Deformations on a Layered Spherical Earth［J］.Earthquake，2012，32（2）：73－87）

［11］賈民育，詹潔暉.中國地震重力監(jiān)測體系的結構與能力［J］.地震學報，2000，22（4）：360－367（Jia Minyu，Zhan Jiehui.The Structure and Ability of the China Seismological Gravity Monitoring System［J］.Acta Siesmologica Sinica，2000，22（4）：360－367）