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

一种基于滑窗检测提取电流信号的断路器故障诊断方法

2022-07-17 00:40:53 来源:中国专利 TAG:


1.本技术涉及高压断路器故障诊断技术领域,具体涉及一种基于滑窗检测提取电流信号的断路器故障诊断方法。


背景技术:

2.高压断路器作为在电力系统中起到保护和控制作用的重要关键设备,诊断其工作运行状态对于维护电力系统的安全运行有着重要意义。高压断路器分合闸线圈的电流信号可以不同程度地反映断路器铁芯卡涩、铁芯空行程过大、线圈匝间短路等工作情况,常用来评估断路器的工作状态。
3.以往在对断路器工作状态的信号提取的研究中,提取断路器分合闸线圈电流波形的关键点特征进行故障分类,可实现波形特征差异较大的识别。一种基于振声联合诊断法利用集合经验模式分解对断路器进行振声信号的特征提取进行故障诊断,取得了良好的效果。而且近年来机器学习算法在故障分类、辨识、预测等领域取得的巨大突破,也逐渐被用于断路器故障分析,如神经网络、随机森林、支持向量机等分类模型。但在上述研究中进行断路器特征提取时不能自适应各类断路器,且传统故障分类模型的准确率不理想,导致诊断精度和效率低下。


技术实现要素:

4.本发明的目的在于,提供一种基于滑窗检测提取电流信号的断路器故障诊断方法,以滑窗检测法提取关键时刻点、电流值,并构建分合闸电流综合特征集,以此建立主成分分析pca—支持向量机svm的故障分类模型,并采用网格搜索法gs优化模型有效提高诊断精度和效率,并且可自适应多型号的断路器故障识别。
5.本发明采取的技术方案是:一种基于滑窗检测提取电流信号的断路器故障诊断方法,包括如下步骤:s1:利用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻点和对应的电流值,所述关键时刻点和对应的电流值包括电磁铁铁芯开始运动时刻t1、电磁铁铁芯开始运动时刻的电流极大值i1、铁芯运动至最大行程时刻t3、铁芯运动至最大行程时刻的电流极小值i3、断路器辅助节点断开时刻t5和断路器辅助节点断开时刻的电流稳态值i5;并结合线圈工作阶段持续时间和以及电流信号的总体时域统计特征,建立分合闸电流综合特征样本集;s2:利用主成分分析法pca对分合闸电流综合特征样本集进行特征降维并依据累计贡献度大小选取主成分,得到降维后的综合特征样本集;s3:将降维后的综合特征样本集作为输入样本集输入pca-svm模型,采用k折交叉验证法将输入样本集划分为训练集和测试集,使用训练集对pca-svm模型进行模型训练,并在模型性训练过程中采用网格搜索法gs对模型的惩罚系数c和核函数参数g进行参数优化,将pca-svm模型优化为pca-gs-svm故障分类模型,最后使用测试集对pca-gs-svm故障分类
模型进行性能验证。
6.进一步地,所述步骤s1中利用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻点和对应电流值的具体步骤为:s101:设定一个长度为d0、高度为h0的检测窗口,窗口中心为断路器的分合闸线圈电流曲线中的点;沿窗口中心变化方向,即断路器的分合闸线圈电流曲线方向以步长l移动检测窗口,形成滑动检测窗口;s102:从分合闸线圈开始带电时刻t0开始以初始步长l0进行极大值点检测,当检测到第一个极大值点a后,将第一个极大值点a对应的时刻和电流极大值分别记作电磁铁铁芯开始运动时刻t1和电磁铁铁芯开始运动时刻的电流极大值i1;s103:从第一个极大值点a开始以第一步长l1,在第一时长范围s1内进行极小值点检测,当检测到第一个极小值点b后,将第一个极小值点b对应的时刻记作tb;s104:从第一个极小值点b开始以第二步长l2,在时长范围sb内进行极大值点检测;若在时长范围sb内检测到极大值点,则将tb记作铁芯撞击弯板时刻t2,并从极大值点开始继续以第二步长l2,在第二时长范围s2中进行极小值点检测,并将检测到的第二个极小值点c对应的时刻和电流极小值分别记作铁芯运动至最大行程时刻t3和铁芯运动至最大行程时刻的电流极小值i3;再从第二个极小值点c开始以第二步长l2,在第三时长范围s3内进行拐点检测,检测是否存在第一拐点d;若未在时长范围sb内检测到极大值点,则将tb和tb时刻对应的电流极大值分别记作铁芯运动至最大行程时刻t3和铁芯运动至最大行程时刻的电流极小值i3;并继续以第二步长l2,在第三时长范围s3内进行拐点检测,检测是否存在第一拐点d;s105:若检测到第一拐点d,则将第一拐点d对应的时刻记作断路器主触头分合闸完成时刻t4,并从第一拐点d开始以第三步长l3进行极大值点检测,直至检测到极大值点e,将极大值点e对应的时刻和电流稳态值分别记作断路器辅助节点断开时刻t5和断路器辅助节点断开时刻的电流稳态值i5;若未检测到第一拐点d,则继续以第三步长l3进行第三阶段检测,直至检测到极大值点e,将极大值点e对应的时刻和电流稳态值分别记作断路器辅助节点断开时刻t5和断路器辅助节点断开时刻的电流稳态值i5。
7.进一步地,所述第一时长范围s1为5ms,时长范围sb为3ms,第二时长范围s2为5ms,第三时长范围s3为30ms;初始步长l0为d0,第一步长l1为0.2ms,第二步长l2为0.3ms,第三步长l3为0.5ms。
8.进一步地,所述电流信号的总体时域统计特征包括电流信号的平均值、电流信号的标准偏差、电流信号分布偏斜方向和程度的度量,即偏度skew和峰值在电流波形中极端程度,即峰值因子cf,具体表达式如下:
其中,xi为第i个电流信号变量的取值,i =1,2,
……
,n,n为随机变量个数。
9.进一步地,所述步骤s2的具体步骤为:s201:对断路器线圈电流信号构建信号样本特征维度向量,所述断路器线圈电流信号有m个样本,每个样本有n个特征值,样本特征维度向量;s202:根据样本特征维度向量构建标准化后的协方差矩阵,其中,为样本的均值向量;s203:求取协方差矩阵rd的特征值λi与对应的特征向量vi,并将特征值λi按照大小顺序进行排列,依次记为λ1,λ2,
……
λ
n-1
,λn,其中λ
1 >λ
2 >
……
>λ
n-1 >λn,对应的特征向量依次记为v1,v2,
……vn-1
,vn,构建特征向量矩阵;s204:当累计贡献率时,则提取特征向量矩阵v的前k个主成分组成一个新的主成分矩阵p;利用主成分矩阵p将样本特征维度向量dm进行降维,转换为降维后的矩阵,降维后的矩阵x为m
×
k的矩阵。
10.进一步地,所述步骤s3的具体步骤为:s301:采用高斯径向基核函数rbf作为pca-gs-svm故障分类模型的核函数进行训练;确定惩罚系数c和核函数参数g的选择集,令,,其中c
strat
表示惩罚系数c待搜索范围的上限,c
end
表示惩罚系数c待搜索范围的下限,g
strat
表示核函数参数g待搜索范围的上限,g
end
表示核函数参数g待搜索范围的下限;通过网格搜索法gs,将惩罚系数c和核函数参数g在搜索范围集内的取值进行排列组合,列出所有可能的组合结果生成二维网格,记某一网格的坐标为(i,j),对应的网格参数组合为;s302:将输入样本集划分为k个子集,任取一个作为测试集,余下的k-1个子集作为训练集,利用训练集对pca-gs-svm故障分类模型进行训练;并使用测试集验证训练好的pca-gs-svm故障分类模型的测试结果,统计当前测试集下测试结果的测试评分sk,即预测准确率,其中,k=1,2,3,
……
,k;
s303:重新选定子集作为测试集,重复步骤s302,直至k个子集中的每一个子集都作为了测试集进行测试,取k个测试评分的均值作为当前网格参数组合下的评分;s304:对网格参数组合进行调整,重复步骤s302~303,直至网格遍历完成,选取评分最优的那组对应的网格参数组合作为最优参数组合,完成pca-gs-svm故障分类模型的训练。
11.本发明的有益效果在于:(1)对于部分断路器,铁芯撞击弯板时刻t2和断路器主触头分合闸完成时刻t4难以观测或是缺少这两个时刻点,导致现有故障诊断方法难以自适应不同型号的断路器;本发明采用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻点和对应的电流值,提取的关键点类型都符合实际波形类型,且误差率较低,具有良好的断路器电流波形关键特征辨识性能,且能自适应不同型号断路器波形;(2)在提升提取关键时刻点和对应的电流值的准确率的基础上,结合线圈工作阶段持续时间和以及电流信号的总体时域统计特征,利用主成分分析pca对构建的综合特征集进行降维,进一步提升本发明的故障识别准确率,同时缩短了故障诊断时间,提升了故障诊断效率;(3)本发明采用k折交叉验证法结合网格搜索法gs对现有的pca-svm模型进行优化,得到本发明所用的pca-gs-svm故障分类模型,与常规的svm模型、pca-svm模型和pca-bpnn模型相比,本发明的pca-gs-svm故障分类模型的识别准确率有显著提升。
附图说明
12.为了更清楚地说明本发明实施例中的技术方案,下面将对实施例中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本技术的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其它的附图。
13.图1为本发明实施例的方法流程图;图2为典型断路器分合闸线圈电流波形;图3为四种常见型号断路器分合闸线圈电流标准波形;图4为本发明实施例滑窗检测法提取的关键时刻点的示意图;图5为本发明实施例滑窗检测法的方法流程图;图6为四种典型故障的电流波形图;图7为本发明实施例网格搜索法gs优化惩罚系数c和核函数参数g的结果示意图;图8为本发明实施例在四种常见型号断路器中关键时刻点的提取值和实际值的对比图;图9为本发明实施例的pca-gs-svm故障分类模型的诊断结果示意图。
具体实施方式
14.为了能够更清楚地理解本发明的上述目的、特征和优点,下面结合附图和具体实施方式对本发明进行进一步的详细描述。在下面的描述中阐述了很多具体细节以便于充分理解本发明,但是,本发明还可以采用其他不同于在此描述的其他方式来实施,因此,本发明并不限于下面公开的具体实施例的限制。
15.除非另作定义,此处使用的技术术语或者科学术语应当为本技术所述领域内具有一般技能的人士所理解的通常意义。本专利申请说明书以及权利要求书中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性,而只是用来区分不同的组成部分。同样,
ꢀ“
一个”或者“一”等类似词语也不表示数量限制,而是表示存在至少一个。
ꢀ“
连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接,而是可以包括电性的连接,不管是直接的还是间接的。
ꢀ“
上”、
ꢀ“
下”、
ꢀ“
左”、
ꢀ“
右”等仅用于表示相对位置关系,当被描述对象的绝对位置改变后,则该相对位置关系也相应地改变。
16.如图1所示,一种基于滑窗检测提取电流信号的断路器故障诊断方法,包括如下步骤:s1:利用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻点和对应的电流值,所述关键时刻点和对应的电流值包括电磁铁铁芯开始运动时刻t1、电磁铁铁芯开始运动时刻的电流极大值i1、铁芯运动至最大行程时刻t3、铁芯运动至最大行程时刻的电流极小值i3、断路器辅助节点断开时刻t5和断路器辅助节点断开时刻的电流稳态值i5;并结合线圈工作阶段持续时间和以及电流信号的总体时域统计特征,建立分合闸电流综合特征样本集。
17.所述电流信号的总体时域统计特征包括电流信号的平均值、电流信号的标准偏差、电流信号分布偏斜方向和程度的度量,即偏度skew和峰值在电流波形中极端程度,即峰值因子cf,具体表达式如下:具体表达式如下:具体表达式如下:具体表达式如下:其中,xi为第i个电流信号变量的取值,i =1,2,
……
,n,n为随机变量个数。
18.典型断路器分合闸线圈电流波形如图2所示,t0为分合闸线圈开始带电时刻;t1表示电磁铁铁芯开始运动时刻,对应电流极大值i1;t2是铁芯撞击弯板时刻,如果电流出现局
部增加,可以明显观察到t2,否则不容易观察到;t3表示铁芯运动至最大行程时刻,对应电流极小值i3;t4为断路器主触头分合闸完成时刻,不容易观察到;t5为断路器辅助节点断开时刻,对应电流稳态值i5。t
0 ~ t1为电磁铁线圈上电到电磁铁铁芯运动前的阶段,由于复位弹簧的反作用力,衔铁或铁芯顶杆并未移动,电流按规律升高;t
1 ~ t3为铁芯运动时刻,即电磁铁吸合时间,包含了铁芯顶杆运动的空程阶段t
1 ~ t2和铁芯顶杆脱扣阶段t
2 ~ t3,当电磁铁吸力较小时可以容易观察到t2;t
3 ~ t5为铁芯停止运动,线圈电流恢复到断路器辅助节点断开的时间。其中时间点t1、t3、t5与相对应的电流节点值i1、i3、t5可以反映线圈电流、电阻和电磁铁运动状态等物理信息,为本发明实施例所要提取的关键特征。
19.以如图3所示的四种市面上常见的不同型号断路器为例,图3中(a)~(d)依次是型号为zf11-252(l)、hpl550b2、lw13a-500y和lw12-500的断路器的分合闸线圈电流标准波形图。从图3中(a)可知,zf11-252(l)断路器的t
0 ~ t5均可较为明显地观察到,与图2所示的典型断路器分合闸线圈电流波形特征点较类似;从图3中(b)可知,hpl550b2断路器的铁芯撞击弯板时刻t2难以观测;从图3中(c)可知,lw13a-500y断路器缺少铁芯撞击弯板时刻t2和断路器主触头分合闸完成时刻t4;从图3中(d)可知,lw12-500断路器缺少断路器主触头分合闸完成时刻t4。对于实际采集到波形中的离散域电流采样点值,可以根据相邻两点斜率的变化情况,检测如极值点和拐点等关键点,但可能出现t2、t3、t4观测点混肴,判断不准确的情况。
20.为了避免混肴t2、t3和t4,本发明实施例采用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻点和对应电流值,具体步骤为:s101:设定一个长度为d0、高度为h0的检测窗口,窗口中心为断路器的分合闸线圈电流曲线中的点;沿窗口中心变化方向,即断路器的分合闸线圈电流曲线方向以步长l移动检测窗口,形成滑动检测窗口。本发明实施例设置检测窗口的长度d
0 = 0.5,高度h
0 = 0.2。
21.s102:从分合闸线圈开始带电时刻t0开始以初始步长l0进行极大值点检测,当检测到第一个极大值点a后,将第一个极大值点a对应的时刻和电流极大值分别记作电磁铁铁芯开始运动时刻t1和电磁铁铁芯开始运动时刻的电流极大值i1。
22.s103:从第一个极大值点a开始以第一步长l1,在第一时长范围s1内进行极小值点检测,当检测到第一个极小值点b后,将第一个极小值点b对应的时刻记作tb。
23.s104:从第一个极小值点b开始以第二步长l2,在时长范围sb内进行极大值点检测;若在时长范围sb内检测到极大值点,则将tb记作铁芯撞击弯板时刻t2,并从极大值点开始继续以第二步长l2,在第二时长范围s2中进行极小值点检测,并将检测到的第二个极小值点c对应的时刻和电流极小值分别记作铁芯运动至最大行程时刻t3和铁芯运动至最大行程时刻的电流极小值i3;再从第二个极小值点c开始以第二步长l2,在第三时长范围s3内进行拐点检测,检测是否存在第一拐点d;若未在时长范围sb内检测到极大值点,则将tb和tb时刻对应的电流极大值分别记作铁芯运动至最大行程时刻t3和铁芯运动至最大行程时刻的电流极小值i3;并继续以第二步长l2,在第三时长范围s3内进行拐点检测,检测是否存在第一拐点d。
24.s105:若检测到第一拐点d,则将第一拐点d对应的时刻记作断路器主触头分合闸完成时刻t4,并从第一拐点d开始以第三步长l3进行极大值点检测,直至检测到极大值点e,将极大值点e对应的时刻和电流稳态值分别记作断路器辅助节点断开时刻t5和断路器辅助
节点断开时刻的电流稳态值i5;若未检测到第一拐点d,则继续以第三步长l3进行第三阶段检测,直至检测到极大值点e,将极大值点e对应的时刻和电流稳态值分别记作断路器辅助节点断开时刻t5和断路器辅助节点断开时刻的电流稳态值i5。
25.由于不同目标检测点之间有范围差异,故本发明实施例在每个检测阶段采用不同步长l来进行分段滑窗检测。断路器分合闸电流信号的一般持续时间为40~120ms。其中,对于绝大部分断路器,如具有铁芯撞击弯板时刻t2,则t2时刻点与下一个靠近的极大值时刻点间隔范围不超过3ms,t
1 ~ t2和t
2 ~ t3不会超过5ms;如不具有铁芯撞击弯板时刻t2,则t
1 ~ t3不会超过5ms;t
3 ~ t4不超过30ms。因此在本发明实施例中,所述第一时长范围s1为5ms,时长范围sb为3ms,第二时长范围s2为5ms,第三时长范围s3为30ms;初始步长l0为d0,d0在[0.1,2]范围内,第一步长l1为0.2ms,第二步长l2为0.3ms,第三步长l3为0.5ms。
[0026]
s2:利用主成分分析法pca对分合闸电流综合特征样本集进行特征降维并依据累计贡献度大小选取主成分,得到降维后的综合特征样本集;具体步骤为:s201:对断路器线圈电流信号构建信号样本特征维度向量,所述断路器线圈电流信号有m个样本,每个样本有n个特征值,样本特征维度向量。
[0027]
s202:根据样本特征维度向量构建标准化后的协方差矩阵,其中,为样本的均值向量。
[0028]
s203:求取协方差矩阵rd的特征值λi与对应的特征向量vi,并将特征值λi按照大小顺序进行排列,依次记为λ1,λ2,
……
λ
n-1
,λn,其中λ
1 >λ
2 >
……
>λ
n-1 >λn,对应的特征向量依次记为v1,v2,
……vn-1
,vn,构建特征向量矩阵。
[0029]
s204:当累计贡献率时,则提取特征向量矩阵v的前k个主成分组成一个新的主成分矩阵p;利用主成分矩阵p将样本特征维度向量dm进行降维,转换为降维后的矩阵,降维后的矩阵x为m
×
k的矩阵。
[0030]
s3:将降维后的综合特征样本集作为输入样本集输入pca-svm模型,采用k折交叉验证法将输入样本集划分为训练集和测试集,使用训练集对pca-svm模型进行模型训练,并在模型性训练过程中采用网格搜索法gs对模型的惩罚系数c和核函数参数g进行参数优化,将pca-svm模型优化为pca-gs-svm故障分类模型,最后使用测试集对pca-gs-svm故障分类模型进行性能验证;具体步骤为:s301:采用高斯径向基核函数rbf作为pca-gs-svm故障分类模型的核函数进行训练;确定惩罚系数c和核函数参数g的选择集,令,,其中c
strat
表示惩罚系数c待搜索范围的上限,c
end
表示惩罚系数c待搜索范围的下限,g
strat
表示核函数参数g待搜索范围的上限,g
end
表示核函数参数g待搜索范围的下限;通过网格搜索法gs,将惩罚系数c和核函数参数g在搜索范围集内的取值进行排列组合,列出所有可能
的组合结果生成二维网格,记某一网格的坐标为(i,j),对应的网格参数组合为。
[0031]
s302:将输入样本集划分为k个子集,任取一个作为测试集,余下的k-1个子集作为训练集,利用训练集对pca-gs-svm故障分类模型进行训练;并使用测试集验证训练好的pca-gs-svm故障分类模型的测试结果,统计当前测试集下测试结果的测试评分sk,即预测准确率,其中,k=1,2,3,
……
,k。
[0032]
s303:重新选定子集作为测试集,重复步骤s302,直至k个子集中的每一个子集都作为了测试集进行测试,取k个测试评分的均值作为当前网格参数组合下的评分。
[0033]
s304:对网格参数组合进行调整,重复步骤s302~303,直至网格遍历完成,选取评分最优的那组对应的网格参数组合作为最优参数组合,完成pca-gs-svm故障分类模型的训练。
[0034]
在本发明实施例中,令c1=[-8,8]、d1=[-8,8],c1和d1的步长均为0.25,惩罚系数c的范围设置为,核函数参数g的范围设置为,k折验证取5折,可得到如图7所示的结果示意图。从图7可知,最高测试评分点所对应最优评分score为0.971,所对应的最优参数惩罚系数c和核函数参数g分别为15.31和2.20。
[0035]
本发明实施例选取4种断路器故障状态进行诊断研究:分别是电源电压过低、铁芯卡涩、铁芯空行程大和线圈匝间短路,每种故障典型故障波形如图6所示,图6中(a)为电源电压过低故障波形,(b)为铁芯卡涩故障波形,(c)为铁芯空行程大故障波形,(d)为线圈匝间短路故障波形。
[0036]
将正常状态、电源电压过低、铁芯卡涩、铁芯空行程大和线圈匝间短路分别设置标签0~4,经过步骤s1得到的关键特征数据维数为6维,结合总体时域统计特征拓展为12维,最终得到的分合闸电流综合特征样本集如表1所示。
[0037]
表1分合闸电流综合特征样本集通过现场收集、实验得到200组不同故障样本特征数据,200组样本特征数据按不同断路器型号划分结果如表2所示,按不同故障样本划分结果如表3所示。
[0038]
表2不同断路器型号划分表
表3 不同故障样本划分表图8为本发明实施例在四种常见型号断路器中关键时刻点的提取值和实际值的对比图,图8中各个关键时刻点对应的数值如表4所示,表4中t为实际值,t' 为采用本发明实施例所述的方法得到的提取值,误差率为,通过误差率可判断基于滑窗检测的波形关键点自适应提取性能。
[0039]
表4 本发明实施例关键时刻点的提取值和实际值对比表
将5
×
12组测试样本输入已建立的pca-gs-svm故障分类模型进行故障诊断,pca-gs-svm故障分类模型的诊断结果如图9所示。同时为了更好地体现本模型故障诊断的优越性,将本发明实施例的pca-gs-svm故障分类模型与传统的断路器故障分类模型进行诊断准确率和诊断时间的综合比较,分别使用svm模型、pca-svm模型、pca-bpnn模型与本发明实施例的pca-gs-svm故障分类模型进行测试,得到如表5所示的诊断结果对比表。
[0040]
表5 四种模型诊断结果对比表从表5中可以看出,本发明实施例的pca-gs-svm故障分类模型在诊断准确率上相较于pca-bpnn模型、pca-svm模型和svm模型具有显著提升,pca-gs-svm故障分类模型的诊断时间虽然略长于pca-svm模型和svm模型,但综合两方面性能,本发明实施例pca-gs-svm故障分类模型相较于pca-bpnn模型、pca-svm模型和svm模型,还是具有明显的优势。
[0041]
本发明实施例采用滑窗检测法提取断路器的分合闸线圈电流曲线中的关键时刻
点和对应的电流值,提取的关键点类型都符合实际波形类型,且误差率较低,具有良好的断路器电流波形关键特征辨识性能,且能自适应不同型号断路器波形;在提升提取关键时刻点和对应的电流值的准确率的基础上,利用主成分分析pca对构建的综合特征集进行降维,进一步提升本发明的故障识别准确率,同时缩短了故障诊断时间,提升了故障诊断效率;并通过k折交叉验证法结合网格搜索法gs对现有的pca-svm模型进行优化,得到本发明所用的pca-gs-svm故障分类模型,与常规的svm模型、pca-svm模型和pca-bpnn模型相比,本发明的pca-gs-svm故障分类模型的识别准确率有显著提升。
[0042]
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
再多了解一些

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

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

相关文献