CN113935246A - 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质 - Google Patents

一种信号鲁棒稀疏时频分析方法、终端设备及存储介质 Download PDF

Info

Publication number
CN113935246A
CN113935246A CN202111270269.8A CN202111270269A CN113935246A CN 113935246 A CN113935246 A CN 113935246A CN 202111270269 A CN202111270269 A CN 202111270269A CN 113935246 A CN113935246 A CN 113935246A
Authority
CN
China
Prior art keywords
signal
component
sparse
iteration
frequency
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111270269.8A
Other languages
English (en)
Inventor
陈颖频
王海光
王灵芝
喻飞
林凡
陈育群
宋建华
何丽
陈悦
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Minnan Normal University
Original Assignee
Minnan Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Minnan Normal University filed Critical Minnan Normal University
Priority to CN202111270269.8A priority Critical patent/CN113935246A/zh
Publication of CN113935246A publication Critical patent/CN113935246A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/10—Noise analysis or noise optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种信号鲁棒稀疏时频分析方法、终端设备及存储介质，该方法中包括：S1：构建基于形态学成分分析、稀疏时频分析和帧小波变换的信号冲击噪声去除模型；S2：针对信号冲击噪声去除模型利用前向分解法，计算冲击噪声去除后信号的纹理成分yT和卡通成分yC；S3：根据信号的纹理成分yT和卡通成分yC，针对信号冲击噪声去除模型利用后向分解法，计算信号冲击噪声去除后的频谱稀疏解xi。本发明引入形态学成分分析算法将待处理信号分解为低频卡通成分、高频纹理成分以及噪声成分，并在形态学成分分析的基础上进一步引入平稳帧小波正则项，充分挖掘信号有效成分与噪声成分的差异性，达到压制噪声和鲁棒时频分析的目的。

Description

一种信号鲁棒稀疏时频分析方法、终端设备及存储介质
技术领域
本发明涉及信号处理领域，尤其涉及一种信号鲁棒稀疏时频分析方法、终端设备及存储介质。
背景技术
时频分析方法能有效反映信号在局部时间的频谱分布，是信号处理领域重要的研究分支，并且广泛应用于地球物理勘探、雷达成像、机械振动信号分析、生物医学信号分析、电力信号分析等众多领域。然而，传统时频分析方法存在分辨率不足、交叉项干扰等问题。
近年来，随着稀疏表示的兴起以及最优化理论的不断成熟和完善，基于稀疏表示的稀疏时频分析方法开始出现，一定程度上解决了时频分析分辨率不足的问题。稀疏时频分析方法虽然较为准确地刻画了有效信号的频谱，但是对噪声干扰，特别是在冲击噪声干扰下，稀疏时频分析则不能较准确地反映有效信号的频谱分布。
发明内容
为了解决上述问题，本发明提出了一种信号鲁棒稀疏时频分析方法、终端设备及存储介质。
具体方案如下：
一种信号鲁棒稀疏时频分析方法，包括以下步骤：
S1：构建基于形态学成分分析、稀疏时频分析和帧小波变换的信号冲击噪声去除模型：
Figure BDA0003327805290000021
其中，MCA(.)表示形态学成分分析算法，yT表示信号yn的纹理成分，yC表示信号yn的卡通成分，yiT表示将信号yT进行分段后的第i个子信号，yiC表示将信号yC进行分段后的第i个子信号，g表示滑动窗口，Θ＝SF-1表示稀疏变换矩阵，S＝[I|O]表示截断矩阵，I表示单位矩阵，O表示所有元素为零的矩阵，F表示傅里叶变换矩阵，
Figure BDA0003327805290000022
表示L2范数，
Figure BDA0003327805290000023
表示Lp伪范数，p为LP伪范数中稀疏变量的稀疏程度的控制参数，xi表示信号频谱的稀疏解，μ表示平衡参数；
S2：针对信号冲击噪声去除模型利用前向分解法，计算冲击噪声去除后信号的纹理成分yT和卡通成分yC；
S3：根据信号的纹理成分yT和卡通成分yC，针对信号冲击噪声去除模型利用后向分解法，计算信号冲击噪声去除后的频谱稀疏解xi。
进一步的，步骤S2具体包括以下步骤：
S21：利用前向分解法得到冲击噪声去除后信号的纹理成分yT和卡通成分yC的求解模型为：
Figure BDA0003327805290000024
其中，α0表示保真项平衡参数，α1和α2表示正则项的平衡参数，
Figure BDA0003327805290000025
Figure BDA0003327805290000026
均表示Lp伪范数，p0、p1、p2均为LP伪范数中稀疏变量的稀疏程度的控制参数，D(.)表示一阶平稳帧小波正变换，m表示掩码向量；
S22：利用交替乘子迭代法，引入中间辅助变量q0、q1、q2对应的拉格朗日乘子
Figure BDA0003327805290000031
二次惩罚项
Figure BDA0003327805290000032
和二次惩罚项系数λ0、λ1、λ2后，将纹理成分yT和卡通成分yC的求解转换为对q0、q1、q2、
Figure BDA0003327805290000033
的求解，即：
Figure BDA0003327805290000034
Figure BDA0003327805290000035
S23：通过迭代训练计算纹理成分yT和卡通成分yC。
进一步的，步骤S23中通过迭代训练计算纹理成分yT和卡通成分yC的具体过程包括：初始设定yC、yT、q0、q1、q2、
Figure BDA0003327805290000036
均为0，设定迭代次数k＝0，迭代次数阈值K；在每次迭代中，通过对q0、q1、q2、
Figure BDA0003327805290000037
的更新来计算更新后的yC和yT，当迭代次数k＝K时，停止迭代，输出此时的yC和yT作为冲击噪声去除后信号后最终的纹理成分yT和卡通成分yC。
进一步的，q0、q1、q2、
Figure BDA0003327805290000038
的更新公式分别为：
Figure BDA0003327805290000039
Figure BDA00033278052900000310
Figure BDA00033278052900000311
Figure BDA00033278052900000312
Figure BDA00033278052900000313
Figure BDA00033278052900000314
其中，上标(k)和(k+1)分别表示第k次和第k+1次迭代，γ表示学习率参数。
进一步的，步骤S3具体包括以下步骤：
S31：利用后向分解法得到的信号冲击噪声去除后的频谱稀疏解xi的求解模型为：
Figure BDA0003327805290000041
S32：将
Figure BDA0003327805290000042
表示为子信号si，引入中间变量z＝xi与其对偶变量
Figure BDA0003327805290000043
则频谱稀疏解xi的求解模型的增广拉格朗日函数为：
Figure BDA0003327805290000044
其中，β表示二次惩罚项的系数；
S33：通过迭代训练对xi、z和
Figure BDA00033278052900000412
进行更新，得到最终的频谱稀疏解xi。
进一步的，步骤S33中迭代训练的具体过程为：初始设定xi、z、
Figure BDA0003327805290000045
均为0，设定迭代次数k＝0，迭代次数阈值K；在每次迭代中，根据迭代前的xi、z、
Figure BDA0003327805290000046
更新迭代后的xi、z、
Figure BDA0003327805290000047
当迭代次数k＝K时，停止迭代，输出此时的xi作为最终的频谱稀疏解xi。
进一步的，xi、z、
Figure BDA0003327805290000048
的更新公式分别为：
Figure BDA0003327805290000049
Figure BDA00033278052900000410
Figure BDA00033278052900000411
其中，shrinkp表示收缩算子，ΘH表示Θ的共轭转置矩阵。
进一步的，每次迭代时，xi按照更新公式更新后，还需要对其进行中心化处理。
一种信号鲁棒稀疏时频分析终端设备，包括处理器、存储器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本发明实施例上述的方法的步骤。
一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现本发明实施例上述的方法的步骤。
本发明采用如上技术方案，并具有有益效果：
(1)采用“先检测，后成像”的策略，利用冲击噪声的幅度特征提前定位冲击噪声的“位置信息”，再结合噪声的“稀疏先验”统计特性将噪声有效去除。
(2)在时频分析算法中引入形态学成分分析，将信号分解为低频卡通成分和高频纹理成分，利用帧波正则化技术在帧小波变换域重构信号的低频和高频成分，在此基础上做稀疏时频成像。
(3)将前向后向分解法和ADMM结合起来有效求解提出模型。
(4)为降低算法复杂度，针对模型中的扁矩阵特有的矩阵结构，提出一种优化算法，加速算法中的求逆运算。
附图说明
图1所示为本发明实施例一中滑动窗加权子信号示意图。
图2所示为该实施例中一阶平稳帧小波变换与逆变换示意图，其中，图2(a)为一阶平稳帧小波正变换示意图，图2(b)为一阶平稳帧小波逆变换示意图。
图3所示为该实施例中各类范数示意图，其中，图3(a)为L2范数等高线示意图，图3(b)为L1范数等高线示意图，图3(c)为Lp伪范数等高线示意图。
图4所示为本发明实施例一的流程图。
图5所示为该实施例中方法在地震信号频谱分析中的应用，其中，图5(a)为未受冲击信号污染的地震信号，图5(b)为受冲击信号污染的地震信号，图5(c)为本实施例方法恢复的地震信号，图5(d)为未受冲击噪声干扰地震信号的稀疏时频分布，图5(e)为受到冲击噪声干扰的地震信号时频分布，图5(f)为本实施例方法得到的鲁棒时频分布。
图6所示为该实施例中方法对单分量信号的测试结果，其中，图6(a)为原始信号，图6(b)为被冲击噪声干扰的信号，图6(c)为理想时频分布，图6(d)为直接对干扰信号做时频分析的结果，图6(e)为采用本实施例方法得到的时频分析结果，图6(f)为采用本实施例方法恢复的信号，图6(g)为采用本实施例方法恢复的低频成分，图6(h)为采用本实施例方法恢复的高频成分。
图7所示为该实施例中方法对多分量信号的测试结果，其中，图7(a)为原始信号，图7(b)为被冲击噪声干扰的信号，图7(c)为理想时频分布，图7(d)为直接对干扰信号做时频分析的结果，图7(e)为采用本实施例方法得到的时频分析结果，图7(f)为采用本实施例方法恢复的信号，图7(g)为采用本实施例方法恢复的低频成分，图7(h)为采用本实施例方法恢复的高频成分。
具体实施方式
为进一步说明各实施例，本发明提供有附图。这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理。配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点。
现结合附图和具体实施方式对本发明进一步说明。
实施例一：
首先介绍与本实施例相关的预备知识。
(1)稀疏时频分析模型
图1展示了短时傅里叶变换利用滑动窗口加权子信号的过程。图1中信号
Figure BDA0003327805290000071
子信号
Figure BDA0003327805290000072
滑动窗口
Figure BDA0003327805290000073
通过短时滑动窗口可将信号分解，从而获得信号局部时间的频谱分布。
图1中，令
Figure BDA0003327805290000074
表示加权后的子信号，稀疏时频分析的出发点是在频域上找到一个频谱的稀疏解
Figure BDA0003327805290000075
使该频谱满足下式：
Figure BDA0003327805290000076
其中，
Figure BDA0003327805290000077
为稀疏变换矩阵，S＝[I|O]为截断矩阵，其作用是获取反演信号
Figure BDA0003327805290000078
的前M个点，
Figure BDA0003327805290000079
为单位矩阵，
Figure BDA00033278052900000710
表示所有元素为零的矩阵，
Figure BDA00033278052900000711
表示傅里叶变换矩阵，
Figure BDA00033278052900000712
Figure BDA00033278052900000713
表示L2范数，
Figure BDA00033278052900000714
表示Lp伪范数，p为LP伪范数中稀疏变量的稀疏程度的控制参数，μ表示平衡参数。
从式(1)可以看出，反演得到的时域信号
Figure BDA0003327805290000081
再经过截断操作
Figure BDA0003327805290000082
表示获取反演时域信号的前M个点使得yi≈Θxi，同时利用稀疏正则项
Figure BDA0003327805290000083
刻画频谱的稀疏性。
(2)形态学成分分析(Morphological Component Analysis，MCA)
形态学成分分析的基本思想是将被处理信号分解成含低频信息的卡通成分、含高频信息的纹理成分和冲击噪声成分。
(3)一阶平稳帧小波变换
一阶平稳帧小波变换及其逆变换示意图如图2所示。图2中，
Figure BDA0003327805290000084
表示低通分析滤波器，
Figure BDA0003327805290000085
和
Figure BDA0003327805290000086
为高通分析滤波器，
Figure BDA0003327805290000087
表示低通综合滤波器，
Figure BDA0003327805290000088
和
Figure BDA0003327805290000089
表示高通综合滤波器。从帧小波变换和逆变换示意图中可以得出，帧小波变换及其逆变换算法没有Mallat算法的下采样和上采样操作，因此不需要数据补零，具有平移不变性，且相比于传统的小波变换，帧小波变换增加了一个高通滤波器，能更加有效地对信号进行分离处理。
(4)Lp伪范数
由于噪声的存在，保真项
Figure BDA00033278052900000810
会受到一定的干扰，如图3中的虚线所示。图3(a)展示了L2范数的等高线与数据保真项的交点，显然该交点并不在坐标轴上，可见L2范数的促稀疏能力不强。如图3(b)所示，基于L1范数的等高线与被扰动的数据保真项交于坐标轴附近，说明L1范数的促稀疏能力强于L2范数。但L1范数依然容易受到噪声干扰。如图3(c)所示，Lp伪范数的等高线则较好地与数据保真项交于坐标轴上的点，可见Lp伪范数继承了L1范数促稀疏优点的同时，对噪声干扰更加鲁棒。
基于上述的预备知识，本发明实施例以地震信号的稀疏时频分析为例，提出了一种信号鲁棒稀疏时频分析方法，如图4所示，所述方法包括以下步骤：
S1：构建基于形态学成分分析、稀疏时频分析和帧小波变换的信号冲击噪声去除模型：
Figure BDA0003327805290000091
其中，MCA(.)表示形态学成分分析算法，yT表示信号yn的纹理成分，yC表示信号yn的卡通成分，yiT表示将信号yT进行分段后的第i个子信号，yiC表示将信号yC进行分段后的第i个子信号。
S2：针对信号冲击噪声去除模型利用前向分解法，计算冲击噪声去除后信号的纹理成分yT和卡通成分yC。
该实施例中步骤S2具体包括以下步骤：
S21：利用前向分解法得到冲击噪声去除后信号的纹理成分yT和卡通成分yC的求解模型为：
Figure BDA0003327805290000092
其中，α0表示保真项平衡参数，α1和α2表示正则项的平衡参数，
Figure BDA0003327805290000093
Figure BDA0003327805290000094
均表示Lp伪范数，p0、p1、p2均为LP伪范数中稀疏变量的稀疏程度的控制参数，D(.)表示一阶平稳帧小波正变换，m表示掩码向量。
S22：利用交替乘子迭代法，引入中间辅助变量q0、q1、q2对应的拉格朗日乘子
Figure BDA0003327805290000095
二次惩罚项
Figure BDA0003327805290000096
和二次惩罚项系数λ0、λ1、λ2后，将纹理成分yT和卡通成分yC的求解转换为对q0、q1、q2、
Figure BDA0003327805290000101
的求解，即将式(3)转换为以下的增广拉格朗日函数：
Figure BDA0003327805290000102
为求解上述增广拉格朗日函数，需要分别求解每个变量的子问题。
(1)yC子问题求解，固定其他变量，则可得：
Figure BDA0003327805290000103
对上式进行配平方，则yC子问题转化为：
Figure BDA0003327805290000104
令
Figure BDA0003327805290000105
得：
Figure BDA0003327805290000106
整理得：
Figure BDA0003327805290000107
又因
Figure BDA0003327805290000108
则上式(8)又可整理为：
Figure BDA0003327805290000109
利用共轭梯度算法(Conjugate Gradient Method，CGM)求解
Figure BDA00033278052900001010
CGM算法如表1中的算法1所示。
表1
Figure BDA0003327805290000111
(2)yT子问题求解：
Figure BDA0003327805290000121
与yC子问题求解类似地，采取配方法可得：
Figure BDA0003327805290000122
同理，令
Figure BDA0003327805290000123
得：
Figure BDA0003327805290000124
利用共轭梯度法求解
Figure BDA0003327805290000125
(3)q0子问题求解：
Figure BDA0003327805290000126
对上式配平方得：
Figure BDA0003327805290000127
得：
Figure BDA0003327805290000128
(4)q1子问题求解：
Figure BDA0003327805290000129
经配平方后得：
Figure BDA0003327805290000131
得：
Figure BDA0003327805290000132
(5)q2子问题求解：
Figure BDA0003327805290000133
经配平方后得：
Figure BDA0003327805290000134
得：
Figure BDA0003327805290000135
(6)
Figure BDA0003327805290000136
子问题求解：
Figure BDA0003327805290000137
利用梯度上升法得：
Figure BDA0003327805290000138
(7)
Figure BDA0003327805290000139
子问题求解：
Figure BDA00033278052900001310
利用梯度上升法得：
Figure BDA00033278052900001311
(8)
Figure BDA00033278052900001312
子问题求解：
Figure BDA0003327805290000141
利用梯度上升法得：
Figure BDA0003327805290000142
S23：通过迭代训练计算纹理成分yT和卡通成分yC。
该实施例中迭代训练过程的整体流程如下：
S231：初始设定yC、yT、q0、q1、q2、
Figure BDA0003327805290000143
均为0；设定迭代次数k＝0，迭代次数阈值K；
S232：根据共轭梯度算法分别求解式(9)和式(12)，更新第k+1次迭代时信号的卡通部分
Figure BDA0003327805290000144
和纹理部分
Figure BDA0003327805290000145
S233：依次根据式(15)、式(18)、式(21)分别更新第k+1次迭代时的中间辅助变量
Figure BDA0003327805290000146
S234：依次根据式(23)、式(25)、式(27)分别更新第k+1次迭代时中间辅助变量
Figure BDA0003327805290000147
对应的拉格朗日乘子
Figure BDA0003327805290000148
S235：判断k＝K是否成立，如果是，进入S236；否则，令k自加1后，返回S232；
S236：输出第k+1次迭代时信号的卡通部分
Figure BDA0003327805290000149
和纹理部分
Figure BDA00033278052900001410
作为冲击噪声去除后信号后最终的纹理成分yT和卡通成分yC。
S3：根据信号的纹理成分yT和卡通成分yC，针对信号冲击噪声去除模型利用后向分解法，计算信号冲击噪声去除后的频谱稀疏解xi。
该实施例中步骤S3具体包括以下步骤：
S31：利用后向分解法得到的信号冲击噪声去除后的频谱稀疏解xi的求解模型为：
Figure BDA0003327805290000151
S32：将
Figure BDA0003327805290000152
表示为子信号si，引入中间变量z＝xi与其对偶变量
Figure BDA0003327805290000153
则频谱稀疏解xi的求解模型的增广拉格朗日函数为：
Figure BDA0003327805290000154
其中，β表示二次惩罚项的系数。
为求解上述增广拉格朗日函数，需要分别求解每个变量的子问题。
(1)固定其他变量，经配方得xi子问题的目标函数为：
Figure BDA0003327805290000155
解得：
Figure BDA0003327805290000156
其中，
Figure BDA0003327805290000157
均为单位矩阵，ΘH表示Θ的共轭转置矩阵。
进一步的，计算的频谱还需进行中心化处理，即：
Figure BDA0003327805290000158
其中，fftshif表示中心化算子，将向量的前半部分和后半部分交换位置。
(2)z子问题的目标函数为：
Figure BDA0003327805290000159
式(33)的闭合解为：
Figure BDA00033278052900001510
其中，
Figure BDA00033278052900001511
表示Lp收缩算子。
(3)
Figure BDA00033278052900001512
子问题的目标函数为：
Figure BDA00033278052900001513
则有：
Figure BDA0003327805290000161
S33：通过迭代训练对xi、z和
Figure BDA0003327805290000162
进行更新，得到最终的频谱稀疏解xi。
该实施例中迭代训练的过程的整体流程如下：
S331：初始设定xi、z、
Figure BDA0003327805290000163
均为0；设定迭代次数k＝0，迭代次数阈值K；
S332：根据式(32)更新第k+1次迭代时的频谱稀疏解
Figure BDA0003327805290000164
S333：根据式(34)更新第k+1次迭代时的中间变量
Figure BDA0003327805290000165
S334：根据式(36)更新第k+1次迭代时的中间变量的对偶变量
Figure BDA0003327805290000166
S335：判断k＝K是否成立，如果是，进入S336；否则，令k自加1后，返回S332；
S336：输出第k+1次迭代时的频谱稀疏解
Figure BDA0003327805290000167
作为最终的频谱稀疏解xi。
实验分析
本实施例中选取实际地震信号数据进行实验以反映提出方法的正确性，实验结果如图5所示。图5(b)受到20％的冲击噪声干扰，从图5(c)可以看到，本实施例提出有效地恢复地震信号的时域信息。观察图5(e)可知，若不对冲击噪声进行压制，得到的时频图将被严重干扰，无法有效获得地震信号的真实频谱，而本实施例提出的鲁棒时频分析方法得到的时频图与未收冲击噪声干扰地震信号的时频图基本一致，可见提出方法有效压制了冲击噪声给时频图带来的干扰频谱。
另外，本实施例还在不同理论信号上进行测试，利用一系列指标将本文提出方法与其他时频分析方法进行客观评价、对比。
(1)单分量信号时频成像
以单分量正弦信号y＝0.8sin(πt)作为测试信号，为验证提出方法的有效性，对测试信号施加了10％的冲击信号对测试信号进行干扰，再用本实施例方法对被干扰的信号进行时频成像。实验如图6所示。
从图6可以看出，当信号被冲击噪声污染时，对污染的信号直接进行时频分析会出现大量的干扰频率。本实施例提出的鲁棒时频分析算法则能够很准确地信号的局部频谱，同时提出的MCA框架能有效地去除冲击噪声的干扰。
(2)多成分信号时频成像
以多分量信号y＝0.5sin(πt)+0.25sin(4πt)作为测试信号，用10％的冲击信号对测试信号进行干扰，再用本实施例提出方法对被干扰的信号进行时频成像。实验如图7所示。
从图7可以看出，提出的MCA框架能较好地分离低频成分和高频成分，在去除冲击噪声的同时恢复时域信号，同时，提出方法能获得稳健的时频分布，与标准时频分布较为接近。
本发明实施例引入形态学成分分析算法将待处理信号分解为低频卡通成分、高频纹理成分以及噪声成分。另外考虑到冲击噪声的频谱也一定程度上存在稀疏性的特征，因此稀疏时频分析方法不能有效区分信号成分和噪声成分。噪声具有随机性，与有效信号相比则缺少良好的结构相似性，因此无法被小波、帧小波等原子组成的字典稀疏表示。考虑到噪声的这个特点，本发明实施例拟在形态学成分分析的基础上进一步引入平稳帧小波正则项，充分挖掘信号有效成分与噪声成分的差异性，达到压制噪声和凸显信号的目的。
本发明实施例采用交替乘子迭代法将提出模型分解成若干个简单的子问题加以求解。最后，通过高斯噪声背景下的信号时频成像来反映提出方法的有效性。从后续实验可以看到，这种处理方式将非常有效压制噪声频谱，并对信号的频率成分加以保护。
实施例二：
本发明还提供一种信号鲁棒稀疏时频分析终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本发明实施例一的上述方法实施例中的步骤。
进一步地，作为一个可执行方案，所述信号鲁棒稀疏时频分析终端设备可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述信号鲁棒稀疏时频分析终端设备可包括，但不仅限于，处理器、存储器。本领域技术人员可以理解，上述信号鲁棒稀疏时频分析终端设备的组成结构仅仅是信号鲁棒稀疏时频分析终端设备的示例，并不构成对信号鲁棒稀疏时频分析终端设备的限定，可以包括比上述更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述信号鲁棒稀疏时频分析终端设备还可以包括输入输出设备、网络接入设备、总线等，本发明实施例对此不做限定。
进一步地，作为一个可执行方案，所称处理器可以是中央处理单元(CentralProcessing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital SignalProcessor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等，所述处理器是所述信号鲁棒稀疏时频分析终端设备的控制中心，利用各种接口和线路连接整个信号鲁棒稀疏时频分析终端设备的各个部分。
所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述信号鲁棒稀疏时频分析终端设备的各种功能。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序；存储数据区可存储根据手机的使用所创建的数据等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。
本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现本发明实施例上述方法的步骤。
所述信号鲁棒稀疏时频分析终端设备集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)以及软件分发介质等。
尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。

Claims (10)

1.一种信号鲁棒稀疏时频分析方法，其特征在于，包括以下步骤：
S1：构建基于形态学成分分析、稀疏时频分析和帧小波变换的信号冲击噪声去除模型：
Figure FDA0003327805280000011
其中，MCA(.)表示形态学成分分析算法，yT表示信号yn的纹理成分，yC表示信号yn的卡通成分，yiT表示将信号yT进行分段后的第i个子信号，yiC表示将信号yC进行分段后的第i个子信号，g表示滑动窗口，Θ＝SF-1表示稀疏变换矩阵，S＝[I|O]表示截断矩阵，I表示单位矩阵，O表示所有元素为零的矩阵，F表示傅里叶变换矩阵，
Figure FDA0003327805280000012
表示L2范数，
Figure FDA0003327805280000013
表示Lp伪范数，p为LP伪范数中稀疏变量的稀疏程度的控制参数，xi表示信号频谱的稀疏解，μ表示平衡参数；
S2：针对信号冲击噪声去除模型利用前向分解法，计算冲击噪声去除后信号的纹理成分yT和卡通成分yC；
S3：根据信号的纹理成分yT和卡通成分yC，针对信号冲击噪声去除模型利用后向分解法，计算信号冲击噪声去除后的频谱稀疏解xi。
2.根据权利要求1所述的信号鲁棒稀疏时频分析方法，其特征在于：步骤S2具体包括以下步骤：
S21：利用前向分解法得到冲击噪声去除后信号的纹理成分yT和卡通成分yC的求解模型为：
Figure FDA0003327805280000014
其中，α0表示保真项平衡参数，α1和α2表示正则项的平衡参数，
Figure FDA0003327805280000015
Figure FDA0003327805280000016
均表示Lp伪范数，p0、p1、p2均为LP伪范数中稀疏变量的稀疏程度的控制参数，D(.)表示一阶平稳帧小波正变换，m表示掩码向量；
S22：利用交替乘子迭代法，引入中间辅助变量q0、q1、q2对应的拉格朗日乘子
Figure FDA0003327805280000021
二次惩罚项
Figure FDA0003327805280000022
和二次惩罚项系数λ0、λ1、λ2后，将纹理成分yT和卡通成分yC的求解转换为对q0、q1、q2、
Figure FDA0003327805280000023
的求解，即：
Figure FDA0003327805280000024
Figure FDA0003327805280000025
S23：通过迭代训练计算纹理成分yT和卡通成分yC。
3.根据权利要求2所述的信号鲁棒稀疏时频分析方法，其特征在于：步骤S23中通过迭代训练计算纹理成分yT和卡通成分yC的具体过程包括：初始设定yC、yT、q0、q1、q2、
Figure FDA0003327805280000026
均为0，设定迭代次数k＝0，迭代次数阈值K；在每次迭代中，通过对q0、q1、q2、
Figure FDA0003327805280000027
的更新来计算更新后的yC和yT，当迭代次数k＝K时，停止迭代，输出此时的yC和yT作为冲击噪声去除后信号后最终的纹理成分yT和卡通成分yC。
4.根据权利要求3所述的信号鲁棒稀疏时频分析方法，其特征在于：q0、q1、q2、
Figure FDA0003327805280000028
的更新公式分别为：
Figure FDA0003327805280000029
Figure FDA00033278052800000210
Figure FDA0003327805280000031
Figure FDA0003327805280000032
Figure FDA0003327805280000033
Figure FDA0003327805280000034
其中，上标(k)和(k+1)分别表示第k次和第k+1次迭代，γ表示学习率参数。
5.根据权利要求1所述的信号鲁棒稀疏时频分析方法，其特征在于：步骤S3具体包括以下步骤：
S31：利用后向分解法得到的信号冲击噪声去除后的频谱稀疏解xi的求解模型为：
Figure FDA0003327805280000035
S32：将
Figure FDA0003327805280000036
表示为子信号si，引入中间变量z＝xi与其对偶变量
Figure FDA0003327805280000037
则频谱稀疏解xi的求解模型的增广拉格朗日函数为：
Figure FDA0003327805280000038
其中，β表示二次惩罚项的系数；
S33：通过迭代训练对xi、z和
Figure FDA0003327805280000039
进行更新，得到最终的频谱稀疏解xi。
6.根据权利要求5所述的信号鲁棒稀疏时频分析方法，其特征在于：步骤S33中迭代训练的具体过程为：初始设定xi、z、
Figure FDA00033278052800000310
均为0，设定迭代次数k＝0，迭代次数阈值K；在每次迭代中，根据迭代前的xi、z、
Figure FDA00033278052800000311
更新迭代后的xi、z、
Figure FDA00033278052800000312
当迭代次数k＝K时，停止迭代，输出此时的xi作为最终的频谱稀疏解xi。
7.根据权利要求6所述的信号鲁棒稀疏时频分析方法，其特征在于：xi、z、
Figure FDA0003327805280000041
的更新公式分别为：
Figure FDA0003327805280000042
Figure FDA0003327805280000043
Figure FDA0003327805280000044
其中，shrinkp表示收缩算子，ΘH表示Θ的共轭转置矩阵。
8.根据权利要求7所述的信号鲁棒稀疏时频分析方法，其特征在于：每次迭代时，xi按照更新公式更新后，还需要对其进行中心化处理。
9.一种信号鲁棒稀疏时频分析终端设备，其特征在于：包括处理器、存储器以及存储在所述存储器中并在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如权利要求1～8中任一所述方法的步骤。
10.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现如权利要求1～8中任一所述方法的步骤。
CN202111270269.8A 2021-10-29 2021-10-29 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质 Pending CN113935246A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111270269.8A CN113935246A (zh) 2021-10-29 2021-10-29 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111270269.8A CN113935246A (zh) 2021-10-29 2021-10-29 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质

Publications (1)

Publication Number Publication Date
CN113935246A true CN113935246A (zh) 2022-01-14

Family

ID=79284970

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111270269.8A Pending CN113935246A (zh) 2021-10-29 2021-10-29 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质

Country Status (1)

Country Link
CN (1) CN113935246A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115086116A (zh) * 2022-06-13 2022-09-20 重庆邮电大学 基于dct和dwt的稀疏贝叶斯电力线信道和脉冲噪声联合估计方法
CN117172135A (zh) * 2023-11-02 2023-12-05 山东省科霖检测有限公司 一种智能噪声监测管理方法与系统

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115086116A (zh) * 2022-06-13 2022-09-20 重庆邮电大学 基于dct和dwt的稀疏贝叶斯电力线信道和脉冲噪声联合估计方法
CN115086116B (zh) * 2022-06-13 2023-05-26 重庆邮电大学 基于dct和dwt的稀疏贝叶斯电力线信道和脉冲噪声联合估计方法
CN117172135A (zh) * 2023-11-02 2023-12-05 山东省科霖检测有限公司 一种智能噪声监测管理方法与系统
CN117172135B (zh) * 2023-11-02 2024-02-06 山东省科霖检测有限公司 一种智能噪声监测管理方法与系统

Similar Documents

Publication Publication Date Title
Chen Fast dictionary learning for noise attenuation of multidimensional seismic data
Yu et al. Monte Carlo data-driven tight frame for seismic data recovery
Wang et al. Fast dictionary learning for high-dimensional seismic reconstruction
CN113935246A (zh) 一种信号鲁棒稀疏时频分析方法、终端设备及存储介质
Clausel et al. The monogenic synchrosqueezed wavelet transform: a tool for the decomposition/demodulation of AM–FM images
Jia et al. A fast rank-reduction algorithm for three-dimensional seismic data interpolation
Meng et al. Estimation of chirp signals with time-varying amplitudes
CN102945548A (zh) 一种基于方向金字塔滤波的图像处理方法及装置
Zheng et al. The surface wave suppression using the second generation curvelet transform
CN110569728A (zh) 一种基于字典训练和正交匹配追踪的核信号提取方法
Zhou et al. A hybrid method for noise suppression using variational mode decomposition and singular spectrum analysis
CN116153329A (zh) 一种基于cwt-lbp的声音信号时频纹理特征提取方法
CN113917540B (zh) 基于稀疏约束的抗假频射线束进行地震数据去噪的方法
CN116699526A (zh) 一种基于稀疏与低秩模型的车载毫米波雷达干扰抑制方法
Liu et al. A dictionary learning method with atom splitting for seismic footprint suppression
Sun et al. Denoising of desert seismic signal based on synchrosqueezing transform and Adaboost algorithm
Jianhua et al. Impulse interference processing for MT data based on a new adaptive wavelet threshold de-noising method
Sun et al. Application of adaptive iterative low-rank algorithm based on transform domain in desert seismic signal analysis
CN114066749A (zh) 相位相关抗噪位移估计方法、设备及存储介质
Li et al. Simple framework for the contrastive learning of visual representations-based data-driven tight frame for seismic denoising and interpolation
Zhou et al. Coherent Noise Attenuation by Kurtosis-Guided Adaptive Dictionary Learning Based on Variational Sparse Representation
CN110687605A (zh) 基于改进的k-svd算法在地震信号处理的算法分析应用
CN112054803A (zh) 一种基于压缩感知的通信信号的分选方法
Wang et al. A Joint Framework for Seismic Signal Denoising Using Total Generalized Variation and Shearlet Transform
Shi et al. Extraction method of weak underwater acoustic signal based on the combination of wavelet transform and empirical mode decomposition

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination