基于JPEG压缩的数字化地质资料篡改探测方法及应用

刘助龙，赵于前，廖 苗，张竣凯，戴塔根

(1. 中南大学 有色金属成矿预测教育部重点实验室，长沙 410083；2. 中南大学 地球科学与信息物理学院，长沙 410083)

地质资料是地质工作形成的重要知识财富，为我国社会经济发展提供了重要的支撑。截止到2006年6月底，全国地质资料馆已收藏了包括区域调查、矿产、石油和天然气、海洋、水工环(灾害地质)、物化遥、信息技术、科研等在内的各类纸质成果地质资料 10万余种；1999年，实施地质大调查计划以来又产生了2000余种成果地质调查资料[1−3]。国家及省级地质资料馆开展了长达10 a的地质资料图文数字化工作，全国地质资料馆已在网上提供了5 000余份资料的全文在线浏览服务。为充分发挥海量地质资料的潜在巨大作用，2009年以来国土资源部全面部署了地质资料信息服务集群化和产业化工作，主要目标是对地质资料信息进行集群开发，通过统一的共享服务平台，为经济社会发展提供系列、权威、集群化的地质资料信息产品。按照统一部署，全国地质资料馆、实物中心、油气中心、部信息中心、经研院等有关单位和上海、山东、青海、湖北、湖南5个省份开展了试点工作。地质资料电子化服务有了新的进展，信息服务面不断拓展，服务质量不断提高，最大限度地满足了社会经济发展对地质资料提出的多层次、多范围和多样化要求，提高了地质资料的利用率。

全新的地质资料服务模式下，已经构建了全国、区域和省的三级地质资料信息服务集群体系，同时开发了共享服务产品。然而，在这些资料通过互联网面向公众服务的过程中，采用Photoshop、Acron等图像编辑软件，就能很容易地对一些数字图像资料进行篡改，但很难肉眼辨别。如果一些被蓄意篡改的地质资料传播出去，可能会对相关工作带来非常严重的后果。如何判断为政府、社会单位和公众提供的地质资料信息产品是权威的、没被篡改过的，同时为使用者提供便捷的图像真实性检测方法，已成为服务过程中面临的一个重要问题。

目前，中国地质调查局发展研究中心组织开展了数字水印技术[4−6]在地质数据成果汇交、分发等过程中实现知识产权保护的应用研究。数字水印技术是一种主动检测技术，由于其必须提前向原始数字图像中嵌入验证信息，所以在一定程度上，主动取证技术会受到应用条件的限制，无法从根本上遏制图像篡改的发展，因而现在更注重被动取证技术[7−10]研究。

本文作者介绍了几种基于JPEG(Joint photographic experts group)格式图像篡改的被动检测算法，并将其应用于地质资料的篡改检测。

1 JPEG压缩原理

数字地质资料图通常是尺寸较大的图像，如果使用无损图像格式保存，将会给资料的储存和传输带来诸多不便。JPEG压缩算法使用一种有损的压缩方式来处理图像，当使用适当的压缩因子压缩图像时，其失真的程度不至于影响图像文字和纹理等信息的传播，即可以同时得到较好的图像质量并占用较小的存储空间，这也就是JPEG图像之所以广泛地应用于地质资料数字化格式的主要原因。本文作者首先简单介绍JPEG压缩的基本原理。

图1 JPEG压缩流程图Fig. 1 Flow chart of JPEG compression

JPEG压缩流程如图1所示，主要包括以下6个步骤：1)图像空间转换；2)缩减采样；3) 8×8分块；4)DCT变换；5)量化；6)熵编码。首先将源图像从RGB(红绿蓝)空间转换到 YUV(亮度色度)空间，然后将色差和色度进行空间缩减采样，再将源图像分成连续不重叠的 8×8像素块，并对每个块进行离散余弦变换(Discrete cosine transform，DCT)，得到64个DCT系数。习惯上将位于左上角的第一个值称为 DC(Direct current)系数，其余 63个系数称为 AC(Alternating current)系数。接着，用由 64个整数值的量化步长组成的8×8的量化表来量化64个DCT系数，并对量化结果进行四舍五入取整。通过使用一个合适的质量因子(Quality factor, QF)，可以对图像质量和压缩比率进行调整，且标准JPEG压缩的QF值为1到100的整数，其中QF=100表示最小的量化步长，因而可以得到最好的图像质量和最小的压缩比率，QF=1则反之。图2(a)和(b)分别表示标准质量因子QF=100和QF=50的量化表。最后，量化后的 DCT(QDCT)系数通过运用熵编码编码成比特流。先使用Z字形编码把QDCT系数变为一维，再使用霍夫曼编码(Huffman coding)生成二进制比特流。

读取图像的过程是解压缩过程，只需要把以上流程反向执行。需要注意的是由于在量化过程中存在着取整，因此量化是不可逆的，这也是JPEG压缩是有损压缩的原因。

2 几种基于 JPEG 压缩的篡改探测算法

2.1 基于BAG的JPEG图像篡改检测算法

JPEG压缩操作在图像中引入了水平和竖直的网格，称为块伪影栅格(Block artifact grids，BAG)[11]。当篡改者对图像进行复制−粘贴操作时，被粘贴过去的图像块的 BAG通常与原来图像中的 BAG不匹配(匹配概念为 1/64)，如图 3所示。图 3(a)和(b)所示是两幅原始图像，图中的每个网格代表一个8×8分块，图 3(a)的虚线方框内区域为被复制区域，图 3(b)所示的虚线方框内区域为被选中的覆盖区域，图 3(c)为将图 3(a)所示的虚线方框内区域复制粘贴至图 3(b)所示的虚线方框内区域后得到的合成图像。图3(a)～(c)中实线为块伪影栅格。从合成图像很容易看出被粘贴区域与原始区域的栅格不匹配。当篡改者对图像进行喷绘[12−13]或者类似于喷绘的操作时，被提取的 BAG图像中常出现大片无栅格区域或者栅格无规则排列区域。因此，可以通过检测 BAG图像中的栅格不匹配或空白栅格现象来检测图像的篡改并对被篡改区域进行定位。该方法的局限性是，当被篡改后的图像再次以JPEG格式保存时，由于引入了新的JPEG栅格而导致无法检测原始图像的BAG，该方法失效。

图2 JPEG 量化表: (a) QF=100; (b) QF=50Fig. 2 JPEG qualification tables: (a) QF=100; (b) QF=50

图3 BAG不匹配示意图: (a) 被复制区域的原始图像栅格; (b) 被粘贴区域的原始图像栅格; (c) 合成图像栅格Fig. 3 Schematic diagrams of BAG mismatching: (a) Original BAG image of copied region; (b) Original BAG image of pasted region; (c) BAG of composite image

BAG图像的提取方法由以下几个步骤组成。

1) 从原始图像中提取弱水平方向的边缘。

图像中的 BAG被视为一种弱边缘，通过计算二次差分绝对值可以提取出这种弱边缘。设 f(x, y)为待检测图像中(x, y)位置的灰度值，二次差分绝对值为

为了减弱图像中强边缘对弱边缘的影响，该方法使用中值滤波器来减少强边缘的干扰，并将所有大于50的d(x, y)值设置为0。接着，对水平方向每33列(前后各两个DCT块)的d(x, y)求和以增强弱水平方向边缘

最后，采用式(3)来均衡化弱水平方向边缘：

式中：右边第二项为增强后弱水平方向边缘的局部均值。

2) 在弱水平方向边缘图像中以 8为周期提取水平方向BAG。

BAG是一种出现在8×8像素块边缘的周期信号，同时它也是一种局部信号，只与周围4个像素块有关。因此，通过使用8为周期的中值滤波器来增强弱边缘图像e，最终得到水平方向BAG：

3) 提取竖直方向BAG以及获得最终BAG结果图。

通过与步骤(1)和(2)类似的方式，可以提取出竖直方向 BAG图像 gv。将 gh和 gv相加，便得到最终的BAG图像

然而，从一幅BAG图像中直接观察BAG不匹配区域的过程十分麻烦且带有很大的主观性。为了使检测结果更直观，对提取到的BAG图像中BAG不匹配区域进行标记。首先，将提取的BAG图像分成8×8不重叠块，对于其中每一个8×8块A，如果它位于一个不匹配位置，那么，在其6×6中央区域将出现错配的栅格线(BAG line)。可以用b作为块BAG位置的标记符号，将块A的6×6中央区域的6横行的每行求和后得到的最大值与其6竖列的每列求和后得到的最大值相加，并减去块A横、竖边缘各自求和后得到的最小值，计算公式如下：

式中：函数max[A{}]和min[A{}]返回矩阵A{}中的最大或最小值。最后，将BAG图像中所有8×8块由式(6)计算得到最终的标记结果图。

2.2 基于BPPM的JPEG图像篡改检测算法

与单压缩不同，JPEG双压缩是指对图像进行连续两次压缩，第一次的压缩因子为QF1，第二次的压缩因子为QF2，且第一次与第二次压缩的栅格是相匹配的。因为双压缩包含两次量化过程，这种过程使JPEG图像离散余弦变化后，DCT系数的直方图产生周期的波峰波谷现象，这种现象被称为双量化效应[14]。为了直观说明这一现象，用Matlab随机生成一组均值为0、标准差为10且包含10 000个数值的正态分布数据集s。图4(a)所示为数据集s的直方图。可以看出，它基本呈现连续分布；图4(b)所示为对该组数据先经过量化步长2量化，再经过量化步长3量化后的直方图。可以看到，有明显的峰谷周期产生；图 4(c)所示为对该组数据先经过量化步长3量化，再经过量化步长2量化后的直方图。可以看到，同样产生了峰谷周期现象。

图 4 双量化效应阐述: (a) 正态分布数据集 s直方图;(b) 数据集s先后经量化步长2和3量化后的直方图; (c) 数据集s先后经量化步长3和2量化后的直方图Fig. 4 Elaboration for double quantization effects:(a) Histogram for set of normal distribution data s;(b) Histogram for data s quantized by quantization step 2 and then quantization step 3; (c) Histogram for data s quantized by quantization step 3 and then quantization step 2

假设原始图像是一幅JPEG图像，将其篡改后再次以JPEG格式保存。所得图像中未篡改区域经历了两次JPEG压缩，因此，可以在其中检测到双量化效应，而篡改区域一般不会出现双量化效应，这是因为：1)篡改区域经喷绘产生或从无压缩图像中复制−粘贴过来的，所以，只经历了一次 JPEG压缩；2)篡改区域来自另一幅JPEG图像，或是从原图中另一区域复制过来，但是被复制区域的栅格与原图像的栅格不匹配(与2.1节情况类似)，此情况不满足产生双压缩的条件。因此，可以通过检测图像中双量化效应的缺失来进行图像篡改检测和定位篡改区域。

该方法由以下几个步骤组成。

1) 提取图像中所有 8×8块同一频率的 DCT系数，并生成64个DCT系数直方图。

2)估计每一个DCT系数直方图的周期。

设s0是直方图中柱(bin)最大值所对应的索引值，smax和 smin分别是直方图中索引的最大值和最小值。对于每个介于1和smax/20的周期p，本文作者可以计算下面的值：式中：imax=[(smax−sp)/p]，imin=[(smin−sp)/p]。H(p)定义的是当周期取p时收集波峰的能力。可以观察到，当p和直方图的周期相等时，能够得到最大的H(p)，此时直方图的周期phist=argmaxpH(p)。

3) 计算编码块后验概率图。

由于原始区域会产生双量化效应而篡改区域不会，因此，属于原始区域的 DCT块会累积到直方图的波峰中去，而篡改区域的块则会随机地累积到直方图的每个部分中。

设一个周期是从第s0个柱开始，在第s0+p−1个柱结束，那么原始区域块累积到第s0+i个柱的概率为

式中：h(k)表示DCT系数直方图中索引为k的值。

同样，篡改区域块在这个周期中累积到第s0+i的个柱概率为

根据贝叶斯估计，如果一个块累积到第s0+i个柱，那么该块是一个篡改块或者原始块的后验概率可以分别表示为

可以通过以上公式计算出图像中每个 DCT块是否篡改的后验概率，得到一个编码块后验概率图(Block posterior probability map, BPPM)[14]，图中的每一个像素点代表的是检测图像中的一个DCT块，像素的值代表的是这个DCT块的累积后验概率。从这个图中就能分辨出篡改区域。

2.3 基于ASAD的JPEG图像篡改检测算法

设被检测图像为f(x, y)，用JPEG质量因子Q对其进行压缩后的图像为fQ(x, y)，可以通过计算被检测图像与压缩后图像所对应像素差的绝对值之和的均值(Averaged sum of absolute difference, ASAD)[15]的方法来检测篡改区域。该方法对于JPEG原始图像进行喷绘，复制−粘贴，填涂等篡改后，并以JPEG格式和无压缩格式保存的图像检测都有效。具体步骤如下。

1) JPEG压缩。

以质量因子QF对待检测图像f进行JPEG压缩，得到压缩图像fQ。

2) 计算该压缩图像与原被检测图像的每个像素的差的绝对值其中：f(x, y)表示点(x, y)的像素值。

3) 由于差的绝对值图像 FQ(x, y)对比度不够明显，以至于很难从差值图像中分辨出篡改区域，因此，以点(x, y)为中心计算大小为(2b+1)×(2b+1)的像素块的差值绝对值之和，并求其均值，公式如下：

通过后面的实验将进一步发现，b值越大，检测到的篡改区域越明显，但是篡改区域的边界将变得越模糊且耗时更多。

4) 如果需要，再以不同的JPEG质量因子Q压缩待检测图像，重复以上步骤。

一般情况下，如果ASAD图像中存在一部分区域明显亮于或者暗于图像其他区域，则认为被检测图像的对应区域是篡改区域。

对于彩色图像，可以对其RGB 3个色彩空间分别应用式(12)和(13)，然后计算其均值来检测篡改区域。因此，式(12)和(13)可以分别用式(14)和(15)代替；

式中：i=1, 2, 3表示RGB 3个色彩空间；fQ(x, y, i)表示以质量因子QF压缩的待检测图像的第i个色彩空间的像素值；f (x, y, i)表示待检测图像的第i个色彩空间的像素值。

3 实验结果与分析

图 5(a)所示为 JPEG 质量因子为 50、大小为617×720的钻孔剖面图的一部分，该部分图中钻孔并没有打到矿体。图 5(b)所示为图 5(a)的篡改图并以无压缩的TIF格式保存，篡改后的图像增加了一处矿体(即图 5(b)中的网格填充区域)，从而使该地区由无矿成为有矿，造成信息失真。图 5(c)所示为对图 5(b)进行 BAG提取后得到的检测结果图。图 5(d)所示为对图 5(c)进行标记的最终结果图。可以看到，篡改区域虽然能被大致检测出来，但同时出现大量噪声影响检测结果。图5(e) 所示为以质量因子50压缩待检测图像，且b=8时计算得到的ASAD图像，图5(f)所示是对图 5(e)进行分割的结果，从图 5(e)中很容易判定待检测图像被篡改过，并定位篡改区域。很明显，ASAD算法的检测结果优于BAG提取算法的检测结果。

图 6(a)所示为矿区激电异常与实际成矿比对图，它是质量因子为70、大小为572×758的JPEG图像的一部分。该部分图中黄色条形区域为经激电中梯扫面形成的两组走向为 70°的 Dη1、Dη2异常带。该异常呈多个椭圆形，包在一起，北东向展布，规模较大，连续性好。因此，圈定此激电异常带为成矿靶区。经钻孔验证，可确定矿体存在的位置与圈定的靶区基本相符。图6(b)所示为采用复制−粘贴篡改方式，将矿体下移，然后将所得图像保存为质量因子为90的JPEG图像。篡改后，矿体物理位置偏离靶区至东南方向300 m处，对后期的勘探开采工作形成极大误导。图6(c)所示为采用BPPM算法得到的检测结果。可以看出，其中对应于篡改区域的位置明显亮于其他区域。图 6(d)所示为图 6(c)的分割结果，但是其中仍存在一些噪声，影响对于篡改区域的判断。图 6(e)所示为以JPEG质量因子70压缩待检测图像且b=8计算得到的ASAD图像，图6(f)所示为对图6(e)进行分割的结果。通过图6(d)和(f)可以对篡改区域作出较准确的判定。也可以看出，ASAD方法的检测结果明显优于BPPM方法的检测结果。

图5 经填涂方式篡改并以TIF格式保存的篡改图像检测: (a) 原始图像; (b) 篡改图像; (c) 对(b)进行BAG提取得到的检测结果图; (d) 对(c)进行标记的结果图; (e) ASAD图; (f) 对(d)进行分割的结果图Fig. 5 Detection of filling tampered image saved at TIF format: (a) Original image; (b) Tampered image; (c) Extracted BAG image form (b); (d) Marked BAG image of (c); (e) ASAD image; (f) Segmentation result of (e)

图6 经复制−粘贴方式篡改并以JPEG格式保存的篡改图像检测: (a) 原始图像; (b) 篡改图像; (c) BPPM图; (d) (c)的分割结果; (e) ASAD图像; (f) 对(e)进行分割的结果图Fig. 6 Detection of copy-paste tampered image saved at JPEG format: (a) Original image; (b) Tampered image; (c) BMMP image;(d) Segmentation result of (c); (e) ASAD image; (f) Segmentation result of (e)

图7 经填涂方式篡改并以JPEG格式保存的篡改图像检测: (a) 原始图像; (b) 篡改图像; (c) BPPM图; (d) 图(c)的分割结果;(e) ASAD图像; (f) 图(e)的分割结果Fig. 7 Detection of filling tampered image saved at JPEG format: (a) Original image; (b) Tampered image; (c) BMMP image;(d) Segmentation result of (c); (e) ASAD image; (f) Segmentation result of (e)

图 7(a)所示为 JPEG 质量因子为 75、大小为693×1 148的某矿区区域地质图图例部分。图中锰矿为沉积变质型矿床，锰矿层主要赋存于震旦下统莲沱组上部地层中，可确定此地层为主要找矿标志。图7(b)所示为对图 7(a)进行填涂方式篡改的区域地质图的图例，且篡改后的图像以JPEG质量因子85保存。将图7(a)中图例Za更改为震旦系上统后，如利用此图进行详查、勘探等工作，就不会将震旦下统莲沱组上部地层作为重点找矿区。图7(c)所示为采用BPPM算法得到的检测结果，图 7(d)所示为图 7(c)的分割结果。图7(e)所示为以质量因子75压缩待检测图像且b=8计算得到的ASAD图像；图7(f)所示为对图7(e)进行分割的结果。可以看出，图7(d)和(f)检测结果理想。

本文作者运用两个衡量指标来评价检测结果。第一个指标是检测出的篡改区域 A1与真实篡改区域 A2的覆盖率(Overlap, OL)：

第二个指标为误检率(Detection error, DE)，定义如下：

式中：W1表示将未篡改区域判定为篡改区域的像素个数；W2表示将篡改区域判定为未篡改区域的像素个数；TR表示真实篡改区域的像素个数。OL值越大、DE值越小，检测结果就越好。对图6和7中两种不同算法的全图检测结果分别用OL和DE进行评价，结果如表1所示。

对于ASAD算法，b值过小计算得到的ASAD图像对比度不明显，而b值过大篡改区域的边界将变得模糊，因此，需通过实验选取最优的b值。在图6和7中，取不同的b值，分别计算OL和DE值，得到的结果如表2所示。可以明显看出，b等于8时检测结果最优。

表1 图6和7的 ASAD算法与BMMP算法检测结果评价比较Table 1 Detection results comparisons between ASAD algorithm and BMMP algorithm for Figs. 6 and 7

表2 图6和7的ASAD算法不同b值计算得到的检测结果评价比较Table 2 Detection results comparisons of ASAD algorithm with different values of b for Figs. 6 and 7

在Intel Pentium PC(2.93GHz CPU,512MB 内存)的Matlab 7.1平台上，对于图6和7，比较以不同b值计算出4幅ASAD图像与计算出BMMP图像所需的时间，比较结果如表3所示。可以很明显看到，随着图像大小与b值的增大，本研究所提出的方法耗时将明显增加。此外，比较两种方法的实验结果和耗时，可以得到如下结论：BMMP算法比 ASAD算法计算更简单更快。

表3 ASAD算法计算得到4幅ASAD图像与BMMP算法的耗时比较Table 3 Detecting time comparison between ASAD algorithm generating 4 images and BMMP algorithm (s)

4 结论

1) 阐述了3种基于JPEG格式的图像篡改检测方法，为数字化成果地质资料在共享服务过程中可能出现的被篡改现象探索了全新的篡改探测方法，对海量数字化成果地质资料版权保护和真实性检测具有极高的理论和实用价值，从而进一步服务于地质资料信息集群化产业化工作。

2) 基于JPEG篡改探测的3种方法能够探测经不同方式篡改后保存为JPEG格式或其他无压缩格式的数字化成果地质资料图像，并能准确定位被篡改区域。

3) 相比数字签名和数字水印鉴别技术，基于JPEG篡改探测方法不需要数字化成果地质资料提供方对图像进行预处理(提取签名或嵌入水印)，操作简单方便，成本较低，应用前景广泛。

[1] 尚 武, 杨东来, 李景朝, 姜作勤. 中国地质信息服务体系的现状、差距及对策[J]. 中国地质, 2007, 34(4): 730−736.

SHANG Wu, YANG Dong-lai, LI Jing-chao, JIANG Zuo-qin.Gap and countermeasures of the geoinformation service system of China [J]. Geology in China, 2007, 34(4): 730−736.

[2] 周进生. 关于成果地质资料社会化服务的理性思考[J]. 资源与产业, 2007, 9(6): 119−121.

ZHOU Jin-sheng. Views on public Service of archived geological data [J]. Resources & Industries, 2007, 9(6):119−121.

[3] 姜作勤, 马智民, 杨东来, 李景朝, 尚 武, 王 群. 地质信息服务体系框架研究[J]. 中国地质, 2007, 34(1): 173−178.

JIANG Zuo-qin, MA Zhi-min, YANG Dong-lai, LI Jing-chao,SHANG Wu, WANG Qun. Framework of the geological information service system [J]. Geology in China, 2007, 34(1):173−178.

[4] ZHANG F, ZHANG X H, ZHANG H B. Digital image watermarking capacity and detection error rate [J]. Pattern Recognition Letters, 2007, 28(1): 1−10.

[5] CELIK M, SHARMA G, SABER E. Hierarchical watermarking for secure image authentication with localization [J]. IEEE Transactions on Image Processing, 2002, 11(6): 585−595.

[6] KUNDUR D, HATZINAKOS D. Digital watermarking for tell-tale tamper proof i ng and authentication [J]. Proc IEEE, 1999,87(7): 1167−1180.

[7] HUANG Y, LU W, SUN W, LONG D. Improved DCT-based detection of copy-move forgery in images [J]. Forensic Science International, 2011, 206(1/3): 178−184.

[8] STAMM M C, LIU K J R. Forensic detection of image manipulation using statistical intrinsic fingerprints [J]. IEEE Transactions on Information Forensics and Security, 2010, 5(3):492−506.

[9] PENG F, NIE Y, LONG M. A complete passive blind image copy-move forensics scheme based on compound statistic features [J]. Forensic Science International, 2011, 212(1/3):21−25.

[10] MAHDIAN B, SAIC S. Using noise inconsistencies for blind image forensics [J]. Image and Vision Computing, 2009, 27(10):1497−1503.

[11] LI Wei-hai, YUAN Yuan, YU Neng-hai. Passive detection of doctored JPEG image via block artifact grid extraction [J].Signal Processing, 2009, 89(9): 1821−1829.

[12] CRIMINISI A, PEREZ P, TOYAMA K. Region filling and object removal by exemplar-based inpainting [J]. IEEE Transactions on Image Processing, 2004, 13(9): 1200−1212.

[13] BERTALMIO M, VESE L, SAPIOR G, OSHER S. Simultaneous structure and texture image inpainting [J]. IEEE Transactions on Image Processing, 2003, 12(8): 882−889.

[14] LIN Z, HE J, TANG X, TANG C K. Fast, automatic and fine-grained tampered JPEG image detection via DCT coefficient analysis [J]. Pattern Recognition, 2009, 42(11):2492−2501.

[15] LIU Z, LI X, ZHAO Y. Passive detection of copy-paste tampering for digital image forensics [C]// Proc Fourth International Conference on Intelligent Computation Technology and Automation. Los Alamitos, Canada: IEEE Computer Society,2011: 649−652.