基于音頻特征MFCC的零水印算法

宋　慧　李　晨　田麗華

(西安交通大學(xué)軟件學(xué)院　陜西 西安 710049)

?

基于音頻特征MFCC的零水印算法

宋慧李晨田麗華*

(西安交通大學(xué)軟件學(xué)院陜西 西安 710049)

利用小波域的低頻部分的奇異值具有比較穩(wěn)定的性能，設(shè)計(jì)一種以音頻特征為基礎(chǔ)的零水印方案。該方案利用原始音頻特征——Mel頻率倒譜系數(shù)(MFCC)，一方面來(lái)構(gòu)造水印，另一方面還被用于音頻內(nèi)容認(rèn)證。先將原始音頻做分段操作，取各分段進(jìn)行DWT變換；選取小波域的低頻部分，并做SVD操作；最后利用奇異值平均值實(shí)現(xiàn)水印的“嵌入”，獲得一個(gè)位置索引向量并據(jù)此提取水印信息。仿真實(shí)驗(yàn)顯示，該方案在完全不影響音頻聽(tīng)覺(jué)質(zhì)量的前提下，在抵抗常見(jiàn)音頻攻擊時(shí)有很強(qiáng)的魯棒性，并能進(jìn)行惡意篡改定位。

MFCC零水印奇異值分解篡改定位

0　引　言

近年來(lái)，由于計(jì)算機(jī)技術(shù)和因特網(wǎng)技術(shù)的飛速發(fā)展，促使圖像、視頻和音頻等各種各樣形式的多媒體的產(chǎn)生、保存及傳輸都變得越來(lái)越容易。但由此帶來(lái)的安全問(wèn)題，如非法復(fù)制和盜版問(wèn)題日益突出。在這種背景下，能夠很好解決版權(quán)問(wèn)題的數(shù)字水印被提出，已經(jīng)變成當(dāng)前的研究熱點(diǎn)之一[1]。

目前提出的算法，多數(shù)是通過(guò)修改原始音頻的時(shí)域或變換域的信息來(lái)嵌入水印，雖可實(shí)現(xiàn)版權(quán)保護(hù)，但都是以犧牲音頻的質(zhì)量為代價(jià)的。為了不影響音頻音質(zhì)，即實(shí)現(xiàn)水印的透明性，在這種應(yīng)用背景下，零水印技術(shù)首次被提出。零水印算法[2,3]是無(wú)需改變?cè)瓉?lái)音頻的任意數(shù)值的算法，因此能夠很好地處理數(shù)字水印的透明性及穩(wěn)健性之間的關(guān)系。

國(guó)內(nèi)已有一些學(xué)者針對(duì)零水印展開(kāi)了相關(guān)研究。文獻(xiàn)[4]所提出的算法是將音頻的小波域特征提取出并對(duì)其通過(guò)Hash變換轉(zhuǎn)為一組二值序列，將其“嵌入”到原來(lái)的音頻中。雖然能夠?qū)崿F(xiàn)數(shù)字音頻的驗(yàn)證，但嵌入的數(shù)字簽名是無(wú)意義的，這樣就不能直觀給出能夠代表版權(quán)的信息。文獻(xiàn)[2]提出了一種基于DWT和SVD的零水印算法，該算法的思路是先將音頻分幀，計(jì)算每幀的短時(shí)能量和過(guò)零率；然后通過(guò)預(yù)先設(shè)定的閾值篩選出特征顯著的音頻幀，并對(duì)選定的每幀進(jìn)行DWT變換，取DWT變換的低頻部分再進(jìn)行SVD；最后通過(guò)分解后得到的最大奇異值來(lái)生成水印信息。該算法雖然對(duì)常見(jiàn)音頻攻擊具有一定的魯棒性，但是由于需要事先確定出短時(shí)能量和過(guò)零率的閾值，這兩個(gè)閾值的選取對(duì)算法效果影響顯著，算法穩(wěn)定性欠佳。文獻(xiàn)[5,6]提出了基于音頻特征MFCC的水印方法，雖都可對(duì)音頻版權(quán)和內(nèi)容進(jìn)行認(rèn)證，但嵌入水印過(guò)程需將原始音頻內(nèi)容進(jìn)行適當(dāng)修改，這樣必然損害音頻的音質(zhì)，對(duì)水印的不可感知性帶來(lái)挑戰(zhàn)。對(duì)于前面提出的難題，參考文獻(xiàn)[4-6]中提出的方法，將小波變換、音頻信號(hào)特征及奇異值分解(SVD)[7-10]的特點(diǎn)相結(jié)合，根據(jù)小波域低頻系數(shù)奇異值的優(yōu)良特性，本文設(shè)計(jì)了音頻自身特性為基礎(chǔ)的零水印算法。該算法提取音頻特征MFCC參數(shù)，一方面和代表版權(quán)信息的二值圖像做異或運(yùn)算，作為最終的水印“嵌入”原始音頻中，得到一個(gè)位置索引向量并據(jù)此提取水印；另一方面，MFCC還被用在音頻內(nèi)容的認(rèn)證。該算法不用更改原來(lái)音頻的任何數(shù)據(jù)，使其有很強(qiáng)的透明性。仿真實(shí)驗(yàn)顯示，該方案在完全不影響音頻聽(tīng)覺(jué)質(zhì)量的前提下，在抵抗常見(jiàn)音頻攻擊時(shí)有很強(qiáng)的魯棒性，并能進(jìn)行惡意篡改定位。

1　理論基礎(chǔ)

1.1奇異值分解

奇異值分解(SVD)經(jīng)常被用在線性代數(shù)里，是特別有效的運(yùn)算工具，可以將矩陣進(jìn)行對(duì)角化操作[11]。B是一個(gè)矩陣，那么對(duì)其進(jìn)行SVD操作，計(jì)算公式如下：

(1)

將S中對(duì)角線的值叫做奇異值，具有以下特點(diǎn)：(1) 穩(wěn)定性很好，對(duì)各種常見(jiàn)的信號(hào)處理，對(duì)角矩陣不會(huì)明顯改變；(2)S體現(xiàn)出來(lái)的特點(diǎn)為矩陣序列具有的專(zhuān)有固性；(3)S中的第一個(gè)數(shù)是S中最大的數(shù)，而且它對(duì)矩陣的作用最大。基于SVD的以上突出特點(diǎn)，對(duì)SVD的研究成為現(xiàn)在的國(guó)內(nèi)外研究熱點(diǎn)。SVD首先在圖像水印中被使用，最近幾年用在音頻水印上的研究也比較多。

1.2MFCC特征參數(shù)提取

MFCC是基于人耳的聽(tīng)覺(jué)機(jī)理，并將其與言語(yǔ)的生成機(jī)制結(jié)合起來(lái)，利用一種非線性的頻率單位(Mel頻率)來(lái)模擬人耳的聽(tīng)覺(jué)系統(tǒng)，能很好反映出人耳的感知特性，當(dāng)前普遍用于音頻剖析、辨別及音頻數(shù)字信號(hào)的建模。MFCC的提取步驟[12]如圖1所示。

圖1　MFCC特征提取流程圖

1) 取音頻X(n)，先做預(yù)加重操作，再做分幀處理。其次對(duì)每幀進(jìn)行加窗操作，便可以獲得預(yù)處理操作的音頻。

2) 對(duì)處理后的每一幀音頻做快速傅里葉操作，則時(shí)域數(shù)據(jù)轉(zhuǎn)為頻域數(shù)據(jù)，進(jìn)而計(jì)算得到功率譜。

3) 設(shè)計(jì)若干帶通濾波器(一般值取24，用M表示)，將得到各幀的功率譜經(jīng)過(guò)濾波器，得到每一個(gè)所對(duì)應(yīng)的對(duì)數(shù)能量譜S(m)。

(2)

S(m)=ln(e(m))0

(3)

式中，X(k)是原始音頻做快速傅里葉操作后的頻域數(shù)據(jù)，Hm(k)為三角濾波器。

4) 將S(m)進(jìn)行離散余弦操作，最終便可以獲得MFCC參數(shù)C(n)：

(4)

每一個(gè)音頻幀將對(duì)應(yīng)幾個(gè)MFCC參數(shù)，比較重要的音頻信息處于MFCC的C1至C12的參數(shù)內(nèi)[13]。通常，前面階數(shù)的MFCC包含更多有用的內(nèi)容，一般選擇前面幾階。在綜合計(jì)算量及水印嵌入量的情況下，本文選擇第一階MFCC。

2　水印生成

2.1MFCC特征選取

為了能夠在水印體系里比較充分地使用音頻的自身特性，本文由提取得到的MFCC來(lái)生成二值特征序列。該算法中，MFCC參數(shù)，一方面與具有實(shí)際版權(quán)意義的二值圖像進(jìn)行異或運(yùn)算，作為水印“嵌入”到原始音頻中；另一方面還被用于音頻內(nèi)容的認(rèn)證。為了更好地證明第一階MFCC針對(duì)原始音頻攻擊非常敏感，對(duì)其進(jìn)行以下攻擊操作，如圖2所示。可以得出，對(duì)音頻的多項(xiàng)處理，MFCC的一階參數(shù)都非常敏感。

圖2　音頻處理后一階MFCC的變化

2.2水印信號(hào)生成

二值特征序列可以看作驗(yàn)證音頻文件的“條形碼”。為了能夠讓此“條形碼”比較直觀地進(jìn)行表現(xiàn)，對(duì)二值特征序列和Arnold置亂處理的有意義二值圖像序列做異或操作，便能獲得一個(gè)既含有音頻自身特征又具有版權(quán)意義的水印序列。具體步驟如下：

1) 選取一個(gè)具有實(shí)際意義的圖像系列，采用Arnold變換[9]對(duì)其實(shí)行置亂處理，這樣可以有效提高該算法的安全性。因?yàn)橐纛l實(shí)際上是一個(gè)一維序列，水印是二維系列，所以加密操作后的水印序列必須進(jìn)行降維操作，變成一個(gè)一維序列，用V表示：

V={V(i),1≤i≤M×N}

(5)

2) 為獲取音頻分段的MFCC，首先將音頻X做分段處理，并根據(jù)前面所述MFCC提取的步驟，得到每個(gè)分段的第一階MFCC；然后將其放于f(i)中，1≤i≤n，i為幀號(hào)，n為總段數(shù)，M×N≤n，并將f(i)作為密鑰k1。

(6)

4) 將二值特征序列w(i)與預(yù)處理操作后的二值圖像序列V做異或操作，便獲取含有音頻特征的序列Wz，并將其作為待嵌入的水印信息。

3　水印的嵌入和提取

由于DWT具有多尺度和多分辨率等優(yōu)點(diǎn)，以及SVD所得的奇異值在信號(hào)處理中具有極強(qiáng)穩(wěn)定性的特性，本文以此為基礎(chǔ)，利用DWT低頻系數(shù)奇異值所具有的穩(wěn)定特性，設(shè)計(jì)了一種以音頻自身特性為基礎(chǔ)的零水印算法。本文算法思路是：先根據(jù)音頻自身特征來(lái)構(gòu)造水印，再將原音頻做分幀處理，并對(duì)每幀實(shí)行DWT變換操作；取DWT變換的低頻部分，進(jìn)行SVD處理；最后使用奇異值平均值實(shí)現(xiàn)“嵌入”過(guò)程，獲得一個(gè)位置索引向量并據(jù)此實(shí)現(xiàn)水印提取。本文算法流程如圖3和圖4所示。

圖3　音頻水印嵌入流程圖

圖4　音頻水印提取流程圖

3.1水印的嵌入過(guò)程

1) 原始音頻分段。按照上節(jié)敘述的水印生成的詳細(xì)步驟，生成二值序列wz(i),其中1≤i≤n(n是水印大小)。然后把音頻X等分為n段，k為音頻幀。

X={X(k),1≤k≤n}

(7)

2) 離散小波變換。對(duì)于每一個(gè)原始音頻幀做H級(jí)小波變換操作，便可以取得低頻部分CAH及高頻部分。

3) 奇異值分解。因?yàn)閷?duì)音頻進(jìn)行DWT操作后，音頻的主要能量集中存儲(chǔ)于CAH中。所以，提取每個(gè)音頻分段的CAH，并將其轉(zhuǎn)換成矩陣Dj，1≤j≤n，再對(duì)矩陣Dj實(shí)行SVD操作。提取各矩陣的奇異值平均值組成數(shù)組S(i)，然后對(duì)S(i)進(jìn)行升序排列，組成升序奇異值平均值數(shù)組M(i)，1

4) 獲取臨界值。得出水印里0和1的數(shù)值并用C0、C1表示，然后得出他們的比值r=C0/C1,則t=?r/(r+1)×n」，n為M(i)里元素的數(shù)目。臨界值T為M(i)中的第t個(gè)元素。將T看成密鑰k2存儲(chǔ)。

5) 分類(lèi)。根據(jù)閾值T對(duì)升序排列的奇異平均值數(shù)組M(i)按照式(8)做分類(lèi)操作，將M(i)分成A類(lèi)和B類(lèi)，然后依次記下A和B類(lèi)中的數(shù)值處于原音頻X(i)中的段號(hào)i。

A：M(i)>TB：M(i)≤T

(8)

6) “嵌入”水印。根據(jù)水印比特值為1或0，將其分別映射到A、B類(lèi)里，并在A、B類(lèi)里選擇相對(duì)應(yīng)數(shù)值；然后記載此數(shù)值處在原來(lái)音頻X中的段號(hào)數(shù)值；最后產(chǎn)生一個(gè)位置索引，并用N表示。將索引向量作為密鑰進(jìn)行存儲(chǔ)，用k3表示，水印嵌入操作完成。在提取水印的時(shí)候，嵌入的索引向量N及臨界值T要被發(fā)送到接收器中。

3.2水印的提取過(guò)程

1) 提取水印時(shí)，首先對(duì)音頻進(jìn)行同嵌入水印時(shí)一樣的操作，獲取各矩陣所對(duì)應(yīng)的奇異值均值，然后組成數(shù)組S。

2) 根據(jù)接收到的嵌入索引向量N(k3)提取數(shù)組S中相應(yīng)位置對(duì)應(yīng)的奇異值均值S(i)，并讓S(i)與接收到的臨界值T(k2)進(jìn)行比較，按照下式進(jìn)行水印比特的提取：

(9)

3) 取出WZ′及特征序列w(i)并實(shí)行異或操作，獲取一維數(shù)組V′(i)。再取其依次進(jìn)行升維操作及Arnold變換，便可以獲得提取的水印圖像。

4　仿真實(shí)驗(yàn)及分析

為測(cè)試本文所提出的水印方案的性能，使用Matlab軟件進(jìn)行了仿真測(cè)試。在測(cè)試過(guò)程中采用采樣頻率是44 100Hz、量化精度為16bit、段長(zhǎng)為7s的音頻信號(hào)進(jìn)行測(cè)試，并選取比較常用的db1小波基對(duì)每個(gè)音頻分段實(shí)行3級(jí)DWT操作。實(shí)驗(yàn)中選用的原來(lái)信號(hào)的波形如圖5所示，有意義的二值圖像如圖6所示，其大小為48×48。

圖5　原始音頻信號(hào)波形圖

圖6　二值水印圖像

4.1魯棒性測(cè)試

為體現(xiàn)版權(quán)保護(hù)的客觀公正，取比特誤碼率(BER)及歸一化相關(guān)系數(shù)(NC)作為本文所提出的算法魯棒性的評(píng)價(jià)指標(biāo)，其計(jì)算公式如下：

(10)

(11)

式中，w代表原水印的信息，w′代表所提取水印的信息，E代表所提取水印錯(cuò)誤的比特?cái)?shù)，M1×M2代表水印大小。

為測(cè)試本文所提出方案的魯棒性，做了如下幾種攻擊操作：(1) 添加信噪比為25dB的高斯噪聲；(2) 添加信噪比為20dB的高斯噪聲；(3) 經(jīng)過(guò)截止頻率是10 000Hz的Butterworth數(shù)字濾波；(4) 上重采樣，速度為88.2kHz；(5) 下重采樣，速度為22.05kHz；(6) 下重量化，量化長(zhǎng)度為8bit；(7) 上重量化，量化長(zhǎng)度為32bit；(8)MP3壓縮，速率為128kb/s。經(jīng)由上述攻擊操作提取的各個(gè)水印如圖7所示，其中(i)是無(wú)攻擊操作時(shí)提取的結(jié)果。

圖7　提取的各水印圖像

從圖7可以較明顯發(fā)現(xiàn)，“嵌入”水印音頻遭到如上攻擊操作時(shí)，提取出的水印清晰可辨。尤其對(duì)于重采樣、加噪、重量化操作得到的水印圖像與原來(lái)水印圖像相比較，肉眼看不出差異之處，從而表明文中提出的水印方法針對(duì)上面的攻擊操作擁有比較不錯(cuò)的魯棒性。

表1　攻擊后的所提取水印的BER和NC

表1中給出了經(jīng)過(guò)上述攻擊后的NC和BER值，并將本文水印算法與文獻(xiàn)[2]和文獻(xiàn)[4]的算法進(jìn)行了比較測(cè)試。通過(guò)比較發(fā)現(xiàn)，本文水印算法在受到前文所列的攻擊時(shí)，NC均大于文獻(xiàn)[2]和文獻(xiàn)[4]算法，BER均低于文獻(xiàn)[2]和文獻(xiàn)[4]算法。主要是由于文獻(xiàn)[2]中，生成水印的閾值是預(yù)先設(shè)定的，具有一定的隨機(jī)性，對(duì)算法性能有較大影響。實(shí)驗(yàn)發(fā)現(xiàn)，有部分音頻段經(jīng)過(guò)攻擊操作后，其短時(shí)能量和過(guò)零率會(huì)發(fā)生變化，采用事先選定的閾值進(jìn)行特征音頻段的篩選會(huì)增加水印信息的誤碼率。文獻(xiàn)[4]是以DWT為基礎(chǔ)的算法，本文算法是結(jié)合了DWT及SVD的算法。由于同時(shí)利用了DWT多尺度、多分辨率等優(yōu)點(diǎn)及SVD所得奇異值均值穩(wěn)定性強(qiáng)的優(yōu)點(diǎn)，使得本文算法提取水印的相似度有所提高，誤碼率相對(duì)較小。實(shí)驗(yàn)結(jié)果也印證了本文算法較之前算法具有更強(qiáng)的魯棒性。

4.2篡改定位測(cè)試

如果音頻的內(nèi)容沒(méi)有發(fā)生惡意的篡改，那么原來(lái)的音頻同接收到的音頻的MFCC值變化不大[6]。因此通過(guò)比較提取的MFCC及接收端獲取到的密鑰k1，來(lái)確認(rèn)音頻有沒(méi)有發(fā)生篡改，還可以定位發(fā)生攻擊的區(qū)域。

采用數(shù)組w1(i)來(lái)存儲(chǔ)原始音頻的MFCC，w2(i)來(lái)存儲(chǔ)接收的信號(hào)的MFCC，m(i)為原音頻X的MFCC同收到的音頻的MFCC做差值并取絕對(duì)值操作的結(jié)果。

(12)

式中：i代表語(yǔ)音第i幀。m(i)的值趨近于0，說(shuō)明兩幀的變化不大，相反則認(rèn)為該幀已經(jīng)發(fā)生改變。通過(guò)m(i)可以確定不同幀所處位置。

為測(cè)試本文算法針對(duì)惡意篡改的檢測(cè)和定位方面的功能，本文中主要進(jìn)行了下列實(shí)驗(yàn)操作。篡改類(lèi)型 1 將音頻的任意采樣值進(jìn)行靜音操作，篡改類(lèi)型 2將測(cè)試音頻中的部分采樣值用其他部分的采樣值進(jìn)行替換操作。圖8(a)為將音頻10 000～30 000之間的采樣點(diǎn)靜音的波形，圖8(b)是把測(cè)試音頻信號(hào)70 001～100 000的采樣點(diǎn)用測(cè)試音頻的100 001～130 000的采樣點(diǎn)替換的波形。經(jīng)過(guò)篡改操作的m值和相應(yīng)的位置定位如圖9所示，測(cè)試表明，本文算法能夠?qū)崿F(xiàn)對(duì)篡改攻擊的音頻段的定位。

圖8　篡改攻擊后的音頻信號(hào)

圖9　經(jīng)過(guò)篡改操作的m值和位置定位

5　結(jié)　語(yǔ)

本文所提出的方案是把離散小波變換同奇異值分解進(jìn)行結(jié)合，以Mel頻率倒譜系數(shù)獲得音頻的特征，并將其特征和二值圖像異或后作為最終的水印。且該算法不更改原來(lái)音頻任何數(shù)據(jù)，具有很強(qiáng)的透明性。仿真實(shí)驗(yàn)表明，該算法在完全不影響音頻聽(tīng)覺(jué)質(zhì)量的前提下，對(duì)多種攻擊操作具備較強(qiáng)的魯棒性，并可以定位惡意篡改。

[1]RossiF,VillaN.Supportvectormachineforfunctionaldataclassification[J].Neurocomputing,2006,69(7-9):730-742.

[2] 任克強(qiáng),李慧,謝斌.基于DWT和SVD的自適應(yīng)數(shù)字音頻零水印算法[J].計(jì)算機(jī)工程與應(yīng)用,2014,50(12):87-90,95.

[3] 李慧.基于變換域的數(shù)字音頻零水印算法研究[D].江西:江西理工大學(xué),2012.

[4] 吳翔,楊曉元,張敏情,等.基于特征的數(shù)字音頻簽名算法[J].計(jì)算機(jī)工程與應(yīng)用,2009,45(10):106-107,171.

[5] 寇占奎,徐江峰.基于音頻特征的半脆弱水印方案[J].計(jì)算機(jī)工程與設(shè)計(jì),2012,33(9):3323-3326,3341.

[6] 林曉丹.用于版權(quán)和內(nèi)容認(rèn)證的音頻水印[J].計(jì)算機(jī)應(yīng)用,2010,30(8):2269-2271.

[7] 劉泉,高慶.基于DWT的自適應(yīng)盲音頻水印算法[J].武漢理工大學(xué)學(xué)報(bào),2008,30(1):95-98.

[8] 雷敏,楊榆.基于DWT-DCT-SVD的音頻盲水印算法[J].北京郵電大學(xué)學(xué)報(bào),2011,34(S1):51-54.

[9] 段歲軍,范九倫.一種基于SVD和DWT的音頻水印算法[J].計(jì)算機(jī)應(yīng)用研究,2014,31(7):2116-2118.

[10] 楊志疆,葉阿勇.多重水印嵌入的DWT-SVD音頻公開(kāi)水印算法[J].系統(tǒng)仿真技術(shù),2014,10(1):27-31,48.

[11]WangJ,HealyR,TimoneyJ.Arobustaudiowatermarkingschemebasedonreducedsingularvaluedecompositionanddistortionremoval[J].SignalProcessing,2011,91(8):1693-1708.

[12] 宋知用.MATLAB在語(yǔ)音信號(hào)分析與合成中的應(yīng)用[M].北京:北京航空航天大學(xué)出版社,2013:37-43.

[13] 甄斌,吳璽宏,劉志敏,等.語(yǔ)音識(shí)別和說(shuō)話人識(shí)別中各倒譜分量的相對(duì)重要性[J].北京大學(xué)學(xué)報(bào):自然科學(xué)版,2001,37(3):371-378.

ZERO-WATERMARKINGALGORITHMBASEDONAUDIOFEATURESOFMFCC

SongHuiLiChenTianLihua*

(SchoolofSoftwareEngineering,Xi’anJiaotongUniversity,Xi’an710049,Shaanxi,China)

Byusingthestabilityofsingularvalueoflowfrequencycoefficientsonwaveletdomain,wedesignedazero-watermarkingschemewhichtakesthefeaturesofaudioasthebasis.Theschemeusesthefeatureoftheoriginalaudio—melfrequencycepstrumcoefficients(MFCC)toconstructthewatermarksequence,inadditionthisfeatureisalsoappliedtoaudiocontentauthentication.Inthealgorithm,originalaudioisdividedintosectionsatfirst,andeachsectioniscarriedoutthediscretewavelettransform.Thenthelowfrequencypartsofwavelettransformareselectedforsingularvaluedecomposition.Finallythemeanofsingularvalueisusedtoimplementthe"embedment"ofwatermarksinoriginalaudio,andthisobtainsapositionindexvector,itcanbeusedtoextractthewatermarkinformation.Simulationexperimentindicatedthattheschemepresentedinthispaperhasquitestrongrobustnessonresistingsomepopularattacksandcanlocatethemalicioustamperingonthepremiseofwithoutlosingauditoryquality.

MFCCZero-watermarkingSingularvaluedecompositionTamperlocalisation

2015-07-01。國(guó)家自然科學(xué)基金項(xiàng)目(61403302)。宋慧，碩士生，主研領(lǐng)域：數(shù)字音頻水印。李晨，講師。田麗華，高工。

TP

ADOI:10.3969/j.issn.1000-386x.2016.10.031