李長文,趙建昕

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

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

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

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

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

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

圖1 無誤差的觀測站A、目標(biāo)B、可觀測點(diǎn) P 的位置圖

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

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

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

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

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

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

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

因此

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

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

因?yàn)?D(t)≥ 0,所以

由式(4)可得

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

將式(4)改寫為

此式兩端的實(shí)部與虛部之比對應(yīng)的等式為

利用這一關(guān)系可以研究目標(biāo)運(yùn)動要素的可觀測性。

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

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

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

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

則Z ≠ 0,且Z與目標(biāo)運(yùn)動要素相互唯一確定。令

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

因?yàn)閆 ≠ 0,不妨假設(shè)z4≠ 0,上述方程組依次相鄰兩式相減可得方程組

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

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

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

因?yàn)?z4≠ 0,所以

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

下面用一個算例質(zhì)疑以上結(jié)論,設(shè)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,則以上計算項(xiàng)目為:

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

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

2 假設(shè)平均聲速為常數(shù)的目標(biāo)運(yùn)動要素估計

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

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

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

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

這個方程組等價于

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

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

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

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

3 運(yùn)動要素估計的經(jīng)典方法

若認(rèn)為 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 的各分量并非相互獨(dú)立同分布的,且X 與方位觀測值有關(guān),所以Y 與Z并非滿足嚴(yán)格線性擬合模型,但在觀測站機(jī)動的條件下,還是可以計算 min || Y – X Z ||2的解

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

4 給定觀測站機(jī)動方式的目標(biāo)運(yùn)動要素估計實(shí)驗(yàn)

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

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

表1 觀測站不機(jī)動的估計實(shí)驗(yàn)

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

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

表2 觀測站機(jī)動的估計實(shí)驗(yàn)

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

5 結(jié)論

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

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

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

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

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

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

使用平均聲速這一概念可使上面的論述比較簡單,然而所得結(jié)論的可用性有待于實(shí)踐檢驗(yàn)。

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

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

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

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

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

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

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

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