针对这一问题，文献[10-11]提出了利用距离-多普勒二维相关后目标距离-多普勒位置的复幅度信息进行角度估计。这种方法虽然通过二维互相关处理提高了目标回波的信噪比，但接收信号快拍数目降为1，传统的超分辨算法不再适用，这对外辐射源雷达目标DOA估计带来了新的挑战。近年来，基于压缩感知(Compressed Sensing，CS)理论的DOA估计方法被应用于外辐射源雷达[12-13]。这种方法依靠构造过完备原子字典，通过迭代求解选择最匹配的原子，能够在单快拍的情况下实现多目标的DOA估计。但在过完备原子字典的构建中，字典矩阵的分割精度影响着DOA估计的精度，若分割的精度提高，则过完备字典矩阵的维数将变大，导致运算量较大。若减小字典矩阵的分割精度存在网格失配的问题，会影响目标角度的估计性能。
针对上述问题，本文将A&M插值迭代算法引入外辐射源雷达DOA估计，将角度估计问题转化为频率估计问题。同时，针对A&M插值迭代算法分辨率受阵元数目及阵列结构的限制，将稀疏阵列引入外辐射源雷达，提出了一种新的SA&M-Relax外辐射源雷达目标DOA估计方法，在减少阵元数目的情况下保证了估计精度。通过仿真实验，验证了本文所提方法的有效性。
1 系统模型 外辐射源雷达监测天线不仅包含了目标回波，还包含了直达波信号、多径杂波信号及噪声信号，则第m个监测通道接收信号可以表示为
 | (1) |
式中：Sd为直达波信号复幅度；m=1, 2, …, M为天线个数；t=1, 2, …, T为采样时间, T为信号的长度；Nm(t)为天线接收到的噪声信号；Aq、τq和fq分别为目标回波信号的复幅度、延时和多普勒频率；Bi和τi分别为监测天线接收到的多径杂波信号的复幅度和延时；q和i分别为目标数目和多径杂波数目；am为信号角度复幅度信息；θq、θi和θd分别为目标回波信号、多径杂波信号和直达波信号的入射角度。
a(θ)为目标接收信号的阵列流形矢量，其表示为
 | (2) |
式中：dm-1为第m个阵元到第1个阵元的距离；λ为信号波长。
由于直达波信号和多径杂波信号的能量强于目标回波信号，通过时域对消处理抑制直达波信号和多径杂波信号得到目标回波信号[14]。各阵元对消后的目标回波信号可以表示为
 | (3) |
由于对消后的目标回波信号能量极低，噪声强度远高于目标回波信号，利用对消后的信号直接进行目标角度估计难以获得角度信息。为了提取出目标的角度信息，利用距离-多普勒二维相关后目标距离-多普勒位置的复幅度信息进行角度估计，此时目标回波信噪比显著增强。目标回波信号和参考信号的互相关函数可以表示为
 | (4) |
式中：Sd(t)为监测通道提纯后的直达波信号；l=0, 1, …, L和k=-K, …, K分别为时延和多普勒频率离散化表示，L和K分别为离散化时延和多普勒频率范围。
考虑一个积累单元内有q个目标同时位于一个距离-多普勒单元[lg, kg]内，则此时位于[lg, kg]内的q个目标的回波信号可以表示为
 | (5) |
式中：A表示目标导向矢量；Xq(lg, kg)表示第q个目标在[lg, kg]距离-多普勒单元的复包络；Nm(lg, kg)表示噪声信号。
通过积累可以保留目标的幅相信息，但得到的接收信号Y为M×1维矩阵。
2 基于SA&M-Relax算法的目标DOA估计方法 2.1 A&M插值迭代算法 A&M插值迭代算法[15]是采用单快拍求解目标波达方向的基础。传统FFT方法求解目标角度时，通过FFT得到其频谱图，根据其频谱图的峰值点位置求解目标的角度。但这种方法估计得到的目标角度性能依赖于天线数量，在天线数量较少时估计得到的角度与实际值有较大差别。对此，A&M插值迭代算法利用FFT粗估计结果，采用插值迭代的方法提升角度估计的精度。
对于阵元间隔为d的一个均匀线阵，当目标入射角度为θ时，接收信号的阵列流形矢量为
 | (6) |
令f=dMsin θ/λ，则对于目标角度的估计可以转化为频率估计问题。假设FFT的峰值点为
 | (7) |
式中：fs为采样频率，在DOA估计中，fs=M。
因此，求解f关键在于找出误差迭代关系：
 | (8) |
式中：
考虑接收信号的离散傅里叶变换系数为
 | (9) |
式中：
将正弦信号表达式代入式(9)得
 | (10) |
式中：Np为噪声的傅里叶变换系数。
当δ-p << M时，则式(10)简化为
 | (11) |
式中：
文献[15]给出了2种误差迭代公式，其迭代关系式为
 | (12) |
 | (13) |
将式(11)代入式(12)可以得到
 | (14) |
可以看出，
综上所述，A&M插值迭代算法的流程如下：
步骤1? 对接收数据做FFT，得到
步骤2? 计算初始频率
步骤3? 初始化δ0=0，i=0。
步骤4? i=i+1, 计算傅里叶系数：
 | (15) |
步骤5 ?通过式(12)或式(13)计算误差
步骤6? 迭代，
步骤7? 通过
2.2 基于稀疏阵列的多源A&M插值迭代算法 通过2.1节分析，A&M插值迭代算法是基于均匀线阵且目标个数为1得到的。由于稀疏阵列是非均匀的，无法直接采用A&M插值迭代算法进行角度估计。对于非均匀阵列，可以等效看作是由多个虚拟均匀阵元叠加而成。此处记为
 | (16) |
式中：
 | (17) |
式中：
采用最小二乘法，可以得到阵列转换矩阵：
 | (18) |
利用得到的转换矩阵，可以将二维相关积累后的接收信号Y转换到虚拟均匀阵列数据，转换后的数据为
 | (19) |
将转化后的信号
 | (20) |
当同一个距离-多普勒单元出现的目标个数大于1时，A&M插值迭代算法不再适用。针对此问题，结合Relax算法对多目标实现角度估计。在对多源信号中某一个信号源角度估计时，可以将其他信号源当作噪声处理。在估计其中一个目标信号的频率点时，减去其余目标在该点的分量(称此分量为泄漏分量)，再采用A&M插值迭代算法估计，通过相互迭代提升估计的精度，得到目标的估计结果。
设第q个目标对第l个目标的泄漏分量为αqΥq[μl+p]，αq为第q个目标的傅里叶系数，μl为离散傅里叶变换后第l个目标在频谱中的位置。减去其余目标对第l个目标点的泄漏后，可得
 | (21) |
 | (22) |
式中：
综上所述，本文基于SA&M-Relax算法的目标DOA估计方法步骤如下：
步骤1 ?通过式(19)将稀疏阵列接收信号转化到虚拟均匀线阵，得到转换后数据
步骤2? 对转换后的数据
步骤3? 初始化μl=0, αl=0, l=1, 2, …, Q。
步骤4? 通过式(21)计算得到
步骤5? 利用式(12)或式(13)计算hjl。
步骤6? 更新μjl，此处加入辛克函数以提高计算精度。
 | (23) |
步骤7 ?更新幅值。
 | (24) |
步骤8? 分别计算每个目标的αjl，并进行迭代，判断是否收敛，未收敛返回步骤4；若收敛计算结束。
2.3 计算复杂度分析 本文方法的计算复杂度分析如下：当阵列结构为均匀阵列，目标个数为1时，2.1节中算法总复杂度为O(K(n2+2n)+nlog n)；当阵列结构为稀疏阵列时，式(18)可以预先计算并存储，故一个目标算法的计算复杂度为O(K(n2+2n)+nlog n+n2)；多个目标时算法的运算复杂度为O(Kq(2n2+2n)+nlog n)；文献[12]中算法复杂度为O(K(5n3+2n))。K表示算法的循环次数，q为目标个数。通过比较可得，本文方法能够有效降低运算复杂度。
3 仿真分析 本节采用地面数字电视多媒体广播信号(Digital Terrestrial Television Multimedia Broadcasting signal，DTMB)[16]对本文方法进行仿真验证，信号带宽为7.56 MHz，波长为0.5 m，采样率设置为2倍的信号带宽。采用的帧体结构为PN420模式，在二维相干积累中以一个帧体时长作为脉冲重复周期，对应的脉冲重复频率为1.8 kHz。设置阵列结构为[0, 0.5λ, 1.5λ, 2λ, 3λ]，直达波信号及多径杂波信号参数如表 1所示。
表 1 实验仿真参数 Table 1 Experimental simulation parameters
| 信号 | 距离/m | 多普勒频率/Hz | 信噪比/dB |
| 直达波 | 0 | 0 | 50 |
| 多径1 | 25 | 0 | 35 |
| 多径2 | 45 | 0 | 30 |
| 多径3 | 65 | 0 | 28 |
| 多径4 | 85 | 0 | 24 |
| 多径5 | 105 | 0 | 20 |
表选项
设定目标速度为20 m/s，距离为2.8 km，设积累时间为1 s，目标信噪比为-20 dB，目标角度在每个实验中分别给出。假设每个阵元接收信号中的直达波与多径杂波已经采用对消算法抑制[14]。通过距离-多普勒二维相关处理，已经获得了目标的距离和速度信息，此时距离-多普勒二维相关结果如图 1所示。可以看出，目标的峰值已经出现。在多普勒频率为0的位置出现的伪峰是由于DTMB信号PN帧头自相关特性而引起的。
 |
| 图 1 目标检测结果 Fig. 1 Target detection results |
| 图选项 |
1) 实验1。验证本文方法的收敛特性。设置稀疏线阵转换到均匀线阵的虚拟阵元数目为21，均匀线阵间距为0.5λ。此时目标对应的频率通过虚拟后的阵元数目与阵元间距来计算。计算表达式为
 | (25) |
单个目标时，设置目标角度为18.7°，通过式(25)计算得到对应的频率值为3.37 Hz。迭代次数设置40次，仿真实验结果如图 2所示。可以看出，仿真大约在第12次开始收敛，且收敛的结果与目标实际的频率基本一致。在实际仿真实验中，迭代次数一般不会超过20。
 |
| 图 2 单个目标时本文方法收敛仿真 Fig. 2 Convergence simulation diagram of proposed method with one target |
| 图选项 |
在实际环境中，多个目标落入同一个距离-多普勒单元的个数一般不超过2个，故本文选取了2个目标进行仿真验证。设置2个目标角度分别为26.6°与3.7°，通过式(25)计算得到对应的频率值分别为4.7 Hz与0.68 Hz。迭代次数设置100次，仿真实验结果如图 3所示。可以看出，仿真大约在第40次左右开始收敛，收敛的结果与目标实际的频率基本一致。在实际仿真实验中，如果迭代50次之后收敛的结果与实际结果仍然不一致，此时增加迭代次数也没有意义了。
 |
| 图 3 两个目标时本文方法收敛仿真 Fig. 3 Convergence simulation diagram of proposed method with two targets |
| 图选项 |
2) 实验2。比较相同阵元数目下的稀疏阵列和均匀阵列估计性能。稀疏阵列为本实验采用阵列，均匀线阵采用相同阵元数目、阵元间距为半波长。设定2个目标角度分别为3.7°与18.4°。在5阵元均匀线阵的情况下, 通过式(25)计算得到对应的频率值分别为0.16 Hz和0.79 Hz，在5阵元稀疏阵列结构下转换后的虚拟阵元个数为21，此时通过式(25)计算得到对应的频率值分别为0.68 Hz和3.31 Hz。图 4为2种阵列的最终收敛结果。其中，S表示稀疏阵列，A表示实际值，E表示估计值，数字代表第几个目标。可以看出，在相同阵元数目的情况下，均匀线阵不能成功分辨2个目标，但稀疏阵列能够成功分辨出2个目标，并估计出目标角度。因此，SA&M-Relax算法相比于A&M-Relax具有更高的分辨率。
 |
| 图 4 不同阵列估计结果 Fig. 4 Different array estimation results |
| 图选项 |
3) 实验3。比较本文方法与压缩感知算法的分辨率性能。设定2个目标角度差不同，分别求2种算法正确估计的概率(估计结果与实际角度差为±1°)，每种情况分别进行100次蒙特卡罗仿真实验。仿真结果如图 5所示。可以看出，当2个目标角度差为12°左右时，本文方法能正确分辨2个目标概率大于0.9。在目标角度差小于12°之前，分辨性能略差于压缩感知算法，这是由于本文方法的分辨性能受到傅里叶分辨率的限制。但若已知2个目标中的一个目标角度，此时本文方法的分辨率将显著提高。图 6为2个目标角度相差5°的情况下，已知其中一个目标的角度，将其作为迭代的初始值，得到的最终迭代结果。可以看出，在已知其中一个目标的角度情况下，本文方法的分辨率将显著提升，最终迭代结果会与目标实际频率基本一致。
 |
| 图 5 两个目标不同角度间隔下正确分辨概率 Fig. 5 Probability of correct resolution between two targets at different angle intervals |
| 图选项 |
 |
| 图 6 一个目标角度已知情况下估计结果 Fig. 6 Estimation results when a target angle is known |
| 图选项 |
4 结论 针对现有的外辐射源雷达目标波达方向估计存在单快拍下的低信噪比角度估计问题。
1) 提出了一种新的SA&M-Relax外辐射源雷达目标DOA估计方法。
2) 解决了A&M插值迭代算法只适用于均匀线阵和单个目标的角度估计问题，提升了估计性能。
3) 相比于压缩感知类算法，能够避免稀疏字典的构建，减小运算量，为外辐射源雷达目标DOA估计问题提供一个新的思路。
参考文献
| [1] | GRIFFITHS H D, BAKER C J. An introduction to passive radar[M]. Norwood: Artech House, 2017: 10. |
| [2] | WANG L, HE Q, BLUM R S, et al. Joint parameter estimation employing coherent passive MIMO radar[J]. The Journal of Engineering, 2019, 2019(20): 6859-6862. DOI:10.1049/joe.2019.0578 |
| [3] | DEBORA P, FABRIZIO S, FEDERICA P, et al. Maritime moving target long time integration for GNSS-based passive bistatic radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 2018, 54(6): 3060-3083. DOI:10.1109/TAES.2018.2840298 |
| [4] | MOHAMMAD F, MOSTAFA D, SHAPOUR K. Rao detector for passive MIMO radar with direct-path interference[J]. IEEE Transactions on Aerospace and Electronic Systems, 2020, 56(4): 2999-3009. DOI:10.1109/TAES.2019.2961841 |
| [5] | FENG Y, SHAN T, LIU S H, et al. Interference suppression using joint spatio-temporal domain filtering in passive radar[C]//2015 IEEE Radar Conference. Piscataway: IEEE Press, 2015: 1156-1160. |
| [6] | SENDALL J L, MAASDORP F D V. Detection state refinement in FM multistatic passive radar[C]//2017 IEEE Radar Conference. Piscataway: IEEE Press, 2017: 717-721. |
| [7] | MAUSSANG F, GARELLO R, SOULAT F, et al. GPS passive bistatic radar system in oceanic environment: Detection performance estimation[C]//2011 IEEE Spain OCEANS. Piscataway: IEEE Press, 2011: 1-5. |
| [8] | 郭帅, 王俊, 陈刚, 等. 外辐射源雷达参考通道多径干扰抑制方法[J]. 西安电子科技大学学报, 2018, 45(3): 18-23. GUO S, WANG J, CHEN G, et al. Method for multipath interference restraint of the reference channel in passive bistatic radar[J]. Journal of Xidian University, 2018, 45(3): 18-23. DOI:10.3969/j.issn.1001-2400.2018.03.004 (in Chinese) |
| [9] | FENG W K, FRIEDT J M, MEROU G G, et al. Passive radar delay and angle of arrival measurements of multiple acoustic delay lines used as passive sensors[J]. IEEE Sensors Journal, 2019, 19(2): 594-602. DOI:10.1109/JSEN.2018.2872867 |
| [10] | HOWLAND P E, MAKSIMIMUK D, REITSMA G. FM radio based bistatic radar[J]. IEE Proceedings-Radar, Sonar and Navigation, 2005, 152(3): 107-115. DOI:10.1049/ip-rsn:20045077 |
| [11] | WANG J, WANG H T, ZHAO Y. Direction finding in frequency-modulated-based passive bistatic radar with a four-element Adcock antenna array[J]. IET Radar, Sonar and Navigation, 2011, 5(8): 807-813. DOI:10.1049/iet-rsn.2010.0262 |
| [12] | 王海涛, 王俊. 基于压缩感知的无源雷达超分辨DOA估计[J]. 电子与信息学报, 2013, 35(4): 877-881. WANG H T, WANG J. Super-resolution DOA estimation in passive radar based on compressed sensing[J]. Journal of Electronics & Information Technology, 2013, 35(4): 877-881. (in Chinese) |
| [13] | 左罗, 王俊, 陈刚, 等. 基于TLS-CS的外辐射源雷达超分辨DOA估计方法[J]. 系统工程与电子技术, 2020, 42(1): 61-66. ZUO L, WANG J, CHEN G, et al. TLS-CS super-resolution DOA estimation method in PBR[J]. Systems Engineering and Electronics, 2020, 42(1): 61-66. (in Chinese) |
| [14] | 张各各, 王俊, 刘玉春. 一种基于数字电视地面广播照射源的外辐射源雷达快速杂波相消算法[J]. 电子与信息学报, 2013, 35(1): 36-40. ZHANG G G, WANG J, LIU Y C. A fast cancellation algorithm for DTTB-based passive radar system[J]. Journal of Electronics & Information Technology, 2013, 35(1): 36-40. (in Chinese) |
| [15] | ABOUTANIOS E, MULGREW B. Iterative frequency estimation by interpolation on Fourier coefficients[J]. IEEE Transactions on Signal Processing, 2005, 53(4): 1237-1242. DOI:10.1109/TSP.2005.843719 |
| [16] | 数字电视地面传输国家标准特别工作组. 数字电视地面广播传输系统帧结构、信道编码和调制: GB 20600-2006[S]. 北京: 中国标准出版社, 2007. National Standards Task Force for Terrestrial Digital Television Transmission. Frame structure, channel coding and modulation for digital television terrestrial broadcasting system: GB 20600-2006[S]. Beijing: China Standard Press, 2007(in Chinese). |
