一种残膜回收机防缠绕挑膜装置的制 一种秧草收获机用电力驱动行走机构

一种新的汽车行驶工况构造方法与流程

2021-11-24 23:48:00 来源:中国专利 TAG:


1.本发明涉及一种数据驱动方法,具体涉及一种新的汽车行驶工况构造方法。


背景技术:

2.汽车行驶工况(driving cycle)又称车辆测试循环,是描述汽车行驶的速度

时间曲线。行驶工况作为车辆能耗/排放测试方法和限值标准的基础,体现汽车道路行驶的运动学特征,是汽车行业的一项重要的、共性基础技术,也是汽车各项性能指标标定优化时的主要基准[1]。本世纪初,我国直接采用新欧洲标准行驶循环(nedc)行驶工况对汽车产品能耗/排放的认证[2]。近年来,随着汽车保有量的快速增长,我国道路交通状况发生很大变化,政府、企业和民众日渐发现以nedc工况为基准所优化标定的汽车,实际油耗与法规认证结果偏差越来越大[3]。欧洲在多年的实践中也发现nedc工况的诸多不足,转而采用世界轻型车测试循环(wltc),但该工况怠速时间比和平均速度这两个最主要的工况特征,与我国实际汽车行驶工况的差异更大[4]。另一方面,我国地域辽广,各个城市的发展程度、气候条件及交通状况的不同,使得各个城市的汽车行驶工况特征存在明显的不同[5]。为了更好地理解建立汽车行驶工况曲线的重要性,深入研究作为车辆开发、评价的基础的依据,制定反映我国实际道路行驶状况的测试工况,显得越来越重要。


技术实现要素:

[0003]
针对上述现有技术存在的问题,由于汽车行驶数据的采集设备直接记录的原始采集数据往往会包含一些不良数据,因此,需设计合理的方法对原始不良数据进行预处理,从而完成运动学片段的提取根据上述经预处理后的数据,利用主成分分析和聚类分析法,生成一条能体现参与数据采集汽车行驶特征的汽车行驶工况曲线,并与原始数据进行比较,归纳检验结果。
[0004]
为了实现上述目的,本发明采用的技术方案是:一种新的汽车行驶工况构造方法,包括如下步骤:
[0005]
步骤一、数据预处理:
[0006]
对数据产生影响的情况进行了数据的分析预处理,具体步骤如下:
[0007]
(1)针对gps信号缺失数据的处理,此时将不连续段的数据变成连续段的数据;为判断时间数据的连续性,通过时间序列分析;直接提取连续的时间片段,从而直接对每个连续的时间片段进行数据分析;
[0008]
(2)针对长时间低速行驶的数据处理,将该数据的速度直接赋值为0m/s,此时的数据仍然是连续的时间段;
[0009]
(3)对加、减速度异常情况的处理,对异常的数据直接去除,此时连续段的时间数据变成不连续段的时间数据;根据每一时刻的加速度公式,计算得出该车每一时刻的加速度,如式(1)所示:
[0010][0011]
式中,a
k,k
‑1为第k秒和第k

1秒之间加速度,v
k
为第k秒速度,t
k
为第k秒时刻,n为连续段内总点数;由于认为一辆汽车从0km/h至100km/h的加速时间大于7s,即一般汽车的最大加速度为3.968m/s2,同时最大减速度为7.5~8m/s2,根据这两个判断标准,对不在该加速度范围内的数据进行剔除;
[0012]
(4)由于步骤(3)中剔除了加减速度异常的数据,此时部分连续时间段数据变成不连续数据;需继续对不连续的数据段重复步骤(1)的处理方法,使之变成连续的数据段;
[0013]
(5)针对怠速数据的处理,对每一个连续时间段,将怠速时间超过180s的数据部分进行删除,保留小于180s的怠速部分;因为该数据的处理过程是在连续片段内进行的,故最后仍可以得出连续的时间片段。
[0014]
步骤二、运动学片段提取
[0015]
从步骤一中处理后的有效数据中提取出若干个运动学片段;
[0016]
运动学片段是指汽车从怠速状态开始至下一个怠速状态开始之间的车速区间,将汽车的行驶过程视为多个运动学片段的组合;
[0017]
运动学片段划分主要有两种方法,一种开始于怠速起点终止于下一个怠速起点;另一种开始于怠速终点终止于下一个怠速终点;
[0018]
步骤三、汽车行驶工况建立:
[0019]
(1)特征参数的选取:通过数据处理后得到的运动学片段数据库,确定9个特征参数来代表,运动学片段信息,并计算完成了每个运动学片段特征参数,得到运动学片段特征参数数据库;
[0020]
9个特征参数分别为平均速度v
m
、平均行驶速度v
r
、平均加速度a
w
、平均减速度a
d
、怠速时间比p0、加速时间比p
w
、减速时间比p
d
、速度标准差σ
v
、加速度标准差σ
a

[0021]
(2)主成分分析:运用主成分分析,在不致损失原变量太多信息的条件下尽可能地减小数据变量个数,降低数据维度、简化计算,得出采用3个主成分变量来代替原来的9个特征参数数据,即,通过spss软件得到的协方差矩阵、特征值和原始变量标准化后的数据计算得出3个主成分变量;
[0022]
(3)k

means聚类分析:
[0023]
聚类分析算法的核心是确定中心点,对于给定的样本集,按照样本之间的距离大小,将样本集划分为k个簇;选取的k值决定了样本的聚类效果;让簇内的点的距离尽量小,而让簇间的距离尽量的大,采用的聚类距离为欧氏距离,如式(14)所示:
[0024][0025]
将步骤三中(2)中主成分分析得到的3个主成分,运用sas软件中的procfastclus过程对这三个主成分进行k

means聚类分析;经主成分分析的数据聚成三类行驶工况,分别为1类为拥堵,2类为畅通,3类为特别拥堵;
[0026]
(4)工况曲线拟合
[0027]
距离最终聚点距离越近的运动学片段的汽车行驶工况曲线,越能代表该类别的总体汽车行驶工况曲线;
[0028]
从聚类后的各个类别中提取一定的比例的片段进行拟合,最终将各个类别的片段组合构成一条长度在规定时间内的汽车行驶工况曲线。
[0029]
进一步的,所述对步骤二中的第一种方案进行运动学片段分析,因此需找出怠速时段的怠速起点和运动学片段的最后一个速度为0m/s的点,进行划分运动学片段;
[0030]
怠速时段的怠速起点,即速度、加速度取值同为0;与运动学片段的最后一个速度为0m/s的点为同一点;由于一个完整的运动学片段中,与怠速起点的速度和加速度特征的相同的点不止一个,而满足怠速终点特征的的点只有一个;因此,只需要找到连续片段最后一个速度为0m/s的点;已知最后一个速度为0m/s的点同时满足如下特征如式(2)所示:
[0031][0032]
进一步的,所述步骤三中(2)主成分分析运用主成分分析,在不致损失原变量太多信息的条件下尽可能地减小数据变量个数,降低数据维度、简化计算;每个主成分的贡献率计算方法如式(13)所示:
[0033][0034]
其中,λ1,λ2,

,λ
p
为第1个主成分到第p个主成分所对应的特征值;
[0035]
运用spss软件对所有个运动学片段的9个参数进行主成分分析,计算得到各主成分的贡献率及累计贡献率,分析出前3个主成分累计贡献率,验证主成分来代替这些参数是合理的;因此得出选择3个主成分变量来代替原来的9个特征参数数据,通过spss软件得到的协方差矩阵、特征值和原始变量标准化后的数据计算得出3个主成分变量。
[0036]
进一步的,所述步骤三中(3)借助excel与spss计算得到分类后各类别整体特征参数特征值,得出1类代表了城市道路拥堵情况,第2类运动学片段代表了城市道路通畅情况,第3类则代表了城市道路特别拥堵的情况。
[0037]
本发明的有益效果是:本发明揭示了一种新的汽车行驶工况构造方法,为建立符合我国实际道路行驶工况提供新的思路;在预处理过程中,借助matlab对庞大的数据集进行处理,通过设置合理的方法对数据进行预处理以剔除不良数据,对数据进行运动学片段提取,并计算运动学片段的特征参数;采用主成分分析对特征参数进行降维处理,在此基础上运用聚类分析,将经主成分分析的数据聚成三类,分别为畅通,较为拥堵,拥堵三类行驶工况。分别从这三类行驶工况中选择运动学片段,构造行驶工况曲线,选取计算9个特征参数,尝试了一种新的加速、减速、怠速和匀速工况的划分标准,得到了一种基于城市典型道路的实际速度、加速度等特征参数的行驶工况;使得到的汽车行驶工况曲线,有助于相关部门制定反映实际道路行驶情况的测试工况,进而提高汽车能耗的管理水平。
附图说明
[0038]
图1为数据处理过程图;
[0039]
图2为运动学片段的解释图;
[0040]
图3为快速聚类散点图;
[0041]
图4为汽车行驶工况拟合曲线图;
[0042]
图5为行驶工况拟合曲线和实际曲线对比图。
具体实施方式
[0043]
为使本发明的目的、技术方案和优点更加清楚明了,下面通过附图及实施例,对本发明进行进一步详细说明。但是应该理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限制本发明的范围。
[0044]
根据同一辆车在不同时间段内所采集的40多万条道路行驶记录,旨在建立一条能体现参与数据采集汽车行驶特征的汽车行驶工况曲线。该曲线所体现的汽车运动特征(如平均速度、平均加速度等)能代表所采集数据源的相应特征,两者间的误差越小,说明所构造的汽车行驶工况的代表性越好。
[0045]
由于汽车行驶数据的采集设备直接记录的原始采集数据往往会包含一些不良数据,因此,需设计合理的方法对原始不良数据进行预处理,从而完成运动学片段的提取。根据上述经预处理后的数据,本文提供一种新的汽车行驶工况构造方法,利用主成分分析和聚类分析法,生成一条能体现参与数据采集汽车行驶特征的汽车行驶工况曲线,并与原始数据进行比较,归纳检验结果;具体内容如下:
[0046]
1车辆行驶工况数据预处理
[0047]
1.1数据误差分析方法
[0048]
由于汽车行驶数据的采集设备直接记录的原始数据往往会包含一些不良数据值,会对采集到的数据试验数据造成一定的误差,包括由于高层建筑覆盖等原因导致gps信号丢失,出现汽车加速度数据异常及长时间怠速。误差产生的原因及处理如下:
[0049]
(1)时间不连续
[0050]
本文考虑造成时间不连续是由于高层建筑覆盖或过隧道等,使gps定位不成功,从而造成原始数据中时间不连续。对于这种误差,可以直接提取连续的时间片段,从而直接对每个连续的时间片段进行数据分析。
[0051]
(2)加、减速度异常
[0052]
一般情况下,汽车速度在0km/h至100km/h内的加速时间大于7s,紧急刹车最大减速度为7.5m/s2至8m/s2。本文通过计算每个时间段的加、减速度的大小,可以筛选出异常的加、减速度值,然后直接去除,降低数据采集过程中产生的误差。
[0053]
(3)长时间怠速
[0054]
汽车长时间怠速的原因有两点,第一由于长时期停车(如停车不熄火等候人、停车熄火了但采集设备仍在运行等);第二长时间堵车、断断续续低速行驶等情况造成了长时间的怠速。本文将对这些长时间的怠速最长时间按180s处理,针对以上两点怠速的原因,对怠速时间不超过180s,则保留原始数据,并对超过180s的怠速段数据看作异常值进行剔除,以减小这些数据对后期结果的影响。
[0055]
依据上文对数据产生的误差分析,本文依次对上述可能对数据产生影响的情况进行了数据的分析处理,主要通过5个步骤对原始数据进行处理,具体步骤如下:
[0056]
(6)针对gps信号缺失数据的处理,此时可以将不连续段的数据变成连续段的数据。为判断时间数据的连续性,通过时间序列分析。直接提取连续的时间片段,从而直接对每个连续的时间片段进行数据分析。
[0057]
(7)针对长时间低速行驶的数据处理,将该数据的速度直接赋值为0m/s,此时的数据仍然是连续的时间段。
[0058]
(8)对加、减速度异常情况的处理,对异常的数据直接去除,此时连续段的时间数据变成不连续段的时间数据。根据每一时刻的加速度公式,计算得出该车每一时刻的加速度,如式(1)所示:
[0059][0060]
式中,a
k,k
‑1为第k秒和第k

1秒之间加速度,v
k
为第k秒速度,t
k
为第k秒时刻,n为连续段内总点数。由于认为一辆汽车从0km/h至100km/h的加速时间大于7s,即一般汽车的最大加速度为3.968m/s2,同时最大减速度为7.5~8m/s2,根据这两个判断标准,对不在该加速度范围内的数据进行剔除。
[0061]
(9)由于步骤(3)中剔除了加减速度异常的数据,此时部分连续时间段数据变成不连续数据。需继续对不连续的数据段重复步骤(1)的处理方法,使之变成连续的数据段。
[0062]
(10)针对怠速数据的处理,对每一个连续时间段,将怠速时间超过180s的数据部分进行删除,保留小于180s的怠速部分。因为该数据的处理过程是在连续片段内进行的,故最后仍可以得出连续的时间片段。
[0063]
为了更清晰直观的展现本文数据处理过程,图1选取了某一连续段的原始数据,通过matlab绘制出如下的三维图,表示每次数据预处理步骤(1)至(5)数据段发生的变化。
[0064]
1.2数据预处理结果
[0065]
利用matlab对搜集的不同时段的三个文件的原始数据,按上述步骤进行处理,输出新的excel文件。其中文件1为12月第3周的采集数据,文件2为11月第1周的采集数据,文件3为12月第1周的采集数据。预处理结果如表1所示:文件1经过数据处理之后剩下的数据点记录数为178771个,文件2经过数据处理之后剩下的数据点记录数为139106个,文件3经过数据处理之后剩下的数据点记录数为158198个,共476075条记录。构建汽车行驶工况曲线需要多个运动学片段,接下来将在上述处理后的有效数据中提取出若干个运动学片段。
[0066]
表1剩余数据点记录数
[0067][0068]
2运动学片段提取
[0069]
2.1提取原理
[0070]
从起始点到目的地,汽车受道路类型和交通状况的影响,会多次经过启动,低速,停止等一系列过程。运动学片段是指汽车从怠速状态开始至下一个怠速状态开始之间的车速区间。如图2是定义运动学片段的解释图,整个运动学片段包含怠速、加速、减速、匀速等情况。本文将汽车的行驶过程视为多个运动学片段的组合。
[0071]
构造汽车行驶工况的模型,经典构造策略将整个过程视为连续,利用统计学对整个过程分析,并按道路区域类型的不同进行工况划分,组建若干个车辆行驶候选工况,最后选取各级候选工况合称为总工况。虽然这种方法满足了构造的行驶工况要求,但过程处理
量过大,且两点与两点间可能隐藏相关性,造成数据信息重复。这不仅加重了计算的负担,而且在构造工况时考虑了重叠无效数据信息,对数据信息的汇总会有诱导及偏差。因此,在数据预处理的基础上,本文采用更为优越的运动学片段方法对数据进行分析。
[0072]
2.2提取方法
[0073]
运动学片段划分主要有两种方法,一种开始于怠速起点终止于下一个怠速起点;另一种开始于怠速终点终止于下一个怠速终点。本文提供第一种方案进行运动学片段分析,因此需找出怠速时段的怠速起点和运动学片段的最后一个速度为0m/s的点,进行划分运动学片段。
[0074]
数据预处理中已经找到各整片运动学片段,且将异常的数据进行了剔除和转化,但不是独立的运动学片段。表2归纳出一个运动学片段中不同运动状态(不同的运动状态见图2)情况下的速度和加速度特点。
[0075]
表2各运动状态速度加速度特点
[0076][0077]
通过分析发现,怠速时段的怠速起点(速度、加速度取值同为0)与运动学片段的最后一个速度为0m/s的点为同一点。由于一个完整的运动学片段中,与怠速起点的速度和加速度特征的相同的点不止一个,而满足怠速终点特征的的点只有一个。因此我们只需要找到连续片段最后一个速度为0m/s的点。已知最后一个速度为0m/s的点同时满足如下特征如式(2)所示:
[0078][0079]
2.3提取结果
[0080]
结合上述怠速点的特征,用matlab循环语句实现对运动学片段进行提取。通过matlab软件遍历循环提取后,我们最终得到文件1的运动学片段数为840个,文件2运动学片段数为502个,文件3的运动学片段为560个,共计得到1902个运动学片段。
[0081]
3汽车行驶工况建立
[0082]
3.1特征参数的选取
[0083]
本文需要通过对某些特征参数对所构造的汽车行驶工况进行合理性分析。合理的汽车运动特征评估体系需要考虑多个因素,其中速度和加速度最为重要,而只用二者无法全面的描述所划分的运动学片段,因此仍需引入其他的特征参数来对运动学片段进行描述,过多的参数会消耗大量计算时间,增加计算难度,过少的参数则不能准确地代表大量数据
[7]
。本文选取了共9个对汽车行驶工况较为重要的特征参数,分别为平均速度v
m
、平均行驶速度v
r
、平均加速度a
w
、平均减速度a
d
、怠速时间比p0、加速时间比p
w
、减速时间比p
d
、速度标准差σ
v
、加速度标准差σ
a

[0084]
特征参数的计算公式如下所示:
[0085]
(1)平均速度v
m
[0086][0087]
式中,t
total
为每个运动学片段总时长。
[0088]
(2)平均行驶速度v
r
[0089][0090]
式中,v
k
为v
k
>0的速度,为汽车非怠速时长。
[0091]
(3)平均加速度a
w
[0092][0093]
式中,a
k
为汽车行驶加速度大于0的值,为汽车的加速总时长。
[0094]
(4)平均减速度a
d
[0095][0096]
式中,a
k

为汽车行驶加速度小于0的值,为汽车的减速总时长。
[0097]
(5)怠速时间比p0[0098][0099]
式中,t0为怠速总时长。
[0100]
(6)加速时间比p
w
[0101][0102]
(7)减速时间比p
d
[0103][0104]
(8)速度标准差σ
v
[0105][0106]
(9)加速度标准差σ
a
[0107][0108]
式中,a
m
为平均加速度,其计算公式如下所示:
[0109][0110]
3.2主成分分析
[0111]
通过数据处理后得到的运动学片段数据库,共确定9个特征参数来代表运动学片段信息,并计算完成了每个运动学片段特征参数,得到运动学片段特征参数数据库。由于运动学片段特征参数数据量较大,若直接对9个特征参数进行分析产生的计算量较大,所以需要通过主成分分析进行数据降维处理
[8]
。运用主成分分析,在不致损失原变量太多信息的条件下尽可能地减小数据变量个数,降低数据维度、简化计算
[9

10]

[0112]
每个主成分的贡献率计算方法如式(13)所示:
[0113][0114]
其中,λ1,λ2,

,λ
p
为第1个主成分到第p个主成分所对应的特征值。
[0115]
本文运用spss软件对1902个运动学片段的9个参数进行主成分分析,可以计算得到各主成分的贡献率及累计贡献率等,从表3可以得出,前3个主成分累计贡献率83.690%,说明主成分来代替这些参数是合理的。
[0116]
表3各主成分贡献率及累计贡献率
[0117][0118]
因此,可以选择3个主成分变量来代替原来的9个特征参数数据,通过spss软件得到的协方差矩阵、特征值和原始变量标准化后的数据可以计算得出3个主成分变量,以便进行下一步聚类分析。
[0119]
3.3 k

means聚类分析
[0120]
聚类分析算法的核心是确定中心点,对于给定的样本集,按照样本之间的距离大小,将样本集划分为k个簇。选取的k值决定了样本的聚类效果。让簇内的点的距离尽量小,而让簇间的距离尽量的大
[11]
。本文采用的聚类距离为欧氏距离,如式(14)所示:
[0121][0122]
将3.2中主成分分析得到的3个主成分,运用sas软件中的procfastclus过程对这三个主成分进行k

means聚类分析。多次聚类分析发现分为3类的效果较好。绘制散点图3如下所示:
[0123]
其中,第一类片段包括805个片段,第二类包括1088个片段,第三类包括9个片段。表4截取了聚类结果的部分运动学片段的聚类结果,运动学片段所有距离结果仅展示部分
下表所示:
[0124]
表4运动学片段聚类结果
[0125][0126][0127]
借助excel与spss计算得到分类后各类别整体特征参数特征值,如下表5所示:
[0128]
表5三类运动学片段的特征参数
[0129][0130]
由表5可知,第2类运动学片段的平均速度为22.64km/h,远高于第1类和第3类,而怠速比例为0.19,比其他两类都要低。可以看出第2类运动学片段代表了城市道路通畅情况。第1类平均加速度和平均减速度远低于第3类。我们可以得到第1类代表了城市道路拥堵情况,而第3类则代表了城市道路特别拥堵的情况。
[0131]
3.4工况曲线拟合
[0132]
在上述聚类分析的基础上,我们可以通过每个运动学片段与最终聚点的距离来考察该运动学片段是否能代表该类别的运动学片段
[12]
。距离最终聚点距离越近的运动学片段,该片段的汽车行驶工况曲线,越能代表该类别的总体汽车行驶工况曲线
[13]

[0133]
从聚类后的各个类别中提取一定的比例的片段进行拟合,最终将各个类别的片段组合构成一条长度在1200

1300s的汽车行驶工况曲线。由于第3类的运动学片段数目较少,为了方便计算,我们仅选取第3类别中的前5个距离最终聚点最近的运动学片段。第1类和第2类的运动学片段,我们各选取了前20个离最终聚点距离最近的运动学片段。通过插值法,挑选出第3类的1运动学片段,第1类的5个运动学片段,第2类的7个运动学片段,共计13个运动学片段。
[0134]
运用matlab软件将以上运动学片段组合,构造出如图4所示长度为1300s的汽车行驶工况拟合曲线,图中[0s,9s]区间的运动学片段为第3类运动学片段,[10s,458s]区间的运动学片段为第1类运动学片段,[459s,1300s]区间的运动学片段为第2类运动学片段。
[0135]
3.5合理性检验
[0136]
3.5.1特征参数误差分析检验
[0137]
根据要求,构造的工况曲线要能代表实际所采集数据源(处理后的数据)的相应特征。构造的工况曲线与实际采集数据曲线误差越小,说明构造的汽车行驶工况代表性越好。为了更直观的显示图4中构造的工况曲线和实际采集数据曲线的误差,通过matlab中的unidrnd函数在1902个运动学片段中随机抽取13个运动学片段绘制出一条随机工况曲线。对比图如图5所示,行驶工况曲线与实际采集数据曲线之间的误差较小,拟合效果较好。
[0138]
特征参数误差分析的基本思想是将行驶工况拟合曲线和随机实际采集数据曲线对比,从3.1所定义的9个特征参数作为分析的指标,分别计算所构造道路行驶工况与该道路行驶工况所对应的整体的相对误差,根据所计算的结果来分析构造出的工况是否有效。在研究中,一般相对误差小于10%都在可以接受的范围内
[14]
。相对误差计算公式如式(15)所示:
[0139][0140]
δ
i
为构造出的道路行驶工况的特征参数,δ
j
为该类别对应总体的特征参数。计算结果如表6所示:
[0141]
表6行驶状况误差分析
[0142][0143]
从表6中,我们可以得出,每类别的平均速度、平均行驶速度、怠速时间比相对误差都在7.04%以内,所以该部分数据拟合的效果较为合理。第2类别的特征参数误差均在7.50%以内。第3类别中的加速度时间比为61.63%,误差较大,但考虑到相对误差较大数据的仅有极个别情况,在可接受范围内,因此所构造的汽车行驶工况较为合理。
[0144]
3.5.2相关系数检验
[0145]
相关系数是一种描述随机变量x与y之间相关关系的变量
[15]
,将每一类的各片段参数特征值序列看成是随机变量x的分布,把各类综合的参数特征值序列看成是随变量y的分布。pearson相关系数计算公式如式(16)所示:
[0146][0147]
通过spss双变量的相关性计算出x与y的相关系数,该数越接近于1,说明该片段的参数特征值与整个片段的参数特征值线性相关性越大,该片段的数据就更有代表性
[16]
。表7列出了三类中选取片段相对应相关系数对比,即每一类的道路行驶工况代表特征参数与
道路行驶工况整体参数的相关系数检验。
[0148]
表7相关系数检验结果
[0149][0150]
表7中我们可以得到三类相关系数均接近于1,说明该片段的参数特征值与整个片段的参数特征值线性相关性较大。因此,我们挑选出来的片段的数据就更有代表性,检验结果合理。
[0151]
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换或改进等,均应包含在本发明的保护范围之内。
再多了解一些

本文用于企业家、创业者技术爱好者查询,结果仅供参考。

发表评论 共有条评论
用户名: 密码:
验证码: 匿名发表

相关文献