基于迭代重建算法的X射線光柵相位CT成像?

戚俊成 陳榮昌 劉賓 陳平 杜國浩 肖體喬

1)(中北大學信息與通信工程學院,太原 030051)

2)(中國科學院上海應用物理研究所,上海 201800)

基于迭代重建算法的X射線光柵相位CT成像?

戚俊成1)?陳榮昌2)劉賓1)陳平1)杜國浩2)肖體喬2)?

1)(中北大學信息與通信工程學院,太原 030051)

2)(中國科學院上海應用物理研究所,上海 201800)

(2016年8月18日收到;2016年12月5日收到修改稿)

基于光柵干涉儀的X射線成像技術可以同時獲得樣品內部的吸收信息、相位信息和散射信息,既保持了傳統X射線衰減成像的優點,又擁有相襯成像和散射成像的優勢.然而基于傳統CT重建算法的X射線光柵成像需要采集大量完整的原始投影數據,數據采集時間過長從而使得物體接受很大的輻射劑量,難以在實際中應用.提出基于傳統代數迭代重建算法的光柵成像技術.該方法利用現有X射線光柵成像系統采集少量原始投影數據,基于傳統代數迭代重建算法,對旋轉變化的相位數據進行CT重構,同時基于傅里葉變換的方法對微分相位數據進行相位恢復.模擬和實驗結果表明,基于少量或不完備的原始投影數據,該方法能夠準確重構成像對象的吸收、相位和散射三維信息,同時還能對微分相位切片進行高信噪比的相位恢復,得到樣品折射率實部衰減率,為X射線光柵成像技術在工業、生物和醫學診斷等領域的應用提供理論和技術支撐.

X射線光柵成像,迭代重建算法,相位恢復

1 引 言

自從X射線發現至今,X射線成像一直受到廣泛的關注.X射線光柵成像技術對X射線源的相干性要求較低,并可以同時獲得樣品的吸收、散射和相位信息[1].更重要的是,基于光柵的X射線相襯成像方法不受視場限制,可以實現大視場成像.因此,該方法被視為最有潛力從實驗室走向實際應用的X射線相襯成像方法,受到了廣泛的的關注,也具有廣闊的實際應用前景.

傳統光柵成像技術的重建算法是以希爾伯特變換為核心的反投影算法,是一種解析重建算法,要求具有完整的投影數據.在數據采集過程中,以步進掃描的方式采集數據需要花費大量時間采集大量投影作為原始數據,使得樣品接收大量的輻射.針對這個問題,國內外科學家采取了很多辦法來減少投影數,從而降低樣品所受輻射劑量.Momose等[2]和Wen等[3]采用一次曝光的方法采集數據,即樣品每旋轉一個角度,采集一張疊柵條紋圖像,然后將疊柵條紋圖像進行傅里葉變換,對比有樣品和沒有樣品情況下經傅里葉變換后的各個參數,就可以得到樣品的各種信息.這種方法的優點是不需要相位步進掃描,時間分辨率較高,但是成像的空間分辨率受到疊柵條紋周期限制,并且對光柵和探測器的要求也很高,要求相位光柵和吸收光柵完全匹配.Zhu等[4]采用反投影的方法,在樣品0?—360?旋轉的同時采集數據,然后對比樣品旋轉角度差為180?時的兩幅圖像,從而得到樣品信息.這種方法在CT過程中省去了步進掃描,使得實驗效率大大提升,但是無法獲得散射信息.Jensen等[5]采用交錯光柵步進的方法采集數據,在樣品旋轉的同時進行光柵掃描.Zanette等[6,7]提出了交錯相位步進掃描法,即在樣品旋轉的同時進行相位步進掃描,從多個樣品旋轉角度中獲取樣品的信息.在同等輻射劑量下,此方法可以獲得較好的結果,但是也存在投影角度不精確的問題,當投影數較少時,遠離旋轉中心的樣品高頻信息會丟失,而且此方法較難實現.

此外,光柵成像在實際應用中經常無法獲得完整的投影數據,主要原因是輻射劑量、經濟成本、無法檢測或很難檢測到完全的投影數據等.對于很大的樣品,要求探測器和光柵面積足夠大才能重建,難以實現;有時投影不是均勻地分布在180?或360?范圍內,這時濾波反投影(filtered back projection,FBP)方法受到了限制.

綜上所述,本文針對X射線光柵成像技術存在的問題,通過傳統代數迭代重建算法(algebraic reconstruction technique,ART),利用少量或者非完整的投影數據,對快速、低劑量光柵成像技術的CT重建和相位恢復算法展開了研究.利用傳統ART對由光柵成像所獲得的吸收、散射和旋轉變化的一階微分相位信息進行了高精度重建.對一階微分相位切片進行了高信噪比的相位恢復,為X射線光柵成像技術在生物醫學、材料和工業等領域的研究應用提供技術支持.

2 原 理

2.1 光柵成像技術

基于準平行光束的X射線光柵成像系統主要由相位光柵G1、吸收光柵G2和探測器組成(圖1)[8].當X射線光束穿過相位光柵G1后,由于泰伯效應(光柵自成像效應)[9],在G1后的泰伯距離處,將會得到光柵G1的自成像.在此像面上放置一個周期與自成像條紋對應的吸收光柵G2,即可得到疊柵干涉條紋.當X射線透過放置于相位光柵前的樣品后,X射線被物體吸收、折射和散射,相位光柵的入射X射線波前發生變化,從而導致相位光柵的泰伯像條紋發生強度變化和形狀扭曲.光柵成像系統的作用之一就是定量測出此泰伯像條紋的扭曲度和光強度的變化,最終得到X射線穿過樣品的吸收、相位和散射信息.

為了從探測器所得到的信息中分離出相位信息和其他信息,采用相位步進技術使一塊光柵以步進的方式掃描另外一塊光柵,獲得每一個像素點上隨著步進數變化的光強位移曲線[10].在忽略高次項的情況下,位移曲線方程可以傅里葉變換展開為[1,11,12]

式中(m,n)為像素點的坐標,k=2π/d,d為光柵G2的周期,yg為相位步進掃描過程中相位光柵所處的位置,ai為振幅系數,?i為相對應的相位系數.

圖1 X射線光柵成像系統基本布局圖Fig.1.Layout for the X-ray grating imaging system.

樣品的吸收信息與每一像素點位移曲線的0階振幅系數a0有關;微分相位信息由位移曲線的相位系數?1表征[11];散射信息主要反映了樣品的局部散射能力,樣品的散射信息可由條紋的可見度V=a1/a0表征[1].

2.2 光柵成像的物理模型分析

X射線光柵成像技術的一大優點就是通過一組原始數據可以同時獲得樣品的吸收、相位和散射信息,而這三種信息的表征形式是不同的.吸收信息和散射信息的表征形式是樣品的吸收系數或者散射因子的線積分,而在整個CT掃描過程中,吸收系數的整體分布隨著樣品一起轉動,即吸收系數都是旋轉不變的;對于大部分成像樣品來講,也可以認為散射因子是旋轉不變的或者近似旋轉不變的.所以,對于吸收信息和旋轉不變的散射信息,只需要通過實驗獲得相應的CT投影數據,利用傳統迭代重建算法就可以對其進行CT重建.

相位信息的數據采集過程如圖2所示,o′x′y′z′為樣品所在坐標系,隨著樣品一起轉動,oxyz是光柵和X射線光源所在坐標系,是靜止不動的.在CT數據采集過程中,樣品折射率實部衰減率δ的整體分布隨著樣品一起轉動,即δ的分布是旋轉不變的.但是光柵所在坐標系相對于樣品坐標系不是靜止的,所以δ沿光柵坐標系y方向的一階偏導的分布是旋轉變化的.因此,微分相位信息不能像吸收信息或者散射信息一樣,直接應用傳統迭代CT重建算法進行CT重構.

圖2 相位信息數據采集過程示意圖Fig.2.Schematic of phase acquisition.

通過對X射線在樣品中傳播過程的分析可知,樣品的相位信息取決于X射線穿過樣品后的折射角,而折射角為?δ/?y沿著X射線傳播方向的積分.在CT數據采集過程中,隨著投影角度?的變化,投影數據、投影角度和δ一階偏導數之間的關系為[13]

代表投影角度為?時恢復得到的相位導數和分別代表初始δ分布沿著垂直于X射線的傳播方向和平行于X射線的傳播方向上的一階偏導數.通過(2)式可以將旋轉變化的的投影信息轉變為旋轉不變的的投影信息,即將旋轉變化的微分相位信息轉變成旋轉不變的微分相位信息.和

2.3 代數迭代CT重建

對于濾波反投影CT重建算法而言,需要完整并且角度分布均勻的投影數據,而迭代重建適合投影數據不完全或投影間隔不均勻等場合的圖像重建[14].傳統代數迭代CT重建方法的基本思想是:首先將圖像離散地分解為N=n×m的矩形區域,每一個區域代表一個像素,同時給定一個初始圖像估計值f(0),然后計算f(0)的投影數據,并用該投影數據與實際投影數據的差來校正圖像,得到一次近似圖像f(1).按照同樣的方法再由f(1)求得二次近似圖像f(2),當投影到最后,得到f(N),即完成一次迭代過程.第二次迭代時,以f(N)為初始值并再次投影,如此繼續,直到滿足預定條件為止,其迭代公式為

式中b為迭代次數,aij為投影矩陣系數,pi為投影數據向量,xj為重建后的圖像的第j個像素值,λ為松弛因子,取值范圍為(0,2).

2.4 相位恢復

利用迭代算法重建得到的微分切片數據,構建新的二維函數f(y′,z′)[15]:

由傅里葉變換的微分特性可得

即[16]

式 中F[f(y′,z′)](k,l)是f(y′,z′)的 傅 里 葉 變 換,(k,l)是傅里葉變換后的坐標,因此相位切片數據可以表示為

式中F和F?1分別為傅里葉變換和傅里葉逆變換.因此,利用δ沿著y′和z′方向的一階偏導數切片數據可以準確恢復得到δ.

綜合以上分析,通過ART可實現光柵成像技術與迭代重建算法的完美結合,并對δ進行高信噪比的相位恢復.

3 結果分析

3.1 模擬結果

在對光柵成像數據進行CT重構時,最需要解決的問題是微分相位數據的CT重建和相位信息的恢復.為了驗證本文算法的準確性,對理想相位模型進行光柵成像模擬并對其進行相應重構和恢復.改進Shepp-Logan模型,得到理想相位模型(圖3(a)和圖3(b)).基于同步輻射X射線光源模型,按照表1所示的參數,模擬光柵成像過程及其相應的相位步進掃描,得到每一個像素點有樣品和無樣品時的位移曲線,通過光柵成像信息恢復算法公式((1)式),得到折射率實部衰減率導數?δ/?y的投影正弦圖(圖3(c)).通過(2)式得到投影正弦圖在垂直和平行于X射線傳播方向上的微分投影信息,即和的投影信息,利用(3)式對其進行CT重構后,得到和的切片(圖4). 最后結合兩個方向上的一階微分信息進行相位恢復,得到圖5所示相位恢復結果.

表1 模擬參數Table 1.Simulation parameters.

圖3 (a)Shepp-Logan模型;(b)(a)中白線位置處的δ分布圖;(c)微分相位正弦圖Fig.3.(a)Shepp-Logan phantom;(b)the profile at the position marked with white line in(a);(c)the sinogram of the differential phase information.

圖4 微分相位信息重建結果 (a)δ沿著y′方向微分切片;(b)δ沿著z′方向微分切片Fig.4.Reconstruction result with iterative reconstruction algorithm:(a)Differential phase image of phantom along y′direction;(b)differential phase image of phantom along z′direction.

圖5 折射率迭代算法重建結果 (a)δ切片圖;(b)(a)中白線位置處的δ分布圖Fig.5.Reconstruction result of the refractive index with iterative reconstruction algorithm:(a)The slice of δ;(b)the profile at the position marked with white line in(a).

圖6 折射率濾波反投影重建結果 (a)δ切片圖;(b)(a)中白線位置處的δ分布圖Fig.6.Reconstruction result of the refractive index with FBP:(a)The slice of δ;(b)the profile at the position marked with white line in(a).

利用旋轉不變的折射率實部衰減率的一階偏導數投影信息,通過傳統迭代重建算法,可以得到圖4所示的微分相位信息的切片數據.由圖4(a)和圖4(b)所示切片信息可以很明顯地看到樣品的邊界信息,且求導方向互相垂直.最后恢復得到的δ信息切片如圖5所示.濾波反投影的重建結果如圖6所示.對比圖5(a)和圖6(a)所示切片可以發現,由于投影數量較少,兩種方法均會出現一定量的偽影,均可依據圖像灰度區分不同δ所對應的區域,但迭代重建算法的圖像質量明顯高于濾波反投影的重建結果.從圖5(b)和圖6(b)所示對應切片中相同位置處的δ分布圖可以看出,迭代重建算法得到切片的定量信息更加準確,信噪比更高,具有重要的實際應用意義.

3.2 實驗結果

實驗在上海光源(SSRF)X射線成像及生物醫學應用光束線站(BL13W1)上進行,該線站采用1.9T的Wiggler光源和液氮冷卻雙晶單色器,能提供8—70 keV的單色X射線[17].成像系統由旋轉樣品臺、相位光柵、吸收光柵、納米位移臺和探測器組成,如圖7所示.光柵成像系統距離光源點約34 m,在此處X射線光束基本可看作準平行光,具有較好的空間相干性[18].相位光柵安裝在PI公司的P-611.ZS型納米精度壓電陶瓷位移平臺上,用以實現納米精度相位步進掃描.吸收光柵安裝在一個精密轉臺上,可根據實驗情況實時同步調整兩光柵刻線之間的夾角以改變疊柵條紋的周期.兩光柵之間的水平距離和垂直高度可分別通過安裝在其底座上的精密位移臺遠程自動控制.在采集圖像過程中,使用Photonic Science公司生產的9μm X射線CCD記錄疊柵條紋變化.

為突出驗證光柵成像和迭代重建算法的優勢,實驗中采用能量為20 keV的X射線對樣品在0?—180?范圍內作CT旋轉掃描,CT掃描間隔為3?;移動吸收光柵,在每一個CT投影角度利用多步相移法對相位光柵進行相位步進掃描,相位步進掃描的步長為0.6μm,即每周期掃描4次,具體實驗參數如表2所示,成像樣品為尼龍6(PA6)包裹的聚甲基丙烯酸甲酯(PMMA)棒材,樣品在20 keV光束能量時,折射率等信息如表3所示.

圖7 上海光源光柵成像裝置圖Fig.7.Grating based X–ray imaging system at SSRF.

通過多步相移法,利用光柵成像系統采集得到原始投影信息,按照(1)式所示的方法進行信息提取后,可以分別得到樣品的吸收、相位和散射信息.吸收和散射信息直接通過傳統代數迭代重構算法重構后,可得到樣品的吸收和散射信息切片(圖8);將旋轉變化的微分相位信息進行變換后,利用迭代重建算法重建,可得到δ的微分切片如圖9(a)和圖9(b)所示.再對圖9(a)和圖9(b)所示微分切片進行相位恢復,最后得到圖9(c)所示的樣品折射率實部衰減信息.

表3 樣品信息列表Table 3.Information of samples.

圖8 重建結果 (a)吸收切片;(b)散射切片Fig.8.Experimental results of reconstruction:(a)Slice of absorption information;(b)slice of scattering information.

圖9 相位信息重建結果 (a)δ沿著y′方向微分切片;(b)δ沿著z′方向微分切片;(c)恢復后的δ信息Fig.9.Reconstructed results of phase information:(a)Differential phase image of the sample along y′direction;(b)differential phase image of the sample along z′direction;(c)the slice after phase retrieval.

由以上重構結果可以看出,在少量投影的情況下,通過傳統代數迭代重建算法可以獲得樣品各種襯度的準確信息.在圖8(a)所示的吸收切片中,樣品為低Z材料,吸收系數很小,只能看到由于同步輻射X射線成像邊緣效應而形成的邊界信息,樣品的內部信息與空氣很難通過圖像予以區別.圖8(b)所示的散射信息中只有各種樣品的邊界信息,這是因為樣品是均勻材料,只有在不同材料的交界處才會對X射線形成散射.因此,通過吸收和散射信息很難將PA6,PMMA和空氣等區分開來.而對于相位信息,雖然兩種材料的δ差別很小(表3),但是相位恢復后由圖9(c)可準確區分PA6和PMMA兩種高分子材料.實驗樣品的成分和結構更加復雜時,不同襯度的切片將會突出展示樣品的不同信息,使得成像結果更加豐富,最終形成吸收、相位、散射三位一體的成像模式.

4 結 論

本文基于上海光源同步輻射平行光模型,利用傳統的代數迭代重建算法系統地研究了少量投影情況下X射線光柵成像的信息重建算法,重點研究了一階微分相位信息的CT重建和樣品折射率實部衰減率δ的恢復問題.通過變換的方式,將隨著投影角度變化而變化的一階微分相位信息轉變為旋轉不變的信息,即沿著光線和垂直于光線兩個方向的微分信息,通過傳統的ART,對其直接進行重建得到了δ的一階微分切片,最后再通過傅里葉變換恢復得到δ的實際分布信息.

理論分析和實驗結果表明,該方法不僅可以獲得樣品的吸收切片、沿著兩個不同方向的微分相位切片和散射切片數據,還可以獲得高信噪比的相位數據,在X射線光柵成像中具有可行性.相對于傳統的反投影相位重建算法,基于ART的相位重建算法在少量投影情況下,在光柵成像物理模型和降低輻射劑量等方面都有所改善;與傳統積分對一階微分信息的恢復算法相比,該方法不需要考慮微分數據的邊界條件,也沒有誤差累積效應,可以恢復得到信噪比非常好的相位數據;由于相位恢復時,使用的兩個方向的微分數據均為同一次實驗結果,也不存在圖像配準等問題.因此,本文提出的快速光柵成像方法既保留了光柵成像和ART重建算法的優點,又可以實現在缺角度或者光束不完全覆蓋樣品等條件下的X射線光柵CT成像,在生物醫學、材料和工業等要求低輻射劑量、快速和高精度的領域具有很好的應用前景.

感謝上海光源BL13W1線站工作人員在實驗工作中的幫助.

[1]Pfeiffer F,Bech M,Bunk O,Kraft P,Eikenberry E F,Br?nnimann C,Grünzweig C,David C 2008Nat.Mat.7 134

[2]Momose A,Yashiro W,Maikusa H,Takeda Y 2009Opt.Express17 12540

[3]Wen H H,Bennett E E,Kopace R,Stein A F,Pai V 2010Opt.Express35 1932

[4]Zhu P P,Zhang K,Wang Z L,Liu Y J,Liu X S,Wu Z Y,McDonald S A,Marone F,Stampanoni M 2010Proc.Natl.Acad.Sci.107 13576

[5]Jensen T H,Bech M,Zanette I,Weitkamp T,David C,Deyhle H,Rutishauser S,Reznikova E,Mohr J,Feidenhans’l R,Pfeiffer F 2010Phys.Rev.B82 214103

[6]Zanette I,Bech M,Pfeiffer F,Weitkamp T 2011Appl.Phys.Lett.98 094101

[7]Zanette I,Bech M,Rack A,Le Duc G,Tafforeau P,David C,Mohr J,Pfeiffer F,Weitkamp T 2012Proc.Natl.Acad.Sci.109 10199

[8]Chen B,Zhu P P,Liu Y J,Wang J Y,Yuan Q X,Huang W X,Ming H,Wu Z Y 2008Acta Phys.Sin.57 1576(in Chinese)[陳博,朱佩平,劉宜晉,王寯越,袁清習,黃萬霞,明海,吳自玉2008物理學報57 1576]

[9]Talbot H F 1936Philos.Mag9 401

[10]Qi J C,Ren Y Q,Du G H,Chen R C,Wang Y D,He Y,Xiao T Q 2013Acta Opt.Sin.33 1034001(in Chinese)[戚俊成,任玉琦,杜國浩,陳榮昌,王玉丹,和友,肖體喬2013光學學報33 1034001]

[11]Bech M,Jensen T H,Bunk O,Donath T,David C,Weitkamp T,Le Duc G,Bravin A,Cloetens P,Pfeiffer F 2010Zeitschrift Fur Medizinische Physik20 7

[12]Momose A,Kawamoto S,Koyama I,Suzuki Y 2004Developments in X-Ray Tomography IV5535 352

[13]Zhu P P,Wang J Y,Yuan Q X,Huang W X,Shu H,Gao B,Hu T D,Wu Z Y 2005Appl.Phys.Lett.87 264101

[14]Yang F Q,Zhang D H,Huang K D Wang K,Xu Z 2014Acta Phys.Sin.63 058701(in Chinese)[楊富強,張定華,黃魁東,王鹍,徐哲2014物理學報63 058701]

[15]Kottler C,David C,Pfeiffer F,Bunk O 2007Opt.Express15 1175

[16]Arnison M R,Larkin K G,Sheppard C J R,Smith N I,Cogswell C J 2004J Microsc.214 7

[17]Xiao T Q,Xie H L,Deng B,Du G H,Chen R C 2014Acta Opt.Sin.34 0100001(in Chinese)[肖體喬,謝紅蘭,鄧彪,杜國浩,陳榮昌2014光學學報34 0100001]

[18]Qi J C,Ye L L,Chen R C,Xie H L,Ren Y Q,Du G H,Deng B,Xiao T Q 2014Acta Phys.Sin.63 104202(in Chinese)[戚俊成,葉琳琳,陳榮昌,謝紅蘭,任玉琦,杜國浩,鄧彪,肖體喬2014物理學報63 104202]

PACS:42.25.Kb,41.60.Ap,07.60.Ly,87.59.–e DOI:10.7498/aps.66.054202

Grating based X-ray phase contrast CT imaging with iterative reconstruction algorithm?

Qi Jun-Cheng1)?Chen Rong-Chang2)Liu Bin1)Chen Ping1)Du Guo-Hao2)Xiao Ti-Qiao2)?

1)(School of Information and Communication Engineering,North University of China,Taiyuan 030051,China)
2)(Shanghai Institute of Applied Physics,Chinese Academy of Sciences,Shanghai 201800,China)

18 August 2016;revised manuscript

5 December 2016)

Grating based X-ray imaging technology is a coherent imaging technique that bears tremendous potential in threedimensional tomographic imaging of weak absorption contrast specimens.Three kinds of contrast information including absorption,phase and scattering can be retrieved separately based on a single set of raw projections.However,the grating based X-ray imaging with the conventional phase-retrieval method using the conventional phase-stepping approach andfiltered back projection(FBP)reconstruction algorithm require large amounts of raw data,so that long exposure time and large amounts of radiation dose is accepted by the sample.According to the traditional grating based X-ray imaging system,we propose a low dose,fast,multi-contrast CT reconstruction approach based on the iterative reconstruction algorithm that optimizes dose efficiency but does not share the main limitations of other reported methods.Prior to reconstruction,firstly,the projections are acquired with the phase stepping approach and multi-contrast projections are retrieved from the raw data by conventional retrieval algorithm.Then the rotational variable differential phase projections are converted to rotational invariable projections by means of decomposing the differential phase projections into the rotational invariable projections in two mutually perpendicular derivative directions via the transformation of coordinates.Finally,the absorption,phase and scattering information are reconstructed with the iterative reconstruction algorithm and the phase is retrieved based on the fast Fourier transform(FFT).We validated and assessed the phase reconstruction approach with a numerical simulation on a phase Shepp-Logan phantom.The experiment was performed at the X-ray imaging and biomedical application beam line(BL-13W)in the Shanghai Synchrotron Radiation Facility(SSRF)where 20 keV X-ray from a Si(111)monochromator is emitted.The X-ray interferometer was positioned at 34 m from the Wiggler source.The images were recorded with a scintillator/lens-coupled CCD camera with 2048 pixel×2048 pixel resolution and an effective pixel size of 9μm.The numerical tests and the experimental results demonstrate that,for the small radiation dose deposited in the sample,the iterative reconstruction algorithm provides phase reconstructions of better quality and higher signal to noise ratio than the conventional FBP reconstruction algorithm,and also provides the multi-contrast 3D images,including absorption image,phase image and scattering image.This development is of particular interest for applications where the samples need inspecting under low dose and high speed conditions,and will play an important role in the nondestructive and quantitative imaging in the industry,biomedical and medical diagnosisfields.

grating based X-ray imaging,iterative reconstruction algorithm,phase retrieval

PACS:42.25.Kb,41.60.Ap,07.60.Ly,87.59.–e

10.7498/aps.66.054202

?國家自然科學基金(批準號:11375257,61301259,U1232205)、中北大學校學科研究基金(批準號:2015110246)和山西省自然科學基金(批準號:2015021099)資助的課題.

?通信作者.E-mail:qijuncheng@nuc.edu.cn

?通信作者.E-mail:tqxiao@sinap.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.11375257,61301259,U1232205),the Foundation of North University of China(Grant No.2015110246),and the Natural Science Foundation of Shanxi Province,China(Grant No.2015021099).

?Corresponding author.E-mail:qijuncheng@nuc.edu.cn

?Corresponding author.E-mail:tqxiao@sinap.ac.cn