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

摘要

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

主权项

一种基于卡尔曼滤波的钢轨裂纹声发射信号提取与去噪方法，其特征在于所述方法步骤如下：步骤一、钢轨裂纹声发射信号到达时间自动识别：1)对含噪声发射信号S进行钢轨裂纹信号检测；2)取检测到的钢轨裂纹声发射信号投影P；3)全局AIC值法对钢轨裂纹声发射信号投影P判断钢轨裂纹信号到达时间T：AIC(p)=ln(σw2)+2p/N;]]>上式为p阶AR模型AIC值计算公式，式中为窗口内信号AR模型x(k)+a1x(k‑1)+a2x(k‑2)…+apx(k‑p)＝w(k)中零均值白噪声模型误差w(k)的方差，N为建模信号长度；全局AIC值定义为：AIC=AICn+AICs=lnσwn2+2pnK-Nn+lnσws2+2psNs-K;]]>式中全局AIC值由噪声信号窗口AICn和钢轨裂纹信号窗口AICs组成，K＝1,2,…,Nn+Ns为噪声信号与钢轨裂纹信号的分界点，将声发射信号投影划分为起始于点Nn的噪声窗口序列{xn(k)}和结束于点Ns的钢轨裂纹信号窗口序列{xs(k)}，和分别为噪声窗口和钢轨裂纹信号窗口的AR模型误差方差，pn和ps分别为噪声窗口和钢轨裂纹信号窗口AR模型的阶数；步骤二、建立有色轮轨接触噪声AR模型及其噪声方程：1)噪声信号序列{n(k)}一阶AR模型：n(k)+a1n(k‑1)＝w1(k)，其中a1为模型参数，模型误差w1(k)是均值为零方差为的高斯白噪声；2)由噪声信号一阶AR模型建立有色轮轨接触噪声方程：Vk＝Ψk,k‑1Vk‑1+ζk‑1，式中转移矩阵Ψk,k‑1＝‑a1，测量噪声系统噪声ζk‑1＝w1(k)；步骤三、建立钢轨裂纹信号时变参数AR模型及钢轨裂纹信号卡尔曼滤波基本方程：1)钢轨裂纹信号{x(k)}的时变参数AR模型：x(k)+a1(k‑1)x(k‑1)+a2(k‑2)x(k‑2)+…+ap(k‑p)x(k‑p)＝wp(k)，其中a1(k),a2(k),…,ap(k)为模型时变参数，模型误差wp(k)是均值为零方差为的高斯白噪声；2)时变参数基函数展开：ai(k)=ai0+ai1k+ai2k22+...+aimkmm!,(i=1,2,...,p),]]>式中令基函数对应系数则3)定义矢量x(t)＝[f0(t)x(t),f1(t)x(t),…,fm(t)x(t)]T，则p阶AR模型可表示为：x(k)+[x‾T(k-1),x‾T(k-2),...,x‾T(k-p)]θ=wp(k),]]>将时变参数估计问题转化为时不变参数θ的估计问题，进而得到非平稳钢轨裂纹信号的时变参数AR模型； 4)由时变参数建立卡尔曼滤波状态方程： 定义状态量：Xk=x1(k)x2(k)...xp-1(k)xp(k)=x(k-p)x(k-p+1)...x(k-2)x(k-1);]]> 状态向量Xk为p维，p维状态分量有关系：x1(k)=x2(k-1)x2(k)=x3(k-1)...xp-1(k)=xp(k-1);]]> 由钢轨裂纹信号时变参数AR模型关系得状态方程： 状态方程的一步转移矩阵： 系统噪声驱动矩阵：Γk,k-1=00...01;]]>则可得卡尔曼滤波状态方程：Xk＝Φk,k‑1Xk‑1+Γk,k‑1Wk‑1，式中系统噪声Wk＝wp(k)；5)测量扩增法建立卡尔曼滤波新的测量方程：Zk*=Hk*Xk+Vk*;]]>原卡尔曼滤波测量方程：Zk＝HkXk+Vk，式中Hk为测量矩阵，Vk为测量噪声，定义新的测量量测量矩阵测量噪声Zk*=Zk-Ψk,k-1Zk-1Hk*=(HkΦk,k-1-Ψk,k-1Hk,k-1)Φk,k-1-1Vk*=HkΓk,k-1Wk-1+ζk-1;]]>其中Zk＝z(k)，{z(k)}为测得含噪声发射信号，为零均值的白噪声，方差为：Rk*=HkΓk,k-1Qk-1Γk,k-1THkT+Rk-1,]]>式中Rk为ζk‑1的方差，Qk为Wk的方差，步骤四、有色轮轨接触噪声卡尔曼滤波估计钢轨裂纹信号：1)有色测量噪声卡尔曼滤波递推算法：状态一步预测X^k,k-1=Φk,k-1X^k-1;]]>状态一步预测误差方差阵Pk,k‑1：Pk,k-1=Φk,k-1Pk-1Φk,k-1T+Γk,k-1Qk-1Γk,k-1T;]]>滤波增益矩阵Kk：Kk=Pk,k-1Hk*T(Hk*Pk,k-1Hk*T+Rk*)-1;]]>状态滤波估计X^k=X^k,k-1+Kk[Zk*-Hk*X^k,k-1]-1;]]>滤波误差方差阵Pk：Pk=[I-KkHk*]Pk,k-1;]]>式中I为单位阵；2)滤波初值确定：X^0=μX0+CX0H0T(H0CX0H0T+R0)-1(Z0-H0μZ0);]]>P0=CX0+CX0H0T(H0CX0H0T+R0)-1H0CX0;]]>式中是X0的均值，是Z0的均值，表示X0的方差。