图像修复是图像处理的一个基本问题，针对图像中遗失或破损的区域，利用周围已知的信息，按照一定的规则来填补，使修复后的图像接近或达到原图的视觉效果。修复技术在文物复原和文物保护方面有着重要和广泛的应用。如从旧字画中移除划痕、复原破损文物等。图像修复的研究方法根据修复的侧重点可以分为3类：基于扩散的修复方法、基于纹理的修复方法和基于稀疏表示的修复方法。基于扩散的修复方法常用于较小区域的修复，如基于偏微分方程的BSCB^[1](由Bertalmio、Sapiro与Casel Ballester提出)修复模型，该修复模型的基本思想是以单个像素为基本单位沿着等照度线按照由外向内的顺序依次进行迭代式修复，修复过程是通过求解一个三阶偏微分方程实现的，同时加入各向异性扩散方程，避免了等照度线的相交，使边缘平滑。基于纹理的修复方法常用于较大区域的修复，Criminisi等人^[2]提出了基于优先度顺序样本填充的修复算法，其核心思想是根据等照度线驱动采样过程，通过搜索最佳相似匹配块进行纹理复制，此方法通过平衡优先权函数中置信项和数据项来同时处理图像中的结构信息和纹理信息。但是该算法容易产生误差累积，同时对非线性结构的信息扩散和深度信息的处理不够理想，且图像中如果没有足够合理的采样块，修复结果通常不佳。基于稀疏表示的修复方法，Zhang^[3]采用解线性规划问题来处理图像修复，分析了稀疏分解在修复方面的一些应用，并给出了详细的理论依据。但是在稀疏表示的字典构造、稀疏表示特性约束与其他统计特性约束的结合等方面还存在较多问题，如容易生成伪影等问题。

分数阶在最近几年引起了广泛的关注，越来越多的学者将分数阶应用于图像处理中。常见应用有图像去噪、图像增强、图像配准、图像边缘检测以及图像修复。如文献[4]将分数阶引入到图像去噪中，在保持图像的弱边缘和纹理细节方面有较好的效果。文献[5]将分数阶应用于图像增强中，从而使图像边缘更加突出、纹理更加清晰且图像平滑区域信息得到较好保留。文献[6]将分数阶应用于图像配准中，由于分数阶微分比整数阶微分具有更好的细节描述能力，该方法能够有效地提高图像配准的精度，能用于包含较多弱纹理和弱边缘信息的医学图像配准。文献[7]将分数阶引入到图像的边缘检测中，将G-L分数阶梯度代替Canny中传统的梯度算子，不但可以增强图像的细节信息，而且可以增强灰度均匀和弱纹理区域的梯度信息，从而提高了边缘检测的精度和稳定性。文献[8]将分数阶引入图像修复中，提出了分数阶与TV模型相结合的方法，该方法对图像细节信息的修复具有较好的效果。

TV模型由于能去除图像中的噪声，近年来学者们将TV模型应用到图像修复中。如Chan等人^[9]提出了将TV模型应用到图像修复中，该方法是基于偏微分方程的图像修复算法，但是其收敛速度较慢，且容易产生过度平滑。针对基于偏微分方程修复算法运行速度较慢的问题，Telea^[10]提出了一种新的基于水平集的快速图像修复方法，该算法解决了修复算法运行速度较慢的问题，但是基于水平集的方法易造成区域模糊并且对边缘信息保持不佳。针对该问题，李开宇等人^[11]提出一种改进方案，在权函数的设计中，引入连续强度来保持边缘信息，并利用等照度线的方向计算出两个像素点的位置关系。当修复单个点时，引入置信因子来加权插值点。该算法在保证运行效率的同时提高了修复精度，但是该方法并没有解决纹理过度平滑问题。文献[8]提出了分数阶与TV模型相结合的修复方法，相对原始的TV模型有了较大提升，能够非线性地保留图像平滑区域的纹理信息，有效地解决纹理过度平滑问题，对图像细节信息具有较好的修复效果。但是该方法对具有弱导数性质的纹理细节等信息的保持仍不够理想，并且在模型最小化过程中导致了计算上的困难，其主要原因是正则项和数据项在零点处不可微。

针对上述问题，本文在文献[8]的基础上提出了一种结合纹理信息与分数阶TV模型的图像修复方法。由于文献[8]对具有弱导数性质的纹理细节等信息的保持仍然不够理想，并且没有充分利用已知区域的纹理先验信息，从而导致修复精度不够理想，本文在求解TV模型时，在梯度计算过程中增加了一个极小值，克服了正则项和数据项在零点处不可微的问题，增加模型的稳定性，从而可以有效地保持纹理区域的弱导数特性；此外在修复过程中充分利用已知区域的纹理先验信息，根据先验信息确定待修复区域的纹理方向，从而更好地保持图像中的纹理细节和弱边缘信息。

目前，分数阶微分主要有3种定义，分别为Grünwald-Letnikov(G-L)、Riemann-Liouville(R-L)和Caputo^[12]，其中G-L定义从整数阶微分的差分定义出发推导，而R-L和Caputo定义则根据整数阶积分的柯西公式计算推导。Caputo定义适用于分数阶微分方程初边值问题的分析，适用于工程等领域；其中G-L和R-L在数值运算时可以转化为卷积积分形式，因此适用于信号处理等领域，在计算中，G-L定义比R-L定义更准确，因此本文选取G-L定义的分数阶微积分。

G-L是比较经典的分数阶定义，它的实质是由整数阶推导出来的，根据整数阶微分的定义，在区间$ t\in [a, b](a<b, a, b\in \bf{R}) $内存在函数$ f\left( t \right) $连续可微，可得该连续函数的1阶微分定义为

$ {f^\prime }(t) = \mathop {\lim }\limits_{h \to 0} \frac{{f(t + h) - f(h)}}{h} $

(1)

式中，$ h(h\in \bf{R}) $表示步长，根据上面的1阶微分定义，可得出连续函数$ f\left( t \right) $的2阶微分，定义为

$ \begin{array}{*{20}{l}} {{f^{\prime \prime }}(t) = \mathop {\lim }\limits_{h \to 0} \frac{{{f^\prime }(t + h) - {f^\prime }(h)}}{h} = }\\ {\mathop {\lim }\limits_{h \to 0} \frac{{f(t + 2h) - 2f(t + h) + f(h)}}{{{h^2}}}} \end{array} $

(2)

以此类推，可将整数阶阶次从一阶推到$ n $阶，从而可得到连续函数$ f\left( t \right) $的$ n $阶微分，定义为

$ {f^n}(t) = \mathop {\lim }\limits_{h \to 0} \frac{1}{{{h^n}}}\sum\limits_{m = 0}^n {{{( - 1)}^m}} \left( {\begin{array}{*{20}{l}} n\\ m \end{array}} \right)f(t - mh) $

(3)

将整数阶$ n $扩展到整个实数范围，引入伽马函数，可以得出任意阶次$ \alpha $的G-L分数阶微分

$ \begin{array}{l} {}_a^G{\rm{D}}_t^\alpha f(t) = \mathop {\lim }\limits_{h \to 0} \frac{1}{{{h^\alpha }}}\sum\limits_{m = 0}^{\left[ {\frac{{t - a}}{h}} \right]} {{{( - 1)}^m}} \times \\ \frac{{\mathit{\Gamma} (\alpha + 1)}}{{m!\mathit{\Gamma} (\alpha - m + 1)}}f(t - mh) \end{array} $

(4)

式中，$ G $表示G-L分数阶定义，$ \alpha $代表任意实数阶次，$ h $为步长，[]表示取整符号，$ \mathit{\Gamma} (n) $为伽马函数，定义为

$ \mathit{\Gamma }(n) = \int_0^\infty {{{\rm{e}}^{ - t}}} {t^{n - 1}}{\rm{d}}t = (n - 1)! $

(5)

对于图像的中低频区域^[13]，当阶次选择在0~1之间时，用分数微分算子处理后得到的值大于整数阶(包括1阶和2阶)微分算子的值。然而在图像高频区域中，分数阶微分算子提升幅度小于1阶和2阶微分算子，阶次越小，提升的幅度越小。因此，可以得出如下结论：分数阶微分具有“弱导数”的性质，它不仅能增强信号的高频部分，而且能非线性地保留中低频。在数字图像中，图像的纹理细节和弱边缘对应的是中低频信息，图像的噪声对应着高频信息。用整数阶微分处理图像，其弱边缘和纹理信息将较大地衰减，噪声信号将被较大地增强。但是，分数阶微分可以用来克服这一缺陷，噪声信号不会被极大地增强，而弱边缘以及纹理信息同时会被非线性地保留。

图像修复方法一般都是由已知区域的边缘区域$\boldsymbol{E}$(邻域)不断地向待修复区域$\boldsymbol{D}$扩散的过程。TV模型图像修复^[9]，即全变分模型^[14-15]。TV模型图像修复方法是一种基于偏微分方程的图像修复方法，其修复模型为

$ {J_\lambda }(\mathit{\boldsymbol{u}}) = \int\limits_{\mathit{\boldsymbol{E}} \cup \mathit{\boldsymbol{D}}} {\left| {\nabla \mathit{\boldsymbol{u}}} \right|{\rm{d}}x{\rm{d}}y} + \frac{\lambda }{2}\int_\mathit{\boldsymbol{E}} {{{\left| {\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^0}} \right|}^2}} {\rm{d}}x{\rm{d}}y $

(6)

式中，第1项为正则项，它的作用是在$ \boldsymbol{D}\cup \boldsymbol{E} $区域内的图像总变分中求取最优。第2项为数据保真项，它的作用是平滑待修复区域。$\boldsymbol{u} $为参考图像，$ {{\boldsymbol{u}}^{0}} $为待修复图像，$ \lambda $为正则化参数，用于平衡正则项和数据项的权重。

由式(6)可知，能量泛函的梯度下降方程为

$ -\operatorname{div}\left(\frac{\nabla \boldsymbol{u}}{|\nabla \boldsymbol{u}|}\right)+\lambda_{e}\left(\boldsymbol{u}-\boldsymbol{u}^{0}\right)=0 $

(7)

式中，${{\lambda }_{e}}=\left\{ \begin{array}{*{35}{l}} \lambda , x\in \boldsymbol{E} \\ 0, x\in \boldsymbol{D} \\ \end{array} \right. $，$ \operatorname{div}\left( \frac{\nabla \boldsymbol{u}}{|\nabla \boldsymbol{u}|} \right) $为扩散项，$ \frac{1}{|\nabla \boldsymbol{u}|} $为扩散系数，其作用是平衡图像在不同梯度下的扩散快慢，以此来避免图像修复后破损区域与其邻域之间平滑过度。由于TV模型中的扩散是各向异性的，所以TV修复模型不仅有较好的修复效果，而且还有优良的保边性能。但是TV修复模型在修复图像时只能有效地逼近分片常数，因此，在修复结果中常常出现“阶梯效应”。

由于图像是离散的，图像中点的邻域如图 1所示。其中$ O $为污染点，邻域$ \boldsymbol{B}\in \{N, S, W, E\} $，半像素邻域$ \boldsymbol{{B}'}\in \left\{ {N}', {S}', {W}', {E}' \right\} $。

因此，将式(6)离散化为

$ \sum\limits_{\begin{array}{*{20}{c}} {X \in \boldsymbol{B}}\\ {x \in \boldsymbol{B'}} \end{array}} {\frac{1}{{\left| {\nabla {\boldsymbol{u}_x}} \right|}}} \left( {{\boldsymbol{u}_o} - {\boldsymbol{u}_X}} \right) + {\lambda _e}(O)\left( {{\boldsymbol{u}_O} - \boldsymbol{u}_O^0} \right) = 0 $

(8)

式中，$ {{\lambda }_{e}}=\left\{ \begin{array}{*{35}{l}} \lambda , x\in \boldsymbol{E} \\ 0, x\in \boldsymbol{D} \\ \end{array} \right. $, $ {{\boldsymbol{u}}_{O}} $为$ O $点修复后的像素，$ \boldsymbol{u}_{O}^{0} $为$ O $点修复前的原始像素，$ {{\boldsymbol{u}}_{x}} $表示点$ O $半邻域的像素，$ {{\boldsymbol{u}}_{X}} $表示点$ O $邻域的像素。

化简得到最终的表达式为

$ {\mathit{\boldsymbol{u}}_O} = \frac{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{{{\mathit{\boldsymbol{u}}_X}}}{{\sqrt {\left| {\nabla {\mathit{\boldsymbol{u}}_x}} \right|} }} + {\lambda _e}(O)\mathit{\boldsymbol{u}}_O^0} }}{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{1}{{\sqrt {\left| {\nabla {\mathit{\boldsymbol{u}}_x}} \right|} }} + {\lambda _e}(O)} }} $

(9)

对上式进行迭代，直到获得满意的修复效果。

分数阶微分的TV模型为如下形式^[8]

$ {J_\lambda }(\mathit{\boldsymbol{u}}) = \int\limits_{\mathit{\boldsymbol{E}} \cup \mathit{\boldsymbol{D}}} {\left| {{\nabla ^\alpha }\mathit{\boldsymbol{u}}} \right|{\rm{d}}x{\rm{d}}y} + \frac{\lambda }{2}\int_\mathit{\boldsymbol{E}} {{{\left| {\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^0}} \right|}^2}} {\rm{d}}x{\rm{d}}y $

(10)

式中，$\boldsymbol{u} $代表参考图像，$ {{\boldsymbol{u}}^{0}} $为待修复图像，$ \lambda $为正则化参数。使用分数阶梯度代替整数阶梯度，并将其离散化，得到

$ -\operatorname{div}\left(\frac{\nabla^{\alpha} \boldsymbol{u}}{\left|\nabla^{\alpha} \boldsymbol{u}\right|}\right)+\lambda_{e}\left(\boldsymbol{u}-\boldsymbol{u}^{0}\right)=0 $

(11)

同式(7)，化简得到最终的表达式为

$ {\mathit{\boldsymbol{u}}_O} = \frac{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{{{\mathit{\boldsymbol{u}}_X}}}{{\sqrt {\left| {{\nabla ^\alpha }{\mathit{\boldsymbol{u}}_x}} \right|} }} + {\lambda _e}(O)\mathit{\boldsymbol{u}}_O^0} }}{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{1}{{\sqrt {\left| {{\nabla ^a}{\mathit{\boldsymbol{u}}_x}} \right|} }} + {\lambda _e}(O)} }} $

(12)

基于TV模型的图像修复会导致纹理过度平滑，文献[8]将G-L分数阶微分引入TV修复模型中，该方法较好地解决了过度平滑的问题。但是该方法对具有弱导数性质的纹理细节等信息的保持仍不够理想，并且在模型最小化过程中导致了计算困难，其主要原因是正则项和数据项在零点处不可微。针对该问题，本文在模型的求解过程中，引入一个极小值，即在式(11)的正则项的分母中增加一个极小值$ p $

$ -\operatorname{div}\left(\frac{\nabla^{\alpha} \boldsymbol{u}}{\left|\nabla^{\alpha} \boldsymbol{u}\right|+p}\right)+\lambda_{e}\left(\boldsymbol{u}-\boldsymbol{u}^{0}\right)=0 $

(13)

同式(7)，进行化简得到如下表达式

$ {\mathit{\boldsymbol{u}}_O} = \frac{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{{{\mathit{\boldsymbol{u}}_X}}}{{\sqrt {\left| {{\nabla ^\alpha }{\mathit{\boldsymbol{u}}_x}} \right| + {p^2}} }} + {\lambda _e}(O)\mathit{\boldsymbol{u}}_O^0} }}{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{1}{{\sqrt {\left| {{\nabla ^a}{\mathit{\boldsymbol{u}}_x}} \right| + {p^2}} }} + {\lambda _e}(O)} }} $

(14)

式中，参数$ \alpha $为分数阶的阶次，当$ \alpha =1 $时，该分数阶TV模型即文献[9]的1阶TV模型，所以分数阶TV模型实际上是整数阶TV模型的一个扩展和延伸。式中的参数$ p $为引入的极小值，其选择原则如下：当$ p $比较小时可保持锐利的边缘；当$ p $较大时可使得扩散程度较高。因此，当待修复区域的面积较小时，选取较小的$ p $值可获得较好的修复效果；但当待修复区域的面积较大时，应选取较大的$ p $值，这是因为$ p $值较小时，扩散度较低，会导致待修复区域中心的初始噪声在边界信息扩散前就被平滑形成锐利边界，从而得到错误的修复结果，而选取较大的$ p $值则可避免该问题。

文献[8]在修复过程中没有充分利用已知区域的纹理先验信息，本文根据图像已知区域的先验来确定待修复区域的纹理方向，从而在图像修复过程中更好地保持了图像中的纹理细节和弱边缘信息。

待修复区域内部可以通过待修复区域的外围区域信息来获取尽可能近似的局部纹理方向。$\boldsymbol{E}$为已知区域，$\boldsymbol{D}$为待修复区域，记$ I\left( i, j \right) $为$ \left( i, j \right)$处的像素值，$\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$为整幅图像区域。若$ \mathit{\boldsymbol{\rho }} = ({\rho _x}, {\rho _y}) \in {{\bf{Z}}^2} $满足：对任意的$ \left( {i, j} \right) \in \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $，$ I\left( {i, j} \right) \approx I(i + {\rho _x}, j + {\rho _y}) $，则称$ \mathit{\boldsymbol{\rho }} $为重复向量。一般情况下图像的重复向量只有一个，寻找纹理的方向就是确定$ \mathit{\boldsymbol{\rho }} $的过程。对区域$\boldsymbol{E}$中任意两点$ ({i_1}, {j_1}) $, $ ({i_2}, {j_2}) $，如果$ ({i_1}, {j_1}) $和$ ({i_2}, {j_2}) $足够接近，则分别以这两点为中心，固定一个3×3像素的匹配窗，记集合

$ \mathit{\boldsymbol{V}} = \left\{ {\begin{array}{*{20}{l}} {\left( {{i_1} + k,{j_1} + h} \right),\left( {{i_2} + k,{j_2} + h} \right) \in \mathit{\boldsymbol{S}}}\\ {(k,h) \in [ - 1,1] \times [ - 1,1]} \end{array}} \right\} $

(15)

若平均差

$ M_{\text{mad}}=\frac{1}{|\boldsymbol{V}|}\left\{\sum\limits_{(k, h) \in \boldsymbol{V}}\left|I\left(i_{1}+k, j_{1}+h\right)-I\left(i_{2}+k, j_{2}+h\right)\right|\right\} $

小于给定的阈值，则把$ ({i_2} - {i_1}, {j_2} - {j_1}) $作为候选方向之一，放入候选集合$ \mathit{\boldsymbol{C}} $。确定$ \mathit{\boldsymbol{C}} $后，采用平均绝对差$ {A_{{\rm{mad}}}} $的方法，从$ \mathit{\boldsymbol{C}} $中选出最优方向。对于某一候选方向$ \mathit{\boldsymbol{c}} = ({c_x}, {c_y}) \notin \mathit{\boldsymbol{C}} $，若$ \mathit{\boldsymbol{\rho }} $满足

$ \mathit{\boldsymbol{\rho }} = \mathop {\arg \min }\limits_{\left( {{c_x},{c_y}} \right) \in \mathit{\boldsymbol{C}}} \frac{1}{{|\mathit{\boldsymbol{R}}|}}\sum\limits_{(u,v) \in \mathit{\boldsymbol{R}},\left( {{u^\prime },{v^\prime }} \right) \in \mathit{\boldsymbol{S}}} {\left| {I(u,v) - I\left( {{u^\prime },{v^\prime }} \right)} \right|} $

(16)

则将$ \mathit{\boldsymbol{\rho }} $作为最佳纹理方向。其中$ \mathit{\boldsymbol{S}} $为已知像素点集合，集合$ \mathit{\boldsymbol{R}} $是区域$\boldsymbol{D}$的外环子集，$ \left( {u, v} \right) \in \mathit{\boldsymbol{R}} $满足：

1) 若关于方向$\mathit{\boldsymbol{c}} = ({c_x}, {c_y}) $，$\left( {u, v} \right) $的两个对应点$ (u + {c_x}, v + {c_y}) $和$ (u - {c_x}, v - {c_y})$至少满足$ \left( {u', v'} \right) = (u + {c_x}, v + {c_y}) \in \mathit{\boldsymbol{S}} $或者$ \left( {u', v'} \right) = (u - {c_x}, v - {c_y}) \in \mathit{\boldsymbol{S}} $；

2) 若$ (u + {c_x}, v + {c_y}) \in \mathit{\boldsymbol{S}} $和$ (u - {c_x}, v - {c_y}) \in \mathit{\boldsymbol{S}} $同时成立，由于在纹理上，$ (u + {c_x}, v + {c_y}) $和$ (u - {c_x}, v - {c_y}) $差别不大，选$ \left( {u', v'} \right) = (u + {c_x}, v + {c_y}) $或$(u - {c_x}, v - {c_y}) $均可。

本文实验中有3个参数$ \lambda $、$ p $和$ \alpha $。正则化参数$ \lambda $值影响着正则项和数据项的权重，根据经验，本文取$ \lambda =0.005 $，参数$ p $初始值取0.05。

下面分析分数阶微分阶次$ \alpha $的选取，对于可积函数$f\left( x \right) $，其傅里叶变换为$ F(u) = \int_{\bf{R}} f (u){{\rm{e}}^{ - vux}}{\rm{d}}x $，根据傅里叶变换的性质，其$ n(n \in {{\rm{Z}}^ + }) $阶导数为$ {D^n}f(x)\mathop \Leftrightarrow \limits^{{\rm{FT}}} {(DF)^n}(u) = {(vu)^n}F(u) = {\hat h^n}(u)F(u) $，把阶次从整数阶扩展到分数阶可得

$ \begin{array}{*{20}{c}} {{D^\alpha }f(x)\mathop \Leftrightarrow \limits^{{\rm{FT}}} {{(DF)}^\alpha }(u) = {{(vu)}^\alpha }F(u) = {\hat h ^\alpha }(u)F(u)}\\ {{{\hat h}^\alpha }(u) = {{\hat b}^\alpha }(u){{\rm{e}}^{v{\theta ^\alpha }(u)}}} \end{array} $

(17)

式中，$ {\hat b^\alpha }(u) = |u{|^\alpha } $，$ {\hat \theta ^\alpha }(u) = \frac{{\alpha {\rm{ \mathsf{ π} }}}}{2}{\mathop{\rm sgn}} (u) $，根据式(17)绘制如图 2所示的分数阶幅频特性曲线图。

从图 2可以看出，频率在0~1之间时，阶次越小的分数阶微分算子对幅值提升越快。当频率大于1时，阶次较小的分数阶微分算子对幅值提升变得缓慢，随着阶次变大，对中低频部分的幅值上升比较快。当阶次小于1时，可作为中低频增强释放和高频压缩的滤波器。而当阶次大于1时，可看做一个高通滤波器。在图像处理中，图像低频部分对应着图像的纹理信息以及图像的弱边缘，图像的轮廓和噪声对应着高频部分。

总之, 根据不同目标选取不同的阶次。分数阶微分的阶次越大，中频和高频的信号越能够被有效地增强。而在图像修复中，中低频对应图像破损或缺失的平滑部分，高频则对应着图像的边界部分。本文目标是解决边界没有得到充分扩散以及过度平滑的修复问题，因此需要保留边缘信息及纹理信息。从幅频特性曲线可以看出，当阶次在0~1之间时，信号可以得到有效增强，故本文的阶次建议取值0~1之间，经过后期的实验也证明了这一点。

1) 根据$ \mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^0} $计算差值图像，确定图像的待修复区域；

2) 根据

$ \mathit{\boldsymbol{\rho }} = \mathop {\arg \min }\limits_{\left( {{c_x},{c_y}} \right) \in \mathit{\boldsymbol{C}}} \frac{1}{{|\mathit{\boldsymbol{R}}|}}\sum\limits_{\begin{array}{*{20}{c}} {(u,v) \in \mathit{\boldsymbol{R}}}\\ {\left( {{u^\prime },{v^\prime }} \right) \in \mathit{\boldsymbol{S}}} \end{array}} {\left| {I(u,v) - I\left( {{u^\prime },{v^\prime }} \right)} \right|} $

确定图像的纹理方向；

3) 为了增加模型的稳定性，根据式(7)，在能量泛函的梯度下降方程引入极小值$ p $，如

$ - {\mathop{\rm div}\nolimits} \left( {\frac{{{\nabla ^\alpha }\mathit{\boldsymbol{u}}}}{{\sqrt {{{\left| {{\nabla ^\alpha }\mathit{\boldsymbol{u}}} \right|}^2} + {p^2}} }}} \right) $

4) 根据式(13)最终化简得到

$ {\mathit{\boldsymbol{u}}_O} = \frac{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{{{\mathit{\boldsymbol{u}}_X}}}{{\sqrt {\left| {{\nabla ^\alpha }{\mathit{\boldsymbol{u}}_x}} \right| + {p^2}} }} + {\lambda _e}(O)\mathit{\boldsymbol{u}}_O^0} }}{{\sum\limits_{\begin{array}{*{20}{c}} {X \in \mathit{\boldsymbol{B}}}\\ {x \in \mathit{\boldsymbol{B'}}} \end{array}} {\frac{1}{{\sqrt {\left| {{\nabla ^a}{\mathit{\boldsymbol{u}}_x}} \right| + {p^2}} }} + {\lambda _e}(O)} }} $

5) 根据经验取$ p = \frac{p}{5} $进行迭代，经过多次循环迭代得到最佳的$ p $值，输出$ {{\boldsymbol{u}}_{O}} $。

实验所使用的计算机环境为：Intel Core i3-2130 CPU，3.40 GHz，内存4 GB，操作系统为64位Windows 7.0，程序使用R2014a版MATLAB实现。为验证本文算法对具有弱纹理、弱边缘特征的图像具有更好的修复性能，实验图像选择标准库中不同程度纹理信息的Barbara和Lena图像(如图 3所示)。为了提高效率，在原图像上的纹理部分截取部分(如图 3中的红色方框所示，截取后的图像如图 4 (a)和4 (d)所示)，待修复图像为在截取的图像上进行人为破坏，如图 4 (b)和4 (e)所示。另外为了显示本文算法更具说服力，本文再选取具有代表性实验图像中的岩石图像(如图 3(c)所示)。为了提高效率，在原图像上的纹理部分截取一小部分(如图 3(c)中的红色方框所示，截取后的图像如图 4 (g)所示)，待修复图像为在截取的图像上进行水印破坏，如图 4 (h)所示。本文从不同参数(阶次及模板大小)以及修复性能验证两个方面进行实验，并对实验结果进行分析和讨论。

为了证明本文方法的有效性，从定性和定量两个方面来评价修复效果。定性的评价即视觉效果，比较修复后的图像和原始图像的差异，通过二者差值图像黑色区域所占大小来评价，定量评价则通过修复后图像和原始图像的灰度均方误差MSE和峰值信噪比PSNR来评价。

灰度均方误差定义为

$ f_{\text{MSE}}=\frac{1}{m \times n} \sum\limits_{i, j}\left(I_{1}(i, j)-I(i, j)\right)^{2} $

(18)

峰值信噪比定义为

$ {f_{{\rm{PSWR}}}} = 10\lg \frac{{\mathop {\max }\limits_{1 \le i \le n,1 \le n} {{\left| {{I_1}(i,j)} \right|}^2}}}{{\frac{1}{{m \times n}}\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{\left[ {{I_1}(i,j) - I(i,j)} \right]}^2}} } }} $

(19)

式中，$ {I_1} $是参考图像的像素点，$ I $是修复后图像的像素点，$ m \times n $为像素个数。理想情况下MSE值应为0，表示两幅图像同一位置上点的灰度值相同。PSNR值越高越好，越高表示修复后图像越接近参考图像。

本实验选用图 4(a)所示的Barbara作为参考图像进行实验，为了证明本文方法的有效性，对Barbara图像进行不同模板和不同阶次的取值，图 5是本文算法取不同阶次和3个不同模板宽度$ k $所对应的修复结果，实验结果如图 5和表 1、表 2所示。

表 1为不同模板和不同阶次之间MSE的实验结果，表 2为不同模板和不同阶次之间峰值信噪比的实验数据。从表中可以看出，MSE最小值不在阶次$ \alpha $等于1时，当阶次从0.1开始逐渐增大，MSE缓慢增加，当$ \alpha $值在0.1处时，MSE值最小、PSNR值最大(如表 1、表 2中加粗字体所示)。为了更直观地表示修复效果，本文绘出均方误差MSE和阶次$ \alpha $的关系图(如图 6)，峰值信噪比PSNR和阶次$ \alpha $、模板$ k $的关系图(如图 7)。图 6展示了均方误差(MSE)和阶次$ \alpha $的关系，MSE值越小表示修复效果越好，图 7展示了峰值信噪比(PSNR)和阶次$ \alpha $、模板宽度$ k $的关系，PSNR值越大表示修复效果越好。由图 6、图 7可以看出，均方误差图中红色线条最低，峰值信噪比图中红色线条最高，并且当阶次$ \alpha $等于0.1时，均方误差最小，同时峰值信噪比最大。即当模板为3×3、阶次$ \alpha $等于0.1时，修复效果最佳。这是因为当模板尺寸过大时，会对图像造成过度平滑，对于纹理比较多的图像，应该选择较小的模板。本文图片纹理信息较多，综合考虑选取3×3大小的模板。对于弱纹理的图像，当阶次$ \alpha $比较大时，会导致高频信号如边缘信息模糊，从而降低修复效果，所以阶次$ \alpha $等于0.1时效果较好。在弱纹理图片中，需要增强中频和低频部分区域的梯度信息。经过多次实验，Barbara图片纹理方向为(-2, 3)时修复效果最好。所以当模板为3×3、阶次$ \alpha $等于0.1时，均方误差最小、峰值信噪比最大，即修复效果最好。

本实验选用图 4(d)所示的Lena图像重复进行上述实验，实验结果如图 8、表 3和表 4所示。

表 3为不同模板和不同阶次之间均方差的实验数据，表 4为不同模板和不同阶次之间峰值信噪比的实验数据。从表中可以看出，MSE最小值不在阶次$ \alpha $等于1时，当阶次从0.1逐渐增大，MSE值急剧下降，$ \alpha $值在0.9处时，MSE值最小、PSNR值最大(如表 3和表 4中加粗字体所示)。为了更直观地看出阶次为0.9时，修复效果最好，绘出均方误差MSE和阶次$ \alpha $的关系图(如图 9)，峰值信噪比PSNR和阶次$ \alpha $、模板$ k $的关系图(如图 10)。

经过多次实验，Lena图像的纹理方向为(-1, 1)时修复效果最好。为了得到最佳的阶次，根据灰度均方误差和峰值信噪比来分析参数$ \alpha $对修复的影响进一步实验。图 9展示了均方误差和阶次$ \alpha $的关系，MSE值越小表示修复效果越好，图 10展示了峰值信噪比和阶次$ \alpha $的关系，PSNR值越大表示修复效果越好。可以看出，均方误差图中蓝色线条最低，峰值信噪比图中蓝色线条最高，并且当阶次$ \alpha $等于0.9时，均方误差最小，同时峰值信噪比最大。即当模板为3×3、阶次$ \alpha $等于0.9时，修复效果最佳。若模板尺寸过小，则会增加算法的计算量，对于纹理不明显的图像，应该选择较大的模板。本图片纹理信息较弱，综合考虑本实验选取5×5大小的模板，经实验验证，当模板为5×5、阶次$ \alpha $等于0.9时，均方误差最小、峰值信噪比最大，即修复效果最好。

本文第3次实验选用图 4(g)所示的岩石图像重复上述实验，原理同上。经过多次实验表明：当$ \alpha $值取0.1，且$ k $等于1时，MSE值最小、PSNR值最大。即当模板为3×3、阶次$ \alpha $等于0.1时，修复效果最佳。这是因为当模板尺寸过大时，会对图像造成过度平滑，对于纹理比较多的图像，应该选择较小的模板。本图片纹理信息较多，综合考虑本实验选取3×3大小的模板。对于弱纹理的图像，当阶次$ \alpha $比较大时，会导致高频信号如边缘信息模糊，从而降低修复的效果，所以阶次$ \alpha $等于0.1时效果较好。在弱纹理图片中，需要增强中频和低频部分区域的梯度信息。经过多次实验，岩石图像纹理方向为(1, 1)时修复效果最好。所以当模板为3×3、阶次$ \alpha $等于0.1时，均方误差最小、峰值信噪比最大，即修复效果最好。

本文实验选取了不同程度的3种不同纹理图片，实验结果显示不同纹理最佳模板和阶次不一样。Barbara图像的纹理比较丰富，Lena图像的纹理较弱，岩石图像的纹理也比较丰富。

通过实验可知，Barbara图像在分数阶阶次$ \alpha $等于0.1，模板宽度$ k $等于1(即为3×3)时，修复效果达到最佳。Lena图像在阶次为0.9时(模板为5×5)修复效果最佳。岩石图像在阶次为0.1时(模板为3×3)修复效果最佳。然而，对于不同的图像，实验参数仍然需要人工选定进行多次实验。同时也可以得出结论：1阶微分不能有效地描绘具有弱导数性质的纹理细节等信息。近年来的研究也表明，分数阶微分在描述纹理信息能力上优于整数阶微分。因此，在TV算法中使用分数微分代替整数阶微分，可以获得更好的修复效果。

为了验证本文算法修复性能，选取了3幅存在不同程度的弱边缘和弱纹理的图像进行实验。将本文算法分别与TV模型、文献[2]算法、文献[8]算法进行对比，并分别进行定性和定量分析。

本次实验仍选用Barbara图像、Lena图像和岩石图像进行实验，如图 11、图 12和图 13所示。由上节实验可知，对于Barbara图像，模板大小为3×3，阶次为0.1时，修复效果最佳；对于Lena图像，模板大小为5×5，阶次为0.9时，修复效果最佳；对于岩石图像，模板大小为3×3，阶次为0.1时，修复效果最佳，本节实验即采用该参数值。图 4(a)-(c)、图 4(d)-(f)以及图 4(g)-(i)分别是参考图像、待修复图像以及掩码图像，图 11(a)-(d)分别是用TV模型、文献[2]算法、文献[8]算法和本文算法的修复结果图像以及差值图像。通过比较可以发现，本文算法的差值图像白色区域最少，这表明本文算法相比其他算法具有更好的修复效果。为了方便观察，TV算法、文献[2]算法、文献[8]算法和本文算法下修复前和修复后的MSE值和PSNR值如表 5和表 6所示。

从表 5和表 6中可以看出，本文算法对3幅破损图像均能得到较好的修复效果。原有的TV模型在图像修复时会导致纹理过度平滑，文献[8]算法加入了分数阶可以很好地保持图像的纹理信息，但并没有充分利用纹理的边缘信息。本文在求解梯度时引入极小化参数来增加模型的稳定性，并且根据已知修复区域的先验来确定纹理的修复方向，可以充分利用纹理的边缘信息。实验结果表明，本文修复效果比其他比较方法更好。此外本文选择的Barbara图像和岩石图像纹理明显比Lena图像的纹理要丰富，但是从实验结果可以看出，对于Lena图像修复的效果比Barbara和岩石修复的效果要好得多。故可以得出如下结论：对于纹理越丰富的弱图像，本文修复效果越佳；对于弱纹理的图像，本文修复效果更佳。

将G-L分数阶微分引入到TV模型中，一方面在TV模型求解梯度时引入极小化参数，解决了正则项和数据项在零点处不可微的问题，增加了模型的稳定性，从而可以有效地保持纹理区域的弱导数特性；另一方面在修复过程中根据图像已知区域的先验来确定待修复区域的纹理方向，充分利用了图像中的纹理细节和弱边缘信息，从而提高了修复的精度。此外，通过实验给出了修复效果和阶次$ \alpha $、模板宽度$ k $的关系，从而为最佳模板参数的选取提供了依据，尽管不同类型的图像其最佳参数是不同的, 但是其最佳修复阶次一般在0~1之间，这是因为图像的平滑部分对应着信号的低频部分，图像的纹理细节对应着信号的中频部分，而TV算法对于弱纹理区域修复效果不理想，增强该区域的梯度信息即需提高低频和中频部分，所以选用阶次0~1之间修复效果更佳。理论分析和实验结果均表明，本文算法能提升弱纹理和弱边缘图像的修复精度，是TV模型的一个重要延伸和扩展。

但是不同图像修复的最佳阶次需要不断测试，因而是比较耗时和费力的。今后可以研究自适应G-L分数阶TV的图像修复算法。此外本文的分数微分掩模是2维的，要完成3维图像的修复, 还需将其扩展到3维空间。