基于机动目标距离-瞬时调频的ISAR成像方法

摘要

本发明公开了一种ISAR瞬时成像的方法。该方法是：首先，对距离压缩后的回波信号进行包络对齐和距离走动校正；然后，利用发明的时间-调频率分布TCD结合Clean算法对各距离单元信号进行处理，得到指定时刻瞬时多普勒谱和瞬时调频率谱；最后，根据得到的瞬时多普勒谱重构目标的距离瞬时调频率图像，根据得到的瞬时调频率谱重构目标的图像，实现对机动目标的动态ISAR成像。本发明不同于以往的方法，是一种新型的瞬时成像方法，可以获得更多的目标信息并且大大改善对机动目标的识别能力，成像效果好。

说明书

技术领域

本发明属于雷达技术领域，涉及成像，具体说的是一种逆合成孔径雷达ISAR成像方法，用于目标检测与识别。 

背景技术

ISAR是一种高分辨小旋转角成像技术，在距离方向，它通过发射大带宽信号获得高分辨；在方位方向，高分辨通过目标相对转动时回波间的相干积累取得。近二十年来，ISAR成像理论和技术有了很大发展，并逐步应用于实际系统，它能对目标进行全天候、全天时和远距离观察，获得目标的高分辨率图像，以对目标进行分类和识别。 

ISAR成像技术已经相对比较成熟。对平稳飞行目标，在成像的各个环节，包络对齐，自聚焦和横向成像，都有一系列卓有成效的快速算法。例如，在包络对齐方面，有基于频域快速相关法，散射重心加近距离单元滑动相关法，基于范数1距离最近的滑动对齐法，复包络相关法，最小熵方法，全局相关，和超分辨方法。在自聚焦方面，有单特显点法，多特显点综合法，散射重心方法，多普勒中心跟踪法。在成像方面，主要用横向傅氏变换，对转角较小不满足横向分辨要求时，用Burg外推，RELAX等超分辨方法。对于目标机动飞行，当目标不很大，散射点不发生大的距离走动时，平稳飞行的包络对齐方法仍然适合于机动飞行。一些自聚焦方法也适用于机动飞行目标，如散射中心法。 

然而，由于目标运动的非合作性，若目标在观测期间作机动运动，即目标相对于雷达视线的姿态成偏航、侧摆和俯仰三维运动。经过平动补偿后，运动目标等效为非均匀转动的转台目标，此时单个散射点方位回波信号的变化规律一般比较复杂，多普勒是时变的，为调幅-调频信号。因而，机动目标的ISAR成像的方位向处理也就可以归结为多分量调幅-调频信号的参数估计和提取的问题。已有的对平动补偿后的机动目标的动态成像的方法，如Radon-Wigner瞬时成像，联合时频分布的方法，基于自适应Chirplet分解的ISAR成像算法等，这些方法通常都假设单个散射点的方位信号为等幅线性调频信号LFM。实际中回波信号的高次相位项会使信号在时频 面上的分布成为曲线并产生严重的交叉项，造成成像质量下降。 

发明内容

本发明的目的在于克服上述已有技术的不足，提供一种运基于机动目标距离-瞬时 

调频的ISAR成像方法，以抑制回波信号的高次相位项产生的交叉项，有效的提高成像质量。 

为实现上述目的，本发明包括如下步骤： 

(1)对原始数据信号依次进行距离压缩、包络对齐和距离走动校正；并将距离走动校正后的数据划分为若干个距离单元，定义每个距离单元的回波序列总信号为： 

$s\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}{s}_k\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}{\sigma }_k\left(t\right)·\exp \left[j({\Phi }_{k,0}+f_k·t+\frac{1}{2}{\gamma }_k·t^2+\frac{1}{6}{\hat{\gamma }}_k·t^3)\right]$

其中，k表示信号分量的序号，σ_k(t)，Φ_k，O，f_k，γ_k， j，t分别表示成像时刻第k个信号分量的幅度，初相，起始频率，起始调频率，二次调频率，虚数符号和时间变量； 

(2)取其中一个距离单元的数据，计算出该距离单元数据信号的能量；取k＝1； 

(3)估计第一个信号分量s₁(t)的起始调频率γ₁和二次调频率 并由此得到瞬时调频率为 

(4)根据起始调频率γ₁和二次调频率 构造参考信号： 

${\mathrm{sref}}_1\left(t\right)=\exp \left[j(\frac{1}{2}{\gamma }_1{t}^2+\frac{1}{6}{\hat{\gamma }}_1{t}^3)\right];$

(5)利用构造的参考信号对总信号进行解调频处理，并对解调频处理后的信号进行傅里叶变换，得到解调频处理后的信号谱S_d(f)； 

(6)从解调频处理后的信号谱中滤出第一信号分量的信号谱：S₁(f)＝WIN₁(f)·S_d(f)， 

其中 f为频率变量，f₁为初始频率，Δf₁为第一信号分量的带宽； 

(7)从解调频处理后的信号谱S_d(f)中除去第一信号分量的信号谱S₁(f)，得到剩余信号的信号谱S_r(f)； 

(8)对剩余信号的信号谱进行逆傅里叶变换，并与构造的参考信号相乘，得到剩余信号S_r(t)，并计算该剩余信号剩余能量：E＝∑|s_r(t)|²，定义剩余能量比率为ξ＝E/E₀，其中E₀＝∑|s(t)|²，为总信号的能量； 

(9)对剩余信号的信号谱S_r(f)做逆傅里叶变换，得到散射点的幅度σ_k(t)；将总信号更新为此时的剩余信号s(t)＝s_r(t)，并将信号分量的序号k加1，进入下一个信号分量的计算； 

(10)重复执行步骤(3)-(9)，直到剩余能量比率ξ小于预定的门限值0.1； 

(11)重复执行步骤(2)-(10)，直到计算完所有的距离单元； 

(12)将各个距离单元中散射点的幅度σ_k(t)按照其瞬时调频率γ_k(t)的值填入距离-瞬时调频二维平面上的相应位置，得到目标t时刻的距离瞬时调频ISAR图像。 

本发明的优势在于，提出了一种对机动目标的距离-瞬时调频率的ISAR成像方法，通过对时间-调频率分布平面上信号的分析和积累，能够抑制回波信号的高次相位产生的交叉项的影响，有效提取由于目标高机动产生的具有复杂回波相位特性的ISAR信号，并重构信号实现对目标的高分辨瞬时RID和RIC成像，提高对目标的ISAR成像概率。具体表现在：通过上面所述的步骤(3)和(4)有效估计复杂运动目标产生的高次相位项，并提取调频率及其变化率特性；之后，通过步骤(5)和(6)对回波信号的中心频率和幅度能力信息进行有效估计；接着，利用已有的信息对目标回波信号进行重构；最后，通过步骤(12)完成对目标的瞬时调频率成像。 

附图说明

图1是本发明基于机动目标的距离-瞬时调频的ISAR成像流程图； 

图2是本发明在估计起始调频率和二次调频率时所使用的交叉项示意图； 

图3是本发明估计起始调频率和二次调频率时对第一个信号分量进行变换的结果图； 

图4是本发明采样TCD和CLEAN方法估计起始调频率和二次调频率的结果图； 

图5是对舰船目标进行仿真的结果图； 

图6是对舰船目标的仿真数据进行成像的结果图； 

图7是对舰船目标的实测数据进行成像的结果图。 

具体实施方式

参照图1，本发明的ISAR成像过程包括： 

步骤1，原始数据为若干个距离单元。 

(1a)对舰船目标原始数据进行距离压缩； 

(1b)对距离压缩后的数据进行包络对齐； 

(1c)对包络对齐后的数据进行距离走动校正； 

(1d)将距离走动校正后的数据划分为若干个距离单元，并定义每个距离单元的回波序列总信号如下： 

$s\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}{s}_k\left(t\right)=\underset{k=1}{\overset{K}{\Sigma }}{\sigma }_k\left(t\right)·\exp \left[j({\Phi }_{k,0}+f_k·t+\frac{1}{2}{\gamma }_k·t^2+\frac{1}{6}{\hat{\gamma }}_k·t^3)\right]---1)$

其中σ_k(t)，Φ_k，0，f_k，γ_k， j和t分别表示成像时刻第k个信号分量的幅度，初相，起始频率，起始调频率，二次调频率，虚数符号和时间变量。 

步骤2，取距离走动校正后的第一个距离单元的数据，计算出该距离单元数据信号的能量E₀＝∑|s(t)|²。 

步骤3，估计第一个信号分量s₁(t)的起始调频率γ₁和二次调频率 得出瞬时调频率为γ_k(t)。 

(3a)将第一个信号分量，即二次调频率QFM信号S₁(t)表示为如下形式 

$s_1\left(t\right)=\exp \left[\mathrm{j\Phi }\left(t\right)\right]=\exp \left[j({\Phi }_0+f_1·t+\frac{1}{2}{\gamma }_1·t^2+\frac{1}{6}{\hat{\gamma }}_1·t^3)\right]$

2)其中，Φ₀表示初相，f₁为初始频率，γ₁代表调频率， 是频率的二次项。 

(3b)根据时间调频率分布非常适合用来估计瞬时调频率的特性，将二次调频率QFM信号S₁(t)转化到时间-调频平面，其分布表示为： 

TCD_s(t，c)＝∫Z_s(t，τ)exp(-jcτ)dτ 

          ＝∫exp{-jτ[c-γ(t)]}dτ            3) 

          ＝2πδ[c-γ(t)] 

其中$Z_s(t,\tau )=\left(\begin{array}{cc}{s}_1(t+\sqrt{\tau }){s_1}^{*}\left(t\right){s_1}^{*}\left(t\right)s_1(t-\sqrt{\tau }),& \tau \ge 0\\ Z_s^{*}(t,-\tau ),& \tau <0\end{array}\right)$

4) 

Z_s(t，τ)为信号的乘积， 为Z_s(t，τ)的复共轭转置信号，τ为时间延时变量，c为时间-调频平面上的直线， 为瞬时调频率，δ[·]为冲击函数；显然，S₁(t)信号的时间-调频率分布为沿着瞬时调频率γ(t)集聚的一条直线； 

(3d)根据时间-调频率分布TCD有时移不变、频移不变和调频率不变性的性质，对式1)中的二次调频率QFM信号S₁(t)进行时移、频移和调频率偏移变换，得到变换后的二次调频率QFM信号： 

${s_1}^{\prime }\left(t\right)=\exp \left\{j\right[{\Phi }_0+(f+f_0)·(t-t_0)+\frac{1}{2}(\gamma +c_0)·{(t-t_0)}^2+\frac{1}{6}\hat{\gamma }·{(t-t_0)}^3\left]\right\}---5)$

其中，t₀为时间偏移量，f₀为频率偏移量，c₀为调频率偏移量，γ为调频率， 是频率的二次项； 

(3e)将式5)代入式3)，得到信号s₁′(t)时间-调频分布TCD： 

${\mathrm{TCD}}_{s^{\prime }}(t,c)\left(\begin{array}{c}=\int \exp \{-j·\tau [c-c_0-\gamma -\hat{\gamma }·(t-t_0)\left]\right\}\mathrm{d\tau }\\ =2\mathrm{\pi \delta }\{c-c_0-[\gamma +\hat{\gamma }·(t-t_0)\left]\right\}={\mathrm{TCD}}_s(t-t_0,c-c_0)\end{array},-,-,-,6,)\right)$

上式的结果表明，时间-调频分布TCD不会受到二次及以下相位项的影响；因此，利用时间-调频分布TCD对一个距离单元内散射点的回波序列信号进行处理，可以得到散射点的瞬时调频率；信号中的线性及二次相位误差不会影响距离-瞬时调频率RIC图像的质量，最多会引起RIC图像的方位向偏移； 

(3f)以两个散射点P₁和P₂的回波信号为例，取[T_C1，T_C2]时间段内两散射点的回波数据进行时间-调频分布TCD变换，得到的两个散射点回波信号的理想瞬时调频率谱如图2所示；由图2可见，时间-调频分布TCD的自身项受到了交叉项干扰，其中，[T_C1，T_C2]为适合做距离-瞬时调频率RIC成像的时间段，该交叉项是由于一个距离单元内多个散射点回波的多分量的二次调频率QFM信号的时间-调频分布TCD存在而产生的。 

抑制交叉项干扰的方法就是对时间-调频平面进行Radon变换，也就是在时间-调频率平面沿着任意直线对TCD_s(t，c)求线积分得Radon变换后的曲线分布： 

D_s(c₀，c₁)＝∫_LTCD_s(t，c)ds         7) 

其中，L为用于线积分的任意直线，c₀和c₁分别为此曲线分布的横、纵坐标变量，c₀，c₁分别表示第一个信号分量s₁(t)的起始调频率γ₁和二次调频率 的值在曲线分布D_s中的尖峰位置的坐标值，并由此可得瞬时调频率γ_k(t)， 其中γ₁＝c₀， 

图3(b)为对图3(a)中的时间-调频分布TCD进行Radon变化得到的曲线分布图，其中两个尖峰代表了两个散射点P₁和P₂各自的回波信号分量，由图3(b)中尖峰的位置可估计出信号分量的瞬时调频率； 

步骤4，根据初始调频率γ₁和二次调频率 构造参考信号。 

在有多个信号分量的情况下，信号分量中存在弱信号分量，弱信号分量有可能会被强信号分量的交叉项所掩盖，从而导致无法准确的估计出信号分量参数，为解决此问题，需要结合′CLEAN′技术来估计多分量二次调频率信号QFM的瞬时多普勒频率和调频率。 

通过估计(3f)中尖峰位置坐标得到的γ₁和 构造以下参考信号： 

${\mathrm{sref}}_1\left(t\right)=\exp \left[j(\frac{1}{2}{\gamma }_1{t}^2+\frac{1}{6}{\hat{\gamma }}_1{t}^3)\right].---8)$

步骤5，利用参考信号sref₁(t)对式1)所示的总信号s(t)进行解调频处理，可得解调频后的信号： 

$s_d\left(t\right)=s\left(t\right)·{\mathrm{sref}}_1^{*}\left(t\right)=s_1\left(t\right)+s_r\left(t\right)---9)$

其中 

S₁(t)＝σ₁·exp[j(Φ_1，0+f₁·t)] 

10) 

$s_r\left(t\right)=\underset{k=2}{\overset{K}{\Sigma }}{\sigma }_k·\exp \left\{j\right[{\Phi }_{k,0}+f_i·t+\frac{1}{2}({\gamma }_k-{\gamma }_1)·t^2+\frac{1}{6}({\hat{\gamma }}_k-{\hat{\gamma }}_1)·t^3\left]\right\}---11)$

上式中s_r(t)则表剩余信号，Φ_1，0表示第一个信号分量的初相， 表示参考信号sref₁(t)的共轭，可以发现S_d(t)中的所有信号分量的调频率和二次调频都分别减去了γ₁和 使第一个信号分量变为一个单频信号。 

步骤6，对S_d(t)进行傅立叶变换，可得解调频处理后的信号谱为 

S_d(f)＝FFT[S_d(t)]＝S₁(f)+S_r(f)              12) 

其中FFT[·]表示傅立叶变换，S₁(f)为第一个信号分量的频谱，S_r(f)为剩余信 号分量的频谱。 

从式12)可见，只有第一个信号分量频谱S₁(f)＝FFT[s₁(t)]为窄谱，剩余信号分量的频谱S_r(f)＝FFT[s_r(t)]都是展宽的，因此很容易从S_d(f)的表达式所示的频谱中估计出第一个信号分量的多普勒频率f₁以备后续估计信号的最强幅度。 

步骤7，用窄窗在频域中从解调频后信号频谱S_d(f)中滤出将第一信号分量的窄谱： 

S₁(f)＝WIN₁(f)·S_d(f)                13) 

其中，${\mathrm{WIN}}_1\left(f\right)=\left(\begin{array}{cc}1& f_1-\Delta f_1/2<f<f_1+\Delta f_1/2\\ 0& \mathrm{others}\end{array}\right)---14)$

其中，WIN₁(f)为滤波器函数，Δf₁表示第一信号分量s₁(t)的带宽。利用式13)可从解调频后的信号频谱S_d(f)中滤除第一信号分量的主瓣，并且保证对其余信号分量的影响很小。 

步骤8，从总的信号谱中除去第一信号分量的频谱S₁(f)，得到剩余信号的频谱S_r(f)，并对其进行逆傅里叶变换，对式8)中的参考信号进行匹配，得到剩余信号： 

$s_r\left(t\right)=\mathrm{IFFT}\left[s_r\left(f\right)\right]·{\mathrm{sref}}_1\left(t\right)=\underset{k=2}{\overset{K}{\Sigma }}{\sigma }_k·\exp \left[j({\Phi }_{k,0}+f_{k}t+\frac{1}{2}{\gamma }_k{t}^2+\frac{1}{6}{\hat{\gamma }}_k{t}^3)\right]---15)$

其中，IFFT[·]代表逆傅立叶变换。 

步骤9，利用式15)中剩余信号S_r(t)计算剩余能量E，其中E＝∑|s_r(t)|²，并定义剩余能量比率为ξ＝E/E₀。 

步骤10，将式13)中滤出的第一信号分量频谱S₁(f)做逆傅里叶变换： 

s₁(t)＝IFFT[S₁(f)]·exp(-jf₁·t)＝σ₁(t)·exp(jΦ_1，0)          16) 

其中，σ₁(t)为第一信号分量s₁(t)的幅度； 

步骤11，用式16)中的第一信号分量幅度σ₁(t)更新式1)中总信号幅度σ_k(t)，即σ_k(t)＝σ₁(t)，此时k为1，这是由于式1)中的信号分量序号k是一个随步骤循环执行而递增的值，所以当重复执行以上步骤时能够得出第k个信号分量的幅度，然后再把信号分量的幅度值赋值给σ_k(t)。 

步骤12，更新总信号s(t)，令s(t)＝s_r(t)，并使信号分量序号k增加1，然后重复执行步骤3-步骤11，处理剩余信号s_r(t)直到剩余能量比率ξ小于预定的门限值，例如ξ为0.1。 

步骤13，重复执行步骤2-步骤12，计算完所有的距离单元后，得到距离单元信号中各散射点子回波的幅度和瞬时调频率。 

由于前面所估计的幅度σ_k(t)和瞬时调频率γ_k(t)分别反映了t时刻目标散射点的强度和横向位置，故将各个距离单元中散射点的幅度σ_k(t)按照其瞬时调频率γ_k(t)的值填入距离-方位二维平面上的相应位置，就能够得到目标t时刻的距离-瞬时调频RIC的ISAR图像。 

本发明的效果可以通过以下仿真和实测数据进一步说明： 

仿真1： 

对两散射点P₁和P₂的回波信号变换到时间-调频平面进行仿真，首先把如图3(a)所示的两闪射点P₁和P₂总回波信号的时间-调频分布TCD根据本发明的时间-调频分布TCD结合CLEAN算法进行处理，估计出信号中两个分量的参数，有调频率，二次调频率，幅度；接着利用这些参数准确重构出两散射点P₁和P₂的单个回波信号分量，即从两散射点P₁和P₂的总回波信号中提取出了其各自的回波信号分量；接着把重构出的两个回波信号分量变换到时间-调频平面，其结果如图4所示，其中： 

图4是对两闪射点P₁和P₂总回波信号用本发明时间-调频分布TCD结合CLEAN算法提取的单个回波信号分量的时间-调频分布TCD图； 

从图4与图3(a)相比可以发现，图4中时间-调频分布TCD的交叉项得到了有效的抑制，同时对信号调频率的估计也是相当准确。 

仿真2： 

对船舰目标进行仿真，首先建立目标舰船的原始模型；接着对该模型的仿真回波信号进行距离压缩和包络对齐；接着取出处理经过距离压缩和包络对齐的回波信号中某一距离单元信号并将其变换到的时间-调频平面；接着对变换到时间-调频平面的信号进行Radon变换；接着用本发明的时间-调频分布TCD结合CLEAN算法对经过距离压缩和包络对齐的回波信号中的某一距离单元信号进行处理，估计得到单个信号分量的调频率，二次调频率，瞬时调频率和幅度；接着利用这些参数而重构出来的单个信号分量，这样重复8次就重构出了8个信号分量；最后将此8个信号分量变换到时间-调频平面，其结果如图5所示，其中： 

图5(a)是仿真目标船舰的原始模型； 

图5(b)是对目标船舰的原始模型的仿真回波信号进行距离压缩和包络对齐后，对某一距离单元信号的时间-调频分布TCD图； 

图5(c)是将图5(b)中时间-调频分布TCD进行Radon变换后的曲线分布图； 

图5(d)是对重构的8个信号分量做时间-调频率分布TCD的结果； 

从图5(c)可见，曲线分布图中的尖峰位置不能明显区分，此尖峰位置坐标分别表示各分量信号的起始调频率和二次调频率，这是由于图5(b)中时间-调频平面上的直线使得各信号分量的调频率是时变的，而且时间-调频分布TCD中的交叉项严重影响了自身项； 

从图5(d)与图5(b)中的结果对比可以发现，图5(d)中时间-调频率分布TCD的交叉项得到了明显的抑制。 

仿真3： 

对目标船舰模型的仿真回波数据，分别用已有技术距离瞬时-多普勒RID方法和本发明提出的时间-调频分布TCD结合CLEAN算法进行成像，其结果如图6所示，其中： 

图6(a)是利用传统的距离瞬时-多普勒RID成像算法对图5(a)中目标船舰模型的仿真回波数据在三个不同的成像时刻进行处理得到的三幅距离瞬时-多普勒RID图像； 

图6(b)是在与图6(a)相同的三个成像时刻，用本发明的时间-调频分布TCD结合CLEAN算法对图5(a)中目标船舰模型的仿真回波数据进行处理得到的三幅距离-瞬时调频RIC图像； 

从图6(a)可见，第二幅图像的分辨率要高于另外两个时刻，这是由于此时的多普勒谱宽要大于其它成像时刻。 

从图6(b)与图6(a)相比可见，在三个不同的成像时刻，用本发明时间-调频分布TCD结合CLEAN算法处理得到的三幅距离-瞬时调频RIC图像都具有良好质量。 

仿真4： 

对海面舰船目标的实测回波数据，分别用已有技术距离瞬时-多普勒RID方法和本发明提出的时间-调频分布TCD结合CLEAN算法进行成像，本试验中雷达工作 在X波段，带宽为240MHz，脉冲重复频率为125Hz，雷达距海平面高度约100米左右，雷达视线与舰船之间的相对斜视角约为45度，其结果如图7所示，其中： 

图7(a)是利用已有技术距离瞬时-多普勒RID成像算法对海面舰船目标的实测回波数据进行处理产生的不同时刻的三幅距离瞬时-多普勒RID图像； 

图7(b)是在与图7(a)相同的三个成像时刻，利用本发明提出的时间-调频分布TCD结合CLEAN算法对海面舰船目标的实测回波数据进行处理产生的三幅距离-瞬时调频RIC图像； 

从图7(b)中距离-瞬时调频RIC图像与图7(a)中的瞬时-多普勒RID图像相比可知，利用本发明提出的算法可以在多普勒谱很窄，已有技术距离瞬时-多普勒RID算法无法有效成像的情况下，得到有效的距离-瞬时调频RIC图像。与已有瞬时-多普勒RID成像方法相比，本发明的距离瞬时调频的成像方法显然可以得到目标的更多信息。