基于L统计量和粒子群优化的微多普勒抑制方法与流程

本发明涉及逆合成孔径雷达成像技术领域，具体为微多普勒效应消除后的刚体恢复方法。

背景技术：

逆合成孔径雷达(inversesyntheticapertureradar,isar)成像时，通常情况将目标假设为刚体，通过运动补偿将雷达远场的物体转化为理想转台后，使用距离多普勒算法对其距离和方向维度分别依次进行压缩得到高分辨率的雷达图像。但是，当目标拥有高速旋转的部件时，该高速旋转部件将破坏刚体假设，如果仍然通过距离多普勒算法(range-doppler，r-d)成像的话会导致图像域中产生不同距离单元的条带干扰，严重影响成像结果，这种高速旋转部件产生的效应就是微多普勒效应。所以，若使用r-d算法对有高速自旋部件的目标成像的话，需要在方位压缩前进行去除微多普勒的处理。

消除微多普勒的方法很多，通过对含有微多普勒成分的距离单元进行基于时频分析的处理方法能很好的将非平稳的微多普勒成分和平稳的刚体成分分开。其中黑山大学的l.stankovic教授通过l统计量的时频域排序消除算法，能高效的除去微多普勒分量。纵使l统计量拥有简单快速的优点，但缺点也很明显，即l统计量方法会在消除微多普勒分量的同时也消除部分刚体的能量，导致重建后刚体的能量集中度下降，最终导致雷达成像的结果分辨率下降。

由带有自旋部件的isar回波建模可知，旋转部件对应的多普勒频率为正弦曲线，而刚体部件为恒定频率。对一个单频信号与一个正弦调频信号的叠加信号采用l统计量方法的结果如图1所示。

其中图1(a)为存在微多普勒效应的短时傅里叶变换(short-timefouriertransform，stft)表示的时频图，可以看出单频信号所对应的刚体部分被正弦调频信号的频谱所干扰，不能清晰地分离出单频信号成分。图1(b)为l统计量的排序结果，随后选取适合的门限值，将幅度较大的元素予以剔除，剔除后的结果见图1(c)，可见正弦调频信号被较好的消除，但刚体部件的能量有所损失，这将导致重建后刚体的能量分散，影响成像效果。

为了提高重建后刚体的能量集中度，l.stankovic提出了使用遗传算法(geneticalgorithm，ga)恢复被消除的刚体，但是复杂度极大，运算时间过长，掩盖了l统计量原有优势。本发明基于上述背景提出了在在不显著提高数据处理时间的前提下，对损失的刚体成分进行时频域的高效恢复方法，即基于粒子群优化(particleswarmoptimization，pso)的刚体恢复算法。

技术实现要素：

本发明的目的是：针对现有技术中进行刚体损失恢复效率低的问题，提供一种基于l统计量和粒子群优化的微多普勒抑制方法。

本发明为了解决上述技术问题采取的技术方案是：

基于l统计量和粒子群优化的微多普勒抑制方法，包括以下步骤：

步骤一：获取原始信号的短时傅里叶变换结果，即stft值；

步骤二：对stft值的绝对值进行排序，根据排序结果，在每个频率单元由高到低剔除一定比例的stft值，此比例门限不低于50％；

步骤三：找到stft值最大值点，将该位置的值及其两侧的d个采样单元的值去除；

步骤四：计算剩余能量；

步骤五：判断剩余能量是否大于阈值，若剩余能量值大于阈值，则kk＝kk+1，kk为刚体数目的计数器，并执行步骤三，若剩余能量值不大于阈值，则执行步骤六；

步骤六：判断kk是否大于零，若kk不大于零，结束，若kk大于零，则将原本剩余stft值的中位数作为刚体的估计模值；

步骤七：初始化一系列粒子群，随机从[0,2π)中选择相位；

步骤八：计算每个粒子的适应性函数，根据适应度选择最佳粒子速度和位置，并判断迭代次数是否达到迭代阈值或适应性函数保持不变，若不符合，则根据位置和速度公式继续迭代更新下一代粒子群，若符合，则更新位置和速度信息，所述位置和速度公式为：

vid(t+1)＝u*vid(t)+c1*rand()*

[pid(t)-xid(t)]+c2*rand()*

[pgd(t)-xid(t)]

xid(t+1)＝xid(t)+vid(t+1)，1≤i≤n,1≤d≤d

其中，搜索空间为d维度，粒子总数为n，xi＝(xi1,xi2,...,xid)表示为第i个粒子的位置，vid表示速度，pi＝(pi1,pi2,...,pid)对应于第i个粒子历史经历过的位置中的最佳位置，pg是pi(i＝1,...,n)中全局最优位置，序号为g，c1,c2分别是符号为正的常数，即加速因子，rand()是一个0到1之间的随机数，u为惰性权因子；

步骤九：获取适应度最高的粒子群最为最后的估计结果，然后将所有恢复的样本沿时间维度求和，得到重建的刚体样本的ft，并将其作为最终刚体成分w1的频谱sw1(k)的估计，令kk←kk-1并执行步骤六。

进一步的，所述步骤一中短时傅里叶变换公式如下

其中，窗函数w(i)长度为mw,s(i)是长度为m的信号，m和k分别是离散时间和频率的自变量，i是累加符号中的累加自变量。

进一步的，所述步骤三的详细步骤为：

步骤三一：对剩余的stft值求和，得到重建的傅立叶变换结果sl(k)；

步骤三二：设置一个阈值ε，去掉m-d分量之后的信号能量的5％，并计算刚体反射器的长度，得到|d|≤1+2(m/mw-1)，初始刚体反射器数目kk初始化为1，查找并记录sl(k)对应的最大值，即刚体反射器位置，然后移除该最大值位置处及其两侧各d个采样单元的值。

进一步的，所述步骤三一中重建的傅立叶变换结果sl(k)的公式为：

其中，χk(m)是频率k在移除完最强的m-mp个stft点后剩下的mp个stft值。

进一步的，所述步骤六中，若kk大于零，则将原本剩余stft值的中位数作为刚体的估计模值的详细步骤为：若kk＞0，假设在频率上没有零值，在刚体点的每一侧，形成一个初始粒子数作为第一代，其中每个粒子的位置代表丢失的非零刚体值的相位，每个相位值从集合[0,2π)中随机选择，将每个粒子的适应度计算为恢复完后刚体信号的总能量，其公式为：

其中，e是信号的能量，kr是刚体的位置，d是刚体信号累加的范围。

进一步的，所述初始粒子数为50、100或200个粒子。

进一步的，所述初始粒子数为100个粒子。

进一步的，所述步骤八中迭代阈值为100次迭代。

进一步的，所述步骤二中在每个频率单元由高到低剔除一定比例的stft值，此比例为50％。

本发明的有益效果是：与l统计量方法相比，本发明在不显著增加计算时间的基础上大幅度提高了重建ft的频谱聚集度，与基于遗传算法的刚体恢复方法比又大量节省了计算时间，本发明可以高效率的恢复并重建刚体频谱。

附图说明

图1(a)为存在微多普勒效应的短时傅里叶变换时频图，图1(b)l统计量的排序结果图，图1(c)为l统计量将幅度较大的元素予以剔除后的结果图；

图2为本发明的整体流程图；

图3(a)为模拟信号的stft结果，图3(b)为ft模拟信号的结果；

图4(a)与图4(d)为l统计量方法的结果，图4(b)与图4(e)ga方法的结果，图4(c)与图4(f)为pso方法的结果；

图5(a)为原始信号的stft结果，图5(b)为原始信号的ft结果；

图6(a)为l统计量方法移除后的信号stft图,图6(b)是遗传算法下的stft图，图6(c)是pso算法下的stft图，图6(d)为l统计量方法移除后的信号ft图,图6(e)是遗传算法下的ft图，图6(f)是pso算法下的ft图；

图7为飞机散射点模型；

图8为直接采用rd算法的成像结果；

图9为第103距离门的stft结果；

图10(a)为采用l统计量算法恢复刚体后重建的r-d成像结果，图10(b)为采用ga算法恢复刚体后重建的r-d成像结果,图10(c)为采用pso算法恢复刚体后重建的r-d成像结果；

图11为安-26飞机外形图；

图12为直接rd算法的成像结果；

图13为第150距离门的stft结果；

图14(a)为采用l-统计量算法的成像结果，图14(b)为采用ga算法的成像结果,图14(c)为采用pso算法的成像结果。

具体实施方式

具体实施方式一：参照图2具体说明本实施方式，本实施方式利用粒子群优化(pso)算法恢复刚体并进行重建刚体图像。其具体过程为：

步骤一：计算原始信号的短时傅里叶变换。由此就将原始信号从时域转换到时频域。

短时傅里叶变换公式如下

其中窗函数w(i)长度为mw,s(i)是长度为m的信号。

步骤二：对stft值的绝对值进行排序，并根据每个频率值对应的stft和l统计量移除适当百分比，一般不低于50％(此处为50％，50％为经验值，一般大于等于50％的情况下能除掉较大比例的微多普勒成分，stft的l统计量就是stft沿着每个频率排序的结果)的最高stft样本值。

步骤三：对剩余的stft值求和，得到重建的傅立叶变换(fouriertransform，ft)结果sl(k)

其中，χk(m)是频率k在移除完最强的m-mp个stft点后剩下的mp个stft值。在这里，我们设置一个阈值ε为去掉m-d(m-d是微多普勒的英文micro-doppler的缩写)分量之后的信号能量的5％，并计算刚体反射器的长度，可以得到|d|≤1+2(m/mw-1)。初始刚体反射器数目kk初始化为1。查找并记录sl(k)对应的最大值也就是刚体反射器位置。然后，移除该最大值位置处及其两侧各d个采样单元的值。

步骤四：计算剩余能量，如果能量值仍然高于阈值则令kk＝kk+1(kk是刚体数目的计数器，每检测出一个刚体成分就自动加一)。同样，类似在步骤2中，查找并记录sl(k)的最大值。然后移除位置处的值，并将d样本放在该位置的每一侧。

步骤五：如果kk＞0，说明尚有需恢复的刚体成分，缺失的刚体值估计为剩余值沿时间维度的中位数。在刚体点的每一侧，假设在频率上没有零值，我们可以形成一个初始粒子数，一般为50、100或者200个粒子(此处为100个粒子)作为第一代，其中每个粒子的位置代表丢失的非零刚体值的相位。每个相位值可以从集合[0,2π)中随机选择。将每个粒子的适应度计算为恢复完后刚体信号的总能量，如下所示

为了得到缺失stft值的估计值，我们将它们设置为剩余stft值的中值。如我们所知，在越接近在kr处被零值环绕的峰值意味着更接近重建刚体ft的最佳解。因此，能量较低的粒子具有较高的适应度。

步骤六：迭代次数到达最大迭代次数100次或直到最高适应度保持不变前，根据适应度函数选择最佳粒子速度和位置。搜索空间为d维度，粒子总数为n，xi＝(xi1,xi2,...,xid)表示为第i个粒子的位置，vid表示速度，速度和位置的更新公式如下所示

xid(t+1)＝xid(t)+vid(t+1)，(1≤i≤n,1≤d≤d)(5)

其中pi＝(pi1,pi2,...,pid)对应于第i个粒子历史经历过的位置中的最佳位置，而pg则是pi(i＝1,...,n)中全局最优位置，序号为g。c1,c2分别是符号为正的常数，被称作加速因子。rand()是一个0到1之间的随机数。u为惰性权因子。

步骤七：为了得到丢失刚体相位的最终估计结果，我们从最后一代粒子中选择适应度最高的粒子，这里丢失的刚体值是非零的。最后，我们将所有恢复的样本沿时间维度求和，得到重建的刚体样本的ft，并将其作为最终sw1(k)的估计。令kk←kk-1并进入第5步。

与l统计量方法相比，本发明在不显著增加计算时间的基础上大幅度提高了重建ft的频谱聚集度，与基于遗传算法的刚体恢复方法比又大量节省了计算时间，达到了高效率恢复并重建刚体频谱的效果。

同样作为基于生物演化的多维优化算法，ga算法和pso算法由于都是随机初始化，所以迭代次数不一定，无法定量分析算法复杂程度。通常，对于已经确定种群大小的ga算法和确定了粒子群规模的pso算法，适应度函数被计算的次数为种群规模和迭代数的乘积。迭代的具体次数没办法直接给出。假设这两个算法均为迭代100次停止迭代或者直到连续五个子代最大适应度保持不变停止迭代，尽管代数无法准确确定，但是最坏的情况可以确定，即迭代100次的时候，假设100个个体。下列定性分析中，假设初始化的个体(pso中的粒子群，ga中的种群)数目为100，迭代次数为100。两个算法的适应性函数一样都是能量最低准则。假设刚体数目kk＝3，m＝256,mw＝32。所以两个算法适应度函数的计算次数都是10000.现在来比较一次迭代时两者的计算量。

pso算法：对于pso算法，每个个体大小为11520个相位值。每个个体要更新一次速度和位置，总共包含5次乘除运算和5次加减运算，更新完后直接对数据求适应度函数即可。

ga算法：对于ga算法，相位先得经过三位二进制的编码才能得到个体，个体大小为11520*3＝34560。选择过程中，每次20个个体被选中参与到竞赛中，经过19次比较得到新的个体，所以总共为1900次比较，或者布尔运算。交叉过程中，70％的数据进行了交叉运算，交叉过程中，两个个体的每一位都被随机的分配给新的子代，一次分配相当于两次二进制加法，每位都进行，即进行34560*2*34560*0.7/2＝836075520次二进制加法。变异过程中，有百分之25的概率每50位就有一位数变异，也就是说每次平均执行17280次二进制加法。

最后，计算适应度函数前需要把得到的个体进行二进制解码，解码为十进制浮点数才能计算能量大小。

从上述pso和ga算法的比较中，可以显而易见的看到，ga算法的运算次数和步骤都远大于pso算法，导致实际仿真与实测数据验证中，ga算法的运行时间也远大于pso算法。

实施例一：

本实施例子中使用最简单的信号仿真。信号为一个平稳刚体信号和一个非平稳信号的叠加信号，其中平稳信号代表刚体成分，而非平稳信号为正弦调制的信号，代表微多普勒成分。该模拟信号的时域表达式为

图3(a)与图3(b)为该信号的stft结果与ft的结果，点数n＝512，可见存在严重的微多普勒效应，无法从ft结果中分离出刚体成分。

对该模拟信号分别采用l统计量，ga算法以及本发明提出的pso算法进行处理，结果如图4所示。

图4(a)与图4(d)为l统计量方法的结果，图4(b)与图4(e)ga方法的结果，图4(c)与图4(f)为pso方法的结果。可见经过ga和pso后刚体部分得到了很好的恢复，恢复效果相近，表1给出了两种方法的运算时间对比，可见pso方法在运行时间上具有明显的优势。

表1运行时间对比

实施例二：

本实施例子中使用的信号为三个平稳刚体信号和三个非平稳信号的和信号，其中平稳信号代表刚体成分，而非平稳信号为正弦调制的信号，代表微多普勒成分。该模拟信号的时域表达式为

图5(a)为原始信号的stft结果，可以看到刚体信号的频率不随时间变化，而被随时间呈正弦变化的微多普勒成分部分的遮盖了。图5(b)为原始信号的ft结果，可以看到，在频谱图中不能区分刚体和旋转部件的成分。

图6(a)为l统计量方法移除后的信号stft图，去除每个频率的最大前50％的stft模值后重新排序后得到去除m-d后的stft结果。图6(b)为去除微多普勒后信号的ft，可以看到l统计量能有效的去除微多普勒成分，但是刚体信号本身的聚集性不理想。

图6(b)和图6(e)是遗传算法下的stft和ft图，图6(c)和图6(f)分别是pso算法下的stft和ft图。观察两种算法的ft图可以看到，ga算法和pso算法的刚体恢复作用效果显著。两种算法运行时间如表1所示，由结果可知，pso算法在刚体恢复效果不亚于ga算法的情况下能大幅度减少运算时间。

表2运行时间对比

实施例三：

本实施例子中将针对理想的isar点散射模型进行距离-多普勒成像。图7为点散射模型在空间中的分布。如图，模型为空间中拥有两个螺旋桨的飞机。表3为仿真参数，信噪比为0。

图8为原始回波数据r-d算法后的成像结果，可以看到，和螺旋桨位于同个距离单元的刚体散射点被因微多普勒效应产生的干扰条带遮盖了，可以看到微多普勒效应出现在了第103距离门及其附近距离门、150距离门及其附近距离门。

图9为第103距离门的stft，刚体散射点对应于时频图中的平稳信号，而旋转点对应于不平稳的正弦调频信号。

图10(a)到(c)分别为采用l统计量、ga和pso算法恢复刚体后重建的r-d成像结果，从图中可以看出，刚体恢复将明显的改善最后图像的质量，而ga算法和pso算法的结果在图像域差距并不明显。表4给出了不同信噪比下的三种微多普勒去除方法运算时间和图像熵(图像质量指标)的量化对比。

表3

表4图像熵和运行时间对比结果

实施例四：

本实施例中采用来自安-26飞机的实测数据，安-26拥有两个螺旋桨，所以其回波将带有微多普勒成分。安-26飞机的三视外形图如图11所示。

图12是不经过微多普勒消除处理直接利用r-d算法的成像结果。从图像可以看出，两条干扰条带分别对应着an-26的两个螺旋桨，而这些条带遮挡了对应距离门的飞机刚体部分图像。

图13为含有微多普勒成分的第150个距离门的stft图，由3图可知，虽然刚体成分能明显的识别出来，但是微多普勒成分因为两个螺旋桨之间的互相干扰以及频谱混叠等原因导致其正弦调频的形式无法辨别。

图14(a)到(c)分别为采用l-统计量、ga算法和pso算法的成像结果，从图中可以看出，刚体恢复将明显的改善最后图像的质量，而ga算法和pso算法的结果在图像域效果接近。三种成像方式所消耗时间和图像熵的量化对比如表5所示。可见，对于实测数据，pso与ga的成像结果近似，但pso在运行时间上比ga具有明显的优势。

表5图像熵和运行时间对比结果

需要注意的是，具体实施方式仅仅是对本发明技术方案的解释和说明，不能以此限定权利保护范围。凡根据本发明权利要求书和说明书所做的仅仅是局部改变的，仍应落入本发明的保护范围内。