CN102540189B - 基于复数后向投影的自旋目标三维成像方法 - Google Patents

Abstract

本发明公开了一种基于复数后向投影的自旋目标三维成像方法，其步骤包括：(1)雷达录取ISAR回波；(2)回波相干化；(3)平动补偿；(4)算法选择；(5)二维切片成像；(6)重构回波并更新；(7)判断目标回波能量是否低于能量门限；(8)操作结束，重构目标三维图像。本发明通过回波相干化处理、平动精细补偿以及针对不同脉冲重复频率进行算法选择，克服了广义Radon变换方法无法对目标回波的平动精确补偿，受目标高速运动影响，对应图像旁瓣较高、分辨率低的不足，具有平动补偿精确，消除目标高速运动对图像的影响，提高不同脉冲重复频率条件下的成像率，对应图像最高分辨率达到波长量级，成像质量高的优点。

Description

基于复数后向投影的自旋目标三维成像方法

技术领域

本发明属于信号处理技术领域，更进一步涉及雷达成像领域中的基于复数后向投影的自旋目标三维成像方法。本发明可以有效地对空间目标进行三维成像，并对目标形状、尺寸等特征进行精确描述，为后续的目标识别提供有力保障。

背景技术

当目标高速自旋时，传统的三维成像方法已经不再适用。在观测时间内，散射点相对于成像中心的距离和多普勒为时变量，较少的有效积累数据将导致图像质量下降，甚至使距离-瞬时多普勒算法失效。因此，有必要根据自旋目标的运动特征研究新的成像方法。

Qi Wang等人在文献“High-Resolution Three-Dimensional Radar Imaging forRapidly Spinning Targets”(IEEE trans.on GRS，VOL.46，no.1，pp.22-30，2008)中提出将广义Radon变换(GRT)与CLEAN技术相结合进行目标三维成像的方法。该方法利用旋转点目标回波在距离慢时间域的周期性变化规律，通过检测正弦曲线估计散射点空间位置，并结合修正的CLEAN技术对目标进行三维成像。但该方法存在的不足是，广义Radon变换方法所得图像的主瓣较宽，分辨率低，影响对空间目标散射点分布特征的准确描述。

发明内容

本发明的目的在于克服现有技术的不足，提出一种基于复数后向投影的自旋目标三维成像方法。该方法弥补了广义Radon变换方法具有存在距离像偏移和展宽，方位像散焦，无法对目标平动进行准确补偿，成像旁瓣较高的不足，充分利用自旋散射点回波在距离-慢时间域具有正弦包络和相位这一特征，采用宽带复数后向投影算法沿不同的正弦曲线进行相干积分，并结合修正的CLEAN技术对目标进行三维成像。

实现本发明的基本思路是：首先将自旋目标的ISAR回波进行回波相干化，以补偿目标高速运动对回波产生的影响，接着以图像熵最小为准则，通过建立有效的搜索算法对目标回波的平动量进行精细补偿以得到转台模型下的目标回波然后根据不同尺寸、不同自旋频率目标对脉冲重复频率的要求，进行二维成像算法的选择，最后对目标沿转轴方向采用相应的二维成像方法对所有二维切片进行成像，最终结合修正的CLEAN技术对目标进行三维成像。

本发明的具体步骤如下：

1.基于复数后向投影的自旋目标三维成像方法，包括如下步骤：

(1)雷达录取ISAR回波；

(2)回波相干化

2a)在距离频域对回波距离向的平移和展宽分别补偿相位因子；

2b)在慢时间域对回波补偿相位因子，使得回波完全相干化；

(3)平动补偿

3a)采用傅里叶变换对回波进行距离脉压，获得距离-慢时间域回波；

3b)采用某单次距离-慢时间域回波，分别计算其与其它各次距离-慢时间域回波的滑动相关系数，具有最大滑动相关系数的回波与该回波之间的时间间隔为目标的自旋周期，对自旋周期取倒数并乘以2π，得到自旋角频率；

3c)对应最大滑动相关系数的回波与该次回波之间移动的距离数为一个周期内目标的平动量，补偿平动量实现包络粗对齐；

3d)从第二次距离-慢时间域回波开始，对于单次距离-慢时间域回波，采用广义Radon变换方法获取不同偏移量下该次回波与之前所有回波对应的二维图像并求图像熵，使图像熵最小的距离偏移量为最优距离偏移量，根据最优偏移量对回波包络进行微调，将每次回波进行包络精对齐后，最终得到平动补偿后的目标回波；

(4)算法选择

4a)采用距离-多普勒算法对平动补偿后的回波成像，得到目标的距离-多普勒图像；

4b)若距离-多普勒图像的方位维没有混叠，则选择CIRT算法，执行步骤(5)；

4c)若距离-多普勒图像出现方位维混叠，则选择RIRT算法；

(5)二维切片成像

5a)建立自旋目标的柱面坐标系，沿目标自转轴方向将自旋目标分成多个二维切片；

5b)将二维切片按照步骤(4)选择的算法进行成像，得到二维切片图像；

(6)重构回波并更新

6a)对于所有的二维切片图像，提取图像域中最大峰值对应的三维坐标，记录为目标散射点的位置坐标，并在距离-慢时间域构造散射点回波；

6b)运用最小范数准则对步骤6a)中所构造的散射点回波的幅度进行估计；

6c)将记录的散射点的位置坐标与估计的回波幅度组合构成散射点回波；

6d)从步骤3d)得到的目标回波中减去散射点回波，更新目标回波；

(7)判断目标回波能量是否低于能量门限

7a)若目标回波能量低于门限，则执行步骤(8)；

7b)若目标回波能量高于门限，则执行步骤(5)；

(8)操作结束，根据6a)中记录的目标散射点位置重构目标三维图像。

与现有技术相比，本发明具有以下优点：

第一，本发明通过采用回波相干化处理，克服了现有技术中广义Radon变换方法受目标高速运动影响，存在距离像偏移和展宽，以及方位像散焦的不足，具有消除高速运动目标对图像的影响，获得聚焦良好图像的优点。

第二，本发明通过采用以图像熵最小为准则，建立有效的搜索算法对目标回波的平动量进行精细补偿。克服了现有技术中广义Radon变换方法无法对目标平动进行准确补偿的不足，具有获得平动准确补偿后的自旋目标回波的优点。

第三，本发明通过采用针对不同脉冲重复频率条件下的算法选择，克服了现有技术中广义Radon变换方法成像旁瓣较高、分辨率低的不足，具有提高不同脉冲重复频率条件下目标的成像率，并使图像最高分辨率达到波长量级，成像质量高的优点。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图1，对本发明具体实施方式作进一步的详细描述：

步骤1，获取目标的ISAR回波，雷达以脉冲重复频率发射并接收脉冲，得到以距离为行向量以方位为列向量的ISAR回波；

步骤2，对目标的ISAR回波进行相干化处理

2a)为消除速度所引起的一维距离像平移和展宽，在距离频域对回波乘以相位补偿因子：

${\mathrm{\Φ}}_r\left(t_m\right)=\exp \left(j\frac{4\mathrm{\π}}{c}(\frac{{2R}_{\mathrm{ref}}\left(t_m\right)v}{c}+\frac{f_{c}v}{\mathrm{\γ}})f\right)\exp \left(j\frac{4\mathrm{\π}}{c}(\frac{v}{\mathrm{\γ}}-\frac{v^2}{\mathrm{c\γ}})f^2\right)$

其中，Φ_r(t_m)为距离向相位补偿因子，t_m为慢时间，j为虚数单位，c为光速，R_ref(t_m)为参考距离，v为根据不同慢时间对应的R_ref(t_m)采用多项式拟合方法估计出的目标速度，f_c为载频，γ为调频率，

为距离频率，为快时间。

2b)为使得回波完全相干化，在慢时间域对回波乘以下式的相位补偿因子：

${\mathrm{\Φ}}_a=\exp \left(j\frac{{8\mathrm{\πR}}_{\mathrm{ref}}\left(t_m\right)v}{c^2}{f}_c\right)$

其中，Φ_a为方位向相位补偿因子，j为虚数单位，R_ref(t_m)为参考距离，v为根据不同慢时间对应的R_ref(t_m)采用多项式拟合方法估计的速度，c为光速，f_c为载频。

步骤3，对目标的ISAR回波进行平动补偿

3a)采用傅里叶变换对回波进行距离脉压，即对距离频域-慢时间域回波沿距离向进行傅里叶变换，从而得到目标的距离-慢时间回波；

3b)对于某单次距离-慢时间域回波，分别计算其与其它各次距离-慢时间域回波的滑动相关系数，具有最大滑动相关系数的回波与该回波之间的时间间隔为目标的自旋周期，对自旋周期取倒数并乘以2π，得到自旋角频率；

3c)对应最大滑动相关系数的回波与该次回波之间移动的距离数为一个周期内目标的平动量，补偿平动量实现包络粗对齐；

3d)从第二次距离-慢时间域回波开始，对于单次距离-慢时间域回波，采用广义Radon变换方法获取不同偏移量下该次回波与之前所有回波对应的二维图像并求图像熵，使图像熵最小的距离偏移量为最优距离偏移量，根据最优偏移量对回波包络进行微调，将每次回波进行包络精对齐后，最终得到平动补偿后的目标回波。具体步骤如下：

第一步：设定终止循环的阈值为Δ_Threshold＝ρ_r/8，ρ_r为距离分辨率，设定初始搜索范围Δ_Scope。对于第m次回波(m≥2)，令Δ₁＝Δ_Scope/(N-1)，N为距离单元数，并令搜索阶段数k＝1；

第二步：对于第m次回波的第k个搜索阶段，设定搜索向量Δ_m，k为：

${\mathrm{\Δ}}_{m,k}=[-\frac{N-1}{2}{\mathrm{\Δ}}_k,-(\frac{N-1}{2}-1){\mathrm{\Δ}}_k,...,0,...,(\frac{N-1}{2}-1){\mathrm{\Δ}}_k,\frac{N-1}{2}{\mathrm{\Δ}}_k]+{\widehat{\mathrm{\Δ}}}_{\mathrm{ini},m}^k$

其中，Δ_k为第k个搜索阶段的搜索步长；

表示距离偏移量的初始值并满足

${\widehat{\mathrm{\Δ}}}_{\mathrm{ini},m}^k=\left\{\begin{array}{cc}0,& m=2,k=1\\ {\widehat{\mathrm{\Δ}}}_{\mathrm{opt},m-1},& m>2,k=1\\ {\widehat{\mathrm{\Δ}}}_{\mathrm{opt},m}^{k-1},& m\≥2,k>1\end{array}\right.$

其中，

表示第m-1次回波对应的最优偏移量；

表示第m次回波第k-1个搜索阶段对应的最优偏移量；

第三步：将第m次回波按照Δ_m，k分别进行平移后，用广义Radon变换方法重构前m次回波的图像I_m(x，y)：

$I_m(x,y)={{\∫}_0^{\mathrm{\Θ}}}^{\′}\left|s_m\right|{\mathrm{d\θ}}^{\′}$

其中，Θ′＝ωm/PRF表示前m次回波对应的目标总转角，ω为目标自旋频率，PRF为脉冲重复频率，|s_m|表示第m次回波的实包络，θ′＝ωt′_m，t′_m∈[0，m/PRF]；

第四步：对于每个偏移量，分别计算重构图像熵，使重构图像熵最小的偏移量为第m次回波第k个搜索阶段的最优偏移量

记录

第五步：如果Δ_k＞Δ_Threshold，则令搜索阶段数k＝k+1，并令Δ_k＝Δ_k-1/2，返回第二步；否则根据估计的最优偏移量对第m次回波进行平移后，令m＝m+1，重复第二步到第四步，直至m＝M，其中M为回波的方位单元数。

步骤4，算法选择

4a)采用距离-多普勒算法对补偿后的回波成像，即沿方位向进行傅里叶变换并得到目标的距离-多普勒图像；

4b)若距离-多普勒图像的方位维没有混叠，则选择CIRT算法，执行步骤5，CIRT算法公式如下：

$\hat{I}(x,y,z)={\∫}_0^{\mathrm{\Θ}}s(r,t_m)\exp \left(j\frac{4\mathrm{\π}}{\mathrm{\λ}}r\right)\mathrm{d\θ}$

其中，

为重构图像，x、y、z为所有可能散射点位置的坐标，Θ＝ωT_a为自旋目标在总观测时间T_a中的转角，ω为步骤3b)中得到的自旋频率，s(r，t_m)为平动补偿后距离-慢时间域的目标回波，r＝xsinθ+ycosθ+z为搜索变量，λ为信号载频对应的波长，θ＝ωt_m，t_m∈[0，T_a]。

4c)若得到目标的距离-多普勒图像出现方位维混叠，则选择RIRT算法，RIRT算法公式如下：

第一步，按照下式对平动补偿后的回波实包络沿距离向进行傅里叶变换：

s(ξ，t_m)＝∫|s(r，t_m)|exp(-jξr)dr

其中，s(ξ，t_m)为平动补偿后对回波实包络沿距离向进行傅里叶变换的结果，ξ为距离频率，ξ的支撑区为[-π，π]，t_m为慢时间，|s(r，t_m)|为目标回波包络的模值，r为距离，j为虚数单位；

第二步，将傅里叶变换后的回波s(ξ，t_m)按照下式进行一维滤波与逆傅里叶变换：

$s^{\′}(r,t_m)={\∫}_{-\mathrm{\π}}^{\mathrm{\π}}\left|\mathrm{\ξ}\right|s(\mathrm{\ξ},t_m)\exp \left(\mathrm{j\ξr}\right)\mathrm{d\ξ}$

其中，s′(r，t_m)为频域滤波并进行逆傅里叶变换后的回波，r为距离，t_m为慢时间，|ξ|为对ξ求模，s(ξ，t_m)为回波实包络沿距离向进行傅里叶变换后的结果；

第三步，按照下式进行后向投影，得到重构的回波图像：

$I(x,y,z)={\∫}_0^{\mathrm{\Θ}}{s}^{\′}(r^{\′},t_m){\mathrm{dt}}_m$

其中，I(x，y，z)为重构图像，x、y、z为所有可能的散射点三维坐标，Θ＝ωT_a为旋转目标在总观测时间T_a中的转角，ω为步骤3b)中得到的自旋角频率，s′(r′，t_m)为频域滤波后的回波，r′＝xcos(ωt_m)+ysin(ωt_m)+z为搜索变量，t_m为慢时间。

步骤5，二维切片成像

5a)建立自旋目标的柱面坐标系，沿目标自旋轴方向Z将自旋目标分成多个二维切片；

5b)将所有二维切片按照步骤4选择的算法进行成像，得到二维切片图像。

步骤6，重构回波；

6a)对于所有的二维切片图像，提取图像域中最大峰值对应的三维坐标，记录为目标散射点的位置坐标，并按照下式在距离-慢时间域构造散射点回波X′(r，t_m)：

$X^{\′}(r,t_m)=\sin c\left(\frac{2B}{c}(r-\hat{x}\cos \left({\mathrm{\ωt}}_m\right)-\hat{y}\sin \left({\mathrm{\ωt}}_m\right)-\hat{z})\right)$

$\×\exp (-j\frac{4\mathrm{\π}}{\mathrm{\λ}}(\hat{x}\left({\mathrm{\ωt}}_m\right)+\hat{y}\sin \left({\mathrm{\ωt}}_m\right)+\hat{z}))$

其中，X′(r，t_m)为散射点在距离-慢时间域单位幅度的回波，r为距离，t_m为慢时间，B为信号带宽，c为光速，

为散射点三维坐标的估计值，ω为步骤3b)中估计的目标自旋角频率，j为虚数单位，λ为信号载频对应的波长。

6b)运用最小范数准则对步骤6a)中所构造的散射点回波的幅度进行估计，最小范数准则的表达式为：

$\min \left|\right|E(r,t_m)-\hat{A}{X}^{\′}(r,t_m)\left|\right|=\min \underset{r,t_m}{\mathrm{\Σ}}{|E(r,t_m)-\hat{A}{X}^{\′}(r,t_m)|}^2$

其中，min(·)为求最小值运算，‖·‖为求模2-范数运算，r为距离，t_m为慢时间，E(r，t_m)为目标距离-慢时间域回波，为散射点回波幅度的估计值，X′(r，t_m)为步骤6a)中重构的散射点在距离-慢时间域单位幅度的回波。

6c)按照下式将记录的散射点的位置坐标与估计的回波幅度组合构成散射点回波：

$X(r,t_m)=\hat{A}{X}^{\′}(r,t_m)$

其中，X(r，t_m)表示信号形式，r表示方向矢量，t_m表示观测时间，

为步骤6b)中估计的散射点回波幅度值，X′(r，t_m)为步骤6a)中构造的散射点在距离-慢时间域单位幅度的回波。

6d)从步骤3d)得到的目标回波中减去重构散射点回波X(r，t)，更新目标回波E(r，t_m)。

步骤7，判断目标回波能量是否低于能量门限，能量门限为回波总能量的5％；

7a)若目标回波能量低于门限，则执行步骤8；

7b)若目标回波能量高于门限，则执行步骤5。

步骤8，操作结束，根据6a)中记录的目标散射点位置重构目标三维图像。

下面结合附图2对本发明的效果做进一步说明。

附图2所示的仿真在MATLAB7.0软件下进行的，仿真数据的参数如下：雷达PRF为1000Hz，带宽为2GHz，对应距离分辨率为7.5cm。信号载频为10GHz。数据矩阵为350×1000，观测时间为1s。图2(a)是仿真目标的三维散射点分布图，其中，三维坐标的单位均为米，9个散射点分别位于不同的旋转平面，目标自旋频率为2Hz，自旋轴与雷达视线的夹角为π/4。

图2(b)是采用基于复数后向投影的自旋目标三维成像方法得到的目标三维图像，其中，三维坐标的单位均为米，9个散射点的位置与仿真分布图中相一致，证明了算法的有效性。

Claims (2)

1.基于复数后向投影的自旋目标三维成像方法，包括如下步骤：

(1)雷达录取ISAR回波；

(2)回波相干化

2a)在距离频域对回波距离向的平移和展宽分别补偿相位因子；

2b)在慢时间域对回波补偿相位因子，使得回波完全相干化；

(3)平动补偿

3a)采用傅里叶变换对回波进行距离脉压，获得距离-慢时间域回波；

3b)对于某单次距离-慢时间域回波，分别计算其与其它各次距离-慢时间域回波的滑动相关系数，具有最大滑动相关系数的回波与该回波之间的时间间隔为目标的自旋周期，对自旋周期取倒数并乘以2π，得到自旋角频率；

3c)对应最大滑动相关系数的回波与该次回波之间移动的距离数为一个周期内目标的平动量，补偿平动量实现包络粗对齐；

3d)从第二次距离-慢时间域回波开始，对于单次距离-慢时间域回波，采用广义Radon变换方法获取不同偏移量下该次回波与之前所有回波对应的二维图像并求图像熵，使图像熵最小的距离偏移量为最优距离偏移量，根据最优偏移量对回波包络进行微调，将每次回波进行包络精对齐后，最终得到平动补偿后的目标回波；

(4)算法选择

4a)采用距离-多普勒算法对平动补偿后的回波成像，得到目标的距离-多普勒图像；

4b)若距离-多普勒图像的方位维没有混叠，则选择CIRT算法，执行步骤(5)；所述的CIRT算法公式为：

$\hat{I}(x,y,z)={\∫}_0^{\mathrm{\Θ}}s(r,t_m)\exp \left(j\frac{4\mathrm{\π}}{\mathrm{\λ}}r\right)\mathrm{d\θ}$

其中，

为重构图像，x、y、z为所有可能散射点位置的坐标，Θ＝ωT_a为自旋目标在总观测时间T_a中的转角，ω为步骤3b)中得到的自旋角频率，s(r，t_m)为平动补偿后距离-慢时间域的目标回波，r＝xsinθ+ycosθ+z为搜索变量，θ＝ωt_m，t_m∈[0，T_a]，λ为信号载频对应的波长；

4c)若距离-多普勒图像出现方位维混叠，则选择RIRT算法；所述的RIRT算法公式为：

第一步，按照下式对平动补偿后的回波实包络沿距离域进行傅里叶变换：

$s(\mathrm{\ξ},t_m)=\∫|s(r,t_m)|\exp (-\mathrm{j\ξr})\mathrm{dr}$

其中，s(ξ，t_m)为平动补偿后的回波实包络沿距离域进行傅里叶变换后的结果，ξ为距离频域，支撑区为[-π，π]，t_m为慢时间，|s(r，t_m)|为自旋目标回波的实包络，r为距离，j为虚数单位；

第二步，将s(ξ，t_m)按照下式进行一维滤波与逆傅里叶变换：

$s^{\′}(r,t_m)={\∫}_{-\mathrm{\π}}^{\mathrm{\π}}\left|\mathrm{\ξ}\right|s(\mathrm{\ξ},t_m)\exp \left(\mathrm{j\ξr}\right)\mathrm{d\ξ}$

其中，s′(r，t_m)为滤波并进行逆傅里叶变换后的回波，r为距离，t_m为慢时间，|ξ|为对ξ求模，s(ξ，t_m)为回波实包络沿距离域进行傅里叶变换后的结果；

第三步，按照下式进行后向投影，得到重构的目标图像：

$I(x,y,z)={\∫}_0^{\mathrm{\Θ}}{s}^{\′}(r^{\′},t_m){\mathrm{dt}}_m$

其中，I(x，y，z)为重构的图像，x、y为所有可能的散射点位置坐标，z为当前切片对应的坐标；Θ＝ωT_a为自旋目标在总观测时间T_a中的转角，ω为步骤3b)中得到的自旋角频率，s′(r′，t_m)为频域滤波后的回波，r′＝xcos(ωt_m)+ysin(ωt_m)+z为搜索变量，t_m为慢时间；

(5)二维切片成像

5a)建立自旋目标的柱面坐标系，沿目标自旋轴方向将自旋目标分成多个二维切片；

5b)将所有二维切片按照步骤(4)选择的算法进行成像，得到二维切片图像；

(6)重构回波并更新

6a)对于所有的二维切片图像，提取图像域中最大峰值对应的三维坐标，记录为目标散射点的位置坐标，并在距离-慢时间域构造散射点回波；

6b)运用最小范数准则对步骤6a)中所构造的散射点回波的幅度进行估计；

6c)将记录的散射点的位置坐标与估计的回波幅度组合构成散射点回波；

6d)从步骤3d)得到的目标回波中减去散射点回波，更新目标回波；

(7)判断目标回波能量是否低于能量门限

7a)若目标回波能量低于门限，则执行步骤(8)；

7b)若目标回波能量高于门限，则执行步骤(5)；

(8)操作结束，根据6a)中记录的目标散射点位置重构目标三维图像。

2.根据权利要求1所述的基于复数后向投影的自旋目标三维成像方法，其特征在于，步骤(7)中所述的能量门限为回波总能量的5％。

Images