CN110133648A - 一种选取逆合成孔径雷达船只成像时窗的方法 - Google Patents

Abstract

本发明公开了一种选取逆合成孔径雷达船只成像时窗的方法，属于雷达技术领域，本发明通过理论分析和数值仿真，验证了散射点的瞬时多普勒频率仅在较短的时间内才近似满足线性调频信号的变化规律，推导了上述短时间段长度的计算公式，针对短时线性调频信号的调频斜率和中心频率估计，设计了一种数据外推与压缩感知相结合的参数估计方法，通过相邻短时间段的数据关联，得到完整数据采集时间内主要强散射点的瞬时多普勒频率变化，综合上述信息，选取出这些强散射点的瞬时多普勒频率变化都比较平稳的时间段。本发明仿真和实测数据的处理结果验证了所提出方法的有效性，通过与基于图像对比度的(Image Contrast Based Algorithm，ICBA)方法进行比较，本发明方法取得了更好的成像质量。

Description

一种选取逆合成孔径雷达船只成像时窗的方法

技术领域

本发明属于雷达技术领域，具体涉及一种选取逆合成孔径雷达船只成像时窗的方法。

背景技术

利用逆合成孔径雷达(Inverse Synthetic Aperture Radar，ISAR)对海面船只进行成像，是对船只实现分类和识别的重要技术途径，在遥感和军事等领域具有广泛的应用。在海面上行驶的船只，受海浪等因素的影响，在横滚、俯仰、偏航三个轴向上均存在着转动，使得船只的运动规律异常复杂。距离多普勒(Range Doppler，RD)算法是ISAR系统中最为常用的算法，但该算法要求在成像积累时间内船只的转动必须保持平稳。因此，时间加窗处理成为ISAR船只成像过程中的重要环节之一。利用时窗选取算法寻找到船只转动相对平稳的时间段，从而仅对该时间段内的回波进行RD成像，可获得较好的成像质量。

目前已有的ISAR时窗选取算法，主要可分为以下几类：基于散射点追踪的算法、基于船只转速和多普勒展宽估计的算法、基于图像对比度的算法、基于相位非线性度的算法、适合于合作式船只的算法、基于多普勒频率信息的算法。

基于散射点追踪的算法，通过追踪位于船首、船尾的2个散射点搜索适合得到侧视图的成像时刻，通过追踪位于上层建筑某处、该处在甲板投影位置的2个散射点搜索适合得到俯视图的成像时刻；该算法的缺点是在海情较高的情况下很难准确跟踪上述4个特定散射点。

基于船只转速和多普勒展宽估计的算法，首先进行分时段成像得到多幅距离多普勒图像，然后通过Radon变换等方法在图像上获得船身斜率随时间的变化，进而近似得到船只有效转动矢量的垂直分量随时间变化的曲线；再从多幅距离多普勒图像上提取多普勒展宽信息，得到其随时间变化的曲线；最后，利用上述2条曲线确定适合得到侧视图或俯视图的成像时刻；该算法的缺点是在海情较高的情况下所得的多幅分时段距离多普勒图像的质量可能较低，从而会影响船只有效转动矢量的垂直分量和多普勒展宽信息提取的准确性。

基于图像对比度的算法，一般又称为ICBA(Image Contrast Based Algorithm)算法，该算法利用图像对比度最大的准则，先对时窗的中心时刻进行一维优化，再对时窗的长度进行一维优化；通过该算法获得的图像具有最高的对比度，因而具有较高的成像质量；但该方法也具有计算效率不高的缺点；另外，图像对比度也并不是衡量船只转动平稳性的直接评价指标。

基于相位非线性度的算法，所依据的原理是在船只转动相对平稳的时间段内两个散射点的相位之间具有较好的线性关系，因此通过计算多个强散射体各自所在距离单元的平均相位非线性度，可寻找到计算出的平均相位非线性度低于设定阈值的时间段，并将这些时间段选为适宜成像的时间段；该算法的缺点是有些情况下并不能找到强散射体占据了主导地位的多个距离单元。

一种仅适合于合作式船只的算法，要求目标上必须安装有运动传感器，用于获取目标运动数据，该算法对于非合作目标不适用。

基于多普勒频率信息的算法，通过估计散射点的转动造成的瞬时多普勒频率变化，寻找到频率变化比较平稳的时间段，并将这些时间段选为适宜成像的时间段。该类算法的优点是多普勒频率变化是评估船只转动平稳性的直接指标，因而所得图像具有较高的成像质量。目前，这类算法的主要不足是要么仅考虑了单个散射点的瞬时多普勒频率变化，要么是用船只整体的瞬时多普勒频率变化代表多个散射体的平均多普勒频率变化。实际上，船只上主要强散射点的瞬时多普勒频率变化并不相同，应寻找这些强散射点的瞬时多普勒频率变化都比较平稳的时间段。

针对上述不足，本申请提出一种基于多散射点瞬时多普勒频率估计的算法，用于选取ISAR船只成像过程中的适宜时间窗。通过理论分析和数值仿真，验证了船只转动引起的瞬时多普勒频率仅在较短的时间内近似服从线性调频信号的规律。推导了上述短时间段的时间长度计算公式。设计了数据外推与压缩感知相结合的短时线性调频信号调频斜率和中心频率的估计方法，保证了即使在信号长度较短、从而频率分辨率较低的情况下，依然具有较高的调频斜率和中心频率估计精度。得到各个短时间段内主要强散射点的瞬时多普勒频率变化以后，对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到这些强散射点在完整数据采集时间内的瞬时多普勒频率变化。之后，综合这些强散射点的瞬时多普勒频率变化信息，选取出这些强散射点的瞬时多普勒频率变化都比较平稳的时间段。分别利用仿真和实测数据，对所提出方法的有效性进行了验证，并与ICBA算法的成像效果进行了对比，本申请算法取得了更好的成像质量。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种逆合成孔径雷达船只成像时窗选取方法，设计合理，克服了现有技术的不足，具有良好的效果。

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

一种选取逆合成孔径雷达船只成像时窗的方法，包括如下步骤：

步骤1：对逆合成孔径雷达回波信号进行包括距离压缩、平动补偿在内的预处理；

将完整的信号采集时间划分为若干个短时间段，在每个短时间段内散射点瞬时多普勒频率的变化如公式(8)所示：

其中，t表示方位向时间，f_dk(t)表示第k个散射点的瞬时多普勒频率，i表示时间段的序号，N表示所有时间段的个数，K_rki表示第k个散射点在第i个时间段内的线性调频斜率，t_start表示信号采集时间的起始时刻，T_seg表示每个时间段的时间长度，f_dcki表示第k个散射点在第i个时间段内的中心频率，rect(t/T_seg)为[-T_seg/2,T_seg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数；

在逆合成孔径雷达系统中，方位时间-距离频率域的船只回波信号表示如公式(13)所示：

其中，f_r表示距离向频率，S(t,f_r)表示方位时间-距离频率域的船只回波信号，T_obs为信号的数据采集时间，f_c为雷达载波频率，B_r为雷达信号带宽，C表示光速，R₀(t)表示从雷达的相位中心到船只转动中心的矢量，R₀(t)为R₀(t)的模，i_R0(t)为R₀(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ₀(r)为位于r处散射点的归一化RCS(Radar Cross Section，雷达后向散射截面积)，A(r,t)为考虑发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点RCS的时变性等因素对信号幅度和相位造成影响的复函数；

经过包括距离对齐、相位补偿在内的平动补偿处理后的信号S′(t,f_r)表示如公式(14)所示：

其中，R′₀为成像中心时刻散射点与雷达间的距离；

对式(14)离散化后的信号进行距离向IFFT(Inverse Fast FourierTransformation，快速逆傅里叶变换)，得到方位时间-距离时间域的信号s′(t,τ)，如公式(15)所示：

其中，τ表示距离向时间，K为散射点的数目，A_k(t)为第k个散射点对应的A(r,t)，σ_k为第k个散射点的RCS，τ_k为第k个散射点相对于雷达的双程延迟时间，r_k为从船只转动中心到第k个散射点的矢量，R′_k为r_k在成像中心时刻在雷达视线上的投影长度，相位项包含着散射点k的多普勒频率信息；

将公式(15)进行简化，得到公式(16)：

步骤2：从预处理后的ISAR回波信号中，挑选出L个功率最强的距离单元，并将这些距离单元按功率从大到小从按照1至L编上编号；

步骤3：以T_seg为时间段长度，在方位向上将步骤2中L个距离单元的信号划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；

步骤4：设j为L个功率最强的距离单元中当前单元的编号，并令j的初值为1；

步骤5：设l_j为第j个距离单元在所有距离单元中的序号，表示序号为l_j的距离单元对应的一维信号，表示中第i个时间段的信号，利用经典的Burg算法，对进行数据外推，得到

其中，表示序号为l_j的距离单元中所包含的散射点的数目，表示序号为l_j的距离单元所对应的距离向时间；

将式(8)代入式(17)后进行简化，得到如下近似表达式：

式中，表示复函数A_k(t)中的相位，为式(17)中积分项所对应的初始相位；因此，有：

步骤6：将按PRI(Pulse Repetition Interval，脉冲重复间隔)为采样周期进行离散化，并将离散化后的信号表示成矩阵形式，将该矩阵乘以随机采样过程对应的测量矩阵Φ得到如公式(22)所示：

其中，为的矩阵形式，m为方位向采样点的序号，M为利用经典的Burg算法进行数据外推后对应的方位向采样点数，M的值取偶数，0≤m≤M-1，Φ为随机采样过程对应的测量矩阵，为对进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M；PRI为脉冲重复间隔；

步骤7：利用OMP(orthogonal matched pursuit，正交匹配追踪)重建算法，根据公式(26)，估计出稀疏系数矢量

信号又能表示为：

其中，W_M＝e^-j2π/M，θ为稀疏系数矢量；

将式(24)代入式(22)，得到如下表达式：

其中，T为感知矩阵；

其中，||·||₁表示l₁范数，||·||₂表示l₂范数，ε为与噪声电平相关的常数；

步骤8：估计K_rki和f_dcki的值，其中下标k的取值范围为从1至

将序号为l_j的距离单元在第i个时间段所估计出的散射点数目记为将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点，等于矩阵中有效局部峰值点的个数；

K_rki和f_dcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式如式(27)-式(30)所示：

n_kic＝n_ki-n_kirM(28)；

其中，n_ki为第k个散射点第i个时间段对应的局部峰值系数在矢量θ中的序号，0≤n_ki≤M²-1，floor(·)为取自变量整数部分的函数，n_kir为对应的调频率序号，0≤n_kir≤M-1，n_kic为对应的多普勒中心频率序号，0≤n_kic≤M-1，PRF(Pulse Repetition Frequency)为脉冲重复频率；K_rki为散射点k第i个时间段对应的调频斜率，f_dcki为散射点k第i个时间段对应的多普勒中心频率；

步骤9：判断j是否等于L；

若：判断结果是j＝L；则执行步骤10；

或判断结果是j≠L；则令j＝j+1，然后执行步骤5；

步骤10：判断i是否等于N；

若：判断结果是i＝N；则执行步骤11；

或判断结果是i≠N；则令i＝i+1，然后执行步骤4；计算出L个最强距离单元所有时间段对应的K_rki和f_dcki；

步骤11：从第2个时间段开始，至第N个时间段结束，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联；

其中，f_dcli+K_rli·T_seg/2为第l个散射点在第i个时间段结束时刻的瞬时多普勒频率，1≤i≤N，为序号为l_j的距离单元在第i个时间段估计出的散射点数；为第l′个散射点在第i+1个时间段开始时刻的瞬时多普勒频率， 为序号为l_j的距离单元在第i+1个时间段估计出的散射点数；K_rl′(i+1)为散射点l′第i+1个时间段对应的调频斜率；K_rli为散射点l第i个时间段对应的调频斜率；取1≤i≤N中的最小值为

步骤12：利用式(32)，计算完整采集时间内各个离散时刻t_exp允许的积累时间长度Δt(t_exp)；

其中，t_exp＝m·PRI,0≤m≤M-1，为序号为l_j的距离单元中第l个散射点对应的搜索终止时刻与t_exp时刻之间的时间差，这里的搜索终止时刻记为t_exp时刻之后的第1个具有如下特点的时刻：序号为l_j的距离单元中第l个散射点在该时刻的瞬时多普勒频率与该散射点在t_exp时刻瞬时多普勒频率之差的绝对值大于期望的多普勒分辨率ρ_f；这里需要说明的是，当搜索过程中遇到完整采集时间内的终止时刻t_start+N·T_seg时，搜索过程自然终止，将搜索终止时刻记为t_start+N·T_seg；

步骤13：利用式(33)和式(34)，计算最终的最佳成像起始时刻t_opt,start和时窗长度Δt_opt；

Δt_opt＝Δt(t_opt,start) (34)；

其中，Δt(t_opt,start)表示Δt函数在t_opt,start时刻的函数值；

步骤14：结束。

本发明所带来的有益技术效果：

本发明考虑了船只上主要强散射点的瞬时多普勒频率变化并不相同，通过估计船只上多个强散射点的瞬时多普勒频率变化，综合多散射点的瞬时多普勒频率信息进行ISAR系统的时窗选取，对时窗长度和中心时刻的选取更加准确，从而可取得聚焦效果更好的图像，并进而提高船只分类和识别的准确度。

附图说明

图1为船只坐标系示意图；

图2为船只散射点的分布图；

图2(a)为驱逐舰的散射点分布示意图；

图2(b)为航空母舰的散射点分布示意图；

图3为船只散射点的瞬时多普勒频率变化示意图；

图3(a)为驱逐舰散射点的瞬时多普勒频率变化示意图；

图3(b)为航空母舰散射点的瞬时多普勒频率变化示意图；

图4为本发明方法的流程图；

图5为单散射点条件下所估计出的瞬时多普勒频率曲线示意图；

图5(a)为第1种情况示意图；图5(b)为第2种情况示意图；

图6为多散射点条件下所估计出的瞬时多普勒频率曲线示意图；

图6(a)为第1个强距离单元的2个强散射点的示意图；图6(b)为第2个强距离单元的2个强散射点的示意图；

图7为选中时窗内强散射点的瞬时多普勒频率示意图；

图7(a)为本发明方法选中时窗内强散射点的瞬时多普勒频率示意图；

图7(b)为ICBA算法选中时窗内强散射点的瞬时多普勒频率示意图；

图8为船只仿真成像结果示意图；

图8(a)为本发明方法的船只仿真成像结果示意图；8(b)为ICBA算法的船只仿真成像结果示意图；

图9为实测数据成像结果示意图；

图9(a)为不加窗处理的实测数据成像结果示意图；

图9(b)为利用本发明方法进行加窗的实测数据成像结果示意图；

图9(c)为利用ICBA算法进行加窗的实测数据成像结果示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

1、本文提出一种基于多散射点瞬时多普勒频率估计的算法，用于选取ISAR船只成像过程中的适宜时间窗。通过理论分析和数值仿真，验证了船只转动引起的瞬时多普勒频率仅在较短的时间内近似服从线性调频信号的规律。推导了上述短时间段的时间长度计算公式。设计了数据外推与压缩感知相结合的短时线性调频信号调频斜率和中心频率的估计方法，保证了即使在信号长度较短、从而频率分辨率较低的情况下，依然具有较高的调频斜率和中心频率估计精度。得到各个短时间段内主要强散射点的瞬时多普勒频率变化以后，对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到这些强散射点在完整数据采集时间内的瞬时多普勒频率变化。之后，综合这些强散射点的瞬时多普勒频率变化信息，选取出这些强散射点的瞬时多普勒频率变化都比较平稳的时间段。分别利用仿真和实测数据，对所提出方法的有效性进行了验证，并与ICBA算法的成像效果进行了对比，本申请算法取得了更好的成像质量。

2、船只转动引起的多普勒频率特性分析

在ISAR系统中，目标的平动分量对成像没有贡献，在信号处理过程中由平动补偿环节完成补偿。因此只分析船只转动引起的多普勒频率的特性。

首先建立如图1所示的船只坐标系。该坐标系以零时刻的船只转动中心为原点o，以零时刻的船头方向和船体左舷方向分别为x轴和y轴，z轴的方向满足右手定则；x轴、y轴和z轴分别是船只横滚、俯仰、偏航转动的中心轴，op是零时刻船只转动中心与雷达的连线在xoy平面的投影线。是op与x轴的夹角，ψ是op与船只转动中心-雷达连线的夹角，即掠射角。R₀是雷达在零时刻的位置矢量，r_k0是散射点k在零时刻的位置矢量。

为简化分析，将横滚、俯仰、偏航转动引起的多普勒频率单独进行分析，然后再将三种转动的作用进行综合。当船只仅进行横滚转动时，设t时刻散射点k的坐标从零时刻的位置r_k0＝(x_k0,y_k0,z_k0)^T转动到r_kt＝(x_kt,y_kt,z_kt)^T，r_kt的表达式如公式(1)所示：

其中，r_kt是散射点k在t时刻的位置矢量；θ_r(t)为船只在t时刻的横滚角；

由于横滚运动造成的散射体相对于雷达的径向速度v_kr(t)表达式如公式(2)所示：

v_kr(t)＝[ω_r(t)×r_kt]·i_R0 (2)

其中，ω_r(t)为横滚转动对应的角速度矢量，i_R0为R₀对应的单位矢量，设i_R0＝(m₁,m₂,m₃)^T；

m₁,m₂,m₃分别是i_R0矢量的x坐标、y坐标和z坐标。将式(1)代入式(2)，计算径向速度v_kr(t)引起的多普勒频率f_dkr(t)，如公式(3)所示：

其中，ω_r(t)为矢量ω_r(t)的模，λ为信号波长，(x_k0,y_k0,z_k0)为散射点k在零时刻时在船只坐标系中的三维坐标。

同理，可得俯仰和偏航转动各自引起的多普勒频率f_dkp(t)和f_dky(t)，对应的公式分别如公式(4)和(5)所示：

其中，ω_p(t)为俯仰转动角速度矢量ω_p(t)的模，ω_y(t)为偏航转动角速度矢量ω_y(t)的模，θ_p(t)和θ_y(t)分别为船只在t时刻的俯仰角和偏航角。

三种转动综合作用引起的多普勒频率f_dk(t)是f_dkr(t)、f_dkp(t)f_dky(t)三者之和。

当成像时间相对较短时，θ_r(t)、θ_p(t)和θ_y(t)的值通常都较小，多普勒频率f_dk(t)表达式如公式(6)所示：

根据船只水动力学的相关理论，ω_r(t)、ω_p(t)和ω_y(t)近似服从余弦函数的变化规律，由此可得：

其中，A_r、T_r和φ_r分别是横滚余弦角速度的幅度、周期和初相位，A_p、T_p和φ_p分别是俯仰余弦角速度的幅度、周期和初相位，A_y、T_y和φ_y分别是偏航余弦角速度的幅度、周期和初相位；

需指出的是，式(7)近似成立的前提是成像时间比较短，当成像时间较长时则误差会较大。

从式(7)可以看出，散射点k瞬时多普勒频率的变化可近似为三个余弦函数之和。

显然，要想从ISAR信号中直接提取多个散射点具有上述特点的瞬时多普勒频率是非常困难的。但另一方面，余弦函数为光滑函数，因此散射点瞬时多普勒频率的变化将具有良好的光滑性。根据这一特点，可以将完整的信号采集时间划分为若干个短时间段，在每个短时间段内散射点瞬时多普勒频率的变化可近似为线性变化。因此，可以将式(7)改写为：

其中，i为时间段的序号，N为信号采集时间内所划分的时间段数量，K_rki为散射点k第i个时间段对应的调频斜率，f_dcki为散射点k第i个时间段对应的多普勒中心频率，t_start为信号采集的起点时刻，T_seg为时间段长度，rect(t/T_seg)为[-T_seg/2,T_seg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数。

显然，式(8)的精度与时间段长度T_seg有关，T_seg越长，式(8)的精度越低。

下面，推导T_seg的上限公式。对式(7)进行泰勒展开，令其二次项小于系统所希望取得的多普勒分辨率ρ_f的一半，即有：

其中，ρ_f为系统所希望取得的多普勒分辨率，K为船只散射点的数目，T_obs为信号采集时间，经整理后，可得：

其中：

考虑最不利的情况，可得：

根据表1数据以及式(12)估计T_seg的典型数值。设ρ_f＝2Hz，ψ＝20°，λ＝0.031m。

驱逐舰的散射点坐标设置为：x_k0＝100m，y_k0＝20m，z_k0＝30m；船只转动角速度余弦函数中的相位设置为：

航空母舰的散射点坐标设置为：x_k0＝200m，y_k0＝50m，z_k0＝60m；船只转动角速度余弦函数中的相位设置为：

所计算得出的值，对于驱逐舰为T_seg≈0.5s，对于航空母舰为T_seg≈2s。

表1 5级海情下船只三维转动的典型参数

为直观地表现散射点的瞬时多普勒频率变化特性，开展了仿真工作。其它参数的设置和前面保持不变，散射点的数目由单个增加到多个。驱逐舰的散射点分布如图2(a)所示，航空母舰的散射点分布如图2(b)所示。另外，设平台高度为5000m，平台速度和船只速度均为零。这里将平台速度和船只速度均设为零，相当于观察的是平动补偿后散射点多普勒频率的变化。

图3(a)和图3(b)分别显示了驱逐舰和航空母舰两艘船只各随机选中的5个散射点在-5s至5s间的瞬时多普勒频率变化。这些曲线是根据散射点与雷达间距离的变化通过数值计算得出的，可以认为是代表了多普勒频率的“真值”。从图中可看出，在比较长的时间范围内，瞬时多普勒频率的变化明显偏离线性变化，但由于曲线比较光滑，当时间段的长度比较短时则瞬时多普勒频率的变化呈现近似的线性规律。

3、算法模型和处理流程

在ISAR系统中，方位时间-距离频率域的船只回波信号S(t,f_r)可用公式(13)表示为：

其中，t表示方位向时间，f_r表示距离向频率，rect(·)为矩形函数，T_obs为信号的数据采集时间，f_c为雷达载波频率，B_r为雷达信号带宽，C表示光速，R₀(t)表示从雷达的相位中心到船只转动中心的矢量，R₀(t)为R₀(t)的模，i_R0(t)为R₀(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ₀(r)为位于r处散射点的归一化RCS(Radar Cross Section，雷达后向散射截面积)，A(r,t)为考虑发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点RCS的时变性等因素对信号幅度和相位造成影响的复函数。

经过距离对齐、相位补偿等平动补偿处理后的信号可表示为：

其中，R′₀为成像中心时刻散射点与雷达间的距离。

对式(14)离散化后的信号进行距离向IFFT(Inverse Fast FourierTransformation，快速逆傅里叶变换)，可得：

其中，τ表示距离向时间，K为散射点的数目，A_k(t)为第k个散射点对应的A(r,t)，σ_k为第k个散射点的RCS，τ_k为第k个散射点相对于雷达的双程延迟时间，r_k为从船只转动中心到第k个散射点的矢量，R′_k为r_k在成像中心时刻在雷达视线上的投影长度，相位项包含着散射点k的多普勒信息。为简化起见，可将式(15)改写为：

当船只有效转动矢量的指向和大小均保持不变时，式(16)中的多普勒频率项f_dk(t)近似为常数，这种情况下直接使用RD算法进行成像即可。当在实际的海洋环境中，由于海浪等因素的作用，船只散射点的多普勒频率为时变的。如前所述，船只散射点的多普勒频率近似满足式(8)，在较短的时间段内可近似地视为线性调频信号，在不同时间段内的中心频率和调频斜率并不相同。

从信号s′(t,τ)中挑选出L个功率最强的距离单元，设l_j表示其中第j个距离单元在所有距离单元中的序号，表示该距离单元对应的一维信号，显然有：

其中，表示序号为l_j的距离单元中所包含的散射点的数目，表示该距离单元所对应的距离向时间。

将式(8)代入式(17)后进行简化表达，并考虑到A_k(t)的变化通常比较缓慢，可得到如下近似表达式：

式中，表示复函数A_k(t)中的相位，为式(17)中积分项所对应的初始相位。

设表示中第i个时间段的信号，则有：

从式(21)可看出，是一个多散射点线性调频信号的和。

根据信号在噪声环境下估计A_ki、θ_ki、K_rki、f_dcki的值。

离散调频傅里叶变换是实现上述估计的经典工具，但受栅栏效应的影响，一般采用改进的离散调频傅里叶变换。但由于时间段的长度T_seg比较短，当利用改进的离散调频傅里叶变换进行分析时分辨率较低，若存在散射点多普勒频率比较接近的情况则严重影响参数估计的准确性。考虑到船只图像中散射点的数目相比于整幅图像的像素数而言通常是比较稀疏的，因此本申请选用压缩感知的方法实现A_ki、θ_ki、K_rki、f_dcki的估计。另外，为进一步提高分辨率，在进行压缩感知计算之前首先利用经典的Burg算法进行数据外推，从而扩大数据量。设表示对进行数据外推后得到的信号，为的离散化形式，其中m＝0,1,···,M-1，M为利用经典的Burg算法进行外推后的数据点数，M的值取偶数，PRI为脉冲重复周期(Pulse Repetition Interval)。接下来针对建立稀疏表达模型，实现A_ki、θ_ki、K_rki、f_dcki的估计。

为减少压缩感知计算过程中的数据量，对进行随机采样测量，可得：

其中，为的矩阵形式，Φ为随机采样过程对应的测量矩阵，为对进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M。

测量矩阵Φ可选用随机高斯矩阵、随机伯努利矩阵等多种形式。由于为多个线性调频信号和的形式，因此采用如下的基矩阵Ψ：

其中，W_M＝e^-j2π/M，信号可表示为：

其中，θ为稀疏系数矢量。

由于船只散射点的数量在整副图像像素数中占比较小，θ中只有少量元素的值相对较大，其它元素的值都接近于0。将式(24)代入式(22)，可得：

其中，T为感知矩阵。考虑噪声的影响下，对系数矢量θ的求解可转化为如下优化问题：

其中，||·||₁表示l₁范数，||·||₂表示l₂范数，ε为与噪声电平相关的常数。

设K_sp表示稀疏度，则当感知矩阵T满足K_sp阶RIP(Restricted IsometryProperty，有限等距性质)条件，且时，通过求解上述优化问题，可保证以极大的概率恢复出系数矢量θ。

求解出系数矢量θ后，就可以得到A_ki、θ_ki、K_rki、f_dcki的值。

由于误差等因素的影响，对于同一个距离单元，在不同的时间段所估计出的散射点个数可能不相等，因此将序号为l_j的距离单元在第i个时间段所估计出的散射点数目记为将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点。等于矩阵中有效局部峰值点的个数。

K_rki和f_dcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式为：

n_kic＝n_ki-n_kirM (28)；

其中，n_ki为第k个散射点第i个时间段时对应的局部峰值系数在矢量θ中的序号(0≤n_ki≤M²-1)，floor(·)为取自变量整数部分的函数，n_kir为对应的调频率序号(0≤n_kir≤M-1)，n_kic为对应的中心频率序号(0≤n_kic≤M-1)，PRF为脉冲重复频率(Pulse RepetitionFrequency)。

提取出L个功率最强的距离单元在所有时间段内对应的K_rki和f_dcki后，可得到这些距离单元每个时间段内多个散射点的瞬时多普勒频率变化。接下来需要对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到各散射点在完整数据采集时间内的多普勒频率变化。设序号为l_j的距离单元在第i个时间段估计出个散射点，其中第l个散射点在第i个时间段结束时刻的瞬时多普勒频率为f_dcli+K_rli·T_seg/2。设该距离单元在第i+1个时间段估计出个散射点，其中第l′个散射点在第i+1个时间段开始时刻的瞬时多普勒频率为f_dcl′(i+1)-K_rl′(i+1)·T_seg/2，按如下原则对l和l′进行关联：

取1≤i≤N中的最小值为

完成了L个功率最强的距离单元各个时间段内散射点的关联后，可得到这些距离单元在完整数据采集时间内各散射点的瞬时多普勒频率变化。

接下来，利用L个功率最强的距离单元内各散射点瞬时多普勒频率的变化，可计算出最优成像时间范围。首先，以序号为l_j的距离单元内的第l个散射点为例，先计算针对该散射点而言各时刻的允许积累时间长度。以任意的t_exp时刻为例，计算以该时刻为基准的允许积累时间长度。从t_exp时刻向后进行逐时刻搜索，直到某个时刻的瞬时多普勒频率与t_exp时刻瞬时多普勒频率之差的绝对值超出期望的多普勒分辨率ρ_f，将搜索终止时刻与t_exp时刻之间的时间差记为

这里需要说明的是，当搜索过程中遇到完整采集时间内的终止时刻t_start+N·T_seg时，搜索过程自然终止，将搜索终止时刻记为t_start+N·T_seg。

利用上述过程，可以得到L个功率最强的距离单元各散射点在完整采集时间内各个时刻的允许积累时间长度。

显然，对于完整采集时间内的各个时刻，允许的积累时间长度为：

最终，取允许积累时间长度最长的时刻为选取时窗的起始时刻t_opt,start，最长的允许积累时间为选取时窗的时长Δt_opt，即有：

Δt_opt＝Δt(t_opt,start) (34)；

其中，Δt(t_opt,start)表示Δt函数在t_opt,start时刻的函数值。

图4给出了本发明方法的处理流程，具体如下：

步骤1：对ISAR回波进行距离压缩、平动补偿等预处理；

步骤2：挑选出L个功率最强的距离单元，并对这些距离单元按功率从高到低编上从1至L的编号；

步骤3：对这L个距离单元的信号在方位向上按T_seg为间隔划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；

步骤4：设j为L个功率最强的距离单元中当前单元的编号，令j的初值为1；

步骤5：并设l_j表示该距离单元在所有距离单元中的序号。之后，对序号为l_j的距离单元在第i个时间段内的信号利用经典的Burg算法进行数据外推后得到将表示成矩阵对该矩阵乘以测量矩阵后Φ得到对矩阵利用正交匹配追踪(OMP，orthogonal matched pursuit)等重建算法，估计出稀疏系数矢量得到系数矢量后，利用式(27)-式(30)计算得到K_rki和f_dcki的值，其中下标k的取值范围为从1至利用同样的处理流程，通过二重循环，计算出L个最强距离单元所有时间段对应的K_rki和f_dcki。然后，从第2个时间段开始，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联。接下来，利用式(32)，计算得出完整采集时间内各个时刻允许的积累时间长度Δt(t_exp)。最后，利用式(33)和式(34)，计算得出最终的最佳成像起始时刻t_opt,start、时窗长度Δt_opt。

4、仿真与实测数据处理结果

为验证所提出方法的正确性，开展了一系列的仿真和实测数据处理实验。首先，开展了单散射点的仿真实验，目的是验证估计出的瞬时多普勒频率曲线具有较高的精度。然后，开展了船只的仿真实验，目的包括：(1)验证算法在多散射点条件下的瞬时多普勒频率曲线估计精度，(2)展现所提出方法的成像结果具有较好的成像质量，(3)与ICBA方法的成像结果进行量化比较。最后，开展了实测数据的处理实验，目的是验证算法在实际环境条件下的正确性。

在单散射点的仿真实验中，采用了如表2所示的仿真参数。其中，设置了两种情况下的船只转动参数。第1种情况下，船只仅在偏航轴进行转动。第2种情况下，船只在横滚、俯仰、偏航轴均进行转动，其中横滚轴的转动幅度明显大于另外两个轴的转动。仿真中加入了高斯热噪声的影响。为保证仿真系统的实际等效噪声系数值与表2中设置的值一致，发射信号的峰值功率P_t由下式计算得出：

其中，k_B为波尔兹曼常数，T₀＝290K，F为接收机噪声系数，S_L为系统损耗，R_max为波束照射范围内与雷达的最远距离，PRI为脉冲重复周期，NEσ⁰为归一化等效噪声系数，ρ_az,e为期望的方位分辨率，ρ_rg为地距分辨率，G为天线增益，T_r为脉冲宽度，T_int为脉冲积累时间。代入表2中的相关参数，计算得出的峰值功率约为120w。令σ表示散射点的RCS，并将其设置为

σ＝(NEσ⁰·10⁶)·ρ_az，eρ_rg (36)；

图5给出了单散射点条件下所估计出的瞬时多普勒频率曲线。其中，子图(a)对应第1种情况，子图(b)对应第2种情况。为进行估计误差的评估，图中还同时给出了瞬时多普勒频率“真值”对应的曲线。这里的真值曲线是根据雷达与散射点间距离的变化利用数值计算而得出的。从图中可以看出，瞬时多普勒频率的估计结果与“真值”之间的吻合度是比较好的。经过计算，子图(a)中的最大误差约为0.5Hz，子图(b)中的最大误差约为1.8Hz。

表2 单散射点仿真实验中使用的主要参数

在船只仿真实验中，散射点的位置分布采用了如图2(a)所示的分布。在这些散射点中，将其中4个散射点的RCS设置为其它散射点的10倍。这4个强散射点分布于2个距离单元中，它们零时刻在船只坐标系中的坐标分别为：[-26,12,0]、[80,12,0]、[-26,-12,0]、[80,-12,0]。仿真实验中，设置L＝2，船只转动参数采用单散射点仿真实验第2种情况中设置的参数。其它参数的设置与单散射点仿真实验中的设置相同。图6给出了多散射点条件下所估计出的瞬时多普勒频率曲线。其中，子图(a)对应第1个强距离单元的2个强散射点，子图(b)对应第2个强距离单元的2个强散射点。从图中可以看出，瞬时多普勒频率的估计结果与“真值”之间的吻合度是比较好的。经过计算，子图(a)中的最大误差约为2.5Hz，子图(b)中的最大误差约为1Hz。

分别采用本申请算法与ICBA算法，对适宜的ISAR时窗进行了估计。本申请算法得出的时间段约为[-4.88s,-4.63s]，ICBA算法得出的时间段约为[3.92s,4.12s]。图7给出了上述两个时间段内4个强散射点瞬时多普勒频率真值的变化。其中，子图(a)对应本申请算法，子图(b)对应ICBA算法。从图7中可看出，本申请算法所选出的时间段长于ICBA算法所选出的时间段，并且本申请算法选出的时间段内4个强散射点多普勒频率的变化均小于所设置的期望多普勒分辨率(ρ_f＝2Hz)，而ICBA算法选出的时间段内有2个强散射点多普勒频率的变化超出了设置的ρ_f＝2Hz。图8给出了船只的成像结果。其中，子图(a)对应本申请算法，子图(b)对应ICBA算法。为评估两种算法的成像质量，对4个强散射点对应的多普勒分辨率、方位向峰值旁瓣比、方位向积分旁瓣比进行了计算并求取平均值。表3列出了4个强散射点对应的多普勒分辨率、方位向峰值旁瓣比、方位向积分旁瓣比的平均值。从表3可看出，采用本申请算法取得的多普勒分辨率更佳。

表3.仿真实验的主要技术指标

为验证本申请算法在实际环境条件下的有效性，针对2018年秋季在中国某海域获取的机载ISAR实测数据，开展了成像处理工作。表4给出了该次载机飞行实验过程中雷达平台的主要参数。

表4.实测数据获取过程中雷达平台的主要参数

对实测ISAR回波数据经距离压缩、平动补偿等预处理后，分别采用本申请算法和ICBA算法进行了加窗处理。之后，利用RD算法对加窗后的数据进行成像，得到两种情况下的成像结果图。另外，为进行对比，还给出了不进行加窗处理情况下的成像结果图。图9(a)、(b)、(c)分别为不加窗、采用本申请算法进行加窗处理、采用ICBA算法进行加窗处理三种情况下的成像结果。表5给出了这三种成像结果分别对应的图像对比度和熵。从图9和表5的计算结果可以看出：(1)当不进行加窗处理时，所得图像在方位向的聚焦效果较差，图像对比度不高，图像的熵值较大；(2)当采用本申请算法和ICBA算法加窗处理后得到的图像聚焦效果均较好。两幅图像的对比度值基本相当，而本申请算法所得结果的熵值略小一些。

表5.实测数据处理实验的主要技术指标

5、结论

基于多普勒频率信息的算法是目前选取ISAR船只成像过程中适宜时间窗的一类常用算法。但现有算法要么仅提取单个散射点的瞬时多普勒频率，要么提取的是整艘船只的平均多普勒频率。实际上，船只上主要强散射点的瞬时多普勒频率变化并不相同。针对上述不足，本申请提出一种基于多散射点瞬时多普勒频率估计的算法。通过理论分析和数值仿真，验证了船只转动引起的瞬时多普勒频率仅在较短的时间内近似服从线性调频信号的规律。推导得到上述短时间段的时间长度计算公式。设计了一种数据外推与压缩感知相结合的估计方法，用于估计短时线性调频信号的调频斜率和中心频率。该方法将数据外推处理的高分辨特性与压缩感知计算的低旁瓣特性相结合，从而保证了即使在信号长度较短、从而频率分辨率较低的情况下，依然具有较高的调频斜率和中心频率估计精度。通过上述方法得到各短时间段内散射点的瞬时多普勒频率变化后，对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到主要强散射点在完整数据采集时间内的瞬时多普勒频率变化。最终，综合多个强散射点的瞬时多普勒频率变化信息，选取出这些强散射点的瞬时多普勒频率变化都比较平稳的时间段，将该时间段作为最终选定的时间窗。分别利用仿真和实测数据，对本申请所提出方法的有效性进行了验证，并与ICBA算法的成像效果进行了对比，本申请算法取得了更好的成像质量。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.一种选取逆合成孔径雷达船只成像时窗的方法，其特征在于：包括如下步骤：

步骤1：对逆合成孔径雷达回波信号进行包括距离压缩、平动补偿在内的预处理；

将完整的信号采集时间划分为若干个短时间段，在每个短时间段内散射点瞬时多普勒频率的变化如公式(8)所示：

其中，t表示方位向时间，f_dk(t)表示第k个散射点的瞬时多普勒频率，i表示时间段的序号，N表示所有时间段的个数，K_rki表示第k个散射点在第i个时间段内的线性调频斜率，t_start表示信号采集时间的起始时刻，T_seg表示每个时间段的时间长度，f_dcki表示第k个散射点在第i个时间段内的中心频率，rect(t/T_seg)为[-T_seg/2,T_seg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数；

在逆合成孔径雷达系统中，方位时间-距离频率域的船只回波信号表示如公式(13)所示：

其中，f_r表示距离向频率，S(t,f_r)表示方位时间-距离频率域的船只回波信号，T_obs为信号的数据采集时间，f_c为雷达载波频率，B_r为雷达信号带宽，C表示光速，R₀(t)表示从雷达的相位中心到船只转动中心的矢量，R₀(t)为R₀(t)的模，i_R0(t)为R₀(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ₀(r)为位于r处散射点的归一化雷达后向散射截面积，A(r,t)为考虑发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点雷达后向散射截面积的时变性等因素对信号幅度和相位造成影响的复函数；

经过包括距离对齐、相位补偿在内的平动补偿处理后的信号S′(t,f_r)表示如公式(14)所示：

其中，R′₀为成像中心时刻散射点与雷达间的距离；

对式(14)离散化后的信号进行距离向快速逆傅里叶变换，得到方位时间-距离时间域的信号s′(t,τ)，如公式(15)所示：

其中，τ表示距离向时间，K为散射点的数目，A_k(t)为第k个散射点对应的A(r,t)，σ_k为第k个散射点的雷达后向散射截面积，τ_k为第k个散射点相对于雷达的双程延迟时间，r_k为从船只转动中心到第k个散射点的矢量，R_k′为r_k在成像中心时刻在雷达视线上的投影长度，相位项包含着散射点k的多普勒频率信息；

将公式(15)进行简化，得到公式(16)：

步骤2：从预处理后的ISAR回波信号中，挑选出L个功率最强的距离单元，并将这些距离单元按功率从大到小从按照1至L编上编号；

步骤3：以T_seg为时间段长度，在方位向上将步骤2中L个距离单元的信号划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；

步骤4：设j为L个功率最强的距离单元中当前单元的编号，并令j的初值为1；

步骤5：设l_j为第j个距离单元在所有距离单元中的序号，表示序号为l_j的距离单元对应的一维信号，表示中第i个时间段的信号，利用Burg算法，对进行数据外推，得到

其中，表示序号为l_j的距离单元中所包含的散射点的数目，表示序号为l_j的距离单元所对应的距离向时间；

将式(8)代入式(17)后进行简化，得到如下近似表达式：

式中，表示复函数A_k(t)中的相位，为式(17)中积分项所对应的初始相位；因此，有：

步骤6：将按脉冲重复间隔为采样周期进行离散化，并将离散化后的信号表示成矩阵形式，将该矩阵乘以随机采样过程对应的测量矩阵Φ得到如公式(22)所示：

其中，为的矩阵形式，m为方位向采样点的序号，M为利用Burg算法进行数据外推后对应的方位向采样点数，M的值取偶数，0≤m≤M-1，Φ为随机采样过程对应的测量矩阵，为对进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M；PRI为脉冲重复间隔；

步骤7：利用正交匹配追踪重建算法，根据公式(26)，估计出稀疏系数矢量

信号又能表示为：

其中，W_M＝e^-j2π/M，θ为稀疏系数矢量；

将式(24)代入式(22)，得到如下表达式：

其中，T为感知矩阵；

其中，||·||₁表示l₁范数，||·||₂表示l₂范数，ε为与噪声电平相关的常数；

步骤8：估计K_rki和f_dcki的值，其中下标k的取值范围为从1至

将序号为l_j的距离单元在第i个时间段所估计出的散射点数目记为将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点，等于矩阵中有效局部峰值点的个数；

K_rki和f_dcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式如式(27)-式(30)所示：

n_kic＝n_ki-n_kirM (28)；

其中，n_ki为第k个散射点第i个时间段对应的局部峰值系数在矢量θ中的序号，0≤n_ki≤M²-1，floor(·)为取自变量整数部分的函数，n_kir为对应的调频率序号，0≤n_kir≤M-1，n_kic为对应的多普勒中心频率序号，0≤n_kic≤M-1，PRF为脉冲重复频率；K_rki为散射点k第i个时间段对应的调频斜率，f_dcki为散射点k第i个时间段对应的多普勒中心频率；

步骤9：判断j是否等于L；

若：判断结果是j＝L；则执行步骤10；

或判断结果是j≠L；则令j＝j+1，然后执行步骤5；

步骤10：判断i是否等于N；

若：判断结果是i＝N；则执行步骤11；

或判断结果是i≠N；则令i＝i+1，然后执行步骤4；计算出L个最强距离单元所有时间段对应的K_rki和f_dcki；

步骤11：从第2个时间段开始，至第N个时间段结束，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联；

其中，f_dcli+K_rli·T_seg/2为第l个散射点在第i个时间段结束时刻的瞬时多普勒频率， 为序号为l_j的距离单元在第i个时间段估计出的散射点数；f_dcl′(i+1)-K_rl′(i+1)·T_seg/2为第l′个散射点在第i+1个时间段开始时刻的瞬时多普勒频率， 为序号为l_j的距离单元在第i+1个时间段估计出的散射点数；K_rl′(i+1)为散射点l′第i+1个时间段对应的调频斜率；K_rli为散射点l第i个时间段对应的调频斜率；取中的最小值为

步骤12：利用式(32)，计算完整采集时间内各个离散时刻t_exp允许的积累时间长度Δt(t_exp)；

其中，t_exp＝m·PRI,0≤m≤M-1，为序号为l_j的距离单元中第l个散射点对应的搜索终止时刻与t_exp时刻之间的时间差，这里的搜索终止时刻记为t_exp时刻之后的第1个具有如下特点的时刻：序号为l_j的距离单元中第l个散射点在该时刻的瞬时多普勒频率与该散射点在t_exp时刻瞬时多普勒频率之差的绝对值大于期望的多普勒分辨率ρ_f；当搜索过程中遇到完整采集时间内的终止时刻t_start+N·T_seg时，搜索过程自然终止，将搜索终止时刻记为t_start+N·T_seg；

步骤13：利用式(33)和式(34)，计算最终的最佳成像起始时刻t_opt,start和时窗长度Δt_opt；

Δt_opt＝Δt(t_opt,start) (34)；

其中，Δt(t_opt,start)表示Δt函数在t_opt,start时刻的函数值；

步骤14：结束。