基于模糊分类算法的盾构掘进地质特征识别方法和系统

文档序号:30072309发布日期:2022-05-18 02:22阅读:249来源:国知局
基于模糊分类算法的盾构掘进地质特征识别方法和系统

1.本发明涉及盾构隧道技术领域,更具体地,涉及一种基于模糊分类算法的盾构掘进地质特征识别方法和系统。


背景技术:

2.随着城市化进程的不断加快,城市之间的交流日益频繁。城际轨道交通系统得到了快速的发展。盾构法施工由于其安全、快速、环保等特点,广泛地应用在城市地下隧道的修建中。盾构法采用盾构刀盘切削岩土,经土仓、螺旋机、传送带等设备将土体传送到渣土车,再经渣土车运出至隧道外部;开挖后的隧道由管片拼装机拼装隧道管片并注浆形成圆形的盾构隧道。盾构机在开挖过程中,地质特征是影响盾构掘进速度最重要的因素。在隧道施工前,通常会进行地质勘察绘制施工区间的地质剖面图。但是这种经勘察孔的地质特征直线连接的方式绘制的地质剖面图不能够准确揭示盾构掘进每一环的地质情况。为了盾构的安全、高效地掘进,盾构掘进过程中的地质特征识别至关重要。
3.盾构掘进地层中岩层的比例对盾构刀具磨损、掘进效率和盾构参数设置具有重大影响,根据不同岩层比例设置盾构参数和盾构掘进时间安排,能够保证盾构安全、高效地运行。田超于2019年在《价值工程》上发表的《复杂条件下长距离硬岩掘进施工关键技术研究》中强调了盾构掘进过程中针对复杂地层进行参数优化的重要性。目前,施工区间的地质特征情况是经过地质勘察孔推测而来,且不能给出盾构掘进所在地层的地质特征,无法为盾构在软硬不均地层中掘进提供准确的地质信息。因此,为了获取盾构掘进刀盘前方的地质情况,指导盾构机的掘进安排及刀具更换,有必要提出一种基于模糊分类算法的盾构掘进地质特征识别方法。
4.经过对现有技术文献检索发现,现有技术公开一种基于深度学习的地质特征检测识别方法及系统,该方案利用深度学习的方法对地质特征信息进行分类,并生成地质特征识别信息,可对地质特征进行细化分析,为地质情况分析提供基础。但是,该方法是针对地质构造数据、地层数据、油气资源数据进行特征识别,为油气探测提供数据信息。此外,该方法通过人为的方式确定地质信息类别标记,人为因素对地质特征的分类结果影响较大,不适用于盾构掘进施工。
5.本发明提供了一种基于模糊分类算法的盾构掘进地质特征识别方法,根据盾构法施工过程中盾构机对地质特征的反应及盾构参数的变化情况,对地质特征进行分类标记,精确获取盾构掘进过程中穿越断面的地质特征情况,能够合理地安排掘进时间和对盾构掘进过程进行提前预警。


技术实现要素:

6.本发明的首要目的是提供一种基于模糊分类算法的盾构掘进地质特征识别方法,根据盾构法施工过程中盾构机对地质特征的反应及盾构参数的变化情况,对地质特征进行分类标记,精确获取盾构掘进过程中穿越断面的地质特征情况,能够合理地安排掘进时间
和对盾构掘进过程进行提前预警。
7.本发明的进一步目的是是提供一种基于模糊分类算法的盾构掘进地质特征识别系统。
8.为解决上述技术问题,本发明的技术方案如下:
9.一种基于模糊分类算法的盾构掘进地质特征识别方法,包括以下步骤:
10.s1:收集隧道施工区间地质资料、地质特征识别参数和出土情况,所述地质资料包括地质剖面图,根据出土情况得到初步的地质特征;
11.s2:利用主成分分析方法对所述地质特征识别参数进行特征提取;
12.s3:将特征提取后得到的数据集采用模糊划分系数fpc、轮廓系数s结合地质剖面图确定地质特征最终类别数k;
13.s4:根据确定的地质特征最终类别数k,结合地质剖面图,更新地质特征;
14.s5:利用模糊分类算法,确定每一盾构掘进环所属的地质特征类别,并与盾构掘进记录的地质特征及地质剖面图对比,确定地质特征识别的准确性,若准确性达到阈值以上,则该识别模型可用于地质特征识别;否则,调整模糊分类算法的参数,直到满足要求为止;
15.s6:将新收集的盾构实时参数经步骤s2处理后,输入到调整后的地质特征识别模型中,确定盾构掘进所在的地质特征类别并输出所属地质特征类别的隶属度。
16.优选地,步骤s1中所述地质剖面图由地质勘察单位在隧道施工前对隧道所在区间进行取土试验后绘制得到,所述地质特征识别参数包括盾构机传感器按照盾构施工环号收集反馈的盾构机参数及其二次变换参数,所述出土情况指的是盾构掘进后由螺旋机排出的渣土情况,包括岩石的比例情况。
17.优选地,所述盾构机参数包括盾构机推力f、推进速度ar、刀盘扭矩t、刀盘转速crs、上部土仓压力uep和下部土仓压力lep,所述二次变换参数包括贯入度pr、场切深指数fpi、扭矩切深指数tpi和比能se,由盾构机参数二次变换得到,并由以下公式确定:
18.pr=ar/crs
19.fpi=f/pr
20.tpi=t/pr
[0021][0022]
式中,ar为推进速度,mm/s;crs为刀盘转速,r/s;pr为贯入度,mm/r;fpi为场切深指数,kn
·
r/mm;tpi为扭矩切深指数,kn
·
r;f为盾构推力,kn;t为刀盘扭矩,kn
·
mm;a为刀盘面积,m2;se为比能,mj/m3。
[0023]
优选地,根据岩石的比例情况,将对应环号的地层记录为:
[0024]
当岩石比例小于20%时,记录为软土地层;
[0025]
当岩石比例在20%至80%时,记录为软硬不均地层;
[0026]
当岩石比例大于80%时,记录为硬岩石层。
[0027]
优选地,所述步骤s2中利用主成分分析对所述地质特征识别参数进行特征提取,具体为:
[0028]
利用主成分分析方法对n
×
m条地质识别参数进行处理,提取前k个主成分,组成新的地质特征识别数据集y,其中,m为所述地质特征识别参数的个数,n为所述地质参数按照
环号记录的条数。
[0029]
优选地,所述主成分分析具体包括以下步骤:
[0030]
s2.1:构建n
×
m维的识别参数数据库,数据库矩阵x由以下公式表示:
[0031][0032]
式中,cm为第m个地质特征识别参数的名称;xn为第m个识别参数中的第n行数据;x
ij
为第i个地质特征识别参数第j行数据的值;
[0033]
s2.2:对数据库矩阵x进行标准化处理,标准化后的第i行,第j列数据x
ij
由以下公式表示:
[0034][0035]
式中,为第m列数据的平均值,sm为第m列数据的标准差;
[0036]
s2.3:计算标准化处理后的数据库矩阵的协方差,构建协方差矩阵cov,第j个因子和第k个因子之间的协方差cov
jk
由以下公式确定:
[0037][0038]
式中,x
ik
为标准化处理后的数据库矩阵中第i行第k列数据,为第j个因子和第k个因子的标准化数据的平均值,根据因子之间的协方差,协方差矩阵由以下公式确定:
[0039][0040]
s2.4:计算协方差矩阵的特征值和相应的特征向量,并根据特征向量计算每个因子所占的权重;
[0041]
s2.5:将特征向量按照对应特征值的大小从左到右排列成一列矩阵,取前k列组成一个降维矩阵p;
[0042]
s2.6:由以下公式将原始数据转变为k维数据:
[0043]
y=xp
[0044]
式中,y为地质特征识别数据集。
[0045]
优选地,步骤s3中所述模糊分类算法,具体为:
[0046]
s301:输入分类类别数c、权重指数m和初始隶属度函数u0;
[0047]
s302:计算聚类中心(vi),聚类中心由以下公式确定:
[0048][0049]
式中,u
ik
为隶属度函数,0《u
ik
≤1;n为样本数量;c为分类类别数,1≤k≤n;1≤i≤c;m为权重指数;vi为第i类聚类中心;xk为第i类中的第k个样本;
[0050]
s303:根据下式更新隶属度矩阵:
[0051][0052]
式中,d
ik
=||x
k-vi||为第i类中第k个样本与第i类的聚类中心的欧氏距离;
[0053]
s304:给定一个小值ε》0,判断是否满足如果满足条件,则结束运算过程输出聚类中心、每一类的样本及其隶属度;否则返回s302。
[0054]
优选地,所述步骤s3中的模糊划分系数是判断模糊分类模糊性程度的标准,分类越明确时,模糊划分系数越大,模糊划分系数由以下公式确定:
[0055][0056]
优选地,所述步骤s3中的轮廓系数s是根据划分类别内样本点的最小组内距离和最大组间距离,是确定模糊分类算法中最终类别数k的一种方法,由以下公式确定:
[0057][0058]bi
=min{b
i1
,b
i2
,b
i3
,

,b
ik
}
[0059]
式中,si为单个样本点的轮廓系数;ai为同一类别内第i个样本点到其他样本点的平均距离;b
ik
为第i个样本点到其他类别内所有样本点的平均距离,轮廓系数(s)为所有样本点轮廓系数的平均值。
[0060]
优选地,所述的最终类别数k指的是当模糊划分系数最大,且轮廓系数较大时的类别数k;当模糊划分系数与轮廓系数不同时出现最大值时,可将地质剖面图与两种方法结合使用,确定最终k值。
[0061]
优选地,所述的地质特征名称是根据地质剖面图中,盾构掘进断面内地层的软硬特性命名,地质特征名称可命名为:软土地层、软硬不均地层、硬岩地层、软土为主地层、硬岩为主地层、孤石地层等。
[0062]
一种基于模糊分类算法的盾构掘进地质特征识别系统,包括:
[0063]
收集模块,所述收集模块收集隧道施工区间地质资料、地质特征识别参数和出土情况,所述地质资料包括地质剖面图,根据出土情况得到初步的地质特征;
[0064]
特征提取模块,所述特征提取模块利用主成分分析方法对所述地质特征识别参数进行特征提取;
[0065]
类别数确定模块,所述类别数确定模块将特征提取后得到的数据集采用模糊划分系数fpc、轮廓系数s结合地质剖面图确定地质特征最终类别数k;
[0066]
更新模块,所述更新模块根据确定的地质特征最终类别数k,结合地质剖面图,更新地质特征;
[0067]
训练模块,所述训练模块利用模糊分类算法,确定每一盾构掘进环所属的地质特征类别,并与盾构掘进记录的地质特征及地质剖面图对比,确定地质特征识别的准确性,若准确性达到阈值以上,则该识别模型可用于地质特征识别;否则,调整模糊分类算法的参数,直到满足要求为止;
[0068]
识别模块,所述识别模块将新收集的盾构实时参数经步骤s2处理后,输入到调整后的地质特征识别模型中,确定盾构掘进所在的地质特征类别并输出所属地质特征类别的隶属度。
[0069]
与现有技术相比,本发明技术方案的有益效果是:
[0070]
本发明利用盾构实时掘进参数,利用主成分分析降低数据的维度,采用模糊划分系数、轮廓系数确定最终的地质特征类别,最后采用模糊分类算法实现盾构掘进过程中地质特征的实时确定,现场工程师可以根据实时确定的地质特征,设置对应的盾构掘进参数、安排刀具更换等任务;此外,该方法能够准确识别各种地质状态,若遇到异常地质特征时,根据识别结果能够对施工过程进行预警,保证盾构施工的安全。同时,本方法通过调整地质特征类别及模糊分类算法的参数,能够实现地质特征识别准确率达90%以上,满足盾构隧道的施工要求。本方法操作简单易行,成本低,能够显著地提高盾构的施工效率,保证了盾构掘进的安全性。
附图说明
[0071]
图1为本发明的方法流程示意图。
[0072]
图2为实施例提供的模糊划分系数与轮廓系数变化示意图。
[0073]
图3为实施例提供的地质特征识别结果2维特征空间图。
[0074]
图4为实施例提供的地质特征识别结果与盾构施工效率对比图。
[0075]
图5为本发明的系统模块示意图。
具体实施方式
[0076]
附图仅用于示例性说明,不能理解为对本专利的限制;
[0077]
为了更好说明本实施例,附图某些部件会有省略、放大或缩小,并不代表实际产品的尺寸;
[0078]
对于本领域技术人员来说,附图中某些公知结构及其说明可能省略是可以理解的。
[0079]
下面结合附图和实施例对本发明的技术方案做进一步的说明。
[0080]
实施例1
[0081]
本实施例提供一种基于模糊分类算法的盾构掘进地质特征识别方法,如图1所示,包括以下步骤:
[0082]
s1:收集隧道施工区间地质资料、地质特征识别参数和出土情况,所述地质资料包括地质剖面图,根据出土情况得到初步的地质特征;
[0083]
s2:利用主成分分析方法对所述地质特征识别参数进行特征提取;
[0084]
s3:将特征提取后得到的数据集采用模糊划分系数fpc、轮廓系数s结合地质剖面图确定地质特征最终类别数k;
[0085]
s4:根据确定的地质特征最终类别数k,结合地质剖面图,更新地质特征;
[0086]
s5:利用模糊分类算法,确定每一盾构掘进环所属的地质特征类别,并与盾构掘进记录的地质特征及地质剖面图对比,确定地质特征识别的准确性,若准确性达到阈值以上,则该识别模型可用于地质特征识别;否则,调整模糊分类算法的参数,直到满足要求为止;
[0087]
s6:将新收集的盾构实时参数经步骤s2处理后,输入到调整后的地质特征识别模型中,确定盾构掘进所在的地质特征类别并输出所属地质特征类别的隶属度。
[0088]
步骤s1中所述地质剖面图由地质勘察单位在隧道施工前对隧道所在区间进行取土试验后绘制得到,所述地质特征识别参数包括盾构机传感器按照盾构施工环号收集反馈的盾构机参数及其二次变换参数,所述出土情况指的是盾构掘进后由螺旋机排出的渣土情况,包括岩石的比例情况。
[0089]
所述盾构机参数包括盾构机推力f、推进速度ar、刀盘扭矩t、刀盘转速crs、上部土仓压力uep和下部土仓压力lep,所述二次变换参数包括贯入度pr、场切深指数fpi、扭矩切深指数tpi和比能se,由盾构机参数二次变换得到,并由以下公式确定:
[0090]
pr=ar/crs
[0091]
fpi=f/pr
[0092]
tpi=t/pr
[0093][0094]
式中,ar为推进速度,mm/s;crs为刀盘转速,r/s;pr为贯入度,mm/r;fpi为场切深指数,kn
·
r/mm;tpi为扭矩切深指数,kn
·
r;f为盾构推力,kn;t为刀盘扭矩,kn
·
mm;a为刀盘面积,m2;se为比能,mj/m3。
[0095]
根据岩石的比例情况,将对应环号的地层记录为:
[0096]
当岩石比例小于20%时,记录为软土地层;
[0097]
当岩石比例在20%至80%时,记录为软硬不均地层;
[0098]
当岩石比例大于80%时,记录为硬岩石层。
[0099]
所述步骤s2中利用主成分分析对所述地质特征识别参数进行特征提取,具体为:
[0100]
利用主成分分析方法对n
×
m条地质识别参数进行处理,提取前k个主成分,组成新的地质特征识别数据集y,其中,m为所述地质特征识别参数的个数,n为所述地质参数按照环号记录的条数。
[0101]
所述主成分分析具体包括以下步骤:
[0102]
s2.1:构建n
×
m维的识别参数数据库,数据库矩阵x由以下公式表示:
[0103][0104]
式中,cm为第m个地质特征识别参数的名称;xn为第m个识别参数中的第n行数据;x
ij
为第i个地质特征识别参数第j行数据的值;
[0105]
s2.2:对数据库矩阵x进行标准化处理,标准化后的第i行,第j列数据x
ij
由以下公式表示:
[0106][0107]
式中,为第m列数据的平均值,sm为第m列数据的标准差;
[0108]
s2.3:计算标准化处理后的数据库矩阵的协方差,构建协方差矩阵cov,第j个因子和第k个因子之间的协方差cov
jk
由以下公式确定:
[0109][0110]
式中,x
ik
为标准化处理后的数据库矩阵中第i行第k列数据,为第j个因子和第k个因子的标准化数据的平均值,根据因子之间的协方差,协方差矩阵由以下公式确定:
[0111][0112]
s2.4:计算协方差矩阵的特征值和相应的特征向量,并根据特征向量计算每个因子所占的权重;
[0113]
s2.5:将特征向量按照对应特征值的大小从左到右排列成一列矩阵,取前k列组成一个降维矩阵p;
[0114]
s2.6:由以下公式将原始数据转变为k维数据:
[0115]
y=xp
[0116]
式中,y为地质特征识别数据集。
[0117]
步骤s3中所述模糊分类算法,具体为:
[0118]
s301:输入分类类别数c、权重指数m和初始隶属度函数u0,初始隶属度函数u0由输出的识别数据集y和步骤s303中的隶属度函数表达式确定;
[0119]
s302:计算聚类中心(vi),聚类中心由以下公式确定:
[0120][0121]
式中,u
ik
为隶属度函数,0《u
ik
≤1;n为样本数量;c为分类类别数,1≤k≤n;1≤i≤c;m为权重指数;vi为第i类聚类中心;xk为第i类中的第k个样本;
[0122]
s303:根据下式更新隶属度矩阵:
[0123][0124]
式中,d
ik
=||x
k-vi||为第i类中第k个样本与第i类的聚类中心的欧氏距离;
[0125]
s304:给定一个小值ε》0,判断是否满足如果满足条件,则结束运算过程输出聚类中心、每一类的样本及其隶属度;否则返回s302。
[0126]
所述步骤s3中的模糊划分系数是判断模糊分类模糊性程度的标准,分类越明确时,模糊划分系数越大,模糊划分系数由以下公式确定:
[0127][0128]
所述步骤s3中的轮廓系数s是根据划分类别内样本点的最小组内距离和最大组间距离,是确定模糊分类算法中最终类别数k的一种方法,由以下公式确定:
[0129][0130]bi
=min{b
i1
,b
i2
,b
i3
,

,b
ik
}
[0131]
式中,si为单个样本点的轮廓系数;ai为同一类别内第i个样本点到其他样本点的平均距离;b
ik
为第i个样本点到其他类别内所有样本点的平均距离,轮廓系数(s)为所有样本点轮廓系数的平均值。
[0132]
实施例2
[0133]
本实施例给出实施例1的一个具体实施例,具体为:
[0134]
某双线盾构隧道位于广州市,隧道全长3.1km,盾构隧道管片设计内径8.0m、外径8.8m、管片环宽1.8m。采用两台土压平衡盾构进行施工,刀盘直径9.15m。盾构机主要穿越全风化花岗岩、中风化花岗岩和微风化花岗岩及其混合地层。由于隧道修建于市区,采用地质勘查孔取土试验获取盾构穿越的地质情况,并绘制地质剖面图。但是,由于两勘察孔之间间距较大,勘察孔之间的地质情况由勘察孔连线推测而得,地质情况不明确,无法精确制导盾构施工,本工程采用模糊分类算法的盾构掘进过程地质特征实时识别方法进行施工。
[0135]
步骤s2中将收集的盾构参数和二次变换参数共10
×
1197条数据进行主成分分析,提取前k个主成分,组成新的地质特征识别数据集y,数据库矩阵x由以下公式表示:
[0136][0137]
协方差矩阵为:
[0138][0139]
计算协方差矩阵的特征值和相应的特征向量,并根据特征向量计算每个因子所占的权重。按照特征向量大小将主成分从大到小排列,计算特征向量累积权重,前6个主成分的累计权重达到99%以上,因此可以将10个识别参数转化为6个主成分。
[0140]
将特征向量按照对应特征值的大小从左到右排列成一列矩阵,取前6列组成一个降维矩阵(p),确定的降维矩阵如下所示:
[0141][0142]
所述模糊分类算法具体为:
[0143]
步骤一:采用遍历的方法将分类类别数c分别设置为2~8,设置权重指数m为1.5,由输入的识别数据集y确定初始隶属度函数u0。
[0144]
步骤二:计算聚类中心vi,聚类中心由以下公式确定:
[0145][0146]
式中,u
ik
为隶属度函数,0《u
ik
≤1;n为样本数量(n=1197);c为分类类别数,1≤k≤n;1≤i≤c;m为权重指数;vi为第i类聚类中心;xk为第i类中的第k个样本。
[0147]
步骤三:根据下式更新隶属度矩阵:
[0148][0149]
式中,d
ik
=||x
k-vi||为第i类中第k个样本与第i类的聚类中心的欧氏距离。
[0150]
步骤四:给定一个较小的值(ε》0.0001),判断是否满足如果满足条件,则结束运算过程输出聚类中心、每一类的样本及其隶属度;否则返回步骤二。
[0151]
(2)图2给出了模糊划分系数、轮廓系数与分类类别数之间的变化关系。模糊化分系数由以下公式确定:
[0152][0153]
(3)根据划分类别内样本点的最小组内距离和最大组间距离,由以下公式确定轮廓系数(s):
[0154][0155]bi
=min{b
i1
,b
i2
,b
i3
,

,b
ik
}
[0156]
式中,si为单个样本点的轮廓系数;ai为同一类别内第i个样本点到其他样本点的平均距离;b
ik
为第i个样本点到其他类别内所有样本点的平均距离。轮廓系数(s)为所有样本点轮廓系数的平均值。
[0157]
(4)根据图2所示,当分类类别数为3时,模糊划分系数最大,且轮廓系数最大;结合地质剖面图,确定地质特征最终类别数k为3。
[0158]
根据确定的地质特征最终类别数(k),结合地质剖面图,更新地质特征的名称。具体的,地质特征名称更新为:软土地层、软硬不均地层、硬岩地层。并将软土地层赋值为1,命名为类别1;软硬不均地层赋值为2,命名为类别2;硬岩地层赋值为3,命名为类别3。
[0159]
将地质特征最终类别数(k)、权重指数(m)、初始隶属度函数(u0)和降维后的识别数据集(y)输入到模糊分类算法,确定每一盾构掘进环所属的地质特征类别,并与盾构掘进记录的地质特征及地质剖面图对比,确定地质特征识别的准确性。根据最终类别数(k=3)和权重指数(m=1.5),确定每一环的地质特征类型,如图3所示,根据主成分1和主成分2绘制关于地质特征识别结果的2维特征空间图。将识别结果与记录的地质特征类型对比,准确度达到91.6%,满足盾构隧道施工要求。
[0160]
如图4所示,将新收集的盾构实时参数经第二步的降维处理后,输入到调整后的地质特征识别模型中,确定盾构掘进所在的地质特征类别并输出所属地质特征类别的隶属
度。图4中的地质特征类型与施工效率作对比,在软土地层中,盾构掘进速度较快;在硬岩地层中,盾构掘进速度最慢。现场工程师可根据图4的识别结果对施工情况进行评估,安排施工情况。
[0161]
实施例3
[0162]
本实施例提供一种基于模糊分类算法的盾构掘进地质特征识别系统,如图5所示,包括:
[0163]
收集模块,所述收集模块收集隧道施工区间地质资料、地质特征识别参数和出土情况,所述地质资料包括地质剖面图,根据出土情况得到初步的地质特征;
[0164]
特征提取模块,所述特征提取模块利用主成分分析方法对所述地质特征识别参数进行特征提取;
[0165]
类别数确定模块,所述类别数确定模块将特征提取后得到的数据集采用模糊划分系数fpc、轮廓系数s结合地质剖面图确定地质特征最终类别数k;
[0166]
更新模块,所述更新模块根据确定的地质特征最终类别数k,结合地质剖面图,更新地质特征;
[0167]
训练模块,所述训练模块利用模糊分类算法,确定每一盾构掘进环所属的地质特征类别,并与盾构掘进记录的地质特征及地质剖面图对比,确定地质特征识别的准确性,若准确性达到阈值以上,则该识别模型可用于地质特征识别;否则,调整模糊分类算法的参数,直到满足要求为止;
[0168]
识别模块,所述识别模块将新收集的盾构实时参数经步骤s2处理后,输入到调整后的地质特征识别模型中,确定盾构掘进所在的地质特征类别并输出所属地质特征类别的隶属度。
[0169]
相同或相似的标号对应相同或相似的部件;
[0170]
附图中描述位置关系的用语仅用于示例性说明,不能理解为对本专利的限制;
[0171]
显然,本发明的上述实施例仅仅是为清楚地说明本发明所作的举例,而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说,在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明权利要求的保护范围之内。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1