一種空間匹配的高速飛行目標(biāo)軌跡估計(jì)方法

鄧 步，李弘毅，顧亞平

（1.中國(guó)科學(xué)院聲學(xué)研究所東海研究站，上海 201815；2.中國(guó)科學(xué)院大學(xué)，北京 100049）

0 引言

相對(duì)于光學(xué)和雷達(dá)目標(biāo)探測(cè)，低空目標(biāo)聲探測(cè)采用被動(dòng)探測(cè)方法，受環(huán)境影響較小，因此受到廣泛關(guān)注[1-3]。當(dāng)目標(biāo)離測(cè)點(diǎn)陣列的距離較近時(shí)，可將傳感器陣列所獲得的信號(hào)視為滿足球面波，此時(shí)可將測(cè)得的信號(hào)到各個(gè)陣元的時(shí)間差作為先驗(yàn)信息 進(jìn)行參數(shù)估計(jì)，得到運(yùn)動(dòng)參數(shù)[4]或定位參數(shù)[5-6]。當(dāng)目標(biāo)離測(cè)點(diǎn)陣列的距離較遠(yuǎn)時(shí)，聲信號(hào)視為平面波，進(jìn)而測(cè)得信源方位[7]，再利用陣列網(wǎng)絡(luò)進(jìn)行定位[8]。

文獻(xiàn)[9]提出了一種基于聲多徑傳播模型的目標(biāo)運(yùn)動(dòng)參數(shù)估計(jì)方法。該方法利用不同路徑聲信號(hào)的到達(dá)時(shí)間估計(jì)值，實(shí)現(xiàn)對(duì)目標(biāo)運(yùn)動(dòng)參數(shù)的估計(jì)，該方法要求布設(shè)盡可能多的長(zhǎng)基線聲傳感器陣。

文獻(xiàn)[10]對(duì)瞬時(shí)頻率方法進(jìn)行了拓展，提出了一種基于單聲傳感器陣列的目標(biāo)運(yùn)動(dòng)軌跡估計(jì)方法。該方法首先利用多普勒效應(yīng)估計(jì)目標(biāo)飛行參數(shù)，結(jié)合波達(dá)方向（Direction of Arrival,DOA）估計(jì)目標(biāo)軌跡方向、水平偏置和高度，依據(jù)幾何模型確定其運(yùn)動(dòng)軌跡。該方法主要適用低速飛行目標(biāo)。

文獻(xiàn)[11]提出了基于五元十字陣的低空目標(biāo)聲測(cè)無(wú)源定位算法。在時(shí)延估計(jì)方面，該文獻(xiàn)提出了基于維納加權(quán)的頻域自適應(yīng)時(shí)延估計(jì)算法和基于雙譜的時(shí)延估計(jì)算法。在小尺寸高精度定向方面，該文獻(xiàn)首次將壓差式矢量聲傳感器引入低空目標(biāo)側(cè)向研究中，提出基于矢量傳感器陣的寬帶相干信號(hào)子空間最優(yōu)波束形成算法。

文獻(xiàn)[12]利用布設(shè)于有限空間內(nèi)的聲傳感器陣實(shí)現(xiàn)低空高速目標(biāo)跟蹤為背景，通過(guò)對(duì)某典型低空高速目標(biāo)噪聲波形的分析，發(fā)現(xiàn)其噪聲呈現(xiàn)寬帶低頻譜特征。因此，該文獻(xiàn)主要研究了基于時(shí)延估計(jì)的目標(biāo)定位跟蹤方法，提出了多面交匯的目標(biāo)軌跡估計(jì)方法。

目前對(duì)低空飛行目標(biāo)的研究主要集中在直升機(jī)、無(wú)人機(jī)等低速目標(biāo)，而低空高速飛行目標(biāo)威脅性更高，且更加難以跟蹤。高速飛行目標(biāo)的噪聲是一類瞬態(tài)信號(hào)[10]，難以準(zhǔn)確獲取其信號(hào)特征，給定位跟蹤帶來(lái)了較大困難。

影響高速目標(biāo)軌跡測(cè)量精度的原因很多，但主要因素有：（1）目標(biāo)高速運(yùn)動(dòng)，會(huì)產(chǎn)生多普勒頻移；（2）環(huán)境中的風(fēng)速是時(shí)變的，且風(fēng)向具有不確定性，造成聲音傳播各向異性。兩者均會(huì)引起陣列中陣元的測(cè)量誤差。

因此，在對(duì)高速飛行目標(biāo)進(jìn)行波達(dá)方位估計(jì)時(shí)，不可避免地存在估計(jì)誤差。同時(shí)，當(dāng)目標(biāo)高速運(yùn)動(dòng)的速度接近聲速時(shí)，在某一時(shí)刻上，不同陣列所接收到的信號(hào)并非信源在同一時(shí)刻所發(fā)出，難以在時(shí)間維度上進(jìn)行對(duì)齊。針對(duì)上述問(wèn)題，借鑒點(diǎn)云數(shù)據(jù)處理的思想[13]，本文在探測(cè)范圍內(nèi)建立字典，提出了一種基于空間匹配的軌跡估計(jì)方法。仿真結(jié)果表明，在一定誤差范圍內(nèi)，軌跡估計(jì)結(jié)果依然表現(xiàn)良好，達(dá)到預(yù)期目的。

1 字典建立

假設(shè)目標(biāo)的探測(cè)范圍已知，即目標(biāo)出現(xiàn)在100 m×100 m×200 m 的空間范圍內(nèi)，在探測(cè)區(qū)域頂點(diǎn)處安裝相同的測(cè)向陣列，如圖1 所示。目標(biāo)從任意方向飛過(guò)探測(cè)區(qū)域時(shí)，至少?gòu)膬蓚€(gè)陣列之間穿過(guò)。將陣列1 位置視為坐標(biāo)原點(diǎn)。每個(gè)陣列從字典建立到平面擬合流程一致，本文以陣列1 進(jìn)行推導(dǎo)說(shuō)明。

圖1 目標(biāo)探測(cè)區(qū)域示意圖 Fig.1 Schematic diagram of target detection area

將空間內(nèi)水平和垂直距離以設(shè)定間距d1為間隔進(jìn)行劃分，即若d1=1 m，該空間被劃分為2×106個(gè)具有獨(dú)立坐標(biāo)的確定的空間點(diǎn)，表示為：（xi,yi,zi）,i=1,2,…,N?？臻g點(diǎn)示意圖如圖2 所示。

圖2 空間點(diǎn)示意圖 Fig.2 Space point diagram

每個(gè)點(diǎn)相對(duì)于位于原點(diǎn)的陣列1 的方位角與俯仰角可通過(guò)坐標(biāo)轉(zhuǎn)換公式進(jìn)行轉(zhuǎn)換，投影到極坐標(biāo)下。將N個(gè)空間點(diǎn)的方位角θ與俯仰角?構(gòu)成字典進(jìn)行數(shù)據(jù)存儲(chǔ)。即：

2 軌跡估計(jì)

將接收到的信號(hào)進(jìn)行分幀處理，計(jì)算每一幀數(shù)的波達(dá)方位，設(shè)有K幀數(shù)據(jù)，結(jié)果表示為

計(jì)算每幀數(shù)據(jù)得到的波達(dá)方向與字典中每個(gè)空間點(diǎn)方位角與俯仰角的均方誤差，表示為

單幀數(shù)據(jù)計(jì)算得到的波達(dá)方位與空間點(diǎn)字典產(chǎn)生N個(gè)誤差結(jié)果。再進(jìn)行最小誤差匹配，即將計(jì)算得到的誤差進(jìn)行排序，搜索誤差最小的n個(gè)點(diǎn)進(jìn)行存儲(chǔ)，表示為

當(dāng)n=100 時(shí)，最小誤差分布圖如圖3 所示。從圖3 中可以看到，匹配得到的n個(gè)誤差最小空間點(diǎn)在空間中沿著一定方向進(jìn)行分布且分布密度與空間點(diǎn)和陣列1 的幾何距離成正比。

圖3 最小誤差分布圖 Fig.3 Schematic diagram of the minimum error distribution

2.1 三維空間直線擬合

本文使用整體最小二乘法[14-15]，對(duì)單幀數(shù)據(jù)得到的n個(gè)誤差最小空間點(diǎn)進(jìn)行擬合，以獲得空間直線的幾何參數(shù)，為后續(xù)平面擬合做準(zhǔn)備。

將空間直線的點(diǎn)向式方程變形為

式中：x0、y0、z0為空間直線的已知點(diǎn)；A、B、C為待求參數(shù)。

設(shè)V為擬合誤差，I為單位矩陣，對(duì)存儲(chǔ)的n個(gè)誤差最小空間點(diǎn)進(jìn)行構(gòu)造誤差方程。根據(jù)最小二乘理論中的間接平差模型[16]，令：

則誤差方程表示為

整體最小二乘法是同時(shí)顧及觀測(cè)值和系數(shù)矩陣誤差的一種平差方法，相比于最小二乘法，其參數(shù)估值是統(tǒng)計(jì)最優(yōu)解。本文采用易于編程實(shí)現(xiàn)的整體最小二乘迭代法，解算式（8）中的空間直線的誤差方程，平差準(zhǔn)則為

代入式（8），對(duì)系數(shù)矩陣B和參數(shù)向量X中各元素進(jìn)行求導(dǎo)，獲得迭代方程式：

其中：

設(shè)初始值為X0，使用整體最小二乘迭代法對(duì)誤差最小的n個(gè)點(diǎn)進(jìn)行擬合，結(jié)果如圖4 所示。

圖4 整體最小二乘直線擬合結(jié)果 Fig.4 The result of overall least squares linear fitting

該直線的方向向量可視為目標(biāo)與陣列坐標(biāo)的矢量方向，即可視為真實(shí)目標(biāo)位置估計(jì)點(diǎn)位于該擬合直線上。

同理對(duì)K幀數(shù)據(jù)連續(xù)進(jìn)行上述處理，擬合結(jié)果如圖5 所示。

圖5 對(duì)K 幀數(shù)據(jù)進(jìn)行的整體最小二乘直線擬合結(jié)果 Fig.5 The result of overall least square linear fitting for the K-frames of data

2.2 穩(wěn)健最小二乘平面擬合

由于高速飛行目標(biāo)在探測(cè)范圍內(nèi)為直線運(yùn)動(dòng)，在一個(gè)陣列的條件下，多幀最小二乘擬合直線理想情況下共面，即均在同一個(gè)平面上。理想條件下，兩條共面直線即可確定平面，但在2.1 節(jié)中對(duì)直線進(jìn)行擬合時(shí)存在一定誤差，因此直線共面條件不可能實(shí)現(xiàn)。此時(shí)使用空間直線上的采樣點(diǎn)來(lái)對(duì)平面進(jìn)行擬合，得到由方向直線構(gòu)成的近似平面。本文在每幀數(shù)據(jù)所確定的直線上均勻采樣100 個(gè)點(diǎn)構(gòu)成空間點(diǎn)數(shù)據(jù)集M，該數(shù)據(jù)集中不存在異常點(diǎn)，受文獻(xiàn)[13,17-18]啟發(fā)，采用一種穩(wěn)健最小二乘平面擬合方法。理論上，兩幀數(shù)據(jù)所構(gòu)成的空間點(diǎn)數(shù)據(jù)集即可確定相應(yīng)平面。

使用多項(xiàng)式函數(shù)表示平面：

式中：u、v、w為待求參數(shù)。

設(shè)數(shù)據(jù)集：

其最佳擬合平面應(yīng)滿足：

由式（12），式（13）可以表示為

對(duì)式（14）中的u、v、w求偏導(dǎo)使其為0。整理后設(shè)：

則可得：

可解得平面方程系數(shù)為

設(shè)Di為數(shù)據(jù)集M中點(diǎn)到擬合平面的距離，為數(shù)據(jù)集M中點(diǎn)到擬合平面的平均距離。精度評(píng)定距離標(biāo)準(zhǔn)偏差為

利用式（17）求得的u、v、w對(duì)平面特征向量進(jìn)行初始化，計(jì)算出數(shù)據(jù)集M中每個(gè)點(diǎn)至擬合平面的距離Di。當(dāng)Di≥ 2σ時(shí)，此點(diǎn)被認(rèn)為是異常點(diǎn)并刪除；反之,則保留。再利用所有保留下來(lái)的點(diǎn)通過(guò)式（13）—式（18）重新計(jì)算進(jìn)行迭代，直到剩下所有點(diǎn)的Di都在規(guī)定的閾值之內(nèi)，即Di都小于標(biāo)準(zhǔn)偏差的2 倍時(shí)為止，獲得最佳的u、v、w構(gòu)造平面特征向量，進(jìn)而獲得擬合平面。

對(duì)2.1 節(jié)中所得擬合直線進(jìn)行平面擬合，結(jié)果如圖6 所示。

圖6 擬合直線的平面擬合 Fig.6 Plane fitting by fitted lines

2.3 軌跡估計(jì)

由圖1 可知，目標(biāo)從任意方向進(jìn)入探測(cè)區(qū)域時(shí)，始終位于至少兩個(gè)陣列的中間區(qū)域。高速飛行目標(biāo)在探測(cè)區(qū)域內(nèi)可以視為直線運(yùn)動(dòng)，而空間中兩個(gè)不平行平面即可確定一條相交直線，即為估計(jì)得到的目標(biāo)運(yùn)動(dòng)軌跡，算法流程圖如圖7 所示。

圖7 估計(jì)目標(biāo)飛行軌跡的算法流程圖 Fig.7 Algorithm flowchart for target flying trajectory estimation

設(shè)陣列1 位于坐標(biāo)點(diǎn)（0,0,0），陣列3 位于坐標(biāo)點(diǎn)（0,100,0），對(duì)于單個(gè)陣列，分別使用兩幀數(shù)據(jù)與多幀數(shù)據(jù)進(jìn)行軌跡估計(jì)，目標(biāo)軌跡估計(jì)與顯示如圖8、9 所示。

如圖8（a）、8（b）所示，兩平面交線即為估計(jì)軌跡。如圖9（a）、9（b）所示，紅色點(diǎn)為真實(shí)軌跡點(diǎn)，綠色點(diǎn)為本文算法估計(jì)得到的軌跡點(diǎn)。仿真結(jié)果驗(yàn)證了兩個(gè)陣列分別僅使用兩幀數(shù)據(jù)即可估計(jì)出目標(biāo)軌跡?？梢酝ㄟ^(guò)增加使用的數(shù)據(jù)幀數(shù)提高估計(jì)精度，該方法可較為準(zhǔn)確地對(duì)目標(biāo)進(jìn)行軌跡估計(jì)。

圖8 不同幀數(shù)據(jù)的軌跡估計(jì)結(jié)果 Fig.8 Trajectory estimation results by different frames of data

圖9 真實(shí)軌跡與估計(jì)軌跡 Fig.9 Real trajectory and the estimated trajectory

3 誤差分析

3.1 軌跡誤差

高速飛行目標(biāo)在探測(cè)范圍內(nèi)為直線運(yùn)動(dòng)，軌跡估計(jì)結(jié)果也為空間中的直線方程。由于探測(cè)范圍已知，設(shè)在x軸上[0,100]范圍內(nèi)，同樣以d1為間隔，對(duì)估計(jì)的軌跡ME與真實(shí)軌跡MT進(jìn)行采樣，數(shù)據(jù)點(diǎn)數(shù)為n，分別表示為

軌跡估計(jì)誤差表示為

3.2 數(shù)據(jù)幀數(shù)對(duì)誤差的影響

由于直線擬合與空間擬合均存在誤差，估計(jì)得出的軌跡誤差與數(shù)據(jù)幀數(shù)有直接關(guān)系，如表1 所示。

表1 數(shù)據(jù)幀數(shù)對(duì)軌跡估計(jì)誤差影響 Table 1 The influence of the number of data frames on the trajectory estimation error

3.3 單陣列安裝誤差對(duì)軌跡估計(jì)結(jié)果影響

由2.3 節(jié)可知，使用兩個(gè)陣列（陣列1、陣列2）即可估計(jì)出目標(biāo)軌跡。若忽略陣列1 誤差，陣列2 由于陣列設(shè)備安裝產(chǎn)生方位估計(jì)誤差。設(shè)陣列2 對(duì)每幀數(shù)據(jù)進(jìn)行方位估計(jì)誤差相同，俯仰角與方位角誤差范圍均為[? 3°,3°]，誤差分辨率為0.1°。目標(biāo)飛行高度分別為60.9、90.9、120.9 m，誤差結(jié)果如圖10 所示。

圖10 單陣列波達(dá)方位估計(jì)誤差對(duì)軌跡估計(jì)誤差的影響 Fig.10 Influence of DOA estimation error of a single array on trajectory estimation error

3.4 雙陣列安裝誤差對(duì)軌跡估計(jì)結(jié)果影響

若考慮陣列1、2 由于陣列設(shè)備安裝產(chǎn)生的方位估計(jì)誤差。設(shè)陣列1、2 對(duì)每幀數(shù)據(jù)進(jìn)行方位估計(jì)誤差相同，俯仰角與方位角誤差范圍均為[? 3°,3°]，誤差仿真步長(zhǎng)為0.1°。目標(biāo)飛行高度分別為60.9、90.9、120.9 m，誤差結(jié)果如圖11 所示。

圖11 雙陣列波達(dá)方位估計(jì)誤差對(duì)軌跡估計(jì)誤差的影響 Fig.11 The influence of double-array DOA estimation error on trajectory estimation error

3.5 雙陣列隨機(jī)誤差對(duì)軌跡估計(jì)結(jié)果影響

同時(shí)考慮陣列1、2 產(chǎn)生的波達(dá)方位估計(jì)誤差對(duì)軌跡估計(jì)造成的影響。設(shè)陣列1、2 對(duì)每幀數(shù)據(jù)進(jìn)行波達(dá)方位估計(jì)誤差為隨機(jī)誤差，設(shè)俯仰角、方位角隨 機(jī)誤差 范圍均 為[? 5 °,5 °]、[? 7 °,7 °]、[? 9 °,9 °]，誤差仿真步長(zhǎng)為0.1°，進(jìn)行蒙特卡洛仿真。軌跡估計(jì)誤差數(shù)值分布與軌跡估計(jì)誤差統(tǒng)計(jì)分布如圖12～14 所示，軌跡估計(jì)平均誤差與最大誤 差如表2 所示。

圖12 波達(dá)方位估計(jì)隨機(jī)誤差在[-5°,5°]對(duì)軌跡估計(jì)影響 Fig.12 Influence of the random error of DOA estimation on the trajectory estimation in [-5°,5°]

圖13 波達(dá)方位估計(jì)隨機(jī)誤差在[-7°,7°]對(duì)軌跡估計(jì)影響 Fig.13 Influence of the random error of DOA estimation on the trajectory estimation in [-7°,7°]

圖14 波達(dá)方位估計(jì)隨機(jī)誤差在[-9°,9°]對(duì)軌跡估計(jì)影響 Fig.14 Influence of the random error of DOA estimation on the trajectory estimation in [-9°,9°]

表2 1 000次蒙特卡洛仿真中軌跡估計(jì)的平均誤差與最大誤差 Table 2 Average error and maximum error of trajectory esti-mation in 1 000 times of Monte Carlo simulations

4 結(jié)論

本文針對(duì)高速飛行目標(biāo)計(jì)算波達(dá)方位時(shí)由于目標(biāo)高速運(yùn)動(dòng)，不同陣列所接收到的信號(hào)并非聲源在同一時(shí)刻所發(fā)出，難以在時(shí)間維度上進(jìn)行對(duì)齊等問(wèn)題，不可避免存在誤差。借鑒點(diǎn)云數(shù)據(jù)處理的思想，提出了一種基于空間匹配的高速飛行目標(biāo)軌跡估計(jì)方法。對(duì)于單個(gè)陣列而言，對(duì)每一幀數(shù)據(jù)計(jì)算得到的波達(dá)方向在字典上進(jìn)行匹配，找到誤差最小的N個(gè)空間點(diǎn)進(jìn)行直線擬合。再對(duì)兩幀以上的數(shù)據(jù)所擬合的直線進(jìn)行空間點(diǎn)采樣，得到空間點(diǎn)數(shù)據(jù)集，通過(guò)空間平面擬合，得到估計(jì)平面。兩個(gè)陣列解算得到的平面相交，相交直線即為高速飛行目標(biāo)的軌跡估計(jì)結(jié)果。仿真結(jié)果表明：

（1）由于空間點(diǎn)直線擬合與平面擬合存在誤差，造成軌跡估計(jì)結(jié)果存在誤差。當(dāng)陣列使用的數(shù)據(jù)幀數(shù)增加時(shí)，軌跡估計(jì)的誤差減小。

（2）該算法可以減小單次波達(dá)方向估計(jì)誤差對(duì)軌跡估計(jì)的影響。單陣列與雙陣列的陣列安裝誤差造成的軌跡估計(jì)誤差隨飛行高度的降低而減小。目標(biāo)飛行高度為120.9 m 時(shí)，單陣列與雙陣列的陣列安裝誤差范圍為[? 2 °,2 °]，軌跡估計(jì)誤差小于10 m。對(duì)于隨機(jī)誤差，該算法表現(xiàn)出了較好的魯棒性。

仿真結(jié)果驗(yàn)證了高速飛行目標(biāo)的軌跡估計(jì)算法的有效性。