李長文,趙建昕

(海軍潛艇學院,山東 青島 266042)

用純方位方法對等速直航目標運動要素進行估計,是一個經典問題,其方法為觀測站使用被動聲納觀測聲源發(fā)出的信號并給出聲源的方位,然后利用一定的原理對目標運動要素進行估計或對其運動狀態(tài)進行預測等。根據文獻[1-4],傳統作法一般不考慮聲音在海水中傳播的時間,即認為聲速為無限值,可觀測的目標方位即當前時刻噪聲源的實際方位?？紤]到聲音在海水中傳播總需要一定時間,對于運動目標,可被觀測的目標方位并非目標噪聲源當前時刻的實際方位,所觀測的方位必與目標運動要素及聲速有關,基于無誤差條件下的這個關系,就可以研究純方位條件下目標運動要素的估計問題。

聲源發(fā)出的瞬時噪聲經水聲信道作用后,于觀測站處所檢測到的信號一般為一個非瞬時信號[5-6],因此由檢測裝置給出的觀測到目標的時刻并非嚴格的聲音沿特定途徑到達觀測站的時刻,所以自聲源發(fā)出噪聲信號至觀測站檢測到信號的時間延遲是由比較復雜的物理過程決定的,其所依賴的變量至少應包括聲源深度、觀測站深度、觀測站到聲源水平距離、水聲環(huán)境等,不能簡單地認為是聲音沿直線勻速運動到觀測站的,即不能人為地設想聲音自聲源傳播至觀測點的速度這一概念。因此本文所說的聲音在海水中的傳播速度,是指一個假想的平均速度Vs(t),即 t時刻到達觀測站的聲信號對應的聲源至觀測站的水平距離與信號發(fā)出至被檢測到的延遲時間之比,對于t時刻所觀測到的聲音信號,將被認為以這個假想的平均速度于觀測站所處水平面上自聲源位置點處向外勻速直線運動。若不考慮聲音信號自發(fā)出至被觀測到的時間延遲,則相當于認為Vs(t)= ∞。

1 等速直航目標可被觀測的方位真值β(t)的規(guī)律及運動要素的可觀測性

假設目標作等速直航運動,以 V0、C0表示其速率與航向,以D0、β0表示 0時刻目標相對觀測站的實際距離與方位真值,稱(V0, C0, D0, β0)為等速直航目標的運動要素。

假設目標及觀測站在同一水平面上運動,如圖1,以0時刻觀測站實際位置O 為原點,正北( N )方向為Ox 軸正方向,正東 ( E )方向為Oy 軸正方向,建立一個平面直角坐標系。角的符號符合標準數學教科書的規(guī)定,即自Ox軸旋轉90 °至Oy軸的方向為正方向,此處與順時針為正的航海習慣一致。為了表達式簡潔,平面上的點、向量及其運算可以用復數表示,其中arg(z)表示復數z的[0, 2π)上的輻角主值,i 為虛數單位,Re(z),表示復數 z的實部,Im(z)表示復數 z的虛部。

圖1 無誤差的觀測站A、目標B、可觀測點 P 的位置圖

以A(t)= xA(t)+ i yA(t)表示t時刻觀測站實際位置,忽略觀測站自身定位誤差,則 A(t)對于觀測者是確定的已知函數。t時刻目標實際位置點為

以 D0(t)= | A(t)B(t)|、β0(t)= arg(A(t)B(t))分別表示t時刻目標相對觀測站實際距離及方位真值,其中A(t)B(t)= B(t)– A(t)。

以β(t)表示 t時刻所觀測到的目標方位,對于Vs(t)＜ ∞的情形,β(t)非當前時刻目標的真實方位β0(t),而是目標于當前時刻之前的某個位置點 P(t)所發(fā)出的聲音信號到達觀測站的到達角(到達方向的反方向相對正北方向的角),即

其中 A(t)P(t)= P(t)– A(t),以下稱β(t)為可被觀測的方位真值。若認為 Vs(t)= ∞,則β(t)= β0(t)。

P(t)相當于 t 時刻可被觀測到的聲源位置點,記D(t)= | A(t)P(t)|,則聲音自P(t)到達觀測站的時間為D(t)/ Vs(t),因此

此方程應對應于唯一的(D(t), β(t)),下面研究其計算公式。由式(1)及式(3)可得

記 k(t)= V0/ Vs(t),則 0 ≤ k(t)＜＜ 1,且

因此

計算后整理成關于D(t)的一元二次方程

其中 a(t)= 1 – k2(t)＞ 0,b(t)的計算公式為

因為 D(t)≥ 0,所以

由式(4)可得

這就是利用(t, V0, C0, D0, β0, Vs(t), A(t))計算 D(t)、β(t)的公式。

將式(4)改寫為

此式兩端的實部與虛部之比對應的等式為

利用這一關系可以研究目標運動要素的可觀測性。

若觀測站靜止,則xA(t)= yA(t)= 0,式(7)變?yōu)?/p>

對于Vs(t)為確定的已知函數的情形,若Vs(t)=∞,則 k(t)= 0,式(8)變?yōu)?/p>

當 sin(β(t)– β0)不恒為 0,即觀測站靜止且不位于目標航向線上時,可觀測的參數為(V0/ D0, C0, β0)而非(V0, C0, D0, β0)。當觀測站機動時,基于式(7)(k(t)= 0),(V0, C0, D0, β0)是可觀測的[4]。

設靜止的觀測站不位于目標航向線上,Vs(t)已知且Vs(t)∈ [V0, ∞ ),下面證明運動要素是可觀測的,且只需 4個無誤差方位即可唯一確定(V0, C0, D0,β0)。記 Z = (z1, …, z4)T,其中

則Z ≠ 0,且Z與目標運動要素相互唯一確定。令

利用這些記號,設 0 ≤ t1＜ t2＜ t3＜ t4,式(8)對應的方程組為

因為Z ≠ 0,不妨假設z4≠ 0,上述方程組依次相鄰兩式相減可得方程組

其中(以r 表示矩陣的行標,c 表示列標)

在觀測站不位于目標航向線上的條件下,矩陣G一般非奇異,因此

其中 p = (p1, p2, p3)T= G–1q,將式(10)代入式(9)的第一個方程,可得

因為 z4≠ 0,所以

代入式(10)計算(z1, z2, z3),由此可得方程式(9)的唯一非零解,這個解對應唯一的目標運動要素:

下面用一個算例質疑以上結論,設V0= 15.4 kn,C0= 65.3°,D0= 100.2鏈,β0= 0.1°,A(t)= 0,(t1, t2, t3, t4)= (5, 6.6, 7.1, 9.9)s,Vs(t1, t2, t3, t4)= (1455, 1495, 1505,1499 )m/s,則以上計算項目為:

計算的 V0, C0, D0, β0與設定值誤差均小于 10–7。

這些結論說明,考慮Vs(t)為確定的已知有限函數對應的純方位條件下的目標運動要素的可觀測性優(yōu)于認為Vs(t)= ∞ 的情形。因此,可以期望相應的觀測度應有同樣的性質,即考慮觀測誤差條件下,對應于相同的觀測站機動方式,Vs(t)＜ ∞ 對應的目標運動要素估計應優(yōu)于認為Vs(t)= ∞ 的情形。

2 假設平均聲速為常數的目標運動要素估計

平均聲速 Vs(t)一般不為常數,只有在一定區(qū)域內可以看成常數,下面假設Vs(t)= Vs為一常數,對于Vs已知的情形 θ = (V0, C0, D0, β0)為純方位觀測條件下的待估計參數向量,其取值空間可定為Θ = [1,25]×[0, 2π )×[1852, 37040)×[0, 2 π )。對于 Vs未知的情形,待估計參數向量為θ = (V0, C0, D0, β0, Vs),其取值空間可定為Θ =[1, 25]×[0, 2π )×[1852, 37040)×[0, 2 π )× [1400, 1550]。

設 tn = (n – 1)Δt時刻β(tn )的觀測值為βn,觀測誤 差 為 Δβ(tn),則 :Δβ(tn)=βn–β(tn),n=1,…,N 假 設Δβ(tn),n =1,…,N 為相互獨立同 N(0,σβ2)分布的隨機誤差,其平方和為

對于給定的Δt ＞ 0,N,A(tn), n = 1, …, N,若 Vs已知,則 J為(V0, C0, D0, β0)的函數,若 Vs未知,則 J 為(V0,C0, D0, β0,Vs)的函數,下面只討論這一情形。稱J在Θ上的最小值點 ( V0*, C0*, D0*, β0*, Vs*)為( V0, C0, D0,β0, Vs)的非線性最小二乘法估計,這個估計應為方程組的解,下面研究這個估計的性質。

根據以上假設,(β1, …, βN)的似然函數為

這個方程組等價于

所以,( V0*, C0*, D0*, β0*, Vs*)為的極大似然估計,且σβ2的極大似然估計為

(σβ2)*并非σ2β的無偏估計,類似于標準的線性最小二乘法的性質,應將式(14)內的 N 改為N – 5,若Vs已知,則改為N – 4,數值實驗發(fā)現,這個方法可以對σβ給出比較精確的估計。因為J( V0*, C0*, D0*, β0*,Vs*)為對應于估計值的剩余平方和,所以可以用σβ*檢驗計算的精度,σβ*與 σβ越接近,則觀測值信息利用得越充分。

因為非線性最小二乘法估計的算法不如線性最小二乘法簡單,除考慮計算的正確性、計算速度外,還必須考慮初值的穩(wěn)定性。數值實驗發(fā)現,計算 V0*,C0*, D0*, β0*, Vs*所需要的算法不太依賴于初值的選取。以區(qū)域 Θ 的中點作為搜索起點,可以給出比較精確的數值解。若有比較準確的目標運動要素的先驗知識,則人工輸入初值,可以加快計算速度。

假設已給出( V0*, C0*, D0*, β0*),則可以對任一時刻t目標實際距離D0(t),方位β0(t)進行估計,其公式為

3 運動要素估計的經典方法

若認為 Vs= ∞,則 k(t)= 0,式(7)變?yōu)?/p>

記 X = (tnsinβn, – tncosβn, sinβn, – cosβn)N×4Y =(xA(tn)sinβn– yA(tn)cosβn)N×1,則無誤差條件下 Y = XZ,在有觀測誤差的條件下,Y – XZ 的各分量并非相互獨立同分布的,且X 與方位觀測值有關,所以Y 與Z并非滿足嚴格線性擬合模型,但在觀測站機動的條件下,還是可以計算 min || Y – X Z ||2的解

Z*對應的目標運動要素估計用(V0l, C0l, D0l, β0l)表示,以下稱之為線性估計,這就是不考慮觀測信號時間延遲條件下,估計目標運動要素的傳統方法。即使觀測站機動,XTX 亦可能表現出近似奇異性,實驗程序內用偽逆 pinv( XTX )代替( XTX )–1。計算(V0l, C0l,D0l, β0l)的作用之一是作為非線性算法(V0, C0, D0,β0)的初值,數值實驗發(fā)現,這個初值對于加快計算速度有時是十分有效的。

4 給定觀測站機動方式的目標運動要素估計實驗

取 Vs(t)=1450m/s,σβ=0.5°,Δt=2s,V0=15.4kn,C0=65.3°,D0=100.2 鏈,β0=10.1°,進行數值實驗。

首先對觀測站等速直航情形,設觀測站以速率U = 4kn 沿航向 W = 50°等速直航,則 A(t)= U t eiW,這時傳統方法不能給出(V0, C0, D0, β0)的估計。取 Θ的中點作為搜索起點, 上述非線性方法計算的直到T = 600s的最近10個估計見表1。

表1 觀測站不機動的估計實驗

自表 1可以看出,觀測站等速直航情形,非線性方法可以對σβ給出比較精確的估計,對( V0, C0, D0,β0, Vs)的估計則比較差,下面對觀測站機動情形進行實驗。

假設觀測站機動方式為:初始時刻沿觀測方位對應的航向航行Tl= 120s,之后向方位增加的反方向轉向ΔW = 50°,假設轉向運動為勻速圓周運動,角速度為ω = 0.5° /s,轉向運動完成后直航Tl時間,之后向方位增加的方向轉向ΔW,…,以此類推進行周期為 Tl+ ΔW/ω的機動。其他參數同上,相同項目的一次實驗見表2。

表2 觀測站機動的估計實驗

可以看出,觀測站機動給出的目標運動要素估計明顯優(yōu)于觀測站不機動的情形。

5 結論

對于觀測站不機動目標等速直航的情形,認為Vs= ∞ 將導致對目標方位的觀測含有系統誤差,單就這個誤差看,似乎可以忽略,但所描述的系統的性質根本不同,簡單地說,此時若將目標速率與初始距離乘以相同的倍數,則對模型無任何影響,這與實際不符,兩種方法的簡單比較見表3。

由式(5)(6)可以看出,聲速有限決定了可觀測的目標方位中包含獨立的目標距離信息,認為 Vs= ∞等于放棄了這些信息,除非出于算法簡單的目的外,沒有理由認為Vs= ∞。

表3 非線性方法與線性方法的簡單比較

數值實驗發(fā)現,線性方法與非線性方法均能給出比較精確的β0的估計,關于V0、C0、D0,線性方法給出的估計穩(wěn)定很慢且不能穩(wěn)定于真值附近,非線性方法所給出的估計可以穩(wěn)定于真值附近。

非線性方法能給出σβ、Vs比較準確的估計,且σβ的估計能很快穩(wěn)定,也就是說,即使不是出于攻擊的目的,在對等速直航目標的純方位跟蹤過程中,可以對這兩個參數進行估計。

數值實驗發(fā)現,Vs(t)= Vs已知的情形,目標運動要素的可觀測度明顯優(yōu)于Vs未知的情形。

使用平均聲速這一概念可使上面的論述比較簡單,然而所得結論的可用性有待于實踐檢驗。

以方位觀測誤差平方和作為目標函數可以給出性質較優(yōu)的運動要素的估計,這個目標函數的性質有待于進一步研究。純方位方法目的是為了跟蹤目標當前態(tài)勢,尋求有利于這一目的的目標函數有待于進一步研究。

以上仿真數值實驗中搜索J 的最小值點使用的是Matlab函數lsqnonlin,算法是通用的,雖然較小的σβ有利于提高觀測度,但這個算法對置σβ= 0不能立即給出精確的估計,因此就上述非線性方法開發(fā)有針對性的快速準確算法是必要的。

[1]Daniel H. Wabgner,W. Charles Mylander, Thomas J.Sanders.海軍運籌分析[M](第3版).姜青山,鄭保華譯.,國防工業(yè)出版社,2008.

[2]董志榮.純方位系統定位與跟蹤的本載體最優(yōu)軌線方程及其最優(yōu)軌線[J].指揮控制與仿真,2007.29(1):7-15.

[3]潘志堅,閻福旺,劉孟奄,等.純方位水下目標運動分析方法研究[J].聲學學報,1997,22(1):87-92.

[4]石章松,劉忠. 單站純方位目標跟蹤系統可觀測性分析[J].火力與指揮控制, 2007,32(2):26-29.

[5]Robert J. Urick.水聲原理[M].第 3 版.哈爾冰:哈爾濱船舶工程學院出版社,1990.

[6]Frieden David R.海軍武器系統原理[M].第9版.上海翻譯出版公司, 1991.