鞠陽， 元勇虎， 林蔚， 徐海生

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

0 引言

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

1 勻加速直線運動模型

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

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

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

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

(1)

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

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

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

(2)

式中：

v(t)=v0+at；

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

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

(3)

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

(4)

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

(5)

f2=[f0,v0,a]=

(6)

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

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

2 快速TMA解算方法

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

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

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

對于下列不等式約束問題：

(7)

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

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

(8)

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

(9)

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

(10)

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

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

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

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

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

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

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

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

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

2)由估計值(f0(i-1),v0(i-1),a(i-1))與各可行域，通過內(nèi)點- LM方法對f2進(jìn)行求解，得到估計值(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單位，共并行搜索計算(max(v0)-2)次，選擇其中使目標(biāo)函數(shù)最小值的對應(yīng)參數(shù)作為新的估計值(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))，由估計值(tc(i-1),l(i-1))與各變量的可行域，通過內(nèi)點- LM方法對f1進(jìn)行4次求解，得到(tc,l)的估計值(tc(i-1,1),l(i-1,1))。

5)得到第i次解算的運動參數(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次的初始估計值，返回步驟2. 若求解的近似值收斂于穩(wěn)定值，則判定為收斂。

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

設(shè)觀察方W靜止不動，運動目標(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è)為tc=250 s. 仿真模擬時間共為400 s，每秒采樣5次瞬時頻率信息，共采樣10 s，每間隔5 s進(jìn)行1次TMA解算，仿真實驗設(shè)定如表1所示，仿真環(huán)境為Windows 10 64位操作系統(tǒng)，Mathematica科學(xué)計軟件。

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

表1 目標(biāo)運動參數(shù)

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

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

圖3 初始速度估計Fig.3 Initial speed estimation

圖4 加速度估計Fig.4 Acceleration estimation

圖5 正橫時間估計Fig.5 Transverse time estimation

圖6 正橫距離估計Fig.6 Transverse range estimation

圖7 距離估計Fig.7 Range estimation

4 結(jié)論

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