CN103698765B - 一种isar成像方位定标方法 - Google Patents

Abstract

一种ISAR成像方位定标方法，包括以下步骤：雷达录取回波，得到以距离为列向量，以方位为行向量的回波数据；选取距离单元，组成新回波矩阵；采用频域加窗技术提取散射点回波；估计距离单元回波的相位二次项系数；基于最小二乘估计方法计算ISAR图像方位分辨率；根据ISAR图像方位分辨率，完成对ISAR图像的方位定标。本发明通过采用2阶DPT变换估计距离维回波的相位二次项系数，弥补了传统方法基于LPFT和图像对比度最大化进行相位二次项系数估计时计算量过大，不便于实时应用的缺点，并通过最小二乘方法估计得到目标的转动角速度，从而实现ISAR图像的方位向定标。

Description

一种ISAR成像方位定标方法

技术领域

本发明属于信号处理技术领域，尤其涉及一种雷达信号处理中的逆合成孔径雷达成像方位定标方法。

背景技术

逆合成孔径雷达（Inverse Synthetic Aperture Radar，以下简称ISAR）具有全天时、全天候、远距离、高分辨等优点，在军事和民用领域具有广泛的应用。为了确定目标的真实尺寸，进而对目标进行识别和分类，需要对ISAR图像进行距离和方位定标。距离维的定标仅与雷达发射信号带宽以及采样频率有关，可事先获得，但是方位维的定标处理与目标相对于雷达的转角有关。由于ISAR成像的目标一般为非合作目标，通常无法事先准确得到其相关运动参数，所以需要对ISAR成像目标的转动角速度进行估计，进而通过转动角度最终实现对ISAR成像的方位定标。

李宁等人在文献“一种改进的ISAR图像方位向定标方法”(雷达科学与技术，2012年)中，根据ISAR距离维回波信号二次项系数与距离维尺寸的关系，首先用最大对比度和质量评估准则筛选出多个质量最好的包含特显点的距离单元，然后提取出这些单元内特显点的回波信号并估计各特显点的调频率，最后由特显点的距离坐标和调频率获取目标的转动角速度，从而实现ISAR图像的方位向定标。但是上述方法存在的不足是，在对含有特显点的回波信号进行调频率估计时，采用局部多项式相位变换(LPFT)并结合对比度最大准则的方法，因此要想获得很高的估计精度，需要设置很小的搜索间隔，在解空间内反复搜索，从而带来巨大的计算量，不便于实时应用。

发明内容

针对现有技术存在的不足，本发明的目的是提供一种计算量小、定标精度高、便于实时应用的基于离散多项式变换的ISAR成像方位定标方法。

为了实现上述目的，本发明采取如下的技术解决方案：

一种ISAR成像方位定标方法，包括以下步骤：

步骤1、雷达录取回波，得到以方位为行向量、以距离为列向量的回波数据S_raw，回波数据S_raw的行数为M、列数为N；

步骤2、选取距离单元，组成新回波矩阵S_p，具体步骤如下：

步骤2a、计算回波数据S_raw中距离单元回波的归一化幅度方差：

${\mathrm{\δ}}_m=1-\frac{{\left[E\left({|s}_m\right|)\right]}^2}{E\left({{|s}_m|}^2\right)},$

δ_m为回波数据S_raw中第m行距离单元回波的归一化幅度方差，E(·)表示求期望运算操作，s_m表示回波数据S_raw中第m行距离单元回波的幅度值，m=0,1,…,M-1；

步骤2b、将所有距离单元回波的归一化幅度方差值按照从小到大的顺序进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为归一化幅度方差排序位置矩阵Q₁；

步骤2c、对回波数据S_raw的每行进行FFT变换，将FFT变换之后的矩阵每行按幅度的绝对值从大到小进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为幅度绝对值排序位置矩阵Q₂；

步骤2d、将Q₁和Q₂相加，并按照从小到大的顺序进行排序，选取前P行距离单元回波组成新回波矩阵S_p，并记录新回波矩阵S_p中各行距离单元回波在回波数据S_raw中的原行号，记为Q，新回波矩阵S_p的行数为P、列数为N；

步骤3、采用频域加窗技术提取散射点回波；

步骤3a、对新回波矩阵S_p的每行进行FFT变换，变换到方位频域，得到距离时域方位频域数据矩阵S_PF；

步骤3b、对距离时域方位频域数据矩阵S_PF每行的幅度值进行搜索，确定幅度最大值的位置，以该幅度最大值的位置为中心添加矩形窗，矩形窗以外的元素置为零，提取出散射点频谱，执行下一步；

步骤3c、对加窗之后的距离时域方位频域数据矩阵S_PF的每行进行IFFT变换，变换到方位时域，得到距离时域方位时域数据矩阵S_PD；

步骤4、估计距离单元回波的相位二次项系数：

步骤4a、对步骤3c得到的距离时域方位时域数据矩阵S_PD的每行进行2阶DPT运算，得到数据矩阵Z：

$z_k\left(n\right)={\mathrm{DPT}}_2[s_k\left(n\right),\mathrm{\τ}]=x_k\left(n\right)\×x_k^{*}(n-\mathrm{\τ}),$

其中，z_k(n)表示数据矩阵Z中第k行的第n个元素，DPT₂[·]表示2阶DPT运算操作，s_k(n)表示步骤3c得到的距离时域方位时域数据矩阵S_PD中第k行的第n个元素，τ表示信号时延，(·)^*表示共轭操作，k=0,1,…P-1，n=0,1,…N-1；

步骤4b、对数据矩阵Z的每行进行FFT变换，得到方位频域数据矩阵Z_F，执行下一步；

步骤4c、对方位频域数据矩阵Z_F每行的幅度值进行搜索，确定每行幅度值最大的位置，记录为G=[g₀ g₁ …g_k… g_P-1]^T，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，执行下一步；

步骤4d、计算方位频域数据矩阵Z_F中每行距离单元回波的相位二次项系数：

$a_k=\frac{\mathrm{\π}\×{\mathrm{PRF}}^2\×(g_k-\frac{N}{2}-1)}{N\×\mathrm{\τ}},$

a_k表示方位频域数据矩阵Z_F中第k行距离单元回波的相位二次项系数，PRF表示雷达发射信号的脉冲重复频率，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，N为回波数据S_raw的列数，τ表示信号时延；

步骤5、基于最小二乘估计方法计算ISAR图像方位分辨率P_a；

步骤5a、计算ISAR图像距离分辨率其中，C为电磁波传播速度，F_s为雷达接收信号的采样频率，γ为雷达发射信号的调频率，M为回波数据S_raw的行数；

步骤5b、计算距离维回波的相位二次项系数随距离的变化率e=LSE₁[u,v]，LSE₁[u,v]表示以v为自变量，u为函数值进行一阶线性拟合，其中，u=[a₀ a₁ …a_k… a_P-1]^T构成相位二次项系数的列向量，v=Q·p_r，Q为新回波矩阵S_p中各行距离单元回波分别在回波数据S_raw中的原行号，P_r为ISAR图像距离分辨率；

步骤5c、计算ISAR目标转动角速度其中，e为距离维回波的相位二次项系数随距离的变化率，λ为雷达发射信号的波长；

步骤5d、计算ISAR图像方位分辨率其中，λ为雷达发射信号的波长，PRF表示雷达发射信号的脉冲重复频率，Ω为ISAR目标转动角速度，N为回波数据S_raw的列数；

步骤6、根据ISAR图像方位分辨率P_a，完成对ISAR图像的方位定标。

进一步的，所述步骤1中的回波数据S_raw通过解线调频、包络对齐和方位自聚焦处理得到。

进一步的，所述步骤2b中新回波矩阵S_p的行数P为8～10。

进一步的，所述步骤3b中的矩形窗长度为L+1，L=N/16，N为新回波矩阵S_p的列数。

由以上技术方案可知，本发明基于距离回波的归一化幅度方差和峰值能量值选取含有较强散射点的距离单元，采用频域加窗技术提取散射点回波，然后采用2阶DPT估计选取的距离回波相位二次项系数，最终基于最小二乘法估计得到目标转动角速度，从而实现对ISAR图像方位维定标。与现有技术相比，本发明通过采用2阶DPT变换估计距离维回波的相位二次项系数，弥补了传统方法基于LPFT和图像对比度最大化进行相位二次项系数估计时计算量过大，不便于实时应用的缺点，并通过最小二乘方法估计得到目标的转动角速度，从而实现ISAR图像的方位向定标，对目标转动的角速度估计准确，计算高效。

附图说明

图1为本发明实施例的方法的流程图；

图2为仿真实验的目标模型图；

图3为仿真实验选取的10行距离单元的相位二次系数与其相应的距离位置（m）的最小二乘拟合结果图；

图4为进行100次蒙塔卡洛仿真后的估计相对误差结果图；

图5为使用估计出的角速度对ISAR图像定标的结果；

图6为估计出的转动角速度相对误差随信噪比的变化曲线。

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明实施例，并配合所附图示，做详细说明如下。

本发明方法的基本思路是：首先，雷达录取回波，对录取到的ISAR原始回波数据进行解线调频、包络对齐和方位自聚焦得到距离压缩后的回波数据；其次，计算各个距离单元回波的归一化幅度方差和峰值能量，并进行排序，选取出一组归一化幅度方差较小和峰值能量较大的距离单元回波；然后，采用频域加窗技术提取强散射点回波；之后，对上述距离单元回波进行2阶离散多项式变换（DiscretePolynomial-PhaseTransform，以下简称DPT），估计得到距离单元回波的相位二次项系数；最后，采用最小二乘法估计得到ISAR目标的转动角速度，最终实现对ISAR的方位定标处理。

以上是本发明的核心思想，下面将结合本发明实施例中的附图，对本发明实施例的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下获得的所有其他实施例，都属于本发明保护的范围。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其它不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

参照图1，图1为本发明方法的流程图，本发明方法包括以下步骤：

步骤1、雷达录取回波，通过解线调频、包络对齐和方位自聚焦处理，得到以方位为行向量、以距离为列向量的回波数据S_raw，回波数据S_raw的行数为M、列数为N；

步骤2、选取距离单元，组成新回波矩阵S_p，具体步骤如下：

步骤2a、计算回波数据S_raw中距离单元回波的归一化幅度方差：

${\mathrm{\δ}}_m=1-\frac{{\left[E\left({|s}_m\right|)\right]}^2}{E\left({{|s}_m|}^2\right)},$

δ_m为回波数据S_raw中第m行距离单元回波的归一化幅度方差，E(·)表示求期望运算操作，s_m表示回波数据S_raw中第m行距离单元回波的幅度值，m=0,1,…,M-1；

步骤2b、将所有距离单元回波的归一化幅度方差值按照从小到大的顺序进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为归一化幅度方差排序位置矩阵Q₁；

步骤2c、对回波数据S_raw的每行进行FFT变换，将FFT变换之后的矩阵每行按幅度的绝对值从大到小进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为幅度绝对值排序位置矩阵Q₂；

步骤2d、将Q₁和Q₂相加，并按照从小到大的顺序进行排序，选取前P行距离单元回波组成新回波矩阵S_p，并记录新回波矩阵S_p中各行距离单元回波在步骤1得到的回波数据S_raw中的原行号，记为Q，新回波矩阵S_p的行数为P，列数为N（与回波数据S_raw列数相同）；更具体的，P可为8～10；

步骤3、采用频域加窗技术提取散射点回波；

步骤3a、对新回波矩阵S_p的每行进行FFT变换，变换到方位频域，得到距离时域方位频域数据矩阵S_PF；

步骤3b、对得到距离时域方位频域数据矩阵S_PF每行的幅度值进行搜索，确定幅度最大值的位置，以该幅度最大值的位置为中心，添加窗长度为L+1的矩形窗，矩形窗以外的元素置为零，提取出散射点频谱，执行下一步；更具体的，L=N/16，N为新回波矩阵S_p的列数；

步骤3c、对加窗之后的距离时域方位频域数据矩阵S_PF的每行进行IFFT变换，变换到方位时域，得到距离时域方位时域数据矩阵S_PD；

步骤4、估计距离单元回波的相位二次项系数：

步骤4a、对步骤3c得到的距离时域方位时域数据矩阵S_PD的每行进行2阶DPT运算，得到数据矩阵Z：

$z_k\left(n\right)={\mathrm{DPT}}_2[s_k\left(n\right),\mathrm{\τ}]=x_k\left(n\right)\×x_k^{*}(n-\mathrm{\τ}),$

其中，z_k(n)表示数据矩阵Z中第k行的第n个元素，DPT₂[·]表示2阶DPT运算操作，s_k(n)表示步骤3c得到的距离时域方位时域数据矩阵S_PD中第k行距离单元的第n个元素，τ表示信号时延，(·)^*表示共轭操作，k=0,1,…P-1，n=0,1,…N-1；更具体的，τ=N/2；

步骤4b、对数据矩阵Z的每行进行FFT变换，得到方位频域数据矩阵Z_F，执行下一步；

步骤4c、对方位频域数据矩阵Z_F每行的幅度值进行搜索，确定每行幅度值最大的位置，记录为G=[g₀ g₁ …g_k… g_P-1]^T，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，执行下一步；

步骤4d、计算方位频域数据矩阵Z_F中每行距离单元回波的相位二次项系数：

$a_k=\frac{\mathrm{\π}\×{\mathrm{PRF}}^2\×(g_k-\frac{N}{2}-1)}{N\×\mathrm{\τ}},$

a_k表示方位频域数据矩阵Z_F中第k行距离单元回波的相位二次项系数，PRF表示雷达发射信号的脉冲重复频率，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，N为回波数据S_raw的列数，τ表示信号时延；

步骤5、基于最小二乘估计方法计算ISAR图像方位分辨率P_a；

步骤5a、计算ISAR图像距离分辨率其中，C为电磁波传播速度，F_s为雷达接收信号的采样频率，γ为雷达发射信号的调频率，M为回波数据S_raw的行数；

步骤5b、计算距离维回波的相位二次项系数随距离的变化率e=LSE₁[u,v]，LSE₁[u,v]表示以v为自变量，u为函数值进行一阶线性拟合，其中，u=[a₀ a₁ …a_k… a_P-1]^T构成相位二次项系数的列向量，v=Q·p_r，Q为新回波矩阵S_p中各行距离单元回波分别在回波数据S_raw中的原行号，P_r为ISAR图像距离分辨率；

步骤5c、计算ISAR目标转动角速度其中，e为距离维回波的相位二次项系数随距离的变化率，λ为雷达发射信号的波长；

步骤5d、计算ISAR图像方位分辨率其中，λ为雷达发射信号的波长，PRF表示雷达发射信号的脉冲重复频率，Ω为ISAR目标转动角速度，N为回波数据S_raw的列数；

步骤6、根据ISAR图像方位分辨率P_a，完成对ISAR图像的方位定标。

与现有技术相比，本发明根据ISAR距离压缩后信号相位随方位维的变化特点，通过2阶DPT估计不同距离单元回波信号的相位二次项系数，进而通过最小二乘估计得到ISAR目标转动角速度的方位定标方法，克服了现有方法需要在解空间内反复搜索估计相位二次项系数计算量大的问题，运算效率高；同时，本方利用距离单元回波的归一化幅度方差和频域加窗技术来选取距离单元回波并提取散射点回波，克服了现有方法没有有效的距离单元回波选取准则的缺点，并且该方法能够有效提高目标转动角速度估计精度，提高ISAR方位维定标精度。

本发明的效果可以通过以下的仿真实验进一步说明，仿真时采用MATLAB（R2010b）软件进行仿真。

仿真条件如下：

雷达发射信号为X波段，波长为0.03m，发射信号带宽为1GHz，脉冲宽度为50μs，重复频率为1000Hz，采样频率为10MHz，距离维取512点，方位维取512点。

仿真内容

如图2所示，图2为仿真目标的模型图，对天线长度为27.21m，太阳能电池板长为47.70m的卫星点目标进行仿真，图2中横轴为方位维，纵轴为距离维，单位均为m。仿真时目标等效转动角速度为0.2rad/s。在得到距离压缩后的距离频域方位时域数据之后，将回波加上信噪比为12dB的零均值复高斯白噪声。

选取归一化幅度方差最小和峰值能量最大的10行距离单元回波，使用本发明方法提出的基于DPT的转动角速度估计算法对角速度进行估计，图4所示为进行100次蒙塔卡洛仿真后的估计相对误差结果。图4中横轴为仿真次数，纵轴为估计相对误差，从图4可看出，角速度估计的相对误差均小于0.01，估计精度很高。

用估计出的角速度对目标进行ISAR图像定标，定标结果如图5所示，图5中横轴为方向维，纵轴为距离维，单位均为m。定标之后卫星的天线长度为27.17m，太阳能电池板长为47.21m，误差分别为-0.15%和-1.03%，可看出定标后的ISAR图像其距离维尺度和方位维尺度与仿真点目标模型基本一致。

图6为取自10行距离单元估计出的转动角速度相对误差随信噪比的变化曲线。图6中横轴为信噪比（单位dB），纵轴为估计相对误差。从图6可以看出，当信噪比大于-5dB时，转动角速度的估计相对误差就已经趋于稳定，估计精度很高，说明本发明提出的方法有着一定的稳健性。

上述实验结果验证了本方法提出的基于DPFT的均匀转动目标转动角速度估计方法的有效性和稳定性。

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (4)

1.一种ISAR成像方位定标方法，其特征在于，包括以下步骤：

步骤1、雷达录取回波，得到以距离为列向量，以方位为行向量的回波数据S_raw，行数为M、列数为N；

步骤2、选取距离单元，组成新回波矩阵S_p，具体步骤如下：

步骤2a、计算回波数据S_raw中距离单元回波的归一化幅度方差：

${\mathrm{\δ}}_m=1-\frac{{\left[E\left(\right|s_m\left|\right)\right]}^2}{E\left({\left|s_m\right|}^2\right)},$

δ_m为回波数据S_raw中第m行距离单元回波的归一化幅度方差，E(·)表示求期望运算操作，s_m表示回波数据S_raw中第m行距离单元回波的幅度值，m＝0,1,…,M-1；

步骤2b、将所有距离单元回波的归一化幅度方差值按照从小到大的顺序进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为归一化幅度方差排序位置矩阵Q₁；

步骤2c、对回波数据S_raw的每行进行FFT变换，将FFT变换之后的矩阵每行按幅度的绝对值从大到小进行排序，并记录回波数据S_raw中各行在排序后的矩阵中的位置，记为幅度绝对值排序位置矩阵Q₂；

步骤2d、将Q₁和Q₂相加，并按照从小到大的顺序进行排序，选取前P行距离单元回波组成新回波矩阵S_p，并记录新回波矩阵S_p中各行距离单元回波在回波数据S_raw中的原行号，记为Q，新回波矩阵S_p的行数为P、列数为N；

步骤3、采用频域加窗技术提取散射点回波；

步骤3a、对新回波矩阵S_p的每行进行FFT变换，变换到方位频域，得到距离时域方位频域数据矩阵S_PF；

步骤3b、对距离时域方位频域数据矩阵S_PF每行的幅度值进行搜索，确定幅度最大值的位置，以该幅度最大值的位置为中心添加矩形窗，矩形窗以外的元素置为零，提取出散射点频谱，执行下一步；

步骤3c、对加窗之后的距离时域方位频域数据矩阵S_PF的每行进行IFFT变换，变换到方位时域，得到距离时域方位时域数据矩阵S_PD；

步骤4、估计距离单元回波的相位二次项系数：

步骤4a、对步骤3c得到的距离时域方位时域数据矩阵S_PD的每行进行2阶DPT运算，得到数据矩阵Z：

$z_k\left(n\right)={\mathrm{DPT}}_2[s_k\left(n\right),\mathrm{\τ}]=x_k\left(n\right)\×x_k^{*}(n-\mathrm{\τ}),$

其中，z_k(n)表示数据矩阵Z中第k行的第n个元素，DPT₂[·]表示2阶DPT运算操作，s_k(n)表示步骤3c得到的距离时域方位时域数据矩阵S_PD中第k行的第n个元素，τ表示信号时延，(·)^*表示共轭操作，k＝0,1,…P-1，n＝0,1,…N-1；

步骤4b、对数据矩阵Z的每行进行FFT变换，得到方位频域数据矩阵Z_F，执行下一步；

步骤4c、对方位频域数据矩阵Z_F每行的幅度值进行搜索，确定每行幅度值最大的位置，记录为G＝[g₀ g₁ …g_k… g_P-1]^T，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，执行下一步；

步骤4d、计算方位频域数据矩阵Z_F中每行距离单元回波的相位二次项系数：

$a_k=\frac{\mathrm{\τ}\×{\mathrm{PRF}}^2\×(g_k-\frac{N}{2}-1)}{N\×\mathrm{\τ}},$

a_k表示方位频域数据矩阵Z_F中第k行距离单元回波的相位二次项系数，PRF表示雷达发射信号的脉冲重复频率，g_k为方位频域数据矩阵Z_F中第k行距离单元回波的幅度值最大值位置，N为回波数据S_raw的列数，τ表示信号时延；

步骤5、基于最小二乘估计方法计算ISAR图像方位分辨率P_a；

步骤5a、计算ISAR图像距离分辨率其中，C为电磁波传播速度，F_s为雷达接收信号的采样频率，γ为雷达发射信号的调频率，M为回波数据S_raw的行数；

步骤5b、计算距离维回波的相位二次项系数随距离的变化率e＝LSE₁[u,v]，LSE₁[u,v]表示以v为自变量，u为函数值进行一阶线性拟合，其中，u＝[a₀ a₁ …a_k… a_P-1]^T构成相位二次项系数的列向量，v＝Q·p_r，Q为新回波矩阵S_p中各行距离单元回波分别在回波数据S_raw中的原行号，P_r为ISAR图像距离分辨率；

步骤5c、计算ISAR目标转动角速度其中，e为距离维回波的相位二次项系数随距离的变化率，λ为雷达发射信号的波长；

步骤5d、计算ISAR图像方位分辨率其中，λ为雷达发射信号的波长，PRF表示雷达发射信号的脉冲重复频率，Ω为ISAR目标转动角速度，N为回波数据S_raw的列数；

步骤6、根据ISAR图像方位分辨率P_a，完成对ISAR图像的方位定标。

2.根据权利要求1所述的ISAR成像方位定标方法，其特征在于：所述步骤1中的回波数据S_raw通过解线调频、包络对齐和方位自聚焦处理得到。

3.根据权利要求1或2所述的ISAR成像方位定标方法，其特征在于：所述步骤2d中新回波矩阵S_p的行数P为8～10。

4.根据权利要求1或2所述的ISAR成像方位定标方法，其特征在于：所述步骤3b中的矩形窗长度为L+1，L＝N/16，N为新回波矩阵S_p的列数。