石耀霖，胡才博

中國科學院大學地球與行星科學學院，計算地球動力學重點實驗室，北京 100049

0 引言

中國是一個多地震的國家.在王仁先生1970年代初期轉入地球動力學研究前后，1966年3月8日邢臺M6.7地震造成了8064人遇難，1976年7月28日唐山7.8級地震造成24萬多人遇難，地震災害的預測和預防是一個不可回避的科學任務.王仁先生急國家和人民之所急，對于地震預報研究給予極大的重視，開展了一系列開拓性工作.他認為要充分利用現(xiàn)代力學成果，不限于定性的描述，更要從力學機制出發(fā)，通過深入研究使地震預報研究更具普遍性和預見性.他率先將有限元方法引進到中國的地球動力學領域，提出了板塊運動學和動力學的聯(lián)合反演方法以及地應力場的反演方法，在深入了解地震斷層的破裂機制基礎上，提出了研究地震遷移、地震危險區(qū)預測的思想與方法.

1994年，王仁與國際著名地震學家安藝敬一(K.Aki)在北京共同主持了由國際理論和應用力學聯(lián)合會(IUTAM)主辦的“地球動力學中的力學問題國際研討會”.會后在國際刊物PureandAppliedGeophysics上出版了會議文集.在這次會議上，王仁回顧了世界百余年來地球動力學發(fā)展的歷史以及中國學者的有關工作，有力地促進了中國的地球動力學研究，也促進了國外同行對中國在這個研究領域工作的交流和了解(Wang, 1995; Wang and Aki, 1996).這次會議上，王仁被安藝敬一譽為“中國的地球動力學之父”.

1 王仁先生的地震預報科學思想

什么是地震？王仁先生概括為：“在地應力作用下，地殼內某些脆弱的地帶，當?shù)貧Τ^這里的強度時，就會造成巖層的破裂，斷層的突然錯動，發(fā)生地震.”地震是一種力學現(xiàn)象.而這個力學過程會引起震源周圍相當大區(qū)域內一系列物理化學響應：“在這個應力逐漸增加直到巖層破裂的整個過程中，在震源周圍的廣大地區(qū)內，地殼巖層的物理、化學性質不斷地發(fā)生變化，引起各種前兆現(xiàn)象.總結這些經驗，找出其規(guī)律，就能把地震預報工作做得更好.”(王仁，1977).

近年來地震學界關于地震預報和地震預測的定義、區(qū)別和聯(lián)系有許多討論，陳運泰(2009)在《地震預測：回顧與展望》一文有非常詳細的論述.在王仁先生開展地震力學研究的年代，尚沒有這樣的仔細劃分，因此本文沿用王仁先生地震預報的提法，主要是著重強調他在從傳統(tǒng)的“基于地震活動性和其他多種前兆現(xiàn)象的經驗預報”向“基于對地震孕育發(fā)生物理過程了解的物理預報”的轉變方面開拓性的工作.

地震預報難題的解決，要抓基礎、抓本質，搞清力學過程和機制.王仁先生認為：“地震是一個極端復雜的過程，要進一步把地震預報工作搞好，需要在廣泛深入實踐的基礎上，對地震的本質進行探索，取得規(guī)律性的認識.由于地震表現(xiàn)出來的是一個力學過程，地震預報工作需要從力學方面進行探討.”(王仁，1977).他的思路可以概括為圖1來表達.

圖1 對王仁先生開展地震預報研究的科學思路概括總結示意圖

解決力學難題必須關注兩個方面：應力場和介質.王仁先生認為：“從力學的觀點分析，巖層斷裂的矛盾雙方，是促使巖層破裂的推動力，同巖層抵抗破裂的阻力的矛盾.前者就是地殼在各種外力作用下的應力分布，后者是巖層的強度.地震的孕育和發(fā)生就是這對矛盾對立統(tǒng)一的發(fā)展過程.”(王仁，1977).

王仁先生對地質力學中的一些力學問題，如與巨型緯向構造體系、經向構造體系、歐亞山字型構造等密切相關的全球性的構造應力場，與雁行排列的壓性構造、多字型構造、山字型構造等有關的區(qū)域性構造應力場，褶皺變形問題的波長和流變性等，都有很深入的研究和論述(王仁，1976).

怎樣了解應力場？地球動力學所研究的問題與普通工程中所遇到的力學問題雖然有相同的方面，都是研究固體在外力作用下的應力分布和變形特性，不過條件復雜多了，人們面對的地質體的變形和現(xiàn)今的應力是歷史的產物.“一個辦法是對地應力場進行實測……現(xiàn)場的應力測量，測得一些地點的現(xiàn)有應力狀態(tài)，使我們對接近地表面的應力狀態(tài)有了一些感性知識.可是離形成一個應力場還很遠，而且深度最多只達到2公里左右，深部的應力還沒法直接測量……地震的震源機制進行了各種分析，提出了震源處的主應力方向和地震前后的應力差.大地測量工作者在地表特別在地震活動區(qū)布置了測量網(wǎng)，逐年進行重復測量來觀察地殼的變形，由此也可推算出地殼的受力情況.地質工作者則根據(jù)地質構造運動的遺跡，追究它們的力學成因，對這地區(qū)曾經受到過的力進行推斷.特別是李四光同志運用力學觀點把表面上錯綜復雜的地質現(xiàn)象，看成是統(tǒng)一的構造應力場的產物，歸納出各種類型的構造體系.這些來自實際的資料對于我們分析一個地區(qū)的應力分布都會有很大幫助，但離開掌握地殼內的應力場還有很大的差距.”(王仁，1977).

全面了解三維應力場的“另一個可能的辦法是充分運用高速電子計算機的優(yōu)越性，先參照實際情況選定一個地區(qū)，選擇一些參數(shù)，應用有限單元計算方法進行計算，然后把計算結果和其它實測結果比較，對原選參數(shù)進行修改.經過不斷的修改使計算和盡量多的實測結果相符.這樣可認為應力場逐漸趨近于實際情況.” “分析地殼內的應力分布是一個反序的問題，……，解是不唯一的.”(王仁，1977).“求解驅動機制和變形過程需要處理兩方面的反演：在空間上反演地球內部的過程，時間上反演歷史的過程.另外，材料性質是高度非線性的，變形非常大，在幾何上也是非線性的，因而是一個高度非線性的反演問題.”(王仁，1989).影響計算結果的因素包括邊界條件、結構物性、初始條件，而對于地質體這些都往往是難以確定的，需要綜合分析和多種途徑去進行.

對于介質，尤其是發(fā)生地震的巖石圈地殼，首先要搞清地質構造背景.“要通過地震地質工作把與本地區(qū)地震密切相關的構造體系劃分出來，根據(jù)地質條件選定這個地區(qū)的邊界，和區(qū)域內一些主要的活動性構造，特別是活動性斷層.我們將認為絕大多數(shù)地震是原有斷層再次活動的結果.”(王仁，1977).“然后就是要確定這個地區(qū)內介質的力學性質，包括彈性模量、粘性和塑性的應力-應變關系，以及斷層的抗剪和抗拉等強度條件.它們應該是在發(fā)震深度處的溫度和圍壓條件下的性質.對于這些參數(shù)，目前知道得還較差，選擇的范圍就較廣，要做出較準確的估計，還需要進行大量的實地考察和實驗室工作.”(王仁，1977).“因為巖石多少具備滲透性，孔隙液體在大地構造過程中的作用是很重要的.孔隙壓力增加時，有效應力減小，破壞所需的剪應力減小，……應力與變形就這樣與液體流動耦合了起來.液體的熱膨脹和增壓進一步使應力和流場與溫度耦合.構造過程中水力-熱-力學的完全耦合是一個待研究的重要課題.在水庫誘發(fā)地震分析中某些礦物的水化減弱效應也很重要.”(王仁，1989).

介質性質研究中另一個重要問題是，要研究巖石的破裂過程和規(guī)律.“需要從力學原理上探討非穩(wěn)定材料的應力-應變關系.這也是對力學提出的新課題.”王仁先生十分重視基于野外現(xiàn)場測量和實驗室內實驗的巖石力學性質的研究工作，強調將實驗室內的試驗結果用于實際情況時需要進行空間尺度的外推和時間尺度的外推，這兩種外推需要考慮地球介質變形的物理機制，對單軸壓縮試驗全過程、橄欖石微觀變形機制等有很好的總結(王仁，1983).

前兆異常的分析.“地震前兆異常的種類很多，有些異常可以和巖石的脆性破裂聯(lián)系起來.”(王仁，1977).“其它手段如地面變形、地表傾斜、重力、地殼電阻率、地震波波速等也都有類似現(xiàn)象，它們和巖石破裂過程及地應力的變化應該都是有密切關系的.”(王仁，1977).“我們覺得，如果把這個地區(qū)的應力場變化規(guī)律摸得清楚些，就能對異常形態(tài)做出更好的解釋.”.雖然當時數(shù)字化的記錄才逐漸推廣，還沒有積累現(xiàn)今海量的大數(shù)據(jù)，但王仁先生已經注意到，“需要創(chuàng)造新的途徑，引用新技術，以及用電子計算機對數(shù)據(jù)進行及時的綜合處理.”

2 王仁先生在地震預報中的創(chuàng)新探索

王仁先生把他的思想付諸于實踐，取得了一系列具有鮮明創(chuàng)新思想、緊扣時代前沿的學術成果.

2.1 應力場反演

在應力場反演方面，利用震源機制解、大地測量和應力測量數(shù)據(jù)作為約束條件，研發(fā)了聯(lián)合反演板塊運動學和動力學參數(shù)的方法.首先把速度和應力邊界條件視為一組待求的參數(shù)，對于給定的巖石力學參數(shù)進行計算，并根據(jù)實際觀測的約束條件，采用疊加原理和最小二乘方法把這組參數(shù)反演出來.通過反復調節(jié)巖石力學參數(shù)和相應的反演，直至得到一個與實際觀測結果擬合較滿意的結果.他通過反演得到了北美板塊的速度場和應力場以及相應的板塊邊界的運動速率和其底部的拖曳力(王仁和梁海華，1985).他也對我國華北應力場用計算和模擬多種方法進行了反演(王仁等，1982a)，這個應力場是預測華北地震危險區(qū)的基礎之一.平面分析可以作為初次近似，但地震本是一個空間問題，要進一步的近似和對有傾角的斷層錯動問題，則需進行三維應力場分析，復雜程度會增加很多(王仁，1977).王仁先生認為本構關系和破壞準則可能是大地構造、地震數(shù)值模擬的最重要問題，本構關系和破壞準則隨環(huán)境條件變化且含有時間因素(王仁，1989)，這些觀點至今還發(fā)揮著巨大作用.

地球內部的應力可以分為動態(tài)和靜態(tài)兩種.前者指那些不能積累起來的短周期變化的應力，它們不能推動構造運動，只能起到觸發(fā)作用，例如地震波和地球自由振蕩、固體潮、極移以及短期地球自轉速率變化引起的應力波動，后者指與重力、長期的構造運動、地球內部的熱和孔隙流體壓力等有關的應力.地應力對地震活動性、巖體工程的穩(wěn)定性以及礦產資源的成礦環(huán)境具有重要意義.在構造地質學中，將偏離靜巖壓力(一點各個方向的壓力都相等，類似于靜水壓力)的那部分地應力稱為構造應力，它可以分為歷史上構造變形遺留下來的古應力和現(xiàn)今構造運動產生的應力.由于受到技術、財力和其他種種因素的制約，人們迄今能得到的現(xiàn)場地應力量數(shù)據(jù)還很少，雖然少數(shù)測點深度可達數(shù)千米，但多數(shù)都不超過 1000 m，而且精度較低.盡管可以通過地震震源機制解、構造地質學分析和大地測量手段得到地應力場的大致方向，但很難精確確定其大小.因此，如何獲得地下的應力狀態(tài)，目前仍是地球動力學中難以解決的問題.

2.2 巖石力學實驗

在介質研究的巖石力學實驗方面，王仁先生的研究團隊把光彈、云紋、激光全息干涉和掃描電鏡技術以及有限元方法引入到巖石力學實驗的綜合研究中，利用這些實驗手段監(jiān)測巖石試件在不同加載時刻的變形，記錄破裂演變過程.他們利用掃描電鏡實時研究了具有割縫的大理巖的破壞過程，觀察到試件宏觀破壞是細觀裂紋逐步集中和聯(lián)結所導致的(趙永紅等，1992，1993).他們還將損傷力學引入到巖石力學實驗中，研究了巖石試件表面裂紋的發(fā)展變化及其與外載和預制缺陷之間的關系(趙永紅等，1994).以實驗觀測和分維統(tǒng)計分析結果為依據(jù)，建立了巖石分維損傷本構模型，這些可以為模擬各種巖石材料的實際問題，以及大尺度地球動力學問題的分析計算提供基礎.王仁先生研究團隊利用實驗室試驗和有限元模擬研究了含有裂縫的大理巖的X型剪切破裂過程，對新破裂的擴展方向的力學機制進行了分析(Wang et al., 1987)；利用熱彈塑增量理論和有限元法研究了熱狀態(tài)對地震發(fā)生的影響，認為斷層外部地震的發(fā)生被認為是巖石應變軟化的結果，而斷層內部地震的發(fā)生被認為是斷層熱軟化的結果(蔡永恩等，1992).

2.3 單一強震和地震序列的數(shù)值模擬

王仁先生在地震預報方面的主要貢獻是，在對應力場和巖體介質研究的基礎上，開展地震序列的數(shù)學物理模擬，進行地震危險區(qū)的預測.一次大地震后，還會不會發(fā)生大的余震；一個地震活動區(qū)域，歷史和未來序列的大地震會怎樣發(fā)生遷移，這是人們最為關切的問題，這對于防震減災具有重要意義.1976年唐山地震后，人們迫切關注京津唐地區(qū)近期地震形勢.王仁先生急國家人民之所急，提出了模擬地震遷移與危險區(qū)預測的思想和方法，帶領研究團隊，基于構造地質背景，利用中國歷史地震資料的約束，開展彈塑性有限元方法計算模擬，先后研究了1966年邢臺大地震后到唐山地震期間京津唐地區(qū)大地震的遷移(王仁等，1980)和華北700年間14個7 級以上歷史大地震的遷移過程(王仁等，1982b)，預測了唐山地震后的危險地區(qū)和北京的安全度.從理論上說，如果能知道一個地區(qū)的地應力分布和巖石強度，就可以估計那里地震的危險性.但是在實踐中這二者都很難確切知道.王仁先生創(chuàng)新研究的方法的核心是，根據(jù)構造地質學和地震地質學確定研究區(qū)的地質構造格局和邊界形狀，利用地震波速和巖石力學實驗確定區(qū)城巖石力學參數(shù)，結合現(xiàn)場應力、形變測量和震源機制解確定邊界外力.通過調整邊界外力和巖石力學參數(shù)，調整外力和介質的力學性質各參數(shù)，使地區(qū)內第一個歷史地震震中處于臨破裂狀態(tài)，然后令其單元破裂釋放應力，繼續(xù)計算變化后構造作用下演變的應力場，并再次對參數(shù)進行修改，使第二個歷史地震處臨近破裂.如此反復，復現(xiàn)歷史上該地區(qū)全部地震序列.修改參數(shù)使計算結果接近這個序列將是一個十分繁復的工作.不過也只有這樣，才能夠有把握地用這個應力場的分析來推測今后的地震趨勢.初步嘗試獲得了令人鼓舞的結果.預測了唐山地震后的地震危險區(qū)和北京的地震安全度.研究成果在1979年巴黎召開的地震預報會議和1981年國際地震與地球內部物理學協(xié)會(International Association of Seismology and Physics of the Earth′s Interior, IASPEI)上進行了報告，引起了與會者的極大關注，開創(chuàng)了研究地震遷移與危險區(qū)預測的新方法.之后，他領導的團隊進一步將不連續(xù)界面的拉格朗日不連續(xù)變形分析(Lagrange Discontinous Deformation Analysis, LDDA)方法引入到地震物理過程的研究(Cai et al., 2000)，計算了唐山地震引起的準靜態(tài)和動態(tài)應力變化.

王仁先生提出了利用大地測量和地震學資料進行板塊運動學和動力學的聯(lián)合反演的方法.為了深入了解地震斷層的破裂機制，他將塑性力學和損傷力學引人到巖石破裂機理的研究中.王仁先生重視力學基礎理論的研究，出版了固體力學專著或論文集(王仁等，1979；王仁和黃杰藩，1981；王仁，1982).王仁和殷有泉系統(tǒng)介紹了工程巖石類介質的彈塑性本構關系，對彈性、塑性、黏性、黏彈塑性本構與巖石材料的變形、破壞、蠕變、松弛等的關系進行了深入的解析，對彈塑性問題的全量理論、增量理論和內在時間理論的關鍵公式、相互關系和程序研發(fā)進行了介紹(王仁和殷有泉，1981)，這些是數(shù)值地震預報的奠基性工作.王仁先生十分重視地質材料的非彈性變形性質和破裂特點，關注地下爆破、水庫水位變化引起的地下滲流、軟弱巖層的蠕變，這是地震數(shù)值模擬和預報不可避免的關鍵科學問題(王仁，1986).王仁先生從地球動力學的名稱起源(由著名彈性力學家勒夫在1911年第一次提出)談起，評述了李四光等人將力學應用于地質構造分析的早期工作，對我國學者在20世紀八九十年代在地球動力學方面取得的進展進行了認真的梳理和總結(王仁，1997).

關于震前和震后效應的區(qū)分.一次地震后，地應力重新分布，它會引起周圍地區(qū)出現(xiàn)各種“異常”，有時要很長時間才會逐漸平息.由于地震是一次接一次的，這些異常，究竟是另一次大地震的前兆呢還是前次地震的震后效應呢，在這兩種情況下的應力場上可能會有區(qū)別(王仁，1977).1976 年唐山大地震后，王仁等為了回答這個問題，提出了利用彈塑性理論和庫侖破裂準則研究震后地震危險區(qū)和地震遷移的方法(王仁等，1980，1982a，b).這個方法的基本思想是：地震的發(fā)生不但與巖石的強度有關，還與地震斷層震前的應力狀態(tài)有關.一個地方發(fā)生了地震，震源處積累的彈性應變能被釋放，使地震斷層內部的應力降低，從而導致應力的重新分布.為了研究大地震后的地震危險性問題，王仁提出了地震序列的數(shù)學模擬方法和地震安全度(等于斷層內的剪應力與斷層摩擦強度之差與斷層的摩擦強度的比值)的概念.使所要研究的地震序列中某次地震的斷層震前應力最接近破裂的臨界值，然后通過降低該地斷層的摩擦系數(shù)，模擬地震的發(fā)生(從靜摩擦到動摩擦)，并且使計算的斷層位錯和所釋放的應變能與地震學得到的該次地震的結果一致.在計算下一次地震時，上一次地震的斷層錯距被保留下來，斷層摩擦系數(shù)恢復到原來靜摩擦的狀態(tài).按照這種方式依次對地震序列中的每一個地震進行模擬.如果每次所得到的地震位錯和所釋放的應變能都與實際地震相符合，那么最后得到的這個應力場就被認為是一個比較接近實際的應力場了.這樣就可以由這個應力場和地震安全度預測將來的地震危險區(qū)(王仁等，1980，1982b).

王仁先生研究團隊利用新的LDDA方法模擬了1976年唐山大地震的破裂、錯動和應力釋放的動力學全過程，考慮了靜摩擦和動摩擦的轉化，考慮了接觸非線性的影響，模擬結果與實際觀測一致(蔡永恩等，1999)；利用三維黏彈性有限元模型研究了1976年大地震同震及震后變形，認為華北板塊下方軟流層黏度為7.1×1018Pa·s，上地幔黏度為2.1×1019Pa·s(孫荀英等，1994)；以新豐江水庫區(qū)構造的三維數(shù)值模擬為例研究了深層構造新活動性的影響，認為深層新構造活動對淺層構造活動有重要作用，不可忽視(丁原章等，1992)；利用震前斷層蠕動模型，研究了1976年唐山大地震之前的地溫異常，認為加速蠕動產生的唐山相對玉田的表層溫度累計變化為2.5 °C，與觀測結果一致(蔡永恩等，1987).在叢書《21世紀100個科學難題》中王仁先生有專門章節(jié)提出第37個科學難題“地球構造運動驅動機制的反演”(王仁，1998)，為地震數(shù)值模擬和預測提出了重要發(fā)展方向.

3 王仁先生地震預報思想得到了繼承和發(fā)展

從王仁先生1970年代投身地球動力學和地震預報研究以來，已經過去了半個世紀.科學技術手段有了巨大的發(fā)展，尤其近20年來，許多新方法新技術已經得到了完善和大規(guī)模應用，大大增強了我們對地震物理過程研究的能力，積累了豐富的觀測大數(shù)據(jù).高性能計算技術有了飛躍的發(fā)展，三維復雜介質的非線性耦合物理和力學問題的有限元計算已經成為可能；GNSS/GPS和InSAR等空間大地測量手段已經得到完善和發(fā)展，并初步積累了有用的資料來認識地殼應變場的孕震、同震和震后變化；地震層析成像分辨能力極大增強，噪聲成像把“干擾”變成了信息的來源；人工地震深部探測技術的發(fā)展及在世界各地大規(guī)模的應用增強了我們對巖石圈結構的“透視”能力；中國、美國、日本等國家在應力變化觀測能力上取得了突出的進展；與王仁先生提出的地震造成應力變化影響后續(xù)地震活動性的地震安全度類似的庫侖應力概念已經被普遍接受采用；有關能源開采的技術應用造成的人工觸發(fā)地震，已經積累很多的研究案例，并受到科學界、產業(yè)界和普通民眾的廣泛重視，同時提供了對孕震機理研究的新途徑；斷層本構關系有了更進一步的實驗探索和經驗總結；新技術裝備的地震地質考察揭示了活斷層的分布和活動特征.這些進展都使得王仁先生的許多想法可以推進和實現(xiàn).

3.1 數(shù)值地震預報在國際學術界的發(fā)展和應用

王仁先生關于通過有限元計算推斷地震活動演變趨勢的思想(Wang and Wu, 1983; 王仁，1994)，被證明符合科學發(fā)展的主流.國際上2001年由美國國家航空航天局(NASA)資助了QuakeSim項目，用于對地震斷層系統(tǒng)進行建模.致力于通過各種邊界元、有限元和分析應用程序對地震過程進行建模，構建一個可互操作的工具系統(tǒng)，用以更好地理解活動構造和地震過程(Pierce et al., 2008).QuakeSim的目標是顯著提高地震預報質量，從而減輕這種自然災害帶來的危險，在2012年已經推出了QuakeSim 2.0.

國際知名學者Keilis-Borok于1985年在意大利成立了結構與非線性動力學小組，在1990年，他成立了前蘇聯(lián)(現(xiàn)俄羅斯)科學院的國際地震預測理論和數(shù)學地球物理研究所(Gabrielov et al., 2014).1997年，亞太經合組織地震科學合作項目ACES獲批，至今已運行了20余年，獲得了大量科技成果，研究領域覆蓋了地震孕育發(fā)生過程涉及的微觀、細觀和宏觀地球物理問題(張永仙等，2020).意大利在2009年4月6日拉奎拉地震后發(fā)展了一套可操作的地震預測(operational earthquake forecast,OEF)技術(Jordan et al., 2011).美國加利福尼亞州發(fā)展了“統(tǒng)一的加利福尼亞州地震破裂預測(Uniform California Earthquake Rupture Forecast, UCERF)”系統(tǒng)，在2007年提出的UCERF2基礎上,2014年又提出了3rd Uniform California Earthquake Rupture Forecast(UCERF3)模型(Field et al, 2014).

3.2 數(shù)值地震預報在國內學術界的發(fā)展和應用

劉啟元和吳建春(2003)提出了需要打破長期徘徊在以地震前兆異常監(jiān)測為基礎的經驗性預測局面，把注意力盡快轉向研究以動力學為基礎的數(shù)值預報.中國地震局編制了2007—2020年的《國家地震科學技術發(fā)展綱要》，在“國家地震減災科學計劃”章節(jié)，明確要求安排“地震數(shù)值預測試驗研究”的專項.朱守彪等認為通過基于嚴格數(shù)學、力學原理的有限元方法計算可以模擬斷層由閉鎖到解鎖產生地震的全過程，認為給出未來強震發(fā)生的時間、空間、強度三要素的數(shù)值預報是可行的(朱守彪等，2008；朱守彪，2012).石耀霖等對我國地震數(shù)值預報路線圖提出了更具體的建議，提出了長期時間無關地震概率預測的具體實施方案，發(fā)展了相關計算方法和程序，給出了青藏高原東北緣地震時空遷移和概率預測的研究案例，并提出了中長期時間相關地震概率預報的技術路線圖，對地震短臨概率預報做了一些探索和研究(石耀霖等，2018).

其中，時間無關的長期預報得到進一步的發(fā)展和改進.孫云強等把復雜斷層系統(tǒng)下地震序列的研究，從王仁先生當年的700年延伸到數(shù)萬年，把二維的模型推廣到三維的模型，把彈塑性模型發(fā)展為黏彈塑性有限元模型，模擬了區(qū)域斷層系統(tǒng)的地震循環(huán)和地震序列的時空演化，獲得了萬年時間尺度的人工合成地震目錄(孫云強和羅綱，2018；孫云強等，2019，2020).在模型滿足區(qū)域地球動力學背景的基礎上，根據(jù)模擬的人工合成地震目錄分析了青藏高原東北緣各斷層上不同位置、不同震級的地震復發(fā)特征(圖2).模型預測可以得到數(shù)千年考古地震資料的驗證支持，模型模擬得到的人工地震目錄b值為0.96，也與幾十年來儀器記錄的地震目錄的b值1.01接近.該黏彈塑性有限元模型的特點是，回避初始條件的不確定性，通過對青藏高原東北緣地震活動進行數(shù)萬年的計算模擬，弱化了較短時間(例如幾百年)內初始條件的明顯影響，而側重于討論數(shù)萬年的模擬中，研究地區(qū)整體和各個斷層不同段落表現(xiàn)出的平均活動特征和發(fā)震概率(圖3).

圖2 青藏高原東北緣三維黏彈塑性有限元計算模型的結構、斷層系統(tǒng)和邊界條件示意圖(據(jù)孫云強和羅綱，2018)

圖3 青藏高原東北緣斷層系統(tǒng)根據(jù)黏彈塑性有限元計算得到的長期平均的地震發(fā)生概率(據(jù)孫云強等，2020)

時間無關的長期預報模型的弱點是，僅僅是一個長時間平均結果，無法考慮最近數(shù)十年、幾百年大震序列對地應力場演變和地震活動性發(fā)展的影響.例如汶川地震的復發(fā)時間大約為數(shù)千年，2008年汶川地震釋放了多年積累的應變能后，會有很長時間不會再在原地發(fā)生同等規(guī)模的大地震，但是這種情況在時間無關模型中無法考慮.而考慮歷史地震序列，推斷未來地震活動，正是王仁先生1980年代華北地區(qū)研究工作的科學思路的精髓.時間相關的中長期預報也得到了發(fā)展.胡才博等利用三維黏彈性模型研究了長期構造加載、歷史大地震和中下地殼和上地幔的黏彈性松弛對2008年汶川大地震觸發(fā)的影響(胡才博，2009；Hu et al., 2012).董培育等基于巖石庫侖摩爾破裂準則，利用青藏高原及鄰區(qū)百年歷史范圍內的強震信息，來反演估算該區(qū)域的初始應力場.然后，考慮區(qū)域構造應力加載及強震造成的應力擾動共同作用，重現(xiàn)了青藏高原一百多年歷史強震的發(fā)展過程，并對今后地震活動進行預測(董培育等，2020).與王仁先生1980年代初的工作比較，不僅由于計算技術的發(fā)展，計算模型從80年代的700多個單元和接近400個節(jié)點，發(fā)展到31萬多單元和17萬多節(jié)點規(guī)模上的擴大；斷層的處理也從各向同性物性改變?yōu)闄M向各向同性近似，能夠對逆掩斷層在平面簡化時的應力變化進行更好的模擬；更重要的是，由于有了80年代還不存在的GPS位移資料積累，董培育等的模型，可以利用邊界上的GPS資料插值給出位移速度邊界條件，并把計算結果的速度值與研究區(qū)域內的GPS觀測值比較、把計算的應力張量與震源機制觀測和實際應力測量得到的主壓應力方向進行比較，得到更加吻合實際的模型.在模擬中發(fā)現(xiàn)，簡單二維模型在一些地區(qū)不能與實際觀測吻合，這些地區(qū)可能是由于下地殼大規(guī)模流動對上地殼拖曳作用，使得二維簡化已經不能再成立.然而考慮到計算時調整參數(shù)的巨大工作量，目前還不適宜直接進行三維模擬，因此利用了2.5維的計算方法，將下地殼的拖曳力近似表達為水平二維模型中水平方向的等效體力，使得復雜的三維問題得到了良好的二維近似表達，計算的每年運動矢量與GPS實際觀測值有很好的吻合(圖4).

關于初始應力場的反演估算，分為兩種情況.一是歷史上發(fā)生過大地震的地方，根據(jù)圖4速度矢量可以計算出應變率，再根據(jù)胡克定律可以計算出應力增長率.地震時應力狀態(tài)一定是達到了破裂強度，如果知道巖體和斷裂強度，又知道應力每年增長速率，就可以推算出模型初始時刻的應力.對于沒有發(fā)生過大地震的地方，由于未地震時斷裂系統(tǒng)處于自組織亞臨界狀態(tài)，應力不會離破裂強度太遠，設定合理下限；而從模型初始年份到現(xiàn)在一直沒有發(fā)生大地震，則從應力增長率可以估計其初始應力狀態(tài)上限.初始應力處于上限和下限中的某個狀態(tài)，我們無法唯一確定.因此，采用王仁先生曾經建議過，但限于計算能力當時沒有實施的 Monte Carlo隨機方法(王仁，1977)，進行大量獨立的隨機試驗計算，生成成千上萬個各自不同的區(qū)域初始應力場模型，每個模型都能復現(xiàn)歷史強震有序發(fā)生過程，但未來應力場演化過程和后續(xù)地震序列發(fā)生會不盡相同.最后，將這成千上萬個模型在未來時間段內的危險性預測結果集成，得到不同斷層、不同段落的地震危險性的統(tǒng)計結果.我們對青藏高原1904—2014年26個M≥6.5破壞性地震序列進行復現(xiàn)模擬.該序列中的第一個地震為1904年8月30日發(fā)生在鮮水河上的道孚MW6.8地震，最后一個地震為2014年2月12日發(fā)生在阿爾金斷裂帶上的于田地震.對1000個初始條件不同的模型進行了計算，它們都能夠重復26個歷史地震序列，但是其余部位應力初始條件是在一定應力范圍內隨機產生而各個彼此不同.從它們計算得到的青藏高原未來地震危險性如圖5所示.其中的概率是按1000個隨機模型中有多大百分比產生了破裂來統(tǒng)計計算.從2014年到現(xiàn)在研究區(qū)內發(fā)生了4個M6.5以上地震，其中2017年九寨溝M7.0地震和2020年于田M6.5地震發(fā)生在預測高概率區(qū)臨近部位.

圖4 考慮下地殼流動對上地殼拖曳作用后，青藏高原研究區(qū)域內計算的和GPS實際觀測的運動速度矢量的對比(據(jù)董培育等，2020)

圖5 2015年到2054年青藏高原研究預測區(qū)內計算的發(fā)震高概率區(qū)(據(jù)董培育等，2020)

在上述模型的同震效應的計算上，王仁先生當年在二維模擬中采用達到破裂準則后，減小震源處摩擦系數(shù)μ來計算破裂引起的位移和應力變化(王仁等，1980，1982a，1982b).這種方法也被后續(xù)研究者采用和發(fā)展(楊樹新等，2012).董培育和石耀霖(2013)更進一步指出采用橫向各向同性殺傷單元與三維位錯模型的解答更吻合，并且可以用于走滑斷層和逆掩斷層不同類型斷層影響的二維簡化之下的同震應力場計算.

計算能力的提高和斷層速度相關、狀態(tài)相關破裂本構關系研究的進展，使得人們不但能模擬地震斷層造成的靜態(tài)應力變化，而且嘗試通過自適應調整時間步長，進一步在一個數(shù)值模擬案例中既能模擬孕震的緩慢過程，也能模擬斷層錯動的迅速動態(tài)過程(朱守彪，2012；Zhu and Zhang, 2013；袁杰等，2021).

固體潮汐應力對地震的觸發(fā)方面.日月引力不但可以引起海水的漲落，而且可以引起地球的變形.前者稱為海潮，后者稱為固體潮.導致固體潮的力稱為引潮力，它是日月對地球的引力和地球自轉引起的慣性離心力的合力.引潮力能否觸發(fā)地震一直是地震預報的一個熱門話題.國內外很多人対這個問題進行了研究，目前還沒有一個確定的結論.地震能否被觸發(fā)，不但取決于斷層面上固體潮引起的應力大小，還取決于地震前斷層面上的初始應力與斷層的摩擦滑動強度相差多少.王仁等利用15層快速地球模型得到的固體潮應力場的理論解，提出了研究固體潮觸發(fā)地震的理論方法(王仁和丁中一，1979).這個方法是首先利用固體潮應力場的理論解，計算出在發(fā)震時刻斷層面上的正應力和沿斷層錯動方向的剪應力.假設斷層已經處于破裂的臨界狀態(tài)，然后根據(jù)巖石的庫侖破裂準則判斷斷層面上的應力對地震的觸發(fā)作用.如果斷層面上的應力落在庫侖破裂線的外面，就表示能夠觸發(fā)，反之，就不能觸發(fā).通過對中國以及國外較大地震的研究得知，潮汐應力對淺源走滑型地震的觸發(fā)效應明顯，可以達到50%以上，對傾滑型地震則不明顯.而新的研究表明，不僅月球在地球上造成的固體潮對地震有觸發(fā)作用，而且地球在月球上造成的固體潮應力變化對深源月震有觸發(fā)作用(張貝等，2016).不僅潮汐力可以觸發(fā)地震，而且氣候變化也可以觸發(fā)(Carlson et al., 2020)；不僅自然因素可以觸發(fā)，而且人類活動也可以穩(wěn)定或觸發(fā)地震，例如高壩水庫蓄水(Cheng et al., 2016)、大面積過量開采地下水(Kundu et al., 2015)、地熱和頁巖氣開發(fā)的注水壓裂(Lei et al., 2017)，都可能會誘發(fā)或者觸發(fā)中小強度的天然地震，也會造成人類生命財產的嚴重損失，必須引起足夠重視.

王仁先生關注地震前兆的力學機制.吳忠良等討論了地震前兆檢驗的地球動力學問題，指出需要考慮地震前觀測到的各種前兆資料，“反演”地震的孕育過程，在這一地球動力學過程的背景下，重新審視觀測到的“前兆”異常信息與地震之間的“對應”(吳忠良和蔣長勝，2006；吳忠良等，2009).

3.3 大數(shù)據(jù)和人工智能在地震預報的發(fā)展和應用

王仁先生在世時，大數(shù)據(jù)及人工智能還沒有得到今天這樣的發(fā)展，但王仁先生一貫重視探索新的途徑，引用新的技術，用電子計算機對數(shù)據(jù)進行綜合處理(王仁，1977).現(xiàn)在，我們也嘗試著把機器學習方法引入地震預報研究(石耀霖等，2021；李林芳等，2021)，通過長短期記憶(Long Short-Term Memory，LSTM)神經網(wǎng)絡處理川滇地區(qū)1970—2004年的地震目錄，利用16個地震預報因子和滑動的時空窗口，對川滇部分地區(qū)開展了基于神經網(wǎng)絡的中期(1年)地震預報研究.發(fā)現(xiàn)用1970年到2004年地震目錄作為訓練集，把2005年到2019年為測試集進行回溯性預報檢驗時，實際震級落在預測震級±0.5內的準確率為70.2%，虛報率為18.7%，漏報率為11.1%，可以回溯性預測2008年汶川MS8.0地震及其他6級以上強震(圖6).典型的預測值偏差分布見圖7.通過擴大研究區(qū)域范圍、改變大震級地震在均方差計算中的權重等測試，方法依然表現(xiàn)穩(wěn)定，結果具有可重復性.這顯示地震數(shù)值預報不僅可以探索中長期地震預報，而且有可能在更短時間的地震預報方面開展實用化研究.

圖6 對川滇研究區(qū)域2005年后6級以上地震回溯性預測檢驗(據(jù)石耀霖等，2021)

圖7 典型的震級預測值與實際值偏差的分布柱狀圖

4 展望

我國從1970年代初設立了國家地震局，統(tǒng)一領導和部署地震相關工作和地震科學研究，每年年底進行下一年度地震活動的預報，這在世界上是絕無僅有的.中國和世界迄今在地震預報方面只取得了有限度的很少量的成功，離攻克實用化的地震預報目標還相距遙遠.但是與1966年邢臺地震前相比，我國在臺站建設方面，從幾十個地震臺站發(fā)展到了僅隸屬中國地震臺網(wǎng)中心的固定地震臺站1100多個，并具備了數(shù)千臺流動臺站部署的能力；從幾乎沒有前兆臺站發(fā)展到各種類型前兆臺站3000多個，特別是從僅有傳統(tǒng)大地測量手段到具備了GNSS觀測的260個基準站和3000多個流動站及InSAR等空間大地測量手段，從沒有地應力監(jiān)測能力到具備數(shù)十個優(yōu)質臺站對地應力長期變化、潮汐變化和地震動態(tài)應力波動進行記錄的能力.在基礎地震地質和地球物理勘察方面，從一般的地質調查到對出露和隱伏的地震活斷層的針對性調查，獲得了全國200多條活動構造帶上的2000余個幾何學和運動學參數(shù).從開展少量地震測深觀測到地震測深剖面總長度接近60000 km.

正演計算方面，考慮地殼介質的三維孔隙黏彈塑性和斷層的速度相關、狀態(tài)相關本構關系，自適應調整時間步長的準靜態(tài)和動態(tài)模型，可以全面地模擬地震的孕育、地震的蠕動、慢地震、低頻地震或一般地震失穩(wěn)破裂錯動發(fā)展的過程，研究地震波的輻射傳播和造成的動態(tài)應力影響，定量計算地震的同震彈性應力變化和震后的黏彈性應力調整.反演計算方面，機器自行調整發(fā)生過地震的地方的初始應力參數(shù)，隨機形成其余部位合理范圍內的初始應力，形成能夠模擬重現(xiàn)歷史地震序列及大地震破裂的單個模型，進而對未來地震趨勢進行計算和預測.集成分析方面，要能夠隨機生成大量初始應力狀態(tài)不同，但都能夠重現(xiàn)歷史大地震序列的模型，從成千上萬甚至更多模型的計算中，統(tǒng)計分析未來研究區(qū)域各個部位的地震危險性，給出大地震發(fā)生的概率，提供基于物理模型的長期和中期地震概率預報.

我國還要開展地震孕育模型的計算結果與實際前兆和地震活動性的對比研究，發(fā)現(xiàn)前兆的時空演變特征和地震活動規(guī)律，在物理思想指導下開展大數(shù)據(jù)驅動的中期和短臨地震預報研究，最終實現(xiàn)長中短臨漸進式的地震數(shù)值預報.王仁先生當年前瞻性的科學思想開拓的道路已經越走越寬.我們一定要進一步推進地震數(shù)值預報的研究，不斷取得新的研究成果，作為對王仁先生最好的紀念.

致謝感謝兩位匿名審稿人的中肯意見.