CN115291169A - 一种声源成像方法、系统、设备及存储介质 - Google Patents
一种声源成像方法、系统、设备及存储介质 Download PDFInfo
- Publication number
- CN115291169A CN115291169A CN202210932721.0A CN202210932721A CN115291169A CN 115291169 A CN115291169 A CN 115291169A CN 202210932721 A CN202210932721 A CN 202210932721A CN 115291169 A CN115291169 A CN 115291169A
- Authority
- CN
- China
- Prior art keywords
- sound source
- model
- matrix
- sound
- imaging
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 107
- 238000000034 method Methods 0.000 claims abstract description 87
- 238000005259 measurement Methods 0.000 claims abstract description 27
- 230000006835 compression Effects 0.000 claims abstract description 7
- 238000007906 compression Methods 0.000 claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims description 61
- 239000013598 vector Substances 0.000 claims description 24
- 230000006870 function Effects 0.000 claims description 21
- 238000004590 computer program Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 10
- 238000007476 Maximum Likelihood Methods 0.000 claims description 5
- 238000004088 simulation Methods 0.000 claims description 4
- 238000012546 transfer Methods 0.000 claims description 4
- 238000004458 analytical method Methods 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 14
- 230000008901 benefit Effects 0.000 abstract description 4
- 238000004422 calculation algorithm Methods 0.000 description 10
- 238000010586 diagram Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 4
- 238000013459 approach Methods 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 208000001992 Autosomal Dominant Optic Atrophy Diseases 0.000 description 1
- 206010011906 Death Diseases 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000003491 array Methods 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000835 fiber Substances 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
- G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
- G01S5/22—Position of source determined by co-ordinating a plurality of position lines defined by path-difference measurements
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
本发明涉及声学成像技术领域，尤其涉及一种声源成像方法、系统、设备及存储介质，包括：以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。本发明通过二维稀疏贝叶斯学习压缩波束形成方法求解声源成像问题模型，能够准确地估计声源强度，且具有方法简单、计算效率高、精度高等优点。
Description
技术领域
本发明涉及声学成像技术领域，尤其涉及一种声源成像方法、系统、设备及存储介质。
背景技术
基于平面传声器阵列的二维压缩波束形成声源成像方法因其具有成像空间范围广、对相干和不相干声源均适用且成像清晰等优点备受关注，其通过建立并求解传声器测量声压与目标区域离散化网格点源强间的欠定方程组实现声源成像，凸松弛算法和贪婪算法常被用于求解上述欠定方程组。
然而，基于凸松弛算法的声源成像方法需严格满足约束等距条件且计算复杂度大、计算速度慢，对于大规模问题不适用；而基于贪婪算法的声源成像方法则可能获得局部最优解而非全局最优解，且其易受感知矩阵列相干性影响，性能对网格间距敏感。
发明内容
本发明的目的在于提供一种声源成像方法、系统、设备及存储介质，以同时兼顾声源成像性能稳健性和高计算效率，使其可以适用于大规模问题。
为解决以上技术问题，本发明提供了一种声源成像方法、系统、设备及存储介质。
第一方面，本发明提供了一种声源成像方法，所述方法包括以下步骤：
以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；
将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；
利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
在进一步的实施方案中，所述将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型的步骤包括：
压缩波束形成将平面传声器阵列测量模型中的目标声源区域沿声源仰角和声源方位角两个方向，以预设间隔离散为若干个固定网格点，且所述固定网格点的数量大于传声器的数量；
假设声源位置估计与固定网格点一致，建立声源成像问题模型。
在进一步的实施方案中，所述声源成像问题模型具体为：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝[a(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集，L为快拍总数。
在进一步的实施方案中，所述利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出包括：
假设所述源强矩阵满足复高斯分布，建立源强矩阵先验分布模型和噪声先验模型；
基于全概率公式，根据所述源强矩阵先验分布模型和所述噪声先验模型得到测量声压向量分布；
基于第二类极大似然估计，根据所述测量声压向量分布得到超参数估计目标函数；
对所述超参数估计目标函数微分逐步更新迭代，得到迭代过程中声源位置估计处的声源能量估计；
若噪声信息已知，则得到波束形成的功率输出；若噪声信息未知，则通过迭代计算噪声估计方差，直至达到预定的收敛条件，则停止迭代。
在进一步的实施方案中，所述超参数估计目标函数具体为：
式中，σSg为源强矩阵的第g个元素的方差；p(P★)为测量声压向量的分布；tr(·)为矩阵的迹；det(·)为矩阵的行列式；ΣP∈CQ×Q为传声器声压P★的协方差矩阵。
在进一步的实施方案中，所述方法还包括：利用所述二维稀疏贝叶斯学习压缩波束形成方法进行声源成像性能仿真分析。
第二方面，本发明提供了一种声源成像系统，所述系统包括：
测量模型建立模块，用于以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；
成像问题建立模块，用于将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；
成像问题求解模块，用于利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
在进一步的实施方案中，所述声源成像问题模型具体为：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝[a(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集，L为快拍总数。
第三方面，本发明还提供了一种计算机设备，包括处理器和存储器，所述处理器与所述存储器相连，所述存储器用于存储计算机程序，所述处理器用于执行所述存储器中存储的计算机程序，以使得所述计算机设备执行实现上述方法的步骤。
第四方面，本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质中存储有计算机程序，所述计算机程序被处理器执行时实现上述方法的步骤。
本发明提供了一种声源成像方法、系统、设备及存储介质，所述方法通过将二维稀疏贝叶斯学习压缩波束形成方法应用于求解声源成像问题模型，解决了传统用于求解欠定方程组的算法不仅需要严格满足约束等距条件且计算复杂度大、计算速度慢，对于大规模问题不适用；而且容易获得局部最优解而非全局最优解，易受感知矩阵列相干性影响，性能对网格间距敏感的问题。相比于现有的凸松弛算法和贪婪算法，稀疏贝叶斯学习方法无显式的稀疏约束，而是通过单独地缩放每个网格点源强对应方差来隐含稀疏促进，兼顾了声源成像性能的稳健性和高计算效率。
附图说明
图1是本发明实施例提供的声源成像方法流程示意图；
图2是本发明实施例提供的平面波传导声场中阵列采样声信号模型示意图；
图3是本发明实施例提供的声源成像仿真试验布局示意图；
图4是本发明实施例提供的三种方法在3000Hz和1000Hz下的声源成像云图；
图5是本发明实施例提供的声源成像系统框图；
图6是本发明实施例提供的计算机设备的结构示意图。
具体实施方式
下面结合附图具体阐明本发明的实施方式，实施例的给出仅仅是为了说明目的，并不能理解为对本发明的限定，包括附图仅供参考和说明使用，不构成对本发明专利保护范围的限制，因为在不脱离本发明精神和范围基础上，可以对本发明进行许多改变。
参考图1，本发明实施例提供了一种声源成像方法，如图1所示，该方法包括以下步骤：
S1.以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型。
如图2所示，本实施例以由若干传声器组成的平面阵列的中心为原点、平面阵列为xoy平面，建立笛卡尔坐标系；在笛卡尔坐标系中，本实施例设置“●”代表传声器；l＝1，2，...，L为快拍索引，L为快拍总数；Si为第i个声源；θ代表仰角，其为声源入射方向与z轴的夹角；φ代表方位角，其为声源入射方向在xoy平面上的投影与x轴的夹角，且0°≤θ≤90°，0°≤φ≤360°，得到声源的DOA(声源位置估计)为(θ，φ)。
S2.将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型。
在本实施例中，所述平面阵列中的所有传声器均在xoy平面上，压缩波束形成将平面传声器阵列测量模型中的目标声源区域(xoy平面)沿声源仰角θ和声源方位角φ两个方向，按照预设间隔离散为G个固定网格点第g个固定网格点对应的声源DOA为(θGg,φGg)，令ΩG＝[[θG1,φG1]T,…,[θGg,φGg]T,...,[θGG,φGG]T]∈R2×G，其中，上标T表示转置运算符。
假设声源位置估计与固定网格点重合，将声源成像问题转换为如下的声源成像问题模型：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝[a(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点(潜在声源点)与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集。
其中，
sg＝[qg,1,…,qg,l,…,qg,L]∈C1×L
式中，qg,l∈C为第l快拍下第g个网格点的声源强度，此处，声源强度qg,l通过声源在坐标原点处产生的复声压来表示。
S3.利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
由于传声器数量通常远少于离散的网格点数量，因此，声源成像问题模型是严重欠定的线性方程组，是典型的不适定逆成像问题，不存在解析解；为了兼顾声源成像性能稳健性和高计算效率，本发明实施例将稀疏贝叶斯学习应用于求解上述欠定方程组，提出了二维稀疏贝叶斯学习压缩波束形成方法。
二维稀疏贝叶斯学习压缩波束形成方法(2D-SBL-CB)使用分层的两级贝叶斯推理求解声源成像问题模型；对于二维稀疏贝叶斯学习压缩波束形成方法，其第一级首先通过假设的信号先验模型和数据似然推理源强的后验概率分布；第二级通过最大化观测数据的概率进行超参数自适应学习以获得稀疏估计，其自适应地从观测信号中学习超参数，进而获得稀疏稳健的估计。
在一个实施例中，所述利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出包括：
假设所述源强矩阵满足复高斯分布，建立源强矩阵先验分布模型和噪声先验模型；
基于全概率公式，根据所述源强矩阵先验分布模型和所述噪声先验模型得到测量声压向量分布；
基于第二类极大似然估计，根据所述测量声压向量分布得到超参数估计目标函数；
对所述超参数估计目标函数微分逐步更新迭代，得到迭代过程中声源位置估计处的声源能量估计；
若噪声信息已知，则得到波束形成的功率输出；若噪声信息未知，则通过迭代计算噪声估计方差，直至达到预定的收敛条件，则停止迭代。
具体地，本实施例基于声源成像问题模型，2D-SBL-CB(二维稀疏贝叶斯学习压缩波束形成)方法假设源强S:,l(S:,l为源强矩阵S的第l列)满足均值为0、协方差矩阵为对角阵Γ＝Diag(σS)∈CG×G的复高斯分布，其中，Diag(·)表示将括号中的向量转化为以其元素为对角元素的对角矩阵；σS＝[σS1,...,σSg,...,σSG]，σSg表示源强S:,l的第g个元素的方差，用于控制声源强度分布的稀疏度以及稀疏剖面，是待求解的超参数；多数据快拍时，源强矩阵S中各列向量具有相同的稀疏剖面，且假设各快拍下模型参数独立，因此，声源源强矩阵的先验分布模型可写为：
式中，Π表示连乘运算符。
类似地，假设噪声矩阵也服从均值为0的复高斯分布，且各个传声器或者快拍间噪声独立，则多数据快拍下噪声的先验分布模型可写为：
式中，
是噪声分布的协方差矩阵，其中，I∈CQ×Q为单位矩阵。
因此，相应的数据似然函数为：
基于源强矩阵的先验分布模型和数据似然函数可得，源强矩阵的后验分布为：
式中，分别为源强矩阵S的后验分布的均值和协方差，其中，
式中，ΣP∈CQ×Q为传声器声压P★的协方差矩阵：
其中，E(·)表示期望函数。
2D-SBL-CB方法利用证据来估计超参数
和σS，根据全概率公式可得测量声压向量的分布(证据)p(P★)为：
式中，tr(·)表示矩阵的迹；det(·)表示矩阵的行列式。
通过第二类极大似然估计(即，最大化证据)，可获得超参数σS的超参数估计目标函数：
式中，
为计算得到的网格点声压的方差，是声压幅值的平方。
由于上式所示最大化问题的目标函数非凸，因此，通过对目标函数微分逐步更新迭代获得参数σSg的近似估计：
式中，
表示第γ次迭代过程中(θGg,φGg)方向的声源能量估计(方差)；(·)H表示共轭转置。
若噪声信息已知，则2D-SBL-CB方法的计算步骤至此结束；若噪声信息未知，则需执行迭代计算估计噪声方差
为计算模型参数方差提供基础，具体为：
初始化
在第γ次迭代中，超参数
可通过随机极大似然获得：
式中，
为感知矩阵A(ΩG)中索引为
的列向量构成的矩阵；
为
的伪逆。
更新η(γ)：
式中，η(γ)为第η(γ)次迭代中的常数。
在本实施例中，所述预定的收敛条件包括当γ≥γmax或η(γ)<10-3时，迭代停止。
表1为二维稀疏贝叶斯学习压缩波束形成方法(2D-SBL-CB)的算法流程：
表1
上述输出的超参数
即为波束形成的功率输出，而所有网格点声压幅值可由
获得；相比于凸松弛方法和贪婪算法，稀疏贝叶斯学习方法无显式的稀疏约束，而是通过单独地缩放每个网格点源强对应方差来隐含稀疏促进。
需要说明的是，由于直接得到
还是有很多网格点有小的声源，这些小的声源实际上不是真正的声源，而是噪声干扰引起的虚假声源，因此，需要采用截断去除方法去掉这些小的声源，只留下为数不多的真正的大声源(即，将小的虚假声源源强置0，以使后续步骤无需参与运算)，此时的A变成了
以使声源成像问题方程由欠定变成了超定了，因此，就可以采用
求解获得结果。
在一个实施例中，本实施例还包括选用平面传声器阵列，采用本发明实施例提供的二维稀疏贝叶斯学习压缩波束形成方法进行声源成像仿真分析，成像时所需稀疏度或信噪比均需被准确估计。
为验证本实施例提出的二维稀疏贝叶斯学习压缩波束形成方法，本实施例在半消声室中进行了扬声器声源成像试验，试验布局如图3所示，试验时采用孔径0.65m的36通道
扇形轮阵列；声源由两个扬声器及它们的地面反射形成的镜像声源组成，它们的坐标分别约为(1.25，0，3)m、(1.25，-2.4，3)m、(-1.25，0，3)m及(-1.25，-2.4，3)m，对应的DOA分别约为(22.62°,360°)、(42.05°,297.51°)、(42.05°,242.49°)和(22.62°,180°)，依次标记为声源1至声源4。
采用2D-OMP-CB(二维正交匹配追踪压缩波束形成方法)、2D-IR l1-CB(二维迭代重加权l1最小化压缩波束形成方法)和本实施例提出的2D-SBL-CB(二维稀疏贝叶斯学习压缩波束形成方法)进行声源成像，计算时θ和
方向均以5°间隔进行网格划分，考虑传声器的频响失配，将信噪比估计为20dB，稀疏度估计为4，图4分别展示了三种方法在3000Hz和1000Hz下的声源成像云图，表2为三种方法在3000Hz和1000Hz下分别获得的平均DOA估计偏差MAEDOA，表2如下所示：
表2
从表2中可以看出，当频率为3000Hz时，三种方法均能成功分离四个声源，但其中2D-IRl1-CB方法(图4(c))估计声源为真实声源附近多个网格点，能量分散，且存在较多虚假源干扰真实声源的成像；同时，结合表2所示各方法的声源DOA估计偏差可以看出：2D-OMP-CB(图4(a))和2D-IRl1-CB(图4(c))方法获得的声源DOA估计偏差明显高于2D-SBL-CB(图4(e))方法。
当频率为1000Hz时，2D-OMP-CB(图4(b))无法成功分离声源2和声源3，表明此时声源2和声源3间夹角已小于2D-OMP-CB方法在1000Hz时对应的空间分辨能力；而2D-IRl1-CB(图4(d))方法仍存在能量分散和虚假源干扰情况；2D-SBL-CB(图4(f))方法却仍能分离并准确定位各声源，且其DOA估计偏差仍小于已有的2D-OMP-CB和2D-IRl1-CB方法，这表明，2D-SBL-CB方法具有较高的空间分辨能力，具有优于已有的2D-OMP-CB和2D-IRl1-CB方法的声源成像性能。
本发明实施例提供了一种声源成像方法，所述方法利用二维稀疏贝叶斯学习压缩波束形成方法求解严重欠定的声源成像问题模型，其中，二维稀疏贝叶斯学习压缩波束形成方法使用分层的两级贝叶斯推理求解声源成像问题模型，从而自适应地从感知矩阵中学习超参数，以获得稳健的稀疏估计。本实施例的优点在于，与已有的凸松弛方法和贪婪算法相比，二维稀疏贝叶斯学习压缩波束形成方法具有更好的估计精度和计算效率，能够实现稳健且准确的源强估计。
需要说明的是，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本申请实施例的实施过程构成任何限定。
在一个实施例中，如图5所示，本发明实施例提供了一种声源成像系统，所述系统包括：
测量模型建立模块101，用于以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；
成像问题建立模块102，用于将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；
成像问题求解模块103，用于利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
在一个实施例中，所述声源成像问题模型具体为：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝[a(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集，L为快拍总数。
关于一种声源成像系统的具体限定可以参见上述对于一种声源成像方法的限定，此处不再赘述。本领域普通技术人员可以意识到，结合本申请所公开的实施例描述的各个模块和步骤，能够以硬件、软件或者两者结合来实现。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
本实施例提供了一种声源成像系统，所述系统通过测量模型建立模块和成像问题建立模块将声源成像问题转变为声源成像问题模型；所述成像问题求解模块通过二维稀疏贝叶斯学习压缩波束形成方法对所述声源成像问题模型进行求解，获取源强估计；与现有技术相比，本申请提出的二维稀疏贝叶斯学习压缩波束形成方法兼顾高计算效率、强空间分辨能力和强抗噪声干扰能力，且其受网格间距影响很小，声源成像性能稳健。
图6是本发明实施例提供的一种计算机设备，包括存储器、处理器和收发器，它们之间通过总线连接；存储器用于存储一组计算机程序指令和数据，并可以将存储的数据传输给处理器，处理器可以执行存储器存储的程序指令，以执行上述方法的步骤。
其中，存储器可以包括易失性存储器或非易失性存储器，或可包括易失性和非易失性存储器两者；处理器可以是中央处理器、微处理器、特定应用集成电路、可编程逻辑器件或其组合。通过示例性但不是限制性说明，上述可编程逻辑器件可以是复杂可编程逻辑器件、现场可编程逻辑门阵列、通用阵列逻辑或其任意组合。
另外，存储器可以是物理上独立的单元，也可以与处理器集成在一起。
本领域普通技术人员可以理解，图6中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有相同的部件布置。
在一个实施例中，本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述方法的步骤。
本发明实施例提供的一种声源成像方法、系统、设备及存储介质，其一种声源成像方法将二维稀疏贝叶斯学习压缩波束形成方法应用于求解声源成像问题模型，能够兼顾声源成像性能稳健性和高计算效率，且其受网格间距影响很小，声源成像性能稳健；相比于已有的2D-OMP-CB方法和2D-IR l1-CB方法，本发明实施例提出的二维稀疏贝叶斯学习压缩波束形成方法综合性能最优。
在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用软件实现时，可以全部或部分地以计算机程序产品的形式实现。所述计算机程序产品包括一个或多个计算机指令。在计算机上加载和执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输。所述计算机可读存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，DVD)、或者半导体介质(例如SSD)等。
本领域技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。
以上所述实施例仅表达了本申请的几种优选实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和替换，这些改进和替换也应视为本申请的保护范围。因此，本申请专利的保护范围应以所述权利要求的保护范围为准。
Claims (10)
1.一种声源成像方法，其特征在于，包括以下步骤：
以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；
将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；
利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
2.如权利要求1所述的一种声源成像方法，其特征在于，所述将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型的步骤包括：
压缩波束形成将平面传声器阵列测量模型中的目标声源区域沿声源仰角和声源方位角两个方向，以预设间隔离散为若干个固定网格点，且所述固定网格点的数量大于传声器的数量；
假设声源位置估计与固定网格点重合，建立声源成像问题模型。
3.如权利要求2所述的一种声源成像方法，其特征在于，所述声源成像问题模型具体为：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝[a(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集，L为快拍总数。
4.如权利要求3所述的一种声源成像方法，其特征在于，所述利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出包括：
假设所述源强矩阵满足复高斯分布，建立源强矩阵先验分布模型和噪声先验模型；
基于全概率公式，根据所述源强矩阵先验分布模型和所述噪声先验模型得到测量声压向量分布；
基于第二类极大似然估计，根据所述测量声压向量分布得到超参数估计目标函数；
对所述超参数估计目标函数微分逐步更新迭代，得到迭代过程中声源位置估计处的声源能量估计；
若噪声信息已知，则得到波束形成的功率输出；若噪声信息未知，则通过迭代计算噪声估计方差，直至达到预定的收敛条件，则停止迭代。
5.如权利要求4所述的一种声源成像方法，其特征在于，所述超参数估计目标函数具体为：
式中，σSg为源强矩阵的第g个元素的方差；p(P★)为测量声压向量的分布；tr(·)为矩阵的迹；det(·)为矩阵的行列式；ΣP∈CQ×Q为传声器声压P★的协方差矩阵。
6.如权利要求1所述的一种声源成像方法，其特征在于，所述方法还包括：利用所述二维稀疏贝叶斯学习压缩波束形成方法进行声源成像性能仿真分析。
7.一种声源成像系统，其特征在于，所述系统包括：
测量模型建立模块，用于以由若干传声器组成的平面阵列的中心为原点建立笛卡尔坐标系，得到平面传声器阵列测量模型；
成像问题建立模块，用于将所述平面传声器阵列测量模型中的目标声源区域划分为若干个固定网格点，并建立声源成像问题模型；
成像问题求解模块，用于利用预先建立的二维稀疏贝叶斯学习压缩波束形成方法求解所述声源成像问题模型，得到波束形成的功率输出。
8.如权利要求7所述的一种声源成像系统，其特征在于，所述声源成像问题模型具体为：
P★＝A(ΩG)*S+N
式中，P★∈CQ×L为平面传声器阵列中Q个传声器在所有L个快拍下的测量声压信号构成的矩阵；A(ΩG)＝pa(θG1,φG1),…,a(θGg,φGg),...,a(θGG,φGG)]∈CQ×G为感知矩阵，其中，a(θGg,φGg)为第g个固定网格点与所有传声器之间的传递函数，θ为声源入射方向与z轴的夹角，φ为声源入射方向在xoy平面上的投影与x轴的夹角；
为所有快拍下声源源强分布矢量构成的源强矩阵，其中，sg为第g个固定网格点的所有快拍下的源强构成的源强分布矢量；N∈CQ×L为测量噪声干扰信号构成的噪声矩阵，C为复数集，L为快拍总数。
9.一种计算机设备，其特征在于：包括处理器和存储器，所述处理器与所述存储器相连，所述存储器用于存储计算机程序，所述处理器用于执行所述存储器中存储的计算机程序，以使得所述计算机设备执行如权利要求1至6中任一项所述的方法。
10.一种计算机可读存储介质，其特征在于：所述计算机可读存储介质中存储有计算机程序，当所述计算机程序被运行时，实现如权利要求1至6任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210932721.0A CN115291169A (zh) | 2022-08-04 | 2022-08-04 | 一种声源成像方法、系统、设备及存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210932721.0A CN115291169A (zh) | 2022-08-04 | 2022-08-04 | 一种声源成像方法、系统、设备及存储介质 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN115291169A true CN115291169A (zh) | 2022-11-04 |
Family
ID=83826551
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210932721.0A Pending CN115291169A (zh) | 2022-08-04 | 2022-08-04 | 一种声源成像方法、系统、设备及存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115291169A (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116338574A (zh) * | 2023-04-10 | 2023-06-27 | 哈尔滨工程大学 | 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法 |
| CN116989888A (zh) * | 2023-09-27 | 2023-11-03 | 之江实验室 | 声成像方法、装置、计算机设备和存储介质 |
-
2022
- 2022-08-04 CN CN202210932721.0A patent/CN115291169A/zh active Pending
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116338574A (zh) * | 2023-04-10 | 2023-06-27 | 哈尔滨工程大学 | 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法 |
| CN116338574B (zh) * | 2023-04-10 | 2023-09-19 | 哈尔滨工程大学 | 一种基于匹配波束的稀疏贝叶斯学习水下声源定位方法 |
| CN116989888A (zh) * | 2023-09-27 | 2023-11-03 | 之江实验室 | 声成像方法、装置、计算机设备和存储介质 |
| CN116989888B (zh) * | 2023-09-27 | 2024-03-12 | 之江实验室 | 声成像方法、装置、计算机设备和存储介质 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Huang et al. | Robust adaptive beamforming with a novel interference-plus-noise covariance matrix reconstruction method | |
| CN115291169A (zh) | 一种声源成像方法、系统、设备及存储介质 | |
| CN108845325B (zh) | 拖曳线列阵声纳子阵误差失配估计方法 | |
| CN110208735A (zh) | 一种基于稀疏贝叶斯学习的相干信号doa估计方法 | |
| Liu et al. | Augmented subspace MUSIC method for DOA estimation using acoustic vector sensor array | |
| CN107015066B (zh) | 一种基于稀疏贝叶斯学习的天线阵列故障诊断方法 | |
| CN109901101A (zh) | 基于电磁矢量传感器互质阵列相干信号到达角估计方法 | |
| CN111812581B (zh) | 基于原子范数的球面阵列声源波达方向估计方法 | |
| CN111273229A (zh) | 基于低秩矩阵重建的水声宽频散射源的定位方法 | |
| CN109298382A (zh) | 一种基于期望极大算法的非均匀直线阵波达方向角估计方法 | |
| CN112766304A (zh) | 一种基于稀疏贝叶斯学习的机动阵列方位估计方法 | |
| CN111142062A (zh) | 一种利用Toeplitz特性的无栅格目标波达方向估计方法 | |
| CN108614235B (zh) | 一种多鸽群信息交互的单快拍测向方法 | |
| CN112731280B (zh) | 互质阵列混合噪声环境下的esprit-doa估计方法 | |
| CN112763972B (zh) | 基于稀疏表示的双平行线阵二维doa估计方法及计算设备 | |
| CN111896929A (zh) | 非均匀mimo雷达的dod/doa估计算法 | |
| CN115267673B (zh) | 考虑重建网格偏移的稀疏声源成像方法、系统 | |
| Nascimento et al. | Acoustic imaging using the Kronecker array transform | |
| CN113381793B (zh) | 一种面向相干信源估计的无网格波达方向估计方法 | |
| CN114184999B (zh) | 一种互耦小孔径阵列的生成式模型处理方法 | |
| CN115130504A (zh) | 基于稀疏贝叶斯学习的稳健波束形成方法 | |
| Rajani et al. | Direction of arrival estimation by using artificial neural networks | |
| Tan et al. | An iterative adaptive dictionary learning approach for multiple snapshot DOA estimation | |
| CN110824484B (zh) | 一种基于恒模算法的阵元位置估计方法 | |
| Wu et al. | Direction finding with the sensors' gains suffering bayesian uncertainty—hybrid CRB and MAP estimation |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |