CN106851140B - 一种使用空域平滑滤波的数字照片图像来源识别方法 - Google Patents

Abstract

本发明公开了一种使用空域平滑滤波的数字照片图像来源识别方法，其基于相机的传感器模式噪声和数字照片图像的传感器模式噪声，并考虑JPEG压缩对相机的传感器模式噪声和数字照片图像的传感器模式噪声的质量的影响，充分利用相机的传感器模式噪声和数字照片图像的传感器模式噪声的本质特性，采用高斯白噪声图像对零均值化后的相机的传感器模式噪声和数字照片图像的传感器模式噪声分别进行空域平滑滤波处理，能够有效抑制JPEG压缩噪声；由于本发明方法能有效剔除相机的传感器模式噪声和数字照片图像的传感器模式噪声中的干扰成分，能有效改善所获得的相机的传感器模式噪声和数字照片图像的传感器模式噪声的质量，因此大大提高了本发明方法的识别率。

Description

一种使用空域平滑滤波的数字照片图像来源识别方法

技术领域

本发明涉及一种图像取证技术，尤其是涉及一种使用空域平滑滤波的数字照片图像来源识别方法。

背景技术

随着数码相机、数码摄像机，尤其是带拍照功能的智能手机等各种摄影产品的普及，数字照片图像在人们的日常生活中得到了广泛使用。作为一种重要的信息载体，数字照片图像在新闻报道、科学研究、保险申诉调查、犯罪调查以及情报分析等领域发挥着重要的作用。然而，在这些特殊场合，如何保证作为证据的数字照片图像来源的可靠，一直是一个具有挑战性的研究课题。

对数字照片图像来源进行鉴别时，最直观和最简单的方法就是查看数字照片图像的EXIF(Exchange Image File，可交换图像文件)，但是随着一些功能强大的图像编辑软件的出现，更改EXIF信息已变得非常容易。一些生产商采用在数码相机中嵌入不可见的脆弱水印的方法来保护数字照片图像的真实性和完整性，然而此种做法不仅会影响数字照片图像的质量，而且成本过高，目前大多数品牌的数码相机都没有采用该技术。随着数字照片图像取证技术研究的深入，研究者们把注意力从主动取证转移到从数字照片图像中寻找数码相机内在特性的被动取证上。一般而言，每幅数字照片图像的形成都要经历感光、成像和图像信号处理等一系列硬件和软件操作，虽然这些操作会不可避免地对数字照片图像中的像素点产生负面影响，但是同时也为数字照片图像的来源识别提供了思路。如：MehdiKharrazi,Nasir Memon.Blind Source Camera Identification[C].IEEE InternationalConference on Image Processing 2004:709-712.(梅迪·哈拉奇，纳西尔·麦蒙，盲源相机识别[C]，IEEE国际图像处理会议，2004:709-712)，提出了利用相机镜头径向失真对数字照片图像统计特性的影响进行数字照片图像来源认证。又如：Sevinc Bayram,NasirMemon.Source camera identification based on CFA interpolation[C].IEEEInternational Conference on Image Processing,2005,3:III-69-72.(塞维克·拜拉姆，纳西尔·麦蒙，基于CFA插值特性的源相机识别技术[C]，IEEE国际图像处理会议，2005,3:III-69-72)，提出了利用不同CFA(Color Filter Array，色彩滤波阵列)插值算法引起的差异进行识别。还如：Hany Farid.Blind Inverse Gamma Correction[J].IEEETransactions on Image Proceessing,2001,10:1428-1433(汉尼·福瑞德，盲逆伽玛校正[J],IEEE图像处理，2001，第10期：1428-1433)认为伽马校正过程会在数字照片图像的频域引入高阶相关性，利用其估计每部相机的校正曲线，从而实现来源检测目的。Yu-Feng Hsu,Shi-Fu Chang.Detecting image splicing using geometry invariants and cameracharacteristics consistency[C].International Conference on Multimedia andExpo,2006:549-652.(徐堉峰，张士辅，使用几何不变量和相机特征一致性检测图像拼接[C]，国际多媒体会议，2006:549-652)，提出了利用相机响应函数识别来源的方法。考虑到不同的相机可能采用不同JPEG量化表，Matthew James Sorrell,Digital camera sourceidentification through JPEG quantization[J],Information Science Reference,2008:291-313.(马修·詹姆斯·索罗，基于JPEG压缩的数字相机识别[J]，信息科学参考，2008:291-313)采用量化差异鉴别不同来源的数字照片图像。

良好的数字照片图像来源识别算法应该具有高检测率和强鲁棒性，上述方法虽然在一定程度上实现了识别目的，但是普遍存在计算复杂度高、检测率低的问题。尤其是利用不同CFA插值算法引起的差异进行识别的方法和采用量化差异鉴别不同来源的方法，若遇到采用相同插值方式或相同JPEG量化表的相机时，对应的方法将会失效。经过深入研究，Jan Lukas,Jessica Fridrich,and Miroslav Goljan,Digital camera identificationfrom sensor pattern noise,IEEE Transaction on Information Forensics Security,2006,1(2):205–214.(简·卢卡斯，杰西卡·弗里德里希，米罗斯拉夫·格兰，基于模式噪声的源相机识别[J]，IEEE信息与安全，2006，第1期(2):205-214)，提出了一种非常有效的数字照片图像来源识别方法——基于相机的传感器模式噪声(Sensor Pattern Noise，SPN)的数字照片图像来源检测算法。

相机的传感器模式噪声主要是由传感器的制作工艺不完美和材质不一致所导致的。由于每部相机的传感器模式噪声都具有唯一性，因此可看作相机的指纹。一种经典的基于模式噪声的照片来源识别方法主要包括以下步骤：1)通过对一部相机拍摄的多幅图像的噪声残差进行最大似然估计获取相机的传感器模式噪声R；2)获取待测图像的噪声残差W，作为待测图像的模式噪声；3)用事先设置的阈值对待测图像的模式噪声W与相机的传感器模式噪声R的相关性ρ(W,R)进行判决，实现数字照片图像的来源识别。然而，在实际中由于相机的传感器模式噪声是弱信号，容易受到数字照片图像的内容、CFA插值和JPEG压缩等不利因素干扰，致使识别结果不理想。为了消除不利因素影响，改善相机的传感器模式噪声的质量，提高来源识别率，研究者提出了多种不同的方法。如：Mo Chen、Jessica Fridrich、Miroslav Goljan.Determining Image Origin and Integrity Using Sensor Noise[J].IEEE Transactions on Information and Security,2008,3(1):74-90.(陈莫、杰西卡·弗里德里希、米罗斯拉夫·格兰，基于传感器模式噪声的图像来源与完整性检测[J].IEEE信息与安全.2008，第3期(1):74-90)，其在获取相机的传感器模式噪声时，提出使用内容平滑(如蓝天)的数字照片图像，并在降噪时使用小波变换方法，消除数字照片图像的内容影响。在此基础上，Giovanni Chierchia,Sara Parrilli,Giovanni Poggi,et al.On theinfluence of denoising in PRNU based forgery detection[C].In Processingsecond ACM Workshop Multimedia Forensics,2010:117-122.(乔凡尼·基耶尔基亚，萨娃·帕里利，乔凡尼·波吉，不同降噪方法对基于模式噪声的图像篡改取证的影响[C]，第二届ACM多媒体信息取证会议，2010:117-122)和Kang Xian-gui,Chen Jie.A context-adaptive SPN predictor for trustworthy source camera identification[J],EURASIPJournal on Image and Video Processing,2014,2014(1):1–11.(康显桂，陈杰.上下文自适应传感器模式噪声的预测值得信赖的相机源辨识[J].欧洲信号处理协会之图像与视频处理杂志，2014，第2014期(1)：1-11)，在降噪时分别提出采用BM3D(BlockMatching 3D，三维块匹配)和PCAI(Predictor based on Context AdaptiveInterpolation，基于上下文自适应插值的预测)等降噪方法，也取得了良好的结果。对于CFA插值和JPEG压缩噪声，Chen M,Fridrich J,Goljan M,et al.Determining ImageOrigin and Integrity Using Sensor Noise[J].Information Forensics&SecurityIEEE Transactions on,2008,3(1):74-90.(陈莫、杰西卡·弗里德里希、米罗斯拉夫·格兰，基于传感器模式噪声的图像来源与完整性检测[J].IEEE信息取证与安全.2008，第3期(1):74-90)，提出了零均值化法(zero-mean，ZM)和维纳滤波(Wiener filter，WF)的方法。Lin X,Li C.-T.Preprocessing Reference Sensor Pattern Noise via SpectrumEqualization[J].IEEE Transactions on Information Forensics and Security,2016,11(1):126-140.(林旭、李长纯，基于频谱均衡方法处理参考相机传感器模式噪声[J]，IEEE信息取证与安全，2016，第11期(1):126-140)，提出了只利用绿色通道数据获取相机的传感器模式噪声，以此消除CFA插值的影响，然后采用频谱均衡化法(Spectrum EqualizationAlgorithm，SEA)抑制JPEG压缩噪声。Hu Yong-jian,Yu Bing-hua,Jian Chao,Sourcecamera identification using large components of sensor pattern noise[C].inprocessing second international Conference Compute,2009:1-5.(胡永健，俞兵华，简超，使用模式噪声大分量的源相机识别[C]，第二届国际计算会议，2009:1-5)，其认为相机的传感器模式噪声中的大分量更能代表相机的传感器模式噪声，为消除干扰因素，建议在识别数字照片图像设备来源时将相机的传感器模式噪声中的像素点按像素值的大小进行排序，只利用像素值最大的一部分像素点(如前10％)。考虑到相机的传感器模式噪声是弱信号，Li C.-T.Source camera identification using enhanced sensor pattern noise[J].IEEE Transactions on Information Forensics and Security,2010,5(2):280-287.(李长纯，基于增强模式噪声的源相机识别[J]，IEEE信息取证与安全，2010，第5期(2):280-287)，提出六种模型用于增强数字照片图像的模式噪声。上述这些方法虽然都能改善相机的传感器模式噪声的质量，但是未能完全消除数字照片图像生成过程中的JPEG压缩噪声，致使识别效果提高有限，且对抵抗JPEG二次压缩的鲁棒性较差。

基于以上原因，为了消除JPEG压缩噪声，改善相机的传感器模式噪声的质量，提高来源识别率，有必要研究一种考虑JPEG压缩影响的数字照片图像来源识别方法。

发明内容

本发明所要解决的技术问题是提供一种使用空域平滑滤波的数字照片图像来源识别方法，其能够同时有效消除CFA插值和JPEG压缩噪声的干扰，能够有效改善相机的传感器模式噪声的质量，提高来源识别率，且对抵抗较低的JPEG压缩具有较强的鲁棒性。

本发明解决上述技术问题所采用的技术方案为：一种使用空域平滑滤波的数字照片图像来源识别方法，其特征在于包括以下步骤：

①获取N幅原始数字照片图像，将获取的第n幅原始数字照片图像记为I_org,n；并任意选取一幅数字照片图像作为待识别照片图像，记为I_test；其中，N≥2，1≤n≤N，获取每幅原始数字照片图像所采用的相机与I_test所声明的来源相机为同一部相机，每幅原始数字照片图像和I_test的尺寸大小相同，且宽度为col、高度为row；

②获取每幅原始数字照片图像的降噪图像，将I_org,n的降噪图像记为I'_org,n；然后根据每幅原始数字照片图像及其降噪图像，获取每幅原始数字照片图像的噪声残差图像，将I_org,n的噪声残差图像记为W_org,n，将W_org,n中坐标位置为(i,j)的像素点的像素值记为W_org,n(i,j)，W_org,n(i,j)＝I_org,n(i,j)-I'_org,n(i,j)；接着根据所有原始数字照片图像的噪声残差图像，获取相机的传感器模式噪声，记为R；再消除R中的CFA插值噪声，得到相机的无CFA插值噪声的传感器模式噪声，记为R'；其中，1≤i≤col,1≤j≤row，I_org,n(i,j)表示I_org,n中坐标位置为(i,j)的像素点的像素值，I'_org,n(i,j)表示I'_org,n中坐标位置为(i,j)的像素点的像素值；

同样，获取I_test的降噪图像，记为I'_test；然后根据I_test及I'_test，获取I_test的噪声残差图像，记为W_test，将W_test中坐标位置为(i,j)的像素点的像素值记为W_test(i,j)，W_test(i,j)＝I_test(i,j)-I'_test(i,j)；接着将W_test作为I_test的传感器模式噪声；再消除W_test中的CFA插值噪声，得到I_test的无CFA插值噪声的传感器模式噪声，记为W'_test；其中，I_test(i,j)表示I_test中坐标位置为(i,j)的像素点的像素值，I'_test(i,j)表示I'_test中坐标位置为(i,j)的像素点的像素值；

③随机生成一幅宽度为col且高度为row的标准高斯白噪声图像，记为G；然后根据G在空域对R'进行空域平滑滤波处理，获得相机的最终的传感器模式噪声，记为R^*；同样，根据G在空域对W'_test进行空域平滑滤除处理，获得I_test的最终的传感器模式噪声，记为W^* _test；

④采用SPCE方法计算W^* _test与R^*的相关性，记为SPCE(W^* _test,R^*)；

⑤用事先设置的I_test所声明的来源相机下的相关性判别阈值对SPCE(W^* _test,R^*)进行判决，若SPCE(W^* _test,R^*)大于或等于相关性判别阈值，则判定I_test来自其所声明的来源相机；若SPCE(W^* _test,R^*)小于相关性判别阈值，则判定I_test来自非其所声明的来源相机，即I_test来自其它相机。

所述的步骤②中获取每幅原始数字照片图像的降噪图像和获取I_test的降噪图像所采用的方法为小波降噪处理方法；所述的步骤②中采用最大似然估计法对所有原始数字照片图像的噪声残差图像进行处理，获得R；所述的步骤②中采用零均值化法对R进行处理以消除R中的CFA插值噪声，获得R'；所述的步骤②中采用零均值化法对W_test进行处理以消除W_test中的CFA插值噪声，获得W'_test。

所述的步骤③中根据G在空域对R'进行空域平滑滤波处理获得R^*及根据G在空域对W'_test进行空域平滑滤除处理获得W^* _test的具体过程为：

③_1、对G的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

同样，对R'的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

同样，对W'_test的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

③_2、计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varG，varG的宽度为col且高度为row；

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varR，varR的宽度为col且高度为row；

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varW，varW的宽度为col且高度为row；

然后计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGR，covGR的宽度为col且高度为row；同样，计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGW，covGW的宽度为col且高度为row；

③_3、计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^R(i,j)和b^R(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^R，a^R的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^R，b^R的宽度为col且高度为row；其中，1≤i≤col,1≤j≤row，covGR(i,j)表示covGR中坐标位置为(i,j)的像素点的像素值，varG(i,j)表示varG中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，ξ为惩罚因子；

同样，计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^W(i,j)和b^W(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^W，a^W的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^W，b^W的宽度为col且高度为row；其中，covGW(i,j)表示covGW中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值；

③_4、对a^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；

同样，对b^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；

同样，对a^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；

同样，对b^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；

③_5、获取R^*，将R^*中坐标位置为(i,j)的像素点的像素值记为R^*(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，G(i,j)表示G中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值；

同样，获取W^* _test，将W^* _test中坐标位置为(i,j)的像素点的像素值记为W^* _test(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值。

所述的步骤③_3中取ξ＝0.75。

所述的步骤④中其中，sign()为求符号函数，P_max表示W^* _test与R^*的相关函数图像P中最大的像素值，Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域的尺寸大小，8≤Ne≤15，Ω表示W^* _test与R^*的相关函数图像P中的所有像素点的坐标位置的集合，Ω_Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域中的所有像素点的坐标位置的集合，P(s)表示W^* _test与R^*的相关函数图像P中坐标位置为s的像素点的像素值。

所述的步骤⑤中相关性判别阈值的获取过程为：

⑤_1、使用I_test所声明的来源相机任意拍摄U幅同源图像；然后按照步骤②和步骤③的过程，以相同的方式获得每幅同源图像的最终的传感器模式噪声；其中，U≥100，每幅同源图像的宽度为col且高度为row；

⑤_2、选取与I_test所声明的来源相机不同的V部相机，要求选取的V部相机的分辨率大于或等于I_test所声明的来源相机的分辨率；然后使用选取的V部相机任意拍摄U幅异源图像，共得到U幅异源图像，再将每幅异源图像按照步骤⑤_1获得的同源图像的尺寸大小截取中心区域作为新的异源图像；接着按照步骤②和步骤③的过程，以相同的方式获得每幅新的异源图像的最终的传感器模式噪声；其中，1≤V≤9；

⑤_3、采用SPCE方法计算R^*与每幅同源图像的最终的传感器模式噪声的相关性，共得到U个同源样本相关性值；

同样，采用SPCE方法计算R^*与每幅异源图像的最终的传感器模式噪声的相关性，共得到U个异源样本相关性值；

⑤_4、在区间[0,30]内选取不同的实数值对步骤⑤_3中获得的U个同源相关性值和U个异源相关性值进行判决，获得每个实数值下的假正率(False Positive Rate，FPR)，记为FPR，其中，N_FP表示将异源图像判决为同源的数量，N_TN表示将异源图像判决为异源的数量；

⑤_5、将FPR＝0.01时对应的实数值作为相关性判决阈值。

与现有技术相比，本发明的优点在于：

1)本发明方法基于相机的传感器模式噪声和数字照片图像的传感器模式噪声，并考虑JPEG压缩对相机的传感器模式噪声和数字照片图像的传感器模式噪声的质量的影响，充分利用相机的传感器模式噪声和数字照片图像的传感器模式噪声的本质特性，采用高斯白噪声图像对零均值化后的相机的传感器模式噪声和数字照片图像的传感器模式噪声分别进行空域平滑滤波处理，能够有效抑制JPEG压缩噪声。

2)本发明方法与现有的模式噪声处理方法相比，由于本发明方法能有效剔除相机的传感器模式噪声和数字照片图像的传感器模式噪声中的干扰成分，能有效改善所获得的相机的传感器模式噪声和数字照片图像的传感器模式噪声的质量，因此大大提高了本发明方法的识别率。

3)本发明方法能有效解决JPEG压缩对相机的传感器模式噪声和数字照片图像的传感器模式噪声的负面影响，对较低的JPEG压缩质量带来的破坏也能恢复相机的传感器模式噪声和数字照片图像的传感器模式噪声的原有特性，因此大大提高了本发明方法对抵抗JPEG压缩的鲁棒性。

附图说明

图1为本发明方法的总体实现框图。

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

本发明提出的一种使用空域平滑滤波的数字照片图像来源识别方法，其总体实现框图如图1所示，其包括以下步骤：

①获取N幅原始数字照片图像，将获取的第n幅原始数字照片图像记为I_org,n；并任意选取一幅数字照片图像作为待识别照片图像，记为I_test；其中，N≥2，在本实施例中取N＝60，1≤n≤N，获取每幅原始数字照片图像所采用的相机与I_test所声明的来源相机为同一部相机，每幅原始数字照片图像和I_test的尺寸大小相同，且宽度为col、高度为row。

在本实施例中，在拍摄原始数字照片图像时，通常选择高亮不饱和的蓝天作为获取场景。

②获取每幅原始数字照片图像的降噪图像，将I_org,n的降噪图像记为I'_org,n；然后根据每幅原始数字照片图像及其降噪图像，获取每幅原始数字照片图像的噪声残差图像，将I_org,n的噪声残差图像记为W_org,n，将W_org,n中坐标位置为(i,j)的像素点的像素值记为W_org,n(i,j)，W_org,n(i,j)＝I_org,n(i,j)-I'_org,n(i,j)；接着根据所有原始数字照片图像的噪声残差图像，获取相机的传感器模式噪声，记为R；再消除R中的CFA插值噪声，得到相机的无CFA插值噪声的传感器模式噪声，记为R'；其中，1≤i≤col,1≤j≤row，I_org,n(i,j)表示I_org,n中坐标位置为(i,j)的像素点的像素值，I'_org,n(i,j)表示I'_org,n中坐标位置为(i,j)的像素点的像素值。

同样，获取I_test的降噪图像，记为I'_test；然后根据I_test及I'_test，获取I_test的噪声残差图像，记为W_test，将W_test中坐标位置为(i,j)的像素点的像素值记为W_test(i,j)，W_test(i,j)＝I_test(i,j)-I'_test(i,j)；接着将W_test作为I_test的传感器模式噪声；再消除W_test中的CFA插值噪声，得到I_test的无CFA插值噪声的传感器模式噪声，记为W'_test；其中，I_test(i,j)表示I_test中坐标位置为(i,j)的像素点的像素值，I'_test(i,j)表示I'_test中坐标位置为(i,j)的像素点的像素值。

在此具体实施例中，步骤②中获取每幅原始数字照片图像的降噪图像和获取I_test的降噪图像所采用的方法为小波降噪处理方法；步骤②中采用现有的最大似然估计法对所有原始数字照片图像的噪声残差图像进行处理，获得R；步骤②中采用现有的零均值化法对R进行处理以消除R中的CFA插值噪声，获得R'，采用零均值化法是为了尽量去除R中的CFA插值噪声；步骤②中采用现有的零均值化法对W_test进行处理以消除W_test中的CFA插值噪声，获得W'_test，采用零均值化法是为了尽量去除W_test中的CFA插值噪声。

③随机生成一幅宽度为col且高度为row的标准高斯白噪声图像，记为G；然后根据G在空域对R'进行空域平滑滤波处理，获得相机的最终的传感器模式噪声，记为R^*；同样，根据G在空域对W'_test进行空域平滑滤除处理，获得I_test的最终的传感器模式噪声，记为W^* _test。

在此具体实施例中，步骤③中根据G在空域对R'进行空域平滑滤波处理获得R^*及根据G在空域对W'_test进行空域平滑滤除处理获得W^* _test的具体过程为：

③_1、对G的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对G的上边界和下边界各扩充一行像素点且像素值设定为0，对G的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块。

同样，对R'的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对R'的上边界和下边界各扩充一行像素点且像素值设定为0，对R'的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块。

同样，对W'_test的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对W'_test的上边界和下边界各扩充一行像素点且像素值设定为0，对W'_test的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块。

③_2、计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varG，varG的宽度为col且高度为row。

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varR，varR的宽度为col且高度为row。

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varW，varW的宽度为col且高度为row。

然后计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGR，covGR的宽度为col且高度为row；同样，计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGW，covGW的宽度为col且高度为row。

③_3、计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^R(i,j)和b^R(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^R，a^R的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^R，b^R的宽度为col且高度为row；其中，1≤i≤col,1≤j≤row，covGR(i,j)表示covGR中坐标位置为(i,j)的像素点的像素值，varG(i,j)表示varG中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，ξ为惩罚因子，在本实施例中取ξ＝0.75。

同样，计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^W(i,j)和b^W(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^W，a^W的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^W，b^W的宽度为col且高度为row；其中，covGW(i,j)表示covGW中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值。

③_4、对a^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对a^R的上边界和下边界各扩充一行像素点且像素值设定为0，对a^R的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row。

同样，对b^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对b^R的上边界和下边界各扩充一行像素点且像素值设定为0，对b^R的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row。

同样，对a^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对a^W的上边界和下边界各扩充一行像素点且像素值设定为0，对a^W的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row。

同样，对b^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为即对b^W的上边界和下边界各扩充一行像素点且像素值设定为0，对b^W的左边界和右边界各扩充一列像素点且像素值设定为0；然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为的宽度为col且高度为row。

③_5、获取R^*，将R^*中坐标位置为(i,j)的像素点的像素值记为R^*(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，G(i,j)表示G中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值。

同样，获取W^* _test，将W^* _test中坐标位置为(i,j)的像素点的像素值记为W^* _test(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值。

④采用SPCE(sign peak-to-correlation energy，有符号的峰值相关能量比)方法计算W^* _test与R^*的相关性，记为SPCE(W^* _test,R^*)。

在此具体实施例中，步骤④中其中，sign()为求符号函数，P_max表示W^* _test与R^*的相关函数图像P中最大的像素值，Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域的尺寸大小，8≤Ne≤15，在本实施例中取Ne＝11，Ω表示W^* _test与R^*的相关函数图像P中的所有像素点的坐标位置的集合，Ω_Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域中的所有像素点的坐标位置的集合，P(s)表示W^* _test与R^*的相关函数图像P中坐标位置为s的像素点的像素值。

⑤用事先设置的I_test所声明的来源相机下的相关性判别阈值对SPCE(W^* _test,R^*)进行判决，若SPCE(W^* _test,R^*)大于或等于相关性判别阈值，则判定I_test来自其所声明的来源相机；若SPCE(W^* _test,R^*)小于相关性判别阈值，则判定I_test来自非其所声明的来源相机，即I_test来自其它相机。

在此具体实施例中，关于步骤⑤中事先设置的相关性判别阈值，由于每部相机的相关性阈值可能不同，因此获得的阈值也不相同，阈值获取的具体过程如下：

⑤_1、使用I_test所声明的来源相机任意拍摄U幅同源图像；然后按照步骤②和步骤③的过程，以相同的方式获得每幅同源图像的最终的传感器模式噪声；其中，U≥100，在本实施例中取U＝300，每幅同源图像的宽度为col且高度为row。

⑤_2、选取与I_test所声明的来源相机不同的V部相机，要求选取的V部相机的分辨率大于或等于I_test所声明的来源相机的分辨率；然后使用选取的V部相机任意拍摄U幅异源图像，共得到U幅异源图像，再将每幅异源图像按照步骤⑤_1获得的同源图像的尺寸大小截取中心区域作为新的异源图像；接着按照步骤②和步骤③的过程，以相同的方式获得每幅新的异源图像的最终的传感器模式噪声；其中，1≤V≤9，在本实施例中取V＝9。

⑤_3、采用SPCE方法计算R^*与每幅同源图像的最终的传感器模式噪声的相关性，共得到U个同源样本相关性值。

同样，采用SPCE方法计算R^*与每幅异源图像的最终的传感器模式噪声的相关性，共得到U个异源样本相关性值。

⑤_4、在区间[0,30]内选取不同的实数值对步骤⑤_3中获得的U个同源相关性值和U个异源相关性值进行判决，获得每个实数值下的假正率(False Positive Rate，FPR)，记为FPR，其中，N_FP表示将异源图像判决为同源的数量，N_TN表示将异源图像判决为异源的数量。

⑤_5、将FPR＝0.01时对应的实数值作为相关性判决阈值。

表1给出了十部不同的相机各自拍摄的数字照片图像的数量。

表1 十部智能手机的相机型号、分辨率、数字照片图像的数量及数字照片图像的格式

按照步骤⑤_1至步骤⑤_5对表1中HTC one M7手机的相机进行处理，获得该相机下的相关性判别阈值为11.4031；同样，获得华为荣耀6手机的相机下的相关性判别阈值为11.0204、华为mate 7手机的相机下的相关性判别阈值为8.1747、苹果5s_1手机的相机下的相关性判别阈值为12.1165、苹果5s_2手机的相机下的相关性判别阈值为10.6006、苹果6手机的相机下的相关性判别阈值为10.3048、OPPO R7手机的相机下的相关性判别阈值为7.9123、一加A0001手机的相机下的相关性判别阈值为10.9117、索尼Z1手机的相机下的相关性判别阈值为8.9025、索尼Z3手机的相机下的相关性判别阈值为8.5082。

为了进一步说明本发明方法的可行性和有效性，对本发明方法进行实验验证。

具有高清拍照功能的智能手机在日常生活中得到了广泛普及，但由于存储问题，智能手机的相机多采用较低的JPEG压缩质量，因此对相机的传感器模式噪声破坏更加严重。考虑到本发明方法主要是为了消除JPEG压缩的影响，因此实验时仍选用表1给出的十部智能手机的相机拍摄的数字照片图像建立一个相机照片图像库，如表2所列。在建立相机照片图像库时，为了有别于表1，针对每部智能手机的相机选取的数字照片图像与表1中的数字照片图像互不相同，另外选取了450幅数字照片图像。由于相机的分辨率不尽相同，为便于实验，当遇到维度不一的情况，按照小尺寸的数字照片图像的维度在大尺寸的数字照片图像的中心截取相同维度参与实验。测试环境为Windows 7平台，Intel Core i5-3470处理器，CPU频率为3.2GHz，编译环境为Matlab R2012a。

表2 十部智能手机的相机型号、分辨率、数字照片图像的数量及数字照片图像的格式

为了能客观公正的评价发明方法的性能，采用阈值判决和Kappa统计两种方法进行分析。

方法一：阈值判决方法

将表2中每部相机的传感器模式噪声分别与450幅同源照片的传感器模式噪声和450幅异源照片(对九部异源相机拍摄的4500幅数字照片图像，随机挑选出450幅)的传感器模式噪声进行相关性计算。以第一部智能手机的相机HTC One M7为例，将其传感器模式噪声分别与450幅来自该相机的数字照片图像的传感器模式噪声计算相关性，得到450个同源样本的相关性值；再将其与来自剩余九部相机的450幅数字照片图像的传感器模式噪声计算相关性，得到450个异源样本的相关性值；然后对这900个相关性值进行阈值决策，可得到检测结果的真正率(True Positive Rate，TPR)，其中，N_TP表示将同源样本判决为同源的数量，N_FN表示将同源样本判决为异源的数量。

按照此方法，依次完成其余九部相机的传感器模式噪声与相机照片图像库中的450幅同源照片和450幅异源照片的相关性计算，再根据阈值决策得到各个相机下的检测结果的真正率，然后对十部相机下的检测结果的真正率求均平均值，获得平均真正率。

当需要分析本发明方法对来源识别的影响，可通过比较平均真正率进行评估。

方法二：Kappa统计方法

采用Kappa统计方法时，将表2中每一幅数字照片图像的最终的传感器模式噪声分别与十部相机各自的最终的传感器模式噪声计算相关性值，比较得到的十个相关性值，找出最大相关性值对应的相机，就认为该数字照片图像来自此相机。将表2中的4500幅数字照片图像全部计算完毕，就可得到混淆矩阵，以此获得Kappa系数：其中，K为总样本的数目，在此K＝4500，O表示实际正确预测的样本数目，e为理论正确预测的样本数目。

采用方法[1](Jan Lukas,Jessica Fridrich,and Miroslav Goljan,Digitalcameraidentification from sensor pattern noise,IEEE Transaction onInformation Forensics Security,2006,1(2):205–214.(简·卢卡斯，杰西卡·弗里德里希，米罗斯拉夫·格兰，基于模式噪声的源相机识别[J]，IEEE信息与安全，2006，第1期(2):205-214))、方法[2](Mo Chen、Jessica Fridrich、Miroslav Goljan.Determining ImageOrigin and Integrity Using Sensor Noise[J].IEEE Transactions on Informationand Security,2008,3(1):74-90.(陈莫、杰西卡·弗里德里希、米罗斯拉夫·格兰，基于传感器模式噪声的图像来源与完整性检测[J].IEEE信息与安全.2008，第3期(1):74-90))、方法[3](Lin X,Li C.-T.Preprocessing Reference Sensor Pattern Noise via SpectrumEqualization[J].IEEE Transactions on Information Forensics and Security,2016,11(1):126-140.(林旭、李长纯，基于频谱均衡方法处理参考相机传感器模式噪声[J]，IEEE信息取证与安全，2016，第11期(1):126-140))及本发明方法在表2所给的数字照片图像库上进行实验。

来源检测结果：

对实验方法一得到的每部相机的450个同源相关性值和450个异源相关性值，根据各自的阈值进行判决，得到每部相机的检测结果的真正率，然后再平均得到对所有相机的检测结果的平均真正率。表3列出了方法[1]、方法[2]、方法[3]、本发明方法对表2中的所有相机的检测结果的平均真正率。观察表3所列的结果，四种方法对所有相机的检测结果的平均真正率分别为0.9336、0.9580,0.9513和0.9820，比较可知，本发明方法对所有相机的检测结果的平均真正率明显高于其他三种方法，由此可证明本发明方法的性能优于其他方法。

表3 方法[1]、方法[2]、方法[3]和本发明方法对表2中所有相机的检测结果的平均真正率

方法	方法[1]	方法[2]	方法[3]	本发明方法
					平均真正率	0.9336	0.9580	0.9513	0.9820

为进一步说明本发明方法的优势，采用方法[1]、方法[2]、方法[3]、本发明方法对表2中的4500幅数字照片图像按实验方法二进行实验，获得混淆矩阵和Kappa系数。表4给出了方法[1]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数，表5给出了方法[2]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数，表6给出了方法[3]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数，表7给出了本发明方法对表3中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数。表4、表5、表6和表7中的混淆矩阵中的每一行代表一部相机拍摄的所有数字照片图像的检测结果。以检测结果最差的索尼Z3智能手机的相机为例，针对方法[1]，只有244幅数字照片图像得到正确的检测，剩余206幅数字照片图像被错认为来自其余九部相机；而方法[2]和方法[3]的检测结果也非常不理想，相比之下，本发明方法对450幅数字照片图像的检测精度达到了86.22％。为了衡量不同方法对十部相机拍摄的数字照片图像的检测精度，对表4到表7中的混淆矩阵求解Kappa系数，分别为0.9472、0.9570、0.9536和0.9830。对比表4到表7可知，本发明方法的检测精度最高，进一步证明了本发明方法具有良好的性能。

表4 方法[1]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数

表5 方法[2]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数

表6 方法[3]对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数

表7 本发明方法对表2中的所有数字照片图像的检测结果的混淆矩阵和Kappa系数

鲁棒性分析：

JPEG压缩会降低相机的最终的传感器模式噪声的质量，为评估本发明方法对JPEG压缩的鲁棒性，对表2中所列的所有数字照片图像分别采用90％、70％和50％等三种JPEG压缩质量进行压缩，然后采用不同方法进行检测。

表8列出了方法[1]、方法[2]、方法[3]、本发明方法对表2中的所有数字照片图像分别采用90％、70％、50％的JPEG压缩质量压缩后的检测结果的平均真正率，表9列出了方法[1]、方法[2]、方法[3]、本发明方法对表2中的所有数字照片图像分别采用90％、70％、50％的JPEG压缩质量压缩后的检测结果的Kappa系数。

观察表8所列的平均真正率和表9所列的Kappa系数可知，在相同JPEG压缩质量下，方法[1]、方法[2]和方法[3]的性能明显低于本发明方法；而随着JPEG压缩质量的降低，方法[1]的性能下降最明显，方法[2]和方法[3]也有所下降，而本发明方法表现良好。统计表8中所列的四种方法在三种JPEG压缩质量下最小的平均真正率，将其与表3中未压缩时的平均真正率相比较，分别下降了0.0132、0.0167、0.0095和0.0036，由此可说明本发明方法对抵抗JPEG压缩具有较强的鲁棒性，观察表9所列的Kappa系数也证实了这一点。

表8 方法[1]、方法[2]、方法[3]、本发明方法对表2中的所有数字照片图像分别采用90％、70％、50％的JPEG压缩质量压缩后的检测结果的平均真正率

表9 方法[1]、方法[2]、方法[3]、本发明方法对表2中的所有数字照片图像分别采用90％、70％、50％的JPEG压缩质量压缩后的检测结果的Kappa系数

大量实验结果表明，与现有方法相比，本发明方法能有效改善相机的传感器模式噪声的质量，提高数字照片图像来源的识别率，并对抵抗较低的JPEG压缩具有较强的鲁棒性。

Claims (4)

1.一种使用空域平滑滤波的数字照片图像来源识别方法，其特征在于包括以下步骤：

①获取N幅原始数字照片图像，将获取的第n幅原始数字照片图像记为I_org,n；并任意选取一幅数字照片图像作为待识别照片图像，记为I_test；其中，N≥2，1≤n≤N，获取每幅原始数字照片图像所采用的相机与I_test所声明的来源相机为同一部相机，每幅原始数字照片图像和I_test的尺寸大小相同，且宽度为col、高度为row；

②获取每幅原始数字照片图像的降噪图像，将I_org,n的降噪图像记为I'_org,n；然后根据每幅原始数字照片图像及其降噪图像，获取每幅原始数字照片图像的噪声残差图像，将I_org,n的噪声残差图像记为W_org,n，将W_org,n中坐标位置为(i,j)的像素点的像素值记为W_org,n(i,j)，W_org,n(i,j)＝I_org,n(i,j)-I'_org,n(i,j)；接着根据所有原始数字照片图像的噪声残差图像，获取相机的传感器模式噪声，记为R；再消除R中的CFA插值噪声，得到相机的无CFA插值噪声的传感器模式噪声，记为R'；其中，1≤i≤col,1≤j≤row，I_org,n(i,j)表示I_org,n中坐标位置为(i,j)的像素点的像素值，I'_org,n(i,j)表示I'_org,n中坐标位置为(i,j)的像素点的像素值；

同样，获取I_test的降噪图像，记为I'_test；然后根据I_test及I'_test，获取I_test的噪声残差图像，记为W_test，将W_test中坐标位置为(i,j)的像素点的像素值记为W_test(i,j)，W_test(i,j)＝I_test(i,j)-I'_test(i,j)；接着将W_test作为I_test的传感器模式噪声；再消除W_test中的CFA插值噪声，得到I_test的无CFA插值噪声的传感器模式噪声，记为W'_test；其中，I_test(i,j)表示I_test中坐标位置为(i,j)的像素点的像素值，I'_test(i,j)表示I'_test中坐标位置为(i,j)的像素点的像素值；

③随机生成一幅宽度为col且高度为row的标准高斯白噪声图像，记为G；然后根据G在空域对R'进行空域平滑滤波处理，获得相机的最终的传感器模式噪声，记为R^*；同样，根据G在空域对W'_test进行空域平滑滤除处理，获得I_test的最终的传感器模式噪声，记为W^* _test；

所述的步骤③中根据G在空域对R'进行空域平滑滤波处理获得R^*及根据G在空域对W'_test进行空域平滑滤除处理获得W^* _test的具体过程为：

③_1、对G的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

同样，对R'的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

同样，对W'_test的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；

③_2、计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varG，varG的宽度为col且高度为row；

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varR，varR的宽度为col且高度为row；

同样，计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；计算中的每个区域块中的所有像素点的像素值的方差，得到对应的方差图像，记为varW，varW的宽度为col且高度为row；

然后计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGR，covGR的宽度为col且高度为row；同样，计算中的每个区域块中的所有像素点的像素值与中对应位置的区域块中的所有像素点的像素值之间的协方差，得到对应的协方差图像，记为covGW，covGW的宽度为col且高度为row；

③_3、计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^R(i,j)和b^R(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^R，a^R的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^R，b^R的宽度为col且高度为row；其中，1≤i≤col,1≤j≤row，covGR(i,j)表示covGR中坐标位置为(i,j)的像素点的像素值，varG(i,j)表示varG中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值，ξ为惩罚因子；

同样，计算中的每个区域块与中对应位置的区域块之间的第一系数和第二系数，将中左上角像素点的坐标位置为(i,j)的区域块与中左上角像素点的坐标位置为(i,j)的区域块之间的第一系数和第二系数记为a^W(i,j)和b^W(i,j)，然后将中的所有区域块与中对应位置的区域块之间的第一系数构成第一系数图像，记为a^W，a^W的宽度为col且高度为row；并将中的所有区域块与中对应位置的区域块之间的第二系数构成第二系数图像，记为b^W，b^W的宽度为col且高度为row；其中，covGW(i,j)表示covGW中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值；

③_4、对a^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；

同样，对b^R的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；

同样，对a^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；

同样，对b^W的四周边界进行扩充，且将扩充的像素点的像素值设定为0，将扩充得到的宽度为col+2且高度为row+2的图像记为然后采用尺寸大小为3×3的滑动窗口，在中将横坐标范围为[1,col]且纵坐标范围为[1,row]中的每个像素点对应滑动窗口的左上角进行逐像素点滑动，将划分为col×row个互相重叠的尺寸大小为3×3的区域块；接着计算中的每个区域块中的所有像素点的像素值的均值，得到对应的均值图像，记为 的宽度为col且高度为row；

③_5、获取R^*，将R^*中坐标位置为(i,j)的像素点的像素值记为R^*(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，G(i,j)表示G中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值；

同样，获取W^* _test，将W^* _test中坐标位置为(i,j)的像素点的像素值记为W^* _test(i,j)，其中，表示中坐标位置为(i,j)的像素点的像素值，表示中坐标位置为(i,j)的像素点的像素值；

④采用SPCE方法计算W^* _test与R^*的相关性，记为SPCE(W^* _test,R^*)；

所述的步骤④中其中，sign()为求符号函数，P_max表示W^* _test与R^*的相关函数图像P中最大的像素值，Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域的尺寸大小，8≤Ne≤15，Ω表示W^* _test与R^*的相关函数图像P中的所有像素点的坐标位置的集合，Ω_Ne×Ne表示W^* _test与R^*的相关函数图像P中以像素值为P_max的像素点为中心的窗口区域中的所有像素点的坐标位置的集合，P(s)表示W^* _test与R^*的相关函数图像P中坐标位置为s的像素点的像素值；

⑤用事先设置的I_test所声明的来源相机下的相关性判别阈值对SPCE(W^* _test,R^*)进行判决，若SPCE(W^* _test,R^*)大于或等于相关性判别阈值，则判定I_test来自其所声明的来源相机；若SPCE(W^* _test,R^*)小于相关性判别阈值，则判定I_test来自非其所声明的来源相机，即I_test来自其它相机。

2.根据权利要求1所述的一种使用空域平滑滤波的数字照片图像来源识别方法，其特征在于所述的步骤②中获取每幅原始数字照片图像的降噪图像和获取I_test的降噪图像所采用的方法为小波降噪处理方法；所述的步骤②中采用最大似然估计法对所有原始数字照片图像的噪声残差图像进行处理，获得R；所述的步骤②中采用零均值化法对R进行处理以消除R中的CFA插值噪声，获得R'；所述的步骤②中采用零均值化法对W_test进行处理以消除W_test中的CFA插值噪声，获得W'_test。

3.根据权利要求1所述的一种使用空域平滑滤波的数字照片图像来源识别方法，其特征在于所述的步骤③_3中取ξ＝0.75。

4.根据权利要求1所述的一种使用空域平滑滤波的数字照片图像来源识别方法，其特征在于所述的步骤⑤中相关性判别阈值的获取过程为：

⑤_1、使用I_test所声明的来源相机任意拍摄U幅同源图像；然后按照步骤②和步骤③的过程，以相同的方式获得每幅同源图像的最终的传感器模式噪声；其中，U≥100，每幅同源图像的宽度为col且高度为row；

⑤_2、选取与I_test所声明的来源相机不同的V部相机，要求选取的V部相机的分辨率大于或等于I_test所声明的来源相机的分辨率；然后使用选取的V部相机任意拍摄U幅异源图像，共得到U幅异源图像，再将每幅异源图像按照步骤⑤_1获得的同源图像的尺寸大小截取中心区域作为新的异源图像；接着按照步骤②和步骤③的过程，以相同的方式获得每幅新的异源图像的最终的传感器模式噪声；其中，1≤V≤9；

⑤_3、采用SPCE方法计算R^*与每幅同源图像的最终的传感器模式噪声的相关性，共得到U个同源样本相关性值；

同样，采用SPCE方法计算R^*与每幅异源图像的最终的传感器模式噪声的相关性，共得到U个异源样本相关性值；

⑤_4、在区间[0,30]内选取不同的实数值对步骤⑤_3中获得的U个同源相关性值和U个异源相关性值进行判决，获得每个实数值下的假正率(False Positive Rate，FPR)，记为FPR，其中，N_FP表示将异源图像判决为同源的数量，N_TN表示将异源图像判决为异源的数量；

⑤_5、将FPR＝0.01时对应的实数值作为相关性判决阈值。