水下三维声场仿真与可视化方法研究

孙雪海，笪良龙，李玉阳

(海军潜艇学院 海军水下作战环境研究所，山东 青岛266071)

0 引 言

海洋水声环境是海军作战的主要环境，其涵盖了大量抽象、复杂的信息，如何快速地从如此庞大繁杂的数据中获取对作战有用的信息，并以直观、易于理解的方式呈现，对帮助作战人员感知环境规律和做出正确的指挥决策具有至关重要的作用。图形、图像等视觉信息具有很强的直观性，有利于加深人们对事物、规律的感知和理解。声作为目前用于水下目标探测的最有效手段，水下声场仿真及其可视化是复杂水声环境研究的关键技术。

水下三维声场数据是通过三维声传播模型计算所得。但三维声传播模型的计算量大、计算时间长，给复杂水声战场环境中的快速声场分析带来了严峻的挑战。随着高性能计算技术的发展，基于集群系统的并行计算为三维声场的快速计算提供了可能。文献［1］研究了三维射线－简正波－抛物方程模型在高性能计算平台的并行计算问题。文献［2］实现了水平不变海洋声道中WKBZ 简正波方法的并行计算。

水下三维声场可视化对帮助作战人员感知和理解水下声传播规律具有重要作用。文献［3 －4］研究了基于预积分体渲染技术的水声环境体可视化方法;文献［5 －6］利用基于三维纹理映射直接体绘制算法实现了水下三维声场的体可视化。

本文在前人研究的基础上，主要研究三维抛物方程模型N×2D 弱三维近似方法的并行计算策略及高质量水下三维声场的实时可视化问题。

1 三维声场计算

目前相对成熟的三维声传播模型主要有三维射线模型 (HARPO)、三维耦合简正波模型(CMM3D)和三维抛物方程模型(FOR3D)等［5］。抛物方程方法较射线方法计算速度慢，但在计算精度上高于射线方法;在保证相同精度的前提下，抛物方程方法的计算速度较简正波方法快得多。在求解低频声波随距离变化的声传播问题上抛物方程方法具有独特的优越性。

在柱坐标系(r，φ，z)下，简谐点源在水平变化声道中声场的三维Helmholtz 方程表示为［7］:

考虑不同方位φ 之间耦合关系后的FOR3D 模型，求解过程复杂、计算量大，不利于进行并行计算。但在一般情况下可以忽略方位之间的耦合，采用N ×2D 弱三维近似计算三维声场。此时，式(1)可简化为标准的二维Helmholtz 方程:

式中:p(r，z)为声压;k0=ω/c 为参考波数;n = c0/c为介质折射率。

采用Tappert 方法，令:

其中包络函数ψ(r，z)随距离缓变，Hankel 函数满足Bessel 微分方程:

当k0r ≥1 时，Hankel 函数可用其渐近形式代替:

利用方程(4)表示的Hankel 函数的性质，在远场假设(k0r ≥1 )条件下，将式(3)和式(5)代入式(2)中，即可得到如下简化的椭圆形波动方程:

引入“近轴近似”，小角度近似可表示为:

将近轴近似用于方程(6)，即可得到标准抛物方程:

采用Hardin 和Tappert 提出的分裂－ 步进FFT算法求解方程(8)，可得到某一方位上的二维声场数据，通过计算N 个方位上的二维声场即可近似地组成整个三维空间的声场。图1 为对三维声场进行离散采样的示意图，图中Δr 为距离分辨率，Δφ 为方位分辨率，Δz 为深度分辨率。

图1 三维声场离散采样示意图Fig.1 Sketch of discrete sampling of 3D acoustics energy field

采用N × 2D 弱三维计算三维声场可以进行2层并行。第1 层是不同方位间二维声场的并行计算;第2 层是求解方程内部算法的并行，如FFT算法的并行计算。利用曙光TC4000L 高性能计算平台，结合平台体系结构的特点，采用两级并行策略:将方位扇面分配到不同节点，每个方位扇面内单独执行与其相关的所有运算，计算过程中进程间没有数据交换，直至计算结束将结果合并输出，采用MPI 模型实现;FFT 算法的并行在节点内进行，计算结果存储在共享内存中，采用OpenMP 模型实现。

采用N ×2D 方法计算得到N 个采样方位切片上声传播损失随r － z 的变化，将合并输出的结果采用笛卡尔坐标的形式进行组织，不同方位上的声场数据依次存储到三维纹理中。

数据组织形式如图2 所示，图中纹理s 坐标表示水平距离r，纹理t 坐标表示深度z，纹理r 坐标表示方位φ。

图2 原始声场数据的三维纹理存储示意图Fig.2 Sketch of storage for 3D texture of original acoustics energy field data

采用三维纹理存储声场数据，可以通过设置三维纹理的参数，利用三次线性插值实现对任意方位上声场数据的重采样。由于使用了硬件加速，三次线性插值所带来的额外开销几乎可以忽略不计。

2 三维声场可视化

三维数据场体绘制的算法有很多种，其中基于GPU 的纹理映射体绘制算法利用了图形流水线的并行处理功能，是众多体绘制算法中速度最快的。基于纹理映射的体绘制算法主要包括以下2个阶段(见图3):第1 阶段是纹理生成，将体数据载入系统内存，并将数据转换为纹理数据格式存入显存中;第2 阶段是纹理绘制，根据一系列采样多边形对纹理数据进行重采样并混合生成最终的图像。

原始声场数据的三维纹理中依次存储了不同方位的声传播损失数据，是按照笛卡尔坐标系组织的，但原始声场数据是在柱坐标系下计算得到的，因此，在传统基于纹理映射体绘制算法的纹理绘制阶段，直接对原始声场数据的三维纹理进行重采样无法还原真实的三维声场。针对这一问题，本文在纹理绘制阶段引入了渲染到纹理的步骤，基本思路是:首先定义与帧缓冲区对象(FBO，Frame Buffer Object)绑定的空的三维纹理。然后对原始声场数据的三维纹理进行重采样和坐标转换，还原成真实三维声场的结构，并将结果渲染到与FBO 绑定的三维纹理上。最后再对与FBO 绑定的三维纹理进行颜色映射和光照计算等操作。最终图像生成阶段用于进行重采样的三维纹理是与FBO 绑定的三维纹理。改进后基于三维纹理映射体绘制的算法流程如图4 所示。

图3 基于纹理映射体绘制算法流程Fig.3 Flow of volume rendering algorithm based on texture mapping

图4 基于三维纹理映射体绘制改进算法流程Fig.4 Flow of improved volume rendering algorithm based on 3D texture mapping

2.1 数据重采样

对原始声场数据进行重采样及坐标转换的具体步骤如下:

1)创建与FBO 绑定的三维纹理。为还原柱坐标系下三维声场，与FBO 绑定的三维纹理的大小定义为:

式中:XFBO，YFBO，ZFBO分别为与FBO 绑定的三维纹理的长、宽、高;XOrig，YOrig分别为存储原始声场数据的三维纹理的长、宽。

2)设置视口及定义投影方式。视口大小设置如下:

glViewport (0，0，XFBO，YFBO);

定义二维正视投影矩阵，圆心位于几何中心:

gluOrtho2D (－1，1，－1，1);

3)对存储原始声场数据的三维纹理进行重采样。使用OpenGL 着色语言，在片元着色器中实现对原始声场数据的重采样及坐标转换，并将结果渲染到与FBO 绑定的三维纹理上。

2.2 颜色映射

颜色映射的目的是将声场数据映射为颜色并赋予不同的透明度值，通过传输函数实现。在水下目标探测过程中，当目标辐射噪声的传播损失大于某一临界值时(声呐优质因子，FOM)，将不能对目标实施可靠探测。为了在可视化过程中表达这一现象，采用非线性的颜色、透明度传输函数，如图5所示，当传播损失小于FOM 值时声场可见，反之接近透明(颜色与光强值接近0)。

图5 非线性传输函数示意图Fig.5 Sketch of nonlinear transfer function

2.3 光照计算

为了增强三维声场绘制的真实感，在片元着色器中引入了光照和明暗计算。采用Phong 局部光照模型，它包括环境光(Iambient)、漫反射(Idiffuse)和镜面反射(Isepcular)3 个部分［8］。

式中，Ka，Kd，Ks分别为环境光反射系数、漫反射系数和镜面反射系数;Ia为环境光强度;Ip为照射光源强度;n，l 和v 分别为归一化的表面法向量、指向光源的向量和指向视的向量;h 为v 和l 的中间向量;n为光照指数。

3 仿真实验

实验条件:声场半径120 km;计算范围内最大海深5 000 m;声源为水面目标，深度5 m;声源频率100 Hz。三维声场离散采样的距离分辨率500 m，方位分辨率10°，深度分辨率20 m，声场体数据的大小为241 ×36 ×251。

曙光TC4000L 高性能计算平台分配36 个节点分别计算不同方位的声场数据。同时，为了比较验证并行计算的效率，在高性能计算平台开辟了一个比较节点，用于串行计算。

图形显示计算机硬件性能:CPU:Inter (R)Core (TM)i7 － 4700HQ，8 核;显卡:NVIDIA GeForce GT 750M，4GB 显存。图6 是分别进行5 次并行和串行计算总耗时的结果比较:

水下三维声场数据在直角坐标空间和球面坐标空间的可视化效果分别如图7 和图8 所示。

图6 并行和串行计算结果比较Fig.6 Results compare in parallel and series computation

图7 直角坐标空间水下三维声场体可视化效果图Fig.7 Effect picture of Underwater 3D acoustics energy field in Cartesian space

图8 球面坐标空间水面目标三维声场体可视化效果图Fig.8 Effect picture of surface target 3D acoustics energy field in spherical space

在图形显示计算机中，图7 和图8 所示场景的刷新频率分别大于45 fps 和30 fps。从仿真实验结果看出:在高性能计算平台下对三维声场进行并行计算具有明显的加速效果，平均加速比约为20;改进后的基于三维纹理体绘制算法实现了水下三维声场的实时体绘制，能够生成高质量的水下三维声场图像。

4 结 语

本文实现了水下三维声场的快速计算和实时可视化表达，为作战人员感知和理解水下声传播规律提供了新的方法，拓展了指挥决策的方式。文中三维声场的并行计算主要依托实验室的高性能计算平台，近年来，随着高性能GPU 的飞速发展，基于CUDA 的并行计算越来越受到人们的重视和青睐，使PC 机进行复杂的科学计算成为可能。因此，下一步将重点研究基于CUDA 的三维声场并行计算以及相应的数据组织和体绘制方法。

［1］张林，笪良龙，李玉阳．三维水声传播建模及并行算法研究［J］．系统仿真学报，2013，24(1):54 －57．ZHANG Lin，DA Liang-long，LI Yu-yang．3D Underwater sound propagation modeling and parallel algorithm［J］．Journal of System Simulation，2013，24(1):54 －57．

［2］笪良龙，范培勤，陈冬滨．水平不变海洋声道中WKBZ简正波方法的并行算法［J］．系统仿真学报，2006，18(7):1980 －1982．DA Liang-long，FAN Pei-qin，CHEN Dong-bin．Parallel programme of WKBZ mode approach to sound propagation in horizontally stratified oceans［J］．Journal of System Simulation，2006，18(7):1980 －1982．

［3］李玉阳，笪良龙，宋洁．水声环境仿真及其可视化技术研究［J］．测绘科学技术学报，2009，26(1):8 －11．LI Yu-yang，DA Liang-long，SONG Jie．Research on underwater acoustic environment simulation and visualization［J］．Journal of Geomatics Science and Technology，2009，26(1):8 －11．

［4］王小林，笪良龙，杨廷武，等．基于预积分体渲染技术的水下三维声场体可视化［J］．系统仿真学报，2007，19(21):4970 －4972．WANG Xiao-lin，DA Liang-long，YANG Ting-wu，et al．Underwater acoustic energy fields visualization based on technology of pre-integrated volume rendering［J］．Journal of System Simulation，2007，19(21):4970 － 4972．

［5］杨廷武，王小林，李玉阳，等．基于三维纹理的水下三维声场直接体可视化［J］．系统仿真学报，2006，18(Suppl．1):322－323．YANG Ting-wu，WANG Xiao-lin，LI Yu-yang，et al．Research of direct volume rendering on underwater 3D acoustic energy field［J］．Journal of System Simulation，2006，18(Suppl．1):322 －323．

［6］王小林，笪良龙，杨廷武，等．水下3 维声场体可视化两种切平面实现方法［J］．中国图象图形学报，2008，13(2):312 －315．WANG Xiao-lin，DA Liang-long，YANG Ting-wu，et al．The realization of two clipping plane technologies in underwater acoustic energy fields visualization［J］．Journal of Image and Graphics，2008，13(2):312 －315．

［7］FINN B J，WILLIAM A K，MICHAEL B P，et al．Computation ocean acoustics［M］．Second Edition．New York，USA:Springer，2011:458 －459．

［8］KLAUS E，MARKUS H，JOE M K，et al．Real-time volume graphics［M］．Wellesley，Massachusetts，USA:A K Peters Ltd，2006:115 －116．