CN105629233B - 一种基于isar图像的多散射点微动提取方法 - Google Patents

Abstract

本发明提供一种基于ISAR图像的多散射点微动提取方法，包括以下步骤：步骤S1，选取ISAR图像的同一距离单元内的多个强散射点；步骤S2，用加窗的方式提取每一强散射点的频域信号，并将其变换为时域信号；步骤S3，提取时域信号的相位，根据最小二乘法，用三阶多项式拟合转动所引起的相位变化曲线，利用所述相位变化曲线对时域信号相位进行转动补偿后得到微动相位变化曲线，进一步根据所述相位变化曲线获得转动补偿后的回波数据；步骤S4，对转动补偿后的回波数据，采用余弦相位模型，利用模拟退火算法，估计微动频率。该方法可以快速有效地提取出多散射点的微动频率，具有广泛的应用前景。

Description

一种基于ISAR图像的多散射点微动提取方法

技术领域

本发明涉及一种基于ISAR图像的多散射点微动提取方法，属于信号处理技术领域。

背景技术

微动可以定义为除了目标整体运动(平动)之外的某部件的机械运动。微动在频谱上的体现是微多普勒现象，即目标主体运动所产生的多普勒频率附近会存在旁瓣，称为微多普勒。微动成分中包含目标的细微特征，这些特征(旋转频率、旋转半径等)对目标识别有很大的应用价值，因此存在微动结构的雷达目标特征分析已经引起了研究人员的广泛关注。

逆合成孔径雷达(ISAR)是一种高分辨的成像雷达，通过对目标平动补偿和自聚焦的方法，可以获取非合作目标的距离-多普勒二维图像即ISAR图像。在雷达回波中，微动相关的回波很弱，尤其在微动幅度很小或者目标整体运动补偿效果不好的时候，如何提取散射点微动信息一直是研究难点。时频分析是一种研究微动的重要工具，对某距离单元内的信号作时频变换可以观察信号时频分布特性，但是当实际的微动频率很小时，从各散射点的时频分布特性也无法分辨出微动分量。而且，当同一个距离单元内，存在多个散射点的回波信号时，在时频图上各散射点的时频分布交织在一起更难以分辨。

因此，为了解决ISAR图像中同一个距离单元内存在多个散射点回波的微动分析问题，且微动频率很小(小于10Hz)幅度很小，同时为了使用合理的数学模型和参数估计方法提取微动频率，提出了基于ISAR图像的多散射点微动分析方法。

发明内容

本发明是为了克服现有技术的缺陷，提出了一种基于ISAR图像的多散射点微动提取方法，该方法能够处理ISAR图像中同一距离单元存在多个散射点回波且微动频率很小(小于10Hz)及幅度很小的情况。

实现本发明的技术方案如下：

一种基于ISAR图像的多散射点微动提取方法，包括以下步骤：

步骤S1，选取ISAR图像的同一距离单元内的多个强散射点；

步骤S2，用加窗的方式提取每一强散射点的频域信号，并将其变换为时域信号；

步骤S3，提取时域信号的相位，根据最小二乘法，用三阶多项式拟合转动所引起的相位变化曲线，利用所述相位变化曲线对时域信号相位进行转动补偿后得到微动相位变化曲线，进一步根据所述相位变化曲线获得转动补偿后的回波数据；

步骤S4，对转动补偿后的回波数据，采用余弦相位模型，利用模拟退火算法，估计微动频率。

进一步地，本发明所述步骤S2的具体过程为：

步骤S21，针对第k个强散射点，构造窗函数向量Window_k，窗函数的中心为f_k，主瓣包括区间[f_k-10,f_k+10]，k＝1,2…K，其中K为所提取的强散射点的总数，f_k为第k个强散射点的多普勒频率；

步骤S22，加窗提取第k个强散射点的频域信号即其中，符号.*表示两个向量的对应元素点乘；

步骤S23，将第k个强散射点的频域信号变换为时域信号S_k。

进一步地，本发明所述步骤S3的具体过程为：

步骤S31：提取第k个强散射点的时域信号S_k的相位P_k；

步骤S32：根据最小二乘法，用三阶多项式拟合转动所引起的相位，得到拟合曲线Q_k，

Q_k＝a₀+a₁t+a₂t²+a₃t³，t＝0,PRT,...,(N-1)PRT，PRT为ISAR图像方位向采样间隔；

步骤S33：转动相位补偿得到微动相位变化曲线P_Micro，即P_Micro＝P_k-Q_k，而微动对应的回波数据表示为

进一步地，本发明所述步骤S4的具体过程为：

步骤S41，设置初始温度t₁，温度终止阈值t_threshold，马尔科夫链长L，初始解初始为循环次数i＝1，代价函数C(Γ)为

PRT为ISAR图像方位向的采样间隔，N是方位向采样点数，S_Micro(n)表示S_Micro的第n个元素；

步骤S42，在初始解Γ⁰的邻域内选择新解Γ¹，若C(Γ¹)＞C(Γ⁰)，则将Γ¹的值赋予Γ⁰；若C(Γ¹)＜C(Γ⁰)，则在区间[0,1]选取一个随机数Random，若exp[-(C(Γ⁰)-C(Γ¹))/t_i]＞Random，则将Γ¹的值赋予Γ⁰，否则Γ⁰的值不变；

步骤S43，重复步骤S42的过程L次，并记录L次过程中使得目标函数最优的解；

步骤S43，令外循环的次数i加一，更新温度t_i＝0.95t_i-1，重复步骤S42和S43，直至温度t_i小于所设的阈值t_threshold；

步骤S44，从各温度所对应的最优解中选取最大值Γ_opt，由此得到微动频率估计值f_Micro＝Γ_opt(2)，即Γ_opt的第二个元素。

有益效果：

本发明提出一种基于ISAR图像的多散射点微动提取方法，以加窗的方式提取单个强散射点的回波，补偿掉转动所引起的相位分量后，得到微动相位变化曲线，转动补偿后的数据可视为微动对应的回波数据，用余弦模型近似表示，模拟退火算法估计出模型中的参数。该方法可以快速有效地提取出多散射点的微动频率，具有广泛的应用前景。

附图说明

图1是本发明的流程示意图；

图2是本发明的飞机ISAR图和选取的机头散射点示意图；

图3是本发明的飞机机头散射点的微动相位变化曲线；

图4是本发明的飞机ISAR图和选取的发动机散射点示意图；

图5是本发明的飞机发动机散射点的微动相位变化曲线。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

如图1所示，本发明一种基于ISAR图像的多散射点微动提取方法，包括以下步骤：

步骤S1，选取ISAR图像的某距离单元内的几个强散射点；

所述ISAR图像是目标的距离-多普勒二维图像，可表示为N行M列的矩阵，行表示方位向(即多普勒频域)，列表示距离向，距离单元可视作ISAR图像矩阵的某一列，记作Rangecell，后续的分析即是针对Rangecell。

所述步骤S1的具体过程为：

步骤S11，设定阈值Threshold；

步骤S12，将所选定的距离单元Rangecell的复数数据按幅值排序，设有K个值大于所设阈值Threshold，将Rangecell内幅度大于阈值Threshold的复数数据记为f₁,f₂,...,f_K，其对应K个强散射点。

步骤S2，用加窗的方式提取每一强散射点的频域信号，并将其变换为时域信号；

该步骤的具体过程为：

步骤S21，针对第k个强散射点，构造窗函数向量Window_k，窗函数的中心为f_k，主瓣包括区间[f_k-10,f_k+10]，k＝1,2…K；

步骤S22，加窗提取第k个散射点的频域信号即其中，符号.*表示两个向量的对应元素点乘；

步骤S23，将第k个散射点的频域信号变换到时域信号S_k，即IFFT为逆傅里叶变换。

步骤S3，提取时域信号S_k的相位，根据最小二乘法，用三阶多项式拟合转动所引起的相位变化曲线，利用所述相位变化曲线对时域信号进行转动补偿后得到微动相位变化曲线，进一步根据所述相位变化曲线获得转动补偿后的回波数据；

所述步骤S3的具体过程为：

步骤S31：提取第k个强散射点的时域信号S_k的相位P_k；

步骤S32：根据最小二乘法，用三阶多项式拟合转动所引起的相位，得到拟合曲线Q_k

Q_k＝a₀+a₁t+a₂t²+a₃t³，t＝0,PRT,...,(N-1)PRT，PRT为ISAR图像方位向采样间隔；

步骤S33：转动相位补偿得到微动相位变化曲线，即P_Micro＝P_k-Q_k，而微动对应的回波数据可以表示为

步骤S4，对转动补偿后的回波数据，采用余弦相位模型，利用模拟退火算法，估计微动频率；

所述步骤S4的具体过程为：

步骤S41，设置初始温度t₁，温度终止阈值t_threshold，马尔科夫链长L，初始解初始为循环次数i＝1，代价函数C(Γ)为

PRT为ISAR图像方位向的采样间隔，N是方位向采样点数，S_Micro(n)表示S_Micro的第n个元素；

步骤S42，在初始解Γ⁰的邻域内选择新解Γ¹，将所述Γ⁰和Γ¹分别带入代价函数C(Γ)中，若C(Γ¹)＞C(Γ⁰)，则接受Γ¹为当前解，即将Γ¹的值赋予Γ⁰；若C(Γ¹)＜C(Γ⁰)，则在区间[0,1]选取一个随机数Random，若exp[-(C(Γ⁰)-C(Γ¹))/t_i]＞Random，则依然接受Γ¹为当前解，即将Γ¹的值赋予Γ⁰，否则Γ⁰的值不变；

步骤S43，重复步骤S42的过程L次，并记录L次过程中使得目标函数最优(即最大)的解；

步骤S43，令外循环的次数i加一，更新温度t_i＝0.95t_i-1，重复步骤S42和S43，直至温度t_i小于所设的阈值t_threshold；

步骤S44，从各个温度所对应的最优解中选取最大值Γ_opt，由此得到微动频率估计值f_Micro＝Γ_opt(2)，即Γ_opt的第二个元素。

目标上的微动往往是规律性很强的振动、旋转等周期性运动，如直升飞机旋翼的旋转、引擎引起的表面振动等，这些周期性微动提供了目标运动的细节信息，因此估计微动频率对于目标检测、识别具有重要的意义。

本发明提供的方法既可以在个人计算机、工控机及服务器上以软件的形式安装并执行，也可将方法做成嵌入式芯片以硬件的形式来体现。

下面给出使用上述方法，对实测的飞机回波数据进行微动分析的实施例。

雷达参数：中心频率3GHz，频带宽度300MHz，脉冲宽度26us，脉冲重复频率245Hz，取128帧回波数据作为实验数据，成像目标为波音777。

图2是波音777的ISAR图像，即距离-多普勒二维图像，横向表示距离，纵向表示多普勒。选取图2中第103个距离单元分析圈中所表示的散射点的微动特性，以该散射点为中心进行加窗提取出该散射点对应的频域回波做逆傅里叶变换得到时域回波S_k。

图3是S_k的相位变化曲线P_k和拟合曲线Q_k的差值曲线，即微动相位变化曲线P_Micro，可以看出该曲线是类似余弦变化的。

对转动补偿后的数据，即微动对应的回波数据，用余弦模型估计的结果是6.5Hz，即飞机的机头散射点存在6.5Hz的微动。

图4是波音777的ISAR图像，即距离-多普勒二维图像，横向表示距离，纵向表示多普勒。选取图4中第125个距离单元分析圈中所表示的散射点的微动特性，以该散射点为中心进行加窗提取出该散射点对应的频域回波做逆傅里叶变换得到时域回波S_k。

图5是S_k的相位变化曲线P_k和拟合曲线Q_k的差值曲线，即微动相位变化曲线P_Micro，可以看出该曲线是类似余弦变化的。

对转动补偿后的数据，即微动对应的回波数据，用余弦模型估计的结果是4.3Hz，即飞机的发动机散射点存在4.3Hz的微动。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于ISAR图像的多散射点微动提取方法，其特征在于，包括以下步骤：

步骤S1，选取ISAR图像的同一距离单元内的多个强散射点；

步骤S2，用加窗的方式提取每一强散射点的频域信号，并将其变换为时域信号；

步骤S3，提取时域信号的相位，根据最小二乘法，用三阶多项式拟合转动所引起的相位变化曲线，利用所述相位变化曲线对时域信号相位进行转动补偿后得到微动相位变化曲线，进一步根据所述相位变化曲线获得转动补偿后的回波数据；

步骤S4，对转动补偿后的回波数据，采用余弦相位模型，利用模拟退火算法，估计微动频率。

2.根据权利要求1所述基于ISAR图像的多散射点微动提取方法，其特征在于，所述步骤S2的具体过程为：

步骤S21，针对第k个强散射点，构造窗函数向量Window_k，窗函数的中心为f_k，主瓣包括区间[f_k-10,f_k+10]，k＝1,2…K，其中K为所提取的强散射点的总数，f_k为第k个强散射点的多普勒频率；

步骤S22，加窗提取第k个强散射点的频域信号S_fk，即S_fk＝Rangecell.*Window_k，其中，符号.*表示两个向量的对应元素点乘；

步骤S23，将第k个强散射点的频域信号S_fk变换为时域信号S_k。

3.根据权利要求1所述基于ISAR图像的多散射点微动提取方法，其特征在于，所述步骤S3的具体过程为：

步骤S31：提取第k个强散射点的时域信号S_k的相位P_k；

步骤S32：根据最小二乘法，用三阶多项式拟合转动所引起的相位，得到拟合曲线Q_k，

Q_k＝a₀+a₁t+a₂t²+a₃t³，t＝0,PRT,...,(N-1)PRT，a₀、a₁、a₂、a₃为待定系数，PRT为ISAR图像方位向采样间隔；

步骤S33：转动相位补偿得到微动相位变化曲线P_Micro，即P_Micro＝P_k-Q_k，而微动对应的回波数据表示为

4.根据权利要求1所述基于ISAR图像的多散射点微动提取方法，其特征在于，所述步骤S4的具体过程为：

步骤S41，设置初始温度t₁，温度终止阈值t_threshold，马尔科夫链长L，初始解初始为循环次数i＝1，代价函数C(Γ)为

其中，B为余弦相位模型的幅度，θ为余弦相位模型的初始相位，PRT为ISAR图像方位向的采样间隔，N是方位向采样点数，S_Micro(n)表示回波数据S_Micro的第n个元素；

步骤S42，在初始解Γ⁰的邻域内选择新解Γ¹，若C(Γ¹)＞C(Γ⁰)，则将Γ¹的值赋予Γ⁰；若C(Γ¹)＜C(Γ⁰)，则在区间[0,1]选取一个随机数Random，若exp[-(C(Γ⁰)-C(Γ¹))/t_i]＞Random，则将Γ¹的值赋予Γ⁰，否则Γ⁰的值不变；

步骤S43，重复步骤S42的过程L次，并记录L次过程中使得目标函数最优的解；

步骤S43，令外循环的次数i加一，更新温度t_i＝0.95t_i-1，重复步骤S42和S43，直至温度t_i小于所设的阈值t_threshold；

步骤S44，从各温度所对应的最优解中选取最大值Γ_opt，由此得到微动频率估计值f_Micro＝Γ_opt(2)，即Γ_opt的第二个元素。