本发明涉及一种城市轨道交通线网客流OD动态估计方法。
背景技术：
随着城市轨道交通的快速发展，各大中城市逐步跨入网络化运营。线网结构的复杂化导致乘客出行行为随机性加强，在客流需求分布结构上呈现出动态化、复杂化特点，对轨道交通的运输组织带来很大挑战，亟需利用适当的需求估计模型来研究较短时间范围内客流需求在时间上和空间上的分布规律，以提高轨道交通的动态化运营管理水平和系统应变能力。
自上世纪80年代以来，国内外学者对于OD动态估计开展了大量的研究，形成了一系列OD动态估计模型。但既有的OD估计方法在轨道交通中的应用，主要还存在以下不足：①现有研究主要集中于道路交通领域，没有考虑轨道交通网络结构及客流特征，在轨道交通线网客流OD动态估计方面存在着精度低、运算效率不高的缺点；②既有模型中OD流与采集信息流之间的流量关系大多建立在断面流量容易获取的基础上，但轨道交通中的实时断面客流信息却很难获取，仅能获取进出站客流信息，所以基于断面采集流量的动态流量方程难以适用于轨道交通网络。因此，需要一种新的方法，结合轨道交通线网客流特征和采集数据信息，实现轨道交通线网客流的实时OD估计。
技术实现要素：
发明目的：针对现有技术中存在的不足，本发明提供一种城市轨道交通线网客流OD动态估计方法，该方法基于历史客流出行数据和轨道交通自动售检票系统终端设备实时上传的交易数据来估算当前实时的客流OD分布结构信息。
技术方案：为实现上述目的，本发明采用的技术方案为：
一种城市轨道交通线网客流OD动态估计方法，包括如下步骤：
(1-1)设定时间间隔Δt，按时间间隔Δt对每天的历史客流数据进行分段；统计每天各时间段内的储值票客流数据和单程票客流数据；
(1-2)采用移动平均法对单程票客流数据进行改进，得到改进后的单程票客流数据为：
式中，i≠j；q0ij(t)表示改进后的第t个时段内由站点i进站的乘客中最终从站点j出站的单程票乘客数量；表示由历史客流数据统计获取的第t-a个时段内由站点i进站的乘客中最终从站点j出站的单程票乘客数量；R表示移动平均的时段个数，且R＜t；
(1-3)计算式中，表示由历史客流出行数据统计获取的第t个时段内由站点i进站的乘客中最终从站点j出站的储值票乘客数量，qij(t)表示在第t个时段内由站点i进站的客流去往站点j出站的总乘客数量；
根据qij(t)构建全线网在第t个时间段内进站客流的OD分布矩阵A(t)和第t个时间段内的客流分流率矩阵B(t)：
其中，n为站点的总数；bij(t)为客流分流率，表示时段t内由i站进站的乘客中去往j站的客流占i站总进站客流的比例，且
将客流分流率矩阵B(t)转化成列向量的形式：
B(t)＝[b12(t)，b13(t)，…，b1n(t)，…，b21(t)，…，b2n(t)，…，bn(n-1)(t)]T (1.4)
(1-4)构建客流出站到达系数为：
其中，为客流出站到达系数，表示在第t-m个时段内从站点i出发且以车站j为目的地的OD客流中在时段t内到达目的车站j的比例，t≥muij(t)表示第t个时间段内由i站出发去往j站的乘客平均行程时间，表示第t个时间段内由i站出发去往j站的乘客平均行程时间的标准方差；fij(x)为概率密度函数，表示由i站出发去往j站的客流在时刻x到达j站的概率；
(1-5)基于实时客流数据建立OD客流分配比例和进出站流量间的约束方程：
qij(t-m)＝Ii(t-m)·bij(t-m) (1.7)
式中，Ii(t-m)为i站在第t-m个时段内的总进站乘客数量；qij(t-m)表示第t-m个时段内由站点i进站的客流去往站点j总乘客数量；Oj(t)表示车站j在第t个时段内的出站乘客数量M为线网内任意两个时间段间乘客行程时间最大跨越时段数；Vij(t)为建立流量约束方程时产生的出站量误差；
(1-6)以客流分流率作为状态变量，构建客流OD动态估计状态空间模型，包括状态转移方程1.9与观测方程1.10：
式(1.9)中，B(t)为实际客流分流率bij(t)组成的Rod×1维矩阵，Rod表示OD对的总数，Rod＝n×(n-1)；Bk(t)为由相同客流特征日条件下的前向第k周历史客流分流率组成的Rod×1维矩阵；F(t)与Gk(t)均为状态转移矩阵，表征系统的状态演变特征，是由权重系数γk得到的Rod×Rod维常量矩阵；W(t)为建立系统状态转移方程所产生的误差Wij(t)组成的白噪声矩阵；
式(1.10)中，Oj(t)与Ii(t-m)为实时进出站客流数据；O(t)为n×1维出站客流量矩阵；H(t)为客流出站到达矩阵，其随着时段动态变化，表征状态变量B(t)与观测变量O(t)间的相互关系，为n×Rod维矩阵；是由构建的Rod×1维矩阵，为包含当前时段以及前向M个连续时间段的客流分流率均值；V(t)为建立系统观测方程所产生的误差vij(t)组成的白噪声矩阵；
(1-7)采用卡尔曼滤波方法对客流OD动态估计状态空间模型进行求解，并采用标准化法对OD估计结果进行修正；根据修正后的OD估计结果建立指标，并用该指标检验构建的客流OD动态估计状态空间模型是否正确；若检验结果正确，则判定客流OD动态估计状态空间模型正确，输出客流OD动态估计状态空间模型的估计结果；若检验结果不正确，则重新设置客流OD动态估计状态空间模型的参数值，返回步骤(1-6)；重新设置的参数包括：移动平均的时段数R和权重系数γk。
进一步的，所述步骤(1-6)中构建客流OD动态估计状态空间模型的步骤为：
(2-1)建立相邻时间段间的客流分流率关系：
式中，为由相同客流特征日下前向第k周历史客流数据统计获取的第t个时间段的客流分流率；γk为权重系数，0≤γk≤1，用于衡量前向第k周历史客流信息的可靠性；wij(t)为正态分布高斯白噪声变量，用于表征构建状态转移方程时产生的状态转移误差；
(2-2)将相邻时间段间的客流分流率关系式转换为标准的矩阵形式，得到状态转移方程为：
式中，W(t)为建立系统状态转移方程所产生的误差wij(t)组成的白噪声矩阵，且W(t)～N(0，Q(t))，Q(t)为状态转移误差方差，表示建立状态转移方程时产生的误差方差，Q(t)的无偏估计表达式如下：
式中，Wk(t)表示相同客流特征日下的前向第k周第t个时段的历史状态转移误差，为p天的历史状态转移误差均值；
(2-3)将相邻时段内的客流分流率平均值近似代替为各时段的客流分流率，将式(1-8)的表达式转换为以下形式：
由上式可得状态空间模型的观测方程为：
式中，V(t)为观测方程误差矩阵，且V(t)～N(0，R(t))，R(t)为出站量误差方差矩阵，表示建立观测方程时产生的误差方差，R(t)的无偏估计表达式为：
式中，Vk(t)表示第k天第t个时段的观测方程历史误差矩阵，为连续n天的历史观测误差均值。
进一步的，所述步骤(1-7)中采用卡尔曼滤波方法对客流OD动态估计状态空间模型进行求解并采用标准化法对OD估计结果进行修正，得出最优估计值的步骤为：
(3-1)定义协方差矩阵为P(t)；初始化t＝1；定义P(1)＝[1]n×n；
其中，Bk(1)为由相同客流特征日条件下的前向第k周历史客流分流率组成的Rod×1维矩阵；
(3-2)根据状态转移方程进行先验估计：
式中，表示第t个时段的状态变量B(t)的先验估计值，表示第t-1个时段的状态变量B(t-1)的后验估计值；
(3-3)计算先验估计协方差矩阵；
其中，表示第t个时段的先验估计协方差矩阵；表示第t-1个时段的后验估计协方差矩阵；
(3-4)计算卡尔曼滤波增益；
(3-5)根据卡尔曼滤波增益及估计值与观测值的残差修正先验估计值得到后验估计值即基于状态空间模型的标准卡尔曼滤波方法的解：
(3-6)标准卡尔曼滤波方法估计值约束修正：在OD动态估计过程中，状态变量B(t)需满足等式(1.21)约束，式(1.21)为：
以均方误差最小化为目标函数，可以得到：
其中|| ||表示向量的二范数；为由修正后客流分流率组成的Rod×1维矩阵，它是卡尔曼滤波标准步骤所得估计值的基础上，经过均方误差法调整后的所得修正估计值组成的向量；等式为修正后的状态向量应满足的等式约束方程，Y为n×1维常量矩阵，其元素值均为1；X为n×Rod维矩阵；
为式(1.22)表示的约束问题构建拉格朗日条件函数，可以得到：
其中，Z为构建的拉格朗日条件函数；β为拉格朗日乘子向量；P(B(t)|O(t))为条件概率密度函数；
假设初始系统状态变量B(1)、W(t)、V(t)为联合高斯变量，结合卡尔曼滤波的性质：当B(1)、W(t)、V(t)为联合高斯变量时，卡尔曼滤波估计值是条件为O(t)的B(t)的条件均值，可以得到：
然后分别对式(1.23)中的和β进行一阶求导，可以解得：
(3-7)更新后验估计协方差矩阵；
进一步的，所述步骤(1-7)中根据修正后的OD估计结果建立指标，并用该指标检验构建的客流OD动态估计状态空间模型是否正确的步骤为：
(4-1)构建标准化均方根误差指标：
RMSN值越低，表明估计模型越准确；
(4-2)判断RMSN的值是否小于预设的阈值RMSNmin，若满足RMSN＜RMSNmin，则判定所述客流OD动态估计状态空间模型正确；否则，判定所述客流OD动态估计状态空间模型不正确。
进一步的，所述步骤(1-7)中，重新设置参数的方法为：
计算：
式中，为预先设置的时段数迭代步长，且为整数；τ为预先设置的权重系数迭代步长，τ＜1。
需要说明的是，当k较小时，表明其越靠近实际客流数据，其γk的值在重新设置时，可使其增加，反之，当k较大时，其γk的值在重新设置时，可使其减小。
若通过检验，则确定所建立的线网客流OD动态估计方法是有效的，可以用于实际轨道交通运营管理。
有益效果：本发明提供的城市轨道交通线网客流OD动态估计方法，结合了历史客流的OD分布规律，采用移动平均法对波动性较大的单程票客流进行处理，提取出改进的历史客流分流率矩阵，能够提高模型估计精度；构建了客流出站到达系数，并提出其二重积分计算方法，用于建立OD流与进出站客流间的动态流量关系，解决了轨道交通断面客流难以采集导致流量关系式大大减少的问题；构建了基于卡尔曼滤波算法的状态空间模型，对线网客流OD进行动态估计，采取标准化法对OD估计结果进行修正，进一步提高了OD估计精度；最后提出采用标准化加权均方根误差法对模型的有效性进行检验，保证估计结果的可靠性。该方法能够有效地对城市轨道交通线网客流进行实时OD估计，为城市轨道交通运营管理决策提供数据支持。
附图说明
图1为本发明的方法操作流程示意图。
具体实施方式
如图1所示为一种城市轨道交通线网客流OD动态估计方法操作流程示意图，下面结合实施过程对本发明做出进一步的说明。
本实施例中，首先设定时间间隔Δt，按时间间隔Δt对运营时间进行分段，由于客流数据的实时上传一般为15min间隔，因此，为满足运营管理的需求，可设定时间间隔Δt＝15min。将历史客流出行数据分为储值票客流数据和单程票客流数据，统计储值票和单程票客流数据中的分时段OD分布客流量，考虑到单程票客流的空间分布稳定性较差且乘客到站时间差异性较大，考虑采用移动平均法对单程票OD客流进行改进。在每一个时间段，执行以下步骤：
Step1：计算客流OD分布矩阵：
其中，q0ij(t)的下标0表示单程票，下标i和j为站点编号，q0ij(t)表示通过移动平均法改进后的在第t个时段内由站点i进站的乘客中最终从站点j出站的单程票乘客数量；表示由原始客流数据统计获取的第t-a个时段内由站点i进站的乘客中最终从站点j出站的单程票乘客数量；R表示移动平均的时段个数，且R＜t。
将经过移动平均法改进后的单程票客流量与储值票客流量累加，可以得到全线网当天在第t个时间段内进站客流的OD分布矩阵A(t)为：
其中，n为站点的总数；qij(t)表示在第t个时段内由站点i进站的客流去往站点j出站的总乘客数量，且i≠j；
移动平均的时段数R的取值应根据历史单程票客流数据计算出的相同特征日的客流分布波动性大小而确定。
Step2：计算客流分流率矩阵：基于OD分布矩阵A(t)可得第t个时段内的客流分流率矩阵B(t)为：
其中，bij(t)为客流分流率，表示时段t内由i站进站的乘客中去往j站的客流占i站总进站客流的比例，可知
为便于后续模型的构建，将客流分流率矩阵B(t)转化成列向量的形式，如下：
B(t)＝[b12(t)，b13(t)，…，b1n(t)，…，b21(t)，…，b2n(t)，…，bn(n-1)(t)]T
Step3：构建客流出站到达系数：由于轨道交通AFC系统记录了乘客出行信息，同时，轨道交通乘客行程时间具有可靠性高的特点(地铁区间行驶时间基本固定)。因此，可通过各个OD间的乘客行程时间分布情况，用于刻画进出站客流和OD流间的流量关系。假设乘客在任意时段t内第x分钟由i站出发去往j站的客流行程时间服从正态分布，即x∈N(uij(t)，其中，uij(t)表示时段t内由i站出发去往j站的乘客平均行程时间；表示时段k内，乘客由i站出发去往j站的平均旅行时间的标准方差，uij(t)和可以通过历史客流OD数据进行统计分析得到。
在此基础上，假设在任意时段t内由站点i进站去往站点j的客流在该时段内均匀到达站点i，可以通过各个OD间的乘客行程时间概率密度函数fij(x)，并对第t-m个时段内由i站出发去往j站的客流在时段t内到达j站的概率进行积分，从而计算出客流出站到达系数来刻画进出站客流和OD流间的流量关系，其表达式为：
其中，为客流出站到达系数，表示在第t-m个时段内从站点i出发且以车站j为目的地的OD客流中在时段t(t≥m)内到达目的车站j的比例(某时段采集到的出站客流信息是前期多个时段由其它车站进站客流到达聚集的结果)；
Step4：构建动态流量关系：由于轨道交通断面客流难以采集，难以借鉴道路交通领域以实时采集路段断面客流为基础从而建立动态流量关系，为此，本方法基于时变的进出站采集客流量，建立OD客流分配比例和进出站流量间的约束方程：
qij(t-m)＝Ii(t-m)·bij(t-m)
其中，n表示线网车站总数；Ii(t-m)为车站i在第t-m个时段内的总进站客流量；qij(t-m)表示第t-m个时段内由站点i进站的客流去往站点j出站的总乘客数量；bij(t-m)为第t-m个时段的客流分离率，表示第t-m个时段内由i站出发去往j站的客流占车站i总进站客流的比例；Oj(t)表示车站j在时段t内的出站量；M为路网内任意OD间乘客行程时间最大跨越时段数，最大跨越时段数M的值取决于当地城市轨道交通的线网规模。；为客流出站到达系数；Vij(t)为建立观测方程时产生的出站量误差。
Step5：建立状态空间模型：选取车站客流分流率作为估计变量，建立线网客流OD动态估计状态空间模型。鉴于短时内客流波动不会太大，可得出相邻时段间客流分流率之间的关系满足：
式中，为由相同客流特征日下前向第k周历史客流数据统计获取的第t个时间段的客流分流率；γk为权重系数，0≤γk≤1，用于衡量前向第k周历史客流信息的可靠性；wij(t)为正态分布高斯白噪声变量，用于表征构建状态转移方程时产生的状态转移误差。
将相邻时间段间的客流分流率关系式转换为标准的矩阵形式，得到状态转移方程为：
其式中，W(t)为建立系统状态转移方程所产生的误差wij(t)组成的白噪声矩阵，且W(t)～N(0，Q(t))，Q(t)为状态转移误差方差，表示建立状态转移方程时产生的误差方差，Q(t)的无偏估计表达式如下：
式中，Wk(t)表示相同客流特征日下的前向第k周第t个时段的历史状态转移误差，为p天的历史状态转移误差均值。
由相邻时段间客流分流率之间的关系式可知，估计过程需包含多个时段的状态变量信息，估计第t个时段的客流分流率bij(t)与bij(t-1)，…，bij(t-m)有关，那么，需要将多个时段的变量信息整合为一个时段，从而方便模型构建。为此，假定在一定时间范围内客流波动变化较小，采用相邻时段内的客流分流率均值来将多个时段的bij(t)整合为一个时段，可将相邻时段间客流分流率之间的关系式转化为以下形式：
将上式化为标准的矩阵形式，可得系统观测方程为：
其中，O(t)为n×1维出站量矩阵；H(t)为由确定的客流出站到达矩阵，其随着时段动态变化，表征状态变量与观测变量间的相互关系，为n×Rod维矩阵；是由连续时间段内客流分流率均值构建的Rod×1维矩阵；V(t)为出站量方程系统误差矩阵，且V(t)～N(0，R(t))，R(t)为出站量误差方差矩阵，表示建立观测方程时产生的误差方差，可以根据历史数据中的观测误差方差样本值V(t)进行统计分析得到，其无偏估计表达式如下：
Vk(t)表示第k天第t个时段的历史观测误差，为连续n天的历史观测误差均值。
Step6：卡尔曼滤波算法求解：卡尔曼滤波算法是状态空间模型的经典求解方法，它实际上是一种最优化的自回归数据处理算法，它假设任意时段t的状态向量B(t)存在两种估计值，即先验估计值以及后验估计值卡尔曼滤波算法的基本思想是：任意时段t的OD估计值均是在先验估计的基础上，经过系统观测值O(t)的进一步修正而求得的，而且时段t的先验估计值总是以t-1时段的后验估计值为基础进行改进。其迭代递推具体步骤为：
1)系统初始化：定义协方差矩阵为P(t)；初始化t＝1；定义P(1)＝[1]n×n；其中，初始时段t＝1为运营日的第1个估计时间段，也即[T0，T0+Δt]，T0表示地铁运营日的起始时间；初始OD客流分流率B(1)可以理解为第1个估计时段的后验客流分流率，是卡尔曼滤波迭代算法的基础，但B(1)却难以通过实时采集数据得到，通常利用历史数据中连续k周的初始时段客流分流率矩阵均值代替；Bk(1)为由相同客流特征日条件下的前向第k周历史客流分流率组成的Rod×1维矩阵；
初始协方差矩阵P(1)可以理解为第1个估计时段的后验协方差矩阵，由于P(1)难以通过实时采集数据或历史数据统计分析得到，为此，可将P(1)设为单位矩阵。
2)根据状态转移方程进行先验估计；
其中，表示时段t的状态变量B(t)的先验估计值，表示时段t-1的状态变量B(t-1)的后验估计值；
3)计算先验估计协方差矩阵；
其中，表示时段t的先验估计误差方差；表示表示时段t-1的后验估计误差方差；
4)计算卡尔曼滤波增益；
5)根据卡尔曼滤波增益及估计值与观测值的残差修正先验估计值得到后验估计值即基于状态空间模型的标准卡尔曼滤波方法的解：
6)更新后验估计协方差矩阵；
以上步骤为卡尔曼滤波方法的标准步骤。需要说明的是，在轨道交通线网客流OD动态估计过程中，状态变量B(t)的估计值需满足等式约束为此，可以在步骤6)执行前，采用标准化法对步骤5)得到的OD估计结果进行修正：
式中，为状态变量客流分流率bij(t)的后验估计值，为采用标准化法进行修正后的客流分流率估计值。
得到修正结果的具体步骤为：
以均方误差最小化为目标函数，可以得到：
其中|| ||表示向量的二范数；为由修正后客流分流率组成的Rod×1维矩阵，它是卡尔曼滤波标准步骤所得估计值的基础上，经过均方误差法调整后的所得修正估计值组成的向量；等式为修正后的状态向量应满足的等式约束方程，Y为n×1维常量矩阵，其元素值均为1；X为n×Rod维矩阵，其元素取值如下：
上式中，zeros函数为产生元素值全为0的矩阵。
为目标函数构建拉格朗日条件函数，可以得到：
其中，Z为构建的拉格朗日条件函数；β为拉格朗日乘子向量；P(B(t)|O(t))为条件概率密度函数。
这里，我们假设初始系统状态变量B(1)、W(t)、V(t)为联合高斯变量，结合卡尔曼滤波的性质：当B(1)、W(t)、V(t)为联合高斯变量时，那么卡尔曼滤波估计值是条件为O(t)的B(t)的条件均值，可以得到：
然后分别对式(1.24)中的和β进行一阶求导，可以解得：
于是，通过以上修正步骤就实现了标准卡尔曼滤波算法估计结果的修正。在每计算出某一个时段t的客流分流率估计值之后，则利用以上方法对其进行修正，以保证客流分流率估计值满足等式约束。
Step7：估计方法检验：使用样本数据对Step6得到的估计值进行检验，若未能通过检验，则返回前述步骤中，重新设置移动平均的时段数R、权重系数α等参数的值；若通过检验，则确定所建立的时间序列预测模型是有效的，可以用于实际轨道交通运营管理。
采用标准化加权均方根误差指标(Weighted Root Mean Square Error Normalized，WRMSN)评估估计模型的有效性，WRMSN指标表达式如下：
其中，Rod表示网络OD总数；n为车站总数；bij(t)和分别为时段t内的线网客流分流率真实值和线网客流分流率估计均值；RMSN值越低，表明估计模型越准确。
当根据样本数据实际值与估计值所计算出的RMSN指标数据在允许的范围内(如RMSN≤20％)时则认为方法是可行的，可以应用于实际的轨道交通运营管理；若RMSN值过大则重新设置移动平均的时段数R以及权重系数γk的值，重复建立估计模型的步骤直到模型通过有效性检验。
重新设置参数的方法为：
计算：
式中，为预先设置的时段数迭代步长，且为整数；τ为预先设置的权重系数迭代步长，τ＜1。
需要说明的是，当k较小时，表明其越靠近实际客流数据，其γk的值在重新设置时，可使其增加，反之，当k较大时，其γk的值在重新设置时，可使其减小。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。