基于联合去噪和伪哈密顿量的Duffing振子弱信号检测方法

摘要

本发明请求保护一种利用自相关-小波阈值变换联合去噪和伪哈密顿量的变尺度Duffing振子来检测微弱周期信号的方法，该方法包括对高频工程信号进行频率/时间尺度变换，使其转换为固定角频率1rad/s的信号，从而方便了设置系统相变阈值；搭建相关与小波阈值变换的联合去噪系统，极大程度改善了信噪比，避免了待测信号初始相位和噪声对检测结果的不利影响；构造Duffing系统伪哈密顿量实时地表征系统动力学行为，解决了定量判断系统状态时计算量大，效率低的难题。理论分析及仿真结果表明，较之传统做法，检测过程更方便快速，在实际应用中具有重大意义，因此具有很好的应用前景。

说明书

技术领域

本发明涉及到信号处理领域，具体为一种利用自相关-小波阈值变换联合去 噪和伪哈密顿量的变尺度Duffing振子来检测微弱周期信号的方法。

背景技术

信号检测技术在工程领域中有极其广泛应用，如果有用信号幅度相对于噪 声很微弱，或者有用信号自身幅度非常小，那么它就很容易被噪声淹没，对其 检测会极为困难。传统微弱信号检测方法在检测信噪比极低信号时效果很差， 而Duffing振子混沌系统由于具有对初值极端敏感、对噪声具有较好免疫力等 优点，在检测微弱信号时表现出良好效果。作为一种新的微弱信号检测方法， 混沌振子方法不是消除噪声，而是从噪声背景中提取信号，针对其独特性可将 其应用到实际工程中，包括纳伏级微弱信号检测、地震信号检测、GPS信号捕获、 转子早期故障诊断、齿轮早期疲劳裂纹检测、汽车飞轮壳检测、超声检漏、电 网局部放电窄带干扰等方面。

在混沌系统的策动力为不同的数值时，系统的运动状态在固定点、同宿轨 道、分岔、混沌、临界周期、大尺度周期等各个状态之间转变。基于混沌理论 进行信号检测的原理就是构造一个混沌系统，该系统只对特定信号敏感。设置 该混沌系统参数使其处于临界运动状态，当外界出现系统所敏感的信号时，虽 然信号很微弱，也会使系统的运动状态发生质变。这样，我们就可以根据系统 状态的转变来检测出要检测的微弱信号。

Duffing振子微弱周期信号检测方法只对低率参数信号比较敏感，并且系统 临界阈值会随待测信号频率改变，利用混沌阵列方法检测不同频率信号（郭阳 明,翟正军,姜红梅,等.基于互相关和混沌振子阵列的转子早期碰擦故障检 测[J].航空动力学报,2008,23(12):2219-2223.），这种方法需要较多杜芬 振子组成的阵元，在宽频信号检测及实际工程中不易操作和实现。而且在实际 工程应用中，待测信号大多是不同的高频信号。Duffing振子检测方法相比其它 检测方法有最低检测下限，但当信噪比很低时，系统不稳定，检测率降低，如 果能改善待测信号信噪比将提高Duffing振子的检测能力。采用自相关方法去 噪（李健,周激流,孙涛,等.自相关级联混沌振子法实现随相弱正弦信号的 检测[J].四川大学学报(工程科学版),2010,1:035.），故意失掉待测信号 相位信息同时抑制了部分噪声，但检测信噪比门限并不理想。在判断有无信号 时，相图和时序图的定性分析方法（李月,杨宝俊,石要武.色噪声背景下微 弱正弦信号的混沌检测[J].物理学报,2003,52(3):526-530.）有很强主观 性，且当所能获得的数据量较小时，很容易产生误判（刘海波，吴德伟，戴传 金，毛虎.基于Duffing振子的弱正弦信号检测方法研究[J].电子学报,2013, 41(1):8-12.）。利用Duffing振子的Lyapunov特性指数（Lyapunov Character  Exponent,LCE）定量判断系统状态时计算较复杂，实时性差，不符合工程应 用要求。利用Floquet指数判别系统状态（杨红英，叶昊，王桂增，吕琛.Dufing 振子的Lyapunov指数与Floquet指数研究[J].仪器仪表学报,2008:29(5): 927-932），计算量有所减小，但其在临界状态的特征不明显且受噪声影响较大， 这极大限制了Duffing振子在工程领域的应用。因此迫切需要寻求能在低信噪 比条件下快速判别有无高频工程信号的方法。

鉴于此，本发明提出一种利用自相关-小波阈值变换联合去噪和伪哈密顿量 的变尺度Duffing振子来检测微弱周期信号的方法。

发明内容

针对以上现有技术中的不足，本发明的目的在于提供一种降低检测复杂度； 提高强噪声背景下待测信号信噪比；实时准确的定量判断检测系统输出状态， 判断有无待测信号的子弱信号检测方法，本发明的技术方案如下：一种基于联 合去噪和伪哈密顿量的Duffing振子弱信号检测方法，其包括以下步骤：

101、对待测高频工程信号经过二次采样，即将采样频率为fs的高频工程 信号信号变换为二次采样频率为fs'的低频参数信号，其中低频参数信号 x(t)=s(t)+n(t)，s(t)是周期信号，n(t)是均值为零的高斯白噪声；

102、将步骤101中得到的低频参数信号x(t)=s(t)+n(t)通过自相关器得到输 出信号的函数表达式为：$\left(\begin{array}{c}{R}_{\mathrm{xx}}\left(\tau \right)=\underset{T_0\to \infty }{\lim }\frac{1}{T_0}{\int }_{-T_0/2}^{T_0/2}x\left(t\right)x(t-\tau )\mathrm{dt}\\ =R_{\mathrm{ss}}\left(\tau \right)+R_{\mathrm{sn}}\left(\tau \right)+R_{\mathrm{ns}}\left(\tau \right)+R_{\mathrm{nn}}\left(\tau \right)\end{array}\right),$其中，T₀为待测周期信号周 期，τ为延迟时间，R_ss(τ)和R_nn(τ)分别为周期信号与白噪声的自相关函数，R_sn(τ)为 周期信号与白噪声的互相关函数，R_ns(τ)为白噪声与周期信号的互相关函数。进 行噪声抑制得到自相关信号x(t)'，所述自相关信号x(t)'中去掉了低频参数信 号的相位信息；

103、对步骤102中得到的自相关信号x(t)'进行小波阈值变换，去除自相 关信号x(t)'中剩余的噪声得到信号x(t)''；

104、将步骤103中去噪后得到的信号x(t)''输入Duffing振子检测系统， 设置Duffing振子检测系统的初始状态参数，包括Duffing振子检测系统驱动 力的临界阈值N1、Duffing振子阻尼比及系统初始状态位置；

105、计算步骤104中输入Duffing振子检测系统的平均伪哈密顿量APH， 用T表示其大小，其中，N为动力系统的时间序列 长度，i为系统的第i个状态，x_i，y_i表示Duffing振子在系统的第i个状态下 的相图位置，设定系统伪哈密顿量门限阈值为μ，当T≥μ时，则判断系统处于 大周期状态，且有信号存在；若T<μ时，则判断信号处于混沌状态，且没有信 号存在。

进一步的，步骤103中得到信号x(t)''采用信噪比改善因子SNIR衡量去噪 效果，其计算式如下：

SNIR=SNR_out-SNR_in式中：SNR_in为输入信噪比，SNR_out为输出信噪比。

进一步的，步骤104中检测系统驱动力临界阈值设为0.8257694，Duffing 振子阻尼比设为k=0.5，系统初始状态

进一步的，步骤105中伪哈密顿量门限阈值μ=0.35。

本发明的优点及有益效果如下：

本发明检测强噪声背景下高频工程微弱信号时通过频率/时间尺度变换把

高频信号转换为固定角频率1rad/s的信号，方便了设置系统相变阈值，克服了 传统方法低频参数信号的限制；搭建相关与小波阈值变换的联合去噪系统，极 大程度改善了信噪比，避免了待测信号初始相位和噪声对检测结果的不利影响； 构造Duffing系统伪哈密顿量实时地表征系统动力学行为，解决了定量判断系 统状态时计算量大，效率低的难题。因此本发明在实际应用中具有重大意义。

附图说明

图1混沌阈值与周期策动力频率关系曲线；

图2-20dB待测工程故障信号；

图3自相关去噪后波形；

图4小波阈值变换去噪后波形；

图5不同信噪比下的SNIR；

图6Duffing系统不同状态下LCE随时间演化曲线；

图7Duffing系统PH值分布；

图8策动力幅值r与Duffing系统APH值关系曲线；

图9为本发明实施例流程框图。

具体实施方式

下面结合附图给出一个非限定性的实施例对本发明作进一步的阐述。

现结合附图1-9及对实现本发明提出的基于自相关-小波阈值变换联合去噪 和伪哈密顿量的变尺度Duffing振子来检测微弱周期信号的方案进行具体描述 如下：

高频工程信号经过二次采样变换为低频参数信号，此低频参数信号通过自 相关器，抑制部分噪声，同时去掉了待测信号相位信息，对相关信号进行小波 阈值变换，去除自相关后的残留噪声。去噪后信号输入Duffing振子检测系统， 并设置系统初始状态参数。然后计算Duffing系统平均伪哈密顿量，设定检测 系统伪哈密顿量门限阈值。最后由系统输出的伪哈密顿量判断系统状态及是否 存在待测信号。

结合实例，本发明实现待测微弱信号检测的具体步骤为：

步骤1：高频工程信号进过二次采样变换为低频参数信号。

噪声背景下微弱正弦信号检测，极为关键一点是确定系统混沌阈值r_c和相 变阈值r_d。传统上用梅尔尼科夫（Melnikov）函数进行理论计算得到混沌阈值 r_c的粗略估计值。Duffing方程的Melnikov函数形式如下：

$\left(\begin{array}{c}M\left(t_0\right)={\int }_{-\infty }^{+\infty }\left[y\left(t\right)(-\mathrm{ky}\left(t\right)+r\cos w(t+t_0))\right]\mathrm{dt}\\ =-\frac{4}{3}k\pm \sqrt{2\pi }\mathrm{wr}\mathrm{sech}\left(\frac{\mathrm{\pi w}}{2}\right)\sin w{t}_0\end{array}\right)---\left(8\right)$

根据Melnikov判据，系统出现混沌运动的条件是Melnikov函数存在简单 零点，即系统过鞍点(0，0)的鞍点型不动点的稳定不变流形与不稳定不变流形 横截相交，系统出现横截同宿点，f/k的最小值称为混沌阈值，即：

$\min \left(\frac{f}{k}\right)=\frac{4\cosh \left(\frac{\mathrm{w\pi }}{2}\right)}{3\sqrt{2\pi }w}{r}_c---\left(9\right)$

式（9）说明混沌阈值与周期策动力频率有关，当系统阻尼比k=0.5时，其 关系如附图1所示。可见当系统阻尼比k固定时，在低频段只需要很小幅度的 驱动力就会使系统产生混沌，而在高频段时需要较大的驱动力。另一方面，仿 真研究表明Duffing系统只有在低频参数条件下有较好动态特性和检测效果， 且Duffing振子检测信号时，不同频率待测信号对应的相变阈值也不同，如果 每次检测过程都要搜索相变阈值，将大大增加检测复杂度。

所以，在处理实际工程信号时，本发明对待测信号进行二次采样，即引入 变尺度系数R，对待测信号进行频率/时间尺度变换。这种变换并不改变参与计 算数据的数值，只是将其在时间尺度上进行压缩或放大。若待测信号角频率为w， 其采样频率为f_s，则数值计算的步长为dt=1/f_s。对检测系统引入变尺度系数R 相当于将信号的时间间隔增大了R倍，相应的信号角频率被压缩R倍后变为w/R， 此时数值计算步长变为dt′=Rdt=R/f_s。经过二次采样后的信号在通过Duffing系 统时可被识别，对该信号进行尺度恢复即可得到原采集信号中的特征信号。

-20dB待测工程故障信号如附图2所示，其表达式为i(t)=hcos((w+△w)t+φ)+z_s， 仿真环境下z_s为高斯白噪声，待测信号幅值h=0.000356V，角频率w=200rad/s， 待测信号与系统内置策动力频率差△w=0，采样频率f_s=10000Hz。针对高频参数 待测信号，引入变尺度系数R=200，使得变换后信号角频率w′=1rad/s，则对应的 二次采样频率f_s′=f_s/R=500Hz，数值计算步长dt′=Rdt=0.02。

步骤2：低频参数信号通过自相关器，抑制部分噪声，同时去掉了待测信号 相位信息。

设已知频率待测信号为：x(t)=s(t)+n(t)，s(t)是周期信号，n(t)是均值为零的 高斯白噪声，信号自相关输出为：

$\left(\begin{array}{c}{R}_{\mathrm{xx}}\left(\tau \right)=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{-T_0/2}^{T_0/2}x\left(t\right)x(t-\tau )\mathrm{dt}\\ =R_{\mathrm{ss}}\left(\tau \right)+R_{\mathrm{sn}}\left(\tau \right)+R_{\mathrm{ns}}\left(\tau \right)+R_{\mathrm{nn}}\left(\tau \right)\end{array}\right)---\left(10\right)$

根据互相关函数性质，由于信号s(t)与噪声n(t)不相关，并且噪声平均值为0， 得到R_sn(τ)=R_ns(τ)=0，则有

R_xx(τ)=R_ss(τ)+R_nn(τ)     （11）

随着τ增大，R_nn(τ)趋近于零，则对充分大的τ，可得R_xx(τ)=R_ss(τ)。

工程应用中，当待测信号为x(t)=h·cos(wt+Φ)+n(t)时，有

$R_{\mathrm{xx}}\left(\tau \right)=\frac{a^2}{2}\cos \left(\mathrm{w\tau }\right)+n^{\prime }\left(t\right)---\left(12\right)$

式（12）中，n′(t)是相关信号中混有的噪声。

可见，实际中由于积分时间不可能无限长，噪声只能得到一定程度抑制， 另外，自相关器输出信号相位为0，由式（7）可知，由于丢失了信号相位信息， 反而方便了检测。

步骤2中-20dB待测工程故障信号经过自相关去噪后的波形如附图3所示。

步骤3：对相关信号进行小波阈值变换，去除信号自相关后残留噪声。

小波阈值消噪过程中，信号经过小波变换后，可以认为由信号产生的小波 系数包含有信号的重要信息，其幅值大，但数目较少，而噪声对应的小波系数 幅值小。因此，通过在不同尺度上选取一合适阈值，并将小于该阈值的小波系 数置零，而保留大于该阈值的小波系数，从而使信号中的噪声得到有效抑制。 最后进行逆小波变换，得到去噪后的重构信号。

假设离散含噪待测信号可由式（13）表示为

x(i)=s(i)+n(i)  (i=1…M)     （13）

式中x=(x₁,x₂,…,x_i,…,x_M)为含噪信号，s=(s₁,s₂,…,s_i,…,s_M)为不含噪的真实信号， n=(n₁,n₂,…,n_i,…,n_M)并且n为噪声，σ为噪声标准差，小波阈值降噪具体步骤如下：

(1)选择合适的小波函数对含噪观测信号x进行L尺度小波分解，以获得 相应的尺度系数${\mathrm{Ax}}^l=({\mathrm{Ax}}_1^l,{\mathrm{Ax}}_2^l,...,{\mathrm{Ax}}_i^l,...,{\mathrm{Ax}}_M^l)$与小波系数${\mathrm{Wx}}^l=({\mathrm{Wx}}_1^l,{\mathrm{Wx}}_2^l,...,{\mathrm{Wx}}_i^l,...,{\mathrm{Wx}}_M^l),$l=1,2,…,L。

(2)由信号x，所对应的各尺度上小波系数Wx^l在某些特定位置上有较大值， 这些点对应于原始信号s的奇变位置和重要信息，对于白噪声n，它所对应的各 尺度上小波系数在Wx^l每一尺度上的分布是均匀的，并随着尺度的增加系数Wx^l的 幅值有所减少，因此，消噪办法是寻找一个合适的小波系数的变化函数，把由 噪声引起的小波系数设为0，而对特定位置的较大系数Wxl保留，计算求得相应 的小波系数${\widehat{\mathrm{Wx}}}^l=({\widehat{\mathrm{Wx}}}_1^l,{\widehat{\mathrm{Wx}}}_2^l,...,{\widehat{\mathrm{Wx}}}_i^l,...,{\widehat{\mathrm{Wx}}}_M^l).$

(3)利用尺度系数与经阈值函数变化后的小波系数 重构信号，可得消噪后真实信号s的估计 $\bar{s}=(\overline{s_1},\overline{s_2},...,\overline{s_i},...,\overline{s_M}).$

对步骤3中相关信号进行小波阈值变换去噪后波形如附图4。

本发明算法采用信噪比改善因子SNIR衡量去噪效果，其计算式如下：

SNIR=SNR_out-SNR_in     （14）

式中：SNR_in为输入信噪比，SNR_out为输出信噪比。

附图5为不同输入信噪比下对应的SNIR值，由图可知，信号输入信噪比越 低，两种去噪方法的改善因子越高，通过相关运算可以抑制部分噪声，对相关 后信号进行小波阈值变换后，信噪比又有一定程度改善，且输入信号信噪比越 低，这种改善越明显，证明了本发明方法的有效性。

步骤4：去噪后信号输入Duffing振子检测系统，设置系统初始状态参数。

传统上用LCE确定系统从混沌态跃变到周期态的相变阈值r_d的基本思想是： 最大LCE（L_max）大于零，是系统处于混沌态的标志，L_max等于零时，系统处于临 界混沌状态，当系统L_max由大于零转为小于零，说明系统从混沌态跃变到了大尺 度周期态，L_max符号转变时刻所对应驱动力幅值即为系统相变阈值。如附图6为 Duffing系统处于不同状态时的LCE随时间演化曲线。由此方法得到系统的相变 阈值r_d=0.8257694。

设置检测系统驱动力为临界阈值0.8257694，Duffing振子阻尼比k=0.5， 系统初始状态由步骤1知数值计算步长为0.02，在Matlab/Simulink 仿真环境下采用四阶龙格-库塔方法进行数值求解。

步骤5：计算Duffing系统平均伪哈密顿量，设定系统伪哈密顿量门限阈值。

考虑平面微分动力系统：

$\stackrel{·}{x}=F\left(x\right)---\left(15\right)$

式中，x=(x₁,x₂)^T，F=(f₁,f₂)^T，F∈C^r(r≥1)。

定义：若存在光滑函数H(x₁,x₂)使得即（14）式满足 方程：

${\stackrel{·}{x}}_1=\frac{\partial H}{{\partial H}_2},{\stackrel{·}{x}}_2=-\frac{\partial H}{{\partial H}_1}$

则称（15）式为平面哈密顿系统，其中H(x₁,x₂)成为该系统的哈密顿量。

对于（1）式Duffing方程，不考虑阻尼项和策动力的影响，可以改写为：

$\left(\begin{array}{c}\stackrel{·}{x}=y\\ \stackrel{·}{y}=x-x^3\end{array}\right)---\left(17\right)$

若记则上式的哈密顿量为：

$H(x,y)=\frac{1}{2}{y}^2+G\left(x\right)=\frac{1}{2}{y}^2-\frac{1}{2}{x}^2+\frac{1}{4}{x}^4---\left(18\right)$

从物理上看，上式可用来描述单位质量的质点在外力作用下所做的运动， 因而哈密顿量中代表质点动能，G(x)表示外力对质点所做功的负值，即为势 能，所以H(x,y)代表了一种能量常数。

实际应用中阻尼项和策动力对于系统的哈密顿量有一定影响，但对系统能 量分布几乎没有影响，此时的哈密顿量为伪哈密顿量（Pseudo-Hamiltonian, PH）。附图7为Duffing系统PH值分布情况，灰度图中灰度越高，PH值越大， 可见PH值越靠近系统两个鞍点越低。由此可知，两个鞍点处PH值最低，系统 混沌特性越明显PH值越低，大尺度周期状态时PH值最高。

本发明用下式构造Duffing系统平均哈密顿量（Average  Pseudo-Hamiltonian,APH）。

$T=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}{\mathrm{PH}}_i=\frac{1}{N}\underset{i=1}{\overset{N}{\Sigma }}(\frac{{y_i}^2}{2}-\frac{{x_i}^2}{2}+\frac{{x_i}^2}{4})---\left(19\right)$

式（19）中，N为动力系统的时间序列长度，i为系统的第i个状态。附图 8是策动力为rcos(t)时，策动力幅值变化时Duffing系统APH值变化情况。由图 可知系统在混沌态时APH值相对较小；而大尺度周期状态时APH值明显增大。 这主要是由于PH值是一个变化速率很快的量，只要系统从混沌态转到大尺度周 期态，PH值就迅速地有明显改变。本发明中算法正是利用这种阶跃型跳变设定 阈值来判断有无信号。

数次实验验证后选APH值判决系统状态的门限值为μ=0.35，则有：

步骤6：由系统输出的伪哈密顿量判断系统状态及是否存在待测工程信号。

由附图4看出去噪后的信号虽然仍被噪声掩盖，但通过混沌振子可以检测 出周期信号的存在。为说明本发明算法检测效果，分别把不同信噪比原始信号 经尺度变换后作为三种检测系统的输入，计算系统APH值与L_max以判决系统状态， 进而判断是否有周期信号，同时验证伪哈密顿量检测方法的有效性。实验中取 t=100π时L_max值作为最终系统状态稳定的L_max值，多次实验得到表1至表3结果。

表1混沌检测系统输出

表2相关-混沌检测系统输出

表3自相关与小波变换联合-混沌检测系统输出

从表1至表3数据看出，当T≥0.35时，L_max≤0；T<0.35时，L_max>0，即基于 伪哈密顿量和Lyapunov指数的系统状态判别方法结果是一致的。另外，混沌检 测系统能够检测的信噪比门限为-10.5dB,相关-混沌检测系统能够检测的信噪 比门限为-35.5dB，相关与小波变换联合-混沌检测系统能够检测的信噪比门限 为-39dB，由此可见本发明检测算法的有效性和优越性。

实验得到，利用系统APH值判别状态的平均计算时间为0.62s，利用系统状 态稳定时L_max值判别状态的平均计算时间为6.7s。可见，APH值算法计算效率明 显高于Lyapunov特性指数算法。

以上这些实施例应理解为仅用于说明本发明而不用于限制本发明的保护 范围。在阅读了本发明的记载的内容之后，技术人员可以对本发明作各种改 动或修改，这些等效变化和修饰同样落入本发明方法权利要求所限定的范围。