一种基于相位平均的旋转声源高分辨率定位方法

初 宁， 黄 乾， 余 亮， 宁 岳， 徐建锋， 吴大转

(1. 浙江大学 能源工程学院， 杭州 310027； 2. 浙江上风高科专风实业有限公司， 浙江 绍兴 312300;3. 上海交通大学 机械与动力工程学院， 上海 200240)

风机等通风设备在工业中应用广泛，发挥着送风、换气和排烟等重要作用，但是设备在运行过程中和潜在故障发生时都不可避免地产生噪声。利用麦克风阵列采集风机声源信号，对旋转叶片进行声源定位、绘制声源能量的空间分布云图，能够有效定位异常声源的位置，对风机的噪声品质评价、故障预警和结构优化提供新的技术支撑。

旋转声源定位的研究目前主要集中于提高算法分辨率和定位精度。Lowis等[1]将管道内旋转声源的声压表示成旋转框架下的频域格林函数与声源源强之间的乘积，再利用波束形成识别声源，但分辨率有待提高。Oerlemans等[2-3]提出了通过时域插值法去除多普勒效应的影响，从而实现旋转声源定位，但是该算法在插值运算时破坏了原始采样序列，会产生较大定位误差。王枭等[4]通过在静止声源框架中加入转速，推导得到DAMAS2 (extensions of DAMAS) 的修正算法，有效提升了单极子声源的定位分辨率，但是定位较复杂声源时存在一定的误差。Zhang等[5]提出了一种基于时域等效源方法的时域反演技术，不仅可以准确定位旋转单极子，而且可以定量地预测声源的辐射声场特征，但是该方法仅适用于近场测试。Yu等[6]提出了一种降阶循环维纳滤波器，可以从转子叶片的风洞测试信号中提取出转子叶片信号以及其他干扰，从而纯化波束形成的定位结果。Ralph等[7-8]分别基于虚拟阵列旋转法实现了旋转风机和四叶片风扇的声源定位。

为了提高旋转声源定位算法的分辨率，不少学者做过大量研究工作。由于常规波束形成给出了一个模糊的低分辨率定位成像结果，附带有较多的旁瓣影响，所以可以考虑对常规波束形成的结果进行解卷积得到高分辨率结果。经典的解卷积方法有DAMAS(deconvolution approach for the mapping of acoustic sources)[9]，CLEAN[10]，DAMAS2[11]以及SC-DAMAS(sparsity constrained-DAMAS)[12]等。将解卷积方法应用到旋转声源定位中，即可实现旋转声源高分辨率定位。莫品西[13]基于去多普勒效应给出了旋转声源的高分辨率成像算法，通过直升机旋翼试验验证了方法的可靠性。许丹[14]提出了一种旋转框架技术来消除多普勒效应的影响，再采用改进的DAMAS解卷积算法识别旋转声源，通过仿真和试验验证了该方法的有效性。Ma等[15]提出了一种通过小波压缩计算网格来提高DAMAS效率的新方法，同时还能保证DAMAS算法的高分辨率。

但是，采用时域插值法对旋转声源进行去多普勒效应时，会破坏原始采样序列的结构和相位的连续性，从而在插值时产生无法预知的误差，而且时域方法的分辨率普遍有待提高。而基于频域去多普勒效应的格林函数只在声源运动过程中的某一个瞬时成立，因此只能适应低速移动目标、远距离探测下旋转声源定位。

为此，本文提出一种基于相位平均法[16-17]的旋转声源定位方法。相位平均主要被应用于消除干扰噪声的影响[18-20]，针对干扰噪声的统计平稳特性，采用平均互功率谱消除功率谱中由干扰噪声引起的成分。针对旋转声源定位时，本文借用相位平均法的思想将旋转声信号划分成足够多个时间序列切片，在每一切片内，微小位移的旋转声源可近似等效为静止声源进行定位，并给出了相位平均技术的应用条件。为进一步提高定位分辨率，运用卷积模型改造基于波束形成的能量传播模型，利用解卷积算法消除定位模糊和提高定位精度。其优势在于，卷积模型中的点扩散函数(point spread function ，PSF)能够建立起阵列测量与声源旋转的定量联系，由此通过PSF来表征麦克风阵列定位性能、分析频率对定位分辨率的影响。最后，通过仿真信号、旋转单极子试验、风扇试验验证了方法的有效性和鲁棒性。

1 理论分析

1.1 旋转声源信号传播模型

如图1所示，在介质运动的声场中，声场声压满足下面的非齐次声学波动方程

图1 旋转单极子在有流声场下的传播模型

(1)

根据声场叠加原理，在某一位置x所测量的声压可以表示成所有等效源对于该处声压的贡献值叠加，那么式(1)的解为

(2)

式中：dVy为位置y处声源的体积微元；G(x,t|y,t′)为格林函数,表示位置y处声源在t′时刻发出的脉冲信号在接收端x处和t时刻的声场响应；δ(·)为狄拉克函数;声源按照运动轨迹xs(t′)旋转;s(t′)为声源在t′时刻的声强。

将式(2)中的3个等式联立，再利用狄拉克函数的筛分性质可以解得单极子旋转声源在自由声场下的声压解为

(3)

在有流声场中，运动声源的频域传播模型表示为

(4)

式(4)给出了运动声源在频域上的连续传播模型，基于此可以得出麦克风阵列采集声信号的离散传播模型。假定麦克风阵列由M个麦克风组成，声源平面被划分成N个网格点，其中有K个不相关单极子声源(K

P=GS+E

(5)

式中：P∈M×Z为M个通道麦克风信号做Z点傅里叶变换后的频域矩阵；S=[S1,S2,…,SK]T∈K×Z为K个声源信号的频域形式,(·)T为转置；E∈M×Z为干扰噪声；G∈M×K为频域格林函数，由式(4)可以得出

(6)

1.2 相位平均旋转声源定位算法

式(5)给出了麦克风阵列采集声信号在频域上的正向传播模型，但是这是一个非线性系统，因为声源信号和声源位置均是未知量，因此求解该反问题时，在声源所在平面划出一块扫描区域并划分成多个网格点，每一个网格点可以看成一个等效声源，这里假定有N个网格点，由此将非线性系统问题的求解转化成线性系统求解。常规波束形成(conventional beamforming, CBF)对于该线性反问题的求解为

(7)

(8)

式中：f为需要分析的目标频率；rm,n为第m个麦克风到第n个网格点的距离，n=1,2,…,N;m=1,2,…,M。RCSM(f)为麦克风阵列频域信号的互谱矩阵，具体表示为

RCSM(f)=[P(f)P(f)H]

(9)

式中：P(f)∈M×Z为阵列频域信号；[·]为数学期望，由数理统计的原理可知，样本较大时，均值近似等于数学期望，式(9)可近似表达为

(10)

式中，L为快拍数，表示将原始信号分成L段，相邻两个子段信号重叠率为50%。将式(10)的互谱矩阵估计值代入式(7)即可得到常规波束形成的输出。常规波束形成的算法原理示意图，如图2上半部分[21]所示。对于静止声源的求解，CBF是一种鲁棒性很高的方法，但是当声源旋转时，由图2可知CBF是对整个采样序列求互谱矩阵，所以无法准确给出第n个网格点的功率输出。

而相位平均法则是将声源旋转一周的采样序列划分成足够多个切片，每一个切片声源位移微小，可以近似等效为静止声源求解。为了提高分辨率，可以利用多个周期内声源经过同一弧度的数据进行多次平均。相位平均的原理示意图，如图2下半部分所示。第z个切片下声信号经过相位平均处理的互谱矩阵表示为

r{z}CSM(ma,mb)=1QLp∑Q-1q=0[({pma(I{z}+qD)})× ({pmb(I{z}+qD)})H]I{z}=Im0,n0+(z-1)Lp+II= -Lp2,-Lp2+1,…,Lp-12 Im0,n0=round(fsΔtm0,n0)

(11)

-Lp2,

Lp-12;·

图2 常规波束形成(上)与相位平均波束形成(下)的对比

(12)

式中:fr为声源转速，r/min;R为旋转半径，m，fs为采样频率；grid_div为扫描网格点的边长，m；δCBF为波束形成分辨率，m；c为声速常量，m/s；r为阵列到声源的平均距离,m;f为分析频率，Hz；Da为阵列孔径，m。式(12)给出了声源转速fr，r/min，采样频率fs,Hz,分析频率f,Hz，旋转半径R,Hz，单个切片采样长度Lp，samples，这几个主要变量之间的关系，可以看出这些变量之间是相互制约的。例如，取c=340 m/s,r=2 m,Da=1 m,f=2 500 Hz,fs=50 kHz,fr=1 200 r/min,R=0.5 m，求得Lp≤217。

式(12)第一个不等式限制了每个切片下声源旋转过的弧长不超过一个扫描网格点，这是因为对每一个切片内的旋转声源定位时将其近似看作是静止声源处理，将划分的每一个网格点看作是一个等效点声源，因此为了满足点声源的物理分辨率，需要保证声源在单个切片内运动过的轨迹不超过一个网格点的大小；第二个不等式限制了扫描网格点边长不应超过常规波束形成(conventional beamformer，CBF)算法的理论分辨率(即划分的网格物理分辨率尺寸应当小于算法的理论分辨率)，这是因为算法给出了理论分辨率的极限值，即理论上可以分辨出的两个声源的最小距离，而我们假定划分的每一个扫描网格点是一个等效点声源，在不牺牲算法的理论分辨率要求下，应当限制单个扫描网格点的大小不超过该理论分辨率。

上述分析给出了其他变量给定条件下Lp取值的上限，需要注意的是，采样长度Lp过小会导致傅里叶变换的频率分辨率过低，能量不够集中；除此之外，采样长度过小还会导致样本量太小，误差较大。因此，应用相位平均法时还需要对采样长度Lp的下界进行约束，根据大量的仿真和试验研究，我们发现当Lp<100时，能量几乎无法聚焦，难以实现成像，所以本文所述的基于相位平均的旋转声源定位方法选取的单个切片采样长度均为100或者200。当Lp=100时，可以通过式(12)确定另外几个变量的取值范围，即相位平均法的适用范围。例如，取c=340 m/s,r=2 m,Da=1 m,f=2 500 Hz,fs=50 kHz,Lp=100 samples,R=0.5 m，求得fr<2 597 r/min。当然，根据式(12)确定的只是转速的理论上限，在实际情况下，转速如果略高于理论上限，采用本方法也可以定位，只是会牺牲一定的空间分辨率。同理，当声源转速fr给定后，也可以通过式(12)确定出分析频率、声源旋转半径等其他变量的适用条件，在此不一一赘述。

将式(11)基于相位平均的互谱矩阵代入式(7)即可获得切片z中，第n个网格点基于相位平均处理的波束形成输出功率，具体表示

(13)

为了提高空间分辨率，Capon在波束形成基础上，提出了最小方差无失真响应方位估计波束器(minimum variance distortionless response，MVDR)，极大地抑制了波束形成的旁瓣[21]，基于相位平均处理的MVDR算法的方向向量为

(14)

将式(14)代入式(13)可得旋转框架下的MVDR空间功率谱表达式

(15)

式中，(·)-1为矩阵求逆。为了区别常规MVDR算法，本文中将旋转框架下基于相位平均处理的MVDR算法称作RMVDR(rotating MVDR)。

1.3 点扩散函数和RDAMAS算法

假定各个单极子声源信号互不相关，将式(5)代入式(9)可得

RCSM=G[SSH]GH+

(16)

式中，xn为[SSH]的对角线一个元素，其代表声源能量，将式(16)代入式(7)可得波束形成输出为

(17)

式中，εn1为由噪声引起的误差项。由1.2节中的分析可知，相位平均可以将每一个切片内产生微小位移的旋转声源近似等效为静止声源求解。基于以上近似，各个切片内信号进行处理时不考虑多普勒效应的影响，运动马赫数Mam,k=0，相当于式(6)中的gm,k可以由式(8)来表示，即gn等价于en，则式(17)可以表示为

(18)

式中：n1,n2=1,…,N；yn1为第n1个网格点的波束形成输出；en2为第n2个网格点的波束形成方向向量。将式(18)表示为矩阵形式

(19)

式中：ε=Diag([EEH])∈N×1为模型误差；Diag(·)为将矩阵对角线元素以向量形式表示；H∈N×N称之为传播矩阵；x∈N×1为扫描区域的所有等效声源。

利用式(19)求解x可以看成一个解卷积问题，在这里采用DAMAS算法求解该问题。用DAMAS求解时，假定ε=0，y可以看作是基于相位平均处理后的波束形成输出或者MVDR输出，即式(13)或式(15)的结果，采用高斯-赛德尔迭代法求解，迭代步骤

(20)

式(19)中，由于传播矩阵H是近似的对称块托普利兹(symmetric block toeplitz，SBT)矩阵，可以通过这个矩阵获得一个不变卷积核[22]，将Hx以二维卷积形式表示为

(21)

式中：x0为向量x的矩阵形式，此处即为扫描区域的网格点矩阵;h*∈A×B为不变卷积核，也叫做点扩散函数(point spread function, PSF)；*为同维卷积操作，即卷积结果依然是一个A×B矩阵；[·]′为将矩阵按照列转化为一个列向量。

按照式(19)可以将传播矩阵H中的元素hn1,n2表示

(22)

(23)

通常，式(21)中的PSF尺寸选为A=B=Nm=min {row(x0),column(x0)}，即取扫描网格矩阵的行数和列数中较小的值，PSF表示

h*a,b=h～n1,n2n1=N+12n2=n1+Nm+12-a Nm+Nm+12-b

(24)

2 仿真模拟

2.1 仿真条件

如图1(a)所示，声源扫描平面设置为2 m×2 m的正方形平面，扫描网格大小0.05 m×0.05 m，设置两个仿真声源，声源的初始坐标为(-0.7,0)，(0.7,0)，旋转半径为0.7 m，圆心为坐标原点；右侧的色标显示了声源的声压级大小，单位为dB，动态显示范围设置为16 dB。仿真阵列采用56通道螺旋阵列，示意图如图1(b)所示。阵列参数如表1所示。阵列平面与声源扫描平面相距1 m。声源转速设置为10 Hz，逆时针方向旋转，为了模拟实际旋转工况，仿真在10 Hz的基本转速上又添加了±1 Hz(即±10%)的频率波动作为随机干扰。采样频率设置为50 kHz，单个切片的采样点数为100个采样点，RDAMAS迭代次数300次，均选取RBF的结果作为RDAMAS的输入。

表1 麦克风阵列参数

2.2 仿真结果

本仿真模拟了1 000 Hz中心频率下的声信号定位仿真，计算了当信噪比为5 dB时的定位结果，声源信号分别采用单频信号和高斯白噪声信号进行模拟，噪声为高斯白噪声，将噪声添加到麦克风阵列接收的声信号上；定位结果云图中环形虚线表示声源旋转轨迹，虚线上的*表示的是该切片的中心时刻声源真实所在位置，如图3所示。

(a) 仿真声源初始位置

图4给出了当源信号为1 000 Hz单频信号时的定位结果，分别展示了未经相位平均的常规波束形成算法以及通过相位平均处理的3种算法的定位结果，同时还给出了点扩散函数的图像。可以看出，CBF将声源的整个旋转轨迹定位出来，无法准确定位某一处位置的声源。经过相位平均处理之后，RBF,RMVDR和RDAMAS均可以实现准确定位，其中，RDAMAS算法在RBF的基础进行了解卷积操作，有效提升了分辨率。但是RMVDR算法分辨率不高，这是由于源信号为单频信号，导致麦克风阵列各通道信号相关性较强，使得互谱矩阵的秩较小，在求逆时会出现病态解。

图5展示了信号源改为高斯白噪声信号，其余各项参数和图4仿真保持一致的定位结果。对比图4可以看出，RBF和RDAMAS并没有发生明显变化，但是RMVDR的性能得到了明显提升，说明RMVDR比较适合于高斯白噪声信号定位。因为信号源为宽带高斯信号时，信号之间相关性小，互谱矩阵几乎是满秩矩阵，求逆不会产生病态解，而对于单频等窄带信号由于式(15)中互谱矩阵求逆产生的误差会影响定位精度。

(a) CBF

(a) CBF

除此之外，从仿真的结果可以看出转速波动对本文所述的相位平均法没有造成影响。这是因为相位平均法的基本原理是将原始信号划分成多个切片分块处理，每一个切片内的声源位移微小，可以近似等效为静止声源求解，所以不必关心每个时刻声源的具体转速与运动轨迹，只需采集到每个切片内的声信号就可进行运算，这也是相位平均法与时域去多普勒效应法的本质区别。时域去多普勒效应是在已知声源运动轨迹和转速的条件下，通过时域插值法由麦克风信号重建声源位置信号，所以当转速存在随机波动时，无法准确判断扫描区域的网格点(等效源)在各个采样时刻的位置，插值时会产生较大误差。这也说明相位平均法比时域去多普勒效应法的鲁棒性更高，能够适应转速波动在±10%以内的旋转声源定位。

3 试验验证

3.1 试验参数

为了验证本文提出的方法的有效性，我们进行了旋转蓝牙音箱试验和工业风扇试验，两组试验均在消声室中进行。其中，蓝牙音箱和工业风扇的各项试验参数如表2所示。麦克风阵列采用和仿真相同的56通道螺旋阵列，麦克风采用1/4英寸驻极体麦克风，频率响应20 Hz～20 kHz。旋转蓝牙音箱正对麦克风阵列平面发声，测试距离1.3 m，位置关系如图6所示。工业风扇由于旋转时产生的气流太大，防止麦克风阵列采集信号被过多的流滞噪声干扰，将风扇背对阵列面，距离1.52 m，位置关系如图7所示。

表2 试验参数

图6 旋转蓝牙音箱位置图

图7 工业风扇位置图

3.2 试验结果

3.2.1 蓝牙音箱旋转试验结果

如图8所示为旋转蓝牙音箱在600 r/min转速，2 500 Hz单频信号下的定位结果。采样频率为50 kHz，因此声源旋转一周总采样点数为5 000，将每一圈的麦克风阵列采集到的原始信号划分成25个切片，单个切片的采样点数为200，单个切片声源旋转过的角度为14.4°。图8分别展示了CBF的定位结果，以及基于相位平均处理的RBF，RMVDR和RDAMAS定位结果，采用RBF的结果作为式(20)中RDAMAS迭代算法的输入。CBF，RBF和RDAMAS动态范围16 dB, RMVDR动态范围8 dB。虚线大圈表示声源旋转圆周，每一个切片下声源实际位置用虚线小圈表示(注：由于测量和对准误差，标示位置可能和实际位置略有微小误差)。

(a) CBF

图8(a)的CBF因为对整段采样信号进行波束形成，所以无法准确定位某一处位置声源，只能得到声源运动轨迹。RBF，RMVDR和RDAMAS均可以准确定位对应切片下的声源位置，由于4个蓝牙音箱一致性不太好，所以展示出来的声压幅值大小不一致，在定位云图上表现为对应声源颜色深浅不一样，图8(c)和图8(d)尤为明显。另外对比图8(c)和图8(d)声源旋转过的角度基本等于4个切片长度下声源真实旋转过的角度(14.4°×4=57.6°)。对比图8(c)和图8(g)，RDAMAS在RBF基础上采用迭代法极好的抑制了旁瓣，提高了定位分辨率。图8(e)和图8(f)展示的RMVDR结果虽然也能够定位出声源位置，但成像较模糊，分析原因可能是单频信号的相关性较大，导致互谱矩阵秩较小，是一个病态矩阵，在求逆时容易产生较大误差，这与图4(e)和图4(f)仿真结果一致。

其余试验参数不变，仅将声源信号频率提高到4 000 Hz，结果如图9所示。对比图8(b)和图9(b)可以看出，随着频率升高，PSF变得尖锐，分辨率提高，这与理论分析相吻合。同时，各种算法的空间分辨率也随着频率提高而升高，尤其是RBF的旁瓣明显削弱，性能提升显著。但是，随着频率升高，波长变短，可能会超出阵列的空间采样极限(麦克风间距不得大于半个波长)，因此定位结果中出现了较多的栅瓣影响。

(a) CBF

将图8的声源信号改为高斯白噪声信号，定位中心频率依然为2 500 Hz，其余各项参数保持不变，结果如图10所示。对比图10和图8，RBF和RDAMAS的结果变化不大，但是RMVDR性能提升明显，分辨率显著提高，说明RMVDR算法对高斯白噪声信号适应性更好，这与仿真的结果相一致。

(a) CBF

3.2.2 工业风扇试验结果

工业风扇为三叶片轴流风扇，分别采集了1 130 r/min转速和1 278 r/min转速下的信号，测试距离1.52 m。图12为分析频率2 500 Hz的结果，单个切片采样点数为100，1个大圆虚线框表示叶尖的运行轨迹，4个小圆虚线框表示对应切片下的扇叶实际位置(注：由于测量和对准误差，标示位置可能和实际位置略有微小误差)。CBF依然无法准确定位声源位置，基于相位平均处理的RBF，RMVDR和RDAMAS均可以准确定位该采样切片对应的声源位置。其中，RBF的分辨率较差，这是因为比起蓝牙音箱旋转，工业风扇的气动噪声源声场更加复杂，干扰大，而且声源不是呈现离散点声源分布，而是呈现分布式声源的分布特征(即不满足式(5)中的K≪N)。相对而言，RMVDR的结果比较理想，这是因为气动噪声源的频谱呈现宽带连续谱+低频离散线谱的特征，比较适用于RMVDR算法，具体表现形式如图11所示。低频相邻线谱频率差值约为55.3 Hz，约等于风扇的叶频(1 130/60×3=56.5 Hz)，中高频为宽带连续谱，分析频率选为中高频率。将RMVDR算法的结果作为RDAMAS迭代的输入，对比图12(e)和图12(g)，经过迭代运算之后，定位分辨率显著提升。

(a) 时域波形

RBF，RMVDR和RDAMAS算法定位的声源位置均分布在扇叶尖部，这是因为对于旋转叶片，气流流过叶尖引起的不规则流动产生的涡流噪声是气动噪声源的主要组成部分，而这一部分噪声主要为中高频噪声[24]。随机涡从尾缘脱落，引起与其相反方向的局部环量变化，从而引起桨盘上诱导速度的变换，导致叶尖产生随机升力的脉动，最终产生频带较宽的涡流脱落噪声。涡流脱落噪声位于叶尖，图12的工业风扇试验定位的最大噪声源就位于叶轮尖部，与理论分析相吻合。

(a) CBF

将风扇转速调为高档(1 278 r/min)，分析频率为3 500 Hz，其余试验条件不变，定位结果如图13所示。RBF，RMVDR和RDAMAS算法依旧保持很好的性能，RDAMAS依旧展现出极高的分辨率。随着分析频率增加，从图13(b)可以看出，PSF变得更加尖锐，分辨率提高；定位结果的分辨率也明显提高，尤其是RBF对旁瓣的抑制作用变好。

(a) CBF

4 结 论

本文基于相位平均原理提出了针对旋转声源的定位方法，并给出了相位平均技术所适用的条件。其基本思路是将旋转声信号进行切片平均，在每一个切片内近似等效于静止声源处理，并在此基础上提出了基于解卷积算法的RDAMAS算法，在波束形成基础上显著提高了定位的空间分辨率，通过仿真、旋转蓝牙音箱试验以及工业风扇试验验证了算法的有效性。由于相位平均只对原始采样信号进行切片重排序，不像插值法去多普勒效应那样会改变信号自身原始结构，所以该方法具有很高的鲁棒性，几乎所有的静止高分辨率方法均可以基于相位平均来实现对旋转声源的定位。除此之外，本文还引入了PSF来表征阵列性能以及频率对于定位分辨率的影响，通过PSF图像的尖锐程度反映定位的分辨率大小。

仿真和试验表明，RBF能保持较高的鲁棒性，但分辨率不足；而RMVDR更适用于高斯白噪声信号定位，空间分辨率高于RBF。RDAMAS在RBF或RMVDR基础上进行迭代运算，大幅提高了定位的分辨率。综合1.2节针对式(12)的分析讨论，再结合大量仿真和实验结果，可以大致给出本文所述方法的适用范围：当采样频率保持50 kHz,测试距离2 m时，本方法适用于声源转速≤3 000 r/min，分析频率≥1 500 Hz,旋转半径0.3 m≤R≤1.0 m，声源转速波动在基本转速的±10%的旋转单极子定位成像。针对三叶片工业风扇，本方法可以实现对转速1 278 r/min风扇叶片的叶尖气动噪声源定位，由此本方法可以推广应用到中低转速(≤1 500 r/min)下，叶片数量较少(不超过5片)的风机气动噪声源定位。由于气动噪声呈现出中高频(≥1 500 Hz)宽带谱的特性，所以能够采用本文所述的算法进行定位，而且定位的能量最大的声源位置位于叶轮尖部。