陳越 劉雄英 吳中堂 范藝 任子良 馮久超

1)(華南理工大學(xué)電子與信息學(xué)院,廣州 510641)

2)(廣東技術(shù)師范學(xué)院電子與信息學(xué)院,廣州 510665)

受污染混沌信號(hào)的協(xié)同濾波降噪?

陳越1)劉雄英1)?吳中堂1)范藝2)任子良1)馮久超1)

1)(華南理工大學(xué)電子與信息學(xué)院,廣州 510641)

2)(廣東技術(shù)師范學(xué)院電子與信息學(xué)院,廣州 510665)

(2017年5月19日收到;2017年6月18日收到修改稿)

混沌信號(hào),協(xié)同濾波,噪聲抑制

1 引 言

混沌行為常存在于生物、電氣、機(jī)械等物理系統(tǒng)[1],從觀測(cè)數(shù)據(jù)中檢測(cè)混沌是認(rèn)識(shí)、分析和預(yù)測(cè)這些系統(tǒng)的關(guān)鍵.通常觀測(cè)信號(hào)是被噪聲污染的混沌行為,掩蓋了系統(tǒng)真實(shí)的動(dòng)力學(xué)特性,有效抑制噪聲是分析和研究這些系統(tǒng)的前提.混沌信號(hào)固有的非周期、寬頻譜等特性使得傳統(tǒng)的線性和頻域?yàn)V波方法無(wú)法用于這類(lèi)信號(hào)[2],因此,針對(duì)混沌信號(hào)的非線性動(dòng)力學(xué)行為研究相應(yīng)的降噪方法具有重要意義.

目前,受污染混沌信號(hào)的降噪問(wèn)題正在引起學(xué)術(shù)界的關(guān)注,并提出了一些降噪方法[3?15]:梯度下降法對(duì)混沌信號(hào)降噪能獲得一條比原時(shí)間序列噪聲更小的時(shí)間序列軌道[3],但不能完全重構(gòu)混沌信號(hào);局部投影方法通過(guò)在參考點(diǎn)的鄰域近似吸引子的局部動(dòng)力學(xué)行為來(lái)降噪[4,5],噪聲較強(qiáng)時(shí)由于鄰域半徑擴(kuò)大導(dǎo)致性能下降嚴(yán)重;自然界中多數(shù)信號(hào)的小波變換是稀疏的,即只有少量變換系數(shù)顯著大于零[16],而高斯噪聲的能量則均勻分布在全部變換系數(shù)中,小波閾值法通過(guò)將較小的變換系數(shù)置零實(shí)現(xiàn)噪聲抑制,然后用逆變換重構(gòu)混沌信號(hào)[6,7],這類(lèi)方法的降噪效果在很大程度上取決于變換系數(shù)的稀疏程度,因此,必須選取合適的小波基;經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)是一種數(shù)據(jù)驅(qū)動(dòng)的分解方法,EMD閾值降噪[8]克服了小波閾值必須針對(duì)信號(hào)特征選取合適小波基的問(wèn)題,并在混沌信號(hào)降噪中逐漸得到應(yīng)用[9?11],但這類(lèi)方法仍存在閾值、篩分迭代次數(shù)難以確定的問(wèn)題;文獻(xiàn)[12]使用壓縮感知領(lǐng)域的基追蹤降噪技術(shù)來(lái)處理混沌信號(hào),但為混沌信號(hào)構(gòu)造“字典”是一個(gè)難題,從實(shí)驗(yàn)結(jié)果看降噪效果也不理想;局部曲線擬合法[13]首先對(duì)數(shù)據(jù)進(jìn)行分段,再利用多項(xiàng)式對(duì)分段數(shù)據(jù)做局部近似,最后通過(guò)加權(quán)平滑來(lái)重構(gòu)混沌信號(hào),文獻(xiàn)[14]通過(guò)自動(dòng)搜索最佳擬合窗長(zhǎng)進(jìn)一步完善了這一方法,但由于混沌的高度非線性,局部線性近似難以實(shí)現(xiàn)精確的重構(gòu).已有算法本質(zhì)上都是傳統(tǒng)非線性降噪方法在混沌信號(hào)上的應(yīng)用,而不是針對(duì)混沌的某些特征專(zhuān)門(mén)設(shè)計(jì)的新方法,因此,它們?cè)趯?duì)混沌信號(hào)降噪時(shí)難以確保獲得好的效果.

自然圖像中包含大量的相似塊,如局部紋理、相似邊沿等,基于這一特性,Dabov等[17]提出了協(xié)同濾波的概念和著名的三維塊匹配(block matching 3D,BM3D)算法,該算法被認(rèn)為是目前性能最好的彩色圖像降噪算法[18].在BM3D中,相似的二維(2D)圖像塊堆疊在一起形成3D分組,利用分組內(nèi)圖像塊間的相似性對(duì)3D分組聯(lián)合濾波的過(guò)程被稱(chēng)為協(xié)同濾波.

分形是混沌運(yùn)動(dòng)的基本特征之一,它在相空間體現(xiàn)為軌道的自相似結(jié)構(gòu)[19],這意味著混沌信號(hào)在時(shí)域上也存在一定程度的自相似性.文獻(xiàn)[20]在Mackey-Glass,Lorenz等信號(hào)的時(shí)域波形中觀察到明顯的自相似特征,并將之歸因?yàn)榛煦缥拥淖韵嗨平Y(jié)構(gòu)和自仿射分形特性;文獻(xiàn)[21]發(fā)現(xiàn)中心頻率附近出現(xiàn)的連續(xù)頻段導(dǎo)致了Duffing信號(hào)波形的自相似,并認(rèn)為混沌信號(hào)具有頻率調(diào)制的相似性.本文利用混沌信號(hào)的這種自相似特性,將協(xié)同濾波的方法用于受污染混沌信號(hào)的降噪.在設(shè)計(jì)降噪算法的過(guò)程中,首先在觀測(cè)信號(hào)中搜索相似的片段并建立分組,然后在二維變換域用閾值法對(duì)每個(gè)分組做協(xié)同濾波,再經(jīng)反變換得到分組中各片段的濾波結(jié)果,最后根據(jù)所有分組的濾波結(jié)果重構(gòu)原信號(hào).由于同一分組內(nèi)的片段具有良好的相關(guān)性,與直接對(duì)觀測(cè)信號(hào)做一維變換相比,分組的二維變換能更稀疏的表示原信號(hào).

2 混沌信號(hào)協(xié)同濾波降噪算法

混沌的分形特性表現(xiàn)為混沌信號(hào)在相空間具有自相似結(jié)構(gòu),這意味著其時(shí)域波形也存在一定的自相似性.圖1是一段被噪聲污染的Lorenz混沌信號(hào),可以看到,信號(hào)波形具有自相似性,圖中用不同線型框標(biāo)示了波形中的部分相似片段.

圖1 (網(wǎng)刊彩色)Lorenz信號(hào)的自相似性Fig.1.(color online)Self-similarity property of Lorenz signal.

利用混沌信號(hào)的這種自相似性抑制噪聲的流程如圖2所示,它包括分組、協(xié)同濾波和重構(gòu)三個(gè)步驟:首先,從觀測(cè)信號(hào)中建立若干相似分組;然后,對(duì)每個(gè)分組進(jìn)行協(xié)同濾波;最后,利用各分組的濾波結(jié)果重構(gòu)原信號(hào).

圖2 (網(wǎng)刊彩色)降噪算法流程示意圖 (a)分組;(b)協(xié)同濾波;(c)重構(gòu)Fig.2.(color online)Process of the proposed denoising algorithm:(a)Grouping;(b)collaborative filtering;(c)reconstruction.

2.1 分 組

記參考片段為R,任意片段S與R的相似程度可以用二者的歸一化距離來(lái)度量[17]:

其中,w是片段寬度,該距離越小,S與R的相似程度越高.

以參考片段R為中心在長(zhǎng)度為l(l?w)的窗口內(nèi)搜索,找到與R距離最小的m個(gè)片段形成分組group(R),并保存到m行w列的二維數(shù)組中,即group(R)∈Rm×w.讓參考片段以步長(zhǎng)δ從觀測(cè)信號(hào)的起始位置開(kāi)始向后移動(dòng),搜索窗口也隨之移動(dòng).記錄在每個(gè)位置得到的相似分組,并標(biāo)記分組中每個(gè)片段的位置.

2.2 協(xié)同濾波

協(xié)同濾波以分組為單位進(jìn)行,分組內(nèi)任一片段的濾波將用到同組其他片段的信息,因此可以看作是一個(gè)相互“協(xié)作”的過(guò)程[17].協(xié)同濾波包括以下步驟.

1)對(duì)分組做二維變換:

其中,G=[gi,j]∈Rm×w(i=1,···,m;j=1,···,w) 是變換系數(shù)矩陣,T[·]是二維離散余弦變換(DCT)或二維小波等任意等距線性變換.本文后續(xù)將采用DCT變換.

2)為每個(gè)分組選定一個(gè)閾值λ(R)并在變換域做閾值降噪,將小于閾值的系數(shù)置零來(lái)衰減噪聲,保留大于閾值的系數(shù).閾值函數(shù)HT(·)由下式定義:

閾值λ(R)使用Donoho和Johnstone提出的Visu-Shrink方法[16]確定,該方法在均方誤差意義上是近似最優(yōu)的,同時(shí)重構(gòu)信號(hào)能保持原信號(hào)的光滑性[22],特別適合于處理混沌流信號(hào).VisuShrink閾值定義為

其中,σ是觀測(cè)信號(hào)的噪聲標(biāo)準(zhǔn)差,當(dāng)噪聲強(qiáng)度未知時(shí)它可以由系數(shù)矩陣G(R)的中位數(shù)絕對(duì)偏差來(lái)估計(jì)[16]:

3)對(duì)變換系數(shù)G(R)做閾值后,通過(guò)逆變換得到分組的濾波結(jié)果:

它包含了參考片段R對(duì)應(yīng)分組的每個(gè)片段的降噪估計(jì).

2.3 信號(hào)的重構(gòu)

由2.1節(jié)的描述可知,片段間存在重疊.一般來(lái)說(shuō)一個(gè)信號(hào)點(diǎn)會(huì)同時(shí)屬于多個(gè)不同片段,這些片段在該點(diǎn)的濾波結(jié)果也會(huì)不同.最終的降噪結(jié)果可以逐點(diǎn)計(jì)算,每個(gè)重構(gòu)信號(hào)點(diǎn)由所有包含該點(diǎn)的片段在此處的濾波輸出做算術(shù)平均得到:

3 算法參數(shù)分析

本節(jié)以疊加高斯白噪聲的Lorenz信號(hào)為例討論算法參數(shù)對(duì)性能的影響.Lorenz系統(tǒng)的方程為

其中,參數(shù)α=10,β=28,γ=8/3.方程通過(guò)四階龍格-庫(kù)塔法求解,步長(zhǎng)為0.01.在每次實(shí)驗(yàn)中,方程從混沌區(qū)的隨機(jī)位置開(kāi)始迭代,然后由狀態(tài)變量x產(chǎn)生混沌信號(hào).算法性能通過(guò)輸出信噪比(SNRout)衡量,在相同的輸入信噪比(SNRin)下,SNRout越高表明降噪效果越好.SNR定義為

其中,s(n)為純凈的混沌信號(hào),在計(jì)算SNRin時(shí)x(n)是含噪觀測(cè)信號(hào),在計(jì)算SNRout時(shí)x(n)是降噪算法的輸出.

協(xié)同濾波降噪算法需要指定的參數(shù)有:片段寬度w、搜索窗長(zhǎng)l、搜索窗移動(dòng)步長(zhǎng)δ和分組片段數(shù)m.分組片段數(shù)m對(duì)性能的影響非常有限,在后面的討論中將它固定為m=30.

3.1 片段寬度w

取搜索窗長(zhǎng)l=4000,搜索窗移動(dòng)步長(zhǎng)δ=30,圖3顯示了片段寬度w與降噪效果的關(guān)系.一般,樣點(diǎn)數(shù)量不足會(huì)導(dǎo)致降噪算法的性能下降,當(dāng)w較小時(shí)SNRout偏低,因此,w的取值不應(yīng)小于100.

圖3 不同SNRin下SNRout與w的關(guān)系Fig.3.SNRoutversus w at different SNRin.

3.2 搜索窗長(zhǎng)l

窗長(zhǎng)l越長(zhǎng),分組越能在較大范圍內(nèi)搜集相似片段,算法性能也就越好,圖4(w=120,δ=30)反映了l與降噪性能的這一關(guān)系.考慮到l增大會(huì)增加算法的復(fù)雜度,并且l較大后性能提升并不明顯,折中后選取l為4000左右.如果要獲得最佳性能,可以將整個(gè)觀測(cè)信號(hào)作為搜索范圍.

圖4 不同SNRin下SNRout與l的關(guān)系Fig.4.SNRoutversus l at different SNRin.

3.3 搜索窗移動(dòng)步長(zhǎng)δ

圖5 (w=120,l=4000)顯示了搜索窗移動(dòng)步長(zhǎng)δ與SNRout的關(guān)系.步長(zhǎng)δ越小建立的相似分組就越多,對(duì)觀測(cè)信號(hào)的分析也越精細(xì),因此,SNRout隨步長(zhǎng)δ減小逐漸提高.由于δ減小會(huì)顯著增加算法的復(fù)雜度,一般折中選取δ為片段寬度w的1/2到1/4左右.

圖5 不同SNRin下SNRout與δ的關(guān)系Fig.5.SNRoutversus δ at different SNRin.

4 實(shí)驗(yàn)結(jié)果

本節(jié)通過(guò)仿真實(shí)驗(yàn)評(píng)估協(xié)同濾波降噪算法的性能,算法參數(shù)為w=120,l=4000,δ=w/4,觀測(cè)信號(hào)由Lorenz信號(hào)疊加高斯白噪聲生成.為了顯示算法的優(yōu)越性,還與以下算法進(jìn)行了比較:參考文獻(xiàn)[14]的方法,自動(dòng)搜索最佳擬合窗長(zhǎng)的局部曲線擬合算法[13];小波閾值降噪算法[6],用db8小波作基函數(shù),分解4層,采用軟閾值;EMD迭代區(qū)間閾值(EMD iterative interval thresholding,EMD-IIT)法[8],篩分9次,迭代15次.算法性能通過(guò)SNRout、重構(gòu)的均方誤差(root mean square error,RMSE)和對(duì)相圖的還原情況進(jìn)行比較.RMSE定義為

其中,s(n)是純凈的混沌信號(hào),x(n)是降噪算法的輸出.

取信號(hào)長(zhǎng)度為10000,圖6顯示了不同噪聲強(qiáng)度下各種算法降噪后的SNRout,圖7顯示了重構(gòu)信號(hào)與純凈混沌信號(hào)的均方誤差,結(jié)果由100次隨機(jī)實(shí)驗(yàn)平均得到.可以看到,局部曲線擬合法的性能略好于小波閾值,EMD-IIT法在低信噪比下性能較差.本文算法在實(shí)驗(yàn)中表現(xiàn)出較好的噪聲抑制能力,其SNRout比另外三種方法高2.5—4 dB,從RMSE來(lái)看本文算法的重構(gòu)精度也明顯優(yōu)于其他方法.

圖6 (網(wǎng)刊彩色)不同算法的SNRout比較Fig.6.(color online)SNRoutcomparison among different algorithms.

圖7 (網(wǎng)刊彩色)不同算法的均方誤差比較Fig.7.(color online)RMSE comparison among different algorithms.

對(duì)SNRin=10 dB的Lorenz信號(hào)進(jìn)行降噪實(shí)驗(yàn),幾種算法重構(gòu)的相圖如圖8所示.可以看到,觀測(cè)信號(hào)的相空間被噪聲湮沒(méi),幾乎觀察不到分形特征.如圖8(c)—(e)所示,局部曲線擬合、小波閾值和EMD-IIT算法都能大致還原純凈吸引子的相圖,但由于軌道扭曲變形,使得原本有序的自相似結(jié)構(gòu)顯得雜亂無(wú)章.從圖8(f)可見(jiàn),本文算法重構(gòu)的相圖具有更光滑的軌道,吸引子在相空間完整的分形結(jié)構(gòu)也被更好地保留了下來(lái).

5 結(jié) 論

結(jié)合混沌信號(hào)的自相似分形特性,本文提出了一種基于協(xié)同濾波的降噪算法.該算法的創(chuàng)新之處在于通過(guò)對(duì)相似片段分組將一維信號(hào)的降噪轉(zhuǎn)化為一個(gè)二維聯(lián)合濾波問(wèn)題.由于相似片段間具有良好的相關(guān)性,與直接在一維變換域做閾值降噪相比,對(duì)分組的二維變換能獲得原信號(hào)更稀疏的表示,從而更好地抑制噪聲.仿真結(jié)果表明:本文算法在信噪比提升和重構(gòu)精度上明顯優(yōu)于小波閾值、局部曲線擬合等現(xiàn)有的混沌信號(hào)降噪算法;從對(duì)相圖的還原情況來(lái)看,本文算法重構(gòu)的吸引子具有更清晰的自相似結(jié)構(gòu)和更光滑的軌道,對(duì)純凈混沌吸引子的還原也更加準(zhǔn)確.所以,協(xié)同濾波結(jié)合混沌自相似特性對(duì)受污染混沌信號(hào)進(jìn)行降噪,在重構(gòu)原始混沌信號(hào)方面具有較大的應(yīng)用潛力.

圖8 Lorenz吸引子相圖 (a)無(wú)噪聲污染;(b)噪聲污染;(c)局部曲線擬合降噪;(d)小波閾值降噪;(e)EMD-IIT降噪;(f)本文算法降噪Fig.8.Phase portraits of Lorenz attractor reconstructed from(a)clean signal,(b)noisy signal,(c)signal denoised by local curve fitting,(d)signal denoised by wavelet thresholding,(e)signal denoised by EMD-IIT,(f)signal denoised by the proposed method.

[1]Feng J C 2012Chaotic Signals and Information Processing(Beijing:Tsinghua University Press)pp32–35(in Chinese)[馮久超 2012混沌信號(hào)與信息處理 (北京:清華大學(xué)出版社)第32—35頁(yè)]

[2]Badii R,Broggi G,Derighetti B,Ravani M,Ciliberto S,Politi A,Rubio M A 1988Phys.Rev.Lett.60 979

[3]Liu X Y,Qiu S S,Lau C M 2005J.Syst.Eng.Electron.16 253

[4]Cawley R,Hsu G H 1992Phys.Rev.A46 3057

[5]Leontitsis A,Bountis T,Pange J 2004Chaos14 106

[6]Han M,Liu Y H,Xi J H,Guo W 2007IEEE Signal Proc.Lett.14 62

[7]Constantine W L B,Reinhall P G 2001Int.J.Bifurcat.Chaos11 483

[8]Kopsinis Y,McLaughlin S 2009IEEE Trans.Signal Proc.57 1351

[9]Wang X F,Qu J L,Gao F,Zhou Y P,Zhang Y X 2014Acta Phys.Sin.63 170203(in Chinese)[王小飛,曲建嶺,高峰,周玉平,張翔宇2014物理學(xué)報(bào)63 170203]

[10]Wei X L,Lin R L,Liu S Y,Zhang C H 2016Shock Vib.2016 1

[11]Wang W B,Zhang X D,Wang X L 2013Acta Phys.Sin.62 050201(in Chinese)[王文波,張曉東,汪祥莉 2013物理學(xué)報(bào)62 050201]

[12]Li G M,Lü S X 2015Acta Phys.Sin.64 160502(in Chinese)[李廣明,呂善翔 2015物理學(xué)報(bào) 64 160502]

[13]Tung W W,Gao J B,Hu J,Yang L 2011Phys.Rev.E83 046210

[14]Wang M J,Wu Z T,Feng J C 2015Acta Phys.Sin.64 040503(in Chinese)[王夢(mèng)蛟,吳中堂,馮久超2015物理學(xué)報(bào)64 040503]

[15]Hu J F,Zhang Y X,Yang M,Li H Y,Xia W,Li J 2016Nonlinear Dynam.84 1469

[16]Donoho D L,Johnstone I M 1994Biometrika81 425

[17]Dabov K,Foi A,Katkovnik V,Egiazarian K 2007IEEE Trans.Image Proc.16 2080

[18]Lebrun M 2012Image Proc.On Line2 175

[19]Yu S M 2011Chaotic Systems and Chaotic Circuits(Xi’an:Xidian University Press)pp10–12(in Chinese)[禹思敏2011混沌系統(tǒng)與混沌電路(西安:西安電子科技大學(xué)出版社)第10—12頁(yè)]

[20]He T,Zhou Z O 2007Acta Phys.Sin.56 693(in Chinese)[賀濤,周正歐 2007物理學(xué)報(bào) 56 693]

[21]Tang Y F,Liu S L,Lei N,Jiang R H,Liu Y H 2012Acta Phys.Sin.61 170504(in Chinese)[唐友福,劉樹(shù)林,雷娜,姜銳紅,劉穎慧2012物理學(xué)報(bào)61 170504]

[22]Coifman R R,Donoho D L 1995Lect.Notes Stat.103 125

PACS:05.45.—a,05.40.Ca DOI:10.7498/aps.66.210501

*Project supported by the National Natural Science Foundation of China(Grant No.61372008)and the Science and Technology Planning Project of Guangdong Province,China(Grant Nos.2015B010101006,2014A010103014).

?Corresponding author.E-mail:liuxy@scut.edu.cn

Denoising of contaminated chaotic signals based on collaborative filtering?

Chen Yue1)Liu Xiong-Ying1)?Wu Zhong-Tang1)Fan Yi2)Ren Zi-Liang1)Feng Jiu-Chao1)

1)(School of Electronic and Information Engineering,South China University of Technology,Guangzhou 510641,China)
2)(School of Electronic and Information Engineering,Guangdong Polytechnic Normal University,Guangzhou 510665,China)

d 19 May 2017;revised manuscript

18 June 2017)

Reconstructing chaotic signals from noised data plays a critical role in many areas of science and engineering.However,the inherent features,such as aperiodic property,wide band spectrum,and extreme sensitivity to initial values,present a big challenge of reducing the noises in the contaminated chaotic signals.To address the above issues,a novel noise reduction algorithm based on the collaborative filtering is investigated in this paper.By exploiting the fractal self-similarity nature of chaotic attractors,the contaminated chaotic signals are reconstructed subsequently in three steps,i.e.,grouping,collaborative filtering,and signal reconstruction.Firstly,the fragments of the noised signal are collected and sorted into different groups by mutual similarity.Secondly,each group is tackled with a hard threshold in the two-dimensional(2D)transforming domain to attenuate the noise.Lastly,an inverse transformation is adopted to estimate the noise-free fragments.As the fragments with in a group are closely correlated due to their mutual similarity,the 2D transform of the group should be sparser than the one-dimensional transform of the original signal in the first step,leading to much more effective noise attenuation.The fragments collected in the grouping step may overlap each other,meaning that a signal point could be included in more than one fragment and have different collaborative filtering results.Therefore,the noise-free signal is reconstructed by averaging these collaborative filtering results point by point.The parameters of the proposed algorithm are discussed and a set of recommended parameters is given.In the simulation,the chaotic signal is generated by the Lorenz system and contaminated by addictive white Gaussian noise.The signalto-noise ratio and the root mean square error are introduced to measure the noise reduction performance.As shown in the simulation results,the proposed algorithm has advantages over the existing chaotic signal denoising methods,such as local curve fitting,wavelet thresholding,and empirical mode decomposition iterative interval thresholding methods,in the reconstruction accuracy,improvement of the signal-to-noise ratio,and recovering quality of the phase portraits.

chaotic signal,collaborative filtering,noise reduction

根據(jù)混沌吸引子的自相似分形特性,提出了一種利用協(xié)同濾波重構(gòu)受污染混沌信號(hào)的降噪算法.所設(shè)計(jì)的降噪算法通過(guò)對(duì)相似片段的分組將一維混沌信號(hào)的降噪轉(zhuǎn)化為一個(gè)二維聯(lián)合濾波問(wèn)題;然后,在二維變換域用閾值法衰減噪聲;最后,通過(guò)反變換獲得原始信號(hào)的估計(jì).由于分組中的相似片段具有良好的相關(guān)性,與直接在一維變換域做閾值降噪相比,分組的二維變換能獲得原信號(hào)更稀疏的表示,更好地抑制噪聲.仿真結(jié)果表明,該算法對(duì)原始混沌信號(hào)的重構(gòu)精度和信噪比的提升都優(yōu)于小波閾值、局部曲線擬合等現(xiàn)有的混沌信號(hào)降噪方法,對(duì)相圖的還原質(zhì)量也更好.

10.7498/aps.66.210501

?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):61372008)和廣東省科技計(jì)劃項(xiàng)目(批準(zhǔn)號(hào):2015B010101006,2014A010103014)資助的課題.

?通信作者.E-mail:liuxy@scut.edu.cn

?2017中國(guó)物理學(xué)會(huì)Chinese Physical Society