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

摘要

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

说明书

技术领域

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

背景技术

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

近年来，利用信号在空域的稀疏性，许多基于稀疏表示的DOA估计方法被提 出。最具代表性的为l₁-SVD算法，它利用l₁范数来重构稀疏信号，并且在多快拍 条件下通过奇异值分解(Singular>

发明内容

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

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

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

S1、由阵列接收的K个信号源的数据$x\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}a\left({\theta }_k\right)s_k\left(t\right)+n\left(t\right)=A\left(\theta \right)s\left(t\right)+n\left(t\right),$得 到空间协方差矩阵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依次按列排列，写成向量的形式 $\left(\begin{array}{c}r=\mathrm{vec}\left(R\right)\\ =\mathrm{vec}(A\left(\theta \right)R_s{A}^H\left(\theta \right)+{\hat{\sigma }}^2{I}_M)\\ =G\left(\theta \right)r_s+{\hat{\sigma }}^2{I}_v\end{array}\right),$其中，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、利用一阶泰勒展开对真实的导向矢量进行逼近 $a\left(\sin {\theta }_k\right)\approx a\left(\sin {\tilde{\theta }}_i\right)+a^{\prime }\left(\sin {\tilde{\theta }}_i\right)(\sin {\tilde{\theta }}_i-\sin {\theta }_k),$其中，θ_k为真实的波达方向角，为网格上 离θ_k最近的点，为导向矢量的一阶导数；

S42、记${\delta }_i=\sin {\tilde{\theta }}_i-\sin {\theta }_k,$则$\left(\begin{array}{c}r=\underset{i=1}{\overset{N}{\Sigma }}{\tilde{r}}_i[g\left({\tilde{\theta }}_i\right)+{\delta }_{i}b\left({\tilde{\theta }}_i\right)]+{\hat{\sigma }}^2{I}_v\\ =(G\left(\tilde{\theta }\right)+B\left(\tilde{\theta }\right)\mathrm{diag}\left(\delta \right)){\tilde{r}}_s+{\hat{\sigma }}^2{I}_v\\ =G\left(\tilde{\theta }\right){\tilde{r}}_s+B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta +{\hat{\sigma }}^2{I}_v\end{array}\right),$其中， $b\left({\tilde{\theta }}_i\right)=g^{\prime }\left({\tilde{\theta }}_i\right).$δ＝[δ₁,δ₂,...,δ_N]^T，则δ与联合稀 疏，即非零元素的位置相同；

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

S51、根据S4所述更新后的空间协方差矩阵向量化模型r，求解优化式 $({\tilde{r}}_s,\delta )=\arg \underset{{\tilde{r}}_s\ge 0,\delta }{\min }{\left|\right|{\tilde{r}}_s\left|\right|}_1+\beta {\left|\right|r-{\hat{\sigma }}^2{I}_v-G\left(\tilde{\theta }\right){\tilde{r}}_s-B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta \left|\right|}_2,$其中，符号≥对向量的每个元 素操作，0为全零的列向量，β为正则化参数；

S52、假设初始修正向量δ⁽⁰⁾＝0，求解当前稀疏空间谱 ${\tilde{r}}_s^{(j+1)}=\arg \underset{{\tilde{r}}_s\ge 0}{\min }{\left|\right|{\tilde{r}}_s\left|\right|}_1+\beta {\left|\right|r-{\hat{\sigma }}^2{I}_v-G\left(\tilde{\theta }\right){\tilde{r}}_s-B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta \left|\right|}_2,$其中，上标^(j)表示迭代次数；

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

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

S55、根据S42所述$r=G\left(\tilde{\theta }\right){\tilde{r}}_s+B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta +{\hat{\sigma }}^2{I}_v,$则 $r=G_{\Lambda }\left(\tilde{\theta }\right)r_s+B_{\Lambda }\left(\tilde{\theta }\right)\mathrm{diag}\left(r_s\right)\delta +{\hat{\sigma }}^2{I}_v,$其中，矩阵为列满秩，矩阵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个信号源的数据$x\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}a\left({\theta }_k\right)s_k\left(t\right)+n\left(t\right)=A\left(\theta \right)s\left(t\right)+n\left(t\right),$得 到空间协方差矩阵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依次按列排列，写成向量的形式 $\left(\begin{array}{c}r=\mathrm{vec}\left(R\right)\\ =\mathrm{vec}(A\left(\theta \right)R_s{A}^H\left(\theta \right)+{\hat{\sigma }}^2{I}_M)\\ =G\left(\theta \right)r_s+{\hat{\sigma }}^2{I}_v\end{array}\right),$其中，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、利用一阶泰勒展开对真实的导向矢量进行逼近 $a\left(\sin {\theta }_k\right)\approx a\left(\sin {\tilde{\theta }}_i\right)+a^{\prime }\left(\sin {\tilde{\theta }}_i\right)(\sin {\tilde{\theta }}_i-\sin {\theta }_k),$其中，θ_k为真实的波达方向角，为网格上 离θ_k最近的点，为导向矢量的一阶导数；

S42、记${\delta }_i=\sin {\tilde{\theta }}_i-\sin {\theta }_k,$则$\left(\begin{array}{c}r=\underset{i=1}{\overset{N}{\Sigma }}{\tilde{r}}_i[g\left({\tilde{\theta }}_i\right)+{\delta }_{i}b\left({\tilde{\theta }}_i\right)]+{\hat{\sigma }}^2{I}_v\\ =(G\left(\tilde{\theta }\right)+B\left(\tilde{\theta }\right)\mathrm{diag}\left(\delta \right)){\tilde{r}}_s+{\hat{\sigma }}^2{I}_v\\ =G\left(\tilde{\theta }\right){\tilde{r}}_s+B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta +{\hat{\sigma }}^2{I}_v\end{array}\right),$其中， $b\left({\tilde{\theta }}_i\right)=g^{\prime }\left({\tilde{\theta }}_i\right).$δ＝[δ₁,δ₂,...,δ_N]^T，则δ与联合稀 疏，即非零元素的位置相同；

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

S51、根据S4所述更新后的空间协方差矩阵向量化模型r，求解优化式 $({\tilde{r}}_s,\delta )=\arg \underset{{\tilde{r}}_s\ge 0,\delta }{\min }{\left|\right|{\tilde{r}}_s\left|\right|}_1+\beta {\left|\right|r-{\hat{\sigma }}^2{I}_v-G\left(\tilde{\theta }\right){\tilde{r}}_s-B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta \left|\right|}_2,$其中，符号≥对向量的每个元 素操作，0为全零的列向量，β为正则化参数，β＝0.5；

S52、假设初始修正向量δ⁽⁰⁾＝0，求解当前稀疏空间谱 ${\tilde{r}}_s^{(j+1)}=\arg \underset{{\tilde{r}}_s\ge 0}{\min }{\left|\right|{\tilde{r}}_s\left|\right|}_1+\beta {\left|\right|r-{\hat{\sigma }}^2{I}_v-G\left(\tilde{\theta }\right){\tilde{r}}_s-B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta \left|\right|}_2,$其中，上标^(j)表示迭代次数；

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

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

S55、根据S42所述$r=G\left(\tilde{\theta }\right){\tilde{r}}_s+B\left(\tilde{\theta }\right)\mathrm{diag}\left({\tilde{r}}_s\right)\delta +{\hat{\sigma }}^2{I}_v,$则 $r=G_{\Lambda }\left(\tilde{\theta }\right)r_s+B_{\Lambda }\left(\tilde{\theta }\right)\mathrm{diag}\left(r_s\right)\delta +{\hat{\sigma }}^2{I}_v,$其中，矩阵为列满秩，矩阵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°左右。本发明的估计方法与其他估计方法相比，估计性能有明显的 提高，说明了本发明的估计方法是有效的。