本发明属于海啸灾害评估技术，尤其涉及一种基于蒙特卡洛随机模拟的多重现期海啸灾害评估方法。
背景技术：
进入二十一世纪以来，全球发生了多次地震海啸灾害，直接经济损失超过2600亿美元。其中，2004年印尼苏门答腊9.2级地震海啸和2011年东日本9.0级地震海啸都给当地的居民生活以及经济发展造成了毁灭性的打击。由此可见，海啸已经成为威胁全球沿海居民生命财产安全最严重的自然灾害之一。
开展海啸灾害风险评估的核心是确定某一评估区域的海啸危险性强度，具体来讲是确定该评估区域可能遭受的最大海啸波幅、海啸爬高、淹没深度等。
目前，国际上普遍采用确定性评估技术路线来确定海啸危险性。确定性海啸风险评估方法主要是根据历史事件推断出最具破坏性的地震海啸源，利用海啸模式研究其传播，爬高以及淹没的过程，给出影响评估区域最不利的淹没情景，并据此评估海啸风险。该方法的优点是简单易行，往往通过一个或者几个场景的模拟即可得到结论，且结论形式简单直观。但该方法也存在一定的不足：很难给出这些所谓“最坏”的场景的发生概率，或者评估区域内发生不同程度海啸灾害的可能性(重现期)。
此外，多参数的复合保守性使得这些“最坏”的场景发生概率极小，甚至是不可能的。图1为业内人士采用马尼拉海沟“最坏的”情景评估南中国海地区的海啸风险的示意图。虽然能够直观的看出哪些区域受海啸灾害的影响最大，但是该结果无法给出该场景的重现期，即不同地区受到海啸灾害影响的概率。
由此，如何合理的给出评估区域不同重现期的海啸风险成为当前需要解决的问题。
技术实现要素：
针对现有技术中的问题，本发明提供一种基于蒙特卡洛随机模拟的多重现期海啸灾害评估方法，该方法能够解决现有技术中风险评估结果偏高及无法给出其对应发生概率的问题。
本发明提供一种基于蒙特卡洛随机模拟的多重现期海啸灾害评估方法，包括：
s1、应用多种震源参数估计方法建立潜在海啸源区地震活动性参数逻辑树；
s2、应用线性海啸数值模拟建立海啸单位源格林函数库；
s3、根据建立的潜在海啸源区地震活动性参数逻辑树，应用蒙特卡洛随机模拟生成随机地震事件集；
s4、根据所述海啸单位源格林函数库和所述随机地震事件集，应用随机滑移量模拟生成海啸波幅集；
s5、根据所述海啸波幅集，对多重现期的海啸灾害分布结果统计及对多重现期的海啸灾害分布结果的不确定性分析。
可选地，所述步骤s1包括：
应用多种震源参数估计方法中的最大似然法拟合tgr关系和地震矩守恒定理建立潜在海啸源区的地震活动性参数的逻辑树。
可选地，应用tgr关系获取海啸源区的震级参数，包括：
所述tgr关系为：
其中，m为地震矩，mt是地震目录完整的临界地震矩，mc是转角地震矩，β是tgr关系的斜率，反映了区域的地震活动性；f(m)是地震矩大于mt的地震发生率；
地震矩m与震级m的近似关系为m＝101.5m+c(2)；
m的单位为牛·米，c为常数；
运用最大似然函数计算最优的参数mc和β，最大似然函数为：
其中，n为地震矩mi大于临界地震矩mt的地震个数；
应用地震矩守恒定理的震源参数估计方法，包括：
根据地质参数估算的地震矩表示为：m＝χμwlv(4)；
其中，χ为地震耦合系数，μ为刚度剪切模量，w为地震带的俯冲宽度，l为断层的长度，v为板块的平均速率；
另一方面，根据tgr关系，地震矩m为：
其中，
at为地震矩大于mt的地震事件的年出现率，γ为伽马函数,c1为常数；
结合公式(4)和公式(5)，当mc＞＞mt时：得到转角地震矩mc：
在获得转角地震矩mc后，根据公式(2)推断出转角震级；
所述tgr关系和地震矩守恒定理采用不同的参数得到多个转角震级值，结合全球俯冲带tgr关系，建立潜在海啸源区地震活动性参数逻辑树。
可选地，在执行步骤s2之前，所述方法还包括：
s2a、将海啸源区域划分为若干个大小相等且滑移量为1米的海啸单位源；单位源的断层几何参数通过历史地震震源机制解和相关地质参数估算方法获得。
可选地，所述s3包括：
根据所述潜在海啸源区地震活动性参数逻辑树中的海啸源参数，运用线性海啸数值模拟计算每个海啸源在评估区域的海啸波幅；
其中，线性海啸数值模式表示为：
式中，η为相对于平均海平面的自由表面位移；
由于海啸传播在近岸不满足线性关系，计算得到的海啸波幅需要通过格林法则换算到近岸的海啸波幅，表示为：
其中ac为近岸输出点的海啸波幅，hc为近岸输出点水深，a0为离岸输出点的海啸波幅，h0为沿岸输出点水深。
可选地，所述s3包括：
采用蒙特卡洛随机模拟基于s1中所建立的逻辑树，在海啸源区范围内随机生成多套未来一定时间段内的地震事件集；
对于每一个独立的地震事件集中的地震事件在海啸源区域内随机分布，且每个海啸源区内的地震事件的震级-频率关系满足tgr逻辑树中的其中一支。
可选地，所述s4包括：
针对随机地震事件集中的每一个地震事件，应用震级-破裂范围公式确定并划定计算海啸波幅所需的单位源；震级-破裂范围经验公式表示为：
mw＝4.868+1.392log10(l)(13)
mw＝4.441+0.846log10(a)(14)
式中mw为地震震级，l为破裂长度，a为破裂面积，根据破裂长度和面积选定该地震事件计算所需的海啸单位源；
应用基于vonkarman函数的地震随机破裂模拟技术随机生成所述地震事件在每个单位源上的破裂滑移量；结合海啸单位源数据库计算每个事件在评估区域的最大海啸波幅用以生成最大海啸波幅事件集。
可选地，所述s5包括：
对于某个地震事件集，假设有n个地震海啸源区对沿岸点的海啸风险有贡献，其中第n个潜在地震海啸源区影响近海沿岸场点h海啸波幅t年的超越概率pn(h≥h)，则场点h海啸波幅t年总得超越概率表示为：
由(15)式推出输出点达到指定海啸波幅hh的重现期th为:
对于某一特定的灾害评估输出点o，在特定的重现期t，n个地震目录中都能通过(16)式计算该点上的最大海啸波幅hi，则对于输出点o在重现期t内的最大海啸波幅为：
此外，通过对不同地震目录所产生的最大波幅hi的分组统计得到不同震源参数对结果的影响程度。
为达到上述目的，本发明还提供一种基于蒙特卡洛随机模拟的多重现期海啸灾害评估装置，包括存储器、处理器、总线以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上述方法任意一项的步骤。
为达到上述目的，一种计算机存储介质，其上存储有计算机程序，所述程序被处理器执行时实现如上述方法任意一项的步骤。
本发明具有的有益效果如下：
第一、由于现有地震目录长度不足(最长的仅仅为100余年)，样本数较少，单纯的利用历史地震资料拟合tgr关系会造成其中的转角震级(图4中的mc，主要作用是对tgr关系高震级部分进行修正，tgr关系中大于此震级的地震事件的重现期呈指数形式迅速增长)偏低，进而低估海啸的风险。本发明中结合地震矩转换定律建立潜在海啸源区tgr关系逻辑树，该定律能够通过一段时间的地震年发生率和俯冲带的区域构造参数(俯冲速率，地震耦合因子，剪切模量等等)来计算tgr关系中的转角震级，由此在逻辑树中给出了转角震级的多种可能，使得评估结果更加准确。
第二、任何完善的评估体系都应该存在不确定分析体系，然而目前应用的确定性海啸风险评估方法是不能够提供结果的不确定性的。本发明中利用蒙特卡洛算法(蒙特卡洛随机模拟)，模拟了海量的地震事件。每个地震事件的位置，震级乃至破裂过程在总体满足一定的规律的基础上都是随机的。这就为不确定性分析提供了足够的样本，也能将不确定性融合到评估结果内。
第三、在海啸数值模拟的过程中，海啸初始位移场模型通常是通过okada提出的弹性位错断层模型计算的。该模型中断层的滑移量决定了初始场的海平面起伏量，对数值模拟的结果影响较大。然而，在目前的海啸风险评估应用中，人们在计算断层的滑移量时通常是假定其均匀分布，并根据震级经验公式计算得出。随着地震有限断层反演技术的发展，人们发现实际地震破裂的滑移量分布是不均匀的，且对区域海啸数值模拟结果的影响很大。本发明应用海啸单位源格林函数库来计算每个地震事件所产生的最大海啸波幅，该方法使得模拟海量随机破裂过程的地震事件所产生的海啸波幅成为可能，既可以增加结果的准确性，同时也可以评估随机破裂过程对结果的不确定性。
由此，采用本发明的方法解决了现有技术中海啸评估结果偏高及无法给出结果对应发生概率的问题，将多种不确定性融合到最终的评估结果中，增加了结果的可信度，同时提高了运行效率，方便决策者针对性的做出防灾减灾部署和城市建设规划。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
图1为马尼拉海沟“最坏的”情景评估南中国海地区的海啸风险的示意图；
图2a和图2b分别为本发明实施例提供的方法流程示意图；
图3为马尼拉海沟历史地震分布；
图4为马尼拉海沟tgr关系示意图；
图5为马尼拉海沟tgr参数逻辑树；
图6为马尼拉海沟海啸单位源划分示意图；
图7为南中国海地区重现期为200年的最大海啸波幅分布示意图；
图8为南中国海地区重现期为500年的最大海啸波幅分布示意图；
图9为南中国海地区重现期为2000年的最大海啸波幅分布示意图；
图10(a)至图10(d)分别为南中国海部分城市海啸波高曲线的示意图(依次为高雄，库里马奥，香港和文昌)。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，将本发明用于评估南中国海海啸风险。
在以下的描述中，将描述本发明的多个不同的方面，然而，对于本领域内的普通技术人员而言，可以仅仅利用本发明的一些或者全部结构或者流程来实施本发明。为了解释的明确性而言，阐述了特定的数目、配置和顺序，但是很明显，在没有这些特定细节的情况下也可以实施本发明。在其它情况下，为了不混淆本发明，对于一些众所周知的特征将不再进行详细阐述。
如图2a和图2b所示，图2a和图2b分别示出了本发明一实施例提供的基于蒙特卡洛随机模拟的多重现期海啸灾害评估方法的流程图，本实施例的方法包括下述步骤：
s1、应用多种震源参数估计方法建立潜在海啸源区地震活动性参数逻辑树。
举例来说，本实施例中可根据地质构造理论，南中国海主要的俯冲带为马尼拉海沟，同时也是该区域的主要潜在海啸源。区域范围内的历史地震如图3所示，本实施例中可应用最大似然法拟合该区域的tgr关系，所述tgr关系为：
其中，m为地震矩，mt是地震目录完整的临界地震矩，mc是转角地震矩，β是tgr关系的斜率，反映了区域的地震活动性；
地震矩m与震级m的近似关系为m＝101.5m+c(2)；
m的单位为牛·米，c为常数。
运用最大似然法计算最优的参数mc和β，最大似然函数为：
应用最大似然法求出该区域的tgr参数，其中β＝0.45，转角震级mc的值为7.5，曲线如图4所示。该关系中的转角震级明显低于全球俯冲带的平均值，主要是由于地震目录过短造成的，因此需要通过其他方法进行修正。
应用地震矩守恒定理的震源参数估计方法，可包括：
根据地质参数估算的地震矩可表示为：m＝χμwlv(4)；
其中χ为地震耦合系数，根据前人研究，该区域的地质耦合系数为0.1-0.3之间；μ为刚度剪切模量，取值为3.0-4.9gpa；w为地震带的俯冲宽度，取值98km；l为断层的长度，取值为1100km；v为板块的平均速率，取值90mm/yr；
另一方面，根据tgr关系，地震矩m可以估算为：
其中，
at为地震矩大于mt的地震事件的年出现率，γ为伽马函数,c1为常数；
结合公式(4)和公式(5)，当mc＞＞mt时：
可得到转角地震矩mc：
在获得转角地震矩mc后，根据公式(2)推断出转角震级。带入不同的参数可以得到修正的四个转角震级值7.9，8.1，8.2和8.4，β值可以选马尼拉区域的tgr拟合结果0.45和全球俯冲带的平均结果0.66。由此建立马尼拉海沟tgr关系逻辑树(如图5所示)。
s2、应用线性海啸数值模拟建立海啸单位源格林函数库。
举例来说，在执行s2之前，可将海啸源区域划分为若干个大小相等且滑移量为1米的海啸单位源；单位源的断层几何参数通过历史地震震源机制解和相关地质参数估算方法获得。
例如，结合马尼拉海沟附近现有的地质参数以及历史震源机制解划分马尼拉区域的海啸单位源(如图6所示)，每个单位源长为100km，宽为50km，滑移量为1米。单位源的断层几何参数可以通过历史地震震源机制解及相关的地质参数估算方法获得。
基于上面的海啸源参数(步骤s1中获取的)，建立南中国海区域的海啸数值模型，地形分辨率为2分，在沿岸选取355个输出点，运用线性海啸数值模拟计算每个海啸源在评估区域的海啸波幅，其中线性海啸数值模式可以表示为：
式中，η为相对于平均海平面的自由表面位移；
由于海啸传播在近岸不满足线性关系，计算得到的海啸波幅需要通过格林法则换算到近岸的海啸波幅，可表示为：
其中ac为近岸输出点的海啸波幅，hc为近岸输出点水深，a0为离岸输出点的海啸波幅，h0为沿岸输出点水深。
s3、根据建立的潜在海啸源区地震活动性参数逻辑树，应用蒙特卡洛随机模拟生成随机地震事件集。
采用蒙特卡洛随机模拟(蒙特卡洛算法)基于s1中所建立的逻辑树，在海啸源区范围内随机生成多套未来一定时间段内的地震事件集。共生成地震目录100个，每个地震目录时间为100，000年，共计地震事件166096个，最大震级为9.05。
s4、根据所述海啸单位源格林函数库和所述随机地震事件集，应用随机滑移量模拟生成海啸波幅集。
针对步骤s3中所产生的随机地震事件集中的每一个地震事件，应用震级-破裂范围公式确定并划定计算海啸波幅所需的单位源。震级-破裂范围经验公式可表示为：
mw＝4.868+1.392log10(l)(13)
mw＝4.441+0.846log10(a)(14)
式中mw为地震震级，l为破裂长度，a为破裂面积，根据破裂长度和面积可以确定选定地震事件计算所需的海啸单位源。
应用基于vonkarman函数的地震随机破裂模拟技术随机生上述地震事件在每个单位源上的破裂滑移量。结合海啸单位源数据库计算每个事件在评估区域的最大海啸波幅用以生成最大海啸波幅事件集。
s5、根据所述海啸波幅集，对多重现期的海啸灾害分布结果统计及对多重现期的海啸灾害分布结果的不确定性分析。
例如，对于某个地震事件集，假设有n个地震海啸源区对沿岸点的海啸风险有贡献，其中第n个潜在地震海啸源区影响近海沿岸场点h海啸波幅t年的超越概率pn(h≥h)，则场点h海啸波幅t年总得超越概率表示为：
由(15)式可以推出输出点达到指定海啸波幅hh的重现期th为:
对于某一特定的灾害评估输出点o，在特定的重现期t，n个地震目录中都能通过(16)式计算该点上的最大海啸波幅hi，则对于输出点o在重现期t内的最大海啸波幅为：
根据上述公式计算南中国海区域内200年，500年以及2000年的最大海啸波幅分布如图7，图8或图9所示，统计得到高雄，库里马奥，香港和文昌四个城市的海啸最大波幅曲线如图10所示，图10(a)至图10(d)分别为南中国海部分城市海啸波高曲线的示意图，例如依次为高雄，库里马奥，香港和文昌城市的海啸波高曲线。
本实施例中运用本发明中结合地震矩转换定律建立潜在海啸源区tgr关系逻辑树，该定律能够通过一段时间的地震年发生率和俯冲带的区域构造参数(俯冲速率，地震耦合因子，剪切模量等等)来计算tgr关系中的转角震级，由此在逻辑树中给出了转角震级的多种可能，使得评估结果更加准确。利用蒙特卡洛算法，模拟了海量的地震事件。每个地震事件的位置，震级乃至破裂过程在总体满足一定的规律的基础上都是随机的。这就为不确定性分析提供了足够的样本，也能将不确定性融合到评估结果内。应用海啸单位源格林函数库来计算每个地震事件所产生的最大海啸波幅，该方法使得模拟海量随机破裂过程的地震事件所产生的海啸波幅成为可能，既可以增加结果的准确性，同时也可以评估随机破裂过程对结果的不确定性。
由此，采用本发明的方法给出了南中国海地区多重现期的海啸风险以及部分城市的海啸波幅曲线。
根据本发明实施例的另一方面，本发明还提供了一种基于蒙特卡洛随机模拟的多重现期海啸灾害评估装置，该装置包括：存储器、处理器、总线以及存储在存储器上并可由处理器执行的计算机程序，所述处理器执行所述程序时实现上述图2a和图2b的方法。例如，
应用多种震源参数估计方法建立潜在海啸源区地震活动性参数逻辑树；例如，应用多种震源参数估计方法中的最大似然法拟合tgr关系和地震矩守恒定理建立潜在海啸源区的地震活动性参数的逻辑树。
应用线性海啸数值模拟建立海啸单位源格林函数库；
根据建立的潜在海啸源区地震活动性参数逻辑树，应用蒙特卡洛随机模拟生成随机地震事件集；
根据所述海啸单位源格林函数库和所述随机地震事件集，应用随机滑移量模拟生成海啸波幅集；
根据所述海啸波幅集，对多重现期的海啸灾害分布结果统计及对多重现期的海啸灾害分布结果的不确定性分析。
另外，本发明实施例还提供一种计算机存储介质，其上存储有计算机程序，所述程序被处理器执行时实现如上述方法任意一项的步骤。
需要明确的是，本发明并不局限于上文所描述并在图中示出的特定配置和处理。为了简明起见，这里省略了对已知方法的详细描述。在上述实施例中，描述和示出了若干具体的步骤作为示例。但是，本发明的方法过程并不限于所描述和示出的具体步骤，本领域的技术人员可以在领会本发明的精神后，做出各种改变、修改和添加，或者改变步骤之间的顺序。
还需要说明的是，本发明中提及的示例性实施例，基于一系列的步骤或者装置描述一些方法或系统。但是，本发明不局限于上述步骤的顺序，也就是说，可以按照实施例中提及的顺序执行步骤，也可以不同于实施例中的顺序，或者若干步骤同时执行。
最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。