海啸危险性概率分析的不确定性研究＊

刘 也 任叶飞， 温瑞智 王宏伟

1）中国哈尔滨 150080 中国地震局工程力学研究所

2）中国哈尔滨 150080 中国地震局地震工程与工程振动重点实验室

引言

随着发展海域经济的国家战略提升，沿海地区海啸灾害的防灾减灾工作受到更多的关注.海啸危险性评估是海啸防灾减灾的主要手段之一，评估方法分为概率性方法和确定性方法.确定性海啸危险性分析通常关注最不利的海啸场景或最有可能发生的海啸场景，其只能提供单一场景的海啸危险性信息，不能满足多水准的海啸防灾减灾设防需求，而海啸危险性概率分析（probabilistic tsunami hazard analysis，缩写为PTHA）能够提供不同概率水平海啸危险性信息，被更多地应用于沿海地区海啸危险性评估实践工作中，而在PTHA计算过程中存在许多不确定因素，如何量化这些因素对海啸危险性结果的影响是当前热点问题之一（洪明理等，2014；任鲁川，2014；王培涛，2016）.传统的PTHA，对于具有较大随机性的震源参数通常取中位值或发生概率较大的值，此种取值方式往往忽略了地震的随机性，得到的危险性结果只能代表发生概率较大的一种情况，不能满足实际工程建设对海啸危险性结果的需求.因此，量化PTHA中重要参数不确定性的影响对加强危险性结果的科学性，提高准确性具有重要现实意义.

近些年，研究人员已对PTHA的不确定因素展开了讨论.Gica等（2007）对地震海啸的震源参数敏感性展开了研究，认为滑动角、倾角、震中位置和震源深度在一定的范围内变化时对远场海啸波高的影响是有限的，远场海啸波高对滑移量、破裂面积和走向角较为敏感.Necmioğlu和Özel （2014）以地中海为研究区域分析参数敏感性时则提出倾角变化对近场海啸的波高的影响较为明显.Thio等（2014）在分析加利福尼亚地区的海啸危险性时，把倾角、滑移角的不确定性视为随机不确定性，并假设海啸波高在这些参数影响下服从对数正态分布.任叶飞等（2015）分析了海啸数值模拟结果对海洋水深数据的敏感性，指出在开阔的外海海域，不同数据源之间的水深数据差异对于海啸数值模拟的影响可忽略.Li等（2016）指出在PTHA中采用一致的滑移分布模型将会低估海啸危险性水平.上述研究大多针对PTHA中一个或少数几个影响因素的不确定性展开研究，但PTHA是一个包含多个物理过程的复杂系统，涉及参数众多，且这些参数对PTHA结果的影响具有很强的非线性耦合特征，如何同时考虑多个不确定性因素的影响对PTHA的不确定性分析是一个挑战.

为此，本文在归纳总结PTHA中不确定性来源的基础上，拟分别采用逻辑树和事件树方法量化不同类型的不确定性因素，并以马尼拉海啸潜源为研究对象，展示PTHA不确定性分析的计算流程，初步分析PTHA中不确定性因素的影响，以期为中国沿海地区开展全面且系统的PTHA不确定性分析工作提供技术支持.

1 不确定性来源

不确定性根据其性质可以分为随机不确定性（aleatory uncertainty）和认知不确定性（epistemic uncertainty）（Hoffman，Hammonds，1994）.随机不确定性，也称偶然不确定性，是由自然现象内在的随机性所产生的.例如掷骰子，在每次掷骰子前，我们只能知道出现点数的概率分布情况，但无法准确的预测骰子的点数.这种不确定性代表的是一个物理过程或自然现象内在的规律，是无法消除的.随着观测数据的增多，可以采用合理的模型来描述随机不确定性，量化不确定性的影响.认知不确定性，是指人类掌握的知识还有限，现有的方法、原理、模型等还得不到充分的认知，由此产生的不确定性，可以通过研究的深入而缩小其变化范围.随机不确定性和认知不确定性并不是两个独立的类别，而是不确定性的两个状态.随着研究的深入，认知不确定性的规律被逐步发现，从而转为随机不确定性，甚至转为确定性的.两种不确定性概念的区别列于表1.

表1 随机不确定性与认知不确定的比较Table 1 Comparison of aleatory uncertainty and epistemic uncertainty

PTHA源于地震危险性概率分析（probabilistic seismic hazard analysis，缩写为PSHA），PSHA的不确定性分析起步较早，许多学者开展了大量研究（Budnitzet al，1997；McGuire，2004；Delavaudet al，2012），方法和过程相对较为成熟.Kulkarni等（1984）把随机不确定性和偶然不确定性称为固有不确定性 （inherent uncertainty）和统计不确定性（statistical uncertainty），以便更好地与PSHA结合.McGuire （2004）根据参数性质，对PSHA中不确定性参数进行了分类（表2）.但PTHA与PSHA的关键区别在于PSHA通过衰减关系估计地震动参数的大小，而在PTHA中，地震参数与海啸波高的相关关系不是线性的，而是受多个地震参数的影响，海啸波高变化非常复杂，通常只能采用数值模拟的方式估计海啸波高.由于数值模拟所需的算力远高于衰减关系，导致了PTHA的不确定性分析不能完全借鉴成熟的PSHA不确定性分析方法.

表2 PSHA中不确定性分类（引自McGuire，2004）Table 2 classification of uncertainties in PSHA（after McGuire，2004）

将PTHA视为由两个过程组成：海啸的生成过程和海啸的传播过程.由于地质构造的不明确和地震的随机性，海啸产生和传播过程包含了大量的不确定性，既有随机不确定性，也有认知不确定性.海啸传播过程中的不确定性是相对较小的，主要与海啸的传播数值模型和水深数据有关.海啸生成过程的不确定性是PTHA中不确定性的主要来源，其中包含了对复杂的地震破裂过程的模拟和对随机性较强的地震发生位置、地震重现期等参数的确定.表3给出了PTHA中主要涉及的不确定性参数及模型，同时依据其性质进行分类，并给出建议的处理方法.对于随机不确定性，通常采用统计的方法给出参数的概率分布，可以对不确定性进行量化；对于认知不确定性，通常采用逻辑树的方法进行处理.需要说明的是，表中认知不确定性参数和随机不确定性参数之间的界限是变化的，由于现阶段人们对某些自然现象的认识还不充足，导致出现了不同模型描述同一个物理现象，对这些模型的选择引起了认知不确定性，随着对物理现象认识的不断深入，认知不确定性会逐渐消除，转换为随机不确定性.

应当注意的是，表3中列出的不确定性参数仅为目前所能认知的，随着我们对海啸的认识逐渐深刻，可能会出现更多不确定性参数，这些不确定性参数同样可以按其性质分为认知不确定性和随机不确定性，采用本文的方法进行处理分析.另外，在进行PTHA的不确定性分析时，并不是每一个不确定性参数或模型都需要进行分析，不确定性参数的变化对危险分析的结果可能产生较大的影响，也可忽略不计.在实际工程应用中，参数敏感性分析方法可以帮助我们确定对危险性结果影响较大的参数或模型，对其进行重点处理，减少计算量，提高计算效率.虽然目前已有一些研究对PTHA中涉及的部分参数的敏感性展开了研究（Gicaet al，2007；Necmioğlu，Özel，2014；Thioet al，2014；任叶飞等，2015；Liet al，2016），但一些PTHA参数的敏感性会随研究区域的不同而产生变化，在实际应用中仍需针对具体研究区域开展具体分析（Necmioğlu，Özel，2014）.

表3 PTHA中涉及参数的不确定性类型Table 3 Classification of uncertainties in PTHA

PTHA不确定性分析的主要流程如图1所示，由于参数敏感性分析并非本文工作的重点，因此未对PTHA 中不确定性参数按其敏感性分为主要影响参数和次要影响参数.后文将以PTHA中的某几个参数为例，介绍这些参数的不确定性对危险性结果的影响以及相应的分析方法.

图1 PTHA不确定性分析流程图Fig. 1 Illustration of the process to quantifying uncertainties in PTHA

2 不确定性的逻辑树分析方法

地震海啸是由发生在海底的地震使海床产生垂直位移扰动海水引起的.地震引起海床垂直位移的大小决定了海啸初始规模，因此震级被认为是影响地震海啸规模的主要因素之一.2004年印度洋海啸和2011年日本东北海啸发生前，由于人们对这两个地区的地震震级上限缺乏科学的认知，低估了海啸危险性，造成灾难来临时缺乏足够的应对措施和手段.因此，海啸潜源的震级上限是PTHA中的重要不确定性来源之一.由表2可知，震级上限属于认知不确定性，通常采用逻辑树的方法进行分析处理.本节以马尼拉海啸潜源震级上限不确定性为例，给出逻辑树方法分析不确定性的流程.

马尼拉海啸潜源（图2）是中国南海地区规模最大的潜源，由于中国南海地区复杂的地质条件和目前海底断层探测手段的匮乏，马尼拉海啸潜源的震级上限存在较大的认知不确定性.根据历史地震记录，马尼拉海啸潜源周边区域1900年以前未发生过震级大于M8.5的地震，1900年以后发生的地震未超过M8.0 （周本刚等，2011）.特大地震往往具有很长的历史重现期，而历史记载取决于当地文明的进步程度，有些地区从有地震记录起至今也仅有几百年甚至更短，若使用记录到的历史最大地震作为马尼拉海啸潜源的震级上限，显然会低估马尼拉海啸潜源的危险性.若从地震破裂面几何尺寸的角度反推震级上限则会得到令人意外的结果.假设所有海啸潜源分段联合破裂，依据Papazachos等（2004）提出的俯冲带地震破裂尺寸与震级的经验关系，可得到马尼拉海啸潜源的震级上限为MW9.5.Li等（2017）依据地震构造矩率守恒原理估计马尼拉海啸潜源的震级上限为MW8.7，这种方法估计的震级上限考虑到了更多的物理意义，但许多参数的确定过程仍存在很大的不确定性.因此，确定震级上限的方法有多种，存在各自的优势和局限性.下文先对PTHA中的海啸潜源震级上限进行敏感性分析，随后给出针对震级上限不确定性的分析途径.

图2 马尼拉海啸潜源及示例场点位置Fig. 2 The location of Manila potential tsunami source and example sites

针对马尼拉海啸潜源，采用蒙特卡罗（Monte Carlo）采样技术制定四个设定地震样本集，每个样本集的震级下限均设置为MW7.0，震级上限Mmax分别设置为MW8.6，MW8.8，MW9.0和MW9.2，每个样本集包含600次地震海啸事件.模拟每个样本集中的海啸事件生成和传播过程，记录每次地震海啸事件自发生10个小时之内的目标场点波高分布，绘制危险性曲线，分析不同震级上限对危险性结果的影响.PTHA的计算过程非本文关注重点，此处不再赘述，详细可参见Ren等（2017）.

根据Choi等（2012）的研究，海啸波高服从对数正态分布，但对数正态分布自变量的取值为零到正无穷，而海啸波高受海啸潜源震级上限的限制不可能无限变大.因此，我们采用截尾对数正态分布（truncated lognormal distribution）代替传统的对数正态分布对波高分布进行拟合，将设定地震样本集中的最大震级地震引发的海啸在监测点产生的波高最大值作为截尾正态分布的上限，其概率密度函数为

式中：h为海啸波高；μ和σ分别为lnh的均值和标准差；b为海啸波高上限，根据设定地震样本集中最大震级地震引发的海啸在目标场点产生的最大波高确定；φ和Φ为标准正态分布概率密度函数和累积分布函数.海啸波高超越H的概率为

以中国东南沿海某场点为例绘制危险性曲线（图3），由图3可见：采用截尾对数正态分布的结果当波高较小时与采用对数正态分布拟合的结果非常接近；当波高高于上限值时出现了截尾现象，年发生率快速趋于零.

图3 示例场点的海啸波高年发生率曲线（a）与不同重现期的海啸波高值（b）Fig. 3 Annual occurrence rate of tsunami waves exceeding a given height （a）and tsunami wave heights in terms of different return period （b）for the sample site

由图3可见，随着震级上限的增大，海啸波高的年发生率随之增大，相同重现期的海啸波高值也随之增大，震级上限的改变对2 500年重现期的海啸波高影响大于975年重现期和475年重现期的海啸波高，重现期越长，海啸波高对震级上限越敏感.因此，海啸潜源的震级上限改变将对海啸危险性结果产生显著影响，进行PTHA计算时应重点考虑震级上限引起的不确定性.

以不确定性参数震级上限Mmax为例，构造一个简单的逻辑树模型，如图4所示，说明如何使用逻辑树方法处理认知不确定性.逻辑树的每一个分支代表根据不同模型估计的震级上限的取值，并对每个分支分配相应的权重.通常来说，逻辑树每一个分支的权重应当根据每个分支所代表的模型可信度进行分配.但本节的目的在于展示逻辑树方法的计算流程，给出处理认知不确定性的主体思路，故在此不对权重的取值进行讨论，只是简单的进行等权重分配，最终采取加权的方式计算年发生率：

图4 逻辑树处理震级上限不确定性的示意图Fig. 4 Outline of the logic-tree approach processing uncertainity of magnitude upper-limits

式中：Wi为第i分支的权重，此处设每个分支权重为0.25；n为分支数，此处n＝4；vi（h≥H）为每个分支其波高大于H的年发生率.

经过加权计算后的海啸危险性曲线如图5所示.由图可以看出，经过加权后的海啸危险性结果大于仅考虑震级上限为MW8.6和MW8.8情况下的结果，小于仅考虑震级上限为MW9.2情况下的结果，在波高0—5 m范围内与仅考虑震级上限为MW9.0情况下的结果较为接近，在5—10 m范围内年平均发生概率大于仅考虑震级上限为MW9.0情况下的结果.上述分析可见，这种采用逻辑树的方法实际上综合考虑了目前对物理或经验模型的多种认知可能性，可将不确定性的问题转化为确定性的问题.

图5 通过逻辑树方法确定示例场点的海啸波高年发生率Fig. 5 Annual occurrence rate of tsunami waves exceeding a given height by the logic-tree approach for the sample site

3 不确定性的事件树分析方法

在PSHA中，对于随机不确定性，通常在地震动预测方程中以标准差的形式进行量化考虑；而在PTHA中，海啸波高值都是通过数值模拟得到，不能采用PSHA类似的方法.一些研究人员采用蒙特卡罗随机取样方法分析随机不确定性（Sørensenet al，2012；Renet al，2017）.但是如果每个不确定参数都以此种方式分析，其计算成本巨大，且不容易被实现.针对此问题，本节采用事件树方法对PTHA中的随机不确定性进行分析.

事件树方法分析潜在危险带来的损害，最早应用于核电站，逐渐推广至化工、交通安全性评估等领域（Andrews，Moss，2002）.事件树的实质是对随机不确定性参数概率分布的离散化，事件树的每一支都代表随机不确定性参数可能的取值，每一支的权重代表参数值发生的概率.在事件树中从左到右每一条完整的路线都代表一种情况的发生，根据事件树各分支的权重，通过贝叶斯公式，可以得到每一种情况的发生概率.

逻辑树和事件树在结构上是一样的，但在概念上并不同，事件树节点后每一个分支都代表一种情况，有自己的发生概率，逻辑树每一节点后每一支代表一种假设或是模型，每一支的权重代表的是可信度.虽然二者在计算方法的处理方式上相似，但是概念不同.

以马尼拉海啸潜源为研究对象，考虑破裂面长、宽、滑移量、倾角的随机不确定性，利用事件树的方法计算一定保证率下的海啸波高年发生率，为系统的进行PTHA不确定性分析提供一个可供参考的算例.

由于需要确定事件树各分支的发生概率，对于倾角和滑移角首先需要建立其概率分布模型.采用GCMT提供的1976—2015年的马尼拉海啸潜源历史地震震源机制解数据集，对倾角和滑移角进行统计分析，给出其累积概率分布函数.首先对倾角数据进行统计分析，结果如图6a所示.利用对数正态分布函数对离散观测数据进行回归拟合，结果显示拟合结果与观测结果符合程度较好.

图6 马尼拉潜源统计区内历史地震的破裂面倾角（a）和滑移角（b）的累积概率分布Fig. 6 Cumulative distribution function of dip （a）and rake （b）from historical seismic data in Manila regional PTS statistical area

通常认为逆冲型地震具有更大的概率引发海啸.因而，本文仅对滑移角在0°—180°范围内的数据进行统计拟合，结果如图6b所示.图中显示，滑移角在50°—120°范围内的累积概率分布曲线形状较陡，表明滑移角分布主要集中在这一范围内，总体上符合逆冲型地震的破裂滑动特征.采用指数函数拟合其累积概率分布函数（图6b红线），拟合结果较好.

将倾角和滑移角的累积概率分布函数划分为0%—15%，15%—85%和85%—100%三个区间，利用其15%，50%，85%分位数对应值作为区间代表值，根据累积概率分布函数求得每一区间的取值概率，即事件树每一分支的权重，在这里每一个分支对应权重分别为0.15，0.7，0.15.这样既保证了合理的计算量，又能在事件树各分支上体现差别.

对于马尼拉海啸潜源的地震破裂长度和宽度，可通过震级相关的经验关系的平均值和统计标准差来构建其事件树.考虑到如果破裂面长度取最大值且破裂面积同时取最小值时，破裂宽度和滑移量取值会出现过大或过小，不符合实际情况.因此，这里采用破裂面积代替破裂长度和宽度表示破裂规模尺寸的随机不确定性参数.将破裂面积与震级的经验关系 ［ 式（4）］ （Papazachoset al，2004）估算得到的平均值＋1倍标准差、平均值和平均值−1倍标准差作为事件数的三个分支，对应权重为0.15，0.7，0.15.

采用蒙特卡罗采样技术对震中位置和震级进行随机取样，样本容量为200.针对每个样本地震，其破裂滑移角、倾角和破裂面积的取值各有三种方式，依次构建事件树，共有27支.每一支事件树实际上对应一个样本集，其形式如图7所示.由于破裂滑移角、倾角和破裂面积的概率分布是相互独立的，因此每一分支的权重等同于三个参数的概率相乘，最终得到各分支的权重列于表4.

图7 事件树表现形式和各分支节点的权重（分支线上数值）Fig. 7 Outline of the event-tree approach and the weight of the nodes on of each branch

表4 事件树各个分支的破裂倾角、滑移角和破裂面积取值及该分支的的权重Table 4 The vaule of dip，rake and rupture area for each branch and the the weight of each branch

根据Ren等（2017）提出的PTHA计算过程，对这5 400 （200×27）个样本地震进行海啸生成和传播数值模拟，并计算沿海目标场点针对每个事件树分支的海啸波高年发生率.选取中国东南沿海某示例场点作为计算示例，27个事件树分支计算得到的27条危险性曲线如图8a所示.可以看出，27条危险性曲线差异明显.当波高低于0.01 m时，它们之间的差异很微小；当波高为1 m时，年发生率的最大值和最小值之间存在约10倍的差距，并且这种差距随着波高的增大愈加显著.表明破裂面倾角、滑移角和破裂面积的随机不确定性对海啸危险性分析结果的影响较为显著.

图8 （a）采用事件树方法计算得到的某示例场点的海啸危险性曲线；（b）海啸波高超过1 m的年发生率累计权重分布；（c）不同保证率的危险性曲线Fig. 8 （a）Tsunami hazard curves processed by the event-tree method for the sample site；（b）Cumulative weight ofthe annual rate of tsunami wave exceeding 1 m；（c）Tsunami hazard curves with different guarantee rate

将事件树的计算结果应用于PTHA中，计算年发生率

若要计算一定保证率下的海啸波高年发生率，以图7给出的中国东南沿海某场点为例计算海啸波高大于1 m保证率达到20%和80%时的年发生率，具体过程如下：

1）针对每个分支，计算波高大于1 m的年发生率vi（h≥1），即图8a中绿色虚线对应的纵坐标值.

2）把27个vi（h≥1）从小到大依次排列，表示为（h≥1），计算每一个（h≥1）对应的累积权重：

式中，Pj为第j个（h≥1）对应的权重.

按照上述步骤可得到不同保证率的危险性曲线，此处以20%和80%保证率结果作为参照对比加权平均后的预期年发生率曲线.如图8c所示，加权平均后的预期年发生率整体高于20%保证率的年发生率，略低于80%保证率的年发生率.表明，若使用加权平均后的危险性值作为预期危险性结果作为工程抗海啸设计的输入，除极端情况外，基本可满足工程设防需要.

4 结论

本文对PTHA中涉及的不确定参数进行归纳总结，根据不确定参数的性质，分为随机不确定性参数和认知不确定性参数，针对这两类不确定性参数，提出分别使用逻辑树与事件树方法分析不确定性的影响.以马尼拉海啸潜源为例展示PTHA不确定性分析的计算流程，给出处理不确定性的主体思路，得到以下结论：

1）海啸潜源的震级上限的认知不确定性将对海啸危险性结果产生显著影响，进行PTHA计算时应重点考虑震级上限引起的不确定性.

2）地震破裂面倾角、滑移角和破裂面积的随机不确定性对海啸危险性分析结果的影响较为显著，取值不同时，危险性曲线差异明显.

3）事件树方法加权平均后的预期危险性值保证率远高于20%，略低于80%，基本满足大多数工程抗海啸设计要求.

综上所述海啸危险性分析结果存在较大不确定性，为了给实际工程提供科学合理的海啸危险性输入，满足不同设防水平抗海啸设计对海啸强度指标输入的需求，对PTHA开展了不确定性分析，给出不同保证率水平的危险性结果是十分必要的.本文旨在提出一个有效的方法可综合分析多个参数的不确定性影响，只对部分参数开展了不确定性分析、提供了简单算例，在未来的工作中还需将方法应用于我国沿海PTHA实践工作中，给出合理考虑不确定性的沿海地震海啸危险性分析结果.