②(澳大利亚伍伦贡大学电气计算机与通信工程学院 伍伦贡 999029)
②(School of Electrical, Computer and Telecommunications Engineering, University of Wollongong, Wollongong 999029, Australia)
在现代战争中,随着导弹、歼击机等空中目标的机动性越来越强,雷达检测与追踪它们的难度也日益增加,所以提高雷达的性能就成为了夺取空战胜利中不可或缺的一环。而机载预警雷达将平台架设到飞机上,具有机动性强、覆盖范围广等优点[1],既提高了雷达在战场环境的生存能力,也弥补了地面雷达在探测距离上和检测低空飞行目标上的不足。但机载预警雷达的下视工作状态决定了其较地基雷达面临更加难以处理的地杂波问题,这使得它在目标检测和参数估计问题的处理上变得尤为困难[2]。空时自适应处理(Space-Time Adaptive Processing,STAP)是一种较为可行的机载雷达地杂波抑制方法,但STAP方法在处理雷达目标回波时,假设回波数据在相干处理时间(Coherent Processing Interval,CPI)内的多普勒频率不变。而机载预警雷达在检测机动目标时,其回波会产生多普勒走动[2](其多普勒频率在CPI内发生变化),这样也使得原有的STAP方法对目标能量的相参积累性能大大下降,从而导致雷达对机动目标的检测性能下降[3, 4]。因此研究机动目标检测和参数估计具有重要意义。
当机动目标在一个CPI内做匀加速运动时,目标回波可以看作线性调频(Linear Frequency Modulation,LFM)信号[5]。目前LFM信号检测和参数估计方法主要有最大似然估计[6](Maximum Likelihood,ML)和时频分析[7]等方法。虽然最大似然估计精度很高,对参数的估计方差接近于Cramer-Rao下界,但其过大的运算量不利于该方法在实际情况中进行实时处理和工程实现。而时频分析方法要求有足够的脉冲采样点数,否则目标参数的估计精度难以满足要求[8, 9, 10, 11, 12]。当脉冲重复频率恒定时,脉冲数的增加会使得CPI变大,此时目标和杂波的回波信号可能会产生距离走动,使得后续的处理更加困难。文献[9]利用Radon-傅里叶变换(Radon-Fourier Transform,RFT)校正目标回波的距离走动,但仅适于匀速运动目标。文献[10]利用Radon-分数阶傅里叶变换(Radon-Fractional Fourier Transform,RFRFT)校正目标回波的距离走动和多普勒走动,但在多天线情况下无法直接应用。文献[11]在单天线体制下利用Radon-分数阶模糊函数(Radon-Fractional Ambiguity Function,RFRAF)技术能够有效地校正高机动目标(目标加速度变化,即存在急动度)的距离走动和多普勒走动,同样不能直接应用于多天线体制。另外文献[12]提出的基于重构时间采样的机动目标检测和参数估计算法能够在脉冲数较少的情况下获得精度较高的参数估计,但该方法需要目标方位先验信息已知,当方位信息未知时,该方法无法直接应用。
本文提出了一种基于WVD和重构时间采样的空中机动目标检测和参数估计方法,该方法能够在脉冲数目较少,方位信息未知的情况下,利用重构时间采样和WVD方法对目标回波信号进行有效积累,然后采用变步长3维搜索实现对方位未知机动目标的参数估计。仿真结果验证了所提方法的有效性。
2 数据模型设机载预警雷达为N元均匀线阵,雷达发射脉冲波长为λ,阵元之间的间隔为d=0.5λ。雷达在一个相干处理间隔内发射K个脉冲,xnk 为第n个阵元在第k个脉冲上的采样值,则在一个相干处理间隔内的待检测单元回波数据是一个N×K的矩阵[12]:
$X = \left[ {\begin{array}{*{20}{c}} {{x_{{\rm{11}}}}}&{{x_{{\rm{12}}}}}& \cdots &{{x_{{\rm{1}}K}}}\\ {{x_{{\rm{21}}}}}&{{x_{{\rm{22}}}}}& \cdots &{{x_{{\rm{2}}K}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{x_{N{\rm{1}}}}}&{{x_{N{\rm{2}}}}}& \cdots &{{x_{NK}}} \end{array}} \right]$ | (1) |
把X中的元素拉成NK×1的列向量,记作x=vec(X),其中${\mathop{\rm vec}\nolimits} ( * )$表示矩阵向量化操作算子,这样就得到了一个空时快拍数据,如式(2):
$ \begin{aligned} {x} = & \Big[{{x_{11}} \; {x_{12}} \; \cdots \; {x_{1K}} \; {x_{21}} \; {x_{22}} \; \cdots \; {x_{2K}} \; \cdots } \\ & { {{x_{N1}} \; {x_{N2}} \; \cdots \; {x_{NK}}} \Big]^{\rm{T}}} \end{aligned}$ | (2) |
当待检测单元内只有1个目标时,该单元的一次快拍数据为:
$x = {x_{\rm{s}}} + {x_{\rm{c}}} + {x_{\rm{n}}}$ | (3) |
其中,xs,xc和xn分别表示目标、杂波和噪声的采样数据,其中xs可以表示为:
${x_{\rm{s}}} = {a_{\rm{s}}}s({\varphi _{\rm{s}}},{\varphi _{\rm{t}}}) = {a_{\rm{s}}}s({\varphi _{\rm{s}}}) \otimes {\rm{ }}s({\varphi _{\rm{t}}})$ | (4) |
其中$ \otimes $表示Kronecker积,as表示目标回波复幅度,$s({\varphi _{\rm{s}}},{\varphi _{\rm{t}}})$为目标空时导向矢量,$s({\varphi _{\rm{s}}})$为空间导向矢量,$s({\varphi _{\rm{t}}})$为时间导向矢量[10],可分别表示为:
$\begin{gathered} s({\varphi _{\text{s}}}) = {[1{{\text{e}}^{{\text{j}}{\varphi _{\text{s}}}}}{{\text{e}}^{{\text{j}}2{\varphi _{\text{s}}}}} \cdots {{\text{e}}^{{\text{j}}(N - 1){\varphi _{\text{s}}}}}]^{\text{T}}} \hfill \\ \qquad \qquad \qquad \qquad \quad= [1\quad {{\text{e}}^{{\text{j}}2\pi d\cos {\theta _{\text{s}}}\cos {\gamma _{\text{s}}}/\lambda }}\quad {{\text{e}}^{{\text{j}}{2^2}\pi d\cos {\theta _{\text{s}}}\cos {\gamma _{\text{s}}}/\lambda }} \hfill \\ \qquad \quad \cdots {{\text{e}}^{{\text{j}}2(N - 1)\pi d\cos {\theta _{\text{s}}}\cos {\gamma _{\text{s}}}/\lambda }}{]^{\text{T}}}, \hfill \\ \qquad \qquad [s({\varphi _{\text{t}}}) = {[\begin{array}{*{20}{c}} 1&{{{\text{e}}^{{\text{j}}{\varphi _{\text{t}}}}}}&{{{\text{e}}^{{\text{j}}2{\varphi _{\text{t}}}}}}& \cdots &{{{\text{e}}^{{\text{j}}(K - 1){\varphi _{\text{t}}}}}} \end{array}]^{\text{T}}} \hfill \\ \qquad \qquad \qquad \qquad \qquad \qquad= 1[1{e^{ - j(\frac{{4\pi v}}{{\lambda fp}} + \frac{{2\pi aT}}{{\lambda fp}})cos{\theta _s}\cos {\gamma _s}}}{e^{ - j(2\frac{{4\pi v}}{{\lambda fp}} + \frac{{2\pi aT}}{{\lambda fp}})cos{\theta _s}\cos {\gamma _s}}} \hfill \\ \qquad \qquad \qquad \qquad \cdots {e^{ - j((K - 1)\frac{{4\pi v}}{{\lambda fp}} + {{(K - 1)}^2}\frac{{2\pi aT}}{{\lambda fp}})cos{\theta _s}\cos {\gamma _s}}}{]^T} \hfill \\ \end{gathered} $ | (5) |
其中${\varphi _{\rm{s}}} = 2\pi d\cos {\theta _{\rm{s}}}\cos {\gamma _{\rm{s}}}$表示归一化的空间角频率,${\theta _{\rm{s}}}$表示目标的方位角,${\gamma _{\rm{s}}}$表示目标的俯仰角,由式(5)可以看出,匀速目标(即a=0时)和机动目标的归一化空间角频率和空间导向矢量完全相同,而归一化时间角频率和时间导向矢量并不相同。机动目标的归一化时间角频率${\varphi _{\rm{t}}}$由两部分组成:第1部分为目标初速度对应的归一化时间角频率;第2部分为目标加速度对应的归一化时间角频率,会随时间发生改变,为多普勒走动项[12]。
3 基于WVD的机动目标检测和参数估计 3.1 WVD和重构时间采样Wigner-Ville分布是一种最基本、也是应用最广的时频分布之一[13]。信号x(t)的Wigner-Ville分布为:
${W_{\rm{x}}}(t,\Omega ) = \int_{ - \infty }^{ + \infty } {x\left( {t + \frac{T}{2}} \right){x^*}\left( {t - \frac{T}{2}} \right)\exp ( - {\rm{j}}\Omega T){\rm{d}}{\rm T}} $ | (6) |
式(6)可以看作是信号x(t)的自相关函数${R_{\rm{x}}}(t,\tau ) = x \left(t + \frac{\tau }{2} \right){x^*} \left( t - \frac{\tau}{2} \right)$对t的傅里叶变换。LFM信号时频分布的能量集中在其瞬时频率周围,故WVD对LFM信号有着最理想的时频聚集性[14]。
重构时间采样是利用空间采样重构时间采样,对各个阵元杂波抑制后的回波数据进行相位补偿后首尾拼接,这相当于增加了单个阵元的脉冲点数进行提高参数估计精度[12, 15]。图 1为4个阵元补偿相位前后时频图,图 1(a)为拼接前的时频效果图,图 1(b)为拼接后的时频效果图。由图中看出,只要对每个阵元补偿其相应的相位,就可以达到增加单个阵元脉冲点数的效果。
在此为简化推导过程,下面以2个阵元为例讨论将多阵元首尾拼接时每个阵元所需补偿的相位,这两个阵元杂波抑制后的回波信号分别为:
$\;\left. \begin{gathered} {g_1} = [1{\mkern 1mu} {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot \frac{{2v}}{{\lambda \cdot {f_{\text{p}}}}} + {\text{j}}\pi \cdot 1 \cdot \frac{{2a}}{{\lambda \cdot f_{\text{p}}^{\;2}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot \frac{{2v}}{{\lambda \cdot {f_{\text{p}}}}} + {\text{j}}\pi \cdot {{(K - 1)}^2} \cdot \frac{{2a}}{{\lambda \cdot f_{\text{p}}^{\;2}}}}}]_{(K \times 1)}^{\text{T}} \hfill \\ = [1{\mkern 1mu} {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot 1 \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(K - 1)}^2} \cdot {a_{\text{d}}}}}]_{(K \times 1)}^{\text{T}}n \hfill \\ {g_2} = [{{\text{e}}^{{\text{j}}2\pi \frac{{d\cos \theta }}{\lambda }}}\quad {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot \frac{{2v}}{{\lambda \cdot {f_{\text{r}}}}} + {\text{j}}\pi \cdot 1 \cdot \frac{{2a}}{{\lambda \cdot f_{\text{r}}^2}} + {\text{j}}2\pi \frac{{d\cos \theta }}{\lambda }}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot \frac{{2v}}{{\lambda \cdot {f_{\text{r}}}}} + {\text{j}}\pi {{(K - 1)}^2} \cdot \frac{{2a}}{{\lambda \cdot f_{\text{r}}^2}} + {\text{j}}2\pi \frac{{d\cos \theta }}{\lambda }}}]_{(K \times 1)}^{\text{T}} \hfill \\ \quad = [{{\text{e}}^{{\text{j}}2\pi (d\cos \theta )/\lambda }}\quad {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot 1 \cdot {a_{\text{d}}} + {\text{j}}2\pi (d\cos \theta )/\lambda }}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(K - 1)}^2} \cdot {a_{\text{d}}} + {\text{j}}2\pi (d\cos \theta )/\lambda }}]_{(K \times 1)}^{\text{T}} \hfill \\ \end{gathered} \right\}$ | (7) |
其中,K为脉冲点数,fp为脉冲重复频率,${f_{\rm{d}}} = \frac{{2v}}{{\lambda \cdot {f_{\rm{p}}}}}$为初始频率,${a_{\rm{d}}} = \frac{{2a}}{{\lambda \cdot \mathop f\nolimits_{\rm{p}}^2 }}$为调频率,v为目标的初始速度,a为目标的加速度,λ为雷达工作波长,$\theta $是目标方位。
当第1个阵元的脉冲数由K增加到2K时,第1个阵元接收到的回波数据可表示为:
$\begin{gathered} {{g'}_1} = [\underbrace {1\quad {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot {f_{\text{d}}} + \pi \cdot 1 \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(K - 1)}^2} \cdot {a_{\text{d}}}}}}_k \hfill \\ \underbrace {{{\text{e}}^{{\text{j}}2\pi \cdot K \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {K^2} \cdot {a_{\text{d}}}}}\quad {{\text{e}}^{{\text{j}}2\pi \cdot (K + 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(K + 1)}^2} \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (2K - 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(2K - 1)}^2} \cdot {a_{\text{d}}}}}}_k]_{(2K \times 1)}^{\text{T}} \hfill \\ = \ [\underbrace {1\quad {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot 1 \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot {f_{\text{d}}} + {\text{j}}\backslash \pi \cdot {{(K - 1)}^2} \cdot {a_{\text{d}}}}}}_K\;\underbrace {1\quad {{\text{e}}^{{\text{j}}2\pi \cdot 1 \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot 1 \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi \cdot (K - 1) \cdot {f_{\text{d}}} + {\text{j}}\pi \cdot {{(K - 1)}^2} \cdot {a_{\text{d}}}}}}_K]_{(2K \times 1)}^{\text{T}} \hfill \\ \odot [\underbrace {1\;1\quad \cdots \quad 1}_K\;\underbrace {{{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}}}}\quad {{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}} + {\text{j}}2\pi K \cdot 1 \cdot {a_{\text{d}}}}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}} + {\text{j}}2\pi K \cdot (K - 1) \cdot {a_{\text{d}}}}}}_K]_{(2K \times 1)}^{\text{T}} \hfill \\ = [{g_1}\;{g_2}]_{(2K \times 1)}^{} \odot [\underbrace {1\;1\; \cdots \;1}_K \hfill \\ \underbrace {{{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}} - \frac{{{\text{j}}2\pi d\cos {\theta _0}}}{\lambda }}}\quad {{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}} + {\text{j}}2\pi K \cdot 1 \cdot {a_{\text{d}}} - \frac{{{\text{j}}2\pi d\cos {\theta _0}}}{\lambda }}}\quad \cdots \quad {{\text{e}}^{{\text{j}}2\pi K{f_{\text{d}}} + {\text{j}}\pi {K^2}{a_{\text{d}}} + {\text{j}}2\pi K \cdot (K - 1) \cdot {a_{\text{d}}} - \frac{{{\text{j}}2\pi d\cos {\theta _0}}}{\lambda }}}}_K]_{(2K \times 1)}^{\text{T}} \hfill \\ = {[{g_1}\;{g_2} \odot {{\text{e}}^{\Delta {\varphi _2}}}]_{(2K \times 1)}} \hfill \\ \end{gathered} $ | (8) |
其中$ \odot $表示Hadamard积,由式(8)可见,只要补偿第2个阵元的相位后再进行拼接就相当于将两个阵元的脉冲点拼接到第1个阵元上,由式(8)得到阵元2相应的相位补偿量应为:
$\Delta {\varphi _2} = {\rm{ }}[{\rm{j}}2\pi K{f_{\rm{d}}} + {\rm{j}}\pi {K^2}{a_{\rm{d}}} - {\rm{j}}2\pi (d\cos \theta ){\rm{ /}}\lambda \quad {\rm{j}}2\pi K{f_{\rm{d}}}$ $ + {\rm{j}}\pi {K^2}{a_{\rm{d}}} + {\rm{j}}2\pi K \cdot 1 \cdot {a_{\rm{d}}} - {\rm{j}}2\pi (d\cos \theta ){\rm{ /}}\lambda \; \cdots {\rm{ }}$ ${\rm{j}}2\pi K{f_{\rm{d}}} + {\rm{j}}\pi {K^2}{a_{\rm{d}}} + {\rm{j}}2\pi K \cdot (K - 1) \cdot {a_{\rm{d}}} - {\rm{j}}2\pi (d\cos \theta ){\rm{ /}}\lambda]_{(K \times 1)}^{\rm{T}}$ | (9) |
以此类推,第n个阵元对应的补偿相位为:
$\Delta {\varphi _n} = {\rm{ }}[{\rm{j}}2\pi (n - 1)K{f_{\rm{d}}} + {\rm{j}}\pi {[(n - 1)K]^2}{a_{\rm{d}}} - {\rm{j}}2\pi [d(n - 1)\cos \theta]{\rm{ /}}\lambda $ ${\rm{j}}2\pi (n - 1)K{f_{\rm{d}}} + {\rm{j}}\pi {[(n - 1)K]^2}{a_{\rm{d}}} + {\rm{j}}2\pi (n - 1) \cdot K \cdot 1 \cdot {a_{\rm{d}}} - {\rm{j}}2\pi [d(n - 1)\cos \theta]{\rm{ /}}\lambda \; \cdots $ ${\rm{j}}2\pi K{f_{\rm{d}}} + {\rm{j}}\pi {[(n - 1)K]^2}{a_{\rm{d}}} + {\rm{j}}2\pi (n - 1)K \cdot (K - 1) \cdot {a_{\rm{d}}} - {\rm{j}}2\pi [d(n - 1)\cos \theta]{\rm{ /}}\lambda]_{(K \times 1)}^{\rm{T}}$ | (10) |
其中,n=2,3,$ \cdots $,N。根据式(10)对N个阵元进行相位补偿后,可以将N个阵元的脉冲数据拼接到1个阵元上,起到了利用空域信息来重构时域信息的效果,进行WVD变换后可使目标能量得到有效积累,用于进行目标参数估计。
3.2 构造代价函数由式(10)可知,在重构时间采样时,所需的补偿相位中包含未知的方位角$\theta $和初始频率项fd以及调频率项ad。${f_{\rm{d}}} = \frac{{2v}}{{\lambda \cdot {f_{\rm{p}}}}}$中包含未知参数v,${a_{\rm{d}}} = \frac{{2a}}{{\lambda \cdot \mathop f\nolimits_{\rm{p}}^2 }}$中包含未知参数项a,所以无法直接将多个阵元的回波数据(杂波抑制后)进行拼接。在构造代价函数时需进行速度、加速度、方位角3维搜索。代价函数可表示为:
$(\hat v,\hat a) = \arg \mathop {\max }\limits_{(v,a,\theta )} \left(\frac{1}{{{\large \sum} {\left\{ {{{\left| {{W_{\rm{s}}}\Big({{\Big[{{[{{x}_{{\rm{s,1}}}} \odot {{\rm{e}}^{\Delta {\varphi _1}}}]}^{\mathop{\rm T}\nolimits} }\; {{[{{x}_{{\rm{s,2}}}} \odot {{\rm{e}}^{\Delta {\varphi _2}}}]}^{\mathop{\rm T}\nolimits} }\; \cdots \; {{[{{x}_{{\rm{s,}}n}} \odot {{\rm{e}}^{\Delta {\varphi _n}}}]}^{\mathop{\rm T}\nolimits} }\Big]}^{\mathop{\rm T}\nolimits} }\Big)} \right|}_{NK \times NK}}} \right\}} }}\right)$ | (11) |
其中xs,n为第n个阵元杂波抑制后的回波数据,${\varphi _n}$为第n个阵元所需的补偿相位。Ws($ \cdot $)为WVD变换算子符号,$ {\small \sum} {\left\{ \cdot \right\}} $是对矩阵内所有元素相加。当拼接效果最好时,交叉项成分最少,由于WVD的交叉项成分具有震荡特性,且交叉项的平均值为零[12],故此时WVD变换后能量密度的绝对值之和最小。将WVD变换后能量密度的绝对值之和的倒数的最大值作为代价函数,参数估计的结果可以通过搜索代价函数的最大值得到。
3.3 方法流程本文所提方法利用了重构时间采样(其作用等效于增加单个阵元脉冲采样数)和WVD对LFM信号的频率聚集性,二者结合来提高参数估计精度,通过构造代价函数进行3维搜索实现目标参数估计。图 2给出了该方法的实现流程,随后给出了具体操作步骤。
步骤1 杂波抑制。本文采用文献[16, 17]中的方法对地杂波进行抑制,杂波抑制后的数据可表示为:
$Y = {\widehat {{\rm{ }}R{\rm{t}}}^{ - 1}}X$ | (12) |
其中$\widehat {{\rm{ }}R}$为估计杂波协方差矩阵,利用待检测单元邻近的距离单元训练样本进行估计得到的,即:
$\widehat {{\rm{ }}R} = \frac{1}{L}\mathop \Sigma \limits_{i = 1}^L {\rm{ }}x(i){x^{\rm{H}}}(i)$ | (13) |
其中,L为训练样本的个数,x(i)为第i(i=1,2,$ \cdots $,L)个训练样本的空时快拍数据。
步骤2 确定初值进行搜索。为降低3维搜索的运算量,先确定参数对中$a,v,\,\theta $初值,首先利用延迟自相关解线性调频算法估计加速度初值,设数据xs,n延迟$\tau $后的数据矢量为xds,n,则瞬时自相关函数为:
${r_n} = {x_{{\rm{s}},n}} \odot {\rm{ }}x_{{\rm{ds}},n}^*$ | (14) |
对瞬时自相关函数rn做FFT估计其频率估值fn,对应的目标加速度的估值为an=fn/2$\tau $。对每一个阵元的接收数据估计加速度并求平均,可得粗略估值作为加速度初值表示为:
${a_c} = \frac{1}{N}\mathop {\LARGE \Sigma} \limits_{n = 1}^N {a_n}$ | (15) |
另外,本文分别使用了单脉冲测角法和FFT测频法来估计方位角$\theta $和速度v的初值,估计初速度时,先用每个阵元数据估计目标速度初值,最终初值为各阵元数据估计结果的平均值。以这些初值构造搜索范围进行3维搜索。
步骤3 3维搜索中,先重构时间采样后WVD变换构造代价函数。本文方法在3维搜索中,先对个阵元的数据进行相位拼接完成重构时间采样,WVD变换后得到重构数据的WVD能量谱。根据式(11)构造代价函数,在最优匹配时WVD变换后能量密度值的绝对值之和最小,即在变步长搜索中代价函数取得最大值,此时所对应的参数v0,a0为该目标参数的最优解。
在本文方法中,为降低算法运算量,3维搜索采用了变步长搜索,先使用较大搜索步长进行粗搜索,来确定后续精搜索的初值,取该初值的左右各两个粗搜步长构造精搜索范围,精搜索使用较小的搜索步长来确保搜索的精确度。
4 仿真分析本节将通过仿真实验来验证该方法的有效性。仿真参数设置:天线为阵元间距是d=0.5λ,阵元数为N=16的正侧视理想均匀线阵。载机速度为120 m/s,载机平台高度为10 km,雷达的距离分辨率为20 m,雷达工作波长为λ=0.32 m,相干处理脉冲数为K=64,脉冲重复频率为fp=1500 Hz,输入信噪比SNR=0 dB,杂噪比CNR=50 dB。设机动目标处于待测单元内,处于方位角90°处,初速度为24.01 m/s,加速度为99.9 m/s2。实验估计参数均方根误差进行了200次蒙特卡罗实验。
图 3(a)为总回波的功率谱,因为信杂比较低,回波信号中杂波信号完全覆盖了目标信号。图 3(b)是用杂波协方差矩阵求逆法来抑制杂波后的功率谱,由于未补偿目标加速度,目标在杂波背景下仍然难以检测。图 3(c)为本文方法得出的功率谱,加速度补偿后,目标能量在多普勒域中得以聚集。图 4中蓝线和红线分别为图 3(b)和图 3(c)中$\cos \psi = 0$(其中$\cos \psi = \cos \theta \cos \gamma $)时的功率谱。由图 4可见,未补偿加速度时,由于目标加速度导致的多普勒走动,使得其回波信号功率谱在多普勒域上有着一定的展宽,不利于对机动目标进行检测。而本文方法利用步骤2得出的加速度初值对目标加速度进行补偿后目标能量得到了更为有效的累积,有利于对机动目标的检测。
图 5表示所估参数的均方根误差随信噪比的变化。本文所使用的3维搜索为变步长搜索,具体搜索参数设置为:加速度的搜索范围为初值的±100 m/s2,粗搜索步长为4.0 m/s2,细搜索步长为0.5 m/s2,速度的搜索范围为初值的±60 m/s,粗搜索步长为2.0 m/s,细搜索步长为0.2 m/s,方位角的搜索范围为初值的±10°,粗搜索步长为1°,细搜索步长为0.2°。图中单阵元法为仅利用1个阵元数据做分数阶傅里叶变换后的参数估计结果,非相干积累法是指对所有阵元做分数阶傅里叶变换后的参数估计结果进行统计平均后得到的均值。图 5(a)为初速度均方根误差随信噪比的变化,图 5(b)为加速度均方根误差随信噪比的变化。结果表明,在补偿加速度并进行重构时间采样后,本文方法比单阵元法和非相参积累法参数估计精度有了显著的提高,且当信噪比较低时本文方法估计效果更好。图 6给出了不同搜索间隔条件下的参数估计性能曲线,由此可见,使用较小的步长可以提高参数估计精度。
表 1给出了本文方法中各个步骤所需运算量,其中L为估计协方差矩阵所需样本数,Vcz1,acz1,$\theta $cz1和$\Delta {v_1}$,$\Delta {a_1}$,$\Delta {\theta _1}$分别表示速度、加速度和方位的粗搜范围和步长,Vcz2,acz2,$\theta $cz2和$\Delta {v_2}$,$\Delta {a_2}$,$\Delta {\theta _2}$分别表示速度、加速度和方位的细搜范围和步长。
表 2给出了本文中3种方法进行1次蒙特卡洛实验的运算时间对比,仿真实验采用Matlab R2014a软件进行仿真,计算机参数为Inter(R) Core(TM) i3-3220处理器;主频:3.30 GHz;内存:8.00 GB。本文方法虽然较之单阵元法和非相干积累法运算量有着增大,但其估计精度有着显著的提高。
本文所提出的基于WVD的机动目标检测与参数估计方法能够在目标方位未知,脉冲数目较少的情况下,得到较高精度的参数估计。通过仿真实验表明,该方法利用WVD和重构时间采样对回波信号进行相干积累,进行变步长3维搜索后得到的参数估计精度较传统的单阵元法和非相参积累法更高,仿真实验证明了该方法的有效性。另外,在本文中虽然通过变步长搜索来实现降低算法运算量,但其运算量仍然较大,所以在保证参数估计精度的同时降低运算量是下一步研究和改进的重点。
[1] | Klemm R. Principle of Space-time Adaptive Processing[M]. 3rd Edition, UK, IET Publishers, 2006.(1) |
[2] | 丁鹭飞, 耿富录. 雷达原理[M]. 西安: 电子科技大学出版社, 2002: 1-24.Ding Lu-fei and Geng Fu-lu. The Principle of Radar[M]. Xi'an: Publishing House of University of Electronics Science and Technology, 2002: 1-24.(2) |
[3] | 罗守贵, 金林. 机载预警雷达发展趋势分析[M]. 南京: 南京电子技术研究所, 2008.Luo Shou-gui and Jin Lin. The Analysis of Airborne Early Warning Radar Development Trend[M]. Nanjing: The Fourteenth Institute, CETC, 2008.(1) |
[4] | 王永良, 彭应宁. 空时自适应信号处理[M]. 北京: 清华大学出版社, 2000: 1-57.Wang Yong-liang and Peng Ying-ning. Space-Time Adaptive Processing[M]. Beijing: Publishing House of Tsinghua University, 2000: 1-57.(1) |
[5] | 刘建成. 加速运动目标检测及跟踪技术研究[D]. [博士论文], 国防科技大学, 2007.Liu Jian-cheng. Study on accelerating target detection and tracking[D]. [Ph.D. dissertation], National University of Defense Technology, 2007.(1) |
[6] | James Ward. Maximum likelihood angle and velocity estimation with space-time adaptive processing radar[C]. Conference Record of the Asilomar Conference on Signals, Systems and Computers, USA, 1997: 1265-1267.(1) |
[7] | Rao P and Taylor F J. Estimation of instantaneous frequency using the discrete Wigner distribution[J]. Electronics Letters, 1990, 26(4): 246-248.(1) |
[8] | Boashash B. Estimating and interpreting the instantaneous frequency of a signal[J]. Proceedings of the IEEE, 1992, 80(4): 519-569.(1) |
[9] | Xu Jia, Xia Xiang-gen, Peng Shi-bao, et al.. Radar maneuvering target motion estimation based on generalized Radon-Fourier transform[J]. IEEE Transactions on Signal Processing, 2012, 60(12): 6190-6201.(1) |
[10] | 陈小龙, 刘宁波, 王国庆, 等. 基于Radon-分数阶傅里叶变换的雷达动目标检测方法[J]. 电子学报, 2014, 42(6): 1074-1080.Chen Xiao-long, Liu Ning-bo, Wang Guo-qing, et al.. Radar detection method for moving target based on Radon-Fractional fourier transform[J]. Acta Electronica Sinica, 2014, 42(6): 1074-1080.(3) |
[11] | Chen Xiao-long, Huang Yong, Liu Ning-bao, et al.. Radon-Fractional ambiguity function-based detection method of low-observable maneuvering target[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(2): 815-832.(2) |
[12] | 王小寒. 基于重构时间采样的空中机动目标检测研究[D]. [硕士论文], 中国民航大学, 2012.Wang Xiao-han. Detection of maneuvering targets via virtually enlarging time samples[D]. [Master dissertation], Civil Aviation University of China, 2012.(6) |
[13] | 谢鼎. Wigner_Ville分布在调频信号处理中应用[D]. [硕士论文], 西安电子科技大学, 2010.Xie Ding. The application of Wigner-Ville Distribution in frequency-modulation signal processing[D]. [Master dissertation], Electronics Science and Technology University of Xi'an, 2010.(1) |
[14] | 张华. 低信噪比下线性调频信号的检测与参量估计研究[D]. [硕士论文], 电子科技大学, 2004.Zhang Hua. Study on LFM signal detection and parameters estimation in low SNR[D]. [Master dissertation], University of Electronic Science and Technology, 2004.(1) |
[15] | Ward J. Space-time Adaptive Processing for Airborne Radar[M]. Massachusetts: Massachusetts Institute of Technology Lincoln Laboratory, 1994: 461-464.(1) |
[16] | 吴仁彪, 贾琼琼, 李海. 机载雷达高速空中微弱动目标检测新方法[J]. 电子与信息学报, 2011, 33(6): 1459-1464.Wu Ren-biao, Jia Qiong-qiong, and Li Hai. Detection of fast moving dim targets on airborne radar via STAP[J]. Journal of Electronics & Information Technology, 2011, 33(6): 1459-1464.(1) |
[17] | 郑景忠. 短相干处理时间情况下的机动目标检测与参数估计[D]. [硕士论文], 中国民航大学, 2015.Zheng Jing-zhong. Detection and parameter estimation of maneuvering target under the circumstance of short coherent processing Interval[D]. [Master dissertation], Civil Aviation University of China, 2015."(1) |