平码五不中公式规律
  • / 17
  • 下载费用:30 金币  

实心球阵列三维声源识别的快速反卷积方法.pdf

关 键 ?#21097;?/dt>
实心球 阵列 三维 声源 识别 快速 卷积 方法
  专利查询网所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
摘要
申请专利号:

CN201610878276.9

申请日:

2016.10.08

公开号:

CN106483503A

公开日:

2017.03.08

当前法律状态:

授权

?#34892;?#24615;:

有权

法?#19978;?#24773;: 授权|||实质审查的生效IPC(主分类):G01S 5/18申请日:20161008|||公开
IPC分类号: G01S5/18 主分类号: G01S5/18
申请人: 重庆大学
发明人: 杨洋; 褚志刚; 余立超; 陈涛
地址: 400044 重庆市沙坪坝区沙正街174号
优?#28909;ǎ?/td>
专利代理机构: 重庆大学专利中心 50201 代理人: 唐开平
PDF完整版下载: PDF下载
法律状态
申请(专利)号:

CN201610878276.9

授权公告号:

||||||

法律状态公告日:

2018.10.19|||2017.04.05|||2017.03.08

法律状态类型:

授权|||实质审查的生效|||公开

摘要

本发明公开了一种实心球阵列三维声源识别的快速反卷积方法,包括以下步骤:步骤1、利用声压输出公式计算各个聚焦网格点的声压输出量;步骤2、利用计算得到的各点输出量构造输出矩阵;步骤3、利用PSF函数计算中心聚焦点处声源的PSF矩阵步骤4、迭代求解声源强度分布矩阵本发明在步骤4利用空间转移不变的性质,无需计算庞大的A矩阵,仅计算小维数矩阵由变形得对进行傅里叶变换,基于FFT加速,使计算时间大幅度缩短,提高了计算效?#21097;?#24182;且能够保持较好的空间分辨?#21097;?#20934;确定位各声源。

权利要求书

1.实心球阵列三维声源识别的快速反卷积方法,其特征是,包括以下步骤:
步骤1、利用下面公式计算各个聚焦网格点的声压输出量
<mrow> <mi>W</mi> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>N</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mn>0</mn> </mrow> <mi>N</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <msup> <mi>m</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </mrow> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </munderover> <msup> <mi>&alpha;</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msubsup> <mi>CoY</mi> <mrow> <msup> <mi>nn</mi> <mo>&prime;</mo> </msup> </mrow> <mrow> <msup> <mi>mm</mi> <mo>&prime;</mo> </msup> </mrow> </msubsup> <mo>)</mo> </mrow> <msup> <mi>&alpha;</mi> <mo>*</mo> </msup> <mfrac> <mrow> <msubsup> <mi>Y</mi> <mi>n</mi> <mi>m</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>R</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <mi>k</mi> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <msubsup> <mi>Y</mi> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> <msup> <mi>m</mi> <mo>&prime;</mo> </msup> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>R</mi> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <mi>k</mi> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>*</mo> </msup> </mrow>
式中,三维空间内?#25105;?#20301;置用(r,Ω)表示,r表示所描述位置与原点间的距离,Ω=(θ,
φ)表示所描述位置的方向,?#20219;?#25152;描述位置的仰角、φ为所描述位置的方位角;式中,(rf,
Ωf)为聚焦点的位置坐标;
k=2πf/c为声波的波数,f为声波频?#21097;琧为声速;
α=[α1,α2,Λ,αq,Λ,αQ]T为传声器计权系数组成的列向量;
上标“T”和“*”分别表示转置运算和?#26597;?#36816;算;
n、m、n'、m'均为球谐函数?#29366;危琋为球谐函数?#29366;?#30340;截断长度;
均为Ωf方向的球谐函数;
Rn(krf,ka)、Rn'(krf,ka)均为聚焦径向函数;
为各传声器球谐函数的互谱矩阵;
C为互谱矩阵;
“ο”表示Hadamard积运算;
步骤2、利用计算得到的各点输出量构造矩阵
记gθ=1,2,Λ,G?#21462;φ=1,2,Λ,Gφ分别为聚焦点在仰角、方位角方向的索引,?#24067;蔥小?#21015;
索引,G?#20219;?#32858;焦点在仰角方向的索引最大值、Gφ为聚焦点方位角方向的索引最大值G=GθGφ
为聚焦点总数。则利用步骤1计算得到的各点输出量构造出G?#21462;罣φ维SHB输出矩阵
步骤3、利用PSF函数计算中心聚焦点处声源的PSF矩阵
PSF函数定义为声源识别算法对单位强度单极子点声源的响应,实心球SHB的PSF函数
为:
<mrow> <mi>p</mi> <mi>s</mi> <mi>f</mi> <mrow> <mo>(</mo> <mo>(</mo> <mrow> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> </mrow> <mo>)</mo> <mo>|</mo> <mo>(</mo> <mrow> <msub> <mi>kr</mi> <mn>0</mn> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mn>0</mn> </msub> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>N</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mn>0</mn> </mrow> <mi>N</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <msup> <mi>m</mi> <mo>&prime;</mo> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </mrow> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </munderover> <msup> <mi>&alpha;</mi> <mi>T</mi> </msup> <mrow> <mo>(</mo> <mo>(</mo> <mrow> <mi>t</mi> <mrow> <mo>(</mo> <mrow> <msub> <mi>kr</mi> <mn>0</mn> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mn>0</mn> </msub> </mrow> <mo>)</mo> </mrow> <msup> <mi>t</mi> <mi>H</mi> </msup> <mrow> <mo>(</mo> <mrow> <msub> <mi>kr</mi> <mn>0</mn> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mn>0</mn> </msub> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <msubsup> <mi>oY</mi> <mrow> <msup> <mi>nn</mi> <mo>&prime;</mo> </msup> </mrow> <mrow> <msup> <mi>mm</mi> <mo>&prime;</mo> </msup> </mrow> </msubsup> <mo>)</mo> </mrow> <msup> <mi>&alpha;</mi> <mo>*</mo> </msup> <mfrac> <mrow> <msubsup> <mi>Y</mi> <mi>n</mi> <mi>m</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>R</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <mi>k</mi> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msup> <mrow> <mo>(</mo> <mfrac> <mrow> <msubsup> <mi>Y</mi> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> <msup> <mi>m</mi> <mo>&prime;</mo> </msup> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mi>f</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>R</mi> <msup> <mi>n</mi> <mo>&prime;</mo> </msup> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mi>f</mi> </msub> <mo>,</mo> <mi>k</mi> <mi>a</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>)</mo> </mrow> <mo>*</mo> </msup> </mrow>
式中,t(kr0,Ω0)=[t1(kr0,Ω0),t2(kr0,Ω0),Λ,tq(kr0,Ω0),Λ,tQ(kr0,Ω0)]T为声
场传递函数列向量;上标“H”表示转置?#26597;?#36816;算;
令tq(kr0,Ω0)表示声源到q号传声器的声场传递函数,根据散射理论,
<mrow> <msub> <mi>t</mi> <mi>q</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mn>0</mn> </msub> <mo>,</mo> <msub> <mi>&Omega;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mo>&equiv;</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <mn>0</mn> </mrow> <mi>&infin;</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mo>-</mo> <mi>n</mi> </mrow> <mi>n</mi> </munderover> <msub> <mi>R</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>kr</mi> <mn>0</mn> </msub> <mo>,</mo> <mi>k</mi> <mi>a</mi> <mo>)</mo> </mrow> <msubsup> <mi>Y</mi> <mi>n</mi> <msup> <mi>m</mi> <mo>*</mo> </msup> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <msubsup> <mi>Y</mi> <mi>n</mi> <mi>m</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&Omega;</mi> <mi>q</mi> </msub> <mo>)</mo> </mrow> </mrow>
式中,为Ω0方向的球谐函数,为Ωq方向的球谐函数,Rn(kr0,ka)为声信
号传播径向函数;实际应用中,用截断长度N替代式中的∞来计算PSF;
记为中心聚焦点处声源的PSF矩阵,即:
<mrow> <mover> <mi>A</mi> <mo>~</mo> </mover> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>psf</mi> <mn>11</mn> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mn>12</mn> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mn>13</mn> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>1</mn> <msub> <mi>g</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>1</mn> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>1</mn> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>psf</mi> <mn>21</mn> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mn>22</mn> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mn>23</mn> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>2</mn> <msub> <mi>g</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>2</mn> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mn>2</mn> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>2</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>3</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <msub> <mi>g</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mi>M</mi> </mtd> <mtd> <mi>M</mi> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mn>2</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mn>3</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <msub> <mi>g</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <mo>(</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>1</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>2</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>g</mi> <mi>&theta;</mi> </msub> <mn>3</mn> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <msub> <mi>g</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> <mtd> <mi>&Lambda;</mi> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> <mo>-</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> </msub> </mrow> </mtd> <mtd> <mrow> <msub> <mi>psf</mi> <mrow> <msub> <mi>G</mi> <mi>&theta;</mi> </msub> <msub> <mi>G</mi> <mi>&phi;</mi> </msub> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>
式中,θmin为聚
焦声源面的最小仰角,Δ?#21462;ⅵう?#20998;别为仰角和方位角间隔步骤,Ωc=(θc,φc)为中心聚焦
点的方向;
步骤4、初始化声源强度分布矩阵设定迭代次数L,迭代求解利用第l次迭代
结果进行第l+1次迭代求解的步骤如下:
步骤1)、计算残差矩阵:
式中,均为G?#21462;罣φ维矩阵,由中元素向上、向左分别移动gθc-1、gφc-1?#25442;?br />得,gθc、gφc为中心聚焦点的?#23567;?#21015;索引,“F”、“F-1”分别表示二维傅里叶正、逆变换;
步骤2)、计算矩阵:
<mrow> <msup> <mover> <mi>&beta;</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> </msup> <mo>=</mo> <mo>-</mo> <msup> <mi>F</mi> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <mrow> <mo>(</mo> <mi>F</mi> <mo>(</mo> <msubsup> <mover> <mi>r</mi> <mo>~</mo> </mover> <mi>e</mi> <mrow> <mo>(</mo> <mi>l</mi> <mo>)</mo> </mrow> </msubsup> <mo>)</mo> <mi>o</mi> <msup> <mrow> <mo>(</mo> <mrow> <mi>F</mi> <mrow> <mo>(</mo> <msub> <mover> <mi>A</mi> <mo>~</mo> </mover> <mi>C</mi> </msub> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mo>*</mo> </msup> <mo>)</mo> </mrow> </mrow>
步骤3)、获得矩阵

式中,β(l)(kr0,Ω0)、分别为的元素,S(l)(kr0,Ω0)为的元素;
步骤4)、计算辅助矩阵:

步骤5)、计算步长:
式中,g(l)、分别为中元素形成的向量,“·”表示内积运算;
步骤6)、确定声源强度分布矩阵:
式中,“Ρ+(·)”表示用0替代括号内矩阵的负元素。

说明书

实心球阵列三维声源识别的快速反卷积方法

?#38469;?#39046;域

本发明涉及声源识别的?#38469;?#39046;域。

背景?#38469;?br />

基于传声器阵列测量的波束形成已在声源识别领域占据不可缺少的地位,二维平
面和三维实心球都是普遍采用的阵列形式。二维平面阵列适宜于识别阵列前方特定张角区
域内的声源,常用算法为延迟求和(DAS);三维实心球阵?#24515;?#22815;全方位识别三维腔?#19968;?#22659;内
的声源,常用算法为球谐函数波束形成(SHB)。不论DAS还是SHB,输出结果均可看作各声源
强度与对应点传播函数(point spread function,PSF函数)的乘积的和,实际应用中,传声
器离散采样等因素使两种算法的PSF函数均无法等于理想的δ函数,不仅在真实声源位置输
出宽“主瓣?#20445;?#36824;在非声源位置输出高“旁瓣?#20445;?#26368;终致使声源识别结果空间分辨率差,且?#34892;?br />动态范围低。

“Deconvolution for three-dimensional acoustic source identification
based on spherical harmonics beamforming?#20445;琙.G.Chu,Y.Yang,Y.S.He.Journal of
Sound and Vibration,Volume 344,Pages 484-502,26May 2015(基于球谐函数波束形成
的三维反卷积声源识别,褚志刚,杨洋,贺?#23452;桑琂ournal of Sound and Vibration,344卷,
第484~502页,2015年5月26日)介绍了常见的实心球阵列三维声源识别反卷积方法,如非
负最小二乘(Non-Negative Least Squares,NNLS),NNLS能获得干净明辨的三维声源?#19978;瘢?br />但是它的计算时间很长、效率很低,不能满足波束形成实时?#19978;?#30340;需求。

发明内容

针对现有?#38469;?#20013;存在的?#38469;?#38382;题,本发明所要解决的?#38469;?#38382;题就是提供一种实心
球阵列三维声源识别的快速反卷积方法,它能够减少计算时间、提高计算效?#21097;?#36798;到实时成
像的要求。

本发明所要解决的?#38469;?#38382;题是通过这样的?#38469;?#26041;案实现的,它包括以下步骤:

步骤1、利用下面公式计算各个聚焦网格点的声压输出量


式中,三维空间内?#25105;?#20301;置用(r,Ω)表示,r表示所描述位置与原点间的距离,Ω
=(θ,φ)表示所描述位置的方向,?#20219;?#25152;描述位置的仰角、φ为所描述位置的方位角;式中,
(rf,Ωf)为聚焦点的位置坐标;

k=2πf/c为声波的波数,f为声波频?#21097;琧为声速;

α=[α1,α2,Λ,αq,Λ,αQ]T为传声器计权系数组成的列向量;

上标“T”和“*”分别表示转置运算和?#26597;?#36816;算;

n、m、n'、m'均为球谐函数?#29366;危琋为球谐函数?#29366;?#30340;截断长度;

均为Ωf方向的球谐函数;

Rn(krf,ka)、Rn'(krf,ka)均为聚焦径向函数;

为各传声器球谐函数的互谱矩阵;

C为互谱矩阵;

“ο”表示Hadamard积运算;

步骤2、利用计算得到的各点输出量构造矩阵

记gθ=1,2,Λ,G?#21462;φ=1,2,Λ,Gφ分别为聚焦点在仰角、方位角方向的索引,?#24067;?br />?#23567;?#21015;索引,G?#20219;?#32858;焦点在仰角方向的索引最大值、Gφ为聚焦点方位角方向的索引最大值G
=GθGφ为聚焦点总数。则利用步骤1计算得到的各点输出量构造出G?#21462;罣φ维SHB输出矩阵

步骤3、利用PSF函数计算中心聚焦点处声源的PSF矩阵

PSF函数定义为声源识别算法对单位强度单极子点声源的响应,实心球SHB的PSF
函数为:



式中,t(kr0,Ω0)=[t1(kr0,Ω0),t2(kr0,Ω0),Λ,tq(kr0,Ω0),Λ,tQ(kr0,Ω0)]T
为声场传递函数列向量;上标“H”表示转置?#26597;?#36816;算;

令tq(kr0,Ω0)表示声源到q号传声器的声场传递函数,根据散射理论,


式中,为Ω0方向的球谐函数,为Ωq方向的球谐函数,Rn(kr0,ka)为
声信号传播径向函数。实际应用中,用截断长度N替代式中的∞来计算PSF;

记为中心聚焦点处声源的PSF矩阵,即:


式中,θmin
为聚焦声源面的最小仰角,Δ?#21462;ⅵう?#20998;别为仰角和方位角间隔步骤,Ωc=(θc,φc)为中心
聚焦点的方向;

步骤4、初始化声源强度分布矩阵设定迭代次数L,迭代求解利用第l次
迭代结果进行第l+1次迭代求解的步骤如下:

步骤1)、计算残差矩阵:

式中,均为G?#21462;罣φ维矩阵,由中元素向上、向左分别移动gθc-1、gφc-1
?#25442;?#24471;,gθc、gφc为中心聚焦点的?#23567;?#21015;索引,“F”、“F-1”分别表示二维傅里叶正、逆变换;

步骤2)、计算矩阵:


步骤3)、获得矩阵


式中,β(l)(kr0,Ω0)、分别为的元素,S(l)(kr0,Ω0)为的元
素;

步骤4)、计算辅助矩阵:


步骤5)、计算步长:

式中,g(l)、分别为中元素形成的向量,“·”表示内积运算;

步骤6)、确定声源强度分布矩阵:

式中,“Ρ+(·)”表示用0替代括号内矩阵的负元素。

由于本发明利用PSF空间转移不变性的性质,省去了现有NNLS方法中庞大的A矩阵
的计算和部分乘法运算,仅计算小维数矩阵由变形得对进行傅里叶变换,能基
于FFT加速,使计算时间大幅度减少,提高了计算效?#21097;?#24182;且能够保持较好的空间分辨?#21097;?#20934;
确定位各声源。

附图说明

本发明的附图说明如下:

图1为本发明的球阵列波束形成的坐标系;

图2不同位置声源的PSF?#19978;?#22270;(频?#21097;?000Hz)

图3 PSF标准差?#19978;?#22270;(频?#21097;?000Hz)

图4为聚焦点和扩展结点分布示意图;

图5为声波2000Hz时不同位置处声源的仿真识别?#19978;?#22270;;

图6为应用本发明试验时的声源与传声器实心球阵列的布置图;

图7为声波2000Hz时不同位置处扬声器声源的试验识别?#19978;?#22270;。

具体实施方式

下面结合附图和实施例对本发明作进一步说明:

本发明的原理?#39057;?#36807;程是:

1、确定实心球SHB的PSF函数

本发明的球阵列波束形成的坐标系如图1所示,球原点位于阵列中心,三维空间内
?#25105;?#20301;置坐标(r,Ω),r表示所描述位置与原点间的距离,Ω=(θ,φ)表示所描述位置的方
向,?#21462;ⅵ?#20998;别为其仰角、方位角。

图1中,符号“●”代表嵌在实心球表面的传声器,记a为球半径,Q为传声器总数,q
=1,2,…,Q为传声器索引,(a,Ωq)为q号传声器的位置坐标,球阵列利用这些传声器空间
采样声压信号。记f为声波频?#21097;琧为声速,k=2πf/c为声波的波数,(rf,Ωf)为聚焦点的位置
坐标,根据背景?#38469;?#25152;引用的文献,声压输出量W(krf,Ωf)的表达式为:


式(1)中,上标“T”和“*”分别表示转置运算和?#26597;?#36816;算,n、m、n'、m'均为球谐函数
?#29366;危琋为球谐函数?#29366;?#30340;截断长度,均为Ωf方向的球谐函数,Rn(krf,
ka)、Rn'(krf,ka)均为聚焦径向函数,α=[α1,α2,Λ,αq,Λ,αQ]T为传声器计权系数组成的列
向量,为各传声器球谐函数的互谱矩阵,C表示传声器信号互谱矩阵,“ο”表示Hadamard
积运算。N的取值?#35272;?#20110;波数k、球半径a及阵列设计截断长度ND:


用(r0,Ω0)表示声源的位置坐标,tq(kr0,Ω0)表示声源到q号传声器的声场传递函
数,根据散射理论:


SHB的PSF函数为:



式中,t(kr0,Ω0)=[t1(kr0,Ω0),t2(kr0,Ω0),Λ,tq(kr0,Ω0),Λ,tQ(kr0,Ω0)]T
为声场传递函数列向量,上标“H”表示转置?#26597;?#36816;算。

2、PSF函数空间转移不变性分析

定义函数psfs(krf,kr0,ΔΩ)为:

psfs(krf,kr0,ΔΩ)=psfs(krf,kr0,Ωf-Ωc)=psf((krf,Ωf)(kr0,Ωc)) (5)

式(5)中,Ωc=(θc,φc)为中心聚焦点的方向,若所有声源的PSF函数均满足

psf((krf,Ωf)(kr0,Ω0))=psfs(krf,kr0,Ωf-Ω0) (6)

则称PSF函数空间转移不变,本发明所建立的反卷积方法以该性质为前提。

为分析PSF的空间转移不变性,以图1所示半径为0.0975m的36通道实心球阵列为
例进行仿真模拟。设聚焦声源面为与阵列同心且半径等于1m的球面,仰角从0°到180°,方位
角从0°到360°,间隔均为5°,共37×73个聚焦点,图2给出了2000Hz时不同聚焦点处声源的
PSF?#19978;?#22270;,?#19978;?#37327;为参考2.0×10-5缩放后的dB值。各图中均在真实声源位置输出主瓣,在
非声源位置输出旁瓣,这些主瓣形状及旁瓣分布不全相同,说明式(6)所示空间转移不变性
不绝对成立。进一步,定义PSF标准差为:



其中,K为部分聚焦点位置坐标组成的集合,GK为集合K的元素总数,与Ωf=(θf,
φf)一样,Ωf'=(θf',φf')也表示聚焦方向,二者间存在如下关系:

θf'=θf-θc+θ0 (8)


这里,φf'的取值考虑了方位角的循环性。对图2所示聚焦声源面,K为:

K={(rf,Ωf=(θf,φf))|rf=1m,max(0,θc-θ0)≤θf≤min(π,π+θc-θ0),0≤φf≤2
π} (10)该标准差越小,意味着(r0,Ω0)位置处声源的PSF相对于(r0,Ωc)位置处声源的
PSF的空间转移变化越弱。计算2000Hz时各聚焦点处声源的PSF标准差,所得结果如图3所
示,显见,θc两侧附近区域内,该标准差均较小,可认为式(6)所示空间转移不变性近似成
立。改变声源辐射声波的频?#24335;?#34892;计算,所得结果亦如此。定义最大仰角或最小仰角与中心
仰角的绝对差为仰角张角,根据上述分析,要使反卷积方法的前提(聚焦声源面内所有声源
的PSF均满足式(6)所示空间转移不变性)近似成立,聚焦声源面需采用小仰角张角,对方位
角无要求。

3、定义实心球SHB输出函数

记θmin、θmax分别为聚焦声源面的最小和最大仰角,Δ?#21462;ⅵう?#20998;别为仰角和方位角
间隔,考虑0和2π方位?#24378;?#38388;重合,聚焦声源面的最小和最大方位角分别取为0和2π-Δφ,
沿仰角、方位角方向均等间隔离散聚焦声源面,所得聚焦点分布如图4(a)所示,其中,gθ=
1,2,Λ,G?#21462;φ=1,2,Λ,Gφ分别为聚焦点在仰角、方位角方向的索引,?#24067;蔥小?#21015;索引,G?#20219;?br />聚焦点在仰角方向的索引最大值、Gφ为聚焦点方位角方向的索引最大值,G=GθGφ为聚焦点
总数。假设各聚焦点处均存在声源,?#20918;?#27492;不相干,则互谱矩阵C可写为:


其中,B为声源位置坐标组成的集合,S(kr0,Ω0)为声源强度。再联立式(1)、式(4)
和式(11)得到:


记W为G×1维SHB输出列向量,A为G×G维PSF矩阵,S为G×1维声源强度分布列向
量,根据式(12),输出函数可写为:

W=AS(13)

这里,W、A可分别由式(1)、式(4)计算得出,S为未知量。

4、构建空间转移不变性条件下的输出函数

记为G?#21462;罣φ维声源强度分布矩阵,即:


式(14)中,

记gθc、gφc为中心聚焦点的?#23567;?#21015;索引,图4(a)中的聚焦点向左、右、上、下依次扩展
Gφ-gφc列、gφc-1列、Gθ-gθc?#23567;θc-1行,得到如图4(b)所示的扩散节点分布。记
分别为左侧、?#20063;唷?#19978;方、下方扩展结点上的声源强度分布矩阵,维数分
别为G?#21462;?Gφ-gφc)、G?#21462;?gφc-1)、(Gθ-gθc)×(2Gφ-1)、(gθc-1)×(2Gφ-1),元素布局与类
同。

构建(2Gθ-1)×(2Gφ-1)维矩阵为:


记为中心聚焦点处声源的PSF矩阵,即:



式(16)中,

以为中心180°旋转得对聚焦点(rf,Ωf),通过对齐(rf,Ωf)
位置使与中的元素相对应,对应元素相乘相加得WS(krf,Ωf),即:


式(17)中,B'=B'1YB'2,B'1、B'2分别为部分聚焦点、部分扩展结点位置坐标组成
的集合,“Y”表示并集运算。

记WS为所有聚焦点处输出组成的列向量,若PSF满足式(5)所示空间转移不变性,
那么进一步得到:

W≈WS=ASS (18)

式(18)中,AS称为空间转移不变PSF矩阵,S为待求解的未知量。

5、迭代求解声源强度分布列向量S

根据背景?#38469;?#25152;引用的文献,初始化S(0)=0,NNLS基于第l次迭代结果S(l)进行第l
+1次迭代的具体方案为:

1.计算残差向量:

2.计算函数(“||·||2”表示2范数)关于S的负梯度向量:

3.获得向量


式中,β(l)(kr0,Ω0)为β(l)的元素。

4.计算辅助向量:

5.计算步长:

式中,“·”表示内积运算。

6.确定声源强度分布向量:

式中,“Ρ+(·)”表示用0替代括号内向量的负元素,即向量在非负象限上的欧几
里德投影。

现采用周期边界条件,即

其中,“:”表示从其左侧数?#31181;?#31034;的行(列)到其?#20063;?#25968;?#31181;?#31034;的行(列),若无数字,则表示
所?#34892;?列)。此时,AS是由中元素组成的BCCB矩阵,“Deblurringimages:matrices,
spectra,and filtering?#20445;琍.C.Hansen,J.G.Nagy,D.P.O’Leary,Pages 40-43,Society
for Industrial and Applied Mathematics,2006(“清晰化?#19978;瘢?#30697;阵,波谱和滤波?#20445;?br />P.C.Hansen,J.G.Nagy,D.P.O’Leary,第40~43页,SIAM,2006年)介绍了BCCB矩阵的构成。
故AS的谱分解为:

AS=FHΛF (25)

其中,F为二维酉离散傅里叶变换矩阵,Λ为特征值矩阵,FH为F的转置?#26597;?#30697;阵,
F、FH与?#25105;?#21521;量的乘积可通过二维傅里叶变换获得,无需明确计算F,记X为?#25105;釭?#21462;罣φ维矩
阵,x为X的列堆叠而成的G×1维向量,下列等价关系成立:


上式中表?#38236;?#20215;关系。由式(25)可得:


其中,为AS特征值形成的列向量,第2个“=”是因为联立
式(26)和式(27)可得:


式(28)中,是AS第1列元素形成的矩阵,可由中元素向上、向左分别移动gθc-1、
gφc-1?#25442;?#24471;。基于上述分析,



式(29)中,第2个“=?#24065;?#27714;AS为实矩阵。

上述NNLS迭代的1、2、4中的式(19)、式(20)和式(22)分别改写为:




分别为W、β(l)、g(l)、中元素形成的矩阵。

综合上述原理?#39057;?#36807;程,得本发明的步骤包括:

步骤1、利用下面公式计算各个聚焦网格点的声压输出量


式中,三维空间内?#25105;?#20301;置用(r,Ω)表示,r表示所描述位置与原点间的距离,Ω
=(θ,φ)表示所描述位置的方向,?#20219;?#25152;描述位置的仰角、φ为所描述位置的方位角;式中,
(rf,Ωf)为聚焦点的位置坐标;

k=2πf/c为声波的波数,f为声波频?#21097;琧为声速;

α=[α1,α2,Λ,αq,Λ,αQ]T为传声器计权系数组成的列向量;

上标“T”和“*”分别表示转置运算和?#26597;?#36816;算;

n、m、n'、m'均为球谐函数?#29366;危琋为球谐函数?#29366;?#30340;截断长度;

均为Ωf方向的球谐函数;

Rn(krf,ka)、Rn'(krf,ka)均为聚焦径向函数;

为各传声器球谐函数的互谱矩阵;

C为互谱矩阵;

“ο”表示Hadamard积运算;

步骤2、利用计算得到的各点输出量构造矩阵

记gθ=1,2,Λ,G?#21462;φ=1,2,Λ,Gφ分别为聚焦点在仰角、方位角方向的索引,?#24067;?br />?#23567;?#21015;索引,G?#20219;?#32858;焦点在仰角方向的索引最大值、Gφ为聚焦点方位角方向的索引最大值G
=GθGφ为聚焦点总数。则利用步骤1计算得到的各点输出量构造出G?#21462;罣φ维SHB输出矩阵

步骤3、利用PSF函数计算中心聚焦点处声源的PSF矩阵

PSF函数定义为声源识别算法对单位强度单极子点声源的响应,实心球SHB的PSF
函数为:



式中,t(kr0,Ω0)=[t1(kr0,Ω0),t2(kr0,Ω0),Λ,tq(kr0,Ω0),Λ,tQ(kr0,Ω0)]T
为声场传递函数列向量;上标“H”表示转置?#26597;?#36816;算;

令tq(kr0,Ω0)表示声源到q号传声器的声场传递函数,根据散射理论,


式中,为Ω0方向的球谐函数,为Ωq方向的球谐函数,Rn(kr0,ka)为
声信号传播径向函数;实际应用中,用截断长度N替代式中的∞来计算PSF;

记为中心聚焦点处声源的PSF矩阵,即:


式中,θmin
为聚焦声源面的最小仰角,Δ?#21462;ⅵう?#20998;别为仰角和方位角间隔步骤,Ωc=(θc,φc)为中心
聚焦点的方向;

步骤4、初始化声源强度分布矩阵设定迭代次数L,迭代求解利用第l次
迭代结果进行第l+1次迭代求解的步骤如下:

步骤1)、计算残差矩阵:

式中,均为G?#21462;罣φ维矩阵,由中元素向上、向左分别移动gθc-1、gφc-1
?#25442;?#24471;,gθc、gφc为中心聚焦点的?#23567;?#21015;索引,“F”、“F-1”分别表示二维傅里叶正、逆变换;

步骤2)、计算矩阵:


步骤3)、获得矩阵


式中,β(l)(kr0,Ω0)、分别为的元素,S(l)(kr0,Ω0)为的元
素;

步骤4)、计算辅助矩阵:


步骤5)、计算步长:

式中,g(l)、分别为中元素形成的向量,“·”表示内积运算;

步骤6)、确定声源强度分布矩阵:

式中,“Ρ+(·)”表示用0替代括号内矩阵的负元素。

本方法发明无需计算庞大的A矩阵,仅计算小维数矩阵由变形得对进
行傅里叶变换,能基于FFT(快速傅里叶变换)加速,这是其提高效?#23454;脑?#22240;,为便于区分,将
方法发明称为FFT-NNLS。

实施例1

为验证建立本发明的准确性、对比探究其性能上的提升,进行声源识别仿真模拟。
具体流程为:

1、在特定位置假设具有特定强度辐射特定频率声波的点声源;

2、根据式(3)及式(11)正向计算各传声器接收声压信号的互谱矩阵,这里,用50替
代式(3)中的∞;

3、根据式(2)确定截断长度N;

4、设定聚焦声源面,分别使用NNLS和本发明重构声源强度分布并?#19978;瘛?br />

这里,两种方法被用于识别15°仰角张角区域(仰角:75°~105°,方位角:0°~
360°)内的声源,聚焦声源面被设为与阵列同心且半径等于声源到阵列中心距离的球面,仰
角和方位角间隔均为5°,两种反卷积方法的迭代次数均设为1000。

图3给出了2000Hz时(r0=1m,θ0=90°,φ0=180°)、(r0=1m,θ0=75°,φ0=270°)、
(r0=1m,θ0=90°,φ0=0°)及(r0=1m,θ0=90°,φ0=5°)位置处声源的识别?#19978;?#22270;,为便于
对比,各图中均参考最大输出?#21040;?#34892;dB缩放,显示动态范围均为0~-15dB,同时,每幅图的
上方还列出了以标准声压大小2.0×10-5为参考的最大输出值。图5(a)~(d)中,NNLS在各声
源位置输出窄的主瓣,在非声源位置未输出旁瓣;图5(e)~(h)为本发明的?#19978;?#22270;,从图5
(e)~(h)中对应看出:本发明在各声源位置同样输出窄的主瓣,仅在部分非声源位置输出
极少量的旁瓣,表明本发明同NNLS一样均拥有高的空间分辨率、能?#34892;?#31227;除旁瓣,准确地定
位各声源。

上述仿真模拟在2.5GHz Intel(R)Core(TM)i5-2450M的CPU上进行,完成单频计
算,NNLS的耗时约为43s,本发明仅需约0.7s,可见,与NNLS相比,本发明的计算时间大幅度
减少、计算效率有了巨大的提升。

实施例2

为检验仿真结论的正确性,在消声室内进行实体测试试验。

图6为试验布局图,将稳态信号激励的扬声器作为声源,采用公司、半
径为0.0975m、集成4958型传声器的36通道实心球阵列采样声压信号。各传声器接收的声压
信号经PULSE 3560D型数据采集系统同时采集并传输到PULSE LABSHOP中进行频谱分析,得
声压信号的互谱矩阵,采样频率为16384Hz,信号加?#32791;?#31383;,采用64段平均、66.7%的重叠
?#21097;?#27599;段时长0.25s、对应的频率分辨率为4Hz。进一步,采用MATLAB编制的NNLS和本发明程
序计算各聚焦点的输出量并?#19978;瘛?#36825;里,声源位置、聚焦声源面及反卷积迭代次数的设置均
与仿真模拟时一致。

图7给出了2000Hz时扬声器声源的识别?#19978;?#22270;,图7中的FFT-NNLS为本发明的?#19978;?br />图。其与图5呈现出完全一致的规律,实体测试过程所用计算的时间与仿真结果相同,证明
本方法发明与NNLS相比,计算时间大幅度减少。

关于本文
本文标题:实心球阵列三维声源识别的快速反卷积方法.pdf
链接地址:http://www.pqiex.tw/p-5994937.html
关于我们 - 网站声明 - 网?#38236;?#22270; - 资源地图 - 友情链接 - 网站客服 - 联系我们

[email protected] 2017-2018 zhuanlichaxun.net网站版权所有
经营许可证编号:粤ICP备17046363号-1 
 


收起
展开
平码五不中公式规律 极速快3app 扑克山庄官方网站 竞彩足球混合过关4串1 上海时时乐开奖今天 浙江体彩6+1 今日河南11选5开奖结果 安徽时时彩计划软件手机版下载手机版下载 双色球140期蓝球分析 新浪棋牌官网手机版 网络赌场