一种形状多级描述方法及在多尺度空间数据几何相似性度量中的应用

安晓亚，孙 群，肖 强，严薇，2

1.信息工程大学测绘学院，河南郑州450052；2.61512部队，北京100088

1 引 言

空间数据几何相似性度量模型被广泛应用于地理信息科学的众多领域，根据应用，将其分为两类：第一类用于度量不同目标之间的几何相似性，主要应用于空间数据匹配［1－6］，其实质是从指定数据集中找到一个与待匹配对象在空间位置、方向、大小、长度和形状上最相似的目标，最终识别不同数据集上的同名实体［4］；第二类主要应用于同一目标在数据处理前后变形程度的度量［7－8］。上述应用中，多尺度几何形状相似性度量难度最大，其关键是形状描述子对形状的整体和局部特征都要能刻画，即满足紧致性，如第一类模型要求描述子尽量对形状的整体特征进行描述，而第二类模型要求描述子必须顾及尺度的变化，对局部细节特征更要考虑。目前有代表性的形状描述方法及其特点是：简单描述子对复杂形状的描述明显不足，中心距离对角度函数仅能刻画全局特征，单一的傅里叶描述对噪声敏感，小波描述严重依赖于形状起始点，形状上下文描述仅能刻画局部特征［7－11］。文献［1—6］基于第一类模型进行数据匹配，但很难达到相似度量结果不受尺度变化影响，文献［4］提出以中心距离对弧长的函数来描述几何形状，在同比例尺数据匹配中效果较好。

对此，提出一种具备对形状多级描述能力的新方法，既能以较详细的程度刻画形状，同时通过调节相关参数也能以较概略的程度刻画形状，可以达到通用于两类模型的目标。基于该方法，建立通用多尺度相似性度量模型，并将它们分别应用于多尺度数据匹配和图形化简前后的相似性度量中。

2 形状的多级描述方法

将形状轮廓表示为一组有序点集：A＝｛Pi＝（xi，yi）｜i＝1，2，…，N｝（图1（a）），O为几何中心，从Pi出发，沿逆时针将A按弧长等分为K个弧段，K为偶数，Gt（t＝1，2，…，K－1）是对应等分点。连接Pi与Gt，得到K－1条弦｛PiG1，PiG2，…，PiGK－1｝，lt（Pi）表示Pi对应第t条弦PiGt的弦长。s（Pi）表示O到Pi的中心距离。当Pi在轮廓线上改变时，每一个Pi对应一个lt（Pi）和s（Pi），故lt（Pi）和s（Pi）是Pi的函数。设P0为起始点，则Pi可表示为弧的长度ci。以ci为自变量，lt（Pi）和s（Pi）可表示为ci的函数Lt（ci）和S（ci）。至此，A可由自变量弧长ci、因变量弦长Lt（ci）和S（ci）来描述。对于线目标（图1（b）），如果lAB／cAB≥0.5，寻找一点O，使△OAB为等边三角形，然后以O为几何中心，计算s（Pi），在计算lt（Pi）时，将△OAB看作面目标，然后依据图1（a）的方法计算；如果lAB／cAB＜0.5，直接连接首末点A、B，将其转化为面目标按图1（a）的方法计算lt（Pi）和s（Pi）。

图1 多级弦长Fig.1 Multilevel chord length

一个ci对应弦长函数L1（ci），L2（ci），…，LK－1（ci）有K－1个，其中，ci∈［0，C］，C为轮廓线周长。将有序集合F（ci）＝｛L1（ci），L2（ci），…，LK－1（ci）｝称之为ci对应的多级弦长函数，Lt（ci）表示ci对应的第t级弦长函数。当t＝K／2时，LK／2（ci）是轮廓线1／2弧长对应的弦长，当t大于K／2时，无法将轮廓线整数等分，因此，将F（ci）压缩一半为F（ci）＝｛L1（ci），L2（ci），…，LK／2（ci）｝，只包含K／2级弦长函数。

上述描述方法满足旋转和平移不变性，用C对ci进行归一，然后以S（ci）和Lt（ci）的均值对它们进行归一化后可满足缩放不变性。

3 相似性度量模型建立及应用

3.1 利用多级描述方法度量形状相似性

多级弦长函数描述形状满足紧致性。以图2为例，计算两图当K＝8时的4级弦长函数，并绘制曲线进行对比，图3为1～4级弦长曲线。计算当横坐标取同一数值（从0开始到1，间隔为0.05）时每幅图对应两曲线纵坐标的差值，求这些差值序列对应的标准差s，分别为0.16、0.11、0.08和0.05，标准差越大，说明两形状越不相似，反之则否。经第一级弦长函数描述的两形状最不相似，说明对细节信息能很好地刻画，第2、3级弦长函数的相似性逐渐增大，到第4级弦长函数，两曲线已非常接近，说明对整体形状相似性能有效地描述。因此，可以利用多级描述方法来度量不同尺度空间数据的形状整体或者局部相似性。

图2 整体相似而细节有较大差异的两个形状Fig.2 Two shapes have whole similarity and local difference

图3 图2中两形状的4级弦长曲线Fig.3 4levels chord length curve for Fig.2

之所以可以采用上述方法利用标准差直接度量两形状的相似性，是因为两形状的起始点位置基本一致，且采样相同的轮廓线点数20。但在实际应用中，两形状起始点位置不一定一致，且轮廓线点数也不相等，导致每一级弦长函数的个数不一致，无法进行比较。为此，通过以下方法来改进：

（1）对轮廓线以等弧长间隔进行重采样N′个点c0，c1，…，cN′－1近似表达原始轮廓线，其中，N′＝2n，n是满足2n＞N的最小整数，然后计算Lt（ci）。

（2）对每一个Lt（ci）进行快速傅里叶变换，其形式为

以｜zt（m）｜描述第t级弦长函数，｜zt（0）｜＝1，取N′个系数中的前M 个，构造向量μt＝［｜zt（1）｜｜zt（2）｜…｜zt（M）｜］代替第t级弦长函数描述形状，得到有序集μ＝［μ1μ2…μK／2］Τ。以相同方法对S（ci）进行傅里叶变换，得到系数向量s＝［｜S（1）｜｜S（2）｜…｜S（M）｜］，于是通过μ和s描述整个形状。因此，即使两轮廓线上的点数不一致，但因为对所有弦长和中心距离进行傅里叶变换，在度量相似性时，两个形状都取变换系数中的前M个进行比较，故可以解决点数不一致的问题，同时，μ和s独立于轮廓线的起始点（证明略）。

至此，得到满足旋转、平移、缩放不变性，且满足紧致性和独立于轮廓线起始点的形状描述子μ和s。故可以用μ和s度量两个形状之间的形状相似度：设形状A和形状B的μ、s分别为μA＝和A和B之间的形状差异度定义为

｜sA－sB｜、｜μAi－μBi｜均指向量之间的欧氏距离，权系数w1、w2之和为1，μ、s均做了归一化处理，故d（A，B）shape∈［0，1］，则形状相似度sim（A，B）shape＝1－d（A，B）shape，且sim（A，B）shape∈［0，1］。因为中心距离函数可以刻画整体形状特征［8］，故当K值较小时，ds和dμ均可度量两形状的整体差异性，当K值较大时，dμ可度量细节差异性。

3.2 面向两类应用的几何相似度量模型建立

无论是第一类还是第二类模型，都要度量两目标对应形状的几何相似性。因此，可以将这两类模型统一于一个模型，仅在应用的过程中调节相关参数即可。设待度量的两目标为A和B，度量其几何相似性的总体思路与文献［4］基本一致，都是从两目标空间位置的邻近程度、形状、方向和大小（或长度）相似性程度综合考虑并加权求和，但文献［4］未考虑方向相似程度，且本文提出计算空间位置邻近程度及形状相似性程度均与文献［4］不同。设［a1a2a3a4］Τ、［b1b2b3b4］Τ分别为描述A和B的向量，各分量代表位置、方向、大小（或长度）和形状分量。d（A，B）表示A和B的差异度，sim（A，B）表示A和B的相似度，则A、B的几何相似性可通过下式来计算

采用加权的Minkowski距离来度量d（A，B），则式（3）转化为

式中，p值取2；ωj为权系数，且bj｜（j＝1，2，3，4）是经过归一化后的值，其含义是A和B中第j个分量的差异度，故｜aj－bj｜∈［0，1］，d（A，B）∈［0，1］，sim（A，B）∈［0，1］。

对空间位置的邻近程度度量，文献［4］提出的方法仅能度量面目标，不能度量线目标。对此，利用Hausdorff距离来度量线目标或面目标之间的邻近距离。但传统Hausdorff距离易受目标局部变形程度影响，为此，基于高斯概率统计模型改进Hausdorff距离。基本思路是：先计算A上的每个点与B上所有点的最小距离集合Df（A，B）；然后计算B上的每个点与A上所有点的最小距离集合Db（B，A），最后根据拉依达3σ准则，剔除局部变形较大点，具体步骤为：

（1）计算Df（A，B）和Db（B，A）及其均值和均方差

（2）对Df（A，B）中的任意元素dft，t＝1，2，…，m，m为A上的点数，如果它满足不等式

则认为此元素为A的内点，放入D′f（A，B）中，否则舍弃；对Db（B，A）存在类似操作，形成D′b（B，A）。

（3）令dfmax＝max｛D′f（A，B）｝，dbmax＝max｛D′b（B，A）｝，则A和B的Hausdorff距离为

设Hkmax（A，B）为A和B点集中两点最大距离，则第一个分量｜a1－b1｜＝Hk（A，B）／Hkmax（A，B）是经过归一化处理后的Hausdorff距离值。

空间目标的方向以MABR（minimum area bounding rectangle）的长轴方向角代替，故第二个分量方向差异度为｜a2－b2｜＝｜θA－θB｜／π，θA、θB分别为A、B的方向角；第三个分量面积或周长差异度为｜a3－b3｜＝｜SA－SB｜／max（SA，SB），SA、SB分别为A、B的面积或周长；第四个分量形状差异分量采用3.1节式（2）的计算方法计算，即｜a4－b4｜＝dshape（A，B）。将上述四个分量带入式（4）即可得到目标A、B的总的几何相似度。

利用上述几何相似度量模型进行空间目标匹配的基本步骤是：设A为待匹配的目标，利用式（4）逐一度量A与指定数据集D中的所有目标的几何相似度，然后从D中找到一个与A最相似的目标C，且当sim（A，C）的值大于给定的阈值sim0时，认为A与C匹配，否则不匹配。为提高匹配效率，可以缩小D的范围，方法是在A的最小外接矩形的横向和纵向加一个两数据集几何位置偏移最大Hausdorff距离maxDist，然后形成扩展的最小外接矩形EMBR（enlarge minimum bounding rectangle），处于EMBR内的目标集合构成集合D。

利用几何相似度量模型度量同一目标在化简前后变形程度的基本思路是：设A为原始目标，以不同的化简阈值对A进行化简后得到的新目标序列为B1，B2，…，Bn，然后利用式（4）分别计算sim（A，Bi）（i＝1，2，..，n）的值，进而探寻不同化简算法对形状保持程度的规律。

3.3 相关参数及阈值的确定

匹配过程中的相关参数及阈值的确定：

（1）M和K的确定。因傅里叶变换的高频易受噪声影响，故M不能太高，根据文献［9］，M暂取为10，根据经验，K值暂取为8。后文还会对M，K的取值进行讨论。

（2）其余参数的确定。引入相关反馈技术来确定匹配过程中的各种阈值［12］，匹配前，先选取一定数量的正例样本（已经确定匹配的目标），分别计算ds、dμ、｜aj－bj｜（j＝1，2，3，4）的值，并求上述每个分量的标准差，取倒数并归一化，分别得到式（2）和式（4）中的权值w1、w2和ωj（j＝1，2，3，4），最后根据式（4）计算每对正例样本的相似值，并计算其均值标准差σ，则

图形化简前后相似度计算过程中相关参数及阈值的确定：M的值与匹配过程一致，由于对形状的细节特征更要顾及，因此K值要更大。ωj、w1和w2的取值可根据对目标变形的关注程度调整。

4 实例及分析

4.1 多尺度几何匹配试验

匹配对应两种数据源分别为：现势性截止于2005年的4幅1∶25万数据，现势性截止于2009年的1幅1∶50万数据。经坐标系统一后将两种数据源叠加在一起，以1∶50万水域要素匹配1∶25万水域要素。

4.1.1 面状水库及岛屿的匹配

经统计，1∶25万数据中共有面状水库102个，岛屿78个，1∶50万数据中共有面状水库67个，岛屿49个。匹配的基本过程是：

（1）选取10对正例样本，就可计算出ds、dμ、｜aj－bj｜，如表1所示。

表1 样本相似特征分量及总相似值Tab.1 The swatch similarity characteristic components and general similarity value

得到ds和dμ对应标准差分别为0.017 5和0.069 1，则w1和w2的值分别为0.797 7和0.202 3。得到｜aj－bj｜（j＝1，2，3，4）的标准差分别为0.024 8、0.013 1、0.071 5、0.022 9，则ωj的值分别为0.231 4、0.437 3、0.080 2、0.251 0。得到sim的标准差为0.022 7，均值为0.950 4，则sim0的值为0.882 3。

（2）将上述阈值代入式（4），然后根据3.2节匹配面状水库和岛屿。记n1为匹配结果集中所有实体对的数目；r为n1中正确的匹配数目；n2为两种数据源中同名实体的对数。定义查准率P＝r／n1，查全率R＝r／n2。本例中，水库漏匹配8个，故n1＝67－8＋49＝108，在n1中没有误匹配，故r＝108，查看相似值小于sim0的实体集，发现漏匹配的8个水库中，有3个在1∶25万数据中为常年湖，另外3个是因为形状差异太大（图4（c）），2个是因为在1∶25万数据没有相应的实体。故n2＝67－3＋49＝111，则P＝100%，R＝97.4%。图4是从结果集中选取的岛屿匹配的正例1个（图4（a）），水库匹配的正例1个（图4（b）），水库漏匹配的负例1个（图4（c）），实线包围区域为1∶50万数据，阴影部分为1∶25万数据。

图4 面目标的匹配Fig.4 Area matching

A1对应搜索对象集合包括B1和B2；A2对应搜索对象集合也包括B1和B2；A3仅包括B3；A4仅包括B4。表2是图4的匹配结果。

表2 图4的匹配结果Tab.2 The results of Fig.4

4.1.2 线状河流的匹配

选取10对已匹配样本，得到w1和w2的值分别为0.660 2和0.339 8。ωj的值分别为0.302 5、0.219 4、0.191 0、0.287 1，sim0的值为0.882 4。图5是两种数据叠加在一起的显示效果，实线为1∶50万数据，虚线为1∶25万数据，图的左边部分是局部放大图。先根据河流代码将1∶25万数据进行合并，合并的方法是先进行节点匹配，然后逐段合并1∶25万数据中河流代码相同的目标，例如，B3是由P1P2和P2P3段合并而成。经统计，该区域共有1∶50万河流58条，其匹配结果为：55条正确匹配，2条漏匹配，1条是因为1∶25万数据无此目标，故P＝100%，R＝96.5%。从图5中选取一个正例（左上图）和负例（左下图），A1对应搜索对象集合包括：B1、B2、B3和B4，A2仅包括B5，表3是具体的匹配结果。

图5 线目标的匹配Fig.5 Line matching

表3 图5的匹配结果Tab.3 The results of Fig.5

4.1.3 算法比较

比较本文与文献［2］所提出的利用缓冲区增长方法而不用几何相似性度量方法在匹配的查准率、查全率和效率上的不同，同时与文献［4］所提出的形状相似度量方法进行比较。所利用的数据与前文所述一致，分别比较线目标和面目标匹配结果，具体结果如表4所示。通过分析，文献［2］与本文相比，查准率一样，即没有错误匹配出现，但查全率降低，即出现漏匹配。漏匹配的原因是当两线目标空间位置相邻较远时，文献［2］提出的方法不能匹配这两目标，而利用本文方法，即使两线目标距离较远，只要总体形状具有很大的相似性，仍然可以匹配。从匹配速度上看，文献［2］的方法较低，因为要构建线目标缓冲区，所以耗时较长；文献［4］与本文相比，查准率仍然一样，查全率也降低，原因是文献［4］提出的形状相似度量方法适用于同比例尺目标的匹配，当度量不同比例尺上的同一目标时，相似值差别会很大，因此可能导致漏匹配，如图4（b）所示情况，匹配速度与本文基本相同。

表4 匹配算法的比较Tab.4 The comparison of matching algorithms

4.1.4 K和M值对匹配结果的影响分析

仍以上述数据为试验数据，比较当K分别取2、4、8、16和32和M取1～20时查准率P和查全率R的变化情况，如图6所示。

图6 K和M值对匹配性能的影响Fig.6 The matching results for different Kand M

由图6可知，当K＝8时，无论对面目标或线目标的查全率和查准率，都能达到最大。当M取 9、10、11、12时线目标匹配查准率和查全率能达到最大。

4.2 图形化简前后的相似性度量

以1∶25万水域要素为原始数据，利用Douglas－Peucker（简称D－P）和Li－Openshaw（简称L－O）对目标以不同的阈值化简，然后利用式（4）度量化简前后的相似值，探寻相似值随化简阈值的变化情况。取参数M＝10，K＝16，w1＝w2＝0.5，ωj的值为0.25，表示对化简前后目标在空间位置、方向、大小和形状上变化的关注程度一致。D－P和L－O方法对应化简阈值d的初始值为0.05mm，步长也为0.05mm，最大值为1mm，对图7所示实心水库进行化简，分别得到20个化简结果和相似值。图7是化简阈值分别为0.55mm和1mm时两种化简方法对应的结果。

图7 两种化简方法的对比Fig.7 The comparison of two methods of simplification

图8（a）是总相似值sim（A，Bi）随d而变化的曲线，图8（b）是形状总相似值1－｜a4－b4｜随d而变化的曲线，根据3.1节，1－ds即为整体形状相似值，图8（c）是1－ds随d而变化的曲线。当K＝16时，本文提出的形状多级描述方法可以度量局部相似值，1－dμ是局部形状相似值，图8（d）是1－dμ随d而变化的曲线。分析图8，可得出以下结论：

（1）当0＜d≤0.55时，D－P和L－O对应总相似值具有基本相同的变化规律，用最小二乘线性拟合该区间的相似值得到sim（A，Bi）D－P＝－0.088 1d＋0.997 4，sim（A，Bi）L－O＝－0.070 3d＋0.991 3；当0.55＜d≤1时，D－P对应总相似值sim（A，Bi）随d而变化的变化率要明显比L－O大，用最小二乘线性拟合该区间的相似值得到sim（A，Bi）D－P＝－0.060 8d＋0.982 0，sim（A，Bi）L－O＝－0.003d＋0.959 8。

（2）与总相似值sim（A，Bi）的变化规律类似，形状总相似值、整体形状相似值和局部形状相似值在0＜d≤0.55时，D－P和L－O具有基本相同的变化规律，但当0.55＜d≤1时，L－O对应形状变化比D－P更具有稳定性。由此说明当化简阈值在较小范围变化时，D－P和L－O对形状的保持基本一致，但当化简阈值超出这一范围时，D－P算法对化简前后形状的保持要比L－O算法差。

图8 相似值随化简阈值变化的曲线Fig.8 The relationships between similarity value and simplification threshold

5 总 结

本文主要工作及创新点包括：

（1）组合多级弦长函数和中心距离函数描述空间数据几何形状，达到形状描述的紧致性和唯一性，该方法的最大优点是通过调节K值可以从整体到详细逐级描述形状，满足多尺度相似性度量要求和两类相似性度量模型的应用需求；通过多级弦长和中心距离的离散傅里叶变换，解决两形状轮廓线点数不一致和起始点不一致的问题。

（2）基于形状的多级描述方法建立顾及形状、方向、长度和空间位置的第一类相似性度量模型，并将其应用于不同比例尺空间数据的几何匹配，在此过程中，利用高斯概率统计模型改进传统的Hausdorff距离，提高抗噪水平；引入信息检索中相关反馈技术解决模型中阈值的确定问题。线目标和面目标的匹配试验表明，该模型可有效实现不同比例尺之间水域要素的匹配。

（3）基于第一类相似性度量模型建立了第二类相似度量模型，然后比较两类经典的化简算法Douglas－Peucker算法和Li－Openshaw算法在化简前后形状相似值随化简阈值变化的规律。

［1］ DUEKHAM M，WORBOY F.An Algebraic Approach to Automated Information Fusion［J］.International Journal of Geographical Information Science，2005，19（5）：537－557.

［2］ ZHANG M，SHI W，MENG L Q.A Generic Matching Algorithm for Line Networks of Different Resolutions［C］∥Proceedings of the 8th ICA Workgroup on Generalization and Multiple Representations.LA Coruna：ICA，2005：1－8.

［3］ GOMBOSI M，ZALIK B，KRIVOGRAD S.Comparing Two Sets of Polygons［J］.International Journal of Geographical Information Science，2003，17（5）：431－443.

［4］ HAO Yanling，TANG Wenjing，ZHAO Yuxin，et al.Areal Feature Matching Algorithm Based on Spatial Similarity［J］.Acta Geodaetica et Cartographica Sinica，2008，37（4）：204－209.（郝燕玲，唐文静，赵玉新，等.基于空间相似性的面实体匹配算法研究［J］.测绘学报，2008，37（4）：204－209.）

［5］ TONG Xiaohua，DENG Susu，SHI Wenzhong.A Probabilistic Theory－based Matching Method［J］.Acta Geodaetica et Cartographica Sinica，2007，36（2）：210－217.（童小华，邓愫愫，史文中.基于概率的地图实体匹配方法［J］.测绘学报，2007，36（2）：210－217.）

［6］ FU Z L，WU J H.Entity Matching in Vector Spatial Data［C］∥The International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences.Beijing：ISPRS，2008：146－147.

［7］ BIAN Lihua，YAN Haowen，LIU Jiping，et al.An Approach to the Calculation of Similarity Degree of a Polygon before and after Simplification［J］.Science of Surveying and Mapping，2008，33（6）：207－208.（边丽华，闫浩文，刘纪平，等.多边形化简前后相似度计算的一种方法［J］.测绘科学，2008，33（6）：207－208.）

［8］ TANG Luliang，LI Qingquan，YANG Bisheng.Shape Similarity Measuring for Multi－resolution Transmission of Spatial Datasets over the Internet［J］.Acta Geodaetica et Cartographica Sinica，2009，38（4）：336－340.（唐炉亮，李清泉，杨必胜.空间数据网络多分辨率传输的几何图形相似性度量［J］.测绘学报，2009，38（4）：336－340.）

［9］ ZHANG D S，LU G J.Study and Evaluation of Different Fourier Methods for Image Retrieval［J］.Image and Vision Computing，2005，23（1）：33－49.

［10］ SEO C K，TAE J K.Texture Classification and Segmentation Using Wavelet Packet Frame and Gaussian Mixture Model［J］.Pattern Recognition，2007，40（4）：1207－1221.

［11］ CHEN Shi，MA Tianjun，HUANG Wanhong，et al.Gait Recognition Based on Shape Context Descriptor［J］.Pattern Recognition and Artificial Intelligence，2007，20（6）：794－799.（陈实，马天骏，黄万红，等.基于形状上下文描述子的步态识别［J］.模式识别与人工智能，2007，20（6）：794－799.）

［12］ WU Hong，LU Hanqing，MA Songde.A Survey of Relevance Feedback Techniques in Content－Based Image Retrieval［J］.Chinese Journal of Computers，2005，28（12）：1969－1979.（吴洪，卢汉清，马颂德.基于内容图像检索中相关反馈技术的回顾［J］.计算机学报，2005，28（12）：1969－1979.）