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

一种锂电池等效电路模型参数的辨识方法

2022-11-19 08:52:30 来源:中国专利 TAG:


1.本发明涉及锂电池技术领域,特别涉及一种锂电池等效电路模型参数的辨识方法,该辨识方法利用混沌自适应分数阶粒子群算法对锂电池等效电路模型的参数进行辨识。


背景技术:

2.促进混合动力或纯电动汽车的发展是缓解环境污染、石油短缺等问题的主要战略措施之一。锂电池由于其高能量密度、低放电率和高循环次数,已成为新能源汽车动力系统的核心组成部分之一。然而,锂电池的安全性和功率性能受到许多因素的影响,如温度和荷电状态(soc)、健康状态(soh)和功率状态(sop)估计精度。准确的soc估计支持可用的容量监测和相关的状态反馈控制,这对维持锂电池的安全和高效运行尤为重要。准确的状态估计的前提是要有一个能够准确捕捉电池系统动态特性的模型。当前锂电池模型主要有电化学模型、数据驱动模型和等效电路模型。电化学模型准确度高,但模型过于复杂参数过多,工程实现难度大。数据驱动模型其精度依赖大量数据训练,可移植性差。等效电路模型由其结构简单,便于实现且精度较高,在工程应用中得到了广泛的使用。常见的等效电路模型有pngv模型、rc等效电路模型和分数阶等效电路模型。研究和应用最广泛的是多阶rc等效电路模型和分数阶等电路模型。高阶次的rc等效电路模型精度高,但是模型也更加复杂需要辨识的参数多。分数阶rc等效电路模型拥有更高的精度和不高的复杂度得到了广泛的应用。对于等效电路模型说其模型精度非常依赖于参数辨识。
3.粒子群优化算法通过设计一种无质量的粒子来模拟鸟群中的鸟,粒子仅具有两个属性:速度和位置,速度代表移动的快慢,位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解,并将其记为当前个体极值,并将个体极值与整个粒子群里的其他粒子共享,找到最优的那个个体极值作为整个粒子群的当前全局最优解,粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。传统粒子群优化算法收敛速度慢,种群初始化多采用随机初始化的方法,这样种群初始分布容易扎堆,容易陷入局部最优解而找不到全局最优解。
4.综上所述,目前采用粒子群算法对锂电池建模中参数的辨识还存在以下问题:
5.1.在传统粒子群优化算法中,初始种群一般采用随机生成的策略,因此个体将不可避免地分布不均匀,导致整个搜索过程中的多样性不足。这个基本上限制了搜索性能。
6.2.若采用传统粒子群优化算法对锂电池进行建模,不改进很容易陷入局部最优解,遗落全局最优解,且传统粒子全优化算法收敛速度慢。


技术实现要素:

7.本发明要解决的技术问题为针对现有技术的不足,为了克服现有锂电池建模中参数辨识精度不高的缺陷,提供了一种基于混沌自适应分数阶粒子群优化算法的锂电池参数的辨识方法。具有收敛速度快、辨识精度高、不易陷入局部最优解等优点。提高了锂电池等
效电路模型的精度。
8.为解决发明的技术问题,所采用的技术方案为:一种锂电池等效电路模型参数的辨识方法,所述辨识方法利用混沌自适应分数阶粒子群算法对锂电池等效电路模型的参数进行辨识,然后将得到的最优参数反馈到锂电池等效电路模型中,具体步骤如下:
9.步骤1,获取锂电池等效电路模型的辨识参数
10.所述锂电池电路等效模型的拓扑结构包括一个直流电源、一个欧姆内阻、一个并网电阻和一个分数阶常相位角元件cpe,所述并网电阻和分数阶常相位角元件cpe并联,并记为并联支路,所述欧姆内阻的一端与直流电源的正极相连,另一端与并联支路的一端相连,所述并联支路的另一端与直流电源的负极组成所述锂电池电路等效模型的输出端;
11.将所述直流电源处的电压记为开路电压u
ocv
,所述欧姆内阻的阻值记为欧姆内阻值r0,所述并网电阻的阻值记为并网电阻值r1,所述锂电池等效电路模型的输出端的电压记为端电压u
t

12.将所述分数阶常相位角元件cpe记为元件cpe,元件cpe的阻抗传递函数的表达式为:其中,z(s)为元件cpe的阻抗,c1为元件cpe的电容值,记为电容值c1,s为拉普拉斯算子,β为元件cpe的分数阶阶次,记为阶次β,当β=1时,元件cpe为理想的电容,当β=0时,元件cpe等价为纯电阻;
13.所述锂电池等效电路模型的辨识参数为:欧姆内阻值r0,并网电阻值r1、电容值c1和阶次β,即在混沌自适应分数阶粒子群算法中,目标搜索空间的维数为4,将欧姆内阻值r0,并网电阻值r1、电容值c1和阶次β分别记为辨识参数1、辨识参数2、辨识参数3和辨识参数4;
14.将所述辨识参数中的每一个辨识参数视为粒子个体的一个维度,所述粒子个体有粒子个体位置与粒子个体速度两个属性,粒子个体位置代表移动的方向,粒子个体速度代表移动的快慢,其中,粒子个体位置为辨识参数的值;
15.步骤2,混沌自适应分数阶粒子群算法参数的设置
16.设定目标搜索空间的维数为d,d=4;设一个种群由n个粒子组成;设定最大迭代次数m;设定粒子惯性权重ω;
17.设定粒子个体的位置范围为{x_min,x_max}、粒子个体的速度范围{v_min,v_max}、粒子惯性权重范围为{ω_min,ω_max},其中,x_min为粒子个体位置的最小值,x_max为粒子个体位置的最大值,v_min为粒子个体速度的最小值,v_max为粒子个体速度的最大值,ω_min为粒子惯性权重最小值,ω_max为粒子惯性权重最大值;
18.将n个粒子个体中的任意一个记为粒子个体i,i为种群中任意一个粒子个体的序号,i=1,2,..n,将粒子个体i的位置向量记为个体位置xi,粒子个体i的速度向量记为个体速度vi,其表达式分别如下:
19.xi={x
i1
,x
i2
,...,x
ij
...,x
id
}
[0020]vi
={v
i1
,v
i2
,...v
ij
...,v
id
}
[0021]
其中,个体位置xi为锂电池等效电路模型的一组辨识参数的解,x
ij
为粒子个体i第j维的位置,j=1,2,..d;个体速度vi为锂电池等效电路模型的一组参数的解在粒子搜索解空间中的速度,v
ij
为粒子个体i第j维的速度;
[0022]
将粒子个体i的适应度值记为个体适应度值fi,将粒子个体i搜索到的适应度值最小的位置记为个体最优位置p
best
,p
best
={p
i1
,p
i2
,...p
ij
...,p
id
},其中,p
ij
为粒子个体i的第j维搜索到的最优位置,将整个粒子群搜索到的适应度值最小的位置记为全局最优位置g
best
,g
best
={p
g1
,p
g2
,...,p
gj
...,p
gd
},p
gj
为第j维度整个粒子群搜索到的全局最优位置;
[0023]
步骤3,混沌自适应分数阶粒子群算法的初始化
[0024]
设迭代次数为k,k=1,2,..m;
[0025]
使用logistic混沌映射来初始化粒子群前4次迭代中的粒子个体位置和粒子个体速度;所述logistic混沌映射表达式如下:
[0026]zij
=2z
ij
(1-z
ij
)
[0027]
其中,z
ij
为粒子个体i第j维的0-1之间的随机数,z
ij
为粒子个体i第j维生成的混沌参数;
[0028]
与目标搜索空间的维数d对应,由logistic混沌映射表达式生成混沌序列{z
i1
,z
i2
,...z
ij
...,z
id
},且通过下式映射到粒子群的目标搜索空间:
[0029][0030]
式中,为第k次迭代中的粒子个体i第j维的位置,为第k次迭代中的粒子个体i第j维的速度,x
ij
_max为粒子个体i第j维位置的最大值,x
ij
_min为粒子个体i第j维位置的最小值,v
ij
_max为粒子个体i第j维速度的最大值,v
ij
_min为粒子个体i第j维速度的最小值,k=1,2,3,4;
[0031]
由此得到混沌映射粒子个体初始化位置和混沌映射粒子个体初始化速度随机初始化粒子个体i最优位置随机初始化粒子群全局最优位置k=1,2,3,4,其表达式分别为:
[0032][0033][0034][0035]
在粒子惯性权重范围{ω_min,ω_max}内随机初始化粒子惯性权重k=1,2,3,4;
[0036]
步骤4,第5次迭代-第m次迭代
[0037]
将第k次迭代记为当前迭代,第k-1次迭代记为前一迭代,第k-2次迭代记为前二迭代,第k-3次迭代记为前三迭代,第k-4次迭代记为前四迭代;
[0038]
步骤4.1,个体位置和个体速度的更新
[0039]
通过分数阶速度更新方程进行个体位置和个体速度的更新,得到粒子个体i在当前迭代中的速度和粒子个体i在当前迭代中的位置所述分数阶速度更新方程的表达式如下:
[0040][0041]
其中,为粒子个体i在前一迭代中的速度,为粒子个体i在前二迭代中的速度,为粒子个体i在前三迭代中的速度,为粒子个体i在前四迭代中的速度,为粒子个体i在前一迭代中的位置,r1为加速度权重系数1,r2为加速度权重系数2,为粒子个体i在前一次迭代中的迭代惯性权重系数,α为分数阶阶次,取α=0.5;为粒子个体i在前一次迭代中的最优位置,记为前一个体最优位置为粒子群在前一迭代中的全局最优位置,记为前一全局最优位置
[0042]
步骤4.2,惯性权重的更新
[0043]
记粒子个体i在当前迭代中的惯性权重为当前惯性权重表达式如下:
[0044][0045]
其中,为粒子个体i在前一迭代中的适应度值,为前一迭代中种群所有粒子个体适应度值的平均值,为前一迭代中种群所有粒子个体适应度值中的最小值;
[0046]
步骤4.3,计算粒子个体i在当前迭代中的适应度值,并记为当前个体最优位置将所述前一个体最优位置对应的适应度值记为前一个体最优适应度值前一全局最优位置对应的适应度值记为前一全局最优适应度值进行如下判断:
[0047]
比较和取适应度值最小者的位置为当前个体最优位置并将当前个体最优位置对应的适应度值记为当前个体最优适应度值
[0048]
比较和取适应度值最小者的位置为当前模拟全局最优位置并将当前模拟全局最优位置对应的适应度值记为当前模拟全局最优适应度值为
[0049]
步骤4.4,为避免陷入局部最优,对当前模拟全局最优位置进行混沌搜索,所述混沌搜索的表达式如下:
[0050][0051]
其中,τ为混沌权重,zr是由logistic混沌映射表达式生成的混沌序列,xr为混沌搜索产生的新的粒子位置;
[0052]
通过混沌搜索的表达式对当前迭代模拟全局最优位置进行10次混沌搜索,产生10组粒子个体的位置,计算10组粒子个体的适应度值,并与当前模拟全局最优适应度值为比较,取其中适应度值最小的位置为当前全局最优位置该当前全局最优位置对应一组辨识参数;
[0053]
步骤4.5,进行如下判断:
[0054]
若k=m,结束迭代,输出当前全局最优位置并记为全局最终位置g,该全局最终位置g对应的一组辨识参数即为锂电池等效电路模型的一组最优参数;
[0055]
若k<m,返回步骤4.1进行下一次迭代。
[0056]
优选地,步骤4中所述适应度伯的计算式如下:
[0057][0058][0059][0060]
其中,t为锂电池工况数据集中数据点的个数,m为锂电池工况数据集中数据点的序号,m=1,2,...t,u
t
(m)为锂电池工况数据集第m个数据点对应的端电压;为将前一迭代中得到的锂电池等效电路模型辨识参数带入锂电池等效电路模型中计算出的端电压,u
ocv
(m)为锂电池工况数据集第m个数据点对应的开路电压,i(m)为锂电池工况数据集第m个数据点对应的电流,为前一迭代中的cpe元件的电压值,β
k-1
为前一迭代中元件cpe的分数阶阶次,为β
k-1
阶导数,为前一迭代中的欧姆内阻值,为前一迭代中的并网电阻值、为当前迭代中的电容值,其中,为前一粒子个体i的位置对应的一组辨识参数,即通过第k-1次迭代辨识得到;
[0061]
所述锂电池工况数据集由锂电池动态应力测试dst工况数据构成。
[0062]
本发明辨识方法,能够在锂电池动态应力测试(dst)的电压电流数据下准确的进行分数阶等效电路模型的参数辨识,其有益效果包括:
[0063]
1.通过混沌映射初始化种群,使得初始种群分布更均匀,更加具有多样性,其中包含更多有价值的信息,从而在收敛到最优解时产生更高的成功率。
[0064]
2.引入了分数阶粒子速度更新律,提高粒子群种群的收敛速度。
[0065]
3.通过惯性权重自适应更新和对全局最优点附近进行混沌搜索来避免陷入局部最优。
附图说明
[0066]
图1是本发明方法中建立的锂电池等效电路模型的示意图。
[0067]
图2是本发明辨识方法的流程图。
[0068]
图3是本发明实施例中锂电池等效电路模型端电压同实际值的对比图。
[0069]
图4是本发明实施例中锂电池等效电路模型端电压误差的对比图。
具体实施方式
[0070]
下面结合附图对本发明的技术方案做进一步的详细说明。
[0071]
图2是本发明辨识方法的流程图。由图2可见,本发明提供了一种锂电池等效电路模型参数的辨识方法,所述辨识方法利用混沌自适应分数阶粒子群算法对锂电池等效电路模型的参数进行辨识,然后将得到的最优参数反馈到锂电池等效电路模型中,具体步骤如下:
[0072]
步骤1,获取锂电池等效电路模型的辨识参数
[0073]
图1为本发明中建立的锂电池等效电路模型的示意图,由图1可见,该锂电池等效电路模型为锂电池一阶rc等效电路模型,引入分数阶常相角(cpe)元件构成,具体的,所述锂电池电路等效模型的拓扑结构包括一个直流电源、一个欧姆内阻、一个并网电阻和一个分数阶常相位角元件cpe,所述并网电阻和分数阶常相位角元件cpe并联,并记为并联支路,所述欧姆内阻的一端与直流电源的正极相连,另一端与并联支路的一端相连,所述并联支路的另一端与直流电源的负极组成所述锂电池电路等效模型的输出端;
[0074]
将所述直流电源处的电压记为开路电压u
ocv
,所述欧姆内阻的阻值记为欧姆内阻值r0,所述并网电阻的阻值记为并网电阻值r1,所述锂电池等效电路模型的输出端的电压记为端电压u
t

[0075]
将所述分数阶常相位角元件cpe记为元件cpe,元件cpe的阻抗传递函数的表达式为:其中,z(s)为元件cpe的阻抗,c1为元件cpe的电容值,记为电容值c1,s为拉普拉斯算子,β为元件cpe的分数阶阶次,记为阶次β。
[0076]
所述锂电池等效电路模型的辨识参数为:欧姆内阻值r0,并网电阻值r1、电容值c1和阶次β,即在混沌自适应分数阶粒子群算法中,目标搜索空间的维数为4,将欧姆内阻值r0,并网电阻值r1、电容值c1和阶次β分别记为辨识参数1、辨识参数2、辨识参数3和辨识参数4。
[0077]
将所述辨识参数中的每一个辨识参数视为粒子个体的一个维度,所述粒子个体有粒子个体位置与粒子个体速度两个属性,粒子个体位置代表移动的方向,粒子个体速度代表移动的快慢,其中,粒子个体位置为辨识参数的值。
[0078]
步骤2,混沌自适应分数阶粒子群算法参数的设置
[0079]
设定目标搜索空间的维数为d,d=4;设一个种群由n个粒子组成;设定最大迭代次数m;设定粒子惯性权重ω;
[0080]
本实例中,取n=50,m=150。
[0081]
设定粒子个体的位置范围为{x_min,x_max}、粒子个体的速度范围{v_min,v_max}、粒子惯性权重范围为{ω_min,ω_max},其中,x_min为粒子个体位置的最小值,x_max为粒子个体位置的最大值,v_min为粒子个体速度的最小值,v_max为粒子个体速度的最大值,ω_min为粒子惯性权重最小值,ω_max为粒子惯性权重最大值;
[0082]
本实例中,取x_min={0.2,0.2,200,0.1}、x_max={0.9,2,1500,1}、v_min={0.01,0.1,100,0.01}、v_max={0.02,0.4,200,0.2}、ω_min=0.4、ω_max=0.9。
[0083]
将n个粒子个体中的任意一个记为粒子个体i,i为种群中任意一个粒子个体的序号,i=1,2,..n,将粒子个体i的位置向量记为个体位置xi,粒子个体i的速度向量记为个体速度vi,其表达式分别如下:
[0084]
xi={x
i1
,x
i2
,...,x
ij
...,x
id
}
[0085]vi
={v
i1
,v
i2
,...v
ij
...,v
id
}
[0086]
其中,个体位置xi为锂电池等效电路模型的一组辨识参数的解,x
ij
为粒子个体i第j维的位置,j=1,2,..d;个体速度vi为锂电池等效电路模型的一组参数的解在粒子搜索解空间中的速度,v
ij
为粒子个体i第j维的速度;
[0087]
将粒子个体i的适应度值记为个体适应度值fi,将粒子个体i搜索到的适应度值最小的位置记为个体最优位置p
best
,p
best
={p
i1
,p
i2
,...p
ij
...,p
id
},其中,p
ij
为粒子个体i的第j维搜索到的最优位置,将整个粒子群搜索到的适应度值最小的位置记为全局最优位置g
best
,g
best
={p
g1
,p
g2
,...,p
gj
...,p
gd
},p
gj
为第j维度整个粒子群搜索到的全局最优位置。
[0088]
步骤3,混沌自适应分数阶粒子群算法的初始化
[0089]
设迭代次数为k,k=1,2,..m;
[0090]
使用logistic混沌映射来初始化粒子群前4次迭代中的粒子个体位置和粒子个体速度;所述logistic混沌映射表达式如下:
[0091]zij
=2z
ij
(1-z
ij
)
[0092]
其中,z
ij
为粒子个体i第j维的0-1之间的随机数,z
ij
为粒子个体i第j维生成的混沌参数;
[0093]
与目标搜索空间的维数d对应,由logistic混沌映射表达式生成混沌序列{z
i1
,z
i2
,...z
ij
...,z
id
},且通过下式映射到粒子群的目标搜索空间:
[0094][0095]
式中,为第k次迭代中的粒子个体i第j维的位置,为第k次迭代中的粒子个体i第j维的速度,x
ij
_max为粒子个体i第j维位置的最大值,x
ij
_min为粒子个体i第j维位置的最小值,v
ij
_max为粒子个体i第j维速度的最大值,v
ij
_min为粒子个体i第j维速度的最小值,k=1,2,3,4;
[0096]
由此得到混沌映射粒子个体初始化位置和混沌映射粒子个体初始化速度随机初始化粒子个体i最优位置随机初始化粒子群全局最优位置k=1,2,3,4,其表达式分别为:
[0097][0098][0099][0100]
在粒子惯性权重范围{ω_min,ω_max}内随机初始化粒子惯性权重k=1,2,
3,4。
[0101]
步骤4,第5次迭代-第m次迭代
[0102]
将第k次迭代记为当前迭代,第k-1次迭代记为前一迭代,第k-2次迭代记为前二迭代,第k-3次迭代记为前三迭代,第k-4次迭代记为前四迭代。
[0103]
步骤4.1,个体位置和个体速度的更新
[0104]
通过分数阶速度更新方程进行个体位置和个体速度的更新,得到粒子个体i在当前迭代中的速度和粒子个体i在当前迭代中的位置所述分数阶速度更新方程的表达式如下:
[0105][0106]
其中,为粒子个体i在前一迭代中的速度,为粒子个体i在前二迭代中的速度,为粒子个体i在前三迭代中的速度,为粒子个体i在前四迭代中的速度,为粒子个体i在前一迭代中的位置,r1为加速度权重系数1,r2为加速度权重系数2,为粒子个体i在前一次迭代中的迭代惯性权重系数,α为分数阶阶次,取α=0.5;为粒子个体i在前一次迭代中的最优位置,记为前一个体最优位置为粒子群在前一迭代中的全局最优位置,记为前一全局最优位置
[0107]
步骤4.2,惯性权重的更新
[0108]
记粒子个体i在当前迭代中的惯性权重为当前惯性权重表达式如下:
[0109][0110]
其中,为粒子个体i在前一迭代中的适应度值,为前一迭代中种群所有粒子个体适应度值的平均值,为前一迭代中种群所有粒子个体适应度值中的最小值。
[0111]
步骤4.3,计算粒子个体i在当前迭代中的适应度值并记为当前个体最优位置将所述前一个体最优位置对应的适应度值记为前一个体最优适应度值前一全局最优位置对应的适应度值记为前一全局最优适应度值进行如下判断:
[0112]
比较和取适应度值最小者的位置为当前个体最优位置并将当前个体最优位置对应的适应度值记为当前个体最优适应度值
[0113]
比较和取适应度值最小者的位置为当前模拟全局最优位置并将当前模拟全局最优位置对应的适应度值记为当前模拟全局最优适应度值为
[0114]
步骤4.4,为避免陷入局部最优,对当前模拟全局最优位置进行混沌搜索,所述混沌搜索的表达式如下:
[0115][0116]
其中,τ为混沌权重,zr是由logistic混沌映射表达式生成的混沌序列,xr为混沌搜索产生的新的粒子位置;
[0117]
通过混沌搜索的表达式对当前迭代模拟全局最优位置进行10次混沌搜索,产生10组粒子个体的位置,计算10组粒子个体的适应度值,并与当前模拟全局最优适应度值为比较,取其中适应度值最小的位置为当前全局最优位置该当前全局最优位置对应一组辨识参数。
[0118]
步骤4.5,进行如下判断:
[0119]
若k=m,结束迭代,输出当前全局最优位置并记为全局最终位置g,该全局最终位置g对应的一组辨识参数即为锂电池等效电路模型的一组最优参数;
[0120]
若k<m,返回步骤4.1进行下一次迭代。
[0121]
在上述步骤中,步骤4中所述适应度值的计算式如下:
[0122][0123][0124][0125]
其中,t为锂电池工况数据集中数据点的个数,m为锂电池工况数据集中数据点的序号,m=1,2,...t,u
t
(m)为锂电池工况数据集第m个数据点对应的端电压;为将前一迭代中得到的锂电池等效电路模型辨识参数带入锂电池等效电路模型中计算出的端电压,u
ocv
(m)为锂电池工况数据集第m个数据点对应的开路电压,i(m)为锂电池工况数据集第m个数据点对应的电流,为前一迭代中的cpe元件的电压值,β
k-1
为前一迭代中元件cpe的分数阶阶次,为β
k-1
阶导数,为前一迭代中的欧姆内阻值,为前一迭代中的并网电阻值、为当前迭代中的电容值,其中,为前一粒子个体i的位置对应的一组辨识参数,即通过第k-1次迭代辨识得到。
[0126]
所述锂电池工况数据集由锂电池动态应力测试dst工况数据构成。
[0127]
为了佐证本发明的效果,进行了仿真。
[0128]
在仿真中使用的锂电池为韩国三星inr18650-20r动力锂电池,锂电池工况数据集由该锂电池动态应力测试dst工况数据构成。
[0129]
仿真的实验结果由表1、图3和图4所示,以端电压的均方根误差(rmse)作为验证辨识结果精度的指标。rmse越小,辨识精度越高。其中pso为传统粒子群优化算法,cafpso为本发明中的混沌自适应分数阶粒子群优化算法。
[0130]
表1
[0131] r0r1c1βrmsepso0.753ω0.288ω702.4f0.5780.041vcafpso0.725ω1.865ω983.4f0.9520.009v
[0132]
表1为给出了两种算法中辨识参数的值和rmse。从表中1可以看出cafpso算法的辨识精度要优于传统pso算法。
[0133]
图3是本发明实施例中等效电路模型端电压同实际值的对比图,该对比图给出了pso算法,本发明cafpso算法和实际测量值中端电压的曲线,很明显本发明cafpso算法的结果更加接近于测量值。
[0134]
图4是本发明实施例中锂电池等效电路模型端电压误差ute的对比图,该对比图给出了pso算法,本发明cafpso算法中锂电池等效电路模型端电压误差ute的对比,很明显,本发明建立的锂电池等效模型端电压误差ute更小更加接近真实值。
再多了解一些

本文用于创业者技术爱好者查询,仅供学习研究,如用于商业用途,请联系技术所有人。

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

相关文献