引言

海啸是由海底地震、 海底火山爆发、 海岸和海底山体滑坡、 小行星和彗星溅落大洋以及海底核爆炸等产生的具有超大波长和周期的大洋行波(叶琳等，2005)． 历史上大部分破坏性的海啸都产生于深海地震或地震引发的海底滑坡和塌陷. 例如，2004年印度尼西亚苏门答腊海啸及2011年日本海啸都是由海底强地震触发的． 两次海啸都造成了大量人员伤亡和重大经济损失．

自2004年苏门答腊海啸之后，各临海国家再次对海啸预警研究给予了高度的关注． 海啸数值模拟能够给出海啸波的到达时间、 最大爬高和洪水灾害场景，是评估海啸灾害风险及实时海啸预警的重要工具(姚远等，2007)． 目前常用的数值模型有日本Tohoku大学研发的TUNAMI-N2 模型(Imamura，Imteaz，1995)，美国国家海洋和大气管理局海啸研究中心应用的MOST模型(Titov，Gonzalez，1997)，美国Cornell大学开发的COMCOT(cornell multi-grid coupled tsunami model)模型(Liu et al，1998; Wang，Liu，2007)等． 数值模拟时，地球自转产生的科氏力、 频散效应、 非线性项和底摩擦项的设置等都是影响模拟结果(海啸波高和到时)的重要因素(Shuto，1991; Myers，Baptista，2001; Dao，Tkalich，2007; Grilli et al，2007). 此外，如果是海底地震引发的海啸，数值模拟需要可靠的震源机制解作为海啸生成过程的初始化数据． 然而就目前的地震监测和资料分析水平，欲在可供预警的时间内准确获取潜在地震海啸源参数(地震震级、 地震震中、 震源深度、 发震断层的走向、 倾角、 滑动角、 位错、 剪切模量和震源体尺度等)尚较为困难(赵荣国，李卫平， 1997，1998; 任枭等，2006； 张晁军等，2010)． 这已成为地震海啸波高预测失准的主要原因之一，因此分析这些参数对海啸波高的影响具有重要的意义(Gica et al，2007; Synolakis et al，1997; Seda，Tuncay，2010)．

Titov等(2001)以阿拉斯加阿留申群岛地震带发生的地震而引发的海啸为模拟情景，分析最大海啸波高对震级、 震源深度、 断层走向和倾角的敏感性，得到了最大波高对震级很敏感，而震源深度、 断层走向和倾角并不是影响最大波高的主要因素等结论. Gica等(2007)考虑当日本、 阿留申群岛或智利发生地震海啸时，夏威夷近岸的越洋海啸波高对震中位置、 位错、 断层尺度、 断层各种角度的敏感性. 分析结果显示断层尺度、 走向、 位错的变化会引起海啸波高发生明显的变化，并且这种敏感性并不随震源和观测点距离的增大而减弱. Okal和Synolakis(2008)认为越洋海啸波高对断层的尺度不敏感. 但是Seda和Tuncay(2010)利用TUNAMI-N2海啸数值模型模拟越洋海啸波在地中海东部区域传播时，却发现观测点的海啸波高会随着断层面积的增大而增大，并且当震中的位置发生变化时，观测点的海啸波高会发生明显的变化． 他们认为这可能是由于震中位置发生变化引起海啸波传播路径发生改变所致． 类似的分析还可参见任鲁川等(2009)及Geist(1999)等文章． 目前，人们主要是通过数值模拟观察各震源参数不同取值对最大海啸波高的影响，其分析结果却忽略了参数间的交互作用对波高的影响. 而可供比较各震源参数对最大海啸波高影响程度的全局敏感性定量分析方法也尚未建立．

我国南海具有潜在发生地震海啸的风险(叶琳等，2005)，马尼拉海沟沉降带是最具危险性的潜在海啸源(USGS Tsunami Sources Workshop，2006)． Kirby等(2006)参考横切断裂的位置将马尼拉海沟分成6个假想断裂带F₁-F₆(图 1)． 经数值模拟显示，如果马尼拉海沟的北断层F₁(图 1)发生强震引发海啸，将对我国东南沿海和台湾南部产生严重的威胁(Liu et al，2007)． 据此我们拟在马尼拉海沟的北断层F₁发生M_W8.0地震引发海啸的假设情景下，利用COMCOT数值模型和E-FAST(extended Fourier amplitude sensitivity test)法(Saltelli et al，1999)定量分析最大海啸波高对震级、 震中位置、 震源深度、 断层走向、 倾角和滑动角的敏感性及各参数间的交互效应对最大海啸波高的影响． 观测点选为B₁(20.1°N，119.4°E)、 B₂(18.4°N，118.1°E)和B₃(13.5°N，117.6°E)(图 1)，其水深分别为3185，3868和3958 m． 我们的分析方法同样适用于分析其它沿海地区的最大海啸波高对潜在海啸源参数的敏感性．

图  1  马尼拉海沟假想断层F₁-F₆和数值模拟区域. B₁, B₂和B₃为观测点

Figure  1.  Hypothetical fault segments F₁-F₆ along Manila Trench and numerical simulation domain. B₁, B₂ and B₃ are observation points

下载: 全尺寸图片 幻灯片

1.   E-FAST法

E-FAST是Saltelli等(1999)结合Sobol法(Sobol，1993)的优点，对FAST(Fourier amplitude sensitivity test)法(Cukier et al， 1973，1978)进行改进后得到的一种基于方差的全局敏感性定量分析方法． 该方法将模型的敏感性概括为单个输入参数的敏感性及参数间的交互作用的敏感性. 单个参数独立作用的敏感性用主效应指标衡量，参数总敏感性(独立作用和交互作用)用全效应指标衡量．

Sobol法认为，模型可以分解为单个参数和组合参数的函数. 从而模型输出y的总方差V(y)可以分解为

This page contains the following errors:

Below is a rendering of the page up to the first error.

式中，S_xi表示参数x_i独自对模型输出总方差的直接贡献率，其值落在[0，1]区间. 可以依据S_xi的大小对输入参数的重要度进行排序． 定义x_i，x_j的二阶敏感性指标为

它反映了参数x_i，x_j的交互作用对模型输出的影响． 其它更高阶的敏感性指标的定义与此类似． 参数x_i对模型输出的总影响应该包括其独自对模型输出的直接影响以及与其它参数的交互效应． 因此可定义x_i的全效应指标为

式(3)值越大，说明x_i直接和间接对模型输出的总影响越大． 因此我们可以通过全效应指标S_xi^T与主效应指标S_xi的差异来判断x_i与其它参数是否有交互作用．

主效应指标和全效应指标是我们分析模型输入参数(或各组参数)对输出的直接影响及各个参数(或各组参数)的交互效应常用的两种敏感性指标． E-FAST法是一种可以同时计算这两种敏感性指标的常用的敏感性试验方法． 它具有独立于模型(不要求模型具有线性性或单调性)，并能处理输入参数不同取值范围和不同分布形状对分析结果的影响等优点． 具体的计算方法可参见Cukier等(1973，1978)及Saltelli等(1999)的文章.

2.   最大海啸波高对潜在海啸源参数的敏感性试验

根据Liu等(2007)的研究，马尼拉海沟的F₁-F₃断层若发生M_W8.0地震在南海引发海啸，均可能对我国东南沿海和台湾南部造成影响． 尤其是F₁断层(长160 km)距台湾岛最近，产生的海啸将在20分钟左右到达我国台湾南部，并造成重大的影响． 如果F₁断层发生地震，根据Wells和Coppersmith(1994)经验关系，我们可由断层长度推测地震震级的均值为M_W7.77，标准差为0.28． 查阅美国地质调查局官方网站(http://earthquake.usgs.gov)的地震资料，自1897年1月1日至今，马尼拉海沟俯冲区的历史最大震级为M_W7.9． 为了以马尼拉海沟为潜在海啸源在南海展开地震海啸危险性分析，我们参考Liu等(2007)方法，以F₁断层发生M_W8.0地震在南海引发海啸为假想情景进行以下的敏感性分析试验．

假设马尼拉海沟的北断层F₁(长160 km、 宽35 km)发生强震在南海引发海啸，实际震级为M_W8.0，震中位置为20°12′N、 120°30′E，震源深度为15 km，断层走向为10°，倾角为10°，滑动角为90°(Liu et al，2007)． 以下我们将利用E-FAST法定量分析由COMCOT数值模型计算得到的3小时内观测点B₁，B₂和B₃(图 1)的最大海啸波高对震级、 震中位置、 震源深度、 断层走向、 倾角和滑动角测量偏差的敏感性以及各参数的交互效应对最大海啸波高的影响．

2.1   震源参数的分布及抽样

赵荣国和李卫平(1997，1998)依震级均等率和震级偏差率两个参数考察世界各主要地震机构测定震级的偏差分布情况，定量统计分析显示: 中国地震台网震级与国外地震台网震级相比，内陆地区震级测值一般偏大0.1—0.7，太平洋西南部岛弧和海沟地带测值一般偏小0.1—0.4，面波震级偏差约为0.38． 王淑贞和郭履灿(1988)利用χ²检验验证中国台网临时报出的震级与面波均匀震级的偏差总体服从正态分布． 由此我们取马尼拉海沟本次假想地震震级的测量偏差上限为±0.4，并设震级的测量值z服从正态分布，即z~N(8，0.13)(表 1)．

孟玉梅等(2001)指出，中国速报台网1997年海洋地震速报震中与中国地震局地球物理研究所临时报告震中的平均偏差为44.24 km，并且具有逐年减少的趋势． 假设速报震中与实际震中的经向距离为x，纬向距离为y． 本文假设马尼拉海沟本次假想地震的速报震中与实际震中偏差为±40 km，速报震中与实际震中的偏差(x，y)服从二维正态分布，其边际分布分别为x~N(0，13.3)，y~N(0，13.3). 此外，我们假设震源深度的测量值h和断层各个角度的测量值(走向ω₁、 倾角ω₂、 滑动角ω₃)与实际值略有偏差，偏差上限分别取±10 km和±10°，并假设它们在其不确定性范围内服从正态分布． 结合各参数测量的实际值，我们假设各参数的测量值服从表 1的分布，并利用E-FAST法进行反复抽样，得到各震源参数的抽样样本．

表  1  震源参数服从的概率分布

Table  1.  Probability distribution of the earthquake parameters

下载: 导出CSV 

| 显示表格

2.2   海啸波高的数值模拟

式中，M₀为地震矩； μ为剪切模量，单位为N/m²; D为平均位错，单位为m； L为断层长度，单位为m； W为断层宽度，单位为m. 首先可将抽样得到的矩震级换算为位错样本，将抽样得到的速报震中与宏观震中的经向距离和纬向距离样本换算成速报震中的经纬度样本． 然后将各抽样样本及相应的位错样本、 速报震中的经纬度样本分别输入COMCOT模型，模拟南海区域(0°—30°N，100°—130°E)(图 1)海啸波传播3小时，地形和水深数据取自国家地球物理数据中心(National Geophysical Data Center)的ETOPO2数据库. 网格取2′，时间步长取1 s，得到3小时内在各参数抽样样本下B₁，B₂和B₃(图 1)的海啸波高． 数值模拟结果显示，震源参数的偏差对海啸波到时无明显影响，但给观测点的海啸波高带来较大的不确定性． B₁的最大海啸波高的不确定性范围为[0.1，1.7] m． 最大值1.7 m，与由实际参数模拟得到的最大海啸波高0.3 m偏差1.4 m； 最小值0.1 m，与由实际参数模拟结果偏差0.2 m(图 2). B₂的最大海啸波高的不确定性范围为[0.03，0.51] m. 最大值0.51 m，与由实际参数模拟得到的最大海啸波高0.11 m偏差0.4 m； 最小值0.03 m，与由实际参数模拟结果偏差0.2 m(图 3). B₃的最大海啸波高不确定性范围为 [0.01，0.33] m，最大值0.33 m与由实际参数模拟得到的最大海啸波高0.07 m偏差0.26 m； 最小值0.01 m，与由实际参数模拟的结果偏差0.06 m(图 4).

图  2  B₁点的海啸波振幅

Figure  2.  Tsunami wave height at the site B₁

下载: 全尺寸图片 幻灯片

图  3  B₂点的海啸波振幅

Figure  3.  Tsunami wave height at the site B₂

下载: 全尺寸图片 幻灯片

图  4  B₃点的海啸波振幅

Figure  4.  Tsunami wave height at the site B₃

下载: 全尺寸图片 幻灯片

2.3   敏感性指标的计算

结合各参数的抽样样本和观测点的最大海啸波高样本，计算了各参数的主效应指标和全效应指标(表 2)． 通过这些计算，我们得到以下几点认识.

表  2  观测点B₁, B₂和B₃的最大海啸波高对各震源参数的主效应指标和全效应指标

Table  2.  Total and first-order effects of the maximum tsunami wave heights on the source parameters at the sites B₁, B₂ and B₃

下载: 导出CSV 

| 显示表格

1)3个观测点的最大海啸波高都对震级最敏感． 震级测量偏差对观测点B₁，B₂和B₃的最大海啸波高的方差的直接贡献率分别高达51.6%，73.4%和76.7%． 此外，震中位置、 断层倾角、 走向对上述3个观测点的最大海啸波高有较大的影响，但是对最大海啸波高输出方差的贡献率和重要度排序有所不同． 观测点B₁的最大海啸波高对震级最敏感，其次是震中位置、 断层倾角和走向，其全效应指标分别为0.773，0.450，0.262和0.214； 观测点B₂的最大海啸波高对震级最敏感，其次是断层走向、 震中位置和断层倾角，其全效应指标分别为0.879，0.243，0.230和0.115； 观测点B₃的最大海啸波高对震级最敏感，其次为震中位置、 断层走向和倾角，其全效应指标分别为0.892，0.231，0.217和0.130． 而断层滑动角和震源深度在偏差±10°和±10 km的范围内对观测点的最大海啸波高没有显著影响．

2)从主效应指标与全效应指标的差异(表 2)可以看出，敏感的震源参数在影响上述3个观测点的最大海啸波高时，与其它震源参数产生了较强的交互效应． 其中，断层走向、 倾角和滑动角影响这3个观测点的最大海啸波高的主效应指标都远不及0.1，但是全效应指标都超过或接近0.1. 同样，震中位置在影响观测点B₂和B₃的最大海啸波高的主效应指标也不及0.1，但是全效应指标也超过了0.1. 这说明这些震源参数主要是以与其它参数交互作用的形式对最大海啸波高的输出值产生影响．

3)因缺乏考虑参数间相互作用对模拟结果的影响，主效应指标不能很好地反映各参数对模拟结果的影响程度． 而全效应指标则给出各参数变化对模拟结果不确定性的总贡献，更有利于分析海啸波高对各震源参数的敏感性和各参数间的交互效应．

3.   结论

本文以马尼拉海沟的北断层发生M_W8.0地震在南海引发海啸为假想的模拟情景，结合COMCOT海啸数值模型和基于方差分析的E-FAST全局敏感性分析方法，分析了南海3个观测点的最大海啸波高对震级，震中位置，震源深度和断层走向、 倾角、 滑动角等震源参数的敏感性，以及各震源参数间的交互效应对最大海啸波高的影响，分析结果表明：

1)基于方差的E-FAST全局敏感性分析方法能够定量比较各震源参数对最大海啸波高影响的程度，以及分析各震源参数间的交互效应对最大海啸波高的影响，是一种有效的全局敏感性分析方法． 该方法同样适用于分析其它观测点的海啸波高对潜在海啸源参数的敏感性．

2)通过分析我们所选的观测点的最大海啸波高对各震源参数的敏感性，发现不同观测点的最大海啸波高都对震级十分敏感，对震中位置、 断层走向和倾角比较敏感． 但是各震源参数在影响不同观测点的最大海啸波高时，其重要度排序和对输出方差的贡献率存在差异． 因此我们所筛选出来的敏感的震源参数仅适用于我们所选的观测点． 为了得到更全面的海啸预警信息，我们应尽可能地对所有潜在海啸风险的区域进行敏感性分析．

3)由表 2可以看出，各震源参数在影响观测点的最大海啸波高时存在较强的交互效应，其全效应指标超过0.1的震源参数都大于4个． 因此，为了提高海啸预警信息的准确性，我们不能只关注地震速报的内容(即发震时间、 地点和强度)，而应该把海啸预警研究与地震监测预测研究相结合，尽可能地提高各震源参数测量的精度．

本文是在潜源参数测量偏差服从正态分布的假设下建立起的海啸波高对潜源参数的定量全局敏感性分析方法． 这种假设在一定程度上存在局限性． 随着南海地震监测台站的建设和历史记录的增多，今后可能对马尼拉海沟震源参数测量偏差做有针对性地统计分析． 此外，海啸的生成和传播是很复杂的物理过程，对引起不同观测点对潜源参数有不同敏感性的原因也需要我们做进一步深入研究．