CN105842656A - 基于雅克比旋转联合对角化的空时频方位估计方法 - Google Patents

Abstract

基于雅克比旋转联合对角化的空时频方位估计方法，本发明涉及空时频方位估计方法。本发明是要解决现有的方位估计方法未利用信源的频域信息以及仅利用单个时频分布点信息的问题。该方法是通过一、构造阵列接收信号的空时频分布矩阵；二、得到K个不同阵列接收信号空时频分布矩阵组；三、抽取出元素和四、构造第k个时频点的矢量；五、得到新的矩阵Re(GG^H)；六、求得特征向量v；七、求得雅克比旋转角度θ和φ；八、构造新的矩阵组九、得到D′_XX(t_k,f_k)；十、得到阵列流型矩阵A(θ)；十一、提取联合特征值；十二、求取阵列输出功率；十三、确定声源来波方向等步骤实现的。本发明应用于空时频方位估计领域。

Description

基于雅克比旋转联合对角化的空时频方位估计方法

技术领域

本发明涉及空时频方位估计方法，特别涉及基于雅克比旋转联合对角化的空时频方位估计方法。

背景技术

空间信源方位估计(Direction Of Arrival：DOA)是阵列信号处理技术研究的重要分支，广泛应用于雷达、声纳、医学成像等国民经济及国防建设领域。经过近50年的发展，广大研究学者提出了诸多经典的方位估计方法，包括最大熵谱法(Maximum Entropy：ME)、最大似然法(Maximum Likelihood：ML)以及特征子空间法等(H.Krim and M.Viberg.Twodecades of array signal processing research.IEEE signal processingmagazine.1996,13(4):67-94.)。随着阵列信号处理技术研究的不断深入，基于空时相关矩阵组和高阶累积量组联合对角化的方位估计方法受到专家学者的广泛关注，该类方法利用二阶统计量或高阶统计量信息，对相关矩阵组进行联合对角化处理，在一定程度上改善了方位估计精度和性能。(参见CNKI检索情况，主题词：联合对角化方位估计)([1]宋海岩，朴胜春.基于高阶累积量矩阵组正交联合对角化的高分辨方位估计方法.电子与信息学报，2010年04期.[2]蒋飚，朱埜，孙长瑜.基于空时相关阵联合对角化的高分辨方位估计.哈尔滨工程大学学报，2005年06期.[3]赵龙龙，宋海岩.联合对角化技术在空间谱估计中的应用.鱼雷技术，2012年05期.[4]宋海岩.具有高稳健性的浅海目标方位估计方法研究.哈尔滨工程大学博士论文，2011.[5]余华，幸高翔，刘忠.时空相关矩阵联合对角化对谱分辨的影响.控制工程，2009年05期)。

然而，值得注意的是，现有的基于空时相关矩阵组和高阶累积量组联合对角化方位估计方法仅利用空域-时域二维统计信息，并未利用信源的频域信息。如若能同时有效利用空域-时域-频域三维信息，可有望进一步改善现有方法的方位估计性能。

1998年，Belouchrani和Amin首先提出了空时频分布(Spatial Time-FrequencyDistribution：STFD)矩阵的概念，并将其引入非平稳信号盲源分离技术领域(A.Belouchrani and M.G.Amin.Blind source separation based on time-frequencysignal representations.IEEE Trans.Signal Processing.1998,46(11):2888-2897.)。随后，该研究团队又进一步将STFD矩阵拓展并应用于方位估计领域，结合信号子空间与噪声子空间正交的性质，提出了时频分析多重信号子空间正交方法(Time-FrequencyMusic)，开启了时频分析类方位估计算法研究的先河([1]A.Belouchrani andM.G.Amin.Time-frequency MUSIC.IEEE Signal Processing Lett..1999,6(5):109-110.[2]Adel Belouchrani，Moeness G.Amin，Nadège Thirion-Moreau，and YiminD.Zhang.Source Separation and Localization Using Time-FrequencyDistributions.IEEE SIGNAL PROCESSING MAGAZINE，2013，30(6)：97-107)。近年来，国内外科研工作者也开展了卓有成效的研究，取得了一定的成果。在国内，哈尔滨工程大学李秀坤研究团队将时频分析方位估计方法应用于水下矢量传感器阵列，联合利用声压和振速信息，对水下目标进行有效方位估计(尹世梅.矢量传感器阵时频联合方位估计.哈尔滨工程大学硕士论文，2009；李秀坤，尹世梅，李婷婷.矢量水听器阵时频MUSIC算法研究.声学技术，2010年04期)。西安电子科技大学以及哈尔滨工业大学等研究团队分别对时频子空间(MUSIC)类算法进行仿真研究和性能分析，克服了传统子空间测向方法的缺陷，提高了测向能力。([1]应武，杨绍全.基于时频子空间的DOA估计.航天电子对抗，2004年03期；[2]谢宝钢.基于时频-MUSIC的波达角估计研究.哈尔滨工业大学硕士论文，2008；[3]黄克骥.时频分析方法在阵列信号处理中的应用.电子科技大学博士论文，2004.[4]袁泉，石昭祥.一种基于空间时频分布的波达方向估计方法.探测与控制学报，2007年02期)。在国外，Gershman等研究学者将STFD矩阵的应用由窄带信号拓展到宽带信号，结合相干信号处理方法，实现宽带调频信号的有效方位估计(A.B.Gershman and M.G.Amin.Coherent wideband DOAestimation of multiple FM signals using spatial time-frequencydistributions.IEEE International Conference on Acoustics,Speech,and SignalProcessing.2000,5:3065-3068.)。Ghofrani研究团队利用匹配跟踪(Matching Pursuit：MP)技术对STFD矩阵进行有效估计和修正，改善了现有算法的方位估计性能(S.Ghofrani.Matching pursuit decomposition for high-resolution direction ofarrival.Multidimensional systems and signal processing.2015,26(3):693-716.)。Khodja等专家学者针对阵列误差以及噪声干扰的影响，对时频分析类方位估计算法的性能进行了详细的分析，为该类算法在实际工程中的应用奠定了理论基础(M.Khodja,A.Belouchrani,and K.Abed-Meraim.Performance analysis for time-frequency MUSICalgorithm in presence of both additive noise and array calibrationerrors.EURASIP Journal on advances in signal processing.2012,94.)；

然而，纵观整个时频分析类方位估计算法的研究现状，现有方法均采用多重信号子空间正交法(子空间类方法)作为处理器，在计算过程中需要估计信号子空间和噪声子空间，不可避免由于矩阵特征值分解而带来的计算误差较大以及运算效率较低等缺点；此外，现有方法大多数仅利用信号的单个时频脊点构造STFD矩阵，以代替传统的阵列相关矩阵，无法融合多个时频脊点的信息，存在估计精度低、旁瓣起伏大等缺点。

发明内容

本发明的目的是为了解决现有的方位估计方法大多数仅利用空域-时域二维统计信息，并未利用信源的频域信息以及利用空时频分布矩阵进行方位估计的方法大多数仅利用单个时频分布点信息的问题，而提出的基于雅克比旋转联合对角化的空时频方位估计方法。

上述的发明目的是通过以下技术方案实现的：

步骤一、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)；

其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，

$D_{XX}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}X(t+m+l)X^{*}(t+m-l)e^{-j4\mathrm{\π}fl}---\left(2\right)$

X^*(·)表示X(·)的复共轭矩阵；l为中间变量；

步骤二、选择第k个时频点(t_k,f_k)，根据公式(2)构造第k个阵列接收信号空时频分布矩阵D_XX(t_k,f_k)；从而得到K个不同阵列接收信号空时频分布矩阵组；其中，k＝1,2,3,…,K；K为时频点的总数；

步骤三、选择正整数p和q作为空时频分布矩阵组D_XX(t_k,f_k)的坐标索引号，抽取出坐标索引号对应的坐标分别为(p,p)、(p,q)、(q,p)和(q,q)；根据坐标(p,p)确定元素坐标(p,q)确定元素坐标(q,p)确定元素以及坐标(q,q)确定元素其中，1≤p≤N；1≤q≤N；

步骤四、根据从空时频分布矩阵组D_XX(t_k,f_k)抽取出的元素和构造第k个时频点(t_k,f_k)的矢量进而构造一个3×K维矩阵G＝[g₁,…,g_k,…g_K]；

步骤五、将矩阵G与矩阵G复共轭转置矩阵G^H相乘，然后取实部运算，得到新的矩阵Re(GG^H)；其中，(·)^H表示复共轭转置，Re(·)表示取实部运算；

步骤六、对Re(GG^H)进行特征值分解，求得最大的特征值和最大特征值对应的特征向量v；

步骤七、由v与雅克比旋转角度θ和φ的关系v＝[cos2θ,-sin2θcosφ,-sin2θsinφ]求得雅克比旋转角度θ和φ；其中，θ定义为幅度雅克比旋转角度，φ定义为相位雅克比旋转角度；

步骤八、根据雅克比旋转角度θ和φ构造复余弦-正弦矩阵该矩阵g即为雅克比旋转矩阵；根据从空时频分布矩阵组D_XX(t_k,f_k)，抽取出的元素 和构造2×2维矩阵组利用复余弦-正弦矩阵g对矩阵组中的每一个矩阵进行变换，即：

$\left[\begin{array}{cc}{b}_k^{\left(pp\right)}& b_k^{\left(pq\right)}\\ b_k^{\left(qp\right)}& b_k^{\left(qq\right)}\end{array}\right]=g^T\left[\begin{array}{cc}{a}_k^{\left(pp\right)}& a_k^{\left(pq\right)}\\ a_k^{\left(qp\right)}& a_k^{\left(qq\right)}\end{array}\right]g---\left(4\right)$

进而构造新的矩阵组

步骤九、利用代替矩阵D_XX(t_k,f_k)中第p行第p列的元素利用代替矩阵D_XX(t_k,f_k)中第p行第q列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第p列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第q列的元素得到D′_XX(t_k,f_k)如公式(8)所示：

步骤十、将p＝1重复步骤三～九时，将q＝1重复步骤三～九，直至q＝N为止；

将p＝2重复步骤三～九时，将q＝2重复步骤三～九，直至q＝N为止；

将p＝3重复步骤三～九时，将q＝3重复步骤三～九，直至q＝N为止；

.

.

.

将p＝N-1重复步骤三～九时，将q＝N重复步骤三～九，直至q＝N为止；

将p＝N重复步骤三～九时，将q＝N重复步骤三～九；从而得到p＝1～N循环计算后的最终的D′_XX(t_k,f_k)；再将所有雅克比旋转矩阵g相乘即得到阵列流型矩阵A(θ)；

步骤十一、由步骤十最终得到的D′_XX(t_k,f_k)即为联合对角化后的对角矩阵，记为D_SS(t_k,f_k)；将所有对角矩阵D_SS(t_k,f_k)累加求和并取平均，构造新的N×N维对角矩阵D，提取对角矩阵D中的对角元素λ_n即为联合特征值；

新的N×N维对角矩阵D具体为：

$D=\frac{1}{K}\underset{k=1}{\overset{K}{\Σ}}{D}_{SS}(t_k,f_k)---\left(5\right)$

其中，n＝1,2,…,N；λ_n为第n个对角元素；

步骤十二、预设入射方位角度θ，生成导向矢量a(θ)，并由步骤十得到的阵列导向矢量A(θ)和步骤得到的联合特征值λ_n，求取阵列输出功率

$P_{MVDR}^{TF-JD}\left(\mathrm{\θ}\right)=\frac{1}{a^H\left(\mathrm{\θ}\right)\underset{n=1}{\overset{N}{\Σ}}\left(\frac{1}{{\mathrm{\λ}}_n}{A}_n\left(\mathrm{\θ}\right)A_n^H\left(\mathrm{\θ}\right)\right)a\left(\mathrm{\θ}\right)}---\left(6\right)$

其中，A_n(θ)表示A(θ)的第n列；

步骤十三、设置合适的方位角扫描步长，重复步骤十二，将扫描角度θ和对应的阵列输出功率绘制空间谱图；通过比较输出功率谱图，由功率谱图谱峰位置确定声源来波方向。

发明效果

本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。

本发明涉及一种基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得稳健的高分辨方位估计结果。

纵观整个时频分析类方位估计算法的研究现状，本发明与现有技术相比较具以下优点：(1)本发明首次采用雅克比旋转技术对STFD矩阵进行联合对角化处理，而以往技术大多数仅利用信号的单个时频脊点构造STFD矩阵。相比之下，本发明内容具有估计结果精确、收敛速度快(量化数据)等优点。(从图2可以看出，在信噪比5dB条件下，本发明方法的主峰宽度较窄(约为1°)，旁瓣较低(约为-27dB)；而传统方法的主峰宽度较宽(约为5°)，旁瓣较高(约为-15dB)，同时还伴随有伪峰的出现。从图3可以看出，在信噪比5dB条件下，本发明方法具有明显的两个主峰，且宽度较窄(约为3°)，主峰之间的凹槽较低(约为-8dB)；而传统方法的主峰宽度较宽(约为5°)，旁瓣较高(约为-5dB)，同时还伴随有伪峰的出现。)

(2)本发明首次将最小方差无畸变处理器(MVDR)应用于时频分析类方位估计领域，而以往技术均采用多重信号子空间正交法(子空间类方法)作为处理器。相比之下，本发明内容勿需估计信号子空间和噪声子空间，避免了由于矩阵特征值分解而带来的误差。

以上阐述的本发明与现有技术的不同和区别也即本发明内容的创新点所在。综上所述，本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。该方法可有效提高现有时频分布类方位估计算法在导向矢量误差、低信噪比、小快拍等条件下的稳健性能，进而获得高分辨率的方位估计空间谱图，正确反映空间信源的波达方位信息。

附图说明

图1为具体实施方式二提出的阵列信号模型示意图；

图2为具体实施方式一提出的单信源方位估计结果示意图；

图3为具体实施方式一提出的双相干信源方位估计结果。

具体实施方式

具体实施方式一：本实施方式的基于雅克比旋转联合对角化的空时频方位估计方法，具体是按照以下步骤制备的：

步骤一、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)；

其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，利用两个下角标X的表达形式则是与空时频分布矩阵D_XX(t,f)的表达式(2)有关，表达式(2)中用到了信号X和X的共轭X^*两组数据，其中X^*又可以由X得到，所以用D_XX(t,f)表示空时频分布矩阵

$D_{XX}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}X(t+m+l)X^{*}(t+m-l)e^{-j4\mathrm{\π}fl}---\left(2\right)$

X^*(·)表示X(·)的复共轭矩阵；l为中间变量；

步骤二、选择第k个时频点(t_k,f_k)，根据公式(2)构造第k个阵列接收信号空时频分布矩阵D_XX(t_k,f_k)；从而得到K个不同阵列接收信号空时频分布矩阵组；其中，k＝1,2,3,…,K；K为时频点的总数；

步骤三、选择正整数p和q作为空时频分布矩阵组D_XX(t_k,f_k)的坐标索引号，抽取出坐标索引号对应的坐标分别为(p,p)、(p,q)、(q,p)和(q,q)；根据坐标(p,p)确定元素坐标(p,q)确定元素坐标(q,p)确定元素以及坐标(q,q)确定元素其中，1≤p≤N；1≤q≤N；

步骤四、根据从空时频分布矩阵组D_XX(t_k,f_k)抽取出的元素和构造第k个时频点(t_k,f_k)的矢量进而构造一个3×K维矩阵G＝[g₁,…,g_k,…g_K]；

步骤五、将矩阵G与矩阵G复共轭转置矩阵G^H相乘，然后取实部运算(一般复数表示为a+jb，a表示实部，b表示虚部，取实部运算就是把a取出来，具体表示为Re[a+jb]＝a，)，得到新的矩阵Re(GG^H)；其中，(·)^H表示复共轭转置，Re(·)表示取实部运算；

步骤六、对Re(GG^H)进行特征值分解，求得最大的特征值和最大特征值对应的特征向量v；(特征值分解是线性代数及矩阵论中的重要内容，基本教材都有涉及；在matlab仿真软件中由特定的函数eig直接求解可得最大的特征值和最大特征值对应的特征向量v)

步骤七、由v与雅克比旋转角度θ和φ的关系v＝[cos2θ,-sin2θcosφ,-sin2θsinφ]求得雅克比旋转角度θ和φ；其中，θ定义为幅度雅克比旋转角度，φ定义为相位雅克比旋转角度；

步骤八、根据雅克比旋转角度θ和φ构造复余弦-正弦矩阵该矩阵g即为雅克比旋转矩阵；根据从空时频分布矩阵组D_XX(t_k,f_k)，抽取出的元素 和构造2×2维矩阵组利用复余弦-正弦矩阵g对矩阵组中的每一个矩阵进行变换，即：

$\left[\begin{array}{cc}{b}_k^{\left(pp\right)}& b_k^{\left(pq\right)}\\ b_k^{\left(qp\right)}& b_k^{\left(qq\right)}\end{array}\right]=g^T\left[\begin{array}{cc}{a}_k^{\left(pp\right)}& a_k^{\left(pq\right)}\\ a_k^{\left(qp\right)}& a_k^{\left(qq\right)}\end{array}\right]g---\left(4\right)$

进而构造新的矩阵组即矩阵组中的每一个矩阵均为对角矩阵；

步骤九、利用代替矩阵D_XX(t_k,f_k)中第p行第p列的元素利用代替矩阵D_XX(t_k,f_k)中第p行第q列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第p列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第q列的元素得到D′_XX(t_k,f_k)如公式(8)所示：

步骤十、将p＝1重复步骤三～九时，将q＝1重复步骤三～九，直至q＝N为止；

将p＝2重复步骤三～九时，将q＝2重复步骤三～九，直至q＝N为止；

将p＝3重复步骤三～九时，将q＝3重复步骤三～九，直至q＝N为止；

.

.

.

将p＝N-1重复步骤三～九时，将q＝N重复步骤三～九，直至q＝N为止；

将p＝N重复步骤三～九时，将q＝N重复步骤三～九；从而得到p＝1～N循环计算后的最终的D′_XX(t_k,f_k)；再将所有雅克比旋转矩阵g相乘即得到联合对角化器U；其中，联合对角化器U即为阵列流型矩阵A(θ)；

步骤十一、由步骤十最终得到的D′_XX(t_k,f_k)即为联合对角化后的对角矩阵，记为D_SS(t_k,f_k)；将所有对角矩阵D_SS(t_k,f_k)累加求和并取平均，构造新的N×N维对角矩阵D，提取对角矩阵D中的对角元素λ_n即为联合特征值；

新的N×N维对角矩阵D具体为：

$D=\frac{1}{K}\underset{k=1}{\overset{K}{\Σ}}{D}_{SS}(t_k,f_k)---\left(5\right)$

其中，n＝1,2,…,N；λ_n为第n个对角元素；

步骤十二、预设入射方位角度θ，生成导向矢量a(θ)，并由步骤十得到的阵列导向矢量A(θ)和步骤得到的联合特征值λ_n，求取阵列输出功率

$P_{MVDR}^{TF-JD}\left(\mathrm{\θ}\right)=\frac{1}{a^H\left(\mathrm{\θ}\right)\underset{n=1}{\overset{N}{\Σ}}\left(\frac{1}{{\mathrm{\λ}}_n}{A}_n\left(\mathrm{\θ}\right)A_n^H\left(\mathrm{\θ}\right)\right)a\left(\mathrm{\θ}\right)}---\left(6\right)$

其中，A_n(θ)表示A(θ)的第n列；

(TF：Time Frequency，表示时频；JD：Joint Diagonalization，表示联合对角化；MVDR：Minimum Variance Distortionless Response，表示最小方差无畸变处理器)

步骤十三、设置合适的方位角扫描步长，重复步骤十二，将扫描角度θ和对应的阵列输出功率绘制空间谱图(根据步骤十二表达式(6)可以看出，每一个扫描角度θ对应一个阵列输出功率两者构成一一对应关系，利用matlab软件即可绘图)；通过比较输出功率谱图即空间谱图，由功率谱图谱峰位置确定声源来波方向。

本实施方式效果：

本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。

本发明涉及一种基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得稳健的高分辨方位估计结果。

纵观整个时频分析类方位估计算法的研究现状，本发明与现有技术具以下优点：(1)本发明首次采用雅克比旋转技术对STFD矩阵进行联合对角化处理，而以往技术大多数仅利用信号的单个时频脊点构造STFD矩阵。相比之下，本发明内容具有估计结果精确、收敛速度快(量化数据)等优点。(从图2可以看出，在信噪比5dB条件下，本发明方法的主峰宽度较窄(约为1°)，旁瓣较低(约为-27dB)；而传统方法的主峰宽度较宽(约为5°)，旁瓣较高(约为-15dB)，同时还伴随有伪峰的出现。从图3可以看出，在信噪比5dB条件下，本发明方法具有明显的两个主峰，且宽度较窄(约为3°)，主峰之间的凹槽较低(约为-8dB)；而传统方法的主峰宽度较宽(约为5°)，旁瓣较高(约为-5dB)，同时还伴随有伪峰的出现。)

(2)本发明首次将最小方差无畸变处理器(MVDR)应用于时频分析类方位估计领域，而以往技术均采用多重信号子空间正交法(子空间类方法)作为处理器。相比之下，本发明内容勿需估计信号子空间和噪声子空间，避免了由于矩阵特征值分解而带来的误差。

以上阐述的本发明与现有技术的不同和区别也即本发明内容的创新点所在。综上所述，本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。该方法可有效提高现有时频分布类方位估计算法在导向矢量误差、低信噪比、小快拍等条件下的稳健性能，进而获得高分辨率的方位估计空间谱图，正确反映空间信源的波达方位信息。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤一中根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)具体为：

步骤一一、如图1所示，假设空间存在一均匀直线阵，阵元个数为N，阵元间距为d；在均匀直线阵远场处存在M个窄带信源信号，且第m个窄带信源的入射角度为θ_m，则第m个窄带信号相对应的导向矢量a(θ_m)为：

$a\left({\mathrm{\θ}}_m\right)={\[1,e^{j2{\mathrm{\πf}}_m{\mathrm{dsin\θ}}_m/c},...,e^{j2{\mathrm{\πf}}_m(N-1){\mathrm{dsin\θ}}_m/c}\]}^T;$

其中，f_m表示第m个窄带信号的频率，c表示空间信号传播速度，[·]^T表示矩阵转置，j表示虚数单位；m＝1,…,M；e是自然常数；M为窄带信源信号的总个数；

由此，阵列接收信号X(t)表示为：

X(t)＝A(θ)S(t) (1)

其中，X(t)＝[x₁(t),x₂(t),…,x_n(t)…,x_N(t)]^T为阵列接收信号，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_m)…,a(θ_M)]为各窄带信号对应的导向矢量形成的阵列流型矩阵，S(t)＝[s₁(t),s₂(t),…,s_m(t)，…,s_M(t)]^T为M个窄带信源信号形成的源信号矢量；x_n(t)为X(t)中第n个阵元的阵列接收信号；s_m(t)M个窄带信源信号中第m个窄带信源信号；

步骤一二、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)：

$D_{XX}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}X(t+m+l)X^{*}(t+m-l)e^{-j4\mathrm{\π}fl}---\left(2\right)$

其中，X^*(·)表示X(·)的复共轭矩阵；l为中间变量；

步骤一三、将式(1)代入式(2)，得：

D_XX(t,f)＝A(θ)D_SS(t,f)A^H(θ) (3)

其中，D_SS(t,f)为信源信号的空时频分布矩阵，(·)^H表示复共轭转置；A^H(·)表示A(·)复共轭转置；

$D_{SS}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}S(t+m+l)S^{*}(t+m-l)e^{-j4\mathrm{\π}fl}$

其中，S^*(·)为S(·)的复共轭矩阵。其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤三中抽取出坐标索引号对应的坐标分别为(p,p)、(p,q)、(q,p)和(q,q)；根据坐标(p,p)的元素坐标(p,q)的元素坐标(q,p)的元素以及坐标(q,q)的元素具体为：

公式(7)为第k个空时频分布矩阵D_XX(t_k,f_k)，令矩阵D_XX(t_k,f_k)中的元素用符号a_k来表示，则该矩阵中第p行第p列的元素为(即坐标(p,p)的元素为)，第p行第q列的元素表示为(即坐标(p,q)的元素为)，第q行第p列的元素表示为(即坐标(q,p)的元素为)，第q行第q列的元素表示为(即坐标(q,q)的元素为)；

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤二所述D_XX(t_k,f_k)＝A(θ)D_SS(t_k,f_k)A^H(θ)。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤二所述其它步骤及参数与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤十三所述步长可根据实际情况自行设定，因此方位角空间扫描角度一般是-90°至90°，步长为-90°：1°：90°或-90°：0.01°：90°。其它步骤及参数与具体实施方式一至五之一相同。

具体实施方式七：本实施方式与具体实施方式一至六之一不同的是：步骤十三所述功率谱图为利用matlab软件绘出的图即阵列输出功率谱图。其它步骤及参数与具体实施方式一至六之一相同。

采用以下实施例验证本发明的有益效果：

实施例一：

本实施例基于雅克比旋转联合对角化的空时频方位估计方法，具体是按照以下步骤制备的：

单信源方位估计分析

基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得稳健的高分辨方位估计结果，下面对仿真实例进行分析。

实例参数设置如下：设空间存在一均匀直线阵，阵元个数为8且阵元间距为半波长。空间信源频率为20Hz，入射方位角为10°。接收阵元采样快拍数为1024，信噪比为5dB。

图2给出单信源条件下，常规MVDR处理器(简记为MVDR)、基于空时频分布矩阵的MVDR处理器(简记为STFD-MVDR)，以及基于雅克比旋转联合对角化的空时频方位估计方法(STFD-JD-MVDR)的方位估计结果图。

其中，STFD-JD-MVDR进行方位估计过程中，在步骤十一中生成的对角矩阵D如下：

$D=\left[\begin{array}{cccccccc}1002.09& 0& 0& 0& 0& 0& 0& 0\\ 0& 1025.80& 0& 0& 0& 0& 0& 0\\ 0& 0& 1012.38& 0& 0& 0& 0& 0\\ 0& 0& 0& 1001.66& 0& 0& 0& 0\\ 0& 0& 0& 0& 1031.53& 0& 0& 0\\ 0& 0& 0& 0& 0& 1030.51& 0& 0\\ 0& 0& 0& 0& 0& 0& 1017.29& 0\\ 0& 0& 0& 0& 0& 0& 0& 996.40\end{array}\right]$

分析图2可知：MVDR方法的主瓣较宽，并且旁瓣级较高；STFD-MVDR和STFD-JD-MVDR具有较尖锐的主瓣，但STFD-MVDR的旁瓣起伏较大，甚至出现伪峰(false peaks)。由此可见，在单信源入射的条件下，本发明内容所提出的方法STFD-JD-MVDR具有更好的方位估计性能。

实例二：双相干信源方位估计分析

为了更好的表明本发明方法STFD-JD-MVDR的优良性能，接下来对空间双相干信源进行方位估计分析。

实例参数设置如下：设空间存在一均匀直线阵，阵元个数为8且阵元间距为半波长。空间存在两相干信源，频率均为20Hz，入射方位角分别为3°和10°。接收阵元采样快拍数为1024，信噪比为10dB。

图3给出双相干信源条件下，MVDR、STFD-MVDR、以及STFD-JD-MVDR的方位估计结果图。

其中，STFD-JD-MVDR进行方位估计过程中，在步骤十一中生成的对角矩阵D如下：

$D=\left[\begin{array}{ccccccc}2211.96& 0& 0& 0& 0& 0& 0\\ 0& 2376.11& 0& 0& 0& 0& 0\\ 0& 0& 2453.43& 0& 0& 0& 0\\ 0& 0& 0& 2446.63& 0& 0& 0\\ 0& 0& 0& 0& 2453.43& 0& 0\\ 0& 0& 0& 0& 0& 2376.12& 0\\ 0& 0& 0& 0& 0& 0& 2211.96\end{array}\right]$

分析图3可知：在此仿真条件下，MVDR已不能有效分辨出空间双相干源，而STFD-MVDR和STFD-JD-MVDR均能分辨出空间双相干源。但相比之下STFD-JD-MVDR具有更尖锐的谱峰和更低的旁瓣，而STFD-MVDR旁瓣起伏较大，甚至出现伪峰(false peaks)。由此可见，在双相干信源入射的条件下，本发明内容所提出的方STFD-JD-MVDR具有更好的方位估计性能。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (9)

1.基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于，该方法具体是按照以下步骤进行的：

步骤一、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)；

其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，

$D_{XX}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}X(t+m+l)X^{*}(t+m-l)e^{-j4\mathrm{\π}fl}---\left(2\right)$

X^*(·)表示X(·)的复共轭矩阵；l为中间变量；

步骤二、选择第k个时频点(t_k,f_k)，根据公式(2)构造第k个阵列接收信号空时频分布矩阵D_XX(t_k,f_k)；从而得到K个不同阵列接收信号空时频分布矩阵组；其中，k＝1,2,3,…,K；K为时频点的总数；

步骤三、选择正整数p和q作为空时频分布矩阵组D_XX(t_k,f_k)的坐标索引号，抽取出坐标索引号对应的坐标分别为(p,p)、(p,q)、(q,p)和(q,q)；根据坐标(p,p)确定元素坐标(p,q)确定元素坐标(q,p)确定元素以及坐标(q,q)确定元素其中，1≤p≤N；1≤q≤N；

步骤四、根据从空时频分布矩阵组D_XX(t_k,f_k)抽取出的元素和构造第k个时频点(t_k,f_k)的矢量进而构造一个3×K维矩阵G＝[g₁,…,g_k,…g_K]；

步骤五、将矩阵G与矩阵G复共轭转置矩阵G^H相乘，然后取实部运算，得到新的矩阵Re(GG^H)；其中，(·)^H表示复共轭转置，Re(·)表示取实部运算；

步骤六、对Re(GG^H)进行特征值分解，求得最大的特征值和最大特征值对应的特征向量v；

步骤七、由v与雅克比旋转角度θ和φ的关系v＝[cos2θ,-sin2θcosφ,-sin2θsinφ]求得雅克比旋转角度θ和φ；其中，θ定义为幅度雅克比旋转角度，φ定义为相位雅克比旋转角度；

步骤八、根据雅克比旋转角度θ和φ构造复余弦-正弦矩阵该矩阵g即为雅克比旋转矩阵；根据从空时频分布矩阵组D_XX(t_k,f_k)，抽取出的元素 和构造2×2维矩阵组利用复余弦-正弦矩阵g对矩阵组中的每一个矩阵进行变换，即：

$\left[\begin{array}{cc}{b}_k^{\left(pp\right)}& b_k^{\left(pq\right)}\\ b_k^{\left(qp\right)}& b_k^{\left(qq\right)}\end{array}\right]=g^T\left[\begin{array}{cc}{a}_k^{\left(pp\right)}& a_k^{\left(pq\right)}\\ a_k^{\left(qp\right)}& a_k^{\left(qq\right)}\end{array}\right]g---\left(4\right)$

进而构造新的矩阵组

步骤九、利用代替矩阵D_XX(t_k,f_k)中第p行第p列的元素利用代替矩阵DXX(t_k,f_k)中第p行第q列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第p列的元素利用代替矩阵D_XX(t_k,f_k)中第q行第q列的元素得到D′_XX(t_k,f_k)如公式(8)所示：

步骤十、将p＝1重复步骤三～九时，将q＝1重复步骤三～九，直至q＝N为止；

将p＝2重复步骤三～九时，将q＝2重复步骤三～九，直至q＝N为止；

将p＝3重复步骤三～九时，将q＝3重复步骤三～九，直至q＝N为止；

.

.

.

将p＝N-1重复步骤三～九时，将q＝N重复步骤三～九，直至q＝N为止；

将p＝N重复步骤三～九时，将q＝N重复步骤三～九；从而得到p＝1～N循环计算后的最终的D′_XX(t_k,f_k)；再将所有雅克比旋转矩阵g相乘即得到阵列流型矩阵A(θ)；

步骤十一、由步骤十最终得到的D′_XX(t_k,f_k)即为联合对角化后的对角矩阵，记为D_SS(t_k,f_k)；将所有对角矩阵D_SS(t_k,f_k)累加求和并取平均，构造新的N×N维对角矩阵D，提取对角矩阵D中的对角元素λ_n即为联合特征值；

新的N×N维对角矩阵D具体为：

$D=\frac{1}{K}\underset{k=1}{\overset{K}{\Σ}}{D}_{SS}(t_k,f_k)---\left(5\right)$

其中，n＝1,2,…,N；λ_n为第n个对角元素；

步骤十二、预设入射方位角度θ，生成导向矢量a(θ)，并由步骤十得到的阵列导向矢量A(θ)和步骤得到的联合特征值λ_n，求取阵列输出功率

$P_{MVDR}^{TF-JD}\left(\mathrm{\θ}\right)=\frac{1}{a^H\left(\mathrm{\θ}\right)\underset{n=1}{\overset{N}{\Σ}}\left(\frac{1}{{\mathrm{\λ}}_n}{A}_n\left(\mathrm{\θ}\right)A_n^H\left(\mathrm{\θ}\right)\right)a\left(\mathrm{\θ}\right)}---\left(6\right)$

其中，A_n(θ)表示A(θ)的第n列；

步骤十三、设置合适的方位角扫描步长，重复步骤十二，将扫描角度θ和对应的阵列输出功率绘制空间谱图；通过比较输出功率谱图，由功率谱图谱峰位置确定声源来波方向。

2.根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤一中根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)具体为：

步骤一一、假设空间存在一均匀直线阵，阵元个数为N，阵元间距为d；在均匀直线阵远场处存在M个窄带信源信号，且第m个窄带信源的入射角度为θ_m，则第m个窄带信号相对应的导向矢量a(θ_m)为：

$a\left({\mathrm{\θ}}_m\right)={\[1,e^{j2{\mathrm{\πf}}_m{\mathrm{dsin\θ}}_m/c},...,e^{j2{\mathrm{\πf}}_m(N-1){\mathrm{dsin\θ}}_m/c}\]}^T;$

其中，f_m表示第m个窄带信号的频率，c表示空间信号传播速度，j表示虚数单位；m＝1,…,M；e是自然常数；M为窄带信源信号的总个数；

由此，阵列接收信号X(t)表示为：

X(t)＝A(θ)S(t) (1)

其中，X(t)＝[x₁(t),x₂(t),…,x_n(t)…,x_N(t)]^T为阵列接收信号，A(θ)＝[a(θ₁),a(θ₂),…,a(θ_m)…,a(θ_M)]为各窄带信号对应的导向矢量形成的阵列流型矩阵，S(t)＝[s₁(t),s₂(t),…,s_m(t)，…,s_M(t)]^T为M个窄带信源信号形成的源信号矢量；s_m(t)M个窄带信源信号中第m个窄带信源信号；

步骤一二、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵D_XX(t,f)：

$D_{XX}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}X(t+m+l)X^{*}(t+m-l)e^{-j4\mathrm{\π}fl}---\left(2\right)$

其中，X^*(·)表示X(·)的复共轭矩阵；l为中间变量；

步骤一三、将式(1)代入式(2)，得：

D_XX(t,f)＝A(θ)D_SS(t,f)A^H(θ) (3)

其中，D_SS(t,f)为信源信号的空时频分布矩阵，(·)^H表示复共轭转置；A^H(·)表示A(·)复共轭转置；

$D_{SS}(t,f)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}S(t+m+l)S^{*}(t+m-l)e^{-j4\mathrm{\π}fl}$

其中，S^*(·)为S(·)的复共轭矩阵。

3.根据权利要求1或2所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤三中抽取出坐标索引号对应的坐标分别为(p,p)、(p,q)、(q,p)和(q,q)；根据坐标(p,p)的元素坐标(p,q)的元素坐标(q,p)的元素以及坐标(q,q)的元素具体为：

公式(7)为第k个空时频分布矩阵D_XX(t_k,f_k)，则该矩阵中第p行第p列的元素为第p行第q列的元素表示为第q行第p列的元素表示为第q行第q列的元素表示为

4.根据权利要求1或2所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤二所述D_XX(t_k,f_k)＝A(θ)D_SS(t_k,f_k)A^H(θ)。

5.根据权利要求3所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤二所述D_XX(t_k,f_k)＝A(θ)D_SS(t_k,f_k)A^H(θ)。

6.根据权利要求1、2或5所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤二所述 $D_{SS}(t_k,f_k)=\underset{l=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\underset{m=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}S(t_k+m+l)S^{*}(t_k+m-l)e^{-j4{\mathrm{\πf}}_{k}l}.$

7.根据权利要求4所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤二所述

8.根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤十三所述步长为-90°：1°：90°或-90°：0.01°：90°。

9.根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于：步骤十三所述功率谱图为利用matlab软件绘出的图即阵列输出功率谱图。