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

摘要

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

主权项

一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，其特征在于所述方法步骤如下：步骤一、钢轨裂纹声发射信号到达时间自动识别：1)对含噪声发射信号S进行钢轨裂纹信号检测；2)取检测到的钢轨裂纹声发射信号投影P；3)全局AIC值法对钢轨裂纹声发射信号投影P判断钢轨裂纹信号到达时间T：$AIC\left(p\right)=ln\left({\sigma }_w^2\right)+2p/N;$上式为p阶AR模型AIC值计算公式，式中为窗口内信号AR模型x(k)+a₁x(k11)+a₂x(k‑2)…+a_px(k‑p)＝w(k)中零均值白噪声模型误差w(k)的方差，N为建模信号长度；全局AIC值定义为：$AIC={\mathrm{AIC}}_n+{\mathrm{AIC}}_s={\mathrm{ln\sigma }}_{wn}^2+\frac{2p_n}{K-N_n}+{\mathrm{ln\sigma }}_{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(k11)＝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),$式中令基函数$\left\{\begin{array}{c}{f}_0\left(k\right)=1\\ f_1\left(k\right)=k\\ f_2\left(k\right)=\frac{k^2}{2}\\ .\\ .\\ .\\ f_m\left(k\right)=\frac{k^m}{m!}\end{array}\right.,$对应系数$\left\{\begin{array}{c}{a}_{i0}\\ a_{i1}\\ a_{i2}\\ .\\ .\\ .\\ a_{im}\end{array}\right.,$则$a_i\left(k\right)=\underset{j=0}{\overset{m}{\Sigma }}{a}_{ij}{f}_j\left(k\right);$3)定义矢量$\bar{x}\left(t\right)={[f_0\left(t\right)x\left(t\right),f_1\left(t\right)x\left(t\right),...,f_m\left(t\right)x\left(t\right)]}^T,$则p阶AR模型可表示为：$x\left(k\right)+[{\bar{x}}^T(k-1),{\bar{x}}^T(k-2),...,{\bar{x}}^T(k-p)]\theta =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模型关系得状态方程：状态方程的一步转移矩阵：系统噪声驱动矩阵：${\Gamma }_{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-{\Psi }_{k,k-1}{Z}_{k-1}\\ H_k^{*}=(H_k{\Phi }_{k,k-1}-{\Psi }_{k,k-1}{H}_{k,k-1}){\Phi }_{k,k-1}^{-1}\\ V_k^{*}=H_k{\Gamma }_{k,k-1}{W}_{k-1}+{\zeta }_{k-1}\end{array}\right.;$其中Z_k＝z(k)，{z(k)}为测得含噪声发射信号，为零均值的白噪声，方差为：$R_k^{*}=H_k{\Gamma }_{k,k-1}{Q}_{k-1}{\Gamma }_{k,k-1}^T{H}_k^T+R_{k-1},$式中R_k为ζ_k‑1的方差，Q_k为W_k的方差，步骤四、有色轮轨接触噪声卡尔曼滤波估计钢轨裂纹信号：1)有色测量噪声卡尔曼滤波递推算法：状态一步预测${\hat{X}}_{k,k-1}={\Phi }_{k,k-1}{\hat{X}}_{k-1};$状态一步预测误差方差阵P_k，k‑1：$P_{k,k-1}={\Phi }_{k,k-1}{P}_{k-1}{\Phi }_{k,k-1}^T+{\Gamma }_{k,k-1}{Q}_{k-1}{\Gamma }_{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={\mu }_{X_0}+C_{X_0}{H}_0^T{(H_0{C}_{X_0}{H}_0^T+R_0)}^{-1}(Z_0-H_0{\mu }_{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₀的均值，$C_{X_0}=E\left[X_0{X}_0^T\right]$表示X₀的方差。