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

权利要求书

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₀的方差。