CN116720059A - CEEMDAN-Grey-SSA弱信号提取方法 - Google Patents

Abstract

本发明公开了一种CEEMDAN‑Grey‑SSA弱信号提取方法，首先对获取得到的车体偏移数据进行完全自适应噪声集合经验模态分解CEEMDAN分解，得到K个IMF分量信号，将IMF₁(t)视为主要包含噪声的分量，因此计算剩余IMF分量与IMF₁(t)的灰色关联度表征IMF分量的噪声水平，并进行降序排列；搭建灰色模型，然后进行软阈值降噪处理，然后与有效信号占主导的IMF分量进行重构；基于奇异谱分析SSA进行二次降噪，得到采集的低信噪比列车偏移量中由于轨道水平不平顺导致的车体横移量。本发明通过对信号进行时频分析，可对有效信号进行准确提取，以便在低信噪比的情况下进行弱信号的提取。

Description

CEEMDAN-Grey-SSA弱信号提取方法

技术领域

本发明属于轨道检测技术领域，具体涉及一种CEEMDAN-Grey-SSA弱信号提取方法。

背景技术

轨道作为铁路的重要组成部分之一，对列车的安全舒适运行具有重要影响。因此需要保持良好的轨道几何状态。轨道不平顺包括高低、轨向、轨距、水平和扭曲不平顺等，都会对乘车质量和交通安全产生重大影响。由于目前主要依靠轨检车及便携检查仪进行轨道不平顺检测，考虑到其检测周期与检测成本之间的矛盾，因此依靠在役列车结合其车体运行姿态数据进行轨道不平顺的快速检测具有重大意义。

以水平不平顺为例，列车行驶过程中，受水平不平顺影响，车身姿态会发生轻微变化，产生车体的动态偏移量。然而列车在运行过程中不仅会受到传感器自身测量误差的影响，还会受到环境中各类噪声的影响。而轨道不平顺激扰对车体姿态的影响极其微弱，因此会出现有效信息被噪声淹没的状况，使车体动态偏移量呈现低信噪比的特点，导致轨道不平顺检测精度受到严重影响。因此在进行轨道不平顺检测之前，首先对车体运行姿态数据进行有效信号提取。综上所述，对有效信号的准确提取就显得尤为重要。

结合噪声信息远远大于有效信号的特点，针对轨道不平顺检测研究所面临的问题，提出一种CEEMDAN-Grey-SSA弱信号提取方法，该方法不仅能够提高轨道不平顺检测的精度，而且对于该学科和行业的进步起着关键性的作用。

发明内容

本发明的目的是提供一种CEEMDAN-Grey-SSA弱信号提取方法，通过对信号进行时频分析，可对有效信号进行准确提取，以便在低信噪比的情况下进行弱信号的提取。

本发明所采用的技术方案是，CEEMDAN-Grey-SSA弱信号提取方法，具体按照以下步骤实施：

步骤1、对获取得到的车体偏移数据X(t)，t＝1,2,...,T，T表示采样点的总数，进行完全自适应噪声集合经验模态分解CEEMDAN分解，得到K个IMF分量信号，依次命名为IMF₁(t),IMF₂(t),…,IMF_k(t),…,IMF_K(t)，其中k＝1,2,...,K，K为CEEMDAN分解得到IMF分量的总数；

步骤2、针对步骤1得到K组IMF分量，将IMF₁(t)视为主要包含噪声的分量，因此计算剩余IMF分量与IMF₁(t)的灰色关联度表征IMF分量的噪声水平，并进行降序排列；

步骤3、搭建灰色模型，根据步骤2中按噪声水平降序排序的IMF分量，自适应筛选出噪声信号占主导的IMF分量和有效信号占主导的IMF分量；

步骤4、根据步骤3选择的噪声信号占主导的IMF分量，进行软阈值降噪处理，然后与有效信号占主导的IMF分量进行重构；

步骤5、根据步骤4重构后的信号，基于奇异谱分析SSA进行二次降噪，得到采集的低信噪比列车偏移量中由于轨道水平不平顺导致的车体横移量。

本发明的特点还在于，

步骤1具体按照以下步骤实施：

步骤1.1、随机生成N组白噪声序列wⁱ(t)，i＝1,2,…,N，并将白噪声序列加入原始的车体偏移量x(t)中，获得N组含噪信号，则第i组含噪信号xⁱ(t)表示如下：

xⁱ(t)＝x(t)+λ_iwⁱ(t) (1)

式中，λ_i为第i组高斯白噪声的权值系数；

步骤1.2、对于步骤1.1产生的含噪信号，利用EMD算法分别进行分解，将每组含噪信号通过EMD算法分解，得到N组一阶IMF分量，将N组一阶IMF分量求平均，即得到x(t)的一阶IMF分量，记为IMF₁(t)，计算如下所示：

式中IMF₁ ⁱ(t)表示第i组一阶IMF分量；

步骤1.3、根据步骤1.2得到的x(t)的一阶IMF分量，计算一阶IMF分量的残余项r₁(t)，如下所示：

r₁(t)＝x(t)-IMF₁(t) (3)

步骤1.4、根据步骤1.3中的残余项r₁(t)，给残余项r₁(t)中继续添加白噪声得到r₁(t)+λ₁E₁(wⁱ(t))，获得N组信号，通过EMD算法分解求平均，同步骤1.1～1.2，获得x(t)的二阶分量IMF₂(t)，如下所示：

式中E₁(·)表示经EMD分解后的第1个IMF分量，λ₁E₁(wⁱ(t))是添加白噪声；

步骤1.5、当k＝2,3,…,K，计算第k-1个残余分量r_k-1(t)，如式所示：

r_k-1(t)＝r_k-2(t)-IMF_k-1(t) (5)

式中r_k-2(t)是第k-2个残余分量；IMF_k-1(t)是第k-1阶IMF分量；

步骤1.6、计算r_k-1(t)+λ_k-1E_k-1(wⁱ(t))中的IMF分量并平均，获得IMF_k(t)，如下所示：

E_k-1(·)表示经EMD分解后的第k-1个IMF分量；r_k-1(t)是第k-1个残余分量；λ_k-1E_k-1(wⁱ(t))是添加白噪声；

步骤1.7、循环执行步骤1.5～步骤1.6，直到残差最多有一个极值无法分解，最终得到第K阶分量IMF_K(t)与残余项r_K，如下所示：

则原始信号x(t)表示为：

步骤1中完全自适应噪声集合经验模态分解CEEMDAN算法中，每个分解阶段均通过噪声系数λ选择信噪比，噪声系数λ为0.2。

步骤2具体如下：

步骤2.1、根据步骤1中分解的各个IMF分量，将其分别表征为序列：

IMF_k(t)＝[IMF_k(1),IMF_k(2),...,IMF_k(j),...,IMF_k(T)] (9)

其中:IMF_k是第k个IMF分量的原始序列；IMF_k(j)为IMF_k的j个元素，j＝1,2,...,T；T是序列的长度；

步骤2.2、根据步骤2.1得到的各个分量的原始序列，对其进行极性统一和平均，获得极性统一序列IMF_k ⁽⁰⁾和平均序列IMF_k ⁽¹⁾：

IMF_k ⁽⁰⁾(j)＝IMF_k(j)+|min(IMF_k(j)) (10)

式中，IMF_k ⁽⁰⁾(j)和IMF_k ⁽¹⁾(j)分别代表IMF_k ⁽⁰⁾和IMF_k ⁽¹⁾序列中的第j个元素；

步骤2.3、计算IMF₁与其余IMF分量之间的灰色关联系数序列，公式如下所示：

式中，γ_k为第k个IMF分量与IMF₁之间的灰色关联系数序列；γ_k(j)表

示γ_k的第j个元素；ρ为分辨系数；

步骤2.4、计算IMF₁与其余IMF分量之间的灰色关联度ξ_k：

式中，w为权重系数；ξ_k为第k个IMF分量与IMF₁之间的灰色关联度，取值范围为ξ_k∈[0,1]，灰色关联度值越大，代表IMF分量噪声水平越高；

步骤2.5、根据步骤2.4得到的各IMF分量灰色关联度，对各IMF分量灰色关联度进行降序排列。

步骤2.3中分辨系数ρ取值为0.5；所述步骤2.4中权重系数w取直为1。

步骤3具体如下：

步骤3.1、通过去趋势波动分析DFA算法得到原始车体偏移量x(t)的参数α(t)；

步骤3.2、建立大量的样本数据，利用原始车体偏移量序列的α(t)和η搭建灰色模型；

步骤3.3、根据建立的灰色模型及x(t)的α(t)获取到噪声信号占主导分量的比例。

步骤3.2中搭建灰色模型的具体步骤如下：

步骤3.2.1、分别获取原始车体偏移量序列的去趋势波动分析DFA参数α(t)的原始序列和η(t)的原始序列：

α(t)＝(α(1),α(2),...α(t),...,α(T)) (14)

η(t)＝(η(1),η(2),...η(t),...,η(T)) (15)

步骤3.2.2、对步骤3.2.1得到的α(t)和η(t)原始序列进行归一化，得到归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)：

α⁽¹⁾是灰色模型输入原始序列的第1个元素；η⁽¹⁾是灰色模型输出原始序列的第1个元素。

步骤3.2.3、对步骤3.2.2得到的归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)进行累计生成操作，获得生成序列α⁽¹⁾(t)和η⁽¹⁾(t)：

α⁽⁰⁾(q)表示归一序列α⁽⁰⁾的第q个元素；η⁽⁰⁾(q)表示归一序列η⁽⁰⁾的第q个元素；

步骤3.2.4、计算通过生成序列η⁽¹⁾(t)计算相邻平均序列z(t)，公式如下：

步骤3.2.5、构建数据矩阵B和数据向量y_T：

式中z(2)表示相邻平均序列z(t)的第2个元素；z(3)表示相邻平均序列z(t)的第3个元素；z(T)表示相邻平均序列z(t)的第T个元素；α⁽¹⁾(2)表示生成序列α⁽¹⁾(t)的第2个元素；α⁽¹⁾(3)表示生成序列α⁽¹⁾(t)的第3个元素；α⁽¹⁾(T)表示生成序列α⁽¹⁾(t)的第T个元素；η⁽⁰⁾(2)表示原始序列η⁽⁰⁾(t)的第2个元素；η⁽⁰⁾(3)表示原始序列η⁽⁰⁾(t)的第3个元素；η⁽⁰⁾(T)表示原始序列η⁽⁰⁾(t)的第T个元素；

步骤3.2.6、根据步骤3.2.5构建数据矩阵B和数据向量y_T，计算参数向量P_T，P_T表示如下：

式中a，b为灰色模型参数；B^T表示矩阵B的转置矩阵。

步骤3.2.7、建立归一化序列模型：

式中为归一化序列的预测值；/>

步骤3.2.8、最终样本数据得到原始序列的模型为：

式中α为输入信号的DFA计算结果；为输入信号中以噪声为主的IMF分量的比例；η(t)为样本数据，其中t由α决定，t＝argminα-α(t)-1；

步骤3.2.9、利用步骤3.2.8中计算出的噪声为主的IMF分量比例计算出噪声占主导的IMF分量的个数：

M₀表示噪声占主导的IMF分量的个数，K为CEEMDAN分解得到IMF分量的总数。

从而得到CEEMDAN分解的IMF分量中，噪声占主导的IMF分量。

步骤4具体如下：

步骤4.1、使用软阈值的方法对步骤3中最后得到的噪声占主导的IMF分量处理，具体如下：

式中表示降噪后的IMF分量；sign(·)返回值为整型变量的函数，若返回值大于0，则sign(·)返回为1，若返回值等于0，则sign()返回为0，若返回值小于0，则sign(·)返回为-1，T₀为根据信号特征获得的阈值，如下式(26)所示：

步骤4.2、最后，得到车体横向偏移量的重构信号如下式(27)所示：

式中表示去噪后的噪声主导信号分量；IMF_s1(t)表示有效信号主导的信号分量；r(t)表示残余分量；

最终得到经过完全自适应噪声集合经验模态分解-灰色模型CEEMDAN-Grey降噪后的车体横向偏移量

步骤5具体如下：

步骤5.1、相空间重构：对步骤4中最后得到的CEEMDAN-Grey降噪后的车体横向偏移量通过选择合适的滑窗长度L，以滑动窗口的方式将序列截取为以L₁为大小的序列片段L个，其中L₁＝T-L+1，将序列片段组合，转换为二维轨迹矩阵Y，如下式(28)所示：

y₁表示序列的第一个元素；y₂表示/>序列的第二个元素；y₃表示/>序列的第三个元素；y_L1表示/>序列的第L₁个元素；y_L1+1表示/>序列的第L₁+1个元素；y_L表示/>序列的第L个元素；y_L+1表示/>序列的第L+1个元素；y_T表示/>序列的第T个元素；

步骤5.2、奇异值分解，定义矩阵S＝YY^T，其中Y^T为二维轨迹矩阵的转置，对矩阵S进行奇异值分解，获得矩阵X的特征值λ和对应的特征向量矩阵U、V，则矩阵S表示为：

S＝UΛV^T (29)

式中U、V分别为正交矩阵，是矩阵Y的左右奇异矩阵；Λ为特征值对角矩阵。

对特征值进行降序排列，表示成向量形式如下：

λ＝[λ₁，λ₂，...，λ_L],|λ₁|≥|λ₂|≥…≥|λ_L|≥0

L表示矩阵X非零奇异点的个数，同时，矩阵S的奇异值向量表示成如下形式：

矩阵S的奇异值向量σ也是轨迹矩阵Y的奇异谱。

轨迹矩阵Y为：

式中σ_r表示第r个奇异值；U_r、V_r表示第r个奇异值对应的特征向量矩阵；

步骤5.3、奇异值分组：基于奇异值差分谱法寻找奇异值突变点，当两个奇异值之间存在较大的突变，代表信号中的主导分量发生变化，降序后的奇异值σ＝[σ₁,σ₂,...,σ_L]，令b_r＝σ_r-σ_r+1，则奇异值的差分值为b＝[b₁,b₂,...b_L-1]；当差分值序列在r处出现波峰时，代表前r个奇异值对应有效信号分量，则矩阵S表达为公式(32)

S＝UΛV^T＝U_w1Λ_w1V_w1 ^T+U_w2Λ_w2V_w2 ^T＝W₁+W₂ (32)

式中，U_w1表示矩阵Y的左奇异矩阵的有效信号矩阵；Λ_w1表示矩阵Y的特征值对角矩阵的有效信号矩阵；V_w1 ^T表示矩阵Y的右奇异矩阵的有效信号矩阵；U_w2表示矩阵Y的左奇异矩阵的噪声信号矩阵；Λ_w2表示矩阵Y的特征值对角矩阵的噪声信号矩阵；V_w2 ^T表示矩阵Y的右奇异矩阵的噪声信号矩阵；W₁表示有效信号矩阵；W₂表示噪声数据矩阵；

步骤5.4、对角平均化：步骤5.3通过奇异值差分谱将轨迹矩阵分为有效信号矩阵和噪声信号矩阵，采用对角平均化的方法使有效信号矩阵W₁重构为长度为T的一维车体动态偏移量序列数据X(t′)，其中有效信号矩阵W₁为L

行L₁列的矩阵：

式中表示矩阵W₁的第d行，第t′-d+1列的元素，L^*＝min(L,L₁)，T＝L+L₁-1；

最终得到经过完全自适应噪声集合经验模态分解-灰色模型-奇异谱分析CEEMDAN-Grey-SSA算法有效提取的车体横向偏移量X(t′)。

本发明的有益效果是，一种CEEMDAN-Grey-SSA弱信号提取方法，通过CEEMDAN方法将噪声信号分解成一系列本征模态函数，利用灰色关联度与灰色模型实现噪声主导分量的自适应筛选，然后对噪声占主导的IMF分量进行软阈值处理，并与有效信号主导的IMF分量进行重构。然后针对重构信号低频区域残留的噪声信息采用SSA二次降噪。通过实验仿真验证了所提方法在弱信号提取方面的准确性。

附图说明

图1为本发明一种CEEMDAN-Grey-SSA弱信号提取方法总体流程图；

图2为本发明一种CEEMDAN-Grey-SSA弱信号提取方法中CEEMDAN分解流程图；

图3(a)是传统SSA的去噪结果图；图3(b)是CEEMDAN-Grey去噪结果图；图3(c)是CEEMDAN-Grey结合SSA去噪结果图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于CEEMDAN-Grey结合SSA的双层自适应弱信号提取方法，流程图如图1所示。首先使用CEEMDAN算法对原始信号进行分解，并通过灰色关联法对IMF分量的噪声水平进行估计，然后利用灰色模型自适应筛选噪声主导分量，最后对其进行噪声处理并重构，从而实现信号高频噪声的有效去除；针对CEEMDAN-Grey重构信号低频去噪不完全问题，接着使用SSA实现重构信号的二次降噪，压制重构信号中剩余的随机噪声，完成信号低频噪声的有效去除，实现低信噪比情况下的弱信号提取。

本发明CEEMDAN-Grey-SSA弱信号提取方法，具体按照以下步骤实施：

步骤1、对获取得到的车体横移偏移数据X(t)，t＝1,2,...,T，T表示采样点的总数，进行完全自适应噪声集合经验模态分解CEEMDAN分解，得到K个IMF分量信号，依次命名为IMF₁(t),IMF₂(t),…,IMF_k(t),…,IMF_K(t)，其中k＝1,2,...,K，K为CEEMDAN分解得到IMF分量的总数；

结合图1、图2，步骤1具体按照以下步骤实施：

步骤1.1、随机生成N组白噪声序列wⁱ(t)，i＝1,2,…,N，并将白噪声序列加入原始的车体偏移量x(t)中，获得N组含噪信号，则第i组含噪信号xⁱ(t)表示如下：

xⁱ(t)＝x(t)+λ_iwⁱ(t) (1)

式中，λ_i为第i组高斯白噪声的权值系数；

步骤1.2、对于步骤1.1产生的含噪信号，利用EMD算法分别进行分解，将每组含噪信号通过EMD算法分解，得到N组一阶IMF分量，将N组一阶IMF分量求平均，即得到x(t)的一阶IMF分量，记为IMF₁(t)，计算如下所示：

式中IMF₁ ⁱ(t)表示第i组一阶IMF分量；

步骤1.3、根据步骤1.2得到的x(t)的一阶IMF分量，计算一阶IMF分量的残余项r₁(t)，如下所示：

r₁(t)＝x(t)-IMF₁(t) (3)

步骤1.4、根据步骤1.3中的残余项r₁(t)，给残余项r₁(t)中继续添加白噪声得到r₁(t)+λ₁E₁(wⁱ(t))，获得N组信号，通过EMD算法分解求平均，同步骤1.1～1.2，获得x(t)的二阶分量IMF₂(t)，如下所示：

式中E₁(·)表示经EMD分解后的第1个IMF分量，λ₁E₁(wⁱ(t))是添加白噪声；

步骤1.5、当k＝2,3,…,K，计算第k-1个残余分量r_k-1(t)，如式所示：

r_k-1(t)＝r_k-2(t)-IMF_k-1(t) (5)

式中r_k-2(t)是第k-2个残余分量；IMF_k-1(t)是第k-1阶IMF分量；

步骤1.6、计算r_k-1(t)+λ_k-1E_k-1(wⁱ(t))中的IMF分量并平均，获得IMF_k(t)，如下所示：

E_k-1(·)表示经EMD分解后的第k-1个IMF分量；r_k-1(t)是第k-1个残余分量；λ_k-1E_k-1(wⁱ(t))是添加白噪声；

步骤1.7、循环执行步骤1.5～步骤1.6，直到残差最多有一个极值无法分解，最终得到第K阶分量IMF_K(t)与残余项r_K，如下所示：

则原始信号x(t)表示为：

。

步骤1中完全自适应噪声集合经验模态分解CEEMDAN算法中，每个分解阶段均通过噪声系数λ选择信噪比，噪声系数λ为0.2。

步骤2、针对步骤1得到K组IMF分量，将IMF₁(t)视为主要包含噪声的分量，因此计算剩余IMF分量与IMF₁(t)的灰色关联度表征IMF分量的噪声水平，并进行降序排列；

步骤2具体如下：

步骤2.1、根据步骤1中分解的各个IMF分量，将其分别表征为序列：

IMF_k(t)＝[IMF_k(1),IMF_k(2),...,IMF_k(j),...,IMF_k(T)] (9)

其中:IMF_k是第k个IMF分量的原始序列；IMF_k(j)为IMF_k的j个元素，j＝1,2,...,T；T是序列的长度；

步骤2.2、根据步骤2.1得到的各个分量的原始序列，对其进行极性统一和平均，获得极性统一序列IMF_k ⁽⁰⁾和平均序列IMF_k ⁽¹⁾：

IMF_k ⁽⁰⁾(j)＝IMF_k(j)+|min(IMF_k(j))| (10)

式中，IMF_k ⁽⁰⁾(j)和IMF_k ⁽¹⁾(j)分别代表IMF_k ⁽⁰⁾和IMF_k ⁽¹⁾序列中的第j个元素；

步骤2.3、计算IMF₁与其余IMF分量之间的灰色关联系数序列，公式如下所示：

式中，γ_k为第k个IMF分量与IMF₁之间的灰色关联系数序列；γ_k(j)表示γ_k的第j个元素；ρ为分辨系数；

步骤2.3中分辨系数ρ取值为0.5；所述步骤2.4中权重系数w取直为1。

步骤2.4、计算IMF₁与其余IMF分量之间的灰色关联度ξ_k：

式中，w为权重系数；ξ_k为第k个IMF分量与IMF₁之间的灰色关联度，取值范围为ξ_k∈[0,1]，灰色关联度值越大，代表IMF分量噪声水平越高；

步骤2.5、根据步骤2.4得到的各IMF分量灰色关联度，对各IMF分量灰色关联度进行降序排列。

步骤3、搭建灰色模型，根据步骤2中按噪声水平降序排序的IMF分量，自适应筛选出噪声信号占主导的IMF分量和有效信号占主导的IMF分量；

步骤3具体如下：

步骤3.1、通过去趋势波动分析DFA算法得到原始车体偏移量x(t)的参数α(t)；

步骤3.2、建立大量的样本数据，利用原始车体偏移量序列的α(t)和η搭建灰色模型；

步骤3.2中搭建灰色模型的具体步骤如下：

步骤3.2.1、分别获取原始车体偏移量序列的去趋势波动分析DFA参数α(t)的原始序列和η(t)的原始序列：

α(t)＝(α(1),α(2),...α(t),...,α(T)) (14)

η(t)＝(η(1),η(2),...η(t),...,η(T)) (15)

步骤3.2.2、对步骤3.2.1得到的α(t)和η(t)原始序列进行归一化，得到归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)：

α(1)是灰色模型输入原始序列的第1个元素；η(1)是灰色模型输出原始序列的第1个元素。

步骤3.2.3、对步骤3.2.2得到的归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)进行累计生成操作，获得生成序列α⁽¹⁾(t)和η⁽¹⁾(t)：

α⁽⁰⁾(q)表示归一序列α⁽⁰⁾的第q个元素；η⁽⁰⁾(q)表示归一序列η⁽⁰⁾的第q个元素；

步骤3.2.4、计算通过生成序列η⁽¹⁾(t)计算相邻平均序列z(t)，公式如下：

步骤3.2.5、构建数据矩阵B和数据向量y_T：

式中z(2)表示相邻平均序列z(t)的第2个元素；z(3)表示相邻平均序列z(t)的第3个元素；z(T)表示相邻平均序列z(t)的第T个元素；α⁽¹⁾(2)表示生成序列α⁽¹⁾(t)的第2个元素；α⁽¹⁾(3)表示生成序列α⁽¹⁾(t)的第3个元素；α⁽¹⁾(T)表示生成序列α⁽¹⁾(t)的第T个元素；η⁽⁰⁾(2)表示原始序列η⁽⁰⁾(t)的第2个元素；η⁽⁰⁾(3)表示原始序列η⁽⁰⁾(t)的第3个元素；η⁽⁰⁾(T)表示原始序列η⁽⁰⁾(t)的第T个元素；

步骤3.2.6、根据步骤3.2.5构建数据矩阵B和数据向量y_T，计算参数向量P_T，P_T表示如下：

式中a，b为灰色模型参数；B^T表示矩阵B的转置矩阵。

步骤3.2.7、建立归一化序列模型：

式中为归一化序列的预测值；/>

步骤3.2.8、最终样本数据得到原始序列的模型为：

式中α为输入信号的DFA计算结果；为输入信号中以噪声为主的IMF分量的比例；η(t)为样本数据，其中t由α决定，t＝argmin|α-α(t)|-1；

步骤3.2.9、利用步骤3.2.8中计算出的噪声为主的IMF分量比例计算出噪声占主导的IMF分量的个数：

M₀表示噪声占主导的IMF分量的个数，K为CEEMDAN分解得到IMF分量的总数。

从而得到CEEMDAN分解的IMF分量中，噪声占主导的IMF分量。

步骤3.3、根据建立的灰色模型及x(t)的α(t)获取到噪声信号占主导分量的比例。

步骤4、根据步骤3选择的噪声信号占主导的IMF分量，进行软阈值降噪处理，然后与有效信号占主导的IMF分量进行重构；

步骤4具体如下：

步骤4.1、使用软阈值的方法对步骤3中最后得到的噪声占主导的IMF分量处理，具体如下：

式中表示降噪后的IMF分量；sign(·)返回值为整型变量的函数，若返回值大于0，则sign(·)返回为1，若返回值等于0，则sign()返回为0，若返回值小于0，则sign(·)返回为-1，T₀为根据信号特征获得的阈值，如下式(26)所示：

步骤4.2、最后，得到车体横向偏移量的重构信号如下式(27)所示：

式中表示去噪后的噪声主导信号分量；IMF_s1(t)表示有效信号主导的信号分量；r(t)表示残余分量；

最终得到经过完全自适应噪声集合经验模态分解-灰色模型CEEMDAN-Grey降噪后的车体横向偏移量

步骤5、根据步骤4重构后的信号，基于奇异谱分析SSA(Singular SpectrumAnalysis)进行二次降噪，得到采集的低信噪比列车偏移量中由于轨道水平不平顺导致的车体横移量。

步骤5具体如下：

步骤5.1、相空间重构：对步骤4中最后得到的CEEMDAN-Grey降噪后的车体横向偏移量通过选择合适的滑窗长度L，以滑动窗口的方式将序列截取为以L₁为大小的序列片段L个，其中L₁＝T-L+1，将序列片段组合，转换为二维轨迹矩阵Y，如下式(28)所示：

y₁表示序列的第一个元素；y₂表示/>序列的第二个元素；y₃表示/>序列的第三个元素；y_L1表示/>序列的第L₁个元素；/>表示/>序列的第L₁+1个元素；y_L表示/>序列的第L个元素；y_L+1表示/>序列的第L+1个元素；y_T表示/>序列的第T个元素；

步骤5.2、奇异值分解，定义矩阵S＝YY^T，其中Y^T为二维轨迹矩阵的转置，对矩阵S进行奇异值分解，获得矩阵X的特征值λ和对应的特征向量矩阵U、V，则矩阵S表示为：

S＝UΛV^T (29)

式中U、V分别为正交矩阵，是矩阵Y的左右奇异矩阵；Λ为特征值对角矩阵。

对特征值进行降序排列，表示成向量形式如下：

λ＝[λ₁，λ₂，...，λ_L],|λ₁|≥|λ₂|≥…≥|λ_L|≥0 (30)

L表示矩阵X非零奇异点的个数，同时，矩阵S的奇异值向量表示成如下形式：

矩阵S的奇异值向量σ也是轨迹矩阵Y的奇异谱。

轨迹矩阵Y为：

式中σ_r表示第r个奇异值；U_r、V_r表示第r个奇异值对应的特征向量矩阵；

步骤5.3、奇异值分组：基于奇异值差分谱法寻找奇异值突变点，当两个奇异值之间存在较大的突变，代表信号中的主导分量发生变化，降序后的奇异值σ＝[σ₁,σ₂,...,σ_L]，令b_r＝σ_r-σ_r+1，则奇异值的差分值为b＝[b₁,b₂,...b_L-1]；当差分值序列在r处出现波峰时，代表前r个奇异值对应有效信号分量，则矩阵S表达为公式(33)

S＝UΛV^T＝U_w1Λ_w1V_w1 ^T+U_w2Λ_w2V_w2 ^T＝W₁+W₂ (33)

式中，U_w1表示矩阵Y的左奇异矩阵的有效信号矩阵；Λ_w1表示矩阵Y的特征值对角矩阵的有效信号矩阵；V_w1 ^T表示矩阵Y的右奇异矩阵的有效信号矩阵；U_w2表示矩阵Y的左奇异矩阵的噪声信号矩阵；Λ_w2表示矩阵Y的特征值对角矩阵的噪声信号矩阵；V_w2 ^T表示矩阵Y的右奇异矩阵的噪声信号矩阵；W₁表示有效信号矩阵；W₂表示噪声数据矩阵；

步骤5.4、对角平均化：步骤5.3通过奇异值差分谱将轨迹矩阵分为有效信号矩阵和噪声信号矩阵，采用对角平均化的方法使有效信号矩阵W₁重构为长度为T的一维车体动态偏移量序列数据X(t′)，其中有效信号矩阵W₁为L

行L₁列的矩阵：

式中表示矩阵W₁的第d行，第t′-d+1列的元素，L^*＝min(L,L₁)，/>T＝L+L₁-1。

最终得到经过完全自适应噪声集合经验模态分解-灰色模型-奇异谱分析CEEMDAN-Grey-SSA算法有效提取的车体横向偏移量X(t′)。

本次实验以列车横向振动数据为基础，为方便计算，选择里程1024m的列车横移振动数据进行分析，数据的空间采样间隔为0.25米采样一次，共计采样点4096个，列车正常运行时速度基本保持不变，速度保持在200km/h。

基于以上数据，采用本发明的方法与SAA和CEEMDAN-Grey方法进行单变量突发故障预警，表1表明各方法去噪效果比较。

表1单变量对比实验

由表1所示的结果可以看出，本发明所提方法相比较于CEEMDAN-Grey和SSA方法信噪比有明显提高，均方根误差最小。所以本发明所提方法对于低信噪比数据有效提取其内部隐藏的弱信号有着明显的效果。

为了更加清晰地描述实验结果，将以上2个仿真结果可视化，结果见图3。通过观察图3(a)～图3(c)可知，图3(a)是原始车体偏移量与经过SSA降噪后车体偏移量的对比图，从图中可以看出SSA去噪后的曲线比较接近真实波形，但其在个别波峰处与真实值相差较大，并且由于信噪比过低，实验过程中容易过度去噪，将部分有效信息去除；图3(b)是原始车体偏移量与经过CEEMDAN-Grey自适应筛选噪声主导的IMF分量并对其进行去噪后车体偏移量的对比图，从图中可以看出在信噪比较低的情况下，去除噪声后的车体偏移量在低频区段依旧残留部分噪声信息，去噪不彻底；图3(c)是原始车体偏移量与本文所提方法CEEMDAN-Grey-SSA去噪后的车体偏移量的对比图，可以看出本文所提方法结合前两者的优点，不仅去除了低信噪比下的列车车体横移振动数据的低频噪声，同时避免了有效信息的缺失，去噪后的曲线更加趋近于理想的车体横移数据。

综上所述，所提CEEMDAN-Grey-SSA算法的去噪效果比SSA与CEEMDAN-Grey的去噪效果更好。通过以上仿真实验的结果对比，验证了使用所提方法CEEMDAN-Grey-SSA算法来提取弱信号的有效性和可行性。

Claims (9)

1.CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，具体按照以下步骤实施：

步骤1、对获取得到的车体偏移数据X(t)进行完全自适应噪声集合经验模态分解CEEMDAN分解，t＝1，2，...，TT表示采样点的总数，得到K个IMF分量信号，依次命名为IMF₁(t),IMF₂(t),…,IMF_k(t),…,IMF_K(t)，其中k＝1,2,…,K，K为CEEMDAN分解得到IMF分量的总数；

步骤2、针对步骤1得到K组IMF分量，将IMF₁(t)视为主要包含噪声的分量，因此计算剩余IMF分量与IMF₁(t)的灰色关联度表征IMF分量的噪声水平，并进行降序排列；

步骤3、搭建灰色模型，根据步骤2中按噪声水平降序排序的IMF分量，自适应筛选出噪声信号占主导的IMF分量和有效信号占主导的IMF分量；

步骤4、根据步骤3选择的噪声信号占主导的IMF分量，进行软阈值降噪处理，然后与有效信号占主导的IMF分量进行重构；

步骤5、根据步骤4重构后的信号，基于奇异谱分析SSA进行二次降噪，得到采集的低信噪比列车偏移量中由于轨道水平不平顺导致的车体横移量。

2.根据权利要求1所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤1具体按照以下步骤实施：

步骤1.1、随机生成N组白噪声序列wⁱ(t)，i＝1,2,…,N，并将白噪声序列加入原始的车体偏移量x(t)中，获得N组含噪信号，则第i组含噪信号xⁱ(t)表示如下：

xⁱ(t)＝x(t)+λ_iwⁱ(t) (1)

式中，λ_i为第i组高斯白噪声的权值系数；

步骤1.2、对于步骤1.1产生的含噪信号，利用EMD算法分别进行分解，将每组含噪信号通过EMD算法分解，得到N组一阶IMF分量，将N组一阶IMF分量求平均，即得到x(t)的一阶IMF分量，记为IMF₁(t)，计算如下所示：

式中表示第i组一阶IMF分量；

步骤1.3、根据步骤1.2得到的x(t)的一阶IMF分量，计算一阶IMF分量的残余项r₁(t)，如下所示：

r₁(t)＝x(t)-IMF₁(t) (3)

步骤1.4、根据步骤1.3中的残余项r₁(t)，给残余项r₁(t)中继续添加白噪声得到r₁(t)+λ₁E₁(wⁱ(t))，获得N组信号，通过EMD算法分解求平均，同步骤1.1～1.2，获得x(t)的二阶分量IMF₂(t)，如下所示：

式中E₁(·)表示经EMD分解后的第1个IMF分量，λ₁E₁(wⁱ(t))是添加白噪声；

步骤1.5、当k＝2,3,…,K，计算第k-1个残余分量r_k-1(t)，如式所示：

r_k-1(t)＝r_k-2(t)-IMF_k-1(t) (5)

式中r_k-2(t)是第k-2个残余分量；IMF_k-1(t)是第k-1阶IMF分量；

步骤1.6、计算r_k-1(t)+λ_k-1E_k-1(wⁱ(t))中的IMF分量并平均，获得IMF_k(t)，如下所示：

E_k-1(·)表示经EMD分解后的第k-1个IMF分量；r_k-1(t)是第k-1个残余分量；λ_k-1E_k-1(wⁱ(t))是添加白噪声；

步骤1.7、循环执行步骤1.5～步骤1.6，直到残差最多有一个极值无法分解，最终得到第K阶分量IMF_K(t)与残余项r_K，如下所示：

则原始信号x(t)表示为：

。

3.根据权利要求2所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤1中完全自适应噪声集合经验模态分解CEEMDAN算法中，每个分解阶段均通过噪声系数λ选择信噪比，噪声系数λ为0.2。

4.根据权利要求2所述的基于CEEMDAN-Grey-SSA的弱信号提取方法，其特征在于，所述步骤2具体如下：

步骤2.1、根据步骤1中分解的各个IMF分量，将其分别表征为序列：

IMF_k(t)＝[IMF_k(1),IMF_k(2),...,IMF_k(j),...,IMF_k(T)] (9)

其中:IMF_k是第k个IMF分量的原始序列；IMF_k(j)为IMF_k的j个元素，j＝1,2,...,T；T是序列的长度；

步骤2.2、根据步骤2.1得到的各个分量的原始序列，对其进行极性统一和平均，获得极性统一序列IMF_k ⁽⁰⁾和平均序列IMF_k ⁽¹⁾：

IMF_k ⁽⁰⁾(j)＝IMF_k(j)+|min(IMF_k(j))| (10)

式中，IMF_k ⁽⁰⁾(j)和IMF_k ⁽¹⁾(j)分别代表IMF_k ⁽⁰⁾和IMF_k ⁽¹⁾序列中的第j个元素；

步骤2.3、计算IMF₁与其余IMF分量之间的灰色关联系数序列，公式如下所示：

式中，γ_k为第k个IMF分量与IMF₁之间的灰色关联系数序列；γ_k(j)表示γ_k的第j个元素；ρ为分辨系数；

步骤2.4、计算IMF₁与其余IMF分量之间的灰色关联度ξ_k：

式中，w为权重系数；ξ_k为第k个IMF分量与IMF₁之间的灰色关联度，取值范围为ξ_k∈[0,1]，灰色关联度值越大，代表IMF分量噪声水平越高；

步骤2.5、根据步骤2.4得到的各IMF分量灰色关联度，对各IMF分量灰色关联度进行降序排列。

5.根据权利要求4所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤2.3中分辨系数ρ取值为0.5；所述步骤2.4中权重系数w取直为1。

6.根据权利要求4所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤3具体如下：

步骤3.1、通过去趋势波动分析DFA算法得到原始车体偏移量x(t)的参数α(t)；

步骤3.2、建立大量的样本数据，利用原始车体偏移量序列的α(t)和η搭建灰色模型；

步骤3.3、根据建立的灰色模型及x(t)的α(t)获取到噪声信号占主导分量的比例。

7.根据权利要求6所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤3.2中搭建灰色模型的具体步骤如下：

步骤3.2.1、分别获取原始车体偏移量序列的去趋势波动分析DFA参数α(t)的原始序列和η(t)的原始序列：

α(t)＝(α(1),α(2),...α(t),...,α(T)) (14)

η(t)＝(η(1),η(2),...η(t),...,η(T)) (15)

步骤3.2.2、对步骤3.2.1得到的α(t)和η(t)原始序列进行归一化，得到归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)：

α(1)是灰色模型输入原始序列的第1个元素；η(1)是灰色模型输出原始序列的第1个元素；

步骤3.2.3、对步骤3.2.2得到的归一化序列α⁽⁰⁾(t)和η⁽⁰⁾(t)进行累计生成操作，获得生成序列α⁽¹⁾(t)和η⁽¹⁾(t)：

α⁽⁰⁾(q)表示归一序列α⁽⁰⁾的第q个元素；η⁽⁰⁾(q)表示归一序列η⁽⁰⁾的第q个元素；

步骤3.2.4、计算通过生成序列η⁽¹⁾(t)计算相邻平均序列z(t)，公式如下：

步骤3.2.5、构建数据矩阵B和数据向量y_T：

式中z(2)表示相邻平均序列z(t)的第2个元素；z(3)表示相邻平均序列z(t)的第3个元素；z(T)表示相邻平均序列z(t)的第T个元素；α⁽¹⁾(2)表示生成序列α⁽¹⁾(t)的第2个元素；α⁽¹⁾(3)表示生成序列α⁽¹⁾(t)的第3个元素；α⁽¹⁾(T)表示生成序列α⁽¹⁾(t)的第T个元素；η⁽⁰⁾(2)表示原始序列η⁽⁰⁾(t)的第2个元素；η⁽⁰⁾(3)表示原始序列η⁽⁰⁾(t)的第3个元素；η⁽⁰⁾(T)表示原始序列η⁽⁰⁾(t)的第T个元素；

步骤3.2.6、根据步骤3.2.5构建数据矩阵B和数据向量y_T，计算参数向量P_T，P_T表示如下：

式中a，b为灰色模型参数；B^T表示矩阵B的转置矩阵；

步骤3.2.7、建立归一化序列模型：

式中为归一化序列的预测值；/>

步骤3.2.8、最终样本数据得到原始序列的模型为：

式中α为输入信号的DFA计算结果；为输入信号中以噪声为主的IMF分量的比例；η(t)为样本数据，其中t由α决定，t＝argmin|α-α(t)|-1；

步骤3.2.9、利用步骤3.2.8中计算出的噪声为主的IMF分量比例计算出噪声占主导的IMF分量的个数：

M₀表示噪声占主导的IMF分量的个数，K为CEEMDAN分解得到IMF分量的总数；

从而得到CEEMDAN分解的IMF分量中，噪声占主导的IMF分量。

8.根据权利要求6所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤4具体如下：

步骤4.1、使用软阈值的方法对步骤3中最后得到的噪声占主导的IMF分量处理，具体如下：

式中表示降噪后的IMF分量；sign(·)返回值为整型变量的函数，若返回值大于0，则sign(·)返回为1，若返回值等于0，则sign()返回为0，若返回值小于0，则sign(·)返回为-1，T₀为根据信号特征获得的阈值，如下式(26)所示：

步骤4.2、最后，得到车体横向偏移量的重构信号如下式(27)所示：

式中表示去噪后的噪声主导信号分量；IMF_s1(t)表示有效信号主导的信号分量；r(t)表示残余分量；

最终得到经过完全自适应噪声集合经验模态分解-灰色模型CEEMDAN-Grey降噪后的车体横向偏移量

9.根据权利要求8所述的CEEMDAN-Grey-SSA弱信号提取方法，其特征在于，所述步骤5具体如下：

步骤5.1、相空间重构：对步骤4中最后得到的CEEMDAN-Grey降噪后的车体横向偏移量通过选择合适的滑窗长度L，以滑动窗口的方式将序列截取为以L₁为大小的序列片段L个，其中L₁＝T-L+1，将序列片段组合，转换为二维轨迹矩阵Y，如下式(28)所示：

y₁表示序列的第一个元素；y₂表示/>序列的第二个元素；y₃表示/>序列的第三个元素；/>表示/>序列的第L₁个元素；/>表示/>序列的第L₁+1个元素；y_L表示/>序列的第L个元素；y_L+1表示/>序列的第L+1个元素；y_T表示/>序列的第T个元素；

步骤5.2、奇异值分解，定义矩阵S＝YY^T，其中Y^T为二维轨迹矩阵的转置，对矩阵S进行奇异值分解，获得矩阵X的特征值λ和对应的特征向量矩阵U、V，则矩阵S表示为：

S＝UΛV^T (29)

式中U、V分别为正交矩阵，是矩阵Y的左右奇异矩阵；Λ为特征值对角矩阵；

对特征值进行降序排列，表示成向量形式如下：

λ＝[λ₁，λ₂，...，λ_L],|λ₁|≥|λ₂|≥…≥|λ_L|≥0 (30)

L表示矩阵X非零奇异点的个数，同时，矩阵S的奇异值向量表示成如下形式：

矩阵S的奇异值向量σ也是轨迹矩阵Y的奇异谱；

轨迹矩阵Y为：

式中σ_r表示第r个奇异值；U_r、V_r表示第r个奇异值对应的特征向量矩阵；

步骤5.3、奇异值分组：基于奇异值差分谱法寻找奇异值突变点，当两个奇异值之间存在较大的突变，代表信号中的主导分量发生变化，降序后的奇异值σ＝[σ₁,σ₂,...,σ_L]，令b_r＝σ_r-σ_r+1，则奇异值的差分值为b＝[b₁,b₂,...b_L-1]；当差分值序列在r处出现波峰时，代表前r个奇异值对应有效信号分量，则矩阵S表达为公式(33)

S＝UΛV^T＝U_w1Λ_w1V_w1 ^T+U_w2Λ_w2V_w2 ^T＝W₁+W₂ (33)

式中，U_w1表示矩阵Y的左奇异矩阵的有效信号矩阵；Λ_w1表示矩阵Y的特征值对角矩阵的有效信号矩阵；V_w1 ^T表示矩阵Y的右奇异矩阵的有效信号矩阵；U_w2表示矩阵Y的左奇异矩阵的噪声信号矩阵；Λ_w2表示矩阵Y的特征值对角矩阵的噪声信号矩阵；V_w2 ^T表示矩阵Y的右奇异矩阵的噪声信号矩阵；W₁表示有效信号矩阵；W₂表示噪声数据矩阵；

步骤5.4、对角平均化：步骤5.3通过奇异值差分谱将轨迹矩阵分为有效信号矩阵和噪声信号矩阵，采用对角平均化的方法使有效信号矩阵W₁重构为长度为T的一维车体动态偏移量序列数据X(t′)，其中有效信号矩阵W₁为L

行L₁列的矩阵：

式中表示矩阵W₁的第d行，第t′-d+1列的元素，L^*＝min(L,L₁)，/>T＝L+L₁-1；

最终得到经过完全自适应噪声集合经验模态分解-灰色模型-奇异谱分析CEEMDAN-Grey-SSA算法有效提取的车体横向偏移量X(t′)。