点群选取是制图综合的重要内容，目的是在点数减少的情况下，尽可能正确地表达点群的空间信息；难点是在语义信息缺乏时，如何仅仅依赖点群的几何信息进行选取。对此，许多学者进行了研究^[1-5]，基本思想大都是依据几何信息构造一定的模型，对点的重要性进行判断，从而对点群进行取舍。在这些方法中，点群的聚类及Voronoi图是综合过程的重要手段，在制图综合中具有广泛的应用^[6-7]。然而现有方法在对聚类算法本身进行深入考虑方面较为欠缺。如文献[4]通过k-means算法提取聚类中心，并用层次Voronoi图逐步细化地表达聚类中心点，实现点群的综合，为点群的多尺度表达提供了一种有效的思路，但存在以下两个方面的问题：一是k-means算法需要事先确定k值及人为指定初始聚类中心，且k-means算法将每个点指派到距离最近的聚类中心，因而不能检测非球面分布的点群数据；二是单纯的聚类只考虑数据本身的分布特点，若用同一层级的聚类中心对点群进行粗略表达，不能顾及点群选取必须满足的密度对比要求。因此，该方法不能实现点群的全自动选取，选取结果也存在一定的不确定性和不合理性。

针对上述不足，本文选用更加合理的聚类算法实现点群的自动聚类，进而提出层次聚类的概念，在此基础上提出一种顾及密度对比的点群选取方法。

1 基于改进密度峰值聚类算法的点群自动聚类

聚类是根据数据对象及其关系对数据进行分组，从而分析数据的潜在结构，判断数据的自然簇属性并压缩数据存储容量的数据分析方法^[8-10]。现有的聚类方法主要分为基于划分的聚类^[11]、基于层次的聚类^[12]、基于密度的聚类^{[9, 13]}、基于网格的聚类^[14]及基于模型的聚类^[15]等，其中，基于密度的聚类方法以高密度区域作为判断依据。这种非参数方法与传统方法相比，不仅适用于处理任何形状的数据集，而且无需预先设定簇的数量，更易于聚类的自动实现^[9]。因此，本文基于Rodriguez等^[16]提出的一种密度峰值聚类算法，实现空间点群的自动聚类。

1.1 密度峰值聚类算法及分析

密度峰值聚类算法原理是基于聚类中心的两个特点提出的，一是聚类中心本身密度大；二是聚类中心与其他密度更大的点之间的距离相对更大。其主要过程是，对于数据集中的每个点，计算局部密度和相对距离，其中表示数据集中到的距离小于或等于的个数，称为截断距离（cutoff distance），需要用户自行确定。文献[16]利用cut-off核（cut-off kernel）计算：

$ {\rho _i} = \sum\limits_{j \in {I_S}\backslash \left\{ i \right\}} {\chi \left( {{d_{ij}} - {d_c}} \right)} $

(1)

式中，d_ij表示点p_i与p_j间的欧氏距离；I_S={1, 2…n}（下同），j∈I_S\{i}表示j∈I_S且j≠i；函数χ(x)计算方法为：

$ \chi \left( x \right) = \left\{ \begin{array}{l} 1,x < 0\\ 0,x \ge 0 \end{array} \right. $

(2)

cut-off核计算的局部密度为离散值，有学者提出利用高斯核（Gaussian kernel）计算局部密度^[17]，得到连续值，减小聚类产生冲突的概率，计算公式为：

$ {\rho _i} = \sum\limits_{j \in {I_S}\backslash \left\{ i \right\}} {{{\rm{e}}^{ - {{\left( {{d_{ij}}/{d_c}} \right)}^2}}}} $

(3)

相对距离δ_i的计算公式为：

$ {\delta _i} = \left\{ \begin{array}{l} \mathop {\min }\limits_{j \in I_S^i} \left\{ {{d_{ij}}} \right\},I_S^i \ne \emptyset \\ \mathop {\max }\limits_{j \in I_S^i} \left\{ {{d_{ij}}} \right\},I_S^i = \emptyset \end{array} \right. $

(4)

式中，I_Sⁱ为指标集，

$ I_S^i = \left\{ {k \in {I_S}:{\rho _k} > {\rho _i}} \right\} $

(5)

计算出每个点的ρ_i和δ_i后，以二元对(ρ_i, δ_i)的形式绘制于二维坐标系中，称之为决策图，如图 1所示。从决策图中选取横纵坐标都比较大的点作为密度峰值点，即聚类中心。根据聚类中心及密度边界阈值，将剩余点分到各簇中，完成聚类过程。

密度峰值聚类算法虽有诸多优点，但将该方法应用于空间点群的聚类还存在两个缺点。一是需要凭经验选取截断距离，降低了算法的科学性和普适性；二是需要人工根据决策图确定聚类中心，对于同样的决策图，不同的人可能得出不同结果。

如图 1（b）中，短虚线框和长虚线框分别选取了不同数量的密度峰值点，则会得到完全不同的聚类结果。

1.2 密度峰值聚类算法的改进 1.2.1 截断距离自适应确定

从式（1）、式（3）可以看出，ρ_i由d_c决定，在文献[16]中，d_c是根据若干数据集的经验值设定的，其选取的标准是使每个数据点的平均邻居个数约为数据点总数的2%。但是对于空间数据，不同数据集在点数及分布模式上都存在差异，因此更科学的方法是根据数据集本身的特点自适应地设置d_c的大小。

文献[18]指出数据的势能与数据的不纯度之间存在一定的关系, 即数据的势能分布较均匀时，数据的不纯度较大；数据的势能分布不均匀时，数据的不纯度较小。数据的不纯度最小时，数据的势能差别最大，更易于聚类。对于空间点群集合P={p_i}_i=1ⁿ，每个点的势能计算^[18]公式为：

$ \varepsilon \left( x \right) = \sum\limits_{i = 1}^n {{{\rm{e}}^{ - {{\left( {{d_{xi}}/\sigma } \right)}^2}}}} $

(6)

式中，σ为影响因子。

数据的不纯度可以由基尼指数计算^[18]：

$ G = 1 - \sum\limits_{i = 1}^n {{{\left( {{\varepsilon _i}/Z} \right)}^2}} $

(7)

式中，$Z = \sum\limits_{i = 1}^n {{\varepsilon _i}} $，为数据域的总势能；显然G是σ的函数，即G=G(σ)。

对比点的势能计算式（6）和局部密度计算式（3）发现二者完全等价。因此，为了计算合理的局部密度，最优的截断距离d_c等价于取得最优的影响因子σ。

根据上述分析，以基尼指数值最小时的影响因子作为截断距离的值，能够达到最优聚类的效果，即：

$ {d_c} = \arg \left( {\min G\left( \sigma \right)} \right) $

(8)

1.2.2 聚类中心的自动确定

为了形象说明聚类中心的确定方法，考虑图 2（a）中的例子，其中包含32个离散点，将二元对{(ρ_i, δ_i)}_i=1³²在二维直角坐标系中画出，如图 2（b）所示。

从图 2（b）中可以看出，点20、12、3由于同时具有较大的局部密度ρ和相对距离δ，位于坐标系的右上方，并从其他数据点中脱颖而出，图 2（a）中这3个点正好可能是原始点群的3个聚类中心。

为了自动探测聚类中心，用ρ和δ归一化的乘积评估点间的差异度，根据差异度的统计特征和变化规律确定聚类中心。定义簇中心权值为：

$ {\gamma _i} = \frac{{{\rho _i} - {\rho _{\min }}}}{{{\rho _{\max }} - {\rho _{\min }}}} \cdot \frac{{{\delta _i} - {\delta _{\min }}}}{{{\delta _{\max }} - {\delta _{\min }}}},i \in {I_S} $

(9)

式中，ρ_min、ρ_max、δ_min、δ_max分别是ρ_i、δ_i的最小值、最大值。

将γ_i从大到小降序排列，以序号N为横坐标，γ_i为纵坐标，将排序结果在二维直角坐标系中画出，如图 3所示。从图 3中可以看出，γ总体呈下降趋势，但是下降程度各阶段有所不同：聚类中心的γ值下降比较剧烈，非聚类中心的γ值下降比较平滑，而从聚类中心到非聚类中心过渡时γ值有一个明显的跳跃。为了确定曲线下降的拐点，用两点连线的斜率表示簇中心权值的下降趋势，即：

$ k_i^m = \left( {{\delta _{i + m}} - {\delta _i}} \right)/m $

(10)

式中，k_i^m表示在区间[i, i+m]上簇中心权值的平均变化率，该参数描述了某一区间内γ的总体变化趋势，则γ值变化的拐点满足：

$ x = \arg \left( {\max \left( {k_1^{i - 1}/k_i^{n - i}} \right)} \right) $

(11)

根据式（10），k₁^i-1表示第1个点到第i个点的斜率，k_i^n-i表示第i个点到n个点的斜率。因此，拐点可以认为是簇中心权值总体变化最快的临界点。在二维坐标系中绘出坐标对(i, k₁^i-1/k_i^n-i)，如图 4所示。显而易见，序号为4的点为拐点，对应图 3中点号为17的点为拐点，从而自动探测出聚类中心点为点20、12、3。

2 顾及密度对比的多层次聚类点群综合 2.1 点群多层次聚类

定义1   多层次聚类。对点群$P = \left\{ {{p_i}} \right\}_{i = 1}^n$的聚类结果P=$\left\{ {{\rm{cluster}}\left( {{c_1}} \right),{\rm{cluster}}\left( {{c_2}} \right) \cdots {\rm{cluster}}\left( {{c_k}} \right)} \right\}$，分别进行再次聚类，其任意一个点簇${\rm{cluster}}\left( {{c_j}} \right)\left( {j = 1,2 \cdots k} \right)$的聚类结果为${\rm{cluster}}\left( {{c_j}} \right) = \left\{ {{\rm{cluster}}\left( {c_j^1} \right),{\rm{cluster}}\left( {c_j^2} \right) \cdots {\rm{cluster}}\left( {c_j^{k'}} \right)} \right\}$，其中${\rm{cluster}}\left( {c_j^l} \right)\left( {j = 1,2 \cdots k,l = 1,2 \cdots k'} \right)$表示点簇，c_j^l为各自聚类中心。重复上述步骤，直到所有簇为单个点的过程，称为点群的多层次聚类。从定义1可以看出，多层次聚类是将上一层级的一个簇作为下一层级的原始聚类点群，进行重复迭代聚类，目的是将点群划分成不同等级的簇，并确定不同等级的聚类中心。点群的多层次聚类结果可以存储为一种层次树结构，且由于点群分布的不均匀性，这种层次树结构表现为非平衡多叉树。将上述算例中的数据进行多层次迭代聚类，得出如图 5所示的层次树结构，其中灰色背景标注的点为各级别簇的聚类中心。

2.2 顾及密度对比的点群选取方法

为了使点群选取满足密度对比要求，本文采用以下选取策略。一是根据各级别簇中点数从上而下分配选取的点数循环，直到分配到各簇的选取数量为1；二是只选取聚类中心。通过上述两条策略，能够根据选取点数定位到选取的点。具体步骤如下：

1）对点群进行多层次迭代聚类，生成点群层次树；

2）根据制图综合需要，依据开方根定律计算需要选取的数量N_select；

3）在第1次聚类的结果上，按照各簇点数比例分配需要选取的数量N_select；

4）以此类推，逐层往下分配需要选取的数量，直到需要选取的数量均分配到聚类中心点为止。

按照上述步骤，即可获得满足数量和密度对比要求的选取结果。需要说明的是，在步骤3）中将选取点数分配到各簇时，由于会出现选取数量为非整数的情况，为了简化分配过程，可按照四舍五入取整的方法获得各类簇需要选取的数量，最终选取数量与开方根定律计算的结果可能会有少量偏差。当然也可以根据小数点后面的数字大小分配开方根计算的选取数量，得到严格满足开方根定律的选取结果。

以第1层级聚类为例，设簇数量为k，各簇点数为n_j(j=1, 2…k)，则各簇应选取的数量为：

$ {N_j} = {\mathop{\rm int}} \left( {{N_{{\rm{select}}}} \cdot {n_j}/n} \right) $

(12)

式中，int(·)表示四舍五入取整。同理对于其他级别点簇，可得到相应的选取数量。

继续上述算例，假设需要选取18个点，为了便于描述，用二元对(L-N)表示簇的编号，其中L=1, 2, 3…表示聚类的层级，N=1, 2, 3…表示簇的序号。则第1次聚类结果中，簇（1-1）、（1-2）、（1-3）应选取点数为4、6、8；将3个簇中应选取的点数向下一级别簇分配，以簇（1-1）为例，簇（2-1）、（2-2）、（2-3）应选取点数为2、1、1；以此类推，簇（1-1）最终选取结果点编号为3、2、4、8。同理，原始点群的最终选取结果点编号为3、2、4、8、9、11、12、15、31、25、19、17、20、21、23、30、27、29，如图 6所示。

3 实验与分析 3.1 不同分布模式点群选取实验

实验由点群的自动聚类、层次树构建及点群的自动选取3部分组成，整个过程在Visual Studio 2010开发环境中基于ArcGIS Engine 10.2开发实现。采用Vmap1格式数据，比例尺为1:25万，从可以公开下载的数据中选取美国3个地区不同分布模式的3类数据，分别是聚集分布的工业抽象点数据（对应图 7中结果1~3）、均匀分布的高程点数据和随机分布的城市点数据，具体见表 1。分析发现属性项中没有完整表征3类数据中各点重要性的等级信息，因而不能简单地依据属性信息进行选取。运用本文方法，分别进行1:50万、1:100万、1:400万制图综合操作，选取点数见表 2，原始点群及选取结果如图 7所示（未按比例尺绘制）。从综合结果来看，不同数量的选取结果形成了从详细到概略的多层次表达，3种分布类型的数据都能按照指定的数量获得相应的选取结果，且均能较好地保持原始点群的空间分布形态。为了对选取结果进行定量评价，利用3类数据初次聚类中心对图幅进行Voronoi剖分，分别划分出14、27、16个综合区，按图 8所示对综合区进行编号。统计各多边形综合区内各比例尺数据点数，以综合区编号为横轴、所在综合区点数为纵轴，绘制3类数据各综合区各次选取点数折线图，如图 9所示。

从图 9可以看出，各综合区不同比例尺综合结果与原始点群在数量上基本保持一定的比例关系。对于各比例尺的综合结果，密集的地方依然密集，稀疏的地方依然稀疏，说明本文方法较好地保持了密度对比要求。从图 9中还可看出，原始点群、1:50万、1:100万及1:400万综合结果4条折线没有出现交叉，表明利用本文方法进行点群选取没有出现违反密度对比要求的情况。

3.2 与现有方法的比较分析

按照文献[19]中基于信息传递的点群综合评价方法，本文方法直接考虑了点群的统计信息（按开方根定律确定选取数量）和度量信息（保持密度对比），间接考虑了点群的拓扑信息（只选取聚类中心点），未直接考虑专题信息（语义信息），可见本文方法主要针对语义信息缺乏的情况。在语义不完备的情况下，即等级信息不完整或者部分要素没有语义信息时，可将所有要素一起聚类，按照本文方法构建多层次树，对于等级高的要素，无论其在层次树的哪个等级，直接标注选取，然后在各级别簇中选取剩余数量的点，实现选取。本文相比文献[4]的优势主要是可选取任意数量的点，且满足制图综合中的密度对比要求。

4 结语

本文提出了一种顾及密度对比的多层次聚类点群选取方法，通过利用聚类中心表达相应层级的点簇，保证选取结果与原始点群空间分布形态上的整体一致；根据点群层次树结构逐层向下分配待选取点数，保证满足密度对比要求。实验验证了本文方法的可行性和普适性。本文方法的优点一是无需依靠点群的属性信息即可实现点群的多尺度表达；二是选取数量完全可控。本文的点群选取方法本质上只利用了点群间的距离关系，如何综合考虑点群间的方向关系和拓扑关系，并通过聚类的方法实现更科学合理的选取是下一步需要研究的内容。