本发明属于信号处理
技术领域：
，特别涉及一种基于时间反演的低仰角波达方向估计方法，适用于低仰角条件下相干目标的波达方向估计，进而实现目标定位与跟踪。
背景技术：
：波达方向DOA估计是雷达、声纳信号处理中的一个重要研究方向，低空目标探测中的低仰角DOA估计问题一直是DOA估计领域的研究热点；因为在低仰角条件下，复杂的地形会产生大量的多径信号，多径信号与直达波信号混叠在一起，使得目标信息提取和DOA估计变得很困难。国内外，现有的以直达波为有效信息的DOA估计模型，一般采用通道等效或者滤波技术抑制多径信号影响，以提高DOA估计效果。但是在低空目标探测问题中，由于存在大量的多径信号，这种处理方法通常无法取得理想的效果。因此，低空目标探测，即低仰角的DOA估计是DOA估计领域的一个难题。目前常用的地仰角DOA估计方法主要有以下几种：(1)基于精确信号模型的最大似然方法：该方法利用天线高度、目标高度、目标和雷达之间的距离、地球曲率、大气折射等先验信息，充分考虑多径信号的影响，建立精确的多径信号反射模型，以减少待估计参数个数，提高参数估计精度；由于实际情况的复杂性、反射角和反射系数往往不能精确、甚至无法得到，因此该方法对环境适应能力很差，进而大大限制了该方法的应用。(2)基于地形匹配的等效DOA估计模型方法：具体来说，该方法首先对地形信心进行探索，然后建立相应的地形信息表；测量前，线进行地形匹配操作，根据地形匹配信息建立相应的等效DOA估计模型，从而进行DOA估计；陈伯孝、朱伟在CircuitsSystemSignalProcess期刊的2013年第32卷第2期第647–662页发表的“AltitudemeasurementbasedonterrainmatchinginVHFarrayradar”详细论述了地形匹配算法，并将该方法应用到了米波雷达测高领域；然而该方法有一个缺陷，即实际雷达目标探测过程中，地形信息很复杂，并不完全已知，因而无法保证地形的完全匹配，进而使得该方法在复杂环境中的应用受到了限制。(3)基于直达波信号的反射信号综合建模DOA估计方法：该算法对直达波信号的反射信号综合建模，并将二者的导向矢量进行合成得到DOA估计模型，然后进行DOA估计；该方法模型的准确性与入射角、反射角的准确性紧密相连；复杂地形中，多径信息很复杂，通常无法准确知道反射信号信息，必然导致一定程度上的模型失配，进而对DOA的估计性能产生了影响。技术实现要素：针对上述现有技术存在的不足，本发明的目的在于提出一种基于时间反演的低仰角波达方向估计方法，该种基于时间反演的低仰角波达方向估计方法基于反演宽带聚焦空间平滑MUSIC算法TRF-SSMUSIC，并以探索复杂环境中低仰角条件下相干目标的波达方向估计，提高波达方向DOA的估计精度，进而实现目标定位与跟踪。本发明的技术思路是:确定天线阵列A，所述天线阵列A包含M个阵元，并且所述天线阵列A发送前向探测信号，所述前向探测信号经过目标散射后传播至天线阵列A，得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)，并计算天线阵列A中第n个阵元发射的时间反演探测信号zn(t)，并计算得到天线阵列A中第n个阵元的时间反演接收信号进而计算天线阵列A中第n个阵元的采样频率点ωq处时间反演频域接收信号的协方差矩阵对进行宽带聚焦处理，依次计算第n个阵元的宽带聚焦频率ω0处的聚焦协方差矩阵和第n个阵元的宽带聚焦频率ω0处的聚焦平滑矩阵并计算天线阵列A的空间谱函数S(θ)，进而计算天线阵列A的空间谱函数S(θ)的H个极大值，每一个极大值对应的角度为波达方向。为达到上述技术目的，本发明采用如下技术方案予以实现。一种基于时间反演的低仰角波达方向估计方法，包括以下步骤：步骤1，确定天线阵列A，所述天线阵列A包含M个阵元，并且所述天线阵列A发送前向探测信号，所述前向探测信号经过目标散射后传播至天线阵列A，得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)；t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数，n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数；所述第n个阵元为天线阵列A中的任意一个阵元；步骤2，对天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)进行时间反演操作，得到天线阵列A中第n个阵元发射的时间反演探测信号zn(t)；n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数；步骤3，根据天线阵列A中第n个阵元发射的时间反演探测信号zn(t)，计算得到天线阵列A中第n个阵元的时间反演接收信号n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数；步骤4，根据天线阵列A中第n个阵元的时间反演接收信号计算得到第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵n∈{1,2,..,M}，q∈{1,…,Q}，ωq表示第q个采样点的采样频率，Q表示天线阵列A中第n个阵元的时间反演频域接收信号进行带内采样的采样点个数，M表示天线阵列A包含的阵元个数；步骤5，对第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵进行宽带聚焦处理，依次计算得到第n个阵元的宽带聚焦频率ω0和第n个阵元的宽带聚焦频率ω0处的聚焦协方差矩阵步骤6，对第n个阵元的宽带聚焦频率ω0处的聚焦协方差矩阵进行解相干处理，计算得到第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵步骤7，根据第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵计算得到天线阵列A中第n个阵元的空间谱函数S(θ)，进而计算天线阵列A中第n个阵元的空间谱函数S(θ)的H个极大值，每一个极大值对应的角度为波达方向；其中，H＝P'，P'表示对第n个阵元的宽带聚焦频率ω0处的聚焦协方差矩阵划分的子阵个数。本发明与现有低仰角DOA估计方法相比具有以下优点：1)更强的稳健性：传统方法将多径信息作为DOA估计的不利信息处理，致力于消除或者抑制回波信号中的多径信号；在低角探测环境中，多径信息复杂程度增加，因此多径信号的抑制效果会退化，因而引起整个DOA估计性能的恶化。本发明利用多径信号与时间反演操作提取信号传播媒介信息，这种自适应操作增强了模型的稳健性，进而提高了DOA估计的稳健性。2)更好的测角性能：现有的低仰角DOA估计算法，在多径信号严重的情况下由于直达波信号被多径信号干扰，DOA估计效果会严重退化甚至失效；而本发明采用时间反演技术，可以同时充分利用直达波和多径信息，在相同条件下可以提供更高的信噪比SNR和有效孔径，因而使得DOA估计拥有更高的测角分辨率和测角精确度，拥有更好的测角性能。附图说明下面结合附图和具体实施方式对本发明作进一步详细说明。图1是本发明的一种基于时间反演的低仰角波达方向估计方法流程图；图2是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图；图3是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图；图4是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图；图5是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图；图6是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图；图7是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图；图8是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图；图9是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图；图10是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图。具体实施方式参照图1，为本发明的一种基于时间反演的低仰角波达方向估计方法流程图；所述基于时间反演的低仰角波达方向估计方法，包括以下步骤：步骤1，确定天线阵列A，所述天线阵列A包含M个阵元，并且所述天线阵列A发送前向探测信号，所述前向探测信号经过目标散射后传播至天线阵列A，得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)；t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数，n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，所述第n个阵元为天线阵列A中的任意一个阵元。步骤1的子步骤为：1a)将所述天线阵列A发射前向探测信号记为f(t)，其表达式为：f(t)=f^(t)ejωct]]>其中，表示前向探测信号f(t)的包络，ωc表示前向探测信号f(t)的载波频率，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数。1b)前向探测信号f(t)通过媒介传播至目标处，经过目标散射再传播至天线阵列A，天线阵列A接收经过目标散射传播回来的前向回波时域信号；其中，将天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号记为r(n,k)(t)，其表达式为：r(n,k)(t)=Σi=1PX(i)f(t-τ(n,1,i)-Δτ(k,i))+v(n,k)(t)]]>其中，n∈{1,2,..,M}，k∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，f(t)表示天线阵列A发射的前向探测信号，i∈{1,2,..,P}，i表示第i条多径传播，P表示前向探测信号经过目标散射后传播至天线阵列A时包含的多径传播总条数，τ(n,1,i)表示由天线阵列A的第n个阵元发射的信号经过目标散射通过第i条多径传播回到天线阵列A的第1个阵元的传播时间，△τ(k,i)表示天线阵列A中第k个阵元经过目标散射通过第i条多径传播的传播延时，v(n,k)(t)表示天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号的前向通道噪声，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数，X(i)表示第i条多径传播的能量衰减因子。1c)根据天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号r(n,k)(t)，分别得到天线阵列A中第n个阵元发射前向探测信号、第1个阵元接收的前向回波时域信号r(n,1)(t)到天线阵列A中第n个阵元发射前向探测信号、第M个阵元接收的前向回波时域信号r(n,M)(t)，并将所述天线阵列A中第n个阵元发射的前向探测信号、第1个阵元接收的前向回波时域信号r(n,1)(t)到天线阵列A中第n个阵元发射的前向探测信号、第M个阵元接收的前向回波时域信号r(n,M)(t)按照阵列天线A中的阵元顺序排列成向量，得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)，其表达式为；Yn(t)＝[r(n,1)(t),…,r(n,M)(t)]T。其中，[·]T表示转置，n∈{1,2,..,M}，k∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数。步骤2，对天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)进行时间反演操作，得到天线阵列A中第n个阵元发射的时间反演探测信号zn(t)；n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数。具体地，在文献IEEETRANSACTIONONSIGNALPROCESSING，第59卷，第6期，2655-2668中发表的“TimeReversalBasedActiveArraySourceLocalization”指出，时间反演操作直接在时域进行过程复杂，一种简便实用的方法是：首先将待处理信号通过傅里叶变换，变换到频域，这样时域的时间反演操作就简化为频域的共轭操作，然后再进行能量归一化，大大简化了处理流程；本发明方法基于时间反演宽带聚焦空间平滑MUSIC算法TRF-SSMUSIC，并在频域进行时间反演操作。2a)对天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波时域信号Yn(t)进行快速傅里叶变换，其中天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号r(n,k)(t)经过快速傅里叶变换后，得到天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波频域信号R(n,k)(ω)，其表达式为：R(n,k)(ω)=FFT[r(n,k)(t)]=Σi=1PX(i)·F(ω)·e-jωτ(n,1,i)·e-jωΔτ(k,i)+V(n,k)(ω),]]>F(ω)＝FFT[f(t)],V(n,k)(ω)＝FFT[v(n,k)(t)]其中，FFT[·]表示快速傅里叶变换，F(ω)表示天线阵列A发射的前向探测信号f(t)的频域表达形式，v(n,k)(t)表示天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号的前向通道噪声，V(n,k)(ω)表示天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号的前向通道噪声v(n,k)(t)的频域表达形式，X(i)表示第i条多径传播的能量衰减因子，ω表示前向回波时域信号的角频率，i∈{1,2,..,P}，i表示第i条多径传播，P表示前向探测信号经过目标散射后传播至天线阵列A时包含的多径传播总条数，τ(n,1,i)表示由天线阵列A的第n个阵元发射的信号经过目标散射通过第i条多径传播回到阵列A的第1个阵元的传播时间，△τ(k,i)表示天线阵列A中第k个阵元经过目标散射通过第i条多径传播的传播延时，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数。2b)将天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波频域信号R(n,k)(ω)，分别得到天线阵列A中第n个阵元发射前向探测信号、第1个阵元接收的前向回波频域信号R(n,1)(ω)到天线阵列A中第n个阵元发射前向探测信号、第M个阵元接收的前向回波频域信号R(n,M)(ω)，然后将所述天线阵列A中第n个阵元发射的前向探测信号、第1个阵元接收的前向回波频域信号R(n,1)(ω)到天线阵列A中第n个阵元发射的前向探测信号、第M个阵元接收的前向回波频域信号R(n,M)(ω)按照天线阵列A中的阵元顺序排列成向量，得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波频域信号Yn(ω)，其表达式为：Yn(ω)=FFT[Yn(t)]=[FFT[r(n,1)(t)],...,FFT[r(n,k)(t)],...,FFT[r(n,M)(t)]T=[R(n,1)(ω),...,R(n,k)(ω),...,R(n,M)(ω)]T,]]>进而得到天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波频域信号Yn(ω)的矩阵形式Y^n(ω)=A(Θ)XΓn(ω)F(ω)+ϵn(ω)]]>其中，A(Θ)表示目标导向矢量，θ1,...θp分别第1条多径传播到第P条多径传播各自对应的目标角度，ω表示前向回波时域信号的角频率，d表示天线阵列A的阵元间距，Ω(1)表示天线阵列A中第2个阵元的波数，Ω(M-1)表示天线阵列A中第M个阵元的波数，c表示前向探测信号的传播速度，X表示目标传播衰减因子对角阵，X(1)表示第1条多径传播对应的目标传播衰减因子，X(P)表示第P条多径传播对应的目标传播衰减因子，τ(1,n,1)表示第n个阵元发射的探测信号经过第1条多径传播反射到达第1个阵元所经过的空间传播延时，Γn(ω)表示第n个阵元经过目标的空间传播延时，εn(ω)表示天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波频域信号Yn(ω)的前向通道传播噪声；i∈{1,2,..,P}，i表示第i条多径传播，P表示前向探测信号经过目标散射后传播至天线阵列A时包含的多径传播总条数，V(n,1)(ω)表示天线阵列A中第n个阵元发射前向探测信号、第1个阵元接收的前向回波时域信号的前向通道噪声v(n,1)(t)的频域表达形式，V(n,M)(ω)表示天线阵列A中第n个阵元发射前向探测信号、第M个阵元接收的前向回波时域信号的前向通道噪声v(n,M)(t)的频域表达形式，r(n,k)(t)表示天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波时域信号，R(n,k)(ω)表示天线阵列A中第n个阵元发射前向探测信号、第k个阵元接收的前向回波频域信号。2c)根据天线阵列A发射的前向探测信号f(t)的频域表达形式F(ω)和天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波频域信号Yn(ω)的矩阵形式计算得到天线阵列A中第n个阵元的能量归一化因子gn，并根据天线阵列A中第n个阵元的能量归一化因子gn对天线阵列A中第n个阵元接收经过目标散射传播回来的前向回波频域信号Yn(ω)的矩阵形式取共轭，得到天线阵列A中第n个阵元发射的时间反演探测信号zn(t)，其表达式分别为：gn=(||F(ω)||2)/(||Yn(ω)||2)]]>zn(ω)＝gnY*n(ω)zn(t)＝IFFT[zn(ω)]其中，||·||2表示2范数运算，[·]*表示共轭，IFFT[·]表示逆快速傅里叶变换。步骤3，根据天线阵列A中第n个阵元发射的时间反演探测信号zn(t)，计算得到天线阵列A中第n个阵元的时间反演接收信号n∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数。步骤3的子步骤为：3a)将天线阵列A中第n个阵元发射的时间反演探测信号zn(t)经过目标散射、且由天线阵列A中第k个阵元接收的时间反演信号记为其表达式为：其中，上标TR表示所述信号为时间反演信号，n∈{1,2,..,M}，k∈{1,2,..,M}，M表示天线阵列A包含的阵元个数，表示时间反演后向传播噪声，t∈{1,2,…,S}，S表示前向回波时域信号包含的时域快拍数，i∈{1,2,..,P}，i表示第i条多径传播，P表示前向探测信号经过目标散射后传播至天线阵列A时包含的多径传播总条数。3b)将天线阵列A中第n个阵元发射的时间反演探测信号zn(t)经过目标散射、且由天线阵列A中第k个阵元接收的时间反演信号按照天线阵列A的阵元顺序排列成向量，得到天线阵列A中第n个阵元的时间反演接收信号其表达式为：ynTR(t)=[r(n,1)TR(t),...,r(n,M)TR(t)]T]]>步骤4，根据天线阵列A中第n个阵元的时间反演接收信号计算得到第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵n∈{1,2,..,M}，q∈{1,…,Q}，ωq表示第q个采样点的采样频率，Q表示天线阵列A中第n个阵元的时间反演频域接收信号进行带内采样的采样点个数，M表示天线阵列A包含的阵元个数。步骤4的子步骤为：4a)对天线阵列A中第n个阵元的时间反演接收信号进行快速傅里叶变换，得到天线阵列A中第n个阵元的时间反演频域接收信号其表达式为：YnTR(ω)=FFT[ynTR(t)]=Σm=1MA(Θ)XΓm(ω)zm(ω)+ζn(ω)]]>其中，A(Θ)表示目标导向矢量，X表示目标传播衰减因子对角阵，Γm(ω)表示天线阵列A中第m个阵元发送前向探测信号的传播延时，τ(1,m,1)表示第m个阵元发射的探测信号经过第1条多径传播反射到达第1个阵元所经过的空间传播延时，τ(P,m,1)表示第m个阵元发射的探测信号经过第P条多径传播反射到达第1个阵元所经过的空间传播延时，所述传播延时只与目标有关，因此有Γ1(ω)＝Γn(ω)，zm(ω)表示天线阵列A中第m个阵元的发射信号，ζn(ω)表示天线阵列A中第n个阵元的时间反演后向传播噪声的频域表达形式，表示第n个阵元发射信号经过目标反射被第1个阵元接收过程中的噪声，表示第n个阵元发射信号经过目标反射被第M个阵元接收过程中的噪声，n∈{1,2,..,M}，m∈{1,2,..,M}，M表示天线阵列A包含的阵元个数。4b)对天线阵列A中第n个阵元的时间反演频域接收信号进行Q点带内采样，得到天线阵列A中第n个阵元的采样频率点ωq处的采样值并利用天线阵列A中第n个阵元的采样频率点ωq处的采样值计算得到第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵其表达式为：其中，q∈{1,…,Q}，ωq表示第q个采样点的采样频率，Q表示天线阵列A中第n个阵元的时间反演频域接收信号进行带内采样的采样点个数，表示数学期望运算，[·]H表示共轭转置。步骤5，对第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵进行宽带聚焦处理，依次计算得到第n个阵元的宽带聚焦频率ω0和第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵ωq表示第q个采样点的采样频率，q∈{1,…,Q}，Q表示天线阵列A中第n个阵元的时间反演频域接收信号进行带内采样的采样点个数。具体地，清华大学2004年11月出版的《空间谱估计理论与算法》指出，目前对宽带信号处理算法的研究重要分为两类，第一类为基于不相干信号的处理方法(IncoherentSignal-subspaceMethod,ISM)，第二类为基于相干信号的处理方法(CoherentSignal-SubspaceMethod,CSM)；由于信号模型是基于低角探测中的相干信源建立的，通常本实施例采用CSM处理方法，CSM处理方法分为：相干信号子空间算法(CSS)、双边相关变换法(TCT)和总体最小二乘法(TLS)等；本实施例采用基于协方差矩阵的TCT处理方法。5a)对第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵进行奇异值分解，得到第p条多径传播中采样频率点ωq点处协方差矩阵的奇异值p∈{1,2,...,P}，P表示前向探测信号经过目标散射后传播至天线阵列A时包含的多径传播总条数。5b)对第p条多径传播中采样频率点ωq点处协方差矩阵的奇异值运用最小误差准则，计算得到第n个阵元的宽带聚焦频率ω0：其中，表示取最小值操作，ω0表示第n个阵元的宽带聚焦频率。5c)将第n个阵元中宽带聚焦频率ω0处的特征向量矩阵记为U(ω0)，然后对第n个阵元中采样频率点ωq处时间反演频域接收信号的协方差矩阵进行特征值分解，得到第n个阵元中采样频率ωq处的特征向量矩阵U(ωq)，利用第n个阵元中宽带聚焦频率ω0处的特征向量矩阵U(ω0)和第n个阵元中采样频率ωq处的特征向量矩阵U(ωq)，计算得到第n个阵元中采样频率ωq处的变换矩阵Tq(ωq)，其表达式为：Tq(ωq)＝U(ω0)UH(ωq)5d)利用第n个阵元中采样频率ωq处时间反演频域接收信号的协方差矩阵和第n个阵元中采样频率ωq处的变换矩阵Tq(ωq)，计算得到第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵其表达式为：其中，q∈{1,…,Q}，Q表示天线阵列A中第n个阵元的时间反演频域接收信号进行带内采样的采样点个数，上标H表示共轭转置。步骤6，采用空间平滑算法对第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵进行解相干处理，计算得到第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵具体地，现有技术中常用的解相干方法是前后相空间平滑、托普列兹矩阵重构算法，本实施例中采用前后相空间平滑技术，其具体子步骤为：6a)将第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵划分成P'个子阵，得到P'个子阵各自在第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵，其中第l个子阵在第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵为并根据P'个子阵各自在第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵，计算得到P'个子阵在第n个阵元中宽带聚焦频率ω0处的前向平滑矩阵其表达式分别为：其中，l∈{1,2,...,P'}，表示取出·中第l行到第l+M-P'行，和·中第l列到第l+M-P'列而构成的矩阵，上标f表示前向平滑矩阵，P'表示对第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵划分的子阵个数。6b)将第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵进行上下倒序重排处理，得到上下倒序重排处理后第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵再对上下倒序重排处理后第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵进行左右倒序重排处理，得到左右倒序重排处理后第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵然后对左右倒序重排处理后第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵取共轭，得到重排后第n个阵元中宽带聚焦频率ω0处聚焦协方差共轭矩阵将重排后第n个阵元中宽带聚焦频率ω0处的聚焦协方差共轭矩阵划分为l'个子阵，得到l'个子阵各自在第n个阵元中宽带聚焦频率ω0处的聚焦协方差共轭矩阵，其中第l”个子阵在第n个阵元中宽带聚焦频率ω0处的聚焦协方差共轭矩阵为l”∈{1,…,l'}，根据l'个子阵各自在第n个阵元中宽带聚焦频率ω0处的聚焦协方差共轭矩阵，计算得到l'个子阵在第n个阵元中宽带聚焦频率ω0处的后向平滑矩阵其表达式分别为：其中，flipup[·]表示·的上下倒序重排操作，fliplr[·]表示·矩阵的左右倒序重排操作，[·]*表示共轭操作，表示取出·中第l”行到第l”+M-l'行，且取出·中第l”列到第l”+M-l'列而构成的矩阵；上标B表示后向平滑矩阵。6c)根据P'个子阵在第n个阵元中宽带聚焦频率ω0处的前向平滑矩阵和l'个子阵在第n个阵元中宽带聚焦频率ω0处的后向平滑矩阵计算得到第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵其表达式为：步骤7，根据第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵计算得到天线阵列A中第n个阵元的空间谱函数S(θ)，进而计算天线阵列A中第n个阵元的空间谱函数S(θ)的H个极大值，每一个极大值对应的角度为波达方向；其中，H＝P'，P'表示对第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵划分的子阵个数。具体地，现有技术中利用宽带聚焦频率ω0处的聚焦平滑矩阵估计出空间谱函数S(θ)有最大似然算法ML、加权子空间拟合法WSF、波束成形法DBF、基于特征空间的多重信号分类法MUSIC和旋转因子不变法ESPRIT；本实施例中采用多重信号分类法MUSIC，其具体子步骤为：7a)对第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵进行特征值分解，得到第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵的M'个特征值：λ1≥λ2≥λ3…≥λJ>λJ+1＝…λM'＝δ其中，λm”表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵的第m”个特征值，Vm”表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵对应于特征值λm”的特征矢量，λm'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵的第m'个特征值，Vm'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵对应于特征值λm'的特征矢量，m”∈{1,2,...,J}，m'∈{1,2,...,M'}，J表示将M'个特征值进行大小排序后大特征值的个数，且J＝P'，P'表示对第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵划分的子阵个数，δ表示M'个特征值的噪声功率，I表示M'×M'维单位矩阵，M'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵包含的特征值个数。7b)根据宽带聚焦频率ω0处的聚焦平滑矩阵的M'个特征值，计算得到天线阵列A中第n个阵元的空间谱函数S(θ)：其中，m”'∈{J+1,J+2,...,M'}，J表示将M'个特征值进行大小排序后大特征值的个数，M'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵包含的特征值个数，a(θ)表示入射角度为θ的天线阵列A入射信号的导向矢量，a(θ)＝[e-jΩ(0)sinθ,e-jΩ(1)sinθ,...,e-jΩ(k-1)sinθ,...,e-jΩ(M-1)sinθ]T,e-jΩ(k-1)sinθ表示入射角度为θ的天线阵列A入射信号在第k个阵元处的旋转因子，k∈{1,…,M}，θ表示天线阵列A的入射角度，Vm”'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵对应于特征值λm”'的特征矢量，λm”'表示第n个阵元中宽带聚焦频率ω0处的聚焦平滑矩阵的第m”'个特征值，d表示天线阵列A的阵元间距，c表示前向探测信号的传播速度，ω0表示第n个阵元的宽带聚焦频率。7c)根据天线阵列A中第n个阵元的空间谱函数S(θ)，计算得到天线阵列A中第n个阵元的空间谱函数S(θ)的H个极大值，每一个极大值对应的角度为波达方向；其中，H＝P'，P'表示对第n个阵元中宽带聚焦频率ω0处的聚焦协方差矩阵划分的子阵个数。通过以下计算仿真对本发明效果作进一步验证说明。本发明所有仿真共用以下公共条件：天线阵列A包含的阵元总数M＝16，频率采样点数Q＝10,直达波延时τ＝2ms，探测信号为中心频率为200MHz，带宽为20MHz的线性调频信号，脉宽为50μs，即探测信号ωc＝200MHz,τ0＝50μs,μ＝B/τ0,B＝20MHz其中，t表示时间变量。仿真1：对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法TRF-SSMUSIC四种算法的DOA估计谱图进行仿真。1.1)仿真条件：多径数目分别为2条，3条，4条，信噪比SNR＝5dB，探测角度θ范围，-30°≤θ≤30°，频域快拍数G＝20,1000次Monte-Carlo实验，对应参数配置见表1。表11.2)仿真内容与结果1.2a)在上述1.1仿真条件下，存在2条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法各自得到的谱图进行仿真，结果如图2所示，图2是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图。1.2b)在上述1.1仿真条件下，存在3条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法各自的谱图进行仿真，结果如图3所示，图3是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图。1.2c)在上述1.1仿真条件下，存在4条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法各自的谱图进行仿真，结果如图4所示，图4是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计谱图。从图2、图3和图4可以看出，当信噪比和阵元数相同时，使用本发明方法在2条多径、3条多径和4条多径条件下都具有最窄的主瓣以及最低的旁瓣，而且都能正确分辨出目标，说明本发明方法具有更高的角度分辨率。仿真2：对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法各自得到的波达方向DOA估计误差分布进行仿真。2.1)仿真条件：多径数目分别为2条、3条、4条，信噪比SNR＝5dB,探测角度θ范围为-30°≤θ≤30°，频域快拍数G＝20,1000次Monte-Carlo实验，其他仿真条件同仿真1。2.2)仿真内容与结果：2.2a)在上述2.1仿真条件下，存在2条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的误差分布进行仿真，结果如图5所示，图5是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图。2.2b)在上述2.1仿真条件下，存在3条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的误差分布进行仿真，结果如图6所示，图6是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图。2.2c)在上述2.1仿真条件下，存在4条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的误差分布进行仿真，结果如图7所示，图7是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计误差分布统计图。从图5、图6和图7可以看出，当信噪比和阵元数相同时，使用本发明方法分别在2条多径、3条多径、4条多径条件下本发明提出的时间反演宽带聚焦空间平滑MUSIC算法TRF-SSMUSIC都是无偏估计，而且具有最小的估计误差。仿真3：对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法、本发明方法TRF-SSMUSIC四种算法的DOA估计的均方根误差根误差与信噪比关系进行仿真。3.1)仿真条件：信噪比-6dB～20dB，频域快拍数G＝20,探测角度θ范围-30°≤θ≤30°，1000次Monte-Carlo实验，其他仿真条件同仿真1。3.2)仿真内容与结果：3.2a)在上述3.1仿真条件下，存在2条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的的均方根误差根误差与信噪比关系进行仿真，结果如图8所示，图8是2条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图。3.2b)在上述3.1仿真条件下，存在3条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的均方根误差根误差与信噪比关系进行仿真，结果如图9所示，图9是3条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图。3.2c)在上述3.1仿真条件下，存在4条多径，对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、时间反演宽带聚焦空间平滑MUSIC算法和本发明方法TRF-SSMUSIC四种算法进行波达方向DOA估计的均方根误差根误差与信噪比关系进行仿真，结果如图10所示，图10是4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差与信噪比关系图。从图8、图9和图10可以看出，在2条多径、3条多径和4条多径条件下分别使用传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC、本发明方法TRF-SSMUSIC得到的波达方向DOA估计均方根误差随着信噪比的增加而降低，并且本发明方法分别在,本发明提出的时间反演宽带聚焦空间平滑MUSIC算法TRF-SSMUSIC算法在2条多径、3条多径和4条多径条件下都具有最小的均方根误差，说明本发明具有更高的估计精度和更好的测角性能。综上，相比于现有对传统Capon算法ConvSSCapon、时间反演Capon算法TRSSCapon、传统宽带聚焦空间平滑MUSIC算法CF-SSMUSIC，本发明方法具有更好的波达方向DOA估计测角性能。以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本
技术领域：
的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。当前第1页1 2 3