联合变化检测与子带对消技术的SAR图像干扰抑制方法

李 宁,吕宗森,郭拯危

(1.河南大学计算机与信息工程学院,河南 开封 475004;2.河南大学河南省大数据分析与处理重点实验室,河南 开封 475004;3.河南大学河南省智能技术与应用工程技术研究中心,河南 开封 475004)

0 引 言

合成孔径雷达(synthetic aperture radar,SAR)是一种主动微波遥感设备,能提供全天时、全天候、高分辨率对地观测影像,在现代遥感中具有举足轻重的地位[1-3]。但是,随着电磁环境的日益恶化,SAR经常遭受来自其他电磁设备的射频干扰(radio frequency interference,RFI)。这些干扰通常以直达波或散射波的方式进入到SAR系统中,扰乱了SAR信号的脉冲响应,加剧了SAR图像解译难度[4]。因此,有效地抑制干扰是发挥SAR应用效能的前提,具有重大研究意义。

如何有效检测和抑制各种形式的干扰,一直是SAR领域的重要研究课题。自20世纪90年代以来,多种SAR干扰抑制方法被相继提出,主要分为参数化抑制方法、非参数化抑制方法以及半参数化抑制方法[5]。参数化方法认为干扰可在时域建模,通过对干扰模型中信号的幅度、频率和相位进行估计,解算干扰模型参数,实现干扰抑制。文献[6]指出了最大似然估计算法的优越性,并将其运用于干扰抑制。文献[7]提出了一种基于参数化最大似然的干扰抑制方法。文献[8]提出了一种基于迭代自适应谱估计技术的干扰抑制方法。然而,干扰信号复杂多变,其模型不容易被精确建立,参数化方法往往出现模型失配的情况,从而影响抑制效果。不同于参数化方法,非参数化方法主要利用干扰信号与回波信号在不同变换域中的能量差异进行抑制操作。陷波法是一种经典的非参数化干扰抑制方法[9]。在距离-频域,通过陷波操作消除强窄带干扰。然而,陷波法缺失了部分的频域信号,导致最终成像结果出现异常的点目标响应。文献[10]提出两步陷波法,基于线性预测模型,对缺失的频谱进行补偿。不同于陷波类方法,干扰抑制问题亦可等效为盲信号分离问题,通过诸如特征值分解[11]、复数经验模态分解[12]、独立空间分析[13]等技术,实现干扰信号和回波信号的有效分离,达到干扰抑制的目的。然而,复杂的矩阵分解运算带来了庞大的计算负担,而且在干扰能量较弱的情况下,抑制效果往往不佳。近年来,相关学者基于干扰信号“低秩+稀疏”[14]的特性,研究出了一些半参数化干扰抑制方法。文献[15]首次从稀疏重构的理论角度出发,将干扰抑制问题转换为稀疏优化问题,实现干扰信号的提取。由于优化问题求解的计算复杂度较高,文献[16]结合干扰信号的低秩特性,在提高半参数化方法抗干扰性能的同时,大幅降低了计算复杂度。随后,基于张量理论和快速矩阵分解技术,文献[17]提出了超参数估计干扰抑制算法。然而,现有半参数化方法的算法复杂度依然较高,阻碍了其在实际工程中的大规模应用。

发展到目前,绝大多数干扰抑制方法都是在回波域实现干扰的建模、提取、抑制等操作。但是现代SAR系统,尤其是星载SAR系统,往往不提供原始回波数据产品,只提供单视复数据(single look complex,SLC)图像及以上等级的产品,导致上述算法无法直接应用。文献[18]提出了基于逆聚焦与陷波技术的后验陷波法,该方法对SAR SLC图像数据进行二维逆聚焦操作,基于陷波技术抑制窄带干扰。然而,与传统频域陷波法的缺点相同,陷波操作缺失的有用信号将使图像中部分目标的点目标响应异常,严重影响图像质量。文献[19]提出了子带对消(subband spectral cancellation,SSC)干扰抑制方法。该方法基于SAR SLC图像数据,通过距离子带图像之间的对消操作,实现干扰的提取和抑制,具有鲁棒性高、计算复杂度低、易于工程实现等优点。然而,受相干斑噪声、发射信号频率差异等客观因素的影响,即使是无干扰情况下的子带图像,其强度也不严格相等,导致SSC方法在实际应用中会损失一部分有用的图像信号,造成部分目标场景的纹理信息丢失等问题,限制了其推广应用。

针对上述问题,本文在已有SSC方法的基础上,联合变化检测技术,提出了一种的SAR图像干扰抑制方法。首先,所提方法将含干扰的SAR SLC图像变换至距离频域,进行频谱矫正、频域干扰检测和子带划分处理。然后,对消距离子带图像得到初步干扰图像,利用改进的对数比变化检测算子,实现干扰的二次检测,消除初步干扰图像中的虚警信息,得到精确的干扰位置信息,从而实现干扰的图像域精确定位与抑制。

1 干扰信号模型与频域特性

1.1 干扰信号模型

在复杂电磁环境下,SAR系统容易接收到来自同频段其他电磁设备带来的RFI[20]。干扰的存在会直接影响后续的SAR信号处理以及图像解译。经过正交解调和数字采样后,含干扰的SAR回波数据可表示为

R(τ,η)=S(τ,η)+N(τ,η)+I(τ,η)

(1)

式中：S(τ,η)、N(τ,η)、I(τ,η)分别表示有用回波信号、高斯白噪声及干扰信号;τ和η分别表示距离向快时间和方位向慢时间。

干扰信号通常可以分为窄带干扰和宽带干扰两类。如文献[5]所述,窄带干扰可看作一系列正弦信号的叠加,其模型可表示为

(2)

式中：N表示干扰信号的个数;An(η)、fn和φn分别表示第n个干扰信号的幅度、频率和初相。窄带干扰一般没有复杂的频率调制,在频谱上占用带宽较小,表现为尖峰状。相比之下,宽带干扰的频率调制较为复杂,且带宽较大。不失一般性,宽带干扰可建模为线性频率调制和正弦频率调制两种调制方式[21]。其信号模型可表示为

(3)

式中:相位项φ分为线性频率调制和正弦频率调制两种相位调制项φCM和φSM,分别表示为

φCM=2πfnτ+πKnτ2

(4)

φSM=βncos(2πfnτ+φn)

(5)

式中:Kn表示第n个干扰信号的线性调频率;βn表示正弦调制系数。

1.2 干扰频域特性

SAR系统为开放式系统,不可避免地会与一些电磁设备共享频谱,该特性使其遭受干扰信号的影响,尤其是低频段SAR(如P、L和C波段),干扰现象更为严重[4]。在大多数情况下,干扰信号来自于地面电磁设备的直达波,其干扰功率明显强于雷达回波、杂波和系统噪声。在距离时域SAR回波中,干扰信号的强度与雷达照射场景中强点目标的回波值类似,难以检测和定位干扰。然而,在距离频域SAR回波中,干扰信号具有较为固定的带宽,且其频谱能量明显强于SAR回波信号,这使得干扰信号在频域更易于检测与识别。图1为某机载L波段SAR系统原始数据在距离频域-方位时域的强度图像。

图1(a)中框出了A～F共6个含干扰的感兴趣区域,图1(b)分别放大显示了框A～F的细节。根据干扰信号带宽占SAR信号总带宽的比例划分,绿色区域A、B和C为窄带干扰,红色区域D、E和F则为宽带干扰。由图1可以看出,上述SAR回波信号受到了较为严重的复合型干扰。干扰A与干扰C频率特性相同,强度特性类似,推测为同一种窄带干扰,在方位向具有一定周期性。干扰B带宽较小,其频率具有时变性,特性较为复杂,不易分析。干扰D、E和F的带宽较大,且辐射功率较强,具有一定的能量压制性。

图1 机载L波段SAR原始数据距离频域-方位时域强度图像Fig.1 Range-frequency azimuth-time domain intensity image of L-band airborne SAR original data

2 SSC方法

SSC方法基于SAR成像的全息性、回波频谱的对称性以及RFI的频域特性,通过距离子带图像对消提取干扰,实现图像域中干扰信号的检测和消除[18]。图2给出了SSC方法的流程图,以两个子带为例给出算法的具体实现步骤,对多个子带的情况也同样适用。

图2 SSC方法流程图Fig.2 Fow chart of SSC method

由前文可知,含干扰的SAR回波数据可表示为式(1)。在不考虑噪声影响的情况下,成像后的图像信号总功率可表示为

|Sp|2=|Sr|2+|Si|2

(6)

式中:Sp，Sr，Si分别表示总图像信号、不含干扰的图像信号和含干扰图像信号。

对图像Sp进行距离向快速傅里叶变换(fast Fourier transform,FFT),并进行频谱矫正处理[22]。采用平均距离谱法[10],在距离频域对干扰进行检测,并将频谱划分为含干扰信号的距离子带,和不含干扰信号的距离子带。通过距离向逆向FFT(inverse FFT,IFFT),可得到两个子图像Sp1和Sp2,分别为

Sp1=Sr1+Si

(7)

Sp2=Sr2

(8)

式中:Sr1和Sr2分别表示两个不含干扰的子图像信号;Si表示干扰图像信号。子图像的信号功率满足:

|Sp1|2≈|Sr1|2+|Si|2

(9)

|Sp2|2=|Sr2|2

(10)

将子图像的信号功率做差运算,即将式(9)和式(10)相减,得

|Sp1|2-|Sp2|2≈|Sr2|2-|Sr1|2+|Si|2

(11)

经过频谱矫正处理后,子图像信号的幅度近似相同,即

|Sr2|2-|Sr1|2≈0

(12)

将式(12)代入式(11)可得到干扰图像信号Si为

|Sp2|2-|Sp1|2≈|Si|2

(13)

联合式(6),可以得到不含干扰的图像Sr,有

|Sr|2=|Sp|2-|Si|2=|Sp|2-||Sp2|2-|Sp1|2|

(14)

需要说明的是,为便于理解,以上推导过程以两个子带的情况为例,对于多子带的情况,上述结论同样成立。

经由上述理论推导可知,通过子图像之间的对消操作,干扰信号可以被有效提取。但是,受相干斑噪声、发射信号频率差异等客观因素的影响,子图像间的信号强度并不严格相等,也即式(12)和式(13)近似成立的原因。该近似导致部分有用信号被误判为干扰信号,提高了干扰检测的虚警率。付出的代价则是在成功检测和抑制干扰的同时,损失了部分图像的细节信息,严重的甚至造成场景中目标的丢失。

3 改进的SSC干扰抑制方法

针对上述SSC方法存在的问题,本文联合子带图像变化检测技术,提出了一种改进的SSC干扰抑制方法,具体处理流程如图3所示。相较于SSC方法,本文所提方法增加了干扰二次检测处理步骤,如图3中黑色框所示。

图3 所提方法流程图Fig.3 Flow chart of the proposed method

3.1 干扰初步检测

按照图3中SSC方法的步骤,依次进行距离向FFT、频谱矫正、干扰检测、子带划分、SSC(即差运算)操作,获得干扰初步检测图像。

3.2 干扰二次检测

初步干扰检测图像由子带图像的差分运算生成。从图像处理的角度来看,干扰初步检测即为最简单的图像变化检测过程,其检测效果容易受噪声等因素的影响,虚警率较高。

现有的SAR图像变化检测方法可分为两大类,即基于目标级别的分类比较法和基于像素级别的图像代数法[23-27]。图像分类比较法一般面向特定的目标或地物类型展开研究,需要研究对象的先验知识。图像代数法一般基于差值、比值等算术算子进行逐像素分析,提取变化信息,具有计算复杂度低、易于实现的优点[23]。由于缺少干扰在图像中的先验特征信息,针对本文所述的干扰检测工作,图像代数法具有明显的应用优势。

针对初步干扰检测图像虚警率较高的问题,本文提出了一套变化检测流程,对初步干扰检测图像中的干扰与虚警进行区分,即找出干扰对应的位置信息,消除其他因素导致的不相关变化信息,具体操作流程如下。

步骤 1图像预处理

SAR图像变化检测的预处理步骤包括几何配准、辐射校正、相干斑滤波等[23]。然而,本文所述的距离多子带图像来源于同一景主图像,图像之间不存在几何和辐射误差,无需进行几何配准和辐射校正处理。然而,子图像中固有的相干斑噪声作为一种乘性噪声,将导致图像分布不均匀,对后续检测造成不利影响。中值滤波[28]是一种简单、易于工程实现的相干斑抑制方法,能够在有效抑制相干斑噪声的同时尽可能地保留图像边缘信息。从便于工程实现和应用的角度出发,本文选取中值滤波来抑制相干斑噪声。

步骤 2图像代数运算

基于像素级的图像代数变化检测方法可以高效地逐像素提取图像之间的变化信息。对数比算子是一类常用的图像变化检测算子,具有高效鲁棒的特点,本文选取一种改进的对数比变化检测算子[24],进行子图像间的变化检测操作。

改进的对数比变化检测算子,其表达式为

(15)

式中：Sp1和Sp2分别表示中值滤波处理后的子图像;Cp表示干扰变化信息图像。相比于传统对数比算子,式(15)中常数1的增加避免了检测结果出现无穷量的情况。

步骤 3图像后处理

受残余噪声和子带图像频率差异的影响,变化图像Cp的精度有限,需要结合图像后处理技术,提升干扰位置检测精度。具体处理流程如下。

步骤 3.1采用最大类间方差(OTSU)法[29],对图像Cp进行自适应阈值分割,获得潜在的干扰位置信息,即获得干扰位置二值图，此时部分虚景信息已被消除。

步骤 3.2基于形态学开操作滤波方法[30],消除干扰位置二值图中的颗粒状斑点噪声,同时将干扰位置划分为不同的块状区域。

步骤 3.3基于块状区域的面积属性,设置合理的面积阈值,进一步消除干扰位置二值图中的虚警信息,提升干扰位置检测精度。

3.3 干扰匹配与抑制

在完成第3.2节中的步骤后,本文联合干扰初步检测图像的幅值信息和干扰变化检测图像的位置信息,获取的干扰二次检测图像,表示为

(16)

根据式(14),可得到干扰抑制后的功率图像为

(17)

需要指出的是,上述处理流程可推广应用于多个子带的情况。

4 实测数据处理及结果分析

为了验证所提方法的有效性,本文以欧空局Sentinel-1A卫星干涉宽幅模式SAR SLC图像为数据源,开展相关实验,主要雷达系统参数如表1所示。

表1 主要系统参数Table 1 Main system parameters

图4为含干扰SAR SLC距离频谱分析结果。图4(a)为距离频域-方位时域图像,可以看到干扰位于距离频域左侧并占据一定带宽。图4(b)为沿方位向求和后的距离频谱曲线,可以看出该频谱曲线对称性好,无需进行频谱矫正处理。

图4 含干扰SLC数据的距离频谱分析Fig.4 Range frequency spectrum analysis with interfered SLC data

图5为实验数据图像和干扰抑制结果。图5(a)为含干扰SAR图像,可以看到该区域地物类型丰富且地形较为复杂,涵盖城镇、湖泊、河流、农田等地物以及山地、丘陵等地貌。受干扰影响,场景中出现周期性交错状明亮条纹,并掩盖了部分真实地物及地貌信息。图5(b)～图5(d)分别为后验陷波法[18]、SSC方法[19]以及本文所提方法干扰的抑制结果。从干扰抑制效果来看,3种方法均能有效抑制图像中的干扰信号。图5(e)和图5(f)分别为SSC方法和本文所提方法提取的干扰图像。从图5(e)中可以看到,SSC方法提取的干扰图像中,除了有效的干扰信息,还包含了大量的目标地物信息,这将导致目标细节信息的丢失,降低图像的质量。如图5(f)所示,本文方法在有效去除干扰的同时,有效地保留了目标信息。

图5 实验数据及RFI抑制结果Fig.5 Experimental data and RFI suppression results

为了进一步定性分析,本文在图5(a)中选取了4景典型区域图像,包括城镇、湖泊、池塘和山丘,进行详细的对比分析,结果如图6所示。从图6中可以看到,后验陷波法在频域进行陷波处理,导致目标响应异常,造成图像模糊。传统SSC方法损失了一部分地物目标信息,导致图像的细节纹理信息没有原图像丰富。相比之下,本文方法在实现干扰二次检测的基础上,有效保留了地物目标信息,图像信息与原图基本保持一致。

图6 实验数据及RFI抑制结果细节对比Fig.6 Details comparison of experimental data and RFI suppression results

为了定量评估本文方法的性能,选取图像熵[21]和平均梯度[31]作为评价指标,对抑制结果进行分析。定量评价结果如表2所示,从表中可以看到,本文方法的平均梯度值和熵值最大,表明本文方法干扰抑制后的图像信息量最大,信息保留效果最好,进一步验证了所提方法的有效性。

表2 3种方法的性能评估Table 2 Performance evaluation of three methods

5 结 论

针对传统SSC方法的不足,本文提出了一种改进的SAR图像干扰抑制方法。基于实测星载SAR数据,开展了实验验证和定量分析。实验结果表明,相比已有方法,所提方法可以在抑制干扰的同时,有效地保留目标信息。同时,该方法具有易于工程实现,且无需雷达原始回波的优势。对于提供SLC级别数据的国内外SAR传感器,具有较高的实用价值。