CN103544703B - 数字图像拼接检测方法 - Google Patents

Abstract

提出了一种检测数字图像是否经过拼接处理的方法，该方法首先将原始图像进行下采样分割成四个子图，利用相机图像色彩滤波阵列特点，采用非负最小二乘法获得每个子图内各像素点的预测值和对应的八维预测加权系数矢量，统计各子图归一化预测误差并据此对子图排序；其次将图像分割为多个宏块，对每个宏块估计其各色彩通道的照度，按照度比例构建照度图；然后对照度图进行树图分割，利用分割结果确定边界宏块和内部宏块，对归一化预测误差最小的子图分别计算域间、域内统计特征量；对归一化预测误差最大的子图在三个方向上构建共生矩阵，提取图像纹理特征；最后将获得的特征向量送入支持矢量机，判决图像是否进行过拼接。

Description

数字图像拼接检测方法

技术领域

检测数字图像是否经过拼接合成处理，应用于数字图像取证领域。

背景技术

随着CCD等硬件技术的发展，数码相机、手机的大规模应用，来自于它们的数字图像取代了传统胶片图像成为图像的最大来源。伴随着计算机软件的普及，这些图像的修改变得很容易，以往只有专业摄影师和修片师才能完成的工作，如今几乎人人都能借助图像处理软件完成。图像目前在法律上允许作为证据，在新闻报道中也大量采用，由此产生了问题：我们所见的图像究竟是否真实？是否被篡改过？鉴别图像的原始性即为图像取证，目前数字图像取证技术主要分二大类：主动取证和被动取证。

主动取证是在图像采集时，通过在数据中加入特殊的数据，如果最终图像中的这些特殊数据不符合一定的规律或者遭到破坏，说明图像被篡改了，数字水印技术是其中的一个代表。但主动取证的问题是它要求图像在生成的时刻就被嵌入特殊信息，而这一要求在实际应用中受到很大限制。

被动取证则不受以上限制，仅仅根据提供的图像本身的某些特性就可判别图像是否被篡改过。被动取证的基本依据是：任何形式的篡改都不可避免地引起数据统计特征的变化。目前取证研究可以采用基于像素的检测方法，通过检测像素统计异常进行判断；对图像篡改中伴随的缩放、旋转和拉伸等操作，可以看成是图像信号上、下采样，即图像发生多重采样，通过检测周围像素之间产生周期性的相关性进行判断；针对JPEG格式图像的取证，可以通过分析离散余弦变换(DCT)系数的直方图的特点，估计首次压缩时所使用的量化系数，利用二次压缩DCT系数直方图上存在的周期性噪声和直方图的傅里叶变换估计质量因子，利用Markov随机过程描述差值矩阵，根据矩阵中元素相关值的分布确定是否发生二次JPEG压缩。一个通用的取证方式是基于成像设备特性，受工作原理和物理特性的影响，数码相机的镜头、成像传感器和数字信号后处理会在成像过程中留下特有的设备痕迹和噪声，通过检查设备痕迹和噪声的一致性可判断图像是否进行了篡改，彩色图像需要经过色彩滤镜阵列计算得到，这些插值运算在图像的各个色彩通道内、像素间引入特殊的周期性的统计相关性，检测插值像素的周期相关性是否被破坏以判断图像原始性，检查图像内相机响应函数的一致性也可判别图像真伪。图像的拼接操作会破坏图像的光照一致性，可根据图像场景中光照模型的一致性鉴定图像真伪。

发明内容

在本发明用于检测数字图像是否被拼接篡改了，是一种被动取证的方法。

发明主要基于以下二个思路：

(1)、色彩滤波阵列的一致性。受成本、技术的限制，CCD在每个像素的位置上通常只有一个色彩通道，基于人眼对绿光的敏感性，色彩R、G、B(红、蓝、绿)三个分量在2×2区域内，一般分配2个绿色通道、一个红色通道和一个蓝色通道。这4个通道位置不重合，构成的排列称之为色彩滤镜阵列(CFA)，图1为CFA示意图。不同品牌相机CFA不尽相同。在实际成像中，每个像素的值由R、G、B三个值组成，那么在该位置上缺少的色彩信息由其他位置上CCD获得的该色彩值内插而成，不同品牌的相机内插公式不尽相同。

当图像中包含来自其它图像的拼接数据时，此时会出现二个问题：(a)由于拼接部分和原图部分其所用的CFA阵列、内插算法不同，当将原图区域的内插算法套用到拼接区域时，得到的结果与实际拼接部分的结果差异较大；(b)即使当拼贴来源为同一品牌相机时问题(a)不存在，另一个问题也无法忽视：拼接部分的所有边缘区域只有恰好与被覆盖原图的边缘区域的CFA排列完全一致、交接覆盖处像素值完全一致，则按照原图内插估计得到的数值才能与拼接区域内部边缘处实际值相吻合，但实际中这种情况出现的概率微乎其微。

有时为了使拼接区域与原图在视觉上更好地融合，需要对拼接部分进行平滑、旋转、伸缩变形等处理，这些处理也使得CFA内插预测得到的值与实际值之间产生差异。以上这些由于图像拼接造成的内插异常为图像取证提供了证据。

(2)、施加在物体上的照度与光源的性质、角度等因素有关，图像拼接部分与其它部分照明度难以一致。

本方法首先将原始图像进行2：1下采样分割成四个图像子图，在红色或蓝色通道内采用非负最小二乘法对每个子图各像素点进行预测，得到预测值以及对应的预测加权系数矢量，统计各子图归一化预测误差并根据结果对子图进行排序；其次将原始图像分割为多个宏块，依据色彩反转强度空间一致性在每个宏块内独立进行各色彩通道的照度估计，并依据各色彩照度比例构建照度图；对照度图进行基于树图的分割，根据分割结果确定边界宏块和内部宏块，对照度图进行下采样；对归一化预测误差最小的图像子图，求其分割域内、域间特征统计量；对归一化预测误差最大的子图求三个方向的共生矩阵，提取纹理特征统计量；将所有提取的特征送入经过训练的支持矢量机，得到鉴定结果。

附图说明

包含三个附图。

图1为CCD色彩滤镜阵列示意图。

图2为对图像2：1下采样示例。

图3为图像分割为宏块示例。

具体实施方法

(一)将原图进行2：1下采样，得到的四个图像子图，图2为下采样示例。由于绿色通道的数据相关性较大，预测误差小，难以区分，因此要求采用红色或蓝色通道的数据，对蓝色或红色通道中任选一个通道的数据进行如下处理。对蓝色、红色通道而言，4个图像子图中必有一个的像素值是CCD采集直接得到的，其他三个图像子图中所有像素值是对前者内插得到的，因此它们内部各像素值之间相关性比前者强。

(二)非负最小二乘法预测，获得内插线性系数矢量和预测误差。每个图像子图在所选色彩通道内，采用非负最小二乘法对像素点由其相邻的8个点的值对其进行预测，记录每个像素对应加权系数组成的8×1维矢量。

设为图像子图中的像素点值，它的8个相邻像素点值分别为，…,，根据8个相邻像素值可以得到对的预测值＝，其中为加权系数。在一个M×M的区域内，这里M为奇数，一般取值5、7、9，共有个像素，构成的矢量为p＝[]^T,对应的非负加权系数矢量,以每个像素的8个邻像素为一行、个像素共构成×8的邻像素矩阵G，则采用非负最小二乘的方法求得令取得最小值的矢量u。

上述预测公式是线性的。在实际CFA内插操作中，内插算法不一定是线性的，但非线性公式可以根据泰勒公式展开，取其一阶导数项，舍去非线性展开的部分，将其近似为线性。理想情况下，若图像摄入过程无噪声、图像未被篡改，则同一图像子图中各像素点对应的8×1维加权系数矢量应该相同。若图像摄入有噪声，则在图像未被篡改的条件下，同一子图内各像素点对应的8×1维加权矢量非常接近。当图像被篡改后，由于内插方式的改变、被打断，在同一子图的某些像素对应的8×1维矢量与其它像素点对应的矢量有显著差异。

(三)像素点预测误差为：，分别计算各子图内预测误差之和，除以各像素值之和，得到各图像子图归一化后的预测误差。将四副图像子图归一化预测误差按升序排序，排序中的最后一个图像子图以A表示，排序中的第一个图像子图以D表示，第二个称为C，第三个称为B。

(四)将原始图像分割成若干宏块，在各宏块内部依据色彩反转强度空间一致性进行照度估计，细节如下：

●将原始图像分成大小为N×N(N为偶数，一般N取32、64或128)的宏块，各宏块相连、互不重叠，如图3所示。宏块内每个像素点对应红、绿、蓝色彩值分别用、和表示。

●对宏块内各像素分别求：

和，

计算，其中c分别表示红、绿、蓝色彩

●求该宏块内所有、中的最大值，分别用、表示

●选择宏块内照明充分的像素，所谓照明充分指该像素对应的、同时满足：

和

●在宏块内部的每一个色彩通道，对每个照明充分的像素做变换，计算，(c分别代表红、蓝、绿三个色彩)，分别以x,为横、纵坐标，将原始图像由r、g、b空间映射到x,构成的色彩反转强度空间

●将色彩反转强度空间用霍夫变换映射到霍夫空间，横轴用(c分别代表r、g、b色彩)表示，纵轴用p表示

●在霍夫空间进行投票表决：如果霍夫空间的一个点属于一条直线，则在以该点对应横坐标为索引的数组中加1；若一点属于m条直线的交汇点，则以该点所对应横坐标为索引的数组中增加m。历遍霍夫空间所有点。

●找出以为索引的数组中的最大值，该值相应的即为此宏块对色彩c的照度估计。理论上得到的三色照度估计之和应为1，即：，但实际会有少量偏差，和不一定为1。

(五)构建图像照度图，以宏块为单位组成，每个宏块内部色彩相同，各宏块内部像素值统一用：(c分别代表红、蓝、绿色彩通道)表示。得到的照度图分辨率与原始图像相同。

(六)对照度图进行树图分割。即使对原始图像，得到的各宏块照度估计也不尽相同，究其原因照度估计是很容易受到干扰的，需要对照度图进行分割，具体细节如下：

●以一个宏块为一个顶点构建树图，所有顶点集合用V表示，E为树图所有边的集合，两个相邻顶点、之间的连线构成一条边(，)E。

●此时构建的树图为q＝0时的树图。计算每条边的权重，其中分别为宏块i、宏块j的照度值矢量，为欧拉范数，权重反应了顶点间的差异程度，该值越大说明二宏块照度差异越大。

●对所有边的权重按照升序进行排列，假设树图共有r条边，令q依次为1，2，……，r,重复下列(a)-(c)步骤r次，直到q＝r为止。

(a)根据第q-1次的分割结果，迭代获得第q次的分割结果：设在第q次中处理的是排序中第q条边连接的二个顶点v_i、v_j，若v_i、v_j仍分属不同的区域(用和表示)，那么分别计算区域各自内部照度矢量平均值并更新经q-1次合并后该边的权重

(b)计算和各自区域内部的不相似性，该值为各自区域内部所有边的权重最大值表示，此处用表示。

(c)比较和

●若满足：

表明迭代后二个区域间不相似度大于任意一个区域内部不相似度的函数，说明二个区域严重不相似，则保持和区域各自独立；

式中Δ(z_i,z_j)＝k/dis(z_i,z_j)/|z_i|

dis(z_i,z_j)为二个区域中心点的距离，|z_i|为区域z_i的面积，这里指该区域包含的顶点数目，k为常数，取值范围为400-600之间的整数，省缺值k＝500。当二个区域距离远时dis(z_i,z_j)值大，不等式右侧变小，因此不等式成立的条件容易满足，这可解释为当图像中二个物体距离较远时，它们的照度估计差别较大是可以理解的。当|z_i|较小时，不等式右侧较大，不等式较难满足，即对小区域要求更强的边界。

●否则合并和区域，重新计算合并后新组成区域的平均照度矢量、区域内部边权重的最大值、区域面积，更新区域与其相邻各区域的边权重。

(七)根据照度图分割结果，确定照度图中的边界宏块：若宏块的相邻宏块中存在与其不属于同一区域的，则该宏块称为照度边界宏块。否则宏块称为照度内部宏块。

(八)将照度图采用与步骤(一)相同的2：1下采样方式，得到4个照度子图，分别与4个图像子图相对应。4个照度子图中对应于照度图的宏块称为照度子宏块，则照度子宏块大小为(N/2)×(N/2)。照度图中照度边界宏块对应于照度子图中的照度边界子宏块，其余宏块对应的子宏块称为照度内部子宏块。在图像子图中，与照度边界子宏块对应的区域称为图像边界子宏块，否则称为图像内部子宏块。

(九)进行域内特征统计，得到四个8×1维特征统计量：对步骤(三)确定的图像子图D，根据其相应照度子图，计算每个区域内部所有图像内部子宏块中像素点对应的8×1维加权系数矢量统计特征，矢量中第i个元素u_i(i＝1,2,…8)的特征如下：

●归一化均值：

●归一化方差：

其中TM为该子图中，此内部区域中包含的总像素数目，为(N/2)×(N/2)的整数倍。

若照度子图中形成多个区域，则在对应图像子图中分别求出这些区域内部μ_i-I、σ_i-I，设当前区域为j，则计算求出所有区域中L最小的那个区域，其对应的σ_i-I和μ_i-I分别记为min(σ_i-I)和min(μ_i-I)。同理L最大的区域，其对应的σ_i-I和μ_i-I分别记为max(σ_i-I)和max(μ_i-I)，如此得到4个特征量，每个特征统计量是一个8×1维矢量。

(十)域间特征统计，得到二个特征统计量。对步骤(三)确定的图像子图D，其中二二相邻、不属于同一区域的图像边界子宏块构成图像边界子宏块对。不妨设其中一个图像边界子宏块属于区域e，另一个属于区域f，对二宏块中各像素所对应的8×1维加权系数矢量，计算矢量中第i个元素u_i(i＝1,2,…8)的特征如下：

●域间归一化均值：

●域间归一化方差对比度：

其中分别为分属区域e、f的图像边界子宏块各自内部像素对应u_i归一化均值，分别各自内部归一化方差。历遍D中所有图像边界子宏块对，求出为最大值的图像边界子宏块，它们所对应的μ^(e,f)、S^(e,f)为本步骤得到的二个特征量，每个特征量为8×1维矢量。

(十一)提取纹理特征参考向量。以上步骤提取的特征值与图像纹理特征也有关系，当图像纹理特征丰富时，原始图像上述特征与纹理特征较少的伪造图像类似，因此在评判图像是否伪造时需要将纹理特征因素考虑进去。将步骤(三)中归一化预测误差最大的图像子图A转化为灰度图像，分别在0°、45°和90°3个方向上求图像的共生矩阵，为减少数据量，采用量化系数16进行量化，这样得到3个16×16的共生矩阵，计算每个共生矩阵的10个特征量：

●加权平均：μ＝∑_i∑_jic_ij，其中i和j为共生矩阵的坐标，c_ij为坐标(i,j)对应的共生矩阵中的元素。

●加权方差：σ²＝∑_i∑_j(i-μ)²c_ij

●矩阵能量：

●熵:-∑_i∑_jc_ijlogc_ij,若c_ij为0，则将c_ijlogc_ij设为0。

●对比度:∑_i,j(i-j)²c_ij

●相关性：

●聚类影：∑_i,j((i-μ)+(j-μ))³c_ij

●二阶矩：

●差分倒矩

●聚类凸值：∑_i,j((i-μ)+(j-μ))⁴c_ij

这10个特征量称为纹理特征量。三个矩阵共得到30个特征量。

(十二)构建建立训练集合，进行训练。从若干已知的真实图像和拼接图像中提取上述特征量，哥伦比亚大学、达特茅斯大学等均有此类图像数据库可供使用，也可自行构建图像数据库。采用以径向基函数为内核的支持矢量机进行训练，注意本步骤在训练完成后即可跳过。

(十三)用支持矢量机进行鉴定。将步骤(九)得到的32个数据(4个8×1维矢量＝4×8×1＝32)、步骤(十)得到的16个数据(2个8×1维矢量＝2×8×1＝16)、以及步骤(十一)得到的30个纹理特征量，组成共78(32+16+30＝78)个特征量，送入训练完毕的支持矢量机进行判决。

用以上步骤提取的特征向量，能够准确检测出图像是否经过了拼接篡改。

Claims (2)

1.一种检测数字图像是否经过拼接合成的方法，其特征在于包含以下几个步骤：

A. 将原图像进行2：1下采样，分割成四个图像子图；

B. 利用相机图像色彩滤镜阵列特点，在红色或蓝色通道内对每个图像子图各像素点进行预测，采用非负最小二乘法由其八个相邻点像素值得到该像素点预测值以及对应的预测加权系数矢量；

C. 统计各图像子图归一化预测误差，并据此对四个图像子图进行排序；

D. 将原始图像分割成多个互不重叠的矩形宏块，依据色彩反转强度空间一致性在每个宏块内独立进行各色彩通道的色彩照度估计，并依据各色彩照度比例构建照度图；

E. 对照度图进行基于树图的分割，构建树图以一个区域为一个树图顶点，相邻区域构成的二个顶点由边连接、边的权重为二个区域各自内部平均照度矢量的欧拉范数，并将所有边的权重排序后依次处理，初始时每个宏块作为一个区域，采用迭代更新的方式，每次历遍当前状态下所有二二相邻的区域，判断是否需要将二个相邻区域合并成一个新区域，合并准则为这二个相邻区域边的权重小于门限则合并，反之则保持各自区域独立，门限取决于这二个相邻区域各自内部所有边权重的最大值、二个区域中心点间距离、二区域各自的面积；

F. 依据步骤E分割结果判决照度图中各宏块性质，若一宏块的所有相邻宏块中有与其不属于相同区域的，则该宏块称为照度边界宏块，否则该宏块称为照度内部宏块；

G. 对步骤D得到的照度图进行与步骤A中相同的下采样得到四个照度子图，如此得到的四个照度子图分别与步骤A得到的四个图像子图相对应，照度图中宏块对应照度子图中的照度子宏块，照度图中照度边界宏块对应于照度子图中的照度边界子宏块，照度图中照度内部宏块对应于照度子图中的照度内部子宏块;

H. 在与照度子图相对应的图像子图中，与照度边界子宏块对应的范围称为图像边界子宏块，与照度内部子宏块相对应的范围称为图像内部子宏块；

I. 对步骤C中归一化预测误差最小的图像子图，历遍其所有区域，对每个区域由该区域包含的所有图像内部子宏块中像素点的预测加权系数矢量计算出该区域的域内归一化方差矢量，分别找出域内归一化方差矢量中各元素之和为最小、最大的二个区域，以这二个区域的域内归一化均值矢量和域内归一化方差矢量为特征量；

J. 对步骤C中归一化预测误差最小的图像子图，任何二个相邻、且不属于同一区域的图像边界子宏块构成一个图像边界子宏块对，历遍该子图所有图像边界子宏块对，由一个图像边界子宏块对所包含各像素点对应的预测加权系数矢量，计算该图像边界子宏块对的域间归一化均值矢量和域间归一化方差对比度矢量，找出域间归一化方差对比度矢量中各元素之和最大的那个图像边界子宏块对，以该图像边界子宏块对的域间归一化均值矢量和域间归一化方差对比度矢量作为特征量；

K. 对步骤C中归一化预测误差最大的图像子图，求0度、45度和90度三个方向的共生矩阵，计算矩阵纹理统计特征值作为特征量；

L. 将步骤I、J、K所求得的特征量送入经过训练的支持矢量机，得到鉴定结果。

2.根据权利要求1所述方法，在步骤E中在判决是否需要合并二个相邻区域时，其门限计算方法如下，首先分别找出这二个区域各自内部所有边的权重之最大值，然后分别将各自的最大值与各自修正值相加求和，选取二个和值中较小的和值做为门限，其中二个区域各自修正值的计算方法为常数除以这二个区域中心点间距离、再除以本区域的面积，这里常数值取值范围为400-600间整数，缺省值设为500。