CN102096921A - 基于邻域对数比值及各向异性扩散的sar图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于邻域对数比值及各向异性扩散的SAR图像变化检测方法，它涉及遥感图像处理领域，主要解决SAR图像变化检测的差异图构造受SAR图像斑点噪声影响严重的问题，其实现步骤为：(1)根据邻域对数比值方法，构造两幅不同时间相同地域图像I₁和I₂的差异影像图I_L；(2)对差异图I_L进行自适应窗口各项异性扩散滤波处理，得到最终的差异图滤波结果图

(3)利用OSTU阈值算法，对最终差异图滤波结果图

进行阈值分割，得到变化检测结果图

用邻域对数比值方法构造差异影像图，可对差异影像图的直方图进行压缩，从而有效去除变化检测结果图中的杂点，自适应窗口各向异性扩散滤波的边缘保持能力良好，不会模糊图像边缘，因此得到的变化检测结果图比较精细。

Description

基于邻域对数比值及各向异性扩散的SAR图像变化检测方法

技术领域

本发明属于遥感图像处理领域，涉及遥感图像的变化检测，具体提供一种基于邻域对数比值及各向异性扩散的SAR图像变化检测方法，可用于对SAR图像变化的检测。

背景技术

随着合成孔径雷达(synthetic aperture radar，SAR)技术的快速发展，SAR系统可以全天候，全天时获取图像数据，是较好的变化检测图像源，遥感图像变化检测是通过对同一地区的不同时期的两幅遥感图像进行比较分析，根据图像之间的差异得到所需的地物变化信息。在近十几年来，SAR图像变化检测技术在很多方面得到了应用，例如，森林覆盖变化，地震区域的定位和灾害评估等领域具有广泛的应用前景。

随着计算机和遥感成像技术的不断发展，SAR图像变化检测技术也在不断提高和改进中，在各国学者的努力下提出了以下几种经典的变化检测算法：

图像差值法：其主要过程是将两幅不同时间相同地域的遥感图像对应像素的灰度值相减得到差异影像图。理论上，在得到的差值图像上，差值为0或接近0的认为是不变区域，不为0的认为是变化区域。这种方法的优点在于理论相对简单，容易理解和掌握。但是缺点也比较明显，由于该方法主要通过分析地物光谱值改变的灰度差值图像来实现变化检测，但在某些情况下，仅仅利用光谱特征的差值图像难以实现地物的变化检测，易产生“伪变化”信息。

图像比值法：其主要过程是计算两幅不同时间相同地域的遥感图像对应像素灰度值的比值以得到差异影像图，其中像素比值为1或近似为1的认为是未发生变化的区域，反之，远大于或远小于1的认为是发生变化的区域。该方法可以克服乘性噪声的干扰，但该方法是假设比值图像呈正态分布的，在很多实际问题中该假设并不总是成立的，缺少适应性。

分类后比较法：其主要过程是采用相同的分类体系分别对两幅不同时间相同地域的遥感图像进行独立的分类，然后比较分类结果，从中提取变化信息。如果对应像素的分类类别相同，则认为该像素没有发生变化，反之，发生了变化。该方法可以回避所用多时相遥感图像数据应获取环境条件不同和传感器不同所带来的辐射归一化问题，并可直接获取变化的类型，数量，位置。但该方法工作量很大，并且存在精度方面的缺陷。

在SAR图像变化检测领域，差异影像图的斑点噪声去除是一个重要的环节，它极大地影响着后继处理的结果。虽然对差异影像图进行斑噪抑制处理可有效减少斑点噪声对变化检测结果的影响。但是传统的图像去噪方法，如中值滤波，自适应加权中值滤波，Lee滤波，kuan滤波，fronst滤波等，细节保持能力较差，常会模糊图像的边缘，会影响最后变化检测结果的正确性。因此变化检测中图像去噪算法的选取十分重要。

发明内容

本发明要解决的技术问题在于：克服上述现有技术的不足，提高SAR图像变化检测的正确率，根据SAR图像变化检测自身的特点，提出一种基于邻域对数比值及自适应窗口各向异性扩散滤波的SAR图像变化检测方法，可抑制SAR图像的部分斑点噪声，提高变化检测结果的正确率。

本发明的技术方案是：利用邻域信息，构建了一个由对数比值算法演化而来的邻域对数比值方法，根据邻域对数比值方法构造两幅不同时间相同地域图像的差异影像图I_L；再通过对差异影像图I_L进行基于自适应窗口的各向异性扩散滤波处理得到滤波后结果图

最后对滤波后结果图

进行OSTU阈值分割得到变化结果图。其具体实现步骤如下：

(1)根据邻域对数比值原理，按照如下步骤构造两幅不同时间相同地域的图像I₁，I₂的差异影像图I_L；

1a)分别取得两时像SAR图像I₁，I₂在同一位置x上的像素值邻域集合S₁(x)和S₂(x)，其大小为N×N，N∈{3，5，7，9}；

1b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差异图I_L在位置x上的像素灰度值I_L(x)：

$I_L\left(x\right)=\frac{1}{N\×M}\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\min (\frac{\log \left(S_1\left(x_i\right)\right)}{\log \left(S_2\left(x_i\right)\right)},\frac{\log \left(S_2\left(x_i\right)\right)}{\log \left(S_1\left(x_i\right)\right)})$

其中，S₁(x_i)和S₂(x_i)分别表示图像I₁，I₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，I_L(x)值越小，说明图像I₁，I₂在位置x上的差别越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；

1c)对图像I₁，I₂的每一个像素点从左到右，从上到下重复步骤1a)和1b)，得到差异影像图I_L；

(2)对差异图I_L进行自适应窗口各项异性扩散滤波处理，得到滤波后的差异图NI_L，所述的自适应窗口各向异性扩散滤波处理包括如下步骤：

2a)根据比值边缘检测原理，以差异图像的位置x为中心取一个M×M大小的窗口，M∈{3，5，7，9}，分别沿着位置x的4个直线方向，即0°，45°，90°，145°这四个方向，将窗口分为两个不重叠区域，并计算直线两侧不重叠区域R1和R2内各N点的灰度均值

和

然后求两均值之比

其中i∈{1，2，3，4}，分别表示4个方向的两侧不重叠区域均值之比。若窗口位于图像中的均匀区域时，则

与

相近，d_i趋于1，反之，当窗口中心位于不同区域的交界区时，由于两个区域的统计特性不同，d_i将小于1，d_i越小，说明区域差异越大，窗口位于边缘的可能性越高。设定一个阈值T，当4个方向中的均值比最小值min(d_i)＜T时，则差异图的位置x被判定为边缘点，同时以差异图像位置x为中心的自适应窗口方向取min(d_i)方向。反之，差异图的位置x被判定为处于均匀区域，窗口设定为一个正方形窗口；

2b)如果差异图的位置x被判定为为边缘点，则自适应窗口的高度和宽度由以下两个公式决定，其中窗口的高度为

宽度为

k为尺度参数，反之，若差异图的位置x处于均匀区域，则窗口取为M×M大小的正方形窗口；

2c)在确定的自适应窗口中计算差异图的位置x的扩散系数q，其中q为：

$q=\sqrt{\mathrm{Var}\left(I_i\right)/{\stackrel{\‾}{I}}^2}=\sqrt{\frac{1}{\mathrm{\η}}\mathrm{\Σ}{(I_i-\stackrel{\‾}{I})}^2/{\stackrel{\‾}{I}}^2}$

这里，η为计算参数q的邻域数量，大小为h*w，

为窗口内的像素均值，Var为方差算子，I_i为以差异图像位置x为中心的自适应窗口的第i个元素；

2d)差异图的位置x的扩散系数方程为

其中噪声误差系数

t为滤波算法迭代次数；

2e)利用计算得到的扩散系数方程c(q)，根据各向异性扩散方程原理得到差异影像滤波结果图在位置x上的值NI_L(x)：

${\mathrm{NI}}_L\left(x\right)=I_L\left(x\right)+\mathrm{\λ}\mathrm{div}[c\left(q\right)*\▿I_L\left(x\right)]$

其中λ为稳定系数，取常数，

是梯度算子，div是散度算子；

2f)对I_L从左到右，从上到下重复步骤2e)，得到差异图的滤波结果图NI_L；

2g)根据预先设定的滤波算法迭代次数t，重复步骤2f)t次，得到最终的差异图滤波结果图

(3)利用OSTU大津阈值算法，对最终的差异图滤波结果图

差异图滤波结果图

进行阈值分割，得到变化检测结果图

本发明与现有的技术相比具有如下有益效果：

1、本发明采用邻域对数比值方法构造差异影像图，可对差异影像图的直方图进行压缩，有效减少了处于变化类和非变化类之间的区域对变化检测结果的影响，对于变化检测结果图中杂点的去除效果非常显著；

2、本发明针对邻域对数比值方法构造的差异影像图，创建的自适应窗口各向异性扩散滤波处理方法，可有效减少差异影像图中的噪声影响，同时该滤波基于各向异性扩散原理，与传统的滤波算法相比，边缘保持能力良好，不会模糊图像边缘，对图像的空间分辨率的影响非常小；

3、仿真结果表明，本发明采用的邻域对数比值方法构造差异影像图方法和自适应窗口各向异性扩散滤波处理方法与比值R算子，对数比值LR算子构造差异影像图方法和Lee滤波，SRAD滤波处理方法相比，正确检测率高，错检和漏检个数低。

附图说明

图1是本发明的主要流程图；

图2是本发明设计的4个自适应滤波窗口方向；

图3是本发明设计的4个自适应滤波窗口模型；

图4是第一组实验仿真图，图4(a)和图4(b)的拍摄时间分别为1997.05和1997.08，大小均为290×350，图4(c)为参考图；

图5是对应图4的不同方法构造的差异影像图，其中图5(a)，图5(b)和图5(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图；

图6是对应图5的OSTU阈值分割变化检测结果图，其中图6(a)，图6(b)图6(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图在经过自适应窗口各向异性扩散滤波后的OSTU阈值分割结果图；

图7是对应图5的差异影像图的直方图分布，其中图7(a)，图7(b)，图7(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图直方图分布；

图8是第二组实验仿真图，图8错误！未找到引用源。(a)和图8(b)的拍摄时间分别为1995.9和1996.7，大小均为412×300，图8(c)为参考图；

图9是对应图8的不同滤波算法得到的差异影像图滤波后结果图，其中图9(a)，图9(b)和图9(c)分别表示由Lee滤波，SRAD滤波和本发明方法AWSRAD滤波处理后的滤波结果图；

图10是对应图9的OSTU阈值分割变化检测结果图，其中图10(a)，图10(b)和图10(c)分别表示由Lee滤波，SRAD滤波和本发明方法处理后得到的OSTU阈值分割变化检测结果图。

具体实施方式

参照图1，本发明的基于邻域对数比值及各向异性扩散的SAR图像变化检测方法包括如下步骤：

步骤一：根据邻域对数比值原理，按照如下步骤构造两幅不同时间相同地域的图像I₁，I₂的差异影像图I_L。

1a)分别取得两时像SAR图像I₁，I₂在同一位置x上的像素值邻域集合S₁(x)和S₂(x)，其大小为N×N，N∈{3，5，7，9}；

1b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差异图I_L在位置x上的像素灰度值I_L(x)：

$I_L\left(x\right)=\frac{1}{N\×M}\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\min (\frac{\log \left(S_1\left(x_i\right)\right)}{\log \left(S_2\left(x_i\right)\right)},\frac{\log \left(S_2\left(x_i\right)\right)}{\log \left(S_1\left(x_i\right)\right)})$

其中，S₁(x_i)和S₂(x_i)分别表示图像I₁，I₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，I_L(x)值越小，说明图像I₁，I₂在位置x上的差别越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；

1c)对图像I₁，I₂的每一个像素点从左到右，从上到下重复步骤1a)和1b)，得到差异影像图I_L；

步骤二：对差异图I_L进行自适应窗口各项异性扩散滤波处理，得到最终的差异图滤波结果图

对差异图I_L作自适应窗口各项异性扩散滤波处理的步骤如下：

2a)根据比值边缘检测原理，以差异图像的位置x为中心取一个M×M大小的窗口，M∈{3，5，7，9}，分别沿着位置x的4个直线方向，如图2(a)，(b)，(c)，(d)所示，将窗口分为两个不重叠区域，并计算直线两侧不重叠区域R1和R2内各N点的灰度均值

和

然后求两均值之比其中i∈{1，2，3，4}，分别表示4个方向的两侧不重叠区域均值之比。若窗口位于图像中的均匀区域时，则

与

相近，d_i趋于1，反之，当窗口中心位于不同区域的交界区时，由于两个区域的统计特性不同，d_i将小于1，d_i越小，说明区域差异越大，窗口位于边缘的可能性越高。设定一个阈值T，当4个方向中的均值比最小值min(d_i)＜T时，则差异图的位置x被判定为边缘点，同时以差异图像位置x为中心的自适应窗口方向取min(d_i)方向。反之，差异图的位置x被判定为处于均匀区域，窗口设定为一个正方形窗口；

2b)如果差异图的位置x被判定为为边缘点，则自适应窗口的高度和宽度由以下两个公式决定，其中窗口的高度为宽度为

k为尺度参数，其中，4个窗口的模型如图3(a)，(b)，(c)，(d)所示，反之，若差异图的位置x处于均匀区域，则窗口取为M×M大小的正方形窗口；

2c)在确定的自适应窗口中计算差异图的位置x的扩散系数q，其中q为：

$q=\sqrt{\mathrm{Var}\left(I_i\right)/{\stackrel{\‾}{I}}^2}=\sqrt{\frac{1}{\mathrm{\η}}\mathrm{\Σ}{(I_i-\stackrel{\‾}{I})}^2/{\stackrel{\‾}{I}}^2}$

这里，η为计算参数q的邻域数量，大小为h*w，

为窗口内的像素均值，Var为方差算子，I_i为以差异图像位置x为中心的自适应窗口的第i个元素；

2d)差异图的位置x的扩散系数方程为其中噪声误差系数

t为滤波算法迭代次数；

2e)利用计算得到的扩散系数方程c(q)，根据各向异性扩散方程原理得到差异影像滤波结果图在位置x上的值M_L(x)：

${\mathrm{NI}}_L\left(x\right)=I_L\left(x\right)+\mathrm{\λ}\mathrm{div}[c\left(q\right)*\▿I_L\left(x\right)]$

其中λ为稳定系数，取常数，

是梯度算子，div是散度算子；

2f)对I_L从左到右，从上到下重复步骤2e)，得到差异图的滤波结果图M_L；

2g)根据预先设定的滤波算法迭代次数t，重复步骤2f)t次，得到最终的差异图滤波结果图

步骤三：利用OSTU大津阈值算法，对最终差异图滤波结果图

进行阈值分割，得到变化检测结果图

3a)设差异图滤波结果图的灰度值为0到m-1级，灰度值i的像素为n_i，此时我们得到的像素总数为

各值概率为

然后用阈值T将其分成两组C₀＝{0～T-1}和C₁＝{T～m-1}；

3b)各组产生的概率如下：

C₀产生概率

C₁产生概率

C₀的平均值C₁的平均值

其中，是滤波结果图

的灰度平均值，

是阈值为T时灰度平均值，所以全部采样的灰度平均值为μ＝w₀μ₀+w₁μ₁；

3c)两组间的方差用下式求出：

${\mathrm{\δ}}^2\left(T\right)=w_0{({\mathrm{\μ}}_0-\mathrm{\μ})}^2+w_1{({\mathrm{\μ}}_1-\mathrm{\μ})}^2=\frac{{[\mathrm{\μ}\·w\left(T\right)-\mathrm{\μ}\left(T\right)]}^2}{w\left(T\right)[1-w\left(T\right)]}$

从1～m-1之间改变T，当δ²(T)为最大值时的T值，即为阈值，对滤波结果图

按照阈值T进行分割，小于T设为像素0，大于T设为像素255，得到变化检测结果图。

本发明的效果可以通过以下仿真进一步说明：

1、仿真参数

对于具有参考图的实验仿真图，可进行定量的变化检测结果分析：

①计算漏检个数：统计实验结果图中发生变化区域的像素个数，与参考图中变化区域的像素个数进行对比，把参考图中发生变化但实验结果图中检测为未变化的像素个数，称为漏检个数FN；

②计算错检个数：统计实验结果图中未发生变化区域的像素个数，与参考图中未变化区域的像素个数进行对比，把参考图中未发生变化但实验结果图中检测为变化的像素个数，称为错检个数FP；

③总检测错误数OE：等于漏检个数和错检个数的和；

2、仿真内容

①将本发明方法邻域对数比值NLR算子构造的差异影像图，与比值R算子，对数比值LR算子构造的差异影像图的直方图数据分布对比，以及变化检测结果对比和定量分析通过第一组实验仿真图完成；

②对差异影像图的Lee滤波，SRAD滤波处理和自适应窗口各向异性扩散AWSRAD滤波处理对应的变化检测效果对比和定量分析，通过第二组实验仿真图完成。

3、仿真实验结果及分析

①反应Ottawa地区水灾的SAR图像如图4所示，图4(a)和图4(b)的拍摄时间分别为1997.05和1997.08，大小均为290×350，图4(c)为参考图。

用不同算子构造的差异影像图如图5所示，其中图5(a)，图5(b)和图5(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图。图5差异影像图的直方图分布如图7所示，其中图7(a)，图7(b)，图7(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图直方图分布，从图7(c)可以看出本发明方法构造的差异影像图直方图分布在经过数据压缩后，变化区域与非变化区域的分布相互独立，介于变化类和非变化类之间的区域较少，这样可有效的去除杂点。用不同算子构造的差异影像图在经过自适应窗口各向异性扩散滤波后的OSTU阈值分割结果如图6所示，其中图6(a)，图6(b)和图6(c)分别表示由R算子，LR算子和本发明方法构造的差异影像图在经过自适应窗口各向异性扩散滤波后的OSTU阈值分割结果图，从图6可以看出，本发明方法构造的差异影像图的OSTU阈值结果图杂点较少，变化检测比较好。用不同算子构造差异影像图的变化检测结果分析如表1所示。

表1不同算子构造差异影像图的变化检测结果分析

从表1中可以看出，本发明方法在变化检测总错误数上，比其他对比实验方法都要少，而且在误检数与漏检数上相对都比较均衡。同时与参考图像图(4)对比可看出，本发明对杂点的去除效果比其他两个方法都要好。

②意大利撒丁岛区域的SAR图像如图8所示，图8错误！未找到引用源。(a)和图8(b)的拍摄时间分别为1995.9和1996.7，大小均为412×300，图8(c)为参考图。用NLR算子构造的差异影像图如图8(d)所示。差异影像图经过不同的滤波处理形成的滤波结果图如图(9)所示，其中图9(a)，图9(b)和图9(c)分别表示由Lee滤波，SRAD滤波和本发明方法AWSRAD滤波处理后的滤波结果图。由于本发明方法AWSRAD滤波是一个迭代算法，这里设定迭代次数t为3。从滤波结果图可以看出，由于SRAD滤波和本发明方法都是基于各向异性异性扩散原理的，所以滤波结果图的边缘保持的非常好，没有出现明显的边缘模糊现象，同时图像的分辨率比Lee滤波结果图要高得多，并且本发明采用自适应窗口估算扩散系数，估算结果更加精确，所以视觉效果较SRAD滤波更为清晰。用不同滤波算法得到的OSTU阈值结果图如图10所示，其中图10(a)，图10(b)和图10(c)分别表示由Lee滤波，SRAD滤波和本发明方法处理后得到的OSTU阈值分割变化检测结果图。从图10可以看出，本发明方法对差异影像图进行滤波处理得到的OSTU阈值结果图较其他两个方法杂点明显减少，同时与参考图像图8(c)进行对比可以看出，本发明方法得到的变化检测结果图对细小边缘部分的检测结果非常好，能有效的检测出细小边缘。用不同滤波算法得到的变化检测结果分析如表2所示。

表2不同滤波算法得到的变化检测结果分析

从表2可以看出本发明滤波方法AWSRAD滤波得到的滤波结果图在经过OSTU阈值分割后，变化检测总错误数比其他对比滤波方法都要少。并且漏检数FN和错检数FP相对较均衡，得到的变化检测结果图比较稳定。

Claims (3)

1.一种基于邻域对数比值及各向异性扩散的SAR图像变化检测方法，其特征在于如下步骤：

(1)根据邻域对数比值NLR方法，构造两幅不同时间相同地域图像I₁和I₂的差异影像图I_L；

(2)对差异图I_L进行自适应窗口各项异性扩散滤波处理，得到最终的差异图滤波结果图

(3)利用OSTU大津阈值算法，对最终差异图滤波结果图

进行阈值分割，得到变化检测结果图

2.根据权利要求1所述的基于邻域对数比值及各向异性扩散的SAR图像变化检测方法，其特征在于：所述步骤(1)构造差异影像图I_L的步骤如下：

2a)分别取得两时像SAR图像I₁，I₂在同一位置x上的像素值邻域集合S₁(x)和S₂(x)，其大小为N×N，N∈{3，5，7，9}；

2b)比较两个邻域集合S₁(x)和S₂(x)的相似性，得到差异图I_L在位置x上的像素灰度值I_L(x)：

$I_L\left(x\right)=\frac{1}{N\×N}\underset{i=1}{\overset{N\×N}{\mathrm{\Σ}}}\min (\frac{\log \left(S_1\left(x_i\right)\right)}{\log \left(S_2\left(x_i\right)\right)},\frac{\log \left(S_2\left(x_i\right)\right)}{\log \left(S_1\left(x_i\right)\right)})$

其中，S₁(x_L)和S₂(x_i)分别表示图像I₁和I₂在位置x上邻域集合S₁(x)和S₂(x)的第i个元素，I_L(x)值越小，说明图像I₁和I₂在位置x上的差别越大，则位置x上的像素点属于变化区域的可能性也越大，反之，该位置像素点属于非变化区域的可能性越大；

2c)对图像I₁和I₂的每一个像素点从左到右，从上到下重复步骤2a)和2b)，得到差异影像图I_L。

3.根据权利要求1所述的基于邻域对数比值及各向异性扩散的SAR图像变化检测方法，其特征在于：所述步骤(2)对差异图I_L进行自适应窗口各项异性扩散滤波处理的步骤如下：

3a)根据比值边缘检测原理，以差异图像的位置x为中心取一个M×M大小的窗口，M∈{3，5，7，9}，分别沿着位置x的四个直线方向即0°，45°，90°，145°，将窗口分为两个不重叠区域，并计算直线两侧不重叠区域R1和R2内各N点的灰度均值

和

然后求两均值之比

其中i∈{1，2，3，4}，分别表示四个方向的两侧不重叠区域均值之比；设定一个阈值T，当四个方向中的均值比最小值min(d_i)＜T时，则差异图的位置x被判定为边缘点，同时以差异图像位置x为中心的自适应窗口方向取min(d_i)方向，反之，差异图的位置x被判定为处于均匀区域，窗口设定为一个正方形窗口；

3b)如果差异图的位置x被判定为边缘点，则自适应窗口的高度和宽度由以下两个公式决定，其中窗口的高度为

宽度为

k为尺度参数，反之，若差异图的位置x处于均匀区域，则窗口取为M×M大小的正方形窗口；

3c)在确定的自适应窗口中计算差异图的位置x的扩散系数q：

$q=\sqrt{\mathrm{Var}\left(I_i\right)/{\stackrel{\‾}{I}}^2}=\sqrt{\frac{1}{\mathrm{\η}}\mathrm{\Σ}{(I_i-\stackrel{\‾}{I})}^2/{\stackrel{\‾}{I}}^2}$

式中，η为计算参数q的邻域数量，大小为h*w，为窗口内的像素均值，Var为方差算子，I_i为以差异图像位置x为中心的自适应窗口的第i个元素；

3d)差异图的位置x的扩散系数方程为其中噪声误差系数t为滤波算法迭代次数；

3e)利用计算得到的扩散系数方程c(q)，根据各向异性扩散方程得到差异影像滤波结果图在位置x上的值NI_L(x)：

${\mathrm{NI}}_L\left(x\right)=I_L\left(x\right)+\mathrm{\λ}\mathrm{div}[c\left(q\right)*\▿I_L\left(x\right)]$

其中λ为稳定系数，取常数，

是梯度算子，div是散度算子；

3f)对I_L从左到右，从上到下重复步骤3e)，得到差异图的滤波结果图NI_L；

3g)根据预先设定的滤波算法迭代次数t，重复步骤3f)t次，得到最终的差异图滤波结果图

Images