CN104330799B - 一种基于粒子群滤波优化的isar成像方法 - Google Patents

Abstract

本发明属于ISAR成像技术领域，特别涉及一种基于粒子群滤波优化的ISAR成像方法，本发明将ISAR目标的平动分量进行多项式建模，通过ISAR图像熵最小化，利用PSO算法对多项式系数进行全局优化求解，最后补偿目标平动分量，得到高分辨率的二维成像结果，同时该技术能够实现目标平动多项式阶数的自适应估计。本发明能够有效避免局部最优问题，且求解精度高，对噪声适应性强。

Description

一种基于粒子群滤波优化的ISAR成像方法

技术领域

本发明属于ISAR(Inverse Synthetic Aperture Radar)成像技术领域，特别涉及一种基于粒子群滤波优化的ISAR成像方法。本发明涉及一种基于粒子群滤波优化(Particle Swarm Optimization，PSO)的低信噪比下逆合成孔径雷达(Inverse SyntheticAperture Radar，ISAR)成像方法，本发明将ISAR目标的平动分量进行多项式建模，通过ISAR图像熵最小化，利用PSO算法对多项式系数进行全局优化求解，最后补偿目标平动分量，得到高分辨率的二维成像结果，同时该技术能够实现目标平动多项式阶数的自适应估计。本发明能够有效避免局部最优问题，且求解精度高，对噪声适应性强。

背景技术

逆合成孔径雷达由于其全天时、全天候对观测目标成像监视的能力，在军事和民用领域都获得了广泛的应用。ISAR成像通过发射大带宽信号从而在距离维上得到高分辨率，方位维的高分辨率则取决于发射信号波长以及目标与雷达的相对转角。由于ISAR目标一般是非合作的，在相干积累时间内，目标的运动通常可分解为平动和转动，转动是ISAR方位成像的基础，而平动分量会引起距离压缩后的包络偏移，并影响回波之间的相干性，是成像时需要估计和补偿掉的。由于受到远观测距离、系统和环境噪声的影响，ISAR回波信号的信噪比往往不高。而在低信噪比下用传统的相邻相关法进行包络对齐无法达到所需效果，因此需要研究低信噪比下的ISAR平动补偿和成像方法。

杨磊等人在文献“应用联合自聚焦实现低信噪比ISAR成像平动补偿”(西安电子科技大学学报(自然科学版)，2012年)中，将目标的平动分量进行多项式建模，利用二维ISAR图像的熵作为目标函数，结合阻尼牛顿法进行最优化求解，在低信噪比条件下完成对平动分量的精确估计，最终得到目标聚焦良好的图像。但是上述方法存在两点不足：上述方法将目标平动用多项式进行建模，但是并没有给出一种有效的确定目标多项式阶数的准则和技术；由于目标函数与多项式系数之间的关系曲线不是凸函数，且存在大量的局部最优点，因此，上述求解方法容易陷入局部最优，从而不能估计出真实的目标平动，最终也就得不到聚焦度最好的ISAR图像。

发明内容

本发明的目的在于提出一种基于粒子群滤波优化的ISAR成像方法，本发明首先将目标平动用多项式进行建模，然后利用粒子群滤波优化(PSO)方法对以二维图像熵为目标函数的全局优化问题进行求解，同时实现对多项式阶数的确定和多项式系数的估计，其次根据估计出的多项式阶数和系数构造目标平动，补偿平动分量，最终得到聚焦良好的二维高分辨ISAR图像。

为实现上述技术目的，本发明采用如下技术方案予以实现。

一种基于粒子群滤波优化的ISAR成像方法包括以下步骤：

步骤1，利用逆合成孔径雷达发射线性调频信号S₀，利用逆合成孔径雷达接收原始回波数据，将原始回波数据沿距离向进行快速傅里叶变换，得到距离频域的数据矩阵S₁，对距离频域的数据矩阵S₁进行匹配滤波，得到匹配滤波之后的数据矩阵S_M；构建长度为M的时间列向量t，t＝[-M/2,-M/2+1,…,M/2-1]^T，上标T表示矩阵或向量的转置，M为距离频域的数据矩阵S₁的列数；设定目标平动分量多项式阶数P，P＝1,2,3…，当P＝1时，跳至步骤2；

步骤2，构建大小为M×P的时间矩阵D_P，时间矩阵D_P的第g列为列向量t^g，列向量t^g表示对时间列向量t的每个元素取g次方得出的新的向量，g取1至P；构建目标平动分量多项式系数向量b_P，当P＝1时，目标平动分量多项式系数向量b_P为0；当P>1时，目标平动分量多项式系数向量b_P是长度为P的列向量，目标平动分量多项式系数向量b_P中的元素全为0；令u＝1,2,3…，当u＝1时，跳至步骤3；

步骤3，设定粒子群的粒子总数为K，令k＝1,2,…K，设定粒子群中第k个粒子的初始位置向量粒子群中第k个粒子的初始速度向量粒子群中第k个粒子的初始最优位置向量以及粒子群第k个粒子所对应的初始图像熵值根据第u-1次迭代后粒子群中第k个粒子的速度向量第u-1次迭代后粒子群中第k个粒子的位置向量以及第u-1次迭代后粒子群中第k个粒子的最优位置列向量得到第u次迭代后粒子群中第k个粒子的速度向量得出第u次迭代后粒子群中第k个粒子的位置向量 $z_k^{P,u}=z_k^{P,u-1}+{\mathrm{\υ}}_k^{P,u};$

步骤4，得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵的图像熵值

步骤5，得出第u次迭代后粒子群所对应的图像最小熵值E^P,u,best、以及将第u次迭代后粒子群的最优位置向量z^P,u,best；

步骤6，用u_max表示设定的迭代次数u的经验迭代门限，如果u>u_max或得出目标平动分量多项式阶数为P时的图像最小熵E_P、以及目标平动分量多项式阶数为P时的目标平动分量多项式系数向量b_P，E_P＝E^P,u,best，b_P＝z^P,u,best，|·|₁表示1范数，ξ为设定的粒子速度平均经验值门限；否则，令u的值自增1，返回至步骤3；

步骤7，令目标平动分量多项式阶数P的值自增1，重复执行步骤2至步骤6；然后，判断E_P和E_P-1的大小关系，如果E_P≤E_P-1，则返回至步骤2；如果E_P>E_P-1，则得到最终图像最小熵、以及最终图像最小熵对应的目标平动多项式系数向量，所述最终图像最最小熵是目标平动分量多项式阶数为P-1时的图像最小熵E_P-1，所述最终图像最小熵对应的目标平动多项式系数向量是：目标平动分量多项式阶数为P-1时的目标平动分量多项式系数向量b_P-1；

步骤8，将时间矩阵D_P-1和最终图像最小熵对应的目标平动多项式系数向量的转置相乘得到目标平动分量R；用目标平动分量R对匹配滤波之后的数据矩阵S_M进行平动补偿得到目标平动分量补偿后数据矩阵A，对目标平动分量补偿后数据矩阵A做基于最终图像最小熵的自聚焦，并对自聚焦形成的图像在方位向做快速傅里叶变换，得到最终ISAR图像。

本发明的有益效果为：1)本发明利用粒子群滤波优化方法，在全局优化问题求解中能够有效避免局部收敛，具有优良的全局搜索性能，且具有模型简单、便于实现的特点，本发明对以二维ISAR图像熵为目标函数的优化问题通过粒子群滤波优化方法进行求解，克服了现有方法容易陷入局部最优解而造成的得不到聚焦良好的图像的问题。2)本发明根据目标平动的估计阶数与实际运动的阶数最接近时熵最小的原则，提出一种对多项式阶数和多项式系数的联合估计算法，克服了传统方法需要事先给定多项式阶数造成的目标平动估计不准或冗余的问题。

附图说明

图1为本发明的一种基于粒子群滤波优化的ISAR成像方法的流程图；

图2a为仿真实验中信噪比为0dB时的仿真数据距离包络示意图；

图2b为仿真实验中信噪比为0dB时仿真数据的理想ISAR成像结果示意图；

图2c为仿真实验中信噪比为0dB时对仿真数据采用本发明进行ISAR成像得出的ISAR成像结果示意图；

图2d为仿真实验中信噪比为-5dB时的仿真数据距离包络示意图；

图2e为仿真实验中信噪比为-5dB时仿真数据的理想ISAR成像结果示意图；

图2f为仿真实验中信噪比为-5dB时对仿真数据采用本发明进行ISAR成像得出的ISAR成像结果示意图；

图3a为仿真实验中信噪比为0dB时采用本发明得出的最终ISAR图像熵值和目标平动分量多项式阶数的关系示意图；

图3b为仿真实验中信噪比为0dB且目标平动分量多项式阶数为3时得出的各阶多项式系数和迭代次数u的关系示意图；

图3c为仿真实验中信噪比为-5dB时采用本发明得出的最终ISAR图像熵值和目标平动分量多项式阶数的关系示意图；

图3d为仿真实验中信噪比为-5dB且目标平动分量多项式阶数为3时得出的各阶多项式系数和迭代次数u的关系示意图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的一种基于粒子群滤波优化的ISAR成像方法的流程图。该基于粒子群滤波优化的ISAR成像方法的流程图包括以下步骤：

步骤1，利用逆合成孔径雷达发射线性调频信号S₀，利用逆合成孔径雷达接收原始回波数据，将原始回波数据沿距离向进行快速傅里叶变换，得到距离频域的数据矩阵S₁，对距离频域的数据矩阵S₁进行匹配滤波，得到匹配滤波之后的数据矩阵S_M；构建长度为M的时间列向量t，t＝[-M/2,-M/2+1,…,M/2-1]^T，上标T表示矩阵或向量的转置，M为距离频域的数据矩阵S₁的列数；设定目标平动分量多项式阶数P，P＝1,2,3…，当P＝1时，跳至步骤2。

其具体子步骤为：

(1.1)利用逆合成孔径雷达发射线性调频信号S₀，利用逆合成孔径雷达接收原始回波数据；对逆合成孔径雷达发射的线性调频信号S₀进行距离向快速傅里叶变换，得到线性调频信号S₀的频谱矩阵S₂。

(1.2)对原始回波数据在距离向进行快速傅里叶变换，得到距离频域的数据矩阵S₁，矩阵S₁的列数表示为M，行数表示为N。

(1.3)对距离频域的数据矩阵S₁进行匹配滤波，得到匹配滤波之后的数据矩阵S_M。对距离频域的数据矩阵S₁进行匹配滤波的过程可表示为下式：

S_M＝S₁·conj(S₂)

其中，·表示矩阵的点乘(两个矩阵对应元素相乘)，conj(S₂)表示将矩阵S₂的各元素取共轭得出的矩阵。

(1.4)构建长度为M的时间列向量t，t＝[-M/2,-M/2+1,…,M/2-1]^T，上标T表示矩阵或向量的转置，M为距离频域的数据矩阵S₁的列数；设定目标平动分量多项式阶数P，P＝1,2,3…，当P＝1时，跳至步骤2。

步骤2，构建目标平动分量多项式系数向量b_P，当P＝1时，目标平动分量多项式系数向量b_P为实数0(此时，b_P是长度为1的向量)；当P>1时，目标平动分量多项式系数向量b_P是长度为P的列向量，目标平动分量多项式系数向量b_P中的元素全为0；构建大小为M×P的时间矩阵D_P，D_P＝[t¹,t²,…,t^P]，其中，时间矩阵D_P的第g列为列向量t^g，列向量t^g表示对时间列向量t的每个元素取g次方得出的新的向量，g取1至P。

令u＝1,2,3…，当u＝1时，跳至步骤3。

步骤3，设定粒子群的粒子总数为K，令k＝1,2,…K，设定粒子群中第k个粒子的初始位置向量粒子群中第k个粒子的初始速度向量粒子群中第k个粒子的初始最优位置向量以及粒子群第k个粒子所对应的初始图像熵值

根据第u-1次迭代后粒子群中第k个粒子的速度向量第u-1次迭代后粒子群中第k个粒子的位置向量以及第u-1次迭代后粒子群中第k个粒子的最优位置列向量得到第u次迭代后粒子群中第k个粒子的速度向量得出第u次迭代后粒子群中第k个粒子的位置向量 $z_k^{P,u}=z_k^{P,u-1}+{\mathrm{\υ}}_k^{P,u}.$

其具体子步骤为：

(3.1)设定粒子群的粒子总数为K，令k＝1,2,…K，设定粒子群中第k个粒子的初始位置向量

$z_k^{P,0}=b_p+4\×[\mathrm{rand}\left(\mathrm{0,1}\right)-0.5]$

其中，k表示粒子编号，rand(0,1)表示产生一个0到1之间均匀分布的随机数；本发明实施例中，粒子群的粒子总数K越多，粒子搜索范围越大，越容易达到全局最优，但算法运行的时间也就越长，对大多数问题，粒子总数一般取20～40个，实验表明，对于大多数问题，30个粒子就够用了。因此，经验性地将K取为30。本发明实施例中，向量与数的加法描述为：将向量中每个元素与加上该数，得到新的向量。乘法符号“×”表示矩阵相乘，[rand(0,1)-0.5]表示标量为1×1的矩阵。

设定粒子群中第k个粒子的初始速度向量

${\mathrm{\υ}}_k^{P,0}={\mathrm{\υ}}_{\max }\×\mathrm{rand}\left(\mathrm{0,1}\right)$

其中，υ_max为设定的粒子最大飞翔速度向量，粒子的速度是随机改变的，但是我们不希望粒子速度不受控制的增大，从而使粒子搜索范围毫无意义地发散，因此用υ_max对粒子速度进行限制，通常粒子最大飞翔速度υ_max设定为粒子变化范围的10％～20％。进一步的，将粒子最大飞翔速度υ_max设成：长度为P且每个元素都为10的列向量。

将粒子群中第k个粒子的初始最优位置向量设为粒子群中第k个粒子的初始位置向量

利用粒子群的初始位置向量求出粒子群第k个粒子所对应的初始图像熵值求出粒子群第k个粒子所对应的初始图像熵值的过程为：

得出粒子群中第k个粒子所对应的初始目标平动分量 DP表示步骤2得出的时间矩阵。

得出粒子群中第k个粒子所对应的初始目标二维图像矩阵目标二维图像矩阵和距离频域的数据矩阵S₁具有相同的大小；目标二维图像矩阵的第q行第l列的元素为：

$I_k^{P,0}(q,l)=\frac{1}{\mathrm{MN}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{S}_M(n,m)\exp \left\{j\frac{4\mathrm{\π}((n-1){\mathrm{\Δf}}_r+f_c)}{C}{X}_k^{P,0}\left(m\right)\right\}\exp \left\{j\frac{2\mathrm{\π}}{N}(n-1)(q-1)\right\}\exp \left\{j\frac{2\mathrm{\π}}{M}(m-1)(l-1)\right\}$

其中，q取1到N，l从1到M，M为距离频域的数据矩阵S₁的列数，N为距离频域的数据矩阵S₁的行数；S_M(n,m)为数据矩阵S₁第n行第m列的元素，C是光速，Δf_r为逆合成孔径雷达的距离频谱分辨率，f_c为逆合成孔径雷达发射信号的载频；为粒子群中第k个粒子所对应的初始目标平动分量的第m个元素。

得出粒子群中第k个粒子所对应的初始图像矩阵的总能量

$S_k^{P,0}=\underset{l=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{\left|I_k^{P,0}(q,l)\right|}^2$

其中，q取1到N，l从1到M，为粒子群中第k个粒子所对应的初始目标二维图像矩阵中第q行第l列的元素，|·|表示取模值。

计算得出粒子群中第k个粒子所对应的初始目标二维图像矩阵的图像熵值

$E_k^{P,0}=-\underset{l=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|I_k^{P,0}(q,l)\right|}^2}{S_k^{P,0}}\ln \frac{{\left|I_k^{P,0}(q,l)\right|}^2}{S_k^{P,0}}$

其中，|·|表示取模值。

将粒子群的初始最优位置向量z^P,0,best设为目标平动分量多项式系数向量b_P。

(3.2)通过下列公式计算得到第u次迭代后粒子群中第k个粒子的速度向量

${\mathrm{\υ}}_k^{P,u}={\mathrm{\φ\υ}}_k^{P,u-1}+c_1{r}_1(z_k^{P,u-1,\mathrm{best}}-z_k^{P,u-1})+c_2{r}_2(z^{P,u-1,\mathrm{best}}-z_k^{P,u-1})$

其中，φ为设定的惯性系数，该值越大，粒子探索新空间的可能性就越大，该值越小，粒子遵循原来的搜索路径的可能性就越大。c₁为设定的每个粒子的个体学习因子，c₂为设定的每个粒子的社会学习因子，如果粒子群没有共享信息，即个体学习因子c₁为0，只有自身经验，及社会学习因子c₂不为0，粒子达到最优位置的几率非常小，如果粒子没有自身经验，只有社会经验，收敛速度可能会很快，但容易越过最优位置。r₁表示一个0到1之间随机数(服从均匀分布)，r₂表示一个0到1之间随机数(服从均匀分布)。本发明实施例中，将φ经验性地设为0.8，将c₁设为2，将c₂设为2。

根据以下公式得出第u次迭代后粒子群中第k个粒子的位置向量

$z_k^{P,u}=z_k^{P,u-1}+{\mathrm{\υ}}_k^{P,u}$

步骤4，得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵的图像熵值

其具体子步骤为：

(4.1)得出第u次迭代后粒子群中第k个粒子所对应的目标平动分量 D_P表示步骤2得出的时间矩阵。

(4.2)得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵目标二维图像矩阵和距离频域的数据矩阵S₁具有相同的维度；目标二维图像矩阵的第q行第l列的元素为：

$I_k^{P,u}(q,l)=\frac{1}{\mathrm{MN}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{S}_M(n,m)\×\exp \left\{j\frac{4\mathrm{\π}({\mathrm{n\Δf}}_r+f_c)}{C}{X}_k^{P,u}\left(m\right)\right\}\×\exp \left\{j\frac{2\mathrm{\π}}{N}(n-1)(q-1)\right\}\×\exp \left\{j\frac{2\mathrm{\π}}{M}(m-1)(l-1)\right\}$

其中，q取1到N，l从1到M，M为距离频域的数据矩阵S₁的列数，N为距离频域的数据矩阵S₁的行数；S_M(n,m)为数据矩阵S₁第n行第m列的元素，C是光速，Δf_r为逆合成孔径雷达的距离频谱分辨率，f_c为逆合成孔径雷达发射信号的载频；为目标平动分量的第m个元素。

(4.3)得出第u次迭代后粒子群中第k个粒子所对应的图像矩阵的总能量

$S_k^{P,u}=\underset{l=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}{\left|I_k^{P,u}(q,l)\right|}^2$

其中，q取1到N，l从1到M，为第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵中第q行第l列的元素，|·|表示取模值。

(4.4)计算得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵的图像熵值

$E_k^{P,u}=-\underset{l=1}{\overset{M}{\mathrm{\Σ}}}\underset{q=1}{\overset{N}{\mathrm{\Σ}}}\frac{{\left|I_k^{P,u}(q,l)\right|}^2}{S_k^{P,u}}\·\ln \frac{{\left|I_k^{P,u}(q,l)\right|}^2}{S_k^{P,u}}$

其中，|·|表示取模值。

步骤5，得出第u次迭代后粒子群所对应的图像最小熵值E^P,u,best、以及将第u次迭代后粒子群的最优位置向量z^P,u,best。

其具体子步骤为：

(5.1)对第u-1次迭代后粒子群中第k个粒子所对应的图像熵值和第u次迭代后粒子群中第k个粒子所对应的图像熵值进行比较，如果则第u次迭代后粒子群中第k个粒子所对应的图像最小熵值第u次迭代后粒子群中第k个粒子的最优位置向量如果则第u次迭代后粒子群中第k个粒子所对应的图像最小熵值第u次迭代后粒子群中第k个粒子的最优位置向量 $z_k^{P,u,\mathrm{best}}=z_k^{P,u-1}.$

(5.2)令k依次取1至K，得出第u次迭代后所有粒子所对应的图像最小熵值得出第u次迭代后所有粒子的最优位置向量K表示设定的粒子群的粒子总数。

(5.3)对第u次迭代后所有粒子所对应的图像最小熵值进行从大到小的排序，将第u次迭代后所有粒子所对应的图像最小熵值的最小值作为第u次迭代后粒子群所对应的图像最小熵值E^P,u,best，将第u次迭代后所有粒子所对应的图像最小熵值的最小值对应的粒子的最优位置向量作为第u次迭代后粒子群的最优位置向量z^P,u,best。

步骤6，判断迭代次数u是否大于迭代次数u的经验迭代门限u_max或者所有粒子速度的平均值是否小于所有粒子速度平均经验值门限ξ；

其具体子步骤为：

设定u_max为迭代次数u的经验迭代门限，如果u_max太小，单个粒子所对应的图像最小熵达不到最小就迭代终止，以至于得不到与真实值最接近的目标平动分量多项式系数向量b_P；如果u_max太大，迭代次数太多，凭空增加计算量。本发明实施例中，u_max取为100。

设定ξ为所有粒子速度平均经验值门限，如果粒子位置已达到所能达到的最优位置，即图像熵最小，则粒子的速度就会下降，为了避免增加计算量，在所有粒子速度平均值小于一定值时，认为粒子位置已达到最优位置，ξ一般取0～0.1之间的数。本发明实施例中，经验性地将ξ取为0.05。

如果u>u_max或者第u次迭代后所有粒子速度的平均值得出目标平动分量多项式阶数为P时的图像最小熵E_P、以及目标平动分量多项式阶数为P时的目标平动分量多项式系数向量b_P，E_P＝E^P,u,best，b_P＝z^P,u,best，|·|₁表示1范数；否则，如果u≤u_max且第u次迭代后所有粒子速度的平均值令u的值自增1，返回至步骤3，重复执行步骤3至步骤6；k＝1,2,…K，K表示设定的粒子群的粒子总数。

步骤7，令目标平动分量多项式阶数P的值自增1，重复执行步骤2至步骤6；然后，判断E_P和E_P-1的大小关系，如果E_P≤E_P-1，则返回至步骤2；如果E_P>E_P-1，则得到最终图像最小熵、以及最终图像最小熵对应的目标平动多项式系数向量，所述最终图像最最小熵是目标平动分量多项式阶数为P-1时的图像最小熵E_P-1，所述最终图像最小熵对应的目标平动多项式系数向量是：目标平动分量多项式阶数为P-1时的目标平动分量多项式系数向量b_P-1；

步骤8，将时间矩阵D_P-1和最终图像最小熵对应的目标平动多项式系数向量的转置相乘得到目标平动分量R；用目标平动分量R对距离频域匹配滤波之后的数据矩阵S_M进行平动补偿得到目标平动分量补偿后数据矩阵A，对目标平动分量补偿后数据矩阵A做基于最终图像最小熵的自聚焦，并对自聚焦形成的图像在方位向做快速傅里叶变换，得到最终ISAR图像。

其具体子步骤为：

(8.1)得出目标平动分量R：R＝D_P-1×b_P-1，其中，D_P-1为步骤2得出的时间矩阵。

(8.2)得出目标平动分量补偿后数据矩阵A：

$A=S_M\exp \left(j\frac{4\mathrm{\π}(f_c+f_r)}{C}R\right)$

其中，f_c为逆合成孔径雷达发射信号的载频，f_r为逆合成孔径雷达的距离频率， $f_r=[-\frac{{\mathrm{N\Δf}}_r}{2},-\frac{(N-1){\mathrm{\Δf}}_r}{2},...,\frac{(N-1){\mathrm{\Δf}}_r}{2}],$ Δf_r为逆合成孔径雷达的距离频谱分辨率，N为距离频域的数据矩阵S₁的行数，C是光速。

(8.3)对目标平动分量补偿后数据矩阵A做基于最终图像最最小熵的自聚焦，并对自聚焦形成的图像在方位向做快速傅里叶变换，得到最终ISAR图像。

本发明的效果可以通过以下仿真实验进行说明。

1)仿真条件

原始回波数据为卫星点目标的仿真数据，逆合成孔径雷达发射信号为X波段，带宽为1GHz，脉冲重复频率为1000Hz，相干积累时间为0.128s，距离和方位的采样点数分别为512和128点。

2)仿真内容

首先，针对原始回波数据采用本发明进行ISAR成像，在ISAR成像的过程中，使用3阶多项式形式的平动(对原始回波数据进行距离压缩，然后给压缩后数据加上3阶多项式形式的平动)，归一化的多项式系数为b＝[25.6 -1.6384 -0.2097]。

具体地，对仿真数据分别加上信噪比为0dB、-5dB的复高斯白噪声，然后针对加上噪声的仿真数据采用本发明进行ISAR成像。参照图2a，为仿真实验中信噪比为0dB时的仿真数据距离包络示意图，参照图2b，为仿真实验中信噪比为0dB时仿真数据的理想ISAR成像结果示意图，参照图2c，为仿真实验中信噪比为0dB时对仿真数据采用本发明进行ISAR成像得出的ISAR成像结果示意图；参照图2d，为仿真实验中信噪比为-5dB时的仿真数据距离包络示意图，参照图2e，为仿真实验中信噪比为-5dB时仿真数据的理想ISAR成像结果示意图，参照图2f，为仿真实验中信噪比为-5dB时对仿真数据采用本发明进行ISAR成像得出的ISAR成像结果示意图。图2a至图2f中，横轴表示方位维，纵轴表示距离维。从图2a和图2d中看出，仿真的原始回波数据的信噪比比较低，目标回波已经被噪声淹没，因此用传统的相邻相关方法进行包络对齐会失效，导致后续的自聚焦方法就无从应用，从而得到聚焦良好的成像结果。而本发明充分利用距离压缩和方位压缩的二维信噪比增益，通过估计目标平动进行包络对齐和自聚焦联合校正，大大降低了ISAR成像所需的信噪比门限。对比图2b和图2c，并对比图2e和图2f，可以发现本发明得出ISAR成像结果与理想的SAR成像结果相差很小。

参照图3a，为仿真实验中信噪比为0dB时采用本发明得出的最终ISAR图像熵值和目标平动分量多项式阶数的关系示意图，图3a中，横轴表示目标平动分量多项式阶数，纵轴表示最终ISAR成像结果的熵值。参照图3b为仿真实验中信噪比为0dB且目标平动分量多项式阶数为3时得出的各阶多项式系数和迭代次数u的关系示意图，图3b中，横轴表示迭代次数，纵轴表示多项式系数，b1表示目标平动分量多项式阶数为3的第1阶多项式系数(平方项)、b2表示目标平动分量多项式阶数为3的第2阶多项式系数(一次项)、b3表示目标平动分量多项式阶数为3的第3阶多项式系数(常数项)。参照图3c，为仿真实验中信噪比为-5dB时采用本发明得出的最终ISAR图像熵值和目标平动分量多项式阶数的关系示意图，图3c中，横轴表示目标平动分量多项式阶数，纵轴表示最终ISAR成像结果的熵值。参照图3d为仿真实验中信噪比为-5dB且目标平动分量多项式阶数为3时得出的各阶多项式系数和迭代次数u的关系示意图，图3d中，横轴表示迭代次数，纵轴表示多项式系数，b1表示目标平动分量多项式阶数为3的第1阶多项式系数(平方项)、b2表示目标平动分量多项式阶数为3的第2阶多项式系数(一次项)、b3表示目标平动分量多项式阶数为3的第3阶多项式系数(常数项)。

从图3a和图3c中看出，熵值变化曲线均在3处(目标平动分量多项式阶数为3)存在一个拐点，即在阶数为3时对目标平动的估计最准确，最终的成像结果熵值最小。从图3b和图3d中看出，当目标平动分量多项式阶数为3时，在不同的信噪比下，本发明采用的多项式系数都能很快收敛到最优解上，并有着较高的估计精度。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (6)

1.一种基于粒子群滤波优化的ISAR成像方法，其特征在于，包括以下步骤：

步骤1，利用逆合成孔径雷达发射线性调频信号S₀，利用逆合成孔径雷达接收原始回波数据，将原始回波数据沿距离向进行快速傅里叶变换，得到距离频域的数据矩阵S₁，对距离频域的数据矩阵S₁进行匹配滤波，得到匹配滤波之后的数据矩阵S_M；构建长度为M的时间列向量t，t＝[-M/2，-M/2+1，…，M/2-1]^T，上标T表示矩阵或向量的转置，M为距离频域的数据矩阵S₁的列数；设定目标平动分量多项式阶数P，P＝1，2，3...，当P＝1时，跳至步骤2；

步骤2，构建大小为M×P的时间矩阵D_P，时间矩阵D_P的第g列为列向量t^g，列向量t^g表示对时间列向量t的每个元素取g次方得出的新的向量，g取1至P；构建目标平动分量多项式系数向量b_P，当P＝i时，目标平动分量多项式系数向量b_P为0；当P＞1时，目标平动分量多项式系数向量b_P是长度为P的列向量，目标平动分量多项式系数向量b_P中的元素全为0；令迭代次数u＝1，2，3...，当u＝1时，跳至步骤3；

步骤3，设定粒子群的粒子总数为K，令k＝1，2，…K，设定粒子群中第k个粒子的初始位置向量粒子群中第k个粒子的初始速度向量粒子群中第k个粒子的初始最优位置向量以及粒子群第k个粒子所对应的初始图像熵值根据第u-1次迭代后粒子群中第k个粒子的速度向量第u-1次迭代后粒子群中第k个粒子的位置向量以及第u-1次迭代后粒子群中第k个粒子的最优位置列向量得到第u次迭代后粒子群中第k个粒子的速度向量得出第u次迭代后粒子群中第k个粒子的位置向量

步骤4，得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵的图像熵值

步骤5，得出第u次迭代后粒子群所对应的图像最小熵值E^P，u，best、以及第u次迭代后粒子群的最优位置向量z^P，u，best；

步骤6，用u_max表示设定的迭代次数u的经验迭代门限，如果u＞u_max或得出目标平动分量多项式阶数为P时的图像最小熵E_P、以及目标平动分量多项式阶数为P时的目标平动分量多项式系数向量b_P，E_P＝E^P，u，best，b_P＝z^P，u，best，|·|₁表示1范数，ξ为设定的粒子速度平均经验值门限；否则，令u的值自增1，返回至步骤3；

步骤7，令目标平动分量多项式阶数P的值自增1，重复执行步骤2至步骤6；然后，判断E_P和E_P-1的大小关系，如果E_P≤E_P-1，则返回至步骤2；如果E_P＞E_P-1，则得到最终图像最小熵、以及最终图像最小熵对应的目标平动多项式系数向量，所述最终图像最小熵是目标平动分量多项式阶数为P-1时的图像最小熵E_P-1，所述最终图像最小熵对应的目标平动多项式系数向量是：目标平动分量多项式阶数为P-1时的目标平动分量多项式系数向量b_P-1；

步骤8，将时间矩阵D_P-1和最终图像最小熵对应的目标平动多项式系数向量的转置相乘得到目标平动分量R；用目标平动分量R对匹配滤波之后的数据矩阵S_M进行平动补偿得到目标平动分量补偿后数据矩阵A，对目标平动分量补偿后数据矩阵A做基于最终图像最小熵的自聚焦，并对自聚焦形成的图像在方位向做快速傅里叶变换，得到最终ISAR图像。

2.如权利要求1所述的一种基于粒子群滤波优化的ISAR成像方法，其特征在于，所述步骤1的具体子步骤为：

(1.1)利用逆合成孔径雷达发射线性调频信号S₀，利用逆合成孔径雷达接收原始回波数据；对逆合成孔径雷达发射的线性调频信号S₀进行距离向快速傅里叶变换，得到线性调频信号S₀的频谱矩阵S₂；

(1.2)对原始回波数据在距离向进行快速傅里叶变换，得到距离频域的数据矩阵S₁，矩阵S₁的列数表示为M，行数表示为N；

(1.3)对距离频域的数据矩阵S₁进行匹配滤波，得到匹配滤波之后的数据矩阵S_M；

(1.4)构建长度为M的时间列向量t，t＝[-M/2，-M/2+1，…，M/2-1]^T，上标T表示矩阵或向量的转置，M为距离频域的数据矩阵S₁的列数；设定目标平动分量多项式阶数P，P＝I，2，3...，当P＝1时，跳至步骤2。

3.如权利要求1所述的一种基于粒子群滤波优化的ISAR成像方法，其特征在于，所述步骤3的具体子步骤为：

(3.1)设定粒子群的粒子总数为K，令k＝1，2，…K，设定粒子群中第k个粒子的初始位置向量

$z_k^{P,0}=b_p+4\×\[rand(0,1)-0.5\]$

其中，rand(0，1)表示产生一个0到1之间均匀分布的随机数；

设定粒子群中第k个粒子的初始速度向量

$v_k^{P,0}=v_{max}\×rand(0,1)$

其中，υ_max为设定的粒子最大飞翔速度向量；

将粒子群中第k个粒子的初始最优位置向量设为粒子群中第k个粒子的初始位置向量设定粒子群第k个粒子所对应的初始图像熵值将粒子群的初始最优位置向量设为目标平动分量多项式系数向量b_P；

(3.2)通过下列公式计算得到第u次迭代后粒子群中第k个粒子的速度向量

$v_k^{P,u}={\mathrm{\φv}}_k^{P,u-1}+c_1{r}_1(z_k^{P,u-1,best}-z_k^{P,u-1})+c_2{r}_2(z^{P,u-1,best}-z_k^{P,u-1})$

其中，φ为设定的惯性系数，c₁为设定的每个粒子的个体学习因子，c₂为设定的每个粒子的社会学习因子，r₁表示一个0到1之间随机数，r₂表示一个0到1之间随机数；

根据以下公式得出第u次迭代后粒子群中第k个粒子的位置向量

$z_k^{P,u}=z_k^{P,u-1}+v_k^{P,u}.$

4.如权利要求1所述的一种基于粒子群滤波优化的ISAR成像方法，其特征在于，所述步骤4的具体子步骤为：

(4.1)得出第u次迭代后粒子群中第k个粒子所对应的目标平动分量 D_P表示步骤2得出的时间矩阵；

(4.2)得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵目标二维图像矩阵和距离频域的数据矩阵S₁具有相同的维度；目标二维图像矩阵的第q行第l列的元素为：

$I_k^{P,u}(q,l)=\frac{1}{MN}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{S}_M(n,m)\×\exp \left\{j\frac{4\mathrm{\π}({\mathrm{n\Δf}}_r+f_c)}{C}{X}_k^{P,u}\right(m\left)\right\}\×\exp \left\{j\frac{2\mathrm{\π}}{N}(n-1)(q-1)\right\}\×\exp \left\{j\frac{2\mathrm{\π}}{M}(m-1)(l-1)\right\}$

其中，q取1到N，l从1到M，M为距离频域的数据矩阵S₁的列数，N为距离频域的数据矩阵S₁的行数；S_M(n，m)为数据矩阵S₁第n行第m列的元素，C是光速，Δf_r为逆合成孔径雷达的距离频谱分辨率，f_c为逆合成孔径雷达发射信号的载频；为目标平动分量的第m个元素；

(4.3)得出第u次迭代后粒子群中第k个粒子所对应的图像矩阵的总能量

$S_k^{P,u}=\underset{l=1}{\overset{M}{\Σ}}\underset{q=1}{\overset{N}{\Σ}}|I_k^{P,u}(q,l){|}^2$

其中，q取1到N，l从1到M，|·|表示取模值；

(4.4)计算得出第u次迭代后粒子群中第k个粒子所对应的目标二维图像矩阵的图像熵值

$E_k^{P,u}=-\underset{l=1}{\overset{M}{\Σ}}\underset{q=1}{\overset{N}{\Σ}}\frac{|I_k^{P,u}(q,l){|}^2}{S_k^{P,u}}\·ln\frac{|I_k^{P,u}(q,l){|}^2}{S_k^{P,u}}.$

5.如权利要求1所述的一种基于粒子群滤波优化的ISAR成像方法，其特征在于，所述步骤5的具体子步骤为：

(5.1)对第u-1次迭代后粒子群中第k个粒子所对应的图像熵值和第u次迭代后粒子群中第k个粒子所对应的图像熵值进行比较，如果则第u次迭代后粒子群中第k个粒子所对应的图像最小熵值第u次迭代后粒子群中第k个粒子的最优位置向量如果则第u次迭代后粒子群中第k个粒子所对应的图像最小熵值第u次迭代后粒子群中第k个粒子的最优位置向量

(5.2)令k依次取1至K，得出第u次迭代后所有粒子所对应的图像最小熵值得出第u次迭代后所有粒子的最优位置向量K表示设定的粒子群的粒子总数；

(5.3)将第u次迭代后所有粒子所对应的图像最小熵值的最小值作为第u次迭代后粒子群所对应的图像最小熵值E^P，u，best，将第u次迭代后所有粒子所对应的图像最小熵值的最小值对应的粒子的最优位置向量作为第u次迭代后粒子群的最优位置向量z^P，u，best。

6.如权利要求1所述的一种基于粒子群滤波优化的ISAR成像方法，其特征在于，所述步骤8的具体子步骤为：

(8.1)得出目标平动分量R：R＝D_P-1×b_P-1，其中，D_P-1为步骤2得出的时间矩阵；

(8.2)得出目标平动分量补偿后数据矩阵A：

$A=S_M\exp \left(j\frac{4\mathrm{\π}(f_c+f_r)}{C}R\right)$

其中，f_c为逆合成孔径雷达发射信号的载频，f_r为逆合成孔径雷达的距离频率，f_r与f_c相加的含义是：f_r的每一个元素与f_c相加后得到的向量，Δf_r为逆合成孔径雷达的距离频谱分辨率，N为距离频域的数据矩阵S₁的行数，C是光速；

(8.3)对目标平动分量补偿后数据矩阵A做基于最终图像最最小熵的自聚焦，并对自聚焦形成的图像在方位向做快速傅里叶变换，得到最终ISAR图像。