2.845

2023影响因子

(CJCR)

  • 中文核心
  • EI
  • 中国科技核心
  • Scopus
  • CSCD
  • 英国科学文摘

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于相关向量机的高速列车牵引系统剩余寿命预测

王秀丽 姜斌 陆宁云

王秀丽, 姜斌, 陆宁云. 基于相关向量机的高速列车牵引系统剩余寿命预测. 自动化学报, 2019, 45(12): 2303−2311 doi: 10.16383/j.aas.c190204
引用本文: 王秀丽, 姜斌, 陆宁云. 基于相关向量机的高速列车牵引系统剩余寿命预测. 自动化学报, 2019, 45(12): 2303−2311 doi: 10.16383/j.aas.c190204
Wang Xiu-Li, Jiang Bin, Lu Ning-Yun. Relevance vector machine based remaining useful life prediction for traction systems of high-speed trains. Acta Automatica Sinica, 2019, 45(12): 2303−2311 doi: 10.16383/j.aas.c190204
Citation: Wang Xiu-Li, Jiang Bin, Lu Ning-Yun. Relevance vector machine based remaining useful life prediction for traction systems of high-speed trains. Acta Automatica Sinica, 2019, 45(12): 2303−2311 doi: 10.16383/j.aas.c190204

基于相关向量机的高速列车牵引系统剩余寿命预测

doi: 10.16383/j.aas.c190204
基金项目: 国家自然科学基金(61490703, 61873122, 61922042), 江苏高校优势学科建设工程资助项目, 南京航空航天大学博士生短期访学项目(180401DF03)资助
详细信息
    作者简介:

    王秀丽:南京航空航天大学自动化学院博士研究生. 主要研究方向为基于数据驱动的故障预测及其应用. E-mail: xiuliwang@nuaa.edu.cn

    姜斌:南京航空航天大学自动化学院教授. 主要研究方向为智能故障诊断与容错控制及其应用. 本文通信作者. E-mail: binjiang@nuaa.edu.cn

    陆宁云:南京航空航天大学自动化学院教授. 主要研究方向为基于数据驱动的故障诊断与预测及其应用. E-mail: luningyun@nuaa.edu.cn

Relevance Vector Machine Based Remaining Useful Life Prediction for Traction Systems of High-speed Trains

Funds: Supported by National Natural Science Foundation of China (61490703, 61873122, 61922042), Priority Academic Program Development of Jiangsu Higher Education Institutions, and Doctoral Student Short-term Visit Project of Nanjing University of Aeronautics and Astronautics (180401DF03)
  • 摘要: 高速列车牵引系统在运行过程中总是受到诸多不确定因素的影响, 例如, 由于列车的负载、运行环境及元器件的老化引起的不确定性, 不确定因素不可避免地影响牵引系统剩余寿命的预测精度. 为了提高不确定情景下剩余寿命预测的准确性, 本文首先采用改进的相关向量机(Relevance vector machine, RVM)方法, 建立鲁棒性能良好的多步回归模型, 由于t分布比常用的高斯分布更具有鲁棒性, 通过权重和随机误差服从t分布而非高斯分布, 改进了相关向量机回归模型, 随后将超参数的先验一并融入似然函数, 通过最大化似然函数估计未知的超参数, 此外, 利用首达时间方法从概率角度对剩余寿命进行了预测, 最后通过牵引系统中电容器退化的案例, 与传统的相关向量机方法、自回归方法和支持向量机方法进行对比, 验证了所提算法的有效性.
  • 牵引系统为高速列车的运行提供主要动力, 正如心脏通过血管为身体提供氧气, 因此, 牵引系统的剩余寿命(Remaining useful life, RUL)决定着高速列车的剩余寿命[1]. 事实上, 牵引系统的元器件无时无刻不因受到老化磨损等影响而退化, 这类退化的特征在退化初期是极其微小的[2-3], 以致很难察觉, 但是, 如若放任不管, 各类故障或失效会相继发生, 更有甚者, 引发重大灾难性事故, 造成人员伤亡及财产损失. 例如, 由于调度问题, $ 2011 $年甬台温铁路列车追尾事故, 造成直接经济损失约$ 1.94 $亿元及$ 40 $人员伤亡; $ 2018 $年发生的$1\cdot25 $G281次列车火灾事故, 因电气设备故障引起整节车厢着火; 1998年, 发生于德国的ICE列车脱轨事故至今仍触目惊心, 共造成$ 101 $人死亡, $ 88 $人受伤. 除此之外, 故障引发的高速列车紧急制动、列车晚点等比比皆是, 鉴于此, 为保障高速列车运行的安全性及可靠性, 提前预知牵引系统的剩余寿命显得尤为重要[4].

    牵引系统剩余寿命预测旨在利用历史及当前退化趋势数据推测当前时刻到停机时刻的有效运行时间. 近年来, 基于数据驱动的剩余寿命预测受到国内外专家学者的广泛关注, 基于数据驱动的算法可进一步分为: 基于状况的方法、基于统计过程的方法及基于机器学习的方法. 基于状况的剩余寿命预测方法利用观测特征推断隐藏的退化状况和潜在的退化趋势, 建立观测模型和系统模型, 随后估计剩余寿命, 其中, 隐Markov模型及粒子滤波是两种具有代表性的方法[5-6]. 基于统计的剩余寿命预测方法通过建立统计过程模型从概率角度对剩余寿命进行预测, 常用的过程模型有: Gamma模型、Wiener模型及逆Gaussian模型. 其中, Gamma过程模型主要适用于退化趋势单向且单调的退化过程[7]; Wiener过程模型则适合随着时间的推移退化过程随高斯噪声双向变化的退化过程[8]; 而逆Gaussian过程相较于其他统计方法更加灵活, 适用于较多的统计过程[9]. 机器学习方法近年来发展迅速, 自回归(Autoregressive, AR)方法是一种较为经典的基于时间序列的机器学习预测方法, 其利用同一变量的历史数据以预测该变量未来的趋势[10], 自回归方法仅适用于受自身历史因素影响较大的系统, 为了弥补自回归方法的不足, 向量自回归模型得以蓬勃发展[11-12]. 支持向量机(Support vector machine, SVM)方法也是一种典型的机器学习方法, 并已广泛用于故障预测中[13-15], 相关向量机(Relevance vector machine, RVM)方法是一种形式上类似于支持向量机的稀疏贝叶斯学习方法, 它不仅继承了支持向量机的优点, 还具有良好的稀疏性及泛化能力, 并且不再受Mercer定理的限制.

    相关向量机于2001年由Tipping提出[16], 随后广泛应用于故障及剩余寿命预测中. 相关向量机模型可用于预测后续车头时距[17], 主要利用蕴含驶过站点的车头时距、运行时间及乘客客流量的时间序列建立相关向量机模型, 相关向量机还应用于预测质子交换膜燃料的剩余寿命中[18], 与粒子滤波结合也是一种常见的方法[19]. 大部分文献研究相关向量机时均推导出预测值服从高斯分布[20-22], 但是, 高速列车在运行过程中, 总是受到诸多不确定性因素的影响, 例如, 经过不同的地段受到的阻力不同, 经过山地或者丘陵时爬坡较多, 并且不同站点载客量也会随机变化[23]. 此时若不寻找鲁棒性能好的模型, 预测精度就会大大降低, 而高斯模型的鲁棒性较差, 因此, 怎样利用受不确定影响的数据得到鲁棒性较好的预测模型成为亟待解决的问题.

    不确定因素影响相关向量机中权重及随机误差的方差项, 而一般涉及相关向量机的算法总是假设方差服从均匀分布[16]. 为了使得预测算法对不确定因素具有良好的鲁棒性, 本文考虑相关向量机中权重及随机误差的精度 (方差的逆)服从$ \Gamma $分布, 进而权重及随机误差便服从$ t $分布, $ t $分布是一种比高斯分布鲁棒性强的分布形式, 这意味着$ t $分布对于不确定因素造成的一些异常值的灵敏度远低于高斯分布.

    本文的主要安排如下: 第1节对牵引系统的构造进行简要描述, 并分析中间直流环节支撑电容器性能退化的退化特征; 第2节利用相关向量机对退化过程进行建模及参数优化; 在第3节中, 采用首达时间从概率角度对剩余寿命进行预测; 第4节, 用实验数据对所提出的方法进行仿真验证, 该仿真实验将本文提出的方法与传统相关向量机、自回归方法及支持向量机方法进行对比; 第5节进行概括性总结及展望.

    本节首先阐述牵引系统运行机理, 并针对支撑电容器性能退化阐述支撑电容器剩余寿命预测的重要性, 随后针对支撑电容器性能退化提取特征量以表征其退化趋势.

    牵引系统是高速列车的核心部件, 其通过一系列设备将电能转换为机械能, 为列车的运行提供驱动力. 首先, 高压电通过受电弓给变压器提供单相交流电, 变压器降低电压后, 单相交流电输送给变流器(包括整流器、中间直流环节和逆变器), 其中单相交流电通过整流器转换为直流电, 直流电再通过中间直流环节输送到逆变器, 然后, 逆变器向电动机提供三相交流电, 最后, 电动机的扭矩和速度传递到车轮通过变速箱驱动列车.

    高速列车的牵引传动方式是“交 − 直 − 交”, 也就是从“交流电”到“直流电”, 再到“交流电”. 为什么中间端要增加直流, 而不是直接“交 − 交”呢? 这主要是技术上的一个问题, 虽然转换器的输入和输出都是交流电, 但交流电不能直接从整流器流向逆变器, 因为从整流器直接到逆变器的“交流 − 交流”之间功率因数太低, 谐波太大, 对电气设备会造成一定的危害. 所以, 在“交流 − 交流”之间建立一个“缓冲带”, 即中间直流环节.

    支撑电容器是中间直流环节中用于维持牵引系统中电压稳定的重要元器件. 当牵引系统的负载突然加大时, 整流器不能立即提供大量能量, 引起电动机转速下降, 此时, 支撑电容器在直流环节中起到能量补充的作用. 当牵引系统突然变为再生制动时, 整流器不能瞬间吸收大量能量, 导致中间直流电压上升, 此时, 支撑电容器起到吸收能量的作用, 从而缓解电压上升. 在列车实际运行中, 中间直流环节是变流器的薄弱环节, 尤其是支撑电容器易受恶劣环境的侵害, 且时刻受到牵引系统内部大电流的冲击, 有更高的机率发生性能退化, 另一方面, 支撑电容器造价昂贵. 因此, 基于电容器性能退化数据预测电容器的剩余寿命显得尤为重要.

    中间直流环节支撑电容器主要由上支撑电容器及下支撑电容器组成, 如图1所示, 支撑电容器在性能退化过程中会导致上下电压波动, 若其中一个出现性能退化, 该支撑电容器会将其储存的电压全部释放直至失效, 此时对应的端电压将不断减小直至中间直流环节停机.

    图 1  中间直流环节结构简图
    Fig. 1  The structure diagram of intermediate DC link

    在支撑电容器性能退化过程中, 上下端电压的幅值会随着时间的推移而增减, 端电压包含着丰富的支撑电容器性能退化信息, 因此, 本文选取上下端电压作为特征量. 支撑电容器的稳定性极易受到诸多不确定因素的影响: 不同牵引系统连续工作的环境温度会使得系统内部温度不同, 内部温度的高低直接影响着支撑电容器的剩余寿命; 不同电容器间的容许耐压有所不同, 当施加过电压时, 电容器内部短路时刻有所差异; 不同电容器间受到的高次谐波的冲击也不同. 可针对多批次电容器多次采集上下端电压来提取特征量, 提取的特征量包含两个维度: 时间维度及批次维度. 时间维度的数据通过时间的变化反映退化趋势, 不同批次间的数据蕴含不确定性因素对特征量的影响.

    不确定因素使得不同批次间的退化数据存在一定的差异性, 每个批次数据的退化趋势随着时间的变化退化率逐渐变化, 基于退化数据不同批次间的差异性及随时间的动态特性, 首先利用相关向量机建立动态模型, 由于不确定性的影响, 对模型中的超参数引入$ \Gamma $分布, 然后基于超参数的似然函数及先验概率对超参数进行优化.

    已知批次 − 时间维度数据

    $ X_{M\times N} = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1N} \\ x_{21}& x_{22} & \cdots & x_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ x_{M1} & x_{M2} & \cdots & x_{MN} \end{bmatrix}\in {\bf{R}}^{M\times N}$

    其中, $ M $表示批次总数; $ N $代表随时间的采样总数. 基于批次对矩阵$ X_{M\times N} $的每行求均值得到 $[x_1,$$x_2,\cdots, x_N]\in{\bf{R}}^{\textit{1}\times N} ,$其中, $ x_i = \dfrac{x_{1i}+\cdots+x_{Mi}}{M},\; i = $$1, 2, \cdots, N $. 接下来, 运用时刻$ \{t_1, t_2, \cdots, t_N\} $的退化数据$ X_{1:N} = $$ \{x_1, x_2, \cdots, x_N\} $, 基于相关向量机及$t_n,\;1\leq n< $$ N $ 时刻的数据$ x_n $建立l -$\!\,$步动态回归模型

    $ x_{n+l} = {{ \phi}}(x_n){{ \omega}}+{\rm e} $

    (1)

    其中, $ x_{n+l} $$ x_n $$l $-$\!\,$步预测值, 且$ 1< n+l\leq N $; $ {{ \phi}}(x_i) \!=\! [1, {\cal K}(x_i,x_{n-l+1}), \cdots, {\cal K}(x_i,x_n)]\in{\bf{R}}^{1\times (l+1)} $为基向量, $ {\cal K}(x_i,x_j),\;j = n-l+1, \cdots, n $是核函数; ${{ \omega}} = [\omega_0, \omega_1, \cdots, \omega_l]^{\rm{T}}\in{\bf{R}}^{(l+1)\times 1} $是对应基向量的权重向量; 随机误差$ e $关于批次服从均值为 $ 0 $、方差为$ \beta^{-1} $(为了后续计算方便, 采用精度$ \beta $)的高斯分布. 因此, 第$ l $-$\!\,$步预测值$ x_{n+l} $也服从高斯分布, 并且均值为$ y(x_n) = {{\phi}}(x_n){{\omega }}$, 方差为$ \beta^{-1} $.

    预测值$ x_{n+l} $关于$ {{\omega }}$$ \beta $的概率密度函数可进一步表示为

    $ p(x_{n+l}|{{\omega}}, \beta) = \left(\dfrac{\beta}{2{{\text{π}}}} \right)^\frac{1}{2}{\rm e}^{\frac{-\beta}{2}(x_{n+l}-\phi \omega )^2} $

    (2)

    我们的目标是寻找合适的参数$ {{\omega }}$$ \beta $使得模型 (1)成立. 如果直接对参数${{ \omega }}$$ \beta $进行估计, 会导致过拟合现象, 尤其是在不确定因素影响下. 为了避免过拟合现象的发生, 假设权重$ {{\omega }}$服从均值为$ 0 $的先验高斯分布, 则其概率密度函数为

    $ \begin{split} p({{\omega}} |{{\alpha}} ) =& \prod\limits_{i = 0}^l {\sqrt {\dfrac{{{\alpha _i}}}{{2{\text{π}} }}} } \exp \left( { - \dfrac{{{{\omega}} _i^2{\alpha _i}}}{2}} \right)=\\ & {(2{\text{π}} )^{ - \frac{{l + 1}}{2}}}|A{|^{\frac{1}{2}}}{{\rm e}^{ - \frac{1}{2}{{{\omega}} ^{\rm{T}}}A{{\omega}} }} \end{split} $

    (3)

    其中, $ {{\alpha}} = [\alpha_0, \alpha_1, \cdots, \alpha_l]^{\rm{T}}{\in{\bf{R}}^{(l+1)\times 1}} ;$矩阵 $ A = $${\rm{diag}} \{\alpha_0, \alpha_1, \cdots, \alpha_l\}\in{\bf{R}}^{(l+1)\times (l+1)} $是权重的精度矩阵. 值得一提的是$ {{\omega }}$中的每个分量是独立同分布的, 则每个超参数$ \alpha_i $分别对应着相应的权重分量$ \omega_i $. 在相关向量机回归模型中, 如果超参数中的某些值$ \alpha_{i0} $, $ \alpha_{i1}\cdots $, $ \alpha_{ik} $$ (0\leq i0, \cdots, ik\leq l) $趋于$ 0 $, 那么, 相对应的权重$ \omega_{i0}, \omega_{i1}, \cdots, \omega_{ik} $$ (0\leq i0, \cdots, ik\leq l) $也趋于0. 事实上, 非零权重对应的样本称之为相关向量.

    涉及相关向量机的诸多文献总是假设超参数$ \alpha $$ \beta $服从均匀分布, 但是不确定因素导致退化数据存在一定的异常点, 若超参数服从均匀分布, 则可得到权重及随机误差服从高斯分布, 事实上, 高斯回归模型对异常点格外敏感. 若对超参数$ \alpha $$ \beta $引入超先验分布 — $ \Gamma $分布, 此时, 权重及随机误差则服从$ t $分布, $ t $分布较高斯分布具有更强的鲁棒性. 为了寻找鲁棒性强的回归模型, 首先给出超参数$ \alpha $$ \beta $的先验概率密度函数

    $ p(\alpha ) = \prod\limits_{i = 0}^l \Gamma ({\alpha _i}|a,b) $

    (4)

    $ p(\beta) = \Gamma(\beta|c,d) \hspace{20pt}$

    (5)

    其中, $ \Gamma(\alpha|a,b) \!=\! \Gamma(a)^{-1}b^a\alpha^{a-1}e^{-b\alpha} , $$ \Gamma(a) \!=\! \int_{0}^{\infty}t^{a-1} ×\!$$e^{-t}{\rm{d}}t $. 此时, 权重分量$ \omega_i $及随机误差$ e $的先验概率分布函数为

    $ \begin{split} p({\omega _i}) =& \int p ({\omega _i}|{\alpha _i})p({\alpha _i}){\rm{d}}{\alpha _i}=\\ & \dfrac{{{b^a}\Gamma \left(a + \dfrac{1}{2}\right)}}{{{{(2{\text{π}} )}^{\frac{1}{2}}}\Gamma (a)}}{\left(b + \dfrac{{\omega _i^2}}{2}\right)^{ - \left(a + \frac{1}{2}\right)}}=\\ & \dfrac{{\Gamma \left(\dfrac{\nu }{2} + \dfrac{1}{2}\right)}}{{\Gamma \left(\dfrac{\nu }{2}\right)}}{\left(\dfrac{\lambda }{{{\text{π}} \nu }}\right)^{\frac{1}{2}}}{\left(1 + \dfrac{{\lambda \omega _i^2}}{\nu }\right)^{ - \left(\frac{\nu }{2} + \frac{1}{2}\right)}}=\\ & t({\omega _i}) \end{split} $

    (6)

    $ \begin{split} p(e) =& \int p (e|\beta )p(\beta ){\rm{d}}\beta= \\ & \dfrac{{{d^c}\Gamma \left(c + \dfrac{1}{2}\right)}}{{{{(2{\text{π}} )}^{\frac{1}{2}}}\Gamma (c)}}{\left(d + \dfrac{{{e^2}}}{2}\right)^{ - \left(c + \frac{1}{2}\right)}}=\\ & \dfrac{{\Gamma \left(\dfrac{{\nu '}}{2} + \dfrac{1}{2}\right)}}{{\Gamma \left(\dfrac{{\nu '}}{2}\right)}}{\left(\dfrac{{\lambda '}}{{{\text{π}} \nu '}}\right)^{\frac{1}{2}}}{\left(1 + \dfrac{{\lambda '{e^2}}}{{\nu '}}\right)^{ - \left(\frac{{\nu '}}{2} + \frac{1}{2}\right)}}=\\ & t(e) \end{split}$

    (7)

    其中, 参数$ \lambda = {a/b} $$ \nu = 2a $分别为分布$ t(\omega_i) $的精度和自由度, 参数 $ \lambda' = {c/d} $$ \nu' = 2c $分别为分布$ t(e) $的精度和自由度.

    通过式(6)和式(7), 观察到权重$ \omega_i,\;i = 1,2,\cdots, $$l$及随机误差$ e $均服从$ t $分布, 事实上, $ t $分布通常有着比高斯分布更长的“尾巴” , 而回归模型本身不具有鲁棒性, 通过让回归模型系数基于长尾的概率分布可得到一个更加鲁棒的模型, 后续将通过仿真实例进一步验证.

    为了得到回归模型(1), 需要对超参数$ \alpha_i,\;i = 0, $$ 1, \cdots, l$$ \beta $进行估计. 根据式(2)和式(3)及卷积公式, 预测值$ x_{n+l} $关于超参数$ \beta $$ \alpha$的条件概率密度函数为

    $ \begin{split}& p({x_{n + l}}|\alpha ,\beta ) =\int_{ - \infty }^{ + \infty } p ({x_{n + l}}|\omega ,\beta )p(\omega |\alpha ){\rm{d}}\omega= \\ &\qquad{\left( {\frac{\beta }{{2{\text{π}} }}} \right)^{\frac{{l + 1}}{2}}}{\left( {\frac{1}{\beta } + \phi {A^{ - 1}}{\phi ^{\rm{T}}}} \right)^{ - \frac{1}{2}}}{{\rm e}^{ - \frac{1}{2}(\frac{1}{\beta } + \phi {A^{ - 1}}\phi )x_{n + l}^2}}\end{split} $

    (8)

    在式(8)中, 超参数$ \alpha_i $$ \beta $的先验概率密度函数由式(4)和式(5)给出, 接下来, 通过最大化边缘似然函数及$ \alpha_i $, $ \beta $的先验的乘积对超参数$ \alpha_i $, $ \beta $进行优化. 以便简化计算量, 利用对数函数的单调性, 对目标似然函数求对数, 则目标似然函数的对数函数为

    $ \begin{split} {\cal L} =\;& - \frac{1}{2}\Big[\ln |\frac{1}{\beta }I + \phi {A^{ - 1}}\phi^{\rm T} | + {x_{n + l}}\Big(\frac{1}{\beta }I+ \\ & \phi {A^{ - 1}}\phi^{\rm T} \Big){^{ - 1}}{x_{n + l}}\Big] + \sum\limits_{i = 0}^l {(a - \ln {\alpha _i} - b{\alpha _i})}+ \\& c\ln \beta - d\beta \end{split} $

    (9)

    然后, $ {\cal L} $关于$ \ln \alpha_i $$ \ln\beta $求一阶偏导数

    $ \frac{\partial{\cal L}}{\partial \ln \alpha_i} = \frac{1}{2}\left[1-\alpha_i\left(\mu_i^2+\Sigma_{ii}\right)\right]+a-b\alpha_i \hspace{22pt}$

    (10)

    $ \begin{split} \dfrac{{\partial L}}{{\partial \ln \beta }} =\;& \dfrac{1}{2}\left[\dfrac{n}{\beta } - {\left({x_{n + l}} - \phi \mu \right)^2} - {\rm{tr}}\left(\Sigma {\phi ^{\rm{T}}}\phi \right)\right]+\\& c - d\beta \end{split}$

    (11)

    其中,

    $ \Sigma = \left(\beta\phi^{\rm{T}}\phi+A\right)^{-1} $

    (12)

    $ \mu = \beta\Sigma\phi^{\rm{T}}x_{n+l} \hspace{20pt}$

    (13)

    分别为权重$ {{\omega}} $的后验协方差及均值; $ \mu_i $是第$ i $个权重的后验均值; $ \Sigma_{ii} $$ \Sigma $的第$ i $项对角元素.

    $ \dfrac{\partial{\cal L}}{\partial \ln \alpha_i} = 0 $$ \dfrac{\partial{\cal L}}{\partial \ln \beta} = 0 $, 则超参数可优化为

    $ \alpha_i = \frac{\gamma_i+2a}{\mu_i^2+2b} \hspace{41pt} $

    (14)

    $ \beta = \frac{n-\Sigma_{ii}\mu_i+2c}{(x_{n+l}-\phi \mu)^2+2d} $

    (15)

    其中, $ \gamma_i = 1-\alpha_i\Sigma_{ii} $.

    式(14)及式(15)中的超参数$ \alpha_i $$ \beta $与变量$ \gamma_i $, $ \mu_i $$ \Sigma_{ii} $有关. 通过式(12)和式(13)可以看出, 变量$ \gamma_i $, $ \mu_i $$ \Sigma_{i} $又反包含未知超参数$ \alpha_i $$ \beta $, 所以, 式(12)和式(13)仅是超参数$ \alpha_i $$ \beta $的隐式表达形式. 为了得到超参数$ \alpha_i $$ \beta $的具体数值, 首先根据经验[16], 令$ a = b = c = d = 10^{-4} $, 再根据式(4)和式(5), 求出$ \alpha_i $$ \beta $的先验概率的均值, 然后代入式(12)和式(13), 得到$ \gamma_i $, $ \mu_i $$ \Sigma_{i} $的值, 接着代入式(14)和式(15)得到超参数$ \alpha_i $$ \beta $, 依次迭代直至$ \alpha_i $$ \beta $收敛.

    按照常规思维, 超参数$ a, b, c, d $需要进行优化而非根据经验得到, 但是, 若想通过式(4)和式(5)优化$ a, b, c, d $, 需对$ \alpha $$ \beta $进行采样, 事实上, $ \alpha $$ \beta $分别是权重$ {{\omega}} $和随机误差$ e $的超参数, 实际操作中无法得到$ \alpha $$ \beta $的采样数据, 故而无法对$ a, b, c, d $利用极大似然估计等算法进行优化, 只能通过经验进行赋值. 我们的目标是求得$ \alpha $$ \beta $, 对$ a, b, c, d $赋值相当于对$ \alpha $$ \beta $的先验概率密度进行赋值, 由于式(14)和式(15)的多次迭代, 可以得到渐近收敛的$ \alpha $$ \beta $的估计值, 故$ \alpha $$ \beta $先验概率密度函数不会影响$ \alpha $$ \beta $的估计.

    通过改进的相关向量机建立退化模型后, 接下来, 将利用首达时间定义寿命, 随后, 鉴于不确定性对剩余寿命的影响, 从概率角度出发预测出剩余寿命.

    已知测量数据$ x_1, x_2, \cdots, x_n $, 在时刻$ t_n $的寿命$ L $可通过首达时间定义为

    $ L = \inf\{t_l: x({t_n}+t_l)\equiv x_{{n}+l}\in\Omega|X_{1:{n}}\} $

    (16)

    其中, $ \inf\{\cdot\} $表示变量的下确界; $ t_l $$ l $步之后的时间; $ X_{1:n} $代表从时刻$ t_1 $$ t_n $的历史测量数据; $ x_{n+l} $表示时刻$ t_{n}+t_l $的预测值, 可通过式(1)得到; $ \Omega $是预设的失效阈值, 通常根据实际应用经验设定, 本文中的失效阈值由中车株洲所提供. 首达时间的提出是从实际应用中的安全性考虑, 当退化变量一旦超过失效阈值边界$ D\in\Omega $就认定为失效, 接下来我们将估计首次到达失效阈值边界$ D $的时间$ L $的概率密度.

    基于测量数据$ X_{1:n} = \{x_1, x_2, \cdots, x_n\} $, 寿命$ L $的条件累积分布为

    $ \begin{split} P\left\{ L \le t|{X_{1:n}}\right\} =\;& P\left\{ {x_{n + l}} \ge D|{X_{1:n}}\right\}= \\ & 1 - P\left\{ {x_{n + l}} \le D|{X_{1:n}}\right\}= \\ & 1 - P\left\{ z \le \dfrac{{D - {{\tilde \mu }_{n + l}}}}{{\sqrt {\tilde \sigma _{n + l}^2} }}\right\}= \\ & P\left\{ z \ge \dfrac{{D - {{\tilde \mu }_{n + l}}}}{{\sqrt {\tilde \sigma _{n + l}^2} }}\right\}= \\ & \Phi ({g_{n + l}}) \end{split} $

    (17)

    其中, $ z $为服从标准分布的随机变量; $ \Phi(\cdot) $是标准累积正态分布; $ g_{n+l} = {\tilde{\mu}_{n+l}-D}/{\sqrt{\tilde{\sigma}^2_{n+l}}} $, $ \tilde{\mu}_{n+l} $表示预测值$ x_{n+l} $的均值, $ \tilde{\sigma}_{n+l}^2 $$ x_{n+l} $的方差, 它们均可以通过式(8)得到. 事实上, 总有$ L\geq 0 $, 所以引入寿命$ L $的截尾条件累积分布

    $ \begin{split} P\{ L \le t|{X_{1:n}},L \ge 0\} =& \dfrac{{P\{ 0 \le L \le t|{X_{1:n}}\} }}{{P\{ L \ge 0|{X_{1:n}}\} }}=\\ & \dfrac{{\Phi ({g_{n + l}}) - \Phi ({g_n})}}{{1 - \Phi ({g_n})}} \end{split} $

    (18)

    得到$ L $的截尾条件概率后, 通过对$P\{L\leq t|X_{1:n}, L\geq $$ 0\} $求导数, 可进一步求得$ L $的截尾条件概率密度.

    $ p_{\{L|X_{1:n}, L\geq 0\}} = \frac{\phi(g_{n+l})\Delta g(n+l)}{1-\Phi(g_n)} $

    (19)

    其中, $ \phi(\cdot) $表示标准正态分布的概率密度. 随后, 可得到剩余寿命的期望

    $ {\rm{E}}_{\{L|X_{1:n}, L\geq 0\}} = \frac{\int_{0}^{+\infty}{l}p_{\{l|X_{1:n}, l\geq 0\}}{\rm{d}}l}{\int_{0}^{+\infty}p_{\{l|X_{1:n}, l\geq 0\}}{\rm{d}}l} $

    (20)

    高速列车牵引系统作为复杂系统, 在负载工况、振动及元器件老化等诸多不确定性因素的影响下, 不同批次间的电容器性能退化存在一定的差异性. 为充分考虑不确定因素对RUL预测的影响, 可选定置信水平利用Neyman置信区间方法给出RUL预测值的置信上下限. 具体实现如下:

    1)利用方差与均值间的关系$ \sigma_L^2 = E_{L^2}-$$ (E_L)^2 $及式(20)求出剩余寿命的标准差$ \sigma_L $;

    2)通过剩余寿命分布的均值和标准差, 构造枢轴量$ \dfrac{E_L-\mu_L}{\sigma_L/\sqrt{M}} $, 其中, $ \mu_L $是剩余分布的已知真实均值; $ M $是批次数据;

    3)给定置信水平$ 1-\alpha $, 利用$ P\Big\{-z_{\alpha/2} < $$\dfrac{E_L-\mu_L} {\sigma_L/\sqrt{M}} < z_{\alpha/2}\Big\} = 1-\alpha $, 求得置信区间为$\Big(E_L- $$ \dfrac{\sigma_L}{\sqrt{M}}z_{\alpha/2}, E_L+\dfrac{\sigma_L}{\sqrt{M}}z_{\alpha/2}\Big). $

    基于中车株洲电力机车研究所开发的CRH2A型高速列车dSPACE平台, 中南大学模拟了直流环节支撑电容器性能退化过程[24-25], 利用此平台采集到的数据, 首先对退化特征进行了提取, 然后针对建立的模型验证其合理性, 随后, 进行了剩余寿命预测验证, 并与传统的相关向量机方法、自回归方法和支持向量机方法进行了对比, 最后通过图形解释了$ t $分布比高斯分布具有较好的鲁棒性.

    我们在株洲研究所的dSPACE平台采集了$ 13 $组支撑电容器退化数据, 进行支撑电容器的剩余寿命预测验证, 退化数据包含中间直流环节上端电压及下端电压, 图2给出上支撑电容器出现性能退化时的端电压从平稳过程到失效停机的退化过程.

    图 2  上下端电压从平稳过程到退化过程直至停机的演变趋势
    Fig. 2  The voltages evolution from the stationary process to the degradation process

    图2可以看出, 中间直流环节支撑电容器在正常工作阶段上下端电压均稳定在2 700 V左右, 当上支撑电容器发生故障时, 上端电压逐渐降低, 下端电压逐渐上升, 直到上下端电压超出设定失效阈值触发牵引系统停机. 由于$ 13 $组退化支撑电容器的停机时间不同, 为了得到上下端电压的失效阈值, 对$ 13 $组数据求触发停机的端电压的均值, 经过计算, 上下端电压的失效阈值分别为2 165 V和3 205 V.

    为了预测支撑电容器的剩余寿命, 首先需要选择特征量以表征退化随时间演变的趋势. 在支撑电容器退化过程中, 上下端电压的幅度随着运行时间呈下降和上升的趋势, 因此, 选取上下端电压作为支撑电容器性能退化过程的特征信号, 然后对$ 13 $组支撑电容器求平均振幅以表征任意高速列车牵引系统支撑电容器的性能退化趋势. 图3给出了退化过程中的上下端电压的平均幅值.

    图 3  上下端电压在支撑电容器退化过程中的平均幅值
    Fig. 3  The mean amplitude of the upper and lower voltages during the capacitors' degradation

    为了建立模型(1), 选取$ 800\sim1\;700\;{\rm{min}} $的平均幅值作为输入数据, $ l $-步的预测步长设为$20 \;{\rm{min}}$, 则目标数据是$820\sim$$1\;720\;{\rm{min}}$的端电压, 其服从的分布可通过$ 13 $组批次数据得到, 随后利用式(14)和式(15)来训练超参数, 使训练输出值与目标数据匹配良好. 为了验证建立的模型的合理性, 将预测的输出电压与实际值进行比较, 如图4所示.

    图 4  建模过程中的预测值与真实值进行比较
    Fig. 4  The comparison between predicted values and the actual values in the modeling process

    图4中, 粗实线表示通过$ 13 $组端电压获得的实际平均幅值, 虚线是以$ 13 $组数据为输出训练超参数后, 代入模型(1)得到的预测平均幅值, 由图4可知, 预测值与相应的目标值能够很好匹配.

    首先利用训练好的模型 (1)对$1\;701\sim2\;400 \;{\rm{min}}$端电压进行趋势预测, 而后利用第$ 3 $节提到的首达时间方法对剩余寿命进行预测. 为了证明本文提出的方法的优越性, 利用自回归模型 (AR)、传统相关向量机方法和支持向量机 (SVM)方法与本文提出的方法进行比较.

    自回归模型是一种非概率线性预测方法, 在本文中, 首先针对退化趋势建立模型: $x_{n+l} =c+ $$ \displaystyle\sum\nolimits_{i = 1}^{p}\phi_i x_{n-i}+e, $其中, $ c $是常数项; $ p $是阶数; $ \phi_i $是方程系数; $ e $是均值为0、标准差为$ \sigma $的随机误差值. 然后利用首达时间对剩余寿命进行预测. 支持向量机模型的退化趋势方程为$ x_{n+l} = \phi(x_n)\omega+e, $其中, $ \phi(\cdot) $是基函数; $ \omega $是权重常数; $ e $是均值为0、标准差为$ \sigma $的随机误差值. 传统相关向量机模型与支持向量机模型形式类似, 但是相关向量机模型建立在稀疏贝叶斯理论框架下, 权重$ \omega $满足式(3), 而本文中改进的相关向量机方法不仅满足式(3), 还满足式(4)和式(5). 接下来针对下端电压趋势模型给出模型中的参数, 为了更有效地对比各算法, 模型中的$ \sigma $均设为$ 0.3, $自回归模型中的阶数$ p $设为2, 系数分别为$ c = 0 ,$$ \phi_1 = 1.022 $$ \phi_2 = 0.994. $支持向量机、相关向量机及改进的相关向量机中的参数选择如表1.

    表 1  退化趋势模型中的参数取值
    Table 1  The parameter values in the degradation model
    SVMRVM改进RVM
    核函数径向基高斯高斯
    $ \alpha_i $/$ 1.0\times 10^{12} $$ 1.2\times 10^{-6} $
    $ \beta $/$ 2.47\times 10^{-4} $$ 6.62\times 10^{-4} $
    下载: 导出CSV 
    | 显示表格

    根据参数建立退化趋势模型后, 利用首达时间进行剩余寿命预测, 寿命预测的对比图如图5所示.

    图 5  剩余寿命预测值与真实值对比图
    Fig. 5  The RUL comparison between the predicted values and the actual values

    图5中的实线均表示预测得到的剩余寿命, 虚线是实际的剩余寿命值, 图5(c)5(d)中的两条细虚线代表置信区间的上下限, 本文置信水平设为95 %, 由于自回归方法及支持向量机方法并不是从概率角度出发的回归方法, 所以图5(a)5(b)不存在置信区间上下限. 进一步观察4个子图中的预测值, 可以看出, 图5(a)出现误差累计, 图5(b)5(c)波动性较大, 且图5(c)$ 1\;820\; {\rm{min}}$附近即超出置信区间, 图5(d)利用改进的RVM方法, 预测得到的剩余寿命几乎均在置信区间内. 为了更加形象地验证本文提出方法的准确性, 表2对比了4种方法的预测值与真实值间的均方根误差(Root mean squared error, RMSE)及相对平均偏差(Relative average deviation, RAD). 从表2明显看出本文提出的方法无论从均方根误差还是相对平均偏差上均优于其他3种方法.

    表 2  剩余寿命均方根误差(RMSE)及相对平均偏差(RAD)比较
    Table 2  The RMSE and RAD comparison of RUL between different methods
    ARSVMRVM改进RVM
    RMSE (min)1.68603.56542.39051.4586
    RAD (%)29.680235.890020.346811.3533
    下载: 导出CSV 
    | 显示表格

    利用高斯分布和$ t $分布拟合若干数据的直方图, 可以看出, 与高斯分布相比$ t $分布更具有鲁棒性. 图6列出了服从高斯分布的数据以及含有少量异常点的情况, 其中, 图6(a)中的直方图根据随机产生的均值为$ 0 $、方差为$ 1 $$ 500 $个随机数绘制而成, 图6(b)中的直方图在图6(a)的基础上, 在区间$ [8,10] $上添加了$ 30 $个随机异常点, 实线均表示使用$ t $分布对相应的直方图进行的拟合, 虚线表示使用高斯分布进行的拟合.

    图 6  高斯分布与$ t $分布拟合效果对比图
    Fig. 6  The fitting comparison between the Gaussian distribution and the $ t $ distribution

    由于高斯分布是$ t $分布的一种特例, 因此, 图6(a)中的数据既可以用高斯分布又可以用$ t $分布拟合. 但是当多了$ 30 $个异常点时, 从图6(b)可以看出, 高斯分布会受到异常点强烈干扰, 而$ t $分布相对不受影响, 所以$ t $分布比高斯分布鲁棒性好.

    考虑高速列车牵引系统运行中遭受的不确定因素干扰, 本文提出了一种鲁棒的剩余寿命预测算法. 基于稀疏贝叶斯理论建立了具有鲁棒性的相关向量机多步退化模型, 模型中权重及随机误差的精度不再服从均匀分布而是服从$ \Gamma $分布, 进而得到权重与随机误差服从$ t $分布, $ t $分布对不确定因素造成的异常点具有良好的鲁棒性. 此外, 利用极大化似然函数对精度进行估计时, 充分地将精度的先验概率密度函数考虑到似然函数中.

    然而, 剩余寿命预测未来仍面临诸多挑战. 建立的预测模型是否精确关系到剩余寿命的准确度. 如何建立模型检验的评估机制是一项亟待解决的问题. 预测模型中的参数一般都是利用批次数据训练得到, 而剩余寿命预测考虑时间序列维度, 寻找自适应方法优化参数也是提高剩余寿命预测精度需要克服的重要难题之一.

  • 图  1  中间直流环节结构简图

    Fig.  1  The structure diagram of intermediate DC link

    图  2  上下端电压从平稳过程到退化过程直至停机的演变趋势

    Fig.  2  The voltages evolution from the stationary process to the degradation process

    图  3  上下端电压在支撑电容器退化过程中的平均幅值

    Fig.  3  The mean amplitude of the upper and lower voltages during the capacitors' degradation

    图  4  建模过程中的预测值与真实值进行比较

    Fig.  4  The comparison between predicted values and the actual values in the modeling process

    图  5  剩余寿命预测值与真实值对比图

    Fig.  5  The RUL comparison between the predicted values and the actual values

    图  6  高斯分布与$ t $分布拟合效果对比图

    Fig.  6  The fitting comparison between the Gaussian distribution and the $ t $ distribution

    表  1  退化趋势模型中的参数取值

    Table  1  The parameter values in the degradation model

    SVMRVM改进RVM
    核函数径向基高斯高斯
    $ \alpha_i $/$ 1.0\times 10^{12} $$ 1.2\times 10^{-6} $
    $ \beta $/$ 2.47\times 10^{-4} $$ 6.62\times 10^{-4} $
    下载: 导出CSV

    表  2  剩余寿命均方根误差(RMSE)及相对平均偏差(RAD)比较

    Table  2  The RMSE and RAD comparison of RUL between different methods

    ARSVMRVM改进RVM
    RMSE (min)1.68603.56542.39051.4586
    RAD (%)29.680235.890020.346811.3533
    下载: 导出CSV
  • [1] 姜斌, 吴云凯, 陆宁云, 冒泽慧. 高速列车牵引系统故障诊断与预测技术综述. 控制与决策, 2018, 33(5): 841−855

    1 Jiang Bin, Wu Yun-Kai, Lu Ning-Yun, Mao Ze-Hui. Review of fault diagnosis and prognosis techniques for high-speed railway traction system. Control and Decision, 2018, 33(5): 841−855
    [2] 2 Shang J, Chen M Y, Ji H Q, Zhou D H. Recursive transformed component statistical analysis for incipient fault detection. Automatica, 2017, 80: 313−327 doi: 10.1016/j.automatica.2017.02.028
    [3] 3 Lei Y G, Qiao Z J, Xu X F, Liu J, Niu S T. An underdamped stochastic resonance method with stable-state matching for incipient fault diagnosis of rolling element bearings. Mechanical Systems and Signal Processing, 2017, 94: 148−164 doi: 10.1016/j.ymssp.2017.02.041
    [4] 周东华, 纪洪泉, 何潇. 高速列车信息控制系统的故障诊断技术. 自动化学报, 2018, 44(7): 1153−1164

    4 Zhou Dong-Hua, Ji Hong-Quan, He Xiao. Fault diagnosis techniques for the information control system of high-speed trains. Acta Automatica Sinica, 2018, 44(7): 1153−1164
    [5] 5 Tobon-Mejia D A, Medjaher K, Zerhouni N, Tripot, G. A data-driven failure prognostics method based on mixture of Gaussians hidden Markov models. IEEE Transactions on Reliability, 2012, 61(2): 491−503 doi: 10.1109/TR.2012.2194177
    [6] 6 Haque M S, Choi S, Baek J. Auxiliary particle filtering-based estimation of remaining useful life of IGBT. IEEE Transactions on Industrial Electronics, 2018, 65(3): 2693−2703 doi: 10.1109/TIE.2017.2740856
    [7] 7 Tseng S T, Balakrishnan N, Tsai C C. Optimal step-stress accelerated degradation test plan for gamma degradation processes. IEEE Transactions on Reliability, 2009, 58(4): 611−618 doi: 10.1109/TR.2009.2033734
    [8] 8 Zhai Q, Ye Z S. RUL prediction of deteriorating products using an adaptive wiener process model. IEEE Transactions on Industrial Informatics, 2017, 13(6): 2911−2921 doi: 10.1109/TII.2017.2684821
    [9] 9 Si X S, Wang W B, Hu C H, Zhou D H. Remaining useful life estimation-a review on the statistical data driven approaches. European Journal of Operational Research, 2011, 213(1): 1−14 doi: 10.1016/j.ejor.2010.11.018
    [10] 10 Lu N Y, Yao Y, Gao F R. Two-dimensional dynamic PCA for batch process monitoring. AIChE Journal, 2005, 51(12): 3300−3304 doi: 10.1002/aic.10568
    [11] 11 Zhao C H, Gao F R. Online fault prognosis with relative deviation analysis and vector autoregressive modeling. Chemical Engineering Science, 2015, 138(22): 531−543
    [12] 12 Zhao C H, Gao F R. Critical-to-fault-degradation variable analysis and direction extraction for online fault prognostic. IEEE Transactions on Control Systems Technology, 2017, 25(3): 842−854 doi: 10.1109/TCST.2016.2576018
    [13] 13 Li Y, Lu N Y, Wang X L, Jiang B. Islanding fault detection based on data-driven approach with active developed reactive power variation. Neurocomputing, 2019, 337(14): 97−109
    [14] 14 Pham H T, Yang B S, Nguyen T T. Machine performance degradation assessment and remaining useful life prediction using proportional hazard model and support vector machine. Mechanical Systems and Signal Processing, 2012, 32: 320−330 doi: 10.1016/j.ymssp.2012.02.015
    [15] 15 Huang H Z, Wang H K, Li Y F, Zhang L L, Liu Z L. Support vector machine based estimation of remaining useful life: current research status and future trends. Journal of Mechanical Science and Technology, 2015, 29(1): 151−163 doi: 10.1007/s12206-014-1222-z
    [16] 16 Tipping M E. Sparse bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 2001, 1(3): 211−44
    [17] 17 Yu H Y, Wu Z H, Chen D W, Ma X L. Probabilistic prediction of bus headway using relevance vector machine regression. IEEE Transactions on Intelligent Transportation Systems, 2017, 18(7): 1772−1781 doi: 10.1109/TITS.2016.2620483
    [18] 18 Wu Y M, Breaz E, Gao F, Miraoui A. A modified relevance vector machine for PEM fuel-cell stack aging prediction. IEEE Transactions on Industry Applications, 2016, 52(3): 2573−2581 doi: 10.1109/TIA.2016.2524402
    [19] 19 Saha B, Goebel K, Poll S, Christophersen J. Prognostics methods for battery health monitoring using a Bayesian framework. IEEE Transactions on Instrumentation and Measurement, 2009, 58(2): 291−296 doi: 10.1109/TIM.2008.2005965
    [20] 20 Widodo A, Kim, E Y, Son J D, Yang B S, Tan A C, Gu D S, Choid B K, Mathew J. Fault diagnosis of low speed bearing based on relevance vector machine and support vector machine. Expert Systems with Applications, 2009, 36(3): 7252−7261 doi: 10.1016/j.eswa.2008.09.033
    [21] 21 Zio E, Di Maio F. Fatigue crack growth estimation by relevance vector machine. Expert Systems with Applications, 2012, 39(12): 10681−10692 doi: 10.1016/j.eswa.2012.02.199
    [22] 22 Widodo A, Yang B S. Application of relevance vector machine and survival probability to machine degradation assessment. Expert Systems with Applications, 2011, 38(3): 2592−2599 doi: 10.1016/j.eswa.2010.08.049
    [23] Wang X L, Jiang B, Lu N Y. Adaptive relevant vector machine based RUL prediction under uncertain conditions. ISA Transactions, 2018, DOI: 10.1016/j.isatra. 2018.11.024
    [24] 24 Yang C H, Yang C, Peng T, Yang X Y, Gui W H. A fault-injection strategy for traction drive control systems. IEEE Transactions on Industrial Electronics, 2017, 64(7): 5719−5727 doi: 10.1109/TIE.2017.2674610
    [25] 25 Yang X Y, Yang C H, Peng T, Chen Z W, Liu B, Gui W H. Hardware-in-the-loop fault injection for traction control system. IEEE Journal of Emerging and Selected Topics in Power Electronics, 2018, 6(2): 696−706 doi: 10.1109/JESTPE.2018.2794339
  • 期刊类型引用(7)

    1. 葛磊蛟,金雨含,李小平,陈艳波. 人工智能赋能电力电容器运行状态感知及优化控制技术综述. 电气工程学报. 2024(04): 1-12 . 百度学术
    2. 王力,石立超. 基于改进相关向量机的模拟电路故障预测. 计算机应用与软件. 2023(03): 52-60 . 百度学术
    3. 汪翔,李小波,吴浩,吴竑霖. 基于非线性维纳过程的电解电容剩余寿命预测. 电子元件与材料. 2023(03): 334-340+346 . 百度学术
    4. 钟麦英,王钦,彭涛,席霄鹏,杨超,薛婷. 高速列车牵引传动系统运行状态监测技术综述. 山东科技大学学报(自然科学版). 2023(02): 88-97 . 百度学术
    5. 李天梅,司小胜,张建勋. 多源传感监测线性退化设备数模联动的剩余寿命预测方法. 航空学报. 2023(08): 94-112 . 百度学术
    6. 成庶,刘嘉文,伍珣,于天剑,向超群,袁东辉. 基于特征提取与误差补偿的金属化薄膜电容器剩余寿命预测. 中国电机工程学报. 2022(07): 2672-2681 . 百度学术
    7. 陆宁云,陈闯,姜斌,邢尹. 复杂系统维护策略最新研究进展:从视情维护到预测性维护. 自动化学报. 2021(01): 1-17 . 本站查看

    其他类型引用(12)

  • 加载中
图(6) / 表(2)
计量
  • 文章访问数:  2184
  • HTML全文浏览量:  638
  • PDF下载量:  311
  • 被引次数: 19
出版历程
  • 收稿日期:  2019-03-20
  • 录用日期:  2019-09-02
  • 刊出日期:  2019-12-01

目录

/

返回文章
返回