CN105334488A - 基于信源数估计的栅格偏移优化目标到达角估计方法 - Google Patents

Abstract

本发明公开了一种基于信源数估计的栅格偏移优化目标到达角估计方法，主要解决现有技术在信源数未知、小样本和相干源情况下到达角估计性能降低的问题。其实现过程为：利用阵元去噪对接收数据进行去噪处理获得数据观测矢量；对计算获得的恢复矢量求峰值进而估计信源数和第一次目标到达角；计算当前峰值和当前代价函数；通过计算并比较左移代价函数值和右移代价函数值来更新当前栅格参数；直到可调栅格步长满足栅格精度且动态信源数大于所估计的信源数时停止角度估计，获得第二次目标到达角估计矢量。本发明实现了在信源数未知、小样本和相干源情况下获得高精度稀疏多目标的到达角估计。

Description

基于信源数估计的栅格偏移优化目标到达角估计方法

技术领域

本发明属于信号处理技术领域，更进一步涉及阵列信号处理和压缩感知技术领域中的一种基于信源数估计的栅格偏移优化目标到达角估计方法。本发明可以有效解决在信源数未知、小样本和相干源情况下到达角估计性能降低的问题，实现稀疏多目标的高精度到达角估计。

背景技术

到达角(DOA：DirectionofArrival)估计是阵列信号处理领域所研究的主要问题之一，随着科技的不断进步，对目标到达角的估计精确度和分辨率要求越来越高。DOA估计技术是在波束形成和时域谱的基础上发展起来的，在无线电通信、雷达、电子侦察等领域有着广泛的应用。

目前，已经提出了很多目标到达角的估计方法，比如最典型的MUSIC(MultipleSignalIdentificationandClassification)算法，该类算法的出现使得DOA估计突破了天线孔径瑞利限，实现了波束内多目标DOA超分辨。但MUSIC方法的DOA估计效果依赖于协方差矩阵的估计准确性，理论上协方差矩阵估计所需样本数大于阵元数的2倍，协方差矩阵估计误差小于3dB，实际中所需样本数通常要求大于阵元数的8倍以上。然而，工程上由于阵列的高动态或目标信号的快运动造成平稳样本数不足，导致MUSIC方法的性能下降甚至完全失效。阵列接收的信号除直达波外，还有由于环境的反射产生的多径延迟信号，这些延迟信号与直达波是相干信号。相干信号的存在使信号协方差矩阵发生秩亏，导致协方差矩阵特征分解获得的信号子空间不能准确表征直达波与多径延迟信号对应的导向矢量，利用空间正交性的子空间类到达角估计方法不再有效。

稀疏到达角估计，利用目标在空域的稀疏性，在小样本情况下可以实现相干/非相干信号的到达角高分辨估计，但是该类方法的估计精度取决于字典的栅格，当栅格太密时，造成字典列与列相关性太高，导致稀疏恢复性能下降甚至失效而无法准确估计目标到达角。通常稀疏恢复DOA估计方法需要准确已知信源数，然而实际中信源数往往预先未知，需要寻找能够准确估计信源数的稀疏信号恢复方法。

中国科学院声学研究所提出的专利申请“一种信源数估计方法及其波达方向估计方法”(申请号：200810008058.5，公开号：101272168)中公开了一种波达方向估计中基于特征空间的信源数估计方法。该方法将信号协方差矩阵进行特征值分解，得到特征矢量矩阵，再由特征矢量求信源数判据，进而求比例值来估计信源数。虽然该方法能够准确估计信源数并能够为波达方向估计信号处理过程省去大量的计算量，减少硬件开销，但是仍然存在的不足之处是，该方法在小样本和相干源的情况下由于协方差矩阵估计不准，则会导致信源数和波达方向估计性能下降。

李磊，李国林，刘润杰在其发表的论文“基于相干积累矩阵重构的波达方向估计新方法”(《雷达学报》2015,4(2)：178-184)中公开了一种基于相干积累矩阵重构的快速解相干方法。该方法主要针对短时小样本条件下的相干信号波达方向估计问题。该方法首先利用相干积累技术对阵列接收快拍进行处理，得到累积快拍矢量，再根据该矢量构造一个非降维等效协方差矩阵，来实现解相干。虽然该方法避免了阵列孔径损失，估计精度高、计算量小，但是仍然存在的不足之处是，该方法需要准确已知信源数，而实际中信源数往往预先未知或难以估准。

桂林电子科技大学所提出的专利申请“一种未知信源数高精度波达方向估计方法”(申请号：200910325063.4，公开号：104076324)中公开了一种未知信源数高精度波达方向估计方法。该方法通过设定K个不同的时间平滑间隔去对天线阵列进行时间平滑处理，以构建K个空时自相关矩阵，再计算最终的组合空时自相关矩阵和波达方向估计的空间谱函数，最后通过逐步改变搜索方向对空间谱函数进行谱峰搜索进而估计信源数和来波方向。该方法虽然可以在未知信源数情况下对目标到达角方向进行估计，但是仍然存在的不足之处是，该方法在小样本和相干源情况下估计性能降低甚至失效。

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种基于信源数估计的栅格偏移优化目标到达角估计方法。本发明对阵列接收数据进行阵元去噪处理来构建数据观测矢量，将信源数进行动态递增来估计信源数，利用硬门限参数化反馈零空间调整获得第一次目标到达角估计矢量和恢复矢量，根据恢复矢量来寻找峰值位置矢量，计算定义为多目标峰值功率和与其它栅格功率和之比的代价函数，比较计算获得的三种代价函数大小进而确定待更新的栅格偏移参数，迭代提升多目标到达角的估计精度。本发明方法无需已知信源数，适用于小样本、相干源情况。为实现上述目的，本发明的方法包括如下步骤：

本发明的具体实施步骤如下：

(1)接收小样本回波数据：

利用多路数据采集器获取天线阵列所接收到的小样本回波数据；

(2)阵元去噪：

按照下式，对天线阵列所接收到的小样本回波数据进行阵元去噪，得到数据观测矢量：

$b=\frac{1}{Q}\underset{q=1}{\overset{Q}{\Σ}}x\left(t_q\right)x_1^{*}\left(t_q\right)$

其中，b表示N×1维的数据观测矢量，N表示天线阵列中的阵元个数，Q表示采样快拍数，q表示采样序号，x(t_q)表示t_q时刻天线阵列所接收到的N×1维小样本回波数据，x₁(t_q)表示t_q时刻天线阵列中第一个天线阵元所接收到的小样本回波数据，*表示共轭操作；

(3)计算粗栅格字典：

(3a)按照下式，计算导向矢量：

$a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{ccccccc}1& e^{\frac{j2\mathrm{\π}(d_2-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& e^{\frac{j2\mathrm{\π}(d_3-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_n-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_N-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}\end{array}\right]}^T$

其中，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量，N表示天线阵列中的阵元个数，θ_i∈Θ，∈表示属于符号，Θ表示角度搜索范围，i表示搜索角度θ_i在角度搜索范围Θ中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N，T表示转置操作；

(3b)按照下式，计算粗栅格字典：

A＝[a(θ₁)a(θ₂)…a(θ_i)…a(θ_L)]

其中，A表示N×L维的粗栅格字典，N表示天线阵列中的阵元个数，L表示总的粗栅格个数，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量；

(4)计算初始恢复矢量：

(4a)按照下式，计算粗栅格字典A的右逆矩阵：

P_a＝A^H(A·A^H)^-1

其中，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典，H表示共轭转置操作，-1表示求逆操作；

(4b)按照下式，计算粗栅格字典A的正交投影矩阵：

P＝I-P_aA

其中，P表示粗栅格字典A的正交投影矩阵，I表示单位矩阵，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典；

(4c)按照下式，计算初始恢复矢量：

${\tilde{r}}^0=P_a\·b$

其中，表示初始恢复矢量，P_a表示粗栅格字典A的右逆矩阵，b表示数据观测矢量；

(5)估计信源数：

(5a)将动态信源数初始化为1；

(5b)按照下式，对第k次内循环中的恢复矢量进行排序操作：

$\[{\tilde{r}}^l,T\]=sort\left(\right|{\tilde{r}}^k\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示恢复矢量取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|,'descend')表示降序排列操作；

(5c)按照下式，计算第k+1次内循环中的恢复矢量：

${\tilde{r}}^{k+1}={\tilde{r}}^k+P(u^k-{\tilde{r}}^k)$

其中，表示第k+1次内循环中的恢复矢量，表示第k次内循环中的恢复矢量，u^k表示第k次内循环中的中间辅助矢量，P表示粗栅格字典A的正交投影矩阵；

(5d)按照下式，计算内循环相对误差值：

$H_2=\frac{\left|\right|{\tilde{r}}^{k+1}-{\tilde{r}}^k|{|}_2}{\left|\right|{\tilde{r}}^k|{|}_2}$

其中，H₂表示内循环相对误差值，和分别表示内循环次数为k+1和k时的恢复矢量，||·||₂表示取2范数操作；

(5e)判断内循环相对误差值H₂是否大于10^‐3，若是，则执行步骤(5b)，否则，执行步骤(5f)；

(5f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；

(5g)按照下式，计算失配相对误差：

${\mathrm{\γ}}^{l+1}=\frac{\left|\right|{\mathrm{Au}}^k-b|{|}_2}{\left|\right|b|{|}_2}$

其中，γ^l+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示粗栅格字典，u^k表示第k次内循环中的中间辅助矢量，b表示数据观测矢量，||·||₂表示取2范数操作；

(5h)按照下式，计算外循环相对误差值：

H₁＝|γ^l+1-γ^l|

其中，H₁表示外循环相对误差值，γ^l+1和γ^l分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；

(5i)判断外循环相对误差值H₁是否大于0.05，若是，则执行步骤(5b)，否则，执行步骤(5j)；

(5j)将外循环结束时的动态信源数值作为信源数的估计值；

(6)第一次估计目标到达角：

(6a)按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(6b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；

(6c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量中；

(7)按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(8)按照下式，计算当前代价函数：

$F^{\left(M\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V\left(s\right)}{\underset{i\∉p_I(1:\hat{P})}{\Σ}\left|\tilde{r}\left(i\right)\right|}$

其中，F^(M)表示当前代价函数，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示恢复矢量中的第i个元素， 表示不属于符号，表示峰值位置矢量p_I的前个元素，表示信源数的估计值，|·|表示取模值操作；

(9)按照下式，计算当前峰值：

$\left\{\begin{array}{c}{p}_V^{\left(M\right)}=p_V\left(s\right)\\ p_I^{\left(M\right)}=p_I\left(s\right)\end{array}\right.$

其中，表示当前峰值，表示当前峰值对应的位置下标，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，p_I(s)表示峰值位置矢量p_I中的第s个元素；

(10)计算左移代价函数：

(10a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；

(10b)按照下式，计算左移栅格字典：

其中，A^(L)表示左移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(L)，∈表示属于符号，Θ^(L)表示左移角度搜索范围，i表示搜索角度θ_i在左移角度搜索范围Θ^(L)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N；

(10c)按照下式，计算左移栅格字典的正交投影矩阵：

P^(L)＝I-(A^(L))^H(A^(L)(A^(L))^H)^-1A^(L)

其中，P^(L)表示左移栅格字典A^(L)的正交投影矩阵，I表示单位矩阵，A^(L)表示左移栅格字典，H表示共轭转置操作，-1表示求逆操作；

(10d)按照下式，计算左移恢复矢量：

${\tilde{r}}^{\left(L\right)}={\tilde{r}}^{\left(L\right)}+P^{\left(L\right)}\left(u^{\left(L\right)}-{\tilde{r}}^{\left(L\right)}\right)$

其中，表示左移恢复矢量，u^(L)表示左移辅助矢量，P^(L)表示左移栅格字典A^(L)的正交投影矩阵；

(10e)按照下式，寻找左移峰值位置矢量：

$\[p_V^{\left(L\right)},p_I^{\left(L\right)}\]=findpeaks\left(\left|{\tilde{r}}^{\left(L\right)}\right|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示左移恢复矢量取模值后按降序重排的左移峰值矢量，表示左移峰值矢量的元素在左移恢复矢量模值中对应的下标组成的左移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(10f)按照下式，计算左移代价函数：

$F^{\left(L\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(L\right)}\left(s\right)}{\underset{i\∉p_I^{\left(L\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(L\right)}\left(i\right)\right|}$

其中，F^(L)表示左移代价函数，表示左移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示左移恢复矢量中的第i个元素， 表示不属于符号，表示左移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(11)计算右移代价函数：

(11a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；

(11b)按照下式，计算右移栅格字典：

其中，A^(R)为右移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(R)，∈表示属于符号，Θ^(R)表示右移角度搜索范围，i表示搜索角度θ_i在右移角度搜索范围Θ^(R)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N；

(11c)按照下式，计算右移栅格字典的正交投影矩阵：

P^(R)＝I-(A^(R))^H(A^(R)(A^(R))^H)^-1A^(R)

其中，P^(R)表示右移栅格字典A^(R)的正交投影矩阵，I表示单位矩阵，A^(R)表示右移栅格字典，H表示共轭转置操作，-1表示求逆操作；

(11d)按照下式，计算右移恢复矢量：

${\tilde{r}}^{\left(R\right)}={\tilde{r}}^{\left(R\right)}+P^{\left(R\right)}\left(u^{\left(R\right)}-{\tilde{r}}^{\left(R\right)}\right)$

其中，表示右移恢复矢量，u^(R)表示右移辅助矢量，P^(R)表示右移栅格字典A^(R)的正交投影矩阵；

(11e)按照下式，寻找右移峰值位置矢量：

$\[p_V^{\left(R\right)},p_I^{\left(R\right)}\]=findpeaks\left(\left|{\tilde{r}}^{\left(R\right)}\right|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示右移恢复矢量取模值后按降序重排的右移峰值矢量，表示右移峰值矢量的元素在右移恢复矢量模值中对应的下标组成的右移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(11f)按照下式，计算右移代价函数：

$F^{\left(R\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(R\right)}\left(s\right)}{\underset{i\∉p_i^{\left(R\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(R\right)}\left(i\right)\right|}$

其中，F^(R)表示右移代价函数，表示右移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示右移恢复矢量中的第i个元素， 表示不属于符号，表示右移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(12)更新当前栅格参数：

(12a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(12b)，否则，执行步骤(12e)；

(12b)将步骤(8)的当前代价函数的值更新为左移代价函数的值；

(12c)将恢复矢量的元素值更新为左移恢复矢量的元素值；

(12d)按照下式，更新目标到达角矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(L\right)}\left(p_I^{\left(L\right)}\left(s\right)\right)$

其中，表示更新后的目标到达角矢量，Θ^(L)表示左移角度搜索范围，表示左移峰值位置矢量的第s个元素，s表示动态信源数；

(12e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(12f)，否则，执行步骤(13)；

(12f)将当前代价函数的值更新为右移代价函数的值；

(12g)将恢复矢量的元素值更新为右移恢复矢量的元素值；

(12h)按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(R\right)}\left(p_I^{\left(R\right)}\left(s\right)\right)$

其中，表示更新后的目标到达角估计矢量，Θ^(R)表示右移角度搜索范围，表示右移峰值位置矢量的第s个元素，s表示动态信源数；

(13)判断可调栅格步长Δ是否大于栅格优化精度ξ_θ，若是，则执行步骤(10)，否则，执行步骤(14)；

(14)按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}\left(s\right)={\widehat{\mathrm{\θ}}}^{\left(1\right)}\left({p_I}^{\left(M\right)}\right)$

其中，表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，表示更新后的目标到达角估计矢量，表示当前峰值对应的位置下标；

(15)判断动态信源数是否小于信源数，若是，则执行步骤(7)，否则，执行步骤(16)；

(16)获得第二次目标到达角估计矢量：

将步骤(14)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。

本发明与现有技术相比具有如下优点：

第一，由于本发明采用对天线阵列所接收到的小样本回波数据进行阵元去噪，得到数据观测矢量，因而不用估计协方差矩阵，不仅克服了现有技术中由于样本数不足造成协方差矩阵估计不准导致目标到达角估计性能下降甚至失效的问题，而且克服了存在相干信号时协方差矩阵秩亏的问题，使得本发明可以在小样本与相干/独立源并存情况下仍然可以准确估计目标到达角。

第二，由于本发明通过迭代更新动态信源数来估计信源数，克服了现有技术的稀疏恢复DOA估计方法中需要准确已知信源数，而实际中因为信源数往往未知或难以估准导致目标到达角估计不准的问题，使得本发明在信源数未知的情况下，仍然可以准确估计目标到达角。

第三，由于本发明通过计算代价函数进而更新当前栅格参数，克服了现有技术的稀疏恢复DOA估计方法中由于栅格太密造成字典列与列相关性太高，导致稀疏恢复性能下降而无法准确估计目标到达角的问题，使得本发明可以在栅格较密时仍然有较好的稀疏恢复性能，从而准确估计目标到达角。

附图说明

图1是本发明的流程图；

图2是本发明的仿真图；

图3是本发明对6路数据采集器采集的天线阵列接收的实测回波数据进行目标到达角估计的仿真试验结果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明的具体实施步骤如下。

步骤1，接收小样本回波数据。

利用多路数据采集器获取天线阵列所接收到的小样本回波数据，其中，小样本回波数据是指在4到512个样本数之间所选取的回波数据。

步骤2，阵元去噪。

按照下式，对天线阵列所接收到的小样本回波数据进行阵元去噪，构建数据观测矢量：

$b=\frac{1}{Q}\underset{q=1}{\overset{Q}{\Σ}}x\left(t_q\right)x_1^{*}\left(t_q\right)$

其中，b表示N×1维的数据观测矢量，N表示天线阵列中的阵元个数，Q表示采样快拍数，q表示采样序号，x(t_q)表示t_q时刻天线阵列所接收到的N×1维小样本回波数据，x₁(t_q)表示t_q时刻天线阵列中第一个天线阵元所接收到的小样本回波数据，*表示共轭操作。

步骤3，计算粗栅格字典。

按照下式，计算导向矢量：

$a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{ccccccc}1& e^{\frac{j2\mathrm{\π}(d_2-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& e^{\frac{j2\mathrm{\π}(d_3-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_n-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_N-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}\end{array}\right]}^T$

其中，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量，N表示天线阵列中的阵元个数，θ_i∈Θ，∈表示属于符号，Θ表示角度搜索范围，i表示搜索角度θ_i在角度搜索范围Θ中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N，T表示转置操作。

按照下式，计算粗栅格字典：

A＝[a(θ₁)a(θ₂)…a(θ_i)…a(θ_L)]

其中，A表示N×L维的粗栅格字典，N表示天线阵列中的阵元个数，L表示总的粗栅格个数，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量。

步骤4，计算初始恢复矢量。

按照下式，计算粗栅格字典A的右逆矩阵：

P_a＝A^H(A·A^H)^-1

其中，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典，H表示共轭转置操作，-1表示求逆操作。

按照下式，计算粗栅格字典A的正交投影矩阵：

P＝I-P_aA

其中，P表示粗栅格字典A的正交投影矩阵，I表示单位矩阵，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典。

按照下式，计算初始恢复矢量：

${\tilde{r}}^0=P_a\·b$

其中，表示初始恢复矢量，P_a表示粗栅格字典A的右逆矩阵，b表示数据观测矢量。

步骤5，估计信源数。

(5a)初始化动态信源数为1。

(5b)按照下式，对第k次内循环中的恢复矢量进行排序操作：

$\[{\tilde{r}}^l,T\]=sort\left(\right|{\tilde{r}}^k\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示恢复矢量取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|,'descend')表示降序排列操作。

(5c)按照下式，计算第k+1次内循环中的恢复矢量：

${\tilde{r}}^{k+1}={\tilde{r}}^k+P(u^k-{\tilde{r}}^k)$

其中，表示第k+1次内循环中的恢复矢量，表示第k次内循环中的恢复矢量，u^k表示第k次内循环中的中间辅助矢量，P表示粗栅格字典A的正交投影矩阵。

第k次内循环中的中间辅助矢量由下式计算得到：

$\left\{\begin{array}{c}{u}_{T_k}^k={\tilde{r}}_{T_k}^k+{\mathrm{\χ}}^k{\left(A_{T_k}\right)}^H\left(A_{T_k}\right){\tilde{r}}_{T_k^c}^k\\ u_{T_k^c}^k=0\\ u^k=u_{T_k}^k\∪u_{T_k^c}^k\end{array}\right.$

其中，u^k表示第k次内循环中的中间辅助矢量，表示第k次内循环中的中间辅助矢量u对应下标子集T_k的元素构成的矢量，T_k表示索引指标集T中的前s个元素构成的下标子集，s表示动态信源数，表示第k次内循环时恢复矢量中对应下标子集T_k的元素构成的矢量，χ^k表示第k次内循环时的反馈参数，表示粗栅格字典A中对应下标子集T_k的列矢量构成的矩阵，H表示共轭转置操作，表示恢复矢量中对应下标子集的元素构成的矢量，表示索引指标集T中下标子集T_k的补集，表示第k次内循环时中间辅助矢量u中对应下标子集的元素构成的矢量。

(5d)按照下式，计算内循环相对误差值：

$H_2=\frac{\left|\right|{\tilde{r}}^{k+1}-{\tilde{r}}^k|{|}_2}{\left|\right|{\tilde{r}}^k|{|}_2}$

其中，H₂表示内循环相对误差值，和分别表示内循环次数为k+1和k时的恢复矢量，||·||₂表示取2范数操作。

(5e)判断内循环相对误差值H₂是否大于10^‐3，若是，则执行步骤(5b)，否则，执行步骤(5f)。

(5f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数。

(5g)按照下式，计算失配相对误差：

${\mathrm{\γ}}^{l+1}=\frac{\left|\right|{\mathrm{Au}}^k-b|{|}_2}{\left|\right|b|{|}_2}$

其中，γ^l+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示粗栅格字典，u^k表示第k次内循环中的中间辅助矢量，b表示数据观测矢量，||·||₂表示取2范数操作。

(5h)按照下式，计算外循环相对误差值：

H₁＝|γ^l+1-γ^l|

其中，H₁表示外循环相对误差值，γ^l+1和γ^l分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作。

(5i)判断外循环相对误差值H₁是否大于0.05，若是，则执行步骤(5b)，否则，执行步骤(5j)。

(5j)将外循环结束时的动态信源数值作为信源数的估计值。

步骤6，第一次估计目标到达角。

按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作。

将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量。

取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量中。

步骤7，按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作。

步骤8，按照下式，计算当前代价函数：

$F^{\left(M\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V\left(s\right)}{\underset{i\∉p_I(1:\hat{P})}{\Σ}\left|\tilde{r}\left(i\right)\right|}$

其中，F^(M)表示当前代价函数，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示恢复矢量中的第i个元素， 表示不属于符号，表示峰值位置矢量p_I的前个元素，表示信源数的估计值，|·|表示取模值操作。

步骤9，按照下式，计算当前峰值：

$\left\{\begin{array}{c}{p}_V^{\left(M\right)}=p_V\left(s\right)\\ p_I^{\left(M\right)}=p_I\left(s\right)\end{array}\right.$

其中，表示当前峰值，表示当前峰值对应的位置下标，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，p_I(s)表示峰值位置矢量p_I中的第s个元素。

步骤10，计算左移代价函数：

用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围。

按照下式，计算左移栅格字典：

其中，A^(L)为左移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(L)，∈表示属于符号，Θ^(L)表示左移角度搜索范围，i表示搜索角度θ_i在左移角度搜索范围Θ^(L)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N。

按照下式，计算左移栅格字典的正交投影矩阵：

P^(L)＝I-(A^(L))^H(A^(L)(A^(L))^H)^-1A^(L)

其中，P^(L)表示左移栅格字典A^(L)的正交投影矩阵，I表示单位矩阵，A^(L)表示左移栅格字典，H表示共轭转置操作，-1表示求逆操作。

按照下式，计算左移恢复矢量：

${\tilde{r}}^{\left(L\right)}={\tilde{r}}^{\left(L\right)}+P^{\left(L\right)}\left(u^{\left(L\right)}-{\tilde{r}}^{\left(L\right)}\right)$

其中，表示左移恢复矢量，u^(L)表示左移辅助矢量，P^(L)表示左移栅格字典A^(L)的正交投影矩阵。

按照下式，寻找左移峰值位置矢量：

$\[p_V^{\left(L\right)},p_I^{\left(L\right)}\]=findpeaks\left(\left|{\tilde{r}}^{\left(L\right)}\right|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示左移恢复矢量取模值后按降序重排的左移峰值矢量，表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作。

按照下式，计算左移代价函数：

$F^{\left(L\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(L\right)}\left(s\right)}{\underset{i\∉p_I^{\left(L\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(L\right)}\left(i\right)\right|}$

其中，F^(L)表示左移代价函数，表示左移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示左移恢复矢量中的第i个元素， 表示不属于符号，表示左移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作。

步骤11，计算右移代价函数。

用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围。

按照下式，计算右移栅格字典：

其中，A^(R)为右移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(R)，∈表示属于符号，Θ^(R)表示右移角度搜索范围，i表示搜索角度θ_i在右移角度搜索范围Θ^(R)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N。

按照下式，计算右移栅格字典的正交投影矩阵：

P^(R)＝I-(A^(R))^H(A^(R)(A^(R))^H)^-1A^(R)

其中，P^(R)表示右移栅格字典A^(R)的正交投影矩阵，I表示单位矩阵，A^(R)表示右移栅格字典，H表示共轭转置操作，-1表示求逆操作。

按照下式，计算右移恢复矢量：

${\tilde{r}}^{\left(R\right)}={\tilde{r}}^{\left(R\right)}+P^{\left(R\right)}\left(u^{\left(R\right)}-{\tilde{r}}^{\left(R\right)}\right)$

其中，表示右移恢复矢量，u^(R)表示右移辅助矢量，P^(R)表示右移栅格字典A^(R)的正交投影矩阵。

按照下式，寻找右移峰值位置矢量：

$\[p_V^{\left(R\right)},p_I^{\left(R\right)}\]=findpeaks\left(\left|{\tilde{r}}^{\left(R\right)}\right|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示右移恢复矢量取模值后按降序重排的右移峰值矢量，表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作。

按照下式，计算右移代价函数：

$F^{\left(R\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(R\right)}\left(s\right)}{\underset{i\∉p_i^{\left(R\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(R\right)}\left(i\right)\right|}$

其中，F^(R)表示右移代价函数，表示右移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示右移恢复矢量中的第i个元素， 表示不属于符号，表示右移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作。

步骤12，更新当前栅格参数。

(12a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(12b)，否则，执行步骤(12e)。

(12b)将步骤8的当前代价函数的值更新为左移代价函数的值。

(12c)将恢复矢量的元素值更新为左移恢复矢量的元素值。

(12d)按照下式，更新目标到达角矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(L\right)}\left(p_I^{\left(L\right)}\left(s\right)\right)$

其中，表示更新后的目标到达角矢量，Θ^(L)表示左移角度搜索范围，表示左移峰值位置矢量的第s个元素，s表示动态信源数。

(12e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(12f)，否则，执行步骤13。

(12f)将当前代价函数的值更新为右移代价函数的值。

(12g)将恢复矢量的元素值更新为右移恢复矢量的元素值。

(12h)按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(R\right)}\left(p_I^{\left(R\right)}\left(s\right)\right)$

其中，表示更新后的目标到达角估计矢量，Θ^(R)表示右移角度搜索范围，表示右移峰值位置矢量的第s个元素，s表示动态信源数。

步骤13，判断可调栅格步长Δ是否大于栅格优化精度ξ_θ，若是，则执行步骤10，否则，执行步骤14。

可调栅格步长从第二次内循环开始由上一次循环中可调栅格步长的0.5倍来更新，可调栅格步长的初始值由下式计算得到：

$\mathrm{\Δ}=\frac{180}{(5\×(N-1))}$

其中，Δ表示可调栅格步长，N表示天线阵列中的阵元个数。

栅格优化精度由下式计算得到：

${\mathrm{\ξ}}_{\mathrm{\θ}}=\frac{9}{(5\×(N-1))}$

其中，ξ_θ表示栅格优化精度，N表示天线阵列中的阵元个数。

步骤14，更新目标到达角估计矢量。

按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}\left(s\right)={\widehat{\mathrm{\θ}}}^{\left(1\right)}\left({p_I}^{\left(M\right)}\right)$

其中，表示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，表示更新后的目标到达角估计矢量，表示当前峰值对应的位置下标。

步骤15，判断动态信源数是否小于信源数，若是，则执行步骤7，否则，执行步骤16。

步骤16，获得第二次目标到达角估计矢量。

将步骤14中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。

下面结合附图2对本发明在小样本情况下对目标到达角估计的效果做进一步的描述。

1.仿真条件：

本发明的仿真运行系统为Intel(R)Core(TM)2i5‐3570CPU3.40GHz，64位Windows操作系统，仿真软件采用MATLABR(2013a)。

仿真参数设置如下表所示：

参数	参数值
		载频	1.5GHz
相控阵阵元个数	32
		阵元间距	0.1m
目标个数	3
		目标中频	5.3MHz，5.6MHz，5.9MHz
目标角度	0°,30°,45°
		信噪比	10dB，5dB，10dB
快拍数	8

2.仿真结果分析：

图2(a)表示采用现有技术的MUSIC方法形成的波达方向估计曲线，图2(a)中横坐标表示角度，纵坐标表示归一化幅度。图2(a)中的曲线表示在各个角度上的归一化幅度值大小情况，曲线中的两个尖峰表示2个真实目标处的归一化幅度值。图2(b)表示采用本发明形成的波达方向估计曲线，图2(b)中横坐标表示角度，纵坐标表示归一化幅度。图2(b)中曲线表示在各个角度上的归一化幅度值大小情况，曲线中的三个尖峰表示3个真实目标处的归一化幅度值。

由图2(a)中现有技术的MUSIC方法形成的波达方向估计曲线可见，现有技术的MUSIC方法形成的波达方向估计曲线只出现了2个代表目标的尖峰，与仿真条件中的目标个数相比，显然没有分辨出仿真条件中的3个目标。由图2(b)中本发明形成的波达方向估计曲线可见，本发明出现了3个代表目标的尖峰，成功的在小样本的情况下分辨出3个目标。下表具体说明三种方法的到达角方向估计结果。

	目标1估计值	目标2估计值	目标3估计值
				MUSIC方法	‐0.1161°	\	44.83°
本发明	‐0.03629°	29.79	44.96°

由上表可知，本发明成功的估计出了3个目标的到达角方向，且估计误差较小，而现有技术的MUSIC方法只估计出了2个目标的到达角方向，且估计误差较大。

下面结合附图3对本发明基于6阵元实测数据在信源数未知且存在相干源的情况下对目标到达角估计的效果做进一步的描述。

1.仿真试验条件：

本发明的仿真试验中的天线阵列采用6元等距线阵，阵元间距为7.5cm，接收信号频段S频段。试验中在2.25GHz处设定2个相干源，2.26GHz处设定1个独立源，经测量3个目标信号分别从-30°，-12°和10°三个方位发射。6阵元接收到回波信号后，下变频本振频率2.15GHz，采样频率1GHz，获得6阵元数字量化实测数据。

2.仿真试验结果分析：

图3(a)表示对天线阵列接收的回波数据做傅里叶变换后得到的频谱图，图3(a)中横坐标表示频率，纵坐标表示归一化幅度，曲线表示在各个频率上归一化幅度大小情况，曲线中的两个尖峰表示3个真实目标的归一化幅度，其中两个尖峰中较高的尖峰值为2个相干的真实目标的归一化幅度。图3(b)表示采用现有技术的MUSIC方法处理天线阵列接收的回波数据形成的到达角估计曲线，图3(b)中横坐标表示角度，纵坐标表示归一化幅度，图3(b)中的曲线表示在各个角度上的归一化幅度值大小情况。曲线中的一个尖峰表示1个真实目标处的归一化幅度。图3(c)表示采用本发明处理天线阵列接收的回波数据形成的到达角估计曲线，图3(c)中横坐标表示角度，纵坐标表示归一化幅度，图3(c)中的曲线表示在各个角度上的归一化幅度值大小情况，3个尖峰表示3个真实目标处的归一化幅度。

由图3(a)中对天线阵列接收的回波数据做傅里叶变换后得到的频谱图可见，频谱中有两个尖峰，分别位于99.61MHz和110.4MHz，由于中频信号频率为目标射频频率减去下变频本振信号频率，因此3个真实目标中有2个真实目标的频率相同，为相干源。由图3(b)中现有技术的MUSIC方法处理天线阵列接收的回波数据形成的到达角估计曲线可见，现有技术的MUSIC方法只出现1个代表目标的尖峰，说明只测出了一个独立源目标的到达角，与本发明所估计出的目标数相比，丢失了2个目标。由图3中(c)本发明处理天线阵列接收的回波数据形成的到达角估计曲线可见，本发明出现了3个代表目标的尖峰，成功的在未知信源数且存在相干信号的情况下分辨出3个目标。下表具体说明了两种方法的到达角方向估计结果。

	目标1	目标2	目标3
				目标方位真实值	‐30°	‐12°	10°
本发明估计值	‐29.92°	‐11.02°	9.9°
				MUSIC方法估计值	\	‐10.8°	\

由上表可知，本发明成功的估计出了3个目标的到达角方向，且估计误差较小，而现有技术的MUSIC方法只估计出了1个独立源目标的到达角方向。

Claims (7)

1.一种基于信源数估计的栅格偏移优化目标到达角估计方法，包括如下步骤：

(1)接收小样本回波数据：

利用多路数据采集器获取天线阵列所接收到的小样本回波数据；

(2)阵元去噪：

按照下式，对天线阵列所接收到的小样本回波数据进行阵元去噪，得到数据观测矢量：

$b=\frac{1}{Q}\underset{q=1}{\overset{Q}{\Σ}}x\left(t_q\right)x_1^{*}\left(t_q\right)$

其中，b表示N×1维的数据观测矢量，N表示天线阵列中的阵元个数，Q表示采样快拍数，q表示采样序号，x(t_q)表示t_q时刻天线阵列所接收到的N×1维小样本回波数据，x₁(t_q)表示t_q时刻天线阵列中第一个天线阵元所接收到的小样本回波数据，*表示共轭操作；

(3)计算粗栅格字典：

(3a)按照下式，计算导向矢量：

$a\left({\mathrm{\θ}}_i\right)={\left[\begin{array}{ccccccc}1& e^{\frac{j2\mathrm{\π}(d_2-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& e^{\frac{j2\mathrm{\π}(d_3-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_n-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}& ...& e^{\frac{j2\mathrm{\π}(d_N-d_1){\mathrm{sin\θ}}_i}{\mathrm{\λ}}}\end{array}\right]}^T$

其中，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量，N表示天线阵列中的阵元个数，θ_i∈Θ，∈表示属于符号，Θ表示角度搜索范围，i表示搜索角度θ_i在角度搜索范围Θ中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N，T表示转置操作；

(3b)按照下式，计算粗栅格字典：

A＝[a(θ₁)a(θ₂)…a(θ_i)…a(θ_L)]

其中，A表示N×L维的粗栅格字典，N表示天线阵列中的阵元个数，L表示总的粗栅格个数，a(θ_i)表示搜索角度θ_i处的N×1维导向矢量；

(4)计算初始恢复矢量：

(4a)按照下式，计算粗栅格字典A的右逆矩阵：

P_a＝A^H(A·A^H)^-1

其中，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典，H表示共轭转置操作，-1表示求逆操作；

(4b)按照下式，计算粗栅格字典A的正交投影矩阵：

P＝I-P_aA

其中，P表示粗栅格字典A的正交投影矩阵，I表示单位矩阵，P_a表示粗栅格字典A的右逆矩阵，A表示粗栅格字典；

(4c)按照下式，计算初始恢复矢量：

${\tilde{r}}^0=P_a\·b$

其中，表示初始恢复矢量，P_a表示粗栅格字典A的右逆矩阵，b表示数据观测矢量；

(5)估计信源数：

(5a)将动态信源数初始化为1；

(5b)按照下式，对第k次内循环中的恢复矢量进行排序操作：

$\[{\tilde{r}}^l,T\]=sort\left(\right|{\tilde{r}}^k\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示恢复矢量取模值后按降序重排的矢量，l表示外循环次数，k表示内循环次数，T表示排序操作后记录中每一个元素在恢复矢量中对应元素的下标组成的索引指标集，|·|表示取模值操作，sort(|·|,'descend')表示降序排列操作；

(5c)按照下式，计算第k+1次内循环中的恢复矢量：

${\tilde{r}}^{k+1}={\tilde{r}}^k+P(u^k-{\tilde{r}}^k)$

其中，表示第k+1次内循环中的恢复矢量，表示第k次内循环中的恢复矢量，u^k表示第k次内循环中的中间辅助矢量，P表示粗栅格字典A的正交投影矩阵；

(5d)按照下式，计算内循环相对误差值：

$H_2=\frac{\left|\right|{\tilde{r}}^{k+1}-{\tilde{r}}^k|{|}_2}{\left|\right|{\tilde{r}}^k|{|}_2}$

其中，H₂表示内循环相对误差值，和分别表示内循环次数为k+1和k时的恢复矢量，||·||₂表示取2范数操作；

(5e)判断内循环相对误差值H₂是否大于10^-3，若是，则执行步骤(5b)，否则，执行步骤(5f)；

(5f)将第l次外循环中的动态信源数加1，用加1后的动态信源数作为下一次外循环中的动态信源数；

(5g)按照下式，计算失配相对误差：

${\mathrm{\γ}}^{l+1}=\frac{\left|\right|{\mathrm{Au}}^k-b|{|}_2}{\left|\right|b|{|}_2}$

其中，γ^l+1表示第l+1次外循环中的失配相对误差，l表示外循环次数，A表示粗栅格字典，u^k表示第k次内循环中的中间辅助矢量，b表示数据观测矢量，||·||₂表示取2范数操作；

(5h)按照下式，计算外循环相对误差值：

H₁＝|γ^l+1-γ^l|

其中，H₁表示外循环相对误差值，γ^l+1和γ^l分别表示外循环次数为l+1和l时的失配相对误差，|·|表示取模值操作；

(5i)判断外循环相对误差值H₁是否大于0.05，若是，则执行步骤(5b)，否则，执行步骤(5j)；

(5j)将外循环结束时的动态信源数值作为信源数的估计值；

(6)第一次估计目标到达角：

(6a)按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(6b)将峰值位置矢量中的第一个元素值，放入信源位置矢量的第一个位置，将峰值位置矢量中的下一个元素值依次放入信源位置矢量的第二个位置，直至放入信源位置矢量的元素个数与估计信源数的值相同时停止取值，得到最终信源位置矢量；

(6c)取出角度搜索范围Θ中与最终信源位置矢量中元素值相同的下标值对应的角度值，将取出的角度值放入第一次估计的目标到达角矢量中；

(7)按照下式，寻找峰值位置矢量：

$\[p_V,p_I\]=findpeaks\left(\right|\tilde{r}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，p_V表示恢复矢量取模值后按降序重排的峰值矢量，p_I表示恢复矢量元素取模值降序重排后，恢复矢量原下标值经过重新排列得到的峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(8)按照下式，计算当前代价函数：

$F^{\left(M\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V\left(s\right)}{\underset{i\∉p_I(1:\hat{P})}{\Σ}\left|\tilde{r}\left(i\right)\right|}$

其中，F^(M)表示当前代价函数，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示恢复矢量中的第i个元素， 表示不属于符号，表示峰值位置矢量p_I的前个元素，表示信源数的估计值，|·|表示取模值操作；

(9)按照下式，计算当前峰值：

$\left\{\begin{array}{c}{p}_V^{\left(M\right)}=p_V\left(s\right)\\ p_I^{\left(M\right)}=p_I\left(s\right)\end{array}\right.$

其中，表示当前峰值，表示当前峰值对应的位置下标，p_V(s)表示峰值矢量p_V中的第s个元素，s表示动态信源数，p_I(s)表示峰值位置矢量p_I中的第s个元素；

(10)计算左移代价函数：

(10a)用角度搜索范围中当前位置下标对应的角度减去可调栅格步长后得到的搜索角度范围作为左移角度搜索范围；

(10b)按照下式，计算左移栅格字典：

其中，A^(L)表示左移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(L)，∈表示属于符号，Θ^(L)表示左移角度搜索范围，i表示搜索角度θ_i在左移角度搜索范围Θ^(L)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N；

(10c)按照下式，计算左移栅格字典的正交投影矩阵：

P^(L)＝I-(A^(L))^H(A^(L)(A^(L))^H)^-1A^(L)

其中，P^(L)表示左移栅格字典A^(L)的正交投影矩阵，I表示单位矩阵，A^(L)表示左移栅格字典，H表示共轭转置操作，-1表示求逆操作；

(10d)按照下式，计算左移恢复矢量：

${\tilde{r}}^{\left(L\right)}={\tilde{r}}^{\left(L\right)}+P^{\left(L\right)}(u^{\left(L\right)}-{\tilde{r}}^{\left(L\right)})$

其中，表示左移恢复矢量，u^(L)表示左移辅助矢量，P^(L)表示左移栅格字典A^(L)的正交投影矩阵；

(10e)按照下式，寻找左移峰值位置矢量：

$\[p_V^{\left(L\right)},p_I^{\left(L\right)}\]=findpeaks\left(\right|{\tilde{r}}^{\left(L\right)}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示左移恢复矢量取模值后按降序重排的左移峰值矢量，表示左移恢复矢量元素取模值降序重排后，左移恢复矢量原下标值经过重新排列得到的左移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(10f)按照下式，计算左移代价函数：

$F^{\left(L\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(L\right)}\left(s\right)}{\underset{i\∉p_I^{\left(L\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(L\right)}\left(i\right)\right|}$

其中，F^(L)表示左移代价函数，表示左移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示左移恢复矢量中的第i个元素， 表示不属于符号，表示左移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(11)计算右移代价函数：

(11a)用角度搜索范围中当前位置下标对应的角度加上可调栅格步长后得到的搜索角度范围作为右移角度搜索范围；

(11b)按照下式，计算右移栅格字典：

其中，A^(R)为右移栅格字典，N表示天线阵列中的阵元个数，θ_i∈Θ^(R)，∈表示属于符号，Θ^(R)表示右移角度搜索范围，i表示搜索角度θ_i在右移角度搜索范围Θ^(R)中的序号，e表示以自然常数为底的指数操作，j表示虚数单位，λ表示天线阵列的工作波长，d_n表示天线阵列中第n个天线阵元的坐标值，n＝1,2,…,N；

(11c)按照下式，计算右移栅格字典的正交投影矩阵：

P^(R)＝I-(A^(R))^H(A^(R)(A^(R))^H)^-1A^(R)

其中，P^(R)表示右移栅格字典A^(R)的正交投影矩阵，I表示单位矩阵，A^(R)表示右移栅格字典，H表示共轭转置操作，-1表示求逆操作；

(11d)按照下式，计算右移恢复矢量：

${\tilde{r}}^{\left(R\right)}={\tilde{r}}^{\left(R\right)}+P^{\left(R\right)}(u^{\left(R\right)}-{\tilde{r}}^{\left(R\right)})$

其中，表示右移恢复矢量，u^(R)表示右移辅助矢量，P^(R)表示右移栅格字典A^(R)的正交投影矩阵；

(11e)按照下式，寻找右移峰值位置矢量：

$\[p_V^{\left(R\right)},p_I^{\left(R\right)}\]=findpeaks\left(\right|{\tilde{r}}^{\left(R\right)}\left|{,}^{\′}{\mathrm{descend}}^{\′}\right)$

其中，表示右移恢复矢量取模值后按降序重排的右移峰值矢量，表示右移峰值矢量元素取模值降序重排后，右移恢复矢量原下标值经过重新排列得到的右移峰值位置矢量，findpeaks(·,'descend')表示寻找局部峰值并将局部峰值按降序排列，|·|表示取模值操作；

(11f)按照下式，计算右移代价函数：

$F^{\left(R\right)}=\frac{\underset{s=1}{\overset{\hat{P}}{\Σ}}{p}_V^{\left(R\right)}\left(s\right)}{\underset{i\∉p_I^{\left(R\right)}(1:\hat{P})}{\Σ}\left|{\tilde{r}}^{\left(R\right)}\left(i\right)\right|}$

其中，F^(R)表示右移代价函数，表示右移峰值矢量中的第s个元素，s表示动态信源数，Σ(·)表示求和，表示右移恢复矢量中的第i个元素， 表示不属于符号，表示右移峰值位置矢量的前个元素，表示信源数的估计值，|·|表示取模值操作；

(12)更新当前栅格参数：

(12a)判断左移代价函数是否同时大于当前代价函数和右移代价函数，若是，则执行步骤(12b)，否则，执行步骤(12e)；

(12b)将步骤(8)的当前代价函数的值更新为左移代价函数的值；

(12c)将恢复矢量的元素值更新为左移恢复矢量的元素值；

(12d)按照下式，更新目标到达角矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(L\right)}\left(p_I^{\left(L\right)}\right(s\left)\right)$

其中，表示更新后的目标到达角矢量，Θ^(L)表示左移角度搜索范围，表示左移峰值位置矢量的第s个元素，s表示动态信源数；

(12e)判断右移代价函数是否同时大于当前代价函数和左移代价函数，若是，则执行步骤(12f)，否则，执行步骤(13)；

(12f)将当前代价函数的值更新为右移代价函数的值；

(12g)将恢复矢量的元素值更新为右移恢复矢量的元素值；

(12h)按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}={\mathrm{\Θ}}^{\left(R\right)}\left(p_I^{\left(R\right)}\right(s\left)\right)$

其中，表示更新后的目标到达角估计矢量，Θ^(R)表示右移角度搜索范围，表示右移峰值位置矢量的第s个元素，s表示动态信源数；

(13)判断可调栅格步长Δ是否大于栅格优化精度ξ_θ，若是，则执行步骤(10)，否则，执行步骤(14)；

(14)按照下式，更新目标到达角估计矢量：

${\widehat{\mathrm{\θ}}}^{\left(1\right)}\left(s\right)={\widehat{\mathrm{\θ}}}^{\left(1\right)}\left({p_I}^{\left(M\right)}\right)$

其中，示更新后的目标到达角估计矢量的第s个元素值，s表示动态信源数，表示更新后的目标到达角估计矢量，表示当前峰值对应的位置下标；

(15)判断动态信源数是否小于信源数，若是，则执行步骤(7)，否则，执行步骤(16)；

(16)获得第二次目标到达角估计矢量：

将步骤(14)中更新后的目标到达角矢量作为第二次的目标到达角估计矢量。

2.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(1)中所述的小样本回波数据是指在4到512个样本数之间所选取的回波数据。

3.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(5b)中所述的外循环次数l的初始值为1，l的终值是外循环相对误差值小于等于外循环误差门限时的外循环次数。

4.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(5b)中所述的内循环次数k的初始值为0，k的终值是内循环相对误差值小于等于内循环误差门限时的内循环次数。

5.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(5c)中所述的第k次内循环中的中间辅助矢量由下式计算得到：

$\left\{\begin{array}{c}{u}_{T_k}^k={\tilde{r}}_{T_k}^k+{\mathrm{\χ}}^k{\left(A_{T_k}\right)}^H\left(A_{T_k}\right){\tilde{r}}_{T_k^c}^k\\ u_{T_k^c}^k=0\\ u^k=u_{T_k}^k\∪u_{T_k^c}^k\end{array}\right.$

其中，u^k表示第k次内循环中的中间辅助矢量，表示第k次内循环中的中间辅助矢量u对应下标子集T_k的元素构成的矢量，T_k表示索引指标集T中的前s个元素构成的下标子集，s表示动态信源数，表示第k次内循环时恢复矢量中对应下标子集T_k的元素构成的矢量，χ^k表示第k次内循环时的反馈参数，表示粗栅格字典A中对应下标子集T_k的列矢量构成的矩阵，H表示共轭转置操作，表示恢复矢量中对应下标子集的元素构成的矢量，表示索引指标集T中下标子集T_k的补集，表示第k次内循环时中间辅助矢量u中对应下标子集的元素构成的矢量。

6.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(13)中所述的可调栅格步长从第二次内循环开始由上一次内循环中可调栅格步长的0.5倍来更新，可调栅格步长的初始值由下式计算得到：

$\mathrm{\Δ}=\frac{180}{(5\×(N-1\left)\right)}$

其中，Δ表示可调栅格步长，N表示天线阵列中的阵元个数。

7.根据权利要求1所述的基于信源数估计的栅格偏移优化目标到达角估计方法，其特征在于：步骤(13)中所述的栅格优化精度由下式计算得到：

${\mathrm{\ξ}}_{\mathrm{\θ}}=\frac{9}{(5\×(N-1\left)\right)}$

其中，ξ_θ表示栅格优化精度，N表示天线阵列中的阵元个数。