一种基于时空数据嵌入的交通流量预测方法

文档序号:29216461发布日期:2022-03-12 11:31阅读:489来源:国知局
一种基于时空数据嵌入的交通流量预测方法

1.本技术属于交通流量预测技术领域,具体涉及一种基于时空数据嵌入的交通流量预测方法。


背景技术:

2.交通流量预测是时空预测中的一项典型任务,其目的是根据历史的交通流量来预测未来的交通流量。交通流量预测不仅可以预测潜在的道路拥堵来帮助管理者及时引导交通,还可以帮助出行者提前计划或调整出行路线。因此,实现准确、稳定的交通流量预测势在必行。然而,由于交通流量具有复杂的时空依赖性,获取精确的交通流量预测结果仍然是一个巨大的挑战。
3.目前,深度学习已经被广泛地应用于各个领域,因为它可以结合简单但非线性的模块来生成原始输入数据的多层次深度表示。最近,基于卷积神经网络(cnn)或图卷积神经网络(gcn)的高级深度学习模型已经被成功地应用于交通预测中。具体而言,基于cnn的模型将交通网络视为一幅图像,并在欧几里得空间(如二维矩阵或规则网格)中提取空间特征。基于gcn的模型将交通网络视为一个图结构,其中节点代表传感器(电警设备),权重代表传感器之间的相关性,并通过对交通网络的非欧几里得拓扑结构进行编码来提取空间特征。与基于图像结构的cnn相比,gcn能够更好地利用交通网络固有的拓扑结构和传感器之间的异质相关性,从而获得更好的预测结果。
4.然而,cnn模型中简单的图像无法准确地表示交通网络的真实结构,因为交通网络具有不规则的非欧几里得拓扑结构。因此,传统的cnn无法有效地提取交通网络复杂的空间特征。其次,尽管大多数现有的基于gcn的研究能够获得较好的预测结果,但其只是构建了一个包含固定权重的静态图,无法准确地反映随时间变化而发生改变的传感器之间的相关性。因此,现有的关于深度学习的交通流量预测研究仍然存在一些不足。


技术实现要素:

5.本技术的目的在于提供一种基于时空数据嵌入的交通流量预测方法,实现准确稳定的交通流量预测。
6.为实现上述目的,本技术所采取的技术方案为:
7.一种基于时空数据嵌入的交通流量预测方法,包括:
8.步骤1、获取历史交通流量数据:采集n个电警设备m个时间点的交通流量作为历史交通流量数据;
9.步骤2、基于历史交通流量数据执行时空数据嵌入,包括:
10.步骤2.1、交通流量的区间表示:
11.步骤2.1.1、取历史交通流量数据中交通流量的最大值和最小值作为交通流量范围的上限和下限,并将交通流量范围等距划分为p个区间,区间符号di代表第i个区间,1≤i≤p;
12.步骤2.1.2、确定历史交通流量数据中的每个交通流量的所属区间,并将每个交通流量表示为其所属区间的区间符号以转化为对应的交通流量区间;
13.步骤2.2、交通流量向量的生成:
14.步骤2.2.1、取历史交通流量数据中预设的时间步长内所有的交通流量区间作为输入数据,记为其中t表示时间步长,采用word2vec模型将输入数据转换为嵌入数据其中z表示交通流量向量的嵌入维度,即得到对应的交通流量向量;
15.步骤3、基于时空数据嵌入后得到的交通流量向量提取时间特征得到节点特征矩阵;
16.步骤4、基于时空数据嵌入后得到的交通流量向量提取电警设备之间的相关性得到动态关联图;
17.步骤5、将所述节点特征矩阵和动态关联图输入图卷积神经网络gcn中,得到图卷积神经网络gcn输出的预测结果,所述预测结果为所有电警设备未来t

个时间点的交通流量。
18.以下还提供了若干可选方式,但并不作为对上述总体方案的额外限定,仅仅是进一步的增补或优选,在没有技术或逻辑矛盾的前提下,各可选方式可单独针对上述总体方案进行组合,还可以是多个可选方式之间进行组合。
19.优选地,所述区间di的下界vr
i,l
和上界vr
i,u
计算如下:
[0020][0021][0022]
式中,vr
min
和vr
max
分别表示交通流量范围的下限和上限。
[0023]
优选地,所述基于时空数据嵌入后得到的交通流量向量提取时间特征得到节点特征矩阵,包括:
[0024]
步骤3.1、以不同电警设备对嵌入数据进行划分,划分得到n个嵌入矩阵
[0025]
步骤3.2、取滤波器组,该滤波器组中包含高度为h1,h2,

,hg的g种不同高度,且每种高度具有r个滤波器;
[0026]
步骤3.3、采用不同高度的滤波器对每个嵌入矩阵进行时间特征提取得到多个特征图,将所有特征图连接并展开后得到高层次的时间特征其中lt表示时间特征的长度;
[0027]
步骤3.4、将所有嵌入矩阵对应的高层次的时间特征进行连接得到节点特征矩阵
[0028]
优选地,所述采用不同高度的滤波器对每个嵌入矩阵进行时间特征提取得到多个特征图,包括:
[0029]
步骤3.3.1、选择高度为hg的r个滤波器,g=1,2,

,g,对第n个嵌入矩阵
进行时间特征提取,n=1,2,

,n,取嵌入矩阵wn中第d行到第e行的元素得到的嵌入矩阵wn的子矩阵
[0030]
步骤3.3.2、设置滤波器的宽度为交通流量向量的嵌入维度z,则高度为hg的滤波器可表示为
[0031]
步骤3.3.3、对子矩阵分别应用其一滤波器以获得相应的特征图fm,应用公式如下:
[0032]
fms=f(θ

wn[s:s+h
g-1]+bs),s=1,2,

,t-hg+1
[0033][0034]
式中,fms表示特征图fm的第s个组成部分,f(
·
)表示激活函数,

表示子矩阵和滤波器之间的点积,bs表示偏置项。
[0035]
优选地,所述基于时空数据嵌入后得到的交通流量向量提取电警设备之间的相关性得到动态关联图,包括:
[0036]
步骤4.1、计算电警设备的平均交通流量向量:
[0037][0038]
式中,表示第x个电警设备在时间步长t内的平均交通流量向量;t表示时间步长;表示在时间步长t内第x个电警设备的第t个交通流量向量;
[0039]
步骤4.2、计算第x个电警设备和第y个电警设备的相似度:
[0040][0041]
式中,sim(x,y)为第x个电警设备和第y个电警设备的相似度,y=1,2,

,n,和分别表示平均交通流量向量和的第z个分量;
[0042]
步骤4.3、计算第x个电警设备和第y个电警设备之间的距离dis
x,y

[0043]
步骤4.4、基于第x个电警设备和第y个电警设备的相似度和距离得到动态相关性:
[0044][0045]
式中,a
x,y
表示第x个电警设备和第y个电警设备之间的动态相关性,α和β表示权值系数;σ1表示相似度的标准差;σ2表示距离的标准差;ε表示阈值;
[0046]
步骤4.5、获得动态关联图,所述动态关联图为不同时间段计算得到的加权邻接矩阵a的组合,时间段的时间跨度对应时间步长,所述加权邻接矩阵a为以a
x,y
为元素的n
×
n维矩阵。
[0047]
本技术提供的基于时空数据嵌入的交通流量预测方法,将原始的交通流量数据映
射到向量空间中,并将每个交通流量数据转化为相应的向量表示以量化和度量交通流量之间隐含的相关性。此外,充分利用交通流量向量中包含的关联信息以获取高层次的时间特征,并且为不同的时间段生成不同的动态关联图,从而对电警设备之间的动态相关性进行有效建模,进而实现非欧几里得空间特征的深度提取,以获取准确、稳定的交通流量预测结果。
附图说明
[0048]
图1为本技术的基于时空数据嵌入的交通流量预测方法的流程图;
[0049]
图2为本技术交通流量的区间表示流程图;
[0050]
图3为本技术提取时间特征的流程图;
[0051]
图4为本技术实验中cltfp、gcn、stgcn和stde-dgcn模型在24小时中的预测误差示意图。
具体实施方式
[0052]
下面将结合本技术实施例中的附图,对本技术实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅是本技术一部分实施例,而不是全部的实施例。基于本技术中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本技术保护的范围。
[0053]
除非另有定义,本文所使用的所有的技术和科学术语与属于本技术的技术领域的技术人员通常理解的含义相同。本文中在本技术的说明书中所使用的术语只是为了描述具体的实施例的目的,不是在于限制本技术。
[0054]
其中一个实施例中,针对现有大部分研究都忽略了交通流量之间隐含的相关性,难以有效地提取深层次的时空依赖特征。且大多数研究仅仅构建一个由固定权重组成的静态图来表示传感器之间的相关性,忽略了电警设备之间的相关性会随时间变化而发生改变的事实,本实施例提出一种基于时空数据嵌入的交通流量预测方法。
[0055]
交通流量预测的任务是根据交通网络中观测到的历史交通流量数据来预测未来的交通流量。为了形式化地表示交通流量预测问题,本实施例首先定义了一些关键的概念:
[0056]
定义1:交通网络图。使用加权无向图g=(v,e,a)来表示交通网络,其中v代表一组图节点;e代表一组边,表示传感器之间的连通性;代表一个加权邻接矩阵,表示传感器之间的相关性。
[0057]
定义2:特征矩阵。在图g上观测到的交通流量由特征矩阵表示,其中r代表每个节点包含的特征数量。特别地,被用来表示t时刻的特征矩阵。
[0058]
因此根据上述定义,可以将交通预测问题看作:根据交通网络图g和特征矩阵x来学习映射函数f,并通过公式(1)来预测未来t

个时刻所有传感器的交通流量。
[0059]
[y
t+1
,y
t+2
,

,y
t+t

]=f(g;[x
t-t+1
,x
t-t+2
,

,x
t
])
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0060]
其中,t表示历史时间段的序列长度,即时间步长,表示模型在t+1时刻预测的所有传感器的交通流量。
[0061]
具体的,如图1所示,本实施例的基于时空数据嵌入的交通流量预测方法,包括以
下步骤:
[0062]
步骤1、获取历史交通流量数据:采集n个电警设备m个时间点的交通流量作为历史交通流量数据。
[0063]
步骤2、基于历史交通流量数据执行时空数据嵌入。
[0064]
在时空数据嵌入中,首先以区间的形式表示每个交通流量的值,并提出了一种std2vec方法来生成包含一系列交通流量向量的向量空间。
[0065]
步骤2.1、交通流量的区间表示。
[0066]
由于观测到的交通流量中存在一些出现频率较低的值,这不仅会产生较高的计算成本,而且会降低嵌入的交通流量向量的质量。因此,为了降低嵌入表示的复杂度和获得高质量的交通流量向量,本实施例将数值接近的交通流量分组在一起并用统一的区间符号来表示。如图2所示,具体包括以下步骤:
[0067]
步骤2.1.1、取历史交通流量数据中交通流量的最大值和最小值作为交通流量范围的上限和下限,并将交通流量范围等距划分为p个区间,区间符号di代表第i个区间,1≤i≤p,其中区间di的下界vr
i,l
和上界vr
i,u
计算如下:
[0068][0069][0070]
式中,vr
min
和vr
max
分别表示交通流量范围的下限和上限。
[0071]
步骤2.1.2、确定历史交通流量数据中的每个交通流量的所属区间,并将每个交通流量表示为其所属区间的区间符号以转化为对应的交通流量区间。
[0072]
步骤2.2、交通流量向量的生成。
[0073]
步骤2.2.1、取历史交通流量数据中预设的时间步长内所有的交通流量区间作为输入数据,记为其中t表示时间步长,采用word2vec模型将输入数据转换为嵌入数据其中z表示交通流量向量的嵌入维度,即得到对应的交通流量向量。
[0074]
一个文档包含了一系列单词,而一个传感器包含了一系列不同历史时间点的交通流量区间。根据这种类比关系,本实施例的std2vec方法将每个时间点的交通流量区间视为一个单词,并将每个包含交通流量区间的传感器视为一个文档。然后,std2vec方法采用word2vec模型以生成一个向量空间,该向量空间包含了交通流量区间的关联信息。最后,根据所生成的向量空间,每个交通流量区间被转换成相应的交通流量向量。使用时空数据嵌入中生成的向量空间将输入数据从二维(即传感器和时间步长)转换为三维(即传感器、时间步长和向量维度)。如以每5分钟获取一次,则一天包含288个时间点,即m=288,时间步长可根据实际需求确定,如为12。
[0075]
word2vec模型是将单词转换为单词向量的最有效的方法之一,该方法通过生成一个嵌入空间来使语义相似的单词在该空间中距离相近,然后学习各单词的向量表示。word2vec可以由两种不同的模型实现:cbow模型和skip-gram模型。本实施例采用了skip-gram模型,因为当训练数据足够时,skip-gram模型通常会表现出更好的性能。
[0076]
在获得交通流量向量之后,两个交通流量之间的相关性可以通过它们的向量来计
算。此外,可以通过类比两个词向量之间的高相关性来说明两个交通流量向量之间的高相关性。具体而言,两个高度相关的词向量是指这两个词在文档中经常相邻出现,或者它们周围的词高度相关。类似地,两个交通流量向量之间的高度相关性可以由两种解释来说明。第一种是这两个交通流量通常沿着时间维度相邻出现,这表示它们之间存在短期时间相关性。第二种是这两个交通流量之间的时间跨度较大,但它们周围的交通流量在时间维度上高度相关,从而使这两个交通流量高度相关,这表示它们之间存在长期时间相关性。
[0077]
因此,与原始的交通流量相比,采用本实施例std2vec方法所生成的交通流量向量同时包含了短期相关性和长期时间相关性,有利于深度学习模型进一步提取深层次的时空依赖特征。
[0078]
步骤3、基于时空数据嵌入后得到的交通流量向量提取时间特征得到节点特征矩阵。
[0079]
实施例提取时间特征包括以下步骤:
[0080]
步骤3.1、以不同电警设备对嵌入数据进行划分,划分得到n个嵌入矩阵
[0081]
步骤3.2、取滤波器组,该滤波器组中包含高度为h1,h2,

,hg的g种不同高度,且每种高度具有r个滤波器;
[0082]
步骤3.3、采用不同高度的滤波器对每个嵌入矩阵进行时间特征提取得到多个特征图,将所有特征图连接并展开后得到高层次的时间特征其中lt表示时间特征的长度。时间特征的维度是1*x,x即为这里的lt,为时间特征的长度;
[0083]
步骤3.3.1、选择高度为hg的r个滤波器,g=1,2,

,g,对第n个嵌入矩阵进行时间特征提取,n=1,2,

,n,取嵌入矩阵wn中第d行到第e行的元素得到的嵌入矩阵wn的子矩阵
[0084]
步骤3.3.2、设置滤波器的宽度为交通流量向量的嵌入维度z,则高度为hg的滤波器可表示为
[0085]
步骤3.3.3、对子矩阵分别应用其一滤波器以获得相应的特征图fm,应用公式如下:
[0086]
fms=f(θ

wn[s:s+h
g-1]+bs),s=1,2,

,t-hg+1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0087][0088]
式中,fms表示特征图fm的第s个组成部分,f(
·
)表示激活函数,

表示子矩阵和滤波器之间的点积,bs表示偏置项。
[0089]
步骤3.4、将所有嵌入矩阵对应的高层次的时间特征进行连接得到节点特征矩阵
[0090]
步骤4、基于时空数据嵌入后得到的交通流量向量提取电警设备之间的相关性得到动态关联图。
[0091]
图生成方法对于gcn的特征提取能力至关重要,现有的研究主要使用传感器之间的距离来构建一个具有固定权重的静态图,忽略了传感器之间的相关性会随时间变化而发
生改变的事实。因此本实施例提出了一种新的图生成方法来为不同的时间段生成不同的动态关联图,从而帮助gcn对传感器之间的动态相关性进行有效建模。
[0092]
步骤4.1、计算电警设备的平均交通流量向量:
[0093][0094]
式中,表示第x个电警设备在时间步长t内的平均交通流量向量;t表示时间步长;表示在时间步长t内第x个电警设备的第t个交通流量向量;
[0095]
步骤4.2、计算第x个电警设备和第y个电警设备的相似度:
[0096][0097]
式中,sim(x,y)为第x个电警设备和第y个电警设备的相似度,y=1,2,

,n,和分别表示平均交通流量向量和的第z个分量;
[0098]
步骤4.3、计算第x个电警设备和第y个电警设备之间的距离dis
x,y
。此处距离dis
x,y
根据传感器(电警设备)的地理位置计算,即直接根据电警设备之间的实际地理位置坐标计算得到;
[0099]
步骤4.4、使用基于阈值的高斯核来计算动态关联图对应的加权邻接矩阵a,即根据第x个电警设备和第y个电警设备的相似度和距离得到动态相关性:
[0100][0101]
式中,a
x,y
表示第x个电警设备和第y个电警设备之间的动态相关性,也代表加权邻接矩阵a中的元素,α和β表示权值系数;σ1表示相似度的标准差;σ2表示距离的标准差;ε表示阈值;
[0102]
步骤4.5、获得动态关联图,所述动态关联图为不同时间段计算得到的加权邻接矩阵a的组合,时间段的时间跨度对应时间步长,所述加权邻接矩阵a为以a
x,y
为元素的n
×
n维矩阵。
[0103]
步骤5、将所述节点特征矩阵和动态关联图输入图卷积神经网络gcn中,得到图卷积神经网络gcn输出的预测结果,所述预测结果为所有电警设备未来t

个时间点的交通流量。t

可以根据实际情况确定。
[0104]
本实施例中的gcn采用频谱公式来编码图结构和节点特征,为现有技术,在此不再赘述。归一化的邻接矩阵的定义如公式(9)所示:
[0105][0106]
式中,代表度矩阵。本文中的图卷积的整个运算过程可以使用公式(10)来
表示:
[0107][0108]
式中,f
(q)
、w
1(q)
、w
2(q)
和lm
(q)
分别表示第q层的隐藏特征,两种权重参数和归一化邻接矩阵。代表单位矩阵,u表示第一层的输入特征,relu(nair and hinton 2010)表示激活函数。关于gcn的详细描述可以参考论文kipf和welling(2016)。
[0109]
在另一个实施例中,通过实验进一步验证本实施例提出的流量预测方法的有效性:
[0110]
1、数据描述和预处理
[0111]
本实施例在两组真实的交通数据集pemsd7和pemsd8上验证所提出模型的性能,这两组数据集来自加利福尼亚运输机构性能测量系统(pems)(http://pems.dot.ca.gov/),该系统通过美国加利福尼亚州的39000多个传感器连续实时地收集交通数据。在本数据集中,每五分钟检索一次交通数据,因此每个传感器每天包含288个时间点。
[0112]
数据集pemsd7包含来自洛杉矶县的交通数据。本实施例选择了150个传感器并提取了2017年6月1日至2017年6月30日这一个月的交通流量进行实验,并将20天的数据作为训练集,5天的数据作为验证集,5天的数据作为测试集。
[0113]
数据集pemsd8包含来自圣贝纳迪诺的交通数据。本实施例选择了102个传感器并提取了2017年6月1日至2017年7月31日这两个月的交通流量进行实验,并将41天的数据作为训练集,10天的数据作为验证集,10天的数据作为测试集。
[0114]
2、评估指标
[0115]
为了评估所提出方法的性能,本实施例采用了三种广泛使用的统计指标,包括平均绝对误差(mae),平均方根误差(rmse)和平均绝对百分比误差(mape)。它们的公式定义如下:
[0116][0117][0118][0119]
式中,和y
t
分别表示时间点t的预测值和真实值,n表示测试样本的大小。
[0120]
3、比较对象
[0121]
为了验证本技术所提出的基于时空数据嵌入的交通流量预测方法(简称stde-dgcn模型)的优越性,本实施例选择了以下几种基线模型进行比较:
[0122]
1)ha:历史平均法,利用12个历史时间点的平均交通流量来预测下一个时间点的流量。
[0123]
2)svr(wu等,2004):支持向量回归(svr)是一种广泛使用的机器学习方法。本文将惩罚参数和容忍度设置为0.1和0.03.
[0124]
3)lstm(ma等,2015):长短时记忆网络(lstm)包含一个lstm层,该层具有64个神经元。最终的预测结果通过一个全连接层获得。
[0125]
4)image-cnn(ma等,2017):它用图像(二维矩阵)来表示时空数据,并利用多个卷积层提取时空特征。
[0126]
5)convlstm(shi等,2015):它是一个扩展的全连接lstm,且具有嵌入式卷积层,可以同时提取空间特征和时间特征。
[0127]
6)cltfp(wu和tan,2016):它分别使用cnn和lsym提取空间和时间特征,并将它们融合以获得高级的时空特征。
[0128]
7)gcn(kipf和welling,2016):它利用传感器之间的距离生成gcn的静态交通网络图,并将每个传感器的历史交通流量作为相应的节点特征。
[0129]
8)stgcn(yu等,2018):它利用传感器之间的距离生成gcn的静态交通网络图,并通过具有门控机制的完整卷积结构提取时空特征。
[0130]
4、参数设置
[0131]
在时空数据嵌入中,交通流量区间的数量为200,交通流量区间的下界和上界分别设置为0和1000。word2vec的窗口大小为5,向量维度为100。在时间特征提取中,构建了四种类型的滤波器,其高度分别为3、5、7和9且每种高度的滤波器都有32个。在空间特征提取中,基于阈值的高斯核α、β和ε均被设置为0.5;gcn的层数为4,前三层包含64个神经元,最后一层的神经元个数等于预测范围的长度。滤波器和gcn的激活函数都是relu。采用rmsprop优化器来训练stde-dgcn模型,且采用均方差数作为目标函数。
[0132]
为了确保实验比较的公平性,我们为基线模型和stde-dgcn模型设置了相同的实验参数:将历史时间段的序列长度(时间步长)设置为12,训练代数设置为100,批量大小设置为64,初始学习率设置为0.001,并采用早停法来防止模型过拟合。
[0133]
5、实验结果
[0134]
以数据集pemsd7和pemsd8为历史数据基础,记录stde-dgcn模型和所有基线模型在pemsd7和pemsd8数据集上的预测结果如表1所示。
[0135]
表1两组数据集中不同模型的预测结果
[0136][0137]
由表1可以看出,预测范围分别为15分钟(3个时间步长)、30分钟(6个时间步长)和1小时(12个时间步长)。显然,深度学习模型lstm、image-cnn、convlstm、cltfp、stgcn和stde-dgcn的性能优于具有简单架构的模型,例如ha和svr,这说明深度学习模型更适合处理复杂的时空数据。
[0138]
与convlstm、cltfp、stgcn和stde-dgcn相比,lstm和image-cnn的预测误差更高,因为它们只侧重于提取时间特征或空间特征,这说明交通流量的时间特征和空间特征都会影响交通流量预测的准确性。
[0139]
与convlstm和cltfp相比,基于gcn的模型gcn、stgcn和stde-dgcn性能更好,这说明交通网络的非欧几里得拓扑结构对模型的预测性能是至关重要的。
[0140]
在基于gcn的模型中,stde-dgcn在所有评价指标上都优于gcn和stgcn,因为只有stde-dgcn可以对传感器之间的动态相关性进行有效建模,并全面地利用交通流量之间隐含的相关性。
[0141]
另外,图4展示了表1中四个最佳模型(cltfp、gcn、stgcn和stde-dgcn)在24小时中的预测误差。结果表明,stde-dgcn获得的mae、rmse和mape均低于cltfp、gcn和stgcn,特别是在高峰期(即8:00和17:00)。例如,对于pemsd7,cltfp在17:00的mae约为40,而stde-dgcn在17:00的mae仅约为23。对于pemsd8,stgcn在8:00的rmse约为21,而stde-dgcn在8:00的rmse约为18。此外,stde-dgcn的预测误差在所有时间点的波动都较小,这表明stde-dgcn模型具有更好的准确性和稳定性。
[0142]
以上所述实施例的各技术特征可以进行任意的组合,为使描述简洁,未对上述实施例中的各技术特征所有可能的组合都进行描述,然而,只要这些技术特征的组合不存在矛盾,都应当认为是本说明书记载的范围。
[0143]
以上所述实施例仅表达了本技术的几种实施方式,其描述较为具体和详细,但并
不能因此而理解为对发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本技术构思的前提下,还可以做出若干变形和改进,这些都属于本技术的保护范围。因此,本技术专利的保护范围应以所附权利要求为准。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1