博文
地震数值预报
精选
|||
地震数值预报
石耀霖 张贝 张斯奇 张怀
(中国科学院大学中国科学院计算动力学重点实验室北京 100049)
Numerical earthquake prediction
SHI Yao-Lin† ZHANG Bei ZHANG Si-Qi ZHANG Huai
(University of the Chinese Academy of Sciences, Laboratory of Computational Geodynamics
of Chinese Academy of Sciences, Beijing 100049, China)
摘要 地震能否预报、地震预报遇到的种种困难如何克服？在思考这一问题时，不妨借助于气象预报的经验。气象预报从经验预报发展到数值预报，虽然用了约半个世纪的时间，但获得了成功。中国地震学界，对于地震预报的发展，也必须从经验预报发展为物理预报，对此几乎没有异议；但是物理预报的基本特点就是要基于定量的物理规律，进行数值预报，却鲜被提及。文章讨论了开展地震数值预报的五个主要环节的现状，认为地震数值预报现在应该提到规划的日程上来，设计我国地震数值预报的路线图，按照它的需要部署地震台站观测工作，开展地震数值预报的理论探讨。
关键词地震预报，地球动力学，数值方法，非线性动力学，天气预报
Abstract Can earthquakes be predicted? How should people overcome the difficulties encountered in the study of earthquake prediction? This issue can take inspiration from the experiences of weather forecast. Although weather forecasting took a period of about half a century to advance from empirical to numerical forecast, it has achieved significant success. A consensus has been reached among the Chinese seismological community that earthquake prediction must also develop from empirical forecasting to physical prediction. However, it is seldom mentioned that physical prediction is characterized by quantitatively numerical predictions based on physical laws. This article discusses five key components for numerical earthquake prediction and their current status. We conclude that numerical earthquake prediction should now be put on the planning agenda and its roadmap designed, seismic stations should be deployed and observations made according to the needs of numerical prediction, and theoretical research should be carried out.
Keywords earthquake prediction, geodynamics, numerical method, nonlinear dyamics, weather forecast
1 地震预报的历史回顾
自然灾害的预测和预防，一直是人类关注的重要问题。128 年前，即1884 年，Science 杂志上的一篇题为《天气预报》的文章[1],指出了当时的天气经验预报“几乎都是基于观测事实而不是基于已确立的科学推理。这是不可避免的，因为大气运动是十分复杂的，而大气科学还没有进展到足以解释成功预报它们所需要的的细节”。90 年前，即1922 年，英国的LewisFry Richardson 提出了数值天气预报的思想[2]，试图用数值方法解天气预报的方程组，尽管当时由于计算能力不够而失败了，但是他当年发表的“Weather Predictionby Numerical Process”一文使他成为天气数值预报的先驱。媒体期待着许多计算人员能在一个天气工厂内，像在一个交响乐队指挥下同步进行计算，从而实现世界的天气预报(见图1)。
图1 当时媒体设想的Richardson的全球计算员组成的计算工厂
在科学家的不倦探索下，1950 年，Charney 等终于利用ENIAC 电子计算机实现了第一次天气数值预报。在20 世纪70到80 年代，天气数值预报逐渐发展为能够实用的技术，现在已经成为先进国家天气预报的主要手段。
在103 年前，即在1909 年的Science杂志上也发表了美国地质学家Gilbert 的一篇题为《地震预报》的文章[3]，文中指出“曾经有这样的时代，天气被视为神所操纵。暴风雨和干旱是对人的惩罚和复仇，人们通过祭祀、祈祷来逃避灾难，祭司是神的中介。现在天气被视为自然现象，祭司为气象局所替代，没有祈祷，没有安抚，也没有试图对自然的控制，只是通过科技预测着未来的变化，使服务对象提前得到警告、做好准备：农作物在雨前收获，牛群在洪水到来之前从低地撤离，船舶在风暴之前回到港湾。人们吟唱的圣歌是对科学的赞美。地震也同样曾是笼罩在神秘之中，只有占星家和甲骨文对它有神秘的预测；而现在神秘的阴影让位于知识之光，文明世界的人们期待着地震学家能兴奋地宣布，科学预测地震的时代已经到来。展望地震预测就是我今天讨论的主题。”
遗憾的是，地震预报的进展远远落后于天气预报，虽然100 多年来人们做过不少地震预报的探讨，但是与气象预报已经从经验预报发展为数值预报不同，在很长一段时间内，地震预报没有被科学家认为具备现实可能性，例如，在1935 年的Science 杂志上，当时著名的地震学家Wood 和Gutenberg 指出[4]：“不时的，通过新闻或其他方式，总是有人宣称他们可以预测地震而吸引公众的注意：其中有些宣称者显然是不合格的；有些人的思路的确是合理的，但没有充分考虑到涉及的困难；有些人可能是哗众取宠者；有些人则无疑是真诚的探索者。”他们还提到了我们现在常提到的地震预报三要素：“因为大量地震在地震活动区内经常发生，因此任何含糊的地震预报都会有不错的机会在表观上似乎得到验证，这就要求地震预报应该限定给出面积不太大的地点、确切或接近的时间、以及地震震级。”
地震预报一直在坎坷起伏中发展。1960 年智利Mw9.5 级大地震、1964 年美国阿拉斯加Mw9.2 级大地震后，20 世纪60—70 年代在美国、日本、前苏联等国家曾经掀起了探索地震短期预报的热潮[5,6]。1966 年，著名美国地球物理学家Press 和Brace 写道[5]：“几年前，地震预测仅是占星家、被误导的业余爱好者、追求出名之徒和末日哲学宗教教派的议题，罕见哪个科学家冒着被同行鄙视的风险来就这个问题发表意见。在过去的3年里，这种情况已经发生了巨大的变化，日本最重要的地球科学家提出了一个研究方案，认为再‘经过十多年的数据积累，应该可以开展地震预报’。……在美国科学与技术办公室特设的地震预测调研小组，经过15 个月的审议，认为具备足够的可能性开展一项10 年的研究计划”。岩石在接近破裂时，微裂隙急剧增多引起的扩容现象和水的渗流，以及所导致的地震波速变化，使美国科学家认为找到了有理论基础的前兆手段[7]，乐观的气氛一时占主导。日本从1965 年起就制定了第一个国家“地震预报研究项目”，并在1978 年就已经以不断的监测和评估作为地震预测的正式目标，应对预计不久的将来会发生在本州岛中部的东海大地震[8]。
我国是一个多地震的国家。解放后为适应国民经济的发展，1953 年在中国科学院设立了“中国科学院地震工作委员会”，由李四光、竺可桢兼任正、副主任，1954 年利用中国漫长历史中丰富的地震资料制定了全国地震烈度区划图。1956 年著名地球物理学家傅承义在“科技长远规划”中率先提出在我国开展地震预报研究工作的规划，并提出实现地震预报的科学途径和实施方法。1963 年，傅承义在《科学通报》上发表《有关地震预报的几个问题》的文章[9]。1962 年M6.2级新丰江水库地震和1966 年3 月8 日M6.8级邢台地震之后，地震预报被提到了议事日程，特别是邢台地震造成了8000 余人死亡，周恩来总理多次视察灾区，向地震工作者提出了地震预报的研究课题，地震预报成为我国地震科学的重要研究内容，并建立了国家地震局，成为把地震预报作为日常工作的唯一国家。1971 年，傅承义先生在《地震战线》上发表的文章《关于地震发生的几点认识》[10]，认为“预告的最直接标志就是前兆，寻找前兆一直是研究地震预报的一条重要途径”。文中列举了一些可能前兆，如前震和地下微弱震动、地倾斜和地形变、地磁要素、地震波速度、地下水位、地温、地电、生物，以及月相、气象等要素的变化，并指出：“地震预告是一个极复杂的科学问题”。从1966 年开展地震预报研究以来，在科学工作者的艰苦努力下，我国成功预报了1975年海城地震。
但是地震预报研究遇到的困难远远超出人们的预期。1976 年中国唐山地震短临预报失败，24 万人遇难[11]。美国西海岸大多数地震沿太平洋板块和北美板块边界的圣安德里斯断层发生，特别是在Parkfield附近20 多公里的段落，自有记录的1857年来，每22 年左右重复发生一次6 级左右的地震已经6 次，当时最近一次地震是在1966 年，因为这种极显眼的规律性特点，故被选为地震预报实验场所[12]，并于1985年发布了预期5 年内发生6 级地震的预测，在1992 年甚至发出过72 小时内可能发生4 级地震的短临预报，但直到预报延期到的1993 年底地震一直没有来[13]，直到2004 年(上一次6 级地震之后38 年，预报之后的12 年)才发生6 级地震[14]。而加州1989 年旧金山南面100 多公里Loma prieta的Mw6.9 级地震和1992 年洛杉矶南面200公里的Landers Mw7.3 级地震，也都没有观测到前兆信息。这些失败导致悲观的论点兴盛，其中以Geller 等人1997 年在Nature上发表的《地震不能预报》(Earthquakecannot be predicted)一文[15]为典型代表。他认为[16]，“地震预报的研究已经进行了100 多年并没有取得明显成果，一些宣称获得的突破都未能经得起推敲，广泛的搜索并没有找到可靠的前兆。理论工作表明，断层是一个非线性的过程，它对庞大地球内的不可测量的精微细节都高度敏感，而不只是局限在震源附近。因此，任何微小的地震都有一定的概率级联成一个大事件。对即将发生大地震发布可靠的警报看来实际上是不可能的”。
然而，地震学界并没有停止对地震预报的探讨。尽管更加强调基础研究，但美国、日本和欧洲科学家均在继续进行包括中长期的和短期的预报[17—20]，中国地震科学家更是从来没有放弃地震预报的探讨。不过从20 世纪60 年代至今，地震预报研究已经近半个世纪，虽然地震科学整体研究取得了许多可喜的进展，但在地震预报问题上至今没有取得突破性进展，仍然在原来的水平上徘徊。2008 年的汶川地震，中长期和和短临都没有能实现预报，日本2011 年的东北大地震也没有任何预报，因此地震预报成为科技界和广大民众十分关心的问题。本文将就地震预报研究的发展方向提出一些看法。
2 地震预报的有关概念
在进一步讨论之前，先审视一下地震能否预报的问题。地震震级M以上地震发生总数目N 和M之间存在地震学中称为Gutenberg—Richter 的经验关系，即：log N = a－bM 。Bak 等人[21，22]设计了著名的沙堆模型，在正方形分割的细胞自动机方格内，随机地一粒一粒抛入沙子，代表应力的增加；当一个方格内沙粒数超过规定的临界值时崩塌，令其沙粒数目减少4 个，模拟地震发生和地震时在震源单元产生的应力下降；而相邻四个方格沙粒各增加1 粒，模拟地震引起应力传递使临近区域应力增加。一个方格的崩塌代表小地震，一个崩塌引起一系列方格连锁反应崩塌代表大地震。他们发现沙堆模型同样可以生成这样的幂律崩塌，系统会发展而处于一种高度敏感的不稳定状态，称为临界状态。由于它是在沙子堆积过程中自己逐渐形成的，Bak 等将它称之为自组织的临界状态(SOC)。这时随机地扔沙子，无法估计丢下沙子会导致什么结果：可能停下，可能触发小规模崩塌，但也可能引起极大规模崩塌。在这种状态下，任何规模的崩塌都有可能发生，是完全不可预测的。Barriere 和Turcotte[23]进一步发展三维的由自相似晶格组成的细胞自动机模型，他们不但可以生成Gutenberg—Richter型的分布，而且模拟了小地震有时可以触发大地震，即作为前震；大地震总是可以触发许多小震，即余震。他们也宣称在大地震前观察不到系统性的“前兆”。这些被Geller 等[15]认为是单个地震三要素不能预报的重要根据。
在以上推断中有一个似是而非的陷阱，模拟实验(不管是真实的还是在计算机上虚拟的)固然是认识复杂事物的一种有效方法，但前提是要满足相似条件。沙堆实验可以有助于理解地震数目幂律分布的形成，但在沙堆实验中扔沙子是随机的，不能和地震孕育发生的实际状况相比拟。地质构造板块运动是有规律的，各种研究表明，全球定位系统(GPS)观测得到的板块运动速度在数十年尺度上是均匀的，在更长的地质历史时间尺度上怎样呢？地球科学家发现有一些从地球深部核幔边界升起的热柱，在地表表现为隆起的火山，叫做“热点”。热点的位置大体是固定的，而板块是运动的。太平洋中夏威夷的“热点”在太平洋上产生了帝国—夏威夷火山岛链，从这些岛屿火山岩的年龄和现今离开夏威夷岛的距离，可以知道各个地质历史时期板块运动的速度，它们表明板块运动在数千万年的尺度上也可以是大体均匀的。这样板块运动造成的地块受到的加荷是有规律的，在边界作用下，地块内部应力的变化也是有规律的。尽管板块运动在运动中空间和时间上会有随机起伏，但与完全随机的运动完全是两回事。设想在自由空间中举枪随便乱打，即初始条件中枪的指向是随机的，子弹的落点也将是随机的，无法预测；但是如果是瞄准靶子射击，尽管会有种种随机因素影响，是否能准确击中靶心也难以预测，但是着弹点围绕靶心会形成正态分布，并不是完全不可预测的。
在以上子弹射击成正态分布的比喻中，隐含的条件是处于一个线性系统中。但是人们发现许多实际过程发生于非线性动力学系统中，而非线性动力学系统对初始条件有高度敏感性。即洛伦兹(Lorenz)在用计算机求解仿真地球大气的非线性方程组中发现的如下规律：初始条件很小的差异，可以造成最后结果极大偏差。1979 年12月，洛伦兹在华盛顿一次讲演中提出：“一只蝴蝶在巴西扇动翅膀，有可能会在美国的德克萨斯引起一场龙卷风”的比喻，把对初始值的极端不稳定性形象表达为“蝴蝶效应”。基于非线性系统对初始条件的极端敏感性，系统长时间演化是不能确切预测的。当然，非线性系统的发展演变并不是没有规律的，有时会最终异途同归，演变为一种或有限几种最终状态，似乎状态空间中存在一些“吸引子”；许多情况下更有可能演变成最终状态有无穷多的可能性，但是这些可能状态不是没有规律的，而是在状态空间中“吸引”到无穷多个特殊的点，这些点的空间分布可能具有自相似性，被称为奇异吸引子，达到所谓“混沌”状态。预测非线性系统最终出现的混沌状态是另一种类型的预测。但是即使是非线性动力学系统，在一定的时限内(时限的长短取决于系统的具体特征)，短期预报仍然是可能的。人们可以用数值模式做几天的气象预报已经为实践所证实，人们可以对多时间段内的地震进行什么程度的预报是需要具体工作回答的问题，而不能武断地宣称任何时间段内的地震预报都是不可能的。
目前国外一些学者建议对forecast 和prediction 加以区分，但由于这两个词在英语中本来就是同义词，在气象预报等学科中也被应用，因此在地震预测中究竟如何定义二者的差别没有统一的概念[24,25]。一种理解是prediction 指对一个地震指定时间、地点和震级做的确定性的预报，而forecast是对一定时空窗口内一定震级地震发生概率的估计。也有的理解是prediction是较精确的预报，而forecast指较为不精确的预报。虽然有的学者认为地震不能predict，但也认为地震是可以forecast 的，并开展了地震概率预报研究[26]。笔者认为，一切预报严格来说，都应该以概率预报形式提出，因为罕有预报能有100%的发生概率。其实，现在中国地震学家往往给出确定性的预报，不是因为确定性很高，而是因为他们无法给出事件发生的概率，因此只能以确定性预报形式给出。对于我国地震预报的成功率，不同学者给出的估计从不到10%至30%不等，也就是说预报者虽然给出的是确定性预报，但发生的概率也就是10%左右。应该指出，一定时空窗口内一定震级的地震发生具有一定的背景概率，在地震活动区，背景概率会相当高，因此只有预报概率显著高于背景概率的预报，即获得较高概率增益时，才是有实际意义的预报。而forecast 的概率预报如果能够在将来把时空窗口(特别是短期预报的时空窗口)缩小，概率提高，那么与现在的prediction 又有多少区别呢？纠缠在forecast 和prediction 两个名词上没有太多的意义。重点应该是对预报效能做出科学的评估，行政上可以对一次预报的社会效益做出评价，但是科学上不能以一时一事论英雄，必须对一个单位或个人，一种手段或方法，长期大量预报做出统计和评价[27，28]。
3 地震预报现状和评价
上个世纪60 年代，人们的一种期望是地震是有前兆的。也就是说“有异常”随后就会有地震，“没有异常”随后就没有地震(图2)。那时的认识是，预报不了地震，是因为我们台站不够多，手段不够多，没有能够检测到这些前兆。因此，我国建设了更多的台站，中国大陆现有测震观测站1300 多个。除了测震台网还要建设更多的前兆台网，中国大陆现有形变观测站358 个、电磁观测站294 个、流体观测站619 个，前兆合计观测站1271 个。尽管有那么多的台站，我们还是漏报了汶川地震。我们可以说，我国面积辽阔，台网仍然不够密，需要建设更多的台站，是也需要反思是否仅仅多建台站就可以取得突破进展。
图2 20 世纪六、七十年代人们对地震前兆预报地震的设想(据Turner, Ralph H的《地震预报和公共政策》一文(此文刊登在1976 年出版的《美国国家科学院报告提要》一书第179—202 页上)
陈运泰院士认为地震预报有三个难点[29]：一是地球内部的不可入性；二是大地震事件的非频发性；三是地震孕育发生的复杂性。因为从我国是实际情况来看，固然有这样的情况：“有异常，随之有地震”，比如海城地震的成功预报；但是也有这样的情况：“有异常，随之没有地震”，造成虚报；或是“无明显异常，却随后发生大地震”，造成漏报。
也许，我们没有找到最好的观测和发现地震前兆的手段？人们一直希望找到一种少漏报、少虚报的可靠手段方法，我国地震界也支持和尝试了不少新手段方法，不时也听到有些研究者宣称自己的方法最灵最好，但是遗憾的是，迄今没有得到地震界公认的“灵丹妙药”。也许，虽然任何单一手段都有局限性，但是我们依靠“综合”预报的方法提高预报水平？最简单的一种方法是把各种异常总数随时间变化与大地震的发生找出联系，图3 开始似乎显示了联系，可是，Mw8.0 级的汶川地震却没有遵从经验规律。
图3 1998 年到2008 年我国西部地区上报异常总数随时间的变化，2006 年前异常数目增多似乎与我国或邻近区域一些大地震存在对应关系，但2008 年5 月汶川Mw8 大地震前异常总数反而是下降的(据中国地震局资料)
有的赌徒愿意相信神秘的数字韵律预报？比如，图4 显示了德国“本来应该”获得2010 年世界杯冠军，可惜被章鱼保罗搅了。好在保罗已死，你相信2014年是乌拉圭获得冠军吗？我们在地震预报中恐怕不能寄希望于这样的神秘数字韵律吧。
图4 一些球迷们发现的预测世界杯冠军的神秘数字韵律
有一种说法，短期天气预报容易，长期预报困难；短期地震预报困难，但中长期预报容易。可是图5 显示日本2010 年最新的30 年地震概率预报图，地震学家们一直期待着的是东海、东南海和南海的大地震，突如其来的2011 年3 月的东北大地震使地震学家们面临尴尬。日本从1979年起每年都发布地震概率分布图，此期间发生的一些其他造成十人以上死亡的地震也在图5 中标明，可以看出它们也大多数发生于低概率区。地震预报徘徊不前的现状，迫使我们要从不同的视角重新审视地震预报问题。
图5 日本2010年起30年内地震烈度超过7度的概率预测图，灰色为期待的可能发生大地震的破裂面在地表的投影，棕色为2011年东北大地震实际断裂面在地表的投影；1979年以来发生的一些其他大地震也在图中标明(引自文献[30])
人们往往把地震预报依据预报时间长短分为长期预报、中期预报、短期预报和临震预报。但与气象预报比较，会发现不论预报时段如何，就预报使用的方法，可以分为三种：经验预报、统计预报和数值预报。它们又可以归纳为两大类：第一类为经验预报和统计预报，主要是基于前兆现象的，尽管人们也试图理解前兆的物理成因和规律，但基本依据的是经验，或是对资料的数理统计挖掘和把经验预报以更严格的数理统计形式表达出来。第二类为数值预报，是基于定量化的物理规律的。我国目前的地震预报仅仅依靠基于前兆现象的经验预报和统计预报，没有基于物理规律的地震数值预报。
我国地震预报方法和实践可以参阅文献[31，32]，总的来说，中长期预报更加重视地震活动性、地震地质、空间大地测量等结果；短临预报除了考虑上述资料外，还依赖地形变、电磁、地下水等前兆手段观测结果，一些新的方法、如卫星热红外、次声等也不断有人提出和研究。我国的研究也许受到中医的影响，首先希望“望闻问切”能抓住前兆，预报地震。西方则更希望在布设前兆观测站之前要有理论根据，在企图为地震前兆预测提供坚实理论基础的探讨中，提出许多类型的机制。例如美国的扩容—扩散模型[33]认为，从地壳中岩石的弹性应变积累能量作为开始阶段，随后的第二阶段是岩石中强烈应变部分出现微小的裂隙张开，岩石扩容成为主导因素，标志着真正前兆现象的开始，因为裂缝岩石改变了岩石的物理性质：它使地震波速度(纵波速度与横波速度的比值)下降；岩石是干燥的，则电阻率增加，岩石是湿的，则电阻率下降；通过岩石渗流的水增加(因此有更多的氡气从岩石进入水中)，随着扩容膨胀区的体积增加，小地震的数量有可能会减少，因为在这个阶段，裂隙数目增加会导致水成为不饱和状态，孔隙应力减小会导致滑动摩擦增加，抑制断层发展。第三阶段水会扩散到膨胀(因此水不饱和)的区域，其主要效应是增加地震波速度和提高裂隙中的孔隙水压力，削弱岩石强度，导致小地震的数量增加，随后发生主震。地震发生时造成应力释放，地壳岩石恢复其原有的特性。但是，在大地震出现前的实际观测中，前兆并不是如同理论所设想的那样出现。
经验需要更完善的总结和表达，因此数理统计的方法广泛应用于地震预报。不仅是简单地在前兆和后续地震之间找出一些统计关系，更要通过数理统计模型挖掘出复杂地震活动图像后面隐藏的物理规律。俄国科学家Kelis-Brook 等发展了一系列算法[17] ，利用地震目录进行地震预报，新西兰的David Vere-Jones 等[34]发展了应力释放模型、ETAS模型等描述主震和余震序列的统计特征。
大地震经验预报遇到的重要困难在于大地震的稀少，与气象预报日日月月都积累新经验不同，大地震一人一生在一个地方难得体验几次，小样本统计规律难求。经验本来就少，老经验还经常碰上新问题。
地震预报虽然很复杂，但如果我们在探索一些复杂事物时，能返璞归真，反而有助于理清我们的思路。比如我们掰一根筷子(图6)，作为常识，我们知道它会在上方中央断裂。这是因为，当弯曲时筷子上方中央出现最大张应力，下表面中央出现最大压应力，但一般材料抗压强度比较大，抗张强度比较小，应力超过筷子的强度时就会断裂，因此断裂一般发生在上方。如果我们预先测定了筷子的强度和知道了弯曲的速率，就可以大致预测断裂的位置和时间。
图6 筷子的断裂：(d)图中红色为计算出的张应力最大的部位((b)图中的A)，蓝色为计算出的压应力最大的部位((b)图中的B)
这个简单的例子表明，对于破裂的预测，应力是一个关键量，材料的强度是另一个关键量。应力超过强度就会破坏，在应力增强过程中，还有一些间接的表现，例如小的破裂(微震活动)，材料如果具有磁性和导电性，并有孔隙和含水，那么也可以期待它们随着应力增大在破裂前会有电磁、孔隙水方面的反映，即有多种物理前兆。我们既要关心间接的派生的前兆，更应该关心直接的因素——应力，如果我们能够计算出在特定的加力方式和加力速率下应力的分布，就有可能去预测破裂的发生。
汶川地震预测的失败，引起了中国地震界的科学总结与反思，在总结报告中提出要“……加强地震预报科学研究，不断积累和奠定地震预报的理论基础，逐步从经验预报向物理预报发展。”可以说这是中国地震界的共识。但是地震物理预报的核心是什么，实现地震物理预报需要一个怎样的路线图，对此我们还缺乏广泛的讨论和清晰的认识。
4 基于物理原理的数值预报
数值天气预报是在给定的初始条件和边界条件下，通过对大气流体力学、热力学的一系列偏微分方程组进行数值计算求解，从而得到未来时刻大气的变化和气象要素的分布。它以物理理论为基础，以计算机数学和高速电子计算为实现手段，用地面台站、气球、飞机和卫星来获得特定预报所需要的三维边界条件和初始条件。气象预报的思想在20 世纪20 年代就已经被提出，但当时由于缺乏计算能力而失败了，只有到50 年代开始利用电子计算机后气象预报的研究才取得进展，现在数十个先进国家数值预报已经成为日常天气预报的手段。数值天气预报应该成为数值地震预报的借鉴。
地震的物理预报不能仅仅局限于对地震孕育发生物理过程的定性的了解，物理科学是精确定量的科学，只有开展基于物理规律的数值地震预报研究，才是真正的走上物理预报的途径。如果人们能够建立地下结构和物性的模型，应用连续介质力学、热力学方程和基于岩体破裂准则或断层本构关系，在了解区域边界条件和三维初始应力的条件下，也可以通过高性能计算，了解应力的演变，预测应力由于超过岩体强度而发生地震的最可能位置和破裂类型，预测高应变能积累区的体积大小和未来地震最大的可能震级，根据现有的应力大小和计算的增长速率预测可能发生地震的时段。
实现地震数值预报必须解决5 个关键环节：(1)对物理机制的认识并通过数学公式和数理方程对物理机制进行定量描述；(2)解这些方程的计算能力；(3)对于特定的预报，还要了解所研究区域地下的结构、物性以建立模型；(4)边界条件及其随时间的变化；(5)初始条件。这5 个环节现状如何呢？
图7 显示地震孕育和发生的几个主要过程。地震孕育、地震破裂发生和发展、震后调整及进入新的地震孕育。地震孕育是一个缓慢的过程，需要数十年乃至数千年的时间，其中在临近大地震发生的时段，一般认为岩体可能会出现非线性变形和前兆现象。地震发生时从微小的破裂发展成为大地震，这是一个迅速的过程，需要数秒到数十秒或上百秒时间，但近年发现有所谓的慢地震或静地震[35]，其破裂速度很低，仅有长波地震波辐射或没有地震波的辐射。震后往往还有一个应力调整阶段，在此阶段，断层或周围介质将继续变形，在几个月到几年内，应力仍会有很大的变化，有必要单独考虑，然后才进入下一次大地震的孕育阶段。
图7 地震孕育和发生主要过程图示
对地震孕育破坏的物理规律的了解现状如何？连续介质力学(弹性、粘性、粘弹性、塑性)发展比较成熟，并在地震学中(如地震波传播)和其他科学和工程问题中都已得到成功的应用。对临震和破裂最初发生的许多过程目前了解得还不够，这也是人们对前兆不能很好把握的原因。但对岩体破坏准则，对破裂过程中断层速度和状态相关的本构关系，人们的了解却有相当的进展。因此，可以进行一些初步模拟的工作。
高性能计算能力如何？高性能计算能力近年来不断飞速发展，我们国家的超级计算机和各个科研单位的计算设备也都有相应的迅速发展。国内外都发展了计算所需要的方法和软件，解连续介质的偏微分方程有一种把连续体划分为由若干离散单元组成的近似计算数值解法，过去计算能力有限，只能划分到数万单元，计算精度不高，现在有能力开展百万单元(乃至更大规模)的粘弹性介质、孔隙弹性介质等复杂条件下地球动力学过程的三维初步模拟。
对地下结构物性了解如何？地震学和岩石实验对此问题提供了一定的了解。近年来，我国开展了地壳深部探测专项研究，其成果被评选为2011 年我国十大科技进展新闻之一，反映了我国在这方面投入大幅度加强，资料在迅速积累，表明人们对地球深部的结构、岩性、物理状态和岩体物理力学性质有一定了解，可以构建一些粗线条的模型。
所研究区域的岩石圈侧面和底面边界条件如何？这是一个困难的问题。边界条件，在力学上可以主要分为应力边界条件或位移(速率)边界条件，当然也还有更复杂的边界条件。
应力随深度变化规律比较复杂，岩石圈表层地壳岩石为脆性，可以承受较大的应力，应力随深度增加而增加；但深部在高温高压下，岩石会产生柔性变形，因此能够承受的应力减小；由于地壳复杂分层，有的地方还可能出现随深度而变化的脆性—柔性—脆性的复杂夹层结构，我们很难给出应力随深度变化的边界条件。目前世界最深的钻孔在俄罗斯克拉半岛，深度也仅有12 公里，其次是德国的钻孔，深度为9公里多，其中地应力测量仅做到7公里多，人们对广大区域深部应力现状直接测量了解甚少；对于更深处的应力，人们可以用地震学方法对地震震源机制进行探索，但只能确定主应力方向，不能完整地确定其大小。对应力随时间的变化我们就知之更少。因此应力边界条件比较难以确定。空间大地测量技术近几十年的飞速发展(例如GPS等测量手段)提供了现今地壳构造运动速度的宝贵资料。地质学研究还表明，位移速度在漫长的地质时期(几千万年)可以存在稳定性。例如太平洋中的帝国—夏威夷岛链，是板块运动经过相对固定的热点时火山活动留下的痕迹，火山年龄和它们到现今夏威夷大岛的距离基本为直线关系，反映了运动速度在7000 万年中方向曾经有过一次转折，而速率大体恒定。现今GPS 测量表明，在数十年尺度上，板块运动也存在稳定性，例如印度板块和青藏高原板块之间的相对位置缩短的速率，虽然在苏门答腊大地震后略有增大，但在大多数时间内大体恒定。海洋板块一般可以视为刚性板块，其地表位移可以外推到深部；大陆板块可以发生柔性变形，将地表位移外推到深部时需要更慎重，但在初步的模拟中，不妨假定边界上深部与地表具有相同的速度。
总之，人们需要开展对岩石圈内孕震过程的研究，在研究过程中，虽然研究者对边界条件(特别是深部边界条件)的了解是很差的，但是不排除开展试探研究的可能。
地下应力的初始条件如何，是人们了解最差的问题。我们缺乏对包括大小和方向在内的深部三维应力状况的了解，因此，我们目前无法进行数值地震预报。在短期内，人们也还没有对岩石圈三维应力的大小进行完全探测的能力，地震数值预报也不是短期内能够实现的任务。但是在建立数学方程、计算数学方程的能力和研究地下结构物性三个环节上，人们已经取得一定的进展；在边界条件上也能有一些初步的约束。因此现在开展一些数值地震预报的理论和方法的探讨已经可以提上议事日程。
下面一节仅仅介绍一下地震的孕育、发生和地震带来的影响方面已有的一些工作，其目的是说明，基于数学物理方程模拟计算构造运动造成的地应力变化，有助于我们得到过去定性研讨得不到的新认识，因而数值预报的尝试是有益的。
5 地震孕育、发生和震后调整过程研究中的一些典型事例
关于断层孕育过程的研究，可以以汶川地震为例。一般认为，汶川地震是逆冲、右旋、挤压型断层地震，是由于印度板块向亚洲板块俯冲，造成青藏高原快速隆升，高原物质向东缓慢流动，在高原东缘沿龙门山构造带向东挤压，遇到四川盆地之下刚性地块的顽强阻挡，造成构造应力能量的长期积累，最终在龙门山北川—映秀地区突然释放而产生。但是，为什么汶川地震发生在龙门山之下，而不是更东或更西？为什么汶川地震发生在十余公里深的上地壳底部，而不是更深或更浅？为什么汶川地震西南段以逆冲为主、而东北端为右旋走滑？要回答这一系列问题，需要定量的数值模拟。
为此，首先要构建该区域的三维地质模型。地震层析成像工作提供了该区的深部构造[36](见图8 和图9)，野外和实验室研究使得我们对深部不同层位的岩石岩性和弹性性质有大致的估计，根据地表地热流量密度测量和深部地震波速(与温度有关)，地球科学家可以对深部温度进行估算，根据得到的温度和基于实验室岩石高温高压实验资料，得到不同地区不同深度的粘滞系数分布[37]。应该指出，在高温高压下，岩石不再是弹性体，而是粘弹性体，在短期高频载荷下，表现为固体；但在长期载荷作用下，表现为流体(例如冰看起来是固体，但在长时间尺度下，冰川的冰可以缓慢地像流体一样的流动)。岩石的粘滞系数可以根据实验室高温高压实验确定，在有大震震后变形观测的地区，可以用野外观测资料来验证根据实验室结果计算的岩石粘滞系数是否与观测吻合。
图8 川滇地区的地下结构。图中红色表示该部分地震波速低于这一深度平均值，对于同类岩石，意味着它的温度可能高于该深度平均值；蓝色表示该部分地震波速高于这一深度平均值，对于同类岩石，意味着它的温度可能低于该深度平均值。可以看出四川盆地的下地壳和上地幔为高速区，青藏高原地壳较厚，且下地壳为低速区。(引自文献[36])
图9 川滇地区的地下速度结构剖面图。上图为通过震中垂直于龙门山的剖面，下图为通过震中东西向剖面。可以看出四川盆地的下地壳和上地幔为高速区，青藏高原地壳较厚，且下地壳为低速区。(引自文献[36])
图10 为构建的川西三维粘弹性有限元计算模型，杨氏模量随深度增加而增加，但粘滞系数状况较为复杂，随深度变化为：上地壳粘滞系数高，下地壳粘滞系数低，Moho 面下地幔粘滞系数又较高，更深部的地幔粘滞系数随温度升高而降低。横向高原和盆地结构物性存在显著差异，特征是高原下地壳粘滞系数最低。初始应力状况未知，因此我们只能简单假定初始应力为零，观测应力的增长积累特征。边界条件由GPS 实际观测到的边界点上的位移速度插值确定，假定不同深度的位移速率相同，底面垂直位移为0，水平可以自由滑动。
图10 模拟区域(右下角内红框内部分)的结构分层，注意右下角图显示四川盆地一侧无大地震，高原一侧有较大地震，但最大的汶川地震发生在二者交界的龙门山下。(引自文献[38])
我们习惯于将常温常压下的岩石看作是理想弹性体，因此，先看看假定全部岩石层岩石均为弹性体时候的计算结果(见图11)。这时垂直于断层走向的正应力和沿断层走向右旋的水平剪应力集中部位与图10 的右下角图显示的震中分布没有什么关系。理想弹性模型无法解释汶川地震的孕育成因。
图11 弹性模型的计算结果。图中左边为垂直于断层走向的压应力，右图为沿断层走向的剪应力。红色为应力集中部位，蓝色为低应力部位，弹性模型计算出的应力集中部位与实际地震活动性毫无关系。(引自文献[39])
既然高温高压下岩石是粘弹性的，会产生柔性变形。因此下面展示粘弹性模型的计算结果。由于初始应力状况未知，因此从零应力状态开始，在侧边界上，加上一个从GPS 实际观测插值得到的位移边界条件，并假定底面是水平的，可以自由运动，垂直方向无位移。在恒定边界速度下，粘弹性介质的变形是随时间变化的，对于岩石圈物质的粘滞系数而言，在经过数百年的暂态变化期间后，进入接近稳态的应力增长期，图12 给出的是这时的应力增长率分布情况。图中的红色(应力增长率最高)的区域集中在龙门山断裂带上地壳的底部(10—20 公里深度)，这里发生了汶川8 级大地震；图中的橙色(应力增长率较高)区域位于高原一侧，这里有不少地震发生，其中松潘地震震级可达7 级；图中的浅蓝色为四川盆地上地壳(应力增长率较低)，这里很少有地震活动；图中的深蓝色为地幔(应力增长率最低)，该区域在模拟区域内，地幔柔性变形占主导，基本没有地震发生。注意，图12 左图显示计算的压应力在龙门山断层西南侧最高，与这里断层错动主要为挤压逆掩一致；右图显示计算的剪应力在断层东北段最大，与这里断层右旋走滑发育实际情况一致。
图12 粘弹性模型在进入准稳态后的应力增长率分布图。左边为垂直于断层走向的正应力增长率，右图为沿断层走向的剪应力增长率。(引自文献[39])
考虑岩石在高温高压下具有柔性，因此采用粘弹性模型(高原下地壳粘滞系数最低)，这样计算模拟的结果，可以解释汶川地震多种特征，包括位置、深度、错动方式等，这是过去仅仅依赖定性的讨论无法认识的。在以上模拟中，我们输入的是基于观测的结构模型数据，基于观测和实验室实验的物性数据，以及基于地表GPS 观测的边界条件。不知道的是深部边界条件，只能做一些简化的假定，也不知道初始应力条件，只能从假定的零应力出发计算，分析在数百年后脱离暂态影响后的准稳态应力变化率。得到的是该应力增长率与地震活动性的相关关系。得不到的是应力的绝对值，也无法预测究竟哪里最危险。因为应力增长率高的地方，模型开始计算时刻的初始应力未必一定最高。只知道增量，不知道原来的初始值，就不知道应力的绝对值，无法确定后续地震最可能发生的地点和规模。如果我们对该区地应力状况有充分的了解，就有可能对更具体的地震危险性作出估计。
关于断层错动过程的研究，可以以美国圣安德列斯断层上的Parkfield 段落的地震研究为例。Barbot 等[40]依据速率和状态相关的本构关系研究了该区地震复发周期问题。物理学熟知静摩擦和动摩擦不相等，而对岩石断面的摩擦实验研究进一步表明，断层错动时摩擦力的大小与速度有关，在从一个速度变化为另一个速度时，还有与状态有关的暂态效应。根据1999 年至2010 期间的干涉合成孔径雷达和GPS 数据的实际观测，他们反演构建了Parkfield 段落的断层模型，希望模拟断层上地震破裂孕育、发生、发展的真实过程，实现地震复发的可预测模型。典型结果见图13。
图13 美国圣安德列斯断层上Parkfield 段落地震破裂孕育、发生和传播过程的模拟，试图建立该段落地震复发周期的可预测模型。在1999 年至2010 期间的干涉合成孔径雷达和GPS数据的反演：(a)为断层横剖面不同深度错断的历史，蓝线代表1 年时间间隔内断层两侧缓慢蠕变的累计错断量，红线代表地震发生时1 秒时间间隔内断层地震破裂时的累计错动量,在5—10km深度，有3 次1 秒钟内大幅度的错动，代表了模拟的3 次大地震，红色星号代表震源(破裂起始点)位置；(b)为沿着圣安得利斯断层走向所做的剖面，在7.5 km深度，可以看见一系列的11 个模拟破裂。可以看出有的地震前有小的前震，有时破裂从北端向南端发展，有时破裂发生于南端然后向北发展。6 级地震复发周期为20 年左右(15 到23 年内起伏)，与实际观测吻合。(引自文献[39])
震后应力调整的计算也十分重要，因为这时的应力变化量仍然可以相当于多年的应力积累量。我们计算了日本2011 年Mw9.0 东北大地震对中国东北、华北的影响。对这样大地震造成的位移和应力场变化，必须考虑地球的球体形状，我们利用结构化网格对全球进行划分，对需要重点注意的层位可以对网格加密，结合自适应网格划分技术，对断层附近进行网格加密。网格总数达到400 万个，断层局部可以加密到0.2km，华北、东北地区网格可加密到20—40km(见图14)。图15 给出计算的震位移(即地震时发生的位移，不包含震后断层或岩石蠕变产生的所谓震后变形)变化，和实际观测吻合得很好，达到厘米量级。图16显示中国东部应力同震变化。
图14 有限元计算的网格划分(a)结构化全球网格；(b) 在全球网格基础上，在计算日本东北大地震及其对我国东部影响时采用的自适应加密网格
图15 日本东北大地震在东亚部分地区造成的位移(图中白色为GPS 水平分量观测值，蓝色为本文计算值)
图16 远场位移必须考虑地球球形和分层及横向不均匀性(a)有限元计算的东西向正应力变化；(b) 水平最大主应力和最小主应力变化量
华北地区大体处于NEE—SWW向的挤压应力场之下，尽管存在一些区域差异，但1966 年邢台地震、1967 年河间地震、1969 年渤海地震、1975 年海城地震、1976 年唐山地震，震源机制解均反映NEE向主压应力。日本大地震造成了华北应力场的减小，特别是东西向压应力变化为张应力，因此造成东西向压应力的减小。华北地区年应变率变化在10-9数量级，年东西向压应力增长率约为0.25kPa。而计算的华北东西向压应力从西向东减小了1 到5kPa 不等。换句话说，日本东北大地震造成的华北东西向压应力减小量，大约相当4—20 年构造挤压应力的积累量。这有利于华北在数年到二十年内不再发生与唐山地震等震源机制类似的大地震。过去人们一般根据经验找日本地震与我国华北、东北地震活动性之间的统计关系，而通过计算使我们可以从应力场的变化来进行物理上的讨论。
注意，不能简单地说东西向压应力的减小就一定使华北不发生大地震。应力场的变化是促进地震发生，还是减小了地震危险性，取决于原来的应力场和断层的走向、倾角及可能滑动方向。我们知道驱使断层产生滑动需要的摩擦剪应力τ = μσn，断层上的挤压正应力σn减小,促使滑动的剪应力增大，都会有利于地震发生，反之则会减少地震发生危险性。因此，人们计算被称为库仑应力变化量的参量ΔCFS = Δτ+μΔσn，它为正值并且越大时，表明地震危险性增大，它为负值时，表明地震危险性减小。在计算库仑应力变化量时，必须知道断层走向、倾角和可能的滑动方向(或应力场主压应力方位)。同一个地点的不同走向、不同性质的断层，在临近的同一个大地震后，有的可能变得更危险，有的可能变得更安全。图17 给出了在假定华北主压应力为NEE方向时(对近几十年来华北主要大地震来说，这一假定成立；但不排除局部地区主压应力偏离这一方向)，一些主要断层的库仑应力变化量。可以看出，在这种假定下，日本东北大地震造成绝大多数华北断层更加安全。
图17 假定华北构造应力主压应力方向为NEE向时，不同断层上的库仑应力变化量。1为乌拉山北断层，2 为榆林—府谷断层，3 为山西断陷，4 为五台山前断层，5 为蔚县—延庆断层，6为尚义—怀安断层，7为张家口—北票西断层，8为张家口—北票中断层，9为张家口—北票东断层，10为紫荆关断层，11为太行山前北断层，12为沧东断层，13为宝坻断层，14为唐山—宁河断层，15为益都断层，16为郯庐断
以上例子表明，在地震孕育、发生和震后调整过程中，数值模拟都可以发挥作用，提供我们靠定性分析不能获得的认识。但目前无法进行数值地震预报，其瓶颈是我们不知道初始应力状态，因此虽然可以计算某些情况下应力的变化量，但无法计算应力绝对值，也就无法对地震风险作出评估。初始应必须来自观测，必须进行应力绝对测量和随时间变化量的相对测量，进行地表浅和深部的三维测量，才能推进数值预报的探讨研究。应力观测系统的改进必须有国家级的部署。
6 结束语和展望
和天气预报一样，地震预报如果想取得较大的突破，改变数十年来徘徊的状况，需要在科学思路上有较大的转变，从基于前兆的经验预报、统计预报，发展到基于对地震发生物理基础理解基础上的数值预报[41]。现在应该把数值地震预报研究的规划工作纳入议事日程。应该加强地震数值预报的理论探讨，进行实验基础研究，在继续开展空间大地测量的同时，大力提高和开拓地应力测量方法的研究，并在全国广泛进行地应力绝对测量和布设地应力台站，观测地应力的变化，也对前兆现象的物理成因进行探讨，制定发展地震数值预报长期规划路线图。
即使有了地震数值预报，并不意味着排斥经验预报和数据统计预报。不断有新的技术方法和观测，因此不断有新的感性经验，不断有新的总结，不断有新的理论。理论也有助于我们更深入地了解前兆机制，提高前兆预报的成功率。但是，仅仅依靠积累经验，特别是对于发生频度很低的大地震积累经验是不够的，因为积累经验是缓慢而漫长的过程，获得的经验也充满了不确定性。要改变过去半个世纪难以取得重大突破的状况，必须要有新的科学思想和战略部署。不能仅仅把物理预报限于对地震物理机制的定性研究，而是要开展定量研究并开展数值预报的探讨；不能仅仅泛泛谈谈要逐渐从经验预报向数值预报过渡，而是要有明确的科学思路和切实的路线图及规划；要突破原有的专业分割局限性，注重具有地质、地球物理基础和数值模拟能力的创新人才的培养。
数值预报不是虚幻的梦想，而是现实可行的理想。气象数值预报从提出设想到初步实现走过了大约半个世纪的历程，从地震经验预报到地震数值预报，其道路只会更艰巨和更漫长，因为地球深部物质运动和应力状况的初始条件和边界条件更加难以获取。地震经验预报已经走过了半个世纪，距地震数值预报的目标可能还有世纪之遥，但千里之行始于足下，对于地震数值预报的理论和方法，现在就可以进行探索，不能畏难止步，而要知难而进。虽然现在远不能吹起地震数值预报的“冲锋号”，但我们这一代应该吹响地震数值预报的“起床号”
致谢 感谢中国科学院计算动力学重点实验室许多师生的帮助。谨以此文纪念齐佳老师。
参考文献
[1] Science editor. Science，26 December 1884，568
[2] Richardson L F. Weather Prediction by Numerical Process.
Cartibridge：Cartibridge Univ. Press，1922
[3] Gilbert G K. Science，1909，29(734)：121
[4] Wood H O，Gutenberg B. Science，6 September 1935，
219
[5] Press F，BraceWF. Science，17 June 1966，1575
[6] Hagiwara T，Rikitake T. Science，18 August 1967：761
[7] Scholz C，Sykes L R，Aggarwal Y P. Science，31 August
1973：803
[8] Mogi K. Science，18 July 1986，324
[9] 傅承义. 科学通报，1963，(3)：30
[10] 傅承义. 地震战线，1971，(8)：35
[11] 国家地震局编辑组编. 一九七六年唐山地震. 北京：
地震出版社，1982
[12] Bakun W H，Lindh A G. Science，16 August 1985，
619
[13] Kerr R A. Science，19 February 1993，1120 [14] Bakun W H，Aagaard B，Dost B et al. Nature，2005，
437：969
[15] Geller R J，Jackson D D，Kagan Y Y et al. Science，
1997，275 (5306)：1616
[16] Geller R J. Geophysical Journal International，1997，
131(3)：425
[17] Keilis-Borok V. Annu. Rev. Earth Planet. Sci.，2002，
30：1
[18] Lombardi M，Marzocchi W. Geophys. Res. Lett.，
2009，36：L21302
[19] Uyeda S，Kamogawa M. Eos，Transactions，American
Geophysical Union，2008，89(39)：363
[20] Chadha R K，Pandey A P，Kuempel H J. Geophys.
Res. Lett.，2003，30(7)：1416
[21] Bak P，Tang C，Wiesenfeld K. Physical Review Letters，
1987，59(4)：381
[22] Per Bak，chao Tang. Journal of Geophysical Research，
1989，94(B1l)：15635
[23] Barriere B，Turcotte D L. Phys.Rev.E，1994，49(2)：1151
[24] Kerr R A. Science，30 October 1992，743
[25] Marzocchi W，Zechar J D. Seismological Research Letters，2011，82
(3)：442
[26] Kagan YY，Jackson D D. Geophys. J. Int.，2000，143：438
[27] Shi Y L，Liu J，Zhang G M.Journal of Applied Probability，2001，
38A：222
[28] 石耀霖.中国地震，1992，8(2)：23
[29] 陈运泰.科技导报，2008，26(10)：19
[30] Geller R J. Nature，2011，472：407
[31] 梅世蓉、冯德益、张国民等著.中国地震预报概论.北京：地震出版
社，1993
[32] 张国民，傅征祥，桂燮泰等．地震预报概论．北京：科学出版社，
2001
[33] Press F. Scientific American，1975，232 (5)：14
[34] David Vere-Jones. Journal of Forecasting，1995，11：503
[35] 张晁军，石耀霖. 中国科学院研究生院学报，2005，22(3)：258
[36] 韦伟，孙若昧，石耀霖. 中国科学(地球科学)，2010，40(7)：831
[37] 石耀霖，曹建玲. 地学前缘，2008，15(3)：82
[38] 柳畅，朱伯靖，石耀霖. 地质学报，2012，86(1)：157
[39] 石耀霖. 科学中国人，2012，(11)：20
[40] Barbot S，Lapusta N，Avouac J-P. Science，2012，336：707
[41] 刘启元，吴建春. 地学前缘，2003,10(Suppl)：217
4.20雅安地震
https://blog.sciencenet.cn/blog-1249-682749.html
上一篇：关于地震
下一篇：北京新型冠状病毒感染的肺炎发展趋势分析
