重力衛星觀測俯沖帶特大地震同震重力場變化

張蘇祥 盛書中 周 新

1 東華理工大學地球物理與測控技術學院，南昌市廣蘭大道418號，330013 2 應急管理部國家自然災害防治研究院，北京市安寧莊路1號，100085

因為同震信號會影響震后松弛衰減信號的估計，從重力衛星觀測時間序列中準確提取同震重力場變化，既是利用重力衛星數據研究震源機制的基礎，又是準確確定地幔粘滯度的基礎。有學者采用差分[1]、震前震后平均值差值[2]、堆疊[3]等方法從GRACE的Level 2數據(球諧系數)中提取同震重力變化，但這些方法沒有考慮震后松弛信號對同震信號的污染。本文采用美國得克薩斯空間中心(CSR)發布的RL06數據，通過時間序列分析方法提取2004年蘇門答臘MW9.3、2010年智利MW8.8和2011年日本東北MW9.0三次俯沖帶特大地震的同震重力變化、大地水準面變化、垂線偏差變化及重力梯度變化，并將所得結果與前人根據球體位錯理論模型計算的結果[1-8]進行對比。

1 GRACE時間序列分析

研究表明[9]，RL06模型數據的精度較RL05模型有明顯的提升，而不同的RL06模型數據中，CSR發布的RL06模型數據的階方差最小，其給出了完全規格化的球諧系數，最高階數為60。本文采用2003-01～2016-09 CSR發布的Level2 RL06 GSM重力場模型數據，受某些因素影響，缺失了16個月的數據。由于GRACE的軌道形狀對系數C20項不敏感，該項精度相對較低，本文采用人衛激光測距得到的C20項作為替代[2]。

GRACE重力衛星的軌道高度約為450 km，考慮到重力信號隨衛星與地球距離的增大而衰減，重力衛星僅獲取了重力場的中長波部分信號。由于受高頻噪聲的影響，使用GRACE Level2數據直接計算得到的重力場變化呈南北條帶效應，很難獲得有用的重力變化信息，需要對數據進行平滑處理。Werth等[10]的研究表明，考慮GRACE衛星軌道誤差的去相關性濾波器在水文研究中的濾波效果最佳，因此本文采用去相關性的DDK3濾波器進行平滑處理。

GRACE衛星重力場時間序列中包含了趨勢、季節、同震和震后等多種信號，需要從這些信息中提取出同震重力變化。本文采用最小二乘擬合時間序列分析方法提取同震信號，對于去除背景場后的GRACE月重力變化時間序列，可用式(1)表示：

(1)

式中，a0為常數項，a1為線性趨勢項，t為相對參考歷元時間差，ωk為振幅bk的周期項角頻率，φk為對應時期的相位，H(t-teq)為階躍函數，teq為地震時刻，τ為震后松弛時間(本文取5.0 a)，c和d分別為同震和震后重力變化幅度。利用每個月的觀測數據，使用最小二乘法求解各項系數。

令

(2)

(3)

式(1)的矩陣形式可寫作：

Ax=y

其法方程為：

ATAx=ATy

最小二乘解為：

x=(ATA)-1ATy

式(1)中的c即為同震重力變化。

2 GRACE同震重力場變化

2004-12-26蘇門答臘地震發生在印度洋板塊與緬甸微板塊(南亞板塊中的微板塊)之間；2010-02-27智利地震發生在南極洲板塊與美洲板塊之間；2011-03-11日本東北地震發生在亞歐板塊與太平洋板塊的交界處。利用上述數據及方法，在1°×1°的網格上計算2003-01～2016-09去除背景重力場(GGM03S)后震源區的大地水準面變化、同震重力變化、垂線偏差變化和重力梯度變化的時間序列，并利用時間序列分析提取同震重力場變化信號。

經DDK3平滑處理后，采用最小二乘擬合法分別得到3次地震的同震大地水準面和重力變化。為觀察這3次地震引起的重力場變化時間特征，在重力變化最大處選取A、B點進行時間序列分析。蘇門答臘、智利和日本東北地震的A、B點坐標分別為(2.75°N，94.25°E)和(7.25°N，97.25°E)、(39.25°S，75.75°W)和(35.25°S，69.25°W)及(38.75°N，138.75°E)和(34.75°N，143.75°E)，各點去除背景場后的同震變化時間序列見圖1～3。鑒于地震發生當月信號的復雜性，在時間序列擬合中排除了地震發生當月的數據。根據多項式模型，通過線性最小二乘法擬合了長期趨勢、季節、同震和震后引起的重力場變化，其中模型假定震后松弛時間為5.0 a。由圖1～3可見，GRACE觀測3次地震的結果均呈正負兩極分布，且負變化區(圖1(f)、2(f)、3(e))的同震效應明顯大于正變化區(圖1(e)、2(e)、3(f))。由時間序列可知，震后重力與大地水準面變化有明顯的指數衰減現象，主要是地幔粘彈性松弛效應引起的。蘇門答臘、智利和日本東北地震引起的同震大地水準面變化范圍分別為-5.9～0.8 mm、-3.0～0.8 mm 和-3.2～0.5 mm；同震重力變化范圍分別為-15.5～6.5 μGal、-9.1～2.1 μGal和-11.1～4.2 μGal(表1)。從圖1(e)中可見，除了2004-12，另一個較為顯著的同震信號是2012-04的M8.6和M8.2走滑地震事件引起的，盡管該信號影響2004年震后時間序列的擬合，但并未對2004年同震信號的提取產生干擾。

圖1 蘇門答臘地震同震變化

圖3 日本東北地震同震變化

垂線偏差對斷層滑動模型較敏感，特別是EW向的垂線偏差，可用于約束地震的震源參數[7-8]。經DDK3平滑處理后，采用最小二乘擬合時間序列分析方法提取3次地震的同震垂線偏差變化空間分布。由圖4可知，GRACE衛星可觀測到俯沖帶特大地震引起的同震垂線偏差變化，信號呈負-正-負或正-負-正三極分布。蘇門答臘、智利和日本東北地震事件引起的同震垂線偏差NS向變化范圍分別為-1.2～2.2 mas、-0.9～1.0 mas 和-1.1～1.4 mas；EW向變化范圍分別為-1.8～1.0 mas、-0.8～0.8 mas和-0.7～1.0 mas(表1)。

圖4 同震垂線偏差變化

計算并提取3次地震的同震重力梯度變化，得到圖5。由圖可知，同震重力梯度變化信號呈多極分布，如θλ分量有明顯的2對正負對稱四象限分布信號。相比于其他重力場物理量，重力梯度有更豐富的地震信號，蘇門答臘地震重力梯度rr分量(徑向)變化幅度為-0.6～0.4 mE，rθ分量為-0.5～0.3 mE，rλ分量為-0.2～0.4 mE，θθ分量(NS向)為-0.3～0.4 mE，θλ分量為-0.1～0.2 mE，λλ分量(EW向)為-0.3～0.2 mE；智利地震重力梯度rr分量變化幅度為-0.4～0.1 mE，rθ分量為-0.2～0.2 mE，rλ分量為-0.2～0.2 mE，θθ分量為-0.1～0.2 mE，θλ分量為-0.1～0.1 mE，λλ分量為-0.2～0.1 mE；日本東北地震重力梯度rr分量變化幅度為-0.5～0.3 mE，rθ分量為-0.3～0.3 mE，rλ分量為-0.3～0.2 mE，θθ分量為-0.2～0.3 mE，θλ分量為-0.1～0.1 mE，λλ分量為-0.2～0.2 mE。由此可知，重力梯度中rr分量的同震變化幅度最大，rθ分量次之。

3 討 論

收集理論模型計算的3次地震同震重力場變化，并將其與本文結果進行對比。由表1可見，利用球體位錯理論模型計算的同震重力場變化[1-8]與本文提取的同震重力場分布信號特征一致。對于一個低傾角的俯沖型地震，其上覆板塊的變形特征比俯沖板塊更顯著，即主要變形集中在上覆板塊[11](圖6)。因此，發生在上覆板塊的同震重力場信號減小的幅度遠大于俯沖板塊的重力場增加幅度，即可以由海底面、莫霍面升降及內部物質膨脹的模型來解釋[12]。震后重力與大地水準面變化均隨時間呈非線性增加，這主要是由震后地幔的粘彈性松弛效應引起的。

實線和虛線分別表示埋深為10 m和2 km的斷層

表1 同震重力場變化范圍

本文采用最小二乘擬合時間序列分析方法提取的同震重力變化范圍均大于前人研究結果[1-4]，其中蘇門答臘、智利、日本東北地震的差異分別為7 μGal、4 μGal、6 μGal。導致這些差異的因素可能是由于前人采用了CSR RL04或RL05數據，而本文采用更高精度的RL06數據；部分研究采用的是差分[1]、震前震后平均值差值[2]和堆疊[3]方法，而本文采用多項式擬合時間序列分析方法提取同震重力變化；有的研究采用350km高斯濾波器提取蘇門答臘地震重力變化[1]、采用300 km高斯濾波器提取智利地震同震重力變化[2]、采用扇形濾波器提取日本東北地震重力變化[3]，而本文均采用DDK3濾波器；震后數據長度不同，震后數據越長，利用時間序列分析方法提取的同震重力變化越精確。

由文獻[1-4]給出的球體位錯理論模型計算的同震大地水準面變化、重力變化結果與本文的GRACE觀測結果均呈正負兩極分布。由表1可見，日本東北地震大地水準面理論值[7]與觀測值的差異約為1 mm。在同震重力變化方面，蘇門答臘地震的理論值[1]和觀測值差異最大約為6.4 μGal，智利地震的理論值[2]與觀測值在重力減小區的差異小于 1 μGal，日本東北地震的理論值[4]和觀測值差異約為1 μGal。考慮到海洋、海潮和大氣等模型的不確定性，理論模型計算值與觀測值有1～2 μGal的差異是合理的，其中蘇門答臘地震理論值與觀測值產生較大差異的原因可能是球體位錯理論模型中采用的滑動斷層模型[12]不同。另外，利用位錯理論模型計算同震重力變化需要考慮海水質量重新分布的貢獻，如未作這一改正，會影響理論計算的結果。

智利地震同震垂線偏差變化經DDK3平滑后的觀測結果與文獻[8]在NS向的差異約為0.5 mas，EW向差異約為0.2 mas。日本東北地震同震垂線偏差變化經DDK3平滑后的觀測結果與文獻[7]在NS向的差異為0.7 mas、EW向差異為0.1 mas?？紤]到GRACE類型的重力衛星軌道為近NS向，因此沿NS向的觀測精度應高于EW向[7]。NS向垂線偏差產生的差異大于EW向，這是由于智利地震和日本東北地震在理論值計算時分別采用300 km高斯濾波器[8]和 350 km高斯濾波器[7]，而本文采用DDK3濾波器，因此認為該差異是合理的。

目前很少有計算同震重力梯度變化的理論模型，故無法比較本文時間序列分析方法得到的GRACE同震信號與理論信號的差異。從圖5給出的3個特大俯沖型地震產生的同震重力梯度變化分布來看，其信號特征較大地水準面、重力和垂線偏差更為復雜，即包含更多的信號特征，如θλ分量有明顯的2對正負對稱四象限分布信號，反映地震造成的地殼密度變化有一個明確的底部。因此，重力梯度可能對某些震源參數更敏感，需要有能夠計算同震和震后重力梯度變化的位錯模型及用其約束震源參數等方面的研究。

4 結 語

本文利用DDK3濾波器對GRACE重力衛星的觀測數據進行處理，利用最小二乘擬合時間序列分析方法成功提取了2004年蘇門答臘、2010年智利和2011年日本東北地震產生的同震重力場變化及其空間分布。蘇門答臘地震、智利地震和日本東北地震的同震重力變化范圍分別為-15.5～6.5 μGal、-9.1～2.1 μGal和-11.1～4.2 μGal，同震大地水準面變化范圍分別為-5.9～0.8 mm、-3.0～0.8 mm和-3.2～0.5 mm；垂線偏差NS向變化范圍分別為-1.2～2.2 mas、-0.9～1.0 mas和-1.1～1.4 mas；EW向變化范圍分別為-1.8～1.0 mas、-0.8～0.8 mas和-0.7～1.0 mas；重力梯度各分量中，rr分量的同震變化幅度最大，其次為rθ分量。3次地震的同震信號空間分布均表現為：1)同震大地水準面和重力變化信號呈非對稱兩極分布；2)垂線偏差呈負-正-負或正-負-正三極分布；3)重力梯度變化信號呈復雜的多極分布，如θλ分量有明顯的2對正負對稱四象限分布信號。通過與位錯理論模型計算結果進行比較，驗證了本文觀測結果與理論計算結果具有較好的一致性，為重力衛星數據在震源機制的計算應用方面提供了可靠的信號提取方法。本文提取的同震重力變化、大地水準面變化、垂線偏差變化和重力梯度變化可為利用重力衛星數據約束震源參數提供新的途徑，也可作為GRACE在其他領域(如水文、冰川等)應用的參考。

致謝：感謝中國科學技術大學胡曉輝博士對本文提出寶貴建議。