熔岩流是在火山喷发或者地壳断裂过程中涌出岩浆在地表形成的流体，其从根本上说是重力流，由重力作为驱动力，熔岩发源区域实际地形多变，实际流动状态受多种因素影响^[1]，因此对熔岩流过程的理解描述十分复杂，常常需要考虑诸多因素。

在传统的熔岩流过程描述中，Kerr和Castruccio等使用了经验观察和简化的相关实验^{[2, 3]}，描述了熔岩流过程背后的控制因素，Gjerløw等^[4]以实地数据，遥感影像等数据为基础，模拟了熔岩喷发场景；Zuccarello等^[5]分析了熔岩的喷发速度与其覆盖距离的影响关系；Hakim等^[6]使用了相关的实物3D打印模型模拟汉唐江火山群的史前喷发场景；Chevrel等^[7]使用了统计分布的策略，对熔岩流可能影响的相关区域进行评估，揭示了其熔岩流未来入侵位置的分布概率；Mossoux等^[8]使用Q-LAVHA插件，用于模拟数字高程模型上的一个或多个分布式喷发口的熔岩流淹没概率；高超等^[9]使用了基于计算机相关技术的熔岩流模拟方法，对其过程进行仿真。以上研究模拟部分熔岩流过程，描述了不同的熔岩流活动状态，分析了熔岩流过程中不同角度下的影响因素，但这些研究在模拟熔岩流过程方面的可视化效果不够直观，学习者对于影响熔岩流过程的因素了解存在一些障碍，缺乏学习者反复试验并观察过程的机会，难以取得良好的理解效果。

投影式增强现实技术是真实物理环境和完全虚拟的计算机渲染环境之间的桥梁^[10]，不仅可以用于地理方面地图表示的增强现实导航^[11]，也能够为相关地理现象模拟和可视化提供新的手段。本研究依此技术实现了面向地理教学的熔岩流过程可视化，通过深度相机获得可塑地形的高程数据，使用分层设色法直观渲染并体现可塑地形起伏的可视化效果。之后以流体动力学的浅水方程为熔岩流过程建模，计算并模拟可视化随时间变化的熔岩流过程。学习者可以观察并学习地形和这些参数改变后对熔岩流过程的影响，操作方便，可视化效果直观，简化了熔岩流过程理解方式。

1 基于浅水方程的熔岩流动模型

在模拟流体运动的过程中，多种因素都会对其造成影响^[12]，本研究中设定对熔岩流过程主要的影响因素有熔岩流速度，熔岩流喷发体积以及地形因素，以Navier-Stokes方程为起点，使用基于浅水方程的流体模拟方法，通过Navier-Stokes方程推导出二维浅水方程，依此建立熔岩流过程动力模型。

本文采用描述流体运动的二维浅水方程，只考虑熔岩流在水平方向上的运动，降低求解Navier-Stokes方程和数值解的复杂度。在此假设下，熔岩流动是二维的，其垂直运动忽略不计^[13]，w和z可以认为是0，将粘度常数导致摩擦力纳入为一种流体本身的体积力，那么μ可以被设置为0，且压力p与水深h有关，故$ \frac{\partial p}{\partial x}=\rho g \frac{\partial h}{\partial x} $；f_x体积力与x方向上坡角θ关系为f_x=sin(θ)g。结合以上的推导过程与Navier-Stokes方程，得到本研究在二维平面上的浅水方程形式

$ \frac{\partial u}{\partial t}+\frac{\partial(h u)}{\partial x}+\frac{\partial(h v)}{\partial y}=0 $

(1)

$ \frac{\partial(h u)}{\partial t}+\frac{\partial\left(h u^{2}+0.5 g h^{2}\right)}{\partial x}+\frac{\partial(h u v)}{\partial y}=g S_{x} $

(2)

$ \frac{\partial(h v)}{\partial t}+\frac{\partial\left(h v^{2}+0.5 g h^{2}\right)}{\partial y}+\frac{\partial(h u v)}{\partial x}=g S_{y} $

(3)

其中，(u, v) 是流体在t时刻，在点(x, y) 处的速度分量；h为流体深度；(S_x, S_y)分别为(x, y) 方向的地形参数；g为重力加速度。

从式（1）~式（3）可以看出，流体在(x, y) 处随时间t运动的计算方程与三个参数有关：流体的某时刻速度(u, v)、流体深度h以及地形参数S。

2 系统介绍及地形数据渲染 2.1 熔岩流可视化的系统设计

熔岩流可视化过程基于数据输入输出硬件，以及熔岩流过程计算的，基本构成包括有可塑地形、深度相机传感器、投影仪、计算机以及用于支撑的物理支架，其结构如图 1所示。

2.2 熔岩流过程基础地形数据投影渲染

本文事先设立了一个简单的映射规则，用分层设色的方式将获取的高程数据与不同颜色建立映射表，高处地形使用浅棕色渲染，平原地形使用绿色，低洼处使用蓝色，投影渲染前后对比见图 2。可以看出，不同高程的地形点被投影渲染成了不同的颜色，有助于学习者理解熔岩流地形的高低起伏状态。

3 熔岩流过程可视化

前文通过深度相机与投影仪相对位置的标定等地形数据预处理过程，已经保证系统可以获取到所有准确的地形数据作为熔岩流过程中的地形参数S，因此学习者只需要通过设定合适的初始速度u以及熔岩喷发的体积V（影响熔岩流体深度h），就可以对浅水方程中的变量进行控制，并观察不同初始条件对于熔岩流过程的影响。

本研究使用有限体积法求解浅水方程，即对每一时刻的熔岩流过程，将区域内的熔岩流离散成若干控制体单元，计算通过每个控制体边界输入和输出的流量、动量的通量，然后对每个控制体分别进行流体量和动量的守恒计算，从而得到计算时段末各控制体内的平均水深h和流体速度u，作为下一阶段的求解参数。

在本研究中，某时刻熔岩流的格网控制体单元表示如图 4所示。

针对每一帧图像，首先以此深度图像高度和宽度作为深度相机摄得空间域的全部范围；依据熔岩流的初始产生位置生成熔岩流，以一定的像素范围表示此时刻下熔岩流区域，如图 4中的蓝色格网范围所示；最后根据像素大小将熔岩流表示区域内每一个像素范围作为最小控制体单元。

将二维浅水方程以向量形式重写，则有：

$ U_{t}+(E(U))_{x}+(G(U))_{y}=S(U) $

(4)

式中，U是守恒变量向量；E(U)和G(U)是流量分量；S(U)是源项。

$ U=\left[\begin{array}{l}h \\ h u \\ h v\end{array}\right], E=\left[\begin{array}{l}h u \\ h u^{2}+0.5 g h^{2} \\ h u v\end{array}\right], G=\left[\begin{array}{l}h v \\ h u v \\ h v^{2}+0.5 g h^{2}\end{array}\right] $

其中向量内各变量含义同式（4）~式（6）。

对某一个控制单元积分，则有：

$ \frac{\mathrm{d}}{\mathrm{d} t} \int\limits_{\Omega} U \mathrm{d} x \mathrm{d} y+\oint\limits_{\partial \Omega} F \cdot \mathit{\boldsymbol{n}} d \Delta s=S $

(5)

式中，F=(E, G)，为控制体边界的通量；∂Ω为单元的边界；Ω为单元的面积；n为单元边上外法线单位向量；S为源项积分值。图 5表示了不同控制单元之间内部与边界的通量关系

图 5表示了不同格网之间的熔岩流的数值通量变化关系，其中i, j表示不同格网的位置。图中蓝色箭头表示了中心网格i, j向(i, j+1) 以及(i+1, j) 网格的熔岩数值通量F。

$ F_{1}\left(q_{i+0.5, j}, t\right)=\frac{F_{1}\left(q_{i, j}, t\right)+F_{1}\left(q_{i+1, j}, t\right)}{2} $

(6)

$ F_{2}\left(q_{i, j+0.5}, t\right)=\frac{F_{2}\left(q_{i, j}, t\right)+F_{2}\left(q_{i, j+1}, t\right)}{2} $

(7)

其中不同格网之间边界数值通量F近似可以用上述公式计算。使用有限体积法对熔岩流模拟计算的总过程如图 6。

4 实验结果及分析

实验测试了不同参数影响下的熔岩流过程可视化效果，结果表明，根据不同的参数控制，学习者可以观察到不同的熔岩流过程。

4.1 初始速度

图 7与图 8对比体现了在不同初始熔岩流速度值影响下，熔岩流过程的不同。

图 7与图 8的分别表示了当时间同样发展4秒后，设定初始速度u=0.03m/s的熔岩流过程已经拓展至地形的最低处，而u=0.01m/s的熔岩流过程还在地形的平缓处发展，说明熔岩流的初始速度影响了熔岩流过程的发展情况，且可视化的效果明显。

4.2 初始喷发体积

图 9体现了不同初始熔岩流喷发体积V影响下，熔岩流过程的不同。

图 9分别表现了在相同的地形下，熔岩流过程直接影响了方框内的区域，例如当熔岩喷发体积V设为2 cm³时，熔岩仅在地形设计的“火山口”处停留，而随着设定熔岩喷发体积的增多，熔岩开始随着地形自由流动发展，甚至从“火山口”流出，并逐渐在低洼处汇聚，导致了更大的熔岩流影响范围。

4.3 可塑地形

图 10与图 11体现了在不同地形影响下，地形的起伏会影响地形参数S导致熔岩流过程的路径不同。图 10与图 11分别表示了在不同地地形条件下，其他参数相同的熔岩流随时间发展的情况，在地形相对平缓的条件下，熔岩流的流动缓慢地向四周平均发散流动，而当设计地形为山地时，熔岩流从山地处向四周迅速流动，最后汇聚集在山地周围的低洼处。

由于深度相机的边缘处数据获取可能会存在一定误差，因此本研究还针对投影区域边缘进行了投影误差分析。

如图 12所示，由于深度相机摄得深度图像的边缘存在一定畸变问题，会导致地形边缘处的地形参数S获取出现误差。在这种地形获取出现误差的情况下，根据上文中地形参数S与浅水方程的关系可以看出，浅水方程地形参数S的值不准确，直接影响了熔岩流过程状态的计算，此时计算投影所得的熔岩流并未停留在实际地形的最低洼处，投影效果出现了偏差。当获取地形参数S的值偏大时，会导致熔岩流的流向向着地形参数S更大的方向流动。

4.4 各类因素共同影响

为了综合体现出本研究的三种不同参数下的结果可视化，所以本研究进行了两种典型地形情况下的熔岩流过程综合实验。图 13和图 14体现了在各类因素的共同影响下，低速、喷发体积小、在平原地形下的熔岩流过程以及高速、喷发体积大、山地地形下的熔岩流过程。

图 13与图 14表示了两种情况下，各类参数共同影响熔岩流随时间发展的情况。从两图对比可以看出，在低速，初始喷发体积较小的平原地形下，熔岩流过程的发展较为缓慢，并且影响范围较小。而在山地地形下发展的高速剧烈熔岩流过程，向着低洼处流动，并且最终覆盖影响范围较大。

5 结束语

本文利用各类硬件以及可塑地形实现了熔岩流过程的可视化。选择适用于熔岩流特点的基于浅水方程的流体模型，通过塑造不同的地形并选择熔岩流的初始速度，熔岩流喷发体积，最后实现不同情况下的熔岩流径及影响范围的可视化。本文的熔岩流过程可视化直观展示熔岩流过程，分析其受各类因素影响的结果，有助于提高学习者的理解水平，可以作为地理研究和教学的相关工具。