楊保國，田坦，張殿倫

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

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

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

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

1 SR理論及模型

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

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

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

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

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

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

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

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

利用式(3)，可在已知系統(tǒng)參數(shù)a、b和采樣間隔時，在一定的輸入SNR下，對SR系統(tǒng)求解數(shù)值輸出.而對于雙穩(wěn)態(tài)SR非線性系統(tǒng)來說，此采樣間隔(即系統(tǒng)步長)對系統(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.注意到對高斯白噪聲的每個采樣值又都服從原來參數(shù)——均值和方差的高斯分布.式(2)可寫為

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

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

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

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

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

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

則使用h為歸一化模型(8)的計算步長.

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

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

1.4 歸一化模型的仿真驗證

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

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

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

現(xiàn)在問題如下:若有一單頻信號，幅度約為A2=3，頻率f2=100 Hz;輸入SNR約為－10.457 6 dB，即噪聲強度D2=50.則選擇采樣頻率fs2=50 kHz，如何選擇參數(shù)使這樣的輸入能在系統(tǒng)模型(2)中達到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時，可使其達到 SR狀態(tài).f0=f2/fs2保持不變，則可反推a=a'fs2、b=a3(A/A0)2和h=Ts2，此時可使得模型(2)的時域輸出達到SR狀態(tài).其輸出時域和能量譜密度如圖4所示.此例驗證了歸一化模型(8)和SR系統(tǒng)(2)的等價性.

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

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

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

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

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

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

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

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

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

1)統(tǒng)計仿真f0=0.002條件下，模型(9)的系統(tǒng)增益－輸入SNR－步長的山脊狀圖.

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

2.2 艦船線譜檢測

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

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

表1 f0=0.002時的步長Table 1 Step when f0=0.002

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

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

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

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

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

3 結(jié)束語

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

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

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

［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］胡崗.隨機力與非線性系統(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的隨機共振檢測系統(tǒng)［J］.中國儀器儀表，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é)隨機共振在數(shù)字信號傳輸中的應(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］江毅.水聲目標檢測中隨機共振方法的應(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］張曉兵，董華玉，玄兆林.水聲時變線譜建模與雙穩(wěn)態(tài)隨機共振分析［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］楊定新，胡蔦慶.隨機共振在微弱信號檢測中的數(shù)值仿真［J］.國防科技大學(xué)學(xué)報，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.