多波束海底地形等深线快速生成算法

李海森，陆丹，魏玉阔，么彬

(哈尔滨工程大学 水声技术国家级重点实验室，黑龙江 哈尔滨 150001)

多波束测深声呐作为当前海洋勘测的重要工具之一，被广泛应用于海洋资源调查、海底地形测量、水下打捞等领域.该声呐获得的厘米级分辨率测深数据经滤波处理后，可进行海底地形建模(digital terrain modeling，DTM)，实现海底地形可视化［1-3］.其中等深线图能够有效反映海底地形的高低起伏特征，是一种兼顾海底地形定性与定量信息的高层次表达手段.从海底地形中生成等深线图是多波束测深声呐海底成图软件的最重要功能之一.多波束测深声呐具有宽覆盖、高分辨、大面积扫海测量特点，随着多波束测深数据量的不断增大，如何提高大数据量海底地形的等深线生成速度成为必须面对的问题，也是国产多波束测深声呐必须解决的紧迫问题之一.

海底地形等深线作为一种深度属性等值线，其生成方法与常规等值线生成方法类似.多年来，众多学者对等值线快速生成方法进行了研究［4-8］.在格网或三角网拓扑关系已知条件下，等值线生成过程通常可分为2部分:1)遍历格网或三角网确定等值线搜索起点;2)从起点开始利用拓扑关系进行等值线跟踪直至生成完整的闭合或非闭合等值线.这两部分共同决定了等值线生成效率.研究人员通常采用各种索引结构来提高第一部分中搜索起点确定的效率［4-7］，采用顾及方向的策略简化网格边上等值点存在的判断条件来提高第二部分中等深线跟踪的效率［7］.然而，索引结构虽然能够提高搜索起点的确定效率但其构建和查询比较复杂并且占用空间较大;顾及方向的策略虽然在一定程度上简化了等值点判断条件，但是这种策略只适用于规则格网DTM，对不规则三角网(triangulated irregular network，TIN)结构的DTM并不具适用性.为同时提高这两部分的效率，实现等深线快速生成，本文提出一种基于等深值索引序列的等深线快速生成算法.该算法在提高等深线生成速度的同时更加适合计算机编程实现，并且能同时适用于规则格网和不规则三角网结构DTM的等深线快速生成.目前，该算法已在国产首台便携式多波束测深声呐的海底成图软件中成功实现，并在吉林松花湖以及湖北清江平洛湖水下地形等深线图显示中得到应用.

1 常规等深线生成

常规等深线可以从规则格网和TIN2种结构的DTM中生成.由于多波束测深声呐获得的测深数据呈不规则离散分布，且TIN结构的海底地形模型直接利用原始深度数据，避免了内插的精度损失，更适合建立TIN海底地形模型，因此本文重点研究从TIN海底地形模型生成等深线.文献［3］给出了使用野值剔除后的离散多波束测深数据进行TIN海底地形建模的方法.以此为基础，常规的等深线生成方法按三角形顺序进行搜索，基本过程如下［9］:

1)对给定的等深值h，与所有数据点zi(i=1，2，…，n)进行比较，若zi=h，则将zi作降一微量(如10－4)处理［10-11］.

2)设立三角形标志数组，设置其初始值为0，每一个元素与一个三角形对应，凡处理过的三角形将标志置为1，以后不再处理，直至等深值改变.

3)按顺序判断每一个三角形的三边中的2条边是否有等深线穿过.若三角形一边的2端点为P1(x1，y1，z1)，P2(x2，y2，z2)，则采用如下公式判断该边上是否存在等深点:

若上式成立，则说明该条边上存在等深点.获得搜索起点，通过线性内插计算等深点在该边上的平面坐标.

4)搜索该等深线在该三角形的出口边，也就是相邻三角形的入口边，并内插该等深点的平面坐标.

5)进入相邻三角形，重复第4)步，进行等深线跟踪.若等深点回到起点，则形成闭合等深线;若等深点到达三角网边界，则从搜索起点开始反向跟踪，直至到达另一边界，形成非闭合等深线.

6)当一条等深线全部跟踪完后，将其输出.然后继续三角形的搜索，直至全部三角形处理完成.改变等深值，重复以上过程，直到全部等深线绘制完成为止.

2 基于等深值索引序列的等深线快速生成算法

在常规等深线生成算法中，要初始化三角形标志数组，首先要判断三角形中包含哪些等深值.这就需要查找哪些等深值在三角形顶点的最小和最大深度值范围内.由于三角网中一个点总是被多个三角形共用，显然，常规算法在获取三角形中所含等深值时，同一数据点的深度值会在不同的三角形中被重复计算.另外，由于实测的深度值为浮点型数据，当DTM中数据量很大时，必然需要大量浮点计算，以致等深线生成速度无法满足多波束海底成图系统的实时显示需求.因此，本文提出一种基于等深值索引序列的等深线快速生成算法.算法先将浮点型深度属性三角网映射为整型索引属性三角网，根据索引序列和深度范围的映射关系避免了标志数组初始化中同一数据点在不同三角形中的重复计算过程，然后提出一种基于0/1异或的等深线出口边判断方法，将等深线跟踪过程建立在整型运算的基础上，通过降低计算机运算量来实现等深线快速生成.由于算法运算量的降低与DTM的具体结构无关，该算法同时适用于规则格网和不规则三角网结构DTM的等深线快速生成.

2.1 深度与索引属性的映射

DTM数据深度与索引属性的映射是将绝大部分运算从浮点运算转化为整型运算的关键，也是等值线快速生成算法的基础，主要包括设置整型索引属性和对每个数据点增加索引属性值2部分，处理过程如下:

1)根据等深值设置深度范围索引.假设Zmin和Zmax为DTM的最小和最大深度值，V1，V2，…，Vn(以升序排列)为待显示的n个等深值.显然，Zmin≤V1＜V2＜…≤Zmax，可得到个递增的深度范围，用一个从0到n的整数序列表示:C={0，1，2…，n}.

2)建立等深值与深度范围索引的直接映射关系.如图1所示，建立索引序列C={0，1，2，…，n}与深度范围的直接映射关系:C0=0，映射［Zmin，V1］深度范围;Ci=i(i=1，2，…，n－1)，映射［Vi，Vi+1］深度范围;Cn=n，映射(Vn，Zmax］深度范围.

图1 索引序列与深度范围的映射关系Fig.1 Relation of the index array and depth ranges

3)给DTM中每个数据点添加索引属性，并将浮点型深度属性三角网转化为整型索引属性三角网.例如，图2(a)中最大深度Zmin=22.3，最小深度Zmax= 27.2，待显示等深值为 V1=23.0，V2=25.0，V3= 27.0，则可得索引序列C={0，1，2，3}，转化为整型索引属性三角网后如图2(b)所示.由于索引序列与等深值存在直接映射关系，因此无需任何计算过程即可从索引属性三角网中直接得到任意三角形边上包含的等深值.三角形A的edge1上包含3个等深值分别为V1，V2和V3;edge2上包含2个等深值V2和V3; edge3上包含一个等深值V1.而三角形B的3个顶点索引值相等，因此该三角形中不包含任何等深值.显然，将深度属性三角网转化为整型索引属性三角网后，无需繁琐地对每个三角形顶点的深度值进行比较以判断其中是否包含等深值或包含了哪些等深值，大大减少了计算机运算量和程序复杂度.

图2 深度属性三角网转换为索引属性三角网Fig.2 Turn depth attribute triangulation net to index attribute triangulation net

2.2 基于0/1异或的等深线走向判断方法

在深度属性映射为索引属性的基础上，进一步研究如何快速获取等深线在三角形中走向的问题.

首先考虑最简单的情况，即三角形顶点索引值只为0或1的情况:假设三角形顶点的索引值分别为Cdot1，Cdot2Cdot3，且最小索引值为0，最大索引值为1.显然，只要三角形边端点索引值的异或结果为1，该边上必然存在等深点.因此，在已知等深线入口边条件下，只需判断三角形第三个顶点的索引值与入口边任一端点索引值的异或结果是否为1即可判断等深线出口边.如图3所示，假设等深线入口边为三角形的edge1，若Cdot1∧Cdot3=1，则判断等深线出口边为三角形的edge3;否则，等深线出口边为三角形的edge2.

图3 基于0/1异或的等深线走向判断Fig.3 Contour trend decision based on 0/1 exclusive-OR

然后延伸到一个三角形内通过多条等深线的情况:假设三角形顶点的最小索引值为Cmin=i，最大索引值为Cmax=j，通过索引值与深度范围的直接映射关系可知该三角形内包含Vi+1到Vj的(j－i)个等深值.要获得Vk(i+1≤k＜j)等深线在三角形内的走向，需进行如下处理将顶点索引值转化为0或1:

图4 3条等深线在三角形内的走向判定Fig.4 Trend decision of three contours in a triangle

2.3 等深线快速生成

本文提出的等深线生成算法主要由2部分组成:一是确定三角形中的等深线起点并初始化三角形标志数组;二是等深线跟踪.为提高等深线生成速度，算法先将深度属性映射为索引属性，在索引属性三角网上直接获取等深线起点并初始化标志数组，然后从起点开始使用0/1异或方法快速判断等深线走向进行等深线跟踪.算法基本流程如图5所示.

图5 等深线快速生成算法流程图Fig.5 Flow chart of fast algorithm for contours generation

在整型索引属性三角网上，遍历所有三角形，当三角形顶点的最大最小索引值之差d=Cmax－Cmin≠0时，表示该三角形内含有d个等深值.将当前三角形号、所含等深值数量d、d个等深值索引及其相应的标志按顺序存入标志数组中.然后，以标志数组中对应的三角形为起始三角形进行等深线跟踪.依次遍历当前三角形内的d个等深值，若当前等深值标志为1，表示该等深值已被搜索过，跟踪该三角形的下个等深值;若当前等深值标志为0，表示该等深值未被搜索过，进行等深线跟踪:先在起始三角形内获得起始等深点，使用0/1异或方法判断等深线在三角形内走向，获得出口点，以该出口点作为出口边相邻三角形的新入口点搜索新出口点，如此反复，直至等深线到达边界或返回起点.若等深线到达边界，则从搜索起点进行反向跟踪，直至到达另一边界，形成一条非闭合等深线;若等深线返回起点，则形成一条闭合等深线.三角形标志数组中的所有三角形均处理完毕后，该DTM的等深线生成完成.

3 算法运算量分析与测试结果

3.1 算法运算量分析

假设海底DTM中包含N个数据点，则最多能构成2N－5个三角形［12］，当N较大且数据点均匀分布时，取三角形个数T≈2N，等深值数量为k，k值可根据显示效果需求来选取.

在三角形标志数组初始化阶段，常规算法可分为2个步骤:1)查找数据点深度值zi(i=1，2，…，N)，是否等于等深值hj(j=1，2，…，k)，若相等，对zi作降微量处理.设降微量处理操作的运算量为c1，则该步骤的运算量为Nc1lbk.2)遍历T个三角形，根据每个三角形的最小最大深度值zmin和zmax查找升序等深值序列 {hi}(i=1，2，…)，k中满足zmin＜hj＜zmax条件的hj，其中h=a，a+1，…，b且1≤a≤b≤k，然后标记该三角形所含等深值并初始化.设在k个等深值中查找hj的运算量为c2lbk，初始化该三角形等深值标志的运算量为c3，则该步骤的运算量共为T(c2lbk+c3).取T≈2N，常规算法进行三角形标志数组初始化的总运算量近似为N(c1lbk+2c2lbk+ 2c3).

本文算法的三角形标志数组初始化也可分为2个步骤:1)将N个数据点的深度属性转化为索引属性，即查找升序等深值序列 {hi}(i=1，2，…，k)中满足z≤hj(j=1)，hj＜z≤hj+1(j=1，2，…，k－1)或z＞ hj(j=k)条件的等深值索引j，其中z为当前数据点深度值.设在k个等深值中查找z位置的运算量为c4lbk，则该步骤的运算量为Nc4lbk.2)遍历T个三角形，根据三角形顶点的最大最小索引值直接获取等深值索引并初始化标志数组.设初始化一个三角形等深值标志的运算量为c5，显然，c5=c3，因此该步骤的运算量减少为Tc3.取T≈2N，本文算法的三角形标志数组初始化总运算量近似为N(c4lbk+2c3).

由于c1、c2、c3和c4都大于0，且在升序等深值序列{hi}中查找满足zmin＜hj＜zmax(j=a，a+1，…，b且1≤a≤b≤k)的 hj所需运算量必然大于在{hi}(i=1，2，…，k)中查找满足z≤hj(j=1)，hj＜z≤hj+1(j=1，2，…，k－1，)或(z＞hj(j=k))的j所需运算量，即c2lbk＞c4lbk，因此，N(c1lbk+2c2lbk+ 2c3)大于N(c4lbk+2c3).由此可知，三角形标志数组初始化运算量随数据量的增加呈线性增长，但本文算法运算量的增长斜率小于常规算法，这意味着本文算法的低运算量优势随着数据量的增加而更加突出.

在等深线跟踪阶段，常规算法根据式(1)需要2次浮点数减法，1次浮点数乘法，1次关系运算来判断等深线在三角形内的出口边.而本文算法利用数据点的整型索引属性仅使用2次整型条件表达式操作，1次异或操作即可判断出口边.由于每次出口边判断都能节省部分运算量且计算机处理整型数据的速度高于处理浮点型数据的速度，因此在等深线跟踪阶段，本文算法所需运算量必然小于常规算法，且所节省运算量随数据量的增加而增大.

为验证本文算法正确性与低运算量优势，根据运算量与运行时间的正比关系，使用实际测深数据对两种算法的运行时间进行测试并给出测试结果.

3.2 算法测试结果

算法使用VC++和OpenGL实现，实验的硬件环境为Intel(R)Core(TM)i3，2.93 GHz，3.46 GB内存.实验数据为吉林松花湖区域和湖北清江平洛湖区域的实际水深数据.

图6和图7给出了本文算法对上述区域实测水深数据进行等深线生成得到的等深线图.图6中松花湖区域DTM数据量为12×104，最小与最大深度值为7.13 m和62.67 m，设置等深值范围为8～62 m，等深间隔1 m，等深值数量k=55，等深线生成耗时约112.8 ms.图7中平洛湖区域DTM数据量为14.5×104，最小与最大深度值为 10.3 m 和110.9 m，设置等深值范围为15～110 m，等深间隔5 m，等深值数量 k=20，等深线生成耗时约51.3 ms.由于图6中设置的k值大于图7中的k值，也就是说图6中等深线分布比图7更加密集，因此虽然松花湖区域DTM的数据量小于平洛湖区域，但是其等深线生成耗时大于平洛湖区域.

图6 松花湖部分区域等深线图Fig.6 Contour chart of part area of Songhua Lake

图7 平洛湖部分区域等深线图Fig.7 Contour chart of part area of Pingluo Lake

图8给出了6组不同数据量实验数据的本文算法和常规算法运行时间测试结果.测试数据的深度值在6.64～66.23 m之间，设置最小与最大等深值范围为8 m与66 m，等深间隔为2 m，即等深值数量k=30.

图8 本文算法和常规算法运行时间比较Fig.8 Runtime comparison of the algorithm in this paper and a general algorithm

测试结果表明本文算法在三角形标志数组初始化阶段和等深线跟踪阶段所需运行时间均小于常规算法.从图8(a)可以看出，三角形标志数组初始化运行时间呈线性增长，本文算法运行时间的增长斜率小于常规算法的增长斜率，符合运算量线性增长和斜率c4lbk+2c3＜c1lbk+2c2lbk+2c3的规律.从图8(b)可以看出，等深线跟踪运行时间也随数据量呈线性增长，本文算法的运行时间总是低于常规算法.图8(c)给出的2种算法的等深线生成总运行时间表明:1)本文算法的等深线生成总运行时间低于常规算法;2)等深线生成运行时间随数据量呈线性增长;3)本文算法的运行时间增长斜率小于常规算法.显然，与常规算法相比，本文算法的优势随着DTM数据量的增加而不断增大，尤其适合大数据量海底地形的等深线快速生成.

4 结论

多波束扫海测量数据量已达到百万、千万乃至海量量级，本文针对大数据量海底地形等深线生成算法进行了深入研究，提出一种基于等深值索引序列的等深线快速生成算法.通过实测海底地形数据处理分析，得出以下结论:

1)通过深度属性到索引属性的映射，将等深线生成过程建立于整型索引属性三角网的基础上，降低了算法的复杂度并为基于0/1异或的等深线走向快速判定提供了条件，从而同时降低了等深线搜索起点确认和等深线跟踪的运算量.

2)算法的低运算量优势随着数据量的增大更加突出，因此该算法尤其适用于数据量达到百万级、千万级乃至海量的多波束测深地形等深线快速生成.

3)算法同时适用于规则格网数据结构和不规则三角网数据结构DTM的等深线深快速生成，对DTM数据结构具有更广泛的适用性.

此外，由于海底地形等深线是一种深度属性的等值线，其生成方法与大部分等值线生成方法类似，因此本文算法不仅适用于海底地形的等深线快速生成，同样能够应用于其他属性等值线，如等高线、等压线、等温线等的快速生成，具有广泛的适用性.

［1］李海森，陈宝伟，么彬，等.多子阵高分辨海底地形探测算法及其FPGA和DSP阵列实现［J］.仪器仪表学报，2010，31(2):281-286.

LI Haisen，CHEN Baowei，YAO Bin，et al.Implementation of high resolution sea bottom terrain detection method based on FPGA and DSP array［J］.Chinese Journal of Scientific Instrument，2010，31(2):281-286.

［2］么彬，李海森，周天，等.多子阵超宽覆盖海底地形探测方法试验研究［J］.哈尔滨工程大学学报，2008，29(10): 1076-1081.

YAO Bin，LI Haisen，ZHOU Tian，et al.Test study of a multiple sub-array super-wide-coverage seabed terrain detection method［J］.Journal of Harbin Engineering University，2008，29(10):1076-1081.

［3］LU D，LI H S，WEI Y K，et al.An improved merging algorithm for Delaunay meshing on 3D visualization multibeam bathymetric data［C］//IEEE ICIA 2010.Harbin，China，2010:1170-1176.

［4］KREVELD M V.Efficient methods for isoline extraction from a TIN［J］.International Journal of Geographical Information Science，1996，10(5):523-540.

［5］王涛，刘纪平，毋河海.基于排序预处理的等高线生成算法［J］.测绘学报，2006，35(11):390-394.

WANG Tao，LIU Jiping，WU Hehai.The extraction of contour lines from grid DEM based on sorting［J］.Acta Geodaetica et Cartographica Sinica，2006，35(11):390-394.

［6］顾滨兵，杨兆海，高宇，等.基于网格和端点量化的等值线生成算法［J］.吉林大学学报:信息科学版，2010，28 (1):89-94.

GU Binbing，YANG Zhaohai，GAO Yu，et al.Algorithm of contour lines drawing based on grid using technique of grid scan and point digital［J］.Journal of Jilin University:Information Science Edition，2010，28(1):89-94.

［7］王涛，毋河海，刘纪平.基于区间树索引的等高线生成算法［J］.武汉大学学报:信息科学版，2007，32(2):131-134.

WANG Tao，WU Hehai，LIU Jiping.An algorithm for extracting contour lines based on interval tree from grid DEM［J］.Geomatics and Information Science of Wuhan University，2007，32(2):131-134.

［8］JONES N L，KENNARD M J，ZUNDEL A K.Fast algorithm for generating sorted contour strings［J］.Computers＆Geosciences，2000，26:831-837.

［9］邬伦，刘瑜，张晶，等.地理信息系统——原理、方法和应用［M］.北京:科学出版社，2005:215-216.

［10］迟宝明，李治军，叶勇，等.基于GIS的地下水水位等值线图自动生成算法研究［J］.吉林大学学报:地球科学版，2007，37(2):261-265.

CHI Baoming，LI Zhijun，YE Yong，et al.Study on the arithmetic of automatically drawing isoline of groundwater level based on GIS［J］.Journal of Jilin University:Earth Science Edition，2007，37(2):261-265.

［11］张立华.港口工程水深测量中的不规则区域、任意方向、不等尺度网格法追踪等深线研究［J］.水运工程，2005(4):9-13.

ZHANG Lihua.A study on drawing depth contours using grids based on irregular area，arbitrary direction and different scale in bathymetric survey for port engineering［J］.Port＆Water Engineering，2005(4):9-13.

［12］周培德.计算几何——算法分析与设计［M］.3版.北京:清华大学出版社，2008:89-91.