作为一种主动微波遥感探测手段，合成孔径雷达(synthetic aperture radar,SAR)具有全天时、全天候和远距离高分辨成像等光学遥感设备无法替代的能力。极化和干涉SAR系统在侦查监视、地质测绘、资源勘查和灾害监测评估等领域^[1]获得了广泛应用。

SAR图像质量直接影响其应用效果。RD(range Doppler)、CS(chirp scaling)和ωK等基于加窗Fourier变换的成像算法^[2]通常用于SAR图像重建，对SAR图像的应用往往以此为基础。而在城市居民区、港口和海面舰船等目标区域内常常形成如二面角和三面角等形状构成的孤立强散射点，其旁瓣强度远高于周围像素(见图 1)，降低了图像质量，影响了SAR图像的应用效果。空间变迹法^{[3, 4]}和自适应旁瓣抑制方法^[5]从窗函数演变而来，以牺牲图像分辨率为代价，在较大程度上抑制了图像旁瓣，但这种非线性变换方法可能对图像统计模型和极化散射特征等构成影响^[6]，不利于后续应用。基于小波理论的波形设计的旁瓣抑制方法^[7]需要改变发射信号波形，不易于实现。

图 1 TerraSAR-X单极化SAR图像

图选项

本文提出了一种基于二维外推陷波滤波器的孤立强散射抑制方法，利用等波纹有限冲击响应(finite impulse response,FIR)陷波滤波器实现强散射旁瓣抑制，同时基于一阶自回归(auto regressive,AR)模型和稳态Kalman预测理论推导了运算简单的外推算法，利用外推序列进行滤波，消除了直接卷积滤波引入的滤波误差，避免了有效数据长度缩短问题。该方法在不降低图像分辨率、不改变图像周围像素统计和极化散射特性的情况下，有效地抑制了孤立强散射旁瓣。通过仿真和实测SAR图像处理，该方法的有效性和处理性能得到了验证。

针对孤立强散射点，本文设计了一个二维等波纹线性相位陷波FIR滤波器对相位历史域回波数据进行滤波，以减弱其散射强度。等波纹滤波器具有有限通道起伏的特性，可以将图像幅度误差限制在给定范围内，有限的过渡带范围可以将强散射能量控制在给定的区域内^[8]。

二维等波纹FIR滤波器的设计方法主要包括采样法、加窗法和变换法。前2种方法难以满足复杂的设计指标。本文是在每个方向上，利用PM (Parks-McClellan)法^{[9, 10]}设计出Chebyshev意义下的最优一维等波纹FIR滤波器，再利用变换法设计出各向同性的二维等波纹线性相位FIR滤波器。

以表 1中的参数为例，本文利用PM法求解了满足滤波器参数的最优滤波器长度，同时利用PM法设计了一维等波纹线性相位FIR陷波滤波器，其幅度和相位响应如图 2所示，利用变换法设计出的各向同性二维滤波器的幅度响应如图 3所示。

图 2 一维陷波滤波器幅相响应

图选项

图 3 二维陷波滤波器幅度响应

图选项

线性FIR滤波通常通过卷积来快速计算。由于时域卷积等价于频域乘积，因此滤波器阻带只会衰减信号频域对应于阻带的部分，而对于散布在阻带以外的强散射旁瓣无能为力，因此并未达到阻带信号能量按照设计参数完全滤除的目的。这种非理想的频率特性是由于卷积滤波会引入大量无意义的零元素代替原始信号参与滤波运算造成的，为了改变这一现状，本文结合文^[11]提出的基于AR模型辨识和稳态Kalman预测的最优外推滤波策略进行滤波处理，以实现滤波器对阻带内信号旁瓣能量的滤除。

由于AR模型辨识需要模型定阶和参数估计等复杂运算，对于长信号来说，模型阶数高，运算复杂度高，处理精度低。为此，本文提出了一种基于Chirp-Z变换的AR(1)模型辨识，该算法计算效率高，计算精度高，数值稳定性好。

陷波滤波器阻带内通常只有一个(或者是有限个)强散射信号成分，滤波的目的在于滤除强散射信号成分，因此外推的关键在于对这个强散射信号成分的外推，那么AR建模的关键也是对这个强散射信号成分的建模。根据AR建模复指数信号的合理性分析可知，1个强散射信号成分对应着距离维和方位维这2个AR模型的特征根λ₀和γ₀:

其中： n和m分别为距离和方位采样时间，u₀和v₀为强散射信号成分的距离维和方位维角频率。

因此，AR模型辨识问题就转化为强散射信号角频率估计问题，而这个问题可以通过插值Fourier变换来解决，为了提高估计精度和运算效率，利用Chirp-Z变换代替Fourier变换在频域插值以提高角频率估计精度:

其中： x_n,m为待外推二维信号，Ω(u₀,v₀)表示强散射点所在的邻域。在距离维和方位维形成2个AR(1):

其中： x⁽¹⁾_n 和x⁽²⁾_m 分别为距离维和方位维信号，ε⁽¹⁾_n 和ε⁽²⁾_m 为对应的驱动噪声。

结合基于高阶AR模型的稳态Kalman预测公式^{[11, 12]}，推导AR(1)模型稳态Kalman预测外推公式可得：

其中P为外推步长。通常选择前向外推长度、后向外推长度、外推步长和滤波器阶数一半这4者相等，从而有P=(L－1)/2。执行滤波过程中，将外推后的数据通过滤波器，考虑到线性相位滤波器的固定延时特性，忽略滤波器输出序列的前(L－1)个数据，所得序列与外推前的输入序列等长，且输出采样均为有效数据，消除了直接卷积滤波处理的误差，客观上保持了信号频谱的分辨率。

陷波滤波对强散射信号进行抑制，在不改变图像性质的前提下，达到了削弱强散射旁瓣的目的。但滤波器阻带和过渡带内的图像因滤波引入误差，甚至通带内的图像因通带起伏而引入误差，为了补偿通带和过渡带内图像的幅度(相位未引入误差)，利用等波纹滤波器通带和过渡带幅度响应反比于滤波后图像的通带和过渡带内像素幅度；为了修正阻带强散射幅度和相位误差，可以用滤波前阻带强散射幅度和相位进行替换。

强散射点抑制处理是对SAR相位历史域回波信号进行的处理，对于单视复(single-look complex,SLC)图像，需要首先变换到相位历史域，并补充窗函数。孤立强散射点抑制处理步骤为：

1) 根据幅度检测一个孤立强散射点；

2) 利用式(2)估计强散射点的角频率；

3) 利用式(4)对二维回波进行前后向外推；

4) 设计等波纹FIR陷波滤波器，根据强散射点角频率对滤波器预调制，将陷波窗口对准强散射点位置；

5) 利用FFT执行二维快速滤波，截取信号的有效部分；

6) 利用加窗Fourier变换获得滤波后SLC图像；

7) 根据滤波器设计参数对滤波后图像进行补偿处理；

8) 迭代步骤1—7，可以抑制多个孤立强散射点旁瓣，直到所有孤立强散射点旁瓣得到处理。

为了验证算法的有效性，本文选择仿真二维点散射响应信号和实测SAR相位历史域回波信号进行强散射点旁瓣抑制处理。

将大量二维谐波信号叠加来仿真二维点散射响应。首先随机产生K个范围限制在[1,1 024× 1 024]内的相异整数C_k(k=1,2,…,K)，定义非负整数p_k和q_k为

$\left\{ \begin{array}{l} {x_{n,m}} = s_{n,m}^{(1)} + s_{n,m}^{(2)},\\ s_{n,m}^{(1)} = \sum\limits_{k = 1}^K {\exp \left[ {j2\pi } \right.\left( {\frac{{n{p_k}}}{{4N}} + \frac{{m{q_k}}}{{4M}}} \right)} + \left. {J{\phi _k}} \right],\\ s_{n,m}^{(2)} = 1000\exp \left[ {j2\pi } \right.\left( {\frac{{118.5n}}{N} + \frac{{95.5m}}{M}} \right) + j\left. \phi \right]. \end{array} \right.$

其中·和mod(·)分别为取整和取余操作。二维谐波信号可得：

$\left\{ \begin{array}{l} {p_k} = \left\lfloor {\frac{{{C_k}}}{{1024}}} \right\rfloor ,\\ {q_k} = \bmod ({C_k},1024). \end{array} \right.$

其中: φ_k和φ服从0到2π之间的独立均匀分布；且n=0,1,…,N－1，m=0,1,…,M－1。本文取K=8 192，N=M=256。

仿真二维谐波信号见图 4a，显然强散射点的旁瓣严重影响了周围信号的幅度。二维等波纹FIR陷 波滤波器设计参数见表 1，经过预调制的滤波器

图 4 仿真二维强散射点旁瓣抑制效果

图选项

幅度响应如图 4c所示。由图 4a和4d的强散射点抑制前后图像对比可见，强散射点旁瓣得到很大程度上的抑制，图像分辨率得以完全保持。

实测SAR图像分析选择中科院电子所C波段VV极化SLC SAR图像进行处理，处理前图像如图 5a所示。二维等波纹陷波滤波器的设计参数为阻带截止频率0.002π rad，通带截止频率0.01π rad，通带和阻带起伏上限均为0.043 dB，设计滤波器长度为487。强散射旁瓣抑制效果见图 5b，旁瓣抑制效果明显，而图像分辨率未受影响。

图 5 VV极化SLC图像强散射旁瓣抑制效果

图选项

本文提出一种强散射点旁瓣抑制处理方法，参数估计性能可靠，计算处理简单，能胜任宽幅SAR图像处理。强散射点旁瓣抑制效果明显，且处理过程不会影响周围图像，为后续极化SAR图像应用提供了方便。与空间变迹法和自适应旁瓣抑制方法相比，该方法可以适应任意信号采样率，且处理过程简单，不会引起周围像素性质的改变。下一步，将扩展该方法到任意斜视SAR图像强散射点旁瓣抑制处理中。