受硬件纹理空间、内存容量、计算性能等限制，动态可视化时不可能将TB级海量气象场数据一次性从外存装入。因此面向三维气象场的数据模型需要将复杂的场数据离散化成多层次的数据单元，以实现在视点调度时可视域内的数据快速装载，从而提高动态可视化效率。为此，本文提出了基于体元对象的数据模型(图 1)。

图  1  基于体元对象的数据模型

Figure 1.  Voxel Object-oriented Data Model

该模型的核心是构建多级体元包围时空上离散动态分布的三维场数据点，用体元的空间大小和变化频率描述三维场的时空尺度。体元是连续时间变化下一定的空间范围，是三维场数据的存储、管理和操作单元。体元包括三维场值、空间立体块以及时间点信息三部分。三维场值包括基本的气象场数据及其空间点位数据；空间立体块对应虚拟地球中相邻尺度的等经纬度格网构成的空间范围，依照八叉树全球等经纬度立体剖分方法^[11]生成，并基于椭球八叉树进行分层分块组织；时间点信息按不同时间跨度进行尺度划分，同一尺度时间间隔相等，上下层互为包含关系，构建尺度从大到小的R⁺树结构。空间立体块和时间点信息都是以编码形式存储为体元对象的属性字段，二者的键值共同构成体元的惟一标识。

在上述数据模型的基础上，本文基于Q⁺R树索引和空间信息多级格网索引^[12]设计了面向三维气象场的多级体元时空索引。如图 2所示，多级体元时空索引是一个先时间索引后空间索引的二级索引结构。图 2中，T表示R树的节点(时间)；A、B、C、D表示八叉树节点(空间体元)，他们是大的时间尺度下的粗级体元，时间树带动空间树更新。如果某个时刻时间尺度发生了变化，那么时间节点T1重建生成了细的时间节点，对应的空间树发生变化，变成D₁，D₂，D₃，D₄甚至更多，D₁~D₄这些树在同一视点下空间层级与D是一致的，只是用于表示变化更快的时间尺度下的场数据。

图  2  多级体元时空索引检索机制

Figure 2.  Voxel Based Multi-level Indexing Mechanism

在进行检索时，首先判断时间尺度是否变化，如果发生变化，需要先对索引树的时间结点进行重建。若无变化，则根据时间点索引，并依照视点调度和椭球八叉树去定位最粗一级的体元，然后按空间索引向子层检索，直到满足所需的空间尺度为止。

矢量场可视化方法包括图标法、流线法、流面法、流管法等，其核心^[13]是将矢量数据映射成通用的、可通过图形直观显示的几何数据。其中，图标法是一种典型的逐数据点模拟矢量场的方法，快速高效，能保留数据所有细节信息并且很好地承载场的方向信息。由于传统二维箭头和流线在三维空间中不同视角下存在二义性，本文以能自由旋转的三维箭头模型作为映射方式，长度代表场值大小，箭头指向场的三维方向。

矢量场可视化的关键是场模型的方向在空天空间里应始终平行于该经纬点处球体切线方向。因此，每当时相发生变化时，就需要对当前的场模型的姿态、长度进行调整。场模型的姿态调整是空间坐标系(X_model，Y_model，Z_model)与世界坐标系(X_E，Y_E，Z_E)两个坐标系之间的变换过程。如图 3所示，变换矩阵的求解分为三步。首先将场模型从合方向校正为正北方向，将场模型的Z轴变换为世界坐标Z轴，最后将场模型的X、Y轴校正为世界坐标的X、Y轴。

图  3  场模型姿态自动校正

Figure 3.  Posture Automatic Adjustment of Field Model

基于体绘制的标量场可视化方法包括光线投射法、抛雪球法、错切—变形法和三维纹理映射法等。其中三维纹理映射法由于具备硬件支持及对密集体数据的高效绘制能力而成为当今较实用的体绘制方法^[14]。

基于三维纹理映射体绘制的三维标量场可视化大致流程^[15]如下。首先，从硬盘中把体元数据装载到内存中；然后，通过对体元数据分类和图像特征映射生成纹理缓存，再装入图形硬件的纹理内存生成三维纹理；最后，随着系统视点调动对体元进行实时纹理面重采样，采样面与三维纹理合成结果图像。

该流程中最重要的是纹理面的重采样。本文采用的是图像空间纹理重采样法，采样面垂直于视线方向并等间距分布在体元内。为了方便采样面方程的求解，本文设定初始采样面垂直视线方向并经过体元中心点，其他采样面以初始采样面为基准沿视线和逆视线方向按等间距增量生成，采样面与体元的交点个数N存在3、4、5、6共4种不同情况，如图 4所示。图 4中d₀为采样面之间的间距。

图  4  采样面等距增量生成

Figure 4.  Sampling Surfaces Generate Equidistantly

采样面与体元交点的坐标是由采样面与体元边的方程求解得到的。本文规定交点在体元点上或体元边外视为无交点。在不影响最终结果的前提下，该方法排除了极少量的正常采样面，但也避免采样面为一个点、一条线或体元面的错误情况，在上述条件下，当采样面与某一条体边相交时，其必与该边所在的左右两个体面的边相交。文献^[16]根据这一特性，将如图 5(b)所示的邻边结构作为体边遍历顺序。图 5(a)中数字代表体元点编号，带下划线数字代表体元边编号，图 5(b)中的粗线代表体边0的邻边。规定邻边的遍历为先右平面后左平面，且在左右平面中平行边优先于起始点和终止点的邻接边。

由于该算法存在盲目遍历和重复遍历的不足，本文在其基础上作进一步改进。

改进一 如图 5(c)所示，由于采样面是沿视线方向及其逆方向生成的，假设有一个垂直于视线方向的采样面K随视线方向移动，则按其经过的体元顶点顺序可以推知体边遍历优先顺序。设体元顶点P_j(j=0,1,2，…，7)与体元中心点O组成向量$\overrightarrow{{{P}_{j}}O}$。首先计算$\overrightarrow{{{P}_{j}}O}$与视线向量的点积，点积值越小，表明该顶点越先被采样面K经过，该顶点的邻接体边也就越早被遍历。然后把点积最小和最大的体元顶点的邻接体边分别设为优先度最大和最小，其余顶点的邻接体边也按顶点点积值由小到大依次设置优先度，但其与最大最小优先度的公共体边的优先度保持不变。最后将所有体边按优先度从大到小排列，同一优先度的多条体边按编号大小排列，构成体边遍历优先序列。这样，当采样方向为逆方向时，只需将优先序列倒过来使用即可。在图 5(c)中，当前视线方向下的优先序列为红色体边(2、3、11)＞蓝色体边(6、7)＞黄色体边(8、10)＞紫色体边(0、1)＞绿色体边(4、5、9)。

图  5  体元拓扑结构

Figure 5.  Voxel Topological Structure

改进二 由于不同体边的邻边之间存在很大的交集，会导致有的体边被重复遍历。对此，可以定义一个全局变量用于记录遍历过的体边编号，每次进行求交点算法前，对体边是否遍历过进行判断。另外，该全局变量应排除记录首个与采样面相交的体边，因为该边是判断遍历结束的重要条件，需要多次遍历。

采样的过程中，首先遍历体边找到第一个交点(初始交点)，然后对该交点所在体边的邻边进行遍历求交找到下一个交点，继续遍历新交点的邻边寻找新交点……如此循环，直到求得的新交点等于初始交点为止，此时采样多边形已经闭合，可继续下一个采样面。如果没有交点存在则说明采样已经结束，退出算法。

由于体元的顶点的空间坐标(存储在体元数据缓存中)和纹理坐标(存储在三维纹理内存中)已知，交点的空间、纹理坐标是由体元顶点的坐标线性插值得到的。交点存在情况判断方法及其空间坐标、纹理坐标的求解如图 6所示。

图  6  交点判断及交点空间坐标、纹理坐标求解

Figure 6.  Spatial and Texture Coordinates Calculation of Intersection Points

在图 6中，P₀为当前所求采样面上的点，其位置由体元中心坐标加上其在视线方向上的相对偏移量求得，P₁P₂表示体元空间当前需要判断是否相交的边，Dir为当前视线方向(已单位化的)，P_x为交点，P_j1、P_j2为P₀和P₂在视线方向上的投影。设D等于向量u与Dir的点乘，N为向量v与Dir的点乘，T=N/D，则存在如下几种情况(δ为极小值)。

步骤1：如果D的绝对值小于δ，则此时采样面与线段P₁P₂平行或重合，执行步骤2。否则执行步骤3。

步骤2：如果N的绝对值小于δ，则体元顶点为交点或线段P₁P₂在采样面内。如果N的绝对值大于或者等于δ，则采样面与线段P₁P₂无交点。

步骤3：如果T＞1 或者 T＜0，则交点P_x存在，但不在线段P₁P₂内。否则交点P_x存在且在线段P₁P₂内。

因此，由相似三角形的关系得：

因此，P_x的坐标(pos)、纹理坐标(tex)为：

为了对上述方法进行验证，本文基于开放式虚拟地球集成平台GeoGlobe，对海量气象场数据的可视化效率进行了实验。实验采用的三维气象场数据是由美国国家海洋和大气管理局(NOAA)提供的NCEP/NCAR再分析数据中的1948年至今的月平均，2009~2013年全年日平均和一日4次的等压面通量数据。数据存储在全球144×73个网格点，共17层等压面上。其中矢量场使用的是约16 GB的风场数据，标量场使用的是约7 GB的比湿数据。实验是在Visual Studio 2010环境下用C⁺⁺和Direct 3D编写的。硬件配置为Intel(R)Core(TM)i3-2100双核3.1 GHz CPU，NVIDIA Quadro 600显卡，1 GB显存、8 GB内存。

此外，本文针对矢量场各异性、标量场聚类性的特点分别设计了不同的时空查询方法。

矢量场由于具有方向性，每个场模型代表独一无二的场状态，因此时空查询需要定位到每一个场模型。由于虚拟地球中一次绘制的场模型高达几千个，为了能快速定位到目标，本文采取先定位体元对象，再进行模型拾取的策略。通过视线与球面的交点反算出体元的行列号，再与当前层次、时间点信息构造键值去定位体元，然后对体元内的场模型进行拾取。

标量场可视化时对所有采样点的场值进行了统计分类，并分别赋予了不同的光学属性(颜色、不透明度)，因此标量场时空查询到的是范围值，其结果以可视化图像和统计图的形式显示。本文设计的标量场时空查询包括三种，即根据场值查询其空间分布，根据空间范围查询场值分布，以及根据时相变化查询场值分布变化。

三维场绘制效果和时空检索效果如图 7~10所示。

图  7  矢量场可视化效果和查询效果(黄色模型)

Figure 7.  Visualization and Query Result(Yellow Model) of Vector Field

图  8  标量场的可视化效果及场值分布查询

Figure 8.  Visualization and Distribution Query of Scalar Field

图  9  标量场中指定场值的空间分布查询

Figure 9.  Spatial Distribution Query of Specific Scalar Data

图  10  时相变化下查询场值分布的变化

Figure 10.  Distribution Query of Scalar Field Undertime-varying condition

本文设定时相变化的最小时间尺度为2 s。矢量场在17层约5 000个场模型同时进行时态变换的情况下，动态可视化帧速能达到33~34帧/s。标量场在保证较佳的图像质量(约150~200个纹理采样面)下，动态可视化帧速达到36~38帧/s，漫游流畅度较佳，证明了快速可视化算法的有效性。另外，矢量场标量场的时空检索时间大约为15~16 ms，证明了基于体元时空模型的高效性。此外，本文针对标量场动态可视化中提出了采样面计算改进算法，还与原算法进行效率对比实验，具体见表 1。

由表 1可知，改进后的三维纹理映射体绘制算法减少了约32%的采样面交点计算，提高了约18%的绘制效率。

本文主要从三维气象场的时空特征出发，在全球多尺度空间数据模型基础上拓展时间尺度，形成基于体元对象的数据模型并构建时空一体化的多级索引机制，然后基于这种数据组织方式和索引方式实现了海量三维气象场数据的流畅动态可视化和高效时空检索，并通过实验对方法的有效性进行了验证。下一步研究工作主要包括GPU加速以实现对更高密度场数据的绘制，设计多特征传输函数以实现更好的体绘制效果，以及引入时空分析算法实现三维气象场的时空分析等。