CN104977357B - 一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法 - Google Patents

Abstract

本发明公开了一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，其步骤如下：步骤一：钢轨裂纹声发射信号到达时间自动识别；步骤二：建立有色轮轨接触噪声AR模型及其噪声方程；步骤三：建立钢轨裂纹信号时变参数AR模型及钢轨裂纹信号卡尔曼滤波基本方程；步骤四：有色轮轨接触噪声卡尔曼滤波估计钢轨裂纹信号。本发明具有如下优点：1)在检测到钢轨裂纹信号的基础上，进一步建立钢轨裂纹信号的时变参数AR模型，采用有色测量噪声卡尔曼滤波方法，直接对钢轨裂纹信号进行估计，提取出钢轨裂纹信号；2)将有色轮轨接触噪声情况下的卡尔曼滤波递推公式化简为与一般卡尔曼滤波递推公式一致的一般形式，简化了算法，减小了算法复杂度。

Description

一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法

技术领域

本发明涉及一种高速铁路钢轨声发射技术探伤的去噪方法，具体涉及一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法。

背景技术

目前我国高速铁路蓬勃发展，但在行车过程中由于钢轨长期受到挤压形变和疲劳磨损，会使钢轨表面和内部发生裂纹、折断或其它形式的伤损，由钢轨裂纹扩展而成的钢轨断裂是列车出轨事故主要原因，而高速列车行车速度的提高使其产生裂纹的概率大大增加。现有的大型探伤车和手推式探伤仪占道时间长、作业效率低，不适合高速铁路的钢轨探伤。声发射技术不同于传统的铁路伤损检测技术(超声技术和电磁感应技术)，是一种动态的无损检测方法，具有实时性好、敏感性强等特点，不仅能够检测钢轨表面的裂纹，而且能够感知钢轨内部伤损的发生，所以声发射技术非常适合于钢轨裂纹的在线检测。然而声发射技术由于其敏感性，容易受到外界噪声的干扰，有效的裂纹信号检测的同时伴随着噪声信号。当列车速度较大时，所产生的噪声信号将有效裂纹信号完全淹没，导致有效裂纹信号无法辨别，是阻碍高速情况下声发射技术探伤应用的主要问题。基于噪声信号的高速淹没钢轨裂纹信号检测，由于没有考虑钢轨裂纹信号的特征，所以得不到量化的钢轨裂纹信号。

发明内容

本发明在检测到钢轨裂纹声发射信号的前提下，为了进一步得到钢轨裂纹声发射信号，提出一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法。此方法能够抑制噪声信号，从噪声中提取出不同行车速度下的钢轨裂纹信号，确定裂纹声发射信号的波形，为钢轨裂纹伤损特征提取与分类提供进一步的指导。

本发明的目的是通过以下技术方案实现的：

一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，包括以下步骤：利用全局AIC值法自动判断检测到的钢轨裂纹信号到达时间，划分噪声信号段和钢轨裂纹信号段；建立原始噪声信号的AR模型及有色测量噪声方程；建立钢轨裂纹信号时变参数AR模型及卡尔曼滤波状态方程及新的卡尔曼滤波测量方程；采用有色测量噪声卡尔曼滤波递推算法得到钢轨裂纹信号的估计。如图1所示，具体步骤如下：

步骤一：钢轨裂纹声发射信号到达时间自动识别，流程如图2所示。

1)对含噪声发射信号S进行钢轨裂纹信号检测；

2)取检测到的钢轨裂纹声发射信号投影P；

3)全局AIC值法对钢轨裂纹声发射信号投影P判断钢轨裂纹信号到达时间T：

上式为p阶AR模型AIC值计算公式，式中为窗口内信号AR模型x(k)+a₁x(k-1)+a₂x(k-2)…+a_px(k-p)＝w(k)中零均值白噪声模型误差w(k)的方差，N为建模信号长度；

全局AIC值定义为：

式中K＝1,2,…,N_n+N_s为噪声信号与钢轨裂纹信号的分界点，将声发射信号投影划分为起始于点N_n的噪声窗口序列{x_n(k)}和结束于点N_s的钢轨裂纹信号窗口序列{x_s(k)}，和分别为噪声窗口和钢轨裂纹信号窗口的AR模型误差方差，p_n和p_s分别为噪声窗口和钢轨裂纹信号窗口AR模型的阶数。全局AIC值由噪声信号窗口AIC_n和钢轨裂纹信号窗口AIC_s组成，使全局AIC值达到最小的K点即为钢轨裂纹信号到达时间T。原始声发射信号S中声发射信号到达时间T以前的数据段为噪声信号{n(k)}，检测到的钢轨裂纹声发射信号投影P中到达时间T以后的数据段为钢轨裂纹信号{x(k)}。

步骤二：建立有色轮轨接触噪声AR模型及其噪声方程。

1)噪声信号序列{n(k)}一阶AR模型：

n(k)+a₁n(k-1)＝w₁(k)，

其中a₁为模型参数，模型误差w₁(k)是均值为零方差为的高斯白噪声；

2)由噪声信号一阶AR模型建立有色轮轨接触噪声方程：

V_k＝Ψ_k,k-1V_k-1+ζ_k-1，

式中转移矩阵Ψ_k,k-1＝-a₁，测量噪声系统噪声ζ_k-1＝w₁(k)。

步骤三：建立钢轨裂纹信号时变参数AR模型及钢轨裂纹信号卡尔曼滤波基本方程，流程如图3所示。

1)钢轨裂纹信号{x(k)}的时变参数AR模型：

x(k)+a₁(k-1)x(k-1)+a₂(k-2)x(k-2)+…+a_p(k-p)x(k-p)＝w_p(k)，

其中a₁(k),a₂(k),…,a_p(k)为模型时变参数，模型误差w_p(k)是均值为零方差为的高斯白噪声；

2)时变参数基函数展开：

式中令基函数对应系数则

3)定义矢量则p阶AR模型可表示为：

将时变参数估计问题转化为时不变参数θ的估计问题，进而得到非平稳钢轨裂纹信号的时变参数AR模型；

4)由时变参数建立卡尔曼滤波状态方程：

定义状态量：

状态向量X_k为p维，p维状态分量有关系：

由钢轨裂纹信号时变参数AR模型关系得状态方程：

状态方程的一步转移矩阵：

系统噪声驱动矩阵：

则可得卡尔曼滤波状态方程：

X_k＝Φ_k,k-1X_k-1+Γ_k,k-1W_k-1，

式中系统噪声W_k＝w_p(k)；

5)测量扩增法建立卡尔曼滤波新的测量方程：

原卡尔曼滤波测量方程：

Z_k＝H_kX_k+V_k，

式中H_k为测量矩阵，V_k为测量噪声，

定义新的测量量测量矩阵测量噪声

其中Z_k＝z(k)，{z(k)}为测得含噪声发射信号，为零均值的白噪声，方差为：

式中R_k为ζ_k-1的方差，Q_k为W_k的方差，

步骤四：有色轮轨接触噪声卡尔曼滤波估计钢轨裂纹信号。

由于轮轨接触噪声不是高斯白噪声，所以需要采用有色测量噪声卡尔曼滤波方法对钢轨裂纹信号进行估计，具体方法如下：

1)有色测量噪声卡尔曼滤波递推算法：

状态一步预测

状态一步预测误差方差阵P_k,k-1：

滤波增益矩阵K_k：

状态滤波估计

滤波误差方差阵P_k：

式中I为单位阵，只要给定初值和P₀，就可以根据测量递推计算得到k时刻的估计

2)滤波初值确定：

式中是X₀的均值，是Z₀的均值，表示Z₀的方差。

本发明与现有技术相比，具有如下优点：

1)在检测到钢轨裂纹信号的基础上，进一步建立钢轨裂纹信号的时变参数AR模型，采用有色测量噪声卡尔曼滤波方法，直接对钢轨裂纹信号进行估计，提取出钢轨裂纹信号；

2)将有色轮轨接触噪声情况下的卡尔曼滤波递推公式化简为与一般卡尔曼滤波递推公式一致的一般形式，简化了算法，减小了算法复杂度。

附图说明

图1为本发明的方框图；

图2为钢轨裂纹声发射信号达到时间自动识别流程图；

图3为有色轮轨接触噪声情况下的卡尔曼滤波基本方程建立流程图；

图4为48km/h检测到的钢轨裂纹信号投影；

图5为48km/h速度下含噪声发射信号；

图6为48km/h速度下提取去噪后得到的钢轨裂纹信号投影；

图7为140km/h速度下声发射信号原始数据；

图8为140km/h速度下含噪声发射信号；

图9为140km/h速度下提取去噪后得钢轨裂纹信号。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

本发明提供了一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，包括以下步骤：对各速度下含噪声发射信号截取检测到钢轨裂纹信号段，并对检测后的声发射信号做钢轨裂纹信号到达时间自动识别，分别对噪声部分建立AR模型，对钢轨裂纹信号部分建立时变参数AR模型，进而建立有色测量噪声下的卡尔曼滤波基本方程，再由有色测量噪声卡尔曼滤波递推算法，最终得到钢轨裂纹信号估计。从6km/h～140km/h每隔6km/h在不同速度下分别测得多组实验数据，在有效信号到达以前的声发射信号均为此速度下的噪声信号，取噪声信号长度为10000点，钢轨裂纹信号长度为100000点。

执行步骤一：全局AIC值法判断钢轨裂纹声发射信号到达时间，在噪声窗口{x_n(k)}和钢轨裂纹声发射信号窗口{x_s(k)}内以1为单位变化K的取值，画出随K变化的全局AIC值曲线，曲线最低点时所取得的K即为钢轨裂纹信号达到时间T，以48km/h的钢轨裂纹声发射信号投影为例，到达时间的判断见图4。

执行步骤二：建立噪声信号{n(k)}的一阶AR模型，并由模型参数建立有色测量噪声方程V_k＝Ψ_k,k-1V_k-1+ζ_k-1，为有色测量噪声情况下的卡尔曼滤波做准备。

执行步骤三：基函数展开法建立钢轨裂纹信号的时变参数AR模型，时变参数AR模型阶数p的选取取决于建模精度和计算时间，模型阶数越高模型精度越高最终的估计效果越好，但同时计算量越大计算时间越长，二者需折中；时变参数AR模型的基函数阶数m的选取取决于所需计算时间，一般不宜取过大，导致计算量增大，计算速度下降。进而以钢轨裂纹声发射信号为状态量，建立有色测量噪声卡尔曼滤波状态方程X_k＝Φ_k,k-1X_k-1+Γ_k,k-1W_k-1。

执行步骤四：测量扩增法建立新的卡尔曼滤波测量方程代入有色测量噪声情况下的卡尔曼滤波递推算法，得到钢轨裂纹信号的估计

图5为48km/h速度下含噪声发生信号，图6为48km/h速度下提取去噪后得到的钢轨裂纹信号。48km/h为钢轨裂纹信号没有被噪声淹没情况的典型速度，去噪前噪声信号幅值在20mV，含噪钢轨裂纹信号的幅值为80mV，噪声信号幅值大约占含噪钢轨裂纹信号幅值的25％；而图6去噪后的有效信号中，噪声信号幅值为4mV左右，而钢轨裂纹信号的幅值约为25mV，噪声信号与钢轨裂纹信号幅值比约为16％。去噪前后对比，噪声信号与钢轨裂纹信号幅值比变小，说明基于卡尔曼滤波的信号去噪检测方法可以大幅度减小噪声，增强有效信号，使信噪比增大。

图7为140km/h速度下声发射信号原始数据，其中包括钢轨裂纹信号被完全淹没信号段如8所示，图9为提取去噪后得钢轨裂纹信号，140km/h时钢轨裂纹信号已经完全淹没在噪声信号中无法辨别。图8中140km/h时噪声信号幅值在85mV左右，在时间为4.06s时有一个饱和信号点(幅值显示为100mV)。去噪后可以从被噪声淹没的信号中提取出钢轨裂纹信号，使噪声得到抑制。经提取去噪后的信号中仍然含有噪声，噪声的幅值为4mV，有效信号幅值在9.5mV左右，最终噪声信号与钢轨裂纹信号幅值比约为42.1％。与被淹没情况下相比，基于卡尔曼滤波的提取去噪方案可以达到抑制噪声，提取钢轨裂纹信号的目的。

综合实施例的上述分析，对于高速铁路钢轨探伤中钢轨裂纹声发射信号的提取，本发明采用基于卡尔曼滤波的钢轨裂纹信号提取去噪方法。采用有色测量噪声卡尔曼滤波方法，分别建立噪声信号和钢轨裂纹信号模型，对钢轨裂纹信号直接估计，得到去噪后的钢轨裂纹信号。由于此法是基于钢轨裂纹信号的方法，所以此法去噪效果取决于模型的建模效果，当建模效果较好时，无论有效信号是否淹没在噪声信号中，都能够从噪声信号中提取出钢轨裂纹声发射信号。

Claims (1)

1.一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，其特征在于所述方法步骤如下：

步骤一、钢轨裂纹声发射信号到达时间自动识别：

1)对含噪声发射信号S进行钢轨裂纹信号检测；

2)取检测到的钢轨裂纹声发射信号投影P；

3)全局AIC值法对钢轨裂纹声发射信号投影P判断钢轨裂纹信号到达时间T：

$AIC\left(p\right)=ln\left({\mathrm{\σ}}_w^2\right)+2p/N;$

上式为p阶AR模型AIC值计算公式，式中为窗口内信号AR模型x(k)+a₁x(k-1)+a₂x(k-2)…+a_px(k-p)＝w(k)中零均值白噪声模型误差w(k)的方差，N为建模信号长度；

全局AIC值定义为：

$AIC={\mathrm{AIC}}_n+{\mathrm{AIC}}_s={\mathrm{ln\σ}}_{wn}^2+\frac{2p_n}{K-N_n}+{\mathrm{ln\σ}}_{ws}^2+\frac{2p_s}{N_s-K};$

式中全局AIC值由噪声信号窗口AIC_n和钢轨裂纹信号窗口AIC_s组成，K＝1,2,…,N_n+N_s为噪声信号与钢轨裂纹信号的分界点，将声发射信号投影划分为起始于点N_n的噪声窗口序列{x_n(k)}和结束于点N_s的钢轨裂纹信号窗口序列{x_s(k)}，和分别为噪声窗口和钢轨裂纹信号窗口的AR模型误差方差，p_n和p_s分别为噪声窗口和钢轨裂纹信号窗口AR模型的阶数；

步骤二、建立有色轮轨接触噪声AR模型及其噪声方程：

1)噪声信号序列{n(k)}一阶AR模型：

n(k)+a₁n(k-1)＝w₁(k)，

其中a₁为模型参数，模型误差w₁(k)是均值为零方差为的高斯白噪声；

2)由噪声信号一阶AR模型建立有色轮轨接触噪声方程：

V_k＝Ψ_k,k-1V_k-1+ζ_k-1，

式中转移矩阵Ψ_k,k-1＝-a₁，测量噪声系统噪声ζ_k-1＝w₁(k)；

步骤三、建立钢轨裂纹信号时变参数AR模型及钢轨裂纹信号卡尔曼滤波基本方程：

1)钢轨裂纹信号{x(k)}的时变参数AR模型：

x(k)+a₁(k-1)x(k-1)+a₂(k-2)x(k-2)+…+a_p(k-p)x(k-p)＝w_p(k)，

其中a₁(k),a₂(k),…,a_p(k)为模型时变参数，模型误差w_p(k)是均值为零方差为的高斯白噪声；

2)时变参数基函数展开：

$a_i\left(k\right)=a_{i0}+a_{i1}k+a_{i2}\frac{k^2}{2}+...+a_{im}\frac{k^m}{m!},(i=1,2,...,p),$

式中令基函数对应系数则

3)定义矢量x(t)＝[f₀(t)x(t),f₁(t)x(t),…,f_m(t)x(t)]^T，则p阶AR模型可表示为：

$x\left(k\right)+\[{\stackrel{\‾}{x}}^T(k-1),{\stackrel{\‾}{x}}^T(k-2),...,{\stackrel{\‾}{x}}^T(k-p)\]\mathrm{\θ}=w_p\left(k\right),$

将时变参数估计问题转化为时不变参数θ的估计问题，进而得到非平稳钢轨裂纹信号的时变参数AR模型；

4)由时变参数建立卡尔曼滤波状态方程：

定义状态量：

$X_k=\left[\begin{array}{c}{x}_1\left(k\right)\\ x_2\left(k\right)\\ .\\ .\\ .\\ x_{p-1}\left(k\right)\\ x_p\left(k\right)\end{array}\right]=\left[\begin{array}{c}x(k-p)\\ x(k-p+1)\\ .\\ .\\ .\\ x(k-2)\\ x(k-1)\end{array}\right];$

状态向量X_k为p维，p维状态分量有关系：

$\left\{\begin{array}{c}{x}_1\left(k\right)=x_2(k-1)\\ x_2\left(k\right)=x_3(k-1)\\ .\\ .\\ .\\ x_{p-1}\left(k\right)=x_p(k-1)\end{array}\right.;$

由钢轨裂纹信号时变参数AR模型关系得状态方程：

状态方程的一步转移矩阵：

系统噪声驱动矩阵：

${\mathrm{\Γ}}_{k,k-1}=\left[\begin{array}{c}0\\ 0\\ .\\ .\\ .\\ 0\\ 1\end{array}\right];$

则可得卡尔曼滤波状态方程：

X_k＝Φ_k,k-1X_k-1+Γ_k,k-1W_k-1，

式中系统噪声W_k＝w_p(k)；

5)测量扩增法建立卡尔曼滤波新的测量方程：

$Z_k^{*}=H_k^{*}{X}_k+V_k^{*};$

原卡尔曼滤波测量方程：

Z_k＝H_kX_k+V_k，

式中H_k为测量矩阵，V_k为测量噪声，

定义新的测量量测量矩阵测量噪声

$\left\{\begin{array}{c}{Z}_k^{*}=Z_k-{\mathrm{\Ψ}}_{k,k-1}{Z}_{k-1}\\ H_k^{*}=(H_k{\mathrm{\Φ}}_{k,k-1}-{\mathrm{\Ψ}}_{k,k-1}{H}_{k,k-1}){\mathrm{\Φ}}_{k,k-1}^{-1}\\ V_k^{*}=H_k{\mathrm{\Γ}}_{k,k-1}{W}_{k-1}+{\mathrm{\ζ}}_{k-1}\end{array}\right.;$

其中Z_k＝z(k)，{z(k)}为测得含噪声发射信号，为零均值的白噪声，方差为：

$R_k^{*}=H_k{\mathrm{\Γ}}_{k,k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k,k-1}^T{H}_k^T+R_{k-1},$

式中R_k为ζ_k-1的方差，Q_k为W_k的方差，

步骤四、有色轮轨接触噪声卡尔曼滤波估计钢轨裂纹信号：

1)有色测量噪声卡尔曼滤波递推算法：

状态一步预测

${\hat{X}}_{k,k-1}={\mathrm{\Φ}}_{k,k-1}{\hat{X}}_{k-1};$

状态一步预测误差方差阵P_k,k-1：

$P_{k,k-1}={\mathrm{\Φ}}_{k,k-1}{P}_{k-1}{\mathrm{\Φ}}_{k,k-1}^T+{\mathrm{\Γ}}_{k,k-1}{Q}_{k-1}{\mathrm{\Γ}}_{k,k-1}^T;$

滤波增益矩阵K_k：

$K_k=P_{k,k-1}{H}_k^{*T}{(H_k^{*}{P}_{k,k-1}{H}_k^{*T}+R_k^{*})}^{-1};$

状态滤波估计

${\hat{X}}_k={\hat{X}}_{k,k-1}+K_k{\[Z_k^{*}-H_k^{*}{\hat{X}}_{k,k-1}\]}^{-1};$

滤波误差方差阵P_k：

$P_k=\[I-K_k{H}_k^{*}\]P_{k,k-1};$

式中I为单位阵；

2)滤波初值确定：

${\hat{X}}_0={\mathrm{\μ}}_{X_0}+C_{X_0}{H}_0^T{(H_0{C}_{X_0}{H}_0^T+R_0)}^{-1}(Z_0-H_0{\mathrm{\μ}}_{Z_0});$

$P_0=C_{X_0}+C_{X_0}{H}_0^T{(H_0{C}_{X_0}{H}_0^T+R_0)}^{-1}{H}_0{C}_{X_0};$

式中是X₀的均值，是Z₀的均值，表示X₀的方差。