CN108061653B - 基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法 - Google Patents

Abstract

本发明公开了一种基于谐波‑冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，通过安装在铁轨两侧的正对于列车轮对轴承的麦克风采集列车高速通过时发出的声音信号x(t)，对该检测信号的处理步骤为：(1)构建过完备参数化多普勒调制复数谐波‑冲击复合字典Datom3；(2)使用匹配追踪算法将轨边信号x(t)在构建好的过完备复数复合字典Datom3中进行稀疏分解得到投影字典Datom4及投影系数K；(3)根据轴承共振频带及麦克风到轮对轴承的几何位置关系从字典Datom4中筛选符合要求的原子组成字典Datom5并进行线性组合得到重构故障信号sig。本发明实现了与故障信号更好的时频结构上的匹配，达到更好的稀疏表示与信号重构，声源分离效果得到提升。

Description

基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声 信号分离方法

技术领域

本发明涉及高速列车轮对轴承轨边声学故障诊断技术领域，具体涉及一种基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，用于从轨边信号中分离出轮对轴承故障信号，以消除噪声、提高故障诊断的准确度。

背景技术

列车在高速运行时轮对轴承发出的声音信号中蕴含了与其健康状况密切相关的信息，轨边声学故障诊断具有非接触式监测的特点。然而，轨边信号中掺杂了来自列车其它部件发出的声音信号、轮轨接触声音和空气动力学噪声等，给有效的故障诊断带来了困难。常见的去噪方法是使用数字滤波器，例如巴特沃斯带通滤波器，但这种滤波器对于与轴承共振频带相同的带内噪声无法消除，因此去噪效果不佳，达不到理想的诊断效果。信号稀疏表示是一种通过在给定的超完备字典中用尽可能少的原子来表示信号的方法，通过原子筛选和重构能够实现更为有效的去噪，已经在音频压缩、图像处理等领域取得了广泛应用。本发明是基于稀疏分解的思想，提出了一种基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，本发明是在已有中国发明专利(专利的名称：一种用于列车轴承轨边声学故障检测的声源分离方法和申请号：CN201710555022.8)基础上的再创新，在该已有发明专利中发明人提出了使用多普勒调制谐波原子进行稀疏分解和信号重构，实现轮对轴承信号的声源分离。而本发明通过引入多普勒调制冲击原子，构成谐波-冲击多普勒调制复合字典，使得稀疏表示中使用的原子与信号内在结构更加匹配，提高了信号表示的稀疏性，从而提高了声源分离的有效性和故障诊断结果的精确度。

发明内容

本发明要解决的技术问题为：克服现有技术和方法不足，提供一种基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，具有能够消除带内噪声和克服单一字典原子对信号成分匹配不能达到最优的问题，有效提升了去噪效果和故障诊断的准确性。

本发明解决上述技术问题采用的技术方案为：一种基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，其方法是通过安装在铁轨两侧的麦克风采集列车高速通过时轮对轴承发出的故障声音信号x(t)，实现步骤如下：

步骤1：构建过完备参数化多普勒调制复数谐波-冲击复合字典矩阵Datom3：Datom3＝{S₃(i₃),i₃＝1,2,…,n3}；其中S₃(i₃)为原子集合，i₃为每个原子对应序号；

步骤2：使用匹配追踪算法将轨边信号x(t)在步骤1构建的复合字典矩阵Datom3中进行稀疏分解得到投影字典矩阵Datom4＝{S₄(j),j＝1,2…,m}(其中S₄(j)为分解信号得到的原子集合)、每个投影原子对应的投影矩阵系数C＝{C(j)^opt,j＝1,2,…,m}(其中C(j)^opt投影系数集合，j为投影原子对应的投影系数序号)、每个投影原子对应的参数集合Γ^opt＝{r^j,X₀ ^j,V^j,f_c ^j},j＝1,2,...,m(其中r^j,X₀ ^j,V^j,f_c ^j分别为匹配分解得到符合麦克风距离声源纵向距离集合、麦克风距离声源初始横向距离集合、声源初始速度集合、声源振荡频率集合)，j为投影原子所对应相关参数的序号；

步骤3：根据轴承共振频带以及麦克风与轮对轴承几何距离从步骤2中得到的投影字典矩阵Datom4中筛选符合要求的原子Datom5＝{S₅(k),k＝1,2…N}(其中S₅(k)为符合麦克风到轮对轴承的横向距离和纵向距离要求原子集合，k为筛选后的符合要求的原子对应的序号)，并进行线性叠加后得到重构信号sig：

其中N为从字典矩阵Datom4中筛选符合要求的原子Datom5原子个数，C(k)^opt为投影系数集合，S₅(k)为符合麦克风到轮对轴承的横向距离和纵向距离要求原子集合，real、image分别为集合实部和虚部，k为信号sig重构所对应字典矩阵Datom5原子的序号。

进一步地，所述步骤1中，构建复数复合Datom3的步骤如下：

2-1构建过完备谐波复数字典矩阵Datom1

(A1)设定参数集合：

其中r为麦克风距离声源纵向距离集合，r¹、r²分别为麦克风距离声源的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离声源横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离声源横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长；f_c为声源振动频率范围，f_c ¹、f_c ²声源振动频率最大、最小值，Δf为设置的声源振动频率步长。

(A2)对于步骤(A1)中Γ₁中的第i₁个参数组合：

按以下步骤生成多普勒调制谐波原子：

(A3)首先计算发声幅值序列S_e(n)：

其中f_s为轨边信号的采样频率，t_s(n)＝0,1/f_s,…,(N-1)/f_s为采样时间序列，N为采集到的轨边信号的长度，为声源振动频率集合；

(A4)计算收声时间序列t_r(n)：

其中t_s(n)为采样时间序列(发声时间序列)，c为声速；

(A5)延迟时间序列t_d(i₁)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示声源在起始点与麦克风的距离，计算公式为：

(A6)收声幅值序列S_r(n)计算：

其中M为马赫数，为声源速度；

(A7)以收声时间序列t_r(i₁)为x变量，以收声幅值序列s_r(i₁)为y变量，以延迟时间序列t_d(i₁)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子D^R(n)；

(A8)将重复步骤(A4)-(A7)得到的多普勒调频原子D^I(n)；

(A9)生成参数化多普勒调制谐波原子S₁(i₁)＝D^R(n)+j*D^I(n)；

(A10)改变i₁值，重复(A2)-(A3)，直至遍历(A1)中Γ₁中每组参数组合，最终得到过完备单位复数调频字典矩阵：

Datom1＝{S₁(i₁),i₁＝1,2…n1}。

2-2构建过完备冲击字典矩阵Datom2

(B1)设定参数集：

f_c为声源振动频率范围，f_c ¹、f_c ²声源振动频率最小、最大值，Δf_c为设置的声源振动频率步长；Ws为小波的长度范围，Ws¹、Ws²为小波长度的最小、最大值，ΔWs为设置的小波搜索的步长；a为阻尼比范围，a¹、a²为阻尼比的最小、最大值，Δa为阻尼比的变化步长；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长。

(B2)对于(B1)中Γ中的第i₂个参数：

(B3)首先计算发声幅值序列We(n)：

遍历Γ中的参数生成拉普拉斯小波序列W_e(n)：

其中：a为阻尼比，为声源振动频率集合；

(B4)设定参数集合：

其中r为麦克风距离声源纵向距离集合，r¹、r²分别为麦克风距离声源的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离声源横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离声源横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；W_e为仿真拉普拉斯小波幅值集合，W_e ¹、W_e ²声源振动幅值最小、最大值，W_e ^*是(B3)中产生的介于W_e ¹、W_e ²之间所有幅值；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长。

(B5)对于(B4)中Γ₂中的第i₂个参数组合：

按以下步骤生成多普勒冲击原子：

(B6)由(B4)可知发声幅值序列W_e(n)：

f_s为轨边信号的采样频率，a为阻尼比；

(B7)计算收声时间序列t_w(n)：

其中t_s(n)为采样时间序列(发声时间序列)，c为声速；

(B8)延迟时间序列t_d(i₂)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示声源在起始点与麦克风的距离，计算公式为：

(B9)收声幅值序列w_r(n)计算：

其中 为声源速度。

(B10)以收声时间序列t_w(n)为x变量，以收声幅值序列w_r(n)为y变量，以延迟时间序列t_d(i₂)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子S2(n)；

(B11)改变i₂值，直至遍历(B5)中Γ₂中参数，最终得到过完备单位冲击字典矩阵：

Datom2＝{S2(i₂),i₂＝1,2…n2}。

2-3构建过完备复合复数字典矩阵Datom3

将过完备谐波复数字典矩阵Datom1与过完备冲击字典矩阵Datom2组合得到过完备复合复数字典矩阵Datom3：

Datom3＝{S₃(i₃),i₃＝1,2…n3}，其中n3＝n1+n2；

进一步地，所述步骤2中，稀疏分解的步骤如下：

(C1)初始化迭代次数J＝1；

(C2)将轨边检测信号x(t)与步骤1中得到的过完备多普勒调制过完备复合字典矩阵Datom3中的每个原子进行内积运算，得到投影值数组C(i₃)：

C(i₃)＝x(t)·S₃(i₃)

(C3)计算最优投影向量：

Datom^j＝real(C(j)^opt)*real(S₃(j))+imag(C(j)^opt)*imag(S₃(j))

其中：

C(j)^opt＝max(C(i₃))

S₃(j)为C(j)^opt对应的原子；

(C4)将x(t)减去最优投影向量得到新的x(t)：

x(t)′＝x(t)-Datom^j

(C5)将J的数值加1，重复步骤(C2)-(C4)直至有至少以下条件之一满足：

以上两个公式为中止指标，norm(x(t))为每次迭代后信号的能量，J为迭代次数，σ₁和σ₂为设定的指标阈值，其中σ₁为残值能量阈值，σ₂为迭代次数阈值；

(C6)经过J次迭代得到投影字典矩阵：

Datom4＝{S₄(j),j＝1,2…m}

及每个投影原子对应的投影系数：

C＝{C(j)^opt,j＝1,2,…,m}

及每个投影原子对应的参数集合：

Γ^opt＝{r^j,X₀ ^j,V^j,f_c ^j},j＝1,2,...,m。

进一步地，所述步骤3中，得到Datom5的步骤如下：

(D1)根据列车与麦克风几何关系，确定重构几何参数，麦克风距离声源距离，根据轴承系统确定共振频带，得到筛选参数范围：

其中r为筛选后得到满足要求的麦克风距离声源纵向距离范围，dr为筛选的纵向距离步长，r_s、r_s+dr麦克风距离声源纵向距离最小值、最大值；X₀为筛选后得到满足要求的麦克风距离声源横向距离范围，dx₀为筛选的横向距离步长，X₀ ¹¹、X₀ ¹¹+dx₀为麦克风距离声源横向距离最小值、最大值；f_c为筛选后得到满足要求的声源振荡频率范围，f_c ¹¹、f_c ²²为声源振荡频率最小、最大值；

(D2)遍历Datom4中的每一个原子，如果原子对应的参数满足步骤(D1)设定的筛选参数范围，则保留，最终得到Datom5＝{S₅(i₅),i₅＝1,2…N}。

本发明与现有技术相比的优点在于：

其一，使用过完备参数化多普勒调制复合字典矩阵对轨边信号进行稀疏分解与重构，复合字典矩阵中原子参数与轨边信号内在几何参数和频率参数结构相匹配，在目标故障信号重构过程中同时使用几何参数和频率参数进行原子筛选，与传统方法相比实现了带内消噪；

其次，复合字典矩阵有效提高单一谐波字典矩阵稀疏分解效果，准确分离故障声源，提升了故障诊断的准确性。

再次，本发明可用单一麦克风实现对轨边信号去噪，有效提取故障信息并可进一步判断不同声源位置，对判断轮对轴承故障而言简单可靠具有实际意义。

附图说明

图1为本发明中的列车轴承轨边声学检测的声源分离方法流程图；

图2为轨边声学监测运动学模型；

图3为仿真分析中声源几何位置关系；

图4为轨边多声源声音采集示意图；

图5为轨边多声源声音采集俯视图；

图6(a)为小波调制的谐波信号时域图，图6(b)为其频域图；

图7(a)为多普勒调制的拉普拉斯小波信号时域图，图7(b)为其频域图；

图8(a)为仿真干扰信号时域图，图8(b)为其频域图；

图9(a)为麦克风采集的具有多普勒畸变的轨边信号，图9(b)为其频域图；

图10(a)为麦克风采集到的轨边信号x(t)与利用本发明提出的方法匹配重构得到的信号的对比图，上面为轨边信号,下面为重构轨边信号；图10(b)为两个信号的叠加；图10(c)为图10(b)局部信号放大图；图10(d)为图10(c)两个信号的对比图，上面是轨边信号，下面为重构轨边信号；

图11为利用本发明提出的方法稀疏分解计算得到不同声源位置示意图；

图12为麦克风采集到的信号x(t)经每次迭代分解后的剩余能量残值normx(t)；

图13(a)为仿真目标故障重构信号与利用本发明提出来的方法得到重构故障信号的对比图，上面为目标故障信号，下面为去噪故障信号；图13(b)为两个信号的叠加；图13(c)为两个信号的局部放大图；图13(d)为两个局部信号的对比图，上图为目标故障信号，下图为去噪故障信号。

具体实施方式

下面结合附图及实施例对本发明进行详细说明。

采用高速列车轮对轴承外圈单点故障轨边仿真声音信号进行分析处理，为了验证方法的有效性，在轨边信号同一平面内添加3个不同几何位置的噪声信号，一个故障信号。采样频率4KHz。

图1为发明中基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法流程图。使用安装在铁轨两侧的麦克风采集列车高速通过时轮对轴承发出的故障声音信号，作为检测信号x(t)，对该检测信号的处理步骤为：

(1)构建过完备参数化多普勒复数复合字典Datom3：

Datom3＝{S₃(i₃),i₃＝1,2,…,n3}

(2)将轨边检测信号x(t)在步骤(1)中构建好的过完备复合字典矩阵Datom3中进行稀疏分解得到投影字典矩阵：

Datom4＝{S₄(j),j＝1,2…m}

及每个投影原子对应的投影系数：

C＝{C(j)^opt,i＝1,2,…,m}

及每个投影原子对应的参数集合：

Γ^opt＝{r^j,X₀ ^j,V^j,f_c ^j},j＝1,2,...,m

(3)根据麦克风到轮对轴承的横向距离XX02和纵向距离RR02从步骤(2)中得到的投影字典矩阵Datom4中筛选符合要求的原子：

Datom5＝{S₅(k),k＝1,2…N}

并进行线性叠加后得到重构信号sig：

进一步地，所述步骤(1)中，构建复数复合Datom3的步骤如下：

2-1构建过完备谐波复数字典矩阵Datom1

(A1)设定参数集合：

其中r为麦克风距离声源纵向距离集合，r¹、r²分别为麦克风距离声源的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离声源横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离声源横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长；f_c为声源振动频率范围，f_c ¹、f_c ²声源振动频率最大、最小值，Δf为设置的声源振动频率步长。

(A2)对于步骤(A1)中Γ₁中的第i₁个参数组合：

按以下步骤生成多普勒调制谐波原子：

(A3)首先计算发声幅值序列S_e(n)：

其中f_s为轨边信号的采样频率，t_s(n)＝0,1/f_s,…,(N-1)/f_s为采样时间序列，N为采集到的轨边信号的长度，为声源振动频率集合；

(A4)计算收声时间序列t_r(n)：

其中t_s(n)为采样时间序列(发声时间序列)，c为声速；

(A5)延迟时间序列t_d(i₁)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示声源在起始点与麦克风的距离，计算公式为：

(A6)收声幅值序列S_r(n)计算：

其中M为马赫数，为声源速度；

(A7)以收声时间序列t_r(i₁)为x变量，以收声幅值序列s_r(i₁)为y变量，以延迟时间序列t_d(i₁)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子D^R(n)；

(A8)将重复步骤(A4)-(A7)得到的多普勒调频原子D^I(n)；

(A9)生成参数化多普勒调制谐波原子S₁(i₁)＝D^R(n)+j*D^I(n)；

(A10)改变i₁值，重复(A2)-(A3)，直至遍历(A1)中Γ₁中每组参数组合，最终得到过完备单位复数调频字典矩阵：

Datom1＝{S₁(i₁),i₁＝1,2…n1}。

2-2构建过完备冲击字典矩阵Datom2

(B1)设定参数集：

f_c为声源振动频率范围，f_c ¹、f_c ²声源振动频率最小、最大值，Δf_c为设置的声源振动频率步长；Ws为小波的长度范围，Ws¹、Ws²为小波长度的最小、最大值，ΔWs为设置的小波搜索的步长；a为阻尼比范围，a¹、a²为阻尼比的最小、最大值，Δa为阻尼比的变化步长；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长。

(B2)对于(B1)中Γ中的第i₂个参数：

(B3)首先计算发声幅值序列We(n)：

遍历Γ中的参数生成拉普拉斯小波序列W_e(n)：

其中：a为阻尼比，为声源振动频率集合；

(B4)设定参数集合：

其中r为麦克风距离声源纵向距离集合，r¹、r²分别为麦克风距离声源的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离声源横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离声源横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；W_e为仿真拉普拉斯小波幅值集合，W_e ¹、W_e ²声源振动幅值最小、最大值，W_e ^*是(B3)中产生的介于W_e ¹、W_e ²之间所有幅值；V为初始速度集合，V¹、V²分别为声源速度最小、最大值，ΔV为设置的声源速度变化步长。

(B5)对于(B4)中Γ₂中的第i₂个参数组合：

按以下步骤生成多普勒冲击原子：

(B6)由(B4)可知发声幅值序列W_e(n)：

f_s为轨边信号的采样频率，a为阻尼比；

(B7)计算收声时间序列t_w(n)：

其中t_s(n)为采样时间序列(发声时间序列)，c为声速；

(B8)延迟时间序列t_d(i₂)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示声源在起始点与麦克风的距离，计算公式为：

(B9)收声幅值序列w_r(n)计算：

其中 为声源速度。

(B10)以收声时间序列t_w(n)为x变量，以收声幅值序列w_r(n)为y变量，以延迟时间序列t_d(i₂)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子S2(n)；

(B11)改变i₂值，直至遍历(B5)中Γ₂中参数，最终得到过完备单位冲击字典矩阵：

Datom2＝{S2(i₂),i₂＝1,2…n2}。

2-3构建过完备复合复数字典矩阵Datom3

将过完备谐波复数字典矩阵Datom1与过完备冲击字典矩阵Datom2组合得到过完备复合复数字典矩阵Datom3：

Datom3＝{S₃(i₃),i₃＝1,2…n3}，其中n3＝n1+n2；

进一步地，所述步骤2中，稀疏分解的步骤如下：

(C1)初始化迭代次数J＝1；

(C2)将轨边检测信号x(t)与步骤1中得到的过完备多普勒调制过完备复合字典矩阵Datom3中的每个原子进行内积运算，得到投影值数组C(i₃)：

C(i₃)＝x(t)·S₃(i₃)

(C3)计算最优投影向量：

Datom^j＝real(C(j)^opt)*real(S₃(j))+imag(C(j)^opt)*imag(S₃(j))

其中：

C(j)^opt＝max(|C(i₃)|)

S₃(j)为C(j)^opt对应的原子；

(C4)将x(t)减去最优投影向量得到新的x(t)：

x(t)′＝x(t)-Datom^j

(C5)将J的数值加1，重复步骤(C2)-(C4)直至有至少以下条件之一满足：

以上两个公式为中止指标，norm(x(t))为每次迭代后信号的能量，J为迭代次数，σ₁和σ₂为设定的指标阈值，其中σ₁为残值能量阈值，σ₂为迭代次数阈值；

(C6)经过J次迭代得到投影字典矩阵：

Datom4＝{S₄(j),j＝1,2…m}

及每个投影原子对应的投影系数：

C＝{C(j)^opt,j＝1,2,…,m}

及每个投影原子对应的参数集合：

Γ^opt＝{r^j,X₀ ^j,V^j,f_c ^j},j＝1,2,...,m。

进一步地，所述步骤3中，得到Datom5的步骤如下：

(D1)根据列车与麦克风几何关系，确定重构几何参数,麦克风距离声源距离，根据轴承系统确定共振频带，得到筛选参数范围：

其中r为筛选后得到满足要求的麦克风距离声源纵向距离范围，dr为筛选的纵向距离步长，r_s、r_s+dr麦克风距离声源纵向距离最小值、最大值；X₀为筛选后得到满足要求的麦克风距离声源横向距离范围，dx₀为筛选的横向距离步长，X₀ ¹¹、X₀ ¹¹+dx₀为麦克风距离声源横向距离最小值、最大值；f_c为筛选后得到满足要求的声源振荡频率范围，f_c ¹¹、f_c ²²为声源振荡频率最小、最大值；

(D2)遍历Datom4中的每一个原子，如果原子对应的参数满足步骤(D1)设定的筛选参数范围，则保留，最终得到Datom5＝{S₅(i₅),i₅＝1,2…N}。

图2轨边声学监测运动学模型。图3为仿真示意中声源几何位置关系，麦克风(M)距离轨边纵向距离1m，N1、N2、N3为干扰信号距麦克风纵向距离分别为1.8m、1.6m、1.2m，横向距离2.4m、1.6m、1.2m，O为仿真故障信号距麦克风纵向距离1.4m，横向距离2.0m,N1、N2、N3、O四点在同一平面内(图平行四边形示)。图4轨边多声源声音采集示意图，dx、dy、dz为扬声器之间距离，O点为目标声源，其余点为实验干扰信号，O′为采集轨边信号的麦克风，设定最外带有O点目标声源的平面为研究对象，移动方向自左向右如图4所示。图5为图4轨边多声源声音采集示意图的俯视图，O点目标声源，O′为采集轨边信号的麦克风，dx₀为麦克风横向采集信号搜索距离步长、dr为麦克风纵向采集信号搜索距离步长。

图6(a)仿真冲击信号时域波形，图6(b)为其频域图，从中可以看出谐波信号的中心频率为1000Hz。图7(a)冲击信号经过多普勒调制后的畸变信号(目标故障信号sig)时域图，图7(b)为其频域图。图8(a)为仿真干扰噪信号noise_sig时域图，干扰噪信号noise_sig是由两个不同频率三个不同空间位置的谐波信号经多普勒畸变调制生成，其与仿真目标故障信号对应在同一平面的不同位置，图8(b)为其频域图。图9(a)为轨边信号，即麦克风采集到的信号，包括仿真干扰信号和去噪信号sig，图9(b)为其频域图，从中可以看出信号有三个频率成分。图10(a)利用本发明提出来的方法得到的重构轨边信号和原轨边信号的对比；图10(b)为轨边信号x(t)以及轨边信号x(t)在利用本发明提出来的复合字典Datom3的分解得Datom4线性组合得到重构轨边信号的叠加；图10(c)为取图10(b)270-330点的匹配轨边信号与原始信号的局部放大图，可以看出仿真信号得到很好的匹配分解；图10(d)为图10(c)重叠信号拆分对比图。

图11为利用本发明提出来的方法计算四个声源在空间的具体位置。本发明提出来的构建复合字典矩阵声源距离麦克风纵向距离范围是R＝[1,1.9]，搜索距离为0.1，搜索长度为10；横向搜索范围X＝[1,3]，搜索距离为0.2，搜索长度为11。干扰信号N3(R₃，X₂)、N2位于(R₇，X₄)、N1位于(R₉，X₈)、O位于(R₅，X₆)(R₅表示表示纵向第五个距离1.4m，X₆表示横向第六个距离2.0m，剩余参数以此类推)，采用本专利提出来的基于复合字典的稀疏分解得到的声源位置如图11所示，此结果与事先设置的声源位置参数保持一致，验证了本专利提出方法的可行性。图12为麦克风采集到的信号x(t)经每次迭代分解后的剩余能量norm(x)。图13(a)为目标故障信号与去噪故障信号对比；图13(b)为仿真目标故障信号与利用本发明提出来的方法得到去噪信号的叠加；图13(c)为图13(b)取280-360点的两个信号的局部放大图；图13(d)是图13(c)重叠信号拆分对比图，可以看出利用本发明提出来的方法仿真故障信号得到很好的重构，达到良好的较好的去噪效果。

提供上述实施例仅仅是为了描述本发明的目的，并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (1)

1.一种基于谐波-冲击多普勒调制复合字典的列车轮对轴承轨边声信号分离方法，其特征在于：该方法是通过安装在铁轨两侧的麦克风采集列车高速通过时轮对轴承发出的轨边信号x(t)，对该信号的处理步骤如下：

步骤1：构建过完备参数化多普勒调制复数谐波-冲击复合字典矩阵Datom3：

Datom3＝{S₃(i₃),i₃＝1,2,…,n3}；其中S₃(i₃)为原子集合，i₃为每个原子对应序号；

步骤2：使用匹配追踪算法将轨边信号x(t)在步骤1构建的复合字典矩阵Datom3中进行匹配稀疏分解得到投影字典矩阵Datom4＝{S₄(j),j＝1,2…m}，其中S₄(j)为分解信号得到的投影原子集合、每个投影原子对应的投影矩阵系数C＝{C(j)^opt,j＝1,2,…,m}，其中C(j)^opt为投影系数集合，j为投影原子对应的投影系数序号、每个投影原子对应的参数集合其中r^j,X₀ ^j,V^j,f_c ^j分别为匹配稀疏分解得到符合麦克风距离轮对轴承纵向距离集合、麦克风距离轮对轴承初始横向距离集合、轮对轴承初始速度集合、轮对轴承振动频率集合；

步骤3：根据轴承共振频带以及麦克风与轮对轴承几何距离从步骤2中得到的投影字典矩阵Datom4中筛选符合要求的投影原子Datom5＝{S₅(k),k＝1,2…L}，其中S₅(k)为符合麦克风到轮对轴承的横向距离和纵向距离要求的投影原子集合，k为筛选后的符合要求的投影原子对应的序号，并进行线性叠加后得到重构信号sig：

其中L为从投影字典矩阵Datom4中筛选符合要求的投影原子Datom5原子个数，C(k)^opt为投影系数集合，real、image分别为集合实部和虚部，k为信号sig重构所对应符合要求的投影原子Datom5原子的序号；

所述步骤1中，构建复合字典矩阵Datom3的步骤如下：

2-1 构建过完备谐波复数字典矩阵Datom1

(A1)设定参数集合：

其中r为麦克风距离轮对轴承纵向距离集合，r¹、r²分别为麦克风距离轮对轴承的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离轮对轴承横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离轮对轴承横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；V为初始速度集合，V¹、V²分别为轮对轴承初始速度最小、最大值，ΔV为设置的轮对轴承初始速度变化步长；f_c为轮对轴承振动频率范围，f_c ¹、f_c ²轮对轴承振动频率最小、最大值，Δf为设置的轮对轴承振动频率步长；

(A2)对于步骤(A1)中Γ₁中的第i₁个参数组合：

按以下步骤生成多普勒调制谐波原子：

(A3)首先计算发声幅值序列S_e(n)：

其中f_s为轨边信号的采样频率，t_s(n)＝0,1/f_s,…,(N-1)/f_s为采样时间序列，N为采集到的轨边信号的长度，为轮对轴承振动频率集合；

(A4)计算收声时间序列t_r(n)：

其中t_s(n)为采样时间序列，c为声速；

(A5)延迟时间序列t_d(n)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示轮对轴承在起始点与麦克风的距离，计算公式为：

(A6)收声幅值序列S_r(n)计算：

其中，M为马赫数，为轮对轴承初始速度；

(A7)以收声时间序列t_r(n)为x变量，以收声幅值序列S_r(n)为y变量，以延迟时间序列t_d(n)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子D^R(n)；

(A8)将重复步骤(A4)-(A7)得到的多普勒调频原子D^I(n)；

(A9)生成参数化多普勒调制谐波原子S₁(i₁)＝D^R(n)+j*D^I(n)；

(A10)改变i₁值，在步骤(A9)所示谐波信号的基础上，重复步骤(A2)-(A9)，直至遍历步骤(A1)中Γ₁每组参数组合，最终得到过完备谐波复数字典矩阵：

Datom1＝{S₁(i₁),i₁＝1,2…n1}

2-2构建过完备冲击字典矩阵Datom2

(B1)设定参数集：

f_c为轮对轴承振动频率范围，f_c ¹、f_c ²为轮对轴承振动频率最小、最大值，Δf_c为设置的轮对轴承振动频率步长；Ws为小波的长度范围，Ws¹、Ws²为小波长度的最小、最大值，ΔWs为设置的小波搜索的步长；a为阻尼比范围，a¹、a²为阻尼比的最小、最大值，Δa为阻尼比的变化步长；V为初始速度集合，V¹、V²分别为轮对轴承初始速度最小、最大值，ΔV为设置的轮对轴承初始速度变化步长；

(B2)对于(B1)中Γ中的第i₂个参数：

(B3)首先计算发声幅值序列We(n)：

遍历Γ中的参数生成拉普拉斯小波序列W_e(n)：

其中：a为阻尼比，为轮对轴承振动频率集合；

(B4)设定参数集合：

其中r为麦克风距离轮对轴承纵向距离集合，r¹、r²分别为麦克风距离轮对轴承的最近、最远距离，Δr为设置的纵向距离变化步长；X₀为麦克风距离轮对轴承横向距离集合，X₀ ¹、X₀ ²分别为麦克风距离轮对轴承横向距离最近、最远距离，ΔX₀为设置的横向距离变化步长；W_e为仿真拉普拉斯小波幅值集合，W_e ¹、W_e ²为拉普拉斯小波幅值最小、最大值，W_e ^*是(B3)中产生的介于W_e ¹、W_e ²之间所有幅值；V为初始速度集合，V¹、V²分别为轮对轴承初始速度最小、最大值，ΔV为设置的轮对轴承初始速度变化步长；

(B5)对于步骤(B4)中Γ₂中的第i₂个参数组合：

按以下步骤生成多普勒冲击原子：

(B6)由步骤(B3)可知发声幅值序列W_e(n)：

f_s为轨边信号的采样频率， 为阻尼比；

(B7)计算收声时间序列t_w(n)：

其中t_s(n)为采样时间序列，c为声速；

(B8)延迟时间序列t_d(n)计算，延迟时间序列t_d(n)为最终得到时间序列，其值t_d(n)＝t_s(n)+R(0)/c，其中R(0)表示轮对轴承在起始点与麦克风的距离，计算公式为：

(B9)收声幅值序列w_r(n)计算：

其中 为轮对轴承初始速度；

(B10)以收声时间序列t_w(n)为x变量，以收声幅值序列w_r(n)为y变量，以延迟时间序列t_d(n)为插值x变量，执行三次样条插值重采样处理，并进行能量归一化得到多普勒调频原子S2(n)；

(B11)改变i₂值，直至遍历步骤(B5)中Γ₂中参数，最终得到过完备冲击字典矩阵：

Datom2＝{S2(i₂),i₂＝1,2…n2}

2-3 构建过完备参数化多普勒调制复数谐波-冲击复合字典矩阵Datom3

将过完备谐波复数字典矩阵Datom1与过完备冲击字典矩阵Datom2组合得到过完备参数化多普勒调制复数谐波-冲击复合字典矩阵Datom3：

Datom3＝{S₃(i₃),i₃＝1,2…n3}，其中n3＝n1+n2；

所述步骤2中，匹配稀疏分解的步骤如下：

(C1)初始化迭代次数J＝1；

(C2)将轨边信号x(t)与步骤1中得到的过完备参数化多普勒调制复数谐波-冲击复合字典矩阵Datom3中的每个原子进行内积运算，得到投影值数组C(i₃)：

C(i₃)＝x(t)·S₃(i₃)

(C3)计算最优投影向量：

Datom^J＝real(C(J)^opt)*real(S₃(J))+imag(C(J)^opt)*imag(S₃(J))

其中：

C(J)^opt＝max(|C(i₃)|)

S₃(J)为C(J)^opt对应的原子；

(C4)将x(t)减去最优投影向量得到新的x(t)′：

x(t)′＝x(t)-Datom^J

(C5)将J的数值加1，重复步骤(C2)-(C4)直至有至少以下条件之一满足：

以上两个公式为中止指标，norm(x(t))为每次迭代后信号的能量，J为迭代次数，σ₁和σ₂为设定的指标阈值，其中σ₁为残值能量阈值，σ₂为迭代次数阈值；

(C6)经过m次迭代得到投影字典矩阵：

Datom4＝{S₄(j),j＝1,2…,m}

及每个投影原子对应的投影系数：

C＝{C(j)^opt,j＝1,2,…,m}

及每个投影原子对应的参数集合：

Γ^opt＝{r^j,X₀ ^j,V^j,f_c ^j},j＝1,2,...,m；

所述步骤3中，得到Datom5的步骤如下：

(D1)根据列车与麦克风几何关系，确定重构几何参数，麦克风距离轮对轴承距离，根据轮对轴承确定共振频带，得到筛选参数范围：

其中r为筛选后得到满足要求的麦克风距离轮对轴承纵向距离范围，dr为筛选的纵向距离步长，r_s、r_s+dr为麦克风距离轮对轴承纵向距离最小值、最大值；X₀为筛选后得到满足要求的麦克风距离轮对轴承横向距离范围，dx₀为筛选的横向距离步长，X₀ ¹¹、X₀ ¹¹+dx₀为麦克风距离轮对轴承横向距离最小值、最大值；f_c为筛选后得到满足要求的轮对轴承振动频率范围，f_c ¹¹、f_c ²²为轮对轴承振动频率最小、最大值；

(D2)遍历Datom4中的每一个原子，如果原子对应的参数满足步骤(D1)设定的筛选参数范围，则保留，最终得到Datom5＝{S₅(i₅),i₅＝1,2…N}。