CN101788683A - 一种基于多层次互动的海啸运动预测方法 - Google Patents

Abstract

本发明公开了一种基于多层次互动的海啸运动预测方法，该方法包括以下步骤：1)根据实际海啸震源的位置生成层次式网格，建立网格与球面坐标的对应关系及层次网格的交互影响机制；2)对海啸运动方程进行扇形和正方形相结合地离散，根据实际条件设置初始条件和边界条件，求解波动方程来计算海啸运动流通量变化；3)根据海平面的高度设置初始海面高度场，建立海啸运动时粒子网格互动机制和粒子的更新迭代规则，通过粒子与网格的互动来预测海啸的高度。本发明可以快速、逼真地对海啸传播现象进行预测，解决了现有模拟方法只注重于数值仿真、仿真结果不直观、计算效率过慢等问题，用于灾难预防、海啸预警、虚拟仿真等领域，具有一定的实用价值。

Description

一种基于多层次互动的海啸运动预测方法

技术领域

本发明涉及一种灾难预警仿真应用中基于多层次互动的海啸运动预测方法。

背景技术

海啸是一种非常重要的自然灾害，每年给全世界人民带来巨大的损失。在2004年12月26日，罕见的海啸席卷印度洋多个国家，约30万人丧生，500万人痛失家园与生计。这次海啸，可以引起人们许多深思。海啸具有突发性、影响距离远、危害大的特点，而普遍缺乏必要的预警机制和科研经费、信息传达不畅、民众和政府灾害意识淡薄等，则是造成这次巨大灾难的主要原因。据测定，这次海啸在震后3分钟到达科科斯岛，2.5小时到达印度。因此，如能在海啸可能经过的路径上设立由海底计浪器、浮标、卫星、地面接收站等组成的全天候的自动监测网，采用计算机和相关软件技术计算传播时间，就可以及时预报，尽快通知受灾地区的人们提前紧急疏散，以减少受灾程度。

随着计算机技术的发展，对于海啸的模拟方法逐渐从观察预测转换到计算机仿真预测。鉴于海啸的巨大破坏性，国内外对海啸数值模拟和可视化技术也开始重视起来，近20年来，日本和美国等有辽阔海域的发达国家在探索海啸防灾的道路上走在世界前列，亚洲多数国家在此方面研究甚少，以至于在海啸来临之时，酿成本来可以预防的巨大损失。

1984年李成收集和研究了世界大洋中的地震海啸分布规律。1986年，国家地震局组织了60多位专家编著了《地震对策》(郭增建、陈鑫连任主编)一书，书中有关海啸的章节全面论述了地震海啸的危害、成因机制、产生的条件、地震海啸的成灾条件、地震海啸的世界分布、中国地震海啸的可能性，并全面讨论了地震海啸的对策。至目前为止已经发展了一系列越洋海啸的传播、近岸海啸变形等物理模型，包括海浪生成阶段、风浪成长阶段、海浪耗散消亡阶段，并进行了可视化尝试，但是没有绘制出整个海啸动态场景。

20世纪80年代末期日本学者Imamura博士提出了一种数值计算方法，既节约了计算机内存又节省了计算时间，从而使海啸数值计算在海洋工程应用中得到了广泛的普及，并掀起了越洋海啸数值模拟研究的新高潮。1983年Nihonkai-Chubu以地震引起的海啸为例研究了海底运动和海啸运动的规律，对海啸波的推进做了二维色带的可视化和三维曲面的可视化。Ohmachi等考虑浅水方程中部分非线性项，可成功模拟海浪的有效波高。由于海啸物理模型的复杂，海啸场景的绘制也一个难题，Krikke给出了印度洋海啸俯视场景的二维彩色抽象模拟和三维海啸波传播抽象模拟。Kowalick等于2004年提出了对下一代海啸仿真模型的想法，考虑海面障碍物对海啸波浪变形的影响下，可采用二维截面网格和VOF曲面造型方法。2005年，Roger Edberg等通过求解一个受地形、粘滞力等影响的浅水波方程来对海啸数据快速可视化。阿拉斯加海啸预警中心根据海啸偏微分方程组，用有限差分的方式对2004年12月26日印度洋大海啸进行的数值模拟效果。

综上所述，目前国内外关于海啸模拟的工作主要是基于试验及相关数值模拟方法，对海啸发生时与震源、海底地形及海水交互作用的规则考虑还不够，不能很好地预测和计算海啸发生时的运动过程及海浪高度。实验数据的匮乏是造成不能很好预测的原因，同时由于模型过于复杂和网格过于密集，计算量常常较大。

发明内容

本发明的目的在于克服现有方法注重于试验仿真、计算效率过慢、仿真过程不准确等问题，提供一种基于多层次互动的海啸运动预测方法。

本发明的目的是这样实现的：

一种基于多层次互动的海啸运动预测方法，该方法包括以下具体步骤：

1)基于震源互动的海啸区域网格生成

(a)以海啸实际发生的震源位置为中心，建立双向自适应的层次网格，该双向自适应的层次网格初始为大尺度的径向圆形网格，从震源产生时沿扇形传播，然后根据海啸的传播逐级渐远地细分。

(b)采用九邻域的线性插值将粗网格细分为细网格，同时细网格采用双线性插值对粗网格反馈影响，构建上述层次网格的基于线性插值的交互影响规则。

(c)建立多层次网格与球面坐标的对应关系，计算球面坐标下的每个层次网格点的坐标，将基于球面坐标的海啸波动方程映射到上述层次网格上；在球面坐标系下初始地震海啸的控制方程表达为如下形式：

$\frac{\∂\mathrm{\η}}{\∂t}+\frac{1}{R\cos \mathrm{\φ}}[\frac{\∂M}{\∂\mathrm{\λ}}+\frac{\∂\left(N\cos \mathrm{\φ}\right)}{\∂\mathrm{\φ}}]=0$

$\frac{\∂M}{\∂t}+\frac{\mathrm{gh}}{R\cos \mathrm{\θ}}\frac{\∂\mathrm{\η}}{\∂\mathrm{\λ}}=\mathrm{fN}---\left(1\right)$

$\frac{\∂N}{\∂t}+\frac{\mathrm{gh}}{R}\frac{\∂\mathrm{\η}}{\∂\mathrm{\φ}}=-\mathrm{fM}$

其中，λ，φ是经度和纬度，海浪表面高度h＝η(λ，φ，t)+d(λ，φ)，表示波高和海深之和；M和N为经度方向和纬度方向流通量；g是重力加速度；f是耗散参数，W为地球自转角频率；R是地球半径。

2)与海底地形互动的海啸运动计算

(a)对海啸波动方程进行扇形与正方形网格相结合的有限差分离散：先采用极坐标系下的扇形差分网格离散，当海啸传递到离震源较远的区域时，采用正方形网格离散，为减少计算量，对正方形网格采用蛙跃格式。

(b)依据海岸及岛屿设置刚壁边界的初始法线和初始流通量，以及辐射边界条件的零反射效应；在公式(1)中，设海岸和岛屿的固定边界上，法线方向的速度为零，同时流通量为零，也就是M＝0和N＝0。具体计算时，水边界取为辐射边界条件，使边界反射效应近乎于0，陆边界为刚壁边界，即取法向速度。

(c)依据地形陡缓程度调整地形系数，将其设置为海啸波动方程的边界条件，求解出每个离散网格点处海啸流通量变化，得到海啸运动路径和强度衰减值；

海洋地形地貌是影响海啸传播的重要参数，采用公式：

$H\left(x\right)=\frac{H_0\·\mathrm{Yscale}}{\mathrm{\π}}\·\arctan [\frac{1.0}{(540-\mathrm{Ajust}\·400)}\·(x\·\mathrm{Xscale}-\mathrm{XOFFSET})]-\mathrm{\δ}---\left(2\right)$

在(2)式中，H₀是在深海远处的深度，x为自变量，表示位置。XSCALE、YSCALE分别为X、Y方向的放大系数，以恰当描述海洋的深度和广度。δ来调整基准常数以实现将海洋摆在最合适的观察位置。Ajust为地形调整系数，以灵活地调整地形陡缓程度，来计算海啸在不同海底地形下的运动传播变化，这样求解出方程(1)每个离散点的海啸流通量变化。从而得到海啸运动路径和强度衰减值。

3)粒子与网格结合的互动海啸高度预测

(a)根据基准海平面的高度，采用基于海浪谱的初始波动方程，即某时刻任意网格点上海面高度场的连续函数，设置初始的海面高度场。

根据海平面的高度，设置初始的海面高度场。高度场即t时刻在水平位置x(x，y)上，海面高度Z的连续函数，

X(x，y)＝Z(t)                      (3)

这里x，y是指海面的横，纵坐标，Z为当前点的垂直高度，在没有外力作用下，Z的初值为0。当受到外力作用时，平面X发生水平和垂直方向的畸变(波浪)，并会形成一个新的Z的高度和一个v的速度：

$\frac{{\∂}^{2}z}{{\∂}^{2}x}+\frac{{\∂}^{2}z}{{\∂}^{2}y}=\frac{1}{v^2}\frac{{\∂}^{2}z}{{\∂}^{2}t}---\left(4\right)$

(b)依据海啸流通量变化，采用粒子量化海啸对原始海面高度场的影响，应用自适应算法，在粒子之间行进到一定间距内，在海面对应的区域上动态增补新的粒子。

粒子和网格的交互是由一个叫偏离函数ξ(X，t)来完成的；首先由公式(3)可知ΔZ是所有外力综合叠加的效果，当外力作用于海面时，作用区域会产生一个偏离高度，使该区域和水平位置分离，把偏离高度记为ΔZ，那么在t时刻，满足方程：

其中ξi是某一个外力作用于海面产生的偏离效果，而算法核心思想就是将各种外力的综合效果和粒子所处当前网格结合起来，所以先把所有外力的合力传递给当前粒子，而粒子的作用首先是定位其所在的网格，然后将外界吸收的能量传递给该网格，使当前网格按公式(5)产生一个偏离高度，当粒子进入下一个网格，重复以上过程，最终完成整个模拟海啸涌浪的效果。

邻近的网格偏离高度由方程(6)依次获得：

$\left\{\begin{array}{c}{\mathrm{\ξ}}_i(X,t)=A_i\cos (\frac{\mathrm{\π}}{2}\·\frac{X-X\left(t\right)}{L})\·\mathrm{\σ}\\ {\mathrm{\ξ}}_i(X,t)=0\end{array}\right.\begin{array}{c}|X-X\left(t\right)|\≤L\\ |X-X\left(t\right)|>L\end{array}---\left(6\right)$

A_i为当前粒子i产生的振幅，σ为液体粘度。L₀为粒子影响的范围。

(c)根据粒子所在的网格，计算粒子的叠加偏离高度，将因海啸流通量变化而产生的偏离高度与海面高度场相叠加，计算得到该网格点处准确的海浪高度。

本发明的有益效果：

充分考虑了海啸发生时内外部环境的影响因素，建立了互动的海啸运动模型和层次式计算网格，可以真实反应海啸时波浪高度及传播的特性，模拟结果的可信度高。

在进行海啸模拟时还充分采用了层次化和自适应的思想，建立了多层次的自适应网格，这样由粗及精，既保证了计算的精度，又提高了计算效率。从而快速计算出不同情况下的海啸运动。

总之，应用本发明可以快速有效地预测海啸的运行；在用户交互的方便性、计算的快捷性和预测结果的逼真性上，本发明都有显著提高。

附图说明

图1为本发明建立的扇形网格剖分图

图2为本发明流程示意图

图3为本发明计算出的不同情况下海啸振幅变化图

图4为本发明预测仿真的参数化界面图

图5为利用本发明绘制出的海啸互动真实感效果图

具体实施方式

本发明涉及的基于多层次互动的海啸运动预测方法，包括以下实施步骤：

第一步：基于震源互动的海啸区域网格生成

(1)首先我们根据海啸震源情况设置相关初始参数。引起海啸的地震可以在深海也可以在浅海。本实例中，假设海啸发生在距海岸2000km左右远处的深海。由这一点起，海底地形的计算由公式(2)决定。当海啸传播到这一点的时候海啸的传播速度可以计算得：

$v\left(x\right)=\sqrt{g\×H\left(x\right)}$

那么海啸走过的路程等于

$S\left(x_1\right)={\∫}_{x_0}^{x_1}v\left(x\right)\mathrm{dt}$

(2)建立海啸区域的多层次网格。以海啸实际发生的震源位置为中心，建立双向自适应的层次网格，该双向自适应的层次网格初始为大尺度的径向圆形网格，从震源产生时沿扇形传播，然后根据海啸的传播逐级渐远地细分。

海啸在传播过程中是以机械波的形式来传播的，一种是沿地面传播的地震波，它的速度很快，一种是沿海面传播的海啸波，它的速度稍慢。海啸波是一种重力波，其波长与深度密切相关，在深海位置的波长特别长，运动到近岸处，波长变短，频率变大。这里采用公式：

λ(x)＝Const₁+h²*Const₂

其中，Const₁为基准波长，在深海中通常取为200～400Km，在近岸中取为10～50Km，h为x点处的海洋深度，单位为米，Const₂为调整系数，取Const₂＝1.4*10^-5。

采用九邻域线性插值将粗网格细分为细网格，同时细网格采用双线性插值对粗网格产生反馈，构建上述层次网格的基于线性插值的交互影响规则。

建立多层次网格与球面坐标的对应关系，计算球面坐标下的每个层次网格点的坐标，将基于球面坐标的海啸波动方程映射到上述层次网格上；在球面坐标系下初始地震海啸的偏微分方程表达为如下形式：

$\frac{\∂\mathrm{\η}}{\∂t}+\frac{1}{R\cos \mathrm{\φ}}[\frac{\∂M}{\∂\mathrm{\λ}}+\frac{\∂\left(N\cos \mathrm{\φ}\right)}{\∂\mathrm{\φ}}]=0$

$\frac{\∂M}{\∂t}+\frac{\mathrm{gh}}{R\cos \mathrm{\θ}}\frac{\∂\mathrm{\η}}{\∂\mathrm{\λ}}=\mathrm{fN}$

$\frac{\∂N}{\∂t}+\frac{\mathrm{gh}}{R}\frac{\∂\mathrm{\η}}{\∂\mathrm{\φ}}=-\mathrm{fM}$

其中，λ，φ是经度和纬度，海浪表面高度h＝η(λ，φ，t)+d(λ，φ)，表示波高和海深之和；M和N为经度方向和纬度方向流通量；g是重力加速度；f是耗散参数，W为地球自转角频率；R是地球半径。

第二步：与海底地形互动的海啸运动计算

(1)偏微分方程的离散。上述偏微分方程是不易求解的，需要对方程进行有限差分地离散，这里我们采用扇形网格离散和正方形网格离散相结合的方式。海啸传播初期采用极坐标系下的扇形差分网格离散，当海啸传递到较远区域时，这时的网格已经接近于正方形，为减少计算量，对正方形网格采用蛙跃格式。

(2)设置初始条件。设海岸和岛屿的固定边界上，法线方向的速度为零和流通量为零，也就是M＝0和N＝0。具体计算时，水边界取为辐射边界条件，使边界反射效应近乎于0，陆边界为刚壁边界，即取法向速度。

(3)根据海底地形的情况设置边界条件。设置海底地形的阻挡影响。海洋地形地貌是影响海啸传播的重要参数，采用公式：

$H\left(x\right)=\frac{H_0\·\mathrm{Yscale}}{\mathrm{\π}}\·\arctan [\frac{1.0}{(540-\mathrm{Ajust}\·400)}\·(x\·\mathrm{Xscale}-\mathrm{XOFFSET})]-\mathrm{\δ}$

在上式中，H0是在深海远处的深度，x表示位置。XSCALE、YSCALE分别为x、y方向的放大系数，以恰当描述海洋的深度和广度。δ来调整基准常数以实现将海洋摆在最合适的观察位置。Ajust为地形调整系数，以灵活地调整地形陡缓程度，来计算海啸在不同海底地形下的运动传播变化。求解出方程即可得到每个离散点的海啸流通量的变化。

第三步：基于粒子与网格交互互动规则，预测海啸的涌浪高度

先根据海平面高度，确定初步的海面高度场，然后在海啸传播到的网格区域内产生粒子，建立粒子与网格的互动规则，从而得到海啸实际涌浪高度和互动效果。步骤如下：

(1)基于高度场的波面建模

高度场即t时刻在水平位置X(x，y)上，液面高度Z的连续函数，

X(x，y)＝Z(t)；

这里x，y指海面的横、纵坐标，Z为当前点的垂直高度。

当受到外力作用时，平面X发生水平和垂直方向的畸变(波浪)，并会形成一个新的Z的高度和一个v的速度，

$\frac{{\∂}^{2}z}{{\∂}^{2}x}+\frac{{\∂}^{2}z}{{\∂}^{2}y}=\frac{1}{v^2}\frac{{\∂}^{2}z}{{\∂}^{2}t}$

假设t₀时刻，水面高度为Z₀，在t₀+Δt时刻，新高度为Z′，则不难得到：

Z′(t₀+Δt)＝Z₀(t₀)+ΔZ；

Z₀为t₀时刻的初始高度，ΔZ为增量。

(2)粒子与网格的交互

在本发明中，粒子和网格的交互是由一个叫偏离函数ξ(X，t)来完成的。ΔZ是所有外力综合叠加的效果，当外力作用于液面时，作用区域会产生一个偏离高度，使该区域和水平位置分离，把偏离高度记为ΔZ，那么在t时刻，满足方程：

其中ξi是某一个外力作用于液面产生的偏离效果，而我们算法核心思想就是将各种外力的综合效果和粒子所处当前网格结合起来，所以我们先把所有外力的合力传递给当前粒子，而粒子的作用首先是定位其所在的网格，然后将外界吸收的能量传递给该网格，使当前网格产生一个偏离高度，当粒子进入下一个网格，重复以上过程，最终完成整个模拟过程。

必须指出的是，由于在计算过程中，需要对粒子的位置进行整数化处理，以此来确定每个粒子所处的网格，这样做的结果是导致有些网格内包含了多个粒子，叠加之后会产生较大的高度值，而邻近的网格可能没有粒子落入其内，而只保持一个基准高度，从而产生不连续的情况，为了避免这种情况的发生，为每个粒子设置了一个区域控制的权值，使得它们不仅能影响当前网格，也对邻近网格产生一定的效果，而最终产生一个平滑过渡的曲面。

(3)波浪粒子的运动与迭代更新

由于海啸是不断运动的，粒子也需要不断迭代更新。具体实现过程：将所有的粒子按照顺序放入一个链表，每次粒子更新时系统都会检测相邻两个粒子的间距，当检测到P_i，P_i+1的间距大于某一阈值时，则在两个粒子之间再新加入两个子粒子，新的粒子将继承母粒子的所有属性(影响范围，生命等)，只是速度方向和初始位置稍作变化，我们这里取两个母粒子夹角的三等份区域，分别放入两个粒子，并将新的粒子写入链表中，而等下一次更新的时候，子粒子又重新会变成母粒子，允许新的节点插入，如此往返，从而完成粒子的迭代过程。粒子间距的阈值大小直接与当前粒子的影响范围相关，通过测试证明间距大小取粒子影响范围的1/3时效果较好。

真实海啸波浪在运动的时候，有产生和消亡的过程，最终波浪的所有能量会由于分子之间的摩擦，消失殆尽。具体地，为粒子引入了生命属性，让粒子在行进的过程中生命值也逐渐消亡，其对周围网格产生的影响也会随着粒子生命的衰减而减少。

然后根据粒子所在的网格，计算粒子的叠加偏离高度，将因海啸流通量变化而产生的偏离高度与海面高度场相加，计算得到该网格点处准确的海浪高度。

最后还可以采用OpenGL图形开发库绘制地形、海水、海岸等环境，从而模拟真实感的动态海啸视景。

通过以上步骤，可以对各种不同情况下的海啸场景进行快速逼真地模拟。

表一列出了在本实施例中，本发明与Kowalick等人的工作在2.8GHz CPU和2.0G内存的环境下，算法处理质量及效率的对比。其中计算速度是指求解偏微分方程的速度，绘制速度是指显示海啸的速度。可以看到，在同样的环境下，本发明有显著改进。

表一不同方法在2.8GHz CPU和2.0G内存的环境下所需处理时间

方法名称	与环境交互	计算速度	绘制速度
				Kowalick	无	10分钟	5分钟
本发明	考虑多层次交互	2分钟	10秒

以上列举的仅是本发明的具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (1)

1.一种基于多层次互动的海啸运动预测方法，其特征在于该方法包括以下步骤：

1)基于震源互动的海啸区域网格生成

(a)以海啸实际发生的震源位置为中心，建立双向自适应的层次网格，该双向自适应的层次网格初始为大尺度的径向圆形网格，从震源产生时沿扇形传播，然后根据海啸的传播逐级渐远地细分；

(b)采用九邻域的线性插值将粗网格细分为细网格，同时细网格采用双线性插值对粗网格进行反馈，构建上述层次网格的基于线性插值的交互影响规则；

(c)建立多层次网格与球面坐标的对应关系，计算球面坐标下的每个层次网格点的坐标，将基于球面坐标的海啸波动方程映射到上述层次网格上；

2)与海底地形互动的海啸运动计算

(a)对海啸波动方程进行扇形与正方形网格相结合的有限差分离散：先采用极坐标系下的扇形差分网格离散，当海啸传递到离震源较远的区域时，采用正方形网格离散；

(b)依据海岸及岛屿设置刚壁边界的初始法线和初始流通量，以及辐射边界条件的零反射效应；

(c)依据地形陡缓程度调整地形系数，将其设置为海啸波动方程的边界条件，求解出每个离散网格点处海啸流通量的变化，得到海啸运动路径和强度衰减值；

3)粒子与网格结合的互动海啸高度预测

(a)根据基准海平面的高度，采用基于海浪谱的初始波动方程，即某时刻任意网格点上海面高度场的连续函数，设置初始的海面水面高度场；

(b)依据海啸流通量变化，采用粒子量化海啸对原始海面高度场的影响，应用自适应算法，在粒子之间行进到一定间距内，在海面对应的区域上动态增补新的粒子；

(c)根据粒子所在的网格，计算粒子的叠加偏离高度，将因海啸流通量而增加的高度与原始浪高相叠加，计算得到该网格点处准确的海浪高度。

Images