信号的波达方向(direction of arrival, DOA)估计在阵列信号处理领域占据重要地位，广泛应用于雷达、医学诊疗等领域，近年来学者们进行了大量的研究，取得了非常显著的成果^[1-3]。以MUSIC^[4]和ESPRIT^[5]算法为代表的经典子空间分解类算法，实现简单且分辨率高，并且在此基础上发展出了多种改进算法^[6-8]，但是这类算法运算量大，需要在信噪比高且快拍数多的情况下才能得到较好的估计性能。压缩感知^[9-10]理论是一个新兴起的领域，它能够突破奈奎斯特采样定理的限制，压缩感知理论中的很多重构算法都被用于DOA估计中，其中平滑 ${l_0}$ 范数重构算法^[11-12]、 ${l_1}{\rm{ - SVD}}$ 算法^[13-14]和贪婪算法是广泛应用的方法。平滑 ${l_0}$ 范数重构算法采用最速下降法^[15]和梯度投影的原理，构造一个函数来近似 ${l_0}$ 范数，但是此算法不能保证重构的精确度。 ${l_1}{\rm{ - SVD}}$ 算法利用 ${l_1}$ 范数重构算法对信号进行重构，估计出信号的DOA，能够在低快拍下实现DOA估计；但是该算法重构过程复杂，会导致较高的计算复杂度，并且在求解过程中，需要人工去设置所需参数，参数选取不好将会对最终的结果产生较大的影响。常用的贪婪算法有匹配追踪(matching pursuit, MP)、正交匹配追踪(orthogonal matching pursuit, OMP)^[16]、正则OMP^[17]等，该类算法需要找到测量矩阵中能够提供最大能量的某列原子，然后一直进行迭代，直到满足迭代条件，得到最优解。OMP算法是一种常用的重构算法，研究者们提出了多种改进算法^[18-19]，但是OMP算法运算时间长，DOA估计精度有待进一步提高。

针对稀疏阵列中贪婪算法DOA估计求解问题，本文引入广义正交匹配追踪(generalized orthogonal matching pursuit, GOMP）用于DOA估计，仿真表明直接将GOMP算法用于DOA估计较OMP算法性能有所提高；为了进一步提高DOA估计的精度，引入最优化理论中的最速下降法，提出基于最速下降法的改进GOMP算法。

1 基于稀疏表示的DOA估计模型

假设有一个由 $M$ 个阵元组成的均匀线阵，阵元间距为 $d$ ，有 $K$ 个方向分别为 ${\theta _{{i}}}(i = 1,2, \cdots ,K)$ 的远场窄带信号 ${{{x}}_k}\left( t \right),k \in \left\{ {1,2, \cdots ,K} \right\}$ 。通过网格划分将空间信号稀疏化，整个空间范围等角度均匀划分为 $N$ 份，则所有可能的方向可以表示为 ${\theta _{{i}}} (i = 1,2, \cdots ,N)$ 。将空间信号进行稀疏化，只有 $K$ 个方向上存在真实的目标信号源，其中， $N$ 远远大于实际信号源的数量 $K$ 。则空间中所有的方向集合都与一个导向矢量相对应，导向矢量为 ${{A}} = \left[ {{{a}}\left( {{\theta _1}} \right),{{a}}\left( {{\theta _2}} \right), \cdots ,{{a}}\left( {{\theta _N}} \right)} \right]$ ，导向矢量也称为阵列流型矩阵，则 $t$ 时刻阵列接收的信号为

${{y}}(t) = {{Ax}}(t) + {{n}}(t)$

(1)

${{A}} = \left[ {{{a}}\left( {{\theta _1}} \right),{{a}}\left( {{\theta _2}} \right), \cdots ,{{a}}\left( {{\theta _N}} \right)} \right]$

(2)

${{a}}\left( {{\theta _n}} \right) = {\left[ {1,{{\rm{e}}^{{\rm{j}}{\varphi _1}}}, \cdots ,{{\rm{e}}^{{\rm{j}}(M - 1){\varphi _n}}}} \right]^{\rm{T}}}$

${\varphi _n} = 2{\text{π}}d\sin {\theta _n}/\lambda ,n \in \left\{ {1,2, \cdots ,N} \right\}$

式中： ${{x}}\left( t \right) = {\left[ {{x_1}\left( t \right),{x_2}\left( t \right), \cdots ,{x_N}\left( t \right)} \right]^{\rm{T}}}$ 为空间入射信号矢量，其稀疏度为 $K\left( {K \ll N} \right)$ ， ${{x}}\left( t \right)$ 中只有 $K$ 个非零元素，其余 $N - K$ 个为零元素，非零元素所在的位置即为真实信号源的目标方向； ${{n}}( t ) = {[ {{n_1}( t ),{n_2}( t ), \cdots ,}} {{{\rm{ }}{n_M}( t )} ]^{\rm{T}}}$ 是均值为0，方差为 ${\sigma ^2}$ 的高斯白噪声； ${{y}}\left( t \right) = {\left[ {{y_1}(t),{y_2}(t), \cdots {y_M}(t)} \right]^{\rm{T}}}$ 为阵列接收信号矢量。

稀疏表示的DOA估计就是在已知阵列接收矢量 ${{y}}$ 以及阵列流型矩阵 ${{A}}$ 的条件下，找出 ${{x}}\left( t \right)$ 中所有非零元素所在的位置集合，该位置集合即是 ${\theta _i}(i = 1,2, \cdots ,N)$ 中实际信号源所在的角度。

2 基于改进GOMP算法的DOA估计 2.1 OMP算法

由式(1)、(2)可知， ${{x}}$ 在空间上是 $K$ 阶稀疏的，只有 $K$ 个非零的元素，而且非零元素对应的位置即为 ${\theta _i}(i = 1,2, \cdots ,N)$ 中实际信号源的角度。

为了找出 ${{x}}$ 中非零元素所在的位置，需要确定阵列流型矩阵 ${{A}}$ 中的有效矢量。OMP算法通过贪婪搜索的方式来搜索有效列。在每次迭代过程中，找出 ${{A}}$ 中与阵列接收信号 ${{y}}$ 中剩余部分最相关的一列，剔除该列对于 ${{y}}$ 的贡献，标记该列并更新残差。当完成 $K$ 次迭代后，即可找出参与对 ${{y}}$ 的线性组合中列的标记集合，最后用该标记集合，在集合 ${\theta _i}(i = 1,2, \cdots ,N)$ 中找到信号源所在的真实角度位置。

2.2 GOMP算法

基于OMP算法的DOA估计方法，耗时比较长且估计的精度不高，为了提高DOA估计精度，采用广义正交匹配追踪算法来求解DOA估计。

OMP算法和GOMP算法均是通过贪婪迭代从阵列流型矩阵 ${{A}}$ 中找到能够提供最大能量的某些列原子，代价函数选用 ${{y}}$ 残差最小。2种算法的区别在于进行重构时，OMP算法每次迭代过程中只选择与 ${{y}}$ 中剩余部分最相关的一列，而GOMP算法在每次迭代过程中会选择与 ${{y}}$ 中剩余部分最相关的 $K$ 列，然后在所有选出的可能的坐标集合中进行搜寻，在索引坐标集合中加入使估计能量值最大的坐标。相对于OMP算法，GOMP算法的搜索选择方式能够更好地保证从信号能量相近或者不完全和不精确的阵列接收信号中恢复信号，能更准确地找出信号源所在的真实角度位置。

2.3 最速下降法的相关理论

最速下降法基于最优化理论^[20]，从梯度的角度来求解最优化问题。函数 $f({{x}})$ 的泰勒展开式为

$f({{x}}) = f({{{x}}_k}) + \nabla f{({{{x}}_k})^{\rm{T}}}({{x}} - {{{x}}_k}) + o({{x}})$

式中： $o({{x}})$ 是高阶无穷小量，其在 ${{{x}}_k}$ 的梯度为 ${{{g}}_k} = \nabla f({{{x}}_k})$ 。记 ${{x}} - {{{x}}_k} = {a_k}{{{d}}_k}$ ，其中， ${a_k}$ 是步长，是一个大于零的数值， ${{{d}}_k}$ 为梯度的方向，则泰勒展开式为 $f({{x}}) = f({{{x}}_k}) + {a_k}{{g}}_k^{\rm{T}}{{{d}}_k} + o({{x}})$ 。在 ${{g}}_k^{\rm{T}}{{{d}}_k} < 0$ 时 $f({{x}}) < f({{{x}}_k})$ ，此时的函数值是下降的， ${{{d}}_k}$ 则是下降的方向。

由柯西−施瓦兹不等式可知：

$\left| {{{d}}_k^{\rm{T}}{{{g}}_k}} \right| \leqslant \left| {{{{d}}_k}} \right|\left| {{{{g}}_k}} \right|$

当且仅当 ${{{d}}_k} = {{{g}}_k}$ 时，等号成立， ${{d}}_k^{\rm{T}}{{{g}}_k}$ 最大。所以当且仅当 ${{{d}}_k} = - {{{g}}_k}$ 时， ${{d}}_k^{\rm{T}}{{{g}}_k}$ 最小，以此为搜索方向，此时 $f({{x}})$ 是下降最快的。

最速下降法在每次迭代时选择一个合适的步长，使得目标函数每次减小的程度最大，即每次得到的值和目标函数值之间的差距最小。用 $f({{x}})$ 代表目标函数，迭代步长为

${a_k} = \arg \min f\left( {{{{x}}_k} - a{{{g}}_k}} \right)$

梯度的迭代为

${{{x}}_{k + 1}} = {{{x}}_k} - a{{{g}}_k}$

因此计算步长的公式就是为找到最好的点 ${{{x}}_{k + 1}}$ ，函数 $f({{x}})$ 在这点可以取得最小值，此时 $f\left( {{{{x}}_k} - a{{{g}}_k}} \right)$ 的极小值就是此次迭代的步长。

最速下降法的过程可以表示为在每次迭代中，从 ${{{x}}_k}$ 出发，按照梯度的负方向去寻找最好的步长值，根据寻找到的步长值来确定下一次迭代过程的出发点 ${{{x}}_{k + 1}}$ 。

将函数的类型扩展至二次型，设函数为

$f({{x}}) = \frac{1}{2}{{{x}}^{\rm{T}}}{{Q}}{{x}} - {{{b}}^{\rm{T}}}{{x}}$

式中： ${{Q}}$ 是对称正定阵； $b \in {{\bf{R}}^n}{\text{；}}{{x}} \in {{\bf{R}}^n}$ ；令

${{{g}}_k} = \nabla f({{x}}) = {{Q}}{{x}} - {{b}}$

迭代表示为

${{{x}}_{k + 1}} = {{{x}}_k} - a{{{g}}_k}$

步长为

${a_k} = \arg \min f\left( {{{{x}}_k} - a{{{g}}_k}} \right)$

在 ${{{g}}_k} = 0$ 时迭代过程停止，在不为零时，设 ${\varphi _k}(a) = f({{{x}}_k} - a{{{g}}_k})$ ，由局部最小点的一阶必要条件得到：

${{{\varphi}}'}_k(a) = {({{{x}}_k} - a{{{g}}_k})^{\rm{T}}}{{Q}}( - {{{g}}_k}) - {{{b}}^{\rm{T}}}( - {{{g}}_k})$

当 ${{{\varphi}}'}_k(a) = 0$ 时， $a{{g}}_k^{\rm{T}}{{Q}}{{{g}}_k}{\rm{ = (}}{{x}}_k^{\rm{T}}{{Q}} - {{{b}}^{\rm{T}}}){{{g}}_k}$ ，因为 ${{Q}}$ 是对称正定阵， ${{Q}} = {{{Q}}^{\rm{T}}}$ ， ${{x}}_k^{\rm{T}}{{Q}} - {{{b}}^{\rm{T}}} = {{g}}_k^{\rm{T}}$ ，由此可以得到步长 ${a_k} = \dfrac{{{{g}}_k^{\rm{T}}{{{g}}_k}}}{{{{g}}_k^{\rm{T}}{{Q}}{{{g}}_k}}}$ ， ${a_k}$ 即为最速下降的步长。

2.4 改进GOMP算法

将最速下降法和GOMP相结合得到改进GOMP算法，改进算法在信号重构时不使用计算量大的最小二乘法，而是沿着负梯度方向进行搜索而得到最优解。

改进算法中迭代选择原子及更新索引集的过程和GOMP算法相一致，在恢复信号过程中，GOMP等贪婪算法均使用最小二乘法进行重构，其矩阵公式为 ${{x}} = {({{ A}^{\rm{T}}}{ A})^{ - 1}}{{ A}^{\rm{T}}}{{y}}$ ，由于需要进行矩阵的求逆运算，导致运算量很大，会消耗很长的时间；而最速下降法虽然需要进行多次运算，但都是基本运算，相对而言计算量较小，所要耗费的时间也会更少。因此改进算法使用最速下降法代替最小二乘法来恢复信号。

用 ${{\varPhi}} = \left\{ {{\varphi _k}} \right\}$ 表示测量矩阵，此测量矩阵即为DOA估计模型中的阵列流型矩阵， ${{y}}$ 表示要恢复的信号， ${{{x}}_i}$ 表示第 $i$ 次循环得到的近似信号， ${{{r}}_i}$ 表示残差， ${{{\varLambda}} _i}$ 表示索引集合， ${{\varPsi}} $ 表示选择原子的列向量集合，先得到矩阵的负梯度方向为 ${{{g}}_i} = {{{\varPsi}} _i}^{\rm{T}}{{{r}}_i}$ ，然后计算每次迭代的步长为

$ \begin{array}{c} {a_i} = \dfrac{{\left\langle {{{{g}}_i},{{{g}}_i}} \right\rangle }}{{\left\langle {{{{g}}_i},{{{\varPsi }}_i}^{\rm{T}}{{{\varPsi }}_i}{{{g}}_i}} \right\rangle }} = \dfrac{{\left\langle {{{{\varPsi }}_i}^{\rm{T}}{{{r}}_i},{{{g}}_i}} \right\rangle }}{{\left\langle {{{{g}}_i},{{{\varPsi }}_i}^{\rm{T}}{{{\varPsi }}_i}{{{g}}_i}} \right\rangle }} =\\ \dfrac{{\left\langle {{{{r}}_i},{{{c}}_i}} \right\rangle }}{{\left\langle {{{{\varPsi }}_i}{{{g}}_i},{{{\varPsi }}_i}{{{g}}_i}} \right\rangle }} = \dfrac{{\left\langle {{{{r}}_i},{{{c}}_i}} \right\rangle }}{{\left\| {{{{c}}_i}} \right\|_2^2}} \end{array} $

(3)

式中 ${{{c}}_i} = {{{\varPsi}} _i}{{{g}}_i}$ 。从而得到恢复信号：

${{{x}}_i} = {{{x}}_{i - 1}} + {a_i}{{{g}}_i}$

最后更新残差：

${{{r}}_i} = {{{r}}_{i - 1}} - {a_i}{{{c}}_i}$

(4)

完成循环，判断是否满足迭代停止条件。

改进后的GOMP算法步骤可总结如下：

1)初始化残差 ${{{r}}_0} = {{y}}$ ，初始化索引集 $\varLambda$ 为空集，初始化迭代次数 $i = 1$ ；

2)计算残差和测量矩阵 ${{\varPhi}} $ 的每一列的内积，找出若干个绝对值最大的元素，并将这些值对应的列序号构成集合 ${{P}}$ ；

3)更新索引集 ${\varLambda _i} = {\varLambda _{i - 1}} \cup P$ ，更新列向量集合 ${{{\varPsi}} _i} = {{{\varPsi}} _{i - 1}} \cup {{{\varphi }} _P}$ ；

4)利用 ${{{g}}_i} = {{{\varPsi}} _i}^{\rm{T}}{{{r}}_i}$ 计算负梯度方向 ${{{g}}_i}$ ，利用式(3)计算步长 ${a_i}$ ；

5)利用式(4)更新残差；

6)若满足迭代停止条件，根据索引集 $\varLambda$ 在角度 ${\theta _{\rm{i}}}(i = 1,2, \cdots ,N)$ 中找到信号源所在的真实角度位置；若不满足，则跳转到步骤2)进行下一次循环。

3 仿真结果分析

本文通过MATLAB仿真验证所提算法的可行性，并与 ${l_1}{\rm{ - SVD}}$ 算法、OMP算法、GOMP算法进行实验对比验证其性能。在实验仿真中，采用均匀直线阵列，噪声为高斯白噪声，均匀线阵的阵元间距为半波长，按照等角度划分方式将空间信号等角度划分为181份，每份角度为1°。均方根误差表示为

${\rm{RMSE}} = \sqrt {\frac{1}{{MK}}{{\sum\limits_{m = 1}^M {\sum\limits_{k = 1}^K {\left( {{{\hat \theta }_{km}} - {\theta _{km}}} \right)^2} } }}} $

(5)

式中： $M$ 表示蒙特卡洛实验次数； $K$ 表示信号源的总个数； ${\theta _{km}}$ 表示第 $k$ 个空间信号源的DOA真实值； ${\hat \theta _{km}}$ 表示第 $k$ 个空间信号源的DOA估计值。

实验1　 残差二范数对比实验。设2个信号源以入射角10°、30°入射到均匀直线阵列上，对结合最速下降法的改进GOMP算法和原始GOMP算法进行性能对比，在单快拍条件下，比较2种情况下的残差二范数，进行500次蒙特卡洛实验，仿真结果如图1所示。

通过图1可以看出，改进GOMP算法比GOMP算法的残差值小，随着信噪比的提升，改进算法的残差值下降更快，可见对GOMP算法的优化是有效的。

实验2　不同信噪比下DOA估计的性能对比实验。设2个信号源以入射角10°、30°入射到均匀直线阵列上，均匀线阵阵元数为12，快拍数为8，信噪比以2 dB为步进从−8 ~8 dB变化，在相同的实验条件下进行500次蒙特卡洛实验。当估计偏差小于或等于1°时，认为估计成功，在不同信噪比条件下，统计各个条件下的成功概率和均方根误差，仿真结果如图2、3所示。

在不同信噪比条件下，由图2可以看出，与其他3种算法相比，改进GOMP算法具有更高的估计成功概率，在信噪比为−2 dB时估计成功概率接近100%。由图3可以看出，在高信噪比情况下，4种算法均趋于稳定，改进GOMP算法具有更小的均方根误差，特别是在低信噪比情况下，改进GOMP算法具有明显的优势。

实验3　不同快拍数下的DOA估计的性能对比实验。设2个信号源以入射角10°、30°入射到均匀直线阵列上，均匀线阵阵元数为12，信噪比为−4 dB，在相同的实验条件下进行500次蒙特卡洛实验。当估计偏差小于或等于1°时，认为估计成功，在快拍数不同的情况下，统计各个条件下的成功概率和均方根误差，仿真结果如图4、5所示。

在不同快拍数情况下，由图4可以看出，改进GOMP算法具有更高的估计成功概率，且具有很高的稳定性。在快拍数为12时，已经达到以100%的概率估计出信号的DOA。随着快拍数的增加，其他3种算法的成功率缓慢增加至趋于平稳。由图5可以看出，4种算法的均方根误差均随着快拍数的增加逐渐减小，改进GOMP算法一直保持着较小的均方根误差，特别是在低快拍数情况下，改进GOMP算法具有明显的优势。

实验4　不同快拍数情况下算法运行时间对比实验。由于l₁-SVD算法运行时间相对于贪婪类算法运行时间长，所以只比较另外3种算法的运行时间。由图6可以看出，随着快拍数的增加，3种算法的运行时间均逐步上升，且改进GOMP算法相较于其他2种算法的运行时间相对较短；在相同的条件下，同其他2种算法相比，改进GOMP算法具有更短的运算时间。因此从整体上看改进GOMP算法的运行时间相对偏低，相比于其他2种算法具有运行时间上的优势。

4 结论

本文提出了一种新的贪婪DOA估计方法，该方法可有效提高DOA估计的精度。GOMP是一种改进的贪婪恢复算法，其利用更好的原子选择方式，提高了贪婪算法的稳定性和角度估计精度。

1)利用GOMP算法求解DOA问题，相较于原有贪婪DOA估计方法提高了估计精度和成功率。

2)引入最速下降法对GOMP算法进行优化，并通过仿真验证了改进方法具有更高的估计精度和成功率，降低了算法的运行时间，相对于原有贪婪DOA估计方法的性能有一定的提升。