摘要:熵是热力学中表征物质状态的参量之一, 是体系混乱程度的度量. 一个信号的熵可以用来表示信号的复杂度. Duffing混沌系统从临界混沌状态向大尺度周期状态跃变的阈值是混沌系统分析的一个重要参数, 它的求解方法是混沌理论目前亟待解决的问题之一. 然而传统的实验分析法或者定量分析法存在一定的局限性. 本文在研究中发现, 系统处于混沌态和周期态时输出的多尺度熵值存在较大差异, 且当系统进入周期态后多尺度熵值趋于平稳, 基于这一现象结合遗传算法提出了基于多尺度熵的Duffing混沌系统阈值确定方法. 利用该方法对正弦信号和方波信号的检测系统跃变阈值进行了计算, 结果表明该方法快速准确且计算简单.
关键词: 微弱信号检测/
Duffing系统/
多尺度熵/
遗传算法

English Abstract

全文HTML

--> --> -->

混沌系统的初值敏感性和噪声免疫性使它在微弱信号检测领域得到了广泛的关注, 为微弱信号的检测提供了新的思路. 美国的Birx博士在他的博士论文中最先提出将混沌理论应用于弱信号的检测. 王冠宇等^[1]在国外****研究工作的基础上, 利用Duffing振子对白噪声背景下的微弱正弦信号进行了检测研究, 实现了对频率已知的待测信号幅值的估计, 检测信噪比可以达到–26 dB, 推动了信号检测领域的发展. 李月等^[2-5]研究了强噪声下混沌振子对微弱正弦信号、方波信号的检测, 色噪声背景下混沌振子对微弱正弦信号、方波信号的检测, 指出混沌振子系统对零均值噪声具有很强的免疫力. 之后大量****对Duffing检测系统进行改进, 赖智慧等^[6]对基于Homles-Duffing振子的混沌检测进行了改进, 提出变尺度微弱信号检测方法; 丛超等^[7]提出了一种基于适应步长型间歇混沌振子的信号检测方法, 用一个振子系统可以实现对任意频率任意相位的微弱周期信号的搜索检测; 牛德智等^[8]针对同频微弱信号检测时存在的盲区, 提出了一种策动力移相法予以消除; 陈志光等^[9]利用间歇混沌现象对频率未知信号进行检测并取得较好的效果; 时培明等^[10]提出了一种基于双耦合Duffing混沌振子与变尺度相结合的微弱信号检测新方法, 将其用于检测任意多频微弱信号具有明显优势. 通过大量****的不懈努力, Duffing检测不断趋于成熟.
Duffing混沌系统提供了与传统方法不同的检测途径, 在实际应用中, 混沌系统跃变阈值的确定是利用Duffing系统准确进行微弱信号检测的关键, 但是传统的时序图方法不利于计算机的自动处理. Lyapunov指数法算法复杂, 且容易受到噪声的影响. 近年来有大量****对阈值求解方法进行了研究, 提出了相图分割法^[11]、二分法^[12]、0-1方法^[13,14]、基于庞加莱界面的定量判别方法^[15]等来计算混沌系统跃变阈值, 但仍具有一定的局限性. 本文将熵引入到阈值求解, 熵可以用来表征信号复杂度, 已有不少****利用熵对信号复杂度进行分析, 梁涤青等^[16]利用小波包能量熵来判断混沌序列的复杂度; 杨孝敬等^[17]利用模糊近似熵对磁共振信号复杂度进行分析; 陈祥龙等^[18]利用多尺度样本熵对压缩机工作时的振动信号进行了分析; 王鸿姗等^[19]利用小波包样本熵通过分析声信号的复杂度进行异常音提取. 本文发现Duffing系统在不同状态下多尺度熵的明显差别, 利用这一现象通过分析系统多尺度熵与策动力幅值的变化关系, 提出一种基于多尺度熵的Duffing混沌系统跃变阈值确定方法, 并通过仿真计算说明其可行性.

Duffing振子检测系统模型为:

$\left\{ \begin{aligned}&\dot x = w \times y,\\&\dot y = w \!\times\! [ - ky \!+\! x - {x^3} \!+\! (r \!+\! h)\cos (wt) \!+\! z],\end{aligned} \right.$

其中k是阻尼比; $ x \!-\! {x^3} $是非线性恢复力; $r\cos (wt)$是周期性策动力, r是策动力幅值, h是目标信号的幅值, w是待测信号的频率, z为白噪声.
调节策动力幅值使r从0逐渐增大, 系统将会遍历同宿轨道状态、分叉状态、混沌状态、临界混沌状态, 当r大于某一值${r_{\rm{d}}}$时, 系统由临界混沌状态进入大尺度周期状态, 相轨迹图将发生明显的变化, 如图1所示, 其中${r_{\rm{d}}}$便是本文所要求解的混沌系统跃变阈值. 使用Duffing混沌系统检测微弱信号正是基于这一跃变现象: 将系统的策动力幅值r调为${r_{\rm{d}}}$, 再将待测信号添加到系统, 观察待测信号添加前后系统的相轨迹是否发生跃变, 如果发生了跃变则说明待测信号中有目标信号的存在. 可知利用这一现象检测微弱信号时, 跃变阈值的确定非常关键, 跃变阈值选取过小, 将会严重影响到系统的检测信噪比; 跃变阈值选取过大, 会导致系统状态发生跃变, 造成对检测结果的误判.
图 1 系统状态跃变　(a) 临界混沌状态; (b) 周期状态
Figure1. System state transition: (a) System state transition; (b) periodic state.

现有计算混沌系统跃变阈值的方法存在一定的局限性, 直观的相轨迹分析方法简单直观、判别方便并且不需要复杂的计算, 但是该方法受到诸多因素的影响, 并且在判别过程中会受到实验人员主观因素的影响, 精度低误差大; Melnikov和Shilnikov方法可以计算Duffing系统的混沌范围, 但是都只能说明当策动力幅值处于这个范围时, 系统可能处于混沌状态, 而不能确定系统一定处于混沌状态; Lyapunov指数用来度量动力学性态的规则性程度, 描述了系统轨迹收敛或发散的比率, 但是精确的Lyapunov指数难以得到且容易受到噪声的影响. 为此本文基于多尺度熵提出一种简单有效的计算系统跃变阈值的方法.

熵是系统混乱或无序程度的度量, 系统混乱或者无序的程度越高, 熵值越大, 反之越小. 分析熵值随策动力幅值变化关系后发现, 系统由混沌状态进入周期状态后, 系统熵值明显变小且趋于平稳, 本文利用该现象基于多尺度熵((multi-scale sample entropy, MsEn)来求解系统跃变阈值.
2

-->

3.1.多尺度熵算法

多尺度样本熵(MsEn)用于描述在不同尺度上时间序列的无规则程度, 包含参数$m, s, \tau $, 其中m是嵌入维数, s是相似系数, $\tau $是尺度因子. 多尺度样本熵由两部分构成: 一是对输入的时间序列按照尺度$\tau $进行粗细化操作得到新的时间序列; 一是计算尺度$\tau $下新时间序列的样本熵(sample entropy, SampEn). 具体步骤为^[20]:
步骤1　对原始时间序列按照尺度$\tau $进行粗细化操作
假定有一离散时间序列${x_1}, {x_2}, {x_3}, \cdots, {x_M}$共有M个点, 在尺度为$\tau $时对时间序列进行粗-断点变换得到新的时间序列为:

$y_j^{(\tau )} = \frac{1}{\tau }\sum\limits_{i = (j - 1)\tau + 1}^{j\tau } {x(i)} ,$

经过粗-断点变换后得到的新时间序列长度为$N = M/\tau $.
步骤2　计算尺度$\tau $下新时间序列的样本熵
1)将新的时间序列从${{{Y}}^{(\tau )}}(1)$到${{{Y}}^{(\tau )}}(N - m + 1)$组成一组m维矢量, 第i个矢量${{{Y}}^\tau }(i)$为:

${{{Y}}^\tau }(i) = [{y^\tau }(i),{y^\tau }(i + 1),...,{y^\tau }(i + m - 1)],$

式中$i = 1 \sim N - m + 1$.
2)定义$d[{{{Y}}^\tau }(i), {{{Y}}^\tau }(j)]$为矢量${{{Y}}^\tau }(i)$与矢量${{{Y}}^\tau }(j)$对应元素中差值最大的一个, 即:

$\begin{split}& d[{{{Y}}^\tau }(i),{{{Y}}^\tau }(j)] \\={}& \max \left| {{y^\tau }(i + k - 1) - {y^\tau }(j + k - 1)} \right|, \end{split} $

式中$k = 0 \sim m - 1$, $i, j = 1 \sim N - m + 1$.
3)给出相似系数s, 计算$d[{{{Y}}^\tau }(i), {{{Y}}^\tau }(j)] \leqslant s$的数目${D_i}$, 及该数目与总的矢量数目$N - m$的比值, 记作$D_i^{{\rm{\tau }}, m}\left( s \right)$:

$D_i^{{\rm{\tau }},m}\left( s \right) = \frac{{{D_i}}}{{N - m}}.$

对所有$D_i^{{\rm{\tau }}, m}\left( s \right)$求平均得到${D^{\tau , m}}\left( s \right)$为

${D^{\tau ,m}}\left( s \right) = \frac{1}{{N - m + 1}}\mathop \sum \limits_{i = 1}^{N - m + 1} D_i^{{\rm{\tau }},m}\left( s \right).$

4)将嵌入维数变成$m + 1$, 重复前面的过程, 得到$D_i^{{\rm{\tau }}, m + 1}\left( s \right)$和${D^{\tau , m}}\left( s \right)$:

$D_i^{{\rm{\tau }},m + 1}\left( s \right) = \frac{{{D_i}}}{{N - m - 1}},$

${D^{\tau ,m + 1}}\left( s \right) = \frac{1}{{N - m}}\mathop \sum \limits_{i = 1}^{N - m} D_i^{{\rm{\tau }},m + 1}\left( s \right).$

5)时间序列为有限长时, 得到在尺度$\tau $下的样本熵为

${\rm{SampEn}}(m,s,\tau ) = - \ln [{D^{\tau ,m + 1}}\left( s \right)/{D^{\tau ,m}}\left( s \right)].$

对原始时间序列计算每个尺度$\tau $所对应的样本熵值, 得到多尺度样本熵的计算公式为

$\begin{split}{\rm{MsEn}} \;&= \left\{ \left. \tau \right|{\rm{SampEn}}(m,s,\tau ) \right.\\ &= \left.- \ln [{D^{\tau ,m + 1}}\left( s \right)/{D^{\tau ,m}}\left( s \right)] \right\}.\end{split}$

2 -->

3.2.阈值求解

研究中发现系统混沌态和周期态输出的时间序列多尺度熵值明显不同, 且当系统进入周期态后, 系统输出序列的多尺度熵趋于平稳, 基于此通过描述系统多尺度熵值随策动力幅值变化情况可以确定系统跃变阈值.
3 -->

3.2.1.方波信号检测系统跃变阈值求解

求解方波信号检测系统的跃变阈值, 系统所采用检测方程为:

$ \left\{ \begin{aligned} & \dot x = w \times y,\\ & \dot y = w \times [ky \!+\! x - {x^3} \!+\! (r \!+\! h)g (w \!\times\! t) \!+\! z]. \end{aligned} \right. $

取w = 1 rad/s, k = 0.5, 初值x₀ = (1, 1), 时间序列长度N = 30000, 尺度因子τ = 10, 嵌入维数m = 2, 相似容限s = 0.1 SD_x (时间序列标准差), 分析多尺度熵值随策动力幅值r变化情况如图2所示. 对图2进行分析可知, 系统策动力幅值$r \leqslant 0.605$, 多尺度熵值较大且会随着策动力幅值r波动, $r > 0.605$后多尺度熵值趋于平稳且较之前明显变小, 根据熵的含义及系统变化规律可以判定$r > 0.605$后系统处于周期状态, 所以确定该系统的跃变阈值为${r_{\rm{d}}} = 0.605$. 注意在$r = 0.603$时, 系统多尺度熵值为极小值, 如果只以熵值的大小作为评判标准, 忽略平稳性的评判, 会认为此时系统已处于周期状态, 但从图中可以看到此时系统熵值并未趋于平稳, 在之后系统的熵值仍在波动, 所以仅根据熵值大小来求解阈值并不准确. 对该系统进行仿真实验, 通过分析系统相轨迹图对比该方法计算结果的准确性, 结果如图3所示, 策动力幅值$r = 0.605$时系统处于左图所示的混沌状态, 幅值增加到0.606时系统便跃变至右图所示的周期状态, 通过仿真实验确定的系统跃变阈值为0.605, 多尺度熵方法所求结果与其一致. 可见, 多尺度熵方法可以准确确定方波信号检测系统跃变阈值.
图 2 方波信号检测系统多尺度熵变化情况
Figure2. Variation of multi-scale entropy in square wave signal detection system.

图 3 仿真实验求解方波信号检测系统阈值
Figure3. Simulation experiment to solve the threshold of square wave signal detection system.

3 -->

3.2.2.正弦信号检测系统跃变阈值求解

求解正弦信号检测系统的跃变阈值, 利用Duffing系统(1)式, 取w = 1 rad/s, $k = 0.5$, 初值x₀ = (1, 1), 对于多尺度熵算法取时间序列长度N = 30000, 尺度因子$\tau = 10$, 嵌入维数$m = 2$, 相似容限$s = 0.1 S{D_x}$, 分析该系统多尺度熵值随策动力幅值r变化情况, 结果如图4所示, 可知该系统的跃变阈值为${r_{\rm{d}}} = 0.826$. 通过图5所示的仿真实验可以确定该系统的阈值为0.826, 与多尺度熵方法所求结果一致.
图 4 正弦信号检测系统多尺度熵变化情况
Figure4. Variation of multi-scale entropy in sinusoidal signal detection system.

图 5 仿真实验求解正弦信号检测系统阈值
Figure5. Simulation experiment to solve the threshold of sinusoidal signal detection system.

3 -->

3.2.3.真实信号检测系统跃变阈值求解

为了验证多尺度熵方法对真实信号检测系统的阈值求解效果, 选取一组真实的舰船信号作为样本数据, 真实待检测信号的波形如图6所示, 频域波形如图7所示, 分析可知真实水声信号包含频率为50.27 Hz的正弦信号. 在系统中内置频率为50.27 Hz的正弦信号, 求得该系统阈值为${r_{\rm{d}}} \!=\! 0.825$, 结果如图8所示, 进行仿真实验也得到${r_{\rm{d}}} = 0.825$的阈值, 如图9所示.
图 6 真实水声信号
Figure6. Real underwater acoustic signal.

图 7 真实水声信号频谱
Figure7. Spectrum of real underwater acoustic signals.

图 8 真实信号检测系统多尺度熵值变化情况
Figure8. Changes in multi-scale entropy of real signal detection system.

图 9 仿真实验求解真实水声信号检测系统阈值
Figure9. Simulation experiment to solve the threshold of real underwater acoustic signal detection system.

将策动力幅值r调节至0.825, 也即使系统处于临界混沌状态, 将真实水声信号加入检测系统, 系统相轨迹图由图10(a)所示混沌状态跃变至图10(b)所示周期状态, 成功实现对实测水声信号中目标信号的检测. 所以多尺度熵方法可以很准确地计算系统阈值.
图 10 对真实水声信号的检测　(a) 系统未添加真实信号; (b)系统添加真实信号
Figure10. Detecting real underwater acoustic signals: (a) The system did not add a real signal; (b) system adds real signal

利用系统(1)将频率参数改为5 rad/s, 其他参数不变, 得到多尺度熵变化情况如图11所示, 分析可知多尺度熵方法求得系统阈值为${r_{\rm{d}}} = 0.825$. 进行仿真实验, 仿真结果如图12所示, 可知仿真实验求得的系统阈值为${r_{\rm{d}}} = 0.826$, 与多尺度熵方法求解结果有偏差, 说明多尺度熵方法存在一定问题.
图 11 5 rad/s正弦信号检测系统多尺度熵变化情况
Figure11. Variation of multi-scale entropy in 5 rad/s sinusoidal signal detection system.

图 12 仿真实验求解5 rad/s正弦信号检测系统跃变阈值
Figure12. Simulation experiment to solve the threshold of 5 rad/s sinusoidal signal detection system.

我们分析是时间序列段的随机选取造成了阈值计算的偏差, 第三部分中计算采用的时间序列是在整个Duffing序列中随机选取的长度为30000的子序列段, 在一般状态下随机选取的子序列段可以代表整个时间序列的复杂度, 但是当系统在接近或处于临界混沌状态时会有部分序列段已经处于周期状态, 如图13所示, 该状态下随机选取的时间序列段不能代表整个时间序列的状态.
图 13 系统临界混沌状态
Figure13. Critical chaotic state of the system.

图13表示系统正处于混沌状态, 但右图中a段序列为周期序列, 取用a序列进行多尺度熵计算会得到较小熵值, 误认为系统已经进入周期状态, 导致所求系统跃变阈值偏小, 这便是在前述计算中多尺度熵方法计算结果与仿真实验计算结果存在偏差且偏小的原因. 针对这一现象考虑在整个序列中寻找复杂度最大的子序列及其对应的多尺度熵值, 用复杂度最大的子序列代表整个时间序列, 当系统复杂度最大的子序列对应的熵值处于较小值且趋于平稳时, 则系统必定已经进入周期状态, 可更加精确计算系统阈值. 为此引入遗传算法, 利用遗传算法寻找Duffing序列最大多尺度熵.
2

-->

4.1.遗传算法

遗传算法(genetic algorithm, GA)是基于生物的自然选择和遗传机理而形成的一种全局寻优算法, 其本质是一种基于概率的随机搜索算法. 一般认为GA的计算流程为^[20]:
1) 问题解的遗传表示;
2) 产生初始染色体;
3) 设计适应度函数, 根据适应值对个体进行优劣判断;
4) 进行选择操作, 选出适应值高的染色体, 使它们成为新一代种群中的染色体;
5) 对新种群进行交叉操作, 产生新的染色体;
6) 进行变异操作, 避免算法陷入局部最优解的情况;
7) 对新的种群重复进行选择、交叉、和变异操作;
8) 经过给定次数的迭代或满足给定的条件后, 把最好的染色体作为优化问题的最优解.
2 -->

4.2.最大多尺度熵算法步骤

基于遗传算法的思想, 本文提出一种最大多尺度熵算法, 算法具体步骤为:
1) 利用Duffing系统方程迭代产生一组时间序列${{X}} = [{x_1}, {x_2}, {x_3}, \cdots, {x_m}]$, 将X作为一个种群, 时间序列的长度为m;
2) 将X分成l个染色体, 分别为${{CC}} = [{C_1}, {C_2}, {C_3}, \cdots, {C_l}]$, 每个染色体的长度为n, 染色体的数目$l = m/n$. $CC$中前i个染色体构成的种群${{C}} = [{C_1}, {C_2}, {C_3}, \cdots, {C_i}]$作为算法的初始种群, 循环次数为$N = l - i$;
3) 选取多尺度样本熵作为适应度函数, 通过计算种群中每个染色体的多尺度样本熵值来判断各个染色体的优劣, 作为选择的依据;
4) 对初始种群${{C}} = [{C_1}, {C_2}, {C_3}, \cdots, {C_i}]$进行适应度测试, 计算每个染色体的适应度值, 得到种群中每个染色体所对应的适应度值分别为${{M}} = [{M_1}, {M_2}, {M_3}, \cdots, {M_i}]$;
5) 按照得到的适应度值, 从大到小对染色体进行排序, 重新排序后得到的种群为${{{C}}'} = [C_1', C_2', C_3', \cdots, C_i']$, 对应的适应度值为${{{M}}'} = [M_1', M_2', M_3', \cdots, M_i']$. 对重新排序后的种群${{{C}}'}$按照随机数$w\left( {0 < w < n} \right)$进行交叉操作, 使用参数$g\left(0 < g < i/2 \right)$来控制交叉次数. ${{\rm{C'}}_i}$由${C_{i + 1}}$替换完成淘汰工作;
6) 得到由i个新的染色体组成的种群${{{C}}^{\prime\prime}} = [C_1^{\prime\prime}, C_2^{\prime\prime}, C_3^{\prime\prime}, \cdots, C_i^{\prime\prime}]$, 返回到步骤(4)进行循环迭代计算, 直到满足循环条件;
7) 得到Duffing时间序列中熵值最大的子序列及其对应的熵值.
3 -->

4.2.1.最大多尺度熵算法有效性

第三部分5 rad/s正弦信号检测系统中, 策动力幅值$r = 0.826$, 系统处于混沌状态, 多尺度熵方法求得MsEn = 0.0778, 利用最大多尺度熵算法计算$r = 0.826$时Duffing序列最大多尺度熵值, 取Duffing序列长度为400000, 每个子序列长度为30000, 初始染色体个数为6, 交叉点为随机数, 每次循环淘汰的染色体数目为1, 循环结束后得到6个复杂度较大的子序列段, 其对应的熵值结果如表1所列.

序列	1	2	3	4	5	6
MsEn	0.1188	0.0779	0.0768	0.0796	0.0766	0.0780

表1Duffing子序列段熵值($r = 0.826$)
Table1.Entropy value of the Duffing subsequence segment ($r = 0.826$)

分析表1可知, 程序经过遗传进化找到了4组熵值大于0.0778的混沌子序列, 最大多尺度熵值为0.1188. $r = 0.827$时, 系统处于周期状态, 多尺度熵方法求得MsEn=0.0754, 利用最大多尺度熵算法计算$r = 0.827$时Duffing子序列熵值, 结果如表2所列.

序列	1	2	3	4	5	6
MsEn	0.0820	0.0822	0.0802	0.0790	0.0777	0.0753

表2Duffing子序列段MsEn值($r = 0.827$)
Table2.Entropy value of the Duffing subsequence segment($r = 0.827$)

从表2中看到, 经过遗传进化找到了5组大于0.0754的子序列, 最大MsEn为0.0822, 但因r = 0.827时系统已经处于周期状态, 所求结果与多尺度熵方法所求结果相差不大. 结合表1和表2可以知道, 最大多尺度熵算法可以很有效地寻找到复杂度较大的子序列, 得到最大多尺度熵, 可更加精确地求解系统的跃变阈值.
针对5 rad/s正弦信号检测系统的偏差问题, 现用最大多尺度熵方法求该系统阈值, 分析偏差是否仍旧存在. 取Duffing序列长度为400000, 每个子序列长度为30000, 初始染色体数目为6, 交叉点为随机数, 结果如图14所示, 分析可以确定系统跃变阈值为0.826, 与图12仿真实验所求阈值0.826相同, 所以最大多尺度熵方法可以很好地解决序列段的选取问题, 准确求解系统跃变阈值.
图 14 5 rad/s正弦信号检测系统最大多尺度熵变化情况
Figure14. Variation of maximum multi-scale entropy of 5 rad/s sinusoidal signal detection system.

本文针对Duffing系统跃变阈值难以确定这一问题进行了研究, 首次提出利用Duffing系统周期态和混沌态多尺度熵值差异明显这一现象, 通过分析多尺度熵与策动力幅值的关系对阈值进行求解, 对于该方法存在的问题, 结合遗传算法进行了改进, 通过对正弦信号与方波信号检测系统的阈值求解, 证明该方法可以快速准确得到系统跃变阈值, 解决了阈值难以快速准确得到的问题, 为混沌振子检测低信噪比信号的实际应用奠定了很好的基础.