CN105608679A - 一种融合结构张量与非局域全变分的图像去噪方法 - Google Patents

Abstract

本发明涉及一种融合结构张量与非局域全变分的图像去噪方法，属于图像处理技术领域。该方法包括以下步骤：1)对带噪图像f进行N*N的分块，N的取值可为N＝3,5,7…；2)对图像去噪建模，建立融合结构张量与非局域全变分的目标函数；3)利用分裂Bregman算法求解目标函数得到去噪图像。本发明所述方法具有有效去噪同时保持图像的纹理特征与几何结构特征的优点，算法快速，同时能够减少阶梯效应的产生。

Description

一种融合结构张量与非局域全变分的图像去噪方法

技术领域

本发明属于图像处理技术领域，涉及一种融合结构张量与非局域全变分的图像去噪方法。

背景技术

图像去噪是图像处理领域的一个基础问题，旨在将图像中的噪声去除的同时尽可能的保留原始图像信息，例如边缘、纹理和细小图像结构等。图像去噪的实质是图像退化的逆过程，经典的去噪方法有频域处理的小波法，空域处理的高斯滤波、中值滤波等,这类方法通常在某几何领域内直接对像素点进行操作，取得了较好的平滑效果，但以丢失细节、纹理信息为代价。为了缓解这种矛盾，Yaroslavsky提出邻域滤波，以像素灰度相似性作为邻域间权重依据，双边滤波器进一步加入像素间的位置信息。但是仅考虑局域信息的滤波方法没有兼顾到图像的结构信息，因此容易模糊图像细节部分。Buades等提出的非局域滤波思想，将传统的空间邻域扩展到几何意义下的结构邻域，很好的利用了自然图像的自相似性，成为非局域方法在图像去噪领域的一个里程碑。Singer等认为非局域滤波实质是图像在块空间上的扩散，Milanfar等人提出了基于高阶核函数回归的图像去噪框架，建立起Bilateral滤波、非局部均值等方法的关系。Peter等用结构相似性代替L2距离作为度量因子，并以鲁棒统计权函数代替高斯函数，取得了较好的去噪效果。

近年来，基于PDE(PartialDifferentialEquation)的图像处理方法引起研究者的广泛关注，这类方法通过偏微分方程建立起图像处理模型与数学方程之间的关系，已形成完整的理论体系和数值方法。其中，变分PDE通过目标函数的最优化建立PDE方程，特别适用于不适定问题的求解。目标函数由保真项与正则项构成，正则项的定义通常基于先验知识，Rudin等人提出的经典图像去噪模型就是采用全变分(Totalvariation)正则项，该正则项基于图像分片光滑假设，允许存在不连续点，在抑制高频噪声的同时保持了边缘。之后，对正则项的改进不断推进：Shu等提出了增强梯度域的稀疏性和方向性的TV-L1模型，Chan等提出了高阶全变分模型，有效抑制了“阶梯效应”并很好的刻画了纹理。

发明内容

有鉴于此，本发明的目的在于提供一种融合结构张量与非局域全变分的图像去噪方法，具体包括以下技术方案：

一种融合结构张量与非局域全变分的图像去噪方法，包括以下步骤：

1)对带噪图像f进行N*N的分块，N的取值可为N＝3,5,7...；

2)对图像去噪建模，建立融合结构张量与非局域全变分的目标函数；

3)利用分裂Bregman算法求解目标函数得到去噪图像。

进一步，在步骤1)中，N的取值可为3-2M+1(M为自然数)的奇数。

进一步，在步骤2)中，具体包括：

21)建立基于全变分正则项的图像去噪模型

$E\left(u\right)=\underset{u\∈R^n}{\min }{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}|\▿u|$

其中f是带噪图像，u是待恢复的去噪图像，图像去噪的目的在于从f中恢复出清晰的图像u；式中第一项为数据项，第二项为平滑项，λ是平衡数据项与平滑项的权重因子，Ω为整个图像空间，为n维实数空间，▽为梯度算子，|▽u|也称为全变分正则项；

22)将传统的全变分正则项扩展到非局域空间，得到非局域全变分并作为正则项约束，进而得到基于非局域信息的全变分图像去噪模型

其中非局域算子▽_NLu由图像空间中任意两个像素点(x，y)∈Ω×Ω定义为：

${\▿}_{NL}u(x,y)=\left(u\right(y)-u(x\left)\right)\sqrt{w_{NL}(x,y)}$

其中u(x)，u(y)分别指u在像素点x，y处的值，w_NL：Ω×Ω为像素点x，y之间的边缘权重，且满足对称性，即w_NL(x，y)＝w_NL(y，x)；用于衡量块P_x和P_y之间的相似性,其中G_σ是标准偏差为σ的高斯卷积核，h为尺度参数；

23)融合图像的结构张量与非局域全变分：首先，定义结构张量矩阵为：

$J_p(\▿I_{\mathrm{\σ}})=G_{\mathrm{\ρ}}*(\▿I_{\mathrm{\σ}}\⊗\▿I_{\mathrm{\σ}})\stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]$

其中，I_σ是经过方差为σ的高斯滤波器的图像I，▽是梯度算子，是张量积运算，*代表卷积运算，G_ρ是方差为ρ的高斯滤波器，J₁₁，J₁₂，J₂₁，J₂₂分别为结构张量矩阵的四个元素；

求解对称半正定的结构张量矩阵J_ρ的两个正交特征向量：

$v_1=\left[\begin{array}{c}2J_{12}\\ J_{22}-J_{11}+\sqrt{{(J_{22}-J_{11})}^2+4J_{12}^2}\end{array}\right],v_1=\frac{v_1}{\left|\right|v_1\left|\right|},v_2=v_1^{\⊥}$

与对应的特征值：

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(J_{11}+J_{22}\±\sqrt{{(J_{22}-J_{11})}^2+4J_{12}^2})$

定义局部对比度为s(x)＝|λ₁-λ₂|，则融合结构张量与非局域的权函数w_STNL(x，y)定义为：

$w_{STNL}(x,y)=\exp (-\frac{G_{\mathrm{\σ}}*|P_x-P_y{|}^2}{h^2})\·\exp (-\frac{G_{\mathrm{\σ}}*|P_{s_x}-P_{s_y}{|}^2}{h^2})$

其中和分别代表块P_x和P_y的局部对比度值；

利用新的权函数w_STNL(x，y)修改正则项，得到融合结构变量与非局域全变分的目标函数

其中|▽_STNLu|称为融合结构变量与非局域信息的全变分正则项。

进一步，在步骤3)中，引入变量目标函数改写为：

$\begin{array}{cc}\underset{u,d}{\min }{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}|d|& s.t.d={\▿}_{STNL}u\end{array}$

采用Bregman迭代法对上式进行求解得到：

$\left\{\begin{array}{c}(u^{k+1},d^{k+1})={\mathrm{argmin}}_{u,d}{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}|d|+\frac{\mathrm{\β}}{2}|d-{\▿}_{STNL}u-b^k{|}^2\\ b^{k+1}=b^k+d^{k+1}-{\▿}_{STNL}{u}^{k+1}\end{array}\right.---\left(1\right)$

其中u^k+1，d^k+1，b^k+1分别为变量u，d，b在第k+1次迭代的值，b^k为变量b在第k次迭代的值，β为权重因子；则u^k+1的最小化方程由以下最优性条件给出：

(u^k+1-f)-βdiv_STNL(d^k-▽_STNLu^k+1-b^k)＝0

其中div为散度算子，d^k为变量d在第k次迭代的值；

采用Gauss-Seidel迭代法得到u^k+1在第i点的迭代方程：

$u_i^{k+1,n+1}=\frac{1}{1+\mathrm{\β}{\Σ}_j{w}_{ij}}(\mathrm{\β}{\Σ}_j{w}_{ij}{u}_i^{k+1,n}+f_i-\mathrm{\β}{\Σ}_j\sqrt{w_{ij}}(d_{ij}^{k+1,n}-d_{ji}^{k+1,n}-b_{ij}^{k+1,n}+b_{ji}^{k+1,n}\left)\right)---\left(2\right)$

u^k+1，n＝0＝u^k(3)

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环，上标n代表第n次内循环,n＝0代表迭代的初始化状态；

变量d的迭代方程由软阈值函数给出：

$d^{k+1}=\frac{{\▿}_{STNL}{u}^{k+1}+b^k}{|{\▿}_{STNL}{u}^{k+1}+b^k|}\max \left(\right|{\▿}_{STNL}{u}^{k+1}+b^k|-\mathrm{\λ}/\mathrm{\β},0)$

变量d的离散化表达式：

$d_{ij}^{k+1}=\frac{\sqrt{w_{ij}}(u_j^{k+1}-u_i^{k+1})+b_{ij}^k}{\sqrt{{\Σ}_j{w}_{ij}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{ij}^k}^2}}\max (\sqrt{{\Σ}_j{w}_{ij}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{ij}^k}^2}-\mathrm{\λ}/\mathrm{\β},0)---\left(4\right)$

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环；

由式(1)得到变量b的迭代方程为：

$b_{ij}^{k+1}=b_{ij}^k+d_{ij}^{k+1}-\sqrt{w_{ij}}(u_j^{k+1}-u_i^{k+1})---\left(5\right)$

联合式(2-5)，得到模型求解的变量迭代方程组；迭代停止后，去噪图像为u＝u^R，N，其中K为外循环最大迭代次数，N为内循环最大迭代次数。

本发明的有益效果在于：本发明在基于全变分的图像去噪模型中引入非局域描述子，可以较好捕捉图像的细小纹理，同时考虑图像的结构信息，提出一种融合结构张量与非局域全变分的图像去噪模型。该模型一方面考虑了图像的非局域特性，通过度量patch间的灰度相似性，能够较好的保持图像的纹理特性；另一方面兼顾了几何结构特性，能够调节不同对比度patch之间权重，尤其是提高了低对比度区域的权重值，同时较少阶梯效应的产生。模型求解采用分裂Bregman算法，克服了传统梯度下降法中存在的非线性和不可微性带来的最小化速度缓慢问题，本发明中采用的分裂Bregman算法快速有效。

附图说明

为了使本发明的目的、技术方案和有益效果更加清楚，本发明提供如下附图进行说明：

图1为本发明所述方法的流程示意图。

具体实施方式

下面将结合附图，对本发明的优选实施例进行详细的描述。

图1为本发明所述方法的流程示意图，如图所示，本方法包括以下步骤：1)对带噪图像f进行N*N的分块，N的取值可为N＝3,5,7...，在本实施例中为3-2M+1(M为自然数)的奇数；2)对图像去噪建模，建立融合结构张量与非局域全变分的目标函数；3)利用分裂Bregman算法求解目标函数得到去噪图像。

本实施例中，步骤2)中建立目标函数的详细步骤为：

21)建立基于全变分正则项的图像去噪模型

其中f是带噪图像，u是待恢复的去噪图像，图像去噪的目的就是希望从f中恢复出清晰的图像u。式中第一项为数据项，第二项为平滑项，λ是平衡数据项与平滑项的权重因子，Ω为整个图像空间，为n维实数空间，▽为梯度算子，|▽u|也称为全变分正则项；

22)将传统的全变分正则项扩展到非局域空间，得到非局域全变分并作为正则项约束，进而得到基于非局域信息的全变分图像去噪模型

其中非局域算子▽_NLu由图像空间中任意两个像素点(x，y)∈Ω×Ω定义为

${\▿}_{NL}u(x,y)=\left(u\right(y)-u(x\left)\right)\sqrt{w_{NL}(x,y)}$

其中u(x)，u(y)分别指u在像素点x，y处的值，w_NL：Ω×Ω为像素点x，y之间的边缘权重，且满足对称性，即w_NL(x，y)＝w_NL(y，x)。用于衡量块P_x和P_y之间的相似性,其中G_σ是标准偏差为σ的高斯卷积核，h为尺度参数。与传统的全变分方法相比，非局域正则化方法可以更好的刻画图像中的纹理特征。

23)融合图像的结构张量与非局域全变分。首先，定义结构张量矩阵为

$J_{\mathrm{\ρ}}(\▿I_{\mathrm{\σ}})=G_{\mathrm{\ρ}}*(\▿I_{\mathrm{\σ}}\⊗\▿I_{\mathrm{\σ}})\stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]$

其中，I_σ是经过方差为σ的高斯滤波器的图像I，目的是使边缘检测对尺度小于σ的图像内容不敏感，从而减小求导运算的影响。▽是梯度算子，是张量积运算，*代表卷积运算，G_ρ是方差为ρ的高斯滤波器，与张量积的卷积是为了增强对线型纹理、断裂边缘以及角型等几何结构的刻画能力。J₁₁，J₁₂，J₂₁，J₂₂分别为结构张量矩阵的四个元素。

求解对称半正定的结构张量矩阵J_ρ的两个正交特征向量

$v_1=\left[\begin{array}{c}2J_{12}\\ J_{22}-J_{11}+\sqrt{{(J_{22}-J_{11})}^2+4J_{12}^2}\end{array}\right],v_1=\frac{v_1}{\left|\right|v_1\left|\right|},v_2=v_1^{\⊥}$

与对应的特征值

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(J_{11}+J_{22}\±\sqrt{{(J_{22}-J_{11})}^2+4J_{12}^2})$

特征值描述的是邻域内的特征方向的平均变化和形状信息。当λ₁≈λ₂≈0，对应平坦区域；当λ₁>>λ₂≈0，对应边缘区域；当λ₁≥λ₂＞＞0，对应角型区域或T型结构。定义局部对比度为s(x)＝|λ₁-λ₂|，可以看出局部对比度较好的体现了图像的几何属性，尤其在含有噪声的低对比度区域。因此，融合结构张量与非局域的权函数w_STNL(x，y)定义为

$w_{STNL}(x,y)=\exp (-\frac{G_{\mathrm{\σ}}*|P_x-P_y{|}^2}{h^2})\·\exp (-\frac{G_{\mathrm{\σ}}*|P_{s_x}-P_{s_y}{|}^2}{h^2})$

其中和分别代表块P_x和P_y的局部对比度值。融合了图像非局域信息与结构张量的相似性权重，一方面度量整个空域内patch间的灰度相似性，能够较好的保持图像的纹理特性；另一方面考虑了几何结构特性，能够调节不同对比度patch之间权重，尤其是提高了低对比度区域的权重值，同时较少阶梯效应的产生。

利用新的权函数w_STNL(x，y)修改正则项，得到融合结构变量与非局域全变分的目标函数

其中|▽_STNLu|称为融合结构变量与非局域信息的全变分正则项，y是已知带噪图像，u是待求去噪图像。新的目标函数同时兼顾了图像中的非局域信息与结构信息，可以很好的保留纹理和结构特征，尤其在低对比度区域。

本实施例中步骤3)的详细步骤为：

引入变量目标函数改写为：

$\begin{array}{cc}\underset{u,d}{\min }{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}|d|& s.t.d={\▿}_{STNL}u\end{array}$

采用Bregman迭代法对上式进行求解得到：

$\left\{\begin{array}{c}(u^{k+1},d^{k+1})={\mathrm{argmin}}_{u,d}{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}|d|+\frac{\mathrm{\β}}{2}|d-{\▿}_{STNL}u-b^k{|}^2\\ b^{k+1}=b^k+d^{k+1}-{\▿}_{STNL}{u}^{k+1}\end{array}\right.---\left(1\right)$

其中u^k+1，d^k+1，b^k+1分别为变量u，d，b在第k+1次迭代的值，b^k为变量b在第k次迭代的值，β为权重因子。则u^k+1的最小化方程由以下最优性条件给出：

(u^k+1-f)-βdiv_STNL(d^k-▽_STNLu^k+1-b^k)＝0

其中div为散度算子，d^k为变量d在第k次迭代的值。

采用Gauss-Seidel迭代法得到u^k+1在第i点的迭代方程：

$u_i^{k+1,n+1}=\frac{1}{1+\mathrm{\β}{\Σ}_j{w}_{ij}}(\mathrm{\β}{\Σ}_j{w}_{ij}{u}_i^{k+1,n}+f_i-\mathrm{\β}{\Σ}_j\sqrt{w_{ij}}(d_{ij}^{k+1,n}-d_{ji}^{k+1,n}-b_{ij}^{k+1,n}+b_{ji}^{k+1,n}\left)\right)---\left(2\right)$

u^k+1，n＝0＝u^k(3)

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环，上标n代表第n次内循环,n＝0代表迭代的初始化状态。

变量d的迭代方程由软阈值函数给出：

$d^{k+1}=\frac{{\▿}_{STNL}{u}^{k+1}+b^k}{|{\▿}_{STNL}{u}^{k+1}+b^k|}\max \left(\right|{\▿}_{STNL}{u}^{k+1}+b^k|-\mathrm{\λ}/\mathrm{\β},0)$

变量d的离散化表达式为：

$d_{ij}^{k+1}=\frac{\sqrt{w_{ij}}(u_j^{k+1}-u_i^{k+1})+b_{ij}^k}{\sqrt{{\Σ}_j{w}_{ij}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{ij}^k}^2}}\max (\sqrt{{\Σ}_j{w}_{ij}{(u_j^{k+1}-u_i^{k+1})}^2+{b_{ij}^k}^2}-\mathrm{\λ}/\mathrm{\β},0)---\left(4\right)$

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环。

由式(1)得到变量b的迭代方程为：

$b_{ij}^{k+1}=b_{ij}^k+d_{ij}^{k+1}-\sqrt{w_{ij}}(u_j^{k+1}-u_i^{k+1})---\left(5\right)$

联合式(2-5)，得到模型求解的变量迭代方程组。迭代停止后，去噪图像为u＝u^K，N，其中K为外循环最大迭代次数，N为内循环最大迭代次数。

最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明权利要求书所限定的范围。

Claims (4)

1.一种融合结构张量与非局域全变分的图像去噪方法，其特征在于：包括以下步骤：

1)对带噪图像f进行N*N的分块，N的取值可为N＝3,5,7...；

2)对图像去噪建模，建立融合结构张量与非局域全变分的目标函数；

3)利用分裂Bregman算法求解目标函数得到去噪图像。

2.根据权利要求1所述的一种融合结构张量与非局域全变分的图像去噪方法，其特征在于：在步骤1)中，N的取值可为3-2M+1(M为自然数)的奇数。

3.根据权利要求1所述的一种融合结构张量与非局域全变分的图像去噪方法，其特征在于：在步骤2)中，具体包括：

21)建立基于全变分正则项的图像去噪模型

其中f是带噪图像，u是待恢复的去噪图像，图像去噪的目的在于从f中恢复出清晰的图像u；式中第一项为数据项，第二项为平滑项，λ是平衡数据项与平滑项的权重因子，Ω为整个图像空间，为n维实数空间，为梯度算子，也称为全变分正则项；

22)将传统的全变分正则项扩展到非局域空间，得到非局域全变分并作为正则项约束，进而得到基于非局域信息的全变分图像去噪模型

其中非局域算子由图像空间中任意两个像素点(x，y)∈Ω×Ω定义为：

${\▿}_{NL}u(x,y)=\left(u\right(y)-u(x\left)\right)\sqrt{w_{NL}(x,y)}$

其中u(x)，u(y)分别指u在像素点x，y处的值，w_NL：Ω×Ω为像素点x，y之间的边缘权重，且满足对称性，即w_NL(x，y)＝w_NL(y，x)；用于衡量块P_x和P_y之间的相似性,其中G_σ是标准偏差为σ的高斯卷积核，h为尺度参数；

23)融合图像的结构张量与非局域全变分：首先，定义结构张量矩阵为：

$J_{\mathrm{\ρ}}(\▿I_{\mathrm{\σ}})=G_{\mathrm{\ρ}}*(\▿I_{\mathrm{\σ}}\⊗\▿I_{\mathrm{\σ}})\stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{cc}{J}_{11}& J_{12}\\ J_{21}& J_{22}\end{array}\right]$

其中，I_σ是经过方差为σ的高斯滤波器的图像I，是梯度算子，是张量积运算，*代表卷积运算，G_ρ是方差为ρ的高斯滤波器，J₁₁，J₁₂，J₂₁，J₂₂分别为结构张量矩阵的四个元素；

求解对称半正定的结构张量矩阵J_ρ的两个正交特征向量：

$v_1=\left[\begin{array}{c}2J_{12}\\ J_{22}-J_{11}+\sqrt{{\left(J_{22}-J_{11}\right)}^2+4J_{12}^2}\end{array}\right],v_1=\frac{v_1}{\left|\right|v_1\left|\right|},v_2=v_1^{\⊥}$

与对应的特征值：

${\mathrm{\λ}}_{1,2}=\frac{1}{2}(J_{11}+J_{22}\±\sqrt{{(J_{22}-J_{11})}^2+4J_{12}^2})$

定义局部对比度为s(x)＝|λ₁-λ₂|，则融合结构张量与非局域的权函数w_STNL(x，y)定义为：

$w_{STNL}(x,y)=\exp (-\frac{G_{\mathrm{\σ}}*|P_x-P_y{|}^2}{h^2})\·\exp (-\frac{G_{\mathrm{\σ}}*|P_{s_x}-P_{s_y}{|}^2}{h^2})$

其中和分别代表块P_x和P_y的局部对比度值；

利用新的权函数w_STNL(x，y)修改正则项，得到融合结构变量与非局域全变分的目标函数

其中称为融合结构变量与非局域信息的全变分正则项。

4.根据权利要求3所述的一种融合结构张量与非局域全变分的图像去噪方法，其特征在于：在步骤3)中，引入变量目标函数改写为：

$\begin{array}{ccc}\underset{u,d}{\min }{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}\left|d\right|& s.t.& d={\▿}_{STNL}u\end{array}$

采用Bregman迭代法对上式进行求解得到：

$\left\{\begin{array}{c}\left(u^{k+1},d^{k+1}\right)={\mathrm{argmin}}_{u,d}{\∫}_{\mathrm{\Ω}}\frac{1}{2}\left|\right|u-f|{|}^2+\mathrm{\λ}\left|d\right|+\frac{\mathrm{\β}}{2}|d-{\▿}_{STNL}u-b^k{|}^2\\ b^{k+1}=b^k+d^{k+1}-{\▿}_{STNL}{u}^{k+1}\end{array}\right.---\left(1\right)$

其中u^k+1，d^k+1，b^k+1分别为变量u，d，b在第k+1次迭代的值，b^k为变量b在第k次迭代的值，β为权重因子；则u^k+1的最小化方程由以下最优性条件给出：

$(u^{k+1}-f)-{\mathrm{\βdiv}}_{STNL}(d^k-{\▿}_{STNL}{u}^{k+1}-b^k)=0$

其中div为散度算子，d^k为变量d在第k次迭代的值；

采用Gauss-Seidel迭代法得到u^k+1在第i点的迭代方程：

$u_i^{k+1,n+1}=\frac{1}{1+{\mathrm{\β\Σ}}_j{w}_{ij}}\left({\mathrm{\β\Σ}}_j{w}_{ij}{u}_i^{k+1,n}+f_i-{\mathrm{\β\Σ}}_j\sqrt{w_{ij}}\left(d_{ij}^{k+1,n}-d_{ji}^{k+1,n}-b_{ij}^{k+1,n}+b_{ji}^{k+1,n}\right)\right)---\left(2\right)$

u^k+1，n＝0＝u^k(3)

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环，上标n代表第n次内循环,n＝0代表迭代的初始化状态；

变量d的迭代方程由软阈值函数给出：

$d^{k+1}=\frac{{\▿}_{STNL}{u}^{k+1}+b^k}{|{\▿}_{STNL}{u}^{k+1}+b^k|}\max \left(|{\▿}_{STNL}{u}^{k+1}+b^k|-\mathrm{\λ}/\mathrm{\β},0\right)$

变量d的离散化表达式：

$d_{ij}^{k+1}=\frac{\sqrt{w_{ij}}\left(u_j^{k+1}-u_i^{k+1}\right)+b_{ij}^k}{\sqrt{{\mathrm{\Σ}}_j{w}_{ij}{\left(u_j^{k+1}-u_i^{k+1}\right)}^2+{b_{ij}^k}^2}}\max \left(\sqrt{{\mathrm{\Σ}}_j{w}_{ij}{\left(u_j^{k+1}-u_i^{k+1}\right)}^2+{b_{ij}^k}^2}-\mathrm{\λ}/\mathrm{\β},0\right)---\left(4\right)$

其中下标i，j分别代表第i，j个像素点，上标k代表第k次外循环；

由式(1)得到变量b的迭代方程为：

$b_{ij}^{k+1}=b_{ij}^k+d_{ij}^{k+1}-\sqrt{w_{ij}}(u_j^{k+1}-u_i^{k+1})---\left(5\right)$

联合式(2-5)，得到模型求解的变量迭代方程组；迭代停止后，去噪图像为u＝u^K，N，其中K为外循环最大迭代次数，N为内循环最大迭代次数。