基于对称循环频率的谱相关子空间测向方法与流程

本发明涉及电子信息技术领域，具体来说，涉及基于对称循环频率的谱相关子空间测向方法。

背景技术：

阵列信号处理技术可用于对阵列接收的信号进行处理，从而测定多个信号的波达方向(简称测向)，在雷达、声纳、通信系统以及智能天线等领域有着广泛的应用。

在现有的各种测向方法中，干涉仪测向方法对单个信号具有良好的测向性能，但在同频多信号时方法失效，在应用上受到较大的限制。传统子空间类测向方法对同频多信号具有良好的测向性能，但是测向自由度受天线阵元数的限制，信号数越多需要的天线阵元数越大，无线电接收机体积大，在特定应用上受到限制；基于信号循环平稳特性的谱相关子空间测向方法因对信号具有选择特性，在扩展天线测向自由度的同时提升算法的抗干扰能力，从而获得广泛应用。但是现有谱相关子空间算法在信号二阶循环统计一致的情况下，算法失效。

针对相关技术中的问题，目前尚未提出有效的解决方案。

技术实现要素：

针对相关技术中的上述技术问题，本发明提出基于对称循环频率的谱相关子空间测向方法，包括以下步骤：

步骤1：初始化处理：初始化技术参数并进行存储；

其中，所述技术参数包括阵列的阵元个数m、阵元位置md、信号的传播速度c、信号载波频率f、信号采样频率fs、方位角θ、方位角θ的划分间隔δθ、采样快拍次数t0、信号循环频率a、信号对称循环频率-a、循环频率为a的信号个数k，自相关时延τ0；

步骤2：根据步骤1中的方位角θ、方位角θ的划分间隔δθ确定信号的方位角的离散取值集合及该离散取值集合对应的导向矢量集合；

步骤3：确定阵列接收信号的时域采样矢量：接收装置采样阵列所有阵元接收的入射信号，从而确定阵列接收信号的时域采样矢量：

x(t)＝[x1(t),x2(t),...,xm(t)]

其中：t是模数转换的时刻，即接收信号的时域采样时刻，t＝1,2,...,t0，m是阵元个数，t0是快拍次数；

步骤4：确定阵元接收信号样本数据的循环互相关矩阵，包括：

步骤4-1：根据步骤1中设定的信号循环频率α,计算步骤3中时域采样矢量中的各个矢量xi(t)，i＝1,2,...,m的循环自相关函数

步骤4-2：根据步骤1中设定的自相关时延τ0确定对应的循环自相关值和其中：i＝1,2,...,m；

将各个阵元对应的组成矩阵以及将各个阵元对应组成矩阵

步骤4-3：根据步骤4-2中的确定其对应的循环互相关矩阵x(α)，x(-α)；

步骤4-4：根据步骤4-3中的x(α)，x(-α)构造阵列接收样本数据循环互相关矩阵ψ(α)；

步骤5：确定导向矢量集合中的各个导向矢量与噪声子空间的伪谱，包括：

步骤5-1：分别对步骤4中ψ(α)进行特征值分解确定其噪声子空间；

步骤5-2：确定步骤2中导向矢量集合中的各个导向矢量a(α,θi)与步骤5-1中噪声子空间的伪谱p(α,θi)；

步骤6：确定信号波达方向：分别在步骤5各伪谱p(α,θi)中搜索最大值，每个伪谱中的最大值对应一个导向矢量，取该导向矢量对应的方位角即是测定的信号波达方向。

进一步的，步骤1中通过信号参数估计方法确定信号循环频率a和信号个数k，其中：自相关时延τ0的取值为τ0＝[12345678]。

进一步的，步骤2中：根据步骤1中的方位角θ、方位角θ的划分间隔δθ确定信号的方位角的离散取值集合及该离散取值集合对应的导向矢量集合具体包括：

步骤2-1，根据步骤1中方位角θ划分间隔δθ、将方位角θ均匀划分成nθ个离散取值为θi的集合；

步骤2-2，对应每个离散取值θi,α，确定阵列导向矢量a(α,θi)，其中：i＝1,2,...,nθ。

进一步的，步骤2中：根据每个离散取值θi，确定阵列导向矢量a(α,θi)为：

a(α,θi)＝[a1(α,θi)，a2(α,θi),lam(α,θi)]

其中，τm(i)＝dmsin(θi)/c是信号到达第m个阵元的时间差，c是信号的传播速度，i＝1,2,...,nθ，m＝1,2,...,m，m是阵元个数，α是第k个信号循环频率。

进一步的，步骤4中所述阵元接收信号样本数据的循环自相关函数和的计算方式为：

其中<·>表示求时间平均运算。

进一步的，步骤4中循环互相关矩阵x(α)和x(-α)的计算方式为：

其中：α为信号的循环频率，-α为与α对应的对称循环频率。

进一步的，步骤4中构造阵列接收样本数据循环互相关矩阵ψ(α)的计算方式为：

ψ(α)＝x(α)+i·x(-α)·i

其中：i为m×m交换矩阵。

进一步的，步骤5中确定噪声子空间的方式为：通过对阵列接收样本数据循环互相关矩阵ψ(α)进行特征值分解得到相应的特征值和特征向量，根据进行特征值分解，

其中：uk＝[uk1,uk2,...,ukm]，uk1,uk2,...,ukm为左特征向量，vk1,vk2,...,vkm为右特征向量，σ＝diag{λk1,λk2,...λkm}，λk1＞λk2＞...＞λkm为对应的特征值；

根据子空间理论，λk(k+1),λk(k+2),...,λkm对应的特征向量张成的空间即为噪声子空间ukn＝[uk(k+1),uk(k+2),...,ukm]。

进一步的，步骤5中：确定步骤2中导向矢量集合中的各个导向矢量a(α,θi)与步骤5-1中噪声子空间的伪谱p(α,θi)的计算方式为：

p(α,θi)＝20glg(||a(α,θi)||/||ukna(α,θi)||)

其中||·||表示求模。

本发明的有益效果：基于对称循环频率的谱相关子空间测向方法，依据接收信号的循环平稳特性建立循环自相关函数，利用对称循环频率的循环自相关函数之间的弱相关性，构造循环互相关矩阵，使测向的精度优于现有谱相关子空间测向方法，且具有更高的测向正确概率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是根据本发明基于对称循环频率的谱相关子空间测向方法工作原理流程框图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。

根据本发明的实施例，提供了基于对称循环频率的谱相关子空间测向方法，采用16阵元的均匀线阵，相邻阵元间距为d＝37.5米，选定阵元所在的直线为x轴，从左边数第一个阵元为原点建立参考坐标系。两个载波频率为40mhz，码速率为4mhz/s的bpsk调制信号从1.43°,5.74°方向入射到均匀线阵，接收装置采集所有阵元接收信号。

具体流程包括以下步骤：

步骤1：初始化处理：初始化技术参数并进行存储；

其中，所述技术参数包括阵列的阵元个数m、阵元位置md、信号的传播速度c、信号载波频率f、信号采样频率fs、方位角θ、方位角θ的划分间隔δθ、采样快拍次数t0、信号循环频率a、信号对称循环频率-a、循环频率为a的信号个数k，自相关时延τ0；具体的，通过信号参数估计方法确定信号循环频率a和信号个数k；

本实施例中：初始化阵列的阵元个数m＝16，阵元位置的直角坐标为md，其中d＝37.5,m＝0,1,...,15，信号的传播速度c＝3×10⁸m/s，信号载波频率f＝40mhz，信号采样频率fs＝320mhz，方位角θ∈[-90°,90°]，划分间隔δθ＝0.1°，划分个数nh＝180°/δθ+1＝1801，信号循环频率α＝4e6，其对称循环频率-α＝-4e6，信号个数k＝2，交换矩阵im×m,自相关时延τ0＝[12345678]，快拍次数t0＝2000。

步骤2：根据步骤1中的方位角θ、方位角θ的划分间隔δθ确定信号的方位角的离散取值集合及该离散取值集合对应的导向矢量集合；

本实施例中：根据每个离散取值θi，确定阵列导向矢量a(α,θi)为：

a(α,θi)＝[a1(α,θi)，a2(α,θi),lam(α,θi)]

其中，τm(i)＝dmsin(θi)/c是信号到达第m个阵元的时间差，c是信号的传播速度，i＝1,2,...,nθ，m＝1,2,...,m，m是阵元个数，α是第k个信号循环频率。

本实施例中：首先，根据步骤1中方位角θ划分间隔0.1度将方位角θ均匀划分成1801离散取值为θi＝(i-1)×0.1度的集合，i＝1,2,...,1801；

其次，对应每个离散取值θi,α，确定信号从波达方向方向θi入射时的阵列导向矢量a(α,θi)，阵列导向矢量a(α,θi)的16个元素分别通过下式确定：

其中：m＝0,1,...,15，i＝1,2,...,1801，d为阵元间距，c为信号的传播速度，α为信号循环频率；对应波达方向1.4°，其导向矢量的前8个元素为：

1.0000+0.0000i,0.9969-0.0785i,0.9877-0.1564i,0.9724-0.2334i

0.9511-0.3090i,0.9239-0.3827i,0.8910-0.4540i,0.8526-0.5225i

对应波达方向5.7°，其导向矢量的前8个元素为：

1.0000+0.0000i,0.9511-0.3090i,0.8090-0.5878i,0.5878-0.8090i

0.3090-0.9511i,0.0000-1.0000i,-0.3090-0.9511i,-0.5878-0.8090i

步骤3：确定阵列接收信号的时域采样矢量：接收装置采样阵列所有阵元接收的入射信号，从而确定阵列接收信号的时域采样矢量：

x(t)＝[x1(t),x2(t),...,xm(t)]

其中：t是模数转换的时刻，即接收信号的时域采样时刻，t＝1,2,...,t0，m是阵元个数，t0是快拍次数；

第一个阵元接收信号的时域采样前8个元素分别为：

-0.3878+0.2340i,-0.3972+0.6042i,0.2163-0.1836i,-0.2488+0.1669i,

-0.2183+0.0431i,-0.3023-0.5051i,0.1275+0.1100i,-0.3059+0.0095i

最后一个阵元接收信号的时域采样前8个元素分别为：

-1.6073-0.4407i,1.9497-0.1891i,0.9620-0.3590i,-1.3549+0.0334i,

-1.9182+0.2205i,-1.4811-0.0252i,-0.9528-0.1640i,0.9736+0.3686i

步骤4：确定阵元接收信号样本数据的循环互相关矩阵，包括：

步骤4-1：根据步骤1中设定的信号循环频率α,计算步骤3中时域采样矢量中的各个矢量xi(t)，i＝1,2,...,m的循环自相关函数

步骤4-2：根据步骤1中设定的自相关时延τ0确定对应的循环自相关值和其中：i＝1,2,...,m；

将各个阵元对应的组成矩阵以及将各个阵元对应组成矩阵

步骤4-3：根据步骤4-2中的确定其对应的循环互相关矩阵x(α)，x(-α)；

步骤4-4：根据步骤4-3中的x(α)，x(-α)构造阵列接收样本数据循环互相关矩阵ψ(α)；

本实施例中：步骤4中所述阵元接收信号样本数据的循环自相关函数和的计算方式为：

其中<·>表示求时间平均运算。

本实施例中：步骤4中循环互相关矩阵x(α)和x(-α)的计算方式为：

其中：α为信号的循环频率，-α为与α对应的对称循环频率。

本实施例中：步骤4中构造阵列接收样本数据循环互相关矩阵ψ(α)的计算方式为：

ψ(α)＝x(α)+i·x(-α)·i

其中：i为m×m交换矩阵。

具体的，步骤4中确定阵元接收信号样本数据的循环互相关矩阵：对应的循环互相关矩阵的第一行前8个元素为：

0.0750,0.0715-0.0124i,0.0653-0.0249i,0.0572-0.0332i,

0.0414-0.0389i,0.0343-0.0463i,0.0237-0.0450i,0.0092-0.0365i

对应的循环互相关矩阵的最后一行前8个元素为：

0.0283-0.0029i,0.0197-0.0062i,0.0086-0.0015i,0.0072-0.0011i

-0.0028+0.0092i,-0.0005+0.0130i,0.0024+0.0238i,0.0029+0.0350i

步骤5：确定导向矢量集合中的各个导向矢量与噪声子空间的伪谱，包括：

步骤5-1：分别对步骤4中ψ(α)进行特征值分解确定其噪声子空间；

步骤5-2：确定步骤2中导向矢量集合中的各个导向矢量a(α,θi)与步骤5-1中噪声子空间的伪谱p(α,θi)；

本实施例中：步骤5中确定噪声子空间的方式为：通过对阵列接收样本数据循环互相关矩阵ψ(α)进行特征值分解得到相应的特征值和特征向量，根据进行特征值分解，

其中：uk＝[uk1,uk2,...,ukm]，uk1,uk2,...,ukm为左特征向量，vk1,vk2,...,vkm为右特征向量，σ＝diag{λk1,λk2,...λkm}，λk1＞λk2＞...＞λkm为对应的特征值；

根据子空间理论，λk(k+1),λk(k+2),...,λkm对应的特征向量张成的空间即为噪声子空间ukn＝[uk(k+1),uk(k+2),...,ukm]。

本实施例中：步骤5中：确定步骤2中导向矢量集合中的各个导向矢量a(α,θi)与步骤5-1中噪声子空间的伪谱p(α,θi)的计算方式为：

p(α,θi)＝20glg(||a(α,θi)||/||ukna(α,θi)||)

其中||·||表示求模。

具体的：步骤5中确定导向矢量集合中的各个导向矢量与噪声子空间的伪谱：波达方向-90°,-89.9°,-89.8°,-89.7°对应的伪谱值为0.009685，0.009685，0.009682，0.009679，总共确定1801个伪谱值。

步骤6：确定信号波达方向：分别在步骤5各伪谱p(α,θi)中搜索最大值，每个伪谱中的最大值对应一个导向矢量，取该导向矢量对应的方位角即是测定的信号波达方向；

本实施例中：在步骤5确定的伪谱值p(α,θi)，i＝1,2,...,1801中搜索k＝2个最大峰值，最大峰值等于1.865，最大峰值对应的导向矢量所对应的波达方向为6.1°，即测定的信号的波达方向；

次最大峰值等于1.836，次最大峰值对应的导向矢量所对应的波达方向为1.5°，即测定的信号的波达方向。

本实施例中本发明方法经测定的信号波达方向的均方根误差检验，统计500次独立试验的测定结果，信噪比为-4db～10db，步进2db的情况下，本发明方法以及现有谱相关子空间方法测定的信号波达方向的均方根误差对比如下表1所示；

表1：波达方向估计方法的误差性能比较

本发明方法经测定的信号波达方向的正确概率检验，统计500次独立试验的测定结果，信噪比为-10db～4db，步进2db的情况下，当测向误差绝对值小于2°时为一次正确估计。本发明方法以及现有谱相关子空间方法测定的信号波达方向的正确概率对比如下表2所示；

表2：波达方向估计方法的正确概率性能比较

可见，本发明方法能够精确测定信号波达方向；相比于现有谱相关子空间方法，测定的信号波达方向与实际的信号波达方向之间的均方根误差更小，并且本发明方法在测向误差绝对值小于2°为正确估计的情况下，正确概率高于现有谱相关子空间方法。

由此可见，借助于本发明的上述技术方案，基于对称循环频率的谱相关子空间测向方法，依据接收信号的循环平稳特性建立循环自相关函数，利用对称循环频率的循环自相关函数之间的弱相关性，构造循环互相关矩阵，使测向的精度优于现有谱相关子空间测向方法，且具有更高的测向正确概率。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。