一种采用时频分析的L型阵列二维波达方向估计算法的制作方法

本发明属于阵列信号处理领域，特别涉及一种基于时频分析的二维波达方向估计仿真方法，可解决时频谱混叠和空间邻近信源的方位估计。

背景技术：

二维波达方向(DOA)估计在阵列信号处理领域非常重要，它是通过处理空间中按一定方式排列的阵元上的接收数据，估计多个来波信号的方位信息，如方位角和仰角。二维DOA估计技术可应用于无线通信，雷达，声呐定位等应用，受到国内外学者广泛研究。研究过程中，发展出许多阵列模型，如均匀矩形阵列，双平行均匀线阵，均匀圆阵及L型阵，其中L型阵凭借较为简单的结构和参数设置，易于传统算法的移植，更因其能够提供更大的阵列孔径，使得该类型阵列在二维DOA估计效果上更佳，从而受到大量关注。

目前，已经有大量针对二维DOA估计的算法提出。子空间类算法主要将接收数据协方差矩阵分解为信号子空间和噪声子空间，其中多重信号分类MUSIC算法思路简单易于操作，但需要在二维方向上谱峰搜索；旋转不变ESPRIT算法无需谱峰搜索，但要求额外的角度配对过程。虽然子空间类算法具有超分辨率，但当信源角度邻近或信噪比较低时，估计性能会急速下降。另一类基于传播算子PM思想的算法仅需线性操作即可获得信号子空间，运算复杂度降低。经过不断发展，已出现许多具有高分辨率、抗噪性强、无需谱峰搜索和能自动配对的变形算法，然而这些算法要求信源个数需小于阵元个数，不适用于欠定情况，在工程应用中受限。除此之外，当多个信源DOA在二维空间十分邻近时，估计效果衰退。

实际应用中，多个待估计信源的时频谱混叠且空间邻近，即单个时间采样点上可能对应多个空间邻近信源的频率。在一维DOA估计中，通过引入时频分析，滤除时频域中多信号项和噪声项时频点，仅采用单信源时频点参与角度计算，已能达到优良的估计效果。虽然随着国内外学者对二维DOA估计研究的不断深入，各方面性能都得以大幅提升，但目前为止，仍未提出有效针对该类信源的二维DOA估计算法。此外，现有的大多数二维DOA估计算法在低信噪比和欠定情况下效果不佳。

技术实现要素：

本发明主要针对现有技术存在的问题，提出了一种经短时傅里叶变换(Short-time Fourier transform,STFT)将时域信号变换到时频域，采用二维方向上各信源对应的单信源时频点STFT值形成接收数据，再结合PM法则和阵列的旋转不变性进行二维角度自动配对的仿真方法。本发明着重解决时频谱混叠和空间邻近信源，对噪声具有高鲁棒性，且能够在阵元数小于待估计信源数情况下保持良好的估计效果。

本发明提出的仿真方法主要包括以下步骤：

1)两个均匀线阵正交摆放，形成xz平面上的L型阵列，并假设K个时频谱混叠信源以方位角θ和仰角φ入射到此阵列上。其中每个方向上有M个阵元，阵元间距为d(0＜d≤λ/2，λ为入射窄带信号波长)，θ代表入射信源与x轴正方向的夹角，φ代表入射信源与z轴正方向夹角；

2)对阵列二维方向上的接收数据分别采样，得到输出数据X(t),Z(t)，获得总体接收数据

3)采用时频分析工具STFT,将Y(t)变换到时频域W_y(t,f),W_x(t,f),W_z(t,f)分别是总体输出数据、x方向输出数据和z方向输出数据的STFT值；

4)利用合理的阈值和k-means聚类法选取时频域中每个时间采样点上有且只有一个频率的时频点。所形成的单信源时频点集用Φ^*表示。已知信源数K的情况下，再次利用k-means聚类法将点集Φ^*划分为K类，每类代表每个信源对应的单信源时频点集Φ_i,i＝1,…,K；

5)依据Φ_i,i＝1,…,K,建立单个信源的时频接收数据矩阵X_tf＝Y_tf(1:M,:),Z_tf＝Y_tf(M+1:2M,:)代表x,z方向上各信源的单信源时频点STFT值。D(i:j,:)代表矩阵D的第i到j行数据组成的子阵；

6)计算时频域互相关矩阵:

其中E{}代表期望运算符，H代表共轭转置运算符；

7)根据线性阵列的共轭对称特性，以及均匀线阵的旋转不变性，计算扩大阵列孔径后的新接收矩阵：

其中Y₁＝R_xz(:,1:M-1),Y₂＝R_xz(:,2:M)，Y₄＝J_MY₁^*J_M-1，J_M为副对角线全1的M×M维置换矩阵，*代表共轭运算符，T代表转置运算符；

8)基于步骤7中所获Y，通过传播算子法原则,求出传播算子P，并定义一个扩展的时频传播算子：

P_e＝[I_K,P]^H；

其中I_K为K×K维单位矩阵；

9)定义两个方位角选择矩阵G_x1,G_x2，两个仰角选择矩阵G_z1,G_z2，结合ESPRIT思想建立旋转不变等式，推导出方位角和仰角计算公式。

若信源数目K＝1，顺序执行步骤1)至9)；若K＞1,执行至第4)步后步骤5)-9)重复K次，共获得K对计算结果。

步骤4)所述的求取各信源在二维方向上的单信源时频点，按如下步骤进行：

步骤4a)、选择高能量信号项时频点：

式中，W_Y(t,v)表示总体输出数据在时间点t时，能获得的最大STFT值，且该值出现于频率v处。δ₀定义为高能量时频点选定阈值，取值范围为0.3～0.5，在此范围内取值大小仅影响选取的时频点个数，对最终估计效果无影响。满足条件的点集为Φ，计算Φ中每个点的空间矢量：

分别代表矩阵Y(t)各行，即各阵元接收数据的STFT值。采用k-means聚类算法将计算出的空间矢量分为K₀类，K₀＞K，得到K₀个方向矢量：

其中γ_k是第k类时频点集合，是集合γ_k中时频点的个数；方向矢量构成矩阵

步骤4b)、计算Φ中每个时频点上最可能出现的两个信源所对应的导向矢量a_n1,a_n2：

其中代表指向的噪声子空间的正交投影矩阵,其元素是矩阵A₀中的随机两列，I为单位矩阵，根据导向矢量a_n1,a_n2可以获得对应的两个STFT值：

其中A₂＝[a_n1,a_n2]，代表Moore-Penrose伪逆运算符；

步骤4c)、选取单信源时频点并分类，步骤如下：

单信源时频点上仅有信号项和噪声项，因此W₁,W₂的数量级有明显差异：

阈值₀反映信号与噪声STFT值相差程度，设为5，τ₀反映噪声门限值，其值设定与噪声大小有关；将满足上述条件的时频点集用Φ^*表示，再次使用k-means聚类法将Φ^*中的点分为K类，每一类用集合Φ_i,i＝1,…,K表示；至此，各个信源的单信源时频点被选取出来。

步骤5)所述的依据Φ_i,i＝1,…,K,建立单个信源的时频接收数据矩阵按如下步骤进行：

步骤5a)、建立一个三维矩阵D_y{m}(t,f)，存放2M个阵元的时频分布：

定义索引矢量其各元素代表点集Φ_i中每个单信源时频点在Φ^*上对应的坐标，坐标定义为：例如有一个3×3矩阵，坐标5代表第二列的第二个元素；坐标7代表第三列的第一个元素。N_i为Φ_i中时频点数，时频接收数据矩阵可用如下计算式获得：

Y_tf(m,:)＝D_y{m}(q_i),m＝1,…,2M

步骤5b)、分离两个方向的时频数据：

X_tf＝Y_tf(1:M,:),Z_tf＝Y_tf(M+1:2M,:)；

x和z方向上的接收信号都来自于相同的信号，因此x和z方向上的信号具有相同的时频分布，对总体时频接收数据W_Y(t,f)进行一次单信源分离，即可获得x,z方向上时频坐标相同的匹配时频点。

与现有算法相比，本发明具有下述优点：

1)能够解决实际应用中时频谱混叠和空间邻近信源的二维DOA估计问题。本发明将时域信源转换到时频域中，放大信源特征，利用合理的阈值设定和k-means聚类法，选取信源在二维方向上具有相同坐标的单信源时频点进行计算，消除了多信号时频点和噪声时频点带来的影响。

2)在低信噪比情况下仍具有良好的估计效果。本发明一方面利用噪声在时频域上均匀分布，而信号在时频处有聚集性的特点集中信号强度，另一方面过滤能量低的噪声时频点，提高对噪声的鲁棒性。

3)能够在欠定情况下保持高估计精度。本发明分别对单个信源进行二维DOA估计，则可以放宽各方向阵元数M必须大于信源数K的要求，仅需少数阵元就可以获得较高精度的估计结果。

附图说明

图1是本发明的实现流程图；

图2是本发明的阵列及参数示意图；

图3是单脉冲(AM),线性调频(LFM)，正弦调频(SFM)以及频移键控(FSK)信号在进行单信源时频点选取前后的时频分布图；

图4是本发明与现有两种二维DOA估计算法在不同信噪比条件下均方根误差对比图；

图5是本发明与现有两种二维DOA估计算法在不同仰角角度间隔下估计成功率对比图；

图6是本发明与现有两种二维DOA估计算法在非欠定(a)和欠定(b)情况下估计效果图，其中，图6(a)为M＝10,K＝4时二维DOA估计图，图6(b)为M＝3,K＝4时二维DOA估计图。

具体实施方式

以下参照附图，对本发明的技术方案和效果作进一步的详细说明。

参照图1，本发明的实现步骤如下：

步骤1：利用两个均匀线阵形成xz平面的L型阵列。

在x和z方向上，分别以相同间隔d放置阵元，各方向上有M个阵元，并且规定原点阵元处于z方向上。两均匀线阵正交摆放形成xz平面上的L型阵列。假设有K个频率混叠远场信号(N代表快拍采样个数)以方位角θ_k和仰角φ_k入射到此阵列上，且信号在传播过程中加入均值为0的复高斯白噪声；

步骤2：对x,z方向的均匀线阵接收数据，分别进行时域采样，获得输出数据X(n)，Z(n),n＝1,…,N,构造总体接收数据Y(n)：

其中Y(n)维度为2M×N,X(n)＝Y(1:M,:),Z(n)＝Y(M+1:2M,:)；

步骤3：采用STFT将时域总体接收数据变换到时频域。

总体接收数据Y(n)的STFT值可以用X(n),Z(n)的STFT值表示：

其中W_Y(t,f),W_X(t,f),W_Z(t,f),W_S(t,f)分别表示Y(n),X(n),Z(n),S(n)的短时傅里叶值，A_x,A_z代表x和z阵上形成接收数据的导向矢量阵。至此，接收信号的信息可以通过时频点的分布来表示；

步骤4：选取各信源在二维方向上的单信源时频点。

4a)选择高能量信号项时频点：

式中W_Y(t,v)表示总体输出数据在时间点t时，能获得的最大STFT值，且该值出现于频率v处。δ₀定义为高能量时频点选定阈值，取值范围为0.3～0.5，在此范围内取值大小仅影响选取的时频点个数，对最终估计效果无影响。满足条件的点集为Φ，计算Φ中每个点的空间矢量：

分别代表矩阵Y(t)各行，即各阵元接收数据的STFT值。采用k-means聚类算法将计算出的空间矢量分为K₀类，K₀＞K，得到K₀个方向矢量：

其中γ_k是第k类时频点集合，N_γk是集合γ_k中时频点的个数；方向矢量构成矩阵

4b)计算Φ中每个时频点上最可能出现的两个信源所对应的导向矢量a_n1,a_n2：

其中代表指向的噪声子空间的正交投影矩阵,其元素是矩阵A₀中的随机两列，I为单位矩阵，根据导向矢量a_n1,a_n2可以获得对应的两个STFT值：

其中A₂＝[a_n1,a_n2]，代表Moore-Penrose伪逆运算符；

4c)选取单信源时频点并分类，步骤如下：

单信源时频点上仅有信号项和噪声项，因此W₁,W₂的数量级有明显差异：

阈值λ₀反映信号与噪声STFT值相差程度，设为5，τ₀反映噪声门限值，其值设定与噪声大小有关；将满足上述条件的时频点集用Φ^*表示，再次使用k-means聚类法将Φ^*中的点分为K类，每一类用集合Φ_i,i＝1,…,K表示；至此，各个信源的单信源时频点被选取出来。

步骤5、根据上述选择结果，建立单信源的时频接收数据矩阵。

5a)建立一个三维矩阵D_y{m}(t,f)，存放2M个阵元的时频分布：

定义索引矢量其各元素代表点集Φ_i中每个单信源时频点在Φ^*上对应的坐标，坐标定义为：例如有一个3×3矩阵，坐标5代表第二列的第二个元素；坐标7代表第三列的第一个元素。N_i为Φ_i中时频点数，时频接收数据矩阵可用如下计算式获得：

Y_tf(m,:)=D_y{m}(q_i),m＝1,…,2M

5b)分离两个方向的时频数据：

X_tf＝Y_tf(1:M,:),Z_tf＝Y_tf(M+1:2M,:)；

x和z方向上的接收信号都来自于相同的信号，因此x和z方向上的信号具有相同的时频分布，对总体时频接收数据W_Y(t,f)进行一次单信源分离，即可获得x,z方向上时频坐标相同的匹配时频点。

步骤6：求解时频互相关矩阵：

步骤7：根据线性阵列的共轭对称特性，以及均匀线阵的旋转不变性，计算扩大阵列孔径后的新接收矩阵。

7a)划分时频互相关子阵：

Y₁＝R_xz(:,1:M-1),Y₂＝R_xz(:,2:M)，

由于接收数据的导向矢量阵是范德蒙特矩阵，且线性阵列具有共轭对称的特性，因此建立两个新的互相关子阵：

其中J_M为副对角线全1的M×M维置换矩阵。

7b)构造类似孔径扩大后的新接收矩阵：

新接收矩阵维度扩大为4M×(M-1)。

步骤8、根据传播算子原理，构成Y的新导向矢量阵可以被传播算子P分为两部分，传播算子的计算公式如下：

其中R_yy1,R_yy2为自相关矩阵R_yy=E{YY^H}前K行和后4M-K行，为伪逆运算符。为了后面利用旋转不变性建立子阵关系，这里定义一个扩大的传播算子:

P_e=[I_K,P]^H

步骤9、定义两个方位角选择矩阵和两个仰角选择矩阵，用于建立孔径扩大后子阵间旋转不变关系：

其中0_(M-1)×1,0_M×M代表(M-1)×1维和M×M的全0矩阵。根据旋转不变等式，可以推导出θ_k，φ_k的计算式：

其中Φ_x为对角阵，对角元素由矩阵奇异值分解后的特征值按降序排列构成；Γ由矩阵奇异值分解后的特征矢量阵构成。

若信源数目K＝1，顺序执行步骤1)至9)；若K＞1,执行至第4)步后步骤5)-9)重复K次，共获得K对计算结果。

本发明的效果可通过以下仿真说明：

仿真条件：可能采用的信号为AM,LFM，SFM以及FSK，各信号时频混叠情况如图3上层所示，经过单信源选取后时频分布如图3下层所示。采样点数为512，观测空域角度范围为[-90°,90°]。均方根误差的计算式为：

其中J表示实验次数，φ_k,θ_k分别表示二维角度估计值与实际值。估计成功率表示J次独立实验中检测成功的次数，其中一次实验检测成功的定义为该次所有的DOA估计值均在真实值±1°范围内。

仿真1：信噪比对算法的影响

假设二维各方向M＝6，K＝3，信源为LFM,SFM和FSK信号，入射角度θ＝[98.49°,101.42°,103.08°]，φ＝[91.72°,65.84°,76.16°]。信噪比由-5-15dB变化，本发明与JSVD，Dong三种算法在实验次数J为1000时，均方根误差结果如图4。

图4中横坐标表示信噪比值，纵坐标表示均方根误差。

从图4中可以看出本发明在各信噪比下，甚至是低信噪比时，均方根误差都明显小于其他两种算法。

仿真2：不同仰角角度间隔下的估计成功率

假设二维各方向M＝5，信噪比为10dB。K＝2,分别是SFM,FSK，入射角度为θ＝[55.49°,65.42°]，φ＝[140.72°,140.72°+Δφ]，Δφ由1°-20°变化。本发明与JSVD，Dong三种算法在实验次数J为100时，估计成功率结果如图5。

图5中横坐标表示仰角角度间隔，纵坐标表示估计成功率。

从图5中可以看出，本发明与其他两种算法相比，大幅提高了估计成功率，即便两个信源邻近时，依然能够分辨二维DOA。

仿真3：阵元数对算法的影响

假设信噪比为10dB,K＝4,分别为AM,LFM,SFM和FSK，入射角度为θ＝[30°,35°,40°,41°]，φ＝[104°,70°,75°,80°]。本发明与其他两种算法在M＝10,M＝3(欠定情况，即阵元数小于信源数情况)时，仿真结果如图6(a),6(b)所示。

图6(a)(b)中横坐标表示方位角，纵坐标表示仰角，圆圈代表信源的真实DOA。

从6(a)中可以看出，阵元数大于信源数时，本发明比其他两种算法估计精度更高；

从6(b)中可以看出，本发明同样适用于欠定情况，且保持较高的估计精度，其他两种算法已失去估计功能。