数字图像处理技术的快速发展在给人们带来生活便捷的同时, 也存在着图像本身的真实性和完整性等问题, 目前篡改后的数字图像出现在医学研究、新闻报道或法庭证据等重要场合, 这会对社会造成极大的负面影响, 因此, 对于检测数字图像是否真实得到了人们的密切关注. 当前图像篡改检测算法主要分为两大类: 一是主动检测算法, 另一种是被动检测算法^[1]. 被动检测算法在没有任何先验信息的情况下, 利用图像的本身特性对原图像进行真实性检测, 此类检测算法已经得到了广泛应用^[2].

复制粘贴篡改是生活中常见的图像篡改方式. 复制粘贴篡改图像的检测方法可以分为两类: 基于特征点的检测算法和基于图像块的检测算法, 基于图像块的算法在准确定位到篡改区域上优于基于特征点的算法. Ryu SJ等^[3]提出了一种基于Zernike矩定位复制图像区域的算法, 设计了一种局部敏感散列的新型块匹配, 并通过检查矩的相位来减少误匹配, 但对于较大尺度变换的检测效果不是很好. Dixit R等^[4]利用傅里叶-梅林变换和对数极坐标映射以及使用K均值聚类的基于颜色的分割技术, 具有较好的平移和旋转不变性, 但算法的实时性差. 谷宗运等^[5]提出了一种基于Tchebichef矩的篡改图像检测算法, 比较提取的DWT和Tchebichef矩相邻特征向量的相似性, 实现篡改区域的定位, 但对于多区域篡改区域定位效果较差.

针对上述算法的不足, 本文提出一种低频快速切比雪夫矩算法的篡改图像检测算法, 首先采用非抽样小波变换对图像二维分解, 对低频部分进行重叠分块, 再提取图像块的快速切比雪夫矩特征, 采用PatchMatch算法对特征向量进行匹配, 最后用稠密线性拟合算法消除误匹配以及形态学操作定位篡改区域.

1 构造特征向量 1.1 Tchebichef矩

假设图像 $f\left( {x,y} \right)$ 上一点(p, q), 位于大小为N×N图像块(p, p+N–1)×(q, q+N–1), 其n+m阶Tchebichef矩为^[6]:

$\left\{\begin{aligned} & T_{n,m}^{p,q} = \sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {f\left( {p + x,q + y} \right)} } {t_n}\left( x \right){t_m}\left( y \right) \\ & n,m = 0,1,\cdots,N - 1 \end{aligned} \right.$

(1)

其中, ${t_n}\left( x \right)$ 和 ${t_m}\left( y \right)$ 为n阶和m阶正交切比雪夫多项式, 定义如下:

$ \begin{aligned} {{t}_{n}}\left( x \right)=& \frac{\left( 1-N \right)!}{\sqrt{{}^{\left( N+{\rm n} \right)!}\!\!\diagup\!\!{}_{\left( 2n+1 \right)\left( N-n-1 \right)!}\;}} \\ & \sum\limits_{l=0}^{n}{\frac{{{(-n)}_{l}}{{\left( -x \right)}_{l}}{{\left( 1+n \right)}_{l}}}{{{\left( l! \right)}^{2}}{{\left( 1-N \right)}_{l}}}},\;n,x=0,1,\cdots,N-1 \end{aligned} $

(2)

其中, ${\left( \cdot \right)_l}$ 表示阶乘幂. ${t_m}\left( y \right)$ 定义类似.

正交切比雪夫多项式在x处的 ${t_n}\left( x \right)$ 和在x+a处的 ${t_n}\left( {x + a} \right)$ 有如下关系^[7]:

${t_n}\left( {x + a} \right) = \sum\limits_{l = 0}^n {\sum\limits_{r = 0}^{n - l} {{g_r}\left( {n,l} \right)} } {\left( { - a} \right)_r}{t_l}\left( x \right)$

(3)

当a=1时, 有:

$t\left( {x + 1} \right) = {t_n}\left( x \right) - \sum\limits_{l = 0}^{n - 1} {{g_r}} \left( {n,l} \right){t_l}\left( x \right)$

(4)

其中, ${g_r}\left( {n,l} \right)$ 定义为:

$ {g_r}\left( {n,l} \right) = {\left( { - 1} \right)^n}\frac{{\left( {2l + 1} \right)\left( {N - 1 - l} \right)!}}{{r!\left( {N - 1 - n} \right)!}} \times \sqrt {\frac{{\rho \left( {l,N} \right)}}{{\rho \left( {n,n} \right)}}} \sum\limits_{s = l}^{n - r} {\frac{{{{\left( { - 1} \right)}^s}s!\left( {n + r + s} \right)!\left( {N - 1 - s - r} \right)!}}{{\left( {r + s} \right)!\left( {n - r - s} \right)!\left( {s + l + 1} \right)!\left( {s - l} \right)!\left( {N - 1 - s} \right)!}}} ,\;0 \le l \le n $

(5)

1.2 改进的低频快速Tchebichef矩

对于图像篡改, 如果非平移不变, 在复制和粘贴两个相同的区域会被破坏, 会出现漏检情况. 下采样使得离散小波变换(DWT)不具有平移不变性, 不仅会对 DWT系数产生巨大影响, 还会对篡改区域产生不同的特征向量. 另外, DWT的伪吉布斯现象使得检测边缘和纹理等信号的效果不够理想. 对于DWT存在的不足, 这里采用具有平移不变性的非抽样小波变换(UWT), 由于UWT不包括下采样和小波系数缩减, 因此称为非抽样的^[8].

对图像利用UWT沿行和列进行二维分解, 分别得到4个子带, 即水平高通子带LH, 垂直高通子带HL、对角高通子带HH以及低通子带LL, 每个子带的尺寸不发生变化^[9].

图像进行过UWT后, 由于图像的主要部分是低频部分, 所以提取图像的低频部分. 对提取的低频部分进行检测, 可以大大降低块的个数, 使得计算量仅为原来的1/4, 同时低频对噪声不敏感, 也可以加强特征的鲁棒性. 再对图像进行分块, 假设待测图像大小为M×N, 用a×a像素的滑动窗口每次移动一个像素点对图像进行扫描, 可以得到(M–a+1)×(N–a+1)个重叠块.

对于一个大小为N×N滑动块(p+1, p+N)×(q, q+N–1), 其n+m阶水平方向的Tchebichef矩为:

$\left\{\begin{aligned} & T_{n,m}^{p + 1,q} = \sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {f\left( {p + 1 + x,q + y} \right){t_n}(x){t_m}(y)} } \\ & n,m = 0,1,\cdots,N - 1 \end{aligned} \right.$

(6)

我们可以看出矩 $T_{n,m}^{p + 1,q}$ 计算量的阶数已经上升到(M, M), 其中 $0 \le n,m \le M,M < N$ . $TC_{\rm{m}}^{p,q}$ 和 $TC_{\rm{m}}^{p + N,q}$ 分别表示输出列 $f\left( {p,q + y} \right),0 \le {{y}} \le N - 1$ , m阶切比雪夫矩和输入列 $f\left( {p + N,q + y} \right),0 \le {\rm{y}} \le N - 1$ , m阶切比雪夫矩, 定义如下^[10]:

$TC_m^{p,q} = \sum\limits_{y = 0}^{N - 1} {f\left( {p,q + y} \right)} {t_m}\left( y \right)$

(7)

$TC_m^{p + N,q} = \sum\limits_{y = 0}^{N - 1} {f\left( {p + N,q + y} \right)} {t_m}\left( y \right)$

(8)

相同地, 对于下一个大小为N×N滑动块(p, p+N–1)×(q+1, q+N), 其n+m阶垂直方向的Tchebichef矩为:

$\left\{\begin{aligned} & T_{n,m}^{p,q + 1} = \sum\limits_{x = 0}^{N - 1} {\sum\limits_{y = 0}^{N - 1} {f\left( {p + x,q + 1 + y} \right){t_n}(x){t_m}(y)} } \\ & n,m = 0,1,\cdots,N - 1 \end{aligned} \right.$

(9)

$TR_{\rm{m}}^{p,q}$ 和 $TR_{\rm{m}}^{p + N,q}$ 分别表示输出行 $f\left( {p + x,q} \right),$ $0 \le {\rm{x}} \le $ $ N - 1 $ , m阶切比雪夫矩和输入行 $f\left( {p + x,q + N} \right),$ $0 \le {\rm{x}} \le $ $ N - 1 $ , m阶切比雪夫矩, 定义如下:

$TR_n^{p,q} = \sum\limits_{x = 0}^{N - 1} {f\left( {p + x,q} \right)} {t_n}\left( x \right)$

(10)

$TR_n^{p,q + N} = \sum\limits_{x = 0}^{N - 1} {f\left( {p + x,q + N} \right)} {t_n}\left( x \right)$

(11)

基于上述理论, 对于图像 $f\left( {x,y} \right)$ , 低频快速切比雪夫矩算法如下:

(1)对图像进行UWT二维分解, 提取其低频部分并进行重叠分块;

(2)计算切比雪夫多项式和系数矩阵;

(3)计算所有的行向量. 对于图像 $f\left( {x,y} \right)$ 每行的第一个行向量, 计算其切比雪夫矩, 每行的其余行向量计算其快速切比雪夫矩;

(4)计算所有的列向量的切比雪夫矩;

(5)计算图像 $f\left( {x,y} \right)$ 第一个块的矩特征, 基于步骤(4)的结果, 采用快速切比雪夫矩计算第一行的所有块的矩特征;

(6)基于步骤(2)和步骤(4)的结果, 采用快速切比雪夫矩计算其他块的矩特征.

2 PatchMatch算法匹配

本文采用PatchMatch算法进行特征匹配, 它是一种图像块之间寻找最近邻匹配的快速随机算法, 将正确的偏移量传播并且迭代更新至全部偏移量, 相比较传统的kd-tree算法而言, 不仅可以大大减少处理时间, 还可以提供准确的匹配率^[11]. 算法步骤如下:

(1)初始化. 对于每个像素s, 随机初始化偏移量:

$\delta \left( s \right) = U\left( s \right) - s$

(12)

其中, U(s)是一个二维随机变量, 并且均匀分布在整个图像. 由于我们正在寻找与目标相对较远的匹配, 这里我们不考虑 $\delta \left( s \right) = 0$ , 以及所有偏移量小于给定阈值的情况. 虽然大多数初始随机偏移量很少用到, 但也有可能是最优的或近似最优的.

(2)传播. 在这个阶段, 图像按照从上到下、从左到右的顺序光栅扫描, 则每个像素s偏移量更新为:

$\delta (s) = \arg \mathop {\min }\limits_{\phi \in {\Delta ^P}\left( s \right)} D\left( {f\left( s \right),f\left( {s + \phi } \right)} \right)$

(13)

其中, ${\Delta ^P}\left( s \right) = \left\{ {\delta \left( s \right),\delta \left( {{s^r}} \right),\delta \left( {{s^c}} \right)} \right\}$ , ${s^r}$ 和 ${s^{\rm{c}}}$ 分别表示光栅按行和列扫描的像素s之前的像素. 每个像素判断相邻块的偏移量是否提供了更好的匹配, 如果是, 则采用相邻块的偏移量. 若给定区域像素具有恒定偏移量, 就可获得正确偏移量, 并迅速传播, 填充整个区域的下方和右侧. 每次迭代更新时, 按照反转的光栅顺序(从下到上、从右到左)扫描图像, 以获得更准确的偏移量.

(3)随机搜索. 由于传播过程依赖于随机初始的偏移量, 所以不能达到最优匹配. 为了尽量避免陷入局部最小值, 采用随机搜素, 在修正式(13)后, 对当前偏移量进行随机采样, 设候选偏移量 ${\delta _i}\left( {\rm{s}} \right),i = 1,\cdots,L$ 为:

${\delta _i}\left( {\rm{s}} \right) = \delta \left( s \right) + {R_i}$

(14)

其中, R_i是一个二维随机变量, 并且均匀分布在除去原点的半径为2^i–1的方形区域中. 事实上, 大多数的候选偏移量非常接近 $\delta \left( s \right)$ , 随机搜索后偏移量更新为:

$\delta \left( s \right) = \arg \mathop {\min }\limits_{\phi \in {\Delta ^P}\left( s \right)} D\left( {f\left( s \right),f\left( {s + \phi } \right)} \right)$

(15)

其中, ${\Delta ^R}\left( s \right) = \left\{ {\delta \left( s \right),{\delta _1}\left( s \right),\cdots,{\delta _L}\left( s \right)} \right\}$ .

3 后处理

除了平滑区域之外, 通过特征匹配获得的偏移量应该是细节化的. 但是由于噪声、压缩、几何变换、光照变化和相似区域等原因^[12], PatchMatch算法获得的偏移量很少遵循该情况, 因此后处理阶段需要: (1)对偏移量进行正则化. (2)添加适当的约束.

由于隐式过滤得到的偏移量已经足够规则, 所以需要添加适当的约束. 这里采用稠密线性拟合的方法^[13], 这是因为它复杂性低并且可以快速地得到正确的偏移量.

通过仿射模型, 在像素s的适当N像素邻域内拟合出真实偏移量:

$\hat \delta \left( {{s_i}} \right) = A{s_i},\;i = 1,\cdots,N$

(16)

转换参数A, 设置为平方误差之和的最小值.

${\varepsilon ^2}\left( s \right) = {\sum\limits_{i = 1}^N {\left\| {\delta \left( {{s_i}} \right) - \hat \delta \left( {{s_i}} \right)} \right\|} ^2}$

(17)

虽然偏移量是二维的, 但是参数A可以针对每一维进行独立地优化, 因此, 可以将 $\delta \left( x \right)$ 视为一维偏移量, 问题转化为:

${a^{opt}} = \arg \mathop {\min }\limits_a {\left\| {\delta - Sa} \right\|^2}$

(18)

其中, $\delta = {\left[ {\delta \left( {{s_1}} \right),\delta \left( {{s_2}} \right),\cdots,\delta \left( {{s_N}} \right)} \right]^{\rm{T}}}$ 是匹配阶段的偏移量, $a = {\left[ {{a_0},{a_1},{a_2}} \right]^{\rm T}}$ 是识别仿射模型的参数向量, S是邻域内所有像素齐次坐标的N×3矩阵.

$S = \left[ {\begin{array}{*{20}{c}} 1&{s11}&{s12} \\ 1&{{s_{21}}}&{{s_{22}}} \\ {}&{...}&{} \\ 1&{{s_{N1}}}&{{s_{N2}}} \end{array}} \right]$

(19)

因此真实偏移量为:

$\hat \delta \left( {{s_i}} \right) = {a_0} + {a_1}{s_{i1}} + {a_2}{s_{i2}},\;i = 1,\cdots,N$

(20)

这是一个多元线性回归问题, 解为:

$ {a^{opt}} = {\left( {{S^{\rm T}}S} \right)^{ - 1}}{S^{\rm T}}\delta $

(21)

因此, 相应的平方误差之和为:

$\begin{aligned} {\varepsilon ^2}\left( s \right)& = {\left\| {\delta - S{{\left( {{S^{\rm T}}S} \right)}^{ - 1}}{S^{\rm T}}\delta } \right\|^2} \\ & = {\left\| {\left( {1 - H} \right)\delta } \right\|^2} = {\delta ^{\rm T}}\left( {1 - H} \right)\delta \\ \end{aligned} $

(22)

其中, $H = S{\left( {{S^{\rm T}}S} \right)^{ - 1}}{S^{\rm T}}$ 是对称的, 进一步简化H矩阵为:

$H = Q{Q^{\rm T}},\;Q = \left[ {{q_1},{q_2},{q_3}} \right]$

(23)

其中, ${q^j}$ 为N行列向量, 所以, 有:

${\varepsilon ^2}\left( s \right) = \left( {{\delta ^{\rm T}}\delta } \right) - {\left( {{\delta ^{\rm T}}{q_1}} \right)^2} - {\left( {{\delta ^{\rm T}}{q_2}} \right)^2} - {\left( {{\delta ^{\rm T}}{q_3}} \right)^2}$

(24)

下面给出具体的后处理步骤:

(1)对得到的偏移量进行中值滤波操作.

(2)通过最小二乘线性模型计算圆形邻域内的阈值为 $T_\varepsilon ^2$ 的拟合误差 ${\varepsilon ^2}$ .

(3)去除误匹配对, 包括较匹配区域对的距离像素更接近的匹配对, 或小于匹配区域面积像素的匹配对.

(4)映射检测到的区域.

(5)使用圆形的结构元素进行形态学操作, 实现篡改图像区域的定位完成检测.

为了证明稠密线性拟合方法的有效性, 对篡改图像后处理阶段的进行检测, 如图1所示. 由图1可知, 匹配阶段的偏移量已经足够规则, 但还存在较多误检, 后处理阶段采用稠密线性拟合方法可有效降低误检率.

4 实验分析

实验是在Windows10操作系统下, 采用Matlab软件测试的. 为了验证本文方法的可行性和有效性, 实验采用的数据集为benchmark data^[14], 该数据集包括48幅原图像以及48幅篡改图像, 篡改区域包括单区域和多区域. 同时将本文方法与文献[5]和文献[14]进行对比来证明本文方法的优越性.

性能分析主要包括图像层面和像素层面, 本文从像素层面来衡量篡改检测算法的性能, 采用precision、recall和F这3个指标, 定义如下:

$precision = \frac{{TP}}{{TP + FP}}$

(25)

$recall = \frac{{TP}}{{TP + FN}}$

(26)

$F = 2 \cdot \frac{{precision \cdot recall}}{{precision + recall}}$

(27)

其中, precision为检测准确率, 即检测为篡改图像的数据集中有多少是真正的篡改图像, 包括把篡改图像检测为篡改图像及把真实图像检测为篡改图像; recall为检测召回率, 即数据集中的篡改图像有多少被正确检测了, 包括把篡改图像检测为篡改图像及把篡改图像检测为真实图像; F为综合评价指标. TP为篡改图像中被正确检测到的篡改像素数量, FP为篡改图像中未篡改部分被检测为篡改像素的数量,FN为篡改图像中的篡改部分未被检测到的篡改像素数量. precision和recall这2个评价指标相互制约, F综合考虑了这两个指标, F值越大, 算法的检测性能越好.

本文对篡改区域为单一区域和多区域分别进行了实验, 单区域篡改图像的检测结果见图2, 单区域多次篡改图像的检测结果见图3, 多区域篡改图像的检测结果见图4, 表1为3种算法篡改检测性能的比较.

由图2~图4可知, 文献[5]和文献[14]都会出现漏检测的情况, 尤其对于图4(c)出现两个不同的篡改区域, 而只定位出一个篡改区域, 虽然文献[14]的检测结果优于文献[5]的检测结果, 但是检测到的篡改区域还是会遗漏一些篡改区域的细节信息, 综合来说, 本文算法的检测结果是最好的. 为了更加准确地体现本文算法的优越性, 由表1可知, 本文算法的综合评价标准F是最高的, 其中precision分别提高了8.1%和1.8%, recall分别提高了3.3%和1.2%, F分别提高了3.5%和1.5%, 与上图中的检测结果一致. 这是因为文献[5]采用DWT和Tchebichef矩结合的特征向量, 比较其相似性定位篡改区域, 特征向量没有很好地表示篡改图像信息, 会存在漏检测的情况; 文献[14]对每个重叠分块提取快速切比雪夫矩特征, 特征匹配使采用kd-tree算法, kd-tree算法不能很好地实现篡改图像的匹配; 本文算法采用低频快速切比雪夫矩作为特征向量, PatchMatch算法进行特征匹配, 最后用稠密线性拟合算法消除误匹配, 这样可以较好地描述图像信息, 防止漏检和误检的情况.

我们还比较了本文算法和文献[5]和文献[14]在数据集benchmark data中每幅图片的平均运行时间, 来证明本文算法的实时性, 实验结果如表2所示. 由表2可知, 除了图3(b)之外, 其他图片都是本文算法的运行时间最短, 最后一列是数据集中每幅图片的平均运行时间, 可知本文算法的平均运行时间是最短的. 这是由于本文算法采用UWT实现图像分解, 提取低频部分的快速切比雪夫矩特征, 并且采用PatchMatch算法进行匹配, 可以缩短算法运行时间, 提高本文算法的实时性.

最后我们比较了本文算法和文献[5]和文献[14] 在数据集benchmark data中每幅图片的平均误匹配率(误匹配对/匹配对), 实验结果如表3所示. 由表3可知, 本文算法的匹配对是最多的, 误匹配对是最少的, 误匹配率是最低的. 这是因为相比较文献[5]的相似性匹配及文献[14]的kd-tree匹配, 本文算法采用PatchMatch算法进行匹配得到规则的偏移量, 并且采用稠密线性拟合方法进行后处理降低误匹配, 具有很好的匹配效果.

5 结论与展望

本文提出了一种低频快速切比雪夫矩的篡改图像检测算法. 利用非抽样小波变换分解图像并提取其低频部分, 对于重叠块提取其快速切比雪夫矩特征, PatchMatch算法匹配块特征, 最后剔除误匹配并定位其篡改区域. 相比之前算法, F达到了92.0%, 分别提高了3.5%和1.5%, 本文算法对于单区域、单区域多次以及多区域的篡改均有很好的定位结果, 并且有效地降低了运行时间, 提高了算法的实时性, 最后也降低了误匹配率实现很好的匹配效果. 后续工作将放在对篡改区域进行旋转、加噪、压缩以及模糊等处理的检测及精确定位上.