合成孔径雷达(SAR) 因相干成像系统存在固有缺陷，导致其影像上遍布大量斑点状噪声，又称相干斑噪声。相干斑噪声是典型随机乘性噪声，严重限制了SAR影像解译技术的发展与应用^[1]。因而，抑制相干斑噪声一直是SAR影像处理中的关键步骤之一。

当前，相干斑噪声抑制方法主要包括：基于空域抑制法^[2-4]、基于变换域抑制法^{[1, 5-6]}与基于偏微分方程各向异性扩散抑制法^[7-11]。其中各向异性扩散依靠边缘检测算子，采取非边缘区平滑，边缘区尽量少甚至不平滑的方式，实现相干斑噪声抑制的同时，可保持甚至增强影像边缘与纹理信息，被广泛用于图像去噪、边缘检测和图像分割等领域。

1990年，Perona与Malik^[7]提出了P-M模型，首次采用各向异性扩散抑制光学影像噪声(加性噪声)，并取得了较好效果。但将P-M模型运用到富含乘性噪声的SAR影像中，较难获取理想效果，甚至起相反作用。针对此问题，Yu与Acton^[8]提出SRAD (speckle reducing anisotropic diffusion) 算法，成功将基于各向异性扩散的相干斑抑制方法拓宽到SAR领域。SRAD是P-M模型的发展，其利用瞬时变差系数代替梯度改进了边缘检测算子，有效提高了边缘检测的准确性。边缘检测的准确性是各向异性扩散滤波方法的关键^[12]。Aja-Fernández与Alberola-López^[9]利用更大尺度局部窗改进了SRAD边缘检测算子，提出了抑制相干斑噪声更好的DPAD (detail preserving anisotropic diffusion) 算法。Liu等人^[10]根据影像结构自适应调整窗口的大小和方向，提出了AWAD (adaptive window anisotropic diffusion) 算法，从而能更优检测影像边缘信息。张良培等人^[13]通过对边缘图进行直方图阈值化处理，有效保持了边缘信息。朱磊与程东^[14]利用平行窗通过比率运算生成边缘强度映射，并作为边缘检测算子，有效提升了算法边缘保护能力。陈少波等人^[15]利用算术-几何均值比作为边缘检测算子，在具有大范围边缘、亮线和全向纹理等特性的SAR影像中已获取较好效果。但上述SRAD算法、DPAP算法与AWAD算法以及基于P-M模型与SRAD的诸多算法的边缘检测均需要依赖同质区域的方差与均值的估计，然而同质区受到相干斑噪声干扰，对其方差与均值的准确估计十分困难。李金才等人^[11]为解决上述有问题，提出了一种基于图像熵的各向异性扩散滤波方法(IEAD)，该方法用图像熵构建新的边缘检测算子，可避开均值与方差的估计，有效提高了边缘检测能力。但需要指出，图像熵的准确获取需足够样本，而样本数过大又将增加异质区的可能性，因此样本数的准确选取较难确定。为进一步解决上述问题，需要寻找更有效的边缘检测算子。2013年，张闯等人^[16]利用欧氏距离有效地检测出图像的边缘信息，为构建边缘检测算子提供一个新的思路。

针对欧氏距离在边缘检测中的成功运用，在传统各向异性扩散相干斑抑制方法基础上，提出一种引入欧氏距离的各向异性扩散(EDAD) 相干斑抑制方法。EDAD算法解决了传统各向异性扩散方法中噪声均值与方差较难准确估计与图像熵样本不足等问题，通过比较区域间像素点的欧氏距离，准确寻找异质区，自适应构建新的各向异性扩散系数。利用各向异性扩散系数，完成相干斑较优抑制，并与SRAD算法、DPAD算法、AWAD算法、IEAD算法进行比较，从相干斑抑制效果与边缘保持能力两个方面对本文算法进行验证和评价。

相干斑是一种典型的乘性噪声，严重影响SAR影像观测质量。相干斑噪声对SAR数据的影响可描述为

$ I\left( {x,y} \right) = {I^R}\left( {x,y} \right) \cdot N\left( {x,y} \right) $

(1)

式中，$(x, y)$表示像素坐标；$I$表示观测影像的强度值，即受到相干斑噪声影响的强度信息；${I^R}$表示影像真实强度信息；$N$表示相干斑噪声。

相干斑抑制的目标是去除乘性噪声$N$，还原影像真实像素信息${I^R}$。为实现SAR相干斑抑制目标，Yu与Acton^[8]在P-M模型^[7]基础上，提出了SRAD算法。SRAD的偏微分方程，可用下式表示

$ \begin{array}{*{20}{c}} {\frac{{\partial I\left( {x,y;t} \right)}}{{\partial t}} = div\left[ {c\left( f \right) \cdot \nabla I\left( {x,y;t} \right)} \right]}\\ {I\left( {x,y;0} \right) = {I_0}\left( {x,y} \right)}\\ {\partial I\left( {x,y;t} \right)/\partial \mathit{\boldsymbol{n}}\left| {_{\partial \mathit{\boldsymbol{B}}}} \right. = 0} \end{array} $

(2)

式中，$t$表示时刻，${I_0}$表示初始噪声影像，$div$表示散度算子，Δ表示梯度算子，$\partial \boldsymbol{B}$表示$\boldsymbol{B}$的边界，$\boldsymbol{n}$表示垂直于边界的单位矢量，$c (f)$表示扩散系数。扩散系数是各向异性扩散的核心，决定扩散尺度，是一种非负单调递增函数，可表示为

$ c\left( f \right) = \frac{1}{{1 + \left[ {{f^2}\left( {x,y;t} \right) - f_0^2\left( t \right)} \right]/\mathit{\boldsymbol{T}}}} $

(3)

或

$ c\left( f \right) = \exp \left\{ { - \left[ {{f^2}\left( {x,y;t} \right) - f_0^2\left( t \right)} \right]/\mathit{\boldsymbol{T}}} \right\} $

(4)

式中，$T$表示扩散阈值；${f_0}\left (t \right)$表示尺度函数；$f (x, y; t)$表示瞬时变差系数，计算公式为

$ f\left( {x,y;t} \right) = \sqrt {\frac{{\left( {1/2} \right){{\left( {\left| {\nabla I} \right|/I} \right)}^2} - \left( {1/{4^2}} \right.{{\left( {\left| {{\nabla ^2}I} \right|/I} \right)}^2}}}{{{{\left[ {1 + \left( {1/4} \right)\left( {\left| {{\nabla ^2}I} \right|/I} \right)} \right]}^2}}}} $

(5)

式中，${\nabla ^2}$表示拉普拉斯算子。

观察式(3) 或式(4)，当${f^2}\left ({x, y; t} \right)-{f^2}_0\left (t \right)$远大于$T$时，$c (f)$趋于0，此时扩散停止；当${f^2}\left ({x, y; t} \right)-{f^2}_0\left (t \right)$远小于$T$时，$c (f)$趋于1，此时扩散表现为各向同性。因此，阈值$T$的选取对相干斑抑制效果产生影响，SRAD算法中给出阈值计算公式为

$ T = f_0^2\left( t \right) \cdot \left( {1 + f_0^2\left( t \right)} \right) $

(6)

式中

$ {f_0}\left( t \right) = \frac{{\sqrt {\mathit{var}\left[ {z\left( t \right)} \right]} }}{{mean\left[ {z\left( t \right)} \right]}} $

(7)

式中，$var[z (t)]$与$mean[z (t)]$分别是$t$时刻同质区域的方差与均值。为自动确定${f_0}\left (t \right)$，估计计算公式为

$ {f_0}\left( t \right) \approx {f_0}\exp \left( { - ct} \right) $

(8)

式中，$c$为常数，$N$是SAR强度影像，${f_0}$取$1/\sqrt N $。

Aja-Fernández与Alberola-López^[9]提出了DPAD算法，用Kuan滤波系数代替SRAD算法原扩散系数，即

$ c\left( f \right) = \frac{{1 + 1/{f^2}\left( {x,y;t} \right)}}{{1 + 1/{f^2}\left( {u;t} \right)}} $

(9)

式中

$ {f^2}\left( u \right) = {\sigma ^2}/{u^2} $

(10)

$ {f^2}\left( {x,y} \right) = \mathit{var}\left[ {I\left( {x,y} \right)} \right]/mean\left[ {{I^2}\left( {x,y} \right)} \right] $

(11)

式中，${\sigma ^2}$与$u$分别表示噪声的方差与均值，噪声变差系数${f^2}\left (u \right)$的准确估计对算法结果产生重要影响。

Liu等人^[10]提出了AWAD算法，在SRAD基础上通过自适应改变滤波窗口尺寸提高相干斑抑制结果。滤波窗口的高$h$与宽$w$可利用方向导数获取

$ h = \frac{a}{{{g_1} + 1}} $

(12)

$ w = \frac{a}{{{g_2} + 1}} $

(13)

式中，$a$表示尺度参数，${g_1}$与${g_2}$分别表示像素点最大与最小方向导数。

此外AWAD算法定义了一个新的扩散函数

$ c\left( f \right) = \frac{1}{{\sqrt {1 + {{\left( {f\left( {x,y;t} \right) - T} \right)}^2}} }} $

(14)

李金才等人^[11]提出IEAD算法，在影像边缘检测中引入图像熵概念，将SARD瞬时变差系数改进为

$ \begin{array}{*{20}{c}} {f\left( {x,y;t} \right) = }\\ { - \sum\limits_{t \in \eta \left( {x,y} \right)} {P\left( {x,y;t} \right)\log P\left( {x,y;t} \right)} } \end{array} $

(15)

式中，$\eta (x, y)$表示点$(x, y)$的邻域点，$P (x, y; t)$表示像素点灰度值在滑动窗口中个数的比例。

IEAD算法采用AWAD扩散系数(式(14))，其扩散阈值$T$通过全部滑动窗口取均值得到。

EDAD算法在传统各向异性扩散相干斑抑制的基础上，引入了基于欧氏距离的边缘检测。

欧氏距离被广泛应用于图像处理领域，是两个像素近似程度的数学体现。欧氏距离越小，参与比较的双方相似程度越高；反之，相似程度则越低。此外，当用两个区域间欧氏距离值代替单纯两个像素间欧氏距离值作为相似程度的描述时，像素差异将被放大^[16]，有益于在影像中寻找异质区(即边缘区)，从而为各向异性扩散中引入欧氏距离提供契机。采用区域间欧氏距离代替单纯像素间欧氏距离，即

$ \begin{array}{*{20}{c}} {D = \sum\limits_{\alpha = - h}^h {\sum\limits_{\beta = - w}^w {\left| {{I_k}\left( {{x_k} + \alpha ,{y_k} + \beta } \right)} \right.} } - }\\ {{{\left. {{I_l}\left( {{x_l} + \alpha ,{y_l} + \beta } \right)} \right|}^2}} \end{array} $

(16)

式中，$D$表示两个像素间的欧氏距离，$h$与$w$分别表示参与计算区域窗口高与宽的1/2，$k$与$l$分别表示$k$点与$l$点。

为进一步描述像素点与周围邻域点的近似程度，依据式(16) 逐一计算该像素点与各邻域内像素点的区域欧氏距离，取均值并作为该像素点与周围邻域点近似程度的描述量。将描述量定义为新的边缘检测算子，具体算法如下：

1) 如图 1所示，首先以目标像元$(x, y)$为中心建立$M \times N$大小的处理窗；

2) 在处理窗内，以$(x, y)$为中心建立$m \times n$大小的目标窗($m \times n$小于$M \times N$)、以处理窗口左上起始像元$({x_0}, {y_0})$为中心建立$m \times n$大小的对应窗。根据式(16) 计算目标窗与对应窗的区域欧氏距离；

3) 在处理窗内，从左至右，从上至下依次逐像元移动对应窗，逐一计算目标窗与对应窗的区域欧氏距离；

4) 求取处理窗内获取的全部区域欧氏距离的均值，并将其赋值于目标像元；

5) 等待整幅影像全部像素点(即目标像元) 计算完成，即完成了边缘检测算子的计算工作。

EDAD仍然采用SRAD偏微分方程，与SRAD主要不同之处在于边缘检测算子的改进。EDAD利用欧氏距离度量像素间相似性，进而实现各向异性扩散边缘检测。EDAD边缘检测算子可描述为

$ {f_{x,y}} = \frac{1}{n}\sum\limits_{k = 0}^{n - 1} {{D_k}} $

(17)

式中，$x$与$y$表示当前像素点位置，$n$表示处理窗像素个数，$D_k$表示处理窗内第$k$对应窗与目标窗的区域欧氏距离值。扩散阈值$T$可由整幅影像的欧氏距离均值获取，即

$ T = \frac{1}{{H \times W}}\sum\limits_{i = 0}^{H - 1} {\sum\limits_{j = 0}^{W - 1} {{f_{x,y}}} } $

(18)

式中，$H$与$W$分别表示整幅影响的高与宽。扩散系数为

$ c\left( f \right) = \frac{1}{{\sqrt {1 + {{\left( {{f_{x,y}} - T} \right)}^2}} }} $

(19)

各项异性扩散算法的实质是偏微分方程的求解，求解方法可采用Jacobi迭代算法^[8]。设$\Delta t$与$\Delta h$分别为足够小的时间步长与空间步长。时间与空间的离散坐标从而表示为

$ \begin{array}{l} t = n\Delta t,\;\;\;n = 0,1,2, \cdots \\ x = i\Delta h,\;\;\;i = 0,1,2, \cdots ,H - 1\\ y = j\Delta h,\;\;\;j = 0,1,2, \cdots ,W - 1 \end{array} $

(20)

于是，${I^n}_{x, y}=I\left ({x, y, t} \right)=I\left ({i\Delta h, j\Delta h, n\Delta t} \right)$。EDAD偏微分方程的计算步骤如下：

1) 首先依据式(16) 与式(17) 逐一计算影像像素点的边缘检测算子$f$；

2) 然后依据式(18) 计算各向异性扩散阈值$T$；

3) 接着依据式(19) 逐一计算影像像素点的扩散系数$c (f)$；

4) 再计算式(2) 偏微分方程中$c\left ({{f_{x, y}}} \right) \cdot \Delta I\left ({x, y; t} \right)$的散度，计算公式为

$ \begin{array}{*{20}{c}} {d_{x,y}^n = \frac{1}{{\Delta {h^2}}}\left[ {c_{x + 1,y}^n\left( {I_{x + 1,y}^n - I_{x,y}^n} \right) + c_{x,y}^n\left( {I_{x - 1,y}^n - I_{x,y}^n} \right) + } \right.}\\ {\left. {c_{x,y + 1}^n\left( {I_{x,y + 1}^n - I_{x,y}^n} \right) + c_{x,y}^n\left( {I_{x,y - 1}^n - I_{x,y}^n} \right)} \right]} \end{array} $

(21)

5) 最后计算微分方程的数值，其近似为

$ I_{x,y}^{n + 1} = I_{x,y}^n + \frac{{\Delta t}}{4}d_{x,y}^n $

(22)

为评估EDAD相干斑抑制效果，实验采取两种评价指标进行评价。

1) 等效视数${V_{{\rm{ENL}}}}$。等效视数可反映相干斑抑制程度，计算公式为

$ {V_{{\rm{ENL}}}} = {E^2}/\mathit{var} $

(23)

式中，$E$与$var$分别表示相干斑抑制后影像均匀区域的均值与方差。${V_{{\rm{ENL}}}}$越大表示噪声抑制效果越好。

2) 比值影像(原始影像与相干斑抑制后影像的比) 的均值$V_E$与方差$V_V$。$V_E$可反映原始影像辐射特性的保持能力，理想值为1；$V_V$可反映相干斑抑制程度，理想值为(4-π)/($L$·π)，其中$L$为影像视数。

截取两景StripMap模式下天津测区3 m分辨率的TanDEM-X影像进行实验研究。图 2(a)与图 3(a)分别为裁剪后的两景影像，像素大小均为1 024×1 024。其中图 2(a)多为农作物区域，相干斑较分散，而图 3(a)中存在大量建筑物，相干斑明显且集中在建筑物周边。实验分别采用SRAD算法、DPAD算法、AWAD算法、IEAD算法与EDAD算法对该两景影像进行相干斑抑制。各算法时间步长均设置为0.1，迭代次数60次。在SRAD阈值估计中，$\alpha $取1/6，${f_0}$取1。在DPAD图像变差系数计算中，滑动窗口尺寸为5×5，噪声局部变差系数估计采用均值模型。AWAD依据方向梯度的大小依次选取3×3，5×3和3×5大小的窗口。IEAD计算图像熵采取5×5窗口。EDAD计算边缘检测算子采用9×9窗口，计算区域欧氏距离采用5×5窗口。实验结果如图 2与图 3所示，评价指标如表 1与表 2所示。

从图 2(a)与图 3(a)原始影像中可以直观看出相干斑的存在严重影响了影像质量。针对图 3农作物区域，各类相干斑抑制算法处理结果均可以较好地抑制相干斑，且差异微弱。这是由于相干斑噪声较弱。由此可见，当相干斑噪声不强，实验采取的各类各向异性扩散相干斑抑制方法均有较好抗噪性。但就目视而言，SRAD算法效果最差，EDAD算法的目视效果最好。

整体上看，图 3的相干斑抑制效果比图 2明显，这是由于图 3相干斑噪声较多且集中，与周围像素差异较大，易于判定为异质区，从而相干斑抑制效果明显。从图 3(b)-图 3(f)可以看出，较为集中的相干斑均被有效去除，去除后影像已不存在明显相干斑点。但是各类各向异性扩散相干斑抑制方法仍存在很大差异。图 3(b)中，SRAD算法虽较好平滑了相干斑噪声，但是影像因此变得相对模糊。这是由于算法中瞬时变差系数易受到强噪声的干扰，影响边缘检测的准确性。相比于SRAD算法，DPAD算法对噪声变差系数估计方法进行改进，有效提高边缘保持能力，如图 3(c)所示。对比图 3(b)，图 3(c)的边缘轮廓相对清晰，但是DPAD算法仍然会模糊掉影像细小纹理。AWAD算法通过自适应控制窗口大小和方向，IEAD算法利用图像熵反演出同质区与异质区，二者均可提高边缘检测能力，可获取相对较好的相干斑抑制效果，但局部区域仍然受到强噪声的干扰，如图 3(d)与图 3(e)白色圆形边框标记区域所示。EDAD算法利用区域欧氏距离值反映像素差异性，自适应完成边缘检测，可获取较优边缘检测效果，目视效果最好，如图 3(f)所示。与图 3(b)及图 3(c)相比，图 3(f)有更强的边缘与纹理细节保持能力，而与图 3(d)及图 3(e)相比，图 3(f)有更好的抗强噪声能力。综合图 2与图 3可以看出，EDAD算法在强噪声区域仍具有较好抗噪能力。

此外，通过不同算法相干斑抑制的参数对比，同样可以显示EDAD算法的优势。表 1与表 2中${V_{{\rm{ENL}}}}$计算所选取的同质区为图 2(a)与图 3(a)所示的白色矩形边框标记的区域。可以看出，EDAD算法在区域1与区域2内选取的同质区的${V_{{\rm{ENL}}}}$具有相对较大的数值；比值影像$V_E$与$V_V$均较优。相干斑抑制参数比对结果与分析一致，体现EDAD算法相干斑抑制的优势。

针对SAR影像解译易受相干斑干扰等问题，利用欧氏距离作为边缘检测算子，提出一种引入欧氏距离的各向异性扩散相干斑抑制方法。该方法避开了传统各向异性扩散方法中噪声均值与方差较难准确估计等问题，充分发挥欧氏距离边缘检测算子边缘检测准确性与简易性等优势，较好实现了相干斑的抑制。实验过程中，利用等效视数，相干斑抑制前后比值影像均值、方差3个算法评估参数，验证本文算法较其他传统各向异性扩散相干斑抑制方法具有更优边缘保持与相干斑抑制能力。实验结果证明，该算法对分布较散的弱相干斑区域与分布较集中的强相干斑区域均有较好适用性。此外，算法中各类参数的选取将影响算法质量，最优参数自适应选取将是今后研究工作的重点。