一种海底浅层声学探测数据综合可视化方法

刘加银，刘海行，苏天赟，李新放

（国家海洋局第一海洋研究所，山东 青岛 266061）

目前的海底探测主要依赖于声学探测技术，包括多波束测深、侧扫声呐影像和浅地层剖面探测，其工作原理基本相似，只是探测目标的不同而有所区别（金翔龙，2007）。

多波束声呐系统以获得高密度、高精度的位置信息为主，覆盖广，但成像质量不高，主要用于描绘海底地形地貌；侧扫声呐以成像为主，能直观的反映海底地貌形态、沉积物类型等信息（阳凡林，2003）；浅地层剖面勘测是一种基于水声学原理的连续走航式探测水下浅部地层和构造的地球物理方法，可以有效获得海底以下的浅部地层结构和构造，并分析出海底以下存在的灾害地质因素如埋藏古河道、浅层气、浅部断层、软弱地层和浅部基岩等（李平等，2011）。综合利用多波束测深数据、侧扫声呐数据和浅地层剖面探测数据，对海底地形地貌的认知、水下目标的探测、海底地质结构的鉴别和分类等都具有重要意义（赵建虎，2002）。然而，由于数据获取方式和数据格式的不同，多波束测深数据、侧扫声呐影像和浅地层剖面的数据集相互独立，这给海底地形地貌特征和地层结构的识别带来了困难。数据综合可视化是解释和理解海底浅层声学探测数据的有效手段。

为实现多波束数据和侧扫声呐影像的综合利用，阳凡林（2003） 对多波束和侧扫声呐数据的融合进行了论述，提出基于等深线和轮廓线的同名特征点对不同图像源进行匹配，然后对匹配后的数字信息使用小波变换进行融合可视化。赵建虎等（2013） 提出了ISURF 多波束图像和侧扫声呐图像改进匹配算法，以及适合二者融合的拉普拉斯金字塔法和小波变换法；利用该方法可以获得具有高精度位置信息和高质量图像信息的海床地貌影像。基于二维空间的数据融合可视化方法虽然能获得较高的精度，但并不能较全面的反映海底地质环境，也不能进行交互操作，在具体应用中给信息的进一步提取带来了困难。

姜小俊（2009） 提出了SFBAO 模型对浅地层剖面和侧扫声呐、多波束数据融合后进行GIS 表达。该模型将浅剖数据作为航迹线的一个属性，然后将航迹线与侧扫声呐或多波束数据基于统一的空间参考系进行融合可视化。SFBAO 模型实质上是以数据库的方式综合管理海量的多波束数据、侧扫声呐影像和浅剖数据，并且在此基础上添加可视化模块。由于侧扫声呐影像、多波束和侧扫浅剖数据并没有在同一视图下显示，在相对较小尺度的空间，该方法对海底地质三维环境显示不够直观。陈超等（2012） 提出利用多重纹理对海底地形与海底底质分类纹理进行三维可视化。然而侧扫声呐影像比海底地质分类纹理数据大很多，具体的实现又对纹理尺寸做出了比较严格的限制（Segal et al，2012）。此外，侧扫声呐影像分辨率较高，多波束数据分辨率较低，使用多重纹理时对地形三角网顶点指定纹理坐标容易出现错误。

侧扫声呐影像的分辨率很高，在局部范围内数据量也非常大，一些经典工具往往通过金字塔模型或者纹理切片的方式处理这样超大的纹理数据，如ArcGIS 和OSGEarth。然而，金字塔或者纹理切片的方式处理繁琐，并且在实际应用中需要反复进行纹理数据切换，此时经典方法效率不高，影响视觉效果的流畅性。

海洋信息可视化时，需要注意体现科学数据的可视化和满足视觉感受的可视化的统一，它们的区别在于是否可以从可视化结果中读取数值（苏奋振等，2014）。现有工具难以满足在对海底浅层声学探测数据三维综合可视化的同时，进行海底地质环境信息的提取。

为更加全面、直观、高效地反映海底地质环境和方便进行信息提取，本文提出了一种基于OpenGL 可编程管线的海底浅层声学探测数据三维可视化方法，并利用南海海底某峡谷区域获得的高分辨率多波束测深数据、侧扫声呐影像和浅地层剖面数据对该方法进行实现和应用。

1 GPU 原理

目前的OpenGL 可编程管线包括顶点处理器（Vertex Processor）、细分控制处理器（Tessellation Control Processor）、细分计算处理器（Tessellation Evaluation Processor）、几何处理器（Geometry Processor）、片元处理器（Fragment Processor） 和计算处理器（Compute Processor）。顶点处理器一次只能操作一个顶点，主要用于对图元顶点位置变换和颜色设置；细分控制处理器和细分计算处理器可以在GPU 中对模型对象进行细化，提升图像的渲染质量；几何处理器可以对整个图元进行处理，计算法线和改变原始的拓扑结构，如将三角形转变成点；片元处理器是着色器的最后一个阶段，最终的像素颜色由它来设置；计算处理器比较特殊，它并不是图像绘制管线的一部分，但是却可以对图像和存储缓冲区等OpenGL 对象中的数据进行改变（Kessenich et al，2012）。

处于效率的考虑，OpenGL 提供了缓冲区对象用于在显卡中存储频繁使用的数据。顶点缓冲区对象适合存储顶点数据，一致缓冲区对象以Group 的形式存放uniform 变量，纹理缓冲区对象则可以存储比一维纹理大几万倍的数据（付飞 等，2012）。帧缓冲区对象封装了帧缓冲区的状态，如像素的颜色值、深度值、模板值等。OpenGL 提供了机制，能够将帧缓冲区的状态渲染到纹理，从而可以在内存中获得这些状态值，实现内存与帧缓冲区的交互。

2 绘制流程

根据海底浅层环境建模的需求，本文采用的绘制流程如图1 所示。在内存中对多波束数据进行三角化处理，生成地形三角网后和浅剖数据进行求交运算，得到浅剖数据的实际深度值。将侧扫声呐影像数据转换成色度盘的形式也在内存中完成，并且详细记录侧扫声呐影像的坐标，大小等信息。经内存处理后的数据通过总线传输到显卡进行可视化处理：顶点处理器完成光照计算和颜色设置；侧扫声呐影像数据和多波束数据的三维综合显示在片元处理器中完成。通过OpenGL 提供的接口读取像素的深度值反算出屏幕像素点对应的实际坐标，而计算出的实际坐标可以用来提取地形剖面线和浅地层剖面的相关信息。

图1 本文采用的绘制管线

3 算法描述

3.1 海底地形可视化

海底地形可视化首先需要根据点云数据生成DTM，为增强三维效果，还需要设置颜色、通过法线计算光照。多波束测深点云数据重构海底地形已经有许多相关的探讨和成熟的工具（Joseph，1998；武强等，2011；Zhang et al，2012），一些开源代码也可以用来生成每个顶点的法线（Wang et al，2012；VTK 6.2.0，2014；PCL 2014）。

对海底地形颜色和光照的设置在顶点处理器中完成。先根据地形的高程分布设置多个边界值height [i]，每个边界值分别对应一种颜色color [i]。设顶点的高程值为z，且，则该点的颜色通过公式（1） 进行插值得到。

本文采用Phong 模型对海底地形进行光照渲染，公式如（2） 所示。

式中：La、Ld、Ls分别为环境光、漫反射光和镜面反射光的强度，Ka、Kd、Ks对应地形表面环境光反射率、漫反射光反射率和镜面光反射率，s 为光源方向向量，n 为顶点法线向量，r 为反射光的方向向量，v 为视点方向向量，参数f 可以根据需要进行调节（Wolf，2011）。

地形剖面可以直观地对目标路径的地形变化情况进行表达。交互操作时，在屏幕上指定一条线段，结合模型视点投影矩阵和像素点的位置以及深度值得到线段端点对应的空间三维坐标；然后将指定线段进行细分，得到一系列（x，y） 坐标对，每个坐标对可以生成一条端点坐标为（x，y，10 000）、（x，y，-10 000） 的线段，通过该线段和地形三角网求交可以得到坐标对的深度值。

3.2 侧扫声呐影像与多波束数据综合显示

侧扫声呐影像和多波束测深数据综合显示是通过将侧扫声呐影像作为纹理添加到多波束测深数据所构建的三角网上实现的。侧扫声呐影像数据经过了领域专家的相关处理，具有可靠的位置信息，但此时侧扫声呐影像还只是二维数据，与多波束数据相比，具有分辨率高、数据量大等特点。侧扫声呐影像之间的格式并不一致，有jpg 和Geotif 两种格式，jpg 格式的影像只有四个角点的坐标，对应的区域不规则；侧扫声呐影像对应的地形区域之间也有重叠。传统的空间数据处理软件难以将两种数据进行镶嵌匹配和综合显示。

3.2.1 侧扫声呐影像数据存储

侧扫声呐影像的分辨率很高，在局部范围内数据量非常大，尺寸也不统一。通常情况下OpenGL要求纹理对象（TO） 的大小是2 的n 次幂，并且对TO 的大小作出了限制；而纹理缓冲区对象（TBO）的限制却十分宽松。如，在本文所使用的编程环境中，TO 的最大值为16 384 Byte，一张大小为3 950×9 500 的侧扫声呐影像即使只用1 Byte 记录每个像素的灰度值，也至少需要分割成2 291 份纹理瓦片，同时需要保证每张图像的尺寸为2 的n 次幂；而TBO 的最大值为134 217 728 Byte，可以一次存储3 张以上这样的侧扫声呐影像，避免实际应用中反复切换纹理数据带来的延迟问题，且不必将纹理的尺寸处理成2 的n 次幂。

侧扫声呐影像数据在CPU 中处理后，得到影像的色度盘，像素的颜色索引值存储在TBO 中，这样对于RGB 格式的图像可以节省大约2/3 的显存空间。单个TBO 存储空间有限，无法容纳所有的侧扫声呐影像数据，因此要将侧扫声呐影像进行分组和切割，分别存进多个TBO 中。为在片元处理器中获得正确的侧扫声呐影像数据，需要计算每个TBO 覆盖的范围、单张影像的覆盖范围、在纹理缓冲区的位置、侧扫声呐影像的尺寸等信息；由于相关信息的数量较多，超出了着色器对uniform变量个数的限制，所以将它们和色度盘组织在一起存储到UBO 中。这些数据只会通过总线（Bus） 往GPU 传输一次，因此减轻了带宽的限制。

3.2.2 纹理映射

传统的纹理贴图方法需要为三角网的每个顶点指定纹理坐标，该过程不仅耗时，并且在侧扫声呐影像和地形数据分辨率不一致的情况下，会出现图2 所示的情况。其中，三角形ACD、ABC 代表多波束数据生成的三角网，红色网格代表侧扫声呐影像，黑色虚线代表屏幕上的像素。由于侧扫声呐影像数据的分辨率比多波束数据的分辨率高很多，所以三角形ACD、ABC 只有部分顶点被侧扫声呐影像覆盖。传统方法只指定三角形的各个顶点的纹理坐标，三角形区域内的纹理坐标通过这三个点自动插值得到，然后固定管线依据纹理坐标进行纹理贴图。如果给C、D 点指定纹理坐标，那么四边形CDFE 覆盖范围内的像素点P1 会获得不该有的纹理，三角形AFE 覆盖区域内的像素点P2 的纹理值也会出现错误。若不给CD 指定纹理坐标，那么AFE 区域内的像素就会丢失纹理。通过细化三角形ACD 然后再指定纹理坐标的方法效率较低。

图2 分辨率不一致的纹理贴图

本文提出根据像素对应的实际坐标值求取纹理坐标的方法。在顶点着色器中设置一个传出变量作为顶点的属性，并将该变量初始化为三角网顶点的坐标值；经过图1 的栅格化阶段时，该变量会自动进行插值，片元处理器中获得的就是每个像素对应的三维空间坐标。在进行纹理贴图时，首先判断该像素是否被纹理覆盖，然后才从纹理缓冲区中获取正确的灰度值进行叠加。在小比例尺显示时，每个像素对应的是一块相对较大的区域，覆盖该区域的可能是侧扫声呐影像的多个灰度值，本文的方法只取其中一个代表整个区域。由于海底地质环境的连续性，该方法不会对宏观效果产生影响。在大比例尺显示时，屏幕上每个像素能够对应侧扫声呐影像的一个像素，所以并不影响侧扫声呐影像的分辨率。算法流程如图3 所示。

图3 纹理贴图过程

其中初始化阶段包括侧扫声呐影像灰度数据及坐标信息的加载，记录当前片元覆盖影像张数的变量和累加颜色清零等。侧扫声呐影像数据具有Qmips、.SDF、.JSF、.XTF 等多种格式。在实际应用中，通常是将原始的侧扫声呐数据转换成通用的tif、Geotiff 或JPG 格式（韩春花 等，2011）。计算侧扫声呐影像和片元的对应关系时，因侧扫声呐影像的格式不同而有所区别。本文针对JPG 格式和Geotiff 格式的数据进行阐述。

（1） Geotiff 图像

Geotiff 是TIFF 格式的一种扩展，在TIFF 的基础上定义了一些GeoTag 来对各种坐标系统、椭球基准、投影信息等进行定义和存储（Ritter，2000）。本文利用GDAL 对Geotiff 进行处理。GDAL 使用抽象数据模型来解析它所支持的数据格式，抽象数据模型包括数据集、坐标系统、仿射地理坐标变换、大地控制点、元数据、图像结构域、栅格波段和颜色表等（GDAL，2014）。其中仿射变换包含6 个参数，可以通过GDAL 接口GDALDataset::GetGeoTransform（GT [6]） 获取。在朝北方向的图像中，GT [2] 和GT [4] 的值为0，（GT [0]，GT[3]） 表示图像左上角的横纵坐标值，GT [1] 表示像元宽度，GT [5] 表示像元高度。图像行列号和地理空间之间的变换关系如下:

其中Yline、Xpixel 表示片元在它所在侧扫声呐影像中的行列号，由于处理的是正北朝向的图像，所以很容易判断当前片元是否被侧扫声呐影像覆盖，被覆盖片元对应的侧扫声呐影像像素计算公式为:

式中：floor 为向下取整函数，texelBufferSampler 对应OpenGL 的纹理缓冲区对象，texelFetch 为纹理获取函数，xsize和offset 分别对应图像的宽度和在纹理缓冲区中的起始位置。纹理缓冲区中只存储了颜色的索引值，colArray 是UBO 中存储真实颜色的色度表。Pixel.x 和pixel.y 分别代表当前像素对应的实际地理坐标。

（2） JPG 图像

JPG 图像只提供了四个角点坐标，如图4 所示：ABCD 4 个顶点坐标及其纹理映射关系已知。首先将图像对应的区域分为两个三角形，以三角形ABC 为例，获得任意像素P 的实际坐标后，可得到如下3 个变量：

图4 纹理空间和坐标空间之间的映射

式中：大写字符代表该点对应的三维坐标。若s1，s2，s3 的第3 个分量均不小于零，则当前像素被侧扫声呐影像覆盖。该像元的纹理坐标需要通过插值计算获得，计算公式如（6） 所示。

式中：coordA、coordB、coordc为三角形ABC 顶点的纹理坐标值，本例中分别为（0.0，0.0），（1.0，0.0），（0.0，1.0）；S（△PBC） 为三角形PBC 的面积，依此类推。落入BC 边的点只在B、C 两点之间插值即可。最终灰度值索引在纹理缓冲区中的位置为：

其中，offset 为覆盖侧扫声呐影像在纹理缓冲区中的起始位置，height，width 为影像的高度和宽度，Px、Py分别为点P 纹理坐标的横坐标值和纵坐标值。

3.3 浅地层剖面数据可视化

原始的浅地层剖面图像内容包括干扰图像、多次反射图像、地层界面线、地层层位图像等。这些图像显示具有一定灰度的点状、线状和面状图形，不具有整体的地理方位性和可量测性，一般将浅地层剖面数据作为航迹线的一个属性和其他声学探测数据进行集成（姜小俊，2009）。原始浅地层剖面图像经过领域专家的地质解释之后，得到一系列解释剖面点，格式为（x，y，id，layer1，layer2，layer3，layer4….），其中x，y 为测点的地理坐标，id 为该点的编号，layer1 代表第一层地层界面和海底地形之间的厚度，layer2 代表第二层地层和第一层地层之间的厚度，依此类推；缺失的地层用无效值标记。进行综合可视化处理时，首先作经过浅剖测点的直线和地形求交得到该测点对应的地形坐标，然后和各层厚度累加得到各层界面的实际坐标值。

转换后的数据按测线对地层进行分层渲染得到栅状图，把同一层的离散顶点进行三角化处理得到层面图，通过似三棱柱三维数据建模方法（武强等，2011） 可以得到体模型。

4 实例分析

本文提出的方法在南海海底峡谷区域的海底地质环境综合显示和分析中进行了应用，数据来源于国家海洋局第一海洋研究所2011年的海洋调查成果。其中多波束测深系统的分辨率为20 m，采用EM302 多波束测深系统收集；浅地层剖面数据通过全海洋浅地层剖面仪（Topographic Parametric Sonar，TOPAS） 探测获得，使用的脉冲信号频率为3.3 kHz；侧扫声呐影像使用无缆水下机器人（AUV） 获取。数据覆盖区域的面积为214.15 km2，多波束测深数据生成的DTM 由372 687 个三角形组成，侧扫声呐影像580M，浅地层剖面处理后形成的测点7 万多个。在配置GTX650Ti 显卡的机器上，基于vc++和OSG 实现了系统原型。图5 为实验数据的综合显示效果，在实现海底浅层声学数据综合可视化的同时，可以提取地形剖面和浅地层剖面数据。图中蓝色矩形内的粗线为拉选航迹线的高亮显示，右上角小窗口为该航迹线对应的浅地层剖面数据和地形数据提取结果的可视化界面，通过拉选可以实现海底浅层数据的查询和提取。图6为地形和侧扫声呐影像的贴图效果，从图6（b）可以看出，大比例尺环境下，侧扫声呐影像可以很好的与地形进行匹配显示，分辨率几乎与原始影像一致。

图5 海底地质环境综合显示效果图

图6 侧扫声呐影像和多波束综合可视化效果图

实验结果表明，此方法数据加载量大，视觉效果好，能直观、全面的反映海底浅层环境，满足实际应用中高效运行的需求。

5 结论与展望

与二维表达方式相比，三维GIS 能更加直观、真实的表达海底地形、地貌和地质结构。本文提出了一种基于OpenGL 可编程管线的海底浅层声学探测数据三维综合可视化方法，在满足视觉感受的同时实现了多维数据的交互式提取。使用该方法可以对多种格式的侧扫声呐影像进行纹理贴图，具有数据加载量大的优点，并且不受侧扫声呐影像和多波束数据分辨率不一致的限制。

下一步将根据实际应用需求，在现有工作基础上参考纹理切片和金字塔模型，研究在超大规模场景中多种高分辨率海底地质环境数据的综合可视化。

Estimating Surface Normals in a PointCloud. http://pointclouds.org/documentation/tutorials/normal_estimation.php

GDAL Data Model. http://www.gdal.org/gdal_datamodel.html.

Joseph O R, 1998. Computional Geometry In C. Cambridge University Press.

Kessenich J, Lunar G. The OpenGL Shading Language. 2012.8/2014.5.http://www.opengl.org/registry/doc/GLSLangSpec.4.30.pdf

Niles D R. GeoTIFF Revision 1.0 Specification. 2000/2014.http://www.remotesensing.org/geotiff/spec/geotiffhome.html.

Segal M,Akeley K.The OpenGL Graphic System:A Specification.2012.6/2014.5.http://www.opengl.org/ registry / doc /glspec44. core.pdf VTK 6.2.0 Documentation. 2014.3/2014.5. http://www.vtk.org/doc/nightly/html/index.html.

Wang R, Qian X L. OpenSceneGraph 3 Cookbook. 2012/2014. http: //www.packtpub.com/openscenegrap-3-for-advanced-3d-programming-using-api-cookbook/book.

Wolf D, 2011. Open G L 4.0 Shading Lauguage Cookbook. Packt Publishing Limited.

Zhang L, Bruno J, Zheng F, 2012. 3D Reconstruction of seabed surface through sonar data of AUVs.Indian Journal of Geo-Marine Sciences(0379-5136),41 (6) :509-515

陈超，王文珂，王怀辉，等，2012.一种海底地形与底质的三维融合可视化方法.系统仿真学报，24（9）：1 936-1 944.

付飞，李艳辉译，2012.OpenGL 超级宝典.北京：人民邮电出版社.

韩春花，张俊明，梁建峰，等，2011.侧扫声呐探测数据管理系统设计与实现.海洋通报，30（20）：188-193.

姜小俊，2009.海底浅层声学探测空间数据集成与融合模型及GIS表达研究.浙江：浙江大学.

金翔龙，2007.海洋地球物理研究与海底探测声学技术的发展.地球物理学进展，22（4）：1 243-1 249.

李平，杜军，2011. 浅地层剖面探测综述. 海洋通报，30（3）：344-349.

苏奋振，吴文周，平博，等，2014.海洋地理信息系统研究进展.海洋通报，33 (4）：361-370.

武强，徐华，2011.虚拟地质建模与可视化.北京：科学出版社.

阳凡林，2003.多波束和侧扫声呐数据融合及其在海底底质分类中的应用.武汉：武汉大学.

赵建虎，2002.多波束深度及图像数据处理方法研究.武汉：武汉大学.

赵建虎，王爱军，郭军，2013.多波束与侧扫声呐图像区块信息融合方法研究.武汉大学学报（信息科学版），38（3）：287-289.