李 鑫 郭良浩 章偉裕 徐 鵬 陳連榮 段正辰
(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)
(2 中國科學院大學 北京 100049)
(3 中國人民解放軍92330 部隊 青島 266102)
目標運動分析(Target motion analysis,TMA)一直以來都是水聲領域研究的熱點問題,其中純方位目標運動分析[1?9](Bearings only target motion analysis,BOTMA)方法只需利用目標的方位信息就可以實現對目標位置、運動參數的估計,因而實現相對簡單,在工程上得到了廣泛應用。由于觀測站測量的目標方位是BOTMA 方法唯一可利用的信息,因此觀測站測量的目標方位精度將影響BOTMA 的性能,通常隨著測量方位誤差的增加,BOTMA性能快速下降。
為了解決目標方位估計誤差造成的BOTMA性能下降的問題,提高基陣被動測向精度是一個重要的研究方向。已有相關的科研人員研究了影響基陣被動測向精度的因素,并取得了一定的研究成果。宮在曉等[10]指出參考聲速的選取是否準確對傳統的水平線列陣被動測向結果有很大的影響,此影響與信號的入射方向有關:對目標信號近基陣端射方向入射時產生的影響大于信號近基陣正橫方向入射時的影響,并基于簡正波理論指出由于淺海中聲波是以多模態形式傳播,參考聲速與接收陣處的相速度密切相關,其中相速度大小是由聲源頻率、聲源距離、收發深度和傳播信道等因素決定。由于被動測向時聲源距離未知,因此在對目標角度估計時選取的參考聲速與接收陣處的相速度往往存在偏差,從而影響基陣測向精度[11]。
近年來,已有學者對減小因參考聲速選取偏差導致的測向誤差,從而提高基陣測向精度的方法開展了大量研究。宮在曉等[10]將匹配場方法應用于聲場測向,通過多次計算拷貝場,將與觀測場最相關的拷貝場的聲源角度作為估計角度,無需進行波束形成等信號處理方法進行角度估計,所以不需要考慮參考聲速,避免了因參考聲速選取不準確造成的目標測向誤差,但該方法要求對所在海洋環境參數準確已知,且通常需要對距離、角度進行極大值的搜索,導致了計算量的大量增加。文獻[12–14]提出了聲速無關水下波達方向估計方法,通過使用兩條線陣,對每條線陣的接收信號分別進行處理,利用兩條線陣之間的幾何關系,得到與聲速無關的估計結果。該方法雖然能消除聲速對于估計結果的影響,但需要復雜度很高的匹配算法,分別對兩條線陣的數據進行處理,計算量大,估計速度慢。吳俊楠[11]提出了令水平線陣繞軸向機動,根據機動前后陣列測向結果相對于參考聲速的變化曲線,修正由于參考聲速選取偏差導致的測向誤差的方法。當目標靜止或運動較慢時參考聲速短時間內變化不大,此方法校準獲得的參考聲速可以提高此范圍內測向的精度。但是對于目標運動較快時,參考聲速變化較大,此時校準參考聲速需要陣列多次機動,代價較高。
為了校準參考聲速,提高目標角度估計精度,本文提出了將聲傳播相速度作為估計狀態量,改善了由于參考聲速選取偏差導致的測向誤差的問題,同時提高了BOTMA 的跟蹤性能。文章首先分析了參考聲速對傳統水平線列陣被動測向精度的影響,給出了在確知環境下的兩種理論參考聲速計算方法,而后推導了一種水平線列陣方位-相速度聯合的純方位擴展卡爾曼濾波(Cs bearing only extended Kalman filter,Cs-BO-EKF)方法,仿真結果表明,Cs-BO-EKF 方法在跟蹤過程中實時估計相速度,提高了測向精度,相比純方位擴展卡爾曼濾波(Bearing only extended Kalman filter,BOEKF)算法,降低了受測向誤差的影響,具有更高的跟蹤精度及穩健的跟蹤性能,具有較好的應用前景。
傳統的水平均勻直線陣被動測向法,是基于信號以平面波傳播的假設[15?16]:各陣元接收到同一目標的信號波形相同,僅是時延不同。如圖1 所示,一個由M個陣元組成的均勻直線陣,陣元間距為d,平面波入射方向與基陣法線方向夾角為θ0,其中令沿著基陣正橫方向入射的角度為0?,沿著基陣左端入射的角度為?90?,沿著基陣右端入射的角度為90?。由于各陣元位置不同,使得信號源到達各個陣元具有不同的時延。

圖1 波束形成示意圖Fig.1 Schematic diagram of beamforming
設首陣元為參考陣元,接收到的信號為x1(t),則第m號陣元接收到的信號為
其中,τm(θ0)為第m號陣元相對于參考陣元的時延,
其中,c為參考聲速。聲信號到達不同陣元的時延是確定值,假設對應于信號真實入射角θ0的參考聲速為c0,在實際計算中如果選用參考聲速為c1,根據式(2)計算得到的角度為θ1,則它們滿足下述公式:
由式(4)可知由于選取的參考聲速c1相對于c0產生了?c的偏差,導致測向結果相對于真值產生了?θ的偏差。當信號以不同的角度入射時,選取的參考聲速在不同的偏差的情況下所產生的測向偏差如圖2所示。

圖2 不同入射角度下參考聲速產生的測向誤差Fig.2 Direction finding error caused by reference sound velocity at different incident angles
可以看到水平線列陣的測向誤差與參考聲速和信號的入射角度有關,當信號入射方向靠近基陣的正橫方向時,參考聲速的選取偏差不會導致明顯的測向誤差;當信號入射方向靠近基陣的端射方向時,參考聲速的選取偏差會導致明顯的測向誤差,且參考聲速偏差越大,測向誤差越大。
根據簡正波理論,淺海中聲波是以多模態形式傳播,參考聲速與接收陣處的相速度密切相關,本節分別從簡正波幅度加權和波束形成兩個方面,給出在確知環境下的理論參考聲速的計算方法。
根據簡正波理論,聲場中任一點處的聲壓是聲源激發的多號簡正波疊加的結果,即[17]
其中,pl為聲源激發的第號簡正波,令cpl為第l號簡正波的相速度。將各號簡正波的相對能量作為加權因子對相應號數的相速度進行加權并求和,作為該位置處的理論參考聲速cs,即
已知目標距離rs與方位θs,選取聲速搜索序列,如[c1,c2,···,cn],由于聲信號到達不同陣元間的時延值τ固定,即
遍歷選取的聲速值進行波束形成計算,得到波束形成角度估計序列[θ1,θ2,···,θN],選取與目標角度相對誤差最小的角度估計值對應的聲速作為該目標距離下的參考聲速。改變聲源與接收陣的水平距離,得到不同距離下的理論參考聲速。
BOTMA 方法的基本原理是對運動目標,通過波束形成等方法得到方位觀測序列,并利用該方位序列估計目標的位置和運動參數。為了分析方便,假設目標與觀測站均在二維平面上,目標與觀測站的運動態勢圖如圖3所示。

圖3 目標與觀測站的運動態勢示意圖Fig.3 Movement situation diagram of target and observation station
將連續時間系統表示為離散形式,直角坐標系下,k時刻目標運動狀態XT(k)和觀測站運動狀態XW(k)分別表示為
其中,rxt(k)和rxw(k)分別為k時刻目標和觀測站沿X軸方向距離;ryt(k)和ryw(k)為k時刻目標和觀測站沿Y軸方向的距離;vxt(k)和vxw(k)分別為k時刻目標和觀測站沿X軸方向的速度;vyt(k)和vyw(k)分別為k時刻目標和觀測站沿Y軸方向的速度。
則目標與觀測站的相對運動狀態為XR(k)=XT(k)?XW(k):
其中,rx(k)和vx(k)分別為k時刻目標相對觀測站沿X軸方向的距離和速度;ry(k)和vy(k)分別為k時刻目標相對觀測站沿Y軸方向的距離和速度。
對于勻速直線運動的目標,目標與觀測站的相對運動狀態方程為
T表示采樣時間間隔,w(k)表示過程噪聲,一般假設其滿足零均值高斯分布,協方差矩陣為Q。目標方位的觀測方程為
其中,εβ(k)表示測量噪聲,滿足零均值高斯分布,εβ(k)~N(0,(k))。由于系統的觀測方程非線性,因此對其進行泰勒級數展開取其一階項來進行測量方程的線性化。一階泰勒級數展開:
因此BO-EKF算法的迭代公式如下[2]:
其中,Xk+1|k為k+1時刻狀態向量的預測值,Pk+1|k為k+1時刻狀態向量的協方差矩陣的預測值,Gk+1為k+1時刻的卡爾曼增益,Xk+1|k+1為k+1 時刻狀態向量的估計值,Pk+1|k+1為k+1時刻狀態向量的協方差矩陣的估計值。
對于勻速直線運動的目標,為校準參考聲速,提高測向精度,將相速度cs(k)也作為估計狀態量,目標與觀測站的相對運動狀態方程為
k時刻,選取參考聲速為c0,β0(k)為按照此參考聲速,利用波束形成獲得的舷角估計值,βb(k)為對應相速度cs(k)的舷角真實值,則滿足:
目標方位角與舷角的關系如圖4 所示,其中β(k)為第k時刻目標方位真值,βb(k)為第k時刻目標相對于觀測站的舷角真實值,βv(k)第k時刻觀測站沿x軸方向的速度vxw(k) 與觀測站速度的夾角。

圖4 角度示意圖Fig.4 Angle diagram
由圖4可知,目標方位角與舷角滿足:
因此觀測方程為
可以表示為
其中,εβ(k)表示測量噪聲,滿足零均值高斯分布,εβ(k)~N(0,(k))。由于系統的觀測方程非線性,因此對其進行泰勒級數展開取其一階項來進行測量方程的線性化。一階泰勒級數展開:
Cs-BO-EKF算法的迭代公式如下:
其中,Xk+1|k為k+1時刻狀態向量的預測值,Pk+1|k為k+1時刻狀態向量的協方差矩陣的預測值,Gk+1為k+1時刻的卡爾曼增益,Xk+1|k+1為k+1 時刻狀態向量的估計值,Pk+1|k+1為k+1時刻狀態向量的協方差矩陣的估計值。
以Pekeris 波導為例,海深100 m,水中聲速1500 m/s,海水密度1.0 g/cm3,海底聲速為1650 m/s,海底密度為1.6 g/cm3,海底衰減系數為0.3 dB/λ。接收陣由24 個間隔0.25 m 的陣元組成的水平均勻直線陣,深度為50 m。目標信號頻率為450~550 Hz,深度為25 m,角度為30?。分別利用簡正波幅度加權方法和波束形成方法計算不同距離處的理論參考聲速如圖5所示。

圖5 理論參考聲速Fig.5 Theoretical reference sound velocity
從仿真結果可以看出,兩種方法計算得到的理論參考聲速隨距離變化的趨勢一致,但使用波束形成方法計算的理論參考聲速值在局部區域內存在抖動變化的情況,這是由于在不同距離處不同號數的簡正波可能存在干涉現象;采用簡正波幅度加權方法計算的理論參考聲速值不存在局部抖動變化的情況,這是由于該方法不考慮簡正波干涉情況,只是根據不同號數簡正波的振幅對對應號數的相速度進行加權。
相速度大小是由聲源頻率、聲源距離、收發深度、海洋環境傳播信道等因素決定,為探究該方法對于海洋環境的適應性,分別以不同海洋環境下陣列接收到的有效簡正波號數多少進行討論。
4.2.1 有效簡正波號數相對較多情況
以Pekeris 波導為例,海深100 m,水中聲速1500 m/s,海水密度1.0 g/cm3,海底聲速為1650 m/s,海底密度為1.6 g/cm3,海底衰減系數為0.3 dB/λ。觀測站為24 元水平均勻直線陣,陣元間距為0.25 m,深度為30 m。目標為單位強度的點源,信號頻率為450~550 Hz,深度為30 m,信噪比為0 dB。在當前海洋環境下,計算與觀測站等深度,距離目標8200 m 處接收到的簡正波號數共有60 號,其中高階簡正波(第24~60 號簡正波)在遠距離傳播時由于衰減過多導致其能量遠小于低階簡正波能量,因此只選取前23 號簡正波,計算各號簡正波的相對能量,統計相對能量不低于0.1 倍的最大相對能量(圖6 中紅線)對應的簡正波號數為該環境下的有效簡正波,如圖6所示,該海洋環境下對應的有效簡正波號數共計13號。

圖6 簡正波能量分布Fig.6 Normal wave energy distribution
目標做勻速直線運動,目標速度為10 m/s,航向為170?。觀測站速度為4 m/s,初始航向為90?,目標與觀測站初始距離為12 km。觀測站前300 s做勻速直線運動,然后以0.5?/s 的角速度沿逆時針方向旋轉150 s,之后做勻速直線運動。采樣周期1 s,總觀測時間50 min。目標與觀測站的航跡如圖7所示,目標相對于觀測站的舷角變化如圖8 所示。參考聲速選取觀測站所處深度處的海水聲速即1500 m/s,進行常規波束形成方法得到的目標舷角誤差如圖9所示。由圖8可以看到目標運動的第7~20 min、第35~50 min,目標位于觀測站的近端射方向附近,因此參考聲速的選取偏差會產生較大的測向誤差。

圖7 目標與觀測站的航跡Fig.7 Movement situation of target and observation station

圖8 目標舷角歷程Fig.8 Angle change of target

圖9 常規波束形成法得到的測向誤差Fig.9 Error of direction finding obtained by conventional beamforming method
利用參考聲速為1500 m/s 進行常規波束形成得到的測向結果,結合BO-EKF 算法和Cs-BOEKF 算法估計得到的目標跟蹤結果如圖10 所示,目標距離和速度的跟蹤結果如圖11和圖12所示。

圖10 兩種濾波算法的跟蹤結果Fig.10 Tracking results of two filtering algorithms

圖11 目標距離跟蹤結果Fig.11 Target range tracking results
由圖10~圖12 可以看出,利用BO-EKF 算法得到的目標跟蹤結果在35 min 之后開始發散,從圖8、圖9 中可以看到,由于在35 min 之后目標逐漸近觀測站端射方向入射,采用參考聲速1500 m/s 進行波束形成的測向結果誤差逐漸變大,因此采用其結果直接作為觀測值時,會因為角度估計誤差過大導致BO-EKF 結果不收斂。由于Cs-BO-EKF 算法引入相速度作為狀態估計量,可以迭代出整個跟蹤過程中的相速度估計值,以此對測向結果進行校準,對比看出,Cs-BO-EKF 算法的跟蹤性能更加穩健,且跟蹤結果誤差更低。
圖13 是Cs-BO-EKF 算法對相速度估計的結果,圖中Cp1 是利用簡正波幅度加權方法計算得到的相速度,Cp2 是利用波束形成方法計算得到的相速度。

圖13 相速度估計曲線Fig.13 Phase velocity estimation curve
可以看到,在前5 min 內,由于目標始終位于觀測站基陣近正橫方向,相速度估計值變化對角度估計值變化的影響較小,且此時濾波算法未收斂,因此會產生較大偏差。第5~7.5 min 內,由于觀測站轉動,使得目標相對于觀測站基陣的方向逐漸趨于端射方向,相速度估計值的變化對角度估計值影響較大,因此濾波算法中迭代的相速度值會逐漸收斂到相速度理論值。之后由于目標入射方向逐漸靠近基陣正橫方向,此時參考聲速的選取誤差對測向結果的影響變小,故Cs-BO-EKF 在此區間內對相速度的估計值沒有產生較大的變化,與理論值誤差變大。隨后在35 min 之后,由于目標入射方向逐漸靠近基陣端射方向,此時Cs-BO-EKF 算法對相速度的估計值與理論值又較為接近。
圖14 是利用Cs-BO-EKF 算法估計的相速度值進行常規波束形成得到的測向結果誤差,和參考聲速為1500 m/s 進行常規波束形成得到的測向結果誤差。當目標入射角度位于區間[?58?,?45?]利用相速度估計值得到的測向結果平均誤差為0.530?,利用參考聲速1500 m/s 得到的測向結果平均誤差為1.270?;當目標入射角度位于區間[45?,75?]利用相速度估計值得到的測向結果平均誤差為?0.768?,利用參考聲速1500 m/s 得到的測向結果平均誤差為?2.173?。可以看到利用Cs-BO-EKF算法提高了目標近端射方向入射時測向結果精度。

圖14 測向結果誤差曲線Fig.14 Error curve of direction finding
4.2.2 有效簡正波號數相對較少情況
圖15所示中等水文條件下的聲速剖面,其中海面聲速為1530 m/s,聲速梯度為?0.05 m2/s,其余參數不變。在當前海洋環境下,計算與觀測站等深度,距離目標8200 m 處接收到的簡正波號數共有58 號,其中高階簡正波(第24~58 號簡正波)在遠距離傳播時由于衰減過多導致其能量遠小于低階簡正波能量,因此只選取前23 號簡正波,計算各號簡正波的相對能量,如圖16所示該海洋環境下對應的有效簡正波號數共計6號。

圖15 聲速剖面Fig.15 Sound speed profile

圖16 簡正波能量分布Fig.16 Normal wave energy distribution
參考聲速選取觀測站所處深度處的海水聲速即1528 m/s,進行常規波束形成。BO-EKF 算法和Cs-BO-EKF 算法估計得到的目標跟蹤結果如圖17所示,目標距離和速度的跟蹤結果如圖18 和圖19所示。

圖17 兩種濾波算法的跟蹤結果Fig.17 Tracking results of two filtering algorithms

圖18 目標距離跟蹤結果Fig.18 Target range tracking results

圖19 目標速度跟蹤結果Fig.19 Target speed tracking results
可以看到,CS-BO-EKF 算法對目標的距離,速度的估計的精度都高于BO-EKF 算法,且Cs-BOEKF 算法的跟蹤性能更加穩健。圖20 是Cs-BOEKF算法對相速度估計的結果,可以看到相速度的估計值與理論計算的相速度值在目標入射方向位于觀測站基陣近端射方向時較為吻合。圖21是利用Cs-BO-EKF 算法估計的相速度進行常規波束形成得到的測向結果相對誤差,和參考聲速為1528 m/s進行常規波束形成得到的測向結果相對誤差。當目標入射角度位于區間[?58?,?45?]利用相速度估計值得到的測向結果平均誤差為0.579?,利用參考聲速1528 m/s 得到的測向結果平均誤差為0.997?;當目標入射角度位于區間[45?,75?]利用相速度估計值得到的測向結果平均誤差為?0.235?,利用參考聲速1528 m/s 得到的測向結果平均誤差為?1.547?,可以看到利用Cs-BO-EKF 算法提高了目標近端射方向入射時測向結果精度。

圖20 相速度估計曲線Fig.20 Phase velocity estimation curve

圖21 測向結果誤差曲線Fig.21 Error curve of direction finding
結合圖13和圖20,可以看到在不同的海洋環境下,由Cs-BO-EKF 方法估計的相速度與理論的相速度在目標近端射方向入射時均較為一致。
針對水平線列陣因參考聲速選取偏差造成的測向誤差,從而導致BOTMA 方法在測向誤差較大時跟蹤結果不收斂的問題,本文提出了一種水平線列陣Cs-BO-EKF 方法,引入相速度作為估計狀態量,實時校準參考聲速,從而提高測向精度。淺海傳播條件下的數值仿真結果表明,Cs-BO-EKF 方法相比BO-EKF 方法的跟蹤性能更加穩健,跟蹤精度更高。當目標接近水平線列陣端射方向時,本文所提方法的相速度估計值與理論計算的相速度值較為吻合,從而提高了目標近端射方向入射的測向精度。然而為滿足BOTMA 系統的可觀測性,該方法要求觀測站至少進行一次有效的機動,在實際中不適用于隱蔽跟蹤。