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

基于雅克比旋转联合对角化的空时频方位估计方法
【专利摘要】基于雅克比旋转联合对角化的空时频方位估计方法，本发明涉及空时频方位估计方法。本发明是要解决现有的方位估计方法未利用信源的频域信息以及仅利用单个时频分布点信息的问题。该方法是通过一、构造阵列接收信号的空时频分布矩阵；二、得到K个不同阵列接收信号空时频分布矩阵组；三、抽取出元素和四、构造第k个时频点的矢量；五、得到新的矩阵Re(GGH)；六、求得特征向量v；七、求得雅克比旋转角度θ和φ；八、构造新的矩阵组九、得到D′XX(tk,fk)；十、得到阵列流型矩阵A(θ)；十一、提取联合特征值；十二、求取阵列输出功率；十三、确定声源来波方向等步骤实现的。本发明应用于空时频方位估计领域。
【专利说明】
基于雅克比旋转联合对角化的空时频方位估计方法
技术领域
[0001] 本发明涉及空时频方位估计方法，特别涉及基于雅克比旋转联合对角化的空时频 方位估计方法。
【背景技术】
[0002] 空间信源方位估计(Direction Of Arrival :D0A)是阵列信号处理技术研究的重 要分支，广泛应用于雷达、声纳、医学成像等国民经济及国防建设领域。经过近50年的发展， 广大研究学者提出了诸多经典的方位估计方法，包括最大熵谱法(Maximum Entropy :ME)、 最大似然法(Maximum Likelihood:ML)以及特征子空间法等(H.Krim and M.Viberg.Two decades of array signal processing research. IEEE signal processing magazine. 1996,13(4) :67-94.)。随着阵列信号处理技术研究的不断深入，基于空时相关矩 阵组和高阶累积量组联合对角化的方位估计方法受到专家学者的广泛关注，该类方法利用 二阶统计量或高阶统计量信息，对相关矩阵组进行联合对角化处理，在一定程度上改善了 方位估计精度和性能。（参见CNKI检索情况，主题词:联合对角化方位估计）（[1 ]宋海岩，朴 胜春.基于高阶累积量矩阵组正交联合对角化的高分辨方位估计方法.电子与信息学报， 2010年04期.[2]蒋飚，朱埜，孙长瑜.基于空时相关阵联合对角化的高分辨方位估计.哈尔 滨工程大学学报，2005年06期.[3]赵龙龙，宋海岩.联合对角化技术在空间谱估计中的应 用.鱼雷技术，2012年05期.[4]宋海岩.具有高稳健性的浅海目标方位估计方法研究.哈尔 滨工程大学博士论文，2011. [5]余华，幸高翔，刘忠.时空相关矩阵联合对角化对谱分辨的 影响.控制工程，2009年05期）。
[0003] 然而，值得注意的是，现有的基于空时相关矩阵组和高阶累积量组联合对角化方 位估计方法仅利用空域-时域二维统计信息，并未利用信源的频域信息。如若能同时有效利 用空域-时域-频域三维信息，可有望进一步改善现有方法的方位估计性能。
[0004] 1998年，Belouchrani和Amin首先提出 了空时频分布（Spatial Time-Frequency Distribution: STFD )矩阵的概念，并将其引入非平稳信号盲源分离技术领域 (A.Belouchrani and M.G.Amin.Blind source separation based on time-frequency signal representations. IEEE Trans. Signal Processing.1998,46(11):2888-2897·)〇 随后，该研究团队又进一步将STFD矩阵拓展并应用于方位估计领域，结合信号子空间与噪 声子空间正交的性质，提出了时频分析多重信号子空间正交方法（Time-Frequency Mus i c )，开启了时频分析类方位估计算法研究的先河（[I ] A . Be louchran i and M.G.Amin.Time-frequency MUSIC.IEEE Signal Processing Lett..1999,6(5):109-110. [2]Adel Belouchrani,Moeness G·Amin，Nadege Thirion-Moreau,and Yimin D . Zhang . Source Separation and Localization Using Time-Frequency Distributions.IEEE SIGNAL PROCESSING 1^6厶21陬，2013,30(6):97-107)。近年来，国内 外科研工作者也开展了卓有成效的研究，取得了一定的成果。在国内，哈尔滨工程大学李秀 坤研究团队将时频分析方位估计方法应用于水下矢量传感器阵列，联合利用声压和振速信 息，对水下目标进行有效方位估计(尹世梅.矢量传感器阵时频联合方位估计.哈尔滨工程 大学硕士论文，2009;李秀坤，尹世梅，李婷婷.矢量水听器阵时频MUSIC算法研究.声学技 术，2010年04期）。西安电子科技大学以及哈尔滨工业大学等研究团队分别对时频子空间 (MUSIC)类算法进行仿真研究和性能分析，克服了传统子空间测向方法的缺陷，提高了测向 能力。（[1]应武，杨绍全.基于时频子空间的DOA估计.航天电子对抗，2004年03期；[2]谢宝 钢.基于时频-MUSIC的波达角估计研究.哈尔滨工业大学硕士论文，2008; [3]黄克骥.时频 分析方法在阵列信号处理中的应用.电子科技大学博士论文，2004. [4]袁泉，石昭祥.一种 基于空间时频分布的波达方向估计方法.探测与控制学报，2007年02期）。在国外，Gershman 等研究学者将STro矩阵的应用由窄带信号拓展到宽带信号，结合相干信号处理方法，实现 宽带调频信号的有效方位估计(A.B.Gershman and M.G.Amin.Coherent wideband DOA estimation of multiple FM signals using spatial time-frequency distributions . IEEE International Conference on Acoustics , Speech,and Signal Processing · 2000，5:3065-3068 ·) aGhofrani研究团队利用匹配跟踪(Matching Pursuit: MP )技术对STFD矩阵进行有效估计和修正，改善了现有算法的方位估计性能 (S.Ghofrani.Matching pursuit decomposition for high-resolution direction of arrival.Multidimensional systems and signal processing·2015，26(3):693-716·)〇 Khodja等专家学者针对阵列误差以及噪声干扰的影响，对时频分析类方位估计算法的性能 进行了详细的分析，为该类算法在实际工程中的应用奠定了理论基础（M. Khodja， A.Belouchrani,and K.Abed-Meraim.Performance analysis for time-frequency MUSIC algorithm in presence of both additive noise and array calibration errors.EURASIP Journal on advances in signal processing.2012,94.)；
[0005] 然而，纵观整个时频分析类方位估计算法的研究现状，现有方法均采用多重信号 子空间正交法(子空间类方法)作为处理器，在计算过程中需要估计信号子空间和噪声子空 间，不可避免由于矩阵特征值分解而带来的计算误差较大以及运算效率较低等缺点;此外， 现有方法大多数仅利用信号的单个时频脊点构造 STro矩阵，以代替传统的阵列相关矩阵， 无法融合多个时频脊点的信息，存在估计精度低、旁瓣起伏大等缺点。

【发明内容】

[0006] 本发明的目的是为了解决现有的方位估计方法大多数仅利用空域-时域二维统计 信息，并未利用信源的频域信息以及利用空时频分布矩阵进行方位估计的方法大多数仅利 用单个时频分布点信息的问题，而提出的基于雅克比旋转联合对角化的空时频方位估计方 法。
[0007] 上述的发明目的是通过以下技术方案实现的：
[0008] 步骤一、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵Dxx(t，f);
[0009] 其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，
[0010]
[0011] Xl ·)表示X( ·)的复共辄矩阵;1为中间变量； (2)
[0012]步骤二、选择第k个时频点（tk，fk)，根据公式(2)构造第k个阵列接收信号空时频分 布矩阵Dxx(tk，fk);从而得到K个不同阵列接收信号空时频分布矩阵组;其中，k=l，2,3,…， K;K为时频点的总数；
[0013]步骤三、选择正整数ρ和q作为空时频分布矩阵组Dxx(tk，f k)的坐标索引号，抽取出 坐标索引号对应的坐标分别为化^八化^八^^彡和^^^根据坐标化^彡确定元素^5^.、 坐标(P，q)确定元素、坐标(q，P)确定元素以及坐标(q，q)确定元素其中，Kp <N；1 <q<N；
[0014] 步骤四、根据从空时频分布矩阵组Dxx(tk，fk)抽取出的元I
1和 构造第k个时频点（tk，f k)的矢i，进而 构造一个3XK维矩阵G=[gi，···，gk，'"gK];
[0015] 步骤五、将矩阵G与矩阵G复共辄转置矩阵Gh相乘，然后取实部运算，得到新的矩阵 Re(GGH);其中，（·）H表示复共辄转置，Re( ·)表示取实部运算；
[0016] 步骤六、对Re (GGh)进行特征值分解，求得最大的特征值和最大特征值对应的特征 向量V;
[0017] 步骤七、由V与雅克比旋转角度0和φ的关系V= [cos20，-sin20cos Φ，-sin20sin Φ ]求得雅克比旋转角度9和Φ ;其中，Θ定义为幅度雅克比旋转角度，Φ定义为相位雅克比 旋转角度；
[0018]步骤八、根据雅克比旋转角度Θ和φ构造复余弦-正弦矩阵
该矩阵g即为雅克比旋转矩阵;根据从空时频分布矩阵组Dxx(t k， fk)，抽取出的元素4#、^^、和构造2X2维矩阵纟」
利用复余弦-正
(4) 1 步骤九、利用W代朁矩阵Dxx(tk，fk)中第p行第p列的元素 利用^ ;)代替矩阵 Dxx(tk，fk)中第p行第q列的元素 α[Μ)，利用^^代替矩阵Dxx(tk，fk)中第q行第p列的元素 ，利用#}代替矩阵Dxx(tk，fk)中第q行第q列的元素得到D'xxU^fk)如公式(8)所 示：
[0022]
(B)
[0023] 步骤十、将P = I重复步骤三~九时，将q = l重复步骤三~九，直至q = N为止；
[0024] 将p = 2重复步骤三~九时，将q = 2重复步骤三~九，直至q = N为止；
[0025] 将p = 3重复步骤三~九时，将q = 3重复步骤三~九，直至q = N为止；
[0026].
[0027] .
[0028] .
[0029] 将P = N-I重复步骤三~九时，将q = N重复步骤三~九，直至q = N为止；
[0030] 将p = N重复步骤三~九时，将q = N重复步骤三~九;从而得到p=l~N循环计算后 的最终的DSxdfk);再将所有雅克比旋转矩阵g相乘即得到阵列流型矩阵Α(θ);
[0031]步骤十一、由步骤十最终得到的DSxUhfk)即为联合对角化后的对角矩阵，记为 05办1^1〇;将所有对角矩阵05办1^1〇累加求和并取平均，构造新的"_1对角矩阵0，提取 对角矩阵D中的对角元素 λη即为联合特征值；
[0032]新的N X N维对角矩阵D具体为：
[0033]
(5)
[0034] 其中，η = 1，2，···，Ν; λη为第η个对角元素；
[0035] 步骤十二、预设入射方位角度Θ，生成导向矢量a(0)，并由步骤十得到的阵列导向 矢量Α(θ)和步骤得到的联合特征值λ η，求取阵列输出功率孟卬)： (())
[0036]
[0037] 其中，Αη(θ)表示Α(θ)的第η列；
[0038]步骤十三、设置合适的方位角扫描步长，重复步骤十二，将扫描角度Θ和对应的阵 列输出功率(的绘制空间谱图；通过比较输出功率谱图，由功率谱图谱峰位置确定声源 来波方向。
[0039]发明效果
[0040]本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对 由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变 处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。
[0041] 本发明涉及一种基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得 稳健的高分辨方位估计结果。
[0042] 纵观整个时频分析类方位估计算法的研究现状，本发明与现有技术相比较具以下 优点：（1)本发明首次采用雅克比旋转技术对STro矩阵进行联合对角化处理，而以往技术大 多数仅利用信号的单个时频脊点构造 STFD矩阵。相比之下，本
【发明内容】
具有估计结果精确、 收敛速度快(量化数据)等优点。（从图2可以看出，在信噪比5dB条件下，本发明方法的主峰 宽度较窄（约为1° )，旁瓣较低(约为_27dB);而传统方法的主峰宽度较宽（约为5° )，旁瓣较 高(约为-15dB)，同时还伴随有伪峰的出现。从图3可以看出，在信噪比5dB条件下，本发明方 法具有明显的两个主峰，且宽度较窄（约为3°)，主峰之间的凹槽较低(约为-8dB);而传统方 法的主峰宽度较宽(约为5° )，旁瓣较高(约为-5dB)，同时还伴随有伪峰的出现。）
[0043] (2)本发明首次将最小方差无畸变处理器(MVDR)应用于时频分析类方位估计领 域，而以往技术均采用多重信号子空间正交法(子空间类方法)作为处理器。相比之下，本发 明内容勿需估计信号子空间和噪声子空间，避免了由于矩阵特征值分解而带来的误差。
[0044] 以上阐述的本发明与现有技术的不同和区别也即本
【发明内容】
的创新点所在。综上 所述，本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多 个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理 器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。该方法可有效提 高现有时频分布类方位估计算法在导向矢量误差、低信噪比、小快拍等条件下的稳健性能， 进而获得高分辨率的方位估计空间谱图，正确反映空间信源的波达方位信息。
【附图说明】
[0045] 图1为【具体实施方式】二提出的阵列信号模型示意图；
[0046] 图2为【具体实施方式】一提出的单信源方位估计结果示意图；
[0047] 图3为【具体实施方式】一提出的双相干信源方位估计结果。
【具体实施方式】
【具体实施方式】 [0048] 一:本实施方式的基于雅克比旋转联合对角化的空时频方位估计方 法，具体是按照以下步骤制备的：
[0049] 步骤一、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵Dxx(t，f);
[0050] 其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，利用两个下角 标X的表达形式则是与空时频分布矩阵Dxx(t，f)的表达式(2)有关，表达式(2)中用到了信号 X和X的共辄#两组数据，其中浐又可以由对导到，所以用Dxx(t，f)表示空时频分布矩阵
[00511
(2)
[0052] X*( ·)表示X( ·)的复共辄矩阵;1为中间变量；
[0053]步骤二、选择第k个时频点（tk，fk)，根据公式(2)构造第k个阵列接收信号空时频分 布矩阵Dxx(tk，fk);从而得到K个不同阵列接收信号空时频分布矩阵组;其中，k=l，2,3,…， κ;κ为时频点的总数；
[0054] 步骤三、选择正整数ρ和q作为空时频分布矩阵组Dxx(tk，fk)的坐标索引号，抽取出 坐标索引号对应的坐标分别为化^八化^八^^彡和^^^根据坐标化^彡确定元素^^、 坐标(P，q)确定元素<W)、坐标(q，P)确定元素以及坐标(q，q)确定元素其中， <N；1 <q<N；
[0055] 步骤四、根据从空时频分布矩阵组Dxx(tk，f k)抽取出的元素和 构造第k个时频点（tk，f k)的矢i
，进而 构造一个3XK维矩阵G=[gi，···，gk，'"gK];
[0056] 步骤五、将矩阵G与矩阵G复共辄转置矩阵Gh相乘，然后取实部运算(一般复数表示 为a+jb，a表示实部，b表示虚部，取实部运算就是把a取出来，具体表示为Re [a+jb ]= a，），得 到新的矩阵Re(GGH);其中，（·）H表示复共辄转置，Re( ·)表示取实部运算；
[0057] 步骤六、对Re(GGH)进行特征值分解，求得最大的特征值和最大特征值对应的特征 向量V;(特征值分解是线性代数及矩阵论中的重要内容，基本教材都有涉及;在matlab仿真 软件中由特定的函数eig直接求解可得最大的特征值和最大特征值对应的特征向量V) [0058] 步骤七、由V与雅克比旋转角度Θ和φ的关系V= [cos29，-sin29cos Φ，_sin29sin Φ ]求得雅克比旋转角度9和Φ ;其中，Θ定义为幅度雅克比旋转角度，φ定义为相位雅克比 旋转角度；
[0059]击3悪八、枏据雅克比旋转角度Θ和φ构造复余弦-正弦矩阵
，该矩阵g即为雅克比旋转矩阵;根据从空时频分布矩阵组Dxx(t k， fk)，抽取出的元素纪和构造 2X2维矩阵：
；利用复余弦-正 弦矩阵g对矩阵；
P的每一个矩阵进行变换，BP :
[0060] (4)
[0061 ] 中的每一个矩阵均为对 角矩阵；
[0062]步骤九、利用4〃>代替矩阵Dxx(tk，fk)中第ρ行第ρ列的元素《^，利用^^代替矩阵 Dxx(tk，fk)中第ρ行第q列的元素 af，利用代替矩阵Dxx(tk，fk)中第q行第ρ列的元素 利用#>代替矩阵Dxx(t k，fk)中第q行第q列的元素导到DSx(Wfk)如公式(8)所 示：
[0063]
(8)
[0064] 步骤十、将p = l重复步骤三~九时，将q = l重复步骤三~九，直至q = N为止；
[0065] 将p = 2重复步骤三~九时，将q = 2重复步骤三~九，直至q = N为止；
[0066] 将p = 3重复步骤三~九时，将q = 3重复步骤三~九，直至q = N为止；
[0067] .
[0068] .
[0069] .
[0070] 将P = N-I重复步骤三~九时，将q = N重复步骤三~九，直至q = N为止；
[0071] 将p = N重复步骤三~九时，将q = N重复步骤三~九;从而得到p=l~N循环计算后 的最终的D7 XX(tk，fk);再将所有雅克比旋转矩阵g相乘即得到联合对角化器U;其中，联合对 角化器U即为阵列流型矩阵Α(θ);
[0072]步骤十一、由步骤十最终得到的DSxUifk)即为联合对角化后的对角矩阵，记为 05办1^1〇;将所有对角矩阵05办1^1〇累加求和并取平均，构造新的"_1对角矩阵0，提取 对角矩阵D中的对角元素 λη即为联合特征值；
[0073]新的N X N维对角矩阵D具体为：
[0074]
(5)
[0075] 其中，η = 1，2，···，Ν; λη为第η个对角元素；
[0076] 步骤十二、预设入射方位角度Θ，生成导向矢量a(0)，并由步骤十得到的阵列导向 矢量Α( Θ)和步骤得到的联合特征值λη，求取阵列输出功率孟⑷)::
(6)
[0077]
[0078] 其中，Αη(θ)表示Α(θ)的第n列；
[0079] (TF:Time Frequency，表示时频；JD:Joint Diagonalization，表示联合对角化； MVDR:Minimum Variance Distortionless Response，表不最小方差无畸变处理器）
[0080] 步骤十三、设置合适的方位角扫描步长，重复步骤十二，将扫描角度Θ和对应的阵 列输出功率尸^孟(的绘制空间谱图（根据步骤十二表达式(6)可以看出，每一个扫描角度Θ 对应一个阵列输出功率/^孟⑷)，两者构成一一对应关系，利用matlab软件即可绘图）；通 过比较输出功率谱图即空间谱图，由功率谱图谱峰位置确定声源来波方向。
[0081]本实施方式效果：
[0082]本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对 由多个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变 处理器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。
[0083]本发明涉及一种基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得 稳健的高分辨方位估计结果。
[0084]纵观整个时频分析类方位估计算法的研究现状，本发明与现有技术具以下优点： (1)本发明首次采用雅克比旋转技术对STro矩阵进行联合对角化处理，而以往技术大多数 仅利用信号的单个时频脊点构造 STro矩阵。相比之下，本
【发明内容】
具有估计结果精确、收敛 速度快(量化数据)等优点。（从图2可以看出，在信噪比5dB条件下，本发明方法的主峰宽度 较窄（约为1° )，旁瓣较低(约为-27dB);而传统方法的主峰宽度较宽（约为5° )，旁瓣较高(约 为-15dB)，同时还伴随有伪峰的出现。从图3可以看出，在信噪比5dB条件下，本发明方法具 有明显的两个主峰，且宽度较窄（约为3°)，主峰之间的凹槽较低（约为-8dB);而传统方法的 主峰宽度较宽(约为5° )，旁瓣较高(约为_5dB)，同时还伴随有伪峰的出现。）
[0085] (2)本发明首次将最小方差无畸变处理器(MVDR)应用于时频分析类方位估计领 域，而以往技术均采用多重信号子空间正交法(子空间类方法)作为处理器。相比之下，本发 明内容勿需估计信号子空间和噪声子空间，避免了由于矩阵特征值分解而带来的误差。
[0086] 以上阐述的本发明与现有技术的不同和区别也即本
【发明内容】
的创新点所在。综上 所述，本发明在综合利用空域-时域-频域三维信息的基础上，通过雅克比旋转技术，对由多 个时频分布点构成的空时频分布矩阵组进行联合对角化处理，并结合最小方差无畸变处理 器，提出了一种新的基于雅克比旋转联合对角化的空时频方位估计方法。该方法可有效提 高现有时频分布类方位估计算法在导向矢量误差、低信噪比、小快拍等条件下的稳健性能， 进而获得高分辨率的方位估计空间谱图，正确反映空间信源的波达方位信息。
【具体实施方式】 [0087] 二:本实施方式与一不同的是:步骤一中根据阵列接 收信号X(t)，构造阵列接收信号的空时频分布矩阵D XX(t，f)具体为：
[0088]步骤一一、如图1所示，假设空间存在一均匀直线阵，阵元个数为N，阵元间距为d; 在均匀直线阵远场处存在M个窄带信源信号，且第m个窄带信源的入射角度为0m，则第m个窄 带信号相对应的导向矢量a(0 m)为：
[0089]
[0090] 其中，fm表示第m个窄带信号的频率，c表示空间信号传播速度，[· ]τ表示矩阵转 置，j表示虚数单位;y = …，M;e是自然常数;M为窄带信源信号的总个数；
[0091] 由此，阵列接收信号X(t)表示为：
[0092] X(t)=A(9)S(t) (1)
[0093] 其中，X(t) = [xi(t)，X2(t)，···，xn(t)…，XN(t)]T为阵列接收信号，A(Q) = IIa(Q1)，a (θ2)，…，a(0m)-_，a(0M)]为各窄带信号对应的导向矢量形成的阵列流型矩阵，S(t) = [S1 (t)，s2(t)，…，sm(t)，…，sM(t)]%M个窄带信源信号形成的源信号矢量;χ η⑴为X(t)中第η 个阵元的阵列接收信号;sm(t)M个窄带信源信号中第m个窄带信源信号；
[0094] 步骤一二、根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵Dxx(t， f)：
[0095]
(2)
[0096] 其中，·)表示X( ·)的复共辄矩阵;1为中间变量；
[0097]步骤一三、将式（1)代入式(2)，得：
[0098] Dxx(t,f)=A(0)Dss(t,f)AH(0) (3)
[0099] 其中，Dss(t，f)为信源信号的空时频分布矩阵，（·）H表示复共辄转置;AH( ·)表示 A( ·)复共轭鮮詈，
[0100]
[0101] 其中，Si ·)为s( ·)的复共辄矩阵。其它步骤及参数与【具体实施方式】一相同。
[0102]
【具体实施方式】三:本实施方式与【具体实施方式】一或二不同的是：步骤三中抽取出 坐标索引号对应的坐标分别为化^八化^八^^彡和^^^根据坐标化^彡的元素^:^^坐 标(P，q)的元素坐标(q，P)的元素以及坐标(q，q)的元素具体为：
[0103]公式(7)为第k个空时频分布矩阵D\\(tk，fk)，令矩阵Dxx(tk，fk)中的元素用符号ak 来表示，则该矩阵中第P行第P列的元素为即坐标(P，P)的元素为4W))，第P行第q列的 元素表示为即坐标(p，q)的元素为4 M))，第q行第p列的元素表示为即坐标(q，p) 的元素为)，第q行第q列的元素表示为(即坐标(q，q)的元素为《广>);
[0104]
(7)其它步骤及参 数与【具体实施方式】一或二相同。
【具体实施方式】 [0105] 四：本实施方式与一至三之一不同的是：步骤二所述 Dxx(tk，fk)=A(0)Dss(tk，fk)AH(0)。其它步骤及参数与一至三之一相同。
[0106] 【具体实施方式】五:本实施方式与【具体实施方式】一至四之一不同的是：步骤二所述
,其它步骤及参数 与【具体实施方式】一至四之一相同。
【具体实施方式】 [0107] 六:本实施方式与一至五之一不同的是：步骤十三所 述步长可根据实际情况自行设定，因此方位角空间扫描角度一般是-90°至90°，步长为_ 90° :1° :90°或-90° :0.01° :90°。其它步骤及参数与一至五之一相同。
【具体实施方式】 [0108] 七:本实施方式与一至六之一不同的是：步骤十三所 述功率谱图为利用matlab软件绘出的图即阵列输出功率谱图。其它步骤及参数与具体实施 方式一至六之一相同。
[0109] 采用以下实施例验证本发明的有益效果：
[0110] 实施例一：
[0111] 本实施例基于雅克比旋转联合对角化的空时频方位估计方法，具体是按照以下步 骤制备的：
[0112] 单信源方位估计分析
[0113] 基于雅克比旋转联合对角化的空时频方位估计方法，利用空域-时域-频域三维信 息，通过雅克比旋转技术，对空时频分布矩阵组进行联合对角化，可获得稳健的高分辨方位 估计结果，下面对仿真实例进行分析。
[0114] 实例参数设置如下：设空间存在一均匀直线阵，阵元个数为8且阵元间距为半波 长。空间信源频率为20Hz，入射方位角为10°。接收阵元采样快拍数为1024,信噪比为5dB。
[0115] 图2给出单信源条件下，常规MVDR处理器（简记为MVDR)、基于空时频分布矩阵的 MVDR处理器(简记为SIFD-MVDR)，以及基于雅克比旋转联合对角化的空时频方位估计方法 (STFD-JD-MVDR)的方位估计结果图。
[0116] 其中，STFD-TD-MVDR进行方位估计过稈中，在步骤H中牛成的对角矩阵D如下：
[0118] 分析图2可知:MVDR方法的主瓣较宽，并且旁瓣级较高;STFD-MVDR和STFD-JD-MVDR 具有较尖锐的主瓣，但STro-MVDR的旁瓣起伏较大，甚至出现伪峰（false peaks)。由此可 见，在单信源入射的条件下，本
【发明内容】
所提出的方法STFD-JD-MVDR具有更好的方位估计 性能。
[0119] 实例二:双相干信源方位估计分析
[0120] 为了更好的表明本发明方法STFD-JD-MVDR的优良性能，接下来对空间双相干信源 进行方位估计分析。
[0121] 实例参数设置如下：设空间存在一均匀直线阵，阵元个数为8且阵元间距为半波 长。空间存在两相干信源，频率均为20Hz，入射方位角分别为3°和10°。接收阵元采样快拍数 为1024,信噪比为10dB。
[0122] 图3给出双相干信源条件下，MVDR、STFD-MVDR、以及STFD-JD-MVDR的方位估计结果 图。
[0123] 其中，STFD-JD-MVDR进行方位估计过程中，在步骤^^一中生成的对角矩阵D如下：
[0125] 分析图3可知：在此仿真条件下，MVDR已不能有效分辨出空间双相干源，而STFD-MVDR和STFD-JD-MVDR均能分辨出空间双相干源。但相比之下STFD-JD-MVDR具有更尖锐的谱 峰和更低的旁瓣，而STro-MVDR旁瓣起伏较大，甚至出现伪峰(false peaks)。由此可见，在 双相干信源入射的条件下，本
【发明内容】
所提出的方STFD-JD-MVDR具有更好的方位估计性 能。
[0126] 本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域 技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于 本发明所附的权利要求的保护范围。
【主权项】
1.基于雅克比旋转联合对角化的空时频方位估计方法，其特征在于，该方法具体是按 照W下步骤进行的： 步骤一、根据阵列接收信号x(t)，构造阵列接收信号的空时频分布矩阵Dxx(t，f); 其中，t表示时间变量，f表示频率变量，下角标X表示阵列接收信号，稱 ·)表示x( ·)的复共辆矩阵;1为中间变量； 步骤二、选择第k个时频点（tk，fk)，根据公式(2)构造第k个阵列接收信号空时频分布矩 阵化x(tk，fk);从而得到K个不同阵列接收信号空时频分布矩阵组;其中，k=l，2,3,…，Κ;Κ 为时频点的总数； 步骤Ξ、选择正整数Ρ和q作为空时频分布矩阵组Dxx(tk，fk)的坐标索引号，抽取出坐标 索引号对应的坐标分别为(口，口）、（口，9)、^，口)和^，9);根据坐标(口，口)确定元素4^)、坐标 (P，q)确定元素、坐标(q，P)确定元素 W及坐标(q，q)确定元素皆"1 ;其中，1如卽； 1 <q<N； 步骤四、根据从空时频分布矩阵组化x(tk，fk)抽取出的元素和构 造第k个时频点（tk，fk)的矢i，.进而构造一 个3XK维矩阵G=[gl，···，阱，···gκ]; 步骤五、将矩阵G与矩阵G复共辆转置矩阵GH相乘，然后取实部运算，得到新的矩阵Re (GGH);其中，（.）H表示复共辆转置，Re( ·)表示取实部运算； 步骤六、对Re(GGH)进行特征值分解，求得最大的特征值和最大特征值对应的特征向量 V; 步骤屯、由V与雅克比旋转角度目和Φ的关系v= k〇s2目,-sin2目cos Φ ,-sin2目sin Φ ]求 得雅克比旋转角度θ和Φ ;其中，θ定义为幅度雅克比旋转角度，Φ定义为相位雅克比旋转角 度； 步骤八、根据雅克比旋转角度Θ和Φ构造复余弦-正弦矩降该 矩阵g即为雅克比旋转矩阵；根据从空时频分布矩阵组化x(tk，fk)，抽取出的元素 α严和却心构造2X2维矩阵注：利用复余弦-正弦矩阵g对矩阵组片的每一个矩阵进行变换，即：巧 进而构造新的矩阵组步骤九、利用4?'>代替矩阵Dxx(tk，fk)中第P行第P列的元素皆^，利用戸i代替矩阵DXX (tk，fk)中第P行第q列的元素皆刑用6户代替矩阵Dxx(tk，fk)中第q行第P列的元素命W， 利用6!'">代替矩阵Dxx(tk，fk)中第q行第q列的元素(护\得到iyxx(tk，fk)如公式(8)所示：步骤十、将P=1重复步骤Ξ~九时，将q=l重复步骤Ξ~九，直至q = N为止； 将P = 2重复步骤Ξ~九时，将q = 2重复步骤Ξ~九，直至q = N为止； 将P = 3重复步骤Ξ~九时，将q = 3重复步骤Ξ~九，直至q = N为止； 将p = N-l重复步骤Ξ~九时，将q = N重复步骤Ξ~九，直至q = N为止； 将P = N重复步骤Ξ~九时，将q = N重复步骤Ξ~九;从而得到p = l~N循环计算后的最 终的iyxx(tk，fk);再将所有雅克比旋转矩阵g相乘即得到阵列流型矩阵Α(θ); 步骤十一、由步骤十最终得到的l/xx(tk，fk)即为联合对角化后的对角矩阵，记为化s(tk, fk);将所有对角矩阵Dss(tk，fk)累加求和并取平均，构造新的NXN维对角矩阵D，提取对角矩 阵D中的对角元素 λη即为联合特征值； 新的ΝΧΝ维对角矩阵D具体为：巧 其中，n=l，2，…，N;λn为第n个对角元素； 步骤十二、预设入射方位角度Θ，生成导向矢量a(0)，并由步骤十得到的阵列导向矢量A (9)和步骤得到的联合特征值λη，求取阵列输出功率f芯.基(0);胸 其中，Αη(θ)表示Α(θ)的第η列； 步骤十Ξ、设置合适的方位角扫描步长，重复步骤十二，将扫描角度Θ和对应的阵列输 出功率巧'贯器(0)绘制空间谱图；通过比较输出功率谱图，由功率谱图谱峰位置确定声源来波 方向。2. 根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在 于:步骤一中根据阵列接收信号X(t)，构造阵列接收信号的空时频分布矩阵Dxx(t，f)具体 为： 步骤一一、假设空间存在一均匀直线阵，阵元个数为N，阵元间距为d;在均匀直线阵远 场处存在Μ个窄带信源信号，且第m个窄带信源的入射角度为θ。，则第m个窄带信号相对应的 导向矢量a(0m)为：其中，fm表示第m个窄带信号的频率，C表示空间信号传播速度，j表示虚数单位； /二麥…，M;e是自然常数;Μ为窄带信源信号的总个数； 由此，阵列接收信号X(t)表示为： X(t)=A(0)S(t) (1) 其中，X(t) = [xi(t)，X2(t)，…，Xn(t)…，XN(t) ]τ为阵列接收信号，A(目）=[a(目1)，a (曰2)，…，3(θη)···，3(θΜ)]为各窄带信号对应的导向矢量形成的阵列流型矩阵，s(t) = kl (t)，S2(t)，…，Sm(t)，…，SM(t)]%M个窄带信源信号形成的源信号矢量;Sm(t)M个窄带信源 信号中第m个窄带信源信号； 步骤一二、根据阵列接收信号x(t)，构造阵列接收信号的空时频分布矩阵Dxx(t，f):妈 其中，·)表示x( ·)的复共辆矩阵;1为中间变量； 步骤一 Ξ、将式(1)代入式(2)，得： Dxx(t，f)=A(目）Dss(t，f)AH(0) (3) 其中，Dss(t，f)为信源信号的空时频分布矩阵，（·）H表示复共辆转置;aH( ·)表示A( ·) 复共辆转置；其中，·)为S( ·)的复共辆矩阵。3. 根据权利要求1或2所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征 在于:步骤Ξ中抽取出坐标索引号对应的坐标分别为(口，口）、（口，9)、^，口)和^，9);根据坐 标(P，P)的元素坐标(P，q)的元素坐标(q，P)的元素 Κ及坐标(q，q)的元素 具体为： 公式(7)为第k个空时频分布矩阵Dxx(tk，fk)，则该矩阵中第Ρ行第Ρ列的元素为皆W，第Ρ 行第q列的元素表示为，第q行第P列的元素表示为如W1，第q行第q列的元素表示为冷W1;4. 根据权利要求1或2所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征 在于:步骤二所述Dxx(tk,fk) =A(目)Dss(tk,fk)AH(0)。5. 根据权利要求3所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在 于:步骤二所述 Dxx(tk，fk)=A(目)Dss(tk，fk)AH(0)。6. 根据权利要求1、2或5所述基于雅克比旋转联合对角化的空时频方位估计方法，其特 征在于:步骤二所述7. 根据权利要求4所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在 于:步骤二所述 8. 根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在 于:步骤十Ξ所述步长为-90° : 1° : 90°或-90° : 0.01° : 90°。9. 根据权利要求1所述基于雅克比旋转联合对角化的空时频方位估计方法，其特征在 于:步骤十Ξ所述功率谱图为利用matlab软件绘出的图即阵列输出功率谱图。
【文档编号】G01S3/802GK105842656SQ201610378776
【公开日】2016年8月10日
【申请日】2016年5月31日
【发明人】宋海岩, 秦进平, 刁鸣, 杨昌益, 高洪元, 刘伯胜, 时洁
【申请人】黑龙江工程学院