本发明属于声信号处理算法中的参数识别方法研究，根据声场频率将高阶谐波分量进行正则化，结合贝塞尔函数的加法定理构造转换矩阵，将局部声场系数转换成全局声场系数，从而减小麦克风误差和贝塞尔函数零值对声场系数的影响，提高了声场重构准确度。
背景技术：
在现代战场上，随着电子干扰技术的发展，电磁探测环境越来越恶劣，特别是反辐射导弹的出现和使用，对雷达等有源设备的生存提出了严重的挑战。基于物理原理的声场重构方法旨在一定的空间区域内重构一个与目标声场完全相同或者尽可能接近的物理声场，要求在采集声信号的过程中较好地保留原始声场的空间信息。针对复杂城市战场环境下的多变性和隐蔽性，波域声场重构技术能够为现代战争近距离作战提供一种新的有效探测手段。然而，在声场重构过程中获取声场的高阶圆柱谐波所需要的麦克风数量是相当多的，尤其是在重构大范围高频声场的条件下。
空间声场重构系统用于捕获和重建空间中某一区域内的连续声场。这需要用麦克风阵列在给定的空间区域内采集声场声压，然后准确地对其进行重构。声场重构的关键在于如何获得待重构声场的高阶圆柱谐波成份，然而，麦克风的噪声和贝塞尔函数零值会影响圆柱谐波系数，即声场系数的测量值。因此，减小声场系数的测量误差对于提高声场重构准确度有着十分重要的意义。
相比于全向性麦克风阵列，高阶麦克风阵列具有自适应干扰抑制以及高分辨率波达方向(directionofarrival,doa)估计等诸多优点。高阶麦克风使用刚性阵列设计，可以进一步减小贝塞尔零点问题。使用多个高阶麦克风包围待重构区域，然后利用高阶麦克风采集局部声场的声压，由此计算局部声场的谐波系数，并根据声场频率将不存在的高阶谐波分量进行正则化，可以减小全局声场系数的重构误差。结合贝塞尔函数的加法定理构造转换矩阵，将在每个高阶麦克风处获得的局部声场谐波系数转换为待重构区域的全局声场谐波系数。由于转换矩阵内的贝塞尔函数项的额外变化，使用高阶麦克风大大减小了转换矩阵的条件数，使用半径可变的圆形阵列可以进一步提高阵列的鲁棒性。
综上，本技术利用高阶麦克风采集局部声场声压，根据离散傅里叶逆变换计算局部声场的谐波系数，然后根据声场频率将不存在的高阶谐波分量进行正则化，结合贝塞尔函数的加法定理构造转换矩阵，将局部声场系数转换成全局声场系数，从而减小麦克风误差和贝塞尔函数零值对声场系数的影响，提高了声场重构准确度。
技术实现要素：
要解决的技术问题
为了提高空间声场的重构准确度，本发明在传统的高阶ambisonics(highorderambisonics,hoa)重构方法的基础上，使用高阶麦克风采集局部声场声压，根据离散傅里叶逆变换计算局部声场的谐波系数，然后根据声场频率将不存在的高阶谐波分量进行正则化，结合贝塞尔函数的加法定理构造转换矩阵，将在每个高阶麦克风处获得的局部声场系数转换为待重构区域的全局声场系数，从而有效提高了声场系数的重构准确度。
技术方案
一种基于高阶麦克风阵列的声场重构方法，采用高阶麦克风阵列采集声场信号；其特征在于步骤如下：
步骤1：计算局部声场系数
步骤1a：根据高阶麦克风的不同阵列设计，选择不同的模态强度：
式中：r0为无限长圆柱形刚性挡板的半径，r0≤rm，jv(·)为v阶贝塞尔函数，hv(·)为第一类v阶汉克尔函数，(·)′为求导运算符号，k为波数；
步骤1b：根据声压的圆柱谐波展开式，利用指数函数的正交性，得到声场系数的表达式，即：
式中：rm为高阶麦克风的半径，φqq′为麦克风极坐标下的极角，e为自然指数，p(k,rm,φqq′)为高阶麦克风采集到的局部声场的声压，nq′为白噪声；
步骤1c：对步骤1b中的积分进行空间采样，利用有限的求和来精确逼近连续傅里叶逆变换，得到离散傅里叶逆变换的形式，并且根据香农采样定理，将麦克风数量限制为最小数量，并采用等角度采样方案，即：
式中：q′为麦克风的数量，aq′为第q′个麦克风的采样权重，v为高阶麦克风的阶数；
局部声场系数改写为：
步骤2：局部声场截断
根据高阶麦克风的阶数，确定v阶麦克风的半径，即：
式中：c为声速，fmax为v阶麦克风的最大可重构声场频率；
当声场频率为f时，声场中实际存在的或者已经被激活的圆柱谐波的阶数仅为：
将局部声场系数
比如v为3时，则
步骤3：根据高阶麦克风的阶数和阵列的拓扑结构来构造转换矩阵，即：
式中：
步骤4：计算全局声场系数向量：
式中：
得到
步骤5：利用声场的谐波表达式，得到重构声场的声压分布，即：
式中：
有益效果
本发明改进了传统的hoa声场重构技术测量声场系数的过程，使用高阶麦克风采集局部声场声源并计算局部声场系数，然后结合贝塞尔函数加法定理构造转换矩阵，将在每个高阶麦克风处获得的局部声场系数转换为待重构区域的全局声场系数，从而有效提高了空间声场的重构准确度，并实现了麦克风阵列系统的精简。
本发明通过对局部声场进行截断，舍弃了不存在的高阶谐波分量，减小了白噪声对局部声场系数的测量误差，从而减小了全局声场系数的重构误差。
本发明通过引入高阶麦克风阶数的变化，使得转换矩阵中贝塞尔函数项产生额外的变化，减小了矩阵的条件数，使用半径可变的圆形阵列进一步提高了阵列的鲁棒性。
附图说明
图1为基于高阶麦克风阵列的声场重构方法结构框图
图2为开放阵列的模态强度曲线
图3为刚性阵列的模态强度曲线
图4为三阶麦克风的转换矩阵条件数变化曲线
图5为半径可变的圆形阵列拓扑结构的转换矩阵条件数变化曲线
图6为平面波的目标声场系数、重构声场系数及其重构系数误差曲线
图7为点声源的目标声场系数、重构声场系数及其重构系数误差曲线
图8为平面波的目标声场和重构声场效果图
图9为点声源的目标声场和重构声场效果图
图10为平面波和点声源的声场重构误差
图11为平面波声场重构的归一化均方误差曲线
图12为点声源声场重构的归一化均方误差曲线
图13为平面波声场重构的白噪声增益曲线
图14为点声源声场重构的白噪声增益曲线
具体实施方式
现结合实施例、附图对本发明作进一步描述：
装置实施例：在运用本发明对空间声场进行重构时，首先利用高阶麦克风采集局部声场声压，并根据傅里叶逆变换的离散形式计算局部声场系数，然后根据低阶声场的截断阶数极限将低阶声场中的高阶谐波分量进行正则化，得到正则化后的低阶声场系数。接着根据高阶麦克风的阶数和阵列的拓扑结构计算转换矩阵，将在每个高阶麦克风处获得的局部声场系数转换为待重构区域的全局声场系数。最后，根据声场的谐波表达式得到重构声场的声压分布。
本发明在matlab仿真平台上完成验证，具体的仿真参数为：目标声场分别选用平面波声场和点声源声场，声源频率为1400hz，声速为342m/s，其中点声源的极坐标为(4m，0)，平面波的方向为φo＝π/4，幅值为单位1。重构区域半径为r＝3m，高阶麦克风阵列的半径为rq＝3m。高阶麦克风的拓扑结构为圆形，麦克风等角度分布。
本发明包括以下步骤：
步骤一：计算局部声场系数。该模块的主要功能是利用高阶麦克风采集局部声场的声压，根据麦克风的位置计算局部声场系数
所述步骤一还包括如下子步骤：
步骤a：傅里叶逆变换。根据声压的圆柱谐波展开式，局部声场的表达式如下：
式中：x为测量点的极坐标，x＝(r,φ)，k为波数，k＝2πf/c，v为高阶麦克风的阶数，
利用指数函数的正交性，在上式两边同乘以e-inφ，然后在整个单位圆上进行积分，可以得到声场系数的表达式，即：
式中：rm为高阶麦克风的半径，bv(krm)为模态强度或模态幅度，p(k,rm,φqq′)为高阶麦克风采集到的局部声场的声压。
步骤b：计算模态强度。根据高阶麦克风的不同阵列设计，选择不同的模态强度。
式中：r0为无限长圆柱形刚性挡板的半径，r0≤rm，hv(krm)为第一类v阶汉克尔函数，(·)′为求导运算符号。
步骤c：空间采样。对步骤a中的积分进行空间采样，利用有限的求和来精确逼近连续傅里叶逆变换，可以得到离散傅里叶逆变换的形式，即：
式中：q′为麦克风的数量，aq′为第q′个麦克风的采样权重。
为了确保指数函数的正交性质，因此采样权重aq′需要满足以下条件，即：
步骤d：确定麦克风数量。声场系数是一个周期函数，且周期为2π,所以声场系数的最大角频率为n/2π。根据香农采样定理可知，只要采样频率大于等于两倍的最大频率，就可以利用其离散采样点准确计算出声场系数，因此声场系数的采样角频率为ζs≥n/π。局部声场的截断阶数为v阶，即该声场中的圆柱谐波函数的阶数最高为v阶，所以至少需要2v 1个采样点才能准确计算出声场系数。因此，将麦克风数量限制为最小数量，并采用等角度采样方案，即：
由此，可以得到局部声场系数，即：
附图2所示的是开放阵列的0-3阶模态强度关于kr的变化曲线。从图中可以看出，对于一些特定的kr值和声场截断阶数n，开放式阵列的模态强度的幅值很低，因此，开放式阵列经常会产生病态矩阵。由于实际测量中模态强度的倒数将放大在那些频率点的麦克风接收到的噪声，使得声场重构系统的鲁棒性大大降低，并且麦克风输出信号的信噪比(signaltonoiseratio,snr)会大幅度减小，特别是在输入信号的信噪比比较低的情况下，声场系数的测量值会出现很大的偏差，从而影响整个空间声场的重构精度。
附图3所示的是刚性阵列的0-3阶模态强度关于kr的变化曲线。从图中可以看出，刚性阵列可以很好的避免这一问题。刚性阵列是开放式阵列的改进设计，刚性阵列本质上是一个圆形阵列，且该阵列围绕着一个半径为r0(r0≤rm)的无限长圆柱形刚性挡板。由刚体的散射效应，刚性阵列的贝塞尔零点问题更少，所以相比于开放式阵列，刚性阵列具有更好的稳健性，并且刚性阵列的衍射作用可以提高信号在低频处的信噪比，使阵列在较宽的频率范围内获得很好的效果
步骤二：局部声场截断。该模块的主要功能是根据声场频率计算局部声场的截断阶数，并将不存在的高阶谐波分量正则化。
根据高阶麦克风的阶数，可以确定v阶麦克风的半径，即：
式中：c为声速，fmax为v阶麦克风的最大可重构声场频率。
在一个半径为rm、声场最高频率为fmax的区域内的圆柱谐波函数的阶数最高限制为v，并且当声场频率高于fmax时，声场的谐波阶数将增加到v 1，即：
式中：fact(v 1)为激活第v 1阶圆柱谐波需要的声场频率。
假设v阶麦克风的最大可重构声场频率为fmax，重构声场频率为f(f<fmax)，只有当fact(v)≤f<fact(v 1)时，才能可以激活第v(v＝1,…,v)阶圆柱谐波。当声场频率为f时，声场中实际存在的或者已经被激活的圆柱谐波的阶数仅为：
步骤三：构造转换矩阵。该模块的主要功能是根据高阶麦克风的阶数和阵列的拓扑结构来构造转换矩阵。
贝塞尔函数的加法定理：
根据贝塞尔函数的加法定理可以得到相对于局部原点的声压，即：
由此，可以得到局部声场系数与全局声场系数的关系：
将上式写为矩阵形式，即：
a＝tb
式中：a为局部声场系数矢量，b为全局声场系数矢量，t为系数转换矩阵，其中
转换矩阵t的2-范数条件数定义为：
式中：||·||2为矩阵的l2范数，
附图4所示的是三阶麦克风的转换矩阵条件数关于频率的变化曲线。从图中可以看出，对于三阶麦克风，当阵列的拓扑结构为半径可变的圆形阵列时，转换矩阵的条件数远小于其他两种拓扑结构，说明半径可变的圆形阵列具有非常好的鲁棒性。
附图5所示的是半径可变的圆形阵列拓扑结构的转换矩阵条件数关于频率的变化曲线。从图中可以看出，对于半径可变的圆形阵列，三阶麦克风的转换矩阵条件数远小于一阶和二阶麦克风，说明高阶麦克风的阶数越高，其鲁棒性越好。
步骤四：谐波表示。该模块的主要功能是根据声场的频率和重构区域的半径计算截断阶数，然后用全局声场系数表征重构后的声场。
白噪声条件下的局部声场系数测量值为：
式中：p(k,rm,φqq′)为真实声压信号，nq′为白噪声。
白噪声是一个复数，即：
nq′＝nq′(re) inq′(im)
式中：nq′(re)为白噪声的实部，nq′(im)为白噪声的虚部。
将局部声场系数的测量值分解为信号分量和噪声分量，即：
式中：
由此，可以得到白噪声条件下全局声场系数矢量：
式中：
利用声场的谐波表达式，可以得到重构声场的声压分布，即：
式中：n为全局声场的截断阶数。
附图6所示的是平面波的目标声场系数、重构声场系数及其重构系数误差曲线。附图7所示的是点声源的目标声场系数、重构声场系数及其重构系数误差曲线。从图中可以看出，基于高阶麦克风阵列重构的平面波和点声源的声场系数误差较小，并且两条重构误差曲线中没有出现峰值点，说明高阶麦克风阵列可以减小贝塞尔函数零值的影响，减小系数重构误差。
附图8所示的是平面波的目标声场和重构声场效果图，附图9所示的是点声源的目标声场和重构声场效果图，附图10所示的是平面波和点声源的声场重构误差。从图中可以看出，基于高阶麦克风阵列的声场重构技术可以较好的恢复目标声场的声压分布，并且有效地减小了白噪声的影响和贝塞尔零值问题。
附图11所示的是不同信噪比环境下，平面波声场重构的归一化均方误差与声源信号频率之间的变化关系。附图12所示的是不同信噪比环境下，点声源声场重构的归一化均方误差与声源信号频率之间的变化关系。从图中可以看出，在设计上限频率为1400hz的情况下，当信噪比高于20db时，平面波声场的重构误差基本上低于-10db，当信噪比高于30db时，点声源声场的重构误差基本上低于-10db。误差曲线中仅出现了个别峰值点，表明基于高阶麦克风阵列的声场重构技术具有抑制贝塞尔零值问题的能力。
附图13所示的是平面波重构声场的白噪声增益关于声场频率之间的变化关系。附图14所示的是点声源重构声场的白噪声增益关于声场频率之间的变化关系。从图中可以看出，对于平面波声场，在整个频率范围内白噪声增益的值大约在-10到10db之间，对于点声源声场，白噪声增益的值约在11到27db之间，并且每个圆柱谐波对应的白噪声增益略有不同。此外，在f＝469hz和f＝937hz处，白噪声增益出现了2次幅值下降，这是因为在这两个频率处高阶麦克风对高阶谐波分量作了正则化处理，以降低白噪声增益为代价，提高了声场重构准确度。
综上，本发明改进了传统hoa声场重构技术测量声场系数的过程，利用高阶麦克风采集局部声场声压，然后利用局部声场截断减小全局声场系数的重构误差，高阶麦克风以及半径可变的圆形阵列拓扑结构可以减小转换矩阵的条件数，提高阵列的鲁棒性，并且降低贝塞尔零值的影响，从而提高了空间声场的重构准确度。
1.一种基于高阶麦克风阵列的声场重构方法，采用高阶麦克风阵列采集声场信号；
其特征在于步骤如下：
步骤1：计算局部声场系数
步骤1a：根据高阶麦克风的不同阵列设计，选择不同的模态强度：
式中：r0为无限长圆柱形刚性挡板的半径，r0≤rm，jv(·)为v阶贝塞尔函数，hv(·)为第一类v阶汉克尔函数，(·)′为求导运算符号，k为波数；
步骤1b：根据声压的圆柱谐波展开式，利用指数函数的正交性，得到声场系数的表达式，即：
式中：rm为高阶麦克风的半径，φqq′为麦克风极坐标下的极角，e为自然指数，p(k,rm,φqq′)为高阶麦克风采集到的局部声场的声压，nq′为白噪声；
步骤1c：对步骤1b中的积分进行空间采样，利用有限的求和来精确逼近连续傅里叶逆变换，得到离散傅里叶逆变换的形式，并且根据香农采样定理，将麦克风数量限制为最小数量，并采用等角度采样方案，即：
式中：q′为麦克风的数量，aq′为第q′个麦克风的采样权重，v为高阶麦克风的阶数；
局部声场系数改写为：
步骤2：局部声场截断
根据高阶麦克风的阶数，确定v阶麦克风的半径，即：
式中：c为声速，fmax为v阶麦克风的最大可重构声场频率；
当声场频率为f时，声场中实际存在的或者已经被激活的圆柱谐波的阶数仅为：
将局部声场系数
步骤3：根据高阶麦克风的阶数和阵列的拓扑结构来构造转换矩阵，即：
式中：
步骤4：计算全局声场系数向量：
式中：
得到
步骤5：利用声场的谐波表达式，得到重构声场的声压分布，即：
式中：
