图1 外辐射源探测场景
Fig.1 The observation scene of passive radar with opportunistic illuminator
Reference format: Song Kekang, Feng Wentao. An Efficient Method of Direct Position Determination of Passive Radar with DIRECT Algorithm[J]. Journal of Signal Processing, 2020, 36(1): 149-154. DOI: 10.16798/j.issn.1003- 0530.2020.01.018.
外辐射源雷达通过接收第三方辐射源(比如广播、数字电视等)照射到目标后的散射信号来实现对目标的检测和定位等,具有隐蔽性好、成本低等优势,受到国内学者的广泛研究[1-2]。在外辐射源雷达的目标定位中常常采用的是先估计中间参数,比如到达时间(Time of arrival, TOA)、到达角度(Direction of arrival, DOA)、到达时差(Time difference of arrival, TDOA)或者到达频差(Frequency difference of arrival, FDOA)以及这些参数的组合等,再通过中间参数的解算实现目标的定位,这类方法也称作两步定位方法[3- 4]。针对两步定位法存在信息损失,且在低信噪比时定位精度低的问题,Weiss等人[3]提出目标直接定位(Direct position determination, DPD)方法,该方法将中间参数直接表示为目标位置和速度的形式,只需一步定位计算便可从接收信号中获取目标位置。Weiss等人经过大量的研究和仿真证明[5- 6],DPD方法得益于信息损失少,在低信噪比时定位精度优于传统的两步定位法。但较为遗憾的是,DPD方法中关于位置的目标函数不是凸函数[5],使用传统的凸优化方法较难获得全局最优值,因此网格遍历搜索成为较多文献中常用的方法[6-7],而遍历搜索的效率较低难以满足实际应用。针对该难题,文献[8]通过分布式计算的思路提升网格遍历搜索的效率,但其思路仍是利用网格遍历搜索,因此计算效率没有到较大提升。文献[9]提出利用期望最大值(Expectation maximization, EM)方法迭代求解目标位置的思路,但该方法局限性于TOA体制且需已知发射信号。文献[10]根据矩阵扰动理论,利用牛顿迭代方法对目标位置进行估计,但牛顿迭代难以保证全局最优且存在Hessian矩阵不可逆的问题。还有学者采用遗传优化算法等智能优化算法对目标函数进行全局搜索[11],该类算法仍然面临高维搜索问题。
因此,针对直接定位方法计算效率不高难以满足实际应用需求的问题,本文以一发多收的外辐射源雷达对目标进行直接定位为背景,在最大似然准则下,理论推导了外辐射源雷达含直达波的直接定位目标函数,并提出基于DIRECT(Dividing rectangles)算法实现的直接定位快速计算方法(DIRECT algorithm position determination, DIRECT-DPD)。该方法具备全局收敛性,并且能够实现目标位置的快速估计,相比已有的网格遍历方法计算时间降低2个数量级,相比遗传算法降低1个数量级。
建立一发多收模式的外辐射源雷达探测场景,如图1所示,对静止或慢速目标(短时间内可认为是静止的)进行探测,比如固定目标或慢速海面目标。图1中外辐射源雷达接收站分为回波通道和直达波通道,并采用物理分离的方式进行隔离,实际中也可利用直达波抑制技术[12]等数字信号处理方法来实现,其中回波通道用于接收第三方辐射源发射信号经目标散射后的回波信号,直达波通道用于接收发射信号的直达波。为便于分析,本文假设回波通道与直达波通道位于相同地理位置,实际中也可在不同位置,不影响本文结论。记目标位置矢量为p,第三方辐射源位置矢量为t,第m个接收站位置矢量为rm,其中m=1,…,M,M为接收站个数。考虑窄带连续波信号,则第m个接收站直达波和回波的基带离散信号可分别表示为
图1 外辐射源探测场景
Fig.1 The observation scene of passive radar with opportunistic illuminator
(1)
(2)
式(1)和(2)中,n=1,…,N,N为信号点数,下标d用于表示直达波信号,下标p用于表示回波信号,
表示直达波路径的复衰减系数,
表示回波路径的复衰减系数(包括目标散射系数等),u(n)表示n时刻的发射信号,
其中
和
分别表示第m个接收站的直达波时延和回波时延,Ts=1/fs为采样时间间隔,
和
分别为第m个接收站直达波和回波通道的噪声,且分别服从均值为零,方差为
和
的复高斯分布。进一步,将接收信号写成矢量形式,并令
类似有
和
则第m个接收站的直达波和回波信号矢量可分别表示为
(3)
(4)
式(3)和(4)中,
和
为时延矩阵[13],其中
为单位离散DFT(Discrete Fourier transform)矩阵,并且有H(y)=diag([ej2π(0)/Nsy,…,ej2π(Ns-1)/Nsy])。若将第m个接收站的直达波和回波信号写成扩展列矢量,则有
(5)
式中,m=1,…,M。
本节在第2节外辐射源雷达信号模型基础上,先分析最大似然直接定位原理,给出外辐射源雷达直接定位的目标函数,针对该目标函数为非凸函数,传统凸优化方法难以获得全局最优解的问题,提出采用DIRECT算法对目标位置进行快速求解的方法。
根据式(5)得,所有接收站直达波与回波扩展信号矢量的似然函数为
(6)
式中
表示块对角矩阵,IN表示N×N的单位矩阵,
表示未知参数矢量。由式(6)可得目标位置的最大对数似然估计为
(7)
式中,Ω为目标可能出现的位置区域。对式(6)中的未知参数
和
进行最大似然估计,m=1,…,M,可得估计结果分别为
(8)
(9)
将式(8)和式(9)代入式(7)中,可进一步得到目标位置的最大似然估计为
(10)
根据瑞利商公式[14],可知当u为式(10)中求和项的最大特征值对应的特征值向量时,式(10)中的目标函数取得最大值,此时有
(11)
式中,G(p)=U(p)HU(p),且U(p)的表达式为
U(p)=[Ud,Up(p)]
(12)
式中,
其中Up(p)为关于p的表达式是因为
中包含回波时延,是关于目标位置的表示式,m=1,…,M。此外,式(11)的成立还利用了矩阵特征值的性质,即λmax(UUH)=λmax(UHU)。进一步将G(p)展开可得
(13)
分析式(13)可知,
为各站直达波与直达波互相关部分,不含目标信息,
为直达波与回波互相关部分,对应于直达波与回波的匹配滤波,由其非对角元素可获取回波到达时延,
为回波与回波互相关部分,由其非对角元素可获取不同接收站回波之间的时差。因此,矩阵G(p)既包含直达波与回波之间的匹配滤波也包含回波与回波之间的互相关部分,既可获取到达时延也能获取到达时差,这表明目标函数充分利用了直达波和回波信息。结合文献[5-6,8]知,矩阵G(p)最大特征值关于目标位置p的函数也为非凸函数,采用传统的凸优化方法难以获得其全局最优值,而常用的网格遍历方法运算量大[6-7],对此本文在下一小节中提出采用DIRECT方法对其进行高效求解。
DIRECT算法是求解全局最小值的经典算法,已被证明具有全局收敛性[15]。DIRECT算法采用的策略是计算区间中点而不是区间端点的函数值,因此其适应性广,且适合高维时的运算。DIRECT算法的主要思想是,先将待搜索的可行域变换为单位超立方体,然后计算超立方体中心的函数值,接着超立方体被三等分为多个超矩形,并计算超矩形中心点函数值,根据这些超矩形和函数值大小,寻找潜最优超矩形(具体参考式(14)的定义),并对潜最优超矩形进行划分,然后在划分后的所有超矩形中继续寻找潜最优超矩形直至满足条件(设定精度或步长等)为止。其搜索过程示例如图2所示,图中黑点表示超矩形中心,第1列为开始迭代的状态,第2列为潜最优超矩形的选择,其中阴影部分对应潜最优超矩形,第3列为潜最优超矩形划分后的状态。
图2 DIRECT算法搜索示例[16]
Fig.2 The example of search process of DIRECT algorithm[16]
如何寻找潜最优超矩形是DIRECT算法的关键,其定义如下[15]:将单位超立方体划分为I个超矩形,令ci为第i个超矩形的几何中心,di表示该中心到各顶点的距离,ε为正常数,若对于第h个超矩形,存在
使得
(14)
同时成立,则称该超矩形为潜最优超矩形。式(14)中,g(·)为直接定位目标函数,表达式为
g(p)=-λmax(G(p)),p∈Ω
(15)
式中存在负号是由于DIRECT算法求解的是全局最小值。潜最优超矩形的寻找方法可以通过Graham的conhull算法[16]高效实现。
根据3.2节中DIRECT算法原理的描述,可将DIRECT-DPD的算法步骤总结如下:
步骤1 将搜索位置空间Ω标准化为单位超立方体,将超立方体的中心表示为c1,计算c1处目标函数值g(c1),同时令gmin=g(c1),并记q=1,k=0(k为迭代次数);
步骤2 寻找潜最优超矩形的集合S;
步骤3 选择任一潜最优超矩形s∈S,确定超矩形s的最大边长对应维数的集合I,比如若I={1,2}则表示超矩形的第1和第2维边长均为最大边长,令ξ为最大边长的1/3;
步骤4 计算点cs±ξei处的目标函数值,其中cs为超矩形s的中心,ei为第i维方向上的单位矢量,即除第i个元素为1外,其余元素均为零,i∈I;
步骤5 计算wi=min{g(cs-ξei),g(cs+ξei)},i∈I;
步骤6 对I中所有维度方向的超矩形进行三等分,其顺序按照wi从小到大进行,比如若I={1,2}且w1<w2,则先对第1维上的超矩形进行三等分,然后再对第2维上的超矩形进行三等分;
步骤7 根据步骤5更新gmin,并令q=q+Δq,其中Δq为新采样点数;
步骤8 令S=S-{s},如果S≠∅,返回步骤3,否则进行下一步;
步骤9 令k=k+1,若k等于设定迭代次数,则停止;否则返回步骤2。
本小节分析文中方法与基于网格搜索[5-6]和遗传算法[11]的直接定位方法的理论计算量。仅考虑复数乘法的次数,且DFT采用基2FFT实现[17],其运算量为
通过式(14)中目标函数知,计算量主要集中在最大特征值的计算,记其计算量为E,则
记网格遍历法中每一维的网格划分数量为L,维数为D,则网格遍历法的计算量为LDE。遗传算法计算量为[11]PsRepochmaxE,其中Ps为种群数量,Repochmax为最大繁衍代数。本文算法只计算超矩形中心点的函数值,将超矩形的计算个数记为Q,则本文算法的计算量为QE。三种算法计算量对比如表1所示。
表1 不同方法计算量比较
Tab.1 The comparison of computation of different methods
算法计算量基于网格遍历的DPD方法LDE基于遗传算法的DPD方法PsRepochmaxE本文方法QE
采用4个外辐射源雷达接收站对目标进行定位,考虑二维平面内,其中接收站位置坐标分别为[-20, 10]km、[-10,-20]km、[20, 20]km和[30,-10]km,发射站位置坐标为[-5, 30]km,第三方辐射源信号为16QAM信号,信号带宽为100 kHz,复采样率为100 kHz,直达波的信噪比为20 dB。目标可能出现区域为[-10,-10]km 到[10, 10]km的方形区域,对应的搜索范围为20 km×20 km。目标位置等其余参数在仿真分析中具体给出。
仿真实验1 DIRECT-DPD方法最优超矩形中心点迭代轨迹和目标函数计算次数。在目标位置为[1, 4]km,信噪比为-15 dB时,图3(a)、(b)和(c)中分别给出迭代次数为1次、20次和50次时DIRECT-DPD方法最优超矩形中心的迭代轨迹。从图3中结果可以看出本文方法能够较快的收敛到目标真实位置附近。表2中给出不同信噪比条件下,蒙特卡罗实验200次时DIRECT-DPD方法中目标函数的平均计算次数,以便于计算量的定量分析。
图3 DIRECT-DPD方法中最优超矩形中心点迭代轨迹
Fig.3 The trace of iteration of the optimal hyper rectangle center in DIRECT-DPD method
仿真实验2 不同方法计算时间对比。选取信号点数为212,网格遍历法中网格大小设置为50 m,对应网格点数为400×400,遗传算法中种群大小和最大繁衍代数分别为40和100[11],根据表2中结果,本文DIRECT-DPD方法中超矩形的计算次数选为390次。
表2 不同回波信噪比时目标函数计算次数
Tab.2 The computation times of object function with different SNR
回波信噪比/dB目标函数计算次数/次回波信噪比/dB目标函数计算次数/次-25309.2-15361.3-20344.8-10383.8
根据仿真实验2中参数,结合表1中计算表达式,得到网格遍历方法的计算量为1010量级,遗传算法的计算量为109量级,而本文方法的计算量在108量级。为进一步验证计算复杂度分析的正确性,在相同仿真平台下,图4中给出信噪比-30 dB到-10 dB时,200次蒙特卡罗实验下,三种方法平均计算时间的对比曲线,其中DPD表示原始网格搜索的直接定位方法,GA-DPD表示基于遗传算法实现的直接定位方法。可以看出,图4中结果与表1中计算量分析的结论基本一致,表明了计算复杂度分析的合理性。因此,理论分析和仿真研究表明,本文方法计算效率最高,相比网格遍历提升2个数量级,相比遗传算法提升1个数量级。需要说明的是,本文网格遍历搜索的计算量与文献[11]中网格遍历搜索的计算量不同,原因在于文献[11]中的网格搜索范围为1 km×1 km,网格划分大小为1 m×1 m,而本文网格搜索范围为20 km×20 km,网格划分大小为50 m×50 m。
图4 不同方法平均计算时间随信噪比变化曲线
Fig.4 The average computation time of different method vs. SNR
仿真实验3 不同方法位置估计的均方根误差对比。在信噪比变化范围为-30 dB到-10 dB,蒙特卡罗仿真500次,且其余参数同仿真实验2相同的条件下,图5给出目标位置估计均方根误差随信噪比变化曲线,图中DPD与GA-DPD对应的算法与图4相同。从图中结果可以看出,本文方法定位精度与网格遍历法基本相当,优于遗传算法。因此结合图4和图5中的结果,不难发现在相同条件下,本文方法具有计算速度快且定位精度高的优势。
图5 不同方法位置估计均方根误差随信噪比变化
Fig.5 The root mean squared error of different method vs. SNR
直接定位无需估计中间参数,已被证明在低信噪比时定位精度优于传统两步定位方法。将直接定位用于外辐射源雷达含直达波的探测模型中,理论推导外辐射源雷达最大似然直接定位的目标函数,该目标函数也是Hermit矩阵的最大特征值,且为非凸函数,存在难以快速获取全局最大值的问题。针对该难题,采用DIRECT算法求解非凸目标函数的全局最优解,并对其计算复杂度进行了理论和数值分析。仿真实验结果验证了本文方法的合理性和高效性,具体地,相比常用的网格遍历算法,计算时间降低2个数量级,相比遗传算法计算时间降低1个数量级,同时定位精度在低信噪比时与网格遍历法相当。将本文算法用于多目标直接定位是下一步需要研究的内容。
[1] Batu K C, Braham H. GLRT detector in single frequency multi-static passive radar systems[J]. Signal Processing, 2018, 142: 504-512.
[2] Song Kekang, Yang Yuxiang, Peng Huafeng. Direct Location for Multiple Passive Radars without and with Reference[C]∥Proceedings of 2018 3rd International Conference on Communication, Image and Signal Processing. Sanya: 2018.
[3] 邓兵, 孙正波, 杨乐. 一种简单有效的TDOA-FDOA-AOA目标定位闭式解[J]. 西安电子科技大学学报, 2018, 45(2): 171-180.
Deng Bing, Sun Zhengbo, Yang Le. Efficient Closed-form Estimator for TDOA-FDOA-AOA Localization[J]. Journal of Xidian University, 2018, 45(2): 171-180.(in Chinese)
[4] 曹景敏, 万群, 欧阳鑫信, 等. 观测站有位置误差的多维标度时频差定位算法[J]. 信号处理, 2017, 33(1): 1-9.
Cao Jingming, Wan Qun, Ouyang Xinxin, et al. Multidimensional Scaling-based Passive Emitter Localization from Time Difference of Arrival and Frequency Difference of Arrival Measurements with Sensor Location Uncertainties[J]. Journal of Signal Processing, 2017, 33(1): 1-9.(in Chinese)
[5] Weiss A J. Direct position determination of narrowband radio frequency transmitters[J]. IEEE Signal Processing Letters, 2004, 11(5): 513-516.
[6] Weiss A J. Direct Geolocation of Wideband Emitters[J]. IEEE Transaction on Signal Processing, 2011, 59(6): 2513-2521.
[7] 邓丽娟, 张展, 魏平, 等. 分布式MIMO雷达中仅使用多普勒频移的直接定位技术[J]. 信号处理, 2018, 34(11): 1377-1384.
Deng Lijuan, Zhang Zhan, Wei Ping, et al. Direct Position Determination Using Merely Doppler Frequency Shifts in Distributed MIMO Radar[J]. Journal of Signal Processing, 2018, 34(11): 1377-1384.(in Chinese)
[8] Pourhomayoun M. Distributed Computation for Direct Position Determination Emitter Location[J]. IEEE Transactions on Aerospace and Electronic Systems, 2014, 50(4): 2878-2889.
[9] Tzoreff E, Weiss A J. Expectation-maximization Algorithm for Direct Position Determination[J]. Signal Processing, 2017, 133: 32-39.
[10] 王鼎, 张刚. 一种基于窄带信号多普勒频率测量的运动目标直接定位方法[J]. 电子学报, 2017, 45(3): 591-598.
Wang Ding, Zhang Gang. A Direct Localization Method for Moving Narrowband Source Based on Doppler Frequency Shifts[J]. ACTA Electronica Sinica, 2017, 45(3): 591-598.(in Chinese)
[11] 任衍青, 逯志宇, 巴斌, 等. 锐化遗传直接定位快速估计算法[J]. 西安电子科技大学学报, 2017, 44(4): 144-150.
Ren Yanqing, Lu Zhiyu, Ba Bin, et al. Fast Direct Position Determination Method Based on the Sharpening Function Genetic Algorithm[J]. Journal of Xidian University, 2017, 44(4): 144-150.(in Chinese)
[12] 邵启红. 低频段外辐射源雷达信号处理若干关键技术研究[D]. 武汉: 武汉大学, 2016.
Shao Qihong. Key Techniques of Low Frequency Band Passive Radar Signal Processing[D]. Wuhan: Wuhan University, 2016.(in Chinese)
[13] Hack D E, Patton L K, Braham Himed, et al. Centralized Passive MIMO Radar Detection without Direct-Path Reference Signals[J]. IEEE Transactions on Signal Processing, 2014, 62(11): 3013-3023.
[14] 张贤达. 矩阵分析与应用[M]. 北京: 清华大学出版社, 2013: 57, 442- 443.
Zhang Xianda. Matrix analysis and applications[M]. Beijing: Tsinghua Press, 2013: 57, 442- 443.(in Chinese)
[15] Finkel D E, Kelley C T. Convergence Analysis of the DIRECT Algorithm[J]. AMS, 2004.
[16] Jones D R, Perttunen C D, Stuckman B E. Lipschitzian Optimization Without the Lipschitz Constant[J]. Journal of Optimization Theory and Application, 1993, 79(1): 157-181.
[17] 吴瑛, 张莉, 张冬玲, 等. 数字信号处理[M]. 西安: 西安电子科技大学出版社, 2009: 124-128.
Wu Ying, Zhang Li, Zhang Dongling, et al. Digital Signal Processing[M]. Xi’an: Xidian University Press, 2009: 124-128.(in Chinese)
宋科康 男, 1990年生, 四川自贡人。盲信号处理国家级重点实验室博士研究生, 主要研究方向为微弱目标检测与定位跟踪。
E-mail: skkybyq@163.com
冯文涛 男, 1982年生, 湖北钟祥人。盲信号处理国家级重点实验室助理研究员, 主要研究方向为数据融合、传感器管理和优化算法。
E-mail: wentaofeng@163.com