1. 四川压电与声光技术研究所, 重庆 400060;
2. 西南科技大学信息工程学院, 四川 绵阳 621010;
3. 中国科学院新疆天文台, 新疆 乌鲁木齐 830011
收稿日期: 2020-01-30; 修订日期: 2020-02-28
基金项目: 国家自然科学基金（61771411）资助
Stochastic Resonance in a Bistable System with Multiplicative Square-wave Signal and Multiplicative Dichotomous Noise
1. Sichuan Institute of Piezoelectric and Acoustooptic Technology, Chongqing 400060, China;
2. School of Information Engineering, Southwest University of Science and Technology, Mianyang 621010, China;
3. Xinjiang Astronomical Observatory, Chinese Academy of Sciences, Urumqi 830011, China
随机共振是乘性在随机动力系统中的非线性现象。当系统的噪声、系统输入信号和系统三者间协同作用时，噪声对系统输出性能有所影响。随机共振的概念由文[1]提出，用于对周期性出现的地球冰川时代进行合理的建模分析。此后，随机共振现象得到广泛研究[2]。随机共振也广泛存在于天文领域，比如，太阳耀斑和喷流中，与磁重联相关的随机共振加速，不但产生电子，而且导致3He/4He千倍的增强[3]。小行星Bennu(101955)在行星碰撞和随机共振的联合作用下，经过几百万至几千万年，到达目前的运行轨道[4]。电子与加速圆环面小区域中的湍动磁场相互作用的随机共振加速，可能是人马座A的毫米和短波频谱的产生机制[5]。在射电天文接收机中也满足随机共振条件，即接收机噪声、输入信号(射电天文信号、无线电干扰信号)以及接收机相关参数，如增益、稳定性、动态范围、噪声系数、矩形系数等相互作用，共同决定了整个观测效果。
方波信号可以看作在两个电平间切换的数字信号。随着集成电路的快速发展，方波信号在数字通信系统中起着越来越重要的作用。比如，可以作为数字逻辑系统的输入信号，作为时钟信号准确地触发同步电路。另一方面，双值噪声是工程应用中常见的干扰。比如，随机电报噪声是电子工程领域一种典型的双值噪声。在半导体器件中，闪存产生的随机电报噪声影响信息的正确读取。
系统的输入信号通常具有加性的形式，即以加法的形式作用于系统。然而，在某些情况下，输入信号是状态依赖性的，输入信号是与系统状态变量相乘的形式。乘性信噪作用对射电天文接收机可能存在干扰，特别是在强干扰信号作用下，均以乘性结果输出[6]。目前，射电天文台址的无线电环境测量较多关注干扰信号频点和噪底水平[7-8]，但对强干扰可能导致的接收机内部交调后的噪声变化却没有相关研究报道。同时，该研究结果可以用于基于锁相放大的微弱信号提取[9-10]，而锁相放大后引起的信噪比提升一直是值得研究的问题。
基于此，本文旨在研究乘性方波信号和乘性双值噪声作用下非线性双稳系统中的非线性共振行为。通过仿真定量分析，得到非线性共振行为的能谱分布规律，为强无线电干扰情况下的接收机设计、干扰信号处理提供依据。
1 乘性方波和乘性双值噪声作用下的双稳系统
乘性信号和乘性噪声作用下的双稳系统满足[11-12]：
|
$
\frac{{{\text{d}}x(t)}}{{{\text{d}}t}} = [a + f(t)]x(t) - b{x^3}(t) + \xi (t)x(t) + \eta (t),
$
|
(1) |
其中，a和b为系统参数(a＞0, b＞0)，
|
$
f(t) = r + \varepsilon \zeta (t) + As(t)
$
|
(2) |
为乘性信号；ζ(t)为对称双值噪声，ζ(t)={1, -1}，其相关率为λ0；s(t)为周期为T的方波信号，
|
$
s(t) = \left\{ {\begin{array}{*{20}{r}}
{1,0 \leqslant t < T/2} \\
{ - 1,T/2 \leqslant t < T}
\end{array}} \right.;
$
|
(3) |
r为直流偏置；ε和A分别为乘性双值噪声和乘性方波信号的幅度。若f(t)=0，则(1)式为一双稳系统，一个非稳态和两个稳态${x_ \pm } = \pm \sqrt {a/b} $。限制参数f(t)的取值，使其满足a+f(t)＞0。ξ(t)和η(t)为相关白噪声，其均值为0，相关函数为
|
$
{\langle \xi ({t_1})\xi ({t_2})\rangle = 2D\delta ({t_1} - {t_2}),}
$
|
(4) |
|
$
{\langle \eta ({t_1})\eta ({t_2})\rangle = 2P\delta ({t_1} - {t_2}),}
$
|
(5) |
|
$
{\langle \xi ({t_1})\eta ({t_2})\rangle = \langle \eta ({t_1})\xi ({t_2})\rangle = 2\lambda \sqrt {DP} \delta ({t_1} - {t_2}),}
$
|
(6) |
其中，D和P分别为乘性和加性噪声的强度；λ为两噪声间的耦合强度，0＜λ＜1。
由(1)、(5)和(6)式，(1)式对应的Fokker-Plank方程可以表示为
|
$
{\frac{{\partial \rho (x,t)}}{{\partial t}} = - \frac{\partial }{{\partial x}}[F(x,t)\rho (x,t)] + \frac{{{\partial ^2}}}{{\partial {x^2}}}[G(x,t)\rho (x,t)],}
$
|
(7) |
其中，
|
$
{F(x,t) = ax - b{x^3} + Dx + \lambda \sqrt {DP} + xf(t),}
$
|
(8) |
|
$
{G(x,t) = D{x^2} + 2\lambda \sqrt {DP} x + P,}
$
|
(9) |
设方波信号为一个缓变信号，(1)式在一个周期T的时间内足够达到局域平衡，即满足绝热近似条件[2]，这样可以由(7)~(9)式推导其长时间的渐近分布函数，即
|
$
{\rho _{{\text{st}}}}(x,t) = \frac{M}{{\sqrt {G(x,t)} }}{\text{exp}}[ - \Phi (x,t)],
$
|
(10) |
其中，M为归一化常数；Φ(x, t)为修正的势能函数，
|
$
\varPhi (x,t) = \int_{ - \infty }^x - {\text{d}}{x^\prime }F({x^\prime },t)/G({x^\prime },t).
$
|
(11) |
在绝热近似条件下，粒子在两稳态的转移率可以表示为
|
$
\begin{array}{*{20}{l}}
{{W_ \pm }(\mu ) = \frac{a}{{\sqrt 2 \pi }}{\text{exp}}[\varPhi ({x_0}) - \varPhi ({x_ \pm })]} \\
{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{a}{{\sqrt 2 \pi }}{\text{exp}}\left[ {\frac{a}{{2D}} \mp 2\lambda \sqrt {abDP} /{D^2}} \right]} \\
{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \times {\text{exp}}\left[ {\frac{{a + D + bP/D - 4b{\lambda ^2}P/D}}{{2D}}{\text{ln}}\left( {\frac{P}{{Da/b \pm 2\lambda \sqrt {aDP/b} + P}}} \right)} \right].} \\
{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \times {\text{exp}}\left\{ {\frac{{\mu - \lambda \sqrt {DP} - 2b\lambda P\sqrt {DP} /{D^2}}}{{\sqrt {D[(1 - {\lambda ^2})P} }}} \right.} \\
{{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \times \left. {\left[ {{\text{arctan}}\left( {\frac{{\lambda \sqrt P }}{{D\sqrt {(1 - {\lambda ^2})P} }}} \right) - {\text{arctan}}\left( {\frac{{ \pm \sqrt {a/b} + \lambda \sqrt P /\sqrt D }}{{\sqrt {D(1 - {\lambda ^2})P} }}} \right)} \right]} \right\}}
\end{array}
$
|
(12) |
系统输出的相关函数可以表示为[12]
|
$
G(t) = [{B_1}(r,\varepsilon ,A) + {B_3}(r,\varepsilon ,A)]{\text{exp}}( - {\lambda _0}|t|) + {B_2}(r,\varepsilon ,A)\varphi (t) + C(r,\varepsilon ,A)\delta (t),
$
|
(13) |
其中，
|
$
{{B_1}(r,\varepsilon ,A) = \frac{1}{{16}}{{[w(r + \varepsilon + A) - w(r - \varepsilon - A) + w(r + \varepsilon - A) - w(r - \varepsilon + A)]}^2},}
$
|
(14) |
|
$
{{B_2}(r,\varepsilon ,A) = \frac{1}{{16}}{{[w(r + \varepsilon + A) - w(r - \varepsilon - A) - w(r + \varepsilon - A) + w(r - \varepsilon + A)]}^2},}
$
|
(15) |
|
$
{{B_3}(r,\varepsilon ,A) = \frac{1}{{16}}{{[w(r + \varepsilon + A) + w(r - \varepsilon - A) - w(r + \varepsilon - A) - w(r - \varepsilon + A)]}^2},}
$
|
(16) |
|
$
{C(r,\varepsilon ,A) = \frac{1}{4}[{C_0}(r + \varepsilon + A) + {C_0}(r - \varepsilon - A) + {C_0}(r + \varepsilon - A) + {C_0}(r - \varepsilon + A),}
$
|
(17) |
|
$
{w(\mu ) = \sqrt {\frac{a}{b}} \frac{{{W_ - }(\mu ) - {W_ + }(\mu )}}{{{W_ - }(\mu ) + {W_ + }(\mu )}},}
$
|
(18) |
|
$
{{C_0}(\mu ) = \frac{{8{W_ + }(\mu ){W_ - }(\mu )}}{{{{[{W_ + }(\mu ) + {W_ - }(\mu )]}^3}}},}
$
|
(19) |
|
$
{\varphi (t) = \frac{4}{{{\pi ^2}}}\sum\limits_{j = 0}^\infty {{{(2j + 1)}^{ - 2}}} {\text{exp}}[ - i(2j + 1)\Omega t],}
$
|
(20) |
对(10)式两边进行傅里叶变换，得
|
$
{S(\omega ) = {S_1}(0) + {S_2}(\omega ),}
$
|
(21) |
其中，
|
$
{{S_1}(0) = \frac{2}{\lambda }[{B_1}(r,\varepsilon ,A) + {B_3}(r,\varepsilon ,A)] + C(r,\varepsilon ,A),}
$
|
(22) |
|
$
{{S_2}(\omega ) = {B_2}(r,\varepsilon ,A)\phi (\omega ),}
$
|
(23) |
|
$
{\phi (\omega ) = \frac{8}{\pi }\sum\limits_{i = 0}^\infty {{{(2i + 1)}^{ - 2}}} \delta [\omega - (2i + 1)\Omega ].}
$
|
(24) |
这里S1(0)为零频率处的噪声背景功率；S2(ω)为输出信号的功率谱。系统输出信噪比定义为信号的基波频率处功率与背景噪声功率之比，即
|
$
{R_{{\text{SN}}}} = \frac{8}{\pi }\frac{{{B_2}(r,\varepsilon ,A)}}{{C(r,\varepsilon ,A) + 2[{B_1}(r,\varepsilon ,A) + {B_3}(r,\varepsilon ,A)]/\lambda }}.
$
|
(25) |
2 讨论
基于绝热近似条件推导了双稳系统的渐近分布函数和稳态间转移率。由相关函数得到了输出信号的功率谱，进而求得了系统的输出信噪比。现在基于(25)式，分析信噪比的非单调行为。利用图 1和表 1分析乘性双值噪声对系统信噪比的影响。由图 1和表 1可见，在信噪比与乘性噪声幅度ε的关系曲线上出现了两个极值，一个极大值，一个极小值。随着ε的增大，信噪比先快速增大至最大值，然后缓慢地减小到最小值，最后单调增大。较大的乘性双值噪声幅度可以提升系统的输出性能，由表 1可得，对于b=3.8的情形，ε=3.5时系统输出幅度为ε=0.25时输出幅值的2.5倍；而当b=4.0时，ε=3.5时输出幅度为ε=0.25时输出幅值的2.15倍。研究发现，乘性双值噪声对信噪比的影响与文[12-14]中加性双值噪声对信噪比的影响不同。对于加性噪声情形，在信噪比与双值噪声幅度的曲线上或者仅有一个极值，即最大值[12]，或者虽然有两个极值，但先出现最小值，后出现最大值[13-14]。
表 1 系统参数b和双值噪声幅度ε取不同值时的信噪比(10-9)
Table 1 The SNR for different values of the system parameter b and amplitude ε of the multiplicative dichotomous noise
| ε | b=3.8 | b=4.0 | b=4.2 | b=4.5 |
| 0.25 | 3.41 | 4.39 | 5.48 | 7.30 |
| 0.50 | 5.22 | 6.71 | 8.35 | 11.17 |
| 1.00 | 5.56 | 7.15 | 8.93 | 11.88 |
| 1.50 | 5.37 | 6.89 | 8.59 | 11.42 |
| 2.00 | 5.19 | 6.64 | 8.26 | 10.94 |
| 2.50 | 5.17 | 6.53 | 8.06 | 10.58 |
| 3.00 | 5.72 | 6.95 | 8.32 | 10.60 |
| 3.50 | 8.52 | 9.44 | 10.44 | 12.14 |
由图 2和表 2研究信噪比与乘性方波信号幅值的依赖关系。由图 2和表 2可见，信噪比随A的增大而非单调变化。随着A的增大，信噪比先减小至一个极小值，然后增大至一个极大值，最后单调减小。可见，相对居中的幅值A(如对于λ=0.65取A≈0.65)可以提升系统的输出信号, 而较高或较低的幅值A(对于λ=0.65，取A＞0.7或A≈0.3)则会抑制系统性能。比如，由表 2可得，对于λ=0.65的情形，A=1.0时系统输出幅度比A=0.3时的输出幅值下降约29%；而当λ=0.67时，A=1.0时系统输出幅度比A=0.3时的输出幅值下降约20%。可见，乘性方波信号与加性方波信号对信噪比的影响是不同的。实际上，文[13-14]的研究结果表明，较高的加性方波信号幅值可以提升系统的信噪比，而由图 2和表 2可见，较高的乘性方波信号幅值会抑制系统的输出性能。另外，信噪比随两个白噪声间的耦合强度λ变化也非单调变化，如图 2和表 2，当A＜0.85，信噪比随着λ的增大而减小；当A＞0.9时，信噪比随着λ的增大而增大。
表 2 噪声耦合强度λ和方波信号幅度A取不同值时的信噪比(10-5)
Table 2 The SNR for different values of coupling strength λ between the noises and the amplitude A of the multiplicative square-wave signal
| A | λ=0.65 | λ=0.67 | λ=0.69 | λ=0.72 |
| 0.30 | 3.17 | 2.95 | 2.76 | 2.51 |
| 0.50 | 3.29 | 3.07 | 2.86 | 2.58 |
| 0.65 | 3.36 | 3.19 | 3.02 | 2.77 |
| 0.70 | 3.33 | 3.20 | 3.05 | 2.83 |
| 0.75 | 3.25 | 3.17 | 3.06 | 2.88 |
| 0.80 | 3.13 | 3.10 | 3.03 | 2.90 |
| 0.85 | 2.97 | 2.98 | 2.96 | 2.88 |
| 0.90 | 2.76 | 2.81 | 2.83 | 2.83 |
| 0.95 | 2.51 | 2.60 | 2.67 | 2.72 |
| 1.00 | 2.26 | 2.36 | 2.46 | 2.58 |
系统参数b取不同值时，信噪比与直流偏值r的关系曲线如图 3，由图 3可见，每条曲线上呈现了一个峰值，即广义的随机共振现象。该现象与文[15]中的结果类似。另外，容易看出，信噪比也是参数b的非单调函数。图 4和图 5清晰地表明，随着乘性噪声强度的增大，信噪比可以出现一个共振峰，即随机共振现象，该现象与乘性噪声作用下双稳系统[12-14, 16]中出现的现象相同，而且系统参数a非单调地影响信噪比，如图 4。
图 6表明，信噪比是加性白噪声强度的非单调函数。随着加性噪声强度的增大，系统输出信噪比先单调增大，直到最大值，然后随着加性噪声强度的进一步增大而单调减小。故相对于无噪声情形，适量的加性噪声强度可以提升系统的输出信噪比。
综上，随着乘性方波信号幅度的增大，信噪比可以取得一个最大值和一个最小值(对于λ=0.72情形，在A=0.8时出现最大值为2.9，A=0.35时出现最小值为2.5)，而随着其他参数的变化仅出现一个极值，即最大值。对于加性噪声情形，在信噪比与双值噪声幅度的曲线上或者仅有一个极值，即最大值，或者虽然有两个极值，但先出现最小值，后出现最大值。另外，以前的研究表明，较高的加性方波信号幅值可以提升系统的信噪比，而本文的研究结果发现，较大的乘性方波信号幅度会抑制系统的输出信号，对于λ=0.65的情形，A=1.0时系统输出幅度比A=0.3时的输出幅值下降约29%；反之，较大的乘性双值噪声幅度可以提升系统输出信噪比，对于b=3.8情形，ε=3.5时系统输出幅度为ε=0.25时输出幅值的2.5倍。
3 总结
本文研究了乘性方波信号和乘性双值噪声作用下的双稳系统，观察到了两种随机共振现象，即传统的随机共振与广义的随机共振。当信噪比随着乘性双值噪声强度、乘性白噪声强度和加性白噪声强度变化时，乘性系统出现随机共振现象。广义的随机共振出现在信噪比与系统参数、直流偏置及乘性方波信号幅度的关系曲线上。进一步将研究结果推广，在天文信号处理方面，例如在不同接收机性能、无线电干扰信号注入情况下，对接收机噪声的影响情况等。
