随着科技的发展，麦克风阵列已广泛应用于军事侦查、工业探测、人机交互等领域[1]．大多数基于麦克风阵列的声源定位、语音增强、语音识别等信号处理算法都依赖于阵列的几何结构信息，阵元位置的误差将直接影响算法的性能．传统的麦克风阵列通常有着较大的孔径和固定的几何结构．近年来随着消费电子的发展，无线耳机、助听器等带有双微型麦克风阵列的嵌入式设备越来越普及，其子阵列之间的距离在使用过程中随时可能发生变化，不确定的阵列间距影响了后续语音信号处理的效果，因此对双微型麦克风阵列进行实时、精确的位置校准就显得尤为重要．麦克风阵列的位置校准可以分为主动校准和被动校准两大类[2]．主动校准方法利用位置已知的人工校准信号进行校准，可以达到较高的精度，但实现较为不便；被动校准方法则直接利用位置未知的声音信号进行校准，实现较为容易，但精度相对较低．现有的阵列位置校准算法利用接收信号的到达时间差(TDOA)、到达时间(TOA)、功率谱密度(PSD)、到达方向(DOA)等信息，通过不同的信号处理算法来实现麦克风阵元的定位[3]．文献[4]利用DOA信息结合最大似然算法，对圆形麦克风阵列进行位置校准，但该方法须要已知阵列结构且要求每个节点内包含多个麦克风．文献[5]利用DOA和TDOA信息共同构建目标函数，实现混响条件下的位置校准，但其不适用于麦克风数量较少的情况．文献[6]利用信号与混响的PSD信息构造相干性矩阵，结合MDS算法校准各个麦克风的位置．文献[7]利用声源传播时间(TOF)信息获得距离矩阵，通过最小Frobenius范数求解声源和麦克风位置，但该方法要求声源可控且与麦克风之间完美同步．文献[8]利用TDOA信息，结合几何解算和LS估计求解出阵列间距，但该方法假设模型处于理想环境，忽略了TDOA误差对校准结果的影响．基于TDOA的被动校准方法实现容易，计算量小，适合于硬件资源有限、实时性要求高的无线耳机、助听器等带有双微型麦克风阵列嵌入式设备的应用场景．针对目前基于TDOA的双微型麦克风阵列位置校准的精度受到采样率和阵元间距影响的问题，提出了一种改进的阵列位置校准算法，通过改进算法提高时延估计和几何解算的精度，从而提高阵列位置校准的精度．1 双微型麦克风阵列建模双微型麦克风阵列由两个微型麦克风子阵列组成，子阵列内的阵元个数为N，阵元间距为1~3 cm，两个子阵列之间的间距为10~20 cm．本研究的阵列中，子阵列阵元数N为2，阵元间距L设置为2 cm，定义两个子阵列的麦克风轴中点之间的距离作为子阵列之间的距离d，声源信号的传播方向与麦克风轴中点连线的夹角θ为声源到达方向．本研究的位置校准问题即是利用麦克风接收到的声源信号对子阵列之间的距离d进行估计．在无线耳机、助听器等带有微型麦克风阵列的嵌入式设备的使用场景中，声源与麦克风之间的距离通常远大于麦克风阵元的间距，可以将声波视为平面波，采用远场模型建模，降低系统的计算复杂度．双微型麦克风阵列的远场声源模型如图1所示．10.13245/j.hust.239526.F001图1双微型麦克风阵列远场声源模型2 TDOA算法分析2.1　广义互相关算法原理在麦克风阵列中，任选两个麦克风，将其在k时刻接收到的信号设为y1(k)和y2(k)，则这两路信号可以表示为y1(k)=a1s(k-t)+v1(k)；(1)y2(k)=a2s(k-t-τ21)+v2(k)，(2)式中：a1和a2为由于传播效应导致的衰减因子；s(k)为位置未知的声源信号；t为从声源到麦克风1的传播时间；v1k和v2k为麦克风处的噪声信号；τ21为麦克风2和1之间的到达时间差(TDOA)．两个观测信号y1k和y2k之间的互相关函数R12(p)可以表示为Ry1y2(p)=E[y1ky2(k+p)]．(3)将式(1)和(2)代入式(3)可以推导出Ry1y2(p)=a1a2Rss(p-τ21)+a1Rsv2(p+t)+ a2Rsv1(p-t-τ21)+Rv1v2(p)，(4)式中：Rmm(∙)为信号m的自相关函数；Rmn(∙)为信号m和n的互相关函数．假设声源信号和噪声都不相关，式(4)可简化为Ry1y2(p)=a1a2Rss(p-τ21)．(5)由式(5)可知，当p=τ21时，Ry1y2(p)达到最大值．因此，求出使Ry1y2(p)取得峰值时对应的p就得到了y1(k)和y2(k)之间的相对时延，即τ^=argmaxp Ry1y2(p)．(6)在实际中互相关函数是未知的，须通过近似得出，将式(3)定义的互相关函数Ry1y2(p)用其时间平均估计进行替代，即R^y1y2(p)=1K∑i=0K-p-1y1k+iy2(k+i+p)，(7)式中K为样本序列的长度．将式(1)、(2)代入式(7)可得实际情况下的互相关函数为R^y1y2(p)=a1a2R^ss(p-τ21)+a1R^sv2(p+t)+a2R^sv1(p-t-τ21)+R^v1v2(p)，(8)式中R^sv1(∙)，R^sv2(∙)和R^v1v2(∙)不再满足为0的条件．式(8)等号右边的后三项均作为噪声对R^y1y2(p)产生的影响．为减少噪声对互相关算法的影响，研究者们提出各种不同的频域加权方法，称为广义互相关法[9]．将广义互相关函数定义为Ry1y2GCC(p)=F-1[ψy1y2(f)]=∫-∞+∞ψy1y2(f)ej2πfpdf=∫-∞+∞ϑ(f)φy1y2(f)ej2πfpdf，(9)式中：F-1[∙]表示离散时间傅里叶逆变换；φy1y2(f)=E[Y1(f)Y2*(f)]为互功率谱，且Yn(f)=∑kyn(k)e-j2πfk,n=1,2；ϑ(f)为频域加权函数；ψy1y2(f)=ϑ(f)φy1y2(f)为广义互频谱．频域加权函数ϑ(f)有许多不同的选择，产生各种不同的广义互相关(generalized cross-correlation，GCC)方法，常见的有Roth，SCOT及PHAT加权等，这里采用PHAT加权函数，具体为ϑ(f)=φy1y2(f)-1．(10)2.2　采样率和阵元间距对TDOA的影响由广义互相关算法原理可知，TDOA误差的来源主要有两处：一是采用了相关函数的估计值代替理论值而引入的量化误差；二是噪声的存在影响互相关函数而导致的峰值估计误差．广义互相关算法中定义的互相关函数是严格数学意义上的统计平均或在平稳遍历条件下替代统计平均的无穷时间平均[10]，然而在实际中，严格数学意义上的统计平均或无穷平均是不可能的，只能用有限的时间平均来代替，在实际应用中系统的离散化必然会引入量化误差．将系统的采样率设为fs，则其时间分辨率为Δt=fs-1，对应于距离上的精度为Δd=cΔt=cfs-1．通常而言，为了保证两个阵元之间接收信号的空间信息，阵元间距应该远大于距离精度．然而在微型阵列中，麦克风阵元间距往往较小，不一定能满足上述条件．例如当采样率为48 kHz时，对应的距离精度为0.708 cm，而本研究的双微阵列中的子阵列内阵元间距为2 cm，子阵列间的距离为10~20 cm，当采样率较低时系统的量化误差将对计算结果产生较大影响．对于噪声的影响，由于双微型麦克风阵列中子阵列内的两个麦克风间距很小，接收到的噪声相关性较大[11]，即互相关函数式(8)中的噪声项R^v1v2(p)较大，从而将导致其TDOA估计的误差较大；两个子阵列之间的麦克风间距较大，受噪声相关性的影响较小，因而TDOA估计精度相对较高，同时较大的间距下受到量化误差的影响较小．3 基于TDOA的阵列位置校准算法3.1　改进的广义互相关算法算法流程如图2所示．为了减少采样率对于TDOA的影响，通常采用插值方法来提高原始数据的采样率，对插值后的数据进行广义互相关计算，从而提高时延估计的精度．然而该方法将导致计算量成倍增加，不适用于硬件资源有限的嵌入式设备．提出一种改进的广义互相关算法，在对广义互10.13245/j.hust.239526.F002图2改进的GCC算法流程相关函数进行峰值检测后，取峰值点附近的部分点构造三次样条函数，在小范围内进行三次样条插值，再进行二次峰值检测后得出时延估计结果．在插值区间内构造三次样条方程S(x)=ai+bix+cix2+dix3，根据样条曲线插值的特点可以推导出：ai=yi;bi=yi+1-yihi-hi2mi-hi6(mi+1-mi);ci=mi/2;di=(mi+1-mi)/(6hi),(11)式中：yi为输入待插值点的值；hi=xi+1-xi为子区间的步长；mi为函数的二次微分值Si''(xi)，须通过边界条件进行求解．采用自然边界(natural spline)条件将边界点的二阶导数指定为0，即S''(x0)=S''(xn)=0．可以构造以m为未知数的线性方程组，即ABC∙m0,m1,⋯,mnT=60,y2-y1h1-y1-y0h0,y3-y2h2-y2-y1h1,⋯,yn-yn-1hn-1-yn-1-yn-2hn-2,0T,式中：A=10h02(h0+h1)0h100⋮⋮0000；B=00⋯00h10⋯002(h0+h1)h2⋯00⋯00⋯hn-32(hn-3+hn-2)00⋯0hn-200⋯00；C=0000⋮⋮00hn-202(hn-2+hn-1)hn-101．求解方程组得出mi值后代入式(11)即可得到子区间内的三次样条方程S(x)，再根据所需的上采样倍数进行数据插值．三次样条插值的光滑性较好，同时可以估计出样本范围之外的局部最大值和最小值[12]．采用局部的三次样条插值可以在互相关函数中产生新的峰值点，提高峰值检测的精度，同时由于仅对原有峰值点两侧的少量点数进行插值，计算量较小，适合于嵌入式设备．3.2　改进的几何解算算法由2.2节的分析可知：在双微型麦克风阵列中，子阵列内两个麦克风较小的阵元间距将导致较大的时延估计误差，因此虽然子阵列内的阵元间距为已知量，但直接采用其时延估计值进行几何计算将会引入较大的误差．提出一种改进的几何解算法，结合子阵列之间的时延值与子阵列内部的阵元间距值，推导出阵元间距的解算公式．在远场模型下，任意两个麦克风i,j之间接收声音信号的时间差有如下关系τijc=dijcos θij，(12)式中：τij为由时延估计得到的麦克风i和j之间接收信号的相对时延值；c为空气中声音传播的速度；dij为两个麦克风间的距离；θij为两个麦克风的连线与声源方向的夹角．分别将麦克风m1和m3、m1和m4的时延估计值代入式(12）可得τ13c=d13cos θ13;τ14c=d14cos θ14, (13)式中：θ13为未知的声源到达角度θ；d13为待校准的子阵列间距d，因此有θ=arccos (τ13c/d)．(14)根据图3中的几何关系可以推出：d14=d2+L2；(15)θ14=θ-arctan (L/d)．(16)10.13245/j.hust.239526.F003图3双微型麦克风阵列与声音信号的几何关系图将式(15)和(16)代入式(13)并整理可得τ14=(dcos θ+Lsin θ)/c．(17)将式(14)代入式(17)可得：τ14=τ13c+Lsin [arccos (τ13c/d)]c=τ13+(L/c)1-(τ13c/d)2.整理可得阵元间距的解算公式为d=τ13c1-[(τ14-τ13)c/L]2．相较于文献[8]提出的解算方法，本文算法避免了子阵列内部误差较大的TDOA估计值，采用的是两个子阵列之间误差较小的TDOA，再结合已知的几何结构信息进行推导计算，得到更为准确的阵列校准结果．4 实验与分析在Matlab中搭建双微型麦克风阵列远场模型，设置子阵列内麦克风间距为2 cm，子阵列间距为16 cm，即4个麦克风位置分别为m1(0，0)，m2(0，0.02)，m3(0.16，0)，m4(0.16，0.02)；使用采样率为48 kHz的语音作为声源，每一帧的长度取22 ms，放置于距离原点3 m的随机角度位置，作为远场的未知声源；在麦克风处加入10 dB的高斯白噪声作为干扰噪声，广义互相关加权函数使用PHAT加权．分别采用三种算法实现TDOA估计：文献[13]采用的原始数据+GCC-PHAT方法，文献[14]采用的线性插值+GCC-PHAT方法，本文改进算法为原始数据+GCC-PHAT+局部3次样条插值方法．对数据的插值倍数均设置为10倍，每种方法分别进行200次随机声源位置下麦克风m1和m3之间的时延估计，定义时延估计的均方根误差为στ=1N∑i=1N(τe-τr)2，式中：N为样本数量；τe为估计的时延值；τr为实际的时延值．由表1可见：本文算法的时延估计误差相比于文献[13]算法降低了85%，即本文算法在少量增加计算量的情况下大幅提高了时延估计的精度．相较于时延估计精度相近的文献[14]算法，本文算法的计算量仅为其计算量的14%左右．10.13245/j.hust.239526.T001表1时延估计结果对比算法用时/ms时延估计误差/μs文献[13]0.122 15.699文献[14]1.146 60.806本文0.156 90.861在TDOA估计的基础上进行几何解算，实现阵列位置的校准，仿真结果如表2所示．可以看出：在10~20 cm范围内，随着实际子阵列间距的增大，校准的误差将随之略微增加，其误差主要来源于远场假设引入的误差．10.13245/j.hust.239526.T002表2距离估计结果实际距离/cm距离估计平均值/cm距离估计均方根误差/cm1010.210.651212.200.691414.380.811616.480.941818.480.982020.561.26在现有的基于TDOA的阵列位置校准文献中，文献[8]与本文采用相同的远场建模，可以直接进行对比.将实际子阵列间距设置为16 cm，采用本文算法与文献[8]算法在相同的TDOA估计值下进行阵列校准，仿真结果如图4所示，图中：n为仿真次数；De为距离估计值．10.13245/j.hust.239526.F004图4阵列校准算法对比可以看出：在相同的时延估计值下，本文算法相比于文献[8]的阵列校准算法，距离估计的结果更为稳定，偏移量更小．在200次阵列校准仿真中，文献[8]算法的均方根误差为1.91 cm，而本文算法的均方根误差为0.94 cm，位置校准的误差降低了51%．5 结语针对带有双微型麦克风阵列的嵌入式设备，搭建了双微型麦克风阵列远场模型，分析了TDOA算法的误差来源，提出了一种基于TDOA的阵列位置校准算法，采用改进的广义互相关算法和改进的几何解算算法来减少采样率和阵列间距对双微型麦克风阵列的影响，提高了时延估计和阵列位置校准的精度．实验结果表明：在10 dB信噪比、48 kHz采样率的条件下，本文算法可在计算量相近的情况下将时延估计误差降低85%，位置校准误差降低51%．