基于頻率切片小波變換的時(shí)頻分析與MT信號(hào)去噪

蔡劍華,熊 銳

(1.湖南文理學(xué)院物理與電子科學(xué)學(xué)院,湖南常德415000;2.湖南文理學(xué)院洞庭湖生態(tài)經(jīng)濟(jì)區(qū)建設(shè)與發(fā)展湖南省協(xié)同創(chuàng)新中心,湖南常德415000)

?

基于頻率切片小波變換的時(shí)頻分析與MT信號(hào)去噪

蔡劍華1,2,熊 銳1,2

(1.湖南文理學(xué)院物理與電子科學(xué)學(xué)院,湖南常德415000;2.湖南文理學(xué)院洞庭湖生態(tài)經(jīng)濟(jì)區(qū)建設(shè)與發(fā)展湖南省協(xié)同創(chuàng)新中心,湖南常德415000)

大地電磁(MT)勘測(cè)方法因存在多種人文噪聲使得某些頻段受噪聲影響顯得更為明顯,提出了一種基于頻率切片小波變換時(shí)頻分析的時(shí)頻域MT信號(hào)去噪新方法。給出了方法原理和實(shí)現(xiàn)步驟,先對(duì)MT信號(hào)進(jìn)行頻率切片小波變換,得到全頻帶下的時(shí)頻分布,在此基礎(chǔ)上對(duì)時(shí)頻譜進(jìn)行時(shí)頻閾值濾波和逆變換，重構(gòu)分離出去噪后的MT信號(hào)。仿真實(shí)驗(yàn)和實(shí)測(cè)信號(hào)分析結(jié)果表明,該方法可有效消除大地電磁信號(hào)中的強(qiáng)噪聲;去噪后計(jì)算的響應(yīng)參數(shù)曲線(xiàn)的突變點(diǎn)得到了有效抑制,曲線(xiàn)變得平滑、連續(xù)。

頻率切片小波變換;時(shí)頻分析;大地電磁信號(hào);去噪

大地電磁測(cè)深(Magnetotelluric,MT)方法是油氣田普查勘探、巖石圈結(jié)構(gòu)探測(cè)等領(lǐng)域的一種重要手段。在石油天然氣勘探中,特別是在地震勘探困難區(qū),電磁法勘探技術(shù)已成為首選的非地震勘探技術(shù),被廣泛用于石油綜合地質(zhì)調(diào)查,如勘查斷裂破碎帶、查明盆地邊界、基底埋深和基本構(gòu)造格架等,為油氣田整體評(píng)價(jià)提供科學(xué)依據(jù)[1-2]。然而,大地電磁測(cè)深數(shù)據(jù)由于其信號(hào)自身的特點(diǎn)和日益嚴(yán)重的人文噪聲的干擾,極易受到污染,影響響應(yīng)參數(shù)的穩(wěn)健估計(jì)[3-4]。如何有效消除噪聲是大地電磁工作者關(guān)心的熱點(diǎn)問(wèn)題,國(guó)內(nèi)外學(xué)者提出了一系列大地電磁信號(hào)的去噪方法,如:Robust處理[5-6]、遠(yuǎn)參考技術(shù)[7]、高階統(tǒng)計(jì)量[8]、小波變換[4,9]和基于經(jīng)驗(yàn)?zāi)B(tài)分解[10-12]的濾波方法等。其中,經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法由于具有良好的自適應(yīng)特征,在大地電磁信號(hào)去噪領(lǐng)域得到廣泛應(yīng)用,但也存在模態(tài)混疊和端點(diǎn)效應(yīng)等問(wèn)題[10]。2009年,YAN等[13]提出的時(shí)頻分析新技術(shù)——頻率切片小波變換(Frequency Slice Wavelet Transform,FSWT),通過(guò)引入頻率切片函數(shù)使傳統(tǒng)傅里葉變換實(shí)現(xiàn)了時(shí)頻分析功能。近年來(lái),FSWT被應(yīng)用到模式識(shí)別和爆破振動(dòng)信號(hào)分析中,展示了其較好的時(shí)頻分析和重構(gòu)能力[13-15],但頻率切片小波變換在大地電磁信號(hào)分析領(lǐng)域的應(yīng)用鮮有報(bào)道。我們將頻率切片小波變換的時(shí)頻分析算法引入到大地電磁信號(hào)的時(shí)頻特征分析中,探討基于頻率切片小波變換的大地電磁信號(hào)時(shí)頻域去噪新方法。

1 基本原理

基于頻率切片小波變換的時(shí)頻域去噪方法是先對(duì)信號(hào)進(jìn)行頻率切片小波變換，得到其時(shí)頻分布，在時(shí)頻域內(nèi)進(jìn)行閾值濾波后，再由頻率切片小波變換的逆變換重構(gòu)去噪后的信號(hào)。

1.1 頻率切片小波變換

(1)

設(shè)尺度因子σ=ω/k,k>0,則:

(2)

式中:k為時(shí)頻分辨系數(shù)[14-15],用來(lái)調(diào)整時(shí)頻域?qū)ψ儞Q的響應(yīng)靈敏度,與ω,u無(wú)關(guān)。

通常,用頻率分辨比率η和幅值期望響應(yīng)比率ν這2個(gè)系數(shù)來(lái)評(píng)價(jià)分析信號(hào),以獲得較高的時(shí)頻分辨率。

(3)

對(duì)于f(t)=eiω0t,若|W(t,ω0+Δω,λ,σ)|/|W(t,ω0,λ,σ)|≤ν,則:

(4)

對(duì)于f(t)=δ(t-t0),若|W(t0+Δt,ω,λ,σ)|/|W(t0,ω,λ,σ)|≤ν,則:

(5)

1.2 頻率切片小波變換的逆變換

理論上,頻率切片小波變換的時(shí)頻域分解是冗余的，其逆變換也可以用不同的形式,其中最為簡(jiǎn)潔的一種可表示為:

(6)

公式(6)表明頻率切片小波變換的逆變換只與參數(shù)k有關(guān),而與函數(shù)p(ω)無(wú)關(guān)。當(dāng)k為定值時(shí),公式(6)為傅里葉逆變換。f(t)的頻率切片小波變換為W(t,ω,σ),則在時(shí)頻區(qū)域(t1,t2,ω1,ω2)的信號(hào)分量為:

(7)

顯然,在f(t)的頻率切片小波變換時(shí)頻區(qū)間內(nèi),時(shí)頻區(qū)域(t1,t2,ω1,ω2)即時(shí)頻切片可任意選擇。所以,可根據(jù)選擇的時(shí)頻切片在時(shí)頻空間上隨意提取所需的信號(hào)分量[13,15]。

1.3 去噪方法與步驟

基于頻率切片小波變換時(shí)頻分析的大地電磁信號(hào)分析方法的實(shí)現(xiàn)過(guò)程如下。

1) 采用FSWT進(jìn)行大地電磁信號(hào)分解,得到在全頻帶下的時(shí)頻分布。

(8)

3) 對(duì)濾波后的時(shí)頻區(qū)域進(jìn)行逆變換重構(gòu)分離出的有效MT信號(hào)。

4) 進(jìn)一步對(duì)去噪后的MT信號(hào)進(jìn)行功率譜計(jì)算和響應(yīng)函數(shù)估計(jì)。

2 仿真實(shí)驗(yàn)

為了驗(yàn)證上述時(shí)頻濾波方法的正確性及其在處理MT信號(hào)方面的優(yōu)勢(shì),我們對(duì)一組來(lái)自嚴(yán)家斌的仿真數(shù)據(jù)[9]進(jìn)行了處理,噪聲為類(lèi)似工頻干擾的正弦信號(hào)。在此,仿真信號(hào)表示為s(t),噪聲是n(t),則含噪的仿真MT信號(hào)f(t)可表示為:

f(ti)=s(ti)+n(ti) i=1,2,…,n

(9)

圖1 仿真信號(hào)的去噪a 仿真MT信號(hào); b 加噪的仿真信號(hào); c 加噪后的時(shí)頻譜; d 本文方法去噪后的信號(hào); e EMD時(shí)空濾波法去噪后的信號(hào)

SNR=10lg(PS/Pn)

(10)

(11)

(12)

表1給出了不同信噪比情況下幾個(gè)參數(shù)的對(duì)比結(jié)果。從表1中可以看出,在不同的SNR條件下,不管是本文提出的時(shí)頻域去噪方法還是基于EMD的方法,信號(hào)失真比(SDR)都不大,且前者的SDR更低。就噪聲抑制率來(lái)說(shuō),本文提出的時(shí)域方法也優(yōu)于EMD去噪方法。在較高的SNR環(huán)境下,噪聲很弱,所以NSR也相對(duì)較小(最大值僅是0.2587),對(duì)應(yīng)的SDR也很小(最大值僅是0.0484)。隨著信噪比的增加,NSR也相對(duì)增大(從0.2587增加到0.7950),盡管SDR也增大了(從0.0209增大到0.0420),但是仍舊保持很小的值。這就說(shuō)明本文提出的時(shí)域?yàn)V波方法有很好的去噪和信號(hào)細(xì)節(jié)保持能力。

表1 評(píng)價(jià)參數(shù)比較

3 實(shí)測(cè)信號(hào)分析

3.1 對(duì)工頻干擾的抑制

實(shí)測(cè)數(shù)據(jù)來(lái)自安徽,用EH-4儀器采集,采樣頻率是12kHz。圖2給出了實(shí)測(cè)大地電磁信號(hào)4個(gè)分量的時(shí)域波形。由于勘測(cè)現(xiàn)場(chǎng)周?chē)懈邏弘娋€(xiàn)經(jīng)過(guò),實(shí)測(cè)數(shù)據(jù)受到工頻干擾的影響,其時(shí)域波形顯示信號(hào)存在一個(gè)呈正弦波形態(tài)的干擾,且幅值很大。該50Hz的工頻干擾幾乎淹沒(méi)了電磁場(chǎng)信號(hào),從時(shí)域幾乎很難辨別出真實(shí)的電磁場(chǎng)信號(hào)。

考慮到文章的篇幅,這里僅給出電場(chǎng)信號(hào)Ex分量的詳盡處理過(guò)程。

圖3a右圖展示了去噪后的信號(hào)。圖3b右圖為去噪后信號(hào)的功率譜,可以看出,去噪后信號(hào)變得平滑,50Hz處的尖峰被消除。圖3c右圖給出了重新對(duì)去噪后的Ex分量進(jìn)行頻率切片小波變換得到的時(shí)頻分布圖。從圖3中可以看出,去噪后信號(hào)變得平滑,正弦波干擾被消除,被噪聲淹沒(méi)的有用信號(hào)得到了顯現(xiàn)。消噪后信號(hào)的統(tǒng)計(jì)參數(shù)如下:最大值為896.178mV,最小值為-1098.569mV,均值為-104.007mV,能量為1.706×102mV2。比較圖3 左、右兩邊的3幅圖,可以清楚地看到,50Hz的干擾得到了有效抑制,時(shí)頻譜中50Hz的能量帶也消失了,突出了有用信號(hào)的細(xì)節(jié)。用時(shí)頻域?yàn)V波方法對(duì)受噪的其它3個(gè)分量做同樣的處理,結(jié)果如圖4所示。從圖4可以看出,噪聲均得到了有效的抑制,信號(hào)變得平穩(wěn)。

圖2 受工頻干擾的大地電磁信號(hào)

圖3 大地電磁信號(hào)Ex分量去噪前(左)、后(右)的對(duì)比a 信號(hào)的時(shí)間序列; b 信號(hào)的功率譜; c 信號(hào)的時(shí)頻譜

分別用去噪前、后的MT數(shù)據(jù)計(jì)算了視電阻率曲線(xiàn)和相位曲線(xiàn)(為節(jié)約篇幅這里僅給出視電阻率曲線(xiàn)ρxy及其對(duì)應(yīng)的相位曲線(xiàn)Фxy,如圖5所示)。

圖4 壓制工頻干擾噪聲后的大地電磁信號(hào)

圖5 ρxy曲線(xiàn)(a)和Фxy曲線(xiàn)(b)去噪前、后的對(duì)比

去噪前、后響應(yīng)參數(shù)的平滑度用方差來(lái)衡量，它們的相似程度則用曲線(xiàn)相似度參數(shù)Ncc來(lái)評(píng)價(jià)。

曲線(xiàn)相似性參數(shù)(Ncc)定義[17,19]為:

(13)

式中:f(n),g(n)分別為兩離散序列;Ncc的取值為-1～1,其中,-1代表變換前、后兩數(shù)據(jù)波形反向,0代表兩波形正交,1代表完全相同。

去噪前、后評(píng)價(jià)參數(shù)的對(duì)比結(jié)果如表2所示。

從前文給出的圖和參數(shù)可以看出,這些曲線(xiàn)相似度Ncc都很大。這意味著去噪前、后得到的兩條視電阻率參數(shù)曲線(xiàn)形態(tài)趨勢(shì)一致。但去噪后曲線(xiàn)的方差減小了很多,曲線(xiàn)的奇異點(diǎn)得到了有效抑制(如在ρxy曲線(xiàn)的50Hz頻點(diǎn),視電阻率的值由2280Ω·m降到了286Ω·m)。對(duì)于所有的響應(yīng)參數(shù)曲線(xiàn),誤差棒的平均值減小了,更多的細(xì)節(jié)信息得到了顯現(xiàn),曲線(xiàn)變得平滑和連續(xù),估算的參數(shù)變得穩(wěn)定。

3.2 對(duì)三角波干擾的抑制

單個(gè)三角波干擾或由多個(gè)三角波構(gòu)成的干擾是大地電磁測(cè)深中常見(jiàn)的噪聲，尤其是在礦集區(qū)，井下大功率設(shè)備的開(kāi)停、開(kāi)采帶來(lái)的震動(dòng)干擾等都會(huì)產(chǎn)生此類(lèi)干擾。它們往往只出現(xiàn)在磁場(chǎng)信號(hào)中，且在時(shí)間域表示為明顯的跳變，其形態(tài)為不規(guī)則的三角波。圖6所示是典型的受三角波干擾的MT數(shù)據(jù)。此類(lèi)噪聲通常在磁道觀測(cè)信號(hào)中出現(xiàn),且為磁場(chǎng)相關(guān)噪聲,電道受其影響較小。噪聲局部能量較強(qiáng),將正常磁場(chǎng)信號(hào)完全湮滅。單個(gè)噪聲幅值為正常信號(hào)幅值的數(shù)倍到數(shù)十倍。

表2 去噪前、后評(píng)價(jià)參數(shù)的對(duì)比結(jié)果

這里僅給出磁場(chǎng)信號(hào)Hy分量的詳盡處理過(guò)程。圖7顯示了Hy分量的時(shí)域信號(hào)(圖7a)、FSWT變換后得到的時(shí)頻譜(圖7b)和功率譜(圖7c)。由時(shí)頻譜和功率譜可見(jiàn),含噪Hy分量的能量主要集中在低頻段(100Hz)以?xún)?nèi),在相同的坐標(biāo)軸刻度下,其它的信號(hào)能量幾乎完全被抑制,這不符合正常大地電磁信號(hào)能量分布規(guī)律。用本文提出的時(shí)頻域閾值濾波方法對(duì)圖7b 進(jìn)行濾波,濾波后再進(jìn)行時(shí)頻域頻率切片小波變換的逆變換,重構(gòu)信號(hào),結(jié)果如圖8所示。從圖8可見(jiàn),去噪后的信號(hào),三角波形得到了有效抑制,信號(hào)表現(xiàn)為以0為均值的正態(tài)分布,信號(hào)的能量分布也變得合理。用時(shí)頻域?yàn)V波的方法對(duì)含噪的另一個(gè)磁場(chǎng)分量Hx也做同樣的處理,處理后的MT信號(hào)如圖9 所示。從圖9可以看出,三角波干擾得到了有效的抑制,為后續(xù)的響應(yīng)參數(shù)估算奠定了基礎(chǔ)。

圖6 受三角波干擾的MT數(shù)據(jù)

圖7 Hy分量的時(shí)域信號(hào)(a)、FSWT變換后得到的時(shí)頻譜(b)和功率譜(c)

圖8 Hy分量去噪后的時(shí)域信號(hào)(a)、FSWT變換后得到的時(shí)頻譜(b)和功率譜(c)

圖9 壓制三角波干擾噪聲后的MT數(shù)據(jù)

4 結(jié)束語(yǔ)

大地電磁勘測(cè)中,野外觀測(cè)數(shù)據(jù)不可避免地會(huì)受到各種人文噪聲的影響,且隨著經(jīng)濟(jì)的發(fā)展而日趨嚴(yán)重。基于頻率切片小波變換的時(shí)頻分析能精確地呈現(xiàn)大地電磁信號(hào)的時(shí)頻能量分布特征,在此基礎(chǔ)上再對(duì)時(shí)頻譜進(jìn)行時(shí)頻閾值濾波和逆變換重構(gòu),能分離出去噪后的MT信號(hào),可以有效地濾除實(shí)測(cè)大地電磁信號(hào)中強(qiáng)噪聲,為后續(xù)穩(wěn)健的響應(yīng)參數(shù)估計(jì)奠定了基礎(chǔ)。應(yīng)用本文所提方法對(duì)仿真數(shù)據(jù)和工程實(shí)踐中常見(jiàn)的受工頻干擾和三角波干擾的實(shí)測(cè)MT信號(hào)進(jìn)行了分析,結(jié)果表明,這種方法可有效地消除強(qiáng)干擾;去噪后計(jì)算的響應(yīng)參數(shù)曲線(xiàn)變得平滑、連續(xù),提高了大地電磁測(cè)深數(shù)據(jù)的質(zhì)量。當(dāng)然,大地電磁噪聲的來(lái)源和成分復(fù)雜,噪聲的處理是一項(xiàng)重要而又艱巨的任務(wù),下一步將研究其它類(lèi)型噪聲的時(shí)頻特征及其處理方法,并討論有效的去噪評(píng)價(jià)體系。

[1] 李曉昌.大地電磁測(cè)深法在石油地質(zhì)調(diào)查中的應(yīng)用[J].物探化探計(jì)算技術(shù),2009,31(6):573-578 LI X C.Application of magnetotelluric sounding to petroleum investigation[J].Computing Techniques for Geophysical and Geochemical Exploration,2009,31(6):573-578

[2] 陳孝雄,王友勝,黃潛生.大地電磁測(cè)深在庫(kù)木庫(kù)里盆地結(jié)構(gòu)研究中的應(yīng)用[J].石油天然氣學(xué)報(bào),2007,29(6):78-82 CHEN X X,WANG Y S,HUANG Q S.Application of MT in structural study of kumukuli basin[J].Journal of Oil and Gas Technology,2007,29(6):78-82

[3] 王書(shū)明,王家映.高階統(tǒng)計(jì)量在大地電磁測(cè)深數(shù)據(jù)處理中的應(yīng)用研究[J].地球物理學(xué)報(bào),2004,47(4):928-934 WANG S M,WANG J Y.Application of higher-orderstatistics in magnetotelluric data processing[J].Chinese Journal of Geophysics,2004,47(4):928-934

[4] 嚴(yán)家斌,劉貴忠,柳建新.小波變換在天然電磁場(chǎng)信號(hào)時(shí)間序列處理中的應(yīng)用[J].地質(zhì)與勘探,2008,44(3):75-78 YAN J B,LIU G Z,LIU J X.Application of wavelet transform in processing nature electromagnetic fieldtime series[J].Geology and Prospecting,2008,44(3):75-78

[5] LARSEN J C,MACHIE R L.Robust smooth magnetotelluric transfer functions[J].Geophysical Journal International,1996,124(3):801-819

[6] EGBERT D,BOOKER R.Robust estimation of geomagnetic transfer function[J].Geophysical Journal of the Royal Astronomical Society,1986,87(1):173-194

[7] GAMBLE T M,GOUBAU W M,CLARKE J.Magnetotelluric with a remote magnetic reference[J].Geophysics,1979,44(1):53-68

[8] 蔡劍華,胡惟文,任政勇,等.基于高階統(tǒng)計(jì)量的大地電磁數(shù)據(jù)處理與仿真[J].中南大學(xué)學(xué)報(bào):自然科學(xué)版,2010,41(4):1556-1560 CAI J H,HU W W,Ren Z Y,et al.Magnetotelluric data processing and simulation based on higher-order statistics[J].Journal of Central South University(Science and Technology),2010,41(4):1556-1560

[9] 嚴(yán)家斌.大地電磁信號(hào)處理理論及方法研究[D].長(zhǎng)沙:中南大學(xué),2003 YAN J B.The study on theory and method of Mag-netotellurie signal proeessing[D].Changsha:Central South University,2003

[10] CAI J H.A combinatorial filtering method for magnetotelluric time-series based on Hilbert-Huang transform[J].Exploration Geophysics,2014,45(2):63-73

[11] CAI J H.Magnetotelluric response function estimation based on hilbert-huang transform[J].Pure and Applied Geophysics,2013,170(11):1899-1911

[12] 蔡劍華,王先春,胡惟文.基于經(jīng)驗(yàn)?zāi)B(tài)分解與小波閾值的MT信號(hào)去噪方法[J].石油地球物理勘探2013,48(2):303-307 CAI J H,WANG X C,HU W W.De-noising of MT signal based on empirical mode decomposition and wavelet threshold method[J].Oil Geophysical Prospecting,2013,48(2):303-307

[13] YAN Z H,MIYAMOTO A,JIANG Z.Frequency slice wavelet transform for transient vibration response analysis[J].Mechanical Systems and Signal Processing,2009,23(5):1474-1489

[14] YAN Z H,MIYAMOTO A,JIANG Z.An overall theoretical description of frequency slice wavelet transform[J].Mechanical Systems and Signal Processing,2010,24(2):491-507

[15] YAN Z H,MIYAMOTO A,JIANG Z.Frequency slice algorithm for modal signal separation and damping identification[J].Computers and Structures,2011,89(1/2):14-26

[16] 尚帥,韓立國(guó),胡瑋,等.壓縮小波變換地震譜分解方法應(yīng)用研究[J].石油物探,2015,54(1):51-55 SHANG S,HAN L G,HU W,et al.Applied research of synchrosqueezing wavelet transform in seismic spectral decomposition method[J].Geophysical Prospecting for Petroleum,2015,54(1):51-55

[17] DONOH D L.De-noising by soft-thresholding[J].IEEE Transactions on Information Theory,1995,41(3):613-627

[18] 段晨東,高強(qiáng).基于時(shí)頻切片分析的故障診斷方法及應(yīng)用[J].振動(dòng)與沖擊,2011,30(9):1-5 DUAN C D,GAO Q.Noval fault diagnosis approach using time-frequency slice analysis and its application[J].Journal of Vibration and Shock,2011,30(9):1-5

[19] 賈海青,姜弢,徐學(xué)純,等.基于全方向波束定向的地震記錄信噪比改善方法[J].石油物探,2015,54(5):521-530 JIA H Q,JIANG T,XU X C,et al.The SNR improvement method based on omnidirectional beam-forming for vibrator seismic record[J].Geophysical Prospecting for Petroleum,2015,54(5):521-530

(編輯：顧石慶)

Magnetotelluric data denosing based on time-frequency analysis of the frequency slice wavelet transform

CAI Jianhua1,2,XIONG Rui1,2

(1.DepartmentofPhysicsandElectronics,HunanUniversityofArtsandScience,Changde415000,China;2.CooperativeInnovationCenterforTheConstruction&DevelopmentofDongtingLakeEcologicalEconomicZone,HunanUniversityofArtsandScience,Changde415000,China)

Magnetotelluric (MT) survey methods always affected by a variety of artificial noises,which is more obvious in some frequency bands.In order to solve this problem,an improved denosing method for magnetotelluric data in time-frequency domain was proposed based on time-frequency analysis of the frequency slice wavelet transform (FSWT).The principle and implementation steps of the method were presented.Firstly,the MT signal was processed with the FSWT to get its time-frequency distribution for full frequency bands.Then the derived time-frequency spectrum was applied with the threshold filtering and the denoised signal was reconstructed by the inverse FSWT of the processed time-frequency distribution.The simulation experiment and the measured signal analysis show that the proposed method can effectively eliminate the strong noise of MT signal.After denoising,the mutation points of the calculated response parameters curve were suppressed and parameter curves become smooth and continuous.

frequency slice wavelet transform,time-frequency analysis,magnetotelluric (MT) signal,denosing

2015-11-16;改回日期：2016-04-18。

蔡劍華(1979—),男,博士,副教授,主要從事大地電磁數(shù)據(jù)處理的研究。

國(guó)家自然科學(xué)基金項(xiàng)目(41304098)和湖南省教育廳重點(diǎn)項(xiàng)目(16A146)共同資助。

This research is financially supported by the National Natural Science Foundation of China (Grant No.41304098) and the Key Research Fund of Hunan Provincial Education Department (Grant No.16A146).

P631

A

1000-1441(2016)06-0904-09

10.3969/j.issn.1000-1441.2016.06.016