利用斷層應力積累探討大地震發生的地點
——以日本2011年東北9級大地震為例

謝周敏，蔡永恩，吉岡祥一，阿部大毅

1 國家自然災害防治院，北京 100085 2 北京大學地球與空間科學學院理論與應用地球物理研究所，北京 100871 3 日本神戸大學都市安全研究中心，神戸 657-8501 4 日本神戶大學研究生院行星科學系，神戶 657-8501

0 引言

地震的孕育是震源區一個復雜的物理、化學過程，它的發生是地下巖石突然破裂的結果.即使科學技術高度發達的今天，預報地震的地點、時間和大小(震級)仍然是一個公認的難題.其問題主要在于地球內部的“不可入性”、大地震的“非頻發性”和地球物理過程的復雜性(陳運泰，2009).

王仁(1997)認為：巖層斷裂的矛盾雙方,是促使巖層破裂的推動力同巖層抵抗破裂的阻力的矛盾.前者是地殼在各種外力作用下的應力分布,后者是巖層的強度.地震的孕育和發生就是這對矛盾對立統一的發展過程.因此，從力學的角度看，解決地震預報問題，需要研究震源巖石的強度和應力隨時間的變化過程.前者，一般通過巖石力學實驗可以得到不同深度、不同溫壓條件下的一些認識；后者由于受技術條件的限制，除了地表淺層外，很難進行震源深度處的直接測量，只能通過應力變化引起的各種前兆現象，例如地震活動、地表變形、地震波速異常、大地電磁場、地下水位和溫度場等，反演震源的運動學特征和動力學環境，進行間接預測.

地震發生地點的預測是地震預報首先要解決的問題.一個地震活躍的區域，其歷史地震活動圖像反映了這個區域應變能積累、釋放和轉移的過程.根據板塊學說，在地質時間尺度內，板塊運動的速率是穩定和連續的.板塊間的邊界是最明顯的地震活動帶，在適當長的時間間隔內，地震帶內的大地震的破裂區往往會連接、無重疊地展布于地震帶.如果地震帶內某個地區長期沒有發生大地震，那么，這個地區就被視為未來最有可能發生大地震的地區，通常把這個地區稱為地震空區.茂木(Mogi，1968)認為，空區是地震應變能的積累區，大地震是當震源區積累的應變能達到極限時，導致巖石突然破裂而發生的.

地震空區用來作為預測地震發生的位置，最初是由Fedotov(1965)提出來的.他把發生在日本—庫頁島—勘察加地震帶的淺源大地震的震源區畫出來后,發現在這些破裂區之間存在多年沒有破裂的區域，并且推斷這些區域就是未來大地震可能發生的位置.幾年后，在他所預報的一些地區先后發生了3次7.8級以上的大地震.

如何確定地震空區是預測板間大地震位置的關鍵.Kelleher等(1973)首次提出的一個準則是：空區是走滑或逆沖斷層地震帶的一部分，并且這部分至少30年沒有破裂.但有些板塊邊界地段的歷史地震記錄并沒有逆沖大地震，對于這樣的地區，歷史地震的平靜要么代表一個長期的平穩過程，要么是一種暫時現象，其中大地震的重復時間超過了歷史記錄的長度.因此，需要區分不同類型的地震空區.后來又對這一標準進行了補充，即把板塊邊界某些部位有一次或多次大地震的歷史記錄作為未來大地震在這個部位發生的證據.

預測板塊邊界地震帶大地震位置的地震空區方法，也被用于板塊內部的地震(陳章立等，1981；魏光興等，1981；陸遠忠等，1985；梁春濤等，2018).

利用地震空區預測大地震發生位置，方法雖然簡單，但是茂木(Mogi，1968)發現，當使用80 年的地震記錄時，日本太平洋沿岸對于大于 7 級的事件幾乎沒有空區.這與使用 30 年得到的結果存在顯著差別.因此，判斷利用多長時間的地震記錄得到的地震空區才是可靠的，至關重要.此外，由于板塊邊界存在無震滑動區，所以，空區不一定是未來大地震的發生區(Mogi，1968).

為了克服地震空區預測大地震位置的缺欠，早期地表三角測量方法被用來觀測大地震前的地表變形，并將計算得到的應變的顯著異常區作為可能的大地震位置.1923 年日本關東大地震前幾十年就發現東京灣西部的日本佐賀三原地區的構造應變比通常的大好幾倍，而在其他地區未發現任何顯著的應變(Fujii and Nakane,1979).這個異常是由于亞洲板塊和菲律賓海板塊間的深部界面的無震逆沖滑動與淺部閉鎖段加速下行，在地表產生的擠壓導致的.這個方法也被用于研究阿拉斯加Shumagin地震空區的應變積累(Savage and Lisowski,1986).

三角測量方法的優點是通過直接觀測地表變形，定量地得到地表應變，但是工作量大，測點稀疏、精度低，受地形影響較大.20世紀70年代全球定位系統(GPS)問世，以其全天候、高精度、自動化、高效益等顯著特點，給地表變形監測帶來了革命性的變化，使得人們能利用GPS數據直接研究地震的變形場和斷層的破裂過程.無論是真實的地震空區還是震前地表應變高異常區，都是由斷層閉鎖區引起的，因此利用GPS觀測數據直接反演地表變形和地震斷層的閉鎖區，就可能直接預測大地震的發生位置.

斷層閉鎖程度通常用斷層閉鎖區 “后滑(backslip)”模型(Savage,1983)反演地表變形得到.這個模型基于斷層彈性位錯理論.假設俯沖板塊沿其長度以匯聚速率均勻地向軟流圈滑動，通過在預期的閉鎖區施加與俯沖匯聚速率相反的滑動速率，同時在預期的斷層閉鎖區外的斷層指定零位錯，就可以實現斷層的閉鎖.由于穩定狀態位移解，理論上在地表不產生應變或地形變化，因此，一般模擬俯沖帶斷層閉鎖，只要在預期的地震斷層滑動部位上施加與同震位錯反方向的位錯(正斷層位錯)就可以了.這個反位錯也稱為滑動虧損，其值越大，意味著應變積累越大.這個方法需要預先人為指定反演地震斷層閉鎖區和板塊的運動方向，并且需要對位錯方向進行約束.有研究者利用這個方法和GPS數據反演了2011年日本東北大地震前斷層的閉鎖區(Ito et al.,2000;Nishimura et al.,2004;Suwa et al.,2006;Loveless and Meade,2010),得到滑動虧損深達40～100 km;使用層狀黏彈性模型反演得到的滑動虧欠深度在10～40 km，峰值在35 km.這些模型得到的滑動虧損區和震間蠕滑區(靠近海溝)與這次地震的破裂區明顯不符.

反演斷層閉鎖區方法的優點是由于正演使用了地震位錯理論模型，位移場有理論解，計算格林函數方便、簡單，但是不能很好地反映介質的不均勻性和復雜的斷層幾何形態.

地震應力模型是一種新的直接用來反演斷層應力變化的動力學模型(Xie and Cai,2018)，這個模型由斷層和彈性塊體組成，考慮到地震的破裂區厚度遠遠小于斷層的長度和寬度，斷層被視為兩個彈性塊體之間的接觸面.與斷層僅有剪切位錯運動的地震位錯模型不同，地震應力模型的斷層面上不但有剪應力變化，而且還有正應力變化.應力矢量在斷層接觸面上大小相等，方向相反.斷層面上的應力增加和減少分別代表應變能的積累和釋放；零剪應力線是斷層破裂區的邊界.用地震應力模型反演地表變形和斷層應力變化,不需要像位錯模型反演斷層滑動或閉鎖那樣需要預先人為地指定斷層的滑移區或閉鎖區，也不需要約束應力變化方向.這個方法已經成功地用來反演2011年日本東北9級大地震的同震斷層應力降(圖1a)及其震前7年地表變形和斷層應力的積累區(Xie et al.,2019).反演得到的3個震前應力積累區，最大的就是主震的破裂區，其余兩個分別是主震后30 min內發生在主震破裂區南面和北面的MW7.8和MW7.4大余震區(圖1c).

圖1 反演得到的斷層應力變化(Xie and Cai,2018;Xie et al.,2019)(a)和(b)分別是斷層同震剪應力和正應力變化.主震破裂區南部是主震后30 min內發生的M7.8大余震的破裂區,計算得到的主震和M7.8大余震的矩震級與USGS的一致;白色曲線為零剪應力變化等值線.(c)和(d)分別是震前7年(2004—2010年)內斷層剪應力和正應力變化.顯著的剪應力異常區與同震破裂區一致，紫色線為圖1c中同震破裂區邊界，主震異常區的南部和北部分別是M7.8和M7.4兩個大余震震前的剪應力積累區.Fig.1 Fault stress changes obtained by inversion (Xie and Cai,2018;Xie et al.,2019)(a)and (b)are coseismic shear and normal stress changes on the fault,respectively.The southern part of the main shock rupture zone is the rupture zone of the M7.8 aftershock that occurred within 30 minutes after the main shock.The calculated moment magnitudes of the main shock and the M7.8 aftershock are consistent with those obtained by the USGS.The white curve line is the contour of zero-shear stress change.(c)and (d)are the fault shear and normal stress changes within 7 years (2004—2010)before the main shock,respectively.The significant shear stress anomaly zone is consistent with the coseismic rupture zone.The purple line is the boundary of the coseismic rupture zone in Fig.1c.The shear stress anomaly areas on the south and north sides of the main shock are the shear stress accumulation areas before the two large aftershocks of M7.8 and M7.4,respectively.

本文是2011年日本9級大地震前7年斷層應力變化反演工作的繼續.為了進一步探討這次大地震應力積累過程和檢驗震前7年反演得到的孕震區的穩定性，我們在原來工作的基礎上(2004—2010年)，利用日本GNSS數據，反演震前14年到震前8年(1997—2003年)的斷層應力變化，同時結合反演得到的前7年斷層的應力變化，探討這次大地震斷層應力積累過程和孕震區隨時間演化.

1 反演方法

利用日本1997—2003年的GNSS數據,反演斷層應力變化所使用的方法與我們提出的方法(Xie and Cai,2018)相同，略述如下.

有限元方法用于正演問題的計算.整個模型的斷層面分為有限個子斷層，這些子斷層都是待反演的未知應力變化區.子斷層的數值格林函數解由下面定解問題確定：

Δσji,i=0,

(1)

Δσji=λΔεkkδij+2μΔεji,

(2)

Δεji=(Δuj,i+Δui,j)/2,

(3)

njΔσji|Γ1=Ti,

(4)

Δσji|Γ2=0,

(5)

Δuini|Γ3=0,

(6)

Δσji|Γ3=0,(i≠j).

(7)

上述方程(1)至(3)分別為平衡方程、本構方程和幾何方程，其中Δσji,Δεji和Δui分別是地震應力模型內部的應力、應變和位移變化，λ和μ是Lame常數，δij是Kronecker 符號，下角標相同的項表示求和；方程(4)至(7)為模型的邊界條件；Ti是法線方向為ni(i=1，2，3)的子斷層面上的單位應力分量，Г1,Г2和Г3分別代表模型的斷層面、上表面、側表面和底面邊界.

利用有限元方法可以得到上述邊值問題的每一個子斷層上單位應力分量載荷的位移場，由此可以合成數值格林函數G，進而構造如下的滿足GNSS觀測數據d(dx，dy，dz)的線性方程組：

(8)

β(Wh-Wf)T=0,

(9a)

λLT=0.

(9b)

(8)和(9)式中，T(Tx，Ty，Tz)是待求的由所有子斷層面上應力矢量構成的地震應力模型的斷層面上的總應力變化矢量,下標h和f分別代表斷層的上盤和下盤；(9)式中的W代表斷層面上的法向位移矢量，其中(9a)式和(9b)式分別代表斷層法向位移約束條件和解的平滑約束條件(Masterlack,2003)，其中β是約束參數，λ是平滑系數，L是拉普拉斯微分算子.

斷層面上的總應力變化矢量T，可以利用帶有約束的最小二乘方法，通過構造如下的目標函數(Xie and Cai,2018)：

Sobj(T)=‖dh-GhT‖+‖df+GfT‖

+β2‖(Wh-Wf)T‖+λ2‖LT‖

(10)

求得.

2 結果與討論

由于本文的研究是以前工作的繼續，所以研究區域、地震應力模型和有限元模型與文獻(Xie and Cai,2018)相同，如圖2所示，模型力學參數和反演中使用的平滑因子、斷層法向位移約束參數與我們同震和震前7年反演工作所使用的參數(Xie and Cai,2018;Xie et al.,2019)相同.使用的GNSS數據來自于http:∥www.gsi.go.jp/ENGLISH/page_e30233.html.

圖2 (a)2011年日本東北大地震研究區域;(b)建立在彈性理論和有限元方法上的地震應力模型.圖中紅線是斷層線;(c)斷層的上盤與下盤,Tr,Ts和Tn表示應力分量變化;(d)斷層面子斷層劃分，字母“B”表示較大的子斷層，遠離震源Fig.2 (a)Study area of the 2011 great Tohoku-Oki earthquake in Japan;(b)Earthquake stress model based on elastic theory and finite element method.The red line in the picture denotes the fault trace;(c)Hanging wall and foot wall of the fault,Tr,Ts and Tn represent changes of stress components,respectively;(d)Sub-fault division of the fault plane,the letter “B”indicates a larger sub-fault,far away from the source

為方便起見，我們把1997年到2003年反演得到的新結果稱為后期結果；把2004年到2010年反演得到的結果稱為前期結果.

圖3給出了14年間反演得到的后期和前期上盤斷層面剪應力的每年變化結果和每期的平均結果.從圖中可以看出，震前應力的變化逐年增加，但不是均勻的，其增加區和減小區與2004—2010年應力變化類似.每期平均剪應力增加最顯著的地區(圖3h和圖3p)與圖1a中的同震破裂區基本一致.平均剪應力在主震初始破裂點位置(五角星處)1997—2003年為0.03 MPa，2004—2010年為0.04 MPa.平均剪應力結果表明，剪應力積累速率不是均勻的，有變快的趨勢.2001年后,在主震的M7.8大余震區的剪應力明顯增加,但是在2008年不增反降，與其他年份不同.后來查明，這個地方在2008年發生了一次M6地震.由于此地震的變形場沒有從GNSS數據中剔除，所以上述剪應力降低區正是M6地震的應變能釋放導致的.這個結果無意中證明了反演方法的可靠性.2007年大余震區震前的剪應力為何也有些下降，是否是2008年M6地震的前震影響，目前原因尚不清楚，需要進一步研究.從圖3中剪應力的年平均值可以看出，剪應力明顯下降區大致出現在北緯36°到40°，經度在140°到142°，最淺深度在零剪應力變化線處，離海平面約29 km.剪應力降低意味著應變能釋放，這個剪應力降低區可以用前震、小的重復地震或者20到60 km深處的加速滑動解釋(Wang and Bilek,2014;Mavrommatis et al.,2014;Yokota and Koketsu,2015).

圖3 反演得到的斷層剪應力變化(a)—(g)1997—2003年每年剪應力變化；(h)7年剪應力變化的平均值；(i)—(o)2004—2010年每年剪應力變化；(p)7年剪應力變化的平均值.Fig.3 Shear stress changes obtained on fault by inversion(a)—(g)Yearly change of shear stress from 1997 to 2003;(h)Average of shear stress change in 7 years;(i)—(o)Yearly change of shear stress from 2004 to 2010;(p)Average of shear stress change in 7 years.

圖4給出了反演得到的后期和前期14年間上盤斷層面正應力的每年變化結果.從圖中可以看出，正應力的主要增加區也基本是穩定的，從1997年開始，在主震破裂區總的趨勢是逐漸增加，意味著在這些區域剪切強度增大.主震初始破裂點位置的平均正應力在1997—2003年為0.001 MPa，2004—2010年為0.006.這個結果表明，正應力積累速率也不是均勻的，有變快的趨勢，其增加值比剪應力小近一個量級，并且正應力增加區與剪應力的不一致.近南北的零正應力等值線深度(20 km)要淺于圖3中零剪應力等值線的深度(29 km).這意味著剪應力增加區包含一部分正應力降低區(零正應力變化等值線下方至零剪切應力變化等值線上方的區域).由于剪切強度與正應力成正比，可以推斷在剪應力增加區內存在剪切強度降低區，而在剪應力增加區內正應力增加的區域剪切強度增加.由此可見，零正應力變化等值線的位置就是斷層強度降低和增加的分界線.如果孕震區在整個孕震期間是基本穩定的，其有效摩擦系數基本相同，由于斷層初始破裂點下方區域的正應力隨著時間逐年降低，意味著那里的剪切強度也隨時間降低，這就解釋了為何這次地震在破裂開始的時候先向斷層下方傳播的原因(Ide et al.,2011).值得注意的是零正應力變化等值線不但在主震的初始破裂點(地震波定的震源位置)附近穿過，而且還在主震的兩個大余震的初始破裂點附近穿過，這一點目前還沒有看到有關文獻報道.這似乎不能用巧合來解釋.這可能是由于這些破裂點基本都處在剪應力增加最顯著的地方，在其上方由于正應力增加顯著，導致剪切強度明顯增加，不利于在那里開始破裂；而在其下方雖然正應力降低，但是幅度較小，因此剪切強度降低不明顯，雖然該區處在剪應力增加區，但是其增加不明顯，因此在該區開始破裂的條件不如在零應力變化等值線附近.

圖4 反演得到的斷層正應力變化(a)—(g)1997—2003年每年正應力變化；(h)7年正應力變化的平均值；(i)—(o)2004—2010年每年正應力變化；(p)7年正應力變化的平均值.Fig.4 Normal stress changes obtained on fault by inversion(a)—(g)Yearly change of normal stress from 1997 to 2003;(h)Average of normal stress change in 7 years;(i)—(o)Yearly change of normal stress from 2004 to 2010;(p)Average of normal stress change in 7 years.

圖5給出了使用最佳光滑因子和最佳約束系數都為0.05的地表觀測位移和反演得到的位移的匹配情況.從圖中可以看到，它們符合得很好.這表明反演得到的斷層應力變化是可靠的，能夠滿足觀測數據.

圖5 符合觀測位移的最佳地震應力模型(1997—2003)(a)—(g)水平位移矢量；(h)—(n)垂直位移矢量.紅色和黑色箭頭分別代表反演和觀測得到的位移.Fig.5 Displacements best fitted to the earthquake stress model of the Tohoku-Oki earthquake (1997—2003)(a)—(g)Horizontal displacement vectors;(h)—(n)Vertical displacement vectors.The red and black arrows represent the inverted and observed displacements,respectively.

圖6解釋了俯沖帶逆沖大地震的孕震體與地震空區、地表應變積累區和斷層面應力的關系.圖中紅色A區代表孕震體,應力或應變能增加區，白色B區代表應力減小區或應變能釋放區，黃色C區代表地震空區和由三角網得到的應變積累區，A和B間的粗黑線是孕震體的邊界，即零剪應力變化等值線，白色區域的黑箭頭表示B區的運動方向，藍色區的黑箭頭表示俯沖板塊的運動方向；在零剪應力上方的綠線是零正應力變化等值線，其上的圓圈是震源體破裂點在海底的投影(海底震中).圖中白色箭頭便是剪切位移(與板塊尺寸不成比例)，白箭頭聯線與過斷層的垂直白線之間的夾角是剪應變大小，其與該點剪應力成正比.

圖6 板塊俯沖孕育大地震的示意圖Fig.6 Schematic drawing of a large earthquake bred by plate subduction

3 結論

1997—2010年先后兩期共14年反演得到的每年斷層剪應力顯著增加區與斷層同震的破裂區結果基本一致，表明利用地震應力模型反演得到的剪應力顯著增加區可以用來預測俯沖帶大地震的發生位置.在觀測數據足夠多的地區，預測的可靠性更有保障.反演得到的14年間的主震的孕震體雖然是穩定的，但是孕震體的大小在整個孕震期間，或主震的復發間隔內，是否也是穩定的，本文結果不能回答這個問題，尚需進一步研究.大地震一般孕震時間很長，譬如這次大地震，復發周期長達千年(宍倉正展,2012).如果震源體積是隨時間變化的，并且觀測資料在時間上不足夠長，在空間上不足夠密，就不能保證所預測的震源體大小的可靠性，盡管如此，預測的地震大概位置還是有參考價值的.

地震產生的變形不但與斷層剪應力有關，而且還與正應力有關.對于傾滑或逆沖型地震，正應力的影響尤為重要.利用地震位錯模型得到的閉鎖區(滑動虧損區)明顯與主震破裂區不符的原因，可能是地震位錯模型中沒有考慮斷層面上法向應力變化對閉鎖區深度和大小的影響(Xie et al.,2019)，使得所預測的地震破裂區位置過深.

我們發現零剪應力和零正應力變化的等值線不重合，前者在斷層面上的深度大于后者，這意味著在剪應力增加的區域存在強度降低區.震源初始破裂點容易發生在零正應力變化等值線附近.本文結果表明，正應力變化對俯沖帶大地震的破裂有影響.

與地震空區方法、地表三角測量尋找應變(二維)異常區方法和反演斷層滑動虧損方法預測地震發生地點相比，斷層應力反演方法是最直接的動力學方法，抓住了預測大地震位置的關鍵因素或主要矛盾.隨著空間技術的日新月異，這種方法可望在地震動力學機制研究和地震位置預報中發揮重要的作用.