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

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

背景技术：

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

技术实现要素：

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

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

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

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

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

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

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

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

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

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

p＝{m2m1d0|0≤m1≤(m1-1)}∪{m1m2d0|0≤m2≤(m2-1)}(1)

其中，d0表示半波长，互质阵列的孔径为max(m1(m2-1)d0，m2(m1-1)d0)。

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

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

xi(t)＝ais(t)+ni(t)(2)

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

s(t)＝ψs0(t)(3)

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

子阵i的方向矩阵ai＝[ai(θ1)，ai(θ2)，...，ai(θk)]，其中，ai(θk)是第k个信源入射子阵i的导向矢量，ai(θk)的表达式为：

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

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

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

其中bi1＝aiψ，

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

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

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

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

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

其中，usi是的信号子空间，λs是信号子空间对应的特征矢量，uni是的噪声子空间，λn是噪声子空间对应的特征矢量，信号子空间usi和扩展方向矩阵bi张成相同线性子空间。

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

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

usi＝bit(10)

定义矩阵bai为矩阵bi删除第mi行元素和第2mi行元素后的矩阵，矩阵bbi为矩阵bi删除第1行元素和第mi+1行元素后的矩阵，扩展方向矩阵bi满足旋转不变性，

baiφi＝bbi(11)

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

定义矩阵usai为矩阵usi删除第mi行元素和第2mi行元素后的矩阵，矩阵usbi为矩阵usi删除第1行元素和第mi+1行元素后的矩阵，根据usai＝bait，usbi＝bbit，得到：

usait^-1φit＝usbi(12)

令γi＝t^-1φit，γi的特征值即为φi的对角元素，根据噪声子空间usi得到γi的表达式为：

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

其中，且

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

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

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

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

其中，

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

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

其中，

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

其中，fi(θ)取得极大值时的θ即为子阵i的doa精估计结果，分别在每个初估计值附近计算子阵1和子阵2的降维nc-music谱函数；在第k个初估计值附近计算f1(θ)和f2(θ)，峰值位置分别是和第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是阵元数为m1的均匀线阵，阵元间距为m2λ/2，子阵2是阵元数为m2的均匀线阵，阵元间距为m1λ/2，所述互质线阵的天线阵元总数为t＝m1+m2-1。那么阵元所在位置的集合可以表示为：

p＝{m2m1d0|0≤m1≤(m1-1)}∪{m1m2d0|0≤m2≤(m2-1)}(1)

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

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

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

xi(t)＝ais(t)+ni(t)(2)

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

s(t)＝ψs0(t)(3)

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

子阵i的方向矩阵ai＝[ai(θ1)，ai(θ2)，...，ai(θk)]，其中ai(θk)是第k个信源入射子阵i的导向矢量，ai(θk)的具体形式为：

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

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

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

其中bi1＝aiψ，

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

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

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

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

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

其中usi是的信号子空间，λs是信号子空间对应的特征矢量，uni是的噪声子空间，λn是噪声子空间对应的特征矢量。可以发现信号子空间usi和扩展方向矩阵bi张成相同线性子空间。

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

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

usi＝bit(10)

定义矩阵bai为矩阵bi删除第mi行元素和第2mi行元素后的矩阵，矩阵bbi为矩阵bi删除第1行元素和第mi+1行元素后的矩阵。扩展方向矩阵bi满足旋转不变性，

baiφi＝bbi(11)

其中是包含doa信息的对角矩阵。定义矩阵usai为矩阵usi删除第mi行元素和第2mi行元素后的矩阵，矩阵usbi为矩阵usi删除第1行元素和第mi+1行元素后的矩阵。显然usai＝bait，usbi＝bbit可以推出：

usait^-1φit＝usbi(12)

令γi＝t^-1φit，显然γi的特征值即为φi的对角元素。γi可以由噪声子空间usi求得：

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

其中，且

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

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

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

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

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

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

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

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

fi(θ)取得极大值时的θ即为子阵i的doa精估计结果。分别在每个初估计值附近计算子阵1和子阵2的降维nc-music谱函数。若在第k个初估计值附近计算f1(θ)和f2(θ)，求得峰值位置分别是和则第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)nl)，其中，m表示互质线阵阵元数，n表示全局搜索次数，ni表示局部搜索次数，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阵元数m1＝6，子阵2阵元数m2＝5，信噪比为10db，信源入射角度分别是10°，30°。

图8是不同信噪比下本发明方法与nc-esprit、nc-music、ncripm算法性能对比图，给出了理论下界crb。从图中本发明方法的性能优于nc-esprit算法和nc-ripm算法，nc-music算法非常接近，且低信噪比时优于nc-music算法。其中子阵1阵元数m1＝7，子阵2阵元数m2＝5，快拍数为100，信源入射角度分别是10°，30°。