陳陽，趙安邦，王自娟，惠俊英

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

陣列信號處理誕生幾十年來，已發(fā)展出許多成熟、系統(tǒng)的理論和技術(shù)［1-2］，但由于實(shí)際聲吶環(huán)境的復(fù)雜性，仍存在需要改進(jìn)的地方.本文致力于提高基陣對線譜目標(biāo)的檢測能力，仿真和海試表明改進(jìn)的方法能有效檢測強(qiáng)干擾背景下線譜目標(biāo).

聲吶視野中通常存在多個目標(biāo)，波束主瓣接收了期望目標(biāo)的信號，同時其旁瓣也接收了其他目標(biāo)的信號，成為探測期望目標(biāo)的干擾，稱為“多目標(biāo)干擾”.此時，理論上最佳陣列處理器是最小方差無畸變(minimum variance distortionless response，MVDR)波束形成器［3］.

對于線譜，瞬時頻率方差(variance of instantaneous frequency，VIF)［4］的估計值可作為檢測統(tǒng)計量用于信號檢測，在理想信道白噪聲背景下的檢測性能略低于最佳處理器(匹配濾波器)，但差別很小.瞬時頻率方差檢測器的虛警率與干擾背景的能量起伏無關(guān)，因而在非平穩(wěn)的干擾背景和信道中具有很好的穩(wěn)健性［5］.

水中航行器輻射噪聲中線譜的譜級信噪比通常較連續(xù)譜高許多［6］，如果充分利用這一特性，可以有效提高系統(tǒng)的檢測性能.寬帶信號分頻帶處理進(jìn)行融合時，按各頻帶的信噪比進(jìn)行加權(quán)是最優(yōu)的［7］.然而在實(shí)際應(yīng)用中，線譜頻率未知，各頻帶的信噪比也無法得知，因而一般假設(shè)各頻帶信噪比相同，即各頻帶的權(quán)相同;而由于線譜能量在寬帶信號能量中所占比例很小，所形成的高性能空間譜被其他低信噪比頻帶形成的空間譜淹沒.解決這一問題的一個途徑是將各頻帶獨(dú)立顯示而不進(jìn)行寬帶融合，這就需要四維顯示，但目前仍沒有十分有效的四維顯示方法.注意到，如果目標(biāo)有線譜，該方位波束的瞬時頻率方差很小甚至為零，否則就較大，利用頻率方差對寬帶方位譜進(jìn)行加權(quán)，線譜目標(biāo)方位的譜強(qiáng)度得到增強(qiáng);而且由于不需顯示每個頻帶的結(jié)果，可以通過常規(guī)的三維方位時間歷程瀑布圖顯示.

1 MVDR波束形成

假設(shè)陣列由一組全向陣元組成，陣元的位置為pn(n=0，1，…，M－1，M為陣元數(shù)).陣列接收的快拍由信號和噪聲組成，在頻域表示為X(ω)= XS(ω)+N(ω)，其中信號矢量可以寫為

式中:N(ω)為噪聲矢量;F(ω)是源信號的頻域快拍;v(ω∶κS)=［e－jκTSp0，e－jκTSp1，…，e－jκTSpM－1］T是一個波數(shù)為κS的平面波對應(yīng)的陣列流形矢量，κS=為單位矢量，λ為對應(yīng)于頻率ω的波長.噪聲快拍N(ω)是零均值隨機(jī)矢量，其互譜密度譜矩陣為Sn(ω)=I.陣列輸入的互譜密度矩陣為

式中:Ss(ω)為期望目標(biāo)的互譜密度矩陣，Sc(ω)為多目標(biāo)干擾的互譜密度矩陣.利用權(quán)向量WH(ω)對X(ω)進(jìn)行處理，波束輸出為Y(ω)=WH(ω)X(ω).

對于寬帶信號，分成若干頻帶，各頻點(diǎn)的頻域快拍矢量Xi(ωj)，i=0，…，L－1，j=0，…，J－1，L為快拍數(shù)，J為頻點(diǎn)數(shù).各頻點(diǎn)互譜密度矩陣通過快拍的時域平均得到其估計值:

MVDR波束形成最優(yōu)權(quán)向量為

空間譜為

J個頻點(diǎn)空間譜融合［8］:

2 導(dǎo)向最小方差波束形成

STMV是一種基于稱為導(dǎo)向協(xié)方差矩陣(steered covariance matrix，STCM)的空時統(tǒng)計量的最優(yōu)波束形成方法.STCM的自由度等于時間帶寬積，因而只需要頻點(diǎn)數(shù)與獨(dú)立快拍數(shù)的積不小于陣元數(shù)即可達(dá)到滿秩可逆，對于寬帶信號，其收斂遠(yuǎn)快于MVDR.

指向協(xié)方差矩陣:

式中:

最優(yōu)權(quán)為

式中:I為M×1的1向量.空間譜為

顯然，當(dāng)J=1時，pstmv(θ)=pwmvdr(κs)，即當(dāng)信號為窄帶時，STMV與MVDR等價.

3 瞬時頻率方差

文獻(xiàn)［5］已經(jīng)證明，對于線譜，瞬時頻率方差是十分有效的檢測量.窄帶高斯白噪聲的瞬時頻率方差為有效帶寬的平方，顯著大于線譜信號的瞬時頻率方差.

文獻(xiàn)［5］還比較研究了多種瞬時頻率方差估計方法:希爾伯特變換類的解析信號相位差分瞬時頻率估計方法、過零頻率估計器及其內(nèi)插修正算法、窄帶信號的自適應(yīng)頻率估計.由于寬帶信號陣列處理采用頻域處理，而上述方法均不適合用于頻域處理，因而下面著重分析基于短時傅里葉時頻分析的瞬時頻率方差估計方法.

噪聲中的線譜信號表示為

式中:As和φ0分別為線譜信號的幅度和初相位，信號帶寬為B，噪聲功率為.將信號以采樣率fs采樣、分段，每段長度為T(單位為s)做離散傅里葉變換.每段信號的離散傅里葉變換可以表示為

式中:Ss、Sn分別表示線譜信號和噪聲分量;X、Y分別為Sx(k)的實(shí)部和虛部.頻譜Sx(k)中幅值最大的頻點(diǎn)即為瞬時頻率.根據(jù)文獻(xiàn)［12］時，X滿足正態(tài)分布滿足正態(tài)分布時，X和Y均滿足正態(tài)分布于是，其包絡(luò) z=|Sx(k)|=滿足萊斯分布［11］:

除線譜外其余所有頻點(diǎn)的最大值x=max{|Sx(k)|，k≠k0}的分布函數(shù)為

其概率密度函數(shù)為

y=z－x的分布函數(shù)為

于是z＜x的概率，即線譜頻率幅值最大的概率為

其他頻點(diǎn)最大的概率為

瞬時頻率的均值為

瞬時頻率的方差為

通過蒙特卡洛仿真驗(yàn)證上述理論.帶寬0.5～1 kHz的高斯白噪聲，線譜頻率800 Hz，F(xiàn)FT窗長1 s，20 000次獨(dú)立統(tǒng)計.和瞬時頻率方差的理論值式(18)和蒙特卡洛仿真統(tǒng)計結(jié)果如圖1和圖2，兩者完全一致.由圖中可以看到，當(dāng)線譜譜級信噪比超過一定門限時，瞬時頻率方差為零.

圖1 線譜頻率幅值最大的概率f(k0)Fig.1 Probability of line spectrum to be the peak f(k0)

圖2 短時傅里葉變換估計瞬時頻率方差Fig.2 Variance of instantaneous frequency through STFT

4 頻率方差加權(quán)STMV波束形成檢測器

由上述可知，如果目標(biāo)信號有線譜，只要波束輸出中線譜譜級信噪比超過門限，其瞬時頻率方差很小甚至為零，而如果目標(biāo)信號為連續(xù)譜白噪聲則瞬時頻率方差較大，這種區(qū)別可以用于線譜目標(biāo)檢測.用各方位波束輸出的瞬時頻率方差作為權(quán)對陣列的STMV方位譜進(jìn)行調(diào)整的檢測器稱為瞬時頻率方差加權(quán)STMV波束形成檢測器(VIF-STMV).如圖3所示，陣列數(shù)據(jù)經(jīng)滑動窗分成相互重疊的快拍，STMV波束形成輸出分別進(jìn)行能量積分和周期圖譜估計;能量積分估計空間譜;周期圖譜估計得到波束輸出的功率譜，峰選確定最大譜線的頻率，并統(tǒng)計該頻率的方差，以方差的倒數(shù)作為系數(shù)對空間譜進(jìn)行加權(quán).

圖3 VIF-STMV框圖Fig.3 Diagram of variance of instantaneous frequencybeamforming detecting

不同的加權(quán)方法可以獲得不同的效果，下面給出2種加權(quán)方法，倒數(shù)加權(quán)和指數(shù)加權(quán):

式中:pstmv為STMV波束形成空間譜，σf為瞬時頻率方差，Δ為一小量，它保證當(dāng)瞬時頻率方差為零時倒數(shù)權(quán)存在.

以陣元間距0.3 m的48元直線陣為例，112°方向存在一信噪比為15 dB強(qiáng)干擾(干擾1)，13°和60°方向分別存在信噪比－13 dB、－18 dB兩干擾(干擾2、3)73°方向?yàn)樾旁氡龋?8 dB的目標(biāo).其中干擾3和目標(biāo)均含有1根線譜，分別為650 Hz和700 Hz，其譜級信噪比較連續(xù)譜高15 dB.信號頻帶為0.5～1 kHz，快拍長度為1 s，方位譜估計的積分時間為4s(4個獨(dú)立快拍)，功率譜估計的滑動步長為0.125 s.圖4中實(shí)線為STMV方位譜，可以看到由于分辨力不夠，無法分開干擾3和目標(biāo)，也就是說，干擾3的存在影響了期望目標(biāo)的檢測.虛線為倒數(shù)加權(quán)的結(jié)果，由于干擾3和目標(biāo)含有線譜，該方向方位譜加權(quán)較大，因而能清晰分辨干擾3和目標(biāo)，干擾1和干擾2由于不含線譜，方位譜強(qiáng)度被抑制.點(diǎn)線為指數(shù)加權(quán)的結(jié)果，可以看到干擾1和干擾2方位譜強(qiáng)度被進(jìn)一步抑制，旁瓣更平滑，更穩(wěn)健.

圖4 VIF-STMV仿真Fig.4 VIF-STMV suppress multisource interference

5 海試數(shù)據(jù)處理

為了驗(yàn)證該方法在工程應(yīng)用中的有效性，在東海進(jìn)行了海試，試驗(yàn)采用的接收陣為陣元間距0.3 m的48元拖曳線列陣.目標(biāo)由遠(yuǎn)處以較快速度接近拖船.試驗(yàn)在漁期進(jìn)行，加上試驗(yàn)海區(qū)靠近航道，因而在附近海域存在較多的其他船只的干擾.

圖5 STMV波束形成方位歷程Fig.5 Bearing-time chart of STMV

圖5為STMV波束形成方位歷程圖，從圖中可以清晰看到至少6個目標(biāo)的歷程，期望目標(biāo)和主要的幾個干擾如圖中標(biāo)示.90 s后拖船轉(zhuǎn)向，左右舷的目標(biāo)分別向船艏和船尾方向偏移.前90 s，由于期望目標(biāo)與干擾1方位重疊，目標(biāo)被強(qiáng)干擾淹沒，無法檢測到目標(biāo).拖船轉(zhuǎn)向后目標(biāo)與該強(qiáng)干擾方位分開，可以較為明顯的檢測到目標(biāo).

圖6 倒數(shù)加權(quán)VIF-STMV的方位歷程Fig.6 Bearing-time chart of VIF-STMV (multiplicative inverse)

圖6為倒數(shù)加權(quán)VIF-STMV的方位歷程圖結(jié)果，圖7為指數(shù)加權(quán)VIF-STMV的方位歷程圖結(jié)果.從圖中可以明顯看到目標(biāo)，90 s前位于70°附近，90 s后偏向0°方向.指數(shù)加權(quán)結(jié)果較倒數(shù)加權(quán)結(jié)果背景更干凈.軌跡出現(xiàn)斷裂是由于淺海聲信道濾波效應(yīng)引起的:淺海聲信道是一梳狀濾波器，子通帶隨信道的變化而變化，當(dāng)線譜頻率處于子止帶時線譜信噪比顯著降低使得瞬時頻率方差變大.

圖7 指數(shù)加權(quán)VIF-STMV的方位歷程Fig.7 Bearing-time chart of VIF-STMV (exponential weight)

進(jìn)一步分析各目標(biāo)信號的功率譜.圖8為70°方向波束時頻LOFAR圖，可以看到90 s前的頻譜存在線譜，為了確認(rèn)該線譜屬于目標(biāo)，分析56°方向和90°方向波束時頻LOFAR圖，如圖9和圖10.56°方向波束時頻LOFAR圖出現(xiàn)該線譜，說明90 s之后目標(biāo)進(jìn)入該波束，而90°方向波束時頻LOFAR圖90 s之后頻譜強(qiáng)度變大但沒有線譜，可見是干擾1進(jìn)入該波束，而且干擾1的頻譜在該工作頻段內(nèi)無線譜.圖11為20°方向(干擾2方位)波束的時頻LOFAR圖，頻譜無線譜，因而盡管干擾2輻射噪聲很強(qiáng)，在VIF-STMV方位歷程圖上仍然被抑制.圖12為160°方向(干擾3方位)波束的時頻LOFAR圖，頻譜存在一根線譜，因而干擾3在VIF-STMV方位歷程圖上被增強(qiáng).

圖8 70°方向波束時頻LOFARFig.8 LOFAR at 70°

圖9 56°方向波束時頻LOFARFig.9 LOFAR at 56°

圖10 90°方向波束時頻LOFARFig.10 LOFAR at 90°

圖11 20°方向波束時頻LOFARFig.11 LOFAR at 20°

圖12 160°方向波束時頻LOFARFig.12 LOFAR at 160°

6 結(jié)束語

文中提出了頻率方差加權(quán)導(dǎo)向最小方差波束形成檢測器，利用目標(biāo)輻射噪聲中線譜有較高的強(qiáng)度和穩(wěn)定度這個特征，用每個方位的波束輸出的頻率方差，對方位譜進(jìn)行加權(quán)，得到線譜目標(biāo)的方位-時間歷程圖，線譜目標(biāo)方位波束因輸出信號的瞬時頻率方差較小得到增強(qiáng)，而其他方位波束被抑制.該檢測器可有效抗多目標(biāo)干擾，實(shí)現(xiàn)在強(qiáng)相干干擾中檢測到弱線譜目標(biāo)，且只須三維顯示，避免了通常的線譜檢測器須四維顯示、觀察費(fèi)力的困擾.仿真和海試結(jié)果表明，在多目標(biāo)、強(qiáng)干擾的環(huán)境下，頻率方差加權(quán)導(dǎo)向最小方差波束形成檢測器可遠(yuǎn)程探測線譜目標(biāo)，與方位譜能量檢測器比較，有更好的探測性能.

［1］Van TREES H L.Optimum array processing［J］.New York:John Wiley＆Sons Inc，2003:3-5.

［2］王永良，陳輝彭，應(yīng)寧，萬群.空間譜估計理論與算法［M］.北京:清華大學(xué)出版社，2004:2-8.

WANG Yongliang，CHEN Huipeng，YING Ning，WAN Qun.Theories and algorithms of spectrum estimation［M］.Beijing:Tsinghua University Press，2004:2-8.

［3］CAPON J.High-resolution frequency-wavenumber spectrum analysis［J］.Proceedings of the IEEE，1969，57(8): 1408-1418.

［4］林茂庸，柯有安.雷達(dá)分辨理論［M］.北京:國防工業(yè)出版社，1984:134-141.

LIN Maoyong，KE Youan.Rada resolution theory［M］.Beijing:National Defence Industrial Press，1984:134-141.

［5］梁國龍.回波信號瞬時參數(shù)序列分析及其應(yīng)用研究［D］.哈爾濱:哈爾濱工程大學(xué)，1997:80-82.

LIANG Guolong.A study on instananeous parameter sequence analysis of echoes and its application［D］.Harbin: Harbin Engineering University，1997:134-141.

［6］吳國清，李靖.艦船噪聲識別-線譜穩(wěn)定性和唯一性［J］.聲學(xué)學(xué)報，1999，24(1):6-11.

WU Guqing，LI Jing.Ship radiated-noise recognition-stability and uniqueness of line spectrum［J］.Acta Acustica，1999，24(1):6-11.

［7］HUNG H，KAVEH M.Focussing matrices for coherent signal-subspace processing［J］.IEEE Trans ASSP，1988，36 (8):1272-1281.

［8］KROLIK J，SWINGLER D.focused wide-band array processing by spatial resampling［J］.IEEE Trans ASSP，1990，38(2):356-360.

［9］LIU W L，DING S X.An efficient method to determine the diagonal loading factor using the constant modulus feature［J］.IEEE Trans SP，2008，56(12):6102-6106.

［10］KROLIK J，SWINGLER D.Multiple broadband source location using steered covariance matrices［J］.IEEE Trans ASSP，1989，37(10):1481-1494.

［11］朱華，黃輝寧，李永慶，梅文博.隨機(jī)信號分析［M］.北京:北京理工大學(xué)出版社，1990:317-319.

ZHU Hua，HUANG Huining，LI Youqing，MEI Wenbo.Random signal analysis［M］.Beijing:Beijing Institute of Technology Press，1990:317-319.

［12］WHALEN A D.Signal detecting in noise［M］.Beijing: Science Press，2006:77-78.