基于相關(guān)峰插值的五元十字陣被動(dòng)聲定位算法

馬少春，劉慶華，黃靈鷺

（桂林電子科技大學(xué)信息與通信學(xué)院，廣西 桂林541004）

0 引言

在電子對(duì)抗領(lǐng)域內(nèi)，被動(dòng)聲定位技術(shù)是指利用目標(biāo)聲信號(hào)獲得目標(biāo)位置信息的技術(shù)，具有隱蔽性好、不受通視條件限制等優(yōu)點(diǎn)。目前使用傳聲器陣列對(duì)聲源進(jìn)行定位的方法主要分為三類：時(shí)延估計(jì)、譜估計(jì)和定向波束形成?；跁r(shí)延估計(jì)的方法是指利用聲源發(fā)出的信號(hào)到達(dá)各個(gè)傳聲器的時(shí)間延遲對(duì)聲源位置進(jìn)行估計(jì)，自廣義互相關(guān)算法［1］提出以來，由于其計(jì)算速度快，能達(dá)到實(shí)時(shí)的要求，定位精度較高，而得到廣泛應(yīng)用。

被動(dòng)聲定位在實(shí)際應(yīng)用時(shí)往往受到各種因素的限制，如時(shí)延估計(jì)精度、基陣尺寸、目標(biāo)距離等［2］，特別是對(duì)時(shí)延估計(jì)的精度要求很高，所以準(zhǔn)確地進(jìn)行時(shí)延估計(jì)是被動(dòng)聲定位的關(guān)鍵。常規(guī)的互相關(guān)時(shí)延估計(jì)算法［3］是基于快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）計(jì)算信號(hào)在頻域的互譜進(jìn)而計(jì)算時(shí)間延遲的，計(jì)算得到的時(shí)延間隔為采樣周期的整數(shù)倍，時(shí)延估計(jì)精度為Ts／2，自適應(yīng)時(shí)延估計(jì)算法［4］精度雖優(yōu)于傳統(tǒng)互相關(guān)算法，然其收斂速度較慢計(jì)算量較大不利于實(shí)時(shí)實(shí)現(xiàn)。針對(duì)互相關(guān)時(shí)延估計(jì)方法使用FFT 進(jìn)行快速計(jì)算會(huì)在頻域產(chǎn)生柵欄效應(yīng)影響時(shí)延估計(jì)精度的問題，提出基于相關(guān)峰插值（Fine Interpolation of Correlation Peak，F(xiàn)ICP）方法［5］的五元十字陣［6-7］被動(dòng)聲定位算法。

1 被動(dòng)聲定位原理

對(duì)于空間中的確定聲源，由于聲信號(hào)到達(dá)各點(diǎn)傳聲器時(shí)刻各不相同，因此可以利用時(shí)間延時(shí)計(jì)算出距離差值，進(jìn)而建立幾何關(guān)系求出目標(biāo)位置信息。確定目標(biāo)聲源在三維坐標(biāo)系中的三個(gè)位置參數(shù)需要通過四個(gè)傳聲器獲得三個(gè)獨(dú)立的時(shí)延值，而五元十字陣克服了四元十字陣定位誤差與目標(biāo)方位角有關(guān)的問題且十字陣具有分維特性，陣列冗余度較小，因此系統(tǒng)中陣型結(jié)構(gòu)采用五元十字陣。平面五陣元十字形陣列是由四元十字陣的中心點(diǎn)加一個(gè)聲傳感器構(gòu)成，如圖1所示。

圖1 五元陣示意圖Fig.1 Five-elementcross array schematic diagram

圖中五個(gè)陣元的坐標(biāo)分別為M0（0，0，0），M1（d，0，0），M2（0，d，0），M3（-d，0，0），M4（0，-d，0），目 標(biāo)聲源S的球坐標(biāo)表示為（r，φ，θ），其中r為目標(biāo)S到坐標(biāo)原點(diǎn)的距離，φ 為目標(biāo)s的方位角，即目標(biāo)點(diǎn)S 在五元十字陣平面的投影與x 軸正方向的夾角（x 軸正方向約定為從傳聲器M0到傳聲器M1的方向），θ為目標(biāo)點(diǎn)S 的俯仰角（即五元陣平面法線與傳聲器M0和目標(biāo)點(diǎn)S 連線的夾角），d 為陣列孔徑，C 為空氣中聲音的傳播速率。定義τi=t0-ti（ti為聲源到達(dá)傳聲器i 的時(shí)刻），則聲源到達(dá)各個(gè)傳感器M1，M2，M3，M4與到達(dá)傳感器M0的時(shí)延可分別記為τ1、τ2、τ3、τ4，由 聲 陣 和 目 標(biāo) 的 幾 何 位 置 有 如 下 方程組：可得目標(biāo)的定位方程為：

考慮實(shí)際中目標(biāo)位于遠(yuǎn)場(chǎng)，此時(shí)有

式（2）求出的方位角范圍是-π／2，π／

［］2 ，實(shí)際中可通過判斷φ 值正負(fù)、τ1與τ3的關(guān)系來確定目標(biāo)投影所在象限對(duì)方位角進(jìn)行修正。由觀測(cè)信號(hào)計(jì)算出到達(dá)傳聲器的時(shí)延差值，利用式（2）-式（4）即可求得聲源目標(biāo)在五元十字陣中的球坐標(biāo)參數(shù)。

2 基于相關(guān)峰插值的定位算法原理

假設(shè)兩個(gè)傳聲器接收的信號(hào)模型為：

其中，s（n）為聲源信號(hào)，τ為時(shí)延，ν1（n）、ν2（n）表示加性噪聲。

給 定 一 個(gè) 新 變 換 MCZT（Modified Chirp Z Transform）及逆變換IMCZT 計(jì)算信號(hào)的細(xì)化頻譜：

式中，x（n）序列長(zhǎng)度為N，頻譜X（k）的間隔由N1決定，Δf=fs／N1，fs為采樣頻率。逆變換IMCZT求得的x′（n）≠x（n），其計(jì)算相當(dāng)于對(duì)MCZT求共軛。MCZT 變換形式與FFT 變換相同，可用FFT 算法實(shí)現(xiàn)快速計(jì)算。利用MCZT 計(jì)算兩路信號(hào)細(xì)化的頻譜和互譜，在頻域?qū)プV補(bǔ)零可以提高相關(guān)函數(shù)波形的分辨率。

兩路信號(hào)MCZT 變換為：

式中k=0，1，…，N-1。

根據(jù)相關(guān)定理計(jì)算得兩信號(hào)的互譜

此處R1（k）只是譜的前一部分，通過譜的共軛對(duì)稱性，可以獲得兩信號(hào)N1點(diǎn)的互譜R2（k）：

為提高相關(guān)函數(shù)的分辨率，對(duì)互譜R2（k）補(bǔ)零，擴(kuò)充為N2（N2≥N1）點(diǎn)的序列R（k），最后對(duì)R（k）進(jìn)行逆變換得到兩信號(hào)的相關(guān)函數(shù)r（n）：

其中n=0，1，…，N2-1。

從式（10）可以看到，相關(guān)函數(shù)的分辨率由Δt0=1／fs提高到了Δt0=N1＊fs／N2，fs為信號(hào)采 樣頻率。N1、N2應(yīng)滿足N2≥N1≥2 N-1，且N1=2γ（γ為整數(shù)），增大N2／N1可以提高相關(guān)函數(shù)的分辨率。在實(shí)際聲探測(cè)系統(tǒng)中時(shí)延差值總是處于有限范圍內(nèi)，相關(guān)函數(shù)的主峰處在零點(diǎn)附近，因此，對(duì)于N2點(diǎn)的r（n）只需取其前N 點(diǎn)與后N 點(diǎn)，組合成相關(guān)函數(shù)，可以大大降低計(jì)算量。

3 算法仿真與外場(chǎng)實(shí)驗(yàn)測(cè)試

3.1 FICP算法仿真

實(shí)驗(yàn)中采用信號(hào)x=sin（2πft）exp（-80t），其中f=100Hz，信號(hào)長(zhǎng)度N 取8192點(diǎn)，采樣頻率64 kHz，兩路信號(hào)時(shí)延差為0.001 6s，混入高斯白噪聲，信噪比SNR=8dB，F(xiàn)ICP計(jì)算中N1=64 000，N2／N1=10。

圖2是兩路信號(hào)的時(shí)域波形，圖3采用蒙特卡羅方法分別對(duì)基于FICP計(jì)算和基于FFT計(jì)算的時(shí)延結(jié)果進(jìn)行比較。從圖中可以看出，當(dāng)真實(shí)時(shí)延值處于采樣周期整數(shù)倍之間時(shí)，基于FFT 計(jì)算的互相關(guān)方法估計(jì)的時(shí)延值總是為采樣周期的整數(shù)倍，F(xiàn)ICP方法提高了相關(guān)函數(shù)的分辨率，計(jì)算得到的時(shí)延值精度更高。

圖2 兩路信號(hào)波形Fig.2 Two-channel signal waveforms

圖3 100次計(jì)算結(jié)果比較Fig.3 Comparison of 100times calculation results

3.2 基于FICP方法的五元十字陣定位算法仿真

實(shí)驗(yàn)中采用信號(hào)x=sin（2πft）exp（-50t），其中f=100Hz，信號(hào)長(zhǎng)度N 取8192點(diǎn)，采樣頻率64 kHz，目標(biāo)球坐標(biāo)系坐標(biāo)為（50，180°，85°），混入高斯白噪聲，信噪比SNR=8dB，F(xiàn)ICP計(jì)算中N1=64 000，N2／N1=10。

傳聲器得到的五路信號(hào)時(shí)域波形如圖4所示，圖5、圖6分別為基于FICP 計(jì)算和基于FFT 計(jì)算的信號(hào)相關(guān)波形圖。仿真中，F(xiàn)ICP 計(jì)算結(jié)果為（48.917，179.908，84.934），誤差率為（-2.17%，-0.05%，-0.08%），F(xiàn)FT 計(jì)算結(jié)果為（47.171，179.848，85.135），誤差率為（-5.66%，-0.08%，0.16%），從計(jì)算結(jié)果不難看出，F(xiàn)ICP方法計(jì)算的定位精度優(yōu)于基于FFT 計(jì)算的互相關(guān)方法的精度。

圖4 五路信號(hào)波形圖Fig.4 Five-channel signal waveforms

圖5 五路信號(hào)FICP計(jì)算相關(guān)波形圖Fig.5 Five-channel signal cross-correlation waveforms based on FICP algorithm

圖6 五路信號(hào)FFT 計(jì)算相關(guān)波形圖Fig.6 Five-channel signal cross-correlation waveforms based on FFT algorithm

3.3 外場(chǎng)測(cè)試實(shí)驗(yàn)

如圖7所示，聲源定位系統(tǒng)主要由傳聲器陣列、信號(hào)調(diào)理器、AD采集模塊、信號(hào)處理模塊等組成。其中信號(hào)處理模塊包含DSP核心板、外部存儲(chǔ)器和輸出顯示部分。傳聲器陣列將采集的聲信號(hào)經(jīng)AD轉(zhuǎn)換送給DSP處理，定位結(jié)果由單片機(jī)控制液晶屏完成顯示。

圖7 被動(dòng)聲定位系統(tǒng)框圖Fig.7 Passive acoustic location system block diagram

為了驗(yàn)證五元十字陣定位效果，在空曠樓頂進(jìn)行靜止目標(biāo)的定位實(shí)驗(yàn)。實(shí)驗(yàn)條件如下：樓頂長(zhǎng)寬為50×10m2，播放音效的音響擺放在樓頂天臺(tái)上，天臺(tái)高5m，陣列高度可調(diào)（0.8～1.2m），實(shí)驗(yàn)測(cè)試距離取10～45 m，五個(gè)傳聲器布置支架的十字臂上，陣元間距d=1m，傳聲器型號(hào)為CHZ-213，主要性能指標(biāo)為：全指向性，頻率響應(yīng)20～20kHz，靈敏度56.36mV／Pa，信號(hào)調(diào)理器×10檔輸出，信號(hào)采樣頻率Fs=64kHz，采樣長(zhǎng)度N=16 384點(diǎn)。實(shí)驗(yàn)系統(tǒng)如圖8、圖9所示。

圖8 信號(hào)調(diào)理器和信號(hào)處理模塊Fig.8 Signal disposal and signal processing module

圖9 五元十字陣陣列Fig.9 Five-element cross array

將五個(gè)傳聲器固定好在支架上，連接好系統(tǒng)設(shè)備，確定音響與傳聲器的位置坐標(biāo)，播放音頻信號(hào)，記錄觀測(cè)數(shù)據(jù)。實(shí)際測(cè)量中，音源固定在距聲陣列一定方位作為點(diǎn)聲源，通過改變支架位置和轉(zhuǎn)動(dòng)支架方向來改變目標(biāo)聲源相對(duì)位置，完成對(duì)被測(cè)目標(biāo)的距離、方位角、俯仰角的實(shí)時(shí)測(cè)量。

在外場(chǎng)試驗(yàn)中，實(shí)際測(cè)量得到的方位角和俯仰角誤差較小，但距離的誤差相對(duì)較大。此外，由于外場(chǎng)環(huán)境存在著干擾，測(cè)量的結(jié)果會(huì)出現(xiàn)抖動(dòng)，方位角和俯仰角測(cè)量結(jié)果波動(dòng)較小，距離測(cè)量結(jié)果跳動(dòng)較大，可通過卡爾曼濾波等后置處理方法得到較為準(zhǔn)確的測(cè)量值。表1、表2和表3對(duì)比分析了基于互相關(guān)方法的被動(dòng)聲定位算法的結(jié)果和基于FICP被動(dòng)聲定位算法計(jì)算的目標(biāo)定位結(jié)果，不難看出，基于FICP方法的被動(dòng)聲定位算法性能優(yōu)于基于互相關(guān)方法的定位算法。

表1 距離測(cè)定值Tab.1 The results of distance measurement

表2 方位角測(cè)定值Tab.2 The results of azimuth measurement

表3 俯仰角測(cè)定值Tab.3 The results of pitching angle measurement

4 結(jié)論

本文提出了基于FICP 方法的被動(dòng)聲定位算法，對(duì)算法的原理進(jìn)行了推導(dǎo)說明。FICP通過計(jì)算信號(hào)的細(xì)化頻譜改善了FFT 變換帶來的柵欄效應(yīng)，并且計(jì)算量小可以在嵌入式系統(tǒng)中實(shí)時(shí)實(shí)現(xiàn)。仿真和實(shí)驗(yàn)結(jié)果都表明，相對(duì)于基于互相關(guān)方法的被動(dòng)聲定位算法，采用FICP 方法的被動(dòng)聲定位算法能夠得到更高的時(shí)延估計(jì)精度，從而對(duì)目標(biāo)進(jìn)行更為準(zhǔn)確的定位。外場(chǎng)實(shí)驗(yàn)中也可以看到，采用FICP方法的被動(dòng)聲定位算法對(duì)方位角和俯仰角的定位精度較高，對(duì)距離的測(cè)量精度以及抗干擾能力有待進(jìn)一步改善。

［1］Knapp C H，Carter G C.The generalized correlation method for estimation of time delay［J］.IEEE Transactions on Acoustics，Speech，Signal Processing，1976，24（4）：320-327.

［2］陳華偉，趙俊渭.五元十字陣被動(dòng)聲定位算法及其性能研究［J］.探測(cè)與控制學(xué)報(bào)，2003，25（4）：11-16.

［3］羅柏文，于宏毅.采用合成方法的多路信號(hào)自適應(yīng)時(shí)延聯(lián)合估計(jì)［J］.信號(hào)處理，2013，29（2）：159-164.

［4］黃曉燕，羅建.自適應(yīng)時(shí)延估計(jì)算法在被動(dòng)聲定位系統(tǒng)中應(yīng)用［J］.探測(cè)與控制學(xué)報(bào)，2007，29（2）：73-76.

［5］楊亦春，馬馳州.相關(guān)峰細(xì)化的精確時(shí)延估計(jì)快速算法研究［J］.聲學(xué)學(xué)報(bào)，2003，28（2）：159-166.

［6］向瑾.基于四元十字陣的被動(dòng)聲定位技術(shù)研究［D］.太原：中北大學(xué)，2008.

［7］黃海軍.基于傳聲器陣列的聲源定位系統(tǒng)的初步研究［D］.南京：南京郵電大學(xué)，2013.