最小二乘矩阵滤波器设计与性能分析

徐驰，韩磊，张书第，韩东

（海军大连舰艇学院，辽宁大连 116018）

0 引言

矩阵滤波技术源于短数据滤波。传统的FIR数字滤波器适用于数据长度较大的信号滤波，当信号数据长度较短时，由于FIR滤波器阶数的制约，不能得到期望的滤波效果，矩阵滤波器能针对短数据长度信号设计，可在信号长度较短时得到所需的通阻带响应效果［1－4］。根据矩阵滤波器的特点，它不仅可以用于短数据滤波，还可应用于阵列数据预处理。在阵列数据用于目标定向或定位算法之前滤波，可滤除或保留某一方位或扇面的信号及噪声，能有效提高目标检测能力及定位精度［5－10］。

在矩阵滤波器的设计方法中，Vaccaro［1－2］将滤波器对离散化的通带响应误差、阻带响应和过渡带振幅作为目标函数及约束条件构造凸规划问题设计滤波矩阵;韩东［2］将Vaccaro的方法改进，得到了连续方式的解表达式，不仅提高了设计效率，也减少了计算量;Zhu［3－4］和鄢社锋［6－7］设计了恒定阻带抑制的矩阵滤波器，并将矩阵滤波器用于目标方位估计的预处理;自适应预滤波方法由Hassanien［8］提出，所构造最优化问题的目标函数中包含阵列数据协方差矩阵，也即包含输入信号在各频率的能量信息，利用此信息实现自适应滤波。以上这些矩阵滤波器的设计算法都存在一个共同的问题，即需要利用复杂的软件或迭代算法才能给出矩阵滤波器的最优解，不适用于大阵元的滤波器设计问题。

本文在最小二乘矩阵滤波器设计的基础上，利用奇异值分解给出了滤波器的简化表达式及总体误差，方法简单，设计效率高。通过改变通带阻带位置、通带阻带带宽及阵元数等条件，分析得出这些因素对最小二乘矩阵滤波器性能的影响趋势，为矩阵滤波器的设计和应用提供参考。

1 最小二乘矩阵滤波器设计

考虑等间隔均匀线列阵，接收阵列数据为X（k） =AS（k）+N（k），设计一个M×M维矩阵滤波器H对均匀线列阵数据进行阵元域滤波，滤波输出为:

式中:AH=HA;NH（k）=HN（k）。

对全空间的波达方向进行细分，细分点数为D，该矩阵滤波器对不同方向到达阵元的方向向量响应为

为使该矩阵滤波器保留感兴趣方向的信号，滤除不感兴趣方向的噪声，并且限制过渡带的输出响应，可通过对不同的方向θi设计不同的k（θi），1≤i≤N值实现。

定义由全空间方向向量构成的阵列流形为X=［a（θ1，…，a（θD））］，维数为M×D，期望响应向量构成的期望阵列流形为Y=［k（θ1）a（θ1），…，k（θD）a（θD）］。所以，矩阵滤波器的设计问题就是求矩阵H，使HX=Y。

假设滤波器矩阵为

式中，hk=［hk1，hk2，…，hkM］，1≤k≤M是对应矩阵滤波器H的第k个行向量。令y=［h1，h2，…，hM］T，维数为M2×1，通过求解向量y就可以重构出矩阵滤波器H。

方向向量a（θj）经矩阵滤波器作用后，输出响应Ha（θj）=k（θj）a（θj），1≤j≤D可转化为如下的线性方程组:

通过矩阵变换及所构造的V和b，矩阵滤波器设计问题HX=Y转化为一个线性方程组求解问题，

该线性方程组的线性方程数目为DM，未知数个数为M2。该线性方程组的解y受方程组系数矩阵V的秩决定，当rank（V）=M2且线性方程组行数目DM大于未知数个数M2时，该方程组存在最小二乘解。

求矩阵V的秩。重排矩阵的各行不改变矩阵的秩，对矩阵V重排，可得:

式中，X为Vandermonde矩阵，当D＞M时，该矩阵列满秩，反之则行满秩。即

通常，对接收阵列数据进行矩阵滤波处理，对通带、过渡带、阻带的离散化向量数目总和D远大于接收阵阵元数目M，这说明矩阵V是列满秩的，即线性方程组Vy=b是超定的，可利用误差平方和最小求出该方程组的最小二乘解。

通过以西北师大的研究生新媒体运营为例进行的受众分析研究，我们发现了存在于高校新媒体运营中的一些问题，了解到受众与平台之间的重要联系，只有在明确受众群体的需求上，进行针对受众的内容创新是新媒体运营发展的关键。以上问题及相应的可行性运营建议，不仅适用于西北师大研究生会新媒体运营，对于甘肃各高校研究生新媒体运营平台也具有一定的借鉴意义。

式中:（·）H为矩阵共轭转置;（·）－1为方阵求逆矩阵。（VHV）－1VH是V的左伪逆矩阵，利用求得的y就可以重构出矩阵滤波器H。

由前面构造可知，V的维数为DM×M2，b的维数为DM×1，所以利用V的左伪逆与b的乘积求解y计算量非常大，影响了运算效率，结果的准确性受运算软件的制约。通过化简可得到矩阵滤波器最小二乘解的阵列流形矩阵运算形式。

式中，Bi（1≤i≤M）的维数为M×1，且

利用式（14）构造矩阵W=［B1，B2，…，BM］，维数为M×M。矩阵W与矩阵X，Y的关系为

（·）*表示矩阵共轭。合并式（13）和式（14），有

重排向量y，可得:

即矩阵滤波器H等于期望阵列流形Y与原阵列流形X的右伪逆XH（XXH）－1的乘积。可以通过期望阵列流形中波达方向系数k（θj），1≤j≤D的设置得到不同类型的矩阵滤波器。

2 奇异值分解滤波器误差分析

本节利用奇异值分解给出最小二乘矩阵滤波器的简化运算形式，并分析滤波器响应整体误差。上节给出了最小二乘矩阵滤波器的数值解，即矩阵滤波器H等于期望阵列流形Y与原阵列流形X的右伪逆XH（XXH）－1的乘积。现利用矩阵奇异值分解简化式（15）的滤波器设计。由前述已知rank（X）=M，假设阵列流形的奇异值分解为

其中，U∈CM×M，Q∈CD×D，U和Q都是酉矩阵，ΣX=［Σ1，0M×（D－M）］∈CM×D，Σ1=diag（r1，r2，…，rM）∈RM×M。由于rank（X）=M，因此ri≠0，i=1，…，M，利用奇异值分解可得:

滤波输出总体误差为:

3 仿真分析

最小二乘矩阵滤波器的性能主要与滤波器通带阻带位置、通带阻带带宽及阵元数等因素有关。通过计算机仿真，分析各因素对矩阵滤波器归一化总体均方误差的影响，为矩阵滤波器设计和使用提供参考。归一化总体均方误差值可以度量矩阵滤波器对空间各向信号的滤波效果，定义归一化总体均方误差为:

其中，ΩP和ΩS分别表示通带区域和阻带区域。最小二乘矩阵滤波器的总体误差与通带位置有关。下面通过计算机仿真分析研究最小二乘矩阵滤波器随通带位置变化的响应及归一化总体均方误差。

假设水听器线列阵阵元间隔为半波长且均匀分布，阵元个数为30。下面设计1个单通带矩阵滤波器，通带的带宽固定为30°。通带、阻带和过渡带的离散化采样间隔均设置为1°。给出通带中心位置为10°，15°，30°，45°，60°时矩阵滤波器的空域响应

从图1的仿真结果可知，在固定通带宽度和阵元数的情况下，通带处于正横位置时，最小二乘矩阵滤波器具有最佳的空域响应效果。随着通带中心位置向端首或端尾的外移，通带响应误差及阻带响应误差随之增大。图2给出了不同阵元数情况下，最小二乘矩阵滤波器总体均方误差随通带中心位置变化的趋势。由仿真结果可知，在相同阵元数情况下，随通带外移，归一化总体均方误差增大。在固定通带中心位置情况下，随阵元数的增加，归一化总体均方误差减少。

矩阵滤波器的响应效果不仅与通带位置有关，还与通带宽度有关，假设通带中心位置固定为0°，通带宽度逐渐增大，观察此时的矩阵滤波器响应效果和归一化总体均方误差。图3中，分别给出了阵元数为30时通带宽度为10°，30°，50°，70°和90°的矩阵滤波器响应和误差图。从图3可知，随通带宽度的增加，通带误差和阻带响应随之变化。图4给出了阵元数为20，40，60，80和100时，通带宽度变化的矩阵滤波器归一化总体均方误差，从仿真结果可知，随通带宽度的增加，归一化总体均方误差先增大后减小。

从矩阵滤波器的通带位置及通带宽度仿真可知，在不同阵元数情况下，矩阵滤波器所得的归一化总体均方误差有不同的结果，随阵元数的增加，总体误差都随之减小。现固定通带位置和通带宽度，观察阵元数对最小二乘矩阵滤波器性能影响的大小。假设水听器线列阵阵元间隔为半波长且均匀分布。设置通带扇面［－15°，15°］，阻带扇面为［－90°，－15°）∪（15°，90°］，离散化采样点间隔为0.1°。滤波器响应和误差效果如图5所示。

从图5的仿真结果可知，随阵元数目的增加，最小二乘矩阵滤波器通带和阻带的误差都在下降，当阵元个数从20变化到100时，矩阵滤波器的通带失真逐渐减小，通带误差从－15 dB变化到－20 dB;而阻带衰减也逐渐增大，从－20 dB左右增大到－30 dB左右。从图6的仿真结果可知，随阵元个数从10增加到100时，最小二乘矩阵滤波器的归一化总体均方误差呈抛物线下降。仿真结果说明，对于最小二乘矩阵滤波器的设计，阵元个数越多所得的滤波器响应效果越好，进而具有更好地阵列信号空域预处理能力。

4 结语

本文在最小二乘矩阵滤波器设计的基础上，利用奇异值分解给出了滤波器的简化表达式及总体误差。通过改变通带阻带位置、通带阻带带宽及阵元数等条件，分析得出这些因素对最小二乘矩阵滤波器性能的影响趋势，为矩阵滤波器的设计和应用提供参考。

［1］VACCARO R J，HARRISON B F．Optimal matrix-filter design［J］．IEEE Trans．Signal Processing，1996，44（3）: 705－709．

［2］HAN Dong，ZHANG Xin-hua．Optimal matrix filter design with application to filtering short data records［J］．IEEE Signal Processing Letters，2010，17（5）:521－524．

［3］ZHU Z W，WANG S，LEUNG H，et al．Matrix filter design using semi-infinite programming with application to DOA estimation［J］．IEEE Trans．Signal Processing，2000，48 （1）:267－271．

［4］WANG Shi，ZHUZhi-wen，LEUNGH．Semi-infinite optimization technique for the design of matrix filters［J］．Statistical Signal and Array Processing，1998，9（14）:204－207．

［5］韩东，章新华，康春玉，李军．零点约束矩阵滤波设计［J］．声学学报，2010，35（3）:353－358．

［6］VACCARO R J，CHHETRI A，HARRISON B F．Matrix filter design for passive sonar interference suppression［J］．J Acoustic Soc Am，2004，115（6）:3010－3020．

［7］鄢社锋，侯朝焕，马晓川．矩阵空域预滤波目标方位估计［J］．声学学报，2007，32（2）:151－157．

［8］鄢社锋，马远良．匹配场噪声抑制:广义空域滤波方法［J］．科学通报，2004，49（18）:1909－1912．

［9］MACINNES C S．Sourcelocalizationusingsubspace estimation and spatial filtering［J］．IEEE J．Ocean．Eng．2004，29（2）:488－497．

［10］HASSANIEN A，ELKADER S A，GERSHMAN A B，et al．Convex optimization based beam-space preprocessing with improved robustness against uut-of-sector sources［J］．IEEE Trans．Signal Processing，2006，54（5）:1587－1595．

［11］袁亚湘，孙文瑜．最优化理论与方法［M］．北京:科学出版社，1997．