CN113253192B - 一种用于非圆信号的互质线阵级联doa估计方法 - Google Patents

Abstract

本发明公开了一种用于非圆信号的互质线阵级联DOA估计方法，具体为：设置互质线阵天线阵列，对接收信号采样；利用非圆信号的非圆特性扩展两个子阵的输出数据；分别计算两个子阵扩展后接收数据的协方差矩阵，对两个协方差矩阵特征值分解获得各自的信号子空间和噪声子空间；利用子阵列信号子空间的旋转不变性求得所有模糊角度信息；消除步骤四DOA估计结果的模糊值，获得DOA初估计值；构造降维谱峰搜索函数，在初估计值附近进行DOA精估计。本发明充分利用了非圆信号的非圆特性，进一步扩展了互质阵列的有效孔径，从而提高了DOA估计精度，同时使用级联的方法避免了全局谱峰搜索，又利用降维的方法降低了谱峰搜索的维度，优化DOA估计算法。

Description

一种用于非圆信号的互质线阵级联DOA估计方法

技术领域

本发明属于阵列信号处理，尤其涉及一种用于非圆信号的互质线阵级联DOA估计方法。

背景技术

阵列信号处理具有波束控制灵活、信号增益高、空间分辨率高、抗干扰能力强等优点，因而在近三十年获得了快速发展，在雷达、通信和电子战等领域都有广泛的应用。波达方向估计(Direction of Arrival，DOA)是阵列信号处理的一个主要研究方向。基于旋转不变的信号参数估计(Estimating Signal Parameters via Rotational InvarianceTechniques，ESPRIT)算法和多重信号分类(Multiple Signal Classification，MUSIC)算法等传统DOA估计算法，在阵列规模较小时估计精度较低，阵列规模较大时计算复杂度很高；传统算法直接用于互质阵，会因为阵元间距大于半波长而失效。

发明内容

发明目的：本发明的目的是提供了一种用于非圆信号的互质线阵级联DOA估计方法，在不改变阵元数目情况下扩展阵列孔径，提高DOA估计精度，降低互耦影响，具有更好的角度估计性能。

技术方案：本发明提供了一种用于非圆信号的互质线阵级联DOA估计方法，具体包括如下步骤：

步骤1：设置互质线阵天线阵列，对接收信号采样；

步骤2：利用非圆信号的非圆特性扩展两个子阵的输出数据；

步骤3：分别计算两个子阵扩展后接收数据的协方差矩阵，对两个协方差矩阵特征值分解获得各自的信号子空间和噪声子空间；

步骤4：利用子阵列信号子空间的旋转不变性求得所有模糊角度信息；

步骤5：消除步骤四DOA估计结果的模糊值，获得DOA初估计值；

步骤6：构造降维谱峰搜索函数，在DOA初估计值附近进行DOA精估计。

进一步地，步骤1具体包括：

互质线阵天线阵列包括两个子阵，两个子阵的最左端重合；子阵1是阵元数为M₁的均匀线阵，阵元间距为M₂λ/2，子阵2是阵元数为M₂的均匀线阵，阵元间距为M₁λ/2，互质线阵的天线阵元总数为T＝M₁+M₂-1，阵元所在位置的集合表达式为：

P＝{M₂m₁d₀|0≤m₁≤(M₁-1)}∪{M₁m₂d₀|0≤m₂≤(M₂-1)} (1)

其中，d₀表示半波长，互质阵列的孔径为max(M₁(M₂-1)d₀，M₂(M₁-1)d₀)。

进一步地，步骤2具体包括：

子阵i在采样时刻t的接收信号数据模型表达式：

x_i(t)＝A_is(t)+n_i(t) (2)

其中，子阵i＝1，2；x_i(t)表示子阵i在t时刻的接收数据，t＝1，…，L，L为总采样快拍数，A_i表示子阵i的方向矩阵，n_i(t)表示阵列接收的方差是的零均值加性高斯白噪声，各阵元接收噪声之间互相独立，噪声与信号也互相独立，s(t)表示远场窄带不相干的非圆信源矢量，

s(t)＝Ψs₀(t) (3)

其中，s₀(t)是一个圆信号矢量，Ψ是一个大小为K×K的对角矩阵，K表示信源数，Ψ的第k个对角元素是 表示第k个信号的非圆相位，k＝1，2，...，K；

子阵i的方向矩阵A_i＝[a_i(θ₁)，a_i(θ₂)，...，a_i(θ_K)]，其中，a_i(θ_k)是第k个信源入射子阵i的导向矢量，a_i(θ_k)的表达式为：

其中，θ_k表示第k个信源发出的信号与接收阵列法线的夹角，表示子阵i相邻阵元之间的间距，其中/>且/>

扩展两个子阵的输出数据，将子阵i的接收数据x_i(t)扩展为

其中，J_i是M_i×M_i维的反单位矩阵，其副对角线的元素全为1，其余元素都是0，B_i是子阵i的扩展方向矩阵，

其中B_i1＝A_iΨ，

进一步地，步骤3具体包括：

子阵i扩展接收数据y_i(t)的协方差矩阵为

其中，表示维度为2M_i×2M_i的单位矩阵，其主对角线元素是1，其余元素都是0，是源协方差矩阵，各信源之间独立，R_s是一个对角矩阵；

通过L次快拍估计子阵i扩展接收信号的协方差矩阵表达式为：

对估计的协方差矩阵进行特征值分解，得到：

其中，U_si是的信号子空间，Λ_s是信号子空间对应的特征矢量，U_ni是/>的噪声子空间，Λ_n是噪声子空间对应的特征矢量，信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间。

进一步地，步骤4具体包括：

信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间，存在可逆矩阵T，使得：

U_si＝B_iT (10)

定义矩阵B_ai为矩阵B_i删除第M_i行元素和第2M_i行元素后的矩阵，矩阵B_bi为矩阵B_i删除第1行元素和第M_i+1行元素后的矩阵，扩展方向矩阵B_i满足旋转不变性，

B_aiΦ_i＝B_bi (11)

其中，是包含DOA信息的对角矩阵；

定义矩阵U_sai为矩阵U_si删除第M_i行元素和第2M_i行元素后的矩阵，矩阵U_sbi为矩阵U_si删除第1行元素和第M_i+1行元素后的矩阵，根据U_sai＝B_aiT，U_sbi＝B_biT，得到：

U_saiT^-1Φ_iT＝U_sbi (12)

令Γ_i＝T^-1Φ_iT，Γ_i的特征值即为Φ_i的对角元素，根据噪声子空间U_si得到Γ_i的表达式为：

对Γ_i进行特征值分解，如果分解得到第k个特征值为λ_k，则根据子阵i接收数据估计的第k处信源DOA表达式为：

其中，且/>

进一步地，步骤5具体包括：

和/>在真实值处不严格相等，第k处信源的初估计值/>由/>和/>的实数解中两个最相近值的平均得到：

进一步地，步骤6具体包括：

子阵i满足其中，b_i表示子阵i的扩展导向矢量，据此构造NC-MUSIC谱函数，寻找谱函数峰值位置求得DOA精估计，求解过程表述为：

其中，

构造降维的NC-MUSIC谱函数降低计算复杂度，表达式为：

其中，只与θ有关，/>只与/>有关；将/>中的/>替换，/>表达式为：

其中，

采用拉格朗日乘子法，添加约束条件以消除平凡解，其中，e＝[1，0]^T，得到子阵i降维的NC-MUSIC谱函数的表达式为：

其中，f_i(θ)取得极大值时的θ即为子阵i的DOA精估计结果，分别在每个初估计值附近计算子阵1和子阵2的降维NC-MUSIC谱函数；在第k个初估计值附近计算f₁(θ)和f₂(θ)，峰值位置分别是/>和/>第k个信源的DOA精估计的表达式为：

有益效果：与现有技术相比，本发明具有如下显著的优点：

(1)本发明方法采用互质阵，降低天线之间互耦影响；

(2)本发明方法利用了非圆信号的非圆特性，提高DOA估计精度；

(3)本发明方法只需要一维局部谱峰搜索，降低的运算复杂度。

附图说明

图1是本发明的算法流程图；

图2是本发明中互质线阵示意图；

图3是不同阵元个数下本发明方法与ESPRIT、RD-MUSIC算法复杂度对比图；

图4是不同快拍数下本发明方法与ESPRIT、RD-MUSIC算法复杂度对比图；

图5是不同信噪比下本发明方法与均匀线阵中MUSIC算法性能对比图；

图6是不同信噪比下本发明方法与ULA-MUSIC、CLA-MUSIC和CLA-NC-MUSIC算法性能对比图；

图7是不同快拍数下本发明方法与NC-ESPRIT、NC-MUSIC、NCRIPM算法性能对比图；

图8是不同信噪比下本发明方法与NC-ESPRIT、NC-MUSIC、NCRIPM算法性能对比图。

具体实施方式

符号表示：(·)^T表示矩阵转置运算，(·)^H表示矩阵共轭转置运算，大写字母X表示矩阵，小写字母x(·)表示矢量，e表示自然常数，j表示虚数符号，^*表示取复数共轭运算，angle(·)表示取复数的相角。

如图1所示，本实施例提供了一种用于非圆信号的互质线阵级联DOA估计方法，具体为：

步骤一：设置互质线阵天线阵列，对接收信号采样；

如图2所示，实施例中的互质天线阵列可以分成的两个子阵，配置天线阵列使得两个子阵的最左端重合。子阵1是阵元数为M₁的均匀线阵，阵元间距为M₂λ/2，子阵2是阵元数为M₂的均匀线阵，阵元间距为M₁λ/2，所述互质线阵的天线阵元总数为T＝M₁+M₂-1。那么阵元所在位置的集合可以表示为：

P＝{M₂m₁d₀|0≤m₁≤(M₁-1)}∪{M₁m₂d₀|0≤m₂≤(M₂-1)} (1)

其中，d₀表示半波长，由此可知互质阵列的孔径为max(M₁(M₂-1)d₀，M₂(M₁-1)d₀)，相较于相同阵元数时均匀线阵(M₁+M₂)d₀可知，互质线阵阵列孔径得到了显著提升。

步骤二：利用非圆信号的非圆特性扩展两个子阵的输出数据；

子阵i(i＝1，2)在采样时刻t的接收信号数据模型为：

x_i(t)＝A_is(t)+n_i(t) (2)

其中，x_i(t)表示子阵i在t时刻的接收数据，其中t＝1，…，L，L为总采样快拍数，A_i表示子阵i的方向矩阵，n_i(t)表示阵列接收的方差是的零均值加性高斯白噪声，各阵元接收噪声之间互相独立，噪声与信号也互相独立，s(t)表示远场窄带不相干的非圆信源矢量，可以分解为：

s(t)＝Ψs₀(t) (3)

其中s₀(t)是一个圆信号矢量，Ψ是一个大小为K×K的对角矩阵，其中K表示信源数，Ψ的第k个对角元素是其中/>表示第k个信号的非圆相位，k＝1，2，...，K。

子阵i的方向矩阵A_i＝[a_i(θ₁)，a_i(θ₂)，...，a_i(θ_K)]，其中a_i(θ_k)是第k个信源入射子阵i的导向矢量，a_i(θ_k)的具体形式为：

其中，θ_k表示第k个信源发出的信号与接收阵列法线的夹角，表示子阵i相邻阵元之间的间距，其中/>且/>

扩展两个子阵的输出数据，将子阵i的接收数据x_i(t)扩展为：

其中J_i是M_i×M_i维的反单位矩阵，其副对角线的元素全为1，其余元素都是0，B_i是子阵i的扩展方向矩阵，

其中B_i1＝A_iΨ，

步骤三：分别计算两个子阵扩展后接收数据的协方差矩阵，对两个协方差矩阵特征值分解获得各自的信号子空间和噪声子空间；

子阵i扩展接收数据y_i(t)的协方差矩阵为：

其中表示维度为2M_i×2M_i的单位矩阵，其主对角线元素是1，其余元素都是0，是源协方差矩阵，因为各信源之间独立，R_s是一个对角矩阵。

实际应用中信号快拍数是有限的，可以通过L次快拍估计子阵i扩展接收信号的协方差矩阵。

对估计的协方差矩阵进行特征值分解，

其中U_si是的信号子空间，Λ_s是信号子空间对应的特征矢量，U_ni是/>的噪声子空间，Λ_n是噪声子空间对应的特征矢量。可以发现信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间。

步骤四：利用子阵列信号子空间的旋转不变性求得所有模糊角度信息无需谱峰搜索就可以获得所有模糊的DOA值，具有较低的复杂度；

信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间，则存在可逆矩阵T，使得：

U_si＝B_iT (10)

定义矩阵B_ai为矩阵B_i删除第M_i行元素和第2M_i行元素后的矩阵，矩阵B_bi为矩阵B_i删除第1行元素和第M_i+1行元素后的矩阵。扩展方向矩阵B_i满足旋转不变性，

B_aiΦ_i＝B_bi (11)

其中是包含DOA信息的对角矩阵。定义矩阵U_sai为矩阵U_si删除第M_i行元素和第2M_i行元素后的矩阵，矩阵U_sbi为矩阵U_si删除第1行元素和第M_i+1行元素后的矩阵。显然U_sai＝B_aiT，U_sbi＝B_biT可以推出：

U_saiT^-1Φ_iT＝U_sbi (12)

令Γ_i＝T^-1Φ_iT，显然Γ_i的特征值即为Φ_i的对角元素。Γ_i可以由噪声子空间U_si求得：

对Γ_i进行特征值分解，如果分解得到第k个特征值为λ_k，则根据子阵i接收数据估计的第k处信源DOA可能为：

其中，且/>

步骤五：消除步骤四DOA估计结果的模糊值，获得DOA初估计值；

步骤四得到的共有/>个取值，但其中/>个解是无意义的复数解，首先去掉这些复数解开。/>个实数解中只有1个是真实角度，其余/>个实数解都是模糊角度。由于M₁和M₂是互质的，可以保证无噪声干扰时/>和/>的实数解只在真实值处重合，可以利用这一点消除模糊值。实际情况中接收信号总是含有噪声的，因此/>和/>在真实值处不会严格相等，第k处信源的初估计值/>可以由/>和/>的实数解中两个最相近值的平均得到：

步骤六：构造降维谱峰搜索函数，在初估计值附近进行DOA精估计。

首先利用扩展导向矢量和噪声子空间的正交性构造降维谱峰搜索函数，然后在初估计值附近计算谱峰搜索函数，进行谱峰搜索获得。方法扩展了互质阵列的有效孔径，具有很高的DOA估计精度，同时只需进行一维的局部谱峰搜索，算法的复杂度也比较低。

由于信号与噪声与不相干的，子阵的扩展导向矢量和其噪声子空间都是正交的。对于子阵i满足其中b_i表示子阵i的扩展导向矢量，可以据此构造NC-MUSIC谱函数，寻找谱函数峰值位置可以求得DOA精估计。求解过程可以表述为：

其中，求解过程是一个二维搜索过程，计算复杂度较高。DOA估计问题一般不关心信号非圆角度/>因此可以考虑构造降维的NC-MUSIC谱函数降低计算复杂度。/>可以写作：

其中只与θ有关，/>只与/>有关。将/>中的/>替换，因此/>可以改写为：

其中采用拉格朗日乘子法解决上述问题，添加约束条件/>以消除平凡解，其中e＝[1，0]^T，最终得到子阵i降维的NC-MUSIC谱函数为：

f_i(θ)取得极大值时的θ即为子阵i的DOA精估计结果。分别在每个初估计值附近计算子阵1和子阵2的降维NC-MUSIC谱函数。若在第k个初估计值附近计算f₁(θ)和f₂(θ)，求得峰值位置分别是/>和/>则第k个信源的DOA精估计为：

图3、图4是本发明方法与ESPRIT、RD-MUSIC算法复杂度对比图。其中图3是不同阵元个数下各算法复杂度对比图，图4是不同快拍数下各算法复杂度对比图，从两幅图可以本发明方法的复杂度低于RD-MUSIC算法但是高于NC-ESPRIT算法。其中NC-ESPRIT算法的复杂度是O(8(M-1)K²+13K³+8M³+4M²L)，RD-MUSIC算法复杂度是O(4M²L+8M³+(8M²-4MK)n)，而本发明方法的复杂度是O(8(M-1)K²+13K³+8M²+4M²L(8M²-4MK)n_l)，其中，M表示互质线阵阵元数，n表示全局搜索次数，n_i表示局部搜索次数，L表示快拍数。

图5是不同信噪比下本发明方法与均匀线阵中MUSIC算法性能对比图，给出了理论下界CRB。从图中可以看出本发明方法性能不但优于均匀线阵中的MUSIC算法，而且优于均匀线阵的理论下界CRB。其中快拍数为100，均匀线阵与互质线阵的阵元数相同，信源入射角度分别是为10°，30°。可以看到由于互质阵阵列的孔径更大，互耦影响较小，具有更好的角度估计性能。

图6是不同信噪比下本发明方法与ULA-MUSIC、CLA-MUSIC和CLA-NC-MUSIC算法性能对比图。从图中可以看出本发明方法的角度估计性能优于均匀线阵中MUSIC算法和互质线阵中MUSIC算法，与互质阵中NC-MUSIC算法非常接近。这是因为本发明方法和互质线阵中MUSIC算法利用了信号的非圆特性，并且使用互质线阵减弱了互耦的影响。其中快拍数为200，均匀线阵与互质线阵的阵元数相同，信源入射角度分别是为10°，30°。

图7是不同快拍数下本发明方法与NC-ESPRIT、NC-MUSIC、NCRIPM算法性能对比图，给出了理论下界CRB。从图中可以看出本发明方法性能优于NC-RIPM算法和NC-ESPRIT算法，与NC-MUSIC算法非常接近。其中子阵1阵元数M₁＝6，子阵2阵元数M₂＝5，信噪比为10dB，信源入射角度分别是10°，30°。

图8是不同信噪比下本发明方法与NC-ESPRIT、NC-MUSIC、NCRIPM算法性能对比图，给出了理论下界CRB。从图中本发明方法的性能优于NC-ESPRIT算法和NC-RIPM算法，NC-MUSIC算法非常接近，且低信噪比时优于NC-MUSIC算法。其中子阵1阵元数M₁＝7，子阵2阵元数M₂＝5，快拍数为100，信源入射角度分别是10°，30°。

Claims (6)

1.一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，具体包括如下步骤：

步骤1：设置互质线阵天线阵列，对接收信号采样；

步骤2：利用非圆信号的非圆特性扩展两个子阵的输出数据；

步骤3：分别计算所述两个子阵扩展后接收数据的协方差矩阵，对两个协方差矩阵特征值分解获得各自的信号子空间和噪声子空间；

步骤4：利用子阵列信号子空间的旋转不变性求得所有模糊角度信息；

步骤5：消除步骤四DOA估计结果的模糊值，获得DOA初估计值；

步骤6：构造降维谱峰搜索函数，在所述DOA初估计值附近进行DOA精估计；

其中，所述步骤2具体包括：

子阵i在采样时刻t的接收信号数据模型表达式：

x_i(t)＝A_is(t)+n_i(t) (2)

其中，子阵i＝1,2；x_i(t)表示子阵i在t时刻的接收数据，t＝1,…,L，L为总采样快拍数，A_i表示子阵i的方向矩阵，n_i(t)表示阵列接收的方差是的零均值加性高斯白噪声，各阵元接收噪声之间互相独立，噪声与信号也互相独立，s(t)表示远场窄带不相干的非圆信源矢量，

s(t)＝Ψs₀(t) (3)

其中，s₀(t)是一个圆信号矢量，Ψ是一个大小为K×K的对角矩阵，K表示信源数，Ψ的第k个对角元素是表示第k个信号的非圆相位，k＝1,2,...,K；

子阵i的方向矩阵A_i＝[a_i(θ₁),a_i(θ₂),...,a_i(θ_K)]，其中，a_i(θ_k)是第k个信源入射子阵i的导向矢量，a_i(θ_k)的表达式为：

其中，θ_k表示第k个信源发出的信号与接收阵列法线的夹角，λ表示信号波长，表示子阵/>的阵元数，/>表示子阵i相邻阵元之间的间距，其中/>且/>

扩展两个子阵的输出数据，将子阵i的接收数据x_i(t)扩展为

其中，J_i是M_i×M_i维的反单位矩阵，其副对角线的元素全为1，其余元素都是0，^*表示取复数共轭运算，B_i是子阵i的扩展方向矩阵，

其中B_i1＝A_iΨ，

2.根据权利要求1所述的一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，所述步骤1具体包括：

所述互质线阵天线阵列包括两个子阵，所述两个子阵的最左端重合；子阵1是阵元数为M₁的均匀线阵，阵元间距为M₂λ/2，子阵2是阵元数为M₂的均匀线阵，阵元间距为M₁λ/2，所述互质线阵的天线阵元总数为T＝M₁+M₂-1，阵元所在位置的集合表达式为：

P＝{M₂m₁d₀|0≤m₁≤(M₁-1)}∪{M₁m₂d₀|0≤m₂≤(M₂-1)} (1)

其中，d₀表示半波长，互质阵列的孔径为max(M₁(M₂-1)d₀,M₂(M₁-1)d₀)。

3.根据权利要求1所述的一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，所述步骤3具体包括：

子阵i扩展接收数据y_i(t)的协方差矩阵为

其中，表示维度为2M_i×2M_i的单位矩阵，其主对角线元素是1，其余元素都是0，是源协方差矩阵，各信源之间独立，R_S是一个对角矩阵；

通过L次快拍估计子阵i扩展接收信号的协方差矩阵表达式为：

对估计的协方差矩阵进行特征值分解，得到：

其中，U_si是的信号子空间，Λ_s是信号子空间对应的特征矢量，U_ni是/>的噪声子空间，Λ_n是噪声子空间对应的特征矢量，信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间。

4.根据权利要求1或3所述的一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，所述步骤4具体包括：

信号子空间U_si和扩展方向矩阵B_i张成相同线性子空间，存在可逆矩阵T,使得：

U_si＝B_iT (10)

定义矩阵B_ai为矩阵B_i删除第M_i行元素和第2M_i行元素后的矩阵，矩阵B_bi为矩阵B_i删除第1行元素和第M_i+1行元素后的矩阵，扩展方向矩阵B_i满足旋转不变性，

B_aiΦ_i＝B_bi (11)

其中，是包含DOA信息的对角矩阵；

定义矩阵U_sai为矩阵U_si删除第M_i行元素和第2M_i行元素后的矩阵，矩阵U_sbi为矩阵U_si删除第1行元素和第M_i+1行元素后的矩阵，根据U_sai＝B_aiT，U_sbi＝B_biT，得到：

U_saiT^-1Φ_iT＝U_sbi (12)

令Γ_i＝T^-1Φ_iT，Γ_i的特征值即为Φ_i的对角元素，根据噪声子空间U_si得到Γ_i的表达式为：

对Γ_i进行特征值分解，如果分解得到第k个特征值为λ_k，则根据子阵i接收数据估计的第k处信源DOA表达式为：

其中，且/>

5.根据权利要求4所述的一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，所述步骤5具体包括：

和/>在真实值处不严格相等，第k处信源的初估计值/>由/>和/>的实数解中两个最相近值的平均得到：

6.根据权利要求5所述的一种用于非圆信号的互质线阵级联DOA估计方法，其特征在于，所述步骤6具体包括：

子阵i满足其中，b_i表示子阵i的扩展导向矢量，据此构造NC-MUSIC谱函数，寻找谱函数峰值位置求得DOA精估计，求解过程表述为：

其中，

构造降维的NC-MUSIC谱函数降低计算复杂度，表达式为：

其中，只与θ有关，/>只与/>有关；将/>中的替换，/>表达式为：

其中，

采用拉格朗日乘子法，添加约束条件以消除平凡解，其中，e＝[1,0]^T，得到子阵i降维的NC-MUSIC谱函数的表达式为：

其中，f_i(θ)取得极大值时的θ即为子阵i的DOA精估计结果，分别在每个初估计值附近计算子阵1和子阵2的降维NC-MUSIC谱函数；在第k个初估计值附近计算f₁(θ)和f₂(θ)，峰值位置分别是/>和/>第k个信源的DOA精估计的表达式为：