一种基于特征融合的数字音频篡改自动检测方法与流程

本发明涉及数字音频信号处理技术领域，特别涉及一种基于特征融合的数字音频篡改自动检测方法。

背景技术：

数字音频信号相较于模拟音频信号有着易于传播、复制、存储等优点，但是这些优点在带给人们各种便利的同时也带来了不容忽视的问题。那就是篡改数字音频信号变得同样简单和快捷。数字音频的真实性和完整性受到来自各方面的威胁，且大多数的篡改操作是很难靠人的感官进行判别。数字音频取证技术就是用来验证数字音频信号的真实性和完整性的技术。数字音频被动取证是数字音频取证的一个重要分支，在司法取证、新闻公正和科学发现等领域有着广泛应用。

某一地区的电网频率(electricnetworkfrequency，enf)的波动在相当一段时间内具有稳定性和唯一性。电网频率的非周期性波动，对于接入同一电网的所有设备的影响是相同的，这就使得enf信号可以作为时间戳来使用。当录音设备采用电网供电时，录制的音频文件中就会残留有enf成分(enfcomponent，enfc)，经过带通滤波可以提取出enfc。利用enfc的稳定性和唯一性进行数字音频被动取证一般有两个研究思路，第一是将提取出的enfc与供电部门的电网频率数据库中的数据进行比对，确定音频录制时间与所宣称的是否一致，建立和保存大范围的enf信号数据库难度高、代价高，目前还没有实用价值比较高的enf数据库；第二是提取enf信号中的某些特征，进行一致性或者规律性分析。

grigoras最早提出基于enf的音频篡改检测算法，主要通过待检测音频中enf的波动和参考年份的数据进行比对，从而判断音频是否被篡改过。接着grigoras验证对音频信号加短时窗进行分析，可以与数据库进行更细致、精确的对比。提出建立一个enf准则，以期可以显示不同信号的供电类型，建议将该准则和传统的检测方法共用，为研究者建立标准化的研究分析方法。liu等提出一系列方法对enf数据库中的信号进行进一步处理以方便在语音取证中的使用。rodríguez等在grigoras的研究基础上，提出不需要使用enf标准数据库的方法，使用数字信号处理的方法估计enf信号相位，将enf相位变化的一致性作为特征来检测音频篡改，当特征值大小超过某一边界值即判断待测音频信号为篡改信号。esquef等根据篡改操作会引起篡改点enf瞬时频率突变，使用hilbert变换计算瞬时频率，提出tpsw(two-passsplit-window)方法估计enf背景变化水平，实际瞬时频率变化范围超过背景变化水平的部分的峰值点即为篡改点。

可以看出目前基于enf进行数字音频被动篡改检测的研究存在一些问题：1)没有权威的enf对比数据库。使用待测信号中的enf成分与enf数据库中进行对比来判断语音信号是否经过篡改没有可靠的结果2)部分检测方法是人从直观视觉上判断语音信号是否经过篡改，不能实现自动化。3)目前大多数的检测方法比较简略，没有囊括更多实际检测情况，不能达到数字音频被动检测的要求。

技术实现要素：

本发明的目的在于提出一种基于特征融合的数字音频篡改自动检测方法。该方法通过融合信号中的enf成分基于dft0的相位特征、基于dft1的相位特征和基于hilbert变换的瞬时频率特征，来分辨原始信号和编辑信号，实现了数字音频插入、删除操作的自动检测。该方法通过融合enf成分中具有代表性的相位和瞬时频率特征，并使用支持向量机进行分类，得到可以进行自动检测模型，提高了检测效率，实现了数字音频篡改检测的自动化。

本发明所采用的技术方案是：一种基于特征融合的数字音频篡改自动检测方法，包括以下步骤：

步骤1，对待测信号进行预处理，获得待测信号中的enf成分；

步骤2，对待测信号中的enf成分进行特征提取，包括以下子步骤，

步骤2.1，对enf成分进行基于dft⁰的相位估计，提取相位波动特征f1；

步骤2.2，对enf成分进行基于dft¹的相位估计，提取相位波动特征f2；

步骤2.3，对enf成分进行基于hilbert的瞬时频率估计,提取相位波动特征f3；

步骤3，对提取出的特征进行特征融合，得到特征集，将特征集进行归一化处理；

步骤4，利用优化的支持向量机对特征集进行分类，得到分类准确率。

进一步的所述步骤1的实现方式如下，

步骤1.1，将待测音频信号读入计算机得到离散的时间序列x[n]，其中fs是待测信号采样频率，i代表连续时间序列，n代表离散时间序列；对x[n]进行预处理，包括下采样、去直流分量，得到xd[n]，设下采样频率fd为1000hz或者1200hz；

步骤1.2，将步骤1.1中得到的xd[n]进行零填充，通过中心频率在enf标准频率处的带通滤波器，得到待测信号中的enf成分xenfg[n]。

进一步的，所述步骤2.1的实现方式如下，

首先，对xenfg[n]进行相位估计：对xenfg[n]信号进行分帧、加窗，分帧的长度是nc个标准enf频率周期，帧移为一个标准enf周期长度；应用一个平滑窗w(n)，得到加窗后的信号为xn[n]＝xenfg[n]w(n)；对xn[n]每帧进行n点的dft变换，得到xenfg[n]每帧加窗后的频域分布x(k)；对每帧信号进行dft变换后，令kpeak作为每帧的|x(k)|最大值的整数索引，则每帧enf信号的估计频率为其中ndft是dft变换的点数，fd是下采样频率，得到enf信号的相位序列

然后，使用特征量f1对基于dft⁰得到的enfc相位变化进行描述：设为相应第nb帧的估计相位，表示第nb帧|x(k)|对应最大值的整数索引，令设共有nblock帧，则相位波动特征f1的计算公式为，

其中，2≤nb≤nblock，表示从nb＝2到nblock的平均值。

进一步的，所述步骤2.2的实现方式如下，

首先，对xenfg[n]进行相位估计：计算enf信号xenfc[n]在点n处的近似一阶导数，

x′enfc[n]＝fd(xenfc[n]-xenfc[n-1])；

对一阶导数信号加窗，x′n[n]＝x′enfc[n]w(n)；对每帧x′n[n]进行n点的dft变换，得到每帧加窗后的频域分布x′(k)，找到|x′(k)|的峰值所在点k'peak，将|x′(k)|乘以一个尺度系数f(k)，其中这样可以得到，

dft⁰[k]＝|x(k)|和dft¹[k]＝f(k)|x′(k)|；

则每帧enf信号的估计频率值为得到enf信号的相位序列

其中klow和khigh分别定义为floor[a]代表小于a的最大整数，ceil[b]代表大于b的最小整数，θlow＝arg[x′(klow)]，θhigh＝arg[x′(khigh)]

然后，使用特征量f2对基于dft¹得到的enfc相位变化进行描述：设为相应第nb帧的估计相位，表示第nb帧|x′(k)|对应最大值的整数索引，令设共有nblock帧，则相位波动特征f2的计算公式为，

其中，2≤nb≤nblock，表示从nb＝2到nblock的平均值。

进一步的，所述步骤2.3的实现方式如下，

首先，对信号xenfc[n]进行离散hilbert变换，估计enf信号的瞬时频率f[n]：先利用离散hilbert变换得到xenfc[n]的解析函数x^(a)enfc[x]＝xenfc[x]+i*h{xenfc[x]}，其中h代表hilbert变换算子，瞬时频率是h{xenfc[n]}相角的变化率，设为f[n]，使用五阶椭圆滤波器iir滤波器对f[n]进行低通滤波去除寄生振荡，同时去掉f[n]头尾各2000采样点，最后得到的即为enfc的瞬时频率估计；

然后，使用特征量f3对enfc的瞬时频率波动进行描述：令len为xenfc[n]的长度，len＝length(xenfc[n])，为第n个采样点的瞬时频率，其中2≤n≤len，则相位波动特征f3的计算公式为，

其中，表示从n＝2到len的平均值。

进一步的，步骤4中通过交叉验证和网格搜索法选择优化的支持向量机参数，包括惩罚因子c和径向基函数参数g；

其中交叉验证和网格搜索法结合获得支持向量机参数的基本步骤是：首先输入c值的可能区间[cmin,cmax]和步进cstep，g值的可能区间[gmin,gmax]和步进gstep，并令最优的c值等于cmin，最优的g值等于gmin，最好预测准确率为零；接着使c值在[cmin,cmax]区间内以cstep为步进进行循环；在c值的循环内部，使g值在区间[gmin,gmax]内以gstep为步进进行循环；在g值循环内部，对每一组(c,g)进行支持向量机十折交叉验证，若所得最优分类准确率大于预测准确率则更新最优c值、g值和最优分类准确率；直到c值收敛，算法结束。

与现有技术相比，本发明的优点和有益效果：

(1)本发明融合了enf信号中对信号截断敏感的相位和瞬时频率特征，提高了数字音频被动篡改检测的准确率；

(2)本发明提出基于enf成分一致性的数字音频被动篡改检测的一般框架和enfc特征提取的一般框架，可以应用到其他基于enf信号的篡改检测方法中；

(3)本发明针对复杂环境录音和含噪语音稳定性高，具有很强的鲁棒性。

(4)本发明为数字音频被动篡改检测的准确性和自动化提出了一种广泛性的算法。

附图说明

图1是本发明实施例的流程图；

图2是本发明实施例的dft⁰相位特征提取流程图；

图3是本发明实施例的dft¹相位特征提取流程图；

图4是本发明实施例的hilbert变换瞬时频率特征提取流程图；

图5是本发明实施例的svm参数优化流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种基于特征融合的数字音频篡改自动检测方法，包括以下步骤：

步骤1：对待测信号进行预处理，获得待测信号中的enf成分；

具体实施包括以下子步骤：

步骤1.1：将待测音频信号读入计算机得到离散的时间序列x[n](fs是待测信号采样频率，i代表连续时间序列，n代表离散时间序列)，对x[n]进行预处理，包括下采样、去直流分量，得到xd[n]；

本实施例中考虑到频率混叠效应、信号信息损失和信号的信噪比(过采样可以提高信号的信噪比)的平衡，将信号的下采样频率fd定为1000hz或者1200hz，将标准的enf频率放在ω0＝π/10rad/sample。

步骤1.2：将步骤1.1中经过下采样的信号xd[n]进行零填充，通过中心频率在enf标准频率处的带通滤波器，得到信号中的enf成分xenfg[n]；

使用10000阶的线性零相位fir滤波器进行窄带滤波防止相位延时。中心频率在enf标准频率处，带宽为0.6hz，通带波纹为0.5db，阻带衰减为100db。使用高阶滤波器是为了得到理想的窄带信号。零填充(zeropadding)是指在时域信号的末端加上零点以增加信号长度的做法，在dft之前使用零填充可以提高频率分辨率，帮助更加准确地找到频谱上的峰值点。

步骤2：对待测信号中的enf成分进行特征提取；

具体实现包括以下子步骤：

步骤a1：对xenfg[n]进行基于dft⁰的相位估计，提取相位波动特征f1；

如图2，本实施例中基于dft⁰对xenfg[n]相位进行估计。对滤波后的xenfg[n]进行相位估计，首先对xenfg[n]信号进行分帧、加窗，分帧的长度是nc个标准enf频率周期，帧移为一个标准enf周期长度。应用一个平滑窗w(n)，得到加窗后的信号为xn[n]＝xenfg[n]w(n)，对xn[n]每帧进行n点的dft变换(ndft，dft点数大于enf信号的长度)，得到xenfg[n]每帧加窗后的频域分布x(k)。对每帧信号进行dft变换后，令kpeak作为每帧|x(k)|最大值的整数索引。所以每帧enf信号的估计频率为(ndft是dft变换的点数，fd是下采样频率)，得到enf信号的相位序列

使用特征量f1对基于dft⁰得到的enfc相位变化进行描述。令为相应第nb帧的估计相位，表示第nb帧|x(k)|对应最大值的整数索引，设共有nblock帧，则有2≤nb≤nblock，表示从nb＝2到nblock的平均值，下式则为相位波动特征f1的计算公式：

步骤a2：对xenfc[n]进行基于dft¹的相位估计，提取相位波动特征f2；

如图3，本实施例中基于dft¹对xenfc[n]相位进行估计。计算enf信号xenfc[n]在点n处的近似一阶导数：

x′enfc[n]＝fd(xenfc[n]-xenfc[n-1])

对一阶导数信号加窗：x′n[n]＝x′enfc[n]w(n)。对每帧x′n[n]进行n点的dft变换，得到每帧加窗后的频域分布x′(k)。同样找到|x′(k)|的峰值所在点k'peak(此处的k'peak与前面的kpeak不一定相同)。将|x′(k)|乘以一个尺度系数f(k)，

这样可以得到，

dft⁰[k]＝|x(k)|和dft¹[k]＝f(k)|x′(k)|

估计频率值为(分子中的k'peak和分母中的kpeak是各自的峰值点)。理想的情况下基于dft⁰得到的kpeak和dft¹得到的k'peak应该是相等的，且k'peak应该是最接近的整数，这样才是一个合理的频率值(validsolution)。为了得到一个合理的频率，需要满足以下式子：

如果定义验证条件也可以写作：

下面使用dft¹方法来估计enf信号的相位。考虑一个单频信号的模型stone(n)＝a(n)cos(ω0n+φ0)，n指离散时间变量，ω0是信号的角频率。信号的相位等于φ(n)＝ω0n+φ0，φ0是信号的初相。对于这个值的估计将限制在-π到π之间，φ(n)×n的图形将是一个锯齿状曲线(包裹相位wrappedphase)。这个模型是一个窄带信号，如果a(n)是一个常数，stone(n)将是一个窄带信号。在实践中，假定a(n)随着时间缓慢发展，因此可以认为a(n)在短时间或一帧内是个近似常量a。这个模型不包含任何随机部分或者宽带成分，因此可以应用到本发明的目标问题上。下采样信号经过带通滤波留下一个窄带信号，大部分在enf标准值以外的频率成分被过滤掉。因此每帧enf信号可以写成：xenfc[n]＝acos(ω0n+φ0)，其中ω0＝2πfenfc/fd，而fenfc是enf是实际频率。依据前面计算一阶导数的公式x′enfc[n]＝fs(xenfc[n]-xenfc[n-1])可得：

另外一个余弦信号的一阶导数实际上是另一个具有相同频率的正弦信号，上述公式可以表示为：

c是一个常数，θ是x′enfc的相位。比较上述两个式子可以得出：

ccos(θ)＝cos(φ0)-cos(φ0-ω)

csin(θ)＝sin(φ0)-sin(φ0-ω)

由上等式可以推导出：

上下同除以cos(φ0)，可以得到：

φ0代表xenfc的初始相位，可以使用dft¹方法估计出相位序列，即为：

其中对于θ的值，我们对x′(k)进行线性插值，令klow和khigh定义为floor[a]代表小于a的最大整数，ceil[b]代表大于b的最小整数。

由在(klow,θlow＝arg[x′(klow)])和(khigh,θhigh＝arg[x′(khigh)])进行线性插值可以逼近点求出的值与上式中θ的值保持一致：

对于上述他具有两个可能的值，如果是一个正值，可能在一个二维笛卡尔坐标系的第一或者第三象限，如果是个负值，可能在第二或者第四象限。可以使用作为一个参考，选择中最接近的值。

同样使用步骤a1中的方法计算出基于dft¹估计相位的波动特征。

使用特征量f2对基于dft¹得到的enfc相位变化进行描述。令为相应第nb帧的估计相位，表示第nb帧|x′(k)|对应最大值的整数索引，设共有nblock帧，则有2≤nb≤nblock，表示从nb＝2到nblock的平均值，下式则为相位波动特征f2的计算公式：

步骤a3：对xenfc[n]进行基于hilbert的瞬时频率估计，提取相位波动特征f3；

对信号xenfc[n]进行离散hilbert变换。首先得到xenfc[n]的解析函数：x^(a)enfc[x]＝xenfc[x]+i*h{xenfc[x]}，其中h代表hilbert变换算子。hilbert变换用于计算时域序列的瞬时属性是非常有用的，特别是振幅和频率，瞬时振幅是h{xenfc[n]}(即解析函数的虚部)的振幅，瞬时频率是h{xenfc[n]}相角的变化率。估计enf信号的瞬时频率f[n]。在使用hilbert变换的过程中由于有数值逼近，所以得到的f[n]存在一定的寄生振荡，需要进一步对f[n]进行低通滤波，去除振荡。使用五阶椭圆滤波器iir滤波器。中心频率在enf标准频率处，带宽为20hz，通带波纹为0.5hz，阻带衰减为64hz。由于频率估计的边界效应，去掉f[n]头尾各2000采样点，最后所得即为enfc的瞬时频率估计。

计算特征量f3对enfc的瞬时频率波动进行描述。令len为xenfc[n]的长度，len＝length(xenfc[n])。为第n个采样点的瞬时频率，其中2≤n≤len，表示从n＝2到len的平均值，下式则为相位波动特征f3的计算公式：

步骤3：对提取出的特征进行特征融合，得到特征集，将特征集进行归一化处理；

将步骤2中得到的3个特征进行特征融合，得到特征集。将特征集进行归一化处理，对特征集中的特征向量进行标记，分为原始信号和编辑信号。

步骤4：使用本发明优化支持向量机方法对特征集进行分类，得到分类准确率。

本发明实施例使用交叉验证和网格搜索法选择最佳支持向量机参数，包括惩罚因子c和径向基函数参数g，使用优化的支持向量机对特征集进行分类；

交叉验证是指将数据集平均分为k份，每次将其中k-1份数据作为训练数据，而另外一份数据作为测试数据。这样重复k次，根据k次迭代的平均值估计分类准确度。网格搜索是来确定两个参数即惩罚因子c与径向基函数参数g，基于网格法将c∈[cmin,cmax]，变化步长为cstep，而g∈[gmin,gmax]，变化步长为gstep。这样，针对每对参数(c,g)进行训练，取效果最好的一对参数作为模型参数。

交叉验证和网格搜索法结合获得支持向量机参数的基本步骤是：首先输入c值的可能区间[cmin,cmax]和步进cstep，g值的可能区间[cmin,cmax]和步进gstep。并令最优的c值等于cmin，最优的g值等于gmin，预测准确率为零；接着使c值在[cmin,cmax]区间内以cstep为步进进行循环；在c值的循环内部，使g值在区间[gmin,gmax]内以gstep为步进进行循环；在g值循环内部，对每一组(c,g)进行支持向量机十折交叉验证，若所得最优分类准确率大预测准确率则更新最优c值、g值和最优分类准确率；算法直到c值收敛循环结束，算法结束。算法流程图如图5所示，其中bestc是最优的c值，bestg是最优的g值，acc是预测准确率，bestacc是最优分类准确率。

本发明的实验结果，与进行特征融合之前的单个特征得到的实验结果进行对比，所得到的结果如下表1所示：

表1本发明实施例的融合特征与单个特征或者组合特征的检测结果对比表

本发明的实验结果，使用多个分类器与本发明所使用的优化svm进行对比，所得实验结果如下表2所示：

表2本发明实施例的优化svm分类器与其他分类器检测结果对比表

本发明的实验结果，从表1、2中可以看出，本方法所选用的特征和优化分类器可以达到更好的识别率。且本方法的实验数据来自于三个不同的数据库，因此可以认为本方法在大多数录音条件下都具有较好的检测效果。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。