序貫貝葉斯濾波估計(jì)模函數(shù)和水平波數(shù)

程澤豐，潘 翔

(浙江大學(xué)信息與電子工程學(xué)院，浙江 杭州 310027)

序貫貝葉斯濾波估計(jì)模函數(shù)和水平波數(shù)

程澤豐，潘 翔

(浙江大學(xué)信息與電子工程學(xué)院，浙江 杭州 310027)

介紹了一種在淺海波導(dǎo)環(huán)境中，聯(lián)合聲場(chǎng)模型和測(cè)量場(chǎng)，利用?；盘?hào)處理器解決非線性狀態(tài)-參數(shù)估計(jì)問題的方法.在分析簡(jiǎn)正波模型表征聲場(chǎng)的基礎(chǔ)上，側(cè)重研究了擴(kuò)展卡爾曼濾波和粒子濾波兩種算法，并將其用于?；曅盘?hào)處理進(jìn)行聲場(chǎng)參數(shù)估計(jì).理論分析和數(shù)值仿真表明，兩種算法均能正確估計(jì)模函數(shù)和水平波數(shù)，同時(shí)粒子濾波處理具有更好的寬容性.

序貫貝葉斯；粒子濾波；模基處理；水聲波導(dǎo)

0 引 言

淺海環(huán)境是一個(gè)時(shí)變和空變的波導(dǎo)環(huán)境，水聲信號(hào)在其傳播表現(xiàn)為時(shí)延擴(kuò)展和多普勒擴(kuò)展.在波導(dǎo)環(huán)境中，通常用簡(jiǎn)正波模型來表征聲場(chǎng)，用匹配場(chǎng)處理進(jìn)行聲源定位[1].由于理論模型和實(shí)際聲場(chǎng)存在失配，寬容的匹配場(chǎng)處理成為研究的熱點(diǎn).

基于空間-狀態(tài)模型的序貫貝葉斯濾波[2]適合在不確定的、非平穩(wěn)和非高斯的變化環(huán)境中進(jìn)行信號(hào)分析和處理，可以為模函數(shù)和水平波數(shù)的跟蹤以及環(huán)境參數(shù)的反演提供有力工具.本文依據(jù)實(shí)際測(cè)量的聲壓場(chǎng)和聲速剖面，基于簡(jiǎn)正模式傳播模型，通過序貫貝葉斯濾波算法發(fā)展而來的基于擴(kuò)展卡爾曼濾波[3](Extended Kalman Filtering, EKF)和粒子濾波(Particle Filtering, PF)的模基處理器，估計(jì)對(duì)應(yīng)的模函數(shù)和水平波數(shù).

1 空間-狀態(tài)模型

運(yùn)用序貫貝葉斯濾波處理海洋聲信號(hào)增強(qiáng)問題，首先需要建立空間-狀態(tài)模型.假設(shè)水平分層的海洋深度為h，已知聲源水平距離rs和深度zs.由點(diǎn)聲源發(fā)出的聲能量用亥姆霍茲方程建模，經(jīng)過標(biāo)準(zhǔn)的變量分離方法，獲得只與深度相關(guān)的波方程表示方法

(1)

式中，本征解{Φm(z)}即所謂的模函數(shù)，κz為z方向上的垂直波數(shù).同時(shí)對(duì)應(yīng)的色散關(guān)系為

(2)

將只含深度參量的本征方程轉(zhuǎn)化為狀態(tài)-空間形式，將第m模式定義為

(3)

將噪聲用加性噪聲表示，并將確定的狀態(tài)方程擴(kuò)展為M個(gè)模式Φ(z)=[Φ1(z)|…|ΦM(z)]T以此導(dǎo)出模型的2M維的高斯-馬爾可夫表示下的狀態(tài)和測(cè)量方程

p(rs,z)=CT(rs,zs)Φ(z)+v(z),

(4)

由于陣列在空間上以深度離散采樣聲壓場(chǎng)，可以通過離散化微分方程來提高數(shù)值穩(wěn)定性.假設(shè)傳感器均勻分布，以zl表示第l個(gè)傳感器所在的位置，得到

(5)

其中Δzl=zl+1-zl，將該近似方法代入式(1)，得到

Φm(zl)-2Φm(zl-1)+Φm(zl-2)+Δzl2κz2(m)Φm(zl-1)=0

將離散模態(tài)狀態(tài)向量定義為Φm(zl)=[Φm(zl-2)|Φm(zl-1)]T，得到關(guān)于第m模式的一組差分方程.最后結(jié)合所有可能的模式，對(duì)簡(jiǎn)正模狀態(tài)方程和測(cè)量方程的整體高斯-馬爾可夫表示為

Φ(zl)=A(zl)Φ(zl-1)+ω(zl),

p(rs,zl)=CT(rs,zs)Φ(zl)+v(zl).

(6)

2 ?；幚砥?/h2>

2.1 擴(kuò)展卡爾曼濾波(EKF)

利用有限差分方法來離散化狀態(tài)方程，導(dǎo)出第一節(jié)中所建立模型的離散形式

Φ(zl+1)=[I+ΔzA(zl)]Φ(zl)+Δzω(zl).

(7)

確立了高斯-馬爾可夫模型，對(duì)狀態(tài)向量Φ確定其遞歸(卡爾曼)估計(jì)器

(8)

其中，AΦ(zl)=[I+ΔzA(zl)]，這種一般化表示方法，能夠反映整個(gè)過程的遞歸本質(zhì).

(9)

增益：卡爾曼增益由下式給出

(10)

因此，在經(jīng)過增強(qiáng)的高斯-馬爾可夫模型中，非線性的測(cè)量方程可以寫作

(11)

其中的c[·]的并非聲速剖面c(z).需要得到增強(qiáng)的狀態(tài)向量和相對(duì)應(yīng)的測(cè)量中所有元素的一階導(dǎo)數(shù)，構(gòu)成整個(gè)系統(tǒng)的雅克比矩陣.

2.2 粒子濾波(PF)

擴(kuò)展的卡爾曼濾波器基于泰勒級(jí)數(shù)的一階展開，將非線性問題近似為線性問題.相比之下，粒子濾波則更加適合處理非線性問題.

(12)

(13)

對(duì)應(yīng)的狀態(tài)和測(cè)量方程分別為

Φi(zl,Θ)=A(zl-1)Φi(zl-1)+ωi(zl-1),

p(rs,zl)=c[Φi(zl-1;Θ)]+v(zl).

(14)

3 仿真分析

3.1 仿真參數(shù)設(shè)置

圖1詳細(xì)給出了仿真環(huán)境參數(shù)的設(shè)定.圖1左側(cè)串狀圓圈為接收陣，陣列的布放從水面以下深度10 m處開始，間隔Δzl=2.5 m設(shè)置陣元，一共有N=71個(gè)陣元，最后一個(gè)陣元放置在深度為185 m處的位置.圖1中右側(cè)標(biāo)志為聲源所在位置，離接收陣水平距離rs為0.5 km，深度zs為90 m，與垂直陣列中心位置等深，聲源發(fā)出聲信號(hào)的頻率為50 Hz.將圖1所示環(huán)境參數(shù)輸入到簡(jiǎn)正模計(jì)算軟件Kraken[4]中，所需要的垂直聲速剖面如圖2左圖所示，計(jì)算獲得的聲壓場(chǎng)如圖2右圖所示.

圖1 仿真參數(shù)設(shè)置

圖2 垂直聲速剖面和聲壓場(chǎng)

圖3給出了Kraken計(jì)算的結(jié)果，從左至右分別為1至8號(hào)模.

圖3 Kraken計(jì)算得到的1至8號(hào)模

3.2 結(jié)果分析

通過Kraken計(jì)算，得到8個(gè)模式.圖4給出分別利用EKF和PF?；幚砥鲗?duì)第5號(hào)模模函數(shù)和水平波數(shù)的估計(jì)結(jié)果.圖4(c)、圖4(d)橫軸表示水平波數(shù)，單位為m-1.

(a) EKF?；幚砥鞴烙?jì)模函數(shù)； (b)PF?；幚砥鞴烙?jì)模函數(shù)；(c) EKF?；幚砥鞴烙?jì)水平波數(shù)； (d)PF?；幚砥鞴烙?jì)水平波數(shù)圖4 EKF和PF?；幚砥鞴烙?jì)結(jié)果

圖4中，實(shí)線表示為真值，圖4(a)中虛線表示EKF?；幚砥鲗?duì)模函數(shù)的估計(jì)，圖4(b)中虛線表示PF對(duì)模函數(shù)的估計(jì)結(jié)果.圖4(c)和(d)中“+”號(hào)線分別為EKF和PF?；幚砥鲗?duì)第5號(hào)模對(duì)應(yīng)的水平波數(shù)的估計(jì)值.對(duì)第5號(hào)模而言，基于PF的?；幚砥鲗?duì)模函數(shù)和水平波數(shù)的估計(jì)結(jié)果均優(yōu)于基于EKF模基處理器，且MMSE在10-4數(shù)量級(jí).分析其他模式，也存在類似的結(jié)論.表明了基于PF的?；幚砥髁己玫膶捜菪?

4 結(jié)束語

在獲得聲速剖面和聲壓場(chǎng)數(shù)據(jù)的前提下，本文研究了表征聲場(chǎng)的模函數(shù)和水平波數(shù)的估計(jì)問題.利用擴(kuò)展的卡爾曼濾波器(EKF)和粒子濾波器(PF)構(gòu)建了?；幚砥鳎Y(jié)合測(cè)量聲壓場(chǎng)和聲速剖面聯(lián)合估計(jì)模函數(shù)和水平波數(shù).結(jié)果表明EKF算法在該問題上存在一定的局限性.而PF?；幚砥?，在設(shè)定的條件下，能夠較為精確地估計(jì)模函數(shù)和水平波數(shù)，但其算法實(shí)時(shí)性遠(yuǎn)劣于前者.在今后的研究中，需要進(jìn)一步檢驗(yàn)PF模基處理器的寬容性.

[1]RISTICB,ARULAMPALAMS,GORDONN.BeyondtheKalmanFilter:ParticleFiltersforTrackingApplications[J].IEEETransofAerospace&ElectronicSystems, 2004, 19(7): 37-38.

[2]CANDYJV.BayesianSignalProcessing:Classical,ModernandParticleFilteringMethods[M].Hoboken,NewJersey, 2009: 369-417.

[3]CANDYJV,SULLIVANEJ.Model-basedIdentification:AnAdaptiveApproachtoOcean-acousticProcessing[J].IEEEJournalofOceanicEngineering, 1996, 21(3): 273-289.

[4]PORTERMB.TheKRAKENNormalModeProgram[R].SACLANTUnderseaResearchCentre,May17, 2001: 1-196.

[5]CANDYJV.EnvironmentallyAdaptiveProcessingforShallowOceanApplications:ASequentialBayesianApproach[J].JournaloftheAcousticalSocietyofAmerican, 2015, 138(3): 1268-1281.

Modal Function and Horizontal Wavenumber Estimation Using Sequential Bayesian Filters

CHENG Zefeng, PAN Xiang

(CollegeofInformationScience&ElectronicEngineering,ZhejiangUniversity,HangzhouZhejiang310027,China)

Combining acoustic pressure-field model and measurements based on model-based signal processing in shallow water waveguide, this paper introduces a solution to nonlinear state-parameter estimation problem. It focuses on the two algorithms, the extended Kalman filters and the particle filters, to estimate the parameters of acoustic pressure-field after analyzing the normal mode model representations of acoustic pressure-field. The theoretical analyses and the simulation results both show that the two algorithms are capable of estimating the modal functions and horizontal wavenumbers correctly. Meanwhile, the particle filter algorithm has a better robustness.

sequential Bayesian estimation; particle filter; model-based processing; underwater acoustic waveguide

10.13954/j.cnki.hdu.2017.01.011

2016-09-14

國家自然科學(xué)基金資助項(xiàng)目(61571397)

程澤豐(1993-)，男，浙江舟山人，碩士研究生，電子與通信工程.通信作者：潘翔副教授，E-mail: panxiang@zju.edu.cn.

TN929.3

A

1001-9146(2017)01-0046-05