基于互模糊函數的高亞音速目標聲學跟蹤方法

陳昭男，付 軍，秦亮亮，李 鵬

(1.中國人民解放軍91550部隊,遼寧 大連 116023；2.中國人民解放軍91604部隊,山東 煙臺 265700)

0 引言

低空高速飛行目標的探測和跟蹤是雷達探測技術的盲點之一。 近年來，基于目標聲信號的探測和定位技術作為雷達探測技術的補充方案，得到了廣泛關注。聲探測技術通過被動接收目標發出的聲波以確定目標位置，具有隱蔽性強、不易受電磁干擾和可探測低空目標等優點[1]。

根據目標運動速度是否超音速，可將運動目標分為超聲速和亞聲速兩類。對于超聲速類目標，其聲信號呈現明顯的激波特性，主要通過檢測激波來實現對目標的跟蹤[1]；對于亞聲速類目標，目前研究主要集中在對直升機、無人機等低速目標的跟蹤[2-3]，對接近聲速的高亞音速類目標還未見相關研究。本文主要針對低空高亞音速目標的聲學跟蹤方法展開研究。在低空低速目標的聲學跟蹤方法方面，主要采用的方法有時延估計定位法[4-7]、基于瞬時頻率估計的目標定位方法[8-9]、空間譜估計定向方法[10-12]等。在目前已知的三大類跟蹤方法中時延估計法計算效率高，易于實現，且應用范圍廣，因此主要針對這一方法展開研究，提出了一種利用互模糊函數的高亞音速目標聲學跟蹤方法。

1 運動目標時差量的估計方法

1.1 聲基陣陣型設計

為了提高隱蔽性，采用被動跟蹤體制，即只通過接收運動目標發出的聲音來完成定位及跟蹤。本文針對窄帶聲源進行定位研究，對于寬帶聲源，如果其頻譜中含有明顯的離散譜，本方法也適用。在此采用無指向性的傳感器陣列來實現對目標的定位跟蹤，傳感器陣列由至少三個傳聲器組成，各個傳聲器分別位于不同地點。在采用平面五元十字陣的情況下，目標聲源與傳感器陣列的幾何關系如圖1所示。

圖1 目標與傳感器陣列的幾何關系圖Fig.1 Geometric diagram of targets and sensor array

通過利用不同傳聲器接收目標信號的觀測量，包括到達時間差(Time Difference of Arrive,TDOA)和到達時間的比例尺差(Scale Difference of Arrive,SDOA)對目標的運動參數進行估計。

假設傳感器陣列一共有N個傳聲器，在傳感器陣列接收到的一組信號中，最早到達的一路信號為f0(t)，其對應的傳聲器稱之為基準傳聲器，其他路信號為fi(t)，i=1,2,…,N-1，第i路信號的時長為Ti，數據采集模塊的采樣頻率為fs，對應的采樣周期Ts=1/fs。目標在k時刻的運動參數包括位置坐標[xk,yk,zk]，速度矢量[vx,k,vy,k,vz,k]以及加速度值ak。

由于目標是運動的，不同路信號間的TDOA和SDOA是隨時間變化的。對目標運動軌跡進行分段，每段長度為L(單位為m)，在這一段長度內，認為目標的運動參數是不變的。目標運動軌跡是對應到接收信號的，相應地，對各路接收信號也按照相同的時間間隔進行均勻分段。通過對每一段的目標運動參數分別進行估計，即可得到目標的運動軌跡及相關運動參數。

1.2 TDOA和SDOA值估計

假設對目標的速度值大小v有一個先驗的約束范圍vmin≤v≤vmax，將目標的運動軌跡均勻分為NL段，也就是將對應的各個接收信號按照同樣的時間間隔且時間同步地均勻分為NL段，假設信號采樣周期為Ts，每段對應的目標運動距離為L，相應的采樣點數N0為?L/vmaxTs」。在實際應用中，為保證定位跟蹤精度，每一段信號中的采樣點數目N0不能太少，采樣周期Ts和劃分段數NL的取值使N0=?L/vmaxTs」大于某個預先設定的最小值N0,min。根據后續信號處理算法的實際需求，N0,min的典型值取256或者512。

傳感器陣列接收的各路信號組成接收信號向量[f0(t),f1(t),f2(t),…,fN-1(t)]，其中f0(t)為最先到達的基準信號。將目標的運動軌跡均勻分為NL段，也就是將對應的各個接收信號時間同步后，按照相同時間間隔均分為NL段，對分段后的信號向量，在每一段內取N0個采樣點，對基準信號和其他各路信號的TDOA和SDOA分別進行估計。采用互模糊函數(Cross Ambiguity Function,CAF)方法進行TDOA和SDOA的聯合估計。

假設對于第k段接收信號，目標相對于基準傳聲器和傳聲器i的徑向速度分別為vr0,k和vri,k。徑向速度為目標相對于傳聲器的速度，相向取正，相離取負。用φik表示傳聲器i接收信號的第k段相對于基準傳聲器的SDOA，則SDOA與徑向速度的關系為[13]：

(1)

基準傳聲器的接收信號y0(t)可表達為：

y0(t)=x(t)+n1(t)，

(2)

式中，x(t)為運動目標發出的聲音傳播到基準傳聲器的信號，n1(t)為噪聲。根據式(2)，傳聲器i接收的信號可表示為：

(3)

式中，n2(t)為噪聲，a為傳聲器i相對于基準傳聲器的信號相對幅度增益，τ為傳聲器i相對于基準傳聲器所接收信號之間的時差TDOA，φ為傳聲器i相對于基準傳聲器的時間伸縮因子SDOA，其表達式如式(1)所示。在得到y0(t)和yi(t)兩路信號后，采用互模糊函數法來對τ和φ進行聯合估計。y0(t)和yi(t)的互模糊函數CAF表達式為：

(4)

式中，T為信號時間長度。使CAF取得最大值的τ和φ的組合，分別為TDOA和SDOA的最優估計值，選用其作為后續目標跟蹤過程所使用的TDOA和SDOA值。

當聯合估計TDOA和SDOA時，隨著目標的運動，目標相對于傳聲器的徑向速度是迅速變化的，導致SDOA也是快速變化的，因此將目標的運動軌跡劃分成小段，在每一段內分別進行TDOA和SDOA的估計。目標的運動軌跡是映射到各路接收信號上的，將各路接收信號劃為小段時，采用連續劃分的方式，即各個小段是相連的，同時每一小段內的所有信號都參與估計運算。假設最早到達的一路信號為f0(t)，其與傳聲器i接收信號fi(t)的TDOA和SDOA聯合估計表達式為：

(5)

式中，T1為信號段的長度。對于TDOA和SDOA的估計，具體做法是通過對TDOA和SDOA的二維搜索，尋找使兩個信號的互模糊函數取得最大值的參數組合，即為該時刻的TDOA和SDOA值。

2 目標運動參數估計

2.1 目標運動參數估計模型

在得到觀測參數TDOA和SDOA后，建立目標運動方程和傳感器觀測方程，對目標運動參數進行估計。假設系統的運動參數為位置坐標及速度矢量，首先根據目標前后時刻運動參數的相互關系，建立目標運動方程：

(6)

式中，a為加速度值，|vk|為k時刻速度模值，加速度值通過初始時刻估計值或者先驗值得到，在后續狀態更新過程中作為已知值代入目標運動方程中；nk為系統隨機輸入噪聲向量。式(6)給出的是目標勻加速運動時的運動方程。

根據觀測量與運動參數的相互關系，建立觀測量與運動參數的解析表達式，即為傳感器觀測方程。當利用目標聲音信號到不同聲傳感器的傳播路徑差進行定位時，其傳感器觀測方程為：

i=1,2,…,N-1,

(7)

φik=(c-vr0,k)/(c-vri,k)+χi,

i=1,2,…,N-1，

(8)

式中，[sxi,syi,szi]為聲傳感器i的位置坐標，[sx0,sy0,sz0]為基準聲傳感器的位置坐標，ωi和χi為系統觀測誤差。目標運動方程和傳感器觀測方程中的噪聲量nk，ωi，χi均服從高斯分布。對于速度信息的估計，通過對前后時刻目標位置差分的方法得到。只要建立觀測量與運動參數的解析表達式，就得到了傳感器觀測方程，因此傳感器觀測方程不限于上述兩種。

2.2 模型求解

對于初始時刻運動參數的估計，采用極大似然估計方法。利用傳感器觀測方程分別構建目標的位置坐標[x0,y0,z0]和速度矢量[vx,0,vy,0,vz,0]的似然函數p(τ/[x0,y0,z0])和p(φ/[vx,0,vy,0,vz,0])，再利用極大似然估計方法對目標位置坐標[x0,y0,z0]和目標速度矢量[vx,0,vy,0,vz,0]分別進行估計。對位置坐標進行極大似然估計為例，對其估計過程進行說明。根據式(7)，系統觀測誤差ωi服從高斯分布，對于某路信號的時差測量值τi0，單路信號得到的似然函數p(τ/[x,y,z])可表達為：

(9)

式中，σω為系統觀測誤差變量ωi的方差。采用N路聲傳感器對目標進行測量，在同一時刻，共得到N-1個測量值，總的似然函數是各路信號得到似然函數式(9)的乘積，因此目標位置坐標的似然函數表達式為：

(10)

在無誤差情況下，觀測變量τi0與目標位置坐標之間的關系式為：

(11)

在某一位置坐標[x0,y0,z0]附近的區域，將式(11)用泰勒級數展開的方法，用一次多項式近似，得到的近似表達式為：

(12)

將式(12)代入式(10)，可以得到：

(13)

式中，指數項部分F(x,y,z)的表達式為：

(14)

F(x,y,z)是一個二次多項式，可以寫成如下標準形式：

F(x,y,z)=a1x2+a2y2+a3z2+a4xy+

a5xz+a6yz+a7x+a8y+a9z+a10。

(15)

(16)

(17)

式中，i=0,1,…,N-1。根據式(17)，分別將vri,0和vri,k的表達式代入式(8)，即可得到聲傳感器觀測量SDOA與初始時刻目標速度矢量[vx,0,vy,0,vz,0]的數學關系式，依據此關系式即可得到目標初始速度矢量[vx,0,vy,0,vz,0]的似然函數P(φ/[vx,0,vy,0,vz,0])。針對目標速度矢量的似然函數，利用對目標位置坐標進行極大似然估計的相同步驟，采用式(9)～式(16)，即可完成對目標速度矢量[vx,0,vy,0,vz,0]的極大似然估計。在得到初始目標速度后，通過計算前后兩個時刻的速度變化率可得到目標加速度的估計值。

在通過極大似然估計方法或者利用先驗信息得到初始時刻目標運動參數后，后續時刻的目標運動參數根據目標運動方程和傳感器觀測方程，采用貝葉斯遞推狀態估計算法來實現。

3 算法仿真與誤差分析

3.1 仿真參數設置

為了驗證算法的可行性和準確性，在Matlab仿真環境中，對該方法的性能進行了仿真分析。對于高亞音速飛行目標的模擬，首先從50 Hz聲源中提取部分實測信號，再模擬產生其多普勒效應，最后在仿真環境中疊加高斯白噪聲。目標運動產生的多普勒頻移效應根據目標與傳感器的徑向速度來進行模擬，假設目標運動方向與目標聲波到傳感器傳播方向的夾角為α，則由運動引入的多普勒頻移為vf0cosα/c，v為目標運動速度，f0=50 Hz，c為聲速。在實際測得的50 Hz聲信號基礎上，對于多普勒效應的仿真模擬，采用對信號分段，并對每段信號持續時間(周期)根據多普勒頻移進行調整的方式完成，如原信號每段持續時間為12 ms，頻率為50 Hz，受到多普勒效應影響后的接收信號頻率為60 Hz，則該段信號波形不變，但持續時間調整為10 ms。

3.2 運動參數估計仿真

首先對目標信號特性進行分析。運動目標聲信號的時域波形和靜止聲源的功率譜如圖2和圖3所示。

圖2 亞音速目標聲信號波形Fig.2 Signal waveform of subsonic target

圖3 靜止聲源功率譜Fig.3 PSD of static sound source

圖4 不同信噪比下的目標位置誤差Fig.4 Target location error in different SNRs

圖5 不同信噪比下的目標角度誤差Fig.5 Target angle error in different SNRs

通過圖4和圖5可以看出，無論是目標位置還是角度，在這種情況下，其誤差大小均與信噪比關系不大，距離誤差維持在2 m左右，方位角誤差在5°以內，俯仰角誤差相對大一些。上述仿真結果表明該算法可以實現對低空運動目標的穩定跟蹤。由上述仿真結果，在上述參數設置下，信噪比不是影響定位精度的主要因素，傳聲器采樣率、傳感器陣列結構等因素都對定位精度構成了制約。由圖5可見，俯仰角誤差明顯大于方位角誤差，這與所采用的傳感器陣列結構有關。傳感器陣列都位于同一平面上，各路接收信號對于俯仰角所產生的TDOA不敏感，因此導致其定位誤差較大。

在上述仿真參數設置的基礎上，首先將傳聲器采樣率提高到256 kHz，在不同信噪比下，對位置和角度誤差進行了蒙特卡洛仿真。其對目標定位的位置誤差和角度誤差與接收信號信噪比的關系如圖6和圖7所示。

圖6 高采樣率下的目標位置誤差Fig.6 Target location error in large sample rates

圖7 高采樣率下的目標角度誤差Fig.7 Angle error in large sample rates

由圖6可知，此時目標位置誤差和角度誤差均隨信噪比的增加而迅速下降，目標位置誤差均維持在1.8 m以下，而目標角度誤差除了在0 dB處大于2°以外，其余均在2°以下。整體而言，目標位置誤差和角度誤差與傳聲器采樣率為50 kHz時相比，均出現了明顯下降。該仿真結果驗證了只有傳聲器采樣率達到一定量級時，目標定位誤差才隨著信噪比呈現下降趨勢。

圖8 不同信噪比下的目標位置誤差Fig.8 Target location error in different SNRs

圖9 不同信噪比下的目標角度誤差Fig.9 Target angle error in different SNRs

由該仿真結果可見，目標位置誤差與第一次仿真結果相差不大，而角度誤差，尤其是俯仰角誤差產生了較大程度的下降，俯仰角誤差均維持在2°以內。該仿真結果表明，傳感器陣列在垂直方向上的分布將對俯仰角定位精度產生較大影響，傳感器陣列在垂直方向上分布越均勻時，其俯仰角定位誤差越小。

綜合以上布陣情況，假設基準陣元布設于原點處，設置空間直角坐標系，在x軸正負半軸對稱布設兩個陣元，其相對于基準陣元的時延分別為τ1和τ3(τ1對應x軸正半軸陣元)，同理y軸正負半軸兩個陣元的相對時延分別為τ2和τ4(τ2對應y軸正半軸陣元)，則目標方位角φ與時延估計量的關系為：

(18)

目標俯仰角與時延估計量的關系為：

(19)

通過比較式(18)和式(19)可見，對于分布于水平面的聲陣列來說，其目標方位角的估計精度有時延估計量的精度決定，而與陣列大小及聲速值無關，而其俯仰角估計精度則同時由陣列大小、聲速值以及時延量估計精度共同決定。時延量估計誤差主要由環境噪聲、各測量通道之間的相位不一致性以及聲波在空間中傳播的隨機起伏等因素引起。環境噪聲和傳播隨機起伏可通過改進算法進行抑制，而測量通道的相位不一致性需要通過改進測量設備，提升各通道相位特性來實現。

4 結束語

本文對亞聲速低空運動目標的定位跟蹤問題進行了研究，通過對目標輻射聲信號的TDOA和FDOA值進行測量，提出了一種基于互模糊函數的目標跟蹤方法。該方法基于互模糊函數，首先利用極大似然估計方法或者先驗信息得到初始時刻目標運動參數后，后續時刻的目標運動參數根據目標運動方程和傳感器觀測方程，采用貝葉斯遞推狀態估計算法來實現。仿真結果表明，該方法可以在較小的位置誤差和角度誤差范圍內，實現對低空運動目標的穩定跟蹤。