鞠陽(yáng)， 元勇虎， 林蔚， 徐海生

(1.中國(guó)船舶工業(yè)系統(tǒng)工程研究院, 北京 100094; 2.哈爾濱工程大學(xué) 理學(xué)院, 黑龍江 哈爾濱 150001)

0 引言

對(duì)海洋進(jìn)行探測(cè)及研究對(duì)維護(hù)海洋和建設(shè)海洋有著極其重要的意義。借助于聲探測(cè)器，研究者可以得到運(yùn)動(dòng)目標(biāo)的某些聲信息，目標(biāo)運(yùn)動(dòng)分析(TMA)就是利用傳感器測(cè)得、受噪聲干擾的聲信息來(lái)估計(jì)運(yùn)動(dòng)目標(biāo)的距離、方位及運(yùn)動(dòng)速度等運(yùn)動(dòng)參數(shù)。目標(biāo)運(yùn)動(dòng)參數(shù)解算方法主要有純方位TMA方法[1-2]、匹配場(chǎng)處理技術(shù)[3-4]、多要素聯(lián)合方位的TMA[5-6]等。純方位TMA解算方法僅使用目標(biāo)的方位角通過(guò)運(yùn)動(dòng)模型的幾何關(guān)系建立數(shù)學(xué)模型，由卡爾曼濾波、粒子濾波等算法對(duì)目標(biāo)的運(yùn)動(dòng)參數(shù)進(jìn)行估計(jì)與跟蹤。匹配場(chǎng)處理技術(shù)通過(guò)波導(dǎo)原理，根據(jù)已知的環(huán)境參數(shù)信息與聲強(qiáng)信息求得聲場(chǎng)模型，進(jìn)而解算聲源目標(biāo)的參數(shù)信息。多要素聯(lián)合方位的TMA方法使用方位及頻率、到達(dá)時(shí)間差、倒譜信息等，通過(guò)卡爾曼濾波、粒子濾波等技術(shù)對(duì)目標(biāo)進(jìn)行參數(shù)估計(jì)。在實(shí)際應(yīng)用中，傳統(tǒng)的純方位TMA解算方法通常需要大約20 min左右才能使TMA運(yùn)算收斂，得出目標(biāo)運(yùn)動(dòng)參數(shù)估計(jì)結(jié)果。而匹配場(chǎng)處理技術(shù)不僅需要環(huán)境參數(shù)等先驗(yàn)知識(shí)，還對(duì)水聽(tīng)器陣列等硬件要求較高。由于目標(biāo)發(fā)生機(jī)動(dòng)過(guò)程一般在較短時(shí)間內(nèi)完成，解算時(shí)間過(guò)長(zhǎng)很難對(duì)目標(biāo)機(jī)動(dòng)做出快速反應(yīng)。本文給出基于勻加速直線運(yùn)動(dòng)模型的多普勒頻移數(shù)學(xué)模型與參數(shù)搜索優(yōu)化方法，提出一種基于多普勒頻率信息的被動(dòng)快速TMA解算方法，能夠較快地估計(jì)目標(biāo)運(yùn)動(dòng)參數(shù)信息。

1 勻加速直線運(yùn)動(dòng)模型

移動(dòng)目標(biāo)的TMA問(wèn)題可以轉(zhuǎn)化為多個(gè)運(yùn)動(dòng)參數(shù)的估計(jì)問(wèn)題，多普勒頻移信息中蘊(yùn)含著目標(biāo)的頻率、速度等參數(shù)信息，通過(guò)估計(jì)多普勒頻移可獲得目標(biāo)的運(yùn)動(dòng)參數(shù)信息。

目前，基于多普勒頻移的TMA解算方法中，均假設(shè)目標(biāo)進(jìn)行勻速直線運(yùn)動(dòng)。而在實(shí)際中，目標(biāo)存在著加速的運(yùn)動(dòng)狀態(tài)。

假設(shè)觀察方W靜止不動(dòng)，目標(biāo)M沿著一個(gè)方向做勻加速直線運(yùn)動(dòng)，M0為目標(biāo)初始位置，如圖1所示。

設(shè)目標(biāo)M的發(fā)射頻率為f0，初始速度為v0，加速度為a，目標(biāo)與觀測(cè)方的當(dāng)前方位角為θ，聲在媒介中的傳播速度為c，靜止的被動(dòng)接收器與目標(biāo)運(yùn)動(dòng)方向的垂直距離為l，目標(biāo)M到達(dá)最短距離的時(shí)間為tc，則根據(jù)多普勒頻移[7]，第t時(shí)刻接收到的瞬時(shí)頻率為

(1)

式中：vr為目標(biāo)M在其與接收點(diǎn)方向上的速度。

根據(jù)運(yùn)動(dòng)模型的幾何關(guān)系，目標(biāo)M在其與接收點(diǎn)方向上的速度為

vr(t)=v(t)cosθ,

(2)

式中：

v(t)=v0+at；

則基于勻加速直線運(yùn)動(dòng)模型的多普勒頻移表達(dá)式為

f(f0,v0,a,tc,l;t)=

(3)

若加入?yún)?shù)方位角θ，則表達(dá)式為

(4)

運(yùn)動(dòng)目標(biāo)的TMA問(wèn)題即為多個(gè)運(yùn)動(dòng)參數(shù)的估計(jì)問(wèn)題。若瞬時(shí)頻率為已知條件，則對(duì)表達(dá)式進(jìn)行轉(zhuǎn)化，參數(shù)的估計(jì)問(wèn)題就可以轉(zhuǎn)化為求解函數(shù)最小化問(wèn)題。則兩類運(yùn)動(dòng)模型的最小化問(wèn)題為

(5)

f2=[f0,v0,a]=

(6)

式中：f(t)為第t時(shí)刻的瞬時(shí)頻率；T為采樣頻率；f1和f2表示上述兩種不同的表達(dá)式。

上述模型中，瞬時(shí)頻率f(t)、方位角θ及時(shí)間t是可測(cè)量的，在本文作為已知條件。瞬時(shí)頻率f(t)通常由接收器接收噪聲信號(hào)通過(guò)時(shí)頻變換方法得到，如快速傅里葉變換算法、高分辨譜算法等，本文不做詳細(xì)敘述。

2 快速TMA解算方法

解算TMA最小值目標(biāo)函數(shù)求解各項(xiàng)參數(shù)的過(guò)程中，由于該模型是非線性的，待求參數(shù)較多，且距離較遠(yuǎn)時(shí)角度變化、頻率變化較小，而其值對(duì)角度與頻率較為敏感，使用傳統(tǒng)方法無(wú)法使求解的參數(shù)結(jié)果收斂。因此，本文提出一種聯(lián)合雙模型與參數(shù)搜索優(yōu)化方法的TMA解算方法，使得TMA參數(shù)解算結(jié)果收斂。

2.1 內(nèi)點(diǎn)- Levenberg-Marquard方法

非線性內(nèi)點(diǎn)法[8]具有對(duì)求解問(wèn)題的規(guī)模不敏感、計(jì)算大規(guī)模非線性規(guī)劃問(wèn)題的時(shí)間復(fù)雜度低、尋優(yōu)速度快、魯棒性強(qiáng)等優(yōu)點(diǎn)。其主要思想為在可行域的邊界筑起一道很高的“圍墻”，當(dāng)?shù)c(diǎn)靠近邊界時(shí)，目標(biāo)函數(shù)徒然增大，以示懲戒，防止迭代點(diǎn)穿越已設(shè)邊界，使得最優(yōu)解收斂于可行域之內(nèi)。

對(duì)于下列不等式約束問(wèn)題：

(7)

式中：n為變量維度；ri(x)為第i個(gè)約束條件；k為總的約束個(gè)數(shù)。

構(gòu)建一個(gè)懲罰函數(shù)來(lái)代替約束條件中的不等式約束。即(7)式變?yōu)?/p>

(8)

式中：I(u)為懲罰函數(shù)，未違反約束時(shí)其值趨于0，違反約束時(shí)其值趨于∞，使得最優(yōu)值收斂于可行域內(nèi)。設(shè)懲罰函數(shù)為

(9)

式中：p為懲罰因子,p>0. 將懲罰函數(shù)代入(8)式中，得到近似優(yōu)化問(wèn)題：

(10)

對(duì)于轉(zhuǎn)化得到的近似優(yōu)化問(wèn)題，本文使用高斯- 牛頓方法的 Levenberg-Marquard(LM)方法[9]進(jìn)行求解。LM方法能借由執(zhí)行時(shí)修改參數(shù)達(dá)到結(jié)合高斯- 牛頓算法以及梯度下降法的優(yōu)點(diǎn)，并對(duì)二者之不足作改善，每次迭代尋找一個(gè)合適的阻尼因子參數(shù)λ. 當(dāng)λ很小時(shí)，算法就變成了高斯- 牛頓法的最優(yōu)步長(zhǎng)計(jì)算式；當(dāng)λ很大時(shí)，蛻化為梯度下降法的最優(yōu)步長(zhǎng)計(jì)算式。

內(nèi)點(diǎn)- LM方法的具體步驟為：

步驟1設(shè)立懲罰因子p，懲罰參數(shù)φ，迭代次數(shù)k=1.

步驟2根據(jù)可行域與目標(biāo)函數(shù)，通過(guò)內(nèi)點(diǎn)法轉(zhuǎn)化得到近似優(yōu)化問(wèn)題，利用LM方法進(jìn)行求解。

步驟3檢測(cè)是否滿足收斂條件，滿足則跳出循環(huán)。

步驟4令p=φp，k=k+1，返回步驟2.

2.2 基于雙模型的參數(shù)搜索優(yōu)化方法

為使TMA解算的各參數(shù)收斂，本文提出基于雙模型的參數(shù)搜索優(yōu)化方法，具體步驟可以表述為以下6個(gè)階段：

1) 設(shè)立各個(gè)變量的初始值(f0(0),v0(0),a(0),tc(0),l(0))及各個(gè)變量的可行域，令TMA解算次數(shù)i=1.

2)由估計(jì)值(f0(i-1),v0(i-1),a(i-1))與各可行域，通過(guò)內(nèi)點(diǎn)- LM方法對(duì)f2進(jìn)行求解，得到估計(jì)值(f0(i-1,1),v0(i-1,1),a(i-1,1))。

3)將f0,v0,a的搜索范圍分別設(shè)為[f0(i-1,1)-0.5,f0(i-1,1)+0.5],[j-1,j+1],[0.9a(i-1),1.1a(i-1)]，其中j∈[1,max(v0)-2]。即將發(fā)射頻率與加速度進(jìn)行小幅度搜索，速度v0的搜索空間以2單位為搜索范圍，間隔為1單位，共并行搜索計(jì)算(max(v0)-2)次，選擇其中使目標(biāo)函數(shù)最小值的對(duì)應(yīng)參數(shù)作為新的估計(jì)值(f0(i-1,2),v0(i-1,2),a(i-1,2))。

4)固定變量(f0,v0,a)為(f0(i-1,2),v0(i-1,2),a(i-1,2))，由估計(jì)值(tc(i-1),l(i-1))與各變量的可行域，通過(guò)內(nèi)點(diǎn)- LM方法對(duì)f1進(jìn)行4次求解，得到(tc,l)的估計(jì)值(tc(i-1,1),l(i-1,1))。

5)得到第i次解算的運(yùn)動(dòng)參數(shù)近似值(f0(i),v0(i),a(i),tc(i),l(i))，其值為(f0(i-1,2),v0(i-1,2),a(i-1,2),tc(i-1,1),l(i-1,1))。

6)將第i次解算的近似值(f0(i),v0(i),a(i),tc(i),l(i))作為第i+1次的初始估計(jì)值，返回步驟2. 若求解的近似值收斂于穩(wěn)定值，則判定為收斂。

3 數(shù)值仿真與分析

設(shè)觀察方W靜止不動(dòng)，運(yùn)動(dòng)目標(biāo)M的發(fā)射頻率設(shè)為f0=135 Hz，初始速度設(shè)為v0=4 m/s，加速度設(shè)為a=0.2 m/s2，正橫距離設(shè)為l=4 000 m，正橫時(shí)間設(shè)為tc=250 s. 仿真模擬時(shí)間共為400 s，每秒采樣5次瞬時(shí)頻率信息，共采樣10 s，每間隔5 s進(jìn)行1次TMA解算，仿真實(shí)驗(yàn)設(shè)定如表1所示，仿真環(huán)境為Windows 10 64位操作系統(tǒng)，Mathematica科學(xué)計(jì)軟件。

運(yùn)用本文所提聯(lián)合雙模型與參數(shù)搜索優(yōu)化方法的TMA解算方法，進(jìn)行參數(shù)估計(jì)解算，與傳統(tǒng)方法——牛頓法進(jìn)行TMA解算的估計(jì)結(jié)果進(jìn)行比較，得到發(fā)射頻率f0、速度v0、加速度a、正橫時(shí)間tc、正橫距離l及距離r的估計(jì)值變化曲線，如圖2～圖7所示。由圖2～圖7可看出：發(fā)射頻率、速度、加速度在整個(gè)解算過(guò)程中趨于穩(wěn)定，誤差均小于5%；而正橫時(shí)間在30 s后趨于穩(wěn)定，收斂后誤差小于5%；正橫距離與距離則是20 s后達(dá)到收斂，170 s后趨于穩(wěn)定，穩(wěn)定后誤差小于3%. 在仿真環(huán)境下，每次本文給出的優(yōu)化方法進(jìn)行TMA解算用時(shí)為7 s，解算速度較快。仿真結(jié)果表明，由于解算模型是非線性的，待求參數(shù)較多，且距離較遠(yuǎn)時(shí)角度變化、頻率變化較小，而其值對(duì)角度與頻率較為敏感，傳統(tǒng)方法無(wú)法使求解的各個(gè)參數(shù)結(jié)果收斂。而本文給出的搜索方法與解算方法使結(jié)果精度與計(jì)算效率有了一定提高，使得各參數(shù)較快地收斂于允許誤差范圍內(nèi)，且求解的目標(biāo)參數(shù)信息估計(jì)值較為穩(wěn)定，在設(shè)定真實(shí)值附近作微弱波動(dòng)。當(dāng)水平距離較遠(yuǎn)時(shí)，由于角度變化值非常小，幾乎沒(méi)有變化，使得解算誤差變得較大，甚至不收斂，這是多普勒效應(yīng)的特點(diǎn)所致，即多普勒效應(yīng)適用于中、近距離的TMA解算。

表1 目標(biāo)運(yùn)動(dòng)參數(shù)

注: 表中收斂誤差大部分設(shè)定為10%.

圖2 發(fā)射頻率估計(jì)Fig.2 Transmitting frequency estimation

圖3 初始速度估計(jì)Fig.3 Initial speed estimation

圖4 加速度估計(jì)Fig.4 Acceleration estimation

圖5 正橫時(shí)間估計(jì)Fig.5 Transverse time estimation

圖6 正橫距離估計(jì)Fig.6 Transverse range estimation

圖7 距離估計(jì)Fig.7 Range estimation

4 結(jié)論

本文利用測(cè)量到的方位信息與多普勒頻移的線譜瞬時(shí)頻率信息，通過(guò)建立的勻加速直線運(yùn)動(dòng)的運(yùn)動(dòng)模型，給出了勻加速直線運(yùn)動(dòng)的多普勒頻移表達(dá)式，提出了一種基于內(nèi)點(diǎn)- LM方法的聯(lián)合雙模型與參數(shù)搜索優(yōu)化方法的TMA解算方法。該方法僅需要被動(dòng)聲納，且不需要接收方發(fā)生機(jī)動(dòng)，就能夠快速地獲得一定精度的目標(biāo)運(yùn)動(dòng)參數(shù)估計(jì)值，有效提高了安全性與隱蔽性。數(shù)值仿真結(jié)果表明：發(fā)射頻率、速度、加速度及正橫時(shí)間在本文方法下快速地達(dá)到收斂，正橫距離與距離也在20 s后快速地達(dá)到收斂，誤差均小于5%. 瞬時(shí)頻率的測(cè)量精度與方位角的測(cè)量精度對(duì)解算模型的值較為敏感，對(duì)運(yùn)動(dòng)參數(shù)估計(jì)的精確度影響較大，研究更加精確的瞬時(shí)頻率估計(jì)方法與方位角估計(jì)方法是TMA進(jìn)一步研究的方向之一。