穆森， 李京華， 張恒， 陳果

(1.西北工業(yè)大學(xué) 電子信息學(xué)院， 陜西 西安 710072; 2.陸軍研究院 近地面探測技術(shù)實(shí)驗(yàn)室， 江蘇 無錫 214035;3.宜昌測試技術(shù)研究所， 湖北 宜昌 443003)

0 引言

反潛直升機(jī)對(duì)于水下潛艇具有很大的威脅，為了提高水下潛艇對(duì)抗反潛機(jī)的能力，有必要盡快解決水下對(duì)空探測的技術(shù)問題。反潛直升機(jī)飛行時(shí)，旋翼輻射噪聲的聲波穿越空- 海界面入水后形成的頻譜具有諧波特性，為水下探測空中反潛飛機(jī)提供了可行途徑。

目前關(guān)于被動(dòng)聲探測的研究多集中在同一種介質(zhì)內(nèi)，如空氣介質(zhì)中的聲探測或海水介質(zhì)中的水下聲探測。魏麗萍等[1]給出了諧波集(HS)檢測算法，并成功地對(duì)實(shí)測的空氣中直升機(jī)聲信號(hào)進(jìn)行了HS檢測；馬寧等[2]提出了一種基于小波包和神經(jīng)網(wǎng)絡(luò)的檢測算法，可以較好地檢測識(shí)別不同型號(hào)的空氣中直升機(jī)信號(hào)；Salloum等[3]開發(fā)并建立了一個(gè)聲學(xué)空氣檢測(AAD)系統(tǒng)，實(shí)現(xiàn)了空中目標(biāo)的檢測和跟蹤，并對(duì)小型飛機(jī)、直升機(jī)和超輕飛機(jī)進(jìn)行了分類；周關(guān)林等[4]提出隨機(jī)共振水聲信號(hào)的檢測，為檢測水下微弱信號(hào)提出了新思路。關(guān)于空中飛行目標(biāo)的水下聲探測，即空氣- 水跨界聲探測問題的研究大多集中在空中聲源激發(fā)的水下聲場理論建模方面，而空中目標(biāo)的水下跨界檢測算法則研究很少。例如：以色列研制了潛艇反直升機(jī)和低空慢速飛行目標(biāo)的聲探測浮標(biāo)[5]，用于檢測海面上空目標(biāo)的聲信號(hào)；韓建輝等提出在水下檢測空中聲源的平滑偽維格納·維爾分布- 高度門限(SPWVD-HT)算法[6]、雙譜冪次法[7]，但所檢測聲信號(hào)均采用空中運(yùn)動(dòng)聲源激發(fā)的水下聲場建模仿真信號(hào)。

本文提出基于HS檢測與最強(qiáng)譜線檢測相結(jié)合的檢測算法，根據(jù)背景噪聲進(jìn)行自適應(yīng)門限計(jì)算，并對(duì)功率譜中HS進(jìn)行檢測，若HS不存在，則將最大線譜譜峰與門限進(jìn)行比較，判斷目標(biāo)是否出現(xiàn)，以實(shí)現(xiàn)空中移動(dòng)聲源的水下跨界檢測。應(yīng)用該算法對(duì)水下單水聽器接收的螺旋槳飛機(jī)飛行噪聲實(shí)測信號(hào)進(jìn)行了目標(biāo)檢測仿真，并對(duì)該算法的檢測性能進(jìn)行了分析。

1 飛行目標(biāo)噪聲及入水傳播特性分析

帶有螺旋槳的飛機(jī)在飛行時(shí)會(huì)有較大的噪聲，尤其是直升機(jī)飛行時(shí)，其輻射噪聲主要由主槳、尾槳和發(fā)動(dòng)機(jī)產(chǎn)生。噪聲可分為氣動(dòng)噪聲和機(jī)械噪聲兩大類[8-9]，機(jī)械噪聲屬于高頻噪聲，傳播到空氣中很快被吸收，其噪聲源只對(duì)飛機(jī)近場噪聲起作用；氣動(dòng)噪聲主要集中在頻譜的低頻部分，對(duì)遠(yuǎn)場噪聲起決定性作用，功率較強(qiáng)且持續(xù)不斷。由于直升機(jī)旋翼在結(jié)構(gòu)上有周期轉(zhuǎn)動(dòng)機(jī)制，其遠(yuǎn)場噪聲功率譜是由具有諧波關(guān)系的離散線譜疊加于寬帶連續(xù)譜上構(gòu)成的，呈諧波關(guān)系的頻率分量構(gòu)成一個(gè)HS。例如，直升機(jī)由于有主槳和尾槳，通常由兩個(gè)不同基頻的HS組成，某直升機(jī)的時(shí)域波形和功率譜如圖1所示。圖1(b)中的慢變曲線是分離出的趨勢項(xiàng)波形，其中*和°表示頻譜中離散線譜構(gòu)成的兩個(gè)HS，分別對(duì)應(yīng)主槳槳葉通過頻率(BPF)及其各次諧波頻率，以及尾槳BPF及其各次諧波頻率。

螺旋槳BPF用fBPF表示，它對(duì)于某一種螺旋槳飛機(jī)而言是一個(gè)固定常數(shù)，可由(1)式計(jì)算[10]：

fBPF=iwr/60，

(1)

式中：i為諧波次數(shù)；w為槳葉個(gè)數(shù)；r為轉(zhuǎn)速。當(dāng)飛機(jī)改變速度和狀態(tài)時(shí)，螺旋槳旋轉(zhuǎn)速度不變，fBPF的值也不變。

一般檢測問題中假設(shè)噪聲為平穩(wěn)高斯隨機(jī)過程，然而對(duì)于水下探測，由于目標(biāo)和傳感器周圍的環(huán)境條件(尤其是風(fēng)速和海浪)影響，有可能改變?cè)肼暤慕y(tǒng)計(jì)特性，成為非平穩(wěn)非高斯過程。為便于處理，基于微分法思想，假設(shè)噪聲在一個(gè)非常短的時(shí)間內(nèi)保持平穩(wěn)，可將其視為平穩(wěn)信號(hào)，故信號(hào)處理時(shí)將數(shù)據(jù)分為若干個(gè)短數(shù)據(jù)段進(jìn)行。

空氣中點(diǎn)聲源發(fā)出的聲波傳播入水過程會(huì)有較大衰減。按照波動(dòng)理論[11]，在滿足空氣中點(diǎn)聲源高度要求和聲傳播路徑的相位變化條件基礎(chǔ)上，空氣中高度為h的點(diǎn)聲源可以等效為水中深度為d的點(diǎn)聲源，同時(shí)其聲壓級(jí)衰減20lg(kαd)，其中kα為聲波在空氣中的波數(shù)。圖2所示為空氣中點(diǎn)聲源的水下等效源[11]示意圖，直觀地表明了空氣中點(diǎn)聲源可等效為水下聲源加上空氣聲源到水面的傳播衰減。

因此可將聲源空氣- 水跨界傳播問題轉(zhuǎn)換為水下聲傳播問題，進(jìn)而分析水下聲場。根據(jù)對(duì)波動(dòng)方程解法的不同，本文選取基于拋物方程的RAMGEO模型進(jìn)行仿真計(jì)算，該模型改進(jìn)了拋物方程模型，并與穩(wěn)定的自身初始場結(jié)合[12]，適用于淺海低頻噪聲環(huán)境。本文具體仿真參數(shù)如下：聲源深度10 m，聲源頻率200 Hz，接收水聽器在水下呈垂直線陣排布，海水深度250 m. 利用RAMGEO模型仿真計(jì)算得到的聲源水下傳播衰減趨勢如圖3所示。

由圖3(a)可以看出，聲源在水下環(huán)境水平距離200 m以內(nèi)的范圍大致呈線性衰減，水平距離每增加約10 m，聲壓級(jí)降低2 dB，在水平距離超出200 m范圍后，聲壓級(jí)衰減明顯減小。這表明只要聲源聲壓級(jí)足夠大，聲波就可以實(shí)現(xiàn)在海水介質(zhì)中的遠(yuǎn)距離傳播。聲源在深度范圍上的衰減如圖3(b)所示。由圖3(b)可以看出，在淺海環(huán)境下水深100 m之后，聲波傳播衰減趨勢發(fā)生了較大變化，這可能是因?yàn)楹５茁暡ǚ瓷湟约暗讓幽嗌吵练e物引起的聲能損耗造成的。

由于海上直升機(jī)目標(biāo)往往帶有較大的噪聲聲壓級(jí)，按照等效源理論完全可以滿足聲波的遠(yuǎn)距離傳播，從而為海上戰(zhàn)場聲目標(biāo)的探測提供了有力手段。

2 基于自適應(yīng)門限的HS檢測算法

2.1 算法原理

目標(biāo)檢測基于以下兩種假設(shè)：1)H0-沒有目標(biāo)；2)H1-目標(biāo)出現(xiàn)。檢測算法原理框圖如圖4所示。

由圖4可知，本文提出的檢測算法是一種在頻域中采用HS和最強(qiáng)譜線檢測相結(jié)合的檢測方法，首先需要計(jì)算采集信號(hào)功率譜和背景噪聲功率譜擬合曲線，并根據(jù)擬合曲線計(jì)算檢測門限(DT)，在功率譜中搜尋線譜，將線譜按譜峰高低進(jìn)行降序存儲(chǔ)，若在搜尋的若干個(gè)線譜中有2個(gè)以上存在諧波關(guān)系，則認(rèn)為有HS并判斷為目標(biāo)出現(xiàn)，輸出狀態(tài)為H1；若沒有HS，則將譜線中最強(qiáng)的譜線峰值與DT進(jìn)行比較，該門限由背景噪聲功率譜擬合曲線確定：若大于DT，則判斷為目標(biāo)出現(xiàn)、狀態(tài)為H1，否則判斷為沒有目標(biāo)、狀態(tài)為H0.

2.2 HS檢測

HS檢測時(shí)，必須首先對(duì)水聽器接收信號(hào)時(shí)間序列采用Welch法[13-14]計(jì)算功率譜：

(2)

式中：PSD(f)為接收信號(hào)的功率譜；M為數(shù)據(jù)窗長度；L為數(shù)據(jù)段總數(shù)；xb(n)為第b個(gè)數(shù)據(jù)段中的第n個(gè)數(shù)據(jù)，b=1,2,3，…；d(n)為第n個(gè)數(shù)據(jù)的窗值。

在HS檢測過程中，線譜的判斷是算法中至關(guān)重要的一步，如第1節(jié)所述，螺旋槳飛機(jī)聲信號(hào)的特征頻率(基頻)及其諧波集中在中低頻段，因此線譜檢測過程中首先需要選擇相應(yīng)的頻率范圍，然后進(jìn)行線譜搜索和判定。在選定頻段內(nèi)以滑動(dòng)窗搜索所有譜峰，按譜峰值大小降序排列，選用前幾個(gè)線譜譜峰做HS檢測。線譜的搜索及判定方法如下：

1) 線譜值必須是功率譜中的局部極大值Lmax，該極大值對(duì)應(yīng)的頻點(diǎn)為fk(k=1,2,…,K，K為線譜數(shù)量)，局部極大值判斷過程是：將當(dāng)前檢測點(diǎn)k的功率譜值分別與其左、右相鄰點(diǎn)k-1和k+1的功率譜值求1階差分Δl和Δr，若滿足Δl<0且Δr>0，則Lmax暫定為此鄰域內(nèi)的局部極大值；

2) 若步驟1中Lmax的大小在k點(diǎn)的左、右峰寬Lsize鄰域范圍內(nèi)均為最大值，則Lmax暫定為線譜峰值；

3) 類似于步驟2，若在Lmax點(diǎn)的左、右峰寬Lsize鄰域范圍內(nèi)，分別存在極小值Lminl和Lminr，且比值2×Lmax/(Lminl+Lminr)大于某個(gè)閾值，則可認(rèn)為譜峰Lmax是顯著的，可將fk頻點(diǎn)處的Lmax作為線譜；

4)Lmax譜線的峰高定義為

hs=|Lmax-(Lminl+Lminr)/2|，

(3)

從而線譜的選取兼顧了譜峰左右邊界、峰寬和峰高等特征判據(jù)，對(duì)線譜特征的辨識(shí)度較高。

HS檢測按照以下方程[15]判斷：

|fk/fm-imk|≤ε，

(4)

式中：fm為假定基頻，m為主槳和尾槳基頻個(gè)數(shù)，m∈(1,2)；imk為諧波次數(shù)，imk=1,2,3；ε為選定的閾值。

若檢測到若干線譜對(duì)應(yīng)的頻率與該假定基頻有明顯的諧波關(guān)系，則認(rèn)為檢測到一個(gè)HS，該假定基頻即為HS的基頻。HS構(gòu)建策略如下：

1) 根據(jù)線譜搜索及判定方法進(jìn)行線譜檢測；

2) 根據(jù)線譜峰高h(yuǎn)s的大小，選出若干峰高較大線譜的頻點(diǎn)fk；

3) 將K個(gè)線譜的頻點(diǎn)依次作為假定基頻fm，例如將第m個(gè)頻點(diǎn)作為基頻，其他K-1個(gè)頻點(diǎn)fk用(4)式進(jìn)行諧頻檢測，即只要(4)式中imk取值為2或3時(shí)fk滿足(4)式，即可認(rèn)為fk與fm具有諧頻關(guān)系，它們對(duì)應(yīng)的線譜具有諧波關(guān)系。

在檢測中對(duì)于某個(gè)基頻，有時(shí)會(huì)出現(xiàn)不止一個(gè)HS，這時(shí)需要用權(quán)函數(shù)決定取舍，對(duì)每個(gè)HS內(nèi)所有頻率的能量進(jìn)行計(jì)算并作為權(quán)值，選擇具有最大能量的HS，與其對(duì)應(yīng)的基頻作為提取的螺旋槳飛機(jī)特征頻率。

在HS檢測中，由于諧波信號(hào)特性不同，例如頻率檢測范圍、線譜個(gè)數(shù)、線譜譜峰寬度、HS內(nèi)諧波數(shù)目以及閾值等參數(shù)，需要對(duì)實(shí)際信號(hào)進(jìn)行分析后確定這些參數(shù)的最佳值，從而使諧波檢測效果最佳。

2.3 自適應(yīng)DT的計(jì)算

由2.1節(jié)的算法原理可知，在進(jìn)行目標(biāo)檢測時(shí)首先進(jìn)行HS檢測，若沒有HS則尋找線譜中最強(qiáng)的譜峰，并與DT比較，若目標(biāo)出現(xiàn)，則由目標(biāo)發(fā)出的噪聲級(jí)再加上背景噪聲的噪聲級(jí)就會(huì)大于該DT. DT計(jì)算公式為

TDT=α+P(f)，

(5)

式中：TDT為實(shí)際計(jì)算檢測門限D(zhuǎn)T；α為門限調(diào)節(jié)因子；P(f)為無目標(biāo)時(shí)背景噪聲信號(hào)功率譜的多項(xiàng)式擬合曲線。在實(shí)際環(huán)境中，背景噪聲的大小是變化的，為了使檢測門限能自適應(yīng)地隨之調(diào)整，本文提出采用α來動(dòng)態(tài)調(diào)節(jié)門限，它可以根據(jù)無目標(biāo)時(shí)不同環(huán)境下測得的P(f)來計(jì)算，使DT成為隨環(huán)境變化的自適應(yīng)門限。

一般情況下，實(shí)際背景噪聲發(fā)生變化的主要原因是環(huán)境變化使得某一部分頻率分量的噪聲能量增加或減少，從而造成背景噪聲功率譜曲線[16]上某些頻率分量的功率譜幅度增大或減小，因此這種變化是一種非線性變化，在缺乏先驗(yàn)知識(shí)的條件下，對(duì)這種變化進(jìn)行定量描述是很困難的。因此本文提出采用最小二乘法對(duì)背景噪聲功率譜進(jìn)行擬合，從而對(duì)背景噪聲級(jí)進(jìn)行近似定量描述，α和DT的計(jì)算過程如下：

1) 計(jì)算第l幀背景噪聲nl(t)的功率譜，并取分貝數(shù)后為

Pn(f,l)=10lg(PSDn(f,l))，

(6)

式中：Pn(f,l)為第l幀背景噪聲的功率譜分貝數(shù)；PSDn(f,l)為第l幀背景噪聲的功率譜幅值。

2) 對(duì)Pn(f,l)采用最小二乘法進(jìn)行曲線擬合，設(shè)擬合多項(xiàng)式的階數(shù)為k，第l幀的擬合公式為

P(f,l)=p1fk+p2fk-1+…+pkf+pk+1，

(7)

式中：P(f,l)為第l幀背景噪聲的功率譜擬合曲線；pk為擬合系數(shù)，擬合曲線如圖5中背景噪聲中的連續(xù)實(shí)線。

3) 門限調(diào)節(jié)因子α的計(jì)算示意如圖5所示，圖5中背景噪聲功率譜最大值對(duì)應(yīng)的頻點(diǎn)為fα，則α取值為背景噪聲功率譜的最大值(圖5中黑色實(shí)心圓點(diǎn))與其擬合曲線P(f,l)在fα頻點(diǎn)取值的差值，第l幀門限調(diào)節(jié)因子α的計(jì)算公式為

α(l)=max[Pn(f,l)]-P(fα,l) ，

(8)

式中：max[Pn(f,l)]為第l幀背景噪聲功率譜的最大值；fα為第l幀背景噪聲功率譜最大值對(duì)應(yīng)的頻率值；P(fα,l)為第l幀背景噪聲功率譜頻率值為fα?xí)r對(duì)應(yīng)的擬合曲線值。

4) 根據(jù)第l幀背景噪聲計(jì)算的判決門限D(zhuǎn)T為

TDT(f,l)=α(l)+P(f,l).

(9)

從(9)式可知，這樣計(jì)算出的門限值是一個(gè)隨頻率變化的動(dòng)態(tài)曲線，且隨著時(shí)間的變化而變化。

3 實(shí)測信號(hào)檢測實(shí)驗(yàn)及結(jié)果分析

3.1 實(shí)測信號(hào)檢測實(shí)驗(yàn)

采用實(shí)測數(shù)據(jù)對(duì)本文算法進(jìn)行驗(yàn)證。實(shí)測數(shù)據(jù)的時(shí)域波形和時(shí)頻譜圖如圖6所示。圖6中3組數(shù)據(jù)文件Data1～Data3的記錄時(shí)間長度分別為622 s、974 s、651 s，3組數(shù)據(jù)為水下5 m單水聽器接收的某中型螺旋槳飛機(jī)水面上空飛行時(shí)的飛行噪聲信號(hào)，飛機(jī)勻速直線往返飛行，飛行高度約為150 m，數(shù)據(jù)記錄采樣頻率為10 kHz.

從時(shí)頻譜圖中可以同時(shí)觀察信號(hào)在時(shí)域和頻域上的信息，包括信號(hào)的頻譜特征以及頻率隨時(shí)間的變化。以實(shí)測信號(hào)Data1為例，分析圖6所示時(shí)頻譜圖中在200 s左右目標(biāo)出現(xiàn)的時(shí)刻，縱軸所示信號(hào)的能量急劇增強(qiáng)并在200 Hz上下伴有明顯的線譜特征，線譜之間存在諧頻特征，因此對(duì)該實(shí)測信號(hào)進(jìn)行諧波檢測是完全可行的。

對(duì)實(shí)測信號(hào)做預(yù)處理后進(jìn)行諧波檢測，每次取1 024個(gè)點(diǎn)作為1個(gè)數(shù)據(jù)段，將其分為互相重疊50%的子段，加窗作快速傅里葉變換(FFT)，所得功率譜的頻率分辨率約為1 Hz. 選取子段中前10個(gè)最大的譜峰進(jìn)一步作HS檢測，若出現(xiàn)3個(gè)以上諧波即1個(gè)HS時(shí)，則認(rèn)為成功檢測到目標(biāo)信號(hào)HS. 這里選取實(shí)測信號(hào)Data1的子段信號(hào)，檢測到的目標(biāo)信號(hào)諧波關(guān)系如圖7所示，具有諧波關(guān)系的線譜在圖中用圓點(diǎn)示出。對(duì)3次實(shí)測信號(hào)Data1～Data3結(jié)合自適應(yīng)門限調(diào)節(jié)因子作諧波檢測，發(fā)現(xiàn)目標(biāo)時(shí)檢測到的HS如表1所示。

Tab.1 Harmonic set of measured signals Hz

3.2 實(shí)測信號(hào)檢測結(jié)果分析

以第1次實(shí)測信號(hào)Data1為例進(jìn)行分析，如表1所示，HS檢測的線譜頻點(diǎn)有：192.9、249.0、273.4、397.9、500.5、556.6、1 018.1，其中存在諧波關(guān)系的HS有{192.9，397.9}，{249.0，500.5，1 018.1}，{273.4，556.6}。在實(shí)測過程中存在船舶噪聲干擾，如圖8所示為采集的一段無飛機(jī)的純船舶噪聲段，由其頻譜分析可知此船噪聲主要集中在245 Hz左右，對(duì)其單獨(dú)進(jìn)行HS檢測，檢測到的HS{244,488,967}與Data1中的HS{249.0,500.5,1 018.1}基本一致，從而目標(biāo)檢測時(shí)可用此HS排除船舶噪聲的干擾，并以此檢測識(shí)別飛機(jī)目標(biāo)。

由于實(shí)測信號(hào)受風(fēng)浪、船只等外界因素的影響，對(duì)3次實(shí)測信號(hào)檢測到的諧波線譜不完全相同。但由于是同一架飛機(jī)，HS檢測所得諧波基頻相差不大，表明諧波檢測算法在不斷變化的環(huán)境條件下仍能保持穩(wěn)定的檢測性能。

根據(jù)圖4的檢測算法原理，若某幀信號(hào)檢測的線譜間不存在諧波關(guān)系，則尋找線譜中最強(qiáng)的譜峰，與DT進(jìn)行比較來判斷有無目標(biāo)，若大于DT則輸出H1，否則判斷為H0，繼續(xù)觀測。為了保證檢測效果，若連續(xù)檢測10次都判斷為H1，則認(rèn)為發(fā)現(xiàn)目標(biāo)，其中DT為發(fā)現(xiàn)目標(biāo)前最后一幀信號(hào)計(jì)算所得的門限值，檢測概率的計(jì)算公式為

(10)

使用該方法進(jìn)行檢測的結(jié)果如圖9所示，根據(jù)圖9的檢測結(jié)果可以獲得實(shí)測信號(hào)的相關(guān)預(yù)警信息，在檢測到目標(biāo)后即發(fā)出預(yù)警信息，同時(shí)若檢測到HS，則可利用HS的基頻對(duì)目標(biāo)進(jìn)行輔助識(shí)別，判斷目標(biāo)是否為飛機(jī)。按照實(shí)驗(yàn)中該型號(hào)飛機(jī)的正常巡航速度v=36 m/s計(jì)算，預(yù)警距離為

(11)

式中：t1為發(fā)現(xiàn)目標(biāo)的時(shí)刻；t2為目標(biāo)過頂?shù)臅r(shí)刻。

由(11)式可得預(yù)警距離如表2所示。從圖9和表2可知，從發(fā)現(xiàn)目標(biāo)開始直至目標(biāo)到達(dá)最靠近點(diǎn)(CPA)[17]，可以保持較遠(yuǎn)的預(yù)警距離，并在檢測到目標(biāo)后保持較高的檢測概率。

表2 實(shí)測信號(hào)預(yù)警信息

4 結(jié)論

本文針對(duì)螺旋槳飛機(jī)的水下聲探測，根據(jù)其輻射噪聲信號(hào)的諧波特性，提出了HS檢測與最強(qiáng)譜線檢測相結(jié)合的檢測算法；采用門限調(diào)節(jié)因子α自適應(yīng)地調(diào)整DT，用實(shí)測信號(hào)進(jìn)行了目標(biāo)檢測試驗(yàn)驗(yàn)證。主要得出以下結(jié)論：

1) 螺旋槳飛機(jī)噪聲信號(hào)具有諧波特性，通過基于拋物方程的RAMGEO模型進(jìn)行仿真計(jì)算表明，可將聲源空氣- 水跨界傳播問題轉(zhuǎn)換為水下聲傳播問題，從而按照等效源理論進(jìn)行空中目標(biāo)的水下探測。

2) 在HS檢測過程中，線譜檢測是至關(guān)重要的一個(gè)環(huán)節(jié)，可以通過選取兼顧譜峰左右邊界、峰寬和峰高等特征來提取線譜。

3) 若某幀信號(hào)檢測的線譜間不存在諧波關(guān)系，則可尋找線譜中最強(qiáng)的譜峰，從而利用基于自適應(yīng)門限調(diào)節(jié)因子α的最強(qiáng)譜線檢測方法增強(qiáng)HS檢測方法的穩(wěn)定性，有效提高目標(biāo)的檢測距離。