论文部分内容阅读
摘 要： 与传统波叠加法通过逆问题求解获得源强不同，提出一种基于粒子滤波的波叠加算法来实现正向求解源强。该方法基于波叠加理论建立状态空间模型，根据初始化粒子的先验信息预测并更新状态向量。通过权重计算和重采样，估计等效源的强度及其位置，进而重建三维辐射声场。仿真分析与实验验证表明，该算法可避免逆问题求解中的不适定性难题，实现声场的高精度重建。
关键词： 粒子滤波; 波叠加; 声源定位; 声场重建
文章编号： 2095-2163（2021）03-0080-05 中图分类号：TP391.9; TB52 文献标志码：A
【Abstract】Different from the traditional wave superposition method which is utilized to obtain the source intensity by solving the inverse problem， an innovative wave superposition algorithm based on particle filtering is proposed to realize the forward solving of source intensity. This method presents a state space model based on the wave superposition theory， then updates the state vector according to the prior information prediction of initializing particle. Through weight calculation and resampling， the intensity and the position of the equivalent source can be estimated. Thus， the three-dimensional sound field can be reconstructed. Simulation analysis and experimental verification show that the algorithm can avoid the ill-posed problem of solving the inverse problem and achieve high-precision reconstruction of the sound field.
【Key words】 particle filter; wave superposition; sound source localization; sound field reconstruction
0 引 言
波疊加法由Koopmann等人[1]在1989年提出，是一种求解声辐射问题的积分方程法，其原理为：任何物体辐射的声场均可以等效为该辐射体内部若干个不同源强的等效源所产生声场的叠加。利用波叠加法近场声全息进行声场重建，既能克服基于空间Fourier变换法的近场声全息[2]只能计算规则形状声源的缺点，又能避免基于边界元法的近场声全息[3]所带来的复杂插值运算和奇异积分处理困扰，国内外学者对该方法进行了深入研究[4-8]。
利用波叠加法进行声场重建时，其关键是获取等效源强度和配置等效源位置。通常采用Tikhonov正则化方法[9-10]求解声学逆问题中的病态方程，获取准确的等效源强度，再合理配置等效源位置，进行声场重建。Antoni[11]提出了贝叶斯正则化准则的基本理论，Pereira等人[12]进一步通过理论和实验对比了等效源法和贝叶斯法。Chardon等人[13]针对传统近场声全息逆问题正则化方法的复杂性，提出了一种新的基于压缩采样的正则化方法。Yang等人[14]提出了一种基于卡尔曼滤波的混合局部近场声全息技术，旨在解决等效源配置不确定性问题。
以上所提方法均是从求解声学逆问题的角度出发获得源强，为避免求解逆问题的不适定性，Bai等人[15-16]提出基于等效源模型的自适应滤波方法，对麦克风阵列信号进行处理，该方法利用递归方式估计未知源振幅，属于声学正问题。本文提出一种基于粒子滤波的波叠加方法（Particle Filter based WSM，PF-WSM），该方法以波叠加法为理论基础，建立状态空间模型，利用粒子滤波算法从正向估计等效源强度和位置，根据估计得到参数计算重建面声压值。本文安排如下：首先对PF-WSM算法进行理论描述，再仿真分析声源位置已知和未知两种情况的声场重建结果，最后通过实验验证，证明了该方法的可行性。
1 理论背景
1.1 PF-WSM算法理论
传统波叠加法通过Tikhonov正则化方法求解声学逆问题中的病态方程，获得准确的等效源强度，故被称为Tikhonov正则化波叠加（TR-WSM）算法。与TR-WSM算法反解源强相比，PF-WSM算法正向估计源强，避免了求解逆问题的不适定性。
粒子滤波（Particle Filter，PF）是一种基于蒙特卡洛方法进行递推贝叶斯估计的滤波算法[17-18]。通过采用一组具有相应权值的随机粒子估计状态的前验分布，并通过观测值逼近状态的后验分布，再求出这组粒子的加权求和作为状态的估计值。
在PF-WSM算法中，用等效源强度和位置描述系统状态，源强与传递矩阵计算的声压描述观测状态，由此，状态空间模型表示为：
其中，qn表示源强的系统状态矢量;sn表示源位置的系统状态矢量且不随递归变化;yn表示观测的声压矢量;n表示估计状态递归次数。源强的系统噪声和观测噪声服从均值为零、方差分别为σ21和σ23的正态分布，即v1～N0，σ21和v3～N0，σ23;源位置的系统噪声服从均匀分布，即v2～U-σ2，σ2。G表示等效源到全息面的传递矩阵。 第i个粒子的权值更新方程表示为：
当经过一定次数的估计状态递归后，会存在粒子退化问题，可通过重采样[19]方法来解决。该方法是对后验概率密度的离散近似表示再进行一次采样，繁殖权值较高的样本，淘汰权值较低的样本，重新生成一个新样本集合，以克服样本权值退化问题。
1.2 PF-WSM算法流程
粒子滤波的核心是利用优选随机粒子来实现状态分布的估计，为更好地进行源强分布估计，给定粒子初始状态分布为高斯分布，这是因为根据中心极限定理，随着粒子数量的增加，大多数粒子均值变得相似并越来越接近总体均值。均值为μ和方差为σ2的高斯分布的概率密度函数为：
对于复数域的源强概率密度函数，式（7）可以看作是一个联合概率密度函数的实部和虚部，即：
利用PF-WSM算法在2种情况下进行声场重建对比，根据公式（1）～公式（3）建立状态空间模型。对这2种情况可简述为：
（1）声源位置已知，其观测方程中声压值的变化仅与粒子的源强有关。
（2）声源位置未知，声压值的变化依赖于粒子的源强和源位置。
文中给出了PF-WSM算法的流程如图1所示。具体过程如下：
（1）粒子初始化：第一次递归时，从高斯分布Nμ，σ2中随机产生K个粒子qini=1，2，…，K，并平均分配权重。
（2）系统状态矢量预测：根据式（1）～式（2）所示的状态转移方程，计算粒子qin在n 1次递归的预测粒子qin 1，服从一阶马尔可夫过程。
（3）粒子权重更新：利用式（3）计算n 1次粒子预测值对应的观测声压，同时计算观测声压与n 1次的实际测量声压之间的差值，根据差值更新每个粒子的权重（式（4）～式（5）），再将权重进行归一化。
（4）源强分布估计：利用当前递归次数中每个粒子的权重win 1以及取值qin 1，估计系统状态矢量（源强分布），计算公式为：
（5）重采样：为了快速有效地选择接近真实源强分布的粒子，在进行下一次系统状态估计之前，需要按照粒子的权重对粒子进行重采样。重复上述步骤，即可基于粒子滤波算法构建系统状态估计器，实现源强分布估计。
（6）声场重建：根据上述构建的系统状态估计器，得到合理的估计源强分布，利用源强与重建面之间的传递矩阵求得重建声压：
其中，H表示等效源面到重建面的传递矩阵，Q表示PF-WSM算法估计源强分布。
2 数值仿真
2.1 仿真设置
为验证PF-WSM算法在声场重建中的可行性，将其与Tikhonov正则化波叠加 （TR-WSM）算法进行比较。单极子源的半径为0.01 m，表面脉动速度为0.2 m/s，由25个单极子声源组成的辐射声场布置在xoy面，以声源面中心为坐标原点o，其法面方向为z轴，建立空间直角坐标系。全息面平行于xoy面，距坐标原点0.2 m，在0.8×0.8 m2的平面上均匀布置25个传声器，间距为0.2 m。重建面尺寸为1×1 m2，网格间距为0.1 m，距坐标原点0.05 m。对全息面的测量数据加上随机高斯白噪声，假设声场中的信噪比为20 dB。设置递归次数N为90次，粒子數K为100个。
仿真过程中，相对误差的计算公式如下：
2.2 仿真分析
声源位置已知时，假设25个单极子源均匀分布在声源面上，如图2（a）所示，其源位置与全息面上传声器位置在xoy平面的俯视图上一一对应。根据1.2节的算法流程计算出重建声压，图2（b）比较了100个粒子的PF-WSM算法和TR-WSM算法声压重建误差。源位置均匀分布在网格点上，PF-WSM算法基于源强为高斯分布的先验信息仅估计源强，与TR-WSM算法相比，PF-WSM算法在声压重建中误差较大。
声源准确位置已知的情况下，TR-WSM算法的性能要优于PF-WSM算法。但在实际场景中，通常无法预先知道声源的准确位置，对源位置假设不正确，导致声场重建值偏离实际值，因此假设声源位置已知是一种理想化的情况。由此，假设源位置是未知的，并且存在位置扰动，更符合真实情况。
声源位置未知时，假设25个单极子源分布在声源面上存在扰动，源位置相对于传声器分布位置如图3（a）所示。图3（b）是PF-WSM算法和TR-WSM算法的重建误差对比。在源位置未知且存在扰动的情况下，PF-WSM算法同时估计源强和源位置，基于源强服从高斯分布和源位置服从均匀分布的先验信息，在不断递归之后，源位置越接近声源真实位置，源强越接近真实源强，PF-WSM算法的声场重建结果越接近真实辐射声场。而TR-WSM算法在源位置未知的情况下，由于等效源配置不准确，重建误差增大。图3（b）中PF-WSM算法的重建误差小于TR-WSM算法，说明PF-WSM算法在源强和源位置均未知的情况下，估计效果更好，而TR-WSM算法在未知源位置时的重建效果较差，其应用有一定局限性。图4给出了2种算法在100 Hz时的声场重建效果对比，图4（a）为理论声压，图4（b）为PF-WSM算法的重建声压，图4（c）为TR-WSM算法的重建声压。由图4可以看出，PF-WSM算法的声场重建效果更好，其热点最大值与理论声压相同，说明粒子滤波可以较准确地估计出源强和源位置，而TR-WSM算法的热点最大值小于理论值，这是因为不知道声源位置时，等效源配置不准确，影响声场重建效果。
3 实验验证
为进一步检验算法的可行性，在半消声室进行实验验证。实验设置如图5（a）所示。半消声室尺寸为9.8 m×8.6 m×3.5 m，背景噪声为18 dB（A），截止频率为125 Hz。采用音箱模拟噪声源，频率为300 Hz，音箱纸盆中心位置为（0， 0， 0）。以音箱纸盆中心为坐标原点，音箱平面法向为z轴建立空间直角坐标系。传声器阵列平面平行于声源面，阵列中心距坐标原点0.6 m，离地面高度为1.2 m。传声器阵列如图5（b）所示，采用36个传声器组成3个同心圆的面阵列，3个同心圆半径分别为：0.07 m、0.18 m、0.25 m。重建面大小与传声器阵列面相同，重建面距坐标原点0.05 m。PF-WSM算法中，粒子的先验信息是基于高斯分布，36个虚源均匀分布在0.25 m×0.25 m的正方形平面上，间隔为0.05 m，距坐标原点-0.01 m。将PF-WSM算法和TR-WSM算法重建的声压与传声器测得的声压进行对比。图6为实验验证结果。图6（a）为传声器阵列测得的声压分布，图6（b）和图6（c）分别为PF-WSM算法和TR-WSM算法重建的声压分布。由式（11）计算PF-WSM算法和TR-WSM算法重建误差，分别为5.89%和13.43%，可以得出结果：PF-WSM算法的重建效果优于TR-WSM算法的重建效果。 4 结束语
本文提出一种基于粒子滤波的波叠加算法进行声场重建。该方法在源位置未知且随机的情况下，能够有效利用参数的先验知识，通过粒子滤波估计器，达到求解源强的效果。特别是在声源位置未知且存在扰动的情况下，源位置是非高斯分布的，传统波叠加法由于等效源配置不准确而声场重建效果不好，PF-WSM算法则利用这一特点，可以更好地重建声场。PF-WSM算法不仅可以估计源强，还可以估计源位置，甚至在源位置未知的情况下，估计性能更好，即对源位置的估计越准确，重建的声压就越准确。另一方面，粒子数的增加直接影响估计状态向量的精度，粒子数越多，估计状态越接近真实状态。本文在权衡运算时间和估计精度之后，选择100个粒子进行估计，若增加粒子数（如1 000个粒子），估计值将越接近真实值，但运算时间相应增加，后续可采用改进粒子滤波算法进行优化，降低采样复杂度，提升运算速度。
参考文献
[1]
关键词： 粒子滤波; 波叠加; 声源定位; 声场重建
文章编号： 2095-2163（2021）03-0080-05 中图分类号：TP391.9; TB52 文献标志码：A
【Abstract】Different from the traditional wave superposition method which is utilized to obtain the source intensity by solving the inverse problem， an innovative wave superposition algorithm based on particle filtering is proposed to realize the forward solving of source intensity. This method presents a state space model based on the wave superposition theory， then updates the state vector according to the prior information prediction of initializing particle. Through weight calculation and resampling， the intensity and the position of the equivalent source can be estimated. Thus， the three-dimensional sound field can be reconstructed. Simulation analysis and experimental verification show that the algorithm can avoid the ill-posed problem of solving the inverse problem and achieve high-precision reconstruction of the sound field.
【Key words】 particle filter; wave superposition; sound source localization; sound field reconstruction
0 引 言
波疊加法由Koopmann等人[1]在1989年提出，是一种求解声辐射问题的积分方程法，其原理为：任何物体辐射的声场均可以等效为该辐射体内部若干个不同源强的等效源所产生声场的叠加。利用波叠加法近场声全息进行声场重建，既能克服基于空间Fourier变换法的近场声全息[2]只能计算规则形状声源的缺点，又能避免基于边界元法的近场声全息[3]所带来的复杂插值运算和奇异积分处理困扰，国内外学者对该方法进行了深入研究[4-8]。
利用波叠加法进行声场重建时，其关键是获取等效源强度和配置等效源位置。通常采用Tikhonov正则化方法[9-10]求解声学逆问题中的病态方程，获取准确的等效源强度，再合理配置等效源位置，进行声场重建。Antoni[11]提出了贝叶斯正则化准则的基本理论，Pereira等人[12]进一步通过理论和实验对比了等效源法和贝叶斯法。Chardon等人[13]针对传统近场声全息逆问题正则化方法的复杂性，提出了一种新的基于压缩采样的正则化方法。Yang等人[14]提出了一种基于卡尔曼滤波的混合局部近场声全息技术，旨在解决等效源配置不确定性问题。
以上所提方法均是从求解声学逆问题的角度出发获得源强，为避免求解逆问题的不适定性，Bai等人[15-16]提出基于等效源模型的自适应滤波方法，对麦克风阵列信号进行处理，该方法利用递归方式估计未知源振幅，属于声学正问题。本文提出一种基于粒子滤波的波叠加方法（Particle Filter based WSM，PF-WSM），该方法以波叠加法为理论基础，建立状态空间模型，利用粒子滤波算法从正向估计等效源强度和位置，根据估计得到参数计算重建面声压值。本文安排如下：首先对PF-WSM算法进行理论描述，再仿真分析声源位置已知和未知两种情况的声场重建结果，最后通过实验验证，证明了该方法的可行性。
1 理论背景
1.1 PF-WSM算法理论
传统波叠加法通过Tikhonov正则化方法求解声学逆问题中的病态方程，获得准确的等效源强度，故被称为Tikhonov正则化波叠加（TR-WSM）算法。与TR-WSM算法反解源强相比，PF-WSM算法正向估计源强，避免了求解逆问题的不适定性。
粒子滤波（Particle Filter，PF）是一种基于蒙特卡洛方法进行递推贝叶斯估计的滤波算法[17-18]。通过采用一组具有相应权值的随机粒子估计状态的前验分布，并通过观测值逼近状态的后验分布，再求出这组粒子的加权求和作为状态的估计值。
在PF-WSM算法中，用等效源强度和位置描述系统状态，源强与传递矩阵计算的声压描述观测状态，由此，状态空间模型表示为：
其中，qn表示源强的系统状态矢量;sn表示源位置的系统状态矢量且不随递归变化;yn表示观测的声压矢量;n表示估计状态递归次数。源强的系统噪声和观测噪声服从均值为零、方差分别为σ21和σ23的正态分布，即v1～N0，σ21和v3～N0，σ23;源位置的系统噪声服从均匀分布，即v2～U-σ2，σ2。G表示等效源到全息面的传递矩阵。 第i个粒子的权值更新方程表示为：
当经过一定次数的估计状态递归后，会存在粒子退化问题，可通过重采样[19]方法来解决。该方法是对后验概率密度的离散近似表示再进行一次采样，繁殖权值较高的样本，淘汰权值较低的样本，重新生成一个新样本集合，以克服样本权值退化问题。
1.2 PF-WSM算法流程
粒子滤波的核心是利用优选随机粒子来实现状态分布的估计，为更好地进行源强分布估计，给定粒子初始状态分布为高斯分布，这是因为根据中心极限定理，随着粒子数量的增加，大多数粒子均值变得相似并越来越接近总体均值。均值为μ和方差为σ2的高斯分布的概率密度函数为：
对于复数域的源强概率密度函数，式（7）可以看作是一个联合概率密度函数的实部和虚部，即：
利用PF-WSM算法在2种情况下进行声场重建对比，根据公式（1）～公式（3）建立状态空间模型。对这2种情况可简述为：
（1）声源位置已知，其观测方程中声压值的变化仅与粒子的源强有关。
（2）声源位置未知，声压值的变化依赖于粒子的源强和源位置。
文中给出了PF-WSM算法的流程如图1所示。具体过程如下：
（1）粒子初始化：第一次递归时，从高斯分布Nμ，σ2中随机产生K个粒子qini=1，2，…，K，并平均分配权重。
（2）系统状态矢量预测：根据式（1）～式（2）所示的状态转移方程，计算粒子qin在n 1次递归的预测粒子qin 1，服从一阶马尔可夫过程。
（3）粒子权重更新：利用式（3）计算n 1次粒子预测值对应的观测声压，同时计算观测声压与n 1次的实际测量声压之间的差值，根据差值更新每个粒子的权重（式（4）～式（5）），再将权重进行归一化。
（4）源强分布估计：利用当前递归次数中每个粒子的权重win 1以及取值qin 1，估计系统状态矢量（源强分布），计算公式为：
（5）重采样：为了快速有效地选择接近真实源强分布的粒子，在进行下一次系统状态估计之前，需要按照粒子的权重对粒子进行重采样。重复上述步骤，即可基于粒子滤波算法构建系统状态估计器，实现源强分布估计。
（6）声场重建：根据上述构建的系统状态估计器，得到合理的估计源强分布，利用源强与重建面之间的传递矩阵求得重建声压：
其中，H表示等效源面到重建面的传递矩阵，Q表示PF-WSM算法估计源强分布。
2 数值仿真
2.1 仿真设置
为验证PF-WSM算法在声场重建中的可行性，将其与Tikhonov正则化波叠加 （TR-WSM）算法进行比较。单极子源的半径为0.01 m，表面脉动速度为0.2 m/s，由25个单极子声源组成的辐射声场布置在xoy面，以声源面中心为坐标原点o，其法面方向为z轴，建立空间直角坐标系。全息面平行于xoy面，距坐标原点0.2 m，在0.8×0.8 m2的平面上均匀布置25个传声器，间距为0.2 m。重建面尺寸为1×1 m2，网格间距为0.1 m，距坐标原点0.05 m。对全息面的测量数据加上随机高斯白噪声，假设声场中的信噪比为20 dB。设置递归次数N为90次，粒子數K为100个。
仿真过程中，相对误差的计算公式如下：
2.2 仿真分析
声源位置已知时，假设25个单极子源均匀分布在声源面上，如图2（a）所示，其源位置与全息面上传声器位置在xoy平面的俯视图上一一对应。根据1.2节的算法流程计算出重建声压，图2（b）比较了100个粒子的PF-WSM算法和TR-WSM算法声压重建误差。源位置均匀分布在网格点上，PF-WSM算法基于源强为高斯分布的先验信息仅估计源强，与TR-WSM算法相比，PF-WSM算法在声压重建中误差较大。
声源准确位置已知的情况下，TR-WSM算法的性能要优于PF-WSM算法。但在实际场景中，通常无法预先知道声源的准确位置，对源位置假设不正确，导致声场重建值偏离实际值，因此假设声源位置已知是一种理想化的情况。由此，假设源位置是未知的，并且存在位置扰动，更符合真实情况。
声源位置未知时，假设25个单极子源分布在声源面上存在扰动，源位置相对于传声器分布位置如图3（a）所示。图3（b）是PF-WSM算法和TR-WSM算法的重建误差对比。在源位置未知且存在扰动的情况下，PF-WSM算法同时估计源强和源位置，基于源强服从高斯分布和源位置服从均匀分布的先验信息，在不断递归之后，源位置越接近声源真实位置，源强越接近真实源强，PF-WSM算法的声场重建结果越接近真实辐射声场。而TR-WSM算法在源位置未知的情况下，由于等效源配置不准确，重建误差增大。图3（b）中PF-WSM算法的重建误差小于TR-WSM算法，说明PF-WSM算法在源强和源位置均未知的情况下，估计效果更好，而TR-WSM算法在未知源位置时的重建效果较差，其应用有一定局限性。图4给出了2种算法在100 Hz时的声场重建效果对比，图4（a）为理论声压，图4（b）为PF-WSM算法的重建声压，图4（c）为TR-WSM算法的重建声压。由图4可以看出，PF-WSM算法的声场重建效果更好，其热点最大值与理论声压相同，说明粒子滤波可以较准确地估计出源强和源位置，而TR-WSM算法的热点最大值小于理论值，这是因为不知道声源位置时，等效源配置不准确，影响声场重建效果。
3 实验验证
为进一步检验算法的可行性，在半消声室进行实验验证。实验设置如图5（a）所示。半消声室尺寸为9.8 m×8.6 m×3.5 m，背景噪声为18 dB（A），截止频率为125 Hz。采用音箱模拟噪声源，频率为300 Hz，音箱纸盆中心位置为（0， 0， 0）。以音箱纸盆中心为坐标原点，音箱平面法向为z轴建立空间直角坐标系。传声器阵列平面平行于声源面，阵列中心距坐标原点0.6 m，离地面高度为1.2 m。传声器阵列如图5（b）所示，采用36个传声器组成3个同心圆的面阵列，3个同心圆半径分别为：0.07 m、0.18 m、0.25 m。重建面大小与传声器阵列面相同，重建面距坐标原点0.05 m。PF-WSM算法中，粒子的先验信息是基于高斯分布，36个虚源均匀分布在0.25 m×0.25 m的正方形平面上，间隔为0.05 m，距坐标原点-0.01 m。将PF-WSM算法和TR-WSM算法重建的声压与传声器测得的声压进行对比。图6为实验验证结果。图6（a）为传声器阵列测得的声压分布，图6（b）和图6（c）分别为PF-WSM算法和TR-WSM算法重建的声压分布。由式（11）计算PF-WSM算法和TR-WSM算法重建误差，分别为5.89%和13.43%，可以得出结果：PF-WSM算法的重建效果优于TR-WSM算法的重建效果。 4 结束语
本文提出一种基于粒子滤波的波叠加算法进行声场重建。该方法在源位置未知且随机的情况下，能够有效利用参数的先验知识，通过粒子滤波估计器，达到求解源强的效果。特别是在声源位置未知且存在扰动的情况下，源位置是非高斯分布的，传统波叠加法由于等效源配置不准确而声场重建效果不好，PF-WSM算法则利用这一特点，可以更好地重建声场。PF-WSM算法不仅可以估计源强，还可以估计源位置，甚至在源位置未知的情况下，估计性能更好，即对源位置的估计越准确，重建的声压就越准确。另一方面，粒子数的增加直接影响估计状态向量的精度，粒子数越多，估计状态越接近真实状态。本文在权衡运算时间和估计精度之后，选择100个粒子进行估计，若增加粒子数（如1 000个粒子），估计值将越接近真实值，但运算时间相应增加，后续可采用改进粒子滤波算法进行优化，降低采样复杂度，提升运算速度。
参考文献
[1]