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

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

Info

Publication number
CN110133648B
CN110133648B CN201910395702.7A CN201910395702A CN110133648B CN 110133648 B CN110133648 B CN 110133648B CN 201910395702 A CN201910395702 A CN 201910395702A CN 110133648 B CN110133648 B CN 110133648B
Authority
CN
China
Prior art keywords
time
scattering point
distance
scattering
signal
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910395702.7A
Other languages
English (en)
Other versions
CN110133648A (zh
Inventor
周鹏
吕伟强
隋明妍
薛莲
戴永寿
孙伟峰
万勇
李立刚
曲晓俊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN201910395702.7A priority Critical patent/CN110133648B/zh
Publication of CN110133648A publication Critical patent/CN110133648A/zh
Application granted granted Critical
Publication of CN110133648B publication Critical patent/CN110133648B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/415—Identification of targets based on measurements of movement associated with the target
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects

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)所示：
Figure BDA0002056724830000031
其中，t表示方位向时间，fdk(t)表示第k个散射点的瞬时多普勒频率，i表示时间段的序号，N表示所有时间段的个数，Krki表示第k个散射点在第i个时间段内的线性调频斜率，tstart表示信号采集时间的起始时刻，Tseg表示每个时间段的时间长度，fdcki表示第k个散射点在第i个时间段内的中心频率，rect(t/Tseg)为[-Tseg/2,Tseg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数；
在逆合成孔径雷达系统中，方位时间-距离频率域的船只回波信号表示如公式(13)所示：
Figure BDA0002056724830000032
其中，fr表示距离向频率，S(t,fr)表示方位时间-距离频率域的船只回波信号，
Figure BDA0002056724830000033
Tobs为信号的数据采集时间，fc为雷达载波频率，Br为雷达信号带宽，C表示光速，R0(t)表示从雷达的相位中心到船只转动中心的矢量，R0(t)为R0(t)的模，iR0(t)为R0(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ0(r)为位于r处散射点的归一化RCS(Radar Cross Section，雷达后向散射截面积)，A(r,t)为考虑发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点RCS的时变性等因素对信号幅度和相位造成影响的复函数；
经过包括距离对齐、相位补偿在内的平动补偿处理后的信号S′(t,fr)表示如公式(14)所示：
Figure BDA0002056724830000034
其中，R′0为成像中心时刻散射点与雷达间的距离；
对式(14)离散化后的信号进行距离向IFFT(Inverse Fast FourierTransformation，快速逆傅里叶变换)，得到方位时间-距离时间域的信号s′(t,τ)，如公式(15)所示：
Figure BDA0002056724830000041
其中，τ表示距离向时间，K为散射点的数目，Ak(t)为第k个散射点对应的A(r,t)，σk为第k个散射点的RCS，τk为第k个散射点相对于雷达的双程延迟时间，rk为从船只转动中心到第k个散射点的矢量，R′k为rk在成像中心时刻在雷达视线上的投影长度，相位项
Figure BDA0002056724830000042
包含着散射点k的多普勒频率信息；
将公式(15)进行简化，得到公式(16)：
Figure BDA0002056724830000043
步骤2：从预处理后的ISAR回波信号中，挑选出L个功率最强的距离单元，并将这些距离单元按功率从大到小从按照1至L编上编号；
步骤3：以Tseg为时间段长度，在方位向上将步骤2中L个距离单元的信号划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；
步骤4：设j为L个功率最强的距离单元中当前单元的编号，并令j的初值为1；
步骤5：设lj为第j个距离单元在所有距离单元中的序号，
Figure BDA0002056724830000044
表示序号为lj的距离单元对应的一维信号，
Figure BDA0002056724830000045
表示
Figure BDA0002056724830000046
中第i个时间段的信号，利用经典的Burg算法，对
Figure BDA0002056724830000047
进行数据外推，得到
Figure BDA0002056724830000048
Figure BDA0002056724830000049
其中，
Figure BDA00020567248300000410
表示序号为lj的距离单元中所包含的散射点的数目，
Figure BDA00020567248300000411
表示序号为lj的距离单元所对应的距离向时间；
将式(8)代入式(17)后进行简化，得到如下近似表达式：
Figure BDA00020567248300000412
Figure BDA0002056724830000051
Figure BDA0002056724830000052
式中，
Figure BDA0002056724830000053
表示复函数Ak(t)中的相位，
Figure BDA0002056724830000054
为式(17)中积分项所对应的初始相位；因此，有：
Figure BDA0002056724830000055
步骤6：将
Figure BDA0002056724830000056
按PRI(Pulse Repetition Interval，脉冲重复间隔)为采样周期进行离散化，并将离散化后的信号表示成矩阵
Figure BDA0002056724830000057
形式，将该矩阵乘以随机采样过程对应的测量矩阵Φ得到
Figure BDA0002056724830000058
如公式(22)所示：
Figure BDA0002056724830000059
其中，
Figure BDA00020567248300000510
为
Figure BDA00020567248300000511
的矩阵形式，m为方位向采样点的序号，M为利用经典的Burg算法进行数据外推后对应的方位向采样点数，M的值取偶数，0≤m≤M-1，Φ为随机采样过程对应的测量矩阵，
Figure BDA00020567248300000512
为对
Figure BDA00020567248300000513
进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M；PRI为脉冲重复间隔；
步骤7：利用OMP(orthogonal matched pursuit，正交匹配追踪)重建算法，根据公式(26)，估计出稀疏系数矢量
Figure BDA00020567248300000514
信号
Figure BDA00020567248300000515
又能表示为：
Figure BDA00020567248300000516
其中，WM＝e-j2π/M，θ为稀疏系数矢量；
Figure BDA00020567248300000517
将式(24)代入式(22)，得到如下表达式：
Figure BDA0002056724830000061
其中，T为感知矩阵；
Figure BDA0002056724830000062
其中，||·||1表示l1范数，||·||2表示l2范数，ε为与噪声电平相关的常数；
步骤8：估计
Figure BDA0002056724830000063
Krki和fdcki的值，其中下标k的取值范围为从1至
Figure BDA0002056724830000064
将序号为lj的距离单元在第i个时间段所估计出的散射点数目记为
Figure BDA0002056724830000065
将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点，
Figure BDA0002056724830000066
等于矩阵中有效局部峰值点的个数；
Krki和fdcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式如式(27)-式(30)所示：
Figure BDA0002056724830000067
nkic＝nki-nkirM(28)；
Figure BDA0002056724830000068
Figure BDA0002056724830000069
其中，nki为第k个散射点第i个时间段对应的局部峰值系数在矢量θ中的序号，0≤nki≤M2-1，floor(·)为取自变量整数部分的函数，nkir为对应的调频率序号，0≤nkir≤M-1，nkic为对应的多普勒中心频率序号，0≤nkic≤M-1，PRF(Pulse Repetition Frequency)为脉冲重复频率；Krki为散射点k第i个时间段对应的调频斜率，fdcki为散射点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个最强距离单元所有时间段对应的
Figure BDA0002056724830000071
Krki和fdcki；
步骤11：从第2个时间段开始，至第N个时间段结束，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联；
Figure BDA0002056724830000072
Figure BDA0002056724830000073
其中，fdcli+Krli·Tseg/2为第l个散射点在第i个时间段结束时刻的瞬时多普勒频率，
Figure BDA0002056724830000074
1≤i≤N，
Figure BDA0002056724830000075
为序号为lj的距离单元在第i个时间段估计出的散射点数；
Figure BDA0002056724830000076
为第l′个散射点在第i+1个时间段开始时刻的瞬时多普勒频率，
Figure BDA0002056724830000077
Figure BDA0002056724830000078
为序号为lj的距离单元在第i+1个时间段估计出的散射点数；Krl′(i+1)为散射点l′第i+1个时间段对应的调频斜率；Krli为散射点l第i个时间段对应的调频斜率；取
Figure BDA0002056724830000079
1≤i≤N中的最小值为
Figure BDA00020567248300000710
步骤12：利用式(32)，计算完整采集时间内各个离散时刻texp允许的积累时间长度Δt(texp)；
Figure BDA00020567248300000711
其中，texp＝m·PRI,0≤m≤M-1，
Figure BDA00020567248300000712
为序号为lj的距离单元中第l个散射点对应的搜索终止时刻与texp时刻之间的时间差，这里的搜索终止时刻记为texp时刻之后的第1个具有如下特点的时刻：序号为lj的距离单元中第l个散射点在该时刻的瞬时多普勒频率与该散射点在texp时刻瞬时多普勒频率之差的绝对值大于期望的多普勒分辨率ρf；这里需要说明的是，当搜索过程中遇到完整采集时间内的终止时刻tstart+N·Tseg时，搜索过程自然终止，将搜索终止时刻记为tstart+N·Tseg；
步骤13：利用式(33)和式(34)，计算最终的最佳成像起始时刻topt,start和时窗长度Δtopt；
Figure BDA0002056724830000081
Δtopt＝Δt(topt,start) (34)；
其中，Δt(topt,start)表示Δt函数在topt,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平面的投影线。
Figure BDA0002056724830000091
是op与x轴的夹角，ψ是op与船只转动中心-雷达连线的夹角，即掠射角。R0是雷达在零时刻的位置矢量，rk0是散射点k在零时刻的位置矢量。
为简化分析，将横滚、俯仰、偏航转动引起的多普勒频率单独进行分析，然后再将三种转动的作用进行综合。当船只仅进行横滚转动时，设t时刻散射点k的坐标从零时刻的位置rk0＝(xk0,yk0,zk0)T转动到rkt＝(xkt,ykt,zkt)T，rkt的表达式如公式(1)所示：
Figure BDA0002056724830000101
其中，rkt是散射点k在t时刻的位置矢量；θr(t)为船只在t时刻的横滚角；
由于横滚运动造成的散射体相对于雷达的径向速度vkr(t)表达式如公式(2)所示：
vkr(t)＝[ωr(t)×rkt]·iR0 (2)
其中，ωr(t)为横滚转动对应的角速度矢量，iR0为R0对应的单位矢量，设iR0＝(m1,m2,m3)T；
m1,m2,m3分别是iR0矢量的x坐标、y坐标和z坐标。将式(1)代入式(2)，计算径向速度vkr(t)引起的多普勒频率fdkr(t)，如公式(3)所示：
Figure BDA0002056724830000102
其中，ωr(t)为矢量ωr(t)的模，λ为信号波长，(xk0,yk0,zk0)为散射点k在零时刻时在船只坐标系中的三维坐标。
同理，可得俯仰和偏航转动各自引起的多普勒频率fdkp(t)和fdky(t)，对应的公式分别如公式(4)和(5)所示：
Figure BDA0002056724830000103
Figure BDA0002056724830000104
其中，ωp(t)为俯仰转动角速度矢量ωp(t)的模，ωy(t)为偏航转动角速度矢量ωy(t)的模，θp(t)和θy(t)分别为船只在t时刻的俯仰角和偏航角。
三种转动综合作用引起的多普勒频率fdk(t)是fdkr(t)、fdkp(t)fdky(t)三者之和。
当成像时间相对较短时，θr(t)、θp(t)和θy(t)的值通常都较小，多普勒频率fdk(t)表达式如公式(6)所示：
Figure BDA0002056724830000105
根据船只水动力学的相关理论，ωr(t)、ωp(t)和ωy(t)近似服从余弦函数的变化规律，由此可得：
Figure BDA0002056724830000111
其中，Ar、Tr和φr分别是横滚余弦角速度的幅度、周期和初相位，Ap、Tp和φp分别是俯仰余弦角速度的幅度、周期和初相位，Ay、Ty和φy分别是偏航余弦角速度的幅度、周期和初相位；
需指出的是，式(7)近似成立的前提是成像时间比较短，当成像时间较长时则误差会较大。
从式(7)可以看出，散射点k瞬时多普勒频率的变化可近似为三个余弦函数之和。
显然，要想从ISAR信号中直接提取多个散射点具有上述特点的瞬时多普勒频率是非常困难的。但另一方面，余弦函数为光滑函数，因此散射点瞬时多普勒频率的变化将具有良好的光滑性。根据这一特点，可以将完整的信号采集时间划分为若干个短时间段，在每个短时间段内散射点瞬时多普勒频率的变化可近似为线性变化。因此，可以将式(7)改写为：
Figure BDA0002056724830000112
其中，i为时间段的序号，N为信号采集时间内所划分的时间段数量，Krki为散射点k第i个时间段对应的调频斜率，fdcki为散射点k第i个时间段对应的多普勒中心频率，tstart为信号采集的起点时刻，Tseg为时间段长度，rect(t/Tseg)为[-Tseg/2,Tseg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数。
显然，式(8)的精度与时间段长度Tseg有关，Tseg越长，式(8)的精度越低。
下面，推导Tseg的上限公式。对式(7)进行泰勒展开，令其二次项小于系统所希望取得的多普勒分辨率ρf的一半，即有：
Figure BDA0002056724830000113
其中，ρf为系统所希望取得的多普勒分辨率，K为船只散射点的数目，Tobs为信号采集时间，经整理后，可得：
Figure BDA0002056724830000121
其中：
Figure BDA0002056724830000122
考虑最不利的情况，可得：
Figure BDA0002056724830000123
根据表1数据以及式(12)估计Tseg的典型数值。设ρf＝2Hz，
Figure BDA0002056724830000124
ψ＝20°，λ＝0.031m。
驱逐舰的散射点坐标设置为：xk0＝100m，yk0＝20m，zk0＝30m；船只转动角速度余弦函数中的相位设置为：
Figure BDA0002056724830000125
航空母舰的散射点坐标设置为：xk0＝200m，yk0＝50m，zk0＝60m；船只转动角速度余弦函数中的相位设置为：
Figure BDA0002056724830000126
所计算得出的值，对于驱逐舰为Tseg≈0.5s，对于航空母舰为Tseg≈2s。
表1 5级海情下船只三维转动的典型参数
Figure BDA0002056724830000127
为直观地表现散射点的瞬时多普勒频率变化特性，开展了仿真工作。其它参数的设置和前面保持不变，散射点的数目由单个增加到多个。驱逐舰的散射点分布如图2(a)所示，航空母舰的散射点分布如图2(b)所示。另外，设平台高度为5000m，平台速度和船只速度均为零。这里将平台速度和船只速度均设为零，相当于观察的是平动补偿后散射点多普勒频率的变化。
图3(a)和图3(b)分别显示了驱逐舰和航空母舰两艘船只各随机选中的5个散射点在-5s至5s间的瞬时多普勒频率变化。这些曲线是根据散射点与雷达间距离的变化通过数值计算得出的，可以认为是代表了多普勒频率的“真值”。从图中可看出，在比较长的时间范围内，瞬时多普勒频率的变化明显偏离线性变化，但由于曲线比较光滑，当时间段的长度比较短时则瞬时多普勒频率的变化呈现近似的线性规律。
3、算法模型和处理流程
在ISAR系统中，方位时间-距离频率域的船只回波信号S(t,fr)可用公式(13)表示为：
Figure BDA0002056724830000131
其中，t表示方位向时间，fr表示距离向频率，
Figure BDA0002056724830000132
rect(·)为矩形函数，Tobs为信号的数据采集时间，fc为雷达载波频率，Br为雷达信号带宽，C表示光速，R0(t)表示从雷达的相位中心到船只转动中心的矢量，R0(t)为R0(t)的模，iR0(t)为R0(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ0(r)为位于r处散射点的归一化RCS(Radar Cross Section，雷达后向散射截面积)，A(r,t)为考虑发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点RCS的时变性等因素对信号幅度和相位造成影响的复函数。
经过距离对齐、相位补偿等平动补偿处理后的信号可表示为：
Figure BDA0002056724830000133
其中，R′0为成像中心时刻散射点与雷达间的距离。
对式(14)离散化后的信号进行距离向IFFT(Inverse Fast FourierTransformation，快速逆傅里叶变换)，可得：
Figure BDA0002056724830000134
其中，τ表示距离向时间，K为散射点的数目，Ak(t)为第k个散射点对应的A(r,t)，σk为第k个散射点的RCS，τk为第k个散射点相对于雷达的双程延迟时间，rk为从船只转动中心到第k个散射点的矢量，R′k为rk在成像中心时刻在雷达视线上的投影长度，相位项
Figure BDA0002056724830000135
包含着散射点k的多普勒信息。为简化起见，可将式(15)改写为：
Figure BDA0002056724830000136
当船只有效转动矢量的指向和大小均保持不变时，式(16)中的多普勒频率项fdk(t)近似为常数，这种情况下直接使用RD算法进行成像即可。当在实际的海洋环境中，由于海浪等因素的作用，船只散射点的多普勒频率为时变的。如前所述，船只散射点的多普勒频率近似满足式(8)，在较短的时间段内可近似地视为线性调频信号，在不同时间段内的中心频率和调频斜率并不相同。
从信号s′(t,τ)中挑选出L个功率最强的距离单元，设lj表示其中第j个距离单元在所有距离单元中的序号，
Figure BDA0002056724830000141
表示该距离单元对应的一维信号，显然有：
Figure BDA0002056724830000142
其中，
Figure BDA0002056724830000143
表示序号为lj的距离单元中所包含的散射点的数目，
Figure BDA0002056724830000144
表示该距离单元所对应的距离向时间。
将式(8)代入式(17)后进行简化表达，并考虑到Ak(t)的变化通常比较缓慢，可得到如下近似表达式：
Figure BDA0002056724830000145
Figure BDA0002056724830000146
Figure BDA0002056724830000147
式中，
Figure BDA0002056724830000148
表示复函数Ak(t)中的相位，
Figure BDA0002056724830000149
为式(17)中积分项所对应的初始相位。
设
Figure BDA00020567248300001410
表示
Figure BDA00020567248300001411
中第i个时间段的信号，则有：
Figure BDA00020567248300001412
从式(21)可看出，
Figure BDA00020567248300001413
是一个多散射点线性调频信号的和。
根据信号
Figure BDA00020567248300001414
在噪声环境下估计
Figure BDA00020567248300001415
Aki、θki、Krki、fdcki的值。
离散调频傅里叶变换是实现上述估计的经典工具，但受栅栏效应的影响，一般采用改进的离散调频傅里叶变换。但由于时间段的长度Tseg比较短，当利用改进的离散调频傅里叶变换进行分析时分辨率较低，若存在散射点多普勒频率比较接近的情况则严重影响参数估计的准确性。考虑到船只图像中散射点的数目相比于整幅图像的像素数而言通常是比较稀疏的，因此本申请选用压缩感知的方法实现
Figure BDA0002056724830000151
Aki、θki、Krki、fdcki的估计。另外，为进一步提高分辨率，在进行压缩感知计算之前首先利用经典的Burg算法进行数据外推，从而扩大数据量。设
Figure BDA0002056724830000152
表示对
Figure BDA0002056724830000153
进行数据外推后得到的信号，
Figure BDA0002056724830000154
为
Figure BDA0002056724830000155
的离散化形式，其中m＝0,1,···,M-1，M为利用经典的Burg算法进行外推后的数据点数，M的值取偶数，PRI为脉冲重复周期(Pulse Repetition Interval)。接下来针对
Figure BDA0002056724830000156
建立稀疏表达模型，实现
Figure BDA0002056724830000157
Aki、θki、Krki、fdcki的估计。
为减少压缩感知计算过程中的数据量，对
Figure BDA0002056724830000158
进行随机采样测量，可得：
Figure BDA0002056724830000159
其中，
Figure BDA00020567248300001510
为
Figure BDA00020567248300001511
的矩阵形式，Φ为随机采样过程对应的测量矩阵，
Figure BDA00020567248300001512
为对
Figure BDA00020567248300001513
进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M。
测量矩阵Φ可选用随机高斯矩阵、随机伯努利矩阵等多种形式。由于
Figure BDA00020567248300001514
为多个线性调频信号和的形式，因此采用如下的基矩阵Ψ：
Figure BDA00020567248300001515
其中，WM＝e-j2π/M，信号
Figure BDA00020567248300001516
可表示为：
Figure BDA00020567248300001517
其中，θ为稀疏系数矢量。
由于船只散射点的数量在整副图像像素数中占比较小，θ中只有少量元素的值相对较大，其它元素的值都接近于0。将式(24)代入式(22)，可得：
Figure BDA00020567248300001518
其中，T为感知矩阵。考虑噪声的影响下，对系数矢量θ的求解可转化为如下优化问题：
Figure BDA0002056724830000161
其中，||·||1表示l1范数，||·||2表示l2范数，ε为与噪声电平相关的常数。
设Ksp表示稀疏度，则当感知矩阵T满足Ksp阶RIP(Restricted IsometryProperty，有限等距性质)条件，且
Figure BDA0002056724830000162
时，通过求解上述优化问题，可保证以极大的概率恢复出系数矢量θ。
求解出系数矢量θ后，就可以得到
Figure BDA0002056724830000163
Aki、θki、Krki、fdcki的值。
由于误差等因素的影响，对于同一个距离单元，在不同的时间段所估计出的散射点个数可能不相等，因此将序号为lj的距离单元在第i个时间段所估计出的散射点数目记为
Figure BDA0002056724830000164
将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点。
Figure BDA0002056724830000165
等于矩阵中有效局部峰值点的个数。
Krki和fdcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式为：
Figure BDA0002056724830000166
nkic＝nki-nkirM (28)；
Figure BDA0002056724830000167
Figure BDA0002056724830000168
其中，nki为第k个散射点第i个时间段时对应的局部峰值系数在矢量θ中的序号(0≤nki≤M2-1)，floor(·)为取自变量整数部分的函数，nkir为对应的调频率序号(0≤nkir≤M-1)，nkic为对应的中心频率序号(0≤nkic≤M-1)，PRF为脉冲重复频率(Pulse RepetitionFrequency)。
提取出L个功率最强的距离单元在所有时间段内对应的
Figure BDA0002056724830000169
Krki和fdcki后，可得到这些距离单元每个时间段内多个散射点的瞬时多普勒频率变化。接下来需要对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到各散射点在完整数据采集时间内的多普勒频率变化。设序号为lj的距离单元在第i个时间段估计出
Figure BDA0002056724830000171
个散射点，其中第l个
Figure BDA0002056724830000172
散射点在第i个时间段结束时刻的瞬时多普勒频率为fdcli+Krli·Tseg/2。设该距离单元在第i+1个时间段估计出
Figure BDA0002056724830000173
个散射点，其中第l′个
Figure BDA0002056724830000174
散射点在第i+1个时间段开始时刻的瞬时多普勒频率为fdcl′(i+1)-Krl′(i+1)·Tseg/2，按如下原则对l和l′进行关联：
Figure BDA0002056724830000175
Figure BDA0002056724830000176
取
Figure BDA0002056724830000177
1≤i≤N中的最小值为
Figure BDA0002056724830000178
完成了L个功率最强的距离单元各个时间段内散射点的关联后，可得到这些距离单元在完整数据采集时间内各散射点的瞬时多普勒频率变化。
接下来，利用L个功率最强的距离单元内各散射点瞬时多普勒频率的变化，可计算出最优成像时间范围。首先，以序号为lj的距离单元内的第l个散射点为例，先计算针对该散射点而言各时刻的允许积累时间长度。以任意的texp时刻为例，计算以该时刻为基准的允许积累时间长度。从texp时刻向后进行逐时刻搜索，直到某个时刻的瞬时多普勒频率与texp时刻瞬时多普勒频率之差的绝对值超出期望的多普勒分辨率ρf，将搜索终止时刻与texp时刻之间的时间差记为
Figure BDA0002056724830000179
这里需要说明的是，当搜索过程中遇到完整采集时间内的终止时刻tstart+N·Tseg时，搜索过程自然终止，将搜索终止时刻记为tstart+N·Tseg。
利用上述过程，可以得到L个功率最强的距离单元各散射点在完整采集时间内各个时刻的允许积累时间长度。
显然，对于完整采集时间内的各个时刻，允许的积累时间长度为：
Figure BDA00020567248300001710
最终，取允许积累时间长度最长的时刻为选取时窗的起始时刻topt,start，最长的允许积累时间为选取时窗的时长Δtopt，即有：
Figure BDA0002056724830000181
Δtopt＝Δt(topt,start) (34)；
其中，Δt(topt,start)表示Δt函数在topt,start时刻的函数值。
图4给出了本发明方法的处理流程，具体如下：
步骤1：对ISAR回波进行距离压缩、平动补偿等预处理；
步骤2：挑选出L个功率最强的距离单元，并对这些距离单元按功率从高到低编上从1至L的编号；
步骤3：对这L个距离单元的信号在方位向上按Tseg为间隔划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；
步骤4：设j为L个功率最强的距离单元中当前单元的编号，令j的初值为1；
步骤5：并设lj表示该距离单元在所有距离单元中的序号。之后，对序号为lj的距离单元在第i个时间段内的信号
Figure BDA0002056724830000182
利用经典的Burg算法进行数据外推后得到
Figure BDA0002056724830000183
将
Figure BDA0002056724830000184
表示成矩阵
Figure BDA0002056724830000185
对该矩阵乘以测量矩阵后Φ得到
Figure BDA0002056724830000186
对矩阵
Figure BDA0002056724830000187
利用正交匹配追踪(OMP，orthogonal matched pursuit)等重建算法，估计出稀疏系数矢量
Figure BDA0002056724830000188
得到系数矢量
Figure BDA0002056724830000189
后，利用式(27)-式(30)计算得到
Figure BDA00020567248300001810
Krki和fdcki的值，其中下标k的取值范围为从1至
Figure BDA00020567248300001811
利用同样的处理流程，通过二重循环，计算出L个最强距离单元所有时间段对应的
Figure BDA00020567248300001812
Krki和fdcki。然后，从第2个时间段开始，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联。接下来，利用式(32)，计算得出完整采集时间内各个时刻允许的积累时间长度Δt(texp)。最后，利用式(33)和式(34)，计算得出最终的最佳成像起始时刻topt,start、时窗长度Δtopt。
4、仿真与实测数据处理结果
为验证所提出方法的正确性，开展了一系列的仿真和实测数据处理实验。首先，开展了单散射点的仿真实验，目的是验证估计出的瞬时多普勒频率曲线具有较高的精度。然后，开展了船只的仿真实验，目的包括：(1)验证算法在多散射点条件下的瞬时多普勒频率曲线估计精度，(2)展现所提出方法的成像结果具有较好的成像质量，(3)与ICBA方法的成像结果进行量化比较。最后，开展了实测数据的处理实验，目的是验证算法在实际环境条件下的正确性。
在单散射点的仿真实验中，采用了如表2所示的仿真参数。其中，设置了两种情况下的船只转动参数。第1种情况下，船只仅在偏航轴进行转动。第2种情况下，船只在横滚、俯仰、偏航轴均进行转动，其中横滚轴的转动幅度明显大于另外两个轴的转动。仿真中加入了高斯热噪声的影响。为保证仿真系统的实际等效噪声系数值与表2中设置的值一致，发射信号的峰值功率Pt由下式计算得出：
Figure BDA0002056724830000191
其中，kB为波尔兹曼常数，T0＝290K，F为接收机噪声系数，SL为系统损耗，Rmax为波束照射范围内与雷达的最远距离，PRI为脉冲重复周期，NEσ0为归一化等效噪声系数，ρaz,e为期望的方位分辨率，ρrg为地距分辨率，G为天线增益，Tr为脉冲宽度，Tint为脉冲积累时间。代入表2中的相关参数，计算得出的峰值功率约为120w。令σ表示散射点的RCS，并将其设置为
σ＝(NEσ0·106)·ρaz，eρrg (36)；
图5给出了单散射点条件下所估计出的瞬时多普勒频率曲线。其中，子图(a)对应第1种情况，子图(b)对应第2种情况。为进行估计误差的评估，图中还同时给出了瞬时多普勒频率“真值”对应的曲线。这里的真值曲线是根据雷达与散射点间距离的变化利用数值计算而得出的。从图中可以看出，瞬时多普勒频率的估计结果与“真值”之间的吻合度是比较好的。经过计算，子图(a)中的最大误差约为0.5Hz，子图(b)中的最大误差约为1.8Hz。
表2 单散射点仿真实验中使用的主要参数
Figure BDA0002056724830000192
在船只仿真实验中，散射点的位置分布采用了如图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.仿真实验的主要技术指标
Figure BDA0002056724830000201
为验证本申请算法在实际环境条件下的有效性，针对2018年秋季在中国某海域获取的机载ISAR实测数据，开展了成像处理工作。表4给出了该次载机飞行实验过程中雷达平台的主要参数。
表4.实测数据获取过程中雷达平台的主要参数
Figure BDA0002056724830000202
对实测ISAR回波数据经距离压缩、平动补偿等预处理后，分别采用本申请算法和ICBA算法进行了加窗处理。之后，利用RD算法对加窗后的数据进行成像，得到两种情况下的成像结果图。另外，为进行对比，还给出了不进行加窗处理情况下的成像结果图。图9(a)、(b)、(c)分别为不加窗、采用本申请算法进行加窗处理、采用ICBA算法进行加窗处理三种情况下的成像结果。表5给出了这三种成像结果分别对应的图像对比度和熵。从图9和表5的计算结果可以看出：(1)当不进行加窗处理时，所得图像在方位向的聚焦效果较差，图像对比度不高，图像的熵值较大；(2)当采用本申请算法和ICBA算法加窗处理后得到的图像聚焦效果均较好。两幅图像的对比度值基本相当，而本申请算法所得结果的熵值略小一些。
表5.实测数据处理实验的主要技术指标
Figure BDA0002056724830000211
5、结论
基于多普勒频率信息的算法是目前选取ISAR船只成像过程中适宜时间窗的一类常用算法。但现有算法要么仅提取单个散射点的瞬时多普勒频率，要么提取的是整艘船只的平均多普勒频率。实际上，船只上主要强散射点的瞬时多普勒频率变化并不相同。针对上述不足，本申请提出一种基于多散射点瞬时多普勒频率估计的算法。通过理论分析和数值仿真，验证了船只转动引起的瞬时多普勒频率仅在较短的时间内近似服从线性调频信号的规律。推导得到上述短时间段的时间长度计算公式。设计了一种数据外推与压缩感知相结合的估计方法，用于估计短时线性调频信号的调频斜率和中心频率。该方法将数据外推处理的高分辨特性与压缩感知计算的低旁瓣特性相结合，从而保证了即使在信号长度较短、从而频率分辨率较低的情况下，依然具有较高的调频斜率和中心频率估计精度。通过上述方法得到各短时间段内散射点的瞬时多普勒频率变化后，对相邻短时间段的瞬时多普勒频率变化进行关联，从而得到主要强散射点在完整数据采集时间内的瞬时多普勒频率变化。最终，综合多个强散射点的瞬时多普勒频率变化信息，选取出这些强散射点的瞬时多普勒频率变化都比较平稳的时间段，将该时间段作为最终选定的时间窗。分别利用仿真和实测数据，对本申请所提出方法的有效性进行了验证，并与ICBA算法的成像效果进行了对比，本申请算法取得了更好的成像质量。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.一种选取逆合成孔径雷达船只成像时窗的方法，其特征在于：包括如下步骤：
步骤1：对逆合成孔径雷达回波信号进行包括距离压缩、平动补偿在内的预处理；
将完整的信号采集时间划分为若干个短时间段，在每个短时间段内散射点瞬时多普勒频率的变化如公式(8)所示：
Figure FDA0003943351680000011
其中，t表示方位向时间，fdk(t)表示第k个散射点的瞬时多普勒频率，i表示时间段的序号，N表示所有时间段的个数，Krki表示第k个散射点在第i个时间段内的线性调频斜率，tstart表示信号采集时间的起始时刻，Tseg表示每个时间段的时间长度，fdcki表示第k个散射点在第i个时间段内的中心频率，rect(t/Tseg)为[-Tseg/2,Tseg/2]时间内函数值等于1、其余时间函数值等于0的矩形函数；
在逆合成孔径雷达系统中，方位时间-距离频率域的船只回波信号表示如公式(13)所示：
Figure FDA0003943351680000012
其中，fr表示距离向频率，S(t,fr)表示方位时间-距离频率域的船只回波信号，
Figure FDA0003943351680000013
Tobs为信号的数据采集时间，fc为雷达载波频率，Br为雷达信号带宽，C表示光速，R0(t)表示从雷达的相位中心到船只转动中心的矢量，R 0 (t)为R0(t)的模，iR0(t)为R0(t)的单位矢量，r为从船只转动中心到船上某散射点的矢量，σ0(r)为位于r处散射点的归一化雷达后向散射截面积，A(r,t)为考虑包括发射信号功率、雷达与散射点之间的距离、天线方向图效应、散射点雷达后向散射截面积的时变性在内的因素对信号幅度和相位造成影响的复函数；
经过包括距离对齐、相位补偿在内的平动补偿处理后的信号S′(t,fr)表示如公式(14)所示：
Figure FDA0003943351680000014
其中，R′0为成像中心时刻散射点与雷达间的距离；
对式(14)离散化后的信号进行距离向快速逆傅里叶变换，得到方位时间-距离时间域的信号s′(t,τ)，如公式(15)所示：
Figure FDA0003943351680000021
其中，τ表示距离向时间，K为散射点的数目，Ak(t)为第k个散射点对应的A(r,t)，σk为第k个散射点的雷达后向散射截面积，τk为第k个散射点相对于雷达的双程延迟时间，rk为从船只转动中心到第k个散射点的矢量，R′k为rk在成像中心时刻在雷达视线上的投影长度，相位项
Figure FDA0003943351680000022
包含着散射点k的多普勒频率信息；
将公式(15)进行简化，得到公式(16)：
Figure FDA0003943351680000023
步骤2：从预处理后的ISAR回波信号中，挑选出L个功率最强的距离单元，并将这些距离单元按功率从大到小从按照1至L编上编号；
步骤3：以Tseg为时间段长度，在方位向上将步骤2中L个距离单元的信号划分为N个时间段，设i为当前时间段的序号，并令i的初值为1；
步骤4：设j为L个功率最强的距离单元中当前单元的编号，并令j的初值为1；
步骤5：设lj为第j个距离单元在所有距离单元中的序号，
Figure FDA0003943351680000024
表示序号为lj的距离单元对应的一维信号，
Figure FDA0003943351680000025
表示
Figure FDA0003943351680000026
中第i个时间段的信号，利用Burg算法，对
Figure FDA0003943351680000027
进行数据外推，得到
Figure FDA0003943351680000028
Figure FDA0003943351680000029
其中，
Figure FDA00039433516800000210
表示序号为lj的距离单元中所包含的散射点的数目，
Figure FDA00039433516800000211
表示序号为lj的距离单元所对应的距离向时间；
将式(8)代入式(17)后进行简化，得到如下近似表达式：
Figure FDA0003943351680000031
Figure FDA0003943351680000032
Figure FDA0003943351680000033
式中，
Figure FDA0003943351680000034
表示复函数Ak(t)中的相位，
Figure FDA0003943351680000035
为式(17)中积分项所对应的初始相位；因此，有：
Figure FDA0003943351680000036
步骤6：将
Figure FDA0003943351680000037
按脉冲重复间隔为采样周期进行离散化，并将离散化后的信号表示成矩阵
Figure FDA0003943351680000038
形式，将该矩阵乘以随机采样过程对应的测量矩阵Φ得到
Figure FDA0003943351680000039
如公式(22)所示：
Figure FDA00039433516800000310
其中，
Figure FDA00039433516800000311
为
Figure FDA00039433516800000312
的矩阵形式，m为方位向采样点的序号，M为利用Burg算法进行数据外推后对应的方位向采样点数，M的值取偶数，0≤m≤M-1，Φ为随机采样过程对应的测量矩阵，
Figure FDA00039433516800000313
为对
Figure FDA00039433516800000314
进行随机采样测量后的矩阵，M′为随机采样后的点数且有M′＜＜M；PRI为脉冲重复间隔；
步骤7：利用正交匹配追踪重建算法，根据公式(26)，估计出稀疏系数矢量
Figure FDA00039433516800000315
信号
Figure FDA00039433516800000316
又能表示为：
Figure FDA00039433516800000317
其中，WM＝e-j2π/M，θ为稀疏系数矢量；
Figure FDA0003943351680000041
将式(24)代入式(22)，得到如下表达式：
Figure FDA0003943351680000042
其中，T为感知矩阵；
Figure FDA0003943351680000043
其中，||·||1表示l1范数，||·||2表示l2范数，ε为与噪声电平相关的常数；
步骤8：估计
Figure FDA0003943351680000044
Krki和fdcki的值，其中下标k的取值范围为从1至
Figure FDA0003943351680000045
将序号为lj的距离单元在第i个时间段所估计出的散射点数目记为
Figure FDA0003943351680000046
将L个最强的距离单元每个时间段求出的系数矢量转化为二维矩阵，并寻找矩阵中的局部峰值点，
Figure FDA0003943351680000047
等于矩阵中有效局部峰值点的个数；
Krki和fdcki的值取决于有效局部峰值点对应的系数在系数矢量中的位置，具体计算公式如式(27)-式(30)所示：
Figure FDA0003943351680000048
nkic＝nki-nkirM (28)；
Figure FDA0003943351680000049
Figure FDA00039433516800000410
其中，nki为第k个散射点第i个时间段对应的局部峰值系数在矢量θ中的序号，0≤nki≤M2-1，floor(·)为取自变量整数部分的函数，nkir为对应的调频率序号，0≤nkir≤M-1，nkic为对应的多普勒中心频率序号，0≤nkic≤M-1，PRF为脉冲重复频率；Krki为散射点k第i个时间段对应的调频斜率，fdcki为散射点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个最强距离单元所有时间段对应的
Figure FDA0003943351680000051
Krki和fdcki；
步骤11：从第2个时间段开始，至第N个时间段结束，利用式(31)，逐时间段地完成当前时间段与前一时间段L个最强距离单元内散射点的关联；
Figure FDA0003943351680000052
其中，fdcli+Krli·Tseg/2为第l个散射点在第i个时间段结束时刻的瞬时多普勒频率，
Figure FDA0003943351680000053
Figure FDA0003943351680000054
为序号为lj的距离单元在第i个时间段估计出的散射点数；fdcl′(i+1)-Krl′(i+1)·Tseg/2为第l′个散射点在第i+1个时间段开始时刻的瞬时多普勒频率，
Figure FDA0003943351680000055
Figure FDA0003943351680000056
为序号为lj的距离单元在第i+1个时间段估计出的散射点数；Krl′(i+1)为散射点l′第i+1个时间段对应的调频斜率；Krli为散射点l第i个时间段对应的调频斜率；取
Figure FDA0003943351680000057
中的最小值为
Figure FDA0003943351680000058
步骤12：利用式(32)，计算完整采集时间内各个离散时刻texp允许的积累时间长度Δt(texp)；
Figure FDA0003943351680000059
其中，texp＝m·PRI,0≤m≤M-1，
Figure FDA00039433516800000510
为序号为lj的距离单元中第l个散射点对应的搜索终止时刻与texp时刻之间的时间差，这里的搜索终止时刻记为texp时刻之后的第1个具有如下特点的时刻：序号为lj的距离单元中第l个散射点在该时刻的瞬时多普勒频率与该散射点在texp时刻瞬时多普勒频率之差的绝对值大于期望的多普勒分辨率ρf；当搜索过程中遇到完整采集时间内的终止时刻tstart+N·Tseg时，搜索过程自然终止，将搜索终止时刻记为tstart+N·Tseg；
步骤13：利用式(33)和式(34)，计算最终的最佳成像起始时刻topt,start和时窗长度Δtopt；
Figure FDA0003943351680000061
Δtopt＝Δt(topt,start) (34)；
其中，Δt(topt,start)表示Δt函数在topt,start时刻的函数值；
步骤14：结束。
CN201910395702.7A 2019-05-13 2019-05-13 一种选取逆合成孔径雷达船只成像时窗的方法 Active CN110133648B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910395702.7A CN110133648B (zh) 2019-05-13 2019-05-13 一种选取逆合成孔径雷达船只成像时窗的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910395702.7A CN110133648B (zh) 2019-05-13 2019-05-13 一种选取逆合成孔径雷达船只成像时窗的方法

Publications (2)

Publication Number Publication Date
CN110133648A CN110133648A (zh) 2019-08-16
CN110133648B true CN110133648B (zh) 2023-01-17

Family

ID=67573644

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910395702.7A Active CN110133648B (zh) 2019-05-13 2019-05-13 一种选取逆合成孔径雷达船只成像时窗的方法

Country Status (1)

Country Link
CN (1) CN110133648B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110456351B (zh) * 2019-08-29 2021-06-15 哈尔滨工业大学 基于时变幅值lfm信号参数估计的机动目标isar成像方法
CN111190151B (zh) * 2020-01-14 2021-10-08 中国石油大学(华东) 扫描模式下多模态小卫星sar的系统参数设计及发射功率优化方法
CN114942419B (zh) * 2022-07-26 2022-10-28 中国石油大学(华东) 长积累时间下的船只散射点三自由度微动特征提取方法
CN115407337B (zh) * 2022-11-01 2023-02-03 中国石油大学(华东) 一种基于时间窗二次选择的船只目标三维成像方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102879783A (zh) * 2012-10-12 2013-01-16 西安电子科技大学 基于稀疏探测频率信号的isar成像方法
JPWO2015008554A1 (ja) * 2013-07-19 2017-03-02 国立大学法人東北大学 合成開口処理を伴うセンサ、そのセンサの処理方法、および、プログラム
CN106772375A (zh) * 2016-12-27 2017-05-31 哈尔滨工业大学 基于参数估计的压缩感知成像方法
CN106918811A (zh) * 2017-04-05 2017-07-04 中国石油大学(华东) 一种逆合成孔径雷达船只成像时窗选取方法
CN107132535A (zh) * 2017-04-07 2017-09-05 西安电子科技大学 基于变分贝叶斯学习算法的isar稀疏频带成像方法
CN107390213A (zh) * 2017-07-14 2017-11-24 中南大学 一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法
CN109116352A (zh) * 2018-07-20 2019-01-01 中国石油大学(华东) 一种圆扫isar模式船只超分辨率成像方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9864054B2 (en) * 2014-03-10 2018-01-09 Mitsubishi Electric Research Laboratories, Inc. System and method for 3D SAR imaging using compressive sensing with multi-platform, multi-baseline and multi-PRF data

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102879783A (zh) * 2012-10-12 2013-01-16 西安电子科技大学 基于稀疏探测频率信号的isar成像方法
JPWO2015008554A1 (ja) * 2013-07-19 2017-03-02 国立大学法人東北大学 合成開口処理を伴うセンサ、そのセンサの処理方法、および、プログラム
CN106772375A (zh) * 2016-12-27 2017-05-31 哈尔滨工业大学 基于参数估计的压缩感知成像方法
CN106918811A (zh) * 2017-04-05 2017-07-04 中国石油大学(华东) 一种逆合成孔径雷达船只成像时窗选取方法
CN107132535A (zh) * 2017-04-07 2017-09-05 西安电子科技大学 基于变分贝叶斯学习算法的isar稀疏频带成像方法
CN107390213A (zh) * 2017-07-14 2017-11-24 中南大学 一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法
CN109116352A (zh) * 2018-07-20 2019-01-01 中国石油大学(华东) 一种圆扫isar模式船只超分辨率成像方法

Non-Patent Citations (13)

* Cited by examiner, † Cited by third party
Title
A METHOD OF AUTOMATED EXTRACTION OF GROUND CONTROL POINT UNDER UNDULATE AREAS;Sanhai Ren;《2009 IET International Radar Conference》;20091231;1-4 *
Circular Scan ISAR Mode Super-Resolution Imaging of Ships Based on a Combination of Data Extrapolation and Compressed Sensing;Peng Zhou;《IEEE SENSORS JOURNAL》;20190430;6883-6894 *
ISAR Imaging of Maneuvering Target Based on the Quadratic Frequency Modulated Signal Model With Time-Varying Amplitude;Yong Wang;《IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING》;20170331;1012-1024 *
ISAR舰船目标成像时间段选取;王冉等;《哈尔滨工业大学学报》;20111231(第007期);57-60 *
一种高分辨的稀疏孔径ISAR成像方法;李军等;《西安电子科技大学学报》;20100620(第03期);63-68、75 *
利用压缩感知实现随机变频雷达散射中心估计;杨新锋等;《红外与激光工程》;20160525(第05期);296-302 *
压缩感知机动目标ISAR成像新方法;陈春利等;《红外与激光工程》;20130825(第08期);324-329 *
含旋转部件目标稀疏孔径ISAR成像方法;徐艺萌等;《空军工程大学学报(自然科学版)》;20130825(第04期);61-65 *
基于匹配傅里叶变换的舰船目标成像算法;张云等;《大连海事大学学报》;20090215(第01期);50-55 *
基于局部窗口K分布的快速舰船检测算法;张颢;《计算机应用》;20160310;859-863 *
基于瞬时谱估计的ISAR距离瞬时多普勒成像算法;卢光跃等;《西安电子科技大学学报》;19981020(第05期);54-58 *
基于进动的旋转对称弹头雷达成像方法;刘进等;《信号处理》;20090925(第09期);3-7 *
海上船只目标多角度成像技术;周鹏等;《中国海洋大学学报(自然科学版)》;20170215(第02期);67-74 *

Also Published As

Publication number Publication date
CN110133648A (zh) 2019-08-16

Similar Documents

Publication Publication Date Title
CN110133648B (zh) 一种选取逆合成孔径雷达船只成像时窗的方法
Martorella Novel approach for ISAR image cross-range scaling
Martorella et al. 3D interferometric ISAR imaging of noncooperative targets
CN109116311B (zh) 基于知识辅助稀疏迭代协方差估计的杂波抑制方法
CN104035095B (zh) 基于空时最优处理器的低空风切变风速估计方法
CN107656255B (zh) 基于多径回波的超宽带雷达动目标二维定位方法
CN103869311B (zh) 实波束扫描雷达超分辨成像方法
US8797206B2 (en) Method and apparatus for simultaneous multi-mode processing performing target detection and tracking using along track interferometry (ATI) and space-time adaptive processing (STAP)
CN104898119B (zh) 一种基于相关函数的动目标参数估计方法
CN102707269B (zh) 一种机载雷达距离走动校正方法
CN105974410A (zh) 机载雷达的多舰船目标sar和isar混合成像方法
Li et al. Scaled Radon-Wigner transform imaging and scaling of maneuvering target
CN108469608A (zh) 一种运动平台雷达多普勒质心精确估计方法
CN104280566A (zh) 基于空时幅相估计的低空风切变风速估计方法
CN105738887A (zh) 基于多普勒通道划分的机载雷达杂波功率谱的优化方法
Zhou et al. Time window selection algorithm for ISAR ship imaging based on instantaneous Doppler frequency estimations of multiple scatterers
CN108107427A (zh) 基于超分辨技术的机载/弹载阵列雷达前视成像方法
CN108761417B (zh) 基于知识辅助最大似然的机载雷达杂波抑制方法
Kusk et al. SAR focusing of P-band ice sounding data using back-projection
CN104914421A (zh) 基于和差波束的低空风切变风速估计方法
CN101846741A (zh) 一种逆合成孔径雷达成像数据段选择方法
CN112305541A (zh) 一种基于抽样序列长度约束条件下的sar成像方法
Gao et al. A new cross-range scaling algorithm based on FrFT
CN114942419A (zh) 长积累时间下的船只散射点三自由度微动特征提取方法
CN103996214B (zh) Bp-sar图像重建误差分析方法及bp-sar图像重建方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant