快速检索        
  武汉大学学报·信息科学版  2015,Vol. 40 Issue (10): 1299-1305

文章信息

胡自和,刘坡,龚建华,王群
HU Zihe,LIU Po,GONG Jianhua,WANG Qun
基于虚拟地球的台风多维动态可视化系统的设计与实现
Design and Implementation of Multidimensional and Animated Visualization System for Typhoon on Virtual Globes
武汉大学学报·信息科学版,2015,40(10): 1299-1305
Geomatics and Information Science of Wuhan University,2015,40(10): 1299-1305
http://dx.doi.org/10.13203/j.whugis20130669

文章历史

收稿日期: 2014-04-17

基于虚拟地球的台风多维动态可视化系统的设计与实现
胡自和1,2, 刘坡1,3 , 龚建华3, 王群2    
1. 城市空间信息工程北京市重点实验室,北京,100038;
2. 武汉智慧政通科技有限公司,湖北 武汉,430079;
3. 中国科学院遥感与数字地球研究所,北京,100039
摘要: 台风作为一种热带气旋,是影响我国东南沿海区域的主要灾害性天气。台风可视化是台风研究及预报应用的重要组成部分,在防灾减灾中发挥着重要作用。但是由于数据量巨大,台风动态交互可视化很难在网络环境下实现。虚拟地球技术的出现提供了一个新的网络可视化平台,提高了公众参与的能力,但是其很难支持专业的气象应用。基于开源的虚拟地球平台构建台风动态可视化环境,首先介绍了坐标转换、数据组织和基于GPU的体可视化这些关键技术,然后详细介绍了系统的主要功能,并在World Wind开源平台上实现相关功能的开发,最后通过一个具体的案例证明了系统的有效性和实用性。
关键词: 台风     体可视化     GPU     虚拟地球    
Design and Implementation of Multidimensional and Animated Visualization System for Typhoon on Virtual Globes
HU Zihe1,2, LIU Po1,3 , GONG Jianhua3, WANG Qun2    
1. Beijing Key Laboratory of Urban Spatial Information Engineering,Beijing 100038,China;
2. Wuhan eGOVA Co.,Ltd.,Wuhan 430079,China;
3. State Key Laboratory of Remote Sensing Science,Institute of Remote Sensing and Digital Earth,Chinese Academy of Sciences,Beijing 100039,China
First author: HU Zihe,engineer,specializes in the theories and methods of visualization.E-mail: huzihe06@163.com
Corresponding author: LIU Po,PhD. E-mail: liuposwust@163.com
Foundation support: The National Natural Science Foundation of China,No. 41171351;the Key Knowledge Innovative Project of the Chinese Academy of Sciences,No. KZCX2-EW-318; the National Key Technology R&D Program of China,No. 2014ZX10003002; Foundation of Beijing Key Laboratory of Urban Spatial Information Engineering,No.2014210.
Abstract: Typhoons,a type of tropical cyclone,are the main severe weather event effecting areas on the Southeastern coast of China. Typhoon visualization is a signifiant component in weather research and forecasting applications,and plays an important role in disaster prevention and harm reduction. However,challenges arise when the volume of data is huge,virtual globes can be regarded as a logical platform to visualize such geospatial data over the Internet,but they provide few advanced visualization tools for rendering volumetric data.This paper proposes a virtual globe-based multidimensional and animated visualization system for typhoons. We describe the key technologies,including coordinate transformation,data organization and GPU-based volume rendering. Then,we present the proposed design and implementation in World Wind,introducing the main functions. To demonstrate the capabilities of this system,the data for a simulated typhoon event are rendered on the globe.
Key words: typhoon     volume rendering     GPU     virtual globes    

台风作为一种热带气旋,通常其活动过程会伴随着强风、暴雨以及风暴潮等自然灾害,是影响我国东南沿海区域的主要灾害性天气。台风可视化是台风预报应用的重要组成部分,通过逼真的三维动态模拟和可视化,不仅可以让用户直观地了解台风的运动状态,更能展现台风的运动过程,分析台风的内部结构,从而为台风预报提供支持,在防灾减灾过程中发挥着重要作用。

早期对台风的研究主要集中在台风风场模型及数值预报研究,随着计算机技术的发展,台风可视化越来越引起人们的关注,常用的可视化技术主要有粒子追踪、流线和体可视化等[1]。Van Thu等[2]利用流体模型、粒子追踪和流线来表现台风;朱震毅将二维VBogus风场模型扩展到三维,基于粒子系统来模拟台风的运动过程[3]。体可视化可以观察物体的内部结构,越来越引起人们的关注[4, 5]。Wang等在Vis5D的基础上实现了多维气象数据的粒子、流线和体可视化[6];秦绪佳等在完善传统的光线投射算法的基础上,实现了在虚拟地球环境下的基于GPU台风体可视化[7],但是没有考虑大规模数据量的情况。由于海洋和气象数据的相似性,相关的海洋可视化系统也可以用来实现台风可视化[8]

目前也有一些支持台风可视化的软件,如VAPOR、AVS/Express和Weather 3D等。VAPOR和AVS/Express主要是用作科学可视化,实现了台风的体可视化。Weather 3D则是一款实时气象显示小软件,显示云雾等的三维效果。这些系统主要都是单机系统,很难支持网络可视化[9],对于大规模气象数据的组织和实时调度效果均不理想,而且缺乏公众参与。

自从Gore提出数字地球概念以来[10],越来越多的虚拟地球平台被开发出来,如World Wind (WW)、Google Earth (GE) 等。虚拟地球作为大规模气象数据可视化的理想平台,已经广泛应用于数据共享、分析和空间辅助决策支持[11, 12, 13, 14, 15]。文献[16]基于GE平台采用KML实现了热带气旋数据的共享和查询。这些虚拟地球平台的可视化能力各异,但是很难满足专业的气象可视化的需要,如体可视化[17]。本文采用开源WW平台,扩展其气象可视化能力,实现了台风多维动态数据的交互可视化。

1 关键技术

开源WW平台已经提供了基本的影像、地形和矢量数据的加载功能,并支持常用的点、线、面对象[18],由于其开源特性,用户可以自定义特定的应用功能。本文主要讨论坐标转换、数据组织和GPU的动态体渲染三个关键的技术。首先,虚拟地球平台和其他平台主要的区别是坐标系统的不同,通过坐标转换,可以将现有的功能转移到虚拟地球平台上;其次,虚拟地球很好地支持各种数据服务,但是对体数据的支持能力不强,通过有效的体数据组织,扩展其数据支持能力。此外,体渲染是气象可视化最重要的一个工具,而且需要巨大的计算量。

1.1 坐标转换

虚拟地球平台和传统的渲染平台的主要区别在于其采用球面坐标系,更符合真实的地球形状,传统的数据可以通过坐标系的转换无缝转移到虚拟地球平台中。虚拟地球系统主要包含4个坐标系:笛卡尔坐标系、空间坐标系、球体坐标系、局部坐标系。台风模拟数据一般采用笛卡尔坐标系的网格数据,在可视化的过程中,需要将其转换为球体坐标系和局部坐标系,最后在虚拟地球上渲染出来。不同坐标系之间的转换如图 1所示。

图 1 坐标系转换 Fig. 1 Transformation Between Different Coordinate Systems

1) 笛卡尔坐标与球体坐标的转换。虚拟地球采用WGS84坐标系,需要将其他坐标系的模拟数据转换为球体坐标,详细的转换方法可以参考文献[19],也可以用现成的转换工具,节省开发的时间[20]

2) 球体坐标与空间坐标的转换。在数据地球平台中,通常需要将球体坐标系 (X,Y,Z)转换为空间坐标系(α,β,γ),空间坐标系的原点为地球的中心。球体坐标转换为空间坐标的公式为:

对应的空间坐标转换为球体坐标的公式为:

3) 局部坐标与空间坐标的转换。在渲染过程中,为了提高渲染的精度,采用局部坐标系。如图 1(b)所示,一般取一个几何物体的中心作为坐标原点( X0Y0Z0),在渲染过程中,需要将相关的空间坐标系转换为局部坐标系,它和空间坐标系之间的区别主要在于其系统原点坐标的不同。局部坐标(x,y,z)转换为空间坐标(X,Y,Z)的公式为:

1.2 数据组织

目前的虚拟地球平台已经支持 WFS、WMS和WCS服务,但是这些服务只针对影像、地形和矢量,不支持体数据。气象数据是一种多维多变量体数据,常用的气象数据组织格式为NetCDF和HDF5格式。在文件组织结构上,气象元数据和变量体数据多数是分开存储的,变量读取多依赖于元数据的描述。其庞大的数据量和多时间累计存放的数据存储方式导致很难实现网络的动态可视化,因此需要对原有数据进行有效的组织。针对不同的可视化需求,如用流线来表示风场,只需要生成流线,而对于体可视化,则需要全部的数据,本文采用基于服务的方式,个性化定制不同的功能,从而满足不同的应用需求。

八叉树是一种用于描述三维空间的树状数据结构,能够有效解决体数据的存放问题。本文采用八叉树对气象变量体数据进行预处理,并进行八叉树分块和编码。在构建八叉树之前,需要考虑以下两个问题:①数据在不同的维度差别很大,如高程维只有几十层数据,而在水平维却有几百层数据;②体渲染在纹理边界处不支持双线性或者三线性插值,边界上会存在裂隙[21, 22, 23, 24]。本文提出一种分层八叉树的数据组织策略,如图 2所示,沿高程方向,将原始数据区分为一系列的块,然后对每个块建立多分辨率八叉树,相邻的块重叠一个体素。通过分块,将原始数据分割之后可以满足网络传输和可视化的需要。

图 2 分层八叉树数据组织 Fig. 2 Hieracical Octree Generation and Code

台风数据同时也是多时态的数据,对于N个时次的数据集,会生成N个独立的八叉树,因为数据在相邻时次上存在相关性,在渲染时,使N+1与N时刻中某些相应位置上的子树完全相同,此时只需在N+1处存储指向N中相应子树的指针,相邻时次上的相同子树共享同一个存储空间。分层八叉树不仅可以应用于体可视化,原始的分块数据也可以应用于风场的离散可视化。

1.3 基于GPU的体可视化

体可视化是台风可视化中的一种重要方法,通过它可以观察台风的内部结构,但是传统的体可视化不能直接应用于虚拟地球中。如图 3所示,风场数据以规则的网格数据传输到显存中,但是其实际上存在一定的弯曲,需要进行修正。本文通过计算体数据的范围和采样坐标来修正体可视化结果[25]。假设数据的范围是从最小点 minL(lat0,lon0,H0) 到最大点maxL (lat1,lon1,H1+alt),图 3(b)中的粗黑线表示体数据的范围,由于地球曲率的影响,其高度大小为:

图 3 体可视化 Fig. 3 Volume Rendering on Globe
式中,R0表示地球半径;θ表示数据对应的最大的经度或纬度范围的50%。

采用纹理映射算法渲染体纹理,在GPU中可以得到每个点的局部坐标,先通过式(2)转换为空间坐标,再通过式(3)将空间坐标转换为球面坐标,从而得到每个点的相对高度,通过相对高度可以剔除无效的区域 (A1A2A3)。同时根据经纬度坐标(α,β,γ)计算得到正确的纹理坐标,纹理坐标(u,v,s)的计算公式为:

u、v、s的有效范围为0~1 ,有效范围之外的数据设置为透明。图 4所示为经过坐标转换之后的效果。此外,在多分辨率数据组织的基础上,采用基于视点的多分辨率渲染策略。当体块距离视线比较远时,采用粗略的数据来表现;比较近时,采用比较详细的数据来展示,可以明显地提高渲染的速度。在八叉树构建的过程中,有部分空余的体块,可以采用空值跳跃和早期光线终止策略来加速体渲染的过程[26, 27]

图 4 体渲染效果 Fig. 4 Effect of Volume Rendering
2 系统设计与实现 2.1 软件架构

图 5所示的是台风动态可视化系统的框架,系统主要包括服务器端和客户端。服务器端负责对数据进行读取与预处理,生成XML配置文件,并发布成服务,满足客户端不同的要求。客户端主要负责数据调度和遍历以及常用的数据处理。梯度计算、光照计算和纹理合成等主要的渲染过程都在GPU中完成。

图 5 台风可视化框架 Fig. 5 Framework of Typhoon Visualization
2.2 主要功能

作为一个多维动态可视化平台,其主要目的是模型结果验证和可视化展示,图 6表示台风多维可视化界面,系统包括菜单栏、工具栏、状态栏和主界面窗口。主要功能如下:

图 6 系统交互界面 Fig. 6 Interactive Interface of System

1) 数据存储与预处理,主要包括数据读取、合成、统计和坐标转换等功能,以及八叉树构建和编码,并自动生成相关的XML配置文件,以支持数据集的快速检索以及基于元数据的数据查询与管理。

2) 风场可视化,系统支持台风的二、三维流线的生成和显示,以场线为代表的线表示方式能够很好地揭示流场的走向。

3) 粒子追踪,主要包括二维粒子和三维粒子追踪,通过粒子追踪可以反映台风随时间变化的运动过程。

4) 交互体可视化,用户可以自己设置需要显示的变量,以及该物理量的颜色和透明度,从而分析台风的内部结构等。

5) 剖面分析,主要包括水平、垂直和任意平面剖面等值线图,通过剖面可以得知物理变量在某一高度层的分布情况,用户可以自己定义颜色和剖面分析的位置。

6) 相关的地形、影像的显示。系统支持本地和远程影像服务数据显示,以及相关矢量数据的显示。

2.3 系统实现

实验采用Visual Studio 2008. C# 语言在Windows XP系统下进行开发,采用的显卡是NVIDIA GeForce GT 240,内存为2 G(DDR3 1 333 MHz),处理器为Intel Core i3 530。采用Microsoft DirectX 3D 9.0图形图像处理库,借助High Level Shading Language(HLSL)实现体绘制中的纹理渲染。图 7所示的是气象元素可视化框架,主要包括时间控制和渲染模块,相关的气象元素(如流线、体渲染、粒子)加载到RenderableObjectlists中,优先级别低于矢量数据。利用现有的WW开源平台缩短了开发的时间,降低了开发的难度。

图 7 气象元素可视化框架 Fig. 7 Visualization Framework of Meterological Element
3 案例分析

实验采用2004年8月12日发生在西北太平洋上的台风云娜(Rananim)数据,由欧洲中尺度数值预测中心分析处理的模拟资料。云娜主要的经度范围为85.963 1°~143.963 0°,纬度范围为9.026 2°~45.107 2°,数据格网大小为429×267×19×54,即经度方向划分429行,纬度方向267列,高程上有19层,时间上有54个记录点,模拟起始时间为2004年8月12日12时到14日18时共54个时间步(时间分辨率为1 h),数据大小为6 GB。这里只简单地介绍粒子追踪、流线、体可视化和剖面分析4个常用的功能。

3.1 粒子追踪

粒子追踪一词来源于流体力学,流场中的某个粒子可根据局部速度场来预测粒子未来的位置。Konikow等 [28]提出了一种二维场中的粒子追踪算法,能够预测粒子在二维流场中的未来位置。将粒子追踪算法扩展到三维空间,采用矢量箭头表示风场的特性,箭头的方向代表该点风的方向,箭头长短代表风速,同时用不同的颜色表示不同的风场速度。

图 8(a)中,绿色圆圈标出的区域是台风眼所在的位置,绿色圆圈所在的位置箭头短小且偏白色,说明此处风力很小,绿色圆圈周边箭头长且为深红色,说明此处风力大,箭头符号形成近似于圆形的环流。橙色圆圈标出的位置是与台风伴生的反气旋,其周边符号短小,也形成一个微弱的椭圆形环流,且环流方向与台风处的环流方向相反,符合台风的运动特点。由于巨大的数据量,很难一次表示所有的数据,图 8(b)8(c)8(d)表示第3层在不同分辨率下的可视化结果。可见,随着视点离风场越近,风场表现越详细。

图 8 风场可视化效果图 Fig. 8 Effect Chart of Wind Field Visualized by the Arrows
3.2 流线

台风风场具有流场的特征,目前主要用数值积分法和双流函数法来构造流线,本文采用经典的Runge-Kutta数值积分法[29]。流线可以从一个初始的起点通过一系列小的时间步长“生长”而成。由于台风气旋是向外扩散的,将流线种子点选在台风中心附近,则流线会随着气流向外扩散,反映出台风风场的走势,因此选择台风眼附近的位置作为流线的种子点是进行风场可视化的最佳选择。台风眼通常在台风中心平均直径约为40 km的圆内,数据中通常没有直接的台风中心数据。在离心力的作用下,外面的空气不易进入到台风的中心区内,它里面的空气几乎不旋转,风很微弱,风眼区压强是整个台风区域最低的,因此通过搜索风场范围内压强的最小值即可找到台风眼所在的位置,在台风眼附近给定范围内生成给定个数的随机点,利用数值积分算法,则可以产生相应的三维流线,用流线模拟风场。图 9表示从风眼位置作为种子点开始生成的三维流线效果图,图 9(a)表示流线的初始状态图,台风登陆之后,风场会向四周扩展,图 9(b)表示扩散之后的效果图。此外,由于流线不好表示方向,三棱锥和示踪小球常用来辅助流线展示,图 9(c)9(d)表示对应的效果图,从三棱锥的方向可以看出台风的方向,相对于三维流线本身来说是静止不动的,而示踪小球表示三维流线运动的效果。

图 9 三维流线效果图 Fig. 9 Effect Chart of Three Dimensional Streamline
3.3 体可视化

台风数据主要包含风场、压强、水汽含量等台风因子,不同变量的取值差异很大。直方图反映一个变量的最小值和最大值及变量值的分布情况,在数据分块过程中可以生成体直方图。为了实现多维多变量的可视化,基于直方图自定义传输函数,系统可以选择不同的变量。对于不同时刻的数据,采用时间插值实现动态效果。图 10表示不同时刻的台风可视化效果图,从图中可以清晰地看到一个细小的、轮廓清晰的台风眼,T0表示台风开始模拟时刻,T5T11T17分别表示开始时刻之后5 h、11 h和17 h。

图 10 台风动态可视化效果图 Fig. 10 Effect Chart of Typhoon Dynamic Visualization Varying with Time
3.4 剖面分析

为了分析体数据的特征,系统还提供了一系列的分析功能,图 11表示一个台风的剖面分析图。系统提供沿着经度和纬度的方向做剖面的功能,并可以设置沿高度方向相对地面的高度,不同的颜色表示不同大小的风速。

图 11 台风剖面分析图 Fig. 11 Profile Analysis Diagram of Typhoon
4 结语

本文讨论了虚拟地球平台下台风多维可视化系统的设计和实现,重点论述了坐标转换、数据组织和基于GPU的体可视化三个关键的技术。在开源WW平台已有功能的基础上,扩展其台风动态可视化,并介绍了系统的体系结构和主要功能,通过模拟数据验证了本系统的有效性和实用性。和传统的台风可视化系统比较,系统充分利用了WW已有的功能,在网络环境下实现台风交互可视化,可视化效果更符合真实的自然效果。

虽然系统可以满足大规模台风数据在网络环境下的动态交互可视化,但系统还需要进一步优化。下一步的工作主要集中在以下两个方面:① 系统在可视化方面的功能需要进一步的开发,提高算法的效率;② 由于海洋和气象数据的相似性,系统将进一步扩展,满足海洋数据的需求。

参考文献
[1] Nocke T, Sterzel T, Böttinger M, et al. Visualization of Climate and Climate Change Data: An Overview[C]. Digital Earth Summit on Geoinformatics 2008: Tools for Global Change Research (ISDE'08), Wichmann, Heidelberg, 2008
[2] Van Thu T, Krishnamurti T. Vortex Initialization for Typhoon Track Prediction[J]. Meteorology and Atmospheric Physics, 1992, 47(2/4): 117-126
[3] Zhu Zhenyi. 3D Modeling and Secene Simulaiton of Typhoon[D]. Shanghai: Tongji University, 2008 (朱震毅. 台风三维建模与视景仿真研究[D]. 上海:同济大学, 2008)
[4] Nocke T, Flechsig M, Bohm U. Visual Exploration and Evaluation of Climate-related Simulation Data[C].2007 Winter Simulation Conference, Washington D C, USA, 2007
[5] Doleisch H. SimVis: Interactive Visual Analysis of Large and Time-dependent 3D Simulation Data[C].The 39th Conference on Winter Simulation: 40 Years! The Best is Yet to Come, Washington D C, USA, 2007
[6] Wang H, Lau K H, Chan W M, et al. A PC-Based Visualization System for Coastal Ocean and Atmospheric Modeling[C].Estuarine and Coastal Modeling (1999), New Orleans, Louisiana, USA, 2002
[7] Qin Xujia, Zhang Qinfeng, Chen Jian, et al. GPU Accelerated Typhoon Visualization Method[J]. Journal of Image and Graphic, 2012, 17(2): 293-300(秦绪佳, 张勤锋, 陈坚, 等. GPU加速的台风可视化方法[J]. 中国图像图形学报, 2012, 17(2): 293-300)
[8] Kim C S, Parks K, Park J, et al. Scientific Visualization of Time-varying Oceanographic and Meteorological Data Using VR[C].The 16th IEEE Visualization 2005 (VIS 2005), Minneapolis, MN, USA, 2005
[9] Xu Min, Fang Chaoyang, Zhu Qing, et al. Design and Implementation of Multidimension and Animated Visualization System for Ocean and Atmosphere[J]. Geomatics and Information Science of Wuhan University, 2009, 34(1): 57-59(徐敏, 方朝阳, 朱庆, 等. 海洋大气环境的多维动态可视化系统的设计与实现[J]. 武汉大学学报·信息科学版, 2009, 34(1): 57-59)
[10] Gore A. The Digital Earth: Understanding Our Planet in the 21st Century[J]. Australian Surveyor, 1998, 43(2): 89-91
[11] Blenkinsop T. Visualizing Structural Geology: From Excel to Google Earth[J]. Computers & Geosciences, 2012, 45:52-56
[12] Smith T M, Lakshmanan V. Real-time, Rapidly Updating Severe Weather Products for Virtual Globes[J]. Computers & Geosciences, 2011, 37(1): 3-12
[13] Nguyen Q C, Soon K T. Google Earth as a Tool in 2D Hydrodynamic Modeling[J]. Computers & Geosciences, 2011, 37(1): 38-46
[14] Ballagh L M, Raup B H, Duerr R E, et al. Representing Scientific Data Sets in KML: Methods and Challenges[J]. Computers & Geosciences, 2011, 37(1): 57-64
[15] Bailey J E, Chen A. The Role of Virtual Globes in Geoscience[J]. Computers & Geosciences, 2011, 37(1): 1-2
[16] Joseph T, Jeff H, Kim R, et al. A Tropical Cyclone Application for Virtual Globes[J]. Computers & Geosciences, 2011, 37(1): 13-24
[17] Li J, Wu H, Yang C, et al. Visualizing Dynamic Geosciences Phenomena Using an Octree-based View-dependent LOD Strategy Within Virtual Globes[J]. Computers & Geosciences, 2011, 37(9): 1 295-1 302
[18] Hogan P, Coughlan J. NASA World Wind, Open Source 4D Geospatial Visualization Platform:*. NET & Java*[C]. AGU Fall Meeting Abstracts, San Francisco, USA, 2006
[19] Snyder J P. Map Projections—A Working Manual[R]. United States Government Printing Office, Washington D C, 1987
[20] Evenden G I. Cartographic Projection Procedures for the UNIX Environment: A User's Manual[R]. Open-File Report 90284, United States Geological Survey, Washington D C, USA, 1990
[21] Guthe S, Wand M, Gonser J, et al. Interactive Rendering of Large Volume Data Sets[C].IEEE Visualization, 2002 VIS 2002, Boston, MA, USA, 2002
[22] Kniss J, Mccormick P, Mcpherson A, et al. Interactive Texture-based Volume Rendering for Large Data Sets[J]. IEEE Computer Graphics and Applications, 2001, 21(4): 52-61
[23] Boada I, Navazo I, Scopigno R. Multiresolution Volume Visualization with a Texture-based Octree[J]. The Visual Computer, 2001, 17(3): 185-197
[24] Lamar E, Hamann B, Joy K I. Multiresolution Techniques for Interactive Texture-based Volume Visualization[C].The Conference on Visualization'99: Celebrating Ten Years, San Francisco, California, United States, 1999
[25] Yang C, Wu L. GPU-Based Volume Rendering for 3D Electromagnetic Environment on Virtual Globe[J]. International Journal of Image, Graphics and Signal Processing, 2010, 2(1): 53-60
[26] Li W, Mueller K, Kaufman A. Empty Space Skipping and Occlusion Clipping for Texture-based Volume Rendering[C]. The 14th IEEE Visualization 2003 (VIS 2003), Seattle, Washington D C, USA, 2003
[27] Kruger J, Westermann R. Acceleration Techniques for GPU-based Volume Rendering[C]. The 14th IEEE Visualization 2003 (VIS 2003), Seattle, Washington D C, USA, 2003
[28] Konikow L F, Bredehoeft J D. Computer Model of Two-dimensional Solute Transport and Dispersion in Ground Water[R]. US Government Printing Office, Washington D C, USA, 1978
[29] Butcher U J C. The Numerical Analysis of Ordinary Differential Equations: Runge-Kutta and General Linear Methods[M]. New York:Wiley-Interscience, 1987