CN110749856B - 一种基于零化去噪技术的互质阵欠定测向方法 - Google Patents

Abstract

本发明提供了一种基于零化去噪技术的互质阵欠定测向方法，具体包括以下步骤：互质阵布阵并接收阵列数据；数据预处理；构造均匀虚拟阵列数据向量及模型噪声协方差矩阵；迭代求解零化系数；方程求根并测向。该方法相比现有技术的优势在于：首先，在欠定求解过程中引入了对阵列噪声方差的估计，且用零化去噪的技术实现了对互质阵虚拟阵列中的孔洞进行了插值，充分利用了所有的阵元信息以及阵列自由度；其次，本发明对由有限快拍数引起的模型噪声进行了建模，可提高测向精度；最后，本发明提出的方法属于无网格算法，不存在网格效应，也无需进行空域网格扫描。

Description

一种基于零化去噪技术的互质阵欠定测向方法

技术领域

本发明属于阵列信号处理领域，特别涉及对雷达、通信等电磁波和声学机械波信号来波方向的测定，具体是一种适用于欠定条件下的互质阵测向方法。

背景技术

阵列测向是用传感器阵列来测定信号源的来波方向，广泛应用于雷达、通信、声呐等军用和民用领域。在过去的几十年里，人们提出了大量优秀的测向方法。但这些方法都存在一个共性问题，即只能对个数小于阵元数的信号进行测向。当信号源数大于或等于阵元数时，传统方法会失效，我们把这种条件下的测向称作欠定测向。

为了解决欠定测向问题，互质阵被提出。不同于最小冗余线阵，互质阵提供了一种系统化稀疏阵架构方案，来实现大于阵元数的阵列自由度，从而为欠定测向提供可能，成为了学术界的研究热点。欠定测向需要在阵列的协方差域将物理阵列转化为自由度更高的虚拟阵列，然后在虚拟阵列上实现测向。但互质阵的虚拟阵列存在孔洞，是一个非连续阵列。文献：Pal P,Vaidyanathan P P.Coprime sampling and the music algorithm[C].2011Digital Signal Processing and Signal Processing Education Meeting(DSP/SPE).IEEE,2011中提出的了一种基于空间平滑的多重信号分类的测向方法，简称SS-MUSIC。但该方法只能选取互质阵的虚拟阵列中的连续阵元，造成了阵列信息提取不完整以及自由度利用不充分。文献：Zhou C,Gu Y,Fan X,et al.Direction-of-arrivalestimation for coprime array via virtual array interpolation[J].IEEETransactions on Signal Processing.2018,66(22):5956-5971中提出基于原子范数最小化的方法对互质阵的虚拟阵的孔洞进行插值(简称原子范数)，从而可以利用全部自由度。但该方法未考虑有限快拍数导致的阵列模型噪声，其测向性能仍有提升的空间。文献：YangJ,Liao G,Li J.An efficient off-grid DOA estimation approach for nested arraysignal processing by using sparse Bayesian learning strategies[J].SignalProcessing.2016,128:110-122中利用离网格稀疏贝叶斯的方法(简称离网格贝叶斯)，考虑了各阵元模型噪声及其相关性及非高斯性。但由于欠定条件下阵列噪声方差难以估计，该方法为了避免估计阵列噪声方差舍弃了一个阵元自由度，造成信息丢失，导致测向性能降低。而且该方法本质上仍属于有网格算法，一方面网格效应会影响测向精度，另一方面需要进行空域网格扫描以输出测向结果。

发明内容

针对现有方法的不足，本发明提出一种基于零化去噪技术的互质阵欠定测向方法。该方法可有效利用阵列全部信息及自由度，并且对模型噪声进行建模，属于无网格算法，可明显改善测向性能。

本发明是通过以下技术方案实现的：一种基于零化去噪技术的互质阵欠定测向方法，所述方法包含以下步骤：

(S1)：设置互质正整数M和N，构造两个向量[0,N,2N,...,(M-1)N]，[M,2M,...,(N-1)M]，并将其合并得到向量

在一维坐标轴上以dp中的元素值为位置坐标放置M+N-1个阵元形成互质阵，其中d为控制阵元间距的参数。在远场窄带信号个数为K条件下，其中K≥M+N-1，采集T个快拍的阵列输出数据，形成数据向量

t＝1,2...,T。

(S2)：设置标识矩阵P，使P的第m行第n列的元素满足P_m,n＝m-n。令p_m＝vec(P)，将p_m中的元素从小到大排列，再找出其中所有第一次出现的元素及其对应位置索引，分别形成向量p_v和索引向量q。此时p_v和q中的元素个数均为MN+M+N-2，其中vec(·)表示将矩阵的各列依次堆叠形成向量。

(S3)：构造均匀虚拟阵对应的数据向量z及对应的模型噪声的协方差矩阵Σ：计算阵列协方差矩阵的估计值

并将其向量化得到

其中(·)^H表示共轭转置。设

令M′＝max(M(N-1),N(M-1))，

i从-M′遍历到M′，如果i∈p_v，则令

Σ[:,i]＝Ξ[:,q[i]]，Σ[i,:]＝Ξ[q[i],:]，否则，令z[i]＝0，Σ[:,i]＝0，Σ[i,:]＝0，其中z[i]表示z的第i个元素，Σ[:,i]表示Σ的第i列，Σ[i,:]表示Σ的第i行，其余方法表示类似，

表示克罗内克积。

(S4)：设置迭代计数变量l＝0，初始化模型噪声ε＝0，阵列噪声方差

零化系数h为矩阵L(z)L^H(z)的最小特征值对应的特征向量，其中L(·)为Toeplitz化算子，有

(S5)：解线性等式约束最小二乘问题：

其中

(·)^T表示转置，

为除了正中间一个元素为1其余为0的列向量，ω为第一个元素为1其余元素为0的列向量，T(·)为另一种Toeplitz化算子，有

求解完成后得到差分量Δε，

和Δh。

(S6)：更新ε，

和h，ε＝ε+Δε，

h＝h+Δh；判断迭代计数变量l是否达到上限L或零化系数h是否收敛，若有一个条件满足，则进入(S7)。若两个条件都不满足，迭代计数变量l＝l+1，再用新的ε，

和h更新γ和J，并返回(S5)。

(S7)：求方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0的K个根，为α_k,k＝1,2,...,K，则K个信号的测向结果为θ_k＝arcsin(angle(α_k)λ/(2πd))),k＝1,2,...,K，其中angle(·)表示求复数的幅角，λ为信号波长。

本发明相比现有技术的有益效果为：

首先，在欠定求解过程中引入了对阵列噪声方差的估计，且用零化去噪的技术实现了对互质阵虚拟阵列中的孔洞进行了插值，充分利用了所有的阵元信息以及阵列自由度；其次，本发明对由有限快拍数引起的模型噪声进行了建模，可提高测向精度；最后，本发明提出的方法属于无网格算法，不存在网格效应，也无需进行空域网格扫描。

附图说明

图1为本发明的流程图；

图2为本发明中的互质阵布阵形式；

图3为本发明欠定测向结果的空间谱与离网格贝叶斯方法的对比；

图4为本发明方法与其他方法的欠定测向均方根误差在不同信噪比下的对比；

图5为本发明方法与其他方法的欠定测向均方根误差在不同快拍数下的对比。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步详细说明。参照图1，本发明的具体实施步骤如下：

(S1)：设置互质正整数M和N，构造两个向量[0,N,2N,...,(M-1)N]，[M,2M,...,(N-1)M]，并将其合并得到向量

在一维坐标轴上以dp中的元素值为位置坐标放置M+N-1个阵元形成互质阵，其中d为控制阵元间距的参数，如图2所示。d的取值根据测向的视场角决定，若视场角为范围-90°～90°，则一般取d＝λ/2，λ为信号波长。在远场窄带信号个数为K条件下，其中K≥M+N-1，采集T个快拍的阵列输出数据，形成数据向量

t＝1,2...,T。此时x(t)可建模为x(t)＝As(t)+n(t) (1)

其中，s(t)为包含K个信号的列向量，K个信号互不相关；n(t)为阵列噪声列向量，阵列噪声与来波信号互不相关，各阵元噪声满足独立复高斯分布，即

为噪声方差，I_M+N-1为有M+N-1个对角元素的单位阵；A＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流型，a(θ_k)为阵列流型向量，有

其中，(·)^T表示转置，p_i表示p的第i个元素，

(S2)：设置标识矩阵P，使P的第m行第n列的元素满足P_m,n＝m-n。令p_m＝vec(P)，将p_m中的元素从小到大排列，再找出其中所有第一次出现的元素及其对应位置索引，分别形成向量p_v和索引向量q。此时p_v和q中的元素个数均为MN+M+N-2，其中vec(·)表示将矩阵的各列依次堆叠形成向量。这一步属于数据预处理操作，为下一步骤服务。其中d与p_v的乘积dp_v表示互质阵的虚拟阵列阵元的位置向量，注意这里的p_v中的整数不连续，说明互质阵的虚拟阵列中存在孔洞。

(S3)：构造均匀虚拟阵对应的数据向量z及对应的模型噪声的协方差矩阵Σ：计算阵列协方差矩阵的估计值

并将其向量化得到

其中(·)^H表示共轭转置。设

令M′＝max(M(N-1),N(M-1))，

i从-M′遍历到M′，如果i∈p_v，则令

Σ[:,i]＝Ξ[:,q[i]]，Σ[i,:]＝Ξ[q[i],:]，否则，令z[i]＝0，Σ[:,i]＝0，Σ[i,:]＝0，其中z[i]表示z的第i个元素，Σ[:,i]表示Σ的第i列，Σ[i,:]表示Σ的第i行，其余方法表示类似，

表示克罗内克积。

根据式(1)中的模型，我们可以得到阵列协方差矩阵向量化后的真实值r为

其中，

(·)^*表示求共轭。

中的元素分别表示各个入射信号的功率。可见，由于快拍数有限，r与其估计值

存在误差Δr，即

这里将该误差称为模型噪声，以便与之前的阵列噪声相区别。这里模型噪声满足渐进高斯分布，即可近似表示为

从式(3)中A_v的表达式中可以发现，r中存在重复元素，所以我们需要将其去除。对应的Δr的协方差矩阵也要去除相应的行与列。但去除完重复的元素后，由于互质阵的虚拟阵存在孔洞，对应的p_v向量中存在空缺整数。这里我们将虚拟阵列的孔洞位置假想为存在阵元，并预设阵元输出数据为0，对应的Δr的协方差矩阵的相应行与列也设为0。以上的操作用步骤(S3)中操作即可完成。操作完成后我们将得到孔洞填充后的均匀虚拟阵对应的数据向量z及对应的模型噪声的协方差矩阵Σ，这里模型噪声Δr在上述操作后将会变为下一步的ε。

(S4)：设置迭代计数变量l＝0，初始化模型噪声ε＝0，阵列噪声方差

零化系数h为矩阵L(z)L^H(z)的最小特征值对应的特征向量，其中L(·)为Toeplitz化算子，有

由于存在未知变量ε，

h，我们准备进行迭代求解。迭代求解需要进行变量初始化，这里对零化系数h的初始化采用的是在无噪零化关系L(z)h＝0下采用总体最小二乘实现的。

(S5)：解线性等式约束最小二乘问题：

其中

(·)^T表示转置，

为除了正中间一个元素为1其余为0的列向量，ω为第一个元素为1其余元素为0的列向量，T(·)为另一种Toeplitz化算子，有

求解完成后得到差分量Δε，

和Δh。

在每一步迭代求解ε，

h时，我们采用的是高斯牛顿迭代的方式，求解其差分量Δε，

和Δh。在上述线性等式约束最小二乘问题中，我们最小化的是

而不是

这是由于Σ^-1/2(ε+Δε)是一种白化操作，可以去除变量间的相关性，可使求解更加精确。这里

中的

是一种去噪操作，去除的噪声包括阵列噪声成分

和模型噪声ε。线性等式约束最小二乘问题中第一个等式约束代表零化关系，第二个等式约束是为了防止零解以及保证解的唯一性。解该线性等式约束最小二乘问题可用拉格朗日数乘法，或直接调用现成的相关算法包。

(S6)：更新ε，

和h，ε＝ε+Δε，

h＝h+Δh；判断迭代计数变量l是否达到上限L或零化系数h是否收敛，若有一个条件满足，则进入(S7)。若两个条件都不满足，迭代计数变量l＝l+1，再用新的ε，

和h更新γ和J，并返回(S5)。

在每一步迭代过程中，算法通过求解一个线性等式约束最小二乘问题得到差分量Δε，

和Δh来更新ε，

h，并且判断迭代是否可以结束，若不可以，则继续迭代。

(S7)：求方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0的K个根，为α_k,k＝1,2,...,K，则K个信号的测向结果为θ_k＝arcsin(angle(α_k)λ/(2πd))),k＝1,2,...,K，其中angle(·)表示求复数的幅角，λ为信号波长。

这里最终的测向只需要对迭代求解后的零化系数h构成的方程求根，再加少量的运算即可完成，所以该方法属于无网格算法，没有网格效应的影响，也无需进行空域网格扫描。

为了验证本发明提出的基于零化去噪技术的互质阵欠定测向方法的正确性和相对于现有技术的优越性，做以下仿真实验。

考虑互质阵M＝3，N＝5，d＝λ/2的情形，于是此时共有7个阵元。对于本发明提出的方法的(S6)步骤中，最大迭代次数设为L＝100，当||h-h′||₂/||h′||₂≤10^-8则认为h已收敛，其中h′为上一次迭代得到的零化系数的值。采用的性能对比方法为背景技术中提到的SS-MUSIC、离网格贝叶斯和原子范数三种方法。

实验一：

设置10个信号，其来波方向均匀分布在-55°～55°之间。信噪比为20dB，快拍数为500个。由于只有7个阵元，所以此时为欠定测向。此时的信号源数已经超过了SS-MUSIC算法理论可解的数量，所以未画出其测向结果。原子范数可解此问题，此处也未画出。这里画出了本发明的方法和离网格贝叶斯的测向结果对比，分别如图3(a)、图3(b)。图中虚线表示来波方向真实值，实线表示测向的空间谱。本发明的空间谱由零化系数h求z变换再求倒数得到。可以发现，本发明的空间谱在真实来波方向附近形成了尖锐的峰，准确测到了10个信号的来波方向，而离网格贝叶斯的只正确测到了一小部分，且出现了很多的伪峰。这是由于离网格贝叶斯算法为了避免估计阵列噪声方差，而牺牲了阵元信息和阵列自由度，导致测向错误。

实验二：

设置快拍数为500，信噪比从-10dB扫描至30dB。由于SS-MUSIC此处最大可解信号源数为7个，所以设置7个信号源，来波方向均匀分布在-55°～55°之间。此时信号源数等于阵元数，仍为欠定测向。仿真结果如图4所示，所示结果为500次蒙特卡洛实验的平均。可以发现，本发明的方法测向结果的均方根误差最小。离网格贝叶斯在高信噪比下均方根小于原子范数和SS-MUSIC，但在低信噪比下变差。原子范数在高信噪比下均方根误差大于离网格贝叶斯，这是由于其未考虑模型噪声。SS-MUSIC的均方根误差较大，其原因在与它只用了互质阵虚拟阵列的连续阵元信息。

实验三：

设置快拍数为0dB，快拍数从100个扫描至1000个，信号源设置同上。仿真结果如图5所示，所示结果为500次蒙特卡洛实验的平均。可以发现，本发明的方法测向结果在各个快拍数下均方根误差都最小。离网格贝叶斯在大快拍数下均方根小于原子范数和SS-MUSIC，但在小快拍数下变差。原子范数在大快拍数下均方根误差大于离网格贝叶斯。SS-MUSIC的均方根误差较大。

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

Claims (1)

1.一种基于零化去噪技术的互质阵欠定测向方法，其特征在于，包括以下步骤：

(S1)：设置互质正整数M和N，构造两个向量[0,N,2N,...,(M-1)N]，[M,2M,...,(N-1)M]，并将其合并得到向量

在一维坐标轴上以dp中的元素值为位置坐标放置M+N-1个阵元形成互质阵，其中d为控制阵元间距的参数；在远场窄带信号个数为K条件下，其中K≥M+N-1，采集T个快拍的阵列输出数据，形成数据向量

(S2)：设置标识矩阵P，使P的第m行第n列的元素满足P_m,n＝m-n；令p_m＝vec(P)，将p_m中的元素从小到大排列，再找出其中所有第一次出现的元素及其对应位置索引，分别形成向量p_v和索引向量q；此时p_v和q中的元素个数均为MN+M+N-2，其中vec(·)表示将矩阵的各列依次堆叠形成向量；

(S3)：构造均匀虚拟阵对应的数据向量z及对应的模型噪声的协方差矩阵Σ：计算阵列协方差矩阵的估计值

并将其向量化得到

其中(·)^H表示共轭转置；设

令M′＝max(M(N-1),N(M-1))，

i从-M′遍历到M′，如果i∈p_v，则令

Σ[:,i]＝Ξ[:,q[i]]，Σ[i,:]＝Ξ[q[i],:]，否则，令z[i]＝0，Σ[:,i]＝0，Σ[i,:]＝0，其中z[i]表示z的第i个元素，Σ[:,i]表示Σ的第i列，Σ[i,:]表示Σ的第i行，其余方法表示类似，

表示克罗内克积；

(S4)：设置迭代计数变量l＝0，初始化模型噪声ε＝0，阵列噪声方差

零化系数h为矩阵L(z)L^H(z)的最小特征值对应的特征向量，其中L(·)为Toeplitz化算子，有

(S5)：解线性等式约束最小二乘问题：

ω^T(h+Δh)＝1

其中

(·)^T表示转置，

为除了正中间一个元素为1其余为0的列向量，ω为第一个元素为1其余元素为0的列向量，T(·)为另一种Toeplitz化算子，有

求解完成后得到差分量Δε，

和Δh；

(S6)：更新ε，

和h，ε＝ε+Δε，

h＝h+Δh；判断迭代计数变量l是否达到上限L或零化系数h是否收敛，若有一个条件满足，则进入(S7)；若两个条件都不满足，迭代计数变量l＝l+1，再用新的ε，

和h更新γ和J，并返回(S5)；

(S7)：求方程h[1]α^K+h[2]α^K-1+...+h[K]α+h[K+1]＝0的K个根，为α_k,k＝1,2,...,K，则K个信号的测向结果为θ_k＝arcsin(angle(α_k)λ/(2πd))),k＝1,2,...,K，其中angle(·)表示求复数的幅角，λ为信号波长。

Images