CN103971029B - 一种用于网格失配下的doa估计的交替迭代方法 - Google Patents

Abstract

本发明提供一种高精度远场窄带DOA估计方法。首先，在波达方向在空域上具有稀疏性的基础上，把协方差矩阵写成稀疏表示的形式。然后，利用一阶泰勒展开把协方差矩阵改写成网格匹配下的稀疏表示模型。最后，通过交替迭代的方法求解出稀疏的空间功率谱和角度修正值。本发明方法利用一个凸优化问题和最小二乘问题之间的交替更新，分别求解两个联合稀疏的向量，提高了算法的鲁棒性，可以在粗糙的网格上到达高精度的DOA估计性能。

Description

一种用于网格失配下的DOA估计的交替迭代方法

技术领域

本发明属于阵列信号处理领域，主要涉及远场窄带DOA估计。

背景技术

波达方向(DOA)估计一直是阵列信号处理中一个重要的研究领域，它在雷达、声纳、无线通信及电子对抗和侦查等领域中都有着广泛的应用。如何快速地，高精度地实现DOA估计一直是阵列信号处理不断研究和努力的方向。其中经典的算法有：多重信号分类(Multiple Signal Classification,MUSIC)算法、旋转不变子空间(Estimation ofSignal Parameters via Rotational Invariance Technique,ESPRIT)算法等子空间类算法和最大似然估计类算法(Maximum Likelihood,ML)等。然而，基于子空间理论的DOA估计方法虽然实现了超分辨侧向，但是一旦阵列快拍数不足或者出现相干信号源时，这类方法不能有效地区分信号子空间和噪声子空间，其性能会急剧下降。而最大似然估计类算法由于要进行复杂的多维搜索而不具有实用性。

近年来，利用信号在空域的稀疏性，许多基于稀疏表示的DOA估计方法被提出。最具代表性的为l₁-SVD算法，它利用l₁范数来重构稀疏信号，并且在多快拍条件下通过奇异值分解(Singular Value Decomposition,SVD)来减小数据矩阵的规模以及降低噪声的影响。然而，这些基于稀疏表示的方法通常假设波达方向角位于离散化网格上，一旦真实的波达方向角不在预先设定的网格上，即在网格失配的情况下，其估计性能会急剧恶化。另一方面，虽然更密集的网格理论上可以减小重构误差，但是太过密集的离散化网格会使得过完备化字典原子间高度相关。针对这种情况，Malioutov等人介绍了一种稀疏总体最小二乘(Sparse Total Least Squares,TLS)方法，但是其重构性能并不理想。还有一种模型不确定下的稀疏普拟合(Sparse Spectral Fitting with Modeling Uncertainty,SSFMU)方法，利用混合范数来重构稀疏信号，然而，由于该方法要求解一个联合稀疏的优化问题，计算复杂度太高。另一种DOA估计的思路是，利用稀疏贝叶斯推论(Off-Grid SparseBayesian Inference,OGSBI)来重构稀疏信号，但是这种方法的缺点是估计结果容易受初值的影响。

发明内容

本发明的目的在于提供一种用于网格失配下的DOA估计的交替迭代方法。针对网络失配的情况，基于协方差矩阵的稀疏表示，通过交替迭代的方法分别求解两个联合稀疏的空间谱和角度修正值，提高了估计精度。

本发明的思路是：本发明基于协方差矩阵稀疏表示的模型，针对网格失配的情况，首先，在波达方向在空域上具有稀疏性的基础上，把协方差矩阵写成稀疏表示的形式。然后，利用一阶泰勒展开把协方差矩阵改写成网格匹配下的稀疏表示模型。最后，通过交替迭代的方法求解出稀疏的空间功率谱和角度修正值。

本发明的目的通过如下步骤实现：

S1、由阵列接收的K个信号源的数据得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与该信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流形矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，空间协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)，r_k为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H和(·)^T分别表示矩阵的转置和共轭转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述噪声功率σ²进行估计，噪声功率估计值其中，λ_i为协方差矩阵R中M-K个最小的特征，下面的步骤均用代替σ²；

S3、对S1所述空间协方差矩阵R进行量化操作，并写成稀疏表示的网络匹配模型，具体如下：

S31、把S1所述空间协方差矩阵R依次按列排列，写成向量的形式其中，g(θ_k)＝vec(a(θ_k)a^H(θ_k))，vec(·)表示向量化操作，r_s即为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量，表示维数为M²×K的复矩阵；

S32、把角度在[-90°,90°)的空间范围上过完备化为一个离散的网格则把向量r写成稀疏表示的形式，其中，N＞＞K， 是一个K稀疏向量，其K个非零元素分别等于r₁,r₂,...,r_K；

S4、因为即网格失配，真实的波达方向角不在预先设定好的离散化网格上，估计精度受网格间距的限制，需要对真实的导向矢量进行逼近，以减小空间协方差矩阵向量化的模型误差，具体为：

S41、利用一阶泰勒展开对真实的导向矢量进行逼近 其中，θ_k为真实的波达方向角，为网格上离θ_k最近的点，为导向矢量的一阶导数；

S42、记则其中， δ＝[δ₁,δ₂,...,δ_N]^T，则δ与联合稀疏，即非零元素的位置相同；

S5、通过交替迭代的方法得到稀疏的空间谱和修正向量δ，具体如下：

S51、根据S4所述更新后的空间协方差矩阵向量化模型r，求解优化式其中，符号≥对向量的每个元素操作，0为全零的列向量，β为正则化参数；

S52、假设初始修正向量δ⁽⁰⁾＝0，求解当前稀疏空间谱 其中，上标^(j)表示迭代次数；

S53、根据S52所述当前空间谱找出的最大K个元素在网格上对应的位置，记支撑集Λ，所述最大K个元素即为

S54、根据S53所述支撑集Λ构造矩阵和即矩阵G和在支撑集Λ上对应的列；

S55、根据S42所述则 其中，矩阵为列满秩，矩阵diag(r_s)为满秩的；

S56、则当前支撑集上的修正值其中，[·]^-1和分别代表矩阵的逆和伪逆；

S57、将S56所述在网格上稀疏化，得到当前稀疏修正向量δ^(j+1)，所述当前稀疏修正向量δ^(j+1)在支撑集Λ上对应的元素为

S58、将S57所述δ^(j+1)代入S51所述中，依次循环，当满足迭代停止条件或者达到最大迭代次数时，即可得到稀疏的空间谱和修正向量δ；

S6、通过S5所述稀疏的空间谱和修正向量δ，计算得到修正后的波达方向估计值：设支撑集Λ的元素在网格上的索引值为i₁,...,i_K，则修正后波达方向角的估计值为其中，和分别表示网格上索引值为i_k的对应元素，表示向量δ中索引值为i_k的对应元素。

进一步地，S51所述β＝0.5。

进一步地，S58所述τ＝10^-6。

本发明的有益效果是：

本发明利用一个凸优化问题和最小二乘问题之间的交替更新，分别求解两个联合稀疏的向量，可在粗糙的网格上进行波达方向的精确估计，避免了密集的网格带来的高计算量，提高了估计精度，且交替迭代的方法提高了算法的鲁棒性。

附图说明

图1是本发明方法的流程图。

图2是远场窄带信号接收阵列模型图。

图3是本发明方法的残差和误差随迭代次数变化曲线图。

图4是本发明方法与其他方法DOA估计的均方根误差随信噪比变化曲线图。

图5是本发明方法与其他方法DOA估计的均方根误差随快拍数变化曲线图。

图6是本发明方法与其他方法DOA估计的均方根误差随网格间距变化曲线图。

具体实施方式

下面结合实施例和附图，详细说明本发明的技术方案。

图1是本发明交替迭代方法用于网格失配下的DOA估计的一种具体实施方式流程图。如图1所示，本发明交替迭代方法用于网格失配下的DOA估计包括以下步骤：

S1、由阵列接收的K个信号源的数据得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与该信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流形矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，空间协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)，r_k为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H和(·)^T分别表示矩阵的转置和共轭转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述噪声功率σ²进行估计，噪声功率估计值其中，λ_i为协方差矩阵R中M-K个最小的特征，下面的步骤均用代替σ²；

S3、对S1所述空间协方差矩阵R进行量化操作，并写成稀疏表示的网络匹配模型，具体如下：

S31、把S1所述空间协方差矩阵R依次按列排列，写成向量的形式其中，g(θ_k)＝vec(a(θ_k)a^H(θ_k))，vec(·)表示向量化操作，r_s即为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量，表示维数为M²×K的复矩阵；

S32、把角度在[-90°,90°)的空间范围上过完备化为一个离散的网格则把向量r写成稀疏表示的形式，其中，N＞＞K， 是一个K稀疏向量，其K个非零元素分别等于r₁,r₂,...,r_K；

S4、因为即网格失配，真实的波达方向角不在预先设定好的离散化网格上，估计精度受网格间距的限制，需要对真实的导向矢量进行逼近，以减小空间协方差矩阵向量化的模型误差，具体为：

S41、利用一阶泰勒展开对真实的导向矢量进行逼近 其中，θ_k为真实的波达方向角，为网格上离θ_k最近的点，为导向矢量的一阶导数；

S42、记则其中， δ＝[δ₁,δ₂,...,δ_N]^T，则δ与联合稀疏，即非零元素的位置相同；

S5、通过交替迭代的方法得到稀疏的空间谱和修正向量δ，具体如下：

S51、根据S4所述更新后的空间协方差矩阵向量化模型r，求解优化式其中，符号≥对向量的每个元素操作，0为全零的列向量，β为正则化参数，β＝0.5；

S52、假设初始修正向量δ⁽⁰⁾＝0，求解当前稀疏空间谱 其中，上标^(j)表示迭代次数；

S53、根据S52所述当前空间谱找出的最大K个元素在网格上对应的位置，记支撑集Λ，所述最大K个元素即为

S54、根据S53所述支撑集Λ构造矩阵和即矩阵和在支撑集Λ上对应的列；

S55、根据S42所述则 其中，矩阵为列满秩，矩阵diag(r_s)为满秩的；

S56、则当前支撑集上的修正值其中，[·]^-1和分别代表矩阵的逆和伪逆；

S57、将S56所述在网格上稀疏化，得到当前稀疏修正向量δ^(j+1)，所述当前稀疏修正向量δ^(j+1)在支撑集Λ上对应的元素为

S58、将S57所述δ^(j+1)代入S51所述中，依次循环，当满足迭代停止条件或者达到最大迭代次数时，即可得到稀疏的空间谱和修正向量δ，τ＝10^-6；

S6、通过S5所述稀疏的空间谱和修正向量δ，计算得到修正后的波达方向估计值：设支撑集Λ的元素在网格上的索引值为i₁,...,i_K，则修正后波达方向角的估计值为其中，和分别表示网格上索引值为i_k的对应元素，表示向量δ中索引值为i_k的对应元素。

实施例中用均方根误差(RMSE)来评估各算法的性能，其定义为：其中，Mon为蒙特卡洛实验次数，和θ_k分别代表第m次蒙特卡洛实验估计得到的第k个角度和第k个真实角度。

实施例1

本发明估计方法的残差和误差曲线随迭代次数变化仿真：

实施例1采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°,36.3°]入射到阵列。为了使入射方向角不落在网格上，取离散化网格为{-90°,-88°,...,88°°，间隔2°。参考信噪比固定为10dB，采样快拍数为250。分别定义残差和误差其中，为第j次迭代得到的稀疏空间谱，为真实的稀疏空间谱。单次试验观察其曲线随迭代次数的变化。

实施例1中DOA估计方法包括以下步骤：

根据阵列接收信号x(t)得到协方差矩阵R；

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2；

通过交替迭代的方法求解出稀疏空间谱和修正向量δ；

通过和δ，计算得到修正后的波达方向估计值：并计算均方根误差。

按照本发明的方法估计得到的残差和误差随迭代次数变化的曲线如图3所示。图3可以看到，利用本发明的估计方法，残差曲线和误差曲线都随着迭代次数的增加而下降，且经过4次迭代后就会基本收敛到一个固定值，且该固定的值逼近真实的空间谱，说明了本发明的估计方法可以收敛到一个最优解。

实施例2

本发明估计值的均方根误差随信噪比变化仿真：

实施例2采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°,36.3°]入射到阵列。为了使入射方向角不落在网格上，取离散化网格为{-90°,-88°,...,88°}，间隔2°。采样快拍数为250。参考信噪比SNR从-4dB到20dB变化，间隔为4dB，每个信噪比进行1000次蒙特卡洛实验。

实施例2中DOA估计方法包括以下步骤：

根据不同信噪比下的阵列接收信号x(t)得到协方差矩阵R。

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2。

通过交替迭代的方法求解出稀疏空间谱和修正向量δ。

通过和δ，计算得到修正后的波达方向估计值并计算均方根误差。

按照本发明的方法估计得到的波达方向角的均方根误差随信噪比变化的曲线如图4所示。图4可以看到，利用本发明的估计方法，即使在粗糙的网格且信噪比为0dB时都能达到0.2°以内的估计精度，当信噪比增加到20dB时，均方根误差可以达到0.05°。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

实施例3

本发明估计值的均方根误差随快拍数变化仿真：

实施例3采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°°,36.3°]入射到阵列。为了使入射方向角不落在网格上，取离散化网格为{-90°,-88°,...，88°}，间隔2°。参考信噪比SNR固定为10dB。快拍数从100到400，间隔50，每个快拍数进行1000次蒙特卡洛实验。

实施例3中DOA估计方法包括以下步骤：

根据不同快拍数下的阵列接收信号x(t)得到协方差矩阵R。

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2。

通过交替迭代的方法求解出稀疏空间谱和修正向量δ。

通过和δ，计算得到修正后的波达方向估计值：并计算均方根误差。

按照本发明的方法估计得到的波达方向角的均方根误差随快拍数变化的曲线如图5所示。图5可以看到，利用本发明的估计方法，估计性能随快拍数的增加明显提高，即使在粗糙的网格且快拍数为100时都能达到0.1°以内的估计精度。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

实施例4

本发明估计值的均方根误差随网格间距变化仿真：

实施例4采用的接收阵列为如附图2所示的由8个阵元组成的半波长均匀线性阵列。参考阵元为编号1的阵元天线。两个相同功率的信号源按入射方向[-14.5°°,36.3°]入射到阵列。参考信噪比SNR固定为10dB，快拍数设定为250。网格间距从1°到4°，间隔1°，每个网格间距进行1000次蒙特卡洛实验。

实施例4中DOA估计方法包括以下步骤：

根据阵列接收信号x(t)得到协方差矩阵R。

对R进行特征值分解得到噪声功率的估计值然后对R进行向量化得到r，其中M＝8，K＝2。

在不同的网格间距下，通过交替迭代的方法求解出稀疏空间谱和修正向量δ。

通过和δ，计算得到修正后的波达方向估计值：并计算均方根误差。

按照本发明的方法估计得到的波达方向角的均方根误差随网格间距变化的曲线如图6所示。图6可以看到，利用本发明的估计方法，均方根误差会随网格间距的增加而增加。当网格间距在3°以内时均方根误差都在0.05°左右，一旦网格间距达到4°，均方根误差急剧增加，但是还是在0.16°左右。本发明的估计方法与其他估计方法相比，估计性能有明显的提高，说明了本发明的估计方法是有效的。

Claims (2)

1.一种用于网络失配下的DOA估计的交替迭代方法，其特征在于，包括如下步骤：

S1、由阵列接收的K个信号源的数据得到空间协方差矩阵R＝E[x(t)x^H(t)]＝A(θ)R_sA^H(θ)+σ²I_M，其中，x(t)＝[x₁(t),x₂(t),...,x_M(t)]^T表示各个阵元接收信号构成的矩阵，M为阵元数目，K为远场窄带信号源个数，θ_k为第k个信号源入射到阵列的角度，为第k个信号源的导向矢量，为第k个信号源入射到第m个阵元与该信号源入射到参考阵元的相位差，λ为入射信号的波长，d为相邻两个阵元的间距，A(θ)＝[a(θ₁),a(θ₂),...,a(θ_K)]为阵列流形矩阵，s(t)＝[s₁(t),s₂(t),...,s_K(t)]^T为入射信号，附加噪声n(t)为与各个信号源不相关的加性零均值高斯白噪声，空间协方差矩阵R中，R_s＝diag(r₁,r₂,...,r_K)，r_k为入射信号的功率，σ²为噪声功率，I_M为M阶的单位矩阵，E[·]表示期望，(·)^H和(·)^T分别表示矩阵的共轭转置和转置，k＝1,2,...,K，m＝1,2,...,M；

S2、对S1所述噪声功率σ²进行估计，噪声功率估计值其中，λ_i为协方差矩阵R中K个最小的特征，下面的步骤均用代替σ²；

S3、对S1所述空间协方差矩阵R进行量化操作，并写成稀疏表示的网络匹配模型，具体如下：

S31、把S1所述空间协方差矩阵R依次按列排列，写成向量的形式其中，g(θ_k)＝vec(a(θ_k)a^H(θ_k))，vec(·)表示向量化操作，r_s即为矩阵R_s对角线元素构成的向量，I_v为单位矩阵I_M按列排列得到的向量，表示维数为M²×K的复矩阵；

S32、把角度在[-90°,90°)的空间范围上过完备化为一个离散的网格则把向量r写成稀疏表示的形式，其中，N＞＞K， 是一个K稀疏向量，其K个非零元素分别等于r₁,r₂,…,r_K；

S4、因为即网格失配，真实的波达方向角不在预先设定好的离散化网格上，估计精度受网格间距的限制，需要对真实的导向矢量进行逼近，以减小空间协方差矩阵向量化的模型误差，具体为：

S41、利用一阶泰勒展开对真实的导向矢量进行逼近 其中，θ_k为真实的波达方向角，为网格上离θ_k最近的点，为导向矢量的一阶导数；

S42、记则其中， δ＝[δ₁,δ₂,…,δ_N]^T，则δ与联合稀疏，即非零元素的位置相同；

S5、通过交替迭代的方法得到稀疏的空间谱和修正向量δ，具体如下：

S51、根据S4所述更新后的空间协方差矩阵向量化模型r，求解优化式 其中，符号≥对向量的每个元素操作，0为全零的列向量，β为正则化参数；

S52、假设初始修正向量δ⁽⁰⁾＝0，求解当前稀疏空间谱 其中，上标^(j)表示迭代次数；

S53、根据S52所述当前空间谱找出的最大K个元素在网格上对应的位置，记支撑集Λ，所述最大K个元素即为

S54、根据S53所述支撑集Λ构造矩阵和即矩阵和在支撑集Λ上对应的列；

S55、根据S42所述则 其中，矩阵为列满秩，矩阵diag(r_s)为满秩的；

S56、则当前支撑集上的修正值其中，[·]^-1和分别代表矩阵的逆和伪逆；

S57、将S56所述在网格上稀疏化，得到当前稀疏修正向量δ^(j+1)，所述当前稀疏修正向量δ^(j+1)在支撑集Λ上对应的元素为

S58、将S57所述δ^(j+1)代入S51所述中，依次循环，当满足迭代停止条件或者达到最大迭代次数时，即可得到稀疏的空间谱和修正向量δ，其中，τ＝10^-6；

S6、通过S5所述稀疏的空间谱和修正向量δ，计算得到修正后的波达方向估计值：设支撑集Λ的元素在网格上的索引值为i₁,...,i_K，则修正后波达方向角的估计值为其中，和分别表示网格上索引值为i_k的对应元素，表示向量δ中索引值为i_k的对应元素。

2.根据权利要求1所述的一种用于网络失配下的DOA估计的交替迭代方法，其特征在于：S51所述β＝0.5。