康 健 童风雨 白雨松 丁 翔 冀腾宇 张 柘

①(苏州大学电子信息学院 苏州 215006)

②(西北工业大学数学与统计学院 西安 710072)

③(苏州市微波成像处理与应用技术重点实验室 苏州 215123)

④(苏州空天信息研究院 苏州 215123)

⑤(中国科学院空天信息创新研究院 北京 100190)

1 引言

合成孔径雷达(Synthetic Aperture Radar,SAR)通过主动发射、接收并处理电磁波信号，成为具有全天时、全天候对地观测能力的遥感技术，在国防、自然资源等领域发挥着重要作用。近年来，新的卫星及信号处理技术的飞速发展极大地促进了SAR对地观测技术的进步，使得同时满足高分辨率、大幅宽以及高重访频率的SAR图像获取成为可能，从而提升了其在高精度、大范围地表动态监测等方面的能力[1–3]。通过将不同时间获取的SAR图像进行配准，得到的SAR图像时间序列能提供被观测地区在时间维度的变化信息，进而能更好地服务于农业监测、灾后评估以及城市规划等方面[4]。

SAR成像过程中的相干处理特性使得图像中像素幅度值呈现出随机扰动的现象，这主要是雷达波与像素点内各个散射点相互作用产生的回波相参叠加所导致的。与常见的加性高斯噪声不同，这种相干斑噪声是一种乘性噪声，其非平稳特性对SAR图像的后续解译构成了严重挑战[5,6]。一直以来，SAR图像的相干斑抑制是SAR图像处理的重要任务之一，早期的去噪方法采用简单的低通滤波器对单视SAR图像进行滤波，在噪声被抑制的同时，图像细节信息也存在一定程度的丢失。为了解决此问题，文献[7,8]均利用空间域自适应滤波的方法，在降噪的同时，可以充分保留图像的细节特征。文献[9]引入了全变分(Total Variation,TV)正则化项，且利用对数变换将噪声的乘性模型转化为满足加性模型的数据拟合项，进而提出了相干斑抑制的优化模型并利用最优化方法得到去噪结果。近几年，随着深度学习方法在图像处理领域的广泛应用，一些工作也将其用来抑制SAR图像相干斑，并取得了良好的效果[10–12]。其他单通道SAR图像抑制方法详见综述类文献[6,13,14]。

近年来，世界各国均着眼于SAR卫星星座的建设，如欧空局的哨兵1号(Sentinel-1)卫星星座、中国的高分3号(Gaofen-3)卫星星座，使得具有全球尺度、高重访频率的SAR图像处理成为领域内研究的重点。空间域配准得到的具有高时间分辨率的SAR图像时间序列能充分地反映出观测时间段内地物变化信息。针对时序SAR图像的相干斑抑制，国内外学者也做了很多相关工作。文献[15–17]将单通道SAR图像搜索同质像素点或者图像块进行去噪的思想拓展到时序SAR图像中，即在时空方向上均进行同质像素点或图像块的搜索，进而再进行滤波。文献[18]提出的多基线干涉SAR技术(SqueeSAR)能自适应地将同质分布式散射点进行选取并联合滤波，进而同时实现对多时相SAR幅度及相位进行去噪。文献[19]提出的多时相SAR图像三维块匹配(Multitemporal SAR-BM3D,MSAR-BM3D)方法利用时间维度上的冗余信息以及变换域协同滤波将单通道SAR图像三维块匹配(SAR Block Matching 3D algorithm,SAR-BM3D)方法[20]进行拓展，使其能实现对时序SAR图像的相干斑抑制效果。文献[21]提出的基于比值的多时相SAR图像去噪(Ratio-based Multitemporal SAR,RABASAR)方法通过利用SAR时间序列图像与其在时间维度平均得到的均值图像计算对数比，使SAR图像相干斑的空间非平稳性能被很好地抑制，从而能进一步提升时序SAR图像的去噪效果。

然而，除了空间域SAR图像中原有的相干斑噪声带来的挑战之外，时间序列中存在的变化信息，如出现人造目标、建筑物的结构变化等，均能使得观测区域在时间维度上出现异常的散射强度，如图1所示。由于这些沿时间维度突变信号来源于观测区域地物变化，特别是人造目标，现有的大部分针对SAR图像时间序列进行相干斑抑制的方法并没有对此信号成分进行建模，从而导致去噪结果中含有人为引入的噪点或者出现“虚假散射点”现象。

图1 含有沿时间维度突变信号的SAR幅度时间序列Fig.1 SAR amplitude time series with outliers along the temporal dimension

另外，在一些应用领域，如目标检测等，这些沿时间维度突变信号通常反映了一些人造目标在观测时间内的变化情况，也能为后续的解译工作提供有价值的信息。为此，Baier等人[22]提出了一种非局部低秩信号分解的方法对时序SAR图像进行相干斑抑制(Robust Nonlocal Low-Rank SAR Time Series Despeckling Considering Speckle Correlation by Total Variation Regularization,DespecKS-NLLRTV)，分别将观测数据分解为低秩信号成分、时间相关的相干成分以及稀疏的沿时间维度突变信号成分，在相干斑抑制的同时也能降低突变信号对去噪结果产生的影响。虽然此方法在实际数据中的效果明显，但是由于其加性信号分解是在原始数据域中进行的，并没有将信号成分与相干斑成分进行充分解耦，而且相干斑成分与突变信号成分均被描述为稀疏项并由 L1范数进行建模，也会导致相干斑成分与突变信号成分存在耦合现象，从而导致信号分解结果不准确，并不能真实地反映观测时间内的地物变化所导致的散射强度突变的情况。此外，DespecKS-NLLRTV在对TV邻近算子进行求解时采用了Chambolle-Pock算法[23]，其运算效率不高，进一步限制了该方法在实际数据中的应用。

为了对时序SAR图像进行相干斑抑制的同时，准确地提取出观测时间段内产生的沿时间维度突变散射点，本文提出了一种基于对数域加性信号分解的方法，首先，通过对对数域SAR强度值进行统计建模，得到时序SAR数据的对数似然，再引入针对不同成分的先验知识，提出相应的正则项，进而得到各个成分的后验概率，通过最大后验概率(Maximum A Posterior,MAP)的估计方法，再得出问题对应的优化模型，最后利用交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)[24]进行求解优化。所提方法在仿真以及实际数据中均验证了效果，通过对比试验分析了所提方法在去噪性能以及突变信号成分分离效果上的提升。

2 问题建模与求解方法

本文提出的方法总体框架如图2所示，为了保持时序SAR图像中永久散射点(Persistent Scatterer,PS)的散射特性，首先计算振幅离差指数(Amplitude Dispersion Index,ADI)[25]剔除可能的永久散射点，对于每个非永久散射点，再通过Kolmogorov-Smirnov (K-S)检测[18]将非局部搜索窗口中的同质像素点进行筛选，进而将所有同质点与其排列并构成待分解矩阵，应用所提出的对数域信号分解方法得到相干斑抑制成分以及沿时间维度突变信号成分，最后采用聚合求均值策略得到最后的去噪以及突变信号结果。

图2 本文所提方法流程图Fig.2 Flow chart of the proposed method

2.1 问题建模

SAR图像的强度值服从以下Gamma分布[9,26]：

其中，I ∈R+表示像素强度值，R ∈R+为待恢复的信号强度值，L>0为视数个数，Γ(·)为Gamma函数。在乘性噪声模型下，强度I表示为

其中，n为均值，为E(n)=1的独立同分布的随机变量，进而得到强度的均值与方差为

从式(4)可以看出，噪声在空间方向上是非平稳的，并与信号真实强度值有关。为了方便处理，通常采用对数函数将乘性噪声模型进行转化，得到满足加性噪声模型的Fisher-Tippett分布[26]：

其中，g=lg(I),z=lg(R)。根据上述对数域中的强度似然模型，在MAP框架下，待恢复的信号可以通过求解以下模型来得到：

其中，lgpg(g|z)为对数域中的信号对数似然，R(z)为待恢复信号的先验知识对应的正则化项，λ为正则化项参数。对于时序SAR图像，空间像素点中的时间序列对数域强度值服从以下对数似然：

其中，t为时间索引，C 为常量，在对变量zt进行优化时，可以将其忽略，并不影响优化结果。

上述MAP问题需要对对数域时间序列强度值的先验知识进行建模，本文主要考虑以下先验：

(1) 非局部自相似性：在一定观测范围内，存在地物目标尺寸远大于SAR图像分辨率或者不同目标属于同种地物类别等情况，使得区域内的目标散射特性有相似性，进而可以利用非局部搜索方法将一些同质像素点对应的时间序列找出并进行联合处理。为此，本文利用K-S检测方法将一定搜索范围内(例如，21×21窗口范围)的同质点进行搜索选取，再将候选时间序列构成观测矩阵。

(2) 平滑性：由于观测矩阵是由同质像素点对应的时间序列构成的，其在时间以及空间维度上(表示为矩阵的行和列)具有平滑特性，即散射强度变化不大，本文采用TV范数来对此进行建模。

(3) 稀疏性：对于时序信号中存在的沿时间维度突变信号，其本身具有稀疏特性，本文利用 L1范数对其进行刻画。

根据上述似然以及正则化项，本文提出了如下模型对时序SAR图像进行相干斑抑制并对观测数据中的稳定信号及突变信号成分进行分离：

其中，元素gti所构成的矩阵是将K-S检测方法选取到的同质像素点对应的时序列向量进行按行排列得到的，其时间维度及空间维度的索引分别为t,i。式(8)中X为相干斑抑制后的信号成分，S为沿时间维度突变信号成分，α,β分别为不同正则化项对应的参数。给出矩阵A，其TV范数∥A∥TV定义为

其中，i1,i2分别为矩阵列和行的索引。

2.2 求解方法

为了有效求解所提模型式(8)，本文首先利用变量替换及分解方法[9]将上述模型进行转化，得到：

其中，D(·):=[Dt(·);Di(·)]表示二维1阶导数算子。再将有约束的优化模型转化为无约束优化模型，得到增广拉格朗日函数：

其中，〈·〉表示为内积算子，∥·∥F为矩阵F-范数，T1,T2为优化引入的辅助变量，µ为惩罚项对应的参数。进而利用ADMM方法对其进行优化求解，求解过程如下。

2.2.1Z子问题

类似于文献[26]，上述问题可以通过牛顿法(Newton’s method)进行逐像素多次迭代求解：

2.2.2X子问题

可以通过计算式(14)对于变量X的梯度，并且令梯度为0，得到线性方程：

其中，D∗(·)表示D(·)的伴随算子，利用频域快速求解方法[27]得到式(14)的解。

2.2.3Y,S 子问题

式(16)与式(17)的解可以由软阈值(Soft Thresholding)[28]算子得出。

2.2.4 乘子T1,2更新

按照式(18)，式(19)子问题求解规则，在迭代一定次数并满足收敛条件后，得到的与为当前像素点及其同质点时间序列的相干斑抑制矩阵和突变信号矩阵。

由于上述分解过程是在对数域中进行的，需要将这两种成分转换到原始信号域，令为去噪后的SAR图像像素强度值，根据所提方法的分解模型可得到：

则在原始信号域中，exp()为去噪后的稳定信号成分，exp()(exp()-1)为沿时间维度突变信号成分。

2.3 永久散射点剔除

2.2 节所提方法是针对每个非永久散射点进行的，由于永久散射点几乎不受相干斑噪声的影响，这些散射点的时序散射信息需要被保留。本文采用ADI方法对其进行剔除，通过计算SAR时间序列幅度值的均值与标准差之间的比值，即mA/σA，并将大于一定阈值的像素点进行剔除。

对于每个非永久散射点，所提方法均会找到相对应的同质点时序信号，并对所构成的矩阵进行分解，在遍历所有像素点后，再利用求平均的策略将每个像素点对应的去噪成分以及沿时间维度突变信号成分进行聚合，从而得到最终的相干斑抑制结果以及突变信号分离结果。

3 仿真与实测数据实验

3.1 仿真数据实验结果

本文采用的仿真数据与文献[22]一致，其幅度图像与其随时间变化如图3所示。产生均值为1的高斯随机变量以及值为0.3的指数相参损失用来模拟SAR图像中的相干斑效应。共生成24幅图像作为仿真的时间序列。此外，分别随机选取了0.5%以及1.0%比例的SAR像素，用服从Rice分布的随机数将其替换，用来模拟沿时间维度突变信号干扰现象。利用K-S检测选择同质点的阈值设置为0.25，正则化项参数α,β分别选取为5.0,10.0。为了验证所提方法的有效性，本文选择了以下4种时序SAR图像去噪方法进行对比：(1) SqueeSAR[18];(2) MSAR-BM3D[19];(3) RABASAR[21];(4) DespecKS-NLLRTV[22]，并采用峰值信噪比(Peak Signal-to-Noise Ratio,PSNR)、平均结构相似性指标(Mean Structure Similarity Index Measure,MSSIM)以及平均互相关(Mean Cross Correlation,MCC)来定量验证方法在相干斑抑制上的效果。

图3 本文采用的仿真数据Fig.3 Simulation data used in this paper

图4(a)—图4(f)、图4(g)—图4(l)分别展示了所有对比方法在没有沿时间维度突变信号以及在0.5%像素的突变信号干扰下得到的相干斑抑制结果。相应的定量结果如表1所示。在没有突变信号的情况下，现有的方法如MSAR-BM3D和RABASAR能取得比DespecKS-NLLRTV和本文方法更好的相干斑抑制效果。其中一个原因是MSAR-BM3D和RABASAR均利用了相似图像块对噪声进行抑制，这种方法在例子中给出的同质区域较多的场景中能更好地在空间维度对信号进行平滑处理，然而，在MSSIM指标下，所提出的方法能取得比MSARBM3D更好的效果，主要原因是基于同质像素点的滤波方法能更好地保持场景的空间结构信息。当突变信号的数量开始增加时，DespecKS-NLLRTV和本文方法均能鲁棒地对相干斑噪声进行抑制，其恢复效果很大程度上并不受突变信号的数量影响。由于其他方法并没有考虑时序数据中存在的沿时间维度突变信号，其图像恢复效果随着突变信号的增多而变差。与本文所提出的对数域信号分解方法相比，直接在原始信号域进行信号分解的方法DespecKS-NLLRTV在PSNR等指标上有大约3 dB的差距。由于原始信号域中的信号与相干斑噪声成分存在耦合，直接对观测得到的数据进行信号分解并不能消除耦合产生的影响，使得信号的恢复结果并不理想。此外，对于非局部自相似性构建的时序矩阵，DespecKS-NLLRTV引入了低秩先验进行约束，这会导致恢复得到的时序信号能量有所丢失，也会使得PSNR评价指标降低。与此相比，本文所提出的方法充分考虑到了相干斑噪声的统计特性，并在对数域中将信号成分与其他成分进行分解，能在噪声抑制的同时，很好地保持恢复信号的原始结构。由于本文方法并没有采用低秩约束，并不会导致恢复出的时序信号能量丢失，并且可以大幅提升算法的计算效率。值得注意的是，在沿时间维度突变信号数量较少或不存在的情况下，所提出算法的PSNR值反而较低，而随着突变信号数量增加，PSNR有所增加。其原因在于，沿时间维度突变信号的数量增加使得观测信号的本质特性更加符合所提方法的信号分解模型，从而使得分解得到的稳定信号成分更加接近于真值。如果不存在突变信号，仍然会有一部分信号成分被分解成为稀疏信号项，反而会影响相干斑抑制得到的信号成分恢复精度。

表1 仿真数据相干斑抑制结果的定量分析Tab.1 Quantitative analysis of speckle suppression results in simulation data

图4 仿真数据相干斑抑制结果Fig.4 Speckle suppression results of the simulated data

此外，为了测试提出方法对于不同参数设置下的敏感程度，以PSNR值为例，图5给出了在不同α以及β值下的时序图像恢复效果。可以看出，当α>1,β >2时，所提方法的恢复效果受参数值变化敏感度不高，而且当突变信号增多时，给予 L1范数项更大的权重可以使得恢复效果有进一步提升。为了进一步验证所提方法受图像数目的影响，分别从原始的24张图像中选取3,5,10,15,20张图像构成新的时序数据，用所提方法对其进行相干斑抑制，相应的超参数设置保持不变，即α=5.0,β=10.0，进而计算得到结果的PSNR值如图6所示。可以看出，随着图像数目的增加，所提方法的相干斑抑制结果也随之提升。如果图像数目过少，会影响K-S检测方法对于同质点选取的精度，进而使得所构成的时序矩阵自相似性较差，导致所提方法的去噪效果降低。因此，所提出的方法适用于图像数目较多(多于20张)的时序SAR数据去噪。此外，尽管实验中突变信号的数目增加，所提方法并不受到突变信号数目不同而产生效果上的影响，从而进一步验证了所提方法的鲁棒性。

图5 本文方法的参数敏感性分析Fig.5 Parameter sensitivity analyses of the proposed method

3.2 仿实测数据实验结果

本文分别选取了上海近海以及浦东机场区域的垂直极化干涉宽幅哨兵1号(Sentinel-1 IW VV)数据，两组时间序列获取时间段均为2019-12-02到2021-02-06，每组包含36张图像。其距离向及方位向分辨率分别约为2.33 m和13.94 m。所提方法的参数设置与3.1节仿真数据一致。

图7展示了所有对比方法分别在2020-09-15以及2021-02-06获取的SAR图像上的相干斑抑制结果。可以看出，近海区域的舰船产生的强散射点构成了时间序列中的沿时间维度突变信号。传统的方法如SqueeSAR对于舰船区域的去噪效果差，其平均结果导致了此区域产生模糊的亮斑，对于MSARBM3D，可以发现舰船附近区域的去噪结果存在人为引入的噪声，使得其并不能真实反映出海面散射特性。尽管RABASAR能对大部分同质区域进行很好的相干斑抑制，但对于舰船目标，其散射结果易出现模糊现象，从而损失了真实的船体结构信息。相比于上述几种方法，采用信号分解思想的DespecKS-NLLRTV和本文方法均能通过沿时间维度突变信号分离的策略降低其对于相干斑抑制的影响。相比于DespecKS-NLLRTV，本文方法更好地利用了SAR图像的统计特性，减小了直接在原始图像域中做分解存在的信号与相干斑噪声的耦合效应。另外，由于DespecKS-NLLRTV中的模型对于突变信号以及相干斑成分均采用L1范数进行建模，会导致两种成分分离结果互相干扰。相比之下，本文方法能更加准确地对突变信号进行提取。为定量评估不同方法的相干斑抑制结果，我们分别在两次观测场景中选取了同质区域A1和A2，并根据所有方法得到的去噪结果计算相应区域内的等效视数(Equivalent Number of Looks,ENL)。从表2可以看出，所提方法可以实现良好的去噪效果，可以对于同质区域内的噪声进行有效抑制。虽然MSAR-BM3D在这些同质区域内取得了最好的等效视数，但其在散射特性复杂的区域并不能得到很好的去噪结果，尤其在具有沿时间维度突变散射点附近，其去噪结果通常会引入人为噪声干扰。与上述分析一致，相比于原始信号域信号分解方法DespecKS-NLLRTV，本文提出的信号分解策略能在同质区域内获得更高的ENL。此外，针对得到的突变信号图像，本文进一步计算其香农熵，从而定量地反映了图像的复杂程度。表3列出了两种方法得到的突变信号图像所计算的香农熵，可以看出本文方法提取得到的突变信号图像质量较高，而DespecKS-NLLRTV并不能很好地将相干斑噪声与突变信号成分进行有效的分离，使得提取得到的突变信号中含有噪声干扰，从而导致了其图像熵值较高。

表2 4块同质区域计算得到的等效视数Tab.2 The equivalent apparent number calculated by four homogeneous regions

表3 沿时间维度突变信号的熵值分析Tab.3 Entropy analysis of the outliers along the temporal dimension

图7 上海近海区域哨兵1号数据的相干斑抑制结果Fig.7 Speckle suppression results for Sentinel-1 data in the Shanghai offshore region

图8给出了所有方法在上海浦东机场区域2019-12-02以及2020-09-27获取到的SAR图像相干斑抑制结果。场景中由于飞机目标的出现(放大区域)，其强后向散射特性构成了机场区域内的SAR像素强度突变，使其成为时间序列中的沿时间维度突变信号。在这些区域内，SqueeSAR方法中的简单空间平均模糊掉了强散射目标，从而使得去噪结果中含有虚假亮斑。尽管MSAR-BM3D和RABASAR能对大部分区域进行有效的相干斑抑制，但由于这些方法没有特别考虑时序图像中的沿时间维度突变信号，使得去噪结果中包含人为引入的噪点或者导致突变信号的目标结构被过度平滑。针对选取的同质区域A3和A4，所提方法均能在取得良好去噪效果的同时，实现对于突变目标的提取，并在所提取的图像中完整地保留目标的结构。表3得到的突变信号图像熵值计算结果进一步说明了两种分解方法对于异常目标分离效果差异。

图8 上海浦东机场区域哨兵1号数据的相干斑抑制结果Fig.8 Speckle suppression results for Sentinel-1 data in Shanghai pudong airport region

为了验证不同方法对于图像细节的保持特性，分别选取了近海区域(如图9所示)以及海上舰船区域(如图10所示)，并利用边缘检测技术对所有算法的去噪结果进行图像结构信息提取。可以看出，本文提出的方法能很好地保持近海区域的地物的结构信息，并没有因相干斑噪声抑制导致图像细节信息丢失。与之对比，RABASAR等方法并没有很好地保持图像的结构信息，从而使得边缘提取结果并不准确。在舰船区域结果中，可以看出所提方法很好地将舰船与背景噪声进行分离，使得舰船边缘提取结果并没有受到背景噪声的干扰。与之相比，RABASAR等方法边缘提取结果模糊，不能准确地反映目标的结构特性。

图9 近海区域图像边缘提取结果对比Fig.9 Comparison of image edge extraction results in offshore areas

图10 海上舰船区域图像边缘提取结果对比Fig.10 Comparison of image edge extraction results of ship region

4 结语

本文针对时序SAR图像中的相干斑噪声，提出了一种基于对数域加性信号分解技术的相干斑抑制方法，所提方法不仅能很好地对相干斑进行抑制，并且能对时间序列中的稳定信号成分以及沿时间维度突变信号成分进行有效分离，得到的突变信号能展示出观测时间内由于人造目标出现或地物变化产生的强散射现象，为后续的图像解译提供了有效的参考数据。为了验证方法的有效性，本文选取了现有的4种主流时序SAR图像相干斑抑制方法进行对比，并且用PSNR等几种指标定量地分析了不同方法在仿真数据上的去噪效果。此外，本文选取了上海近海以及浦东机场区域的哨兵1号时序SAR数据，并将上述方法所得结果进一步进行对比分析，均能验证所提方法的有效性。在后续工作中，会进一步验证本文所提方法在观测角度大幅变化情况下获取的SAR时序数据上的相干斑抑制效果。