CN110334861B - 一种基于轨迹数据的城市区域划分方法 - Google Patents

Abstract

本发明公开了一种基于轨迹数据的城市区域划分方法，包括数据预处理、轨迹数据网格划分和城市区域划分；数据预处理主要包括轨迹提取、数据清洗、数据过滤与插值和坐标转换四部分；经过数据预处理后，得到经过平面坐标转换后的轨迹数据集，利用轨迹数据集进行网格划分；在网格划分之后，利用形态学方法对网格划分得到的二值矩阵进行城市区域划分。本发明方法划分得到的城市子区域边界语义含义更加丰富，且可根据需求选择合适的参数进行区域划分，尺度变化多样，适用性更广。

Description

一种基于轨迹数据的城市区域划分方法

技术领域

本发明涉及城市区域划分领域，具体是一种基于轨迹数据的城市区域划分方法。

背景技术

城市为人们提供了生活和工作的场所，城市区域划分是许多与地理空间相关研究的基础和关键，如人类迁移、城市规划、路径规划等。随着城市化进程的加快，解决城市规划和道路规划不合理的问题以为人们提供更好的服务迫在眉睫。人类活动的轨迹数据蕴含了大量的时空信息，能有效地反映人类活动的规律，同时，这些数据为进行交通拥堵分析、交通模式推断、出行热点挖掘、移动模式发现等研究提供了良好的帮助。因此，提取轨迹数据中有价值的信息并应用到城市区域划分中是进行城市区域划分一种有效的解决方案。

现有的城市区域划分方法主要可以分为三类：基于几何图形的区域划分、基于路网的区域划分和基于数据挖掘的区域划分。

基于几何图形的区域划分最简单的即是将数字地图分割成大小相同的网格，其中每个网格都是矩形。这种方法虽然在实现上比较简单，但是它也有不足之处：细分的网格没有其他含义，因此基于这些网格的分析不能为我们提供对道路网络的自然理解。对比基于几何图形划分区域的方法，道路网络的形状和层级为区域划分提供了很好的基础，基于路网的划分方法的城市分区结果比较直观，具有一定的空间语义，但是没有结合实际的出行数据，区域不包含交通信息，且道路网络层级数量少划分的城市区域尺度单一。近几年的基于数据挖掘的区域划分方法，从真实数据入手，对区域进行划分，具有直接的说服力，体现了很强的空间语义。但是，现有的研究主要是基于数据的分析甄别交通热区或者判断土地功能，没有考虑区域边界的交通情况、区域内的联通性以及区域间的可达性等因素。

发明内容

本发明所要解决的技术问题是提供一种基于轨迹数据的城市区域划分方法，利用交通轨迹数据，通过一定的处理手段，使得划分的城市子区域交通便利且具有良好的连通性。

为解决上述技术问题，本发明采用的技术方案是：

一种基于轨迹数据的城市区域划分方法，包括以下步骤：

步骤1：对轨迹数据进行预处理，包括轨迹提取、数据过滤与差值、坐标变换；

轨迹提取：提取能够反映城市交通状况的轨迹数据序列，将这些轨迹数据序列作为城市区域划分的数据源；

数据过滤与插值：数据过滤是将相邻两个轨迹点距离小于一定阈值的轨迹点过滤掉；数据插值是将相邻两个轨迹点距离大于一定阈值的轨迹点之间插入一定数量的伪点；

坐标转换：对于上述数据，已知轨迹的经、纬度坐标分别为lon、lat，则其对应的平面横、纵坐标x、y为：

x＝(lon×20037508.34)/180

y＝log(tan((90+lat)×π/360)/(π/180)×20037508.34/180

步骤2：利用经过平面坐标转换后的轨迹数据集进行网格划分；

1)遍历城市区域边界数据的集合，确定研究区域的范围；

对于区域S，假设其平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，将区域S划分成若干个大小为k×k的单元网格，对每个网格单元内的轨迹点进行计数统计，并存储到二维矩阵中，即得到G_mn；假设水平方向被划分成了m个区间，竖直方向被划分成了n个区间，

2)根据网格单元大小，将轨迹数据映射到相应的网格单元内，得到网格化的数字矩阵；

遍历所有的轨迹点数据，根据映射关系将轨迹点数据映射到所属的网格单元；对每个网格单元的轨迹点数量进行计数，即得到网格化的数字矩阵；

3)矩阵编码，设定网格内轨迹点数目阈值θ，对上述数字矩阵进行重编码；

利用公式

将二维数字矩阵G_mn编码为0-1二值矩阵，其中，i、j表示网格单元的索引，i∈(0,m-1),j∈(0,n-1)；

4)矩阵变换，变换网格化数字矩阵，使得数字矩阵分布与北半球经纬度变化的方向相同；

利用公式G_ij′＝G_jn-1-i将二维矩阵G_mn变换成G_mn′，使得其输出的图像与轨迹数据在坐标系的分布相同；

步骤3：利用形态学方法对网格划分得到的二值矩阵进行城市区域划分；

膨胀操作：利用3×3十字形的卷积核对步骤2中得到的二值矩阵进行膨胀操作，即将目标像素向其4-邻接区域扩张，迭代两次；

细化操作：对膨胀操作后的二值矩阵进行细化操作，抽取二值图像的骨架；

连通区域标记：对处理后的二值矩阵采用3×3十字形的卷积核进行连通区域标记，即判断目标像素的4-邻接区域是否连通，将二值图像中各个连通区域标记出来，得到不同的城市子区域。

进一步的，所述步骤1中还包括对提取的数据进行数据清洗：即对数据轨迹序列，将包含行车速度或相近轨迹点距离超出正常阈值的GPS数据的轨迹序列删除。

进一步的，在步骤1中，设dis为相邻两点的距离，两点间距离的最小阈值、最大阈值设置为min_dis、max_dis；在数据过滤中，当dis小于min_dis时，删除后一个点，并累计计算多个点之间的距离，直至相邻两点之间距离达到最小阈值；在数据插值中，当dis大于max_dis时，线性添加一系列的伪点，通过

计算插入伪点的数量，再依次计算插入轨迹点的经纬度；插入的第n个点的经纬度坐标(lon_n,lat_n)由下式计算得出：

进一步的，在步骤2中，所述映射关系具体为：

已知区域边界平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，网格单元尺度为k，划分的网格个数为m×n，网格划分即为将[x_min,x_max]投影到[0,m-1]、将 [y_min,y_max]投影到[0,n-1]；设区域中的一个轨迹点的坐标为(x,y)，则轨迹点与对应网格单元的所属关系为：

其中，i、j表示网格单元的索引，i∈(0,m-1),j∈(0,n-1)，根据轨迹点与网格单元的所属关系，轨迹点与网格单元索引的映射函数为：

相应的，一个网格单元对应一个二维空间位置，其映射关系为：

x＝(i×k)+x_min,y＝(j×k)+y_min。

进一步的，所述步骤3中还包括细节处理：即针对膨胀操作和细化操作后的二值矩阵，删除区域内部留有的与区域边界不相连的轨迹段以及很小的连通区域。

与现有技术相比，本发明的有益效果是：1)划分得到的城市子区域边界为路网所在的位置，不仅包含了路网信息，还包含了一定的交通信息，语义含义更加丰富。2)具有多尺度划分区域的特点，可以根据需求选择合适的参数进行区域划分，尺度变化多样，适用性更广。

附图说明

图1是本发明整体流程示意图；

图2是本发明轨迹数据网格划分示意图；

图3是本发明轨迹数据网格化数字矩阵的示意图；

图4是本发明中不同网格单元尺度k下城市区域划分的结果图；

图5是本发明中不同矩阵编码阈值θ下城市区域划分的结果图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的说明。本发明方法包括数据预处理、轨迹数据网格划分、城市区域划分。

出租车作为城市公共交通的有效补充，每天都运行在城市道路上，出租车司机掌握了清晰的路况信息，出租车的轨迹数据更能有效地从轨迹数量上反应城市的交通状况，因此，本发明以出租车轨迹数据为研究对象进行说明。

数据预处理主要包括轨迹提取、数据清洗、数据过滤与插值和坐标转换四部分。具体详述如下：

(1)轨迹提取：以出租车为例，出租车轨迹数据包含了载客和空载两种状态的出行轨迹。出租车空载时，需要游离寻客，轨迹的规律性不明显；而出租车在载人的状态下选择的行驶路径通常为最短或最快行驶路径，能够反映城市的交通状况。因此，提取载客状态下的出租车轨迹数据序列作为区域划分的数据源。

(2)数据清洗：对于上述轨迹序列，由于环境等因素的影响少量轨迹序列存在噪声，删除此数据对实验结果影响不大，本发明中将包含行车速度或相近轨迹点距离超出正常阈值的GPS数据的轨迹序列删除掉。即对存在行车速度大于33.33m/s(即120km/h)轨迹段、或者相邻轨迹点距离大于1.5km轨迹段的轨迹序列，执行删除操作。

(3)数据过滤与插值：数据过滤是将相邻两个轨迹点距离小于一定阈值的轨迹点过滤掉，数据插值是将相邻两个轨迹点距离大于一定阈值的轨迹点之间插入一定数量的伪点。当两点的距离dis小于最大阈值max_dis时，删除后一个点，并累计计算多个点之间的距离，直至相邻两点之间的距离达到最小阈值min_dis；当两点的距离dis大于最大阈值max_dis时，线性添加一系列的伪点，使得两点之间的距离小于最大阈值max_dis，通过

计算插入伪点的数量，然后依次计算插入轨迹点的经纬度，插入的第n个点的经纬度坐标(lon_n,lat_n)可以由下式计算得出：

(4)坐标转换：对于上述数据，已知轨迹的经、纬度坐标分别为lon、lat，则其对应的平面坐标x、y为：

x＝(lon×20037508.34)/180

y＝log(tan((90+lat)×π/360)/(π/180)×20037508.34/180

经过数据预处理后，得到经过平面坐标转换后的轨迹数据集，利用该数据集进行网格划分，网格划分详细过程如下：

(1)遍历城市区域边界数据的集合，确定研究区域的范围。

对于区域S，假设其平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，对区域S分别从水平方向和竖直方向两个维度进行划分，把每个维度的取值范围等分为若干个边长为k的小区间，即可将区域S划分成若干个大小为k×k的单元网格，对每个网格单元内的轨迹点进行计数统计，并存储到二维矩阵中，即得到G_mn。假设水平方向被划分成了m个区间，竖直方向被划分成了n个区间，m、n分别由如下公式计算得到，采用向上取整是为了防止区域边界点落入网格外。

(2)根据网格单元大小，将轨迹数据映射到相应的网格单元内，得到网格化的数字矩阵。

已知区域边界平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，网格单元尺度为k，划分的网格个数为m×n，网格划分即为将[x_min,x_max]投影到[0,m-1]、将 [y_min,y_max]投影到[0,n-1]。设区域中的一个轨迹点的坐标为(x,y)，则轨迹点与对应网格单元的所属关系为：

其中，i、j表示网格单元的索引，i∈(0,m-1),j∈(0,n-1)，根据轨迹点与网格单元的所属关系，轨迹点与网格单元索引的映射函数为：

相应的，一个网格单元对应一个二维空间位置，其映射关系为：

x＝(i×k)+x_min,y＝(j×k)+y_min

遍历所有的轨迹点数据，根据映射关系将轨迹点数据映射到所属的网格单元，如图2所示。对每个网格单元的轨迹点数量进行计数，即得到网格化的数字矩阵，如图3所示。

(3)矩阵编码，设定网格内轨迹点数目阈值θ，对上述矩阵进行重编码。

网格单元内轨迹点数目一定程度的反映了该单元网格所在路段车辆行驶的频次，根据需求设定阈值，根据此阈值将轨迹点数量小于给定阈值θ的矩阵元素编码为0，反之，将矩阵元素编码为1。设置阈值θ用来筛选行车数量大于一定数量的道路，利用下面的公式将二维数字矩阵G_mn编码为0-1二值矩阵。

二值化的二维数字矩阵可以看作一个二值图像(Binary Image)。其中，“1”表示包含轨迹点的网格单元，可看作图像前景；“0”表示空白的网格单元，可看作图像背景。因此，此矩阵即可看作为一个二值图像。

(4)矩阵变换，变换网格化数字矩阵，使得数字矩阵分布与北半球经纬度变化的方向相同。

网格化的数字矩阵存储的形式如图3所示，将其输出为图像，与轨迹数据在坐标系中的分布不同，因此，需要进一步进行矩阵变换，利用下面的公式将二维矩阵G_mn变换成G_mn′，使得其输出的图像与轨迹数据在坐标系的分布相同。

G_ij′＝G_jn-1-i

矩阵G_mn′可以看作是基于栅格数据模型表示的轨迹数据。

在网格划分之后，本发明利用形态学方法对网格划分得到的二值矩阵进行城市区域划分。详细过程如下：

(1)膨胀操作

利用3×3十字形的卷积核对上述二值矩阵进行膨胀操作，即将目标像素向其4-邻接区域扩张，迭代两次，消除轨迹数据二值图像中间的噪音点以及孔洞。

(2)细化操作

二值图像的膨胀操作加宽了轨迹点实际网格数，图像细化可以抽取二值图像的骨架，在保证区域间连通性的同时，恢复因膨胀而减小的背景区域大小。利用3×3正方形的卷积核，即判断目标像素的8-邻接区域是否符合一定的条件，对上步操作得到的二值矩阵进行细化操作，抽取二值图像的骨架。

(3)细节处理

经过膨胀操作和细化操作后的图形连通区域标记后存在一些小的孔洞以及游离的轨迹段，因此我们需要在每一步进行一定的细节处理，删除区域内部留有的一些与区域边界不相连的轨迹段以及一些很小的连通区域。在膨胀操作中进行局部填补孔洞操作，将小于给定阈值的连通孔洞的像素值修改为前景像素值。在细化操作中进行轨迹段删除操作，将小于给定阈值的连通轨迹段的像素值修改为背景像素值。

(4)连通区域标记

对上述操作的二值矩阵采用3×3十字形的卷积核进行连通区域标记，即判断目标像素的4-邻接区域是否连通，将二值图像中各个连通区域标记出来，得到不同的城市子区域。

设置不同网格单元尺度k和矩阵编码阈值θ划分的城市结果分别如图4、图 5所示。

除出租车轨迹数据，公交车、私家车、共享单车等由出行而产生的轨迹数据都具备很强的研究价值，也都可以应用于本发明中。根据数据源的不同，划分的城市子区域也可以具有不同语义含义。比如，利用公交车的轨迹数据进行区域划分得到的子区域即表示公共交通便利的城市子区域；利用共享单车轨迹数据划分得到的子区域即表示共享单车出行便利的城市子区域。

Claims (5)

1.一种基于轨迹数据的城市区域划分方法，其特征在于，包括以下步骤：

步骤1：对轨迹数据进行预处理，包括轨迹提取、数据过滤与差值、坐标变换；

轨迹提取：提取能够反映城市交通状况的轨迹数据序列，将这些轨迹数据序列作为城市区域划分的数据源；

数据过滤与插值：数据过滤是将相邻两个轨迹点距离小于一定阈值的轨迹点过滤掉；数据插值是将相邻两个轨迹点距离大于一定阈值的轨迹点之间插入一定数量的伪点；

坐标转换：对于上述数据，已知轨迹的经、纬度坐标分别为lon、lat，则其对应的平面坐标x、y为：

x＝(lon×20037508.34)/180

y＝log(tan((90+lat)×π/360)/(π/180)×20037508.34/180

步骤2：利用经过平面坐标转换后的轨迹数据集进行网格划分；

1)遍历城市区域边界数据的集合，确定研究区域的范围；

对于区域S，假设其平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，将区域S划分成若干个大小为k×k的单元网格，对每个网格单元内的轨迹点进行计数统计，并存储到二维矩阵中，即得到G_mn；假设水平方向被划分成了m个区间，竖直方向被划分成了n个区间，

2)根据网格单元大小，将轨迹数据映射到相应的网格单元内，得到网格化的数字矩阵；

遍历所有的轨迹点数据，根据映射关系将轨迹点数据映射到所属的网格单元；对每个网格单元的轨迹点数量进行计数，即得到网格化的数字矩阵；

3)矩阵编码，设定网格内轨迹点数目阈值θ，对上述数字矩阵进行重编码；

利用公式

将二维数字矩阵G_mn编码为0-1二值矩阵，其中，i、j表示网格单元的索引，i∈(0,m-1),j∈(0,n-1)；

4)矩阵变换，变换网格化数字矩阵，使得数字矩阵分布与北半球经纬度变化的方向相同；

利用公式G_ij′＝G_jn-1-i将二维矩阵G_mn变换成G_mn′，使得其输出的图像与轨迹数据在坐标系的分布相同；

步骤3：利用形态学方法对网格划分得到的二值矩阵进行城市区域划分；

膨胀操作：利用3×3十字形的卷积核对步骤2中得到的二值矩阵进行膨胀操作，即将目标像素向其4-邻接区域扩张，迭代两次；

细化操作：对膨胀操作后的二值矩阵进行细化操作，抽取二值图像的骨架；

连通区域标记：对处理后的二值矩阵采用3×3十字形的卷积核进行连通区域标记，即判断目标像素的4-邻接区域是否连通，将二值图像中各个连通区域标记出来，得到不同的城市子区域。

2.如权利要求1所述的一种基于轨迹数据的城市区域划分方法，其特征在于，所述步骤1中还包括对提取的数据进行数据清洗：即对数据轨迹序列，将包含行车速度或相近轨迹点距离超出正常阈值的GPS数据的轨迹序列删除。

3.如权利要求1所述的一种基于轨迹数据的城市区域划分方法，其特征在于，在步骤1中，设dis为相邻两点的距离，两点间距离的最小阈值、最大阈值设置为min_dis、max_dis；在数据过滤中，当dis小于min_dis时，删除后一个点，并累计计算多个点之间的距离，直至相邻两点之间距离达到最小阈值；在数据插值中，当dis大于max_dis时，线性添加一系列的伪点，通过

计算插入伪点的数量，再依次计算插入轨迹点的经纬度；插入的第n个点的经纬度坐标(lon_n,lat_n)由下式计算得出：

4.如权利要求1所述的一种基于轨迹数据的城市区域划分方法，其特征在于，在步骤2中，所述映射关系具体为：

已知区域边界平面坐标的取值范围为x∈[x_min,x_max]、y∈[y_min,y_max]，网格单元尺度为k，划分的网格个数为m×n，网格划分即为将[x_min,x_max]投影到[0,m-1]、将[y_min,y_max]投影到[0,n-1]；设区域中的一个轨迹点的坐标为(x,y)，则轨迹点与对应网格单元的所属关系为：

其中，i、j表示网格单元的索引，i∈(0,m-1),j∈(0,n-1)，根据轨迹点与网格单元的所属关系，轨迹点与网格单元索引的映射函数为：

相应的，一个网格单元对应一个二维空间位置，其映射关系为：

x＝(i×k)+x_min,y＝(j×k)+y_min。

5.如权利要求1所述的一种基于轨迹数据的城市区域划分方法，其特征在于，所述步骤3中还包括细节处理：即针对膨胀操作和细化操作后的二值矩阵，删除区域内部留有的与区域边界不相连的轨迹段以及很小的连通区域。

Images