随着遥感和地理信息技术的发展，基于遥感影像数据的场景分类和基于地理时空大数据的社会感知被广泛应用于城市功能识别和居民行为模式分析中^[1-2]。其中，场景分类方法是一种面向对象的地物识别方法，基于遥感图像的光谱、形状和纹理等特征，结合概率语义模型，能够有效识别城市地物^[1]。然而，这种方法对于遥感影像的自然属性依赖性强，而城市功能不仅和地物的自然物理属性相关，在很大程度上，更是受到人群活动与社会经济因素的影响^[3-4]。地理空间大数据的兴起，使得大量能体现人群活动特征的数据被应用于城市功能结构的识别中^[5-8]，如百度感兴趣点(point of interests, POIs)数据、地铁公交刷卡数据和出租车轨迹数据等^{[6, 8-9]}。其中，出租车是城市生活的重要传感器，为城市居民行为模式研究和土地利用分类研究提供了强有力的工具^[10-13]；另外，也有学者基于POIs的频率分布对城市功能区进行识别^[10-14]。

出租车轨迹数据对城市功能区识别和土地利用分类作用显著，但之前的研究对于其时序特征的挖掘还不够充分，且缺乏对城市功能层次性的研究。如Pan等利用出租车轨迹数据提取的时间序列对城市功能区进行分类^[15]，该模型虽然效果显著，但分类算法在度量特征之间的相似性时未考虑时间序列的特殊性。此外，城市功能具有明显的层次性，不同层次的城市功能具有较为明显的层级关系^[16]，从系统功能上探讨这种层次性有助于从整体到细部逐步把握城市空间结构的功能变化和分布，但目前针对城市功能层次性的研究还较少。因此，本文引入动态时间规整描述出租车轨迹数据之间的相似性，以聚类算法挖掘其时序特征，并在广州市中心城区开展不同层次下的功能结构研究。本文模型首先对时序出租车出行数据构建出发-到达位置(orientation-destination, OD)的时序密度曲线；然后，基于动态时间规整对研究区进行聚类分析，并采用Silhouette值评估多层次对象；最后，结合OD时序曲线和POIs信息，对不同类别的功能区进行识别分析，在多层次的视角下讨论城市功能的功能属性和空间结构特征。

1 研究区和数据 1.1 研究区

研究区位于广州市中心城区，包括荔湾、越秀、天河和海珠4个区(见图 1)，是珠三角经济区的经济、政治、文化中心。据广州市2015年统计年鉴，研究区总面积为279.63 km²，约占广州市区面积的7.28%。2014年研究区内常住人口为514.38万人，约占广州市区人数的45.16%。

广州市的城市空间结构十分复杂^[17]，中心城区的城市功能分布彼此交错，比如，居民区、商业区等功能区的分布虽然在总体上呈现出一定的规律性，但很多区域的功能并不单一，而是由多种功能相互作用而形成的混合功能区域。这种混合功能区域的出现吻合《马丘比丘宪章》对城市分区综合性和多功能的追求^[18]，本研究也将从多层次的角度探究广州市城市土地混合利用的现状。

1.2 时序出租车轨迹数据

出租车轨迹数据源自广州市交通委员会，为了使研究结果能够反映工作日和休息日的居民出行情况，本文选取2009-05-12(星期二)和2009-05-16(星期六)的轨迹数据进行分析，共包括由14 379辆出租车产生的48 957 065条记录。其中，每条记录包含车辆的唯一ID号、入库时间、GPS信号时间、经纬度、速度、有效标示以及车辆状态，其中车辆状态包括空车、重车、点火、熄火等，GPS发送时间的间隔与车辆的状态有关，间隔4~7 s不等。相关研究表明^[10]，交通分析小区(traffic analysis zone, TAZ)中具有相同或相似的功能及特征，在城市功能结构的研究中以其作为分析单元是合理的。

2 多层次城市功能结构识别

首先从原始的出租车轨迹数据中提取出每一条轨迹的O点和D点，并利用点密度分析以小时为间隔生成工作日和休息日共48 h的O点密度分布和D点密度分布图，之后使用基于动态时间规整的K-MEDOIDS算法对OD点时序密度向量进行聚类分析，在Silhouette值的评估下获取研究区域的多层次聚类对象，同时引入POIs的统计结果进行辅助分析，具体流程如图 2所示。

2.1 特征数据集处理

在原始的出租车轨迹数据中，根据出租车的唯一标识符和轨迹时间顺序提取出每一条轨迹的OD点。对于任意的轨迹数据Tr，按轨迹点生成的先后顺序有轨迹向量T = [t₁…t_i…t_j…t_n]，提取规则为：若Tr中t_i时刻为空车，而t_i+1状态为重车，则t_i+1时刻对应的地理位置为起点；时间后推，遇到的第一个状态为空车的时刻t_j，则t_j时刻对应的地理位置为轨迹终点。由此，得到每一条轨迹的O点和D点，经处理后，有效数据记录为6 687 394条。另外，由于出租车O点和D点主要分布在道路上，而居民在出行时其实际出发点和目的点多分布于道路两侧的居民区、商业区等区域，估计该距离小于1 km，可通过搜索半径为1 km的点密度分析获取TAZ小区内的O、D点密度值。以1 h作为时间间隔，最终生成每个TAZ工作日和休息日共48 h的O点密度变化曲线和D点密度变化曲线。Pan等在利用出租车轨迹上下车点进行土地利用分类时，构造了6种OD点的时序组成方式^[15]：

1) $\left[ {o_w^1 \cdots o_w^i \cdots o_w^{24}, o_r^1 \cdots o_r^i \cdots o_r^{24}} \right] $为48维O点向量，$o_w^i、o_r^i $分别是工作日、休息日第i小时的O点平均数量。

2) $ \left[ {d_w^1 \cdots d_w^i \cdots d_w^{24}, d_r^1 \cdots d_r^i \cdots d_r^{24}} \right]$为48维D点向量，$ d_w^i、d_r^i$分别是工作日、休息日第i小时的D点平均数量。

3) $ \left[ {m_w^1 \cdots m_w^i \cdots m_w^{24}, m_r^1 \cdots m_r^i \cdots m_r^{24}} \right]$为48维O点和D点的差值向量，其中$m_w^i = o_w^i - d_w^i, m_r^i = o_r^i - d_r^i $。

4) $\left[ {r_w^1 \cdots r_w^i \cdots r_w^{24}, r_r^1 \cdots r_r^i \cdots r_r^{24}} \right] $为48维O点和D点的比值向量，$r_{w}^{i}=o_{w}^{i} / d_{w}^{i}, r_{r}^{i}=o_{r}^{i} / d_{r}^{i} $。

5) $ \left[ {{o_w}^1 \cdots {o_w}^i \cdots {o_w}^7} \right]$为7维向量，${{o_w}^i} $是一周第i天O点数量变化。

6) $\left[ {{d_w}^1 \cdots {d_w}^i \cdots {d_w}^7} \right] $为7维向量，${{d_w}^i} $是一周第i天D点数量变化。

在对上述6种特征进行评估时，48维O点向量的判别能力最高；组合多种特征评估时，96维OD点组合向量判别能力最优异，即$\left[ {o_w^1 \cdots o_w^{24}} \right., o_r^1 \cdots o_r^{24}, d_w^1 \cdots d_w^{24}, d_r^1 \cdots d_r^{24}] $。因此，本研究以96维O、D点时序密度向量作为城市功能区识别的输入向量进行聚类分析。

2.2 动态时间规整的K-MEDOIDS聚类

在进行时序数据分析时，需要度量时间序列的相似性，也即判断不同的时间序列是否表现为相似的行为^[19]。目前，时间序列聚类模型在测度相似性时根据目的不同可分为3类，分别依据序列在时间、形状和变化上的相似性进行聚类分析^[20]。出租车轨迹在分析城市居民的出行模式时，着重点在于识别曲线在形状上的变化模式，如以居住功能为主导的区域，早高峰阶段的O点处于高峰、D点处于低谷，而以工作功能为主导的区域，早高峰阶段的O点处于低谷、D点处于高峰。动态时间规整(dynamic time warping, DTW)是基于形状相似性的时间序列度量算法，可以实现时间序列的最佳匹配对齐，挖掘其中灵活的相似性和差异性，是度量时间序列相似性的最普遍使用的定量方法^[19-21]。

给定时间序列P =[p₁…p_i…p_n]和Q =[q₁…q_j…q_m]，构建n×m的矩阵M，其元素(i_th, j_th)是点p_i和点q_j之间的距离。DTW需保证对齐P和Q时差异最小，构建路径W= [w₁…w_k…w_K]，其中max(m, n) ≤K≤m+n－1，且需要满足3个条件：边界条件、连续性和单调性。边界条件需保证起点和终点分别是矩阵中左下角的点(1, 1)和右上角的点(n, m)。连续性和单调性限制w_k的下一个点w_k+1只能在w_k的上方或者右方，不能跨越两个格网，也不能向下方或者左方移动。在所有满足以上3个条件的路径中，选择d_DTW最小的一条作为输出结果，也即度量P和Q中距离最小的路径，其中d_DTW为：

$ d_{\mathrm{DTW}}=\min \frac{\sum\limits_{k=1}^{K} w_{k}}{K} $

(1)

本研究中，DTW在衡量交通小区96维OD点时间序列的差异性时，可消除地块时间变化趋势中的移位等现象，但这种移位不是毫无限制的，故实际度量两地块的相似度时需加以局部限制。针对96维时间序列$ \left[ {o_w^1 \cdots o_w^{24}} \right., o_r^1 \cdots o_r^{24}, d_w^1 \cdots d_w^{24}, d_r^1 \cdots d_r^{24}]$，首先施加第一重限制，在衡量任意两个时间序列时，$ \left( {o_w^1, o_w^1} \right)、\left( {o_w^{24}, o_w^{24}} \right)、\left( {o_r^1, o_r^1} \right)、\left( {o_r^{24}, o_r^{24}} \right)、\left( {d_w^1, d_w^1} \right)、\left( {d_w^{24}, d_w^{24}} \right)、\left( {d_r^1, d_r^1} \right)、\left( {d_r^{{2_1}}, d_r^{24}} \right)$这8个节点在DTW路径中必须匹配，保证工作日和休息日的O、D点匹配至对应的时间范围；其次施加第二重限制，在上述8个节点中任意相邻的两个节点，如$\left( {o_w^1, o_w^1} \right) $和$\left( {o_w^{24}, o_w^{24}} \right) $，需保证每隔6 h的时间点应得以匹配，也即在一天的时间序列中，06:00、12:00、18:00、24:00必须与对应的时间点匹配，以此限制移位的范围。

以DTW衡量任意两对时序曲线的差异后，本文使用K-MEDOIDS算法对TAZ进行聚类分析，同时为了提升算法效率，使用AP (anomalous pattern)方法帮助确定初始化全局聚类中心^[22]。另外，为了评估不同聚类数目结果的可靠性，本文引入Silhouette评估各簇的聚类质量^[23]，在K-MEDOIDS算法的背景下，假定在一个已有的聚类结果中，a(i)指代样本点i与同类中其他所有样本点的DTW距离的中位值，b(i)指代其他类与样本点i之间的最小DTW距离的中位值，则有：

$ s\left( i \right){\rm{ }} = {\rm{ }}\frac{{b\left( i \right) - a\left( i \right)}}{{{\rm{max}}\left\{ {{\rm{}}a\left( i \right), b\left( i \right)} \right\}}} $

(2)

如果s(i)接近于1，意味着样本点i和现有的聚类结果十分匹配；如果s(i)接近于-1，则意味着样本点i应该隶属于其邻类中。所有点的s(i)的均值越高，则代表聚类结果越优。

2.3 POIs数据辅助分析

POIs点包含区域功能的语义信息，在城市功能识别中发挥着重要的作用^[10]，本研究使用POIs作为辅助信息帮助识别广州市中心城区的TAZ功能属性。POIs点的密度反映了该功能区不同种类感兴趣点的绝对量值差异，可用于辅助判断功能区的实际功能属性；但POIs的绝对量值也可能掩盖区域内具备实际主导性的属性信息，为避免TAZ中不同类型POIs的数量差异对识别结果的影响，同时引入POIs富集指数作为辅助判别信息^[24]。

POIs点的密度表示为：

$ {\rm{Densit}}{{\rm{y}}_{{\rm{POI(}}i, j{\rm{)}}}} = \frac{{{\rm{Nu}}{{\rm{m}}_{{\rm{POI(}}i, j{\rm{)}}}}}}{{{\rm{Are}}{{\rm{a}}_j}}} $

(3)

式中，Density_{POI(i, j)}是第j类功能区中第i类POI的密度；Num_{POI(i, j)}是第j类功能区中第i类POI的数量；Area_j是第j类功能区的总面积。

POIs富集指数表示为：

$ {F_{i, j}} = \frac{{{\rm{ }}{n_{i, j}}/{n_j}}}{{{N_i}/N}} $

(4)

式中，F_{i, j}是第j类功能区中第i类POI的富集指数；n_{i, j}是第j类功能区中第i类POI的数量；n_j是第j类功能区中POIs的总数量；N_i是第i类POI的总数量；N是整个研究区中所有POIs的总数量。F_{i, j}越高，表明第j类功能区中第i类POI的富集程度越高。

3 结果和讨论

本研究数据处理和算法实现均采用C++编程实现，使用GDAL、Shark等开源库，在程序执行过程中通过启用OpenMP实现高速并行运算。

为了研究广州市城市功能结构的层次性，本研究针对96维OD点向量进行多次聚类，并从中选取Silhouette值较高的聚类结果进行分析。如图 3所示，聚类数目分别设为k=2~10，Silhouette值都大于0.5，具备较好的解释性，从中选择较具代表性的3个层次进行详细分析。Silhouette值可划分为3段：k=2时，聚类结果最为理想，Silhouette值高达0.87；k=3~4时，聚类结果次之，Silhouette值约为0.76；k=5~10时，聚类结果再次之，Silhouette值分布在0.51~0.62之间。为了使分析结果更具代表性，本研究从3段中各选择Silhouette值最高的聚类结果进行分析，即k=2、3、7。

当聚类数目k=2、3时，分别如图 4、图 5所示。图 4中A1位于研究区域的中心地带，A0呈环状围绕在A1四周，构成简单的环状圈层结构。

图 5中其聚类结果继承了k=2时的环状圈层结构，B0区域几乎包含于A0的范围内，B2区域几乎包含于A1的范围内，B1区域由A0、A1相接部分构成。

图 6为A0和A1的OD点时序变化曲线，在工作日：A0区域在上午时段的O点密度高于同期D点密度，晚上时段的D点密度明显高于同期O点密度，属于典型的居住功能区；A1区域在10:00、15:00、22:00存在明显的峰值，居民的流通高峰较为集中，并且17:00-24:00时，O点密度略高于同期的D点密度，符合工作功能区的居民行为模式。在休息日：A0区域08:00的O、D点密度是较为明显的峰值，符合居民休息日外出游玩的特征，且上午时段O点密度持续高于D点密度；A1区域的O、D点密度曲线不存在明显的峰值，变化幅度较工作日显著减弱，其工作属性特征降低，进一步验证了A0、A1区域的主导功能。此时，研究区的功能体系呈现为以居民区环绕工作区的“职-住”二元环状圈层结构。

图 7为B0、B1、B2 3类功能区的O、D点时序密度曲线，B2区域的O、D点密度曲线变化模式和A1区域几乎一致，但在量值上高于后者，故B2是工作功能区；B0、B1区域的OD点密度曲线分别在14:00—24:00、07:00—14:00, 与A0区域保持一致，推定B0、B1都具备居住区的功能特点。针对B0区域，与A0相比，O、D点密度在上午时段的变化幅度更小，而在17:00—24:00时, O、D点密度差的差异变大，符合居住功能为主导的居民出行特点，故B0是居住功能区；针对B1区域，与A0区域相比，不论工作日和休息日，O、D点密度整体较大，且上午时段的变化幅度增大，而下午时段的O、D点密度差异变小，该曲线模式符合居民的游憩行为，推定B1区域具备游憩功能，同时B1区域的O、D点在工作日10:00存在明显的峰值，符合居民的工作行为，故最终推定B1区域是居住游憩工作混合功能区，此时，混合功能区出现，城市功能的综合性和多元性已逐渐显现。

当聚类数目k=7时，如图 8所示，由C0~C6组成了由外向内的环状圈层结构。从空间继承性上分析，C0区中96.88%的TAZ属于B0，C6区中98.97%的TAZ属于B2，C1由B0、B1相接部分组成，C5由B2、B1相接部分组成，C2、C3、C4基本包含在B1内。如图 9、图 10和表 1所示，考察OD点时序变化曲线、POIs密度及富集指数表，C6的曲线模式与B2一致，且各类型的POIs密度都居于前列，其中公司企业、金融服务及休闲娱乐的POIs具有高富集性，商业配套设施十分完善，故C6是中心商业区；C5的O、D点曲线模式和B1、B2曲线的组合模式相似，如工作日上午时段与B1较为一致，晚上时段与B2较为一致，同时商业性、休闲性的POIs数量丰富，各类POIs的富集指数较为均匀，故C5是商业休闲区。C0的O、D曲线模式与B0一致，涵盖居住区的功能属性，同时公司企业、休闲娱乐等POIs类型具有很高的富集性，故C0是居住游憩混合功能区；C1区域的O、D曲线模式和B0、B1都存在差异，而与A0区域的O、D曲线十分相似，且该区域内地产小区的富集指数极高，故C1是居住功能区。C2、C3、C4区域主要位于B1内，曲线模式和B1都存在相似性，属于居住游憩工作混合区域，在将其进一步细化为3种不同类型的功能区时，三者的城市功能表现为更为复杂的特征，每一种功能区都呈现为“职-住-休”的三元混合状态，但其功能属性各有倾向，如C2中，生活服务的POIs富集性较高，C3中，金融服务和休闲娱乐的POIs富集性较高，C4中，医疗教育的POIs富集性较高。此时，城市功能的综合性和多元性已充分表现出来，研究区的功能体系更为复杂，但依然保留了k=3时以“职-住-休”为基本功能框架的特点，同时继承了k=2、3时的环状圈层结构，功能结构的层次性已十分明显。

总体来说，从k=2到k=3，再到k=7，广州市中心城区的功能结构呈现出明显的层次性：当k=2时，研究区域的功能表现是“职-住”二元结构，空间结构以环状圈层实现了中心城区的二元划分，可为城市居民通勤提供参考；当k=3时，城市功能细化，由“职-住”二元结构向“职-住-休”三元结构发展，功能结构更为完善，环状圈层结构的层次性提升，由外围的居住区逐渐向中心城区的工作区过渡；当k=7时，城市功能在细化的同时逐渐向多元性和混合性发展，虽然依然符合“职-住-休”的三元框架，但土地混合利用成为城市功能的总体特征，环状圈层结构进一步细化，由外围的居住游憩向中心的商业休闲过渡，混合功能区根据POIs的富集和主体特征表现为不同的功能属性倾向。总体来说，本研究将城市功能和空间结构结合起来，有效地阐述了城市功能及空间结构的多层次性特征，这对于从系统性和层次性的特征把握城市土地利用和结构特性十分重要。

4 结语

城市化进程的不断加快，一方面造就了城市的快速扩张^[25]，另一方面导致城市内部土地利用的多样性^[26]，有效识别城市内部功能和空间结构对于城市管理和规划至关重要。本研究构建了基于出租车时序轨迹数据的多层次城市功能识别模型。首先，使用点密度分析构建TAZ的工作日和休息日出租车时序密度曲线；其次，引入动态时间规整和K-MEDOIDS聚类对广州市中心城区的功能属性进行聚类分析，使用Silhouette值评估聚类结果的质量；最后，辅以POIs语义信息，探究广州市中心城区的功能特征和空间结构。从层次性上分析，其功能属性逐渐由“职-住”二元结构向“职-住-休”三元结构转变；从空间结构上分析，功能属性具备明显的环状圈层特征，随着层次的细化，该结构愈加明显；从功能变化上看，中心城区的功能属性明显存在由环状外围向中心过渡的特征，当层次低时，由居住向工作过渡，当层次稍高时，过渡之中增添了游憩功能，当层次再次提高时，功能属性过渡逐渐复杂，并在不同的圈层上呈现出不同的属性倾向。

总体来说，本研究为城市管理和规划提供了多层次分析的视角，对于城市内部土地利用模拟具有重要的意义。此外，本研究尚有不足，如对于城市功能的层次性和环状圈层结构论证还不够充分，出租车作为城市交通流的重要构成对城市交通模式具有很强的代表性^[10]，但公交车、地铁、自驾车等也承载了大量的交通流量，对城市交通模式的挖掘存在影响；同时，时间尺度和空间尺度的选择对于聚类结果会产生一定的影响，在今后的研究中，需要多源融合不同交通流数据在更多的城市进行城市功能结构的论证工作，并尝试着从不同的时空尺度分析中心城区的功能结构，验证研究结果的普遍性。