2. 地球空间信息技术协同创新中心, 湖北 武汉 430079
2. Collaborative Innovation Center of Geospatial Technology, Wuhan 430079, China
海洋覆盖了全球近70%的表面，是虚拟地球表达的重要内容之一。海面可视化即通过对海面进行网格化处理，并使网格变成高低起伏的海浪高度场，同时在每一帧绘制时按一定规律对网格点位置进行更新，进而实现海浪的动态效果。与传统海面可视化相比，虚拟地球中海洋表面可视化，能够突破传统海面可视化在空间尺度上的限制，实现全球多尺度海面的无缝浏览，因此是目前海面可视化系统中最为流行的一种模式，广泛应用于舰船导航、航线规划、渔船监控、虚拟战场等领域。
现有的虚拟地球中的海面动态可视化方法主要采用投影网格[1-2]的方式实现全球海面动态效果，如商业公司Sundog开发的Triton海洋可视化中间件 (WGS-84模式) 即采用这类方法。该方法通过在投影空间构建海面格网，并根据地球形状确定格网的覆盖范围，由于投影空间构建的格网能够跟随视点位置进行自动缩放平移，因此在不构建新的格网情况下实现了海域的无限延伸及尺度缩放。然而这种特殊的网格组织方式无法与虚拟地球离散的空间剖分框架相融合，且所有地区采用同一格网表示，在反映不同海域的差异化特征时存在缺陷，尤其体现在基于投影网格法的精细海陆分界上。
基于全球离散格网的海面可视化方法[3-4]弥补了以上缺陷，它主要基于瓦片四叉树[3]的方式组织海面格网，因此能够与虚拟地球中的多尺度空间数据相匹配，除了根据海底地形高度与海平面高度区分陆地海洋外，对于缺少海底地形数据的地区可以调用相应位置海域分布数据 (如从电子海图中提取的海岸线数据等) 进行更加精细海陆分界。然而与投影网格法相比，现有的离散格网法也存在较为明显的缺陷：首先，全球离散格网法在场景浏览时需要频繁进行海面格网的动态加载及卸载，尽管可以在后台线程中进行，但如果格网构建复杂或内存占用量较大时，仍然会对整体效率造成影响；其次，全球离散格网存在因不同尺度格网拼接导致的格网缝隙问题，而现有的缝隙修补算法[4-6]普遍计算复杂，严重影响场景绘制及加载效率。
综上所述，基于全球离散格网的海面动态可视化方法能够弥补主流方法—投影网格法在反映不同海域差异性特征时存在的缺陷，实现更为精细海陆分界效果。同时，其格网组织方式与虚拟地球离散空间剖分框架相匹配，因此更符合虚拟地球的应用需要。然而现有的全球离散格网法存在绘制效率低、加载速度慢、需要进行格网缝隙修补等问题。这些都限制了全球离散格网法的应用。对此，本文在现有离散格网法基础上进行优化：首先，在数据结构上，本文对传统等经纬度离散格网进行扩展，设计一种面向GPU的多尺度海面网格模型来组织管理海面格网；其次，为了实现风场驱动下的海浪动态绘制，本文在多尺度海面网格模型的基础上提出了一种支持实时风场更新的海浪动态绘制方法；同时，考虑到格网缝隙修补对绘制效率的影响，本文针对海面格网特点并结合GPU技术提出了一种高效的海面格网缝隙修补方法。
1 面向GPU的全球多尺度海面网格模型为了对全球海面格网进行有效组织和管理，同时充分利用GPU并行计算的能力提高海面绘制效率，本文在全球等经纬度离散格网模型[7-8]的基础上进行扩展，设计了一种面向GPU绘制的全球多尺度海面网格模型，如图 1所示。
|
| 图 1 面向GPU的多尺度海面网格模型 Fig. 1 GPU-oriented multi-scale grid model of ocean surface |
海面网格是海面数据组织及绘制的基本单元，按照全球等经纬度网格剖分的规则进行组织与构建。由此，每个海面格网都有唯一的空间编码[8]与全球离散格网相对应，它记录了当前格网所在的级数、行列号。通过空间编码可以检索到当前海域的风速、风向和海水深度、海域分布等信息以及确定网格之间的邻接关系。
海面网格的几何对象描述了海面的几何结构信息，由顶点和顶点索引组成。顶点即格网点，用于表达海面高低起伏的最小单元，除基本的位置信息外，为了实现基于GPU的海面动态可视化，对顶点属性进行扩展，添加了包括经纬度、顶点类型在内的扩展属性。顶点索引记录了离散顶点之间的组织关系，用于实现点向面图元的转换。
为了实现海浪的GPU动态绘制，将海面网格的纹理对象分为位移纹理和表面纹理，位移纹理记录了海面波浪的动态起伏过程，以当前海域的风向、风速为参数再基于海浪建模算法实时生成，最后以三维纹理的方式载入到GPU中。表面纹理则记录了当前海面格网的颜色及透明度信息，可通过海水深度或记录海域分布的掩模数据生成，其中颜色的深浅可以反映因为海水的深度导致的水体明暗效果。而透明度除了反映水体的质感外，还可用于隐去陆表水体，以达到海陆分界效果。
这种模型结构的优点在于将海面格网的几何结构与其具体的形态特征 (记录在位移纹理中) 相互分离，仅在绘制时通过GPU进行关联。海面格网初次加载时只需构建最基本的几何结构，而较为复杂形态特征可通过资源共享的方式被多个海面格网重复利用，因此有效提高了格网的构建速度。
图 2为基于多尺度海面网格模型实现的海面多尺度绘制流程。首先对虚拟地球中的海域进行等经纬度划分，随着视点的变换，动态调度可视海域内相应尺度的空间数据，根据空间数据构建海面网格的几何与纹理。对于同一海域不同尺度的海面格网则以自顶向下四叉树金字塔的方式进行组织，随着视点的递进格网的分辨率逐倍增加，海面波浪的起伏效果也逐渐明显。
|
| 图 2 全球海面多尺度绘制流程 Fig. 2 Drawing Process of Multiscale Global Ocean Surface |
2 风场驱动下的海浪动态绘制方法
海面可视化的真实感绝大程度上取决于海面浪花的动态起伏效果，尤其是风场条件驱动下的海面浪花动态绘制，不仅在视觉上满足用户的需要，同时能够与真实地理环境相匹配，更符合地理信息仿真应用的要求。对此，本文在多尺度海面网格模型基础上，提出了一种风场驱动下的海浪动态绘制方法。
为了实现海面波浪的动态绘制，首先需要对海面波浪的运动规律进行归纳建模，本文采用了文献[9]提出的一种基于海浪谱的海浪建模方法。该方法认为海面波浪高度h(x, t) 是由一系列正弦、余弦波叠加而成，因而可通过二维快速傅里叶逆变换的方式求解, 其公式如下
(1)
式中，m=(x, y) 表示海面点的水平位置；
然而文献[9]这种以高度场序列的形式描述海浪运动规律，仅能对海浪在垂直方向的起伏特征进行模拟，如图 3(a)所示。在真实地理环境中，海浪往往会在风场作用下发生倾斜，如图 3(b)所示。其中z方向表示海面的垂直方向，x表示水平方向，箭头表示风场。
|
| 图 3 水平方向的海浪剖面 Fig. 3 Horizontal Profile of Ocean Wave |
本文对文献[9]的方法进行了改进，根据海浪当前的起伏状态及当前海面风场状况计算每个海浪点在水平方向可能发生的偏移值，从而实现更为逼真的海浪模拟。
主要步骤如下：
首先对式 (1) 求梯度
式中，Gxm, t、Gym, t分别表示高度场横纵方向的梯度值，同样基于二维快速傅里叶逆变化求解，由于梯度指向高度场增长最快的方向，对于海面高度场来说，其梯度矢量的单位向量 (gx, gy) 可看作海浪发生倾斜程度最大的方向。则海浪在水平方向的偏移值为
式中，f为梯度方向 (gx, gy) 与当前风场矢量 (wx, wy) 点积，反映了风场对海浪的影响程度。
由此，建立了以垂直、横、纵3个方向的偏移值描述当前时刻海浪运动的海浪模型，通过将3个方向偏移值作用于海面格网点的初始位置之上即可实现海面的动态绘制。在动态绘制方式上，常规基于CPU实时修改格网点位置的方式，计算量较大，影响绘制效率。而预先构建多时态格网的方式会增加内存的占用量及影响网格的加载效率，不适合大规模的海面动态可视化，且格网构建后修改困难，无法与实时风场进行匹配。对此，本文在多尺度海面网格模型的基础上提出了一种基于GPU的海面动态绘制方法，通过GPU可以并行地对海面格网点位置进行更新，减少CPU计算量，同时将海浪运动的偏移值以位移纹理的方式载入到GPU中供着色器采样，当风场条件发生变化，只需重新载入位移纹理，不需要对构建好的海面格网进行任何修改。其详细算法流程如图 4所示。
|
| 图 4 海面动态绘制流程 Fig. 4 Dynamic drawing process of ocean surface |
(1) 首先遍历虚拟地球瓦片四叉树，确定当前可视海域范围，并构建当前海域的海面格网。
(2) 根据格网的空间编码获取视点最近的海面格网的风向、风速。取其平均值为参数，根据海浪建模算法计算出海浪在横、纵、垂直3个方向的偏移值，其值分别存储在位移纹理的r、g、b 3个通道内。最后将不同时刻的位移纹理进行合并，得到立体纹理后载入顶点着色器中，如果当前场景的风速、风向与上一帧变化不大则不需要进行位移纹理的更新，以减少计算量。
(3) 在顶点着色器中，根据网格顶点的经纬度坐标 (lat，lon) 及当前帧数 (frame) 确定纹理坐标
式中，hres为海浪高度场的经纬度跨度，值为m表示海面高度场对应m°×m°的地球表面区域；totalframe表示一次海浪运动循环的总帧数；%表示取余。
(4) 通过获取的纹理坐标采样当前海域的立体纹理，得到顶点的偏移值D，并将偏移值叠加在原有顶点位置之上
(2)
式中，Dx=D.r·Nx；Dy=D.g·Ny；Dz=D.b·gl_normal；Nx、Ny为顶点所在位置切地球表面时的横纵向量；gl_position为最终输出的顶点坐标；gl_Model View Projection Matrix为顶点变换矩阵gl_vertex为原始顶点坐标；gl_normal为顶点法向量。
3 基于GPU的海面格网缝隙修补方法海浪动态绘制中相邻的不同尺度的海面格网因为分辨率不同，在拼接时，边缘线无法重合而产生缝隙，影响整体视觉效果，如图 5(a)所示。对此，文献[4]预先构建适应不同邻接关系的边缘格网，在绘制时根据需要动态的进行替换，然而其方法需要对相邻格网级差进行约束以减少预构的边缘格网数量。文献[5]在相邻格网之间根据顶点位置实时生成一系列三角形条带进而达到修补缝隙的目的，但其算法复杂，计算量较大，实现较为困难。文献[6]则利用最新几何着色器技术动态的修改边缘几何体的分形方式，从而避免裂缝的产生，但随着相邻格网级差的增加，分形方式会变得十分复杂，进而影响绘制效率。
|
| 图 5 海面裂缝修补 Fig. 5 Repair the crack of ocean grid |
综上所述，现有的海面格网缝隙修补算法，主要是在传统地形格网缝隙修补算法基础上的改进，容错性高，但计算也较为复杂，对格网的绘制效率及加载速度会造成较大的影响。考虑到海面格网相对于地形格网，其几何形态灵活可变，限制条件较少 (全球海平面高度可视为定值，波浪起伏高度不受具体值限制)。因此本文针对海面格网这一特点提出了一种以GPU技术为核心的高效海面格网缝隙修补算法。
本文实现的海面缝隙修补基本思路如图 5(a)和(b)所示，首先确定相邻格网边缘上的重合顶点 (A点)，对重合顶点之间的顶点 (B点，在高级别格网的边缘线上) 通过顶点着色器重新计算其空间坐标。使其刚好落在重合顶点的连线上。其详细算法流程如下。
(1) 场景发生变换时，更新格网的邻接关系，即记录下每个格网的相邻格网的空间编码 (根据空间编码可得到格网的空间范围等信息)。
(2) 在顶点着色器中，首先根据顶点的类型判断是否需要进行坐标修改。如图 5(c)所示，海面格网中的每个顶点根据其位置的不同可分为上边缘点 (a点)、下边缘点 (b点)、左边缘点 (c点)、右边缘点 (d点)、角点 (e点)、内部点 (f点)6种。如果顶点属于内部点则不需要进行修改。如果顶点属于边缘点，但与其相邻的格网级别较高或相等，其顶点坐标同样不需要修改。对于角点，将其当作与相邻格网 (不考虑对角格网) 中级别最大的格网相邻边缘的边缘点处理。
(3) 获取离待修改边缘点最近的左右重叠点的经纬度坐标 (lon1，lat1)(lon2，lat2)，其公式如下
式中，假定每个格网具有m×m个格网点，m为奇数。(lon，lat) 为待修改边缘点的经纬度坐标，(lon0，lat0) 为相邻格网的边缘线上最小经纬度坐标；res0为相邻格网的分辨率 (单个格边的经纬度跨度)；INT表示取整数。
(4) 根据重叠点的经纬度坐标 (lon1，lat1)(lon2，lat2) 及当前时刻 (帧数) 采样位移纹理，分别得到重叠点的位移值D1、D2。根据线性比例计算边缘点的偏移值
对于上下边缘点，a=lon、b=lon1、c=lon2。对于左右边缘点，a=lat、b=lat1、c=lat2。
(5) 将位移值D导入式 (2) 中得到修改后的边缘点坐标。
4 试验为了验证本文提出方法的有效性和可行性，基于开源虚拟地球平台osgEarth进行试验，软件环境为Window 7，Visual Studio 2010和OpenGL。硬件配置为Intel (R) Core (TM) i3-2100双核3.1 GHZ CPU，NVIDIA Quadro 600显卡，1 GB显存8 GB内存。
试验数据中影像采用ESRI发布的ESRI_Imagery_World_2D地图服务。地形采用ReadMap发布的地形瓦片服务，其陆地部分为全球90 m分辨率 (3″×3″) 的STRM数据，海底部分则采用GEBCO (General Bathymetric Chart of the Oceans) 提供全球的深海测量数据 (30″×30″)。风场数据则根据舰船传感器实时获取的风速、风向信息进行模拟插值生成。海域分布数据 (海水与地表分别以不同颜色表示的图像) 则通过电子海图中的海岸线数据编辑生成。试验前根据多尺度海面网格模型，对数据进行多尺度组织，构建四叉树结构的瓦片金字塔。
|
| 图 6 本文方法海面全球视图 Fig. 6 Global view of ocean (the suggested method) |
|
| 图 7 本文方法海面局部视图 Fig. 7 Local view of ocean (the suggested method) |
|
| 图 8 本文方法海面近景视图 Fig. 8 Near view of ocean (the suggested method) |
对比3副图片可以发现，随着视角的递进，海洋表面逐渐清晰，海水色调过渡平滑，宏观海陆分界效果及微观海面波浪效果都能得到有效体现，基本满足全球多尺度海面无缝浏览的应用需要。
图 8和图 9为本文方法与投影网格法实现的海面近景视图的比较。
|
| 图 9 投影网格法海面近景视图 Fig. 9 Near view of ocean (projected grid method) |
从图 8和图 9中可以看出，本文方法在离散格网框架的支持下，能够快速地调用当前区域的海域分布数据，实现精细海陆分界。而投影网格法目前仅能通过RTT技术实时获取的地表高程实现大致的海陆分界效果 (Trition提供的海陆分界方案)，容易出现海水覆盖真实地表等问题。
图 10和图 11为缝隙修补前后海面效果比较。如图 10所示，不同尺度海面格网拼接产生的“T型”缝隙，在海面波浪动态绘制时，变得十分明显，严重影响视觉效果。而采用本文提出缝隙修补方法，可以在不占用过多系统资源情况下，实现海面缝隙快速修补。
|
| 图 10 缝隙修补前的海面 Fig. 10 Ocean surface before repair crevices |
|
| 图 11 缝隙修补后的海面 Fig. 11 Ocean surface after repair crevices |
图 12、图 13为本文方法在不同风场条件下实现的海面波浪效果比较 (图 12风速为5 m/s微风，图 13为24 m/s和风)。风速较高时海面波浪起伏效果更为明显。
|
| 图 12 风速为5 m/s海面波浪效果 Fig. 12 Ocean wave drived by wind (5 m/s) |
|
| 图 13 风速为24 m/s海面波浪效果 Fig. 13 Ocean wave drived by wind (24 m/s) |
为了验证本文方法的绘制及加载效率，在相同试验条件下，将本文方法与基于投影网格的海面可视化方法及现有的基于离散格网的海面可视化方法 (基于文献[3]及文献[4]实现的算法) 进行了比较，试验结果如图 14、图 15所示。
|
| 图 14 3种方法帧率比较 Fig. 14 Frame rate comparison |
|
| 图 15 3种方法加载时间比较 Fig. 15 Loading time comparison |
从图 14可以看出商业软件Triton实现的基于投影网格的海面可视化方法在绘制效率方面具有较为明显的优势，主要原因在于基于离散格网的海面可视化方法，在每一帧绘制时都需要进行海面缝隙修补，因此对帧率有一定影响。而本文实现的改进的离散格网法。由于采用GPU技术进行缝隙修补，在效率上较现有的离散格网法有很大提高，且与投影网格法相比帧率差距能够保持在6帧/秒以内。
图 15可以看出本文方法较其他两种方法加载速度更快。而现有的离散格网法由于需要实时限制相邻格网的级差以减少缝隙修补的复杂性，因此随着格网数量的增加，加载时间呈指数上升的趋势。而投影网格法尽管不需要重新构建新的格网，但其海陆分界方式存在效率问题，最终影响整体的加载速度。相比之下，本文方法采用瓦片四叉树的方式组织管理格网，同级格网之间没有约束、相互独立，因此可以充分利用多线程进行快速地并行加载，90万个网格顶点，开启2个线程的情况下，基本能在5 s左右完成加载。
为了验证本文海面可视化方法具有能够快速反映实时风场的能力，对3种方法在风场条件发生变化时，海浪效果的更新时间进行了比较。如图 16所示。
|
| 图 16 3种方法更新时间比较 Fig. 16 Update time comparison |
从图 16可以看出，由于本文方法在风场条件发生变化时，只需修改位移纹理即可对整个海浪效果进行更新，因此更新时间能够稳定在0.5 s左右，且该部分计算可在后台线程中进行，不会影响当前正在绘制的海浪。而其他两种方法，海浪的更新相当于重新加载一次海面格网，从更新时间及视觉体验上都无法满足实时动态风场更新的需要。
5 结语针对现有的基于全球离散格网的海面动态可视化方法存在绘制效率差，加载速度慢，需要进行格网缝隙修补等问题，本文提出了一种优化方法，主要包含3部分内容：首先，为了能够对全球海面网格进行有效组织和管理，同时充分利用GPU并行计算能力，本文在传统等经纬度离散格网基础上进行扩展，设计了一种面向GPU绘制的全球多尺度海面网格模型；其次，考虑到风场驱动下的海浪起伏效果更符合地理信息仿真应用的需要，本文在多尺度海面网格模型的基础上提出了一种支持风场条件实时更新的海浪动态绘制方法；同时，由于现有的海面格网缝隙修补算法会严重影响海面的绘制及加载效率，本文针对海面格网特点，提出了一种基于GPU的高效海面缝隙修补算法。最后通过试验对本文方法进行验证。试验结果表明，本文方法绘制效率稳定、加载速度快，且具有反映不同海域的差异性特征、实现精细的海陆分界及支持风场条件动态更新等功能，因此更符合虚拟地球的应用需要。下一步将讨论全球海洋表面如何与舰船航行等进行实时交互动态绘制，为面向全球的海洋具体应用奠定基础。
| [1] | JOHANSON C. Real-time Water Rendering:Introducing the Projected Grid Concept[D]. Lund:Lund University, 2004. |
| [2] | 赵欣, 裴炳南. 一种快速的海浪仿真方法[J]. 系统仿真学报, 2012, 24(1): 132–135, 145. ZHAO Xin, PEI Bingnan. Fast Ocean Wave Simulation Method[J]. Journal of System Simulation, 2012, 24(1): 132–135, 145. |
| [3] | 明德烈, 徐秋程, 李向春. 面向全球应用的海洋仿真系统的实现研究[J]. 系统仿真学报, 2012, 24(8): 1741–1745, 1750. MING Delie, XU Qiucheng, LI Xiangchun. On Realization of Visualization System for Global Ocean Simulation[J]. Journal of System Simulation, 2012, 24(8): 1741–1745, 1750. |
| [4] | YANG Xudong, PI Xuexian, ZENG Liang, et al. GPU-Based Real-Time Simulation and Rendering of Unbounded Ocean Surface[C]//Proceedings of the Ninth International Conference on Computer Aided Design and Computer Graphics. Hong Kong, China:IEEE, 2005:6. |
| [5] | LI Bo, WANG Cheng, LI Zhongyi, et al. A Practical Method for Real-Time Ocean Simulation[C]//Proceedings of the 4th International Conference on Computer Science & Education. Nanning:IEEE, 2009:742-747. |
| [6] | PUIG-CENTELLES A, RAMOS F, RIPOLLES O, et al. View-Dependent Tessellation and Simulation of Ocean Surfaces[J]. The Scientific World Journal, 2014, 2014: 979418. |
| [7] | 龚健雅, 陈静, 向隆刚, 等. 开放式虚拟地球集成共享平台GeoGlobe[J]. 测绘学报, 2010, 39(6): 551–553. GONG Jianya, CHEN Jing, XIANG Longgang, et al. GeoGlobe:Geo-spatial Information Sharing Platform as Open Virtual Earth[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(6): 551–553. |
| [8] | 童晓冲. 空间信息剖分组织的全球离散格网理论与方法[J]. 测绘学报, 2011, 40(4): 536. TONG Xiaochong. The Principles and Methods of Discrete Global Grid Systems For Geospatial Information Subdivision Organization[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(4): 536. |
| [9] | TESSENDORF J. Simulating Ocean Water[C]//Simulating Nature:Realistic and Interactive Techniques, SIGGRAPH2001 Course Notes. Los Angeles, CA:Addison Wesley, 2001:47-58. |