董 閣 郭良浩 徐 鵬 閆 超

(1 中國科學(xué)院聲學(xué)研究所 聲場聲信息國家重點實驗室 北京 100190)

(2 中國科學(xué)院大學(xué) 北京 100049)

0 引言

純方位目標運動分析方法僅利用目標方位信息估計目標的距離、速度和航向等運動狀態(tài)參數(shù)信息[1?2]。但根據(jù)純方位目標運動分析的可觀測性理論[3]，若觀測平臺未進行有效機動，目標運動狀態(tài)信息具有不完全可測性。當(dāng)目標進行勻速直線運動時，對于單一靜止的觀測平臺，僅目標航向信息是可觀測的。在實際應(yīng)用中，針對單一靜止的觀測平臺，如何獲取目標航向信息對目標的探測和識別具有重要的意義。

針對勻速直線運動目標，部分學(xué)者[4?6]利用目標方位與目標航向的幾何關(guān)系估計目標航向。根據(jù)目標航向的數(shù)學(xué)模型，張曉勇等[4]提出了最小二乘目標航向估計方法。孫少杰等[5]采用偽線性方法對目標航向的數(shù)學(xué)模型進行處理，提出了偽線性最小二乘目標航向估計方法。陳喆等[6]首先利用線性回歸理論分析了數(shù)據(jù)樣本的線性相關(guān)性，而后利用最小二乘方法進行航向估計。僅利用目標方位的目標航向估計方法收斂時間較長，航向估計精度較低。

考慮到海洋信道對聲傳播的影響，將信號處理方法與水聲傳播規(guī)律相結(jié)合，可以有效地提高聲吶設(shè)備的探測性能。在淺海水平不變遠場條件下，如何利用低頻聲場的聲強在距離-頻率平面上形成的干涉條紋提取目標運動狀態(tài)信息成為一個重要的研究方向。當(dāng)目標的運動軌跡相對于與觀測平臺存在最近通過距離時，余赟等[7]首先對干涉條紋進行Hough 變換，估計波導(dǎo)不變量和目標最近通過距離與目標速度的比值，而后估計目標航向。當(dāng)目標的運動軌跡相對于觀測平臺不存在最近通過距離時，王炳輝等[8]直接建立了目標航向、波導(dǎo)不變量和干涉條紋的關(guān)系模型，提出了利用Hough 變換對目標航向和波導(dǎo)不變量進行聯(lián)合估計的方法。同樣，針對無目標最近通過距離的情況，在淺海Pekeris波導(dǎo)條件下，余赟等[9]首先利用Radon變換提取干涉條紋，假定波導(dǎo)不變量值為1，采用二次測量法對目標進行估計。

相比于僅利用目標方位的航向估計的方法，利用淺海低頻聲場干涉條紋提取目標航向信息的方法有效地提高了目標航向的估計精度，但要求LOFAR譜中具有明顯的干涉條紋，且具有較高的信噪比。同時，以上方法均假定目標是勻速直線運動的，但在實際應(yīng)用中，在觀測時間范圍內(nèi)，很難保證目標持續(xù)進行勻速直線運動。因此，在目標進行變速直線運動時，結(jié)合淺海低頻聲場的聲場特征，如何快速有效地估計目標航向具有重要的研究價值。

近年來，大量國內(nèi)外學(xué)者致力于利用warping變換從簡正波互相關(guān)項中提取與目標距離相關(guān)的信息，再進行目標運動狀態(tài)參數(shù)的估計。根據(jù)簡正波水平波數(shù)差與波導(dǎo)不變量的關(guān)系，戚聿波等[10]提出了一種基于頻域β-warping變換算子的被動測距方法。在之前的研究中，將當(dāng)前時刻與初始時刻目標距離的比值定義為距離特征量，針對已有的利用LOFAR譜圖估計距離特征量的方法[11?12]估計精度不夠的問題，提出了利用頻域β-warping變換的距離特征量估計方法[13]，得到了修正純方位擴展卡爾曼濾波方法。但頻域β-warping變換需要已知波導(dǎo)不變量的值，在實際應(yīng)用中，很難獲取波導(dǎo)不變量的真實值。郭良浩等[14?15]選取特定的頻域β-warping變換算子的系數(shù)，研究表明，當(dāng)波導(dǎo)不變量的估計值為真實值的q倍時，若q、信號帶寬和頻域β-warping變換算子的系數(shù)之間滿足一定的關(guān)系，則簡正波互相關(guān)項通過頻域β-warping后的脈沖時延也為波導(dǎo)不變量取真實值時的q倍。

針對已有目標航向估計方法需要目標進行勻速直線運動的問題，本文將通過頻域β-warping變換得到的距離特征量信息引入到純方位目標航向估計方法中，提出了一種利用頻域β-warping變換的目標航向估計方法。該方法無需目標保持勻速直線運動，首先利用頻域β-warping變換估計距離特征量，而后根據(jù)目標方位和距離特征量利用漸近無偏最小二乘方法估計目標航向?？紤]到波導(dǎo)不變量的估計誤差對頻域β-warping的影響，本文選用文獻[14]提出的頻域β-warping 算子的系數(shù)進行頻域β-warping變換。數(shù)值仿真結(jié)果表明，在淺海Pekeris波導(dǎo)環(huán)境下，當(dāng)目標進行勻速直線運動時，相比于常規(guī)純方位目標航向估計方法，本文方法的收斂時間明顯縮短，航向估計精度更高。當(dāng)目標進行變速直線運動時，該方法依然具有較好的航向估計性能。另外，該方法航向估計性能與觀測時間內(nèi)的平均方位變化率密切相關(guān)。最后，利用一次實際海試數(shù)據(jù)驗證了本文方法的有效性。

1 利用頻域β-warping變換估計距離特征量

在淺海水平不變波導(dǎo)遠場條件下，接收信號自相關(guān)函數(shù)中的簡正波互相關(guān)的部分可以表示為

其中，S(f)為聲源幅度，An(f)為第n階簡正波的幅度，r為接收點和聲源的相對距離，γnm為僅與簡正波號數(shù)有關(guān)的常數(shù)，β為波導(dǎo)不變量的真實值。假設(shè)接收信號的有效頻帶為[f1,f2]，頻域βwarping 變換算子為C為常數(shù)。當(dāng)β1 時，C取內(nèi)的任意值均可保證頻域β-warping 變換重采樣后的頻帶包含接收信號的有效頻帶范圍[10]。假定波導(dǎo)不變量的估計值為且則其中對式(1)進行頻域β-warping變換可得

若所選參數(shù)滿足以下條件[15]：

式(2)可近似為

根據(jù)式(4)可以看出，第n階和第m階簡正波互相關(guān)項通過頻域β-warping變換得到的脈沖時延為

k時刻的距離特征量定義為k時刻的目標距離與初始時刻的目標距離的比值，即

對初始時刻和k時刻的接收信號自相關(guān)函數(shù)分別進行頻域β-warping變換，可得

將式(7)和式(8)代入式(6)可得

在實際應(yīng)用中，如果可以獲得粗略的水體聲速剖面信息，使得波導(dǎo)不變量的估計值在合理的范圍內(nèi)，可以保證式(3)的條件成立，則第n階和第m階簡正波互相關(guān)項通過頻域β-warping后的脈沖時延為波導(dǎo)不變量取真實值時的q倍[15]。由式(9)可以看到，k時刻的距離特征量可以表示為k時刻和初始時刻的第n階和第m階簡正波互相關(guān)項對應(yīng)的脈沖時延的比值。因此，波導(dǎo)不變量的估計誤差不會影響距離特征量的估計精度。

2 目標航向估計方法

2.1 漸近無偏最小二乘方法

目標與觀測平臺在x-y二維平面上的運動態(tài)勢示意圖如圖1所示。假定觀測平臺靜止且位于坐標原點，目標做直航運動，目標航向為?，k時刻的目標方位和距離特征量分別為θk和Mk。根據(jù)距離特征量的定義和正弦定理可得

根據(jù)幾何關(guān)系可以得到

將式(11)和式(12)代入式(10)得

圖1 目標與觀測平臺運動態(tài)勢示意圖Fig.1 The motion state of target and observer platform

假定k時刻的目標方位的真值和估計值分別為和距離特征量的真值和估計值分別為和將目標方位和距離特征量的真值代入式(13)，進一步整理得

其中，X=tan?。定義增廣狀態(tài)向量γ=[X,1]T，將目標方位和距離特征量的估計值代入式(14)可得

因為εθ,k足夠小，所以cosεθ,k ≈1，sinεθ,k ≈εθ,k，則式(16)化簡為

根據(jù)式(17)可以得到，均方誤差為

假設(shè)觀測總時長為K+1，則總的均方誤差為

其中，

令

可得

其中，

為了得到γ的無偏估計，將γTW γ限定為一個常數(shù)同時使達到最小[16?17]。上述問題可以轉(zhuǎn)化為如下條件極值問題：

針對此條件極值問題，可通過拉格朗日乘數(shù)法求解,可得[17]

對γ求偏導(dǎo)，并令偏導(dǎo)值為0，得到

根據(jù)式(23)可以看到，此問題為廣義特征值問題，λ為相對于W的特征值，γ為屬于λ的特征向量，最小特征值所對應(yīng)的特征向量即為該條件極值的解，進而得到目標航向正切值的估計值為

2.2 航向判別方法

目標航向范圍為[0,2π]，但反正切函數(shù)在[0,2π]內(nèi)具有多值性，所以獲得目標航向正切值的估計值后，需要進行航向判別。圖2為目標航向判別示意圖，下面以目標初始時刻位于第一象限為例，進行航向判別。

圖2 目標航向判別示意圖Fig.2 Target course discrimination

當(dāng)目標沿圖2中的虛線運動時，即目標方位不隨時間變化，此時可通過判斷目標與觀測平臺的相對距離的變化來估計目標航向。當(dāng)目標徑向遠離觀測平臺時，=θ0；當(dāng)目標徑向接近觀測平臺時，=θ0+π。

除目標沿圖2中的虛線運動的情形外，其余情形可通過目標方位隨時間的變化判斷目標航向的范圍。根據(jù)初始時刻的目標方位，將目標航向范圍劃分為兩部分。當(dāng)目標由初始位置向第I 部分運動時，目標方位范圍為[0,θ0)∪(θ0+π,2π]，在[0,θ0)和(θ0+π,2π]范圍內(nèi)，目標方位逐漸減小，航向范圍為[0,θ0)∪(θ0+π,2π]。當(dāng)目標由初始位置向第II部分運動時，目標方位范圍為(θ0,θ0+π)，目標方位逐漸增大，航向范圍為(θ0,θ0+π)。

根據(jù)以上的航向判別的劃分方法，若目標初始時刻位于第二象限，當(dāng)目標方位逐漸減小時，航向范圍為[0,θ0)∪(θ0+π,2π]；當(dāng)目標方位逐漸增大時，航向范圍為(θ0,θ0+π)。若目標初始時刻位于第三或第四象限，當(dāng)目標方位逐漸減小時，航向范圍為(θ0?π,θ0)；當(dāng)目標方位逐漸增大時，航向范圍為[0,θ0?π)∪(θ0,2π]。

2.3 性能評價指標

(1)收斂時間

則Kc為航向估計的收斂時間，是用來衡量航向估計方法的解算時間。

(2)均方根誤差

本文利用均方根誤差(Root mean square error,RMSE)來衡量目標航向的估計精度，k時刻的航向估計的均方根誤差定義為

3 仿真研究

采用Pekeris波導(dǎo)模型，海水深度為45 m，海水聲速為1505 m/s，海水密度為1.0 g/cm3；海底聲速為1596 m/s，海底密度為1.6 g/cm3，海底衰減系數(shù)為0.12 dB/λ。觀測平臺為位于水下20 m的64元水平均勻直線陣，陣元間距為0.5 m。水面目標為單位強度的點源，深度為3 m。信號頻率為100～500 Hz，采樣頻率為4000 Hz，單陣元接收信號的工作帶寬內(nèi)信噪比為?15 dB。觀測平臺靜止于坐標原點，目標做勻速直線運動，速度為5 m/s，航向為180?，初始距離為10 km，初始方位為60?。目標進行勻速直線運動時的目標與觀測平臺的運動態(tài)勢如圖3所示。觀測時間為20 min，以1.024 s作為一幀進行數(shù)據(jù)處理。

圖3 目標進行勻速直線運動時的目標與觀測平臺的運動態(tài)勢Fig.3 The motion state of target and observer platform for the target with uniform linear motion

圖4 目標進行勻速直線運動時的目標航向估計結(jié)果Fig.4 Course estimation results for the target with uniform linear motion

圖4(a)為最小方差無失真響應(yīng)波束形成方法得到的目標方位歷程。圖4(b)為頻域β-warping變換后的時域序列，波導(dǎo)不變量取為1，處理頻段為200～400 Hz。通過提取的不同時刻的脈沖序列時延，估計距離特征量，結(jié)果如圖4(c)所示。本文利用文獻[13]提出的距離特征量平均估計誤差來衡量距離特征量的估計精度。距離特征量的平均估計誤差定義為

頻域β-warping變換方法得到的距離特征量平均估計誤差為1.21%。本文將文獻[5]中提出的純方位最小二乘目標航向估計方法作為對比，分析本文提出的利用目標方位和距離特征量的漸進無偏最小二乘目標航向估計方法的性能。根據(jù)以上估計的目標方位和距離特征量，分別利用文獻[5]和本文方法對目標航向進行估計，結(jié)果如圖4(d)所示。文獻[5]和本文方法的收斂時間分別為13.1 min和4.2 min，第20 min，文獻[5]和本文方法的航向估計誤差分別為2.57?和2.30?。

以上仿真為目標進行勻速直線運動的情形，下面分析本文方法在目標進行變速直線運動的條件下的性能。目標初始速度為5 m/s，初始距離為10 km，初始方位為60?，航向為180?。目標先進行勻速直線運動，持續(xù)時間為5 min，然后進行勻加速直線運動，加速度為0.01 m/s2，持續(xù)時間為10 min，之后繼續(xù)進行勻速直線運。其余仿真條件與目標進行勻速直線運動時的仿真條件一致。目標進行變速直線運動時的目標與觀測平臺的運動態(tài)勢如圖5所示。圖6為觀測時間范圍內(nèi)的目標速度。

圖5 目標進行變速直線運動時的目標與觀測平臺的運動態(tài)勢Fig.5 The motion state of target and observer platform for the target with variable velocity linear motion

圖6 目標速度Fig.6 Target velocity

圖7(a)和圖7(b)分別為目標方位歷程圖和頻域β-warping變換后的時域序列，圖7(c)為目標距離特征量的估計結(jié)果，距離特征量的平均估計誤差為1.11%。目標航向估計結(jié)果如圖7(d)所示，本文方法的收斂時間為4.2 min，第20 min，本文方法的航向估計誤差為1.30?。由于文獻[5]的方法要求目標進行勻速直線運動，所以針對目標進行變速直線運動的情形，文獻[5]的方法無法正確估計目標航向。而本文方法由于引入了距離特征量信息，無需目標保持勻速直線運動，在目標進行變速直線運動時，依然具有較好的航向估計性能。

下面分析兩種方法在不同觀測誤差條件下的估計性能。目標與觀測平臺的運動態(tài)勢與目標進行勻速直線運動時的仿真條件相同，目標方位和距離特征量的觀測誤差均服從零均值的高斯分布，蒙特卡洛仿真的次數(shù)為100次。在不同觀測誤差條件下，兩種方法的性能如表1所示。根據(jù)表1可以看出，當(dāng)距離特征量估計誤差一定時，隨著方位估計誤差的增大，兩種方法的收斂時間逐漸增大，航向估計精度逐漸降低，但相比于文獻[5]的方法，本文方法的收斂時間更短，估計精度更高。當(dāng)方位估計誤差一定時，隨著距離特征量估計誤差的增大，本文方法的收斂時間逐漸增大，航向估計精度逐漸降低。本文方法由于引入了距離特征量信息，獲得了更好的航向估計性能，但當(dāng)距離特征量誤差較大時，本文方法的性能下降，收斂時間大于文獻[5]的方法，且估計精度更低。

圖7 目標進行變速直線運動時的目標航向估計結(jié)果Fig.7 Course estimation results for the target with variable velocity linear motion

表1 在不同觀測誤差條件下的兩種方法性能Table1 The estimation performance of the method with different observation error

下面分別研究兩種方法在不同目標航向、不同初始距離和不同目標速度條件下的方法性能，其余仿真條件與目標進行勻速直線運動時的仿真條件相同。目標方位和距離特征量的觀測誤差均服從零均值的高斯分布，方位估計誤差標準差σθ=0.5?，距離特征量估計誤差標準差σm=0.05，蒙特卡洛仿真的次數(shù)為100次。定義θr為觀測時間范圍內(nèi)的平均方位變化率，即

其中，T為采樣時間間隔。在不同目標航向條件下，兩種方法的航向估計均方根誤差如圖8所示，可以看到，本文方法相比于文獻[5]方法的收斂時間更短，航向估計均方根誤差更小。在不同目標航向條件下，觀測時間內(nèi)的平均方位變化率如圖9所示，可以看到，除目標航向為60?和240?兩種情形外，當(dāng)目標航向接近60?或240?時，平均方位變化率較小，對應(yīng)的航向均方根誤差較大。當(dāng)目標平均方位變化率較大時，對應(yīng)的航向均方根誤差較小。當(dāng)目標航向為60?或240?時，即目標沿圖2中的虛線進行運動，則觀測時間內(nèi)的目標方位與不隨時間變化，目標平均方位變化率為0，此時根據(jù)航向判別方法即可確定目標航向，無需利用目標方位和距離特征量進行估計，因此，在平均方位變化率為0時，依然具有較好的航向估計性能。通過以上分析可以看到，兩種方法的航向估計性能與觀測時間內(nèi)的平均方位變化率密切相關(guān)。不同初始距離條件下，兩種方法的航向估計均方根誤差如圖10所示，觀測時間內(nèi)的平均方位變化率如圖11所示。平均方位變化率是隨著初始距離的增大而逐漸減小，兩種方法的航向估計均方根誤差隨著初始距離的增大而逐漸增大。不同目標速度條件下，兩種方法的航向估計均方根誤差如圖12所示，觀測時間內(nèi)的平均方位變化率如圖13所示。平均方位變化率是隨著目標速度的增大而逐漸增大，兩種方法的航向估計均方根誤差隨著目標速度的增大而減小。因此，根據(jù)不同初始距離條件和不同目標速度條件下的航向估計均方根誤差與平均方位變化率的結(jié)果，同樣可以得到，兩種方法的航向估計性能與觀測時間內(nèi)的平均方位變化率密切相關(guān)，當(dāng)平均方位變化率較大時，可以獲得較好的航向估計性能。

圖8 不同目標航向條件下的航向估計均方根誤差Fig.8 The RMSE of course estimation with different target course

圖9 不同目標航向條件下的平均方位變化率Fig.9 The average bearing rate with different target course

圖10 不同初始距離條件下的航向估計均方根誤差Fig.10 The RMSE of course estimation with different initial target range

圖11 不同初始距離條件下的平均方位變化率Fig.11 The average bearing rate with different initial target range

圖12 不同目標速度條件下的航向估計均方根誤差Fig.12 The RMSE of course estimation with different target velocity

圖13 不同目標速度條件下的平均方位變化率Fig.13 The average bearing rate with different target velocity

4 實驗研究

下面利用一次實際海試數(shù)據(jù)對本文方法進行驗證分析。2005年6月，中國科學(xué)院聲學(xué)研究所聲場聲信息國家重點實驗室在黃海海域進行了一次海底水平陣的聲學(xué)實驗，合作目標為一艘小型漁船。水平陣的有效陣元數(shù)為43個，陣元間距約為1.5 m。實驗海區(qū)聲速剖面如圖14所示，海水深度為27 m。采樣頻率為4000 Hz，觀測時間為10 min。實驗過程中，目標近似做勻速直線運動，目標運動航向約為142?，運動態(tài)勢如圖15所示。

目標方位歷程和頻域β-warping變換后的時域序列分別如圖16(a)和圖16(b)所示，距離特征量的估計結(jié)果如圖16(c)所示，距離特征量的平均估計誤差為4.37%。目標航向估計結(jié)果如圖16(d)所示，通過實驗數(shù)據(jù)處理結(jié)果可以看到，本文方法有效地縮短了航向估計的收斂時間，與仿真研究得到的結(jié)論一致。文獻[5]和本文方法的收斂時間分別為6.5 min和3.5 min。

圖14 實驗海區(qū)聲速剖面Fig.14 The sound speed profile of the experiment

圖15 實驗過程中的目標與觀測平臺的運動態(tài)勢Fig.15 The motion state of target and observer platform for the experiment

圖16 實驗數(shù)據(jù)處理的目標航向估計結(jié)果Fig.16 The course estimation results for the experiment

5 結(jié)論

針對單一靜止觀測平臺的目標航向估計問題，常規(guī)純方位目標航向估計方法不僅需要目標保持勻速直線運動，而且航向估計的收斂時間較長，估計精度較低。針對此問題，本文提出了一種利用頻域β-warping變換的淺海目標航向估計方法。該方法首先通過頻域β-warping變換估計距離特征量，根據(jù)幾何關(guān)系，利用目標方位和距離特征量確定目標航向，再利用漸近無偏最小二乘方法對目標航向進行估計。數(shù)值仿真和實驗數(shù)據(jù)處理結(jié)果表明，在淺海水平不變波導(dǎo)遠場條件下，對于單一靜止的觀測平臺，該方法無需目標保持勻速直線運動，對直航運動的目標均可以進行可靠地航向估計。同時，該方法航向估計性能與觀測時間內(nèi)的平均方位變化率密切相關(guān)。復(fù)雜淺海環(huán)境下的方法的適用性是下一步研究的重點。

致謝感謝參加實驗的全體人員，他們的辛勤工作為本文的科學(xué)研究提供了寶貴的實驗數(shù)據(jù)。