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

Abstract

本发明公开了一种实心球阵列三维声源识别的快速反卷积方法，包括以下步骤：步骤1、利用声压输出公式计算各个聚焦网格点的声压输出量；步骤2、利用计算得到的各点输出量构造输出矩阵；步骤3、利用PSF函数计算中心聚焦点处声源的PSF矩阵步骤4、迭代求解声源强度分布矩阵本发明在步骤4利用空间转移不变的性质，无需计算庞大的A矩阵，仅计算小维数矩阵由变形得对

Description

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

技术领域

本发明涉及声源识别的技术领域。

背景技术

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

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

发明内容

针对现有技术中存在的技术问题，本发明所要解决的技术问题就是提供一种实心球阵列三维声源识别的快速反卷积方法，它能够减少计算时间、提高计算效率，达到实时成像的要求。

本发明所要解决的技术问题是通过这样的技术方案实现的，它包括以下步骤：

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

式中，三维空间内任意位置用(r,Ω)表示，r表示所描述位置与原点间的距离，Ω＝(θ,φ)表示所描述位置的方向，θ为所描述位置的仰角、φ为所描述位置的方位角；式中，(r_f,Ω_f)为聚焦点的位置坐标；

k＝2πf/c为声波的波数，f为声波频率，c为声速；

α＝[α₁,α₂,Λ,α_q,Λ,α_Q]^T为传声器计权系数组成的列向量；

上标“_T”和“*”分别表示转置运算和共轭运算；

n、m、n'、m'均为球谐函数阶次，N为球谐函数阶次的截断长度；

均为Ω_f方向的球谐函数；

R_n(kr_f,ka)、R_n'(kr_f,ka)均为聚焦径向函数；

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

C为互谱矩阵；

“ο”表示Hadamard积运算；

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

记g_θ＝1,2,Λ,G_θ、g_φ＝1,2,Λ,G_φ分别为聚焦点在仰角、方位角方向的索引，也即行、列索引，Gθ为聚焦点在仰角方向的索引最大值、G_φ为聚焦点方位角方向的索引最大值G＝G_θG_φ为聚焦点总数。则利用步骤1计算得到的各点输出量构造出G_θ×G_φ维SHB输出矩阵

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

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

式中，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),Λ,t_q(kr₀,Ω₀),Λ,t_Q(kr₀,Ω₀)]^T为声场传递函数列向量；上标“H”表示转置共轭运算；

令t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，根据散射理论，

式中，为Ω₀方向的球谐函数，为Ω_q方向的球谐函数，R_n(kr₀,ka)为声信号传播径向函数。实际应用中，用截断长度N替代式中的∞来计算PSF；

记为中心聚焦点处声源的PSF矩阵，即：

式中，θ_min为聚焦声源面的最小仰角，Δθ、Δφ分别为仰角和方位角间隔步骤，Ω_c＝(θ_c,φ_c)为中心聚焦点的方向；

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

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

式中，均为G_θ×G_φ维矩阵，由中元素向上、向左分别移动g_θc-1、g_φc-1位获得，g_θc、g_φc为中心聚焦点的行、列索引，“F”、“F^-1”分别表示二维傅里叶正、逆变换；

步骤2)、计算矩阵：

步骤3)、获得矩阵

式中，β^(l)(kr₀,Ω₀)、分别为的元素，S^(l)(kr₀,Ω₀)为的元素；

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

步骤5)、计算步长：

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

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

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

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

附图说明

本发明的附图说明如下：

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

图2不同位置声源的PSF成像图(频率：2000Hz)

图3PSF标准差成像图(频率：2000Hz)

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

图5为声波2000Hz时不同位置处声源的仿真识别成像图；

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

图7为声波2000Hz时不同位置处扬声器声源的试验识别成像图。

具体实施方式

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

本发明的原理推导过程是：

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

本发明的球阵列波束形成的坐标系如图1所示，球原点位于阵列中心，三维空间内任意位置坐标(r,Ω)，r表示所描述位置与原点间的距离，Ω＝(θ,φ)表示所描述位置的方向，θ、φ分别为其仰角、方位角。

图1中，符号“●”代表嵌在实心球表面的传声器，记a为球半径，Q为传声器总数，q＝1,2,…,Q为传声器索引，(a,Ω_q)为q号传声器的位置坐标，球阵列利用这些传声器空间采样声压信号。记f为声波频率，c为声速，k＝2πf/c为声波的波数，(r_f,Ω_f)为聚焦点的位置坐标，根据背景技术所引用的文献，声压输出量W(kr_f,Ω_f)的表达式为：

式(1)中，上标“T”和“*”分别表示转置运算和共轭运算，n、m、n'、m'均为球谐函数阶次，N为球谐函数阶次的截断长度，均为Ω_f方向的球谐函数，R_n(kr_f,ka)、R_n'(kr_f,ka)均为聚焦径向函数，α＝[α₁,α₂,Λ,α_q,Λ,α_Q]^T为传声器计权系数组成的列向量，为各传声器球谐函数的互谱矩阵，C表示传声器信号互谱矩阵，“ο”表示Hadamard积运算。N的取值依赖于波数k、球半径a及阵列设计截断长度N_D：

用(r₀,Ω₀)表示声源的位置坐标，t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，根据散射理论：

SHB的PSF函数为：

式中，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),Λ,t_q(kr₀,Ω₀),Λ,t_Q(kr₀,Ω₀)]^T为声场传递函数列向量，上标“H”表示转置共轭运算。

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

定义函数psf_s(kr_f,kr₀,ΔΩ)为：

psf_s(kr_f,kr₀,ΔΩ)＝psf_s(kr_f,kr₀,Ω_f-Ω_c)＝psf((kr_f,Ω_f)(kr₀,Ω_c)) (5)

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

psf((kr_f,Ω_f)(kr₀,Ω₀))＝psf_s(kr_f,kr₀,Ω_f-Ω₀) (6)

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

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

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

θ_f'＝θ_f-θ_c+θ₀ (8)

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

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

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

记θ_min、θ_max分别为聚焦声源面的最小和最大仰角，Δθ、Δφ分别为仰角和方位角间隔，考虑0和2π方位角空间重合，聚焦声源面的最小和最大方位角分别取为0和2π-Δφ，沿仰角、方位角方向均等间隔离散聚焦声源面，所得聚焦点分布如图4(a)所示，其中，g_θ＝1,2,Λ,G_θ、g_φ＝1,2,Λ,G_φ分别为聚焦点在仰角、方位角方向的索引，也即行、列索引，G_θ为聚焦点在仰角方向的索引最大值、G_φ为聚焦点方位角方向的索引最大值，G＝G_θG_φ为聚焦点总数。假设各聚焦点处均存在声源，且彼此不相干，则互谱矩阵C可写为：

其中，B为声源位置坐标组成的集合，S(kr₀,Ω₀)为声源强度。再联立式(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_θ×G_φ维声源强度分布矩阵，即：

式(14)中，

记g_θc、g_φc为中心聚焦点的行、列索引，图4(a)中的聚焦点向左、右、上、下依次扩展G_φ-g_φc列、g_φc-1列、G_θ-g_θc行、g_θc-1行，得到如图4(b)所示的扩散节点分布。记分别为左侧、右侧、上方、下方扩展结点上的声源强度分布矩阵，维数分别为G_θ×(G_φ-g_φc)、G_θ×(g_φc-1)、(G_θ-g_θc)×(2G_φ-1)、(g_θc-1)×(2G_φ-1)，元素布局与类同。

构建(2G_θ-1)×(2G_φ-1)维矩阵为：

记为中心聚焦点处声源的PSF矩阵，即：

式(16)中，

以为中心180°旋转得对聚焦点(r_f,Ω_f)，通过对齐(r_f,Ω_f)位置使与中的元素相对应，对应元素相乘相加得W_S(kr_f,Ω_f)，即：

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

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

W≈W_S＝A_SS (18)

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

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

根据背景技术所引用的文献，初始化S⁽⁰⁾＝0，NNLS基于第l次迭代结果S^(l)进行第l+1次迭代的具体方案为：

1.计算残差向量：

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

3.获得向量

式中，β^(l)(kr₀,Ω₀)为β^(l)的元素。

4.计算辅助向量：

5.计算步长：

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

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

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

现采用周期边界条件，即 其中，“:”表示从其左侧数字指示的行(列)到其右侧数字指示的行(列)，若无数字，则表示所有行(列)。此时，A_S是由中元素组成的BCCB矩阵，“Deblurringimages:matrices,spectra,and filtering”，P.C.Hansen,J.G.Nagy,D.P.O’Leary，Pages 40-43,Society for Industrial and AppliedMathematics,2006(“清晰化成像：矩阵，波谱和滤波”，P.C.Hansen,J.G.Nagy,D.P.O’Leary，第40～43页，SIAM，2006年)介绍了BCCB矩阵的构成。故A_S的谱分解为：

A_S＝F^HΛF (25)

其中，F为二维酉离散傅里叶变换矩阵，Λ为特征值矩阵，F^H为F的转置共轭矩阵，F、F^H与任意向量的乘积可通过二维傅里叶变换获得，无需明确计算F，记X为任意G_θ×G_φ维矩阵，x为X的列堆叠而成的G×1维向量，下列等价关系成立：

上式中表示等价关系。由式(25)可得：

其中，为A_S特征值形成的列向量，第2个“＝”是因为联立式(26)和式(27)可得：

式(28)中，是A_S第1列元素形成的矩阵，可由中元素向上、向左分别移动g_θc-1、g_φc-1位获得。基于上述分析，

式(29)中，第2个“＝”要求A_S为实矩阵。

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

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

综合上述原理推导过程，得本发明的步骤包括：

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

式中，三维空间内任意位置用(r,Ω)表示，r表示所描述位置与原点间的距离，Ω＝(θ,φ)表示所描述位置的方向，θ为所描述位置的仰角、φ为所描述位置的方位角；式中，(r_f,Ω_f)为聚焦点的位置坐标；

k＝2πf/c为声波的波数，f为声波频率，c为声速；

α＝[α₁,α₂,Λ,α_q,Λ,α_Q]^T为传声器计权系数组成的列向量；

上标“T”和“*”分别表示转置运算和共轭运算；

n、m、n'、m'均为球谐函数阶次，N为球谐函数阶次的截断长度；

均为Ω_f方向的球谐函数；

R_n(kr_f,ka)、R_n'(kr_f,ka)均为聚焦径向函数；

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

C为互谱矩阵；

“ο”表示Hadamard积运算；

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

记g_θ＝1,2,Λ,G_θ、g_φ＝1,2,Λ,G_φ分别为聚焦点在仰角、方位角方向的索引，也即行、列索引，G_θ为聚焦点在仰角方向的索引最大值、G_φ为聚焦点方位角方向的索引最大值G＝G_θG_φ为聚焦点总数。则利用步骤1计算得到的各点输出量构造出G_θ×G_φ维SHB输出矩阵

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

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

式中，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),Λ,t_q(kr₀,Ω₀),Λ,t_Q(kr₀,Ω₀)]^T为声场传递函数列向量；上标“H”表示转置共轭运算；

令t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，根据散射理论，

式中，为Ω₀方向的球谐函数，为Ω_q方向的球谐函数，R_n(kr₀,ka)为声信号传播径向函数；实际应用中，用截断长度N替代式中的∞来计算PSF；

记为中心聚焦点处声源的PSF矩阵，即：

式中，θ_min为聚焦声源面的最小仰角，Δθ、Δφ分别为仰角和方位角间隔步骤，Ω_c＝(θ_c,φ_c)为中心聚焦点的方向；

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

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

式中，均为G_θ×G_φ维矩阵，由中元素向上、向左分别移动g_θc-1、g_φc-1位获得，g_θc、g_φc为中心聚焦点的行、列索引，“F”、“F^-1”分别表示二维傅里叶正、逆变换；

步骤2)、计算矩阵：

步骤3)、获得矩阵

式中，β^(l)(kr₀,Ω₀)、分别为的元素，S^(l)(kr₀,Ω₀)为的元素；

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

步骤5)、计算步长：

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

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

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

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

实施例1

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

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

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

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

4、设定聚焦声源面，分别使用NNLS和本发明重构声源强度分布并成像。

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

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

上述仿真模拟在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，信号加汉宁窗，采用64段平均、66.7％的重叠率，每段时长0.25s、对应的频率分辨率为4Hz。进一步，采用MATLAB编制的NNLS和本发明程序计算各聚焦点的输出量并成像。这里，声源位置、聚焦声源面及反卷积迭代次数的设置均与仿真模拟时一致。

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

Claims (1)

1.实心球阵列三维声源识别的快速反卷积方法，其特征是，包括以下步骤：

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

三维空间内任意位置用(r,Ω)表示，r表示所描述位置与原点间的距离，Ω＝(θ,φ)表示所描述位置的方向，θ为所描述位置的仰角、φ为所描述位置的方位角；式中，(r_f,Ω_f)为聚焦点的位置坐标；

k＝2πf/c为声波的波数，f为声波频率，c为声速；

α＝[α₁,α₂,Λ,α_q,Λ,α_Q]^T为传声器计权系数组成的列向量；

上标“T”和“*”分别表示转置运算和共轭运算；

n、m、n'、m'均为球谐函数阶次，N为球谐函数阶次的截断长度；

均为Ω_f方向的球谐函数；

R_n(kr_f,ka)、R_n'(kr_f,ka)均为聚焦径向函数；

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

C为互谱矩阵；

“o”表示Hadamard积运算；

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

记g_θ＝1,2,Λ,G_θ、g_φ＝1,2,Λ,G_φ分别为聚焦点在仰角、方位角方向的索引，也即行、列索引，G_θ为聚焦点在仰角方向的索引最大值、G_φ为聚焦点在方位角方向的索引最大值，G＝G_θG_φ为聚焦点总数，则利用步骤1计算得到的各点输出量构造出G_θ×G_φ维SHB输出矩阵

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

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

式中，t(kr₀,Ω₀)＝[t₁(kr₀,Ω₀),t₂(kr₀,Ω₀),Λ,t_q(kr₀,Ω₀),Λ,t_Q(kr₀,Ω₀)]^T为声场传递函数列向量；上标“H”表示转置共轭运算；

令t_q(kr₀,Ω₀)表示声源到q号传声器的声场传递函数，根据散射理论，

其中，为Ω₀方向的球谐函数，为Ω_q方向的球谐函数，R_n(kr₀,ka)为声信号传播径向函数；实际应用中，用截断长度N替代式中的∞来计算PSF；

记为中心聚焦点处声源的PSF矩阵，即：

其中，θ_min为聚焦声源面的最小仰角，Δθ、Δφ分别为仰角和方位角间隔步骤，Ω_c＝(θ_c,φ_c)为中心聚焦点的方向；

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

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

式中，均为G_θ×G_φ维矩阵，由中元素向上、向左分别移动g_θc-1、g_φc-1位获得，g_θc、g_φc为中心聚焦点的行、列索引，“F”、“F^-1”分别表示二维傅里叶正、逆变换；

步骤2)、计算矩阵：

步骤3)、获得矩阵

式中，β^(l)(kr₀,Ω₀)、分别为的元素，S^(l)(kr₀,Ω₀)为的元素；

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

步骤5)、计算步长：

其中，g^(l)、分别为中元素形成的向量，“·”表示内积运算；

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

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