游波，蔡志明

(海軍工程大學(xué)水聲工程系，湖北武漢430033)

水下目標(biāo)在進(jìn)行一系列戰(zhàn)術(shù)動(dòng)作時(shí)輻射瞬態(tài)信號(hào)，例如魚雷在空投和潛射過程中產(chǎn)生的瞬態(tài)信號(hào)，潛艇在轉(zhuǎn)向、加速時(shí)產(chǎn)生的瞬態(tài)信號(hào)，以及浮標(biāo)聲吶信號(hào)和主動(dòng)聲吶偵察脈沖信號(hào)等。由于潛艇和艦船吸聲降噪技術(shù)的發(fā)展，聲吶對水下目標(biāo)穩(wěn)態(tài)輻射噪聲的檢測愈加困難，而對瞬態(tài)信號(hào)的檢測將彌補(bǔ)穩(wěn)態(tài)信號(hào)檢測的弱點(diǎn)，產(chǎn)生積極的軍事效益。

瞬態(tài)信號(hào)檢測技術(shù)基于多種原理。高階譜分析方法與Power-law檢測器的結(jié)合［1］抓住瞬態(tài)信號(hào)的非高斯特征，對信號(hào)先驗(yàn)知識(shí)要求少?；谛〔ㄗ儞Q的適應(yīng)性時(shí)頻檢測方法［2］利用非高斯小波系數(shù)和噪聲系數(shù)四階矩的差異，進(jìn)一步突出瞬態(tài)信號(hào)，但小波基函數(shù)的選擇對分析結(jié)果有較大影響。希爾伯特－黃變換［3］的特征時(shí)間尺度基于信號(hào)自身特性，經(jīng)驗(yàn)?zāi)B(tài)分解是希爾伯特－黃變換理論的核心，將信號(hào)依據(jù)不同的時(shí)間尺度分解成一系列的固有模態(tài)函數(shù)。但固有模態(tài)函數(shù)截止階數(shù)的選擇會(huì)影響模態(tài)分解結(jié)果。累積和檢驗(yàn)方法(Page Test)［4-7］是一種基于概率特征分析的檢測方法，能敏感捕捉瞬態(tài)信號(hào)出現(xiàn)引起的概率分布函數(shù)變化。研究表明［8］該方法性能優(yōu)于Power-law檢測器、短時(shí)傅里葉變換、小波變換等其他方法。

根據(jù)紐曼皮爾遜準(zhǔn)則，利用檢測性能分析方法計(jì)算門限和相關(guān)參數(shù)是實(shí)現(xiàn)累積和檢驗(yàn)的關(guān)鍵。檢測性能的數(shù)值分析方法有矩陣法和快速傅里葉變換法2種。矩陣法的量化過程由于引入了高維(約幾千階)矩陣的連乘，運(yùn)算量非常大，計(jì)算精度不高?？焖俑道锶~變換法則將2個(gè)時(shí)間序列的卷積轉(zhuǎn)化為頻譜乘積的反傅里葉變換，利用快速傅里葉變換來實(shí)現(xiàn)每一步檢驗(yàn)統(tǒng)計(jì)量概率分布的計(jì)算，在高精度要求下依然具有合理的運(yùn)算量。文獻(xiàn)［9-10］對快速傅里葉變換法的介紹僅限于簡單的原理，缺乏對更新量量化范圍的分析推導(dǎo)，也沒有提及基于馬爾可夫狀態(tài)轉(zhuǎn)移矩陣的檢測概率遞歸計(jì)算方法，而這2個(gè)問題是利用快速傅里葉變換方法進(jìn)行累積和性能分析的關(guān)鍵步驟。本文解決了上述問題，并將檢測概率理論計(jì)算結(jié)果和蒙特卡洛仿真實(shí)驗(yàn)結(jié)果進(jìn)行對比，吻合度好，水池和海試數(shù)據(jù)的檢測結(jié)果進(jìn)一步表明了累積和檢驗(yàn)方法檢測低信噪比水聲瞬態(tài)信號(hào)的有效性。

1 累積和檢驗(yàn)方法性能分析模型

1.1 高斯分布方差變化假設(shè)下累積和檢驗(yàn)統(tǒng)計(jì)量

累積和檢驗(yàn)方法檢測性能與信號(hào)出現(xiàn)前后概率分布模型或參數(shù)的變化密切相關(guān)。由于瞬態(tài)信號(hào)的突發(fā)性和不穩(wěn)定性、聲信道傳播特性未知等因素，概率密度模型復(fù)雜，從一般意義上可將其看作信號(hào)出現(xiàn)前后方差發(fā)生變化的高斯分布。

給定觀測值序列X={x1，x2，…，xN}，假設(shè)瞬態(tài)信號(hào)出現(xiàn)前后，概率分布律從f0變化至f1，即從標(biāo)準(zhǔn)正態(tài)分布變化至方差為σ2的高斯正態(tài)分布:

式中:變化發(fā)生點(diǎn)θ未知。依據(jù)最大似然比準(zhǔn)則［11］得到累積和檢驗(yàn)的對數(shù)似然比統(tǒng)計(jì)量寫成:

式中:g(xn)稱為非線性量或更新量，可寫成g(xn)=x2n－b，b=σ2·lnσ2/σ2－( 1)為偏差。g(xn)要求滿足:累積和檢驗(yàn)實(shí)際上是一種基于最大對數(shù)似然比準(zhǔn)則的序貫檢測方法。設(shè)0和h分別為上、下門限，檢測過程可描述為

1.2 模型原理和量化方法

FFT法利用馬爾可夫過程來關(guān)聯(lián)被檢測序列在不同時(shí)刻狀態(tài)之間統(tǒng)計(jì)聯(lián)系。因此在該序列以及更新量可能值域范圍進(jìn)行量化，建立馬爾可夫狀態(tài)轉(zhuǎn)移矩陣是計(jì)算檢測概率的必要前提條件。

由式(2)得到n時(shí)刻Zn的概率分布律φn(z)等于n－1時(shí)刻Zn－1概率分布律fn－1(z)與更新量g(xn)概率分布律fg(z)的卷積，由于2個(gè)時(shí)間序列的卷積等同于頻譜乘積的反傅里葉變換，上述關(guān)系可表為

將復(fù)雜的卷積運(yùn)算轉(zhuǎn)化為快速離散傅里葉變換的乘積是快速傅里葉變換方法的主要思想。若檢測過程能持續(xù)到n時(shí)刻，則對φn(z)(0＜z＜h)進(jìn)行歸一化操作后，n時(shí)刻Zn的概率分布律fn(z)可表為

式中:fn(z)、fg(z)為矩陣，分別表示檢驗(yàn)統(tǒng)計(jì)量Zn和更新量g在可能值域范圍每個(gè)量化點(diǎn)的概率。將上式(5)進(jìn)行遞歸運(yùn)算，可計(jì)算任意時(shí)刻檢驗(yàn)統(tǒng)計(jì)量概率分布律。

一般可假設(shè)檢測過程從H0假設(shè)下的零狀態(tài)開始。下面以方差發(fā)生變化的高斯分布假設(shè)為例分析更新量g(xn)的量化方法。由契比雪夫不等式［12］有為數(shù)學(xué)期望，ε為任意正數(shù)。令μ=0，ε=10σ，則X以不小于99% 的 概 率 分 布 在 10σ 的 范 圍，令更新量g(xn)=yn=－b，則在H0和H12種假設(shè)下yn分布范圍的上限分別寫為

yn的下限為yn≥－b。設(shè)yn的量化范圍為［hl，hh)，檢測門限為h，下門限可表為hl= max(－h(huán)，－b)，上門限可表為hh=max(h，ynmax)。ynmax在H0和H1假設(shè)下的取值參見式(6)。設(shè)量化階數(shù)為N，量化步長為Δ=(hh－h(huán)l)/N，每一個(gè)量化點(diǎn)可表示為yi=g(xi)=hl+i·Δ，i=1，2，…，N。由雅可比行列式，經(jīng)過計(jì)算得更新量yn=g(xn)服從分布:

將式(7)、(8)代入yn=g(xn)在每一個(gè)量化點(diǎn)的值計(jì)算可得n時(shí)刻更新量等于各個(gè)量化值的概率，即fg(z)=Pr{g(xi)=yi}，組成向量，代入式(4)、(5)即可進(jìn)行遞歸運(yùn)算。文獻(xiàn)［10］的量化范圍直接取為［－h(huán)，h)，這樣做會(huì)導(dǎo)致量化值域分布范圍的不準(zhǔn)確，也會(huì)直接影響更新量概率分布律fg(z)的精度，從而給后面的參數(shù)估計(jì)帶來影響。特別是當(dāng)信號(hào)點(diǎn)數(shù)很少時(shí)，門限給出的量化范圍小于更新量可能的數(shù)值分布范圍，在一定的虛警和檢測概率下無法計(jì)算出門限和方差。

1.3 利用快速傅里葉變換方法計(jì)算檢測概率

H1假設(shè)下的檢測一般考慮2種初始狀態(tài)，一是零初始狀態(tài)，即f0(z)=δ(z)，二是在H0假設(shè)下已經(jīng)進(jìn)入的穩(wěn)定狀態(tài)fss(z)。累積和方法檢測過程中的置零復(fù)位操作(如式(2))意味著下一時(shí)刻檢測初始狀態(tài)的改變，因而置零次數(shù)直接影響更新量概率分布律的計(jì)算。在計(jì)算檢測概率時(shí)，一般將H1假設(shè)下超門限檢測具體分成2種情況:1)檢驗(yàn)統(tǒng)計(jì)量未經(jīng)任何置零復(fù)位操作超門限，置零次數(shù)k為零;2)檢驗(yàn)統(tǒng)計(jì)量經(jīng)過k(k≥1)次置零超門限。

在N個(gè)樣本長度內(nèi)未置零復(fù)位操作而發(fā)生的檢測概率可寫為

另假設(shè)在被檢測的m(1≤m≤N－1)長序列中有k次置零復(fù)位操作，1≤k≤m，記p(k)0(m)為k次置零復(fù)位(檢測判決為無信號(hào)，記為0)的概率。序列中剩余(N－m)個(gè)點(diǎn)從0狀態(tài)開始，檢測結(jié)果為有信號(hào)。因而發(fā)生k≥1次置零操作而超門限檢測概率可寫為

求解k次置零的概率p(k)0(m)，需要用到狀態(tài)轉(zhuǎn)移矩陣的遞歸求解過程。這個(gè)問題在文獻(xiàn)［7-8］中并沒有提及。如圖1所示，設(shè)檢測N長序列，k次置零操作發(fā)生在m點(diǎn)范圍內(nèi)，即1≤k≤m≤N－1。必須指出的是，k次置零操作的最后一次應(yīng)發(fā)生在第m點(diǎn)，否則這k種置零組合方式會(huì)與更小長度序列的k種置零組合方式重合，出現(xiàn)重復(fù)計(jì)算的情況。其余(k－1)次置零操作可能發(fā)生在(m－1)個(gè)點(diǎn)的范圍里，發(fā)生概率p(k－1)0(m－1)是(k－1)次獨(dú)立置零概率乘積的多種組合方式之和。

圖1 檢測概率計(jì)算的遞歸示意圖Fig.1 Illustration of the recursive calculation of Pd

建立H1假設(shè)下穩(wěn)定初始狀態(tài)時(shí)的狀態(tài)轉(zhuǎn)移矩陣C1:

式中:ei表示1·(m－k+1)維，在第i點(diǎn)上為1，其余位置為0的向量。式(13)可實(shí)現(xiàn)圖1所示k次置零概率的求解過程。將式(13)代入式(10)，得到穩(wěn)定初始狀態(tài)時(shí)發(fā)生k≥1次置零操作而超門限檢測概率，再結(jié)合式(9)，最終得到在穩(wěn)定初始狀態(tài)時(shí)瞬態(tài)信號(hào)持續(xù)過程中發(fā)生的檢測概率Pd:

當(dāng)初始狀態(tài)為零時(shí)，瞬態(tài)信號(hào)持續(xù)過程中發(fā)生的檢測概率Pd寫作:

在Pr(k≥1)的運(yùn)算過程中，第1點(diǎn)雖為零，但不計(jì)入置零操作次數(shù)，檢測m個(gè)點(diǎn)發(fā)生置零操作的次數(shù)1≤k＜m≤N－1，其余推導(dǎo)過程類似。瞬態(tài)信號(hào)持續(xù)過程中發(fā)生的檢測概率應(yīng)為式(14)和(15)之和。關(guān)于在瞬態(tài)信號(hào)結(jié)束后出現(xiàn)的延遲檢測概率計(jì)算詳見文獻(xiàn)［7］。

在n時(shí)刻結(jié)束且檢測為無信號(hào)的概率可表為

在n時(shí)刻結(jié)束且檢測為有信號(hào)的概率可表為

式中:pon的初值pon(0)設(shè)為1。式(16)～(18)體現(xiàn)了累積和算法序貫檢測的特征。其中n時(shí)刻Zn的概率分布律φn(z)利用快速傅里葉變換法由式(4)和式(5)遞歸計(jì)算得到。前面所提及H0下的穩(wěn)定狀態(tài)是指在H0假設(shè)下，當(dāng)檢測過程進(jìn)行一段時(shí)間后仍未做出判斷，即可認(rèn)為進(jìn)入了穩(wěn)定狀態(tài)［10］。穩(wěn)態(tài)時(shí)檢驗(yàn)統(tǒng)計(jì)量Zn的歸一化概率分布律可表為

2 理論計(jì)算與仿真實(shí)驗(yàn)結(jié)果比對

圖2為當(dāng)門限為40，方差為1.36時(shí)檢測概率理論計(jì)算結(jié)果和10 000次蒙特卡羅仿真結(jié)果，兩者的結(jié)果非常接近。圖3表示當(dāng)虛警間的平均間隔T－=104s，Pd=0.6時(shí)，在紐曼皮爾遜準(zhǔn)則下利用上述虛警概率和檢測概率的模型，尋優(yōu)搜索計(jì)算得到的檢測門限h和方差σ2隨瞬態(tài)信號(hào)長度的變化圖。

圖2 方差變化的高斯分布下檢測概率理論計(jì)算結(jié)果和蒙特卡羅仿真結(jié)果比較Fig.2 Comparison between the theoretical result and the Monte Carlo simulation result in Gaussian shift-in-variance transient problem

圖3 當(dāng)=104，Pd=0.6時(shí)快速傅里葉變換方法計(jì)算所得檢測門限和方差隨瞬態(tài)信號(hào)長度的變化Fig.3 Illustration of the threshold and variance calculated for different transient signal lengths with =104and Pd=0.6.

3 水池實(shí)驗(yàn)和海試數(shù)據(jù)驗(yàn)證

圖4為在消聲水池中采集到的重物落水瞬態(tài)信號(hào)，信號(hào)主要能量分布在0.2～0.5 s，頻率分布在100～210 Hz范圍內(nèi)，采樣率為800 Hz，瞬態(tài)信號(hào)長度為200個(gè)點(diǎn)。將0.3 s的信號(hào)以－10 dB疊加到高斯白噪聲上，信號(hào)持續(xù)時(shí)間為0.725 s～1.025 s。根據(jù)圖3的計(jì)算結(jié)果，算法參數(shù)設(shè)定為門限47.499，偏差1.120 9。圖5累積和檢驗(yàn)結(jié)果中在0.725～1.025 s處信號(hào)清晰，可見該方法能有效檢測低信噪比水聲瞬態(tài)信號(hào)。

圖4 水池采集信號(hào)濾波后波形Fig.4 The transient signal waveform collected after filtering

圖5 瞬態(tài)信號(hào)按－10 dB重新加入后的波形及累積和檢驗(yàn)統(tǒng)計(jì)量輸出波形Fig.5 The detected waveform with transient signal of－10 dB and the result of the page test

2014年6月舟山外海標(biāo)定聲源拉距試驗(yàn)。圖6是波束形成后能量積分(2 s)檢測和被動(dòng)累積和檢測結(jié)果比對。

圖6 波束形成后能量積分檢測和累積和檢測海試結(jié)果比對Fig.6 Comparison between the results of the energy integral test and the Page Test using the data collected in the sea

聲源船距離15.97 km，聲源是脈寬為2 s的低頻脈沖信號(hào)。試驗(yàn)背景為多艘船只的穩(wěn)態(tài)噪聲干擾。可見累積和檢驗(yàn)算法能抑制穩(wěn)定噪聲干擾，有效檢測瞬態(tài)信號(hào)，信噪比增益明顯。

4 結(jié)論

本文基于馬爾可夫狀態(tài)轉(zhuǎn)移思想，推導(dǎo)了被檢測序列在高斯分布方差變化假設(shè)下快速傅里葉變換法量化值域確定方法、更新量概率分布律的計(jì)算方法和檢測概率的遞歸計(jì)算方法。在驗(yàn)證部分，將基于快速傅里葉變換法的檢測概率理論計(jì)算結(jié)果和蒙特卡洛仿真實(shí)驗(yàn)結(jié)果相比較，吻合度好。水池和海試數(shù)據(jù)的驗(yàn)證結(jié)果表明，利用本文建立的性能分析模型尋優(yōu)搜索出的門限和方差，累積和檢驗(yàn)方法能有效檢測低信噪比的水聲瞬態(tài)信號(hào)。

當(dāng)然對水聲瞬態(tài)信號(hào)而言，高斯分布方差發(fā)生變化的假設(shè)顯然不夠準(zhǔn)確，對各種瞬態(tài)信號(hào)更為準(zhǔn)確的概率密度模型的建立將有助于提高檢測性能。這也是下一步的研究方向。

［1］呂俊軍，吳國清，杜波.非高斯水聲瞬態(tài)信號(hào)Power-Law檢測［J］.聲學(xué)學(xué)報(bào)，2004，29(4):359-362.

LYU Junjun，WU Guoqing，DU Bo.Non-Gaussian underwater transient signals detection using power-law detector［J］.Acta Acoustica，2004，29(4):359-362.

［2］CORNEL L，ANDRE Q.Transient signal detection using overcomplete wavelet transform and high-order statistics［C］//ICASSP2003.［S.l.］，2003.

［3］呂成剛.水下瞬態(tài)信號(hào)特性獲取與分析［D］.哈爾濱:哈爾濱工程大學(xué)，2009:44-51.

LYU Chenggang.Characteristic acquiring and analyzing of underwater transient signal［D］.Harbin:Harbin Engineering University，2009:44-51.

［4］ABRAHAM A D，WILLETT P.Active sonar detection in shallow water using the Page Test［J］.IEEE Trans on Aerospace and Electronic System，2002，33:1225-1229.

［5］WANG Zhen.New approaches to transient detection and signal segmentation［D］.Storrs:University of Connecticut，2002:72-83.

［6］游波，蔡志明.累積和檢驗(yàn)算法中的反饋機(jī)制研究［J］.電子學(xué)報(bào)，2010，38(6):1434-1437.

YOU Bo，CAI Zhiming.On feedback mechanism of the page test［J］.Acta Electronica Sinica，2010，38(6):1434-1437.

［7］游波，蔡志明.一種估計(jì)主動(dòng)聲吶回波擴(kuò)展時(shí)間的有效方法［J］.電子學(xué)報(bào)，2012，40(12):2223-2556.

YOU Bo，CAI Zhiming.An effective method to estimate spreading time of active sonar echoes［J］.Acta Electronica Sinica，2012，40(12):2223-2556.

［8］WANG Zhen，WILLETT P.A performance study of some transient detectors［J］.IEEE Transaction on Signal Processing，2000，48(9):2682-2685.

［9］游波，蔡志明.累積和算法應(yīng)用于主動(dòng)聲吶檢測時(shí)的性能分析［J］.海軍工程大學(xué)學(xué)報(bào)，2009，21(6):80-83.

YOU Bo，CAI Zhiming.Analysis of performance of the Page's Test as used to detect active sonar signals［J］.Journal of Naval University of Engineering，2009，21(6):80-83.

［10］HAN Chunming，WILLETT P.Some methods to evaluate the performance of Page'sTest as used to detect transient signals［J］.IEEE Trans on Signal Processing，1999，47: 2112-2127.

［11］濮曉龍.關(guān)于累積和(CUSUM)檢驗(yàn)的改進(jìn)［J］.應(yīng)用數(shù)學(xué)學(xué)報(bào)，2002，26(2):225-241.

PU Xiaolong.The improvement of the CUSUM test［J］.Acta Mathematicae Applicatae Sinica，2002，26(2):225-241.

［12］盛驟，謝式千.概率論與數(shù)理統(tǒng)計(jì)［M］.北京:高等教育出版社，2008:36-43.

SHENG Ju，XIE Shiqian.Probability and mathematical statistics［M］.Beijing:Higher Education Press，2008:36-43.