王仁先生在地震预报中的开拓性工作

石耀霖，胡才博

中国科学院大学地球与行星科学学院，计算地球动力学重点实验室，北京 100049

0 引言

中国是一个多地震的国家.在王仁先生1970年代初期转入地球动力学研究前后，1966年3月8日邢台M6.7地震造成了8064人遇难，1976年7月28日唐山7.8级地震造成24万多人遇难，地震灾害的预测和预防是一个不可回避的科学任务.王仁先生急国家和人民之所急，对于地震预报研究给予极大的重视，开展了一系列开拓性工作.他认为要充分利用现代力学成果，不限于定性的描述，更要从力学机制出发，通过深入研究使地震预报研究更具普遍性和预见性.他率先将有限元方法引进到中国的地球动力学领域，提出了板块运动学和动力学的联合反演方法以及地应力场的反演方法，在深入了解地震断层的破裂机制基础上，提出了研究地震迁移、地震危险区预测的思想与方法.

1994年，王仁与国际著名地震学家安艺敬一(K.Aki)在北京共同主持了由国际理论和应用力学联合会(IUTAM)主办的“地球动力学中的力学问题国际研讨会”.会后在国际刊物PureandAppliedGeophysics上出版了会议文集.在这次会议上，王仁回顾了世界百余年来地球动力学发展的历史以及中国学者的有关工作，有力地促进了中国的地球动力学研究，也促进了国外同行对中国在这个研究领域工作的交流和了解(Wang, 1995; Wang and Aki, 1996).这次会议上，王仁被安艺敬一誉为“中国的地球动力学之父”.

1 王仁先生的地震预报科学思想

什么是地震？王仁先生概括为：“在地应力作用下，地壳内某些脆弱的地带，当地壳应力超过这里的强度时，就会造成岩层的破裂，断层的突然错动，发生地震.”地震是一种力学现象.而这个力学过程会引起震源周围相当大区域内一系列物理化学响应：“在这个应力逐渐增加直到岩层破裂的整个过程中，在震源周围的广大地区内，地壳岩层的物理、化学性质不断地发生变化，引起各种前兆现象.总结这些经验，找出其规律，就能把地震预报工作做得更好.”(王仁，1977).

近年来地震学界关于地震预报和地震预测的定义、区别和联系有许多讨论，陈运泰(2009)在《地震预测：回顾与展望》一文有非常详细的论述.在王仁先生开展地震力学研究的年代，尚没有这样的仔细划分，因此本文沿用王仁先生地震预报的提法，主要是着重强调他在从传统的“基于地震活动性和其他多种前兆现象的经验预报”向“基于对地震孕育发生物理过程了解的物理预报”的转变方面开拓性的工作.

地震预报难题的解决，要抓基础、抓本质，搞清力学过程和机制.王仁先生认为：“地震是一个极端复杂的过程，要进一步把地震预报工作搞好，需要在广泛深入实践的基础上，对地震的本质进行探索，取得规律性的认识.由于地震表现出来的是一个力学过程，地震预报工作需要从力学方面进行探讨.”(王仁，1977).他的思路可以概括为图1来表达.

图1 对王仁先生开展地震预报研究的科学思路概括总结示意图

解决力学难题必须关注两个方面：应力场和介质.王仁先生认为：“从力学的观点分析，岩层断裂的矛盾双方，是促使岩层破裂的推动力，同岩层抵抗破裂的阻力的矛盾.前者就是地壳在各种外力作用下的应力分布，后者是岩层的强度.地震的孕育和发生就是这对矛盾对立统一的发展过程.”(王仁，1977).

王仁先生对地质力学中的一些力学问题，如与巨型纬向构造体系、经向构造体系、欧亚山字型构造等密切相关的全球性的构造应力场，与雁行排列的压性构造、多字型构造、山字型构造等有关的区域性构造应力场，褶皱变形问题的波长和流变性等，都有很深入的研究和论述(王仁，1976).

怎样了解应力场？地球动力学所研究的问题与普通工程中所遇到的力学问题虽然有相同的方面，都是研究固体在外力作用下的应力分布和变形特性，不过条件复杂多了，人们面对的地质体的变形和现今的应力是历史的产物.“一个办法是对地应力场进行实测……现场的应力测量，测得一些地点的现有应力状态，使我们对接近地表面的应力状态有了一些感性知识.可是离形成一个应力场还很远，而且深度最多只达到2公里左右，深部的应力还没法直接测量……地震的震源机制进行了各种分析，提出了震源处的主应力方向和地震前后的应力差.大地测量工作者在地表特别在地震活动区布置了测量网，逐年进行重复测量来观察地壳的变形，由此也可推算出地壳的受力情况.地质工作者则根据地质构造运动的遗迹，追究它们的力学成因，对这地区曾经受到过的力进行推断.特别是李四光同志运用力学观点把表面上错综复杂的地质现象，看成是统一的构造应力场的产物，归纳出各种类型的构造体系.这些来自实际的资料对于我们分析一个地区的应力分布都会有很大帮助，但离开掌握地壳内的应力场还有很大的差距.”(王仁，1977).

全面了解三维应力场的“另一个可能的办法是充分运用高速电子计算机的优越性，先参照实际情况选定一个地区，选择一些参数，应用有限单元计算方法进行计算，然后把计算结果和其它实测结果比较，对原选参数进行修改.经过不断的修改使计算和尽量多的实测结果相符.这样可认为应力场逐渐趋近于实际情况.” “分析地壳内的应力分布是一个反序的问题，……，解是不唯一的.”(王仁，1977).“求解驱动机制和变形过程需要处理两方面的反演：在空间上反演地球内部的过程，时间上反演历史的过程.另外，材料性质是高度非线性的，变形非常大，在几何上也是非线性的，因而是一个高度非线性的反演问题.”(王仁，1989).影响计算结果的因素包括边界条件、结构物性、初始条件，而对于地质体这些都往往是难以确定的，需要综合分析和多种途径去进行.

对于介质，尤其是发生地震的岩石圈地壳，首先要搞清地质构造背景.“要通过地震地质工作把与本地区地震密切相关的构造体系划分出来，根据地质条件选定这个地区的边界，和区域内一些主要的活动性构造，特别是活动性断层.我们将认为绝大多数地震是原有断层再次活动的结果.”(王仁，1977).“然后就是要确定这个地区内介质的力学性质，包括弹性模量、粘性和塑性的应力-应变关系，以及断层的抗剪和抗拉等强度条件.它们应该是在发震深度处的温度和围压条件下的性质.对于这些参数，目前知道得还较差，选择的范围就较广，要做出较准确的估计，还需要进行大量的实地考察和实验室工作.”(王仁，1977).“因为岩石多少具备渗透性，孔隙液体在大地构造过程中的作用是很重要的.孔隙压力增加时，有效应力减小，破坏所需的剪应力减小，……应力与变形就这样与液体流动耦合了起来.液体的热膨胀和增压进一步使应力和流场与温度耦合.构造过程中水力-热-力学的完全耦合是一个待研究的重要课题.在水库诱发地震分析中某些矿物的水化减弱效应也很重要.”(王仁，1989).

介质性质研究中另一个重要问题是，要研究岩石的破裂过程和规律.“需要从力学原理上探讨非稳定材料的应力-应变关系.这也是对力学提出的新课题.”王仁先生十分重视基于野外现场测量和实验室内实验的岩石力学性质的研究工作，强调将实验室内的试验结果用于实际情况时需要进行空间尺度的外推和时间尺度的外推，这两种外推需要考虑地球介质变形的物理机制，对单轴压缩试验全过程、橄榄石微观变形机制等有很好的总结(王仁，1983).

前兆异常的分析.“地震前兆异常的种类很多，有些异常可以和岩石的脆性破裂联系起来.”(王仁，1977).“其它手段如地面变形、地表倾斜、重力、地壳电阻率、地震波波速等也都有类似现象，它们和岩石破裂过程及地应力的变化应该都是有密切关系的.”(王仁，1977).“我们觉得，如果把这个地区的应力场变化规律摸得清楚些，就能对异常形态做出更好的解释.”.虽然当时数字化的记录才逐渐推广，还没有积累现今海量的大数据，但王仁先生已经注意到，“需要创造新的途径，引用新技术，以及用电子计算机对数据进行及时的综合处理.”

2 王仁先生在地震预报中的创新探索

王仁先生把他的思想付诸于实践，取得了一系列具有鲜明创新思想、紧扣时代前沿的学术成果.

2.1 应力场反演

在应力场反演方面，利用震源机制解、大地测量和应力测量数据作为约束条件，研发了联合反演板块运动学和动力学参数的方法.首先把速度和应力边界条件视为一组待求的参数，对于给定的岩石力学参数进行计算，并根据实际观测的约束条件，采用叠加原理和最小二乘方法把这组参数反演出来.通过反复调节岩石力学参数和相应的反演，直至得到一个与实际观测结果拟合较满意的结果.他通过反演得到了北美板块的速度场和应力场以及相应的板块边界的运动速率和其底部的拖曳力(王仁和梁海华，1985).他也对我国华北应力场用计算和模拟多种方法进行了反演(王仁等，1982a)，这个应力场是预测华北地震危险区的基础之一.平面分析可以作为初次近似，但地震本是一个空间问题，要进一步的近似和对有倾角的断层错动问题，则需进行三维应力场分析，复杂程度会增加很多(王仁，1977).王仁先生认为本构关系和破坏准则可能是大地构造、地震数值模拟的最重要问题，本构关系和破坏准则随环境条件变化且含有时间因素(王仁，1989)，这些观点至今还发挥着巨大作用.

地球内部的应力可以分为动态和静态两种.前者指那些不能积累起来的短周期变化的应力，它们不能推动构造运动，只能起到触发作用，例如地震波和地球自由振荡、固体潮、极移以及短期地球自转速率变化引起的应力波动，后者指与重力、长期的构造运动、地球内部的热和孔隙流体压力等有关的应力.地应力对地震活动性、岩体工程的稳定性以及矿产资源的成矿环境具有重要意义.在构造地质学中，将偏离静岩压力(一点各个方向的压力都相等，类似于静水压力)的那部分地应力称为构造应力，它可以分为历史上构造变形遗留下来的古应力和现今构造运动产生的应力.由于受到技术、财力和其他种种因素的制约，人们迄今能得到的现场地应力量数据还很少，虽然少数测点深度可达数千米，但多数都不超过 1000 m，而且精度较低.尽管可以通过地震震源机制解、构造地质学分析和大地测量手段得到地应力场的大致方向，但很难精确确定其大小.因此，如何获得地下的应力状态，目前仍是地球动力学中难以解决的问题.

2.2 岩石力学实验

在介质研究的岩石力学实验方面，王仁先生的研究团队把光弹、云纹、激光全息干涉和扫描电镜技术以及有限元方法引入到岩石力学实验的综合研究中，利用这些实验手段监测岩石试件在不同加载时刻的变形，记录破裂演变过程.他们利用扫描电镜实时研究了具有割缝的大理岩的破坏过程，观察到试件宏观破坏是细观裂纹逐步集中和联结所导致的(赵永红等，1992，1993).他们还将损伤力学引入到岩石力学实验中，研究了岩石试件表面裂纹的发展变化及其与外载和预制缺陷之间的关系(赵永红等，1994).以实验观测和分维统计分析结果为依据，建立了岩石分维损伤本构模型，这些可以为模拟各种岩石材料的实际问题，以及大尺度地球动力学问题的分析计算提供基础.王仁先生研究团队利用实验室试验和有限元模拟研究了含有裂缝的大理岩的X型剪切破裂过程，对新破裂的扩展方向的力学机制进行了分析(Wang et al., 1987)；利用热弹塑增量理论和有限元法研究了热状态对地震发生的影响，认为断层外部地震的发生被认为是岩石应变软化的结果，而断层内部地震的发生被认为是断层热软化的结果(蔡永恩等，1992).

2.3 单一强震和地震序列的数值模拟

王仁先生在地震预报方面的主要贡献是，在对应力场和岩体介质研究的基础上，开展地震序列的数学物理模拟，进行地震危险区的预测.一次大地震后，还会不会发生大的余震；一个地震活动区域，历史和未来序列的大地震会怎样发生迁移，这是人们最为关切的问题，这对于防震减灾具有重要意义.1976年唐山地震后，人们迫切关注京津唐地区近期地震形势.王仁先生急国家人民之所急，提出了模拟地震迁移与危险区预测的思想和方法，带领研究团队，基于构造地质背景，利用中国历史地震资料的约束，开展弹塑性有限元方法计算模拟，先后研究了1966年邢台大地震后到唐山地震期间京津唐地区大地震的迁移(王仁等，1980)和华北700年间14个7 级以上历史大地震的迁移过程(王仁等，1982b)，预测了唐山地震后的危险地区和北京的安全度.从理论上说，如果能知道一个地区的地应力分布和岩石强度，就可以估计那里地震的危险性.但是在实践中这二者都很难确切知道.王仁先生创新研究的方法的核心是，根据构造地质学和地震地质学确定研究区的地质构造格局和边界形状，利用地震波速和岩石力学实验确定区城岩石力学参数，结合现场应力、形变测量和震源机制解确定边界外力.通过调整边界外力和岩石力学参数，调整外力和介质的力学性质各参数，使地区内第一个历史地震震中处于临破裂状态，然后令其单元破裂释放应力，继续计算变化后构造作用下演变的应力场，并再次对参数进行修改，使第二个历史地震处临近破裂.如此反复，复现历史上该地区全部地震序列.修改参数使计算结果接近这个序列将是一个十分繁复的工作.不过也只有这样，才能够有把握地用这个应力场的分析来推测今后的地震趋势.初步尝试获得了令人鼓舞的结果.预测了唐山地震后的地震危险区和北京的地震安全度.研究成果在1979年巴黎召开的地震预报会议和1981年国际地震与地球内部物理学协会(International Association of Seismology and Physics of the Earth′s Interior, IASPEI)上进行了报告，引起了与会者的极大关注，开创了研究地震迁移与危险区预测的新方法.之后，他领导的团队进一步将不连续界面的拉格朗日不连续变形分析(Lagrange Discontinous Deformation Analysis, LDDA)方法引入到地震物理过程的研究(Cai et al., 2000)，计算了唐山地震引起的准静态和动态应力变化.

王仁先生提出了利用大地测量和地震学资料进行板块运动学和动力学的联合反演的方法.为了深入了解地震断层的破裂机制，他将塑性力学和损伤力学引人到岩石破裂机理的研究中.王仁先生重视力学基础理论的研究，出版了固体力学专著或论文集(王仁等，1979；王仁和黄杰藩，1981；王仁，1982).王仁和殷有泉系统介绍了工程岩石类介质的弹塑性本构关系，对弹性、塑性、黏性、黏弹塑性本构与岩石材料的变形、破坏、蠕变、松弛等的关系进行了深入的解析，对弹塑性问题的全量理论、增量理论和内在时间理论的关键公式、相互关系和程序研发进行了介绍(王仁和殷有泉，1981)，这些是数值地震预报的奠基性工作.王仁先生十分重视地质材料的非弹性变形性质和破裂特点，关注地下爆破、水库水位变化引起的地下渗流、软弱岩层的蠕变，这是地震数值模拟和预报不可避免的关键科学问题(王仁，1986).王仁先生从地球动力学的名称起源(由著名弹性力学家勒夫在1911年第一次提出)谈起，评述了李四光等人将力学应用于地质构造分析的早期工作，对我国学者在20世纪八九十年代在地球动力学方面取得的进展进行了认真的梳理和总结(王仁，1997).

关于震前和震后效应的区分.一次地震后，地应力重新分布，它会引起周围地区出现各种“异常”，有时要很长时间才会逐渐平息.由于地震是一次接一次的，这些异常，究竟是另一次大地震的前兆呢还是前次地震的震后效应呢，在这两种情况下的应力场上可能会有区别(王仁，1977).1976 年唐山大地震后，王仁等为了回答这个问题，提出了利用弹塑性理论和库仑破裂准则研究震后地震危险区和地震迁移的方法(王仁等，1980，1982a，b).这个方法的基本思想是：地震的发生不但与岩石的强度有关，还与地震断层震前的应力状态有关.一个地方发生了地震，震源处积累的弹性应变能被释放，使地震断层内部的应力降低，从而导致应力的重新分布.为了研究大地震后的地震危险性问题，王仁提出了地震序列的数学模拟方法和地震安全度(等于断层内的剪应力与断层摩擦强度之差与断层的摩擦强度的比值)的概念.使所要研究的地震序列中某次地震的断层震前应力最接近破裂的临界值，然后通过降低该地断层的摩擦系数，模拟地震的发生(从静摩擦到动摩擦)，并且使计算的断层位错和所释放的应变能与地震学得到的该次地震的结果一致.在计算下一次地震时，上一次地震的断层错距被保留下来，断层摩擦系数恢复到原来静摩擦的状态.按照这种方式依次对地震序列中的每一个地震进行模拟.如果每次所得到的地震位错和所释放的应变能都与实际地震相符合，那么最后得到的这个应力场就被认为是一个比较接近实际的应力场了.这样就可以由这个应力场和地震安全度预测将来的地震危险区(王仁等，1980，1982b).

王仁先生研究团队利用新的LDDA方法模拟了1976年唐山大地震的破裂、错动和应力释放的动力学全过程，考虑了静摩擦和动摩擦的转化，考虑了接触非线性的影响，模拟结果与实际观测一致(蔡永恩等，1999)；利用三维黏弹性有限元模型研究了1976年大地震同震及震后变形，认为华北板块下方软流层黏度为7.1×1018Pa·s，上地幔黏度为2.1×1019Pa·s(孙荀英等，1994)；以新丰江水库区构造的三维数值模拟为例研究了深层构造新活动性的影响，认为深层新构造活动对浅层构造活动有重要作用，不可忽视(丁原章等，1992)；利用震前断层蠕动模型，研究了1976年唐山大地震之前的地温异常，认为加速蠕动产生的唐山相对玉田的表层温度累计变化为2.5 °C，与观测结果一致(蔡永恩等，1987).在丛书《21世纪100个科学难题》中王仁先生有专门章节提出第37个科学难题“地球构造运动驱动机制的反演”(王仁，1998)，为地震数值模拟和预测提出了重要发展方向.

3 王仁先生地震预报思想得到了继承和发展

从王仁先生1970年代投身地球动力学和地震预报研究以来，已经过去了半个世纪.科学技术手段有了巨大的发展，尤其近20年来，许多新方法新技术已经得到了完善和大规模应用，大大增强了我们对地震物理过程研究的能力，积累了丰富的观测大数据.高性能计算技术有了飞跃的发展，三维复杂介质的非线性耦合物理和力学问题的有限元计算已经成为可能；GNSS/GPS和InSAR等空间大地测量手段已经得到完善和发展，并初步积累了有用的资料来认识地壳应变场的孕震、同震和震后变化；地震层析成像分辨能力极大增强，噪声成像把“干扰”变成了信息的来源；人工地震深部探测技术的发展及在世界各地大规模的应用增强了我们对岩石圈结构的“透视”能力；中国、美国、日本等国家在应力变化观测能力上取得了突出的进展；与王仁先生提出的地震造成应力变化影响后续地震活动性的地震安全度类似的库仑应力概念已经被普遍接受采用；有关能源开采的技术应用造成的人工触发地震，已经积累很多的研究案例，并受到科学界、产业界和普通民众的广泛重视，同时提供了对孕震机理研究的新途径；断层本构关系有了更进一步的实验探索和经验总结；新技术装备的地震地质考察揭示了活断层的分布和活动特征.这些进展都使得王仁先生的许多想法可以推进和实现.

3.1 数值地震预报在国际学术界的发展和应用

王仁先生关于通过有限元计算推断地震活动演变趋势的思想(Wang and Wu, 1983; 王仁，1994)，被证明符合科学发展的主流.国际上2001年由美国国家航空航天局(NASA)资助了QuakeSim项目，用于对地震断层系统进行建模.致力于通过各种边界元、有限元和分析应用程序对地震过程进行建模，构建一个可互操作的工具系统，用以更好地理解活动构造和地震过程(Pierce et al., 2008).QuakeSim的目标是显著提高地震预报质量，从而减轻这种自然灾害带来的危险，在2012年已经推出了QuakeSim 2.0.

国际知名学者Keilis-Borok于1985年在意大利成立了结构与非线性动力学小组，在1990年，他成立了前苏联(现俄罗斯)科学院的国际地震预测理论和数学地球物理研究所(Gabrielov et al., 2014).1997年，亚太经合组织地震科学合作项目ACES获批，至今已运行了20余年，获得了大量科技成果，研究领域覆盖了地震孕育发生过程涉及的微观、细观和宏观地球物理问题(张永仙等，2020).意大利在2009年4月6日拉奎拉地震后发展了一套可操作的地震预测(operational earthquake forecast,OEF)技术(Jordan et al., 2011).美国加利福尼亚州发展了“统一的加利福尼亚州地震破裂预测(Uniform California Earthquake Rupture Forecast, UCERF)”系统，在2007年提出的UCERF2基础上,2014年又提出了3rd Uniform California Earthquake Rupture Forecast(UCERF3)模型(Field et al, 2014).

3.2 数值地震预报在国内学术界的发展和应用

刘启元和吴建春(2003)提出了需要打破长期徘徊在以地震前兆异常监测为基础的经验性预测局面，把注意力尽快转向研究以动力学为基础的数值预报.中国地震局编制了2007—2020年的《国家地震科学技术发展纲要》，在“国家地震减灾科学计划”章节，明确要求安排“地震数值预测试验研究”的专项.朱守彪等认为通过基于严格数学、力学原理的有限元方法计算可以模拟断层由闭锁到解锁产生地震的全过程，认为给出未来强震发生的时间、空间、强度三要素的数值预报是可行的(朱守彪等，2008；朱守彪，2012).石耀霖等对我国地震数值预报路线图提出了更具体的建议，提出了长期时间无关地震概率预测的具体实施方案，发展了相关计算方法和程序，给出了青藏高原东北缘地震时空迁移和概率预测的研究案例，并提出了中长期时间相关地震概率预报的技术路线图，对地震短临概率预报做了一些探索和研究(石耀霖等，2018).

其中，时间无关的长期预报得到进一步的发展和改进.孙云强等把复杂断层系统下地震序列的研究，从王仁先生当年的700年延伸到数万年，把二维的模型推广到三维的模型，把弹塑性模型发展为黏弹塑性有限元模型，模拟了区域断层系统的地震循环和地震序列的时空演化，获得了万年时间尺度的人工合成地震目录(孙云强和罗纲，2018；孙云强等，2019，2020).在模型满足区域地球动力学背景的基础上，根据模拟的人工合成地震目录分析了青藏高原东北缘各断层上不同位置、不同震级的地震复发特征(图2).模型预测可以得到数千年考古地震资料的验证支持，模型模拟得到的人工地震目录b值为0.96，也与几十年来仪器记录的地震目录的b值1.01接近.该黏弹塑性有限元模型的特点是，回避初始条件的不确定性，通过对青藏高原东北缘地震活动进行数万年的计算模拟，弱化了较短时间(例如几百年)内初始条件的明显影响，而侧重于讨论数万年的模拟中，研究地区整体和各个断层不同段落表现出的平均活动特征和发震概率(图3).

图2 青藏高原东北缘三维黏弹塑性有限元计算模型的结构、断层系统和边界条件示意图(据孙云强和罗纲，2018)

图3 青藏高原东北缘断层系统根据黏弹塑性有限元计算得到的长期平均的地震发生概率(据孙云强等，2020)

时间无关的长期预报模型的弱点是，仅仅是一个长时间平均结果，无法考虑最近数十年、几百年大震序列对地应力场演变和地震活动性发展的影响.例如汶川地震的复发时间大约为数千年，2008年汶川地震释放了多年积累的应变能后，会有很长时间不会再在原地发生同等规模的大地震，但是这种情况在时间无关模型中无法考虑.而考虑历史地震序列，推断未来地震活动，正是王仁先生1980年代华北地区研究工作的科学思路的精髓.时间相关的中长期预报也得到了发展.胡才博等利用三维黏弹性模型研究了长期构造加载、历史大地震和中下地壳和上地幔的黏弹性松弛对2008年汶川大地震触发的影响(胡才博，2009；Hu et al., 2012).董培育等基于岩石库仑摩尔破裂准则，利用青藏高原及邻区百年历史范围内的强震信息，来反演估算该区域的初始应力场.然后，考虑区域构造应力加载及强震造成的应力扰动共同作用，重现了青藏高原一百多年历史强震的发展过程，并对今后地震活动进行预测(董培育等，2020).与王仁先生1980年代初的工作比较，不仅由于计算技术的发展，计算模型从80年代的700多个单元和接近400个节点，发展到31万多单元和17万多节点规模上的扩大；断层的处理也从各向同性物性改变为横向各向同性近似，能够对逆掩断层在平面简化时的应力变化进行更好的模拟；更重要的是，由于有了80年代还不存在的GPS位移资料积累，董培育等的模型，可以利用边界上的GPS资料插值给出位移速度边界条件，并把计算结果的速度值与研究区域内的GPS观测值比较、把计算的应力张量与震源机制观测和实际应力测量得到的主压应力方向进行比较，得到更加吻合实际的模型.在模拟中发现，简单二维模型在一些地区不能与实际观测吻合，这些地区可能是由于下地壳大规模流动对上地壳拖曳作用，使得二维简化已经不能再成立.然而考虑到计算时调整参数的巨大工作量，目前还不适宜直接进行三维模拟，因此利用了2.5维的计算方法，将下地壳的拖曳力近似表达为水平二维模型中水平方向的等效体力，使得复杂的三维问题得到了良好的二维近似表达，计算的每年运动矢量与GPS实际观测值有很好的吻合(图4).

关于初始应力场的反演估算，分为两种情况.一是历史上发生过大地震的地方，根据图4速度矢量可以计算出应变率，再根据胡克定律可以计算出应力增长率.地震时应力状态一定是达到了破裂强度，如果知道岩体和断裂强度，又知道应力每年增长速率，就可以推算出模型初始时刻的应力.对于没有发生过大地震的地方，由于未地震时断裂系统处于自组织亚临界状态，应力不会离破裂强度太远，设定合理下限；而从模型初始年份到现在一直没有发生大地震，则从应力增长率可以估计其初始应力状态上限.初始应力处于上限和下限中的某个状态，我们无法唯一确定.因此，采用王仁先生曾经建议过，但限于计算能力当时没有实施的 Monte Carlo随机方法(王仁，1977)，进行大量独立的随机试验计算，生成成千上万个各自不同的区域初始应力场模型，每个模型都能复现历史强震有序发生过程，但未来应力场演化过程和后续地震序列发生会不尽相同.最后，将这成千上万个模型在未来时间段内的危险性预测结果集成，得到不同断层、不同段落的地震危险性的统计结果.我们对青藏高原1904—2014年26个M≥6.5破坏性地震序列进行复现模拟.该序列中的第一个地震为1904年8月30日发生在鲜水河上的道孚MW6.8地震，最后一个地震为2014年2月12日发生在阿尔金断裂带上的于田地震.对1000个初始条件不同的模型进行了计算，它们都能够重复26个历史地震序列，但是其余部位应力初始条件是在一定应力范围内随机产生而各个彼此不同.从它们计算得到的青藏高原未来地震危险性如图5所示.其中的概率是按1000个随机模型中有多大百分比产生了破裂来统计计算.从2014年到现在研究区内发生了4个M6.5以上地震，其中2017年九寨沟M7.0地震和2020年于田M6.5地震发生在预测高概率区临近部位.

图4 考虑下地壳流动对上地壳拖曳作用后，青藏高原研究区域内计算的和GPS实际观测的运动速度矢量的对比(据董培育等，2020)

图5 2015年到2054年青藏高原研究预测区内计算的发震高概率区(据董培育等，2020)

在上述模型的同震效应的计算上，王仁先生当年在二维模拟中采用达到破裂准则后，减小震源处摩擦系数μ来计算破裂引起的位移和应力变化(王仁等，1980，1982a，1982b).这种方法也被后续研究者采用和发展(杨树新等，2012).董培育和石耀霖(2013)更进一步指出采用横向各向同性杀伤单元与三维位错模型的解答更吻合，并且可以用于走滑断层和逆掩断层不同类型断层影响的二维简化之下的同震应力场计算.

计算能力的提高和断层速度相关、状态相关破裂本构关系研究的进展，使得人们不但能模拟地震断层造成的静态应力变化，而且尝试通过自适应调整时间步长，进一步在一个数值模拟案例中既能模拟孕震的缓慢过程，也能模拟断层错动的迅速动态过程(朱守彪，2012；Zhu and Zhang, 2013；袁杰等，2021).

固体潮汐应力对地震的触发方面.日月引力不但可以引起海水的涨落，而且可以引起地球的变形.前者称为海潮，后者称为固体潮.导致固体潮的力称为引潮力，它是日月对地球的引力和地球自转引起的惯性离心力的合力.引潮力能否触发地震一直是地震预报的一个热门话题.国内外很多人対这个问题进行了研究，目前还没有一个确定的结论.地震能否被触发，不但取决于断层面上固体潮引起的应力大小，还取决于地震前断层面上的初始应力与断层的摩擦滑动强度相差多少.王仁等利用15层快速地球模型得到的固体潮应力场的理论解，提出了研究固体潮触发地震的理论方法(王仁和丁中一，1979).这个方法是首先利用固体潮应力场的理论解，计算出在发震时刻断层面上的正应力和沿断层错动方向的剪应力.假设断层已经处于破裂的临界状态，然后根据岩石的库仑破裂准则判断断层面上的应力对地震的触发作用.如果断层面上的应力落在库仑破裂线的外面，就表示能够触发，反之，就不能触发.通过对中国以及国外较大地震的研究得知，潮汐应力对浅源走滑型地震的触发效应明显，可以达到50%以上，对倾滑型地震则不明显.而新的研究表明，不仅月球在地球上造成的固体潮对地震有触发作用，而且地球在月球上造成的固体潮应力变化对深源月震有触发作用(张贝等，2016).不仅潮汐力可以触发地震，而且气候变化也可以触发(Carlson et al., 2020)；不仅自然因素可以触发，而且人类活动也可以稳定或触发地震，例如高坝水库蓄水(Cheng et al., 2016)、大面积过量开采地下水(Kundu et al., 2015)、地热和页岩气开发的注水压裂(Lei et al., 2017)，都可能会诱发或者触发中小强度的天然地震，也会造成人类生命财产的严重损失，必须引起足够重视.

王仁先生关注地震前兆的力学机制.吴忠良等讨论了地震前兆检验的地球动力学问题，指出需要考虑地震前观测到的各种前兆资料，“反演”地震的孕育过程，在这一地球动力学过程的背景下，重新审视观测到的“前兆”异常信息与地震之间的“对应”(吴忠良和蒋长胜，2006；吴忠良等，2009).

3.3 大数据和人工智能在地震预报的发展和应用

王仁先生在世时，大数据及人工智能还没有得到今天这样的发展，但王仁先生一贯重视探索新的途径，引用新的技术，用电子计算机对数据进行综合处理(王仁，1977).现在，我们也尝试着把机器学习方法引入地震预报研究(石耀霖等，2021；李林芳等，2021)，通过长短期记忆(Long Short-Term Memory，LSTM)神经网络处理川滇地区1970—2004年的地震目录，利用16个地震预报因子和滑动的时空窗口，对川滇部分地区开展了基于神经网络的中期(1年)地震预报研究.发现用1970年到2004年地震目录作为训练集，把2005年到2019年为测试集进行回溯性预报检验时，实际震级落在预测震级±0.5内的准确率为70.2%，虚报率为18.7%，漏报率为11.1%，可以回溯性预测2008年汶川MS8.0地震及其他6级以上强震(图6).典型的预测值偏差分布见图7.通过扩大研究区域范围、改变大震级地震在均方差计算中的权重等测试，方法依然表现稳定，结果具有可重复性.这显示地震数值预报不仅可以探索中长期地震预报，而且有可能在更短时间的地震预报方面开展实用化研究.

图6 对川滇研究区域2005年后6级以上地震回溯性预测检验(据石耀霖等，2021)

图7 典型的震级预测值与实际值偏差的分布柱状图

4 展望

我国从1970年代初设立了国家地震局，统一领导和部署地震相关工作和地震科学研究，每年年底进行下一年度地震活动的预报，这在世界上是绝无仅有的.中国和世界迄今在地震预报方面只取得了有限度的很少量的成功，离攻克实用化的地震预报目标还相距遥远.但是与1966年邢台地震前相比，我国在台站建设方面，从几十个地震台站发展到了仅隶属中国地震台网中心的固定地震台站1100多个，并具备了数千台流动台站部署的能力；从几乎没有前兆台站发展到各种类型前兆台站3000多个，特别是从仅有传统大地测量手段到具备了GNSS观测的260个基准站和3000多个流动站及InSAR等空间大地测量手段，从没有地应力监测能力到具备数十个优质台站对地应力长期变化、潮汐变化和地震动态应力波动进行记录的能力.在基础地震地质和地球物理勘察方面，从一般的地质调查到对出露和隐伏的地震活断层的针对性调查，获得了全国200多条活动构造带上的2000余个几何学和运动学参数.从开展少量地震测深观测到地震测深剖面总长度接近60000 km.

正演计算方面，考虑地壳介质的三维孔隙黏弹塑性和断层的速度相关、状态相关本构关系，自适应调整时间步长的准静态和动态模型，可以全面地模拟地震的孕育、地震的蠕动、慢地震、低频地震或一般地震失稳破裂错动发展的过程，研究地震波的辐射传播和造成的动态应力影响，定量计算地震的同震弹性应力变化和震后的黏弹性应力调整.反演计算方面，机器自行调整发生过地震的地方的初始应力参数，随机形成其余部位合理范围内的初始应力，形成能够模拟重现历史地震序列及大地震破裂的单个模型，进而对未来地震趋势进行计算和预测.集成分析方面，要能够随机生成大量初始应力状态不同，但都能够重现历史大地震序列的模型，从成千上万甚至更多模型的计算中，统计分析未来研究区域各个部位的地震危险性，给出大地震发生的概率，提供基于物理模型的长期和中期地震概率预报.

我国还要开展地震孕育模型的计算结果与实际前兆和地震活动性的对比研究，发现前兆的时空演变特征和地震活动规律，在物理思想指导下开展大数据驱动的中期和短临地震预报研究，最终实现长中短临渐进式的地震数值预报.王仁先生当年前瞻性的科学思想开拓的道路已经越走越宽.我们一定要进一步推进地震数值预报的研究，不断取得新的研究成果，作为对王仁先生最好的纪念.

致谢感谢两位匿名审稿人的中肯意见.