戴文舒,陳亞,陳新華,孫長(zhǎng)瑜,余華兵

（1.中國(guó)科學(xué)院聲學(xué)研究所,北京 100190；2.中國(guó)科學(xué)院大學(xué),北京 100190；3.海軍駐無(wú)錫地區(qū)軍事代表室,江蘇無(wú)錫 214061）

差分二次平均修正的頻域相位補(bǔ)償線譜檢測(cè)方法

戴文舒1,2,陳亞3,陳新華1,孫長(zhǎng)瑜1,余華兵1

（1.中國(guó)科學(xué)院聲學(xué)研究所,北京 100190；2.中國(guó)科學(xué)院大學(xué),北京 100190；3.海軍駐無(wú)錫地區(qū)軍事代表室,江蘇無(wú)錫 214061）

針對(duì)常規(guī)平均功率譜方法幅度平方后相位信息丟失造成檢測(cè)性能損失問(wèn)題,通過(guò)對(duì)觀測(cè)數(shù)據(jù)分段進(jìn)行快速傅里葉變換后的統(tǒng)計(jì)特性進(jìn)行分析,提出了高斯噪聲背景下檢測(cè)未知線譜的頻域相位補(bǔ)償方法。采用蒙特卡洛法估計(jì)相位補(bǔ)償因子,先進(jìn)行差分法去除零星野值,然后使用二次平均（TPM）去除連續(xù)野值,最后使用估計(jì)的相位補(bǔ)償因子均值構(gòu)造廣義似然比檢驗(yàn)統(tǒng)計(jì)量實(shí)現(xiàn)檢測(cè)。從理論上對(duì)比了廣義平均周期圖（AVGPR）法、廣義功率譜（GPR）法及文中方法的檢測(cè)性能。仿真結(jié)果表明,文中方法充分利用信號(hào)的相位信息,參數(shù)估計(jì)簡(jiǎn)單,相位補(bǔ)償因子估計(jì)準(zhǔn)確,在一定虛警概率下較AVGPR法的檢測(cè)性能提高了接近5 dB左右。

聲學(xué)；快速傅里葉變換分析；線譜檢測(cè)；頻域相位補(bǔ)償；檢驗(yàn)統(tǒng)計(jì)量；功率譜

0 引言

Urick指出在一定假設(shè)下,海洋背景噪聲服從高斯分布[1]。高斯背景下未知線譜的檢測(cè)近年來(lái)一直被關(guān)注,隨著隱身技術(shù)的不斷提高,水下目標(biāo)輻射線譜信號(hào)越安靜,對(duì)被動(dòng)聲納的檢測(cè)性能要求越苛刻。功率譜分析與低頻分析和記錄（LOFAR）結(jié)合來(lái)檢測(cè)線譜是使用最普遍的方法,然而該方法的檢測(cè)模型中功率譜的計(jì)算只用到了線譜信號(hào)的幅度信息,相位信息在平方觀測(cè)值的離散傅里葉變換（DFT）后丟失。相位信息的丟失會(huì)降低被動(dòng)聲納檢測(cè)的性能。針對(duì)相位信息的補(bǔ)償,國(guó)內(nèi)外的相關(guān)學(xué)者進(jìn)行了一系列研究,大多從時(shí)域平均角度出發(fā)。沈國(guó)際等指出當(dāng)信號(hào)頻率未知或由采樣頻率的關(guān)系無(wú)法實(shí)現(xiàn)對(duì)信號(hào)的整數(shù)倍周期截取,會(huì)導(dǎo)致性能大大下降,甚至反而降低信噪比[2]。陳紹華等提出在頻率已知的條件下,對(duì)每段截取信號(hào)進(jìn)行時(shí)延補(bǔ)償,從而得到無(wú)相位誤差的時(shí)域平均方法[3]。而對(duì)于頻域未知的情況,則補(bǔ)償時(shí)需要進(jìn)行頻率搜索,計(jì)算量大,為此,陳新華等提出了搜索頻率與信號(hào)頻率無(wú)關(guān),僅是截取長(zhǎng)度倒數(shù)的方法,大大減小了計(jì)算量,根據(jù)不同的搜索頻率對(duì)每段截取信號(hào)進(jìn)行相位補(bǔ)償,使時(shí)域平均后的輸出功率到達(dá)最大[4]。Chun等提出的最優(yōu)線譜廣義似然比檢測(cè)（PGLRT）改進(jìn)了廣義平均周期圖（AVGPR）的幅度功率譜檢測(cè)缺陷,充分利用了功率譜的分布特性,但是其檢驗(yàn)統(tǒng)計(jì)量的估計(jì)值參數(shù)遇到了低信噪比時(shí)的虛數(shù)解,導(dǎo)致檢測(cè)結(jié)果發(fā)散[5]。

本文從頻域功率譜角度出發(fā),先對(duì)二元統(tǒng)計(jì)模型下功率譜分布特性進(jìn)行分析,然后研究分段快速傅里葉變換（FFT）后相位因子的存在帶來(lái)的影響及補(bǔ)償方法及估計(jì)方法,進(jìn)而使用估計(jì)得到的相位影響因子使用最大似然準(zhǔn)則構(gòu)造檢驗(yàn)統(tǒng)計(jì)量,從模型出發(fā)理論分析AVGPR、廣義功率譜（GPR）及本文方法的檢測(cè)性能,最后給出了算法的仿真結(jié)果和主要結(jié)論。

1 檢測(cè)問(wèn)題和模型

檢測(cè)問(wèn)題可以描述為二元假設(shè)檢驗(yàn),

式中：s（n）=Aexp（j（w0n+φ））是待檢測(cè)信號(hào),A是信號(hào)的未知幅度,w0是未知?dú)w一化角頻率,φ∈（0, 2π）是初始化相位,未知但是固定；g（n）是方差為的高斯噪聲。將x（n）分成長(zhǎng)度為N的L段,總長(zhǎng)度為Q,R與重疊個(gè)數(shù)有關(guān),第l段表示為

上式表明在最接近信號(hào)的數(shù)字頻率單元觀測(cè)值的功率將最大。由于信號(hào)頻率未知,需要通過(guò)搜索來(lái)確定頻率單元。這顯然計(jì)算量大,不可取。

記信號(hào)部分為ρl（k）,可以看出Δ的存在帶來(lái)的相移,所以需要進(jìn)行補(bǔ)償修正。

對(duì)于噪聲頻率單元,k≠k0,信號(hào)成分很小, Zl（k）服從0均值、方差σ2的高斯分布；對(duì)信號(hào)頻率單元Zl（k）服從μl均值、方差σ2的高斯分布,其中

為使用相位信息,Δ的估計(jì)很重要,從（11）式看到,Δ與2Δ（（l-1）R+N-1）相乘,每段的相移不同,為了充分利用各段的相位信息,由Δ引起的相移需要補(bǔ)償。

2 Δ因子的估計(jì)

補(bǔ)償相移的關(guān)鍵是對(duì)Δ因子進(jìn)行估計(jì)。文獻(xiàn)[5]中指出了用解非線性方程的方法估計(jì)Δ,涉及多解問(wèn)題和初始化猜想。而通過(guò)對(duì)Zl（k）進(jìn)行L點(diǎn)DFT分析：

3 差分二次平均方法去野值

由中心極限定理,滿足任何分布的P個(gè)獨(dú)立觀測(cè)值的均值服從正態(tài)分布[7],所以采用蒙特卡洛仿真500次,由仿真結(jié)果^Δ的估值會(huì)出現(xiàn)野值。為了得到準(zhǔn)確的估計(jì)結(jié)果,要先對(duì)這些數(shù)據(jù)進(jìn)行預(yù)處理。

根據(jù)存在野值的特點(diǎn)分為零星野值和連續(xù)野值,它們的存在會(huì)使得^Δ的均值的估計(jì)值偏離真實(shí)值,所以要進(jìn)行判斷和去除。零星野值點(diǎn)與它左右領(lǐng)域點(diǎn)的差值都很大,如圖1中第i次估計(jì)；若待研究點(diǎn)與其左右領(lǐng)域點(diǎn)的差值只有一個(gè)很大,那么既可以是野值點(diǎn),也可以是臨界點(diǎn),如圖1中第i-1次估計(jì)。

圖1 零星野值和連續(xù)野值Fig.1 Sporadic and continuous outliers

文獻(xiàn)[8]比較了估計(jì)背景均值的4種方法。本文采用廣泛使用的二次平均（TPM）估計(jì)均值,再根據(jù)萊特準(zhǔn)則去除連續(xù)野值。下面方法描述TPM技術(shù),μk是第k個(gè)單元的估計(jì)均值,它是其領(lǐng)域Ωk輸出的函數(shù),位于中心,K=2M+1窗長(zhǎng)。

可見μk是一系列值,計(jì)算μk的均值μ和方差σ2.采用萊斯準(zhǔn)則判斷去零星野值后的3σ,則去除該點(diǎn)。最后用剩余的可信估計(jì)來(lái)估計(jì)Δ.

4 由估計(jì)值^Δ補(bǔ)償各段相移

由（10）式看到在（9）式兩邊乘以因子exp（-jπ· Δ（2（l-1）R+N-1））來(lái)補(bǔ)償每段的相移得到Y(jié)l（k）,其表達(dá)式如（20）式。

在兩種假設(shè)下的聯(lián)合概率密度函數(shù)（PDF）分別為

式中：未知參數(shù)為A、σ2.構(gòu)造對(duì)數(shù)似然比

在假設(shè)H0下,未知參數(shù)為σ2,最大似然估計(jì)^σ20可由（24）式得到。

綜上,基于差分TPM修正的頻域相位補(bǔ)償算法流程圖如圖2所示。

圖2 基于差分TPM修正的頻域相位補(bǔ)償流程圖Fig.2 Flow chart of frequency domain phase compensation using difference and TPM

5 性能分析

文獻(xiàn)[9]分析了平均周期圖法的檢測(cè)性能。以下列出主要公式：

Matlab工具箱中提供了計(jì)算這些概率值的函數(shù),仿真可直接調(diào)用。

6 仿真性能曲線

分別根據(jù)（36）式、（42）式畫出AVGPR、本文方法的受制者工作特性（ROC）檢測(cè)性能曲線,GPR是自由度為2的AVGPR.本文方法記作RGLRT.圖3表示信噪比為-25 dB、L=32下的理論ROC曲線。圖4表示信噪比為-30 dB下的理論ROC曲線。圖5為虛警概率PFA=10-3下的檢測(cè)概率PD隨信噪比SNA的變化規(guī)律。

理論分析RGLRT算法比AVGPR算法在PFA= 10-3下檢測(cè)性能增加了提高了近5 dB左右。而眾所周知,AVGPR算法相比GPR算法的性能提高是由分段累積帶來(lái)的。

7 仿真驗(yàn)證

圖3 信噪比為-25 dB下的ROC性能曲線Fig.3 Detection performance for SNR=-25 dB

圖4 信噪比為-30 dB下的ROC性能曲線Fig.4 Detection performance for SNR=-30 dB

圖5 虛警概率一定下的ROC性能曲線Fig.5 Detection performance at a certain false alarm rate

1）仿真條件：采樣頻率fs=512 Hz,觀測(cè)總數(shù)Q=20 000,分段數(shù)L=39,每段長(zhǎng)度512,蒙特卡羅仿真次數(shù)P=500,取TPM窗長(zhǎng)度K=3,信噪比-29 dB.仿真信號(hào)頻率變化,對(duì)比了真實(shí)Δ與去野值前后估計(jì)值的結(jié)果。

由表1可知,未去野值時(shí)^Δ的估計(jì)均值較去野值后的^Δ的估計(jì)均值與真實(shí)Δ的誤差較大,會(huì)導(dǎo)致各段相位補(bǔ)償?shù)牟粶?zhǔn)確,影響檢測(cè)性能,而且未去野值的的估計(jì)方差要比去野值后的大。差分TPM提高了相移估計(jì)因子的估計(jì)準(zhǔn)確性。

表1 信噪比-29 dB下估計(jì)均值和真實(shí)值Δ的仿真結(jié)果Tab.1 The simulation results of estimated and true values for SNR=-29 dB

表1 信噪比-29 dB下估計(jì)均值和真實(shí)值Δ的仿真結(jié)果Tab.1 The simulation results of estimated and true values for SNR=-29 dB

信號(hào)頻率f*0/Hz真實(shí)Δ值方差1/3（170.6）-6.510 4×10-4-5.974 6×10-4-6.526 6×10-42.346 5×10-45.347 4×10未去野值時(shí)^Δ值去野值時(shí)^Δ值未去野值時(shí)^Δ方差去野值時(shí)^Δ -5 1/4.1（124.8）-2.381 9×10-4-2.072 3×10-4-2.498 0×10-42.396 4×10-44.248 2×10-51/5（102.4）7.812 5×10-44.408 1×10-47.816 7×10-45.369 6×10-41.833 9×10-41/6（85.3）6.510 4×10-45.999 6×10-46.524 9×10-42.418 1×10-48.691 4×10-51/7（73.1）2.790 2×10-41.753 8×10-42.891 2×10-43.905 0×10-47.659 5×10-51/8（64）0-1.201 9×10-603.168 2×10-50 1/9（56.9）-2.170 1×10-4-1.681 7×10-4-2.028 5×10-42.573 9×10-43.939 8×10-5

由第3節(jié)信噪比-30 dB下的ROC理論分析可知,在虛警概率PFA=10-3下AVGPR算法的檢測(cè)概率為0.3,而RGLRT算法的檢測(cè)概率為1,PGLRT算法不能實(shí)現(xiàn)檢測(cè)。

圖6 AVGPR算法下的檢測(cè)結(jié)果Fig.6 Detected result of AVGPR algorithm

8 結(jié)論

圖7 RGLRT算法下的檢測(cè)結(jié)果Fig.7 Detected result of RGLRT algorithm

圖8 PGLRT算法下的檢測(cè)結(jié)果Fig.8 Detected result of PGLRT algorithm

1）對(duì)觀測(cè)數(shù)據(jù)分段進(jìn)行FFT后的功率譜二元統(tǒng)計(jì)模型及分布進(jìn)行研究,得出了各段的AVGPR在沒(méi)有線譜信號(hào)時(shí)服從自由度為2L的中心chi平方分布,而線譜存在時(shí)服從自由度為2L的非中心chi平方分布,非中心參數(shù)為λ=N本文基于頻域相位補(bǔ)償?shù)乃迫槐葯z測(cè)統(tǒng)計(jì)量在沒(méi)有線譜存在時(shí)服從中心F分布,在線譜存在時(shí)服從非中心F分布,非中心參數(shù)為λ=NGPR是自由度為2的AVGPR.

2）相同分段數(shù)下,理論分析本文算法比AVGPR算法在PFA=10-3下檢測(cè)性能提高了近5 dB左右,而且F分布較chi平方分布在達(dá)到同樣檢測(cè)概率時(shí)的檢測(cè)門限低,所以性能更優(yōu)。

3）基于前后向差分去除零星野值后采用TPM估計(jì)均值方法,有效地去除了連續(xù)野值,提高了相位影響因子的估計(jì)準(zhǔn)確性,提高了檢測(cè)統(tǒng)計(jì)量的檢測(cè)性能。

（References）

[1] Urick R J.Models for the amplitude of fluctuations of narrow-band signals and noise in the sea[J].J Acoustic Soc Amer,1977,62： 878-887.

[2] 沈國(guó)際,陶利民,徐永成.時(shí)域同步平均的相位誤差累積效益研究[J].振動(dòng)工程學(xué)報(bào),2007,20（4）：335-339.

SHEN Guo-ji,TAO Li-min,XU Yong-cheng.Research on phase error accumulation effect of time synchronous averaging[J].Journal of Detection and Control,2007,20（4）：335-339.（in Chinese）

[3] 陳紹華,相敬林.一種改進(jìn)的時(shí)域平均法檢測(cè)微弱信號(hào)研究[J].探測(cè)與控制學(xué)報(bào),2003,25（4）：56-59.

CHEN Shao-hua,XIANG Jing-lin.Amodified time averaging method in weak signal detection[J].Journal of Detection and Control,2003,25（4）：56-59.（in Chinese）

[4] 陳新華,孫長(zhǎng)瑜,鮑習(xí)中.基于相位補(bǔ)償?shù)臅r(shí)域平均方法[J].應(yīng)用聲學(xué),2011,30（4）：268-274.

CHEN Xin-hua,SUN Chang-yu,BAO Xi-zhong.Time domain averaging based on phase compensations[J].Applied Acoustics, 2011,30（4）：268-274.（in Chinese）

[5] Chun R W,Joo T G,Hong T C.Optimaltonal detectors based on the power spectrum[J].IEEE Journal of Oceanic Engineering, 2000,25（4）：540-552.

[6] Qing W,Chun R W.Anovel CFAR tonal detector using phase compensation[J].IEEE Journal of Oceanic Engineering,2005, 30（4）：900-911.

[7] 施久玉.概率論和數(shù)理統(tǒng)計(jì)[M].哈爾濱：哈爾濱工程大學(xué)出版社,2005.

SHI Jiu-yu.Probability andmathematics[M].Harbin：Harbin Engineering University Press,2005.（in Chinese）

[8] Willan A S,Edward D L.A performace comparison of four noise background normalization schemes proposed for signal detection systems[J].J Acoustic Soc Amer,1984.76（6）：1738-1742.

[9] Qing W,Chun W,Joo T G.Theoretical performance analysis and simulation of a GLRT tonal detector[C]∥Oceans 2001 MTS/ IEEE ConfErence and Exhibition.Honolulu,HI：IEEE,2001： 1654-1659.

[10] Kay S M.統(tǒng)計(jì)信號(hào)處理基礎(chǔ)：估計(jì)與檢測(cè)理論[M].羅鵬飛,譯.北京：電子工業(yè)出版社,2001.

Kay S M.Fundamentals of statistical signal processing：estimation theory and detection theory[M].LUO Peng-fei,translated. Beijing：Publishing House of Electronics Industry,2001.（in Chinese）

An Algorithm of Line-spectrum Detection Using Frequency Domain Phase Compensation with Difference and TPM Technique

DAI Wen-shu1,2,CHEN Ya3,CHEN Xin-hua1,SUN Chang-yu1,YU Hua-bing1
（1.Institute of Acoustics,Chinese Academy of Sciences,Beijing 100190,China； 2.University of Chinese Academy of Sciences,Beijing 100190,China； 3.Navy Military Representative Office in Wuxi Area,Wuxi 214061,Jiangsu,China）

The common average power spectrum estimation algorithm may cause the loss in detection performance because power spectrum amplitude is simply squared.The statistical properties of FFT of the segmented observations are analyzed,and a novel algorithm using frequency domain phase compensation for detecting line-spectrum in Gaussian noise is proposed.The phase compensation factor is estimated using Monte Carlo method.The sporadic and continuous outliers are removed by using difference and TPM method.Finally the statistics is checked by using generalized likelihood ratio.The detection performances of AVGPR,GPR and the proposed algorithm are theoretically compared.The simulation results show that the detection performance of the proposed algorithm is improved by 5 dB under a certain false alarm rate compared to AVGPR algorithm.It fully utilizes the phase information of the signal and has the advantagesof simple and exact parameter estimation and good robustness.

acoustics；FFT analysis；line-spectrum detection；frequency domain phase compensation； test statistics；power spectrum

TB566

A

1000-1093（2014）10-1630-08

10.3969/j.issn.1000-1093.2014.10.017

2013-08-06

國(guó)家自然科學(xué)基金項(xiàng)目（61372180）

戴文舒（1987—）,女,博士研究生。E-mail：dws.01@163.com；陳亞（1978—）,男,工程師。E-mail：cheya113@sina.com