一种基于振动信号的铁路钢轨几何形变检测方法

摘要

本发明涉及一种基于振动信号的铁路钢轨几何形变检测方法，属于列车轨道故障监控与诊断技术领域。基于列车不同位置几何形变种类的振动数据，对信号进行短时傅里叶变换得到振动频谱；提取频谱特征，采用主成分分析对提取的特征降维，得到主频谱特征；针对每个几何形变种类的支持向量机，建立每个几何形变种类的信度计算模型；在线检测时，根据实时测量振动数据，提取频谱特征，进而获得主频谱特征，采用离线计算得到的几何形变种类的信度模型计算隶属于每类的信度，根据信度得到诊断信度向量。本方法可精确诊断并定位几何形变种类的情况。安装在实时运行的高速列车上即可进行轨道的几何形变种类在线分析与诊断。

主权项

1.一种基于振动信号的铁路钢轨几何形变检测方法，其特征在于该方法包括以下各步骤：（1）设定列车钢轨几何形变种类集合为Θ＝{F₁,F₂,F₃,F₄}，其中，F_j代表几何形变种类集合Θ中的第j个类型的几何形变种类，j＝1,2,3,4，第一种几何形变为钢轨实际的轨道中心线与理想的轨道中心线沿长度方向的垂向几何位置偏差，第二种几何形变为钢轨表面呈波面状磨损，第三种几何形变为钢轨表面局部区域剥落，第四种几何形变为钢轨表面无几何形变；（2）在列车的轴箱、转向架和车厢底部分别设置振动传感器s_k，列车运行时，三个振动传感器分别采集的相应部位的振动信号，得到三个振动信号序列{x_k(m)}，k＝1,2,3，将轴箱上测得的振动信号序列记为{x₁(m)}，将转向架上测得的振动信号序列记为{x₂(m)}，车厢底部测得的振动信号序列记为{x₃(m)}，其中，m为采样时刻，m＝1,...,M，M代表时间序列的长度，对上述三个振动信号序列分别进行短时傅里叶变换，得到三个振动频谱$X_k^m\left(\omega \right)=\underset{n=-\infty }{\overset{\infty }{\Sigma }}\bar{\omega }(n-m)x_k\left(n\right)e^{-\mathrm{j\omega n}}$上式中，n为计算参数，为汉明窗函数，窗函数为：$\bar{\omega }(n-m)=0.54-0.46\cos \left(\frac{\pi (n-m)}{100}\right),$0≤n-m≤200；（3）从振动频谱中，提取11个频谱特征，包括熵、最大幅值、最大幅值角频率、平均均值、均方根值、方差、歪度、峭度、第一频带的平均能量、第二频带的平均能量和第三频带的平均能量，得到频谱特征提取频谱特征的方法如下：（a）将振动频谱的频率范围等分为100个区间，每个频率区间出现的概率密度为l＝1,...,100，则熵$y_{k,1}^m={-p}_{k,l}^m*\log \left(p_{k,l}^m\right);$（b）从振动频谱中直接寻找到最大幅值为i＝1,2,3,...,I，其中X(ω_i)表示与角频率ω_i相对应的幅值，I为角频率ω离散取值的个数，定义ω₁为最低角频率；（c）与上述最大幅值相对应的角频率ω_i为最大幅值角频率（d）平均均值为：$y_{k,4}^m=\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{X}_k^m\left({\omega }_i\right);$（e）均方根值为：$y_{k,5}^m=\sqrt{\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^2};$（f）方差为：$y_{k,6}^m=\frac{1}{I-1}\underset{i=1}{\overset{I}{\Sigma }}{(X_k^m\left({\omega }_i\right)-y_{k,4}^m)}^2,$为（d）求得的平均均值；（g）歪度为：$y_{k,7}^m=\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^3;$（h）峭度为：$y_{k,8}^m=\frac{1}{I}\underset{i=1}{\overset{I}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^4;$（i）第一频带平均能量$y_{k,9}^m=\frac{1}{I_1}\underset{i=1}{\overset{I_1}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^2,$其中（j）第二频带平均能量$y_{k,10}^m=\frac{1}{I_2-I_1}\underset{i=I_1+1}{\overset{I_2}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^2,$（k）第三频带平均能量$y_{k,11}^m=\frac{1}{I-I_2}\underset{i=I_2+1}{\overset{I}{\Sigma }}{X}_k^m{\left({\omega }_i\right)}^2,$（4）对上述频谱特征进行中心化和归一化处理：对中的第p维频谱特征，提取特征后形成第p维频谱特征集合为其中m为采样时刻，m＝1,...,M，M代表时间序列的长度，对第p维频谱特征集合进行中心化和归一化的步骤为：对均值化得到：${\bar{y}}_{k,p}=\underset{m}{\Sigma }{y}_{k,p}^m$对中心化得到：对归一化得到：对$Y_k^m=[y_{k,1}^m,y_{k,2}^m,···,y_{k,11}^m]$中的每维频谱特征进行上述计算，得到${\hat{Y}}_k^m=[{\hat{y}}_{k,1}^m,{\hat{y}}_{k,2}^m,···,{\hat{y}}_{k,11}^m];$（5）传感器s_k对列车进行测量，设定测量的是第j个类型的几何形变F_j，采样点个数为M_j，重复步骤（2）－步骤（4），对测量得到的振动信号{x_k(m_j)}进行处理，得到频谱特征其中，m_j＝1,...,M_j，k＝1,2,3，j＝1,2,3,4，记测量得到的全部振动信号为m'＝1,2,...,M'，其中，M'为总采样点个数并记录步骤（4）的中间变量${\mu }_{k,p}={\bar{y}}_{k,p}^{m^{\prime }}$和计算频谱特征的均值向量μ_j,计算频谱特征的协方差矩阵Σ_j，${\Sigma }_j=\frac{1}{M_j}\underset{m_j=1}{\overset{M_j}{\Sigma }}{({\hat{Y}}_{k,j}^{m_j}-{\mu }_j)}^T({\hat{Y}}_{k,j}^{m_j}-{\mu }_j),$计算频谱特征的总类内离散度矩阵S，（6）对上述总类内离散度矩阵S进行特征值分解，求得S的本征值矩阵Λ和本征向量U：S=UΛU^-1其中，U为11×11维的方阵，方阵中第i列为总类内离散度矩阵S的特征向量u_i，Λ为11×11维的对角矩阵，对角矩阵中对角线上的元素为特征向量u_i相对应的特征值为λ_i，即Λ_ii=λ_i；（7）对特征值λ_i从大到小进行排序，从其中选出前q个最大的特征值，得到与前q个最大的特征值相对应的特征向量u_i'，由q个u_i'组成11×q维的变换矩阵B_k，采用变换矩阵对频谱特征进行变换，得到频谱特征的主频谱特征为：$z_{k,j}^{m_j}=B_k^T{\hat{Y}}_{k,j}^{m_j}$合并总的主频谱特征为m'＝1,2,...,M'，其中，M'为总采样点个数$M^{\prime }=\underset{j}{\Sigma }{M}_j,$本方法取q＝3；（8）根据上述主频谱特征用指示向量y_k,j,m'∈{-1,1}标记几何形变种类F_j的分类，利用带松弛项的线性支持向量机模型，求解该模型，得到判别向量系数w_k,j和判别阈值b_k,j：$\underset{w,b,\xi }{\min }\frac{1}{2}{w_{k,j}}^T{w}_{k,j}+\underset{m=1}{\overset{M}{\Sigma }}{\xi }_{k,j,m^{\prime }}$$\left\{\begin{array}{c}{y}_{k,j,m^{\prime }}({w_{k,j}}^T{z}_k^{m^{\prime }}+b_{k,j})\ge 1-{\xi }_{k,j,m^{\prime }}\\ {\xi }_{k,j,m^{\prime }}\ge 0,m^{\prime }=1,···,M\end{array}\right.$上式中，ξ_k,j,m'为松弛项，w_k,j^T代表向量w_k,j的转置；（9）对传感器s_k在采样时刻n采集的列车振动信号，利用上述步骤（2）—步骤（3）的方法进行处理，得到频谱特征对中的第p维频谱特征，提取特征后形成第p维频谱特征集合为利用步骤(5)得到的μ_k,p和δ_k,p对频谱特征进行中心化和归一化，变换公式为：${\hat{y}}_{k,p}^n=\frac{y_{k,p}^n-{\mu }_{k,p}}{{\delta }_{k,p}}$对中的每维频谱特征进行上述计算，得到频谱特征利用步骤（7）得到的变换矩阵B_k对频谱特征进行变换，得到传感器s_k在采样时刻n的主频谱特征，记为利用下式计算主频谱特征所对应几何形变种类F_j的初始信度$p_{k,j}^n=\left\{\begin{array}{cc}1-\frac{{a\times b}^{-d_{k,j}^n}}{a+1}& d_{k,j}^n\ge 0\\ \frac{b^{d_{k,j}^n}}{a+1}& d_{k,j}^n<0\end{array}\right.$其中，a＝9，b＝3，w_k,j为判别向量系数，b_k,j为判别阈值；（10）利用上述步骤（9）的几何形变种类初始信度得到主频谱特征的初始信度向量其中和由步骤（9）计算得到，为几何形变种类不确定的初始信度，计算方式为：对上述初始信度向量进行归一化，得到主频谱特征的信度向量为：记信度向量为$m_k^n=[m_k^n\left(F_1\right),···,m_k^n\left(F_4\right),$其中，代表与主频谱特征所对应几何形变种类F_j的信度，代表与主频谱特征所对应几何形变种类不确定的信度；（11）使用定义的折扣因子矩阵A，对步骤（10）得到的信度进行折扣处理，k＝1,...,3，j＝1,...,4，采用证据理论进行信度向量的融合，得到铁路钢轨几何形变的检测结果，包括以下步骤：（11-1）定义折扣因子矩阵为：$A=\left[\begin{array}{cccc}{\alpha }_{\mathrm{1,1}}& {\alpha }_{\mathrm{1,2}}& {\alpha }_{\mathrm{1,3}}& {\alpha }_{\mathrm{1,4}}\\ {\alpha }_{\mathrm{2,1}}& {\alpha }_{\mathrm{2,2}}& {\alpha }_{\mathrm{2,3}}& {\alpha }_{\mathrm{2,4}}\\ {\alpha }_{\mathrm{3,1}}& {\alpha }_{\mathrm{3,2}}& {\alpha }_{\mathrm{3,3}}& {\alpha }_{\mathrm{3,4}}\end{array}\right]=\left[\begin{array}{cccc}0.2025& 0.0000& 0.0445& 0.1610\\ 0.2396& 0.1911& 0.0000& 0.0000\\ 0.0000& 0.1122& 0.2481& 0.0127\end{array}\right]$（11－2）对几何形变种类F_j进行检测时，j＝1,...,4，折扣公式为：$m_{k,j}^n\left(F_r\right)=\left\{\begin{array}{cc}(1-{\alpha }_{k,j})·m_k^n\left(F_r\right)& F_r\subseteq \Theta \\ (1-{\alpha }_{k,j})·m_k^n\left(F_r\right)+{\alpha }_{k,j}& F_r=\Theta \end{array}\right.$融合公式为：其中，α_k,j为折扣因子，k＝1,...,3，k'＝1,...,3，F₅为几何形变种类的全集，F₅＝Θ，r＝1,...,5，s＝1,...,5；根据上述融合公式，得到融合结果作为几何形变种类F_j的最终信度；（11－3）遍历几何形变种类，重复步骤（11－2），得到j＝1,...,4，融合之后的最终信度向量为：$m^n=[m_1^n\left(F_1\right),···,m_4^n\left(F_4\right),m_5^n\left(\Theta \right)],$其中，几何形变种类不确定的最终信度$m_5^n\left(\Theta \right)=1-\underset{j=1,...,4}{\max }\left(m_j^n\left(F_j\right)\right)=m_5^n\left(F_5\right);$（11－4）对上述步骤（11—3）的最终信度向量进行归一化处理：得到检测信度向量为：（12）设定一个检测阈值t，根据上述步骤（11）的检测信度向量，对列车钢轨的几何形变种类进行诊断，j＝1,2,3：若mⁿ(F_j)大于设定阈值t，则判定几何形变F_j发生，若mⁿ(F_j)小于或等于设定阈值t，则判定几何形变F_j未发生，t的取值范围为0.8≤t≤1。