本发明涉及外辐射源雷达信号处理技术领域，具体为一种基于radon变换与拉普拉斯算子的微多普勒抑制方法。
背景技术：
逆合成孔径雷达成像时，通常情况将目标假设为刚体，通过运动补偿将雷达远场的物体转化为理想转台后，使用距离多普勒算法对其距离和方向维度分别依次进行压缩得到高分辨率的雷达图像。但是，当目标拥有高速旋转的部件时，该高速旋转部件将破坏刚体假设，如果仍然通过距离多普勒算法成像的话会导致图像域中产生不同距离单元的条带干扰，使原本成像的目标，即目标在空间中的姿态遮挡，对成像结果不利，这种高速旋转部件产生的效应就是所谓的微多普勒效应。所以，欲使用距离多普勒算法(range-doppler算法)对有高速旋转物体成像的话，需要在方位压缩前进行去除微多普勒的处理。
为了区别非平稳的微多普勒以及平稳的刚体成分，通过时频分析工具进行两者的分离是非常普遍的方法。其中，短时傅里叶变换(shorttimefouriertransform)拥有可接受的时频分辨率并且不存在信号之间的交叉项，所以通过处理短时傅里叶变换的时频图可以高效地分离微多普勒和刚体成分。l.stankovic等人提出了l统计量，逆瑞登变换(inverseradontransform)等方法，利用了微多普勒在时频域的特点对微多普勒成分进行了检测和抑制，最终得到了质量较好的雷达图像。由于这类方法太过依赖微多普勒成分在时频域中形式的理想程度，所以，当雷达的回波数据被杂波，噪声污染或者因为采样率不够造成频谱混叠时，上述方法会因为微多普勒的形式被破坏导致方法的效果变差。相对于微多普勒成分，刚体成分的形式在时频域中稳定性高，不会因为上述原因使得信号形式改变。
技术实现要素：
本发明的目的是：针对现有技术中在刚体信号能量较低的情况下，或者信噪比较低的情况时不能有效的检测出刚体信号成分的问题，提出基于radon变换与拉普拉斯算子的微多普勒抑制方法。
本发明为了解决上述技术问题采取的技术方案是：
基于radon变换与拉普拉斯算子的微多普勒抑制方法，包括以下步骤：
步骤一：对原始信号进行短时傅里叶变换；
步骤二：对短时傅里叶变换的时频图进行radon变换，并将变换域中的极值点标记为可能的刚体成分；
步骤三：使用拉普拉斯算子对图像的radon变换域进行处理，得到变换域的陡峭程度信息；
步骤四：确定radon变换域中产生可能的刚体成分的峰值点；
步骤五：根据步骤四中刚体成分产生的峰值点在radon变换域中的坐标位置确定刚体在短时傅里叶时频图中的位置参数；
步骤六：根据步骤五得到的刚体成分参数设计带通滤波器，对原信号进行滤波，滤除微多普勒信号成分，得到纯净的刚体信号成分。
进一步的，所述步骤一中短时傅里叶变换公式如下：
其中，窗函数w(i)长度为mw，s(i)是原信号，m是信号长度，w(i-m)是窗函数。
进一步的，所述mw取信号长度m的八分之一。
进一步的，所述步骤二中radon变换的公式如下：
其中，f(x,y)是信号短时傅里叶变换后的时频谱图，x是时间，y是频率，δ是冲激函数，α是积分投影的角度，β是投影轴上的位置。
进一步的，所述步骤三中拉普拉斯算子表达式如下
进一步的，所述步骤四中确定radon变换域中产生可能的刚体成分的峰值点的方法为：
步骤四一：对变换域的值进行排序，设置经验百分比为1％，取变换域中的函数值为前1％的点；
步骤四二：点函数值大于等于其邻域的函数值即为radon变换域中产生刚体成分的峰值点。
进一步的，所述步骤四中确定radon变换域中产生刚体成分的峰值点通过判断步骤二中可能的刚体峰值点位置经过步骤三拉普拉斯算子处理后在该点还是否为峰值确定，若处理后在该点还是峰值，则此位置为刚体成分产生的峰值点。
进一步的，所述原始信号为雷达回波中所有含有微多普勒成分的距离单元信号。
本发明的有益效果是：
本发明利用了刚体信号在时频域中呈现为直线的特点，利用图像直线检测工具，radon变换进行刚体信号的检测。本发明利用了拉普拉斯算子，对radon变换域中的峰值进行了非刚体成分的排除，即只有足够陡峭的峰值点才被认为是极值点。通过这种方式，可以在刚体信号能量较低的情况下，或者信噪比较低的情况还是能有效的检测出刚体信号成分并相应的进行滤波。
附图说明
图1为本发明流程图；
图2(a)为原信号的变换域图像；
图2(b)为原信号的变换域图像；
图2(c)为原信号的变换域图像；
图3(a)为信号频谱图；
图3(b)为信号频谱图；
图4为飞机散射点模型；
图5为直接采用rd算法的成像结果；
图6(a)为含有微多普勒的距离单元的信号变换域图像；
图6(b)为含有微多普勒的距离单元的信号变换域图像；
图6(c)为含有微多普勒的距离单元的信号变换域图像；
图7为本发明使用rd算法进行成像的结果；
图8为安-26飞机外形结构图；
图9为直接rd算法的成像结果；
图10(a)为含有微多普勒的距离单元的信号变换域图像；
图10(b)为含有微多普勒的距离单元的信号变换域图像；
图10(c)为含有微多普勒的距离单元的信号变换域图像
图11(a)为信号的频谱图；
图11(b)为信号的频谱图；
图12为经过微多普勒去除的成像结果。
具体实施方式
具体实施方式一：参照图1具体说明本实施方式，本实施方式所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，包括以下步骤：
步骤一：计算原始信号的短时傅里叶变换(stft)。由此就将原始信号从时域转换到时频域。短时傅里叶变换公式如下
其中窗函数w(i)长度为mw，mw此处取信号长度m的八分之一，也就是m/8,s(i)是长度为m的信号，即距离-多普勒成像算法中含有微多普勒成分的距离单元。
步骤二：对短时傅里叶变换的时频图进行瑞登变换，并将变换域中的极值点标记为可能的刚体成分。瑞登变换的公式如下
这里，f(x,y)是信号短时傅里叶变换后的时频谱图，x是时间，y是频率，δ是冲激函数，α是积分投影的角度，β是投影轴上的位置。将时频图沿着直线xcosα ysinα＝β积分，相当于得到时频图沿着该直线的投影结果，根据瑞登变换的实现过程可以得到，该变换有聚焦图像中线性成分的能力，图像中的线段成分会在瑞登变换的作用下在变换域中聚焦成峰值点。而通过瑞登变换域中峰值的位置坐标(α，β)能知道，原图像中直线xcosα ysinα＝β位置上有直线成分，在本发明中，原图像为时频图，根据直线位置便能知道刚体信号的频率信息。
步骤三：使用拉普拉斯算子对图像的瑞登变换域进行处理，由此得到变换域的陡峭程度信息，图像的拉普拉斯算子表达式如下
通过对瑞登变换域进行拉普拉斯算子的处理，周围陡峭的像素点会产生极值点，标记这些点，并在下个步骤中，利用此步骤标记的峰值点排除非刚体点。寻找图像峰值点的方法分为以下两个子步骤。由于峰值点为局部的极大值点，虽然不一定要整个变换域中的最大值点，不过需要时局部的最大值，通过这些特点，对变换域中的点进行两部筛选得到峰值。一，对变换域的值进行排序，设置一个经验百分比，此处可取1％，变换域中的函数值为前1％的点才能参与峰值点的讨论。二，该点函数值大于等于其邻域的函数值。通过这两子步骤便能找到图像中的峰值。
步骤四：根据刚体信号在瑞登变换域中的性质，刚体极值产生的峰值点附近比其他非刚体分量产生的峰值点更陡峭，所以若一个瑞登变换域中的峰值在拉普拉斯算子处理过后附近还是存在峰值点的话，便说明该点为刚体信号产生的峰值点，若该峰值点在拉普拉斯算子处理过后附近不存在峰值点，那就说明该点不是刚体成分产生的峰值。由此方法能确定，哪些瑞登变换域的峰值点代表刚体分量。(此步骤中寻找峰值点的方法和上个步骤一致)
步骤五：根据上述确定的刚体成分产生的峰值在瑞登变换域中的坐标位置(α，β)确定刚体在短时傅里叶时频图中的位置参数，由于刚体在整个距离单元的积累时间内都存在，所以根据坐标便能完全确认刚体成分的全部位置信息，即在直线xcosα ysinα＝β，直线长度为整个时间轴的长度
步骤六：根据上述步骤得到的刚体成分参数设计带通滤波器，对原信号进行滤波，滤除微多普勒信号成分，最终得到纯净的刚体信号成分。
在使用距离-多普勒算法时，如果存在微多普勒成分，那么方位压缩后在雷达图像中会产生干扰的条带，在进行方位压缩前对含有微多普勒成分的距离单元进行刚体的检测提取来去除微多普勒成分能有效提高雷达图像的质量。本发明利用了刚体信号在时频域中呈现为直线的特点，利用图像直线检测工具，瑞登变换进行刚体信号的检测。但是单凭瑞登变换识别刚体信号是远远不够的，微多普勒信号局部的线性性质也会在瑞登变换域中产生峰值，所以需要进一步判断。通过观察可发现，刚体信号产生的峰值点附近陡峭程度远大于微多普勒信号产生的峰值点附近陡峭程度，所以本发明利用了拉普拉斯算子，对瑞登变换域中的峰值进行了非刚体成分的排除，即只有足够陡峭的峰值点才被认为是极值点。通过这种方式，可以在刚体信号能量较低的情况下，或者信噪比较低的情况还是能有效的检测出刚体信号成分并相应的进行滤波。
实施例一
本实施例子中使用最简单的信号仿真。信号为三个平稳刚体信号和三个非平稳信号的和信号，其中平稳信号代表刚体成分，而非平稳信号为正弦调制的信号，代表微多普勒成分。该模拟信号的时域表达式为
图2(a)为原始信号的stft，stft使用的是长度n/8的汉明窗，可以看到，刚体信号的频率不随时间变化，而且被随时间呈正弦变化的m-d成分部分的遮盖了。图2(b)为原始信号时频图的瑞登变换结果，可以看到，时频图中的直线被聚焦为了峰值点，同时，除了刚体部分，微多普勒的部分线性性质也会在瑞登变换域中产生峰值。图2(c)为瑞登变换域经过拉普拉斯算子处理的结果。
图3(a)为原信号的频谱图，图3(b)为刚体提取的结果，可以看到，本发明能够成功的提取出刚体成分，去除微多普勒成分。
实施例二
本实施例子中将针对理想的isar点散射模型进行距离-多普勒域的成像。图4为点散射模型在空间中的分布。如图4，模型为空间中拥有两个螺旋桨的飞机。表1为仿真参数，信噪比为0。
表1
图5为原始回波数据r-d算法后的图像效果，可以看到，和螺旋桨位于同个距离单元的刚体散射点被因微多普勒效应产生的干扰条带遮盖了，可以看到微多普勒效应出现在了第103距离门及其附近距离门。
图6(a)为第103距离门的stft，两个刚体散射点对应于时频图中两个平稳信号，而旋转点对应于不平稳的正弦调频信号。图6(b)为图6(a)也就是短时傅里叶频谱的瑞登变换结果。图6(c)为瑞登变换域经过拉普拉斯算子处理后的结果。
图7为最终去除微多普勒成分后的结果，可以看到，本发明能够有效的去除微多普勒成分，提高r-d算法的成像质量。
实施例三
实施例中采用来自an-26飞机的实测数据，an-26拥有两个螺旋桨，所以其回波将带有微多普勒成分。an-26飞机的三视图如图8所示，相关雷达参数可见表2。
表2
图9是不经过m-d消除处理的数据利用r-d算法的成像结果。从图像可以看出，两条干扰条带分别对应着an-26的两个螺旋桨，而这些条带遮挡了对应距离门的飞机刚体部分图像。
图10(a)为含有微多普勒成分的第150个距离门的stft图，由图可知，虽然刚体成分还是以平稳信号的形式存在，但是微多普勒成分因为两个螺旋桨之间的互相干扰以及频谱混叠等原因导致其正弦调频的形式无法辨别。图10(b)为时频图的瑞登变换域。图10(c)为瑞登变换域经过拉普拉斯算子处理的结果。
图11(a)和图11(b)分别为滤波前和滤波后的频谱图，可以看到微多普勒成分已经去除。
图12为滤波后的成像结果，和图8比较能明显看到本发明能有效的保留刚体成分，去除微多普勒分量。
需要注意的是，具体实施方式仅仅是对本发明技术方案的解释和说明，不能以此限定权利保护范围。凡根据本发明权利要求书和说明书所做的仅仅是局部改变的，仍应落入本发明的保护范围内。
1.基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于包括以下步骤：
步骤一：对原始信号进行短时傅里叶变换；
步骤二：对短时傅里叶变换的时频图进行radon变换，并将变换域中的极值点标记为可能的刚体成分；
步骤三：使用拉普拉斯算子对图像的radon变换域进行处理，得到变换域的陡峭程度信息；
步骤四：确定radon变换域中产生可能的刚体成分的峰值点；
步骤五：根据步骤四中刚体成分产生的峰值点在radon变换域中的坐标位置确定刚体在短时傅里叶时频图中的位置参数；
步骤六：根据步骤五得到的刚体成分参数设计带通滤波器，对原信号进行滤波，滤除微多普勒信号成分，得到纯净的刚体信号成分。
2.根据权利要求1所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述步骤一中短时傅里叶变换公式如下：
其中，窗函数w(i)长度为mw，s(i)是原信号，m是信号长度，w(i-m)是窗函数。
3.根据权利要求2所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述mw取信号长度m的八分之一。
4.根据权利要求1所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述步骤二中radon变换的公式如下：
其中，f(x,y)是信号短时傅里叶变换后的时频谱图，x是时间，y是频率，δ是冲激函数，α是积分投影的角度，β是投影轴上的位置。
5.根据权利要求4所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述步骤三中拉普拉斯算子表达式如下
6.根据权利要求1所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述步骤四中确定radon变换域中产生可能的刚体成分的峰值点的方法为：
步骤四一：对变换域的值进行排序，设置经验百分比为1％，取变换域中的函数值为前1％的点；
步骤四二：点函数值大于等于其邻域的函数值即为radon变换域中产生刚体成分的峰值点。
7.根据权利要求1所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述步骤四中确定radon变换域中产生刚体成分的峰值点通过判断步骤二中可能的刚体峰值点位置经过步骤三拉普拉斯算子处理后在该点还是否为峰值确定，若处理后在该点还是峰值，则此位置为刚体成分产生的峰值点。
8.根据权利要求1所述的基于radon变换与拉普拉斯算子的微多普勒抑制方法，其特征在于所述原始信号为雷达回波中所有含有微多普勒成分的距离单元信号。
技术总结