图像成像与传输的过程会受到噪声干扰，导致图像信息的可靠性下降，进而影响到高层次图像信息的处理。对于具有独立同分布统计特性的加性高斯白噪声$ \mathit{\boldsymbol{\eta }} $，灰度图像退化模型可以表示为

$ {\mathit{\boldsymbol{u}}_0} = \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{\eta }} $

(1)

图像去噪的目标是衰减退化图像$ {\mathit{\boldsymbol{u}}_0} $中的噪声$ \mathit{\boldsymbol{\eta }} $，并保护真实图像$ \mathit{\boldsymbol{u}} $的结构。近年来，基于偏微分方程(PDE)的图像去噪方法由于其高度的灵活性和更强的局部自适应性吸引了大量关注。它通过极小化某一表征图像光滑程度的能量泛函$ E\left( \mathit{\boldsymbol{u}} \right) $得到相应的自适应扩散PDE，通过控制扩散的方向与强度使图像的结构在平滑噪声的过程中得以保持。

彩色图像去噪模型大部分源于相关的灰度图像去噪模型。Perona等人关于各向异性扩散的研究在该领域产生了深远影响^[1]。提出的非线性扩散模型为

$ \begin{array}{*{20}{c}} {\frac{{\partial u\left( {x,y} \right)}}{{\partial t}} = {\rm{div(}}g\left( {\left\| {\nabla u\left( {x,y} \right)} \right\|} \right) \times }\\ {\nabla u\left( {x,y} \right))} \end{array} $

(2)

使用梯度模$ \left\| {\nabla u\left( {x, y} \right)} \right\| $表征图像的局部结构，设计扩散系数$ g'\left( {\left\| {\nabla u\left( {x, y} \right)} \right\|} \right) < 0 $来减缓边界处的扩散速率，加快平坦区域内部的平滑。一般地，式(2) 对应的能量泛函模型为

$ E\left( {u\left( {x,y} \right)} \right) = \int_\Omega {f\left( {\left\| {\nabla u\left( {x,y} \right)} \right\|} \right)} {\rm{d}}x{\rm{d}}y $

(3)

通常设计$ f\left( \cdot \right) \ge 0 $为关于$ \left\| {\nabla u\left( {x, y} \right)} \right\| $的增函数，极小化上述能量泛函就相当于光滑了图像。此外，还需要考虑设计扩散系数$ g\left( {\left\| {\nabla u\left( {x, y} \right)} \right\|} \right) = f'\left( {\left\| {\nabla u\left( {x, y} \right)} \right\|} \right)/\left\| {\nabla u\left( {x, y} \right)} \right\| $来减缓边缘处的扩散速率^[2-3]。Rudin等人阐释了能量泛函与PDE扩散方程的紧密联系，提出了著名的ROF模型^[4]

$ \begin{array}{*{20}{c}} {\int_\Omega {\left\| {\nabla u\left( {x,y} \right)} \right\|{\rm{d}}x{\rm{d}}y + } }\\ {\frac{\lambda }{2}\int_\Omega {(u\left( {x,y} \right) - {u_0}{{\left( {x,y} \right)}^2}{\rm{d}}x{\rm{d}}y} } \end{array} $

(4)

使用图像$ u $的全变分(TV)作为正则项来衡量图像的光滑程度，极小化TV则能够在保持图像结构的同时衰减噪声。

上述模型需要求解一个二阶PDE，其结果往往存在“阶梯效应”，会在光滑渐变区域产生虚假边界。Lysaker等人使用二阶导数来度量图像的光滑程度，提出了LLT模型^[5]

$ \int_{\Omega }{\sqrt{\begin{array}{*{35}{l}} u_{xx}^{2}\left( x,y \right)+u_{xy}^{2}\left( x,y \right)+ \\ u_{yx}^{2}\left( x,y \right)+u_{yy}^{2}\left( x,y \right) \\ \end{array}}}\text{d}x\text{d}y $

(5)

极小化上述模型得到分片线性的稳态解，避免了“阶梯效应”的产生。Zhu等人^[6]使用灰度图像曲面的平均曲率(MC)来表征图像的局部结构。图像的几何特征比如角点、对比度等得到了较好的保护，并且也不产生“阶梯效应”。

彩色图像$ \mathit{\boldsymbol{u}}\left( {x, y} \right) = \left( {{\mathit{\boldsymbol{u}}_R}\left( {x, y} \right), {\mathit{\boldsymbol{u}}_G}\left( {x, y} \right), {\mathit{\boldsymbol{u}}_B}\left( {x, y} \right)} \right) $去噪时需要考虑各个单色通道的处理方式。最直接的处理方式是利用RGB颜色模型各个通道相互独立的性质，将每个单色通道视为灰度图像逐个通道进行处理后合成为一幅彩图(CBC-L1)，即考虑的能量泛函为$ \sum\limits_{l = 1}^3 {E\left( {{\mathit{\boldsymbol{u}}_l}} \right)} $，$ l=R $，$ G $，$ B $，$ {\mathit{\boldsymbol{u}}_l} $表示某个单色通道且$ E\left( {{\mathit{\boldsymbol{u}}_l}} \right) \ge 0 $。但RGB颜色模型所描述的图像特征在通道之间具有较强的相关性，孤立地处理各个通道难以保持彩色图像的结构特征。Blomgren和Chan将ROF模型通过CBC-L2的方式耦合起来处理彩色图像^[7]，考虑的能量泛函为$ \sqrt {{{\sum\limits_{l = 1}^3 {E\left( {{\mathit{\boldsymbol{u}}_l}} \right)} }^2}} $。Brito-Loeza等人以类似的方法推广了MC模型，得到了基于全局曲率的模型(GCM)^[8]。但这种全局的耦合方式在通道内部给出了相同的自适应系数，无法刻画彩色图像的局部结构。Moreno等人^[9]使用各通道的梯度信息构造了一个交叉耦合项来处理彩色图像。Ono等人^[10]将RGB图像转换到对色颜色空间中处理，以消除通道间结构信息的相关性，在此基础上提出了去相关的彩色全变分模型(DVTV)，他们还使用局部颜色矩阵的核范数抑制颜色失真^[11]。

黎曼几何也是处理彩图的一种有效方法。Sapiro^[12]在黎曼几何框架下定义了一族彩色图像全变分$ \int_\Omega {f\left( {{\lambda _ + }, {\lambda _-}} \right)} {\rm{d}}\Omega $，其中$ {{\lambda _ + }} $、$ {{\lambda _-}} $分别是图像某点度量张量的最大、最小特征值，选择$ f\left( {{\lambda _ + }, {\lambda _-}} \right) = \sqrt {{\lambda _ + }+{\lambda _-}} $时得到TVF模型。Goldluecke等人^[13]提出使用谱范数$ f\left( {{\lambda _ + }, {\lambda _-}} \right) = \sqrt {{\lambda _ + }} $来度量图像的光滑程度。Holt^[14]则提出极小化图像雅可比矩阵的核范数。Lefkimmiatis等人基于块的雅可比矩阵的Schatten范数归纳了一些彩色图像的TV^[15]。但是目前该框架下的几何模型都是低阶的，会产生灰度图像类似的“阶梯效应”。

本文在黎曼几何的框架下，提出了新的高阶彩色图像去噪模型。首先将几何框架下的低阶模型进行推广，得到高阶矩阵，然后利用其Frobenius范数构造变分模型。为了提高边界保持能力，使用单通道内梯度卷积作为加权函数引导高阶模型的扩散。新模型利用单通道和多通道的一阶和二阶信息进行耦合去噪，既保留了高阶模型的平滑特性，也充分体现了低阶模型的边界保持能力。

在黎曼几何框架下，彩色图像$ \mathit{\boldsymbol{u}}\rm{:}\Omega \subset {{\mathbb{R}}^{2}}\to {{\mathbb{R}}^{3}} $每一空间位置对应矢量终点构成了$ {{\mathbb{R}}^{3}} $空间中的曲面。取曲面$ \mathit{\boldsymbol{u}} $上某点$ P $，其度量张量$ \mathit{\boldsymbol{G}} $的分量可以定义为$ {{g}_{i, j}}={{\mathit{\boldsymbol{g}}}_{i}}\cdot {{\mathit{\boldsymbol{g}}}_{j}}, i, j=1, 2 $，其中$ {{\mathit{\boldsymbol{g}}}_{i}} $和$ {{\mathit{\boldsymbol{g}}}_{j}} $是该点切平面$ {{\mathit{\boldsymbol{T}}}_{p}} $的自然基矢量，即$ {{\mathit{\boldsymbol{g}}}_{1}}=\frac{\partial \mathit{\boldsymbol{u}}}{\partial x}, {{\mathit{\boldsymbol{g}}}_{2}}=\frac{\partial \mathit{\boldsymbol{u}}}{\partial y} $。

$ P $点切平面上的微元可表示为

$ {\rm{d}}\mathit{\boldsymbol{u}}{\rm{ = }}{\mathit{\boldsymbol{g}}_{\rm{1}}}{\rm{d}}x + {\mathit{\boldsymbol{g}}_2}{\rm{d}}y $

(6)

其欧几里得范数的平方

$ \begin{array}{l} {\rm{d}}\mathit{\boldsymbol{u}} \cdot {\rm{d}}\mathit{\boldsymbol{u = }}{\mathit{\boldsymbol{g}}_{\rm{1}}} \cdot {\mathit{\boldsymbol{g}}_1}{\rm{d}}{x^2} + 2{\mathit{\boldsymbol{g}}_1} \cdot {\mathit{\boldsymbol{g}}_2}{\rm{d}}x{\rm{d}}y + \\ \;\;\;\;\;\;{\mathit{\boldsymbol{g}}_2} \cdot {\mathit{\boldsymbol{g}}_2}{\rm{d}}{y^2} = \left[ {{\rm{d}}x,{\rm{d}}y} \right]\mathit{\boldsymbol{G}}{\left[ {{\rm{d}}x,{\rm{d}}y} \right]^{\rm{T}}} \end{array} $

(7)

称为曲面$ \mathit{\boldsymbol{u}} $的第一基本形式。度量张量$ \mathit{\boldsymbol{G}}={{\left( \rm{D}\mathit{\boldsymbol{u}} \right)}^{\rm{T}}}\rm{D}\mathit{\boldsymbol{u}}={{\left[{{\mathit{\boldsymbol{g}}}_{1}}, {{\mathit{\boldsymbol{g}}}_{2}} \right]}^{\rm{T}}}\left[{{\mathit{\boldsymbol{g}}}_{1}}, {{\mathit{\boldsymbol{g}}}_{2}} \right] $包含了点$ P $沿切平面各方向的变化信息。一阶导数矩阵(雅可比矩阵) $ \rm{D}\mathit{\boldsymbol{u}} $的Frobenius范数为

$ \begin{array}{l} {\left\| {\text{D}\mathit{\boldsymbol{u}}} \right\|_{\rm{F}}} = \sqrt {\sum\limits_{i = 1}^2 {{{\left\| {{\mathit{\boldsymbol{g}}_i}} \right\|}^2}} } = \sqrt {{\rm{tr}}\left( \mathit{\boldsymbol{G}} \right)} = \\ \;\;\;\;\;\;\;\;\;\;\;\sqrt {\sigma _1^2 + \sigma _2^2} = \sqrt {{\lambda _ + } + {\lambda _ - }} \end{array} $

(8)

式中，$ {{\sigma }_{1}} $、$ {{\sigma }_{2}} $是$ \rm{D}\mathit{\boldsymbol{u}} $的奇异值；$ {{\lambda }_{+}}=\sigma _{1}^{2} $、$ {{\lambda }_{-}}=\sigma _{2}^{2} $是$ \mathit{\boldsymbol{G}} $的特征值，反映了该点最大、小的变化；$ \rm{tr}\left( \mathit{\boldsymbol{G}} \right) $表示矩阵$ \mathit{\boldsymbol{G}} $的迹。构造能量泛函$ E\left( \mathit{\boldsymbol{u}} \right)={{\int_{\Omega }{\left\| \rm{D}\mathit{\boldsymbol{u}} \right\|}}_{\rm{F}}}\rm{d}\Omega $，其对应的PDE为

$ \frac{{\partial {u_l}\left( {x,y} \right)}}{{\partial t}} = {\rm{div}}\left( {\frac{{\left\| {\nabla {u_l}\left( {x,y} \right)} \right\|}}{{\sqrt {\sum\limits_{l = 1}^3 {{{\left\| {\nabla {u_l}\left( {x,y} \right)} \right\|}^2}} } }}\frac{{\nabla {u_l}\left( {x,y} \right)}}{{\left\| {\nabla {u_l}\left( {x,y} \right)} \right\|}}} \right) $

(9)

记该模型为TVF，在处理灰度图像时退化为TV。Bresson和Chan为此模型提出了一个有效的数值算法^[16]。

在黎曼几何的框架下，现有模型都是基于$ \mathit{\boldsymbol{G}} $的各种矩阵范数构建的，这些低阶模型会产生“阶梯效应”。为了建立高阶模型，对式(6) 进行微分得到

$ \begin{array}{l} {{\rm{d}}^2}\mathit{\boldsymbol{u}}{\rm{ = }}{\mathit{\boldsymbol{u}}_{{\rm{11}}}}{\rm{d}}{x^2} + {\mathit{\boldsymbol{u}}_{12}}{\rm{d}}x{\rm{d}}y + \\ \;\;\;\;\;\;{\mathit{\boldsymbol{u}}_{21}}{\rm{d}}y{\rm{d}}x + {\mathit{\boldsymbol{u}}_{22}}{\rm{d}}{y^2} \end{array} $

(10)

式中，$ {{\mathit{\boldsymbol{u}}}_{11}}=\frac{{{\partial }^{2}}\mathit{\boldsymbol{u}}}{\partial {{x}^{2}}} $，$ {{\mathit{\boldsymbol{u}}}_{12}}=\frac{{{\partial }^{2}}\mathit{\boldsymbol{u}}}{\partial x\partial y} $，$ {{\mathit{\boldsymbol{u}}}_{21}}=\frac{{{\partial }^{2}}\mathit{\boldsymbol{u}}}{\partial y\partial x} $，$ {{\mathit{\boldsymbol{u}}}_{22}}=\frac{{{\partial }^{2}}\mathit{\boldsymbol{u}}}{\partial {{y}^{2}}} $。

式(10) 欧几里得范数的平方为

$ \begin{array}{*{20}{c}} {{{\rm{d}}^2}\mathit{\boldsymbol{u}} \cdot {{\rm{d}}^2}\mathit{\boldsymbol{u}}{\rm{ = }}}\\ {\left[ {{\rm{d}}{x^2},{\rm{d}}x{\rm{d}}y,{\rm{d}}y{\rm{d}}x,{\rm{d}}{y^2}} \right]\mathit{\boldsymbol{T}}{{\left[ {{\rm{d}}{x^2},{\rm{d}}x{\rm{d}}y,{\rm{d}}y{\rm{d}}x,{\rm{d}}{y^2}} \right]}^{\rm{T}}}} \end{array} $

(11)

式中，$ \mathit{\boldsymbol{T}}={{\left[{{\mathit{\boldsymbol{u}}}_{11}}, {{\mathit{\boldsymbol{u}}}_{12}}, {{\mathit{\boldsymbol{u}}}_{21}}, {{\mathit{\boldsymbol{u}}}_{22}} \right]}^{\rm{T}}}\left[{{\mathit{\boldsymbol{u}}}_{11}}, {{\mathit{\boldsymbol{u}}}_{12}}, {{\mathit{\boldsymbol{u}}}_{21}}, {{\mathit{\boldsymbol{u}}}_{22}} \right] $，记$ {{\rm{D}}^{\rm{2}}}\mathit{\boldsymbol{u}}=\left[{{\mathit{\boldsymbol{u}}}_{11}}, {{\mathit{\boldsymbol{u}}}_{12}}, {{\mathit{\boldsymbol{u}}}_{21}}, {{\mathit{\boldsymbol{u}}}_{22}} \right] $，于是可以得到二阶导数矩阵$ {{\rm{D}}^{\rm{2}}}\mathit{\boldsymbol{u}} $的Frobenius范数(以下简称F范数)

$ \begin{array}{l} {\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|_{\rm{F}}} = \sqrt {\sum\limits_{i = 1}^2 {\sum\limits_{j = 1}^2 {{{\left\| {{\mathit{\boldsymbol{u}}_{i,j}}} \right\|}^2}} } } = \\ \;\;\sqrt {{\rm{tr}}\left( \mathit{\boldsymbol{T}} \right)} = \sqrt {\sigma _1^2 + \sigma _2^2 + \sigma _3^2} \end{array} $

(12)

式中，$ {{\sigma }_{1}} $、$ {{\sigma }_{2}} $、$ {{\sigma }_{3}} $是$ {{\rm{D}}^{\rm{2}}}\mathit{\boldsymbol{u}} $的奇异值。构造能量泛函$ E\left( \mathit{\boldsymbol{u}} \right) = {\int_\Omega {\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|} _{\rm{F}}}{\rm{d}}\Omega $，其对应的PDE为

$ \begin{array}{l} \frac{{\partial {u_l}\left( {x,y} \right)}}{{\partial t}} = - ({\left( {\frac{{{u_{l,xx}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{xx}} + {\left( {\frac{{{u_{l,xy}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{xy}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\left( {\frac{{{u_{l,yx}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{yx}} + {\left( {\frac{{{u_{l,yy}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{yy}}) \end{array} $

(13)

式中，$ {u_{l, xx}}\left( {x, y} \right) $、$ {u_{l, xy}}\left( {x, y} \right) $和$ {u_{l, yy}}\left( {x, y} \right) $分别表示$ {u_l} $关于下标的二阶导数。上述模型的扩散与F范数$ {\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|_{\rm{F}}} $相关，F范数包含的二阶导数信息体现了各通道之间的耦合关系。

需要注意的是当输入为灰度图像时，式(13) 退化为标准的LLT模型。但在灰度图像去噪中，LLT模型对边界信息的保持效果不佳，原因在于二阶导数受噪声干扰大，无法提供可靠的边界信息。因此可以预见，使用式(13) 对彩色图像去噪时，边界效果效果可能会有所欠缺。

从引导的观点看，全变分等低阶模型中使用了梯度信息引导扩散，因此避免了过度光滑边界^[17-18]。为此，本文引入加权函数$ \omega \left( {x, y} \right) $控制扩散速率以引导高阶扩散过程，得到如下高阶彩色图像去噪模型

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\partial {u_l}\left( {x,y} \right)}}{{\partial t}} = \\ - ({\left( {\frac{{\omega \left( {x,y} \right){u_{l,xx}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{xx}} + {\left( {\frac{{\omega \left( {x,y} \right){u_{l,xy}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{xy}} + \\ \;\;\;\;{\left( {\frac{{\omega \left( {x,y} \right){u_{l,yx}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{yx}} + {\left( {\frac{{\omega \left( {x,y} \right){u_{l,yy}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}}} \right)_{yy}}) \end{array} $

(14)

式中，$ \omega \left( {x, y} \right) = \omega \left( {\left\| {\nabla {u_{l\sigma }}\left( {x, y} \right)} \right\|} \right) $，$ {u_{l\sigma }}\left( {x, y} \right) = {G_\sigma }\left( {x, y} \right)*{u_l}\left( {x, y} \right) $表示高斯卷积后$ l $通道图像，使用高斯卷积是为了消除噪声对梯度可靠性的影响。$ \sigma $为高斯函数$ {G_\sigma }\left( {x, y} \right) = \frac{1}{{2{\rm{\pi }}{\sigma ^2}}}{{\rm{e}}^{{\rm{-}}\frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}}} $的标准差。

不失一般性，以式(14) 等号右侧第1项为例对模型展开分析。第1项的扩散核可以重写为

$ \begin{array}{l} \frac{{\omega \left( {\left\| {\nabla {u_{l\sigma }}\left( {x,y} \right)} \right\|} \right){u_{l,xx}}\left( {x,y} \right)}}{{{{\left\| {{{\rm{D}}^2}\mathit{\boldsymbol{u}}} \right\|}_{\rm{F}}}}} = \\ \;\;\;\;\omega \left( {\left\| {\nabla {u_{l\sigma }}\left( {x,y} \right)} \right\|} \right){c_l}{d_{l,xx}} \end{array} $

(15)

式中，$ {c_l} = {\left\| {{D^2}{u_l}} \right\|_{\rm{F}}}/{\left\| {{D^2}\mathit{\boldsymbol{u}}} \right\|_{\rm{F}}} $，$ {d_{l, xx}} = {u_{l, xx}}/{\left\| {{{\rm{D}}^2}{u_l}} \right\|_{\rm{F}}} $，$ {\left\| {{{\rm{D}}^2}{u_l}} \right\|_{\rm{F}}} = \sqrt {u_{l, xx}^2\left( {x, y} \right) + u_{l, xy}^2\left( {x, y} \right) + u_{l, yx}^2\left( {x, y} \right) + u_{l, yy}^2\left( {x, y} \right)} $。式(15) 刻画了多个方面对扩散的影响。

1) 单通道内一阶导数信息对扩散的引导作用：$ \omega \left( {\left\| {\nabla {u_{l\sigma }}\left( {x, y} \right)} \right\|} \right) $使用本通道的一阶梯度信息，描述本通道局部边界结构对本通道扩散的影响。为了在扩散时保持边界，通常设计为关于梯度的单调递减函数。

2) 多通道二阶导数信息对单通道扩散的耦合作用：$ {c_l} \le 1 $是单通道二阶导数矩阵F范数与多通道二阶导数矩阵的F范数之比，体现了其他通道信息与本通道扩散的耦合作用。

3) 单通道内部二阶导数信息对某个方向扩散的影响：$ {d_{l, xx}} \ge 1 $是水平方向二阶导数与4个二阶导数之和的比值。线性区域此项取值为0，在此方向上不发生扩散，可以避免“阶梯效应”的产生，体现了高阶模型的特性。

使用$ {u_{l\sigma }}\left( {x, y} \right) $的梯度模来表征图像边界，利用低阶模型保持边界的优势来更好地刻画通道内的结构信息。由于$ \left\| {\nabla {u_{l\sigma }}\left( {x, y} \right)} \right\| $在图像的边缘处很大，而噪声处则相对较小，因此设计加权函数$ \omega \left( {x, y} \right) $关于$ \left\| {\nabla {u_{l\sigma }}\left( {x, y} \right)} \right\| $单调递减就可以减缓边界处的扩散速率。受ROF模型的启发，采用$ \omega \left( {x, y} \right) = 1/\left( {1 + \left\| {\nabla {u_{l\sigma }}\left( {x, y} \right)} \right\|} \right) $，分母加1主要是为了保证算法稳定。

采用峰值信噪比(PSNR)与结构相似性(SSIM)作为评价图像的客观指标。由于人眼视觉对于误差的敏感度并不是绝对的，导致PSNR和SSIM等客观指标和人眼看到的视觉品质并不完全一致，因此在部分实验中选用放大局部图像、展示残差图像等方法对比不同算法的去噪结果差异。对模型式(14) 使用标准的有限差分法进行离散和计算^[1]，使用图 2所示图像集中的6幅图像进行仿真实验，并与文献[8]中以CBC-L2方式耦合的GCM模型，文献[10]中使用对色颜色空间处理的DVTV模型以及文献[16]中TVF模型相比较。实验选取高斯核函数的标准差$ \sigma = 0.1 \sim 0.3 $。

合成图 1只在水平方向有颜色变化，因此只考虑水平剖面的去噪并在图 3给出去噪结果，其中输入图像噪声的标准差$ {\sigma _\eta } = 15 $。图 3(c)中噪声没有被充分平滑，信号的下界小于0，而如果继续演化则边界难以保持造成PSNR下降。图 3(d)和图 3(e)得到的去噪后的信号边缘得到较好保持，但是线性和二次渐变区域的信号并不平滑，出现了不同程度的“阶梯效应”。图 3(f)是本文模型的去噪结果，渐变区域相对平滑，不产生“阶梯效应”，并且很好地保持了信号边界。

图 4给出了合成图像1的水平扫描线在RGB颜色空间中的分布情况。在图 4中，3维坐标系的3个轴分别表示红色、绿色、蓝色的强度。RGB颜色空间中的一个点代表一种颜色。从图 4(a)可以看出，原始信号在RGB颜色空间中映射为点和连续的直线或曲线。原始信号中的平坦区域颜色恒定，因此被映射为RGB颜色空间中的同一个点。RGB空间中的两个孤立点代表两个颜色不同的平坦区域。信号中的颜色渐变在RGB颜色空间中表现为不间断的直线或曲线，其中直线对应信号的线性渐变，曲线对应原始信号中更加光滑的部分。从图 4(b)可以看出，加入噪声后，信号在RGB颜色空间中的分布杂乱无章。图 4(c)—(e)中较好地恢复了孤立点，但无法在RGB颜色空间中得到连续的直线或曲线，表明其去噪结果中存在较为明显的“阶梯效应”。本文模型去噪结果对应的图 4(f)中，孤立点清晰，曲线光滑，表明该模型对颜色恒定区域和渐变区域都进行了良好的恢复。

图 5给出合成图 2的去噪结果与残差，噪声标准差$ {\sigma _\eta } = 15 $。合成图 2是分片常数图像，含有多个方向的边界。图 5(b)中能够观察到大量污点，视觉效果欠佳。图 5(c)—(e)整体效果相近，其中图 5(c)(d)的边界比较粗糙，而本文模型的则更加整齐。图 5(f)—(j)是相应的残差图像，从图 5(g)(i)可以观察到显著的结构信息，图 5(h)中的结构信息不够明显，而本文模型的残差则几乎看不到。这表明本文模型能够有效地保护图像的结构。

图 6给出使用Canny算子对合成图 2的去噪结果进行边缘检测的情况。图 6(c)(e)的平坦区域比图 6(d)(f)的平坦区域存在更多的边界响应，说明这些区域不够平滑，GCM形成了视觉上的斑点而DVTV则整体比较粗糙；图 6(d)与图 6(f)的差别不大，但是结合残差图 5(h)(j)来看，本文模型能够在充分平滑噪声的同时更好地保持边界。

图 7给出彩色Peppers图像的去噪结果，右上和左下是局部放大图像，噪声标准差$ {\sigma _\eta } = 30 $。除了图 7(c)所显示的辣椒表面颜色不均匀，形成了视觉上的斑点外，图 7(d)—(f)均具有较好的视觉效果。然而通过左下角的局部放大图可以看出，图 7(d)的辣椒表面存在明显的“阶梯效应”，图 7(e)也相对粗糙。本文模型给出的图 7(f)中的辣椒表面则更加平滑。观察右上角的局部放大图像可以发现本文模型的边界也得到了较好的保护。表明本文模型能够有效地平滑噪声，消除“阶梯效应”并且保持图像结构。

图 8—图 10分别给出People、Monarch、Fish图像的去噪结果以及相应的局部放大图像，噪声标准差分别为20、25、30。图 8(e)人脸比较粗糙，但是另一方面细节得以保持；嘴唇的颜色明显变暗，说明颜色的对比度有所下降，这一现象在图 9(e)中也有所体现。图 10(d)(e)鱼头存在明显的“阶梯效应”。而本文模型能够很好地保持颜色边界，有效地平滑噪声并且消除“阶梯效应”。

表 1给出了各模型对测试图像集中各图像在不同噪声水平下去噪结果的PSNR与SSIM。图 11给出了各模型在不同噪声水平下PSNR和SSIM增加的均值。只有在噪声标准差$ {\sigma _\eta } = 30 $的情况下，本文模型SSIM的平均增量稍低于TVF，原因在于此时图像的分片线性结构受噪声影响较大。而在其他情况下，本文模型的PSNR与SSIM的平均增量都具有优势。

本文提出一种梯度引导的高阶彩色图像去噪模型。首先将TVF模型推广到高阶，然后利用图像的一阶梯度信息构造加权函数引导高阶扩散，得到新的高阶彩色图像去噪模型。模型结合了高、低阶模型的优点，利用单通道和多通道的信息进行耦合去噪。实验效果表明，与GCM、TVF、DVTV等相关模型相比，本文模型不仅能够有效地消除噪声，还可以较好地保持图像边界，并且有效消除“阶梯效应”。

本文方法适用于去除彩色图像中不同噪声水平的高斯白噪声，与相关模型相比，在分片线性区域恢复方面具有一定的优势。高阶模型的复杂结构影响了计算效率，因此未来工作将关注模型快速求解和高效计算方面的研究。