水下被动目标参数估计的群体蒙特卡罗方法

陈 钊 王海燕 刘郑国，2 申晓红

（1.西北工业大学航海学院电子与通信工程系，陕西 西安710072；2.中国船舶重工集团公司办公厅，北京100097）

引 言

精确定位远场中靠得很近的多个声源是水声信号处理中的关键问题之一。传统的波束形成法受到瑞利限的制约，分辨率很低。最大熵谱法和最小方差法是对常规波束形成法的简单修正，两者在估计性能上没有大的突破。子空间类方法［1］极大地突破了瑞利限的制约，使得高分辨目标方位估计方法向前迈进了一大步，但是在低信噪比，少快拍，信号源相关以及存在系统误差时［2］，此类方法的性能会受到很大的影响。最近十几年，贝叶斯高分辨方位估计方法受到普遍关注并被广泛研究。贝叶斯高分辨方位估计技术是利用信号和噪声的联合后验概率密度函数对信号进行谱估计的一种方法。Lasenby和Fitzgerald最早提出了这种方法［3］，它利用贝叶斯定理把被估计量看作是随机变量，并且引入被估计量的先验知识，从而改善了参数估计的性能。此后不同学者对贝叶斯高分辨方位估计方法进行了探索和改进，提出了许多新的算法［4－5］。

对于水下目标的检测，方位估计和跟踪问题，与主动方式相比，被动方式具有更好的隐蔽性。被动方式利用的是水下目标的辐射噪声，通常涉及宽带信号处理问题。传统的基于子空间类的高分辨方位估计方法通常需要较多的信号快拍数且须在频域内进行处理，涉及频域子带的划分。被广泛使用的相干信号子空间方法需要方位预估，预估的准确性严重影响了目标方位估计的精度。

利用水下被动目标的辐射噪声进行目标方位估计，采用一种时域宽带信号模型对阵列接收目标辐射噪声进行建模［6］。利用贝叶斯定理得到基于贝叶斯模型的目标参数的后验概率密度函数；对得到的后验概率密度函数进行峰值搜索，峰值对应的参数值即所需的目标方位估计值。传统的贝叶斯高分辨方法需要进行多维搜索，随着被估计参数维数的增加，计算量呈指数增长，许多文献提出使用可逆跳变马尔科夫链蒙特卡罗（RJMCMC）方法进行后验概率密度函数的峰值搜索［7－9］。但 RJMCMC算法需要一个样本被丢弃的“烧录”周期，每次迭代只产生一个样本，严重制约了其并行执行的能力。采用一种最新的群体蒙特卡罗（PMC）采样策略执行贝叶斯计算［10］，它在每次迭代时使用重要性采样从目标分布（即后验概率密度函数）中产生一组近似无关的样本，从而显著地增强了算法并行执行的能力。利用这些样本便可进行模型阶数和目标方位的联合估计。同时采用最大后验（MAP）声源恢复方法［6］可以对各个声源进行幅度恢复。

1 水下目标辐射噪声模型

水下目标辐射噪声的建模有许多方法［11－12］。采用与文献［12］类似的方法产生阵列接收宽带辐射噪声来仿真水下被动目标的辐射噪声，它是连续谱成分和线谱成分的叠加，如式（1）所示

式中：x0（t）表示连续谱噪声，由高斯白噪声通过一个低通滤波器来模拟；xr（t）表示多个线谱成分；R为线谱数；Ar、fr和φr分别为第r个线谱成份的幅度、频率和初始相位。

2 基于时域宽带信号模型的贝叶斯高分辨DOA估计方法

2.1 时域宽带信号模型概述

给出用于水下被动目标辐射噪声插值建模的时域宽带信号模型。假定有一M元均匀线列阵（ULA），K个声源sk（t）分别从θk，k＝0，…，K－1方向入射，声源sk（t）为有限带宽的宽带或窄带信号

式中：和分别为声源sk（t）的频率上限和频率下限；Δfk为其带宽。

当声源入射方位角φk≤π/2时，阵元1首先接收到信号，（M－1）τk时延之后阵元M才接收到信号，如图1所示。声源sk（t）在第m个阵元上产生的信号可以表示为

图1 入射角φk≤π/2时的阵列信号接收模型

式中：τk＝d＊sin（θk）/c为第k个声源信号的阵元间延迟，c为水中声速，d为阵元间距；hl（mτk）是对应于时延mτk的插值序列的第l个元素，它采用加权sinc函数，并且将汉明窗作为加权窗。以合适的采样频率fs（对应采样间隔Ts）对sk（t）进行离散化可得离散序列sk（nTs），且省写为sk（n）.所以式（3）可表示为

式中，D＝mτk/Ts为归一化时延。

假设噪声是服从均值为零，标准差为σw的高斯白噪声，且不同阵元和不同快拍上的噪声相互独立，即空域和时域均独立。由式（3）和式（4）可得t时刻阵元接收向量y（n）∈RM的表达为［6］

式中：Hl（τ）∈RM×K是第l个插值矩阵；τ＝［τ0，τ1，…，τK－1］T∈RK为与声源入射方位相对应的时延向量；是n次快拍时刻的声源幅度向量；w（n）∈RM是n次快拍时刻服从N（0，IM）分布的高斯白噪声。最终定义一个向量z（n）∈RM

它可以看作是快拍y（n）与其基于插值逼近之间的误差，可写作为

对于入射方位角φk＞π/2的情况，阵元M首先接收到入射信号，（M－1）τk时延之后阵元1才接收到信号，如图2所示。声源sk（t）在第m个阵元上产生的信号可以表示为

式中，EM为交换矩阵，它交换矩阵Hl（τ）的每一行。重新定义插值矩阵并且重写时域宽带信号模型为

图2 入射角φk＞π/2时的阵列信号接收模型

2.2 贝叶斯后验概率密度函数

对于所有的快拍，全似然函数为

式中k为真实声源数的一个估计。由此通过贝叶斯理论，声源参数的后验概率密度函数可以表示为全似然函数和声源参数先验函数的乘积为

式中：δ2是一个超参数，它与信噪比有密切关系；声源幅度a先验的选取参考文献［6］；τ选为均匀先验；σw的先验选为无信息先验；k的先验选为期望值为Λ的泊松分布，Λ是预期的声源数。由式（14）给出

对于上述参数，我们所感兴趣的是声源数k和时延τ.将式（14）代入式（13）并把“多余”参数a和σw积分掉可以得到［6］

式（15）即参数τ和声源数k的联合后验概率密度函数。通常需要对联合后验概率密度函数进行多维峰值搜索［3］，最大峰值所对应的模型阶数和时延值便是入射声源数和目标方位的估计值。

3 基于PMC算法的贝叶斯高分辨方位估计

目标是在贝叶斯框架下进行联合模型选择和参数估计，即通过最大化式（15）求得模型阶数k和时延τ.后验概率密度函数的最大化通常采用RJMCMC算法。而本文的方法采用 PMC算法［10，13－14］来执行上述贝叶斯计算。PMC算法是新近提出的一种有效的采样技术，是MCMC采样的另一种选择，它在构建转移核时的灵活性使得它允许对不同模型参数空间进行联合采样。相对于较常使用的RJMCMC算法，它具有如下显著优点：首先，与RJMCMC算法相比，它不需要一个样本被丢弃的“烧录”周期；其次，RJMCMC算法一次迭代产生一个样本，而PMC算法每次迭代产生一组近似无关的样本，这大大地增强了算法并行执行的能力。

3.1 PMC算法的样本初始化

从式（14）采样t＝0 时刻模型阶数k（i，0），以k（i，0）为条件采样初始时延估计τ（i，0），其中i＝1，2，…，I表示样本索引，I为样本总数

将k（i，0）和τ（i，0）依次代入式（15），计算每个样本所对应的后验概率，进行归一化得到w（i，0），

依据式（17）和式（18）的结果可以得到t＝0时刻的各个模型阶数所对应的概率和时延估计

3.2 PMC算法的迭代程序

第t次迭代时，PMC算法以t－1次迭代所得到的估计结果为基础进行样本的抽取。PMC算法以一定的概率（d＝1，…，D）从D个转移核中产生一组近似无 关 的 样 本θ（i，t）＝ ｛k（i，t），τ（i，t）｝［10］.使 用与式（18）～（20）类似的方法计算第t次迭代时模型阶数k（t）和时延估计τ（t），算法仅须把上标由0改为t.最后，进行样本重采样并计算下次迭代每个转移核的概率，并进入第t＋1次迭代程序。一个通用的PMC算法的一次迭代可概括为以下几个步骤：

1）从混合转移核产生样本θ（i，t），θ（i，t）表示t时刻第i个样本；

2）计算样本归一化权值w（i，t），w（i，t）表示t时刻第i个样本的权值；

3）基于样本θ（i，t）和权值w（i，t）执行参数估计；

如此迭代下去直到估计结果满足所需的精度或者达到最大迭代次数。

4 计算机仿真

采用如下仿真条件。使用M＝10元ULA，相应波束宽度为WB＝10.09°.阵列快拍数为N＝64.仿真考虑稍复杂的三个声源情况，三个声源具有相同的信噪比，定义见式（21），σs表示目标辐射噪声功率，σw表示阵列加性高斯噪声的标准差。

三个入射声源的方位角分别为θ1＝6.00°，θ2＝12.00°和θ3＝18.00°，相邻的两个声源间的间距接近半个波束宽度；其频谱范围分别为f1∈［600，2 800］Hz，f2∈［500，3 100］Hz和f3∈［500，3 400］Hz；系统采样频率为fs＝8 000Hz；水下声速c＝1 500m/s；阵元间距为d＝c/fs＝0.187 5m；其他参数详见表1.

表1 仿真采用的参数值

将基于PMC算法的贝叶斯高分辨水下被动目标方位估计方法称为PMCMAP（PMC maximum aposterior）算法。图3和图4给出了单次运行时PMCMAP算法的性能。PMC算法采用10次迭代，每次迭代使用100个样本。由图3可以看出：算法收敛相当迅速，到第4次迭代时算法便以很高的精度收敛到了真实的目标方位，到第10次迭代便能很好地估计出目标方位。并且由多次独立仿真可知，大多数情况下算法在首次迭代时就能正确地估计出了模型阶数为3，图3即为此类情况。图4给出了采用PMCMAP算法的声源幅度恢复结果，由图4可知三个声源的幅度都被很好地恢复出来，这为后续的基于频谱估计的目标识别打下基础。

图3 单次运行时PMCMAP算法的DOA估计结果

图5和图6给出了在信噪比为15dB时进行80次独立蒙特卡罗（MC）仿真时PMCMAP算法的性能。在80次独立仿真中有60次算法正确地估计出模型阶数为3；其余20次倾向于低估模型阶数，且估计结果均为2，如表2所示。

图4 单次运行时PMCMAP算法声源幅度恢复结果

图5给出了模型阶数估计结果为2时声源的波达方向（DOA）估计。由图5可以看出：在这20次独立的MC仿真中，算法总是大致地估计出了声源1和3而遗漏了声源2；且声源1的角度估计值倾向于偏大而声源3的角度估计值倾向于偏小。由此可以看出具有相同信噪比的三个声源中处于中间位置的声源的估计结果受到了其两侧声源的干扰，同时两侧声源的DOA估计结果同样受到了中间声源的影响而趋于向其靠近。图5中箭头朝向右下侧的移动趋势清晰地表现了这一结论。

图5 模型阶数估计为2时的方位估计结果

图6给出了正确地估计出模型阶数3时声源的DOA估计，为了表述方便，分别绘出声源1和2、声源2和3的DOA估计结果。由图6可以看出单次运行的DOA估计值（星形）绝大多数都集中地分布在DOA真实值（菱形）附近，其中心位置（三角形）稍有偏离。声源1和2的估计值中心约为［6.1°，11.7°］；而声源2和3的估计值中心约为［11.7°，18.1°］.算法很好地估计出了三个声源的方位，性能较好。此外还可以看出，声源3的估计较声源1更为集中，所以算法对较大的入射角的估计精度要大于对入射角接近于法线方向的声源的估计精度。

图7给出了不同信噪比下使用PMCMAP算法的均方根误差（RMSE）.由图7可见随着RSN的增加方位估计的精度也不断地增加，当RSN≥10dB时，三个声源方位估计的RMSE均小于1°，三个相距约为半束宽的宽带声源能被很好地区分开且被精确地定位。此外还可以看出，三个声源的估计精度大致上依次为声源3＞声源1＞声源2，这和前边得出的结论一致。表2给出了不同信噪比下，算法的模型阶数估计结果和方位估计的均方根误差。

图7 不同信噪比下方位估计的RMSE

表2 不同信噪比下的模型阶数估计结果及DOA估计的均方根误差

5 结 论

将一种基于时域宽带信号模型的贝叶斯高分辨方位估计方法［6］引入水下被动目标的方位估计当中，研究了用插值矩阵表示法来表征阵列接收的宽带水下目标辐射噪声。与传统的处理宽带信号的聚焦类子空间方法相比，该方法具有显著的优势。首先，它可以直接在时域内进行信号处理，避免了聚焦类方法的频域变换及频域子带划分；其次，它在阵列快拍数很少的情况下也能取得很好的估计结果，而相比之下聚焦类方法在少快拍情况下性能会急剧恶化；第三，它可以在进行方位估计的同时进行模型阶数（声源数）的检测和声源幅度的恢复，而计算量增加不大，其中信号幅度的恢复有助于基于谱估计的目标识别。

同时将最新的并行采样算法—PMC算法引入到水下目标方位估计领域，采用PMC算法对求得的目标参数的后验概率密度函数进行峰值搜索。相对于常用的RJMCMC采样算法，PMC算法具有显著的优势。首先，它不需要一个样本被丢弃的“烧录”周期；其次，它避免了传统的MCMC方法一次一个候选的采样方式，PMC算法每次迭代产生一组样本，可以大大地提高采样效率，从而大大增强了算法并行执行的能力。

选取了较复杂的三个目标辐射声源的情况进行计算机仿真。仿真结果显示：采用基于时域宽带信号模型和PMC方法的PMCMAP算法在较少的快拍情况下不仅能够准确地估计出声源数，同时可以精确地估计出目标方位，并且声源幅度也能被很好地恢复出来，取得了很好的估计效果。

［1］杨 鹏，杨 峰，聂在平，等.MUSIC算法在柱面共形天线阵DOA估计中的应用研究［J］.电波科学学报，2008，23（2）：288－291.YANG Peng，YANG Feng，NIE Zaiping，et al.DOA estimation of cylindrical conformal array by MUSIC algorithm［J］.Chinese Journal of Radio Science，2008，23（2）：288－291.（in Chinese）

［2］陶建武，石要武，常文秀.一般阵列误差情况下信号二维方向角估计［J］.电波科学学报，2006，21（4）：606－611.TAO Jianwu，SHI Yaowu，CHANG Wenxiu.Estimation of 2Dangle for signals with general array error［J］.Chinese Journal of Radio Science，2006，21（4）：606－611.（in Chinese）

［3］LASENBY J，FITZGERALD W J.A Bayesian approach to high－resolution beamforming［J］.Radar and Signal Processing，IEE Proceedings－F，1991，138（6）：539－544.

［4］HUANG Jianguo，SUN Yi，LIU Kewei，et al.A new Gibbs sampling DOA estimator based on Bayesian method［C］//IEEE Conference on Acoustics，Speech and Signal Processing，2003，5：229－232.

［5］徐 朴，黄建国.波达方向估计的贝叶斯高分辨方法［J］.电波科学学报，2001，16（2）：149－152.XU Pu，HUANG Jianguo.A Bayesian high－resolution technique for estimation of DOA［J］.Chinese Journal of Radio Science，2001，16（2）：149－152.（in Chinese）

［6］NG W，REILLY J P，KIRUBARAJAN T，et al.Wideband array signal processing using MCMC methods［J］.IEEE Transactions on Signal Processing，2005，53（2）：411－426.

［7］ANDRIEU C，DOUCET A.Joint Bayesian model selection and estimation of noisy sinusoids via reversible jump MCMC［J］.IEEE Transactions on Signal Processing，1999，47（10）：2667－2676.

［8］GREEN P J.Reversible jump Markov Chain Monte Carlo computation and Bayesian model determination［J］.Biometrika，1995，82（4）：711－732.

［9］LAROCQUE J R，REILLY J P.Reversible jump MCMC for joint detection and estimation of sources in colored noise［J］.IEEE Transactions on Signal Processing，2002，50（2）：231－240.

［10］HONG Mingyi，BUGALLO M F，DJURIC P M.Joint model selection and parameter estimation by population Monte Carlo simulation［J］.IEEE Journal of Selected Topics in Signal Processing，2010，4（3）：526－539.

［11］廖宏宇.被动声纳目标/背景建模与仿真［J］.计算机仿真，2006，23（4）：1－4.LIAO Hongyu. Object－background modeling and simulation for passive sonar［J］.Computer Simulation，2006，23（4）：1－4.（in Chinese）

［12］梁民赞，王旭升，李志伟.舰船辐射噪声动态特性建模与重构［J］.舰船电子工程，2008，28（2）：95－97.LIANG Minzan，WANG Xusheng，LI Zhiwei.Dynamic properties modeling of radiating noise from vessels［J］.Ship Electronic Engineering，2008，28（2）：95－97.（in Chinese）

［13］BUGALLO M F，HONG Mingyi，DJURIC P M.Marginalized Population Monte Carlo［C］//Proceedings of the 2009IEEE International Conference on Acoustics，Speech and Signal Processing，2009.2925－2928.

［14］LASKEY K B，MYERS J W.Population Markov Chain Monte Carlo［J］.Machine Learning，2003，50（1/2）：175－196.