CN101963662B - 基于短时分数阶傅里叶域滤波的自聚焦预处理方法 - Google Patents

Abstract

本发明涉及一种基于短时分数阶傅里叶域滤波的自聚焦预处理方法，属于合成孔径雷达(SAR)信号处理领域。所述方法的步骤如下：1)提取单个特显点的相位误差信号；2)判定该相位误差信号的相位是否以二次为主；3)在变换域滤除噪声：如果相位以二次相位误差为主，则使用短时分数阶傅里叶域滤波；否则使用短时傅里叶域滤波；4)将结果时频逆变换到时域，对相位误差进行估计和补偿，实现自聚焦处理。本方法有效解决了低信噪比环境下相位误差，特别是二次相位误差的估计和补偿问题，可以获得相位误差历程的主要部分，将其从原始信号中除去之后可以大大改善成像质量，进而提高了特显点信号的信杂比。

Description

基于短时分数阶傅里叶域滤波的自聚焦预处理方法

技术领域

本发明涉及自聚焦方法，尤其是一种自聚焦预处理方法，属于合成孔径雷达(SAR)信号处理领域。

背景技术

SAR通过距离向的脉冲压缩和方位向合成孔径处理实现二维高分辨率成像，可以全天时、全天候作业，可以对地下和水下场景成像，可以对地面沉降进行精确监控(使用干涉SAR)等，已广泛应用于战场侦查、精确制导、目标识别、资源探测、地形测绘、环境保护和灾情监测等诸多军用和民用领域。SAR的成像模型建立在雷达载荷平稳飞行的假设之上，然而，由于实际中存在各种随机扰动，雷达载荷在运动过程中出现误差是不可避免的。运动误差会造成图像分辨率下降和几何失真，甚至不能成像，因此必须对运动误差进行补偿。运动补偿(Motion Compensation，MoCo)有两种方法：一种是用雷达载荷上配备的惯性导航系统(Inertial Navigation System，INS)或惯性测量单元(InertialMeasurement Unit，IMU)佐以全球定位系统(Global Positioning System，GPS)直接测量雷达的位置和姿态等运动参数，然后根据这些参数进行补偿；另一种是用信号处理的方法从雷达回波中估计出运动误差，然后进行补偿，这种方法称为自聚焦(Autofocus)。随着成像分辨率的逐步提高，对运动补偿精度的要求也不断提高。

利用INS、IMU和GPS等设备测得的雷达运动参数无法完全满足高分辨率成像的要求，这一方面是因为目前测量仪器的精度有限，另一方面是因为除了运动误差之外还有一些误差是仪器无法测量的，如电磁波在大气中的非理想传播和雷达系统自身的噪声等，因此自聚焦是高分辨率成像中不可缺少的步骤。

现代记载SAR中，PGA(Phase Gradient Autofocus)是一种常用的自聚焦算法，这种算法有EICHEL等人在1989年提出。它估计的相位误差不是基于模型，而是巧妙地利用了相位误差的冗余性，故而能够估计任意阶相位误差，具有良好的鲁棒性，适用于多数图像场景，被广泛应用。传统的自聚焦算法中(包括PGA)，一般认为运动误差是非空变的。然而，在测绘带较宽的情况下，运动误差具有距离向的空变性；在方位波束角比较大的情况下，运动误差具有方位向的空变性。也就是说，高分辨率宽测绘SAR成像自聚焦处理必须考虑空变性，否则空变性将导致聚焦质量下降。

解决空变相位误差补偿问题的主要思路是对数据进行分块，划分所得的每一块数据中相位误差的空变性可以忽略。不过，分块处理会带来新的问题。PGA是一种基于特显点的自聚焦算法，分块越多，每块数据中可用的特显点数目就会减少，即运动误差的冗余性越低，自聚焦质量难以保证。为了解决新的问题，本发明提出一种自聚焦预处理方案，引入滤波技术，其目的在于提高SAR图像信噪比，保证较少特显点情况下能完成高质量的自聚焦。

滤波是信号处理领域中的经典问题，按照运算所在的域可将其划分为频域滤波、时域滤波、时频域滤波和空域滤波。频域滤波的效果取决于信号与噪声在频域的分离程度：平稳信号在频域占据很窄的带宽，信号与噪声高度分离，滤波效果很好；但非平稳信号在频域可能占据很大的带宽，很多噪声无法与之分离，因此滤波效果可能会很差。对非平稳信号而言，时频域滤波是更好的选择。因为非平稳信号虽然整体上占据很大的带宽(附图1a)，但是在局部却占据小得多的带宽(附图1b)，所以在时频域可以对信号和噪声进行更为有效地分离。

发明内容

为解决以上问题，实现提高SAR图像信杂比的处理目的，本发明提出基于短时分数阶傅里叶域滤波的自聚焦预处理方法。在预处理阶段使用单个特显点进行高稳健性的相位误差历程估计，获得相位误差历程的主要部分，将其从原始信号中除去之后，可以大大改善成像质量，进而提高特显点信号的信杂比。

所述方法包括如下步骤：

步骤1：从SAR图像数据中提取单个特显点的含有噪声的相位误差信号；

步骤2：判定该相位误差信号的相位是否以二次为主，判定的方法如下：

求出步骤1获得的相位误差信号的分数阶傅里叶变换最优匹配角度α，对最优匹配角度α下的分数阶傅里叶变换结果X_α(u)求模得到I(α，u)＝|X_α(u)|，然后构造如下检验统计量：

$J=\frac{S\{I(\mathrm{\α},u)>\mathrm{\η}\·\underset{\mathrm{\α},u}{\max }I(\mathrm{\α},u)\}}{S\left\{I(\mathrm{\α},u)\right\}},---\left(1\right)$

其中S{·}表示求面积运算符(是求信号与横坐标之间的面积)，系数η一般取0.6-0.8。如果η过高，易造成虚检；如果η过低，则易造成漏检。检验统计量J表示变换结果的幅值I(α，u)超过

的面积与I(α，u)面积之间所占的比例。该比例越小，则说明变换结果的聚集性越高，相位则以二次为主；反之则否。

若检验统计量J小于预设的门限J_th，认为相位以二次为主；若检验统计量J大于门限J_th，认为相位不是以二次为主；

步骤3：对含有噪声的相位误差信号进行时频变换，在变换域滤除噪声，具体方法是：

若步骤2判定相位以二次为主(即该相位误差信号判定为chirp类信号)，则使用短时分数阶傅里叶域滤波技术对含有噪声的相位误差信号进行滤波；

若步骤2判定相位不是以二次为主，则使用短时傅里叶域滤波技术对含有噪声的相位误差信号进行滤波；

进一步地，所述滤波时使用恒虚警率检测技术进行检测，方法为，根据预设的恒虚警率判定准则对相应的短时分数阶傅里叶变换或短时傅里叶变换的结果进行检测，结果中低于检测阈值的，认为是噪声，高于检测阈值的，认为是信号；然后根据检测情况，对相应时频变换结果进行噪声滤除。

也就是说，在进行滤波操作时，时频变换和CFAR是联合使用的。

步骤4：将滤波结果进行相应的时频逆变换到时域(这里进行与步骤3所进行的时频变换相应的逆变换，即短时分数阶傅里叶逆变换或短时傅里叶逆变换)，对相位误差信号进行估计和补偿，实现自聚焦处理。

进一步地，在所述步骤3中，

a)对于相位以二次为主的情况，对含有噪声的相位误差信号

进行短时分数阶傅里叶域滤波时，首先对其进行步骤2中获得的最优匹配角度α的短时分数阶傅里叶变换STFRFT_x，α(t，u)：

${\mathrm{STFRFT}}_{x,\mathrm{\α}}(t,u)={\∫}_{-\∞}^{+\∞}x\left(\mathrm{\τ}\right)g(\mathrm{\τ}-t)K_{\mathrm{\α}}(\mathrm{\τ},u)\mathrm{d\τ}---\left(2\right)$

其中g(t)为窗函数，一般取对称的实函数，K_α(t，u)是分数阶傅里叶变换核函数。当窗函数g(t)取如下形式的高斯函数时，所述短时分数阶傅里叶变换的聚集性可以达到最高：

$g\left(\mathrm{\τ}\right)={\left(\mathrm{\π}\frac{T_x\left|\sin \mathrm{\α}\right|}{B_{x,\mathrm{\α}}}\right)}^{-\frac{1}{4}}\exp (-\frac{B_{x,\mathrm{\α}}{\mathrm{\τ}}^2}{2T_x\left|\sin \mathrm{\α}\right|})---\left(3\right)$

其中，T_x和B_x，α分别是含有噪声的相位误差信号x(t)的时宽和分数阶域带宽。该分数阶傅里叶域带宽可以通过在该短时分数阶傅里叶域对信号最高频率和最低频率相减得到。

对于相位不是以二次为主的情况，则对含有噪声的相位误差信号

进行短时傅里叶域滤波。由于传统的傅里叶变换是分数阶傅里叶变换在旋转角度为α＝π/2时的特例，因此只需要将式(2)和(3)进行旋转角度α＝π/2的操作，就可以得到相对应的短时傅里叶变换结果。

b)对步骤3a的时频变换结果进行检测和滤波。为了获得信号能量在时频域的分布情况，本发明借助恒虚警率(constant false alarm rate，CFAR)检测技术。

所述恒虚警率判定准则如下：

其中σ与含噪声的相位误差信号中的噪声均值mean存在如下关系：

$\mathrm{\σ}=\sqrt{2/\mathrm{\π}}\·\mathrm{mean},$

STR(t，u)是步骤3中获得短时分数阶傅里叶变换或短时傅里叶变换结果，即时频变换结果，即时频表示(Time-frequency representation，TFR)；Th是检测阈值，由预设的虚警概率P_fa决定：

虚警概率P_fa一般设为10^-5～10^-4。

低于该检测阈值Th的，认为是噪声，高于这一阈值Th的，认为是信号，然后根据检测情况，对时频变换结果进行噪声滤除。

对比现有技术，本发明的有益效果在于提出了一种基于时频滤波的自聚焦预处理方案。首先对相位误差类型进行判断，当相位误差信号中的相位主要为二次时，即该相位误差信号为chirp信号，则采用短时分数阶傅里叶域滤波方法进行滤波，提取特显点的实际相位(步骤4)；反之则采用短时傅里叶域滤波。由于短时分数阶傅里叶变换的核函数与chirp类信号的匹配程度更高，短时分数阶傅里叶变换对chirp类信号具有更高的二维聚集性，可以在滤除噪声的同时保存更多的信号能量，提高信噪比。因此，在时频域使用短时分数阶傅里叶域滤波比传统的滤波方法可以得到更好的滤波效果，提高了自聚焦精度。

附图说明

图1非平稳信号的时频域表示示意图；其中图1(a)为频谱，图1(b)为谱图。

图2基于短时分数阶傅里叶域滤波的自聚焦预处理方法流程图。

图3多普勒相位误差历程。

图4短时分数阶傅里叶域滤波前后的单个特显点信号。其中图4(a)为原始信号，图4(b)为经过时频滤波之后。

图5运动误差补偿。其中图5(a)为相位误差估计，图5(b)为经过预处理之后的单个特显点图像

图6SNR＝-10dB时PGA相位误差估计结果。

图7SNR＝-20dB时基于短时分数阶傅里叶域滤波的相位误差估计结果。

图8实测数据处理结果；其中图8(a)为散焦图像，图8(b)为传统SPGA聚焦结果，图8(c)为使用基于时频滤波的预处理的SPGA聚焦结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合具体实例，并参照附图，对本发明实施方式进行详细说明。本发明的仿真是在matlab软件中设计实现的。

本发明提供了一种基于时频域滤波的自聚焦方法，其基本流程如图2所示。结合如下实施例，所述方法具体步骤包括：

步骤1：从SAR二维回波数据矩阵中搜索出特显点所在的最大距离门的位置，并提取单个特显点的含有噪声的相位误差信号。设该相位误差信号模型为：

其中

为多普勒相位误差历程，其形状如图3所示；μ₀＝0.3为调频率；f₀为中心频率，本实施例中令f₀＝0；a(t)为信号包络；w(t)为杂波，其概率分布因场景而异，这里用复高斯白噪声加以近似。采样序列长度为N＝512，归一化的采样间隔为

步骤2：取相位误差信号的不同旋转角度下的分数阶傅里叶变换(仿真中设定旋转角度的取值在-π～π之间)，并找出最优匹配角度α(关于最优匹配角度请参见文献：L.Qi，R.Tao，S.Y.Zhou，and Y.Wang，“Detection and parameterestimation of multicomponent LFM signal based on the fractional Fourier transform，”Science in China(Ser.F.Information Sciences)，vol.47，pp.184-198，Aug.2003.”)。本例中信号模型中已给出相位误差信号调频率参数，可以计算出相对应的最优匹配角度即分数阶傅里叶变换旋转角度为α＝-0.4878。对该旋转角度下的分数阶傅里叶变换结果求模I(α，u)＝|X_α(u)|。

利用检测统计量公式(1)判断信号能量的聚集性。式中取η＝0.5；J_th＝0.01。计算得J＝0.0018＜J_th，因此，可以判定相位以二次为主，即该特显点信号为chirp类信号。为了降低计算量，可以根据先验知识缩小阶次搜索范围，同时适当加大搜索补偿，不过此时应将检测门限J_th做相应地修改。

本发明使用了分数阶傅里叶变换。下面对这一变换进行简要介绍。信号x(t)的分数阶傅里叶变换定义为：

$X_{\mathrm{\α}}\left(u\right)={\∫}_{-\∞}^{+\∞}x\left(t\right)K_{\mathrm{\α}}(t,u)\mathrm{dt}$

$=\left\{\begin{array}{cc}{A}_{\mathrm{\α}}{\∫}_{-\∞}^{+\∞}x\left(t\right)e^{j\frac{(u^2+t^2)}{2}\cot \mathrm{\α}-\mathrm{jut}\csc \mathrm{\α}}& \mathrm{dt\α}\≠\mathrm{n\π}\\ x\left(t\right)& \mathrm{\α}=2\mathrm{n\π}\\ x(-t)& \mathrm{\α}=(2n+1)\mathrm{\π}\end{array}\right.----\left(6\right)$

其中K_α(t，u)是分数阶傅里叶变换核函数，

n为整数，α为分数阶傅里叶变换在时频平面上的旋转角度，且周期为2π，即X_α(u)＝X_α+2π(u)。

步骤3：对于相位以二次项为主的情况，则对含有噪声的相位误差信号进行短时分数阶傅里叶域滤波。滤波的目的是滤除噪声，得到相位误差信号。

根据式(2)和(3)及CFAR检测对信号进行滤波。如图4所示，可以看到在分数阶傅里叶域更好的去除了噪声的影响，保留了信号的能量。结果表明，时频滤波可以有效去除带外和带内噪声，提高信噪比。

步骤4：相位误差补偿。

将滤波后的信号进行相应的分数阶傅里叶逆变换到时域，估计出相位误差。估计结果如图5(a)所示；单个特显点的成像结果如图5(b)所示。可以看出，预处理之后所得相位误差历程估计即可满足要求。对相位误差信号进行估计和补偿的方法采用现有技术即可。

为了做进一步对比，这里给出不同信噪比环境下PGA方法和本发明提出的自聚焦方法估计误差的结果。图6为信噪比-16dB的情况下，使用PGA算法估计的相位误差，估计结果出现散焦现象。图7为信噪比-20dB情况下使用本发明所述的短时分数阶傅里叶域滤波的自聚焦方法，使用本方法之后，可以在很低的信噪比下实现对具有运动误差的图像实现自聚焦处理。

下面给出的是实测数据处理结果。实验所用数据来源于《合成孔径雷达成像——算法与实现》附带的数据光盘，数据由Radarsat International的GordonStaples提供，加拿大航天局版权所有。SAR数据参数详见下表：

参数	符号	值	单位
				采样率	Fr	32.317	MHz
脉冲宽度		30.111	MHz
				脉冲重复频率		0	MHz
距离调频率	Fr	0.72135	MHz/μs
				数据窗开始时间		6.8956	ms
脉宽	Tr	41.74	μs
				复制信号采样数		1349
每回波行采样数		9280
				雷达频率	f0	5.300	GHz
雷达波长	λ	0.05657	M
				脉冲重复频率	Fa	1256.98	Hz
有效雷达速度	Vr	7062	m/s
				方位调频率	Ka	1733	Hz/s
多普勒中心频率	fηc	-6900	Hz

实验结果如图8所示。从实验结果可以看出，用传统的PGA算法聚焦结果尚有一些散焦的痕迹，而采用本发明所述短时分数阶傅里叶滤波算法后可以更好地实现图像聚焦。

本发明包括但不局限于以上的实施例，凡是在本发明的精神和原则之下的任何局部改动和等同替换，都将视为在本发明的保护范围之内。

Claims (7)

1.基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，包括如下步骤：

步骤1：从SAR图像数据中提取单个特显点的含有噪声的相位误差信号；

步骤2：判定该相位误差信号的相位是否以二次为主，判定的方法如下：

求出步骤1获得的相位误差信号的分数阶傅里叶变换最优匹配角度α，对最优匹配角度α下的分数阶傅里叶变换结果X_α(u)求模得到I(α，u)＝|X_α(u)|，然后构造如下检验统计量：

$J=\frac{S\{I(\mathrm{\α},u)>\mathrm{\η}\·\underset{\mathrm{\α},u}{\max }I(\mathrm{\α},u)\}}{S\left\{I(\mathrm{\α},u)\right\}},$

其中运算符S{·}表示求面积；η表示系数；

若检验统计量J小于预设的门限J_th，则认为相位以二次为主；

若检验统计量J大于预设的门限J_th，则认为相位不是以二次为主；

步骤3：对含有噪声的相位误差信号进行时频变换，在变换域滤除噪声，具体方法是：

若步骤2判定相位以二次为主，则使用短时分数阶傅里叶域滤波技术对含有噪声的相位误差信号进行滤波；

若步骤2判定相位不是以二次为主，则使用短时傅里叶域滤波技术对含有噪声的相位误差信号进行滤波；

步骤4：将滤波结果进行与步骤3所进行的时频变换相应的逆变换到时域，对相位误差信号进行估计和补偿，实现自聚焦处理。

2.根据权利要求1所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，步骤2中系数η的取值范围为0.6-0.8。

3.根据权利要求1所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，步骤3所述滤波时使用恒虚警率检测技术进行检测，方法为，根据预设的恒虚警率判定准则对相应的短时分数阶傅里叶变换或短时傅里叶变换的结果进行检测，结果中低于检测阈值的，认为是噪声，高于检测阈值的，认为是信号；然后根据检测情况，对相应时频变换结果进行噪声滤除。

4.根据权利要求1所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，步骤3中，

对于相位以二次为主的情况，对含有噪声的相位误差信号进行短时分数阶傅里叶域滤波时，首先对其进行步骤2中获得的最优匹配角度α的短时分数阶傅里叶变换STFRFT_x，α(t，u)：

${\mathrm{STFRFT}}_{x,\mathrm{\α}}(t,u)={\∫}_{-\∞}^{+\∞}x\left(\mathrm{\τ}\right)g(\mathrm{\τ}-t)K_{\mathrm{\α}}(\mathrm{\τ},u)\mathrm{d\τ}$

其中g(t)为窗函数；K_α(t，u)是分数阶傅里叶变换核函数；

对于相位不是以二次为主的情况，对含有噪声的相位误差信号

进行短时傅里叶域滤波时，进行的旋转角度α＝π/2的短时分数阶傅里叶变换，即可得到相应的短时傅里叶变换结果。

5.根据权利要求4所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，当窗函数g(t)取如下形式的高斯函数时，所述短时分数阶傅里叶变换的聚集性达到最高：

$g\left(\mathrm{\τ}\right)={\left(\mathrm{\π}\frac{T_x\left|\sin \mathrm{\α}\right|}{B_{x,\mathrm{\α}}}\right)}^{-\frac{1}{4}}\exp (-\frac{B_{x,\mathrm{\α}}{\mathrm{\τ}}^2}{{2T}_x\left|\sin \mathrm{\α}\right|})$

其中，T_x和B_x，α分别是含有噪声的相位误差信号x(t)的时宽和分数阶域带宽；该分数阶傅里叶域带宽可以通过在该短时分数阶傅里叶域对信号最高频率和最低频率相减得到。

6.根据权利要求3所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，所述恒虚警率判定准则如下：

其中σ与含噪声的相位误差信号中的噪声均值mean存在如下关系：

$\mathrm{\σ}=\sqrt{2/\mathrm{\π}}\·\mathrm{mean},$

STR(t，u)是步骤3中获得短时分数阶傅里叶变换或短时傅里叶变换结果，即时频变换结果；Th是检测阈值，由预设的虚警概率P_fa决定：

7.根据权利要求6所述基于短时分数阶傅里叶域滤波的自聚焦预处理方法，其特征在于，虚警概率P_fa设为10^-5～10^-4。

Images