王永东，韩尚炜，谢小菊

（1.武汉中迪联创科技有限公司，湖北 武汉 430074；2.中国地质大学（武汉） 地理与信息工程学院， 湖北 武汉 430074）

等值线是由某一数量指标值相等的各点连成的平滑曲线，能反映该指标在一定空间或时间范围内的分布规律，被广泛应用于各种图件中。等值线的绘制一般分为离散数据格网化、等值线追踪和等值线光滑 3个步骤，其中离散数据格网化是将数据以格网的形式重新组织，为等值线追踪做准备，常见的格网类型包括不规则三角网和规则格网；等值线追踪是针对某一指标，利用格网找到若干等值点，并根据格网拓扑关系确定点的连接顺序，从而生成折线型等值线；等值线光滑是对折线型等值线进行光滑处理，使等值线更加符合实际、更加美观。目前，等值线光滑方法主要包括手工勾绘和计算机自动生成两类，手工勾绘方法的图形美观、精度较高，但费时费力，一般多用于复杂地形图的出版；计算机自动生成主要依靠光滑曲线生成算法，通过在已知等值点间添加曲线点的方式实现等值线的光滑。目前的光滑曲线生成算法可以分为两类：①过点的光滑算法，即生成的曲线严格经过等值点，以张力样条函数法为代表[1]，由于曲线必须过 等值点，当等值点较密集时，易产生震荡现象，与实际地形特征不符合；②不过点的光滑算法，即光滑过程中曲线会经过等值线附近，但不一定严格经过等值点[2]，以Bezier函数法为代表[3-4]。

针对Bezier函数法存在的问题，本文采用在拐角处以等值点和切线点共同作为控制点的方法使光滑曲线经过所有等值点[5]，并对比分析了五点法、张力样条函数法、改进后的Bezier函数法在不同类型地形图中的适用性。

1 常用的等值线光滑算法

1.1 五点法

五点法的基本原理为：在两个相邻数据点之间建立一个三次多项式曲线方程，且要求整条曲线具有连续的一阶导数以保证曲线光滑，每个点的一阶导数由该点及其相邻两侧各两点（共5个点）来确定。对于离散点序列中相邻两点(xi,yi)、(xi+1,yi+1)，经过这两点的光滑曲线可表示为：

为保证曲线光滑，要求从(xi-1,yi-1)到(xi,yi)以及从(xi,yi)到(xi+1,yi+1)这两个区间上的曲线在(xi,yi)点处具有相等的一阶导数。假定已知(xi,yi)和(xi+1,yi+1)处的导数，结合已知点(xi,yi)、(xi+1,yi+1)，代入式（1）可求解参数c0、c1、c2、c3，进而得到从离散点(xi,yi)到(xi+1,yi+1)的曲线方程。

若已知曲线上5个点以及曲线经过各点的顺序，可由Akima五点法估算曲线在第三点处的导数值[6]。设平面上5个点为Pi(xi,yi)，i=1,2,3,4,5，则每相邻两点间折线段的斜率为：

第三点处斜率为：

若等值线为开曲线，则首尾处各有两点无法找齐计算所需的5个点，此时采用补点的方法进行处理。具体方法为：对于等值线末端最后3个点(x1,y1)、(x2,y2)、(x3,y3)，补充点(x4,y4)、(x5,y5)满足式（4），且5个点同属于式（5）所表示的抛物线上。

通过式（4）、式（5）可求解抛物线方程，进而求解补点坐标。通过补点可计算得到等值线首尾4个点处的一阶导数值。

1.2 张力样条函数法

1.2.1 张力样条函数

设平面上一组已知点为(x1,y1)、(x2,y2)、…、(xn,yn)，x1＜x2＜…＜xn，常数σ≠0为张力系数，求解一个二阶导数连续的三次样条函数应满足：

式（7）是一个二阶非齐次的常系数线性微分方程，令

求解微分方程可得：

式（9）为通过所有等值点的张力样条函数，式中 未知项为各等值点处的二阶导数值，因此若能确定 各等值点处的二阶导数值，则可确定张力函数。

1.2.2 节点关系式与端点条件

对式（9）进行微分，再结合式（10），可得到节点关系式如式（11）所示。

节点关系式可简记为：

式（12）是含有n个未知量，i=2,3,…,n−1的n-2个方程的线性方程组，要唯一定解，需附加两个方程式，由于等值线会出现开曲线和闭曲线两种情况，需要分开讨论。

对于开曲线，首尾端点处则有：

由式（12）和式（14）组成三对角线性非奇异方程组，可求得唯一解，再代入式（9）可求得张力样条 函数方程。

对于闭曲线，首尾端点处则有：

可得到一个三对角线并带有右上角和左下角系数矩阵的线性方程组，它也是非奇异的，求得唯一解后代入式（9）求得张力样条函数[7]。

1.2.3 张力系数选择

在地形图上采用等高线描绘地形时，特征点之间的曲线取最短为优，为避免相邻光滑曲线相交且要通过等值点，采用规范化的张力系数σ作为标准值，即

σ为规范化张力系数，需要预先进行试验，取最适合当前地图的张力系数，当发生变化时，则由标准值重新计算张力系数[8-9]。

1.3 Bezier函数法

1.3.1 Bezier曲线插值

Bezier曲线是一条需要在起点和终点之间构建插值多项式的混合函数的不规则曲线，可拟合任意数目的控制点，其曲线方程次数比控制点个数少1。通过n+1个顶点定义一个n次多项式，假设n+1个点的位置矢量为Pi，i=0,1n，则Bezier曲线上各点的坐标插值为[10-11]：

式中，Pi为原特征多边形的等值点；Bi,n(t)为n次Bezier 多项式。

Bezier函数法的缺陷是光滑后的曲线不再通过所有等值点，为解决该问题，本文首先对曲线进行分段拟合，每段取相邻3个等值点，拟合得到的曲线将不再通过中间的等值点，再利用切线抹角法进行处理，使光滑后的曲线通过所有等值点。

1.3.2 切线抹角法

如图1a所示，对于折线B-A-C，利用Bezier函数法生成的光滑曲线不再经过A点，为解决该问题，先过点A作线段P1P2，A为P1P2中点（图1b），再以P1、A、P2共同作为切线抹角法的控制点，即可得到经过 点A的光滑曲线，如图1c所示。

图1 切线抹角法原理图

为避免出现多余拐点，切线P1P2应垂直于∠A的平分线，设∠A的平分线与BC相交于D，，，进而求得切线的斜率为：可计算得到

曲线的松紧程度由P1P2两点间的距离决定，如 图2所示，延长DA到E点，取AB与AC的中点F、G，连接FE、GE与切点交于P1、P2，建立系数控制曲线松紧程度。

图2 切线端点计算

设dt已知，则点E坐标的计算公式为：

基于点F、G、E的坐标以及直线P1P2方程，通过直线求交可计算得到点P1与P2的坐标值[5,12]。

2 算法编程实现与结果分析

本文主要针对地形图中生成的等值线（即等高线）结果进行分析，测试数据来源于对一些具有代表性地形的实际测量数据。

2.1 五点法测试

为实现五点法，首先在程序中获取等值线的条数以及每条等值线的点数，然后根据曲线的开闭决定是否需要在首尾各插补两个点，最后循环建立所有相邻两点的三次曲线方程。

五点法测试结果如表1所示，可以看出，五点法的光滑效果主要受等值点密度和等高线密度的影响，当等高距较小、等高线较密集且等值点较稀疏时，光滑处理后的等高线易出现相交的情况。其原因在于如果相邻两条等高线走向相背，进行光滑处理时，生成的曲线会向二者中间凸起，若等值点分布稀疏，则会导致凸起部分相交（图3a），中部的两条等高线向二者中间凸起，若继续缩小等高距，则会导致两条等高线相交（图3b）。

表1 五点法测试结果

图3 不同等高距下五点法结果（局部放大图）

2.2 张力样条函数法测试

张力样条函数法首先采用经验值作为张力系数，再根据等值线的开闭，调用不同的计算函数生成等值线。若光滑效果不佳，则调整张力系数重新进行等值线光滑。在计算等值线的过程中，涉及带有右上角和左下角的三对角线系数矩阵的线性方程组求解时，可采用追赶法（Doolittle分解）将矩阵分解为一个特殊的下三角矩阵和一个特殊的上三角矩阵的内积，再进行求解。

张力样条函数法的光滑效果可通过张力系数进行调节。结果表明，当张力系数增大时，曲线会更加贴近原始折线；当张力系数减小时，曲线弯曲程度会增大。考虑到实际应用场景，当等高距较小、等值点密度不够时，应适当增大张力系数以避免等值线相交，当等高距较大、等值点密度有余时，应适当减小张力系数使地形特征更加明显。张力样条函数法也可通过改变张力系数来调整等值线，有效解决了五点法中等值线相交的问题，但该方法首先需要提供经验值作为张力系数，且需要多次调节张力系数使效果达到最优，因此操作复杂、使用门槛较高。

2.3 改进后的Bezier函数法测试

采用改进后的Bezier函数法生成光滑曲线需分段进行，每次取相邻3个点坐标，首先计算切线点P1、P2的坐标，然后利用3个等值点和2个切线点生成Bezier曲线。若生成曲线效果欠优，可调整相关系数dt后重新进行等值线光滑。

与张力样条函数法类似，改进后的Bezier函数法利用dt来控制曲线与原始折线的偏离程度，当dt趋近于1时，曲线趋向原始折线；当dt较小时，曲线会偏离原始折线。因此，当等高线较密集时，应适当增大dt以避免等值线相交；当等高线较稀疏时，应适当减小dt使地形特征更加明显。与张力样条函数法类似，改进后的Bezier函数法也可通过调节系数值来使最终效果达到最优。同时，相较于张力样条函数法，Bezier曲线不会产生多余拐点，光滑效果更优。

3 算法特点和适用性

五点法的计算速度快，但可能出现等值线相交等情况，适用于原始数据质量高、数据分布平缓、等值线稀疏、结果要求不严格的场合，操作简单，可一键生成；张力样条函数法的等值线严格正确，但需手动调节张力系数，适用于结果要求严格、结果需逐步调整以达到最优状态的场合；改进后的Bezier函数法的等值线严格正确、不产生多余拐点，但需手动调节相关参数，适用于结果要求严格、结果需逐步调整以达到最优状态、不产生多余拐点的场合。