雙穩(wěn)態(tài)隨機(jī)共振系統(tǒng)參數(shù)選擇快速算法及應(yīng)用

楊保國(guó)，田坦，張殿倫

(哈爾濱工程大學(xué) 水聲技術(shù)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001)

隨機(jī)共振(stochastic resonance，SR)理論由意大利學(xué)者Benz［1］提出，用來解釋地球遠(yuǎn)古氣象中出現(xiàn)的冰期與暖氣候期周期交替出現(xiàn)的現(xiàn)象.其理論簡(jiǎn)述［1-2］為:當(dāng)淹沒在噪聲背景下的微弱信號(hào)通過一個(gè)非線性系統(tǒng)，在非線性系統(tǒng)、信號(hào)和噪聲之間達(dá)到某種匹配關(guān)系，將會(huì)發(fā)生強(qiáng)噪聲能量向微弱信號(hào)轉(zhuǎn)移的現(xiàn)象，從而達(dá)到放大信號(hào)和抑制噪聲的目的.

由于受線性響應(yīng)理論以及絕熱近似的限制［3］，工程上SR只適合頻率較低的信號(hào).而對(duì)雙穩(wěn)態(tài)的SR系統(tǒng)進(jìn)行數(shù)值計(jì)算時(shí)大多采用四階龍格庫塔算法.對(duì)于非線性系統(tǒng)，數(shù)值解法的優(yōu)劣都直接決定了系統(tǒng)的性能.使用龍格庫塔算法進(jìn)行仿真，其步長(zhǎng)就會(huì)對(duì)輸出結(jié)果產(chǎn)生影響.從物理意義上講，它應(yīng)該等于采樣頻率的倒數(shù)［5］.于是當(dāng)采用數(shù)值方法時(shí)，影響SR系統(tǒng)輸出的因素除了系統(tǒng)參數(shù)和輸入SNR之外，系統(tǒng)采樣周期也成為一個(gè)重要的因素.隨著對(duì)雙穩(wěn)態(tài)SR系統(tǒng)研究的深入以及利用數(shù)字信號(hào)處理器對(duì)其應(yīng)用的開展，參數(shù)調(diào)節(jié)SR比噪聲調(diào)節(jié)SR變得更受學(xué)者的關(guān)注［6-7］.而絕大部分的研究成果并未注明其參數(shù)是如何選取的［8-13］，即使是SR自適應(yīng)輸出算法［14］，也是改變系統(tǒng)參數(shù)，以試算的辦法來求得輸出信噪比最大.因此研究如何較為快速選取相對(duì)優(yōu)化的參數(shù)，對(duì)實(shí)現(xiàn)雙穩(wěn)態(tài)SR系統(tǒng)的數(shù)字應(yīng)用(如使用DSP器件)具有重要意義［15-19］的.

本文針對(duì)非線性雙穩(wěn)態(tài)SR模型，根據(jù)其歸一化后的形式，提出一種參數(shù)選擇快速算法.根據(jù)仿真和數(shù)據(jù)處理，驗(yàn)證了其正確性.

1 SR理論及模型

1.1 雙穩(wěn)態(tài)SR原理

非線性雙穩(wěn)態(tài)SR系統(tǒng)的勢(shì)函數(shù)［3］為

圖1 SR系統(tǒng)的對(duì)稱雙勢(shì)阱Fig.1 A symmetric double well of SR system

當(dāng)受到周期信號(hào)和高斯噪聲調(diào)制時(shí)，用郎之萬方程表示SR系統(tǒng):

式中:ξ(t)為噪聲強(qiáng)度為D的高斯白噪聲，其自相關(guān)函數(shù)為E(ξ(t)ξ(t+τ))=2Dδ(t－τ);Acos(Ωt)為周期信號(hào)，幅度為A.

1.2 雙穩(wěn)態(tài)系統(tǒng)的數(shù)值解法

郎之萬方程的最廣泛的解法為四階龍格庫塔算法［4］.具體的解法如下:

式中:xn為輸出采樣.Pn=Acos(ωt)+Γ(t)為輸入.且令xn=0.其中h為采樣步長(zhǎng)，其取值實(shí)質(zhì)上為采樣間隔.

利用式(3)，可在已知系統(tǒng)參數(shù)a、b和采樣間隔時(shí)，在一定的輸入SNR下，對(duì)SR系統(tǒng)求解數(shù)值輸出.而對(duì)于雙穩(wěn)態(tài)SR非線性系統(tǒng)來說，此采樣間隔(即系統(tǒng)步長(zhǎng))對(duì)系統(tǒng)輸出SNR增益有較大影響.

1.3 模型的歸一化與分析

由于關(guān)心的是輸出SNR及其系統(tǒng)增益，因此，把式(3)記為

式中，SNRout為輸出信噪比.

令SNRin=10lg(A2/2D)，則式(4)可寫為

將式(2)中t寫為nTs的形式，n為自然數(shù).且Ts=1/fs.令x/Ts=x'，Tsa=a'，=b'，f/fs=f0.注意到對(duì)高斯白噪聲的每個(gè)采樣值又都服從原來參數(shù)——均值和方差的高斯分布.式(2)可寫為

將此式作為新的SR系統(tǒng)，記為

用龍格庫塔算法(3)進(jìn)行數(shù)值求解時(shí)，令h'=1，則Pin(t|t=nTs)=Pin'(n)，x'(n)=kx(t|t=nTs)，k為比例系數(shù).

可以看出，在對(duì)此SR系統(tǒng)進(jìn)行數(shù)值計(jì)算時(shí)，已經(jīng)“分離”了數(shù)值解法步長(zhǎng)和采樣周期關(guān)系.

如果式(8)在某輸入SNR下達(dá)到SR狀態(tài)，只要估算出輸入序列中諧波幅值，則可反換算出式(2)的系統(tǒng)參數(shù)a、b，該參數(shù)正好使輸入達(dá)到SR狀態(tài).同時(shí)對(duì)比式(2)、(6)和(8)，可以得到結(jié)論:當(dāng)式(2)不滿足絕熱理論的條件，信號(hào)頻率較高時(shí)，保持f/fs不變且適當(dāng)?shù)姆糯笙到y(tǒng)參數(shù)a、b即可使其達(dá)到SR狀態(tài).

對(duì)比函數(shù)(5)、(7)和(9)，如果不考慮輸出序列相互的比例系數(shù)而計(jì)算輸出SNR，那么3個(gè)函數(shù)就是等價(jià)的.但是影響輸出的因素卻減少了.針對(duì)函數(shù)(9)，注意到f0完全可以用多速率數(shù)字信號(hào)處理技術(shù)來處理，則系統(tǒng)輸出SNR就變成了關(guān)于輸入SNR和算法步長(zhǎng)h″的二維函數(shù)(當(dāng)然有隱含條件A0略小于Ac).則可根據(jù)計(jì)算得到的系統(tǒng)增益和輸入SNR、算法步長(zhǎng)h″的函數(shù)圖，使用增益最大準(zhǔn)則，來選擇式(8)的參數(shù)，進(jìn)而可計(jì)算式(2)的參數(shù)a、b.

系統(tǒng)增益最大準(zhǔn)則的數(shù)學(xué)描述為:對(duì)于?SNRin(f0)，若?h，使得

則使用h為歸一化模型(8)的計(jì)算步長(zhǎng).

圖2(a)為在f0=0.01條件下系統(tǒng)增益和輸入SNR及其計(jì)算步長(zhǎng)的關(guān)系圖，可以看出其圖為山脊?fàn)?圖2(b)為SNRin=－6.989 7 dB的剖面圖.

圖2 系統(tǒng)增益Fig.2 The gain of the system

1.4 歸一化模型的仿真驗(yàn)證

已知在模型(2)中，輸入信號(hào)周期為低頻f1= 0.01 Hz，幅度為A1=0.3的余弦信號(hào)，加入噪聲強(qiáng)度D1=0.5的噪聲，則

以fs1=5 Hz采樣，通過參數(shù)a1=1、b1=1的系統(tǒng)，可使其達(dá)到SR狀態(tài).其輸入和輸出的時(shí)域波形及能量譜密度如圖3所示.

圖3 在SR系統(tǒng)模型(2)下的輸入輸出Fig.3 The input and output of the SR system model(2)

現(xiàn)在問題如下:若有一單頻信號(hào)，幅度約為A2=3，頻率f2=100 Hz;輸入SNR約為－10.457 6 dB，即噪聲強(qiáng)度D2=50.則選擇采樣頻率fs2=50 kHz，如何選擇參數(shù)使這樣的輸入能在系統(tǒng)模型(2)中達(dá)到SR狀態(tài)是本文要研究的主要問題.

圖4 模型(2)中新參數(shù)下的SR狀態(tài)Fig.4 The SR output in model(2)using new parameters

由已知可知，在模型(8)中，a'=h″=0.2時(shí)，可使其達(dá)到 SR狀態(tài).f0=f2/fs2保持不變，則可反推a=a'fs2、b=a3(A/A0)2和h=Ts2，此時(shí)可使得模型(2)的時(shí)域輸出達(dá)到SR狀態(tài).其輸出時(shí)域和能量譜密度如圖4所示.此例驗(yàn)證了歸一化模型(8)和SR系統(tǒng)(2)的等價(jià)性.

2 參數(shù)選擇快速算法及其應(yīng)用

2.1 參數(shù)選擇快速算法

在離散信號(hào)處理中，對(duì)于線譜檢測(cè)最主要的指標(biāo)是SNR增益.因此采用系統(tǒng)增益最大的次佳準(zhǔn)則來選擇系統(tǒng)參數(shù).選用次佳準(zhǔn)則的好處:1)簡(jiǎn)化參數(shù)選擇模塊的復(fù)雜程度，并保持較優(yōu)的系統(tǒng)增益; 2)當(dāng)信噪比較低時(shí)，估計(jì)諧波幅值將有較大偏差，而此時(shí)使用次佳準(zhǔn)則可使系統(tǒng)保持一定的適應(yīng)性.

次佳準(zhǔn)則的數(shù)學(xué)描述為

于?SNRin(f0)∈(SNRin1SNRin2)，均?h使得G(h)|SNRin(f0)∈(max(G(h)|SNRin(f0))－σmax(G(h) |SNRin(f0)))，則在SNRin(f0)條件下模型(8)采用h為計(jì)算步長(zhǎng).其中σ為表示區(qū)間范圍的實(shí)數(shù).

圖5 系統(tǒng)框圖Fig.5 System diagram

圖6 參數(shù)選擇流程Fig.6 Preferences flow chart

系統(tǒng)框圖如圖5所示.其中參數(shù)選擇模塊中具體流程如圖6.首先利用多速率信號(hào)處理技術(shù)使f0為參考值，計(jì)算fs，其次計(jì)算SNR，利用諧波估計(jì)來估計(jì)周期幅值A(chǔ)，并選擇相應(yīng)的步長(zhǎng)h.選擇步長(zhǎng)時(shí)參考表1，同時(shí)按照式(12)便可求得系統(tǒng)參數(shù)，在此參數(shù)下可令步長(zhǎng)h=1/fs.利用此參數(shù)計(jì)算得到的SNR增益就應(yīng)為近似最大.

其中建立表1的具體步驟如下:

1)統(tǒng)計(jì)仿真f0=0.002條件下，模型(9)的系統(tǒng)增益－輸入SNR－步長(zhǎng)的山脊?fàn)顖D.

2)根據(jù)其剖面圖，采用次佳準(zhǔn)則，即給定輸入信噪比的范圍SNRin(f0)∈(SNRin1SNRin2)，均存在h使得其對(duì)應(yīng)的增益，等于或略小于各剖面的最大增益.則在此輸入信噪比范圍內(nèi)，采用h作為系統(tǒng)步長(zhǎng).

2.2 艦船線譜檢測(cè)

在某目標(biāo)艦船離水聽器距離較近時(shí)，其輻射噪聲的功率譜密度圖為7(a)所示，其中已經(jīng)通過一個(gè)帶通濾波器，濾掉低頻干擾和高頻率成分，并且采樣頻率 f0=10 000 Hz.目標(biāo)噪聲分別在200 Hz和310 Hz附近有很強(qiáng)的線譜.

在目標(biāo)艦船較遠(yuǎn)時(shí)，功率譜密度圖如圖7(b)所示.輻射噪聲較小.并且由于環(huán)境噪聲和自噪聲的影響，目標(biāo)線譜不明顯.

表1 f0=0.002時(shí)的步長(zhǎng)Table 1 Step when f0=0.002

圖7 艦船噪聲Fig.7 Noise of a vessel

利用參數(shù)選擇快速算法來選擇SR系統(tǒng)參數(shù)，來對(duì)其進(jìn)行SR處理，以便能夠增強(qiáng)目標(biāo)線譜的強(qiáng)度.

首先利用多速率信號(hào)處理技術(shù)，使采樣頻率提高10倍，加漢寧窗后利用諧波估計(jì)技術(shù)估計(jì)諧波幅值和頻率，得到A=0.023 7、f=195 Hz.將此頻率諧波作為信號(hào)可估計(jì)SNR約為SNR=1.013 5 dB.因此選擇h=0.2，代入式(12)后計(jì)算出系統(tǒng)參數(shù)a、b.通入SR系統(tǒng)后，經(jīng)過放大，其時(shí)域波形和頻域波形如圖8所示.

可以看出，在f=195 Hz處，其線譜已經(jīng)凸顯了出來.

圖8 SR系統(tǒng)的輸出時(shí)域波形和功率譜密度Fig.8 Time-domain waveform and its spectrum of the out-put of SR system

3 結(jié)束語

本文中將常規(guī)的對(duì)稱雙穩(wěn)態(tài)SR系統(tǒng)方程進(jìn)行歸一化，據(jù)此作為參考模型，并基于龍格庫塔數(shù)值算法，提出一種針對(duì)雙穩(wěn)態(tài)SR系統(tǒng)的參數(shù)選擇快速算法.通過對(duì)模型的仿真和目標(biāo)艦船輻射噪聲的數(shù)據(jù)處理，驗(yàn)證了此算法能夠正確和快速地選擇系統(tǒng)參數(shù).但仍有兩點(diǎn)需要說明:

1)由于SR系統(tǒng)對(duì)噪聲分布的適應(yīng)性和本算法確定步長(zhǎng)時(shí)的非連續(xù)性，使得確定的系統(tǒng)參數(shù)為非最佳的.

2)利用此算法所取得的系統(tǒng)參數(shù)恢復(fù)模擬SR系統(tǒng)時(shí)，得到數(shù)值解應(yīng)逼近真值.但是由于比例系數(shù)k存在，將需要較大的放大器才能觀測(cè)到輸出信號(hào).

［1］BENZI R，SUTERA A，VULPIANI A.The mechanism of stochastic resonance［J］.J Phys A，1981，14(11):453-457.

［2］BENZI R，PARISI G，STERA A，et al.Stochastic resonance in climatic change［J］.Tellus，1982，34(10):10-16.

［3］胡崗.隨機(jī)力與非線性系統(tǒng)［M］.上海:上?？萍冀逃霭嫔纾?994:220-239.

［4］ASDIA S，TEWFIK A H.Detection of weak signals using adaptive stochastic resonance［C］//Proceedings of the 1995 IEEE International Conference on Acoustics Speech and Signal Processing.(s.l.)，1995:1332-1335.

［5］李榮華.偏微分方程數(shù)值解法［M］.北京:高等教育出版社，2005:201-350.

［6］滕萬龍，薛根，羅子健.基于DSP的隨機(jī)共振檢測(cè)系統(tǒng)［J］.中國(guó)儀器儀表，2005(11):72-74.

TENG Wanlong，XUE Gen，LUO Zijian.A stochastic resonance detection system based on DSP［J］.China Instrumentation，2005(11):72-74.

［7］段法兵.參數(shù)調(diào)節(jié)隨機(jī)共振在數(shù)字信號(hào)傳輸中的應(yīng)用［D］.杭州:浙江大學(xué)，2002:50-75.

DUAN Fabing.The application of parameter-turning stochastic resonance in digital signal transmission［D］.Hangzhou:Zhejiang Univercity，2002:50-75.

［8］FAUVE S，HELSLOT F.Stochastic resonance in a bistable system［J］.Phys Lett，1983，16(11):234-239.

［9］XU B，DUAN F，BAO R.Stochastic resonance with tuning system parameters:the application bistable systems in signal processing［J］.Chaos Solitons and Fractals，2002，13 (5):633-644.

［10］MISONO M，KOHMOTO T.Noise-enhanced transmission of information in a bistable system［J］.Phys Rev E，1998，58(10):5602-5607.

［11］江毅.水聲目標(biāo)檢測(cè)中隨機(jī)共振方法的應(yīng)用［D］.南京:東南大學(xué)，2006:105-130.

JIANG Yi.The application of stochastic resonance in the acoustic signal detection［D］.Nanjing:Southeast University，2006:105-130.

［12］張曉兵，董華玉，玄兆林.水聲時(shí)變線譜建模與雙穩(wěn)態(tài)隨機(jī)共振分析［J］.應(yīng)用聲學(xué)，2008，27(4):283-288.

ZHANG Xiaobing，DONG Huayu，XUAN Yaolin.Bistable stochastic resonance analysis for time-varying underwater acoustic line-spectrum detection［J］.Applied Acoustics，2008，27(4):283-288.

［13］楊定新，胡蔦慶.隨機(jī)共振在微弱信號(hào)檢測(cè)中的數(shù)值仿真［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2003，6(25):91-95.

YANG Dinxin，HU Niaoqing.Numerical simulation of stochastic resonance in bistable system for detecting weak signal［J］.Journal of National University of Defense Technology，2003，6(25):91-95.

［14］MITAIM S，KO B K.Adaptive stochastic resonance［J］.Proceedings of the IEEE，1998，86(11):2152-2183.

［15］DUAN F，XU B.Parameter-induced stochastic resonance and baseband binary signals transmission over AWGN channel［J］.Int J Bifurcation and Chaos，2003，13(2): 315-319.

［16］GAMMAITONI L，HANGGI P，JUNG P，MARCHESONI F.Stochastic resonace［J］.Rev Mod Phys，1998，1(70): 223-285.

［17］ZOZOR S，AMBLARD P O.Stochastic resonance in discrete time nonlinear AR(1)models［J］.IEEE Transactions on Signal Processing，1999，49(1):105-111.

［18］LYONS A P，ABRAHAM D A.Statistical characterization of high frequency shallow-water seafloor backscatter［J］.J Acoust Soc Amer，1999，2(15):524-532.

［19］CHEN H，KAY S.Theory of the stochastic resonance effect in signal detection:part 1－fixed detectors［J］.IEEE Transactions on Signal Processing，2007，6:2145-2162.