CN112924924B - 一种基于l1范数稀疏表示的宽带信号DOA估计方法 - Google Patents

Abstract

本发明公开了一种基于l₁范数稀疏表示的宽带信号DOA估计方法，通过优化子空间的正交性确定加权矩阵，以加权l₁范数作为最小化的目标函数来增加信号的稀疏性，进而使l₁范数更接近于l₀范数，解决了L1‑SVD算法在低信噪比时宽带DOA估计精度不高，伪峰较多的问题。本发明的算法主瓣更窄，旁瓣更低，能提高分辨力，并准确稳定地估计出宽带信号波达方向；特别在低信噪比时能有效抑制伪峰，且宽带DOA估计精度更高。

Description

一种基于l1范数稀疏表示的宽带信号DOA估计方法

技术领域

本发明属于电磁探测、无源定位及无线电监测管理技术领域，尤其涉及一种基于l₁范数稀疏表示的宽带信号低信噪比DOA估计方法。

背景技术

DOA估计在电磁环境环境探测领域有着十分重要的作用。近几十年来，DOA估计广泛用于雷达系统、声呐系统、无线通信系统和射电天文系统。随着无线电技术的发展，实际系统使用的频带越来越宽，这给DOA估计带来了新的挑战，因此宽带信号DOA估计研究意义重大。然而，传统的DOA估计算法都是基于窄带信号提出，宽带信号DOA估计会受频率影响。宽带信号DOA估计方法主要有不相关信号子空间方法(ISM)和相干信号子空间方法(CSM)。

目前常用的DOA估计方法是子空间分析方法，这些方法需要多通道数、多快拍及数据量大。近年来，由Donoho与Candes等人提出的压缩感知(CS)理论有潜力解决上述不足。该理论表明，当信号具有稀疏性时，可以用远低于奈奎斯特采样定理的速率对信号进行采样并且能够进行准确的恢复。其中，Malioutov等人提出的L1-SVD算法可用于宽带信号DOA估计，其首先将宽带信号均分为若干个窄带信号，对数据矩阵进行奇异值分解降维，再通过l₁范数构建稀疏信号模型，最后通过二阶锥规划(SOCP)化求解。还有学者提出了宽带协方差矩阵稀疏表示方法估计宽带信号DOA。这些方法在信噪比较高时效果较好，然而，大部分算法用l₁范数代替l₀范数近似求解，在低信噪比时，解的稀疏性变差，空间谱会有较多伪峰出现，随着信噪比的下降，这些方法性能明显降低。针对这一问题，现有技术根据子空间的正交性确定加权矩阵。由Capon谱设计了一个加权范数罚函数，能在一定程度上有效抑制伪峰，并估计出波达方向。这几种方法均对该问题有所，改进提高了低信噪比下宽带信号DOA精度。当信噪比更低时，这些方法效果不好。

本发明通过优化子空间的正交性确定加权矩阵，以加权l₁范数作为最小化的目标函数来增加信号的稀疏性，进而提高低信噪比下宽带DOA估计精度并抑制伪峰。

发明内容

为了解决上述已有技术存在的不足，本发明提出一种基于l₁范数稀疏表示的宽带信号低信噪比DOA估计方法，在低信噪比主瓣更窄，旁瓣更低，能提高分辨力和有效抑制伪峰，并准确稳定地估计出宽带信号波达方向。本发明的具体技术方案如下：

一种基于l₁范数稀疏表示的宽带信号DOA估计方法，包括以下步骤：

S1：接收阵列数据并做FFT，构成宽带信号的J个窄带阵列接收模型；

设P个远场宽带信号入射到一个阵元数为M的均匀阵列(P<M)，将整个观测周期内的信号等分成L段，每段包含K个快拍，每段的阵列接收信号进行时域的快速傅里叶变换后获得信号带宽内的J个窄带阵列接收信号为：

X(f_j)＝A(f_j,θ)S(f_j)+N(f_j)(j＝1,2,…,J) (1)

其中，X(f_j)＝[x₁(f_j),x₂(f_j),…,x_M(f_j)]^T是一个M×L矩阵,f_j是子窄带中心频率；

A(f_j)＝[a(f_j,θ₁),a(f_j,θ₂),…,a(f_j,θ_P)]是频率f_j处的方向矢量，是一个M×P矩阵；

d是两阵元间距，θ_p是入射角；

S(f_j)＝[s₁(f_j),s₂(f_j),…,s_p(f_j)]^T是频率f_j处P×L稀疏信号矢量；

N(f_j)＝[N₁(f_j),N₂(f_j),…,N_M(f_j)]^T是均值为零、方差为的加性高斯白噪声；

为方便描述，将式(1)简化为：

X_j＝A_jS_j+N_j(j＝1,2,…,J) (2)

S2：将频域数据分为J个频段，每个频段的数据矩阵奇异值分解降维，并预处理；

式(2)中，如果X_j能恢复出S_j，则能够根据S_j中的非零行的位置确定信源的DOA估计，即l₀范数问题；但l₀范数的优化是个NP-hard问题，因此通常将l₁范数代替l₀范数近似求解；

其中，是S_j中每一行的行向量的l₂范数组成的列向量，/> λ_j是受噪声影响的正则化参数，通常为一个常量；‖X_j-A_jS_j‖_F是将矩阵X_j-A_jS_j按列拉直并求l₂范数的结果，即：

L1-SVD算法先对接收数据进行奇异值分解，得到M×P的降维矩阵X_svj，设X_j的奇异值分解为：

X_j＝UΛV^H (5)

其中，U为X_j的左奇异向量，V为X_j的右奇异向量，Λ为包含奇异值的对角矩阵，则降维后的稀疏模型为：

X_svj＝X_jVD_P＝A_jS_jVD_P+N_jVD_P＝A_jS_svj+N_svj (6)

其中，S_svj＝S_jVD_P，N_svj＝N_jVDP，D_P＝[I_P 0]^H，I_P是一个P×P的单位矩阵,0为P×(L-P)的零矩阵；

对式(6)预处理；令R_x为原信号矢量X_j的协方差矩阵；式(6)转化为：

利用二阶锥规划即SOCP求解式(7)得：

其中，λ为正则化参数，p和q均为变量；

S3：根据ISM算法思想求出权值W；

由于权值是由噪声子空间与信号方向矢量正交得到：

b^H(θ_i)E_n＝0,i＝1,2,…,N (9)

其中，b(θ_i)是方向矢量，(·)^H是共轭转置，E_n是噪声子空间，N是空间[0°,180°]等间隔划分角度的个数；

权值w_i是将阵列方向矢量投影到噪声子空间的值除以投影到信号子空间的值：

其中，E_s是信号子空间，加权矩阵W为：

W＝diag{w_i} (11)

S4：将步骤S3得到的权值带入步骤S2的稀疏求解模型，分别求出每段的谱估计，谱估计平均即宽带信号DOA；

由式(11)和式(7)得到新的求解模型：

转化为二阶锥规划求解：

求出每段的谱估计，谱估计平均即宽带信号DOA。

本发明的有益效果在于：

1.本发明的方法通过优化子空间的正交性确定加权矩阵，以加权l₁范数作为最小化的目标函数来增加信号的稀疏性，进而提高低信噪比下DOA估计精度并抑制伪峰，解决了L1-SVD算法在低信噪比时宽带DOA估计精度不高，伪峰较多的问题。

2.本发明的算法在低信噪比的情况下主瓣更窄，旁瓣更低，能提高分辨力和有效抑制伪峰，并准确稳定地估计出宽带信号波达方向；特别在低信噪比时能有效抑制伪峰，且宽带DOA估计精度更高。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：

图1是本发明的信号对应的方向矢量与噪声对应的方向矢量的投影图；

图2是本发明实施例1的SNR＝-6dB时归一化谱图；

图3是本发明实施例1的SNR＝10dB时归一化谱图；

图4是本发明实施例1的SNR＝-6dB时4种算法分辨力；

图5是本发明实施例1的SNR与RMSE的关系；

图6是本发明实施例2的室外测量照片；

图7是本发明实施例2的SNR＝0dB户外测试数据谱图，其中，(a)是户外测试数据谱图，(b)是图(a)部分放大图；

图8是本发明实施例2的SNR＝10dB户外测试数据谱图，(a)是户外测试数据谱图，(b)是图(a)部分放大图。

具体实施方式

为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。

宽带阵列信号DOA估计分为分频处理与直接处理，本发明根据分频处理中ISM的思想，利用阵列稀疏表示求宽带DOA估计。即将宽带信号均分为J个频段，分别计算每个频段的DOA估计结果，后将这些DOA估计平均。

一种基于l₁范数稀疏表示的宽带信号低信噪比DOA估计方法，包括以下步骤：

S1：接收阵列数据并做FFT，构成宽带信号的J个窄带阵列接收模型；

设P个远场宽带信号入射到一个阵元数为M的均匀阵列(P<M)，将整个观测周期内的信号等分成L段，每段包含K个快拍，每段的阵列接收信号进行时域的快速傅里叶变换(FFT)后获得信号带宽内的J个窄带阵列接收信号为：

X(f_j)＝A(f_j,θ)S(f_j)+N(f_j)(j＝1,2,…,J) (1)

其中，X(f_j)＝[x₁(f_j),x₂(f_j),…,x_M(f_j)]^T是一个M×L矩阵,f_j是子窄带中心频率；

A(f_j)＝[a(f_j,θ₁),a(f_j,θ₂),…,a(f_j,θ_P)]是频率f_j处的方向矢量，是一个M×P矩阵；

d是两阵元间距，θ_p是入射角；

S(f_j)＝[s₁(f_j),s₂(f_j),…,s_p(f_j)]^T是频率f_j处P×L稀疏信号矢量；

N(f_j)＝[N₁(f_j),N₂(f_j),…,N_M(f_j)]^T是均值为零、方差为的加性高斯白噪声；

为方便描述，将式(1)简化为：

X_j＝A_jS_j+N_j(j＝1,2,…,J) (2)

S2：将频域数据分为J个频段，每个频段的数据矩阵奇异值分解降维，并预处理；

如果X_j能恢复出S_j，则能根据S_j中的非零行的位置来确定信源的DOA估计。式(2)在稀疏表示是l₀范数问题，但l₀范数的优化是个NP-hard问题，因此通常将l₁范数代替l₀范数近似求解；

其中，是S_j中每一行的行向量的l₂范数组成的列向量，/> λ_j是受噪声影响的正则化参数，通常为一个很小的常量；‖X_j-AjS_j‖_F是将矩阵X_j-A_jS_j按列拉直并求l₂范数的结果，即：

当快拍K较大时，必然导致运算量过大。L1-SVD算法先对接收数据进行奇异值分解，得到M×P的降维矩阵X_svj，设X_j的奇异值分解为：

X_j＝UΛ^VH (5)

其中，U为X_j的左奇异向量，V为X_j的右奇异向量，Λ为包含奇异值的对角矩阵，则降维后的稀疏模型为：

X_svj＝X_jVD_P＝A_jS_jVD_P+N_jVD_P＝A_jS_svj+N_svj (6)

其中，S_svj＝S_jVD_P，N_svj＝N_jVD_P，D_P＝[I_P o]^H，I_P是一个P×P的单位矩阵，0为P×(L-P)的零矩阵；；

对式(6)预处理；令R_x为原信号矢量X_j的协方差矩阵；式(6)转化为：

利用二阶锥规划即SOCP求解式(7)得：

其中，λ为正则化参数，p和q均为变量；

由于L1-SVD算法是由l₁范数代替l₀范数，其稀疏性难以保证，特别是低信噪比时存在较大的误差，导致低信噪比时伪峰较多，DOA估计精度较差。在l₁范数下的约束，求解的是模值最小情况下的解，此时稀疏信号S_j中大系数对应的模值大，小系数对应模值小。因此，可对稀疏信号S_j加权，提高解的稀疏性。

S3：根据ISM算法思想求出权值W；

由于权值是由噪声子空间与信号方向矢量正交得到：

b^H(θ_i)E_n＝0,i＝1,2,…,N (9)

其中，b(θ_i)是方向矢量，(·)^H是共轭转置，E_n是噪声子空间，N是空间[0°,180°]等间隔划分角度的个数；但由于噪声、信号处理等影响，实际并不等于0，而是一个很小的数。即信号对应的阵列方向矢量在噪声子空间的投影较小，噪声对应的阵列方向矢量在噪声子空间的投影较大。因此将投影的结果的幅值平方取倒数后，便可以在信号方向得到一个尖峰。

当信噪比变小时，信号子空间和噪声子空间的正交性会变差，谱估计的尖峰值会下降，权值的估计也会不准确，从而导致稀疏性变差。本发明提出权值w_i是将阵列方向矢量投影到噪声子空间的值除以投影到信号子空间的值：

其中，E_s是信号子空间，加权矩阵W为：

W＝diag{w_i} (11)

图1是信号对应的方向矢量与噪声对应的方向矢量的投影，噪声方向矢量在信号子空间E_s的投影N_s比较小，在噪声子空间E_n的投影N_n比较大；反之，信号方向矢量在信号子空间E_s的投影S_s比较大，在噪声子空间的E_n的投影S_n比较小。采用窦慧晶、高立菁、朱子云等提出的采用为加权值并改进奇异值矩阵求得DOA估计(详见窦慧晶,高立菁,朱子云.基于加权l₁范数稀疏信号表示的DOA估计[J].北京工业大学学报,2018 44(10):1297-1302，以下简称MG-L1-SVD)，加权时其峰值比为/>本发明提出算法的峰值比为/>而/>因此加权值增大了峰值比，减小了旁瓣，提高了解的稀疏性。

S4：将步骤S3得到的权值带入步骤S2的稀疏求解模型，分别求出每段的谱估计，谱估计平均即宽带信号DOA；

由式(11)和式(7)得到新的求解模型：

转化为二阶锥规划求解：

求出每段的谱估计，谱估计平均即宽带信号DOA。

为了方便理解本发明的上述技术方案，以下通过具体实施例对本发明的上述技术方案进行详细说明。

实施例1

F.L.Liu等提出采用为加权值求得DOA估计。(详见F.L.Liu,L.Peng,M.Wei1,P.P.Chen,and S.M.Guo.An improved L1-SVD algorithm basedon noise subspace for DOA estimation[J].Progress In Electromagnetics ResearchC,vol.29,pp.109-122,2012.以下简称M-L1-SVD.)

本实施例通过仿真实验对本发明的方法与L1-SVD、MG-L1-SVD、M-L1-SVD在不同信噪比下宽带信号DOA估计性能、分辨力、SNR与RSME的关系进行对比。在仿真实验中，宽带信号是OFDM信号，采用8均匀线阵，阵元间距为0.06m，快拍数为4。

首先，比较本发明的方法与L1-SVD、M-L1-SVD、MG-L1-SVD在不同信噪比下宽带DOA估计性能。一个OFDM信号中心频率为2048MHz，带宽为20.48MHz，来波方向为117°，分别在SNR＝-6dB及SNR＝10dB情况下对比测向性能。图2是SNR＝-6dB时4种算法的归一化谱图，图3是SNR＝10dB时4种算法的归一化谱图。在高信噪比时4种算法都能准确估计宽带DOA，本发明的方法主瓣宽度小，旁瓣最低。当信噪比较低时，本发明的方法效果最好，能抑制伪峰，其主瓣最窄，旁瓣最低。在低信噪比时其他三种算法几乎不能DOA估计。

其次，比较本发明的方法与L1-SVD、M-L1-SVD、MG-L1-SVD在低性噪比下宽带信号的分辨力。2个等幅非相干OFDM信号中心频率为2048MHz，带宽为20.48MHz，2来波方向分别为(84°、100°)，宽带OFDM分为5个窄带，正则化因子λ＝2.7，在SNR＝-6dB情况下对比分辨力。图4是SNR＝-6dB时4种算法分辨力，本发明的方法宽带DOA估计值为：(83,102°)。其他3种方法有许多伪峰且几乎不能分辨2个信号。仿真结果表明，本发明的方法在低信噪比下分辨力好，旁瓣低。

第三，比较本发明的方法与L1-SVD、M-L1-SVD、MG-L1-SVD算法信噪比与角度均方根误差(RMSE)的关系。一个OFDM信号中心频率为2048MHz，带宽为20.48MHz，来波方向为117°，信噪比(SNR)范围为[-6,10]dB，采用50次蒙特卡洛仿真。DOA估计采用均方根误差：

其中，N_c是蒙特卡洛仿真次数，信号真实角度是θ，是信号源第n_c次蒙特卡洛的DOA估计。图5是SNR与RMSE的关系，仿真结果表明，随着信噪比的升高，RMSE越来越小。4种算法在高信噪比时宽带DOA精度相差不大。低信噪比时宽带DOA精度：L1-SVD<MG-L1-SVD<M-L1-SVD<本发明的方法。

实施例2

通过户外实验来验证本发明的方法。图6是在北航沙河校区西北角测量时的照片，左边是发射，距离地面1m，右边是接收，距离地面1.5m。

采用RS SMW200A矢量信号源发射一个中心频率为2.38GHz，带宽为20MHz，功率为35dBm，SNR＝0dB或SNR＝10dB的OFDM信号；玻璃钢天线构成相距为6cm的8天线均匀线阵；接收端采用威视锐2个Y590构成8通道同步接收机，IQ采样率为32MHz；信号源在天线阵的117.2°方向，收发端相距21.9m，此时符合远场条件。

计算时采用10000*100个点，正则化参数λ＝2.7，子窄带频段数J＝5。SNR＝0dB时，DOA估计结果如图7(a)和图7(b)所示，DOA估计值分别为118，118，119，118。SNR＝10dB时，DOA估计结果如图8(a)和图8(b)所示，DOA估计值分别为112，117，119，119。结果表明，利用本发明的方法通过优化子空间的正交性确定加权矩阵增加了信号的稀疏性，产生的谱图主瓣尖锐，旁瓣低,进而提高DOA估计精度并抑制伪峰，能明显提高低信噪比下的正确率，证明了本发明方法的有效性。

本发明提出的基于加权l₁范数的稀疏重构DOA估计算法，通过优化子空间的正交性确定加权矩阵，以加权l₁范数作为最小化的目标函数来增加信号的稀疏性，进而提高低信噪比下DOA估计精度并抑制伪峰。以移动通信中常见的OFDM信号作为仿真对象，仿真与实验验证了算法在低信噪比的情况下主瓣尖锐，旁瓣低，能提高分辨力和有效抑制伪峰，并准确稳定地估计出宽带信号波达方向。综上，本发明的方法对工程有指导作用，可用于反无人机与伪基站定位。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于l₁范数稀疏表示的宽带信号DOA估计方法，其特征在于，包括以下步骤：

S1：接收阵列数据并做FFT，构成宽带信号的J个窄带阵列接收模型；

设P个远场宽带信号入射到一个阵元数为M的均匀阵列，P＜M，将整个观测周期内的信号等分成L段，每段包含K个快拍，每段的阵列接收信号进行时域的快速傅里叶变换后获得信号带宽内的J个窄带阵列接收信号为：

X(f_j)＝A(f_j，θ)S(f_j)+N(f_j)，j＝1，2，…，J (1)

其中，X(f_j)＝[x₁(f_j)，x₂(f_j)，…，x_M(f_j)]^T是一个M×L矩阵，f_j是子窄带中心频率；

A(f_j)＝[a(f_j，θ₁)，a(f_j，θ₂)，…，a(f_j，θ_P)]是频率f_j处的方向矢量，是一个M×P矩阵；

d是两阵元间距，θ_p是入射角；

S(f_j)＝[s₁(f_j)，s₂(f_j)，…，s_p(f_j)]^T是频率f_j处P×L稀疏信号矢量；

N(f_j)＝[N₁(f_j)，N₂(f_j)，…，N_M(f_j)]^T是均值为零、方差为的加性高斯白噪声；

为方便描述，将式(1)简化为：

X_j＝A_jS_j+N_j，j＝1，2，…，J (2)

S2：将频域数据分为J个频段，每个频段的数据矩阵奇异值分解降维，并预处理；

式(2)中，如果X_j能恢复出S_j，则能够根据S_j中的非零行的位置确定信源的DOA估计，即l₀范数问题；但l₀范数的优化是个NP-hard问题，因此将l₁范数代替l₀范数近似求解；

其中，是S_j中每一行的行向量的l₂范数组成的列向量，/>λ_j是受噪声影响的正则化参数，为一个常量；||X_j-A_jS_j||_F是将矩阵X_j-A_jS_j按列拉直并求l₂范数的结果，即：

L1-SVD算法先对接收数据进行奇异值分解，得到M×P的降维矩阵X_svj，设X_j的奇异值分解为：

X_j＝UΛV^H (5)其中，U为X_j的左奇异向量，V为X_j的右奇异向量，Λ为包含奇异值的对角矩阵，则降维后的稀疏模型为：

X_svj＝X_jVD_P＝A_jS_jVD_P+N_jVD_P＝A_jS_svj+N_svj (6)

其中，S_svj＝S_jVD_p，N_svj＝N_jVD_p，D_p＝[I_P 0]^H，I_P是一个P×P的单位矩阵，0为P×(L-P)的零矩阵；

对式(6)预处理；令R_x为原信号矢量X_j的协方差矩阵；式(6)转化为：

利用二阶锥规划即SOCP求解式(7)得：

其中，λ为正则化参数，p和q均为变量；

S3：根据ISM算法思想求出权值W；

由于权值是由噪声子空间与信号方向矢量正交得到：

b^H(θ_i)E_n＝0，i＝1，2，…，N (9)

其中，b(θ_i)是方向矢量，(·)^H是共轭转置，E_n是噪声子空间，N是空间[0°，180°]等间隔划分角度的个数；

权值w_i是将阵列方向矢量投影到噪声子空间的值除以投影到信号子空间的值：

其中，E_s是信号子空间，加权矩阵W为：

W＝diag{w_i} (11)

S4：将步骤S3得到的权值带入步骤S2的稀疏求解模型，分别求出每段的谱估计，谱估计平均即宽带信号DOA；

由式(11)和式(7)得到新的求解模型：

转化为二阶锥规划求解：

求出每段的谱估计，谱估计平均即宽带信号DOA。