建筑物要素是基础地理信息的重要组成部分，其形状相似性度量对建筑物空间模式识别、分类、查询和制图综合具有重要意义^[1-6]。矢量几何图形相似性度量涉及两个关键问题：即几何形状的定量描述方法及形状之间差异性(或相似性)的度量模型^[7]。

几何图形的定量描述方法可以分成基于图形区域、基于图形结构特征和基于图形轮廓特征3种类别^[8]。基于图形区域的方法以图形的整体指标来反映图形结构，通常使用定量的形状参数符(如面积、延展度等)来表征形状的整体结构。典型的方法有网格描述符^[9]和基于矩的形状描述方法^[10]，前者将网格用于描述形状特征的区域统计信息，适用于基于内容的图像表达；后者使用低阶矩构造出由7个不变量所构成的不变矩组，适用于矢量多边形的区域统计信息表达。此类方法主要基于图形区域的整体特征，侧重于反映形状的全局结构，会导致形状轮廓细节信息的丢失。基于图形结构特征的方法通常以图形的骨架特征作为形状的结构化表达，该方法主要应用于非刚性变化的形状相似性度量^[11]。文献[12]通过比较骨架端点之间的测地线距离来匹配骨架图。文献[13]使用不连续骨架的形状相似性度量方法解决传统骨架连接中的不稳定问题。然而，由于尺度变化，地图要素的形状骨架特征并不唯一，该方法主要适用于度量图像数据集(如动物形状)和非刚性变换形状之间的相似性。基于图形轮廓特征的方法，其基本思想是将形状节点转换为几何表达式，以字符串或者函数的形式来近似拟合形状边界。常用方法有傅里叶描述子、链码法、转角函数法、形状上下文、轮廓点分布直方图等。其中，傅里叶描述子通过将空间域的坐标转换为频域的频率来实现形状的定量描述^[14]；链码法是一种基于边界的编码描述法，该方法通过形状起始点的坐标和边界点方向代码来描述图形^[15]；转角函数法将多边形形状转换为函数，在函数空间中分析结构并计算相似性^[16]；形状上下文以两个图形轮廓采样点之间的相互关系(如距离、梯度)作为形状匹配策略^[17]；轮廓点分布直方图在极坐标下对目标形状的轮廓特征点进行描述^[18]。此类方法强调邻近节点之间的相互关系，几何轮廓的局部噪声会影响形状的定量表达。综上，传统的几何形状定量描述方法存在以下不足：基于图形区域的方法以统计特征描述整体信息，容易忽略局部细节；基于图形结构特征的方法，其结构信息提取算法往往复杂度高且提取结果不唯一；基于图形轮廓的方法需要构建复杂的图形编码序列或者形状轮廓描述符，且多数方法并不具备图形的旋转、平移和缩放不变性，不能直接用于形状相似性度量。

在多边形形态特征的定量描述基础上，形状之间的差异通常用形状相似性或相异性(距离)表示，形状相似性越大，则形状之间的距离越短，两个形状越相似^[19]。建筑物的形状相似性度量方法往往基于图形的几何特征，主要包括简单几何特征法、编码法、形状描述符法等。简单几何特征法主要选择单个或多个指标，例如距离相似度(欧氏距离、Hausdorff距离^[20]等)、方向相似度^[21](余弦法、方向关系矩阵^[22]等)、几何相似度^[23](面积比、周长比、面轮廓线相似性^[24]等)，或者使用加权平均法进行综合相似性度量。编码法主要是对比计算编码相似度或构建相似度描述函数^[15]，其计算量相对较小，但对复杂图形的识别准确率不高。形状描述法则通过形状描述子之间的差异来度量形状相似性。例如，傅里叶描述函数将图形的轮廓特征以傅里叶函数表示，通过计算不同阶次系数之间的差异(欧氏距离)来测量形状相似性距离^[25]。转角函数表示了形状边界上切角对弧长之间的关系，形状间的相似性以切角的差异来衡量^[26]。在以上方法中，几何特征法对于相似性指标的选取和权重的确定具有较大的主观性，编码法难以识别复杂的几何图形，而基于形状描述符的方法需进行复杂的傅里叶或转角函数变换，其计算性能较低，且存在傅里叶展开级数和转角起始点难以确定的问题。

针对目前几何图形定量描述方法和形状之间差异性度量模型存在的问题，本文提出了一种基于典型时间规整(canonical time warping，CTW)算法的建筑物形状相似性度量方法，该方法具有以下特点：①以建筑物矢量坐标序列作为形状的定量描述方法，可以直接根据矢量坐标序列进行形状相似性的计算，无须将问题拆分为几何形状的定量描述和相似性度量两个过程。因为直接利用矢量坐标序列，故可同时兼顾几何图形的整体特征与局部特征。②在相似性度量中，该方法可以对齐具有不同节点个数的建筑物形状，对于多尺度地图数据中建筑物图形的旋转、平移和缩放等操作不敏感，具有旋转、平移和缩放不变性，稳健性较强^[27]。③该方法基于降维思想，将二维矢量坐标序列降维至一维序列，以降维后一维序列间的距离作为建筑物形状之间的相似性距离，整个算法不需要人为控制相似性指标权重或阈值，可操作性强。

1 典型时间规整原理

CTW算法计算得到的距离被称为典型规整距离^[28]，该算法结合了典型相关分析(canonical correlation analysis，CCA)和动态时间规整(dynamic time warping，DTW)算法^[27]，主要应用于步态序列识别和多个序列的时间对齐^[29-30]。其中，DTW算法用于时间序列对齐，CCA算法实现多维序列进行特征提取和降维^[31]。

给定时间序列矩阵Q=[q₁, q₂, …, q_i, …, q_m]和C=[c₁, c₂, …, c_j, …, c_n]，其长度分别为m和n。首先构造m×n的距离矩阵，矩阵内的第(i, j)个元素表示q_i和c_j之间的匹配距离，标记为d(q_i, c_j)，通常以欧几里得距离表示^[31]，。之后构造一条从矩阵元素(1, 1)到矩阵元素(m, n)的弯曲路径P(P=p₁, p₂, …, p_k, …, p_K(max(m, n)≤K≤m+n-1))。求解DTW距离即为求解弯曲路径下累计距离和的最小值，即使得下列距离平方和的代价函数最小化

式中，弯曲路径P需要同时满足以下3个条件：①端点对齐。路径P从第一个点对开始，在最后一个点对结束，即p₁=(1, 1)，p_K=(m, n)；②连续性。路径上的任意两个相邻点p_k=(a, b)与p_k-1=(a′, b′)满足0≤|a-a′|≤1，0≤|b-b′|≤1；③单调性。若p_k=(a, b)与p_k-1=(a′, b′)为路径上前后两个点，则需满足a-a′≥0，b-b′≥0。K为满足上述条件后对齐两个时间序列所需的索引数。符合上述条件的弯曲路径可能有多条，其中最短的弯曲路径可以通过动态规划算法求得，公式如下

式中，γ(i, j)表示序列Q[1:i]与序列C[1:j]之间的DTW距离，也就是当前单元格和相邻元素的最小累计距离。

为了方便在DTW算法中添加特征选择机制从而实现DTW算法与CCA算法的结合，式(1)中的代价函数可以修改为如下形式

式中，W_q∈0, 1^K×m和W_c∈0, 1^K×n分别为对齐时间序列Q和C的规整矩阵，对齐后序列Q和C的长度一致。K为式(1)中对齐两个序列所需的索引数。

图 1展示了通过DTW算法进行序列对齐的基本流程。以一维时间序列Q=[4, 1, 3, 3, 5, 2]和C=[4, 1, 3, 5, 4, 5, 2]为例，序列Q的长度为6，序列C的长度为7。首先，构建两个序列之间的距离矩阵，由不同的弯曲路径得到累计距离矩阵，通过路径回溯得到最短弯曲路径(红色路径)。其次，根据路径长度，可知对齐两个序列所需的索引数K为8；根据路径走向(c₃对应于q₃和q₄；q₅对应于c₄、c₅和c₆)，可知序列的匹配关系及对齐序列Q和序列C的规整矩阵W_q和W_c，其中W_q为6×8的矩阵，W_c为7×8的矩阵。最后，将序列Q和C分别与两个规整矩阵相乘，得QW_q^T=[4, 1, 3, 3, 5, 5, 5, 2]，CW_c^T=[4, 1, 3, 3, 5, 4, 5, 2]，二者之间的欧氏距离即为原序列Q和C之间的DTW距离，经计算该距离为1。

经过DTW算法可以对齐两个时间序列，使得两个序列的长度达到一致。为了计算具有多维特征的时间序列之间的差异，CTW算法在DTW算法的基础上添加了线性变换(V_q^T, V_c^T)(CCA算法)，目的是对多维特征序列进行降维。通过求解一组线性变换矩阵，使得变换后的新序列之间的相关性最大，优化模型如下

序列Q和C标准化后，式(4)可以化简为。当V_q^TQQ^TV_q和V_c^TCC^TV_c由于正交约束而被替换为单位矩阵时，式(4)可以转换为如下带有约束条件的优化模型

式(5)可以根据奇异值分解求出原始序列Q和C的最佳线性变换系数V_q和V_c^[32]。由于V_q^T和V_c^T也可以理解为对序列Q和序列C的投影矩阵，因此这种降维方式对图形具有平移、旋转和缩放不变性。在此基础上，CTW算法通过式(6)实现DTW算法与CCA算法的结合

式中，W_q^T和W_c^T通过对齐序列Q和C使得两个序列的长度保持一致，V_q^T和V_c^T实现了多维特征序列在低维特征空间的映射。与CCA算法类似，式(6)需要添加以下约束：①QW_q^TI_K=0，CW_c^TI_K=0；②V_q^TQD_qQ^TV_q=V_c^TCD_cC^TV_c=I；③V_q^TQWC^TV_c为对角矩阵。其中，I_K是大小为K的单位矩阵，D_q=W_q^TW_q，D_c=W_c^TW_c，W=W_q^TW_c。整个算法的输入值为序列Q和C，输出值为V_q、V_c、W_q、W_c。初始化V_q^T=I，V_c^T=I，通过迭代交替计算V_q、V_c、W_q、W_c得到最终的算法收敛结果，该结果即为CTW距离。

2 基于CTW算法的建筑物形状相似性度量

地图中建筑物矢量多边形轮廓由一系列首尾相连的有序坐标点组成，这种有序坐标点串的表达方式与时间序列的表达方式相似。对于两个建筑物多边形形状相似性的比较，可以提取二者的轮廓坐标，将其归一化到相同的参考体系，然后计算二者坐标点位之间的累积距离差值。显然，差值越小则形状越相似，差值越大则形状差异性越大。在概念上，这一思想与时间序列相似性度量方法CTW的原理契合。地图中不同的建筑物多边形具有不同的位置、方向和大小，且多边形之间的节点数目也各不相同。CTW算法具有平移、旋转和缩放不变性，且能对点数各异的矢量建筑物坐标序列进行特征提取，并计算这些序列在特征空间的相似性距离，从而能够度量两个多边形形状之间的相似性。

在具体的应用中，CTW算法通过添加规整矩阵W_q^T和W_c^T对齐具有不同节点数的建筑物形状(DTW算法)，通过加入线性变换V_q^T和V_c^T提取建筑物矢量坐标的低维特征(CCA算法)，不断迭代交替计算线性变换系数和规整矩阵，得到最终的CTW距离。该距离的物理意义为两建筑物多边形轮廓差异程度的定量化描述，表达了两多边形形状的相似性，具有平移、旋转和缩放不变性。显然，CTW距离越小则形状之间的相似度越高。

在形状相似性度量的距离矩阵计算中，假设构成建筑物A的有序节点集合P₁^a=[p₁^a, p₂^a, …p_i^a, …p_n^a, p₁^a]，构成建筑物B的有序节点集合为P^b=[p₁^b, p₂^b, …p_j^b…p_m^b, p₁^b]，序列中的首尾节点坐标一致。其中p_i^a=(x_i^a, y_i^a), i∈n，(x_i^a, y_i^a)为建筑物A中节点p_i^a的坐标；p_j^b=(x_j^b, y_j^b), j∈m，(x_j^b, y_j^b)为建筑物B中节点p_j^b的坐标。在构建(n+1)·(m+1)的距离矩阵时，矩阵内第(i, j)个元素表示节点p_i^a和p_j^b之间的欧氏距离，标记为d(p_i^a, p_j^b)，公式如下

假设DTW算法对齐后的序列长度为K(式(1)中的索引数)，经过CTW算法收敛后，建筑物A可以表示为a=[a₁, a₂, …, a_K]，建筑物B可以表示为b=[b₁, b₂, …, b_K]，序列a和b之间的欧氏距离即为建筑物A和B之间的CTW距离

然而，矢量坐标序列的不同构建顺序及不同起始点的选择会导致不同的节点匹配结果，从而产生不同的CTW距离。本文参考了文献[26]的策略，统一以顺时针方向生成矢量坐标序列，采用移动边界起始点的方法求得所有序列之间CTW距离的最小值，作为最终的相似性计算结果。另外，建筑物形状中常常会出现具有孔洞结构的多边形。根据Gestalt原则，对于有孔洞的复杂多边形来说，人们对于其形状特征的认知以最外围轮廓为主，其内环对于多边形整体形状认知的影响不大。因此，本文不考虑图形内部的孔洞结构对形状认知的影响，统一以图形外围轮廓构成的矢量坐标序列进行相似性度量。

图 2展示了基于CTW算法进行建筑物形状相似性度量的流程，图中以建筑物形状A和形状B作为具体示例。首先，根据建筑物形状可得相应的矢量坐标序列，其中形状A具有4个节点，形状B具有5个节点，形状A和B的矢量坐标序列分别为：P^a=[p₁^a, p₂^a, p₃^a, p₄^a, p₁^a]，P^b=[p₁^b, p₂^b, p₃^b, p₄^b, p₅^b, p₁^b]。然后，初始化低维特征V_a^T=I，V_b^T=I，根据算法调整低维特征和规整矩阵，形状A和形状B经过DTW算法匹配可得到的各个节点之间的对应关系。其中，形状A的节点p₃^a对应于形状B的节点p₃^b和p₄^b，需在p₃^a处插入一个相同点使得序列P^a和P^b之间的长度一致。序列对齐后CCA算法降维结果表现为一一对应关系，即序列a=[a₁, a₂, a₃, a₄, a₅, a₆]，b=[b₁, b₂, b₃, b₄, b₅, b₆]。最后，通过迭代交替计算低维特征和规整矩阵直至算法收敛，可得建筑物形状A和B之间的CTW距离为0.001。这两个形状之间的DTW距离为11.95，远大于CTW距离，说明CTW算法更适合度量建筑物形状之间的相似性。

图 3是L型建筑物(图 3(a))和经过旋转(图 3(b))、平移(图 3(c))和放大(图 3(d))后的形状示例图。为了证实CTW算法的平移、旋转和缩放不变性，本文对比了形状a与形状b、c、d之间的CTW距离与DTW距离。由表 1可以看出，形状a与形状b、c、d之间DTW距离远大于CTW距离，且CTW距离几乎为0，说明平移、旋转和放大后的形状与原始形状非常接近。这一结果符合人类的视觉空间认知，证实了使用CTW算法计算形状相似性的可行性。

3 试验分析 3.1 试数据

深圳是1979年成立的中国第一个经济特区，现在是珠江三角洲的先驱城市^[33]。2019年，深圳城市化率达到99.52%(http://tjj.sz.gov.cn/)。因此，深圳的建筑形态和分布可以被视为中国大城市的典型代表。本次试验以深圳市1∶10 000的矢量建筑物分布图为参考数据。深圳的建筑物主要分布在道路沿线，建筑形态特征相对规则。

本文选出了7种典型的形状：矩形、L形、E形、十字形、C形、T形和Z形，并选取了420个类似的形状作为试验数据，其中每种类型分别包含60个建筑形状。在数据采集和存储的过程中，参考数据中可能会有一些异常节点，如冗余点和共线点^[34]，这些点需要在试验前删除。之后将这些建筑物图形缩放、平移至同一图幅，如图 4所示。另外，本文构建了7个基本形状模板来方便比较形状相似性(图 4中的红色形状)。

3.2 试验方法

试验使用Mathematica 12.1中的Canonical Warping Distance函数来计算420个建筑物形状与7个模板形状之间的相似性距离(CTW距离)^[28]，根据与模板最短距离确定建筑物形状识别结果。本文与快速傅里叶变换(fast Fourier transform，FFT)计算得到的建筑物形状相似性结果进行了对比，对比结果见表 2。其中，TP表示算法识别结果与人类视觉感知一致的建筑物数目，FP表示算法识别结果与人类视觉感知不一致的建筑物数目，FN表示算法识别错误的结果中与人类视觉感知一致的建筑物数目。CTW算法的时间复杂度为O(N³)，试验耗时约617 s(11.2 min)，FFT算法的展开级数为500，时间复杂度为O(Nlog₂N)，试验耗时约40 s。

由表 2可知，基于CTW算法的试验耗时较大，但是其查全率和查准率均高于98%，远高于FFT的识别准确率(52.857%)。同时，针对不同类型的建筑物形状，基于CTW算法的查全率和查准率均达到95%以上，且对于矩形、十字形和T型建筑物形状而言，查全率和查准率均为100%。这是因为CTW算法的本质是基于图形轮廓进行的形状相似性度量，其特点是同时考虑了图形的整体和局部特征，如果图形轮廓存在噪声节点，将会影响算法匹配结果，进而会降低算法的准确率。本次试验数据中的矩形、十字形和T型建筑物形状比较规整和对称，影响形状轮廓的噪声节点数目较少，因此基于CTW算法的准确度较高。对于FFT的算法而言，该算法在识别L型和T型建筑物时的查全率和查准率均低于50%，准确率较低。这是因为FFT算法无法有效顾及建筑物形状中多直角表达的形态特征，该算法相对更适合非多直角表达的复杂形状之间的相似性度量，如面状湖泊^[8]。试验结果表明，与FFT算法相比，CTW算法在建筑物形状相似性度量上具有更好的性能。

3.3 试验分析

针对试验数据中出现的有孔多边形，本文统一以图形外围轮廓构成的矢量坐标序列进行了相似性度量，这些形状在本次试验中都能得到正确的度量结果。表 3展示了本次试验数据中与人类认知结果不一致的图形(建筑物1—6)，及部分几何复杂度较高但是被本文算法正确识别的图形(建筑物7—10)相似度量结果。表 3中的视觉认知结果表示使用CTW算法获得的最相似形状与人类认知结果一致。对比建筑物1、2、3、6和建筑物9，其中建筑物1、2、3和9在人类视觉认知中都属于E型建筑物但不属于标准的E型。由于建筑物9的整体结构较规整，相对模板形状而言影响其形状轮廓的噪声节点较少，因此仍能得到正确的相似性度量结果。建筑物1、2、3和建筑物6，形状轮廓的噪声节点极大影响了形状整体结构，如形成了齿轮状(建筑物1)或圆弧状(建筑物6)的多余结构，这会影响图形与模板坐标序列节点的对齐，从而导致错误的相似性度量结果。建筑物4，该形状具有二义性，可识别成L型或C型，这说明CTW算法可以反映一些具有二义性的建筑物形状。对于建筑物5和建筑物7、8、10，虽然建筑物7、8、10的几何复杂度较高，噪声节点较多，但是这些噪声节点分布较集中，比如建筑物10的噪声节点只分布在L型两个直角拐弯处，该图形仍具有与传统形状类似的整体结构，因此仍能得到正确的相似性度量结果。相对而言，建筑物5噪声节点数量多且分布分散，且该形状与人类视觉认知中传统的Z型建筑物相差过大，因此造成了错误的度量结果。

表 4为本次试验结果与人类认知结果不一致的复杂多边形(表 3中的建筑物1—6)经Douglas-Peucker算法化简后的相似性度量结果。结果表明，复杂多边形经轮廓化简后的形状在CTW算法下均能得到正确的相似性度量结果。因此，针对数量较多且分布分散的噪声节点，可以考虑在相似性度量前使用传统的化简算法进行多余节点的抽稀，从而减少图形的噪声节点，提升算法的准确率。建筑物形状中的少量局部噪声节点一般不会影响建筑物的整体结构，对CTW算法的相似性度量结果影响较小。对于表 3中的所有复杂建筑物及表 4中化简后的建筑物1—5，使用FFT算法均不能得到正确的相似性度量结果，而CTW算法可以正确识别表 3中的建筑物7—10及表 4中的所有形状。

3.4 应用

形状检索是以人类认知的角度从GIS数据库中找到与模板形状匹配的目标图形，建筑物形状检索的关键问题在于形状相似性的度量。本文方法可有效应用建筑物形状检索。表 5记录了与每个模板形状最相似的前6个建筑物以及它们之间的CTW距离，可以看出，这些形状与模板形状在视觉上都较为相似，表明CTW算法可以有效地应用于形状检索场景。其中，对于矩形和十字型建筑物而言，模板形状与前6个形状在视觉上非常相似，其形状之间的CTW距离最小，且矩形形状的长宽比对形状相似性度量结果没有太大影响；对于C型和T型建筑物而言，其模板形状具有对称性，而第4个C型建筑物和第5个T型建筑物并不具备，表明多边形的矢量坐标序列对形状的整体特征具有较好的表示能力。此外，在Z型建筑物的检索结果中，前5个建筑物形状均为Z型模板形状旋转、缩放或平移后的图形，而第6个形状为Z型模板的镜像图形，说明CTW算法不仅具有旋转、缩放、平移不变性，还具有镜像不变性。

4 结论

在传统的度量建筑物形状相似性的算法中，往往需要构建复杂的图形编码序列或形状轮廓描述符。本文提出了一种基于CTW算法的建筑物形状相似性度量方法，该方法可以对齐具有不同节点个数的建筑物坐标序列，能直接根据矢量坐标序列进行形状相似性的计算，降低了算法的复杂性，提升了方法的可操作性。其中建筑物多边形的矢量坐标序列可以顾及图形原始的轮廓特征。试验证明CTW算法具有平移、旋转、缩放和镜像不变性，该算法计算得到的建筑物形状相似性符合人的视觉空间认知。值得注意的是，本文算法并不能处理结构复杂和噪声节点过多的建筑物形状，需要结合地图综合算法进行处理。未来将进一步探究复杂多边形的形状相似性度量，提升CTW算法的适用性。