本发明提供一种基于GPS数据的工业园区物流关键路段识别方法，该方法先获取电子围栏内供应商及工业园区的物流信息、网络拓扑结构数据、北斗GPS数据，通过对上述数据进行计算以获取节点拓扑结构参数、节点物流吸引力、路段拓扑结构参数、路段动态平均速度、路段动态服务水平及路段动态交通量，之后再利用上述数据获得路段的物流吸引力和路段动态通行能力，进而获取路段动态关键度，然后根据关键度排序动态识别物流关键路段。该方法在识别时考虑城市建成环境因素、物流企业分布特征和物流车辆的时空特征，对物流关键路段的识别更具针对性和准确性，而且实时读取的物流车辆的GPS数据保持了较高的实时度。
一种基于GPS数据的工业园区物流关键路段识别方法
技术领域
[0001] 本发明属于城市路网规划领域，涉及关键路段识别，具体涉及一种基于GPS数据的工业园区物流关键路段识别方法。
背景技术
[0002] 在城市路网规划领域，由于城市道路建成环境因素会影响生产物料同城配送；因此，十分需要针对工业园区及其吸引范围内的干线交通进行精细化组织，以保障生产物料
的有效送达和精益化供应。混合交通流条件下，生产服务型物流车辆的同城配送与居民出
行共同对城市交通产生影响，采集分析物流车辆的动态轨迹数据，挖掘物流车辆的出行时
空分布特征，实时识别出物流关键路段，可为工业园区乃至城市交通精细化组织提供有力
决策支持。
[0003] 现有的物流关键路段识别方法主要有三种，第一种是采用出租车、公交车GPS数据挖掘方法识别交通异常[201510060839.9,201610483132.3]、调配车辆[201710131675.3]、
判定公交运行状态[202010081991.6,201811407198.X，201610483132.3]，通过大量个体轨
迹揭示群体的时空行为规律，但利用浮动车的大规模数据识别关键路段，不能体现出物流
车辆在工业园区或城市配送的动态出行规律，无法体现物流关键路段的本质特征，不能直
接用于识别工业园区物流关键路段。
[0004] 第二种是直接采用北斗GPS车辆动态轨迹数据的数据挖掘方法，对物流车辆进行轨迹纠偏[201810675702.8，201711337827.1]或者重型货车交通流量预测
[201911336098.7]，但上述两种方式只考虑重型货车出行偏好和时空特征，都未考虑城市
建成环境因素(如道路路况和物流设施空间分布)对交通出行、道路活力以及对物流车辆运
动轨迹的影响，识别参数单一，缺乏路段相关性数据的综合考量和判断，不适用于工业园区
及其附近的物流关键路段识别。
[0005] 第三种是采用实地勘探和人工测算方法，获得定义的复杂交通网络交通运行数据，基于复杂网络理论和方法，建立加权复杂交通网络仿真模型，通过去除和添加某路段，
仿真评估网络性能变化来判定该路段是否为关键路段，例如一种基于复杂网络的城市关键
道路的发掘方法[201910405988.2]，类似相关的方法，还有一种基于复杂网络相关性的关
键节点识别方法[201911021967.7]，基于复杂网络的多层城市交通网构建及关键节点识别
方法[201811232769.0]，基于复杂网络和手机信令数据的城市交通枢纽点评价方法
[201710616366.5]等。但是这种方法总体上效率较低，所需时间较长，无法对关键路段进行
实时识别，分析结果不能体现不同时序下关键路段的变化，时效性较差，导致分析结果与实
际情况变化不能同步。同时在关键步骤上，该方法依赖特定数据采集器(线圈传感器、录像
设备等)收集车辆在路段或节点的速度数据，且在构建复杂交通网络模型时只关注城市交
通网络的拓扑属性和空间属性，未能结合城市建成环境因素，如道路路况和物流设施空间
分布对物流关键路段的实质性影响，因而，所发掘出的关键路段不具有物流属性，亦缺乏针
对性和有效性。
发明内容
[0006] 鉴于上述技术问题，本发明的目的在于提供一种基于GPS数据的工业园区物流关键路段识别方法，该方法在识别时考虑城市建成环境因素、物流企业分布特征和物流车辆
的时空特征，对物流关键路段的识别更具针对性和准确性，而且通过实时读取的物流车辆
的GPS数据，便可自动识别物流关键路段，保持了较高的实时度。
[0007] 为实现上述目的，本发明是采用如下技术方案实现的：
[0008] 一种基于GPS数据的工业园区物流关键路段识别方法，该方法包括以下步骤：
[0009] 步骤1、根据供应商分布和工业园区的交通影响范围，划定电子围栏；
[0010] 步骤2、获取电子围栏内供应商及工业园区的物流信息，包括供应商的经纬度坐标、供给物流量、工业园区入口的经纬度坐标、需求物流量；
[0011] 步骤3、预处理电子围栏范围内网络拓扑结构数据，使格式统一，所述网络拓扑结构数据包括围栏范围内的节点信息和路段信息；所述节点信息包括节点编号、节点经纬度
坐标、节点各相位绿信比；所述路段信息包括路段编号、路段起止点、路段等级、路段长度、
路段车道数、路段设计通行速度；
[0012] 步骤4、获取北斗GPS数据，并进行预处理，数据包括物流车辆ID信息、物流车辆位置的经纬度信息与相应的时间戳、物流车辆与正北方向夹角、物流车辆所处位置海拔、物流
车辆瞬时速度；
[0013] 步骤5、根据步骤3的节点信息和路段信息计算节点拓扑结构参数，包括节点连接度和节点的通行能力；
[0014] 节点i的连接度的计算公式如下：
[0015]
[0016]
[0017]
[0018]
[0019] 其中，ki为节点i的连接度， 为第i个节点的入度，指与节点i相连且方向指向该节点的边的数目； 为第i个节点的出度，指与节点i相连且由该节点指向其他节点的边
的数目，N为节点数量；
[0020] 节点i的通行能力CAPi的计算公式如下：
[0021]
[0022] 其中，CAPi为节点i的通行能力；CAPi‑a为第a条进口车道的通行能力；Si‑a为第a条进口车道的规划饱和流量，按照国标标准并根据路段等级确定；λi‑a为第a条进口车道对应
的相位绿信比；
[0023] 步骤6、根据步骤2的供应商和工业园区的物流信息以及步骤5的节点通行能力计算节点的物流吸引力，
[0024] 节点i的物流吸引力的计算公式如下：
[0025]
[0026]
[0027]
[0028] 其中，fi‑A为节点i对供应商的吸引力，fi‑Γ为节点i对工业园区入口的吸引力fi‑Γ，fi为节点i的物流吸引力，CAPi为节点i的通行能力，di‑A为节点i到供应商A的距离；di‑Γ为节
点i到工业园区入口Γ的距离；WA为供应商A的供给物流量，当某供应商的供给物流量未知
时，取所有物流量的平均值计算，ω1为对供应商A的节点i的物流吸引力的权重，ω2为对工
业园区入口Γ的节点i的物流吸引力的权重，ω1+ω2＝1；
[0029] 步骤7、根据步骤3的路段信息和步骤5的节点连接度计算路段拓扑结构参数，包括路段连接度、路段权值和路段拓扑值；
[0030] 路段a的连接度Ka的计算公式如下：
[0031] Ka＝ki+kj
[0032] 其中，节点i和节点j分别为路段的起止点，ki和kj分别为两节点的连接度；
[0033] 路段a的拓扑值Ta的计算公式如下：
[0034] Ta＝β1·Ka·Da+β2·La
[0035] 其中，La为路段长度，β1为路段连接度Ka和路段权值Da乘积的权重，β2为路段长度的权重，β1+β2＝1，Da根据路段等级确定；
[0036] 步骤8、根据步骤4的GPS数据计算路段的动态平均速度，之后判断给定粒度的时间窗下是否有多辆车经过路段a，若是，则计算该时间窗下路段a的所有动态平均速度的平均
值，记录为
[0037] 路段a的动态平均速度Va的计算公式如下：
[0038]
[0039] 其中，dij为路段a的节点i到节点j的距离，ti为令x＝1时，GPS数据的第x行的轨迹位置位于节点i处的时间戳，tj为令x＝x+1时，GPS数据的第x行的轨迹位置位于节点j处的
时间戳；
[0040] 步骤9、根据步骤3的路段等级和步骤8路段的所有动态平均速度的平均值，确定路段a的动态服务水平Sa；
[0041] 步骤10、读取步骤4的GPS数据中的全部轨迹点，进行DBSCAN聚类分析,记录聚类分析所得簇数，根据簇中点的经纬度，确定每个簇所包含的路段，之后计算路段a的动态交通
量Qa，计算公式如下：
[0042]
[0043]
[0044] 其中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σa为路段a对应大车*ε
比，D 为每个簇平均邻域密度参数列表， 为第K个簇的平均邻域密度；
[0045] 步骤11、根据步骤6节点的物流吸引力计算路段a的物流吸引力，计算公式如下：
[0046] Fa＝fi+fj
[0047] 其中，节点i和节点j分别为路段a的起止点，fi和fj分别为两节点的物流吸引力；
[0048] 步骤12、根据步骤7的路段拓扑值、步骤9的路段动态服务水平以及步骤10的路段动态交通量，计算路段的动态通行能力，计算公式如下：
[0049]
[0050] 其中，Ea表示路段a的动态通行能力， 表示路段a的拓扑值Ta所占权重，
[0051] 表示路段a的服务水平Sa所占权重， 表示路段动态交通量Qa所占权重，且
[0052] 步骤13、计算路段的动态关键度，并按降序排序，计算公式如下：
[0053]
[0054] 其中，Ra表示路段a的动态关键度，λ1表示路段a的物流吸引力Fa所占权重； 表示路段a物流吸引力的排序等级；λ2表示路段动态通行能力Ea所占权重， 表示路段a的动态
通行能力排序等级，λ1+λ2＝1；
[0055] 步骤14、识别关键路段数量，根据路段动态关键度顺序，动态识别物流关键路段，关键路段数量计算公式如下：
[0056] Lα＝α·L
[0057] 其中，Lα为识别关键路段数量，α为预设关键路段比例，L为网络路段总数。
[0058] 作为本发明的优选，所述路段等级与路段权值及路段动态服务水平关系如下表：
[0059]
[0060] 作为本发明的优选，节点i到供应商a的距离di‑A的计算公式如下：
[0061]
[0062] 其中，节点i的经纬度坐标为(loni，lati)，供应商A的经纬度坐标为(lonA，latA)，R为地球半径，节点i到工业园区入口的距离di‑Γ以及节点i到节点j的距离dij计算过程同上。
[0063] 作为本发明的优选，步骤10路段a的动态交通量Qa的获取方法包括以下步骤：
[0064] 步骤10.1、读取步骤4的GPS数据中的全部轨迹点；
[0065] 步骤10.2、计算样本内每个点与所有点之间的距离dmn；
[0066] 步骤10.3、生成距离分布矩阵Dn×n；
[0067] 步骤10.4、将距离分布矩阵Dn×n每行元素分别按照升序排序，得到矩阵
[0068] 步骤10.5、计算K‑平均最近距离R参数列表，计算公式如下：
[0069]
[0070] 其中，R为邻域半径，DR为R对应的参数列表；K为矩阵 所对应列数， 为对应列中所有元素的均值；
[0071] 步骤10.6、对于给定的R，依次计算邻域密度阈值ε，计算公式如下：
[0072]
[0073] 其中，Pm为第m个点在给定邻域半径R下的邻域对象数量；
[0074] 步骤10.7、生成邻域密度阈值ε参数列表Dε(1≤K≤n)；
[0075] 步骤10.8、给定初始K＝1；
[0076] 步骤10.9、将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0077] 步骤10.10、记录聚类分析所得簇数MK；
[0078] 步骤10.11、判断簇数MK是否连续3次相等，若是，则转步骤10.13，若否，则转步骤10.12；
[0079] 步骤10.12、令K＝K+1，并转步骤10.9；
[0080] 步骤10.13、记录最优簇数M＝MK，并令K＝K+1；
[0081] 步骤10.14、将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0082] 步骤10.15、判断聚类分析所得簇数MK是否等于M，若是，则转步骤10.16，若否，则转步骤10.17；
[0083] 步骤10.16、令K＝K+1，并步骤10.14；
[0084] 步骤10.17、令K＝K‑1，并将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0085] 步骤10.18、根据簇中点的经纬度，确定每个簇所包含的路段；
[0086] 步骤10.19、计算每个簇平均邻域密度参数列表，计算公式如下：
[0087]
[0088] 其中， 为第K个簇的平均邻域密度；
[0089] 步骤10.20、计算路段动态交通量Qa，计算公式如下：
[0090]
[0091] 其中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σa为路段a对应大车比。
[0092] 本发明的优点和有益效果：
[0093] (1)本发明提供的关键路段识别方法只需获取物流车辆的GPS数据，即可自动识别物流关键路段，无需线圈传感器、录像设备等硬件设备，也不必进行传统录像、观测、记录，
即可为交通精细化组织提供决策依据。
[0094] (2)本发明自动更新北斗GPS历史数据并实时读取，通过分析海量数据，实现对物流关键路段的实时识别，保持了较高的实时度和准确性。
[0095] (3)本发明通过动态数据特征提取和自适应调整的DBSCAN聚类算法，自动识别物流关键路段，防止人工调整参数和后期维护产生误差，更具准确性。
[0096] (4)本发明提供的关键路段识别方法考虑了城市建成环境因素、物流企业分布特征和物流车辆的时空特征，对物流关键路段的识别更具针对性和准确性。
附图说明
[0097] 图1为本发明工业园区物流关键路段识别方法流程图。
具体实施方式
[0098] 为使本发明要解决的技术问题、技术方案更加清楚，下面将结合附图及具体实施案例进行详细描述。应当理解，此处所描述的实施例仅用于说明和解释本发明，并不用于限
定本发明。
[0099] 参阅图1，本发明提供的一种基于GPS数据的工业园区物流关键路段识别方法，该方法包括以下步骤：
[0100] 步骤1、根据供应商分布和工业园区的交通影响范围，划定电子围栏；
[0101] 步骤2、获取电子围栏内供应商及工业园区的物流信息，包括供应商的经纬度坐标、供给物流量(通过供应商获得)、工业园区入口的经纬度坐标、需求物流量(通过供应商
获得)；
[0102] 步骤3、预处理电子围栏范围内网络拓扑结构数据，使格式统一(可以统一转换成84坐标系)，所述网络拓扑结构数据包括围栏范围内的节点信息和路段信息；所述节点信息
包括节点编号(根据区域内有多少个十字口确定)、节点经纬度坐标、节点各相位绿信比等；
所述路段信息包括路段编号、路段起止点、路段等级、路段长度、路段车道数、路段设计通行
速度等；
[0103] 步骤4、获取北斗GPS数据，并进行预处理，数据包括物流车辆ID信息，物流车辆位置的经纬度信息与相应的时间戳，物流车辆与正北方向夹角，物流车辆所处位置海拔，物流
车辆瞬时速度；
[0104] 其中，采集到的GPS数据往往存在缺失与误差，这时候需要进行数据筛选和数据补偿，同时还需要根据电子围栏范围，筛选GPS数据，剔除不在电子围栏范围内的GPS轨迹数
据；设置时间窗粒度θ，单位为min，一般取30min，将步骤4的GPS数据按照时间窗粒度划分为
不同的数据集合，之后根据分析需要，选择对应时间窗粒度的北斗GPS数据进行分析；
[0105] 步骤5、根据步骤3的节点信息和路段信息计算节点拓扑结构参数，包括节点连接度和节点的通行能力，具体过程包括步骤5.1至步骤5.7：
[0106] 步骤5.1、根据步骤3的路段起止点信息，确定节点集合的邻接矩阵；
[0107] 步骤5.2、根据步骤5.1的邻接矩阵计算所有节点的入度、出度和连接度；
[0108] 设有向路网中有N个节点，第i个节点的入度 指与节点i相连且方向指向该节点的边的数目，出度 指与节点i相连且由该节点指向其他节点的边的数目，节点i的连接
度ki指该节点入度与出度的和，连接度的计算公式如下：
[0109]
[0110]
[0111]
[0112]
[0113] 步骤5.3、判断节点i与节点j是否可达，若可达，则转步骤5.4，若不可达，则转步骤5.5；
[0114] 步骤5.4、根据步骤3的路段等级确定左转/直行/右转车道数，然后按照国标确定该路段各车道的规划饱和流量；
[0115] 步骤5.5、读取下一个节点j，并转步骤5.3；
[0116] 步骤5.6、判断计算次数是否等于节点i的连接度，若等于，则转步骤5.7，否则，转步骤5.5；
[0117] 步骤5.7、根据步骤3的节点各相位绿信比，计算节点的通行能力，计算公式如下：
[0118]
[0119] 式中，CAPi为节点(信号控制交叉口进口道)通行能力(单位为pcu/h)；CAPi‑a为第a条进口车道的通行能力；Si‑a为第a条进口车道的规划饱和流量；λi‑a为第a条进口车道对应
的相位绿信比；
[0120] 步骤5.8、读取下一个节点i，并转步骤5.3；
[0121] 步骤5.9、判断所有节点是否均参与计算，若是，则转步骤6，否则，转步骤5.8；
[0122] 步骤6、根据步骤2的供应商和工业园区的物流信息以及步骤5的节点通行能力计算节点的物流吸引力，具体过程包括步骤6.1至步骤6.10：
[0123] 步骤6.1、根据步骤2的供应商经纬度坐标，计算节点i到供应商A的距离di‑A和到工业园区Γ入口的距离di‑Γ；
[0124] 设节点i的经纬度坐标为(loni，lati)，供应商A的经纬度坐标为(lonA，latA)，两者间距离di‑A计算公式如下：
[0125]
[0126] 式中，R为地球半径，通常取6378.137(KM)，到工业园区的距离di‑Γ计算过程同上；
[0127] 步骤6.2、对供应商A，计算节点i对供应商的吸引力fi‑A，计算公式如下：
[0128]
[0129] 步骤6.3、对工业园区Γ，计算节点i对工业厂区入口的吸引力fi‑Γ，计算公式如下：
[0130]
[0131] 步骤6.4、判断所有节点是否均参与计算，若是，则转步骤6.6，否则，转步骤6.5；
[0132] 步骤6.5、读取下一个节点i，并转步骤6.1；
[0133] 步骤6.6、判断所有供应商是否均参与计算，若是，则转步骤6.8，否则，转步骤6.7；
[0134] 步骤6.7、读取下一个供应商A，并转步骤6.1；
[0135] 步骤6.8、根据步骤2中的供应商供给物流量，计算节点i的物流吸引力fi，计算公式如下：
[0136]
[0137] 式中，WA为供应商A的供给物流量，当某供应商的供给物流量未知时，取所有物流量的平均值计算，ω1为对供应商A的节点i物流吸引力的权重，ω2为对工业园区Γ的节点i
物流吸引力的权重；且ω1+ω2＝1；一般取ω1＝0.3，ω2＝0.7；
[0138] 步骤6.9、判断所有节点是否均参与计算，若是，则转步骤7，否则，转步骤6.10；
[0139] 步骤6.10、读取下一个节点i，并转步骤6.8；
[0140] 步骤7、根据步骤3的路段信息和步骤5的节点连接度计算路段拓扑结构参数，包括路段连接度、路段权值和路段拓扑值，具体过程包括步骤7.1至步骤7.5：
[0141] 步骤7.1、计算路段a的连接度Ka，计算公式如下：
[0142] Ka＝ki+kj
[0143] 式中，节点i和节点j分别为路段的起止点，ki和kj分别为两节点的连接度；
[0144] 步骤7.2、根据步骤3的路段等级，确定该路段的权值Da，确定方法如下：
[0145] 表1路段权值的确定方式
[0146]路段等级 1 2 3 4
路段名称 快速路 主干道 次干道 支路
路段权值 10 8 5 3
[0147] 步骤7.3、计算路段a的拓扑值Ta，计算公式如下：
[0148] Ta＝β1·Ka·Da+β2·La
[0149] 式中，La为路段长度，β1为路段连接度Ka和路段权值Da乘积的权重，β2为路段长度的权重，且β1+β2＝1，一般取β1＝β2＝0.5；
[0150] 步骤7.4、判断所有路段是否均参与计算，若是，则转步骤8，否则，转步骤7.5；
[0151] 步骤7.5、读取下一个路段a，并转步骤7.1；
[0152] 步骤8、计算路段的动态平均速度，具体过程包括步骤8.1至步骤8.15：
[0153] 步骤8.1、根据步骤4的一组GPS数据，令x＝1，读取GPS数据的第x行信息；
[0154] 步骤8.2、判断第x行的轨迹位置是否位于节点i处，若位于，则转步骤8.3，否则，转步骤8.14，判断方法为：
[0155] 第x行点的经度与节点i的经度差小于等于阈值μ，且第x行点的维度与节点i的维度差小于等于阈值μ，μ一般取0.002；
[0156] 步骤8.3、记录此时时间戳为ti；
[0157] 步骤8.4、令x＝x+1，并读取GPS数据的第x行信息；
[0158] 步骤8.5、判断第x行的轨迹位置是否位于节点i和节点j之间，若是，则转步骤8.6，否则，转步骤8.7，判断方法为：
[0159] 第x行的轨迹位置到节点i和节点j的距离之和与两节点间距离之差是否小于等于0.1，单位为km；
[0160] 步骤8.6、令x＝x+1，并读取GPS数据的第x行信息，并转步骤8.5；
[0161] 步骤8.7、判断第x行的轨迹位置是否位于节点j处，若位于，则转步骤8.8，否则，转步骤8.12；
[0162] 步骤8.8、记录此时时间戳为tj；
[0163] 步骤8.9、计算路段a(节点i和节点j之间)的动态平均速度Va，计算公式如下：
[0164]
[0165] 式中，dij为节点i和节点j之间的距离，计算公式参照节点i到供应商A的距离的计算公式；
[0166] 步骤8.10、判断给定粒度的时间窗下，是否有多辆车经过路段a，若是，则转步骤8.11，若否，则转步骤8.12；
[0167] 步骤8.11、计算该时间窗下，路段a的所有动态平均速度的平均值，记录为
[0168] 步骤8.12、判断所有节点j是否均参与计算，若是，则转步骤8.14，否则，转步骤8.13；
[0169] 步骤8.13、读取下一个节点j，并转步骤8.1；
[0170] 步骤8.14、判断所有节点i是否均参与计算，若是，则转步骤9，否则，转步骤8.15；
[0171] 步骤8.15、读取下一个节点i，并转步骤8.1；
[0172] 步骤9、确定路段的动态服务水平，具体过程包括步骤9.1至步骤9.3：
[0173] 步骤9.1、根据步骤3的路段等级和步骤8.11的路段的所有动态平均速度的平均值，确定路段a的动态服务水平Sa，确定方法如下：
[0174] 表2路动态段服务水平确定方式
[0175]
[0176] 步骤9.2、判断所有路段是否均参与计算，若是，则转步骤10，否则，转步骤9.3；
[0177] 步骤9.3、读取下一个路段a，并转步骤9.1；
[0178] 步骤10、计算路段的动态交通量，具体过程包括步骤10.1至步骤10.20：
[0179] 步骤10.1、读取步骤4的GPS数据中的全部轨迹点；
[0180] 步骤10.2、计算样本内每个点与所有点之间的距离dmn；
[0181] 步骤10.3、生成距离分布矩阵Dn×n；
[0182] 步骤10.4、将距离分布矩阵Dn×n每行元素分别按照升序排序，得到矩阵
[0183] 步骤10.5、计算K‑平均最近距离R参数列表，计算公式如下：
[0184]R
[0185] 式中，R为邻域半径，D为R对应的参数列表；K为矩阵 所对应列数， 为对应列中所有元素的均值；
[0186] 步骤10.6、对于给定的R，依次计算邻域密度阈值ε，计算公式如下：
[0187]
[0188] 式中，Pm为第m个点在给定邻域半径R下的邻域对象数量；
[0189] 步骤10.7、生成邻域密度阈值ε参数列表Dε(1≤K≤n)；
[0190] 步骤10.8、给定初始K＝1；
[0191] 步骤10.9、将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0192] 步骤10.10、记录聚类分析所得簇数MK；
[0193] 步骤10.11、判断簇数MK是否连续3次相等，若是，则转步骤10.13，若否，则转步骤10.12；
[0194] 步骤10.12、令K＝K+1，并转步骤10.9；
[0195] 步骤10.13、记录最优簇数M＝MK，并令K＝K+1；
[0196] 步骤10.14、将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0197] 步骤10.15、判断聚类分析所得簇数MK是否等于M，若是，则转步骤10.16，若否，则转步骤10.17；
[0198] 步骤10.16、令K＝K+1，并转步骤10.14；
[0199] 步骤10.17、令K＝K‑1，并将邻域半径 和邻域密度阈值 作为参数，进行DBSCAN聚类分析；
[0200] 步骤10.18、根据簇中点的经纬度，确定每个簇所包含的路段；
[0201] 步骤10.19、计算每个簇平均邻域密度参数列表，计算公式如下：
[0202]
[0203] 式中， 为第K个簇的平均邻域密度；
[0204] 步骤10.20、计算路段动态交通量Qa，计算公式如下：
[0205]
[0206] 式中，路段a属于第K个簇，节点i和j分别为路段a的起止点，σa为路段a对应大车比；
[0207] 步骤11、计算路段的物流吸引力，具体过程包括步骤11.1至步骤11.3：
[0208] 步骤11.1、计算路段a的物流吸引力Fa，计算公式如下：
[0209] Fa＝fi+fj
[0210] 式中，节点i和节点j分别为路段a的起止点，fi和fj分别为两节点的物流吸引力；
[0211] 步骤11.2、判断所有路段是否均参与计算，若是，则转步骤12，否则，转步骤11.3；
[0212] 步骤11.3、读取下一个路段a，并转步骤11.1；
[0213] 步骤12、计算路段的动态通行能力，具体过程包括步骤12.1至步骤12.3：
[0214] 步骤12.1、计算路段a的动态通行能力Ea，计算公式如下：
[0215]
[0216] 式中，Ea表示路段a的动态通行能力， 表示路段a的拓扑值Ta所占权重， 表示路段a的服务水平Sa所占权重， 表示路段动态交通量Qa所占权重，且 一
般取
[0217] 步骤12.2、判断所有路段是否均参与计算，若是，则转步骤13，否则，转步骤12.3；
[0218] 步骤12.3、读取下一个路段a，并转步骤12.1；
[0219] 步骤13、计算路段的动态关键度，并按降序排序，计算公式如下：
[0220]
[0221] 式中，Ra表示路段a的动态关键度，λ1表示路段a的物流吸引力Fa所占权重； 表示路段a物流吸引力的排序等级(按照降序排列的等级)；λ2表示路段动态通行能力Ea所占权
重， 表示路段a的通行能力排序等级(按照降序排列的等级)，λ1+λ2＝1，一般取λ1＝λ2＝
0.5；
[0222] 步骤14、识别关键路段数量，根据路段动态关键度顺序，动态识别物流关键路段，关键路段数量计算公式如下：
[0223] Lα＝α·L
[0224] 式中，Lα为识别关键路段数量，α为预设关键路段比例，一般取5％，L为网络路段总数。
