等值线的追踪是对离散数据用数学方法插值，将具有相同量值的点变换成图的过程，它将三维信息显示于二维平面，是进行地理要素空间特征分析的强大工具，通过绘制相应的等值线图简单而直观地进行综合分析^[1]。现有的等值线追踪算法在进行出口判断时，会依据最近原则进行出口选择，该原则既能保证追踪算法的顺利运行又能保证算法效率最优^{[2, 3]}。

国内外现有的等值线追踪包括：格网线性内插^[4]和高次曲面内插^[5]。其中，格网线性内插法又分为三角网格法^[6-8]和矩形网格法，两者都是通过网格边实现线性内插，再将得到的等值点进行追踪生成等值线^[9]。

高次曲面内插分为局部曲面和整体曲面内插^[10]。高次曲面内插法的空间变量坐标不易匹配特定的函数表达式，故编程程序复杂且曲面难以契合空间实体^[11-13]。

格网线性内插法和高次曲面内插法都是在格网建好后，在格网边上内插等值点进行等值线追踪^[14]，但是当4条边都存在等值点时，就会出现追踪出口的二义性问题，而现有的常规追踪算法和单个网格内部的改进算法^[15]都无法正确寻找到出口。从而导致在山谷与山脊同时出现的混合地貌下，例如马鞍面地貌，会出现与真实情况严重不符的整体地貌被割裂的错误现象。尤其是在地势凹凸走向非常明显，或者插值网格较粗且取值间距较大的情况下，更容易导致异常凸起或凹陷的等值线小圆出现。

针对现有理论方法上存在的不足，本文提出一种基于8点法的等值线追踪算法，通过外向扩张来探索整体地势的走向，由此选择出正确的等值线追踪出口，避免了整体地貌割裂的错误现象，最终使得等值线可以符合并还原原始真实地貌。

1 构建网格区域

本文的8点追踪方法，首先需要根据等值点构建规则网格区域。

遍历所有的样本点数据，根据每个样本点在网格区域中的坐标和样本数据的位置、网格距，得出样本点在网格区域中的坐标编号，使用行遍历模式，根据式(1)计算所有节点的插值^[16]。

$ u = \sum\limits_{i = 1}^n {{a_i}} {u_i} $

(1)

式中，n表示样本的个数；a_i为距离加权；u_i为样本值。距离权重为网格点到样本点距离平方的倒数，设d_i为距离，则权重系数为^[16]：

$ {a_i} = \frac{{{{\left( {1/{d_i}} \right)}^2}}}{{\sum\limits_{i = 1}^n {{{\left( {1/{d_i}} \right)}^2}} }} $

(2)

然后，利用8点追踪方法在构建好的规则网格中进行等值线追踪，计算两个格网点中间是否有等值点。设追踪值为Z₀，则m、n两点之间等值点存在的判断公式为：

$ \left( {{Z_m} - {Z_0}} \right)\left( {{Z_n} - {Z_0}} \right) \le 0 $

(3)

2 等值线追踪 2.1 现有等值线的追踪方法

等值线追踪进入网格后，只能从网格的另外3个边方向追踪出去。如果不在算法上给予处理，就会出现等值线的交叉和不确定现象，如图 1所示。

现有的研究技术或文献中的解决方案为：假设前一等值点为a₁，当前等值点为a₂，要追踪的下一等值点为a₃，显然可能出现3种情况(a₃₁、a₃₂、a₃₃)，如图 2所示，而实际追踪时只能选择一种。则选择的次序为：

1) 当a₃₁、a₃₃都存在时，选择靠近网格左侧边者为a₃；

2) 当a₃₁、a₃₃靠近左侧边距离相等时，则选择与a₂距离近者为a₃；

3) 当a₃₁、a₃₃只有一个存在时，存在点作为a₃；

4) 当a₃₁、a₃₃都不存在时，对边必存在a₃₂作为a₃。

但是该方法会存在错误，如图 3所示。

根据最近法则，会形成图 3中虚线箭头指向的等值线(正确形式应该为图 3中实线箭头)，而该等值线将跨越山脊、山谷线，从而导致地貌割裂现象，等值线绘制结果不符合地貌特征。

2.2 基于8点法的等值线追踪算法

如图 4所示，以从左边进入为例，中心格的角上的点称为内4点，●表示高于当前高程称为高程点，○表示低于当前高程称为低程点，图形上面c-1、c、c+1、c+2是列号，左面r-1、r、r+1、r+2是行号。

除中心格4个顶点外，还应该考虑虚线圆标出的8个点(外8点)，判断属于哪种情况。用H_{i, j}表示i行j列点的高程，算法流程如下：

1) 当外8点的平均高程大于中心格4角中的黑点(即高程点平均值)时，即H_外>(1/2)(H_{r, c}+H_{r+1, c+1})，应认为是山谷。而当外8点的平均高程小于白点即中心格4角中的低程点平均值时，即H_外 < (1/2)(H_{r+1, c}+H_{r, c+1})时，应认为是山脊。按照等值线不穿过山谷、山脊规则即可判断出口。

2) 当H外介于两者之间，即外8点的平均高程大于内4点中低程点平均值，且小于内4点中高程点的平均值时，可认为中心格类似马鞍面。仍然用前高程值与中心点比较的方法，用待定系数法可知，中心点表达式为：

$ \begin{array}{*{20}{c}} {{H_中心} = A \times \left( {{H_{r, c}} + {H_{r + 1, c}} + 1 + {H_{r + 1, c}} + } \right.}\\ {\left. {{H_{r, c + 1}}} \right)/4 + B \times {H_外}} \end{array} $

(4)

3) 要满足两高程点连接的对角线所示条件，由图 4可知，必有：

若

$ \begin{array}{*{20}{c}} {\left( {{H_{\rm{外}}} = (1/2)\left( {{H_{r, c}} + {H_{r + 1, c + 1}}} \right), } \right.}\\ {{H_中心} = \left( {{H_{r + 1, c}} + {H_{r, c + 1}}} \right)/2} \end{array} $

(5)

若

$ \begin{array}{*{20}{c}} {\left( {{H_{\rm{外}}} = (1/2)\left( {{H_{r + 1, c}} + {H_{r, c + 1}}} \right), } \right.}\\ {{H_{\rm{中心}}} = \left( {{H_{r, c}} + {H_{r + 1, c + 1}}} \right)/2} \end{array} $

(6)

4) 将式(5)、式(6)代入式(4)得A=2, B=-1，最后将A、B代入式(4)得：

$ \begin{array}{*{20}{c}} {{H_{\rm{中心}}} = 2 \times \left( {{H_{r, c}} + {H_{r + 1, c + 1}} + {H_{r + 1, c}} + } \right.}\\ {\left. {{H_{r, c + 1}}} \right)/4 - {H_{\rm{外}}} = 2{H_{\rm{内}}} - {H_{\rm{外}}}} \end{array} $

(7)

式中，H_内表示内4点的平均高程。

5) 求出H_中心后，可与当前高程比较，若当前高程大于H_中心，应向左转；若当前高程小于H_中心，应向右转。

3 验证与对比

为了验证本文的8点法等值线绘制方法，使用文献[16]中涉及的位于中国第一、二阶梯过渡段四川省大渡河流域的数字高程模型数据作为实验数据，当地地形地貌复杂、山高谷深、河网密布，并同时具有丘陵、平坝、沟谷等特征，是理想的实验区域。使用常规等值线追踪方法，绘制的高程值1 800 m等值线如图 5所示。

从图 5中圈出的部分可以看出，等值线受当地地貌连续起伏影响，会形成割裂原始真实地貌的“等值线小圈”现象，且无法绘制出当前整个流域的地形走势，其根本原因即等值线追踪时以最近距离为闭合依据所导致。对比图 6，在采用本文的8点法等值线追踪算法后，整个流域的整体地貌得到了较好的还原。

4 结束语

本文融入了现实中的地貌特征影响因素，在等值线遇到出口选择情况时，进行了8点法扩张判断，在判断出当前区域地貌特征后，再进行追踪的出口取舍。从而避免了等值线追踪到错误出口，使得等值线不割裂整体地貌，不在山谷、山脊地区出现错误小圆，最终等值线可还原真实地貌特征，绘制的等值线图更利于简单而直观地进行综合分析。在水文模型的建立、地形模拟、公路铁路勘测和规划、自然资源管理等地学领域得到了广泛应用，而且显示出了许多优越性，本算法简单明了，便于用计算机进行处理，能够很容易地计算等高线、坡向坡度，自动提取地形。