CN109143154A - 一种应用于l型阵列的信号二维doa与频率联合估计方法 - Google Patents

Abstract

本发明公开了一种应用于L型阵列的信号二维DOA与频率联合估计方法，包括如下步骤：(1)将L型阵列接收信号构建为平行因子三线性模型；(2)对模型进行初等行变换，并根据压缩感知理论进行分块压缩；(3)利用三线性交替最小二乘算法，对压缩后的模型进行循环迭代分解，直到收敛，以得到参数矩阵的估计；(4)构造频率与角度估计的完备字典，并依次利用稀疏恢复方法获得信号的频率和二维DOA估计。本发明的有益效果为：(1)本发明适用于均匀L阵和非均匀L阵，阵列在实际应用中易于实现；(2)本发明算法计算量较小，同时对数据存储的容量要求较低；(3)本发明算法借助于平行因子模型，能够得到配对的角度和频率估计，无需再进行额外配对。

Description

一种应用于L型阵列的信号二维DOA与频率联合估计方法

技术领域

本发明涉及阵列信号处理多参数联合估计技术领域，尤其是一种应用于L型阵列的信号二维DOA与频率联合估计方法。

背景技术

近年来，阵列信号处理在雷达、声呐、无线通信等多个领域得到了迅速的发展，而接收信号多参数联合估计作为阵列信号处理中的一个重要的研究方向，也受到了相关学者们的广泛关注。阵列接收信号通常包含有波达方向DOA、多普勒频率等信息，如何准确的从接收信号中得到这些参数的联合估计，已成为阵列信号处理算法研究的热点。

近年来，国内外学者已经提出了许多关于接收信号DOA与频率联合估计的算法。其中包括传统的多重信号分类(Multiple signal classification，MUSIC)方法、借助旋转不变性估计信号参数(Estimation of signal parameters via rotational invariancetechniques，ESPRIT)的算法和传播算子(Propagator method，PM)算法。

平行因子(Parallel factor，PARAFAC)方法，最早于生理学中提出，用于对多维矩阵的分解，近几年，该方法被逐渐的运用于阵列信号参数估计中。平行因子方法通常需要对接收信号模型进行循环迭代分解，以获得参数矩阵的估计。该过程一般复杂度较高，特别是在信号接收阵列较大或信号采样数较多的情况下。压缩感知(Compressedsensing，CS)是一项用于对高维数据分析的技术，借助于压缩感知思想，可以对高维数据进行压缩，以达到减少采样数据和降低数据存储空间的目的。因此，将压缩感知理论和平行因子模型相结合，可以有效地降低平行因子方法的计算复杂度和对数据存储空间的要求。

发明内容

本发明所要解决的技术问题在于，提供一种应用于L型阵列的信号二维DOA与频率联合估计方法，能够得到性能与传统平行因子方法性能接近的角度和频率估计，同时其计算复杂度比传统的平行因子方法低。

为解决上述技术问题，本发明提供一种应用于L型阵列的信号二维DOA与频率联合估计方法，包括如下步骤：

(1)将L型阵列接收信号构建为平行因子三线性模型；

(2)对模型进行初等行变换，并根据压缩感知理论进行分块压缩；

(3)利用三线性交替最小二乘算法，对压缩后的模型进行循环迭代分解，直到收敛，以得到参数矩阵的估计；

(4)构造频率与角度估计的完备字典，并依次利用稀疏恢复方法获得信号的频率和二维DOA估计。

优选的，步骤(1)具体为：

L型阵列接收端的信号模型可以构建为z₀(t)＝As(t)+n(t)，其中A为整个阵列方向矩阵，n(t)为接收信号噪声，s(t)为信源矢量；为了对接收信号的频率进行估计，对阵列引入多级时延接收，考虑P-1级时延，整个阵列J个快拍时延接收数据可以表示为

其中，F为时延矩阵，表示矩阵的Khatri-Rao积，为接收噪声。

优选的，步骤(2)具体为：

对信号Z进行初等行变换，得到其中，G∈C^{(M+N-1)P×(M+N-1)P}为行变换矩阵，N_I为变换后的噪声；在无噪声环境下，变换后的接收信号Z_I可以表示为平行因子模型；为了降低三线性分解算法的复杂度，将接收数据三维矩阵Z_I ∈C^{(M+N-1)×J×P}压缩为一个较小的矩阵Z_I′∈C^{(M′+N′)×J′×P′}，压缩后的接收信号Z_I′∈C^{(M′+N′)J′×P′}可以表示为

其中，A′、F′、S′和N_I′分别为压缩后的方向矩阵，时延矩阵，信源矩阵和噪声；将Z_I′重排可以得到其他两种形式Z′_II和Z′_III。

优选的，步骤(3)具体为：

根据三线性最小二乘算法，参数矩阵S′、F′、A′的最小乘更新分别为其中Z_I′、Z′_II和Z′_III是压缩后的含噪信号，和为S′、F′、A′的估计值；利用上述更新方法，循环迭代分解得到S′、F′、A′的最终估计值。

优选的，步骤(4)具体为：

定义矩阵的第k列为则有其中，f_k为矩阵F的第k行，n_1k为相应的估计误差，为尺度因子；

利用采样序列构造关于频率估计的完备字典并将上式改写为可以利用l₀范数约束可以得到f_k的估计为

设矩阵和分别为和的第k列，则有 其中，a_xk和a_yk分别为A_x和A_y的第k列，n_2xk和n_2yk为相应的噪声，和缩放因子；

定义设和构造角度估计完备字典 则同频率估计过程类似，u_k和v_k可以如下通过l₀范数获得

最后，获得θ_k和的估计为

本发明的有益效果为：(1)本发明算法适用于均匀L阵和非均匀L阵，阵列在实际应用中易于实现；(2)本发明算法计算量较小，同时对数据存储的容量要求较低；(3)本发明算法借助于平行因子模型，能够得到配对的角度和频率估计，无需再进行额外配对。

附图说明

图1为本发明的L型接收阵列示意图。

图2为本发明的阵列多级时延接收示意图。

图3为本发明的接收信号模型压缩示意图。

图4为本发明与传统平行因子算法复杂度对比示意图。

图5(a)为本发明二维角度和频率估计结果分布图。

图5(b)为本发明方位角和仰角估计结果分布图。

图5(c)为本发明方位角和频率估计结果分布图。

图6(a)为本发明与三种常用算法角度估计性能对比图。

图6(b)为本发明与三种常用算法频率估计性能对比图。

图7(a)为本发明在不同快拍数下的角度估计性能对比图。

图7(b)为本发明在不同快拍数下的频率估计性能对比图。

具体实施方式

本发明提出一种对于L型阵列所接收的非相干信号进行二维DOA和频率联合估计的方法。该方法首先为了对接收信号频率进行估计，对信号进行多级时延接收，并将多级时延输出信号构建为平行因子三线性模型。为了降低对该模型进行分解的复杂度，该算法结合压缩感知理论，对该模型进行分块压缩，然后再利用三线性交替最小二乘算法对该模型分解，以得到信号参数矩阵的估计，并从所获得参数矩阵中得到信号二维DOA和频率估计。该算法能够同时应用于均匀L阵和非均匀L阵，同时在分解过程中，参数矩阵具有相同的列模糊，因此其最终所得到的角度和频率估计是配对的。根据仿真结果，本发明所提出的算法角度和频率估计性能均接近于传统平行因子方法，同时优于一些常用的多参数联合估计算法，同时其复杂度比传统的平行因子方法低。

(1)信号模型

如图1所示的L型阵列，接收到K个信号。该阵列由位于x轴上的M个传感器和位于y轴上的N个传感器组成，两组传感器的参考阵元在原点处重合。x轴上第m个阵元到参考阵元的间距为y轴上第n个阵元到参考阵元的间距为

假设入射信号均为非相干的远场信号，所含噪声为零均值的加性高斯白噪声且与信号独立。第k个信号的中心频率为f_k，波达方向为其中和θ_k分别表示信号的方位角和仰角。x轴和y轴上的传感器接收信号分别可以表示为

x(t)＝A_xs(t)+n_x(t)

y(t)＝A_ys(t)+n_y(t)

其中，s(t)＝[s₁(t),s₂(t),…,s_K(t)]^T为信源矢量， 为x轴上M个阵元的方向矩阵， 为y轴上N-1个阵元(不含参考阵元)的方向矩阵，n_x(t)∈C^M×1和n_y(t)∈C^(N-1)×1为接收噪声，c为信号传播速度。整个阵列t时刻接收信号可以表示为

其中为整个阵列方向矩阵，为接收信号噪声。

为了对接收信号的频率进行估计，我们对阵列引入多级时延接收。如图2所示，考虑P-1级时延，第p级时延时间τ_p＝pτ，其中τ满足关系

对于第p级时延τ_p，整个阵列t时刻的接收信号为

其中，D_P+1(F)表示取矩阵F第p+1行作为对角元素组成对角矩阵。F为时延矩阵，定义为

n_p(t)＝n(t-pt)为时延接收噪声，则第p级时延J个快拍的接收数据为

其中，S^K×J＝[s(t₁),s(t₂),...,s(t_J)]，N_p＝[n_p(t₁),n_p(t₂),...,n_p(t_J)]。整个阵列全部时延接收数据Z^(M+N-1)P×J可以表示为

其中，表示矩阵的Khatri-Rao积，

(2)模型变换与压缩

对接收信号Z进行初等行变换，得到

其中，G∈C^{(M+N-1)P×(M+N-1)P}为行变换矩阵，定义为

在无噪声环境下，变换后的接收信号Z_I可以表示为平行因子模型

其中，F(p,k)表示矩阵F的第(p,k)元素，S(j,k)表示矩阵S的第(j,k)元素，A(m+n,k)表示矩阵A的第(m+n,k)元素。根据平行因子模型的结构特征，我们Z_I进行重排，可以得到另外两个矩阵

为了降低三线性分解算法的复杂度，我们将接收数据三维矩阵Z_I ∈C^{(M+N-1)×J×P}压缩为一个较小的矩阵Z_I′∈C^{(M′+N′)×J′×P′}。其中，M<M′、N′<N、J′<J以及P′<P。图3为压缩过程示意图。

圧缩过程中的四个压缩矩阵为U₁∈C^M×M′(M′<M)、U₂∈C^(N-1)×N′(N′<N-1)、V∈C^P×P′(P′<P)和W∈C^J×J′(J′<J)。上述压缩矩阵需满足有限等距性质可以通过Tucker3分解或者一些随机矩阵如高斯随机矩阵、伯努随机矩阵或局部傅里叶矩阵得到。

Z_I中的接收信号可以写为

其中，分别为x轴和y轴上阵元接收数据。压缩后的接收信号Z_I′∈C^{(M′+N′)J′×P′}可以由下式得到

其中，表示矩阵的Kronecker积，F′＝V^TF、S′＝W^TS。压缩后的噪声为

压缩后的接收信号Z_I′也可以表示为平行因子模型，同时我们可以得到压缩信号另外两种重排矩阵

其中，N′_II和N′_III为压缩后的噪声。

(3)模型分解

根据得到最小二乘拟合为

其中，||·||_F表示矩阵的Forbenius范数。由上式得到S′的最小二乘更新为

其中，表示矩阵的广义逆，和分别为上一次迭代中所得到的A′和F′的估计值。

同样的，根据得到

由上式得到的最小二乘更新为

其中，和分别为上一次迭代中所得到的S′和A′的估计值。

由得到最小二乘拟合为

由上式得到A′的最小二乘更新为

其中，和分别为之前所得到的F′和S′的估计值。

上述已经给出了三线性模型TALS算法分解过程。为了判断上述迭代分解过程是否收敛，我们计算矩阵其中，和为分解过程中A′、F′和S′的估计值。定义第k次迭代中参数矩阵估计的残差平方和(Sum of Squared Residuals，SSR)为

其中，e_ji为矩阵E的第(j,i)元素。定义迭代收敛速率SSR_rate＝(SSR_k-SSR_k-1)/SSR_k-1。根据上述过程，我们不停的迭代更新和的值，当SSR_rate小于某个设定的较小值时，即可认为上述分解过程已经收敛，所得到的和即为A′、F′和S′的最终估计值，表示为

其中，Π为列交换矩阵，Δ₁、Δ₂和Δ₃为尺度模糊矩阵，满足Δ₁Δ₂Δ₃＝I。W₁、W₂和W₃为估计误差。其中尺度模糊问题我们可以通过归一化方法消除，而列模糊对于我们的角度和频率估计是没有影响的。

(4)频率与二维DOA估计

上述我们已通过三线性分解得到压缩后的矩阵F′、A′和S′的估计，而通过稀疏恢复即可从所得到的压缩矩阵中得到信号的二维DOA和频率估计。

a)频率估计

定义矩阵的第k列为根据以及F′=V^TF，可以得到

其中，f_k为矩阵F的第k行，n_1k为相应的估计误差，为尺度因子。

假设为一个关于频率的采样序列，其中R>>K。利用该序列我们可以构造关于频率估计的完备字典其中， 如果满足则存在一个稀疏矩阵e_r∈C^R×1满足其中e_r第r个元素为1，其他元素为0。利用e_r，我们可以将上式改写为

可以利用l₀范数约束可以得到e_r为

上式可以进一步简化为

利用则可以得到f_k的估计为

b)二维DOA估计

根据有

其中，为矩阵的前M′行。为矩阵的后N′行。设矩阵和分别为和的第k列，则有

其中，a_xk和a_yk分别为A_x和A_y的第k列，n_2xk和n_2yk为相应的噪声,和缩放因子。

定义设和为潜在信源位置所组成的序列，其中G>>K、Q>>K。构造两个范德蒙德矩阵其中 和可以看作二维DOA估计的完备字典,如果有则存在两个稀疏矩阵e_g∈C^G×1、e_q∈C^Q×1，分别满足则和可以写为

则同频率估计过程类似，u_k和v_k可以如下通过l₀范数获得

最后，我们可以获得θ_k和的估计为

其中，sin^-1(·)为反正弦函数。

(5)算法复杂度分析与克拉美罗界(CRB)

下面分析本发明所提出的算法的复杂度，以及L型阵列角度和频率估计的克拉美罗界(Cramér-Rao bound，CRB)。

a)复杂度分析

对于所提出的交合和频率估计算法，接收信号压缩过程复杂度为O[P′(M′+N′)P(M+N-1)J+JP′(M′+N′)J′]，三线性分解过程每次迭代的复杂度为O[3K²+K²(J′+M′+N′+P′)+6K²((M′+N′)(P′+J′)+J′P′)+3K(M′+N′)J′P′]，以及信号稀疏恢复过程复杂度为O[K(MG+Q(N-1)+PR)]，因此所提算法总复杂度为O[P′(M′+N′)P(M+N-1)J+JP′(M′+N′)J′+n₁(3K²+K²(J′+M′+N′+P′)+6K²((M′+N′)(P′+J′)+J′P′)+3K(M′+N′)J′P′)+K(MG+Q(N-1)+PR)]，其中n₁为分解过程迭代次数。而对于传统的平行因子算法，复杂度为O(n₂(3K²+K²(J+M+N-1+P)+6K²((M+N-1)(P+J)+JP)+3K(M+N-1)JP)+2K²(M+N-1+P)+9K³)，其中n₂为分解迭代次数。为了直观的比较两种算法的计算复杂度，我们考虑M＝N＝P＝10、K＝3，两种算法分解迭代次数大约为几十次，我们为n₁＝n₂＝30，以及G＝Q＝200、R＝300，M′/M＝N′/N＝J′/J＝0.5。图4给出了传统平行因子算法和本发明所提算法在不同的快拍数下的计算复杂度对比，从图4中可以看出，所提算法由于对接收信号进行了压缩，计算法复杂度大大降低。

b)克拉美罗界

对于L型阵列二维DOA与频率估计，其克拉美罗界为

其中，b_k为矩阵B的第k列。其中1₃为3×3元素全为1的矩阵，

(6)仿真结果

我们假设L型阵列接收到3个远场非相干信号，信号二维DOA和频率分别为J表示信号快拍数，P为时延级数，M和N分别为x轴和y轴上的阵元数。仿真中，我们利用信号角度和频率估计的求根均方误差(Root mean square error,RMSE)来评估所提算法的参数估计性能，其中，角度和频率估计的RMSE分别定义为

其中，θ_k和f_k分别为第k个接收信号的准确方位角、仰角和频率。和分别为第l次蒙特卡洛仿真中θ_k和f_k的估计值。L为总的仿真次数，下面的仿真中我们取L＝1000。

仿真1：图5(a)、图5(b)和图5(c)为发明所提算法在信噪比为10dB时的二维DOA和频率估计结果分布图。定义算法中平行因子模型的大小为(M+N-1)×P×J。仿真1中，接收信号原始模型大小为(10+10-1)×10×200，压缩后变为(5+5)×5×100((M′+N′)×P′×J′)。x轴和y轴上的阵元到参考阵元的距离分别为d^x＝[0,60,100,162,190,265,298,364,410,450]、d^y＝[0,50,95,150,210,256,294,360,410,455]。从图5(a)、图5(b)和图5(c)可以看出，本发明所提算法可以有效用于L阵中的二维DOA和频率估计。

仿真2：图6(a)和图6(b)为所提算法与传统平行因子算法、ESPRIT算法以及PM算法二维DOA和频率估计性能对比图。仿真中，所接收的信号原始模型大小为(10+10-1)×10×200((M+N-1)×P×J)，压缩后为(8+8)×8×160((M′+N′)×P′×J′)。从图6(a)和图6(b)可以看出本章算法角度和频率估计均接近于传统的平行因子算法，同时优于ESPRIT算法和PM算法。

仿真3：图7(a)和图7(b)为所提算法在不同的快拍数J下的角度和频率估计性能。仿真3中，阵列大小为M＝10、N＝10，时延级数P＝10。分别设置快拍数J＝100、J＝200和J＝300。从图7(a)和图7(b)可以看出所提算法角度和频率估计性能会随着J的增大而变好。

本发明提出一种对于L型阵列所接收的非相干信号进行二维DOA和频率联合估计的方法。首先，将L型阵列接收信号构建为平行因子三线性模型，然后对该模型进行变换后再压缩，最后利用三线性交替最小二乘算法对该模型分解，以得到信号参数矩阵的估计，并从所获得参数矩阵中得到信号二维DOA和频率估计。由于上述分解过程所得到的参数矩阵的估计拥有相同的列模糊，因此，最终所获得的角度和频率估计已自动配对。同时该方法对阵列中的阵元间距没有要求，因此该方法对于均匀L阵和非均匀L阵都适用。根据仿真结果，本发明所提出的算法能够得到性能与传统平行因子方法性能接近的角度和频率估计，同时其计算复杂度比传统的平行因子方法低。

Claims (5)

1.一种应用于L型阵列的信号二维DOA与频率联合估计方法，其特征在于，包括如下步骤：

(1)将L型阵列接收信号构建为平行因子三线性模型；

(2)对模型进行初等行变换，并根据压缩感知理论进行分块压缩；

(3)利用三线性交替最小二乘算法，对压缩后的模型进行循环迭代分解，直到收敛，以得到参数矩阵的估计；

(4)构造频率与角度估计的完备字典，并依次利用稀疏恢复方法获得信号的频率和二维DOA估计。

2.如权利要求1所述的应用于L型阵列的信号二维DOA与频率联合估计方法，其特征在于，步骤(1)具体为：

L型阵列接收端的信号模型可以构建为z₀(t)＝As(t)+n(t)，其中A为整个阵列方向矩阵，n(t)为接收信号噪声，s(t)为信源矢量；为了对接收信号的频率进行估计，对阵列引入多级时延接收，考虑P-1级时延，整个阵列J个快拍时延接收数据可以表示为

Z＝[F⊙A]S+N

其中，F为时延矩阵，⊙表示矩阵的Khatri-Rao积，为接收噪声。

3.如权利要求1所述的应用于L型阵列的信号二维DOA与频率联合估计方法，其特征在于，步骤(2)具体为：

对信号Z进行初等行变换，得到其中，G∈C^{(M+N-1)P×(M+N-1)P}为行变换矩阵，N_I为变换后的噪声；在无噪声环境下，变换后的接收信号Z_I可以表示为平行因子模型；为了降低三线性分解算法的复杂度，将接收数据三维矩阵Z_I ∈C^{(M+N-1)×J×P}压缩为一个较小的矩阵Z′_I ∈C^{(M′+N′)×J′×P′}，压缩后的信号Z′_I∈C^{(M′+N′)J′×P′}可以表示为

Z′_I＝[A′⊙F′]S′^T+N′_I

其中，A′、F′、S′和N_I′分别为压缩后的方向矩阵，时延矩阵，信源矩阵和噪声；将Z′_I重排可以得到其他两种形式Z′_II和Z′_III。

4.如权利要求1所述的应用于L型阵列的信号二维DOA与频率联合估计方法，其特征在于，步骤(3)具体为：

根据三线性最小二乘算法，参数矩阵S′、F′、A′的最小乘更新分别为其中Z′_I、Z′_II和Z′_III是压缩后的含噪信号，和为S′、F′、A′的估计值；利用上述更新方法，循环迭代分解得到S′、F′、A′的最终估计值。

5.如权利要求1所述的应用于L型阵列的信号二维DOA与频率联合估计方法，其特征在于，步骤(4)具体为：

定义矩阵的第k列为则有其中，f_k为矩阵F的第k行，n_1k为相应的估计误差，为尺度因子；

利用采样序列构造关于频率估计的完备字典并将上式改写为可以利用l₀范数约束可以得到f_k的估计为

设矩阵和分别为和的第k列，则有 其中，a_xk和a_yk分别为A_x和A_y的第k列，n_2xk和n_2yk为相应的噪声，和缩放因子；

定义设和构造角度估计完备字典 则同频率估计过程类似，u_k和v_k可以如下通过l₀范数获得

最后，获得θ_k和的估计为