一种基于GPS数据的工业园区物流关键路段识别方法

摘要

本发明提供一种基于GPS数据的工业园区物流关键路段识别方法，该方法先获取电子围栏内供应商及工业园区的物流信息、网络拓扑结构数据、北斗GPS数据，通过对上述数据进行计算以获取节点拓扑结构参数、节点物流吸引力、路段拓扑结构参数、路段动态平均速度、路段动态服务水平及路段动态交通量，之后再利用上述数据获得路段的物流吸引力和路段动态通行能力，进而获取路段动态关键度，然后根据关键度排序动态识别物流关键路段。该方法在识别时考虑城市建成环境因素、物流企业分布特征和物流车辆的时空特征，对物流关键路段的识别更具针对性和准确性，而且实时读取的物流车辆的GPS数据保持了较高的实时度。

说明书

技术领域

本发明属于城市路网规划领域，涉及关键路段识别，具体涉及一种基于GPS数据的工业园区物流关键路段识别方法。

背景技术

在城市路网规划领域，由于城市道路建成环境因素会影响生产物料同城配送；因此，十分需要针对工业园区及其吸引范围内的干线交通进行精细化组织，以保障生产物料的有效送达和精益化供应。混合交通流条件下，生产服务型物流车辆的同城配送与居民出行共同对城市交通产生影响，采集分析物流车辆的动态轨迹数据，挖掘物流车辆的出行时空分布特征，实时识别出物流关键路段，可为工业园区乃至城市交通精细化组织提供有力决策支持。

现有的物流关键路段识别方法主要有三种，第一种是采用出租车、公交车GPS数据挖掘方法识别交通异常[201510060839.9,201610483132.3]、调配车辆[201710131675.3]、判定公交运行状态[202010081991.6,201811407198.X，201610483132.3]，通过大量个体轨迹揭示群体的时空行为规律，但利用浮动车的大规模数据识别关键路段，不能体现出物流车辆在工业园区或城市配送的动态出行规律，无法体现物流关键路段的本质特征，不能直接用于识别工业园区物流关键路段。

第二种是直接采用北斗GPS车辆动态轨迹数据的数据挖掘方法，对物流车辆进行轨迹纠偏[201810675702.8，201711337827.1]或者重型货车交通流量预测[201911336098.7]，但上述两种方式只考虑重型货车出行偏好和时空特征，都未考虑城市建成环境因素(如道路路况和物流设施空间分布)对交通出行、道路活力以及对物流车辆运动轨迹的影响，识别参数单一，缺乏路段相关性数据的综合考量和判断，不适用于工业园区及其附近的物流关键路段识别。

第三种是采用实地勘探和人工测算方法，获得定义的复杂交通网络交通运行数据，基于复杂网络理论和方法，建立加权复杂交通网络仿真模型，通过去除和添加某路段，仿真评估网络性能变化来判定该路段是否为关键路段，例如一种基于复杂网络的城市关键道路的发掘方法[201910405988.2]，类似相关的方法，还有一种基于复杂网络相关性的关键节点识别方法[201911021967.7]，基于复杂网络的多层城市交通网构建及关键节点识别方法[201811232769.0]，基于复杂网络和手机信令数据的城市交通枢纽点评价方法[201710616366.5]等。但是这种方法总体上效率较低，所需时间较长，无法对关键路段进行实时识别，分析结果不能体现不同时序下关键路段的变化，时效性较差，导致分析结果与实际情况变化不能同步。同时在关键步骤上，该方法依赖特定数据采集器(线圈传感器、录像设备等)收集车辆在路段或节点的速度数据，且在构建复杂交通网络模型时只关注城市交通网络的拓扑属性和空间属性，未能结合城市建成环境因素，如道路路况和物流设施空间分布对物流关键路段的实质性影响，因而，所发掘出的关键路段不具有物流属性，亦缺乏针对性和有效性。

发明内容

鉴于上述技术问题，本发明的目的在于提供一种基于GPS数据的工业园区物流关键路段识别方法，该方法在识别时考虑城市建成环境因素、物流企业分布特征和物流车辆的时空特征，对物流关键路段的识别更具针对性和准确性，而且通过实时读取的物流车辆的GPS数据，便可自动识别物流关键路段，保持了较高的实时度。

为实现上述目的，本发明是采用如下技术方案实现的：

一种基于GPS数据的工业园区物流关键路段识别方法，该方法包括以下步骤：

步骤1、根据供应商分布和工业园区的交通影响范围，划定电子围栏；

步骤2、获取电子围栏内供应商及工业园区的物流信息，包括供应商的经纬度坐标、供给物流量、工业园区入口的经纬度坐标、需求物流量；

步骤3、预处理电子围栏范围内网络拓扑结构数据，使格式统一，所述网络拓扑结构数据包括围栏范围内的节点信息和路段信息；所述节点信息包括节点编号、节点经纬度坐标、节点各相位绿信比；所述路段信息包括路段编号、路段起止点、路段等级、路段长度、路段车道数、路段设计通行速度；

步骤4、获取北斗GPS数据，并进行预处理，数据包括物流车辆ID信息、物流车辆位置的经纬度信息与相应的时间戳、物流车辆与正北方向夹角、物流车辆所处位置海拔、物流车辆瞬时速度；

步骤5、根据步骤3的节点信息和路段信息计算节点拓扑结构参数，包括节点连接度和节点的通行能力；

节点i的连接度的计算公式如下：

其中，k

节点i的通行能力CAP

其中，CAP

步骤6、根据步骤2的供应商和工业园区的物流信息以及步骤5的节点通行能力计算节点的物流吸引力，

节点i的物流吸引力的计算公式如下：

其中，f

步骤7、根据步骤3的路段信息和步骤5的节点连接度计算路段拓扑结构参数，包括路段连接度、路段权值和路段拓扑值；

路段a的连接度K

K

其中，节点i和节点j分别为路段的起止点，k

路段a的拓扑值T

T

其中，L

步骤8、根据步骤4的GPS数据计算路段的动态平均速度，之后判断给定粒度的时间窗下是否有多辆车经过路段a，若是，则计算该时间窗下路段a的所有动态平均速度的平均值，记录为

路段a的动态平均速度V

其中，d

步骤9、根据步骤3的路段等级和步骤8路段的所有动态平均速度的平均值，确定路段a的动态服务水平S

步骤10、读取步骤4的GPS数据中的全部轨迹点，进行DBSCAN聚类分析,记录聚类分析所得簇数，根据簇中点的经纬度，确定每个簇所包含的路段，之后计算路段a的动态交通量Q

其中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σ

步骤11、根据步骤6节点的物流吸引力计算路段a的物流吸引力，计算公式如下：

F

其中，节点i和节点j分别为路段a的起止点，f

步骤12、根据步骤7的路段拓扑值、步骤9的路段动态服务水平以及步骤10的路段动态交通量，计算路段的动态通行能力，计算公式如下：

其中，E

步骤13、计算路段的动态关键度，并按降序排序，计算公式如下：

其中，R

步骤14、识别关键路段数量，根据路段动态关键度顺序，动态识别物流关键路段，关键路段数量计算公式如下：

L

其中，L

作为本发明的优选，所述路段等级与路段权值及路段动态服务水平关系如下表：

作为本发明的优选，节点i到供应商a的距离d

其中，节点i的经纬度坐标为(lon

作为本发明的优选，步骤10路段a的动态交通量Q

步骤10.1、读取步骤4的GPS数据中的全部轨迹点；

步骤10.2、计算样本内每个点与所有点之间的距离d

步骤10.3、生成距离分布矩阵D

步骤10.4、将距离分布矩阵D

步骤10.5、计算K-平均最近距离R参数列表，计算公式如下：

其中，R为邻域半径，D

步骤10.6、对于给定的R，依次计算邻域密度阈值ε，计算公式如下：

其中，P

步骤10.7、生成邻域密度阈值ε参数列表D

步骤10.8、给定初始K＝1；

步骤10.9、将邻域半径

步骤10.10、记录聚类分析所得簇数M

步骤10.11、判断簇数M

步骤10.12、令K＝K+1，并转步骤10.9；

步骤10.13、记录最优簇数M＝M

步骤10.14、将邻域半径

步骤10.15、判断聚类分析所得簇数M

步骤10.16、令K＝K+1，并步骤10.14；

步骤10.17、令K＝K-1，并将邻域半径

步骤10.18、根据簇中点的经纬度，确定每个簇所包含的路段；

步骤10.19、计算每个簇平均邻域密度参数列表，计算公式如下：

其中，

步骤10.20、计算路段动态交通量Q

其中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σ

本发明的优点和有益效果：

(1)本发明提供的关键路段识别方法只需获取物流车辆的GPS数据，即可自动识别物流关键路段，无需线圈传感器、录像设备等硬件设备，也不必进行传统录像、观测、记录，即可为交通精细化组织提供决策依据。

(2)本发明自动更新北斗GPS历史数据并实时读取，通过分析海量数据，实现对物流关键路段的实时识别，保持了较高的实时度和准确性。

(3)本发明通过动态数据特征提取和自适应调整的DBSCAN聚类算法，自动识别物流关键路段，防止人工调整参数和后期维护产生误差，更具准确性。

(4)本发明提供的关键路段识别方法考虑了城市建成环境因素、物流企业分布特征和物流车辆的时空特征，对物流关键路段的识别更具针对性和准确性。

附图说明

图1为本发明工业园区物流关键路段识别方法流程图。

具体实施方式

为使本发明要解决的技术问题、技术方案更加清楚，下面将结合附图及具体实施案例进行详细描述。应当理解，此处所描述的实施例仅用于说明和解释本发明，并不用于限定本发明。

参阅图1，本发明提供的一种基于GPS数据的工业园区物流关键路段识别方法，该方法包括以下步骤：

步骤1、根据供应商分布和工业园区的交通影响范围，划定电子围栏；

步骤2、获取电子围栏内供应商及工业园区的物流信息，包括供应商的经纬度坐标、供给物流量(通过供应商获得)、工业园区入口的经纬度坐标、需求物流量(通过供应商获得)；

步骤3、预处理电子围栏范围内网络拓扑结构数据，使格式统一(可以统一转换成84坐标系)，所述网络拓扑结构数据包括围栏范围内的节点信息和路段信息；所述节点信息包括节点编号(根据区域内有多少个十字口确定)、节点经纬度坐标、节点各相位绿信比等；所述路段信息包括路段编号、路段起止点、路段等级、路段长度、路段车道数、路段设计通行速度等；

步骤4、获取北斗GPS数据，并进行预处理，数据包括物流车辆ID信息，物流车辆位置的经纬度信息与相应的时间戳，物流车辆与正北方向夹角，物流车辆所处位置海拔，物流车辆瞬时速度；

其中，采集到的GPS数据往往存在缺失与误差，这时候需要进行数据筛选和数据补偿，同时还需要根据电子围栏范围，筛选GPS数据，剔除不在电子围栏范围内的GPS轨迹数据；设置时间窗粒度θ，单位为min，一般取30min，将步骤4的GPS数据按照时间窗粒度划分为不同的数据集合，之后根据分析需要，选择对应时间窗粒度的北斗GPS数据进行分析；

步骤5、根据步骤3的节点信息和路段信息计算节点拓扑结构参数，包括节点连接度和节点的通行能力，具体过程包括步骤5.1至步骤5.7：

步骤5.1、根据步骤3的路段起止点信息，确定节点集合的邻接矩阵；

步骤5.2、根据步骤5.1的邻接矩阵计算所有节点的入度、出度和连接度；

设有向路网中有N个节点，第i个节点的入度

步骤5.3、判断节点i与节点j是否可达，若可达，则转步骤5.4，若不可达，则转步骤5.5；

步骤5.4、根据步骤3的路段等级确定左转/直行/右转车道数，然后按照国标确定该路段各车道的规划饱和流量；

步骤5.5、读取下一个节点j，并转步骤5.3；

步骤5.6、判断计算次数是否等于节点i的连接度，若等于，则转步骤5.7，否则，转步骤5.5；

步骤5.7、根据步骤3的节点各相位绿信比，计算节点的通行能力，计算公式如下：

式中，CAP

步骤5.8、读取下一个节点i，并转步骤5.3；

步骤5.9、判断所有节点是否均参与计算，若是，则转步骤6，否则，转步骤5.8；

步骤6、根据步骤2的供应商和工业园区的物流信息以及步骤5的节点通行能力计算节点的物流吸引力，具体过程包括步骤6.1至步骤6.10：

步骤6.1、根据步骤2的供应商经纬度坐标，计算节点i到供应商A的距离d

设节点i的经纬度坐标为(lon

式中，R为地球半径，通常取6378.137(KM)，到工业园区的距离d

步骤6.2、对供应商A，计算节点i对供应商的吸引力f

步骤6.3、对工业园区Γ，计算节点i对工业厂区入口的吸引力f

步骤6.4、判断所有节点是否均参与计算，若是，则转步骤6.6，否则，转步骤6.5；

步骤6.5、读取下一个节点i，并转步骤6.1；

步骤6.6、判断所有供应商是否均参与计算，若是，则转步骤6.8，否则，转步骤6.7；

步骤6.7、读取下一个供应商A，并转步骤6.1；

步骤6.8、根据步骤2中的供应商供给物流量，计算节点i的物流吸引力f

式中，W

步骤6.9、判断所有节点是否均参与计算，若是，则转步骤7，否则，转步骤6.10；

步骤6.10、读取下一个节点i，并转步骤6.8；

步骤7、根据步骤3的路段信息和步骤5的节点连接度计算路段拓扑结构参数，包括路段连接度、路段权值和路段拓扑值，具体过程包括步骤7.1至步骤7.5：

步骤7.1、计算路段a的连接度K

K

式中，节点i和节点j分别为路段的起止点，k

步骤7.2、根据步骤3的路段等级，确定该路段的权值D

表1路段权值的确定方式

步骤7.3、计算路段a的拓扑值T

T

式中，L

步骤7.4、判断所有路段是否均参与计算，若是，则转步骤8，否则，转步骤7.5；

步骤7.5、读取下一个路段a，并转步骤7.1；

步骤8、计算路段的动态平均速度，具体过程包括步骤8.1至步骤8.15：

步骤8.1、根据步骤4的一组GPS数据，令x＝1，读取GPS数据的第x行信息；

步骤8.2、判断第x行的轨迹位置是否位于节点i处，若位于，则转步骤8.3，否则，转步骤8.14，判断方法为：

第x行点的经度与节点i的经度差小于等于阈值μ，且第x行点的维度与节点i的维度差小于等于阈值μ，μ一般取0.002；

步骤8.3、记录此时时间戳为t

步骤8.4、令x＝x+1，并读取GPS数据的第x行信息；

步骤8.5、判断第x行的轨迹位置是否位于节点i和节点j之间，若是，则转步骤8.6，否则，转步骤8.7，判断方法为：

第x行的轨迹位置到节点i和节点j的距离之和与两节点间距离之差是否小于等于0.1，单位为km；

步骤8.6、令x＝x+1，并读取GPS数据的第x行信息，并转步骤8.5；

步骤8.7、判断第x行的轨迹位置是否位于节点j处，若位于，则转步骤8.8，否则，转步骤8.12；

步骤8.8、记录此时时间戳为t

步骤8.9、计算路段a(节点i和节点j之间)的动态平均速度V

式中，d

步骤8.10、判断给定粒度的时间窗下，是否有多辆车经过路段a，若是，则转步骤8.11，若否，则转步骤8.12；

步骤8.11、计算该时间窗下，路段a的所有动态平均速度的平均值，记录为

步骤8.12、判断所有节点j是否均参与计算，若是，则转步骤8.14，否则，转步骤8.13；

步骤8.13、读取下一个节点j，并转步骤8.1；

步骤8.14、判断所有节点i是否均参与计算，若是，则转步骤9，否则，转步骤8.15；

步骤8.15、读取下一个节点i，并转步骤8.1；

步骤9、确定路段的动态服务水平，具体过程包括步骤9.1至步骤9.3：

步骤9.1、根据步骤3的路段等级和步骤8.11的路段的所有动态平均速度的平均值，确定路段a的动态服务水平S

表2路动态段服务水平确定方式

步骤9.2、判断所有路段是否均参与计算，若是，则转步骤10，否则，转步骤9.3；

步骤9.3、读取下一个路段a，并转步骤9.1；

步骤10、计算路段的动态交通量，具体过程包括步骤10.1至步骤10.20：

步骤10.1、读取步骤4的GPS数据中的全部轨迹点；

步骤10.2、计算样本内每个点与所有点之间的距离d

步骤10.3、生成距离分布矩阵D

步骤10.4、将距离分布矩阵D

步骤10.5、计算K-平均最近距离R参数列表，计算公式如下：

式中，R为邻域半径，D

步骤10.6、对于给定的R，依次计算邻域密度阈值ε，计算公式如下：

式中，P

步骤10.7、生成邻域密度阈值ε参数列表D

步骤10.8、给定初始K＝1；

步骤10.9、将邻域半径

步骤10.10、记录聚类分析所得簇数M

步骤10.11、判断簇数M

步骤10.12、令K＝K+1，并转步骤10.9；

步骤10.13、记录最优簇数M＝M

步骤10.14、将邻域半径

步骤10.15、判断聚类分析所得簇数M

步骤10.16、令K＝K+1，并转步骤10.14；

步骤10.17、令K＝K-1，并将邻域半径

步骤10.18、根据簇中点的经纬度，确定每个簇所包含的路段；

步骤10.19、计算每个簇平均邻域密度参数列表，计算公式如下：

式中，

步骤10.20、计算路段动态交通量Q

式中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σ

步骤11、计算路段的物流吸引力，具体过程包括步骤11.1至步骤11.3：

步骤11.1、计算路段a的物流吸引力F

F

式中，节点i和节点j分别为路段a的起止点，f

步骤11.2、判断所有路段是否均参与计算，若是，则转步骤12，否则，转步骤11.3；

步骤11.3、读取下一个路段a，并转步骤11.1；

步骤12、计算路段的动态通行能力，具体过程包括步骤12.1至步骤12.3：

步骤12.1、计算路段a的动态通行能力E

式中，E

步骤12.2、判断所有路段是否均参与计算，若是，则转步骤13，否则，转步骤12.3；

步骤12.3、读取下一个路段a，并转步骤12.1；

步骤13、计算路段的动态关键度，并按降序排序，计算公式如下：

式中，R

步骤14、识别关键路段数量，根据路段动态关键度顺序，动态识别物流关键路段，关键路段数量计算公式如下：

L

式中，L