㊀周晨静２һ㊀荣㊀建１,３㊀李㊀巍４㊀李㊀琪４
(１．北京工业大学北京市城市交通运行保障工程技术研究中心㊀北京１００１２４;
２．北京建筑大学土木与交通工程学院㊀北京１０２６１６;
３．北京工业大学北京市交通工程重点实验室㊀北京１００１２４;４．北京市市政工程设计研究总院㊀北京１０００８２
)摘㊀要：微观交通仿真中的参数标定是科学应用仿真技术的前提.为了减少参数标定的维度，提高仿真手段工程应用效率，将微观仿真参数分为基础参数和模型参数２个部分，并设计道路工程试验方案，对部分基础参数进行工程实测研究，构建实际运动场景下车辆运行加速度取值特征及加速度 速度的关系探究.根据车辆运行过程，将车辆加速过程分为车辆起步㊁提速㊁稳定３个阶段，并依据道路实测数据，应用数据包络的方法确定加速度与速度的函数关系及取值范围.在车辆起步阶段(０~７．２ k m /h )，需要较大动力完成起步动作，速度低加速度高，加速度 速度斜率k ＝０．１８３２；车辆提
速阶段(７．２~３５．４k m /h )，需要确定周边运行环境处于安全状态，加速度有一定程度减少然后逐渐上升，加速度 速度关系整体斜率k ＝０．０２５９；车辆稳定阶段(３５．４k m /h 以上)，速度相对稳定，加速度 速度斜率约等于０.由此确定了基础参数中加速度的阈值范围，对仿真软件中默认加速度 速度关系进行优化.
关键词：交通工程；微观交通仿真；基础参数；参数标定；工程实测；加速度
中图分类号:U ４９１㊀㊀文献标志码:A㊀㊀d o i :１０．３９６３/j
．i s s n ．１６７４Ｇ４８６１．２０１９．０２．０１５A nA n a l y s i s o fA c c e l e r a t i o n s p e e dP a r a m e t e rT h r e s h o l dB a s e d o n A c t u a l E n g i n e e r i n g M
e a s u r e m e n t L I NZ i h e １,
３㊀Z H O UC h e n j i n g ２һ㊀R
O N GJ i a n １,
３㊀L IW e i ４㊀L I Q i ４(１．B e i j i n g E n g i n e e r i n g R e s e a r c hC e n t e r o f U r b a nT r a n s p o r tO p e r a t i o nG u a r a n t e e ,B e i j i n g U n i v e r s i t y o f T e c h n o l o g y ,B e i j i n g １
００１２４,C h i n a ;２．S c h o o l o f C i v i l a n dT r a n s p o r t a t i o nE n g i n e e r i n g ,B e i j i n g U n i v e r s i t y o f C i v i lE n g i n e e r i n g a
n d A r c h i t e c t u r e ,B e i j i n g １０２６１６,C h i n a ;３．B e i j i n g K e y L a b o r a t o r y o f T r a f f i cE n g i n e e r i n g ,B e i j i n g U n i v e r s i t y o f T e c h n o l o g y ,
B e i j i n g １
００１２４,C h i n a ;４．B e i j i n g G e n e r a lM u n i c i p a l E n g i n e e r i n g D e s i g n&R e s e a r c h I n s t i t u t eC o ．,L t d ．,B e i j i n g １
０００８２,C h i n a )A b s t r a c t :P a r a m e t e r c a l i b r a t i o n i n t r a f f i cm i c r o s i m u l a t i o n i sa p r e m i s eo f s c i e n t i f i ca p p
l i c a t i o n f o r s i m u l a t i o nt e c h Ｇn o l o g y ．I n o r d e r t o r e d u c e d i m e n s i o n o f p a r a m e t e r c a l i b r a t i o n a n d i m p r o v e e f f i c i e n c y o f s i m u l a t i o n t o o l s i n e n g i n e e r i n g a p
Ｇp l i c a t i o n ,m i c r o s i m u l a t i o n p a r a m e t e r s a r e d i v i d e d i n t o t w o p a r t s :b a s i c p a r a m e t e r s a n dm o d e l p
a r a m e t e r s ．A n a c t u a l r o a d t e s t i s d e s i g n e d f o r e n g i n e e r i n g m e a s u r e m e n t o nt h e
b a s i
c p a r a m e t e r s ．R e l a t i o n s h i p
sb e t w e e na c c e l e r a t i o nc h a r a c t e r i s t i c s a n d a c c e l e r a t i o n s p e e da r e a n a l y z e d ．A c c o r d i n g t ov e h i c l e o p e r a t i o n ,a na c c e l e r a t i o n p r o c e s s i sd i v i d e d i n t o t h r e e p h a s e s :s t a r t i n g ,i n c r e a s i n g ,a n d s t a b i l i z i n g ．A c c o r d i n g t o a c t u a l d a t a ,am e t h o d o f d a t a e n v e l o p e i s a p p
l i e d t o d e t e r m i n e r e l a t i o n Ｇs h i p b e t w e e na c c e l e r a t i o n a n d s p e e d ,a n d i t s r a n g e o f v a l u e s ．I n s t a r t i n gp h a s e (０－７．２k m /h ),v e h i c l e n e e d s l a r g e p o w e r t o s t a r t r u n n i n g ．I t s s p e e d i s l o wa n da c c e l e r a t i o n i sh i g h ,t h es l o p eo f a c c e l e r a t i o n s p e e d i s k ＝０．１８３２．I n i n c r e a s i n g
p h a s e (７．２－３５．４k m /h ),i n o r d e r t o e n s u r e s a f e o p e r a t i o n ,i t s a c c e l e r a t i o n r e d u c e s a l i t t l e ,t h e n g r a d u a l l y i n c r e a s e s ,t h e o v e r a l l s l o p e o f a c
c e l e r a t i o n s p e e d i s k ＝０．０２５９．I ns t a b i l i z i n gp h a s e (a b o v e３５．４k m /h ),i t s s p e e d i s r e l a t i v e l y s
t a b l e ,a n d t h e s l o p e o f a c c e l e r a t i o n s p e e d i s a p p r o x i m a t e l y z e r o ．I t d e t e r m i n e s a f i n e r e l a t i o n s h i p b e t w e e n a c c e l e r a t i o n a n d s p e e d i nb a s i c p a r a m e t e r s ,a n do p t i m i z e s a r e l a t i o n s h i p b e t w e e nd e f a u l t a c c e l e r a t i o na n d s p e e d i n s i m u l a t i o n s o f t w a r e ．K e y w
o r d s :t r a f f i ce n g i n e e r i n g ;t r a f f i c m i c r os i m u l a t i o n ;b a s i c p a r a m e t e r s ;p a r a m e t e r sc a l i b r a t i o n ;e n g i n e e r i n g m e a s u r e m e n t ;a c c e l e r a t i o n
７
０１基于工程实测的加速度
速度参数阈值分析 林子赫㊀周晨静㊀荣㊀建㊀李㊀巍㊀李㊀琪
收稿日期:２０１８Ｇ１０Ｇ２０
㊀∗北京市科技计划项目(D １７１１００００３９１７００２
)资助㊀㊀第一作者简介：林子赫(１９９４ )，硕士研究生．研究方向：驾驶行为㊁车道边通行能力．E Ｇm a i l :l i n s a n i t i e s t ＠１６３．c o m ㊀һ通信作者：周晨静(１９８７ )，博士，讲师．研究方向：驾驶行为㊁交通安全．E Ｇm a i l :１３４３９０９９７８１＠１３９．c o m
0㊀引㊀言
现实交通场景是一个多维度㊁多因素相互作
用的复杂系统，加上系统在时间与空间上的变化特性多种多样，用传统的数学方法难以描述现实场景运行的特征与过程.微观仿真技术以事件和时间为推动机制，演化交通系统运行过程，其克服了传统数学的缺点，可以描述高维度的复杂交通系统，因而得到大范围的应用，因此，如何保证微观仿真系统的真实性成为该领域的核心研究内容.
随着对交通仿真技术的深入研究，仿真标定的方法逐渐趋于成熟.交通仿真模型参数标定工作可以提高模型的精确度，为仿真过程中再现实
际驾驶员行为与交通特征提供可能[
１]
.国外学者B e n e k o h a l [２]
早在１９９１年就开始对微观交通仿真
模型参数进行校正.H e l i n g a [３]
阐述的仿真校正工作内容有３个部分：前期分析准备；仿真模型初
始化校正；分析结果评价.B r o c k f e l d 等[４]
利用
G P S 采集车辆运行记录，对１０种跟车模型参数进行标定，标定结果显示在一定误差范围内各个
模型均可再现车辆运动过程.B a r t i n 等[５]通过实
际案例研究，介绍了基于多源数据建立大规模交通仿真模型的标定过程.国外仿真技术逐渐成熟，而国内学者对于仿真标定方法的研究起步较
晚.荣建等[６]在２００１年最早利用交通仿真工具
求解道路通行能力，通过对实测数据及模拟数据的速度Ｇ
流量曲线差异控制，对仿真模型进行标定.张智勇[７]
利用G P S 获取跟驰模型数据,
对C A 跟车模型进行标定.王晓原[８]
利用五轮仪测量前后相邻实验车辆运行过程中车位置㊁速度及加速度数据，对测量运行指标进行标定工作.
随着遗传算法等智能优化技术的发展，仿真标定可以开展全部参数的系统性寻优工作，避免
了实验设计方法中以偏概全的缺点.李志明等[
９]应用遗传算法开展了基于V i s s i m 仿真平台的模
型参数标定工作.张长春等[１０]
通过正交实验法
对V i s s i m 仿真模型参数进行了标定.张清华[１
１]通过改进的遗传算法标定V i s s i m 仿真模型中的最
大减速度由原来的－３m /s ２标定为－５．８６m /s
２
.唐少虎和刘小明[１２Ｇ１３]针对人工萤火虫算法的不
足，提出了基于I A G S O 算法的北京市信号交叉口仿真参数自动校正流程和方法，并标定出车辆
的最大减速度为－５．５m /s ２
.G a l l e l l i 等[１
４]指出尽管仿真校准过程非常重要，但由于现场数据收集的困难和现有仿真指南的不足，部分仿真参数通常在在默认值下进行.从整个交通仿真技术的发展进程中看，多数文献采用智能算法和实验设计的方式来标定参数，用物理实测标定参数的文献非常少，但物理实测的方法反而会得到更精准
的参数取值[１５]
.E s s a 等[１６]提出V i s s i m 部分参
数的取值可以在指定区域内进行选择，文章还提
到与使用默认的V i s s i m 参数相比，在现场实测的动态参数与仿真之间提供了更好的相关性.
丝比论坛经过近３０年的发展,
微观仿真模型标定流程逐渐成熟，但对于仿真参数阈值分析有待深入研究.仿真过程中通常选用软件默认参数开展仿真模型标定工作，然而仿真参数标定后的值却存在较大波动.当前的微观仿真参数标定是将所有的参数进行数据优化.若模型参数过多，则存在不止一个最优解，那么最优解的特性分析是模糊的，无法判断参数之间的相互影响.虽然有大量的参数寻优方法，但寻优结果如果没有物理实测数据的约束，会失去自身的物理意义.因此本文根据参数标定特性，提出将微观仿真参数分为基础参数和模型参数２个部分.其中基础参数从工程实践中得到，模型参数由数学方法标定得到.基础参数作用于所有的车辆运动模型，包括车辆最大加速度㊁可接受加速度㊁最大减速度㊁可接受减速度及期望车速等，利用物理实测工程试验方法获取，而非全部采用数据标定；模型参数仅作用于特定驾驶行为，如换道间隙指标仅仅作用于换道行为过程，与其他车辆运动行为无关.本文利用工程实测的方法，确定基础参数的具体取值范围，突破以往微观仿真模型参数标定工作纬度高㊁难具有工程解释能力的局限.
1㊀车辆动力性能与运行特性分析
在进行车辆性能测量测评时，通常以车辆由
０~１００k m /h 极限加速所用的时间(
百公里加速时间)作为主要指标进行分析.通常２．０T 的中
８
０１交通信息与安全㊀２０１９年２期㊀第３７卷㊀总２１７期
型轿车百公里加速时间分布在７~８s(对应加速度为３．４７m/s２和３．９６m/s２),１．６L紧凑型轿车的百公里加速时间分布在１１~１４s(对应加速度为２．５２m/s２和１．９８m/s２)，但在现实中很少有驾驶员会以极限加速度来驾驶.从汽车网站上寻市面上销量较高家用车(售价１０万~２５万元)：伊兰特㊁桑塔纳和帕萨特,３款车型均属于经济型家用轿车.结合官方公布百公里加速时间，用对应车型的实测数据计算出车辆加速度各指标，统计数据见表１.各车之间标准差差值较大，可见加速度差异较大，以单一趋势对加速度进行标定并不准确.
表１㊀车辆加速度统计表
T a b．１㊀S t a t i s t i c s o f t h e v e h i c l e s a c c e l e r a t i o n
m/s２车辆系列最大值最小值平均值标准差伊兰特２０１７款２．７６１．９８２．３７０．３１８桑塔纳２０１７款３．４４２．９９２．３７０．１８３帕萨特２０１７款４．０３２．５１３．２７０．６２１整体均值３．４１２．４９２．９５０．３７４㊀㊀期望车速是指驾驶员在行驶过程中希望达到的最高 安全 行驶车速，车辆行驶不受其他车辆影响.车辆在加速至期望车速前通常分为３个
加速阶段：车辆起步阶段㊁车辆提速阶段和车辆稳定阶段，各个阶段加速度感受程度是不一样的.笔者选取调查车辆的８５％位车速作为期望车速[１７].由于惯性的存在，车辆在行驶过程中的加(减)速状态下会导致车内人员有不同程度的后仰或前倾的感受.乘客在行驶过程中的加㊁减速未有前倾或后仰的感觉，称为适应可接受加(减)速度；若乘客有后仰或前倾的感觉，但不需要借助身体的其他部分完成自身平衡，如手㊁脚等，称为不适应可接受加(减)速度；乘客若有明显的后仰或前倾的感觉，并需要手㊁脚力量保持平衡的加(减)速，称为不适应不可接受加(减)速度.在城市道路中，驾驶员通常采取适应可接受和不适应可接受加(减)速度.车辆在达到期望车速前，驾驶员往往会选取适应可接受加速度或不适应可接受加速度加速运行；车辆达到期望车速后，驾驶员会选取适应可接受加速度，使车辆保持一种期望运行状态.在紧急情况下，为保证行驶安全不与前方车辆碰撞，车辆选取最大减速度减速运行，车辆最大减速度可达１０m/s２.车辆运动的期望加(减)速度是指驾驶员加(减)速行驶过程中选取的加(减)速度的期望数值，该参数与车辆自身动力性能和驾驶员驾驶感受有关.
在当前微观仿真模型标定工作中，对车辆期望车速开展调查的较多，而对车辆运行加速度往往采用软件默认参数，以V i s s i m为例，车辆默认启动加速度为３．５m/s２㊁减速度为－２．８m/s２，加(减)速度的默认值在起步阶段达到了赛车在比赛
中的加速水平，有强烈的推背感和不适感，属于不适应不可接受加速度，而在现实中车辆运动性能表现远远达不到车辆测试性能，默认加速度随速度增加而降低的趋势与实际加速度相差甚远，见图１，
可见仿真软件默认参数与车辆自身实际性能有较大差异，与城市道路严重不符
.
图１㊀系统默认期望加速度与实际道路加速度对比
F i g．１㊀C o m p a r i s o nb e t w e e n s y s t e md e f a u l t
e x p e c t a t i o na c c e l e r a t i o na n da c t u a l r o a da c c e l e r a t i o n 2㊀数据调查与分析
２．１㊀实验调查内容与方案设计
微观仿真模型参数计算与仿真车辆自身速度密切相关.目前微观仿真模型标定工作注重仿真参数的寻优，却缺少车辆基本运动特性实测数据，使得仿真模型标定结果不具有实际物理意义，降低了标定结果的可移植性.笔者以工程实测的方法，开展对基础参数具体取值范围的研究.此次实验采用V i d e oV B O X P r o产品进行测试，其大量应用于运动车辆性能测试[１８Ｇ１９].该设备可以以２０H z的更新率更新所有G P S参数，并且能记录车辆运行环境㊁判断车辆运动变化，及时便捷的获取车辆运动数据.同时，设备搭载视频图像模块可以实时采集车辆运动数字及图像数据.１)调查内容.实验以获取实际驾驶行为特性为目的，主要调查内容包括车辆运动最大加速度㊁最大减速度㊁可接受加速度㊁可接受减速度㊁速度与加速度变化关系等.实验使用道路最常见的１．６L自动挡车辆，并搭载V i d e oV B O X P r o高精度车辆行驶记录仪完成.
２)方案设计.该实验在社会道路上完成，选取由三环快速路㊁南磨房路㊁西大望路㊁建国路辅
９０１
基于工程实测的加速度 速度参数阈值分析 林子赫㊀周晨静㊀荣㊀建㊀李㊀巍㊀李㊀琪
路围成的环状路径作为实验道路，见图２.该路径包括城市快速路㊁主干路和次干路３个层次路网，以便研究不同道路等级的期望车速研究.道路具体几何参数特性见表２
.
图２㊀调查路线设计
F i g ．２㊀T h e r o u t e i n t h e s t u d y a r e a 表２㊀调查道路几何特性表
T a b ．２㊀G e o m e t r y c
h a r a c t e r i s t i c s o f t h e r o a d 道路名称道路等级限制速度
/(k m /h
)长度
/k m 车道数
(单向)
东三环快速路８０２．７４车道南磨房路次干路５０１．３３车道西大望路主干路６０
２．７３车道建国路
呜莎测试
快速路辅路
１．３
３车道
㊀㊀实验选取２０名有熟练驾驶经验的驾驶员,
并记录驾驶员性别㊁年龄㊁驾龄㊁驾驶时状态等特征.每位驾驶员驾车行驶１个循环，由V i d e oV B O X P r o 设备记录实验过程中的车辆运行数据.该实验研究的重点是加(减)速度的大小及车内人员感受的加减速状态，故选择有较多信号灯的路线，其中该路线有１７个信号灯，保证了驾驶员在驾驶过程中有足够加㊁减速的数量.本实验的目的是探究加(减)速度和速度之间的关系，与驾驶员操作特性无关，故驾驶员无需做重复实验.
２．２㊀调查数据分析
应用G P S 设备获取车辆运行数据，通常会存在大量数据噪声，利用V i d e oV B O X P r o 设备配套软件P e r f o r m a n c e T o o l s 对调查数据进行优化.此外,P e r f o r m a n c eT o o l s 可以依据用户设置速度变化区间，对数据进行加减速区间的截取，十分便捷地获取车辆加减速运行数据.２０名驾驶员开展实验数据采集过程中，由于设备及其他随
机因素影响，共收集有效驾驶数据１７份.
２．２．１㊀确定期望车速
在道路实际调查中，车辆运行受较多外界因素干扰，驾驶员减速让行行为较多.当车辆运行速度高于３０k m /h 时达到较为稳定运行状态(道路设计速度５０k m /h ).对实验路段中车辆运行高于３０k m /h 的运行状态进行统计分析，统计得出各路段不同运行速度累计表，如表３所示.由于期望车速受驾驶员自身特性影响很大，故选取８５％位车速作为道路期望车速,
以避免过多的实验数据难以达到期望运行状态[２０
]，因此，选取各个驾驶员的８５％位车速均值作为路段期望车速,
以此作为车辆加(减)速度分析的区间范围.
表３㊀实测路段速度统计表
T a b ．３㊀M e a s u r e d s p
e e d s t a t i s t i c s 速度统计/(k m /h )百分位速度
５０％８５％９０％M a x
南磨房路速度统计４１．９５０．５５１．８５５．０西大望路速度统计４５．３５５．７５７．２６２．９三环路速度统计
６４．３
７１．１
７２．３
７７．４２．２．２㊀加速度数据分析
当车辆运行未达到期望运行状态时，驾驶员将加速行驶以逼近期望车速.实测过程中，车辆加速运行包括自由启动加速及跟车启动加速两种类型.其中自由启动加速是本文研究的重点内容，跟车启动加速过程是车辆相互作用的结果，故不予考虑.研究针对南磨房路㊁大望路及三环路所获取样本的８５％位车速，分别设置加速数据获取区间为０~５０．５,０~５５．７,０~７１．１k m /h ，结合视频信息获取持续车辆持续加速数据.分析时发现,
由于南磨房路为城市次干道，信号交叉口较多，很难获取自由加速启动过程；而当车辆进入三环道路
isia
时，已经是高速运动状态.所以，以西大望路车辆自由启动过程为重点分析对象.大望路部分实验车辆运动运行速度及加速关系散点图见图３.２．３㊀加速度包络区区间分析以往数据分类，通常采用聚类分析方法开展，然而车辆启动加速过程中，加速度的选取受很多因素的影响，不同速度下加速度变化差异性很大.单纯依靠数据进行聚类分析，分类结果并不具有物理意义.依据实际驾驶经验，车辆在加速过程中可以分为挂档起步轻踩油门㊁适当增加油门踩踏力度提升车速㊁适当减少油门踩踏力度保持车速３个阶段.同样，依据数据表现特性可依据车辆自由加速过程，将其划分为３个阶段：起步加速０
１１交通信息与安全㊀２０１９年２期㊀第３７卷㊀总２１７期
图３㊀车辆运行速度与加速度关系样本
F i g ．３㊀S a m p l e o f r e l a t i o n s h i p b
e t w e e n ㊀㊀㊀v e h i c l e s p
e e da n da c c e l e r a t i o n 阶段，该阶段车辆以较小的加速度运行，车辆具有一定速度且在确保行车安全后进入稳定加速阶段，该阶段车辆按照适应可接受加速度运行，迅速提升车速；在接近期望车速时，车辆加速度逐渐减小，车辆进入速度稳定阶段，车速趋于稳定.为此，需要确定各个阶段分界速度阈值与加速度的关系，才可确定整体加速过程速度与加速度变化关系，因此，对以上数据进行汇总，以峰值作为临界点将数据划分为３个阶段，见图４.由车辆运行速度与加速度关系散点图可以看出，速度与加速度并没有确切的数学关系，但在一定的速度条件下，加速度可以在一定区域内随机选取
.
图４㊀车辆运行速度与加速度关系样本汇总F i g ．４㊀S u mo f r e l a t i o n s h i p b
e t w e e nv e h i c l e s p
e e da n da c c e l e r a t i o n 然而实测数据往往包含大量随机因素的影响致使数据分析难度很大，在同一速度条件下也往往对应多个加速度值.为能给出确定性的速度及加速度变化关系，对相同速度条件下加速度值做均值化处理，以减少随机因素的影响，数据处理结果见图５.
由于数据点较为分散，单一方程求解的R ２
德庆公社值均小于０．４，故打算用包络线的方式对数据进行分类划分.在几何学中，某个曲线族的包络线,
是跟该曲线族的每条线都有至少一点相切的一条
图５㊀均值化后车辆运行速度与加速度关系样本汇总
F i g ．５㊀S u mo f r e l a t i o n s h i p b e t w e e nv e h i c l e s p
e e d ㊀㊀㊀a n da c c e l e r a t i o na
f t e r a v e r a
g e p r o c e s s i n g
曲线.因此用处理后的数据以数据包络线的形
江苏广播网式，对加速度取值区间进行描述.根据图５可得阶段一与阶段二的速度临界值为７．２k m /h ，阶段二与阶段三的速度临界值为３５．４k m /h ，从而分别确定了车辆起步㊁提速㊁稳定阶段的加速度与速度关系.利用线性回归分析方法对数据进行拟合，并按照相同拟合斜率求解各阶段数据分布包络线，得到第一阶段加速度Ｇ
速度拟合方程组:y
云芝糖肽
１＝０．１８３２x ＋b １,b ɪ(－０．３９８９,０．１８１１)，具体结果见图６.按照相同数据处理方法，对第二㊁第三阶段数据进行拟合分析，分别得到相应的加速度Ｇ速度拟合方程组:y ２＝０．０２５９x ＋b ２,b ɪ(－０．０８０３,０．９９３１),y ３＝－０．０６４４x ＋b ３,b ɪ(３．１１７２,４．０８０２)，分析结果见图７~８.图６㊀第一阶段数据回归分析
F i g ．６㊀R e g r e s s i o na n a l y
s i s o f t h e f i r s t p h a s
e 图７㊀第二阶段数据回归分析
F i g ．７㊀R e g r e s s i o na n a l y
s i s o f t h e s e c o n d p h a s e 不同道路等级在各阶段分界阈值划分上有所
差别.期望车速越高，第二㊁第三阶段的速度临界阈值越大.通过实验可以确定第二㊁第三阶段速度临界阈值为３５．２~４６．０k m /h 之间.相对于加速度，车辆减速度并没有明显运动阶段的特征表现,
１
１１基于工程实测的加速度
速度参数阈值分析 林子赫㊀周晨静㊀荣㊀建㊀李㊀巍㊀李㊀琪
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议