CN111179595A - 一种交通流时间序列分解方法 - Google Patents

Abstract

本发明公开了一种交通流时间序列分解方法，它包括以下步骤：(1)采集所需要分析道路断面的交通流历史数据；(2)将历史交通流序列分解为趋势、周期和残差3个成分；(3)将新观测的交通流序列分解为趋势、周期和残差3个成分。该方法充分考虑了交通流序列的周期性，可以剥离出随机噪声的影响，提取交通流周期波动以及趋势变化的规律。该方法适用于具有周期性的交通流状态序列，分解结果可用于交通状态分析、检测与预测，辅助智能交通系统的研究。

Description

一种交通流时间序列分解方法

技术领域

本发明涉及智能交通控制技术领域，具体为一种交通流时间序列分解方法。

背景技术

交通流序列是通过道路检测设备获得的第一手交通状态数据，对交通流序列分析一直以来都是交通工程领域的研究热点。运用科学的方法与技术手段，研究分析交通流序列的变化特征，可以挖掘出交通流的增长趋势、时段瓶颈和异常状态等，一方面可以为道路使用者提供合理的路径规划方案，另一方面，可以为政府道路交通管理部门制定管理与控制方案提供决策依据。

现阶段研究普遍认为现实生活中的交通流序列是由多种信号源叠加产生的复合结果，一方面，交通流序列在短期内不稳定、波动较大，而另一方面，交通流在长期内又呈现周期变化的规律。而时间序列分解技术则是通过特定的方法把时间序列分解为多个简单的、具有一定规律的时间序列分量。根据交通流周期变化的特征，对交通流序列进行时间序列分解，可以剥离出随机噪声的影响，提取交通流周期波动规律。因此，对交通流序列分解可以分析交通流的长期变化趋势、周期波动、异常状态等，辅助交通流运行状况监测、交通流预测等智能交通技术的实施。

发明内容

为了可以将交通流序列分解为趋势、周期和残差3个成分；并提供一种动态分解方法，基于已建立的分解模型，对新观测的交通流数据进行动态分解，本发明提供一种交通流时间序列分解方法，可以将交通流序列分解为趋势、周期和残差3个成分；并提供了一种动态分解方法，基于已建立的分解模型，对新观测的交通流数据实现动态分解。该方法可以用于分析交通流序列的整体趋势、周期波动以及异常状况，为达此目的，本发明提供一种交通流时间序列分解方法，其特征在于：所述方法可以将交通流序列分解为趋势、周期和残差3个成分，并附有一种动态分解方法，基于已建立的分解模型，对新观测的交通流数据进行动态分解，

(1)采集所需要分析的道路断面的交通流历史数据；

(2)将历史交通流序列分解为趋势、周期和残差3个成分；

(3)将新观测的交通流序列动态分解为趋势、周期和残差3个成分。

作为本发明进一步改进，所述交通流历史数据为交通流量、车辆平均速度或其他交通状态指标所组成的时间序列。

作为本发明进一步改进，步骤二将历史交通流序列分解为趋势、周期和残差3个成分具有以下步骤：

(1)定义所采集的交通流序列为Y(t)，周期长度为C，周期数量为M，则交通流序列Y(t)的样本数量为MC，初始化趋势分量T(t)＝0；

(2)去趋势，用原始交通流序列Y(t)减去趋势分量T(t)，得到去趋势的时间序列

即：

(3)周期子序列平滑，对每个周期子序列作K为M/2的LOESS，并分别向前和向后各延展一个周期；

即表示对于自变量为X＝[1,2,3,4,..,M]^Τ，因变量为周期子序列Y＝[Y(c),Y(C+c),Y(2C+c),...,Y((M-1)C+c))]^T，基于K＝M/2的LOESS求X＝[0,1,2,3,4,..,M,M+1]^Τ时对应的估计值

再把所有周期子序列估计值

按照时间先后重组成临时时间序列

上述周期子序列S_c(t)＝{Y(c),Y(C+c),Y(2C+c),...,Y((M-1)C+c))},c＝1,2,...,C表示由每个周期同一位置的数据所组成的时间序列；

(4)周期子序列的低通量过滤，对临时时间序列

依次做长度为C、C和3的移动平均，然后再做K＝C的LOESS，得到时间序列L(t)，t＝1,2,...,MC；

(5)去平滑周期子序列趋势，用

减去L(t)，即：

计算周期分量，依据F(t)，首先计算单个周期的周期分量Q(t),t＝1,2,...,C：

然后再把将Q(t)扩展M个周期，得到周期分量P(t)，t＝1,2,...,MC

(6)去周期，用原始交通流序列Y(t)减去周期分量P(t)，得到去周期的时间序列

即：

计算趋势分量，对去周期的时间序列

做K＝C的LOESS，得到趋势分量T(t)；

(8)终止条件判断，若满足终止条件，则执行下一步骤(9)，否则继续执行先前步骤(2)至(7)；

(9)计算残差分量，用原始序列Y(t)减去趋势分量T(t)和周期分量P(t)，得到残差分量R(t)，即：

R(t)＝Y(t)-T(t)-P(t)。

作为本发明进一步改进，步骤三将新观测的交通流序列动态分解为趋势、周期和残差3个成分，其具有以下步骤：

首先定义新观测数据Z(t),t＝MC+1,MC+2,...

(1)按照时间先后顺序遍历Z(t)中的所有的数据z_t，执行步骤(2)至(5)，直到所有的观测数据都完成了分解；

(2)计算周期分量，根据z_t所在周期的位置和权利要求3中步骤(5)得到的单周期分量Q(t)，计算其周期分量p_t；

(3)去周期，用原始观测值z_t减去周期分量p_t得到去周期数值

(4)计算趋势分量，把

附加

的尾部，并做K＝C的LOESS求该点的估计值，得到趋势分量t_t；

(5)计算残差分量，用原始观测序列z_t减去趋势分量t_t和周期分量p_t，得到残差分量r_t，即：

r_t＝z_t-t_t-p_t。

与现有技术相比，本发明的有益效果是：

本发明充分考虑了交通流序列的周期性，可以将交通流序列分解为趋势、周期和残差3个成分，每一成分都有具体的实际含义和分析价值，趋势成分反映了交通流状态的趋势变化，周期成分反映了交通状态的周期波动，残差成分可以监测交通状态的异常变化，分解结果可用于交通状态分析、监测与预测，辅助智能交通系统的研究。

本发明充分考虑了交通流的动态性，提供了动态分解方法，在已建立的分解模型的基础上，可以将实时观测数据动态分解为趋势、周期和残差3个成分，能够满足智能交通系统实时响应和动态更新的需求。

附图说明

图1为本发明的总体流程图；

图2为历史交通流序列分解的流程图；

图3为观测交通流序列动态分解的流程图；

图4为LOESS的示意图；

图5为周期子序列平滑和周期子序列重组的示意图；

图6为实施例1的分解结果；

图7为实施例2的分解结果。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

本发明提供一种交通流时间序列分解方法，可以将交通流序列分解为趋势、周期和残差3个成分；并提供了一种动态分解方法，基于已建立的分解模型，对新观测的交通流数据实现动态分解。该方法可以用于分析交通流序列的整体趋势、周期波动以及异常状况。

具体技术方案如下：

一种交通流时间序列分解方法，其包括以下步骤，具体总体流程图如图1所示：

(1)选择需要预测的道路断面，采集道路断面的交通流历史数据；

(2)通过时间序列分解方法，将历史交通流序列分解为趋势、周期和残差3个成分；

(3)通过动态分解方法，将新观测的交通流序列分解为趋势、周期和残差3个成分。

进一步的，步骤(1)中所述交通流历史数据可以为交通流量、平均速度或其他交通状态指标所组成的时间序列。

进一步的，步骤(2)中所述交通流时间序列分解方法，分解流程如附图2所示，具体步骤如下：

(a)定义步骤(1)中所采集的交通流序列为Y(t)，周期长度为C，周期数量为M，则交通流序列Y(t)的样本数量为MC，初始化趋势分量T(t)＝0，t＝1,2,...,MC；

(b)去趋势，用原始交通流序列Y(t)减去趋势分量T(t)，得到去趋势的时间序列

即：

(c)周期子序列平滑，如附图5所示，对每个周期子序列作K为M/2的LOESS回归，并分别向前和向后各延展一个周期，即：

对自变量为X＝[1,2,3,4,..,M]^T，因变量为周期子序列

Y＝[Y(c),Y(C+c),Y(2C+c),...,Y((M-1)C+c))]^T，基于LOESS求

X＝[0,1,2,3,4,..,M,M+1]^Τ时对应的估计值

再把所有周期子序列的LOESS估计值

按照时间先后重组成临时时间序列

其中，周期子序列S_c(t)＝{Y(c),Y(C+c),Y(2C+c),...,Y((M-1)C+c))},c＝1,2,...,C表示由每个周期同一位置的数据所组成的时间序列。

其中，LOESS英文全称为Locally Weighted Regression，中文翻译为局部加权回归，是一种非参数回归方法；如附图4所示，x_i的估计值

由距离x_i最近的K个点加权得到，K为LOESS需要指定的参数：

上式中w_k表示邻近点y_k的权重，计算公式如下：

上式中，D(x)为Epanechnikov核函数，函数形式如下：

λ为Epanechnikov核函数的核宽度，表示x_i与其第K近的点的x_[K]距离：

λ＝|x_i-x_[K]|

(d)周期子序列的低通量过滤，对步骤(c)得到的

依次做长度为C、C和3的移动平均，然后再做K＝C的LOESS回归，得到时间序列L(t)，t＝1,2,...,MC；

(e)去平滑周期子序列趋势，用步骤(c)得到的

减去步骤(d)得到的L(t)，即：

计算周期分量，依据步骤(e)得到的F(t)，首先计算单个周期的周期分量

Q(t),t＝1,2,...,C

然后把将Q(t)扩展C个周期，得到周期分量P(t)，t＝1,2,...,MC

(f)去周期，用原始交通流序列Y(t)减去步骤(e)得到周期分量P(t)，得到去周期的时间序列

即：

(g)计算趋势分量，对步骤(f)得到的去周期的时间序列

做K＝C的LOESS回归，得到趋势分量T(t)；

(h)终止条件判断，若满足终止条件，则执行下一步骤(i)，否则继续执行先前步骤(b)至(h)；

(i)计算残差分量，用原始序列Y(t)减去趋势分量T(t)和周期分量P(t)，得到残差分量R(t)，即：

R(t)＝Y(t)-T(t)-P(t)

通过步骤(a)至(i)，历史交通流序列可以Y(t)被分解为趋势序列T(t)、周期序列P(t)和残差序列R(t)3个成分相加。

进一步的，步骤(3)中对新观测数据Z(t),t＝MC+1,MC+2...的动态分解方法，具体步骤如下：

(a)按照时间先后顺序遍历Z(t)中的所有的数据z_t，执行步骤(b)至(e)，直到所有的观测数据都完成了分解。

(b)计算周期分量，根据z_t所在周期的位置和单周期分量Q(t)，计算其周期分量p_t；

(c)去周期，用原始观测值z_t减去周期分量p_t得到去周期数值

(d)计算趋势分量，把

附加到

的尾部，并做LOESS求该点的估计值，得到趋势分量t_t；

(e)计算残差分量，用原始观测序列z_t减去趋势分量t_t和周期分量p_t，得到残差分量r_t，即：

r_t＝z_t-t_t-p_t

通过步骤(a)至(e)，交通流的实时观测值z_t可以被分解为趋势分量t_t和周期分量p_t，得到残差分量r_t 3个部分相加。

实施例1：

(1)选择需要分析的道路断面，采集道路断面的交通流历史数据；

本实施例所采集的交通数据来源美国明尼苏达大学公开的明尼苏达州高速公路交通流数据集，采集的监测点为S195，采集的交通流序列为5分钟间隔的交通流量，采集的时间为2019年4月29日至5月6日工作日(每周的周一到周五)共20天的全天24小时交通流量。采集样本序列共5760个，其中前3周(共4320条)作为历史数据用于时间序列分解，后1周(共1440条)作为实时观测数据用于动态分解。所采集的样本如附图6中原始序列所示。

(2)通过时间序列分解方法，将历史交通流序列分解为趋势、周期和残差3个成分，实施流程如附图2所示，具体实施步骤如下：

(a)将步骤(1)中所采集的前4320条5分钟交通流量按照时间先后顺序排列组成时间序列Y(t)，t＝1,2,...,4320，则周期长度为C＝288，周期数量为M＝15，交通流序列Y(t)的样本数量为4320，初始化趋势分量T(t)＝0，t＝1,2,...,4320；

(b)去趋势，用原始交通流序列Y(t)减去趋势分量T(t)，得到去趋势的时间序列

即：

(c)周期子序列平滑，对每个周期子序列作K＝8(K＝M/2＝7.5，并四舍五入取整)的LOESS回归，并分别向前和向后各延展一个周期。

即对于自变量X＝[1,2,3,4,..,15]^Τ，因变量为周期子序列

Y＝[Y(c),Y(288+c),Y(2×288+c),...,Y(14×288+c))]^T，基于LOESS求

X＝[0,1,2,3,4,..,15,16]^Τ时对应的估计值

再把所有周期子序列的LOESS估计值

按照时间先后重组成临时时间序列

其中，周期子序列

S_c(t)＝{Y(c),Y(288+c),Y(2×288+c),...,Y(14×288+c))},c＝1,2,...,288指在每个周期同一位置组成的时间序列。

(d)周期子序列的低通量过滤，对步骤(c)得到的

依次做长度为288、288和3的移动平均，然后再做LOESS回归，得到时间序列L(t)，t＝1,2,...,4320；

(e)去平滑周期子序列趋势，用步骤(c)得到的

减去步骤(d)得到的L(t)，即：

计算周期分量，依据步骤(e)得到的F(t)，首先计算单个周期的周期分量Q(t),t＝1,2,...,288

然后把将Q(t)扩展15个周期，得到周期分量P(t)，t＝1,2,...,4320

(f)去周期，用原始交通流序列Y(t)减去周期分量P(t)，得到去周期的时间序列

即：

(g)计算趋势分量，对步骤(e)得到的去周期的时间序列

做K＝288的LOESS回归，得到趋势分量T(t)，t＝1,2,...,4320；

(h)循环终止条件判断，设置循环次数为3，若循环次数大于3，则执行下一步骤(i)，否则继续执行先前步骤(b)至(h)；

(i)计算残差分量，用原始序列Y(t)减去趋势分量T(t)和周期分量P(t)，得到残差分量R(t)，即：

R(t)＝Y(t)-T(t)-P(t)，t＝1,2,...,4320

(3)通过动态分解方法，将新观测的交通流序列分解为趋势、周期和残差3个成分，实施流程如附图3所示，具体步骤如下：首先记步骤(1)中所述后1周(共1440条)实时观测数据为Z(t),t＝4321,4322,...,5760

(a)按照时间先后顺序遍历Z(t)中的所有的数据z_t，t＝4321,4322,...,5760，执行步骤(b)至(e)，直到所有的观测数据都完成了分解。

(b)计算周期分量，根据z_t所在周期的位置和步骤(2)中子步骤(e)计算得到的单周期分量Q(t)，计算其周期分量p_t；

(c)去周期，用原始观测值z_t减去周期分量p_t得到去周期数值

(d)计算趋势分量，把

附加

的尾部，并做K＝288的LOESS求该点的估计值，得到趋势分量t_t；

(e)计算残差分量，用原始观测序列z_t减去趋势分量t_t和周期分量p_t，得到残差分量r_t，即：

r_t＝z_t-t_t-p_t

分解结果如附图6中趋势序列、周期序列和残差序列所示，其中附图6中实线表示对历史时间序列的分解，虚线表示对实时观测数据的分解。从附图6中趋势序列可以看出，采集时段内的交通流量在一定的范围内波动，没有明显的增加或减少的趋势；附图6中的周期序列表示交通流量在每一天内的周期变化，可以看出该交通流量序列具有很强的周期性，存在明显的早高峰与晚高峰，且早高峰交通流量高于晚高峰；附图6的残差序列表示噪声干扰或者是其他因素对交通流量的影响，当残差序列出现过大或过小的波动时，可以认为交通流量出现异常，并可以结合其他监测数据进行验证。

实施例1结束。

实施例2：

(1)选择需要预测的道路断面，采集道路断面的交通流历史数据；

本实施例所采集的交通数据来源美国明尼苏达大学公开的明尼苏达州高速公路交通流数据集，采集的监测点为S1924，采集的交通流序列为10分钟间隔的车辆平均速度，采集的时间为2019年7月22日至8月16日工作日(每周的周一到周五)共20天的上午时段(5:00至12：00)的车辆平均速度。采集样本序列共840个，其中前3周(共630条)作为历史数据用于时间序列分解，后1周(共210条)作为实时观测数据用于动态分解。所采集的样本如附图7中原始序列所示。

(2)通过时间序列分解方法，将历史交通流序列分解为趋势、周期和残差3个成分，实施流程如附图2所示，具体实施步骤如下：

(a)将步骤(1)中所采集的前630条10分钟车辆平均速度按照时间先后顺序排列组成时间序列Y(t)，t＝1,2,...,630，则周期长度为C＝42，周期数量为M＝15，则交通流序列Y(t)的样本数量为MC＝630，初始化趋势分量T(t)＝0，t＝1,2,...,630；

(b)去趋势，用原始交通流序列Y(t)减去趋势分量T(t)，得到去趋势的时间序列

即：

(c)周期子序列平滑，对每个周期子序列作K＝8(K＝M/2＝7.5，并四舍五入取整)的LOESS回归，并分别向前和向后各延展一个周期。

即对于自变量X＝[1,2,3,4,..,15]^Τ，因变量

Y＝[Y(c),Y(42+c),Y(2×42+c),...,Y(14×42+c))]^T，基于LOESS求

X＝[0,1,2,3,4,..,15,16]^Τ时对应的估计值

再把所有周期子序列的LOESS估计值

按照时间先后重组成临时时间序列

其中，周期子序列S_c(t)＝{Y(c),Y(42+c),Y(2×42+c),...,Y(14×42+c))},c＝1,2,...,42指在每个周期同一位置组成的时间序列，

(d)周期子序列的低通量过滤，对步骤(c)得到的

依次做长度为42、42和3的移动平均，然后再做LOESS回归，得到时间序列L(t)，t＝1,2,...,630；

(e)去平滑周期子序列趋势，用步骤(c)得到的

减去步骤(d)得到的L(t)，即：

计算周期分量，依据步骤(e)得到的F(t)，首先计算单个周期的周期分量

Q(t),t＝1,2,...,42

然后把将Q(t)扩展15个周期，得到周期分量P(t)，t＝1,2,...,630

(f)去周期，用原始交通流序列Y(t)减去周期分量P(t)，得到去周期的时间序列

即：

(g)计算趋势分量，对步骤(e)得到的去周期的时间序列

做K＝42的LOESS回归，得到趋势分量T(t)，t＝1,2,...,630；

(h)循环终止条件判断，设置循环次数为3，若循环次数大于3，则执行下一步骤(i)，否则继续执行先前步骤(b)至(h)；

(i)计算残差分量，用原始序列Y(t)减去趋势分量T(t)和周期分量P(t)，得到残差分量R(t)，即：

R(t)＝Y(t)-T(t)-P(t)，t＝1,2,...,630

(3)通过动态分解方法，将新观测的交通流序列分解为趋势、周期和残差3个成分，实施流程如附图3所示，具体步骤如下：首先记步骤(1)中所述后1周(共1440条)实时观测数据为Z(t),t＝631,632,...,840

(a)按照时间先后顺序遍历Z(t)中的所有的数据z_t，t＝631,632,...,840，执行步骤(b)至(e)，直到所有的观测数据都完成了分解。

(b)计算周期分量，根据z_t所在周期的位置和和步骤(2)中子步骤(e)计算得到的单周期分量Q(t)，计算其周期分量p_t；

(c)去周期，用原始观测值z_t减去周期分量p_t得到去周期数值

(d)计算趋势分量，把

附加到

的尾部，并做K＝64的LOESS求该点的估计值，得到趋势分量t_t；

(e)计算残差分量，用原始观测序列z_t减去趋势分量t_t和周期分量p_t，得到残差分量r_t，即：

r_t＝z_t-t_t-p_t

分解结果如附图7中趋势序列、周期序列和残差序列所示，其中附图7中实线表示对历史平均车速序列的分解，虚线表示对实时观测数据的分解。从附图7中趋势序列可以看出，每天上午(早上5:00至12:00)的车辆平均速度的增长趋势在72至75英里/小时之间波动；附图7中的周期序列表示车辆平均速度在每天上午的周期变化，可以看出每天上午的车辆速度随时间递减，接近中午时车辆速度最低，这与道路交通流量的变化相关；附图7的残差序列表示噪声干扰或者是其他因素对道路平均速度的影响，当残差序列出现过大或过小的波动时，可以认为道路平均速度出现异常，道路通行出现瓶颈，并可以结合其他监测数据进行验证。

实施例2结束。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。

Claims (4)

1.一种交通流时间序列分解方法，其特征在于：所述方法可以将交通流序列分解为趋势、周期和残差3个成分，并附有一种动态分解方法，基于已建立的分解模型，对新观测的交通流数据进行动态分解，

(1)采集所需要分析的道路断面的交通流历史数据；

(2)将历史交通流序列分解为趋势、周期和残差3个成分；

(3)将新观测的交通流序列动态分解为趋势、周期和残差3个成分。

2.根据权利要求1所述的一种交通流时间序列分解方法，其特征在于：所述交通流历史数据为交通流量、车辆平均速度或其他交通状态指标所组成的时间序列。

3.根据权利要求1所述的一种交通流时间序列分解方法，其特征在于：步骤二将历史交通流序列分解为趋势、周期和残差3个成分具有以下步骤：

(1)定义所采集的交通流序列为Y(t)，周期长度为C，周期数量为M，则交通流序列Y(t)的样本数量为MC，初始化趋势分量T(t)＝0；

(2)去趋势，用原始交通流序列Y(t)减去趋势分量T(t)，得到去趋势的时间序列

即：

(3)周期子序列平滑，对每个周期子序列作K为M/2的LOESS，并分别向前和向后各延展一个周期；

即表示对于自变量为X＝[1，2，3，4，..，M]^T，因变量为周期子序列Y＝[Y(c)，Y(C+c)，Y(2C+c)，...，Y((M-1)C+c))]^T，基于K＝M/2的LOESS求X＝[0，1，2，3，4，..，M，M+1]^T时对应的估计值

再把所有周期子序列估计值

按照时间先后重组成临时时间序列

上述周期子序列S_c(t)＝{Y(c)，Y(C+c)，Y(2C+c)，...，Y((M-1)C+c))}，c＝1，2，...，C表示由每个周期同一位置的数据所组成的时间序列；

(4)周期子序列的低通量过滤，对临时时间序列

依次做长度为C、C和3的移动平均，然后再做K＝C的LOESS，得到时间序列L(t)，t＝1，2，...，MC；

(5)去平滑周期子序列趋势，用

减去L(t)，即：

计算周期分量，依据F(t)，首先计算单个周期的周期分量Q(t)，t＝1，2，...，C：

然后再把将Q(t)扩展M个周期，得到周期分量P(t)，t＝1，2，...，MC

(6)去周期，用原始交通流序列Y(t)减去周期分量P(t)，得到去周期的时间序列

即：

计算趋势分量，对去周期的时间序列

做K＝C的LOESS，得到趋势分量T(t)；

(8)终止条件判断，若满足终止条件，则执行下一步骤(9)，否则继续执行先前步骤(2)至(7)；

(9)计算残差分量，用原始序列Y(t)减去趋势分量T(t)和周期分量P(t)，得到残差分量R(t)，即：

R(t)＝Y(t)-T(t)-P(t)。

4.根据权利要求1所述的一种交通流时间序列分解方法，其特征在于：步骤三将新观测的交通流序列动态分解为趋势、周期和残差3个成分，其具有以下步骤：

首先定义新观测数据Z(t)，t＝MC+1，MC+2，...

(1)按照时间先后顺序遍历Z(t)中的所有的数据z_t，执行步骤(2)至(5)，直到所有的观测数据都完成了分解；

(2)计算周期分量，根据z_t所在周期的位置和权利要求3中步骤(5)得到的单周期分量Q(t)，计算其周期分量p_t；

(3)去周期，用原始观测值z_t减去周期分量p_t得到去周期数值

(4)计算趋势分量，把

附加

的尾部，并做K＝C的LOESS求该点的估计值，得到趋势分量t_t；

(5)计算残差分量，用原始观测序列z_t减去趋势分量t_t和周期分量p_t，得到残差分量r_t，即：

r_t＝z_t-t_t-p_t。

Images