一种基于网格重心的等值线追踪算法

王俊骄 陈 晴 张姗姗

(1.杭州市气象局，浙江 杭州 310008;2.浙江省气象信息网络中心，浙江 杭州 310017)

0 引 言

在气象要素图形显示中经常用到等值线图，它能较直观地给出气象要素的分布情况。不管是网格化的数值预报产品分析，还是在离散站点数据网格化应用中，等值线图的应用都非常广泛。等值线图是以相等数值点的连线表示连续分布且逐渐变化数量特征的一种图型，等值线生成常用算法有标准网格法和三角网格法两种[1]。而标准网格法利用区域内的高程离散点,建立逼近的曲面模型,然后通过内插得到一组规则的矩形网格数据,再在矩形网格的基础上线性插值追踪得到各条等值线。

目前，国内对标准网格等值线追踪的算法很多，如按一定方向直接对每个网格各边的等值点进行等值线追踪[2]、通过对角线转轴法进行追踪[3]、利用数据关联表生成矩形网格等值线[4]等等，对追踪出来的等值线再用一些插值拟合方法进行平滑处理,如Bezier、B样条拟合等[5]，或为了使等值线图更直观,还需要对等值线进行颜色填充[7]。

以上等值线追踪算法的设计均旨在提高整体的搜索效率，提升速度，在搜索过程中只考虑小网格的进边和出边，未分析等值线在小网格内的走向，在小网格尺寸很大时尤其明显。因此，本文提出一种基于网格重心的等值线追踪算法，在等值线追踪过程中，引入重心点分析每个网格，通过增加等值点序列中点的数量，不仅可更真实地反应等值线整体逐渐变化的特征，而且为后续的等值线光滑提供更多的控制点，同时由于本文算法在判断出边时仅需通过重心判断，避免了网格出口边判断复杂的问题。

1 等值线的原理及算法

对于标准网格的等值线生成，首先需获得网格数据，在气象应用中，网格数据一般通过内插方法将离散数据转换成网格数据，也可从气象应用中常用的Micaps四类格式、Grib格式和Netcdf文件获取。

在整个大网格中，分析每个小网格获取等值点，再判断走向，将它们依次用线连接起来，最后光滑等值线，从而构成等值线。本文重点讨论分析网格数据等值线追踪[8]。

1.1 等值点的判断

在几何中，对于由点A和点B构成的一条边AB，其中XA，YA分别表示点A的横坐标和纵坐标，XB，YB分别表示点B的横坐标和纵坐标，VA，VB分别代表点A和点B的高程值。判断在该边上是否存在满足高程值为W的等值点R的条件如下：

(VA-W)(VB-W)<0

(1)

若边AB满足(1)式，则通过内插可知等值点R的横坐标XR和纵坐标YR：

(2)

假设规则网格是由m×n个网格数据点构成,用Di,j表示纵向第i行，横向第j列的网格点，则用Vi,j表示该点的高程值，Di,jDi,j+1表示由点Di,j和Di,j+1所构成的横向边，Di,jDi+1,j表示由点Di,j和Di+1,j所构成的纵向边，其中1≤i≤m，1≤j≤n。

定义1对于任意一条边Di,jDi+1,j，若满足(1)式，则边Di,jDi+1,j定义为等值边，而等值点R的横向位置x和纵向位置y：

(3)

同理对于任意一条横边Di,jDi,j+1为等值边，可得到等值点相应的位置。

定理1由等值点的满足条件(1)式可知，在任意一个小网格中，若某边为等值边，即存在等值点，则该网格的另外3条边中至少有1条边同时为等值边。

另外，在网格数据等值线追踪过程中，须对每次搜索到的等值点进行判别，分析该点是否为终点，以停止搜索。而根据等值线特性可知，在搜索某高程值等值线时，终点需满足以下条件：该点为边界上的点或者该点为搜索起始点。

1.2 传统网格追踪

一条等值线的传统追踪算法描述如下：首先，依照自上而下，自左向右的顺序搜索每个小网格，将第一个找到等值边上的等值点作为搜索起始点，以该等值边为进边分析该网格，计算该网格上另外一个等值点和相应的等值边，连接两个等值点，即完成第一个网格的追踪。然后，以新的等值边作为进边，得到下一个网格，依此类推，直至等值点为终点。

然而，在传统追踪过程中，当网格内不止存在一条出边时(见图1)，就需判断等值线的走向。不仅降低了追踪的效率，而且可能由于不同的判断方法导致产生不同的等值线。当网格尺寸较大时，虽然能加快等值线搜索速度，但无法更真实地反映网格内部和整体的变化特征。

图1 传统追踪示意图(括号中代表该点高程值)

2 本文算法

2.1 基于小网格重心的追踪

本文在传统的追踪算法基础上，提出利用网格的重心来控制等值线的走向，通过重心引导，避免了等值线走向问题，同时使等值线在网格内的走向更符合等值线的实际特征(见图2中b，c，d)。

图2 等值点追踪时小网格内部的走向

在一个小网格中，已知四个点A,B,C,D，括号内为该点的高程值。则重心点E的横坐标Ex、纵坐标Ey和高程值VE为：

基于小网格重心的追踪算法可描述为：当某一边进入小网格时，分析以该网格重心与进入边构成的三角形。由于在该三角形中另外两条边中有且仅有一条为等值边，则找到该等值边及相应的等值点。继续分析下一个三角形，直到找到的等值点在该网格的四条边上，即结束该小网格的搜索。一般该追踪算法可出现图2中b，c，d 3种情况。

以图2中的c为例，追踪高程值为5的等值线。首先得到网格ABCD的重心E，以边AB作为进入边，则分析三角形ABE得到新的等值边AE，由于AE不是网格的边，则继续分析三角形AED，CDE，直到找到网格的出边CD，才完成整个小网格的追踪。在该网格的追踪中，共增加了2个新点，将这2个点存入到整条等值线点的序列中，可为后续的光滑，增加新控制点。

2.2 算法

等值线追踪到的所有等值点坐标按顺序连接就形成了等值线。然而该线条为一条折线，故需对这些等值点进行曲线平滑处理，将处理后的坐标连接起来即可绘出连续平滑曲线。曲线平滑常用的方法有Bezier方法、B样条方法、三次样条方法及最小二乘法，同时对等值线进行填充颜色，边界裁剪等后续加工处理。

为避免在追踪过程中，对某些等值点重复追踪，造成算法冗余，在追踪顺序时，本文先依序搜索网格最外面4条边，根据等值线要么自封闭,要么边界封闭的特性，先将边界封闭的等值线追踪完毕，再追踪自封闭的等值线。因此，本文的算法流程图如下。

图3 本文算法的流程框图

3 本文算法的应用

为验证本文算法的可用性，选取2016年4月5日08时欧洲中心细网格的小时降水预报，在数值预报中绘制高程值为8的等值线，为能清楚地展示网格内部等值线追踪的效果，截取了中心10×10网格数据,用于展示等值线追逐效果，并通过传统和本文的算法进行比较(结果见图4和图5)。

图4 传统追踪的示范图

图5 本文基于重心追踪的示范图

由图4和图5可知：采用传统的追踪算法比基于重心追踪的等值点少，所以等值线会相对比较粗糙，基于重心的追踪法更能反映等值线在小网格内的走向。传统的追踪方法因为需对网格出口边判断，从而导致整体生成的等值线会不同，而等于重心追踪在小网格中不涉及出口边判断问题，确保了等值线的唯一性，可以应用到实际业务中。

4 结 语

本文在分析传统矩形网格等值线追踪算法的基础上,引入小网格重心的算法加以改进,实现寻点与追踪同步。与传统方法比较,本文的方法有以下改进：1)可更真实地反应等值线整体逐渐变化的特征；2)在等值线点序列中增加了节点，为后续的等值线光滑提供更多的控制点；3)在追踪过程中，不涉及出口边判断问题，确保了等值线的唯一性；4)可利用同样原理扩展到三角网格的等值线追踪中。

本文算法通过编码，已成功应用到气象数值预报产品、Micaps四类格式的资料处理中，并应用到实际业务中。

[1] 林毅,金烨，马登哲等.正规网格等值线的虚路径扫描算法[J].计算机工程及应用，2001,37(13):92-94.

[2] 余卫东,李湘阁,王靖.气象场等值线自动绘制[J].气象科技，2002,30(4):222-225.

[3] Dayhoff M O. A contour map program for X-ray crystallogra-phy[J]. Co mmunications of the Association for Computing Machinery,1963,6(10):620-62.

[4] 孙科峰,孙根正,李洁.一种新的矩形网格生成等值线算法[J].东华大学学报(自然科学版)，2005,31(4):66-69.

[5] 韩丽娜,张红祥,张群会.Bezier曲线修改的一种分割算法[J].计算机工程与科学，2006(7):77-79.

[6] 韩家新,王家华.一种等值线填充的连通区域搜索算法[J].计算机应用与软件，2003(10)：5-6.

[7] 陈正旭,封秀燕,王亚云.多岛屿地图上绘制气象要素等值线色块的自适应方法[J].气象科技，2009,37(3):356-359.

[8] 王亚强.等值线相关算法类库的开发与应用[J].气象科技，2010(4):478-482.

[9] 张红祥,车鹏飞.一种多段Bezier曲线光顺拟合方法[J].科技信息,2008(13):48-49.

[10] 孙桂茹,马亮,路登平等.等值线生成与图形填充算法[J].天津大学学报,2000,33(6):816-818.