Page 19 - 《应用声学》2022年第1期

第 41 卷 第 1 期 高玥等： 采用球谐分解的二范数广义逆波束形成 15

min ∥q∥ subject to
¯ p − A¯ q 谐函数的谱，表示为
¯
x∈C N 1 2 ∫ 2π ∫ π

m
∗

N
p nm = p(θ, ϕ)Y (θ, ϕ) sin θdθdϕ
∑ n

=
¯ p − ¯ a n ¯q n
6 ¯ε, (12) 0 0

Q
n=1 2 ∑
m
∗
其中， ¯ p 和 ¯ a n 是由 p 和 a n 归一化之后得到的向量， ≈ α nm p(θ q , ϕ q )Y (θ q , ϕ q ) , (14)
q
n
q=1
即 ¯ p = p/ ∥p∥ ， ¯ a n = a n / ∥a n ∥ ；且 ¯ε = ε/ ∥p∥ ， ¯ x
2 2 2 其中，p(θ q , ϕ q ) 为第 q 个阵元的声压分布，α nm 是由
的一个分量 ¯q n = (∥a n ∥ 2 /∥p∥ 2 )q n 。 q
积分转换到累加和形式的近似，对第 q 个阵元的 m
借助数学工具中凸优化算法建模软件包对
阶 n 次项的加权系数。由于球形阵列传声器个数 Q
式 (12) 进行求解。由于它的稀疏性，即使在单频情
是有限的，因此 p(θ q , ϕ q ) 可以通过截断 P 阶后的逆
况下，此算法也能提供高分辨率的声成像结果。然
球傅里叶变换(Inverse spherical Fourier transform,
而，当测量数据受到噪声污染时，稀疏解的非零分量
ISFT)近似表示为 [25]
会出现在与实际源位置不同的位置。
P n
∑ ∑ m
1.2.2 L2范数广义逆波束形成 p(θ q , ϕ q ) = p nm Y (θ q , ϕ q ). (15)
n
n=0 m=−n
L2 范 数 的 广 义 逆 波 束 形 成 算 法 (L2-
则可将式(15)转换成矩阵形式表示为
Generalized inverse beamforming, L2-GIBF) 是用
Tikhonov 正 则 化 [24] 修 正 问 题 的 方 法 来 求 解 公 p = Y P p nm , (16)
式 (9) 欠定方程的解。通过添加残差范数项和正 其中，p 为各个阵元采集声压数据经过短时傅里叶
则化参数λ 修正等式后，得到正则化的解为 变换的频谱，p nm 为球谐域谱，Y P 矩阵的阶数为
2
2
2 2 2 M × (P + 1) ，包含在每个阵元上计算得到 P 阶的
ˆ q(λ) = arg min{∥p − Aq∥ + λ ∥q∥ }
m
球谐函数Y (θ q , ϕ q )。
H
H
2
= A (AA + λ I Q ) −1 p n
对式 (16) 中的线性系统进行反演，得到了 SFT
2
2
∗
= V (S + λ I Q ) −1 SU p
的数值解为式(17)：
Q 2
∗
∑ s u p
= i i v i , (13) p nm = T P p, (17)
2
s + λ 2 s i
i=1 i
其中，T P 为Y P 的逆矩阵。
其中，A = USV 为奇异值分解结果，U ∈ C Q×Q
∗
短时傅里叶变换后得到的信号频谱先经过
是酉矩阵，S ∈ R Q×Q 是一个对角矩阵，其中包
ISFT 转换到模态域，然后通过 SFT 转换回来，则实
含的奇异值按降序 s 1 > s 2 > · · · > s Q > 0 排列，
现了球谐分解步骤 [25] 。但是SFT和ISFT只能用于
V ∈ C N×Q 满足V V = I Q ，正则化系数λ 定义了
2
∗
单球阵阵列条件下，为了克服这一限制，本文提出
滤波器的强度。
将球谐系数估计法引入此步球谐分解操作中。因此
1.3 基于球谐分解的L2范数广义逆波束形成 公式(16)、公式(17)可以用公式(1)、公式(8)的形式
基 于 球 谐 分 解 的 L2 范 数 广 义 逆 波 表示，则可在任意阵型下实现球谐分解预处理步骤。
束 形 成 (Spherical harmonics decomposition-L2- 将预处理后的数据再进行本文 1.2.2 节中叙述的二
Generalized inverse beamforming, SHD-L2-GIBF) 范数广义逆波束形成算法，实现基于球谐分解的 L2
算法的原理是先通过球谐分解 (SHD) 的预处理来 范数广义逆波束形成算法。
稳定后续的求逆数值问题，尤其在中低频范围内效
1.3.2 正则化参数选择
果更好，而后通过广义逆波束形成 (GIBF) 处理截
对于公式 (13) 中正则化参数的选取方法，本文
断后的阵元接收数据进行声成像。
采用的为 L-曲线准则 [26] ，在二维上，以 lg-lg 为尺
1.3.1 算法原理 度，(∥p − Aq∥ , ∥q∥) 的值构成了如图 1 所示形状类
p nm 为复数球谐波系数，是由球傅里叶变换 似于 L 的曲线，故被称为 L 曲线。为了平衡欠正则
(Spherical Fourier transform, SFT) 计算得到的球 化和过正则化，在L曲线曲率最大的点(即由竖直到