徐 涛，曹枝东

(1. 中国民航大学计算机科学与技术学院 天津 东丽区 300300; 2. 中国民航大学中国民航信息技术科研基地 天津 东丽区 300300)

随着民航业的发展以及人们环境意识的逐步提高，机场噪声问题日益突出。合理规划使用机场周围的土地是解决机场噪声问题的有效手段，而机场噪声等值线图则是确定机场噪声对居民的影响范围、控制机场噪声以及合理规划机场周围土地使用的重要依据。因此，绘制机场噪声等值线图是机场噪声控制工作和机场规划设计的重要环节[1]。

等值线图是数字高程模型的表示形式之一[2]，其绘制一般分为以下步骤：离散数据点网格化、开放等值线起始网格单元的确定及追踪、封闭等值线起始网格单元的确定及追踪、等值线光滑、等值线填充[3]。文献[4]介绍了一种栅格右关联标志的规则网格数据等值线绘制算法，该算法虽避免了等值线交叉，但在等值点追踪上对某个特定网格点需查看与其相邻的4个网格单元，在存储空间和搜索方式上具有一定的缺陷，算法效率不高。文献[5]提出了一种基于搜索圆法的等值线追踪算法，该算法具有快速扫描、追踪效率高等优点，但仍需要先追踪开放等值线，后追踪封闭等值线，如网格边界不规则，该方法效率会下降。文献[6]针对数字地形图中等值线的提取提出了一种改进的摩尔相邻追踪算法，可以完整地追踪到整个地形图等值线，但对同一个网格单元的重复追踪需回溯到前一个网格单元，一旦很多网格单元被回溯时，算法效率急剧降低。机场噪声问题过去几十年受到极少关注，传统的环境噪声等值线图大多结合现有软件Surfer加载噪声数据绘制而成[7]，而机场噪声具有影响范围广以及因航迹覆盖范围不同而噪声分布差异性大等特点[8]。

基于上述原因，本文提出了一种基于路径栅格的机场噪声等值线追踪算法，通过构造有效网格、生成路径栅格，利用入口方向和路径栅格顶点编码和准确判断等值线走向趋势，唯一确定下一个等值点，避免封闭等值线和开放等值线[9]的分别追踪，可以快速生成机场噪声等值线图。

1 基本概念

1.1 等值线基本性质

等值线需满足以下要求[10]：

1) 通常为一条光滑连续曲线；

2) 给定高程值相应域上的等值线不限于一条；

3) 等值线可以是封闭或者开放的；

4) 不能相互交叉。

从这些特点可以看出，对于特定高程值的等值线(不限于一条)，与特定网格单元的交点个数可能为0、2、4[11]，与特定网格单元的一条边上的交点个数为0或1。

1.2 顶点编码

不失一般性，以规则地形网格为例，一个规则地形网格水平和垂直方向分别由M和N个等距离排列的点组成，其中整个网格用G表示，单个网格单元用cells[i,j]表示，网格单元顶点用pts[i,j]表示，其中网格单元的坐标为该网格单元右下角顶点的坐标，如图1所示。当某条等值线的高程值与网格单元的顶点属性值相等时，等值线刚好经过该点。为避免这种情况对后续算法造成不便，沿袭一贯做法，在该网格单元顶点属性值上加一个微小的数e，使等值线不直接经过网格单元顶点。

顶点编码是特定高程值的等值线绕过网格单元的顶点而得到约定的固定数值。任意网格单元只有4条边(不包括顶点)与等值线相交，以网格单元的4个顶点为参考，如图2a～图2d所示。当网格单元顶点A、B、C、D的相邻两边与同一条等值线相交时，约定该网格单元的顶点编码为1、2、4、8。当一个网格单元的两个顶点被特定高程值的等值线绕过时，将网格单元的顶点编码相加得到网格单元的路径栅格顶点编码和。

图1 规则网格

1.3 路径栅格顶点编码和

相对于某个特定高程值，任意网格单元与等值线的交叉关系只有3种情况：

1) 没有等值线经过：约定路径栅格顶点编码和PointCodeSum= 0。

2) 有一条等值线经过：每个网格单元中只有一个顶点被等值线环绕，如图2a～图2d所示，于是网格单元的路径栅格顶点编码和分别为1、2、4、8。

3) 有两条等值线(不同连续等值线上的片段)经过[11]：网格单元分别有两个顶点被等值线环绕，如图2e、图2f所示，网格单元的路径栅格顶点编码和分别为1+4=5和2+8=10，在图2g、图2h状况下，分别约定以AB点和BC点为参考，网格单元的路径栅格顶点编码和分别为1+2=3和2+4=6。

图2 顶点编码及路径栅格顶点编码和

2 基于路径栅格的机场噪声等值线追踪算法

2.1 关键步骤

2.1.1 有效网格的构造

机场及其附近区域的噪声都是由航班的飞行引起，然而机场四周并非所有地方都有航班经过，距航迹较远的地方噪声较小(可以不考虑)。图3为一个100*100的网格(其中浅色部分无噪声，深色部分有噪声)，U、V、W、X、Y、Z6个区域无噪声，剔除这些区域的网格，于是有效网格为S=(G-U-V-W-X-Y-Z)。

图3 机场噪声有效网格

2.1.2 路径栅格的生成

路径栅格(相对特定高程值)是在有效网格中剔除没有等值线经过的网格单元后剩余的网格部分[12]，其中每个网格单元都有路径栅格顶点编码和(不为0)，它是进行等值线追踪的依据。路径栅格的生成步骤如下。

1) 从高程值列表中选取一个高程值Hi，依次扫描有效网格，对每个网格单元做如下处理：分别计算Hi与网格单元cells[i,j]4个顶点(A、B、C、D)属性值Attribute的差值，记为s1、s2、s3、s4，若有为0的情况，在差值上加一个微小的数e，并令F=s1s2s3s4。

①s1、s2、s3、s4都大于0或都小于0时，该网格单元的路径栅格顶点编码和cells[i,j].PointCodeSum=0。

②F<0(如图2a～图2d情况)，做如下判断：

③F>0和s1s3>0(如图2e、图2f情况)，为避免两条等值线交叉情况的出现，这里判断等值线的走向需取网格单元的中心点P，令P点属性值为centerAttr，则：

2) 遍历完毕有效网格，高程值Hi的路径栅格形成。

2.1.3 等值线的追踪

基于生成的路径栅格，构造网格单元出口方向判断表，如表1所示。通过表1确定网格单元的出口方向。在等值线追踪判断出口方向时需要入口方向和路径栅格顶点编码和，确定了起始追踪网格单元，可以得到路径栅格顶点编码和。如该网格单元的两边或者四边有等值线穿过，可任取其中一个方向作为入口方向。

表1 网格单元出口方向判断表

结合表1，给出一个等值线追踪的示例。给定高程值Hi=30，已知起始追踪网格单元cells[1,10]，如图4所示，该网格单元PointCodeSum=3，上下两条边有等值线穿过。1) 取inDirection=上，由表1知，该网格单元outDirection=下，等值线进入网格单元cells[2,10]中，该网格单元PointCodeSum=10，inDirection=上，由表1知，该网格单元outDirection=右，LL，等值线进入网格单元cells[1,13]中，该网格单元PointCodeSum=3，inDirection=下，由表1知，该网格单元outDirection=上。等值线追踪到边界，Hi=30的一条等值线追踪完毕。2) 取inDirection=下，则outDirection=上，等值线追踪到达边界，而此时开放的等值线并没有追踪完毕(两端顶点都为边界点的等值线为开放等值线)，需记录等值线起始追踪网格单元，记Begin[i,j]=cells[1,10]，当追踪没完成，返回起始追踪网格单元，更换进入方向继续追踪，直至到达边界。

图4 等值线追踪示例

在等值线追踪过程中，需确定等值线与网格单元边的交点。已知网格单元边的两个顶点(P1、P2)和特定高程值Hi，可以通过线性插值[11]来确定交点Pcross坐标，追踪完毕顺次连接交点即可绘制需要的等值线。

2.2 算法描述

综上讨论，对于规则的机场噪声数据，可以通过路径栅格等值线追踪算法进行绘制，算法完整描述如下：

1) 构造有效网格。读取机场规则噪声数据，将网格单元顶点保存为二维数组pts，网格单元保存为二维数组cells。剔除明显没有等值线经过的网格区域，形成初始遍历的有效网格S。

2) 生成路径栅格。对于特定的等值线高程值Hi，遍历S区域网格单元，确定每个网格单元的路径栅格顶点编码和。若PointCodeSum>0，将网格单元插入路径栅格链表CellsArray末尾。

3) 确定等值点。取CellsArray首元素，以此为等值线起始追踪网格单元，利用路径栅格顶点编码和任取一个可能的入口方向，记录起始网格单元(Begin[i,j])和起始入口方向，依次确定下一个等值点，通过线性插值计算等值线与网格单元边的交点并记录，顺次存入等值点链表Path中(将起始追踪网格单元的入口方向边的等值点作为Path首元素)。

4) 追踪判断等值点。若正在追踪的等值点到达边界并且起始等值点Path[0]为边界点，说明等值线为开放等值线，当前等值线追踪完毕；若正在追踪的等值点到达边界但Path[0]不为边界点，则重新以Begin[i,j]为起始追踪网格单元，更换入口方向转入步骤3)继续追踪，直至追踪的等值点到达边界；若正在追踪的等值点与Path[0]相同，说明等值线为封闭等值线，当前等值线追踪完毕。在追踪过程中，每次遍历一个网格单元，取当前被等值线绕过的顶点编码number，计算PointCodeSum= PointCodeSum-number，若PointCodeSum=0，则从CellsArray中剔除网格单元。

5) 依次连接Path中相邻两点，高程值Hi的一条等值线绘制完毕。

6) 重复步骤3)～步骤5)，直至路径栅格链表CellsArray长度为0(高程值Hi的所有等值线绘制完毕)。

7) 重复步骤2)～步骤6)，直至预先设定的等值线高程值Hi集合遍历完毕。

3 实验结果及分析

实验使用两个机场的噪声数据集，共选取6组实验数据，分别是规模为100*100、200*200、300*300、300、400*400、1000*1000、1 989*1 989的规则地形网格数据。选取机场噪声分贝值层数分别为10、30、66。其中10层的噪声分贝值数据为：{40,45,50,55,60,65,70,75,80, 85}，30层是从30 dB到88 dB范围以2 dB为间隔的噪声分贝值数据构成的集合，66层则是从20 dB到85 dB范围以1 dB为间隔的噪声分贝值数据构成的集合。实验机器配置：Intel® Core™2 Duo CPU E4600 2.40 GHz，内存2.00 GB，操作系统为Microsoft Windows XP Professional 2002 Service Pack 3。

本文算法绘制的等值线图实验结果如图5、图6所示，图5是数据规模为400*400层数为10的某单跑道机场噪声等值线图，图6是数据规模为100*100层数为30的某双跑道机场噪声等值线图。从图中可以看出，即使在噪声分贝值很密集的情况下，绘制的等值线图仍能够很好地符合等值线的基本特点，特别是没有等值线交叉的情况出现。图5中心的短粗线部分为机场跑道，沿着跑道两端的圆点曲线为该机场主要航迹，从中可以看出，单跑道机场的噪声等值线以跑道为中心，依附航迹分散开来，任意两条相同分贝值间隔的相邻等值线，距离航迹越近的区域等值线越密集，噪声分贝值也越大。双跑道的情况与单跑道情况类似，只是介于两条跑道航迹之间的区域噪声分贝值会被叠加，其值大于单跑道航迹上该地的噪声分贝值(如图6中的A点)。

此外，本文选取3种算法(传统的等值线传播算法，文献[5]中的路径栅格边界追踪算法和本文算法)对上述6组实验数据3种层次的噪声分贝值进行算法时间效率对比，其结果如表2所示。从表中可以看出，数据规模为100*100、200*200、300*300、400*400时，本文算法时间效率总体上优于传统的等值线传播算法，明显优于路径栅格边界追踪算法；当数据规模增大至1 000*1 000和1 989*1 989时，与传统的等值线传播算法和路径栅格边界追踪算法相比，本文算法效率优势更为明显。

图5 某单跑道机场噪声等值线图

图6 某双跑道机场噪声等值线图

表2 3种算法效率比较

本文算法时间效率较高的一个重要原因就是路径栅格链表的使用。算法在构造路径栅格时，将特定高程值等值线经过的所有网格单元以链表存储，形成路径栅格链表，如图7所示，图中深色部分为等值线经过的网格单元，等值线追踪时只需从路径栅格链表中取出首元素(某个深色网格单元)开始追踪即可，从而避免常规的先扫描边界网格单元追踪开放等值线，再扫描内部网格单元追踪封闭等值线。因此通过预先构造路径栅格，建立路径栅格链表，避免了开放等值线和封闭等值线的分开追踪，提高了等值线追踪效率。

图7 路径栅格

4 总 结

等值线追踪算法很多[13-14]，对于不同的应用环境需要考虑不同的具体情形。基于路径栅格的机场噪声等值线追踪算法首先根据机场噪声数据的特殊情况建立有效网格，扫描有效网格确定网格单元的路径栅格顶点编码和，从而建立路径栅格，在初始追踪网格单元和入口方向确定的情况下，通过路径栅格顶点编码和唯一确定下一个等值点，并且路径栅格的提前建立避免了传统的先开放后封闭的等值线追踪方式。对比实验可以看出，本文算法在时间效率上具有明显优势，绘制的等值线图很理想，没有两条等值线交叉情况出现。

[1] 李冉. 机场航空噪声预测及其影响因素研究[D]. 天津:中国民航大学, 2008.LI Ran. Airport air noise prediction and influencing factors study [D]. Tianjin: Civil Aviation University of China, 2008.

[2] KIDNERA D B. High-order interpolation of regular grid digital elevation models[J]. International Journal of Remote Sensing, 2003, 21(14): 2981-2987.

[3] 苗润忠. 光滑的等值线生成算法[J]. 长春理工大学学报,2004, 27(1): 16-18.MIAO Run-zhong. A new smoothed contour lines generating algorithm for quadrilateral meshes[J]. Journal of Changchun University of Science and Technology, 2004,27(1): 16-18.

[4] 王结臣, 钱晨晖, 芮一康. 栅格数据生成等值线的一种实用方法[J]. 测绘科学, 2007, 32(6): 88-90.WANG Jie-chen, QIAN Chen-hui, RUI Yi-kang. A practical method for contour generation based on raster data[J].Science of Surveying and Mapping, 2007, 32(6): 208-282.

[5] 常会, 柴华彬, 邹友峰. 基于搜索圆法的等值线追踪技术[J]. 测绘科学, 2009, 34(2): 119-121.CHANG Hui, CHAI Hua-bin, ZHOU You-feng. Isoline tracing technique based on search-circle[J]. Science of Surveying and Mapping, 2009, 34(2): 119-121.

[6] PRADHAN R, KUMAR S, AGARWAL R, et al. Contour line tracing algorithm for digital topographic maps[J].International Journal of Image Processing, 2010, 4(2):156-163.

[7] 过春燕, 张邦俊. 基于Surfer的机场噪声等值线计算机绘制方法[J]. 中国环境科学, 2003, 23(6): 631-634.GUO Chun-yan, ZHANG Bang-jun. Method of drawing airport noise contours on computer based on Surfer[J].China Environmental Science, 2003, 23(6): 631-634.

[8] BERNARD M S, VAN P, BAARSMA B E. Using happiness surveys to value intangibles: the case of airport noise[J]. The Economic Journal, 2005, 115(500): 224-246.

[9] 孙桂茹, 马亮, 路登平, 等. 等值线生成与图形填充算法[J]. 天津大学学报, 2000, 33(6): 816-818.SUN Gui-ru, MA Liang, LU Deng-ping, et al. Investigation on the algorithm of making and filling isoline[J]. Journal of Tianjin University, 2000, 33(6): 816-818.

[10] 余明辉, 万远扬, 余飞. 一种绘制等值线图的新方法[J].武汉大学学报(工学版), 2006, 39(3): 52-54.YU Ming-hui, WAN Yuan-yang, YU Fei. A new method of drawing isoline map[J]. Engineering Journal of Wuhan University, 2006, 39(3): 52-54.

[11] 张显全, 刘忠平. 基于格网模型的等高线算法[J]. 计算机科学, 2005, 32(9): 199-201.ZHANG Xian-quan, LIU Zhong-ping. An algorithm of contour lines based on regular grid[J]. Computer Science,2005, 32(9): 199-201.

[12] JONES N L, KENNARD M J, ZUNDEL A K. Fast algorithm for generating sorted contour strings[J].Computers and Geosciences, 2000, 26(7): 831-837.

[13] CHANG F, CHEN C J, LU C J. A linear time component labeling algorithm using contour tracing technique[J].Computer Vision and Image Understanding, 2004, 93(2):206-220.

[14] ANGELOVA D, MIHAYLOVA L. Contour extraction from ultrasound images viewed as a tracking problem[C]//Proceedings of the 12th International Conference on Information Fusion. [S.l.]: [s.n.], 2009: 284-291.