合成孔徑序列波束形成軸向運動補償方法

孟曉輝 理 華 王君琳

(中國科學院聲學研究所，北京 100190)

?

合成孔徑序列波束形成軸向運動補償方法

孟曉輝 理 華*王君琳

(中國科學院聲學研究所，北京 100190)

合成孔徑序列波束方法是一種新穎的醫(yī)學超聲成像方法，采用兩個階段的波束形成，在傳統(tǒng)的超聲成像系統(tǒng)中實現(xiàn)合成孔徑成像，在不需要存儲和傳輸大量射頻回波數(shù)據(jù)的情況下，提高醫(yī)學超聲圖像的分辨率。該方法的前提是假設成像目標靜止不動，而通過仿真分析發(fā)現(xiàn)運動會造成成像目標位置錯誤。針對此問題，提出一種合成孔徑序列波束形成運動估計和補償方法：首先通過在同一位置連續(xù)發(fā)射接收兩次以獲取用于運動估計的數(shù)據(jù)，然后采用互相關方法對第一階段波束形成得到的低分辨率圖像進行運動估計和補償，再對其做第二次波束形成得到高分辨率圖像。Field II仿真結果顯示，所提出的運動補償方法可以得到正確的目標位置。對于運動速度為0.1 m/s的點目標，運動補償后的平均橫向分辨率與靜止點目標相比僅降低2.03%，對比度降低2.7 dB。對于運動速度為0.2 m/s的囊目標,運動補償后圖像的對比度分辨率較靜止情況僅降低8.53%，進一步說明所提出的運動補償方法有效。

合成孔徑；合成孔徑序列波束形成；運動補償

引言

醫(yī)學超聲成像與同位素掃描、核磁共振成像(MRI)以及計算機 X 射線層析成像(CT)構成現(xiàn)代醫(yī)學影像診斷的4項主要檢查方法，并且較其他幾種成像方式而言，具有無創(chuàng)、實時、診斷成本低等特點，目前已經(jīng)成為許多臨床疾病診斷的首選方法[1]。合成孔徑技術于上個世紀60年代后期被應用到超聲成像領域[2]。90年代以前，合成孔徑技術多數(shù)是用于無損檢測領域，僅有少數(shù)用于醫(yī)學超聲的研究[3-4]。Nikolov 的博士論文專門論述了合成孔徑醫(yī)學超聲成像[5]。 目前在醫(yī)學超聲成像中，較為常用的合成孔徑技術為合成發(fā)射孔徑(STA)方法[4]，即通過單陣元發(fā)射、全孔徑接收，經(jīng)波束形成得到低分辨率圖像(LRI)，再通過將所有發(fā)射陣元得到的LRI相干相加得到高分辨率圖像(HRI)。由于HRI對發(fā)射和接收都實現(xiàn)了動態(tài)聚焦，所以其較傳統(tǒng)的超聲成像方法就有更高的分辨率[6]。但是，這種發(fā)射合成孔徑成像方法面臨3個主要問題：一是單陣元發(fā)射造成的穿透性差和回波信號信噪比低；二是被測組織的運動會降低低分辨率圖像間的相干性，從而影響成像質(zhì)量；三是系統(tǒng)實現(xiàn)過程中需要存儲和傳輸大量的射頻回波(RF)線，系統(tǒng)復雜度很高。

對于上述3個問題，目前有多種改進方法,比如：為了改善單陣元發(fā)射造成的穿透性差和回波信號信噪比低的問題，可以采用多陣元模擬發(fā)射球面波和編碼發(fā)射方法來解決[7-8]；對于低分辨率圖像間運動影響成像質(zhì)量的問題，也提出了多種運動補償方案[9-10]；稀疏合成孔徑是降低系統(tǒng)復雜度的一種方法[11-12]。此外，合成孔徑序列波束形成[13]，也是一種改善系統(tǒng)復雜性高問題的重要方法。

合成孔徑序列波束形成方法是一種分兩個階段進行的波束形成成像方法，可以采用傳統(tǒng)的線陣列探頭實現(xiàn)B模式成像。第一個階段是采用傳統(tǒng)的線性掃描方式，多陣元單焦點聚焦發(fā)射。對于多個陣元接收的RF數(shù)據(jù)進行單焦點波束形成，得到一條掃描線，發(fā)射和接收的焦點相同。這些焦點可以理解為虛擬發(fā)射源和虛擬接收源，構成一個虛擬陣列。第二個階段的波束形成以第一階段得到的掃描線為輸入，利用虛擬陣列進行延時相加做第二次波束形成。合成孔徑序列波束形成方法由于不需要存儲和傳輸大量的單陣元接收的RF線數(shù)據(jù)，其實現(xiàn)復雜度較合成發(fā)射孔徑方法大大降低。此外，由于是多陣元發(fā)射，所以也不存在合成發(fā)射孔徑方法單陣元發(fā)射造成的穿透性差和回波信號信噪比低的問題，因此這種方法近年來得到很多學者的重視。但是，作為一種合成孔徑方法，其中仍然存在被測組織運動影響成像質(zhì)量的問題。

組織位移可以分解為橫向分量和軸向分量，有研究表明影響合成孔徑成像質(zhì)量的主要是軸向位移分量[14]。筆者提出一種合成孔徑序列波束形成的軸向運動補償方法，并對該方法的效果進行了仿真分析，結果表明該方法可以有效地進行運動補償，從而提高了成像質(zhì)量。

1 方法

1.1 合成孔徑序列波束形成方法的原理

合成孔徑序列波束方法的原理如圖1所示。

圖1 合成孔徑序列波束方法原理Fig.1 Model of SASB

1)第一階段波束形成。采用傳統(tǒng)的線性掃描方式，每次發(fā)射和接收采用單焦點，且發(fā)射與接收焦點相同。在成像掃描過程中，每次聚焦發(fā)射接收均可看作是一個虛擬陣元VSi發(fā)射了一定夾角的球面波，形成了一條掃描線li。多個虛擬陣元構成一個虛擬陣列。

2)第二階段波束形成。在傳統(tǒng)成像中，線li上的點P2只利用以VSi為發(fā)射焦點得到的發(fā)射聲場來獲得。實際上，其鄰近掃描線(如li-1和li+1)上也包含P2點的信息。

例如，由li-1提取P2點信息，計算時延值，有

(1)

式中：c為聲速，Zv為虛擬源深度，|P2VSi-1|)為P2點到線li-1上的虛擬源VSi-1的距離；“+”表示P2點深度大于Zv，“-”表示深度小于Zv。

由于P2點與P1點都在以VSi-1為圓心的弧上，所以式(1)的含義為P1點含有P2點的信息。

隨著目標點與虛擬源距離的增加，更多掃描線包含目標點的信息，第二次波束形成使用到的虛擬源數(shù)目K也隨之增多，實現(xiàn)了動態(tài)孔徑，從而實現(xiàn)與深度無關的橫向成像分辨率[13]。K可直接由圖1的幾何關系來計算，有

(2)

式中，α=2arctan(LA/(2Zv))為虛擬源的開角，L(z)為深度z處的聲場寬度，Δ為虛擬源間距。

在依次對所有掃描li進行第二階段波束形成來合成高分辨圖像時，各圖像點可計算為

(3)

式中，x為圖像的橫坐標，z為深度坐標，lk(td(k))表示第一階段所得的掃描線lk上延時td(k)對應的值，w為加權函數(shù)，隨深度z和掃描線的水平位置xk而改變。

式(3)的含義是將K(z)條掃描線上的對應值加權疊加，得到高分辨率的成像點。

1.2 運動對合成孔徑序列波束形成方法的影響

對于成像目標靜止的情況，上述合成孔徑序列波束形成方法能夠不增加硬件系統(tǒng)的復雜性，有效提高橫向分辨率。仿真實驗采用Field II[15]軟件進行，其參數(shù)設置如表1所示。

表1 仿真實驗參數(shù)設置

如圖2所示，在組織靜止情況下，采用合成孔徑序列波束形成方法，得到兩個階段波束形成。

圖2 靜止情況下合成孔徑序列波束形成(顯示動態(tài)范圍60 dB)。(a)第一階段波束形成；(b)第二階段波束形成Fig.2 Images of static points simulation(The dynamic range is 60 dB). (a)The first stage beamforming；(b)The second stage beamforming

對于被測組織存在運動的情況，以組織沿軸向向下以0.2 m/s的速度勻速運動為例，按照表1的參數(shù)進行仿真，采用合成孔徑序列波束形成方法，得到兩個階段的波束形成，如圖3所示。

由圖3(a)可以看出，由于組織存在向下的軸向運動，第一階段成像較靜止情況(見圖2(a))產(chǎn)生了變化。在此情況下，仍按照靜止情況的采樣點位置關系進行第二階段波束形成，會造成成像點位置錯誤，如圖3(b)所示。

圖3 運動情況下合成孔徑序列波束形成(顯示動態(tài)范圍60 dB)。(a)第一階段波束形成；(b)第二階段波束形成Fig.3 Images of points motion simulation (The dynamic range is 60 dB). (a)The first stage beamforming；(b)The second stage beamforming

1.3 運動補償方法

為了解決運動對合成孔徑序列波束形成方法造成的影響，筆者在本研究中提出一種運動補償方法。

首先使用互相關方法進行運動估計。由于相鄰兩條掃描線的發(fā)射陣元不同，其相關性較差，所以不能直接采用第一階段波束形成的相鄰兩條掃描線li與li+1做互相關來估計這兩條掃描之間的運動。由于兩條掃描線的物理位置非常接近，所以用同一位置連續(xù)兩次發(fā)射間的運動來估計近似相鄰兩條掃描線間的運動。具體的數(shù)據(jù)獲取和處理過程如圖4所示，具體實現(xiàn)流程如下：

(4)

(5)

如果所有掃描線以第一條掃描線為基準進行補償，即dcomp(1)=0，則第i條掃描線的補償值為

(6)

為了提高運動估計精度，可以先對掃描線進行插值，也可以對時延估計結果進行插值[11]，即

(7)

則

(8)

采用式(8)進行運動補償，得到補償后的低分辨率掃描線cli。其中，Round()為四舍五入取整函數(shù)。

3)第二階段波束形成。以運動補償后得到的掃描線為輸入，利用式(3)做波束形成，得到高分辨率掃描線，即

(9)

式中，w(xk，z)采用漢明窗。

圖4 合成孔徑序列波束形成運動補償方法Fig.4 Illustration of SASB motion compensation method

1.4 仿真測試方法

為了驗證所提出方法的有效性，用Field II軟件進行了仿真實驗，與表1不同或增加的仿真參數(shù)如表2所示。

表2 仿真實驗參數(shù)設置

仿真中的陣元總數(shù)為128，每次發(fā)射和接收的陣元數(shù)都為64，具體成像過程為：

2)按照本文1.3節(jié)中描述的方法，以第一根掃描線為基準進行運動估計和補償。

3)第二階段波束形成，以運動補償后得到的掃描線為輸入，利用式(3)做波束形成，得到高分辨率掃描線cli(i=1,2,…,65)。

為了驗證本方法對于點目標的運動補償效果，

設空間有8個點，深度分別為10～80 mm，間隔為10 mm。用橫向-6 dB寬度表征橫向分辨率，用第一旁瓣高度表征對比度。由于發(fā)射和接收陣元數(shù)確定后，發(fā)射和接收焦點深度以及組織的軸向運動速度對本方法的成像效果影響較大，所以分別對焦點深度為25、35、45 mm和點目標靜止及豎直向下運動速度為0.1、0.15、0.2 m/s共12種情況進行仿真實驗。

為了驗證本方法對于囊目標的運動補償效果，在空間設置兩個囊目標，深度分別為40和60 mm，半徑都為4 mm，豎直向下運動，運動速度為0.2 m/s，焦點深度為25 mm。定義對比度分辨率[8]為

(10)

式中，CRref為參考圖像的對比度分辨率，CRobj為測試圖像的對比度分辨率。對比度分辨率CR[16]可定義為

(11)

式中，Ic為囊腫內(nèi)的強度，Is為背景的強度。

文中取囊腫在圖像中的橫向對稱位置為參考圖像。

2 結果

2.1 點目標仿真結果

圖5 點目標仿真結果。(a)未做運動補償；(b)靜止點目標；(c)做運動補償Fig.5 Image results of point simulation. (a)Uncompensated；(b)Stationary points；(c)Compensated

圖5為合成孔徑序列波束形成方法對點目標的仿真結果，其中焦點深度為25 mm，運動速度為0.1 m/s，顯示動態(tài)范圍為40 dB；(a)為未做運動補償?shù)某上窠Y果，(b)為靜止情況下的成像結果，(c)為用本方法以第一條掃描線為基準做運動補償后的成像結果。由圖像可以看出，運動對于合成孔徑序列波束形成方法的成像效果影響較嚴重，未做運動補償?shù)某上窠Y果嚴重失真。所提出的運動補償方法可以有效避免運動對合成孔徑序列波束形成方法成像效果的影響，成像結果與靜止點目標相似。

圖6 不同焦點深度情況下的橫向-6 dB寬度。(a)Zv=25 mm；(b)Zv=35 mm；(c)Zv=45 mmFig.6 Lateral width at -6 dB for different focus depth. (a)Zv=25 mm；(b)Zv=35 mm；(c)Zv=45 mm

圖6為不同焦點深度點目標靜止和軸向運動速度分別為0.1、0.15、0.2 m/s時運動補償后的橫向-6 dB寬度，其中圖6(a)～(c)對應的焦點深度分別為25、35、45 mm?？梢钥闯?，當運動速度為0.2 m/s時，多數(shù)點目標對應的橫向-6 dB寬度都有所增加，而運動速度0.1和0.15 m/s時的結果與靜止時的情況類似。

進一步地，表3給出了焦點深度為25 mm時不同運動速度下8個散射點橫向-6 dB寬度均值，以及相對于靜止情況的增加百分比。當運動速度為0.1 m/s時，橫向-6 dB寬度均值僅增加2.03%，隨著速度加大，橫向-6 dB寬度均值增加幅度增大。

表3 Zv=25 mm時8個散射點橫向-6 dB寬度均值

圖7為不同焦點深度點目標靜止和軸向運動速度分別為0.1、0.15、0.2 m/s時運動補償后的第一旁瓣高度，其中圖7(a)～(c)對應的焦點深度分別為25、35、45 mm。可以看出，運動速度為0.1 m/s時，補償后的第一旁邊高度與靜止情況時的類似，隨著速度的增加，第一旁瓣高度有所升高，尤其是運動速度達到0.2 m/s時升高比較明顯。

圖7 不同焦點深度情況下的第一旁瓣高度。(a)Zv=25 mm；(b)Zv=35 mm；(c)Zv=45 mmFig.7 Height of the first sideband for different focus depth. (a)Zv=25 mm；(b)Zv=35 mm；(c)Zv=45 mm

表4給出了焦點深度為35 mm時8個散射點第一旁瓣高度均值的增加情況。其中，速度為0.1 m/s時補償后的第一旁瓣高度均值較靜止情況增加2.7 dB，速度為0.15 m/s時增加5.4 dB，速度為0.2 m/s時增加10.3 dB。

表4 焦點深度為35 mm時8個散射點第一旁瓣高度均值

2.2 囊目標仿真結果

圖8為合成孔徑序列波束形成方法對囊目標的仿真結果，顯示動態(tài)范圍為40 dB，其中(a)為未做運動補償?shù)某上窠Y果，(b)為靜止情況下的成像結果，(c)為用本方法以第一條掃描線為基準做運動補償后的成像結果。由圖像可以看出，運動對于合成孔徑序列波束形成方法的影響較大，未做運動補償?shù)某上窠Y果嚴重失真，深處的囊目標有一部分已經(jīng)不在成像區(qū)域內(nèi)。所提出的運動補償方法對于囊目標可以有效避免運動對合成孔徑序列波束形成方法成像效果的影響，成像結果與靜止囊目標相似。

圖8 囊目標仿真結果。(a)未做運動補償；(b)靜止囊目標；(c)做運動補償Fig.8 Image results of cysts motion simulation. (a)Uncompensated；(b)Stationary cysts；(c) Compensated

對于上述設定的囊目標，根據(jù)式(10)計算運動補償后的圖像(見圖8(c))較囊腫靜止情況的圖像(見圖8(b))的對比度分辨率僅降低8.53%。

3 討論和結論

由于不需要存儲和傳輸大量單陣元接收的RF線數(shù)據(jù)，合成孔徑序列波束方法的實現(xiàn)復雜度較合成發(fā)射孔徑方法的實現(xiàn)復雜度大大降低。合成孔徑成像的前提是假設成像目標靜止不動，對于傳統(tǒng)的合成發(fā)射孔徑方法，成像目標的運動會造成低分辨率圖像間產(chǎn)生某種錯位，從而造成高分辨率圖像對比度的嚴重下降[10]。而從本文對組織運動對合成孔徑序列波束方法影響的分析中可以看出(見圖3)，成像目標運動最主要的影響是會造成最終成像結果中目標點的位置錯誤。在采用所提出的運動估計和補償方法后，從點目標和囊目標的仿真結果看(見圖5、8)，運動造成的位置錯誤問題都得到了解決。但是，由于采樣率等因素對運動估計和補償精度的限制，比較靜止情況，從點目標的仿真結果看，運動補償后的成像結果的橫向分辨率和對比度均有下降，運動速度越快，下降越明顯：當運動速度為0.1 m/s時，橫向-6 dB寬度均值和第一旁瓣高度均值較靜止情況分別增加2.03%和2.7 dB；而運動速度為0.2 m/s時，則分別增加21.62%和10.3 dB。從囊目標的仿真結果看，速度為0.2 m/s時，運動補償后的成像結果的對比度分辨率下降8.53%。值得注意的是，影響運動補償后成像效果實際是兩次脈沖發(fā)射間組織的位移，所以與彩色血流成像類似，可以依據(jù)具體情況通過調(diào)整系統(tǒng)的脈沖重復頻率，使運動補償保持較好的效果。此外，提高采樣率或者信號處理過程中插值的倍數(shù)，可以降低運動估計和補償誤差，從而進一步提高存在組織運動情況下合成孔徑序列波束形成方法的成像分辨率和對比度。

為了獲得用于運動估計的數(shù)據(jù)，本研究用運動補償方法在同一位置連續(xù)發(fā)射接收2次，這樣做的優(yōu)點：一是用相同陣元獲取的數(shù)據(jù)做運動估計，相關性高、魯棒性強，否則速度的快慢變化容易引起數(shù)據(jù)相關性差而造成估計失?。欢强蓪崿F(xiàn)性強，硬件修改小或者不用修改，有利于現(xiàn)有設備的技術改造；三是成像的效果較好，可以有效地避免組織運動造成的成像點位置錯誤問題。但此方法也存在一些缺點：一是數(shù)據(jù)計算量增大，波束形成計算量是不補償?shù)?倍，另外還增加了互相關算法的計算量,但隨著硬件水平的增加和算法的優(yōu)化，計算量大的問題將會得到解決；二是由于要在每個位置發(fā)射2次，成像的幀頻率降低1倍，因此限制了算法在高幀頻場景的應用。

本研究首先分析了組織運動對合成孔徑序列波束形成方法成像效果的影響，發(fā)現(xiàn)組織的軸向運動造成成像點的位置錯誤；然后針對這一問題提出了一種合成孔徑序列波束形成方法的運動補償算法，用在第一階段波束形成過程中，同一位置連續(xù)發(fā)射接收2次以獲取用于運動估計的數(shù)據(jù)。仿真實驗結果顯示：當組織運動速度不大時，如0.1 m/s，所提出的運動補償方法對于運動點目標可以得到與靜止點目標相似的橫向分辨率；在對運動速度為0.2 m/s的囊目標運動補償后，圖像的對比度分辨率較靜止情況僅降低8.53%，說明該運動補償方法有效。

下一步研究的重點，是如何在不降低成像幀頻的前提下做有效的運動補償。此外，筆者提出的運動補償方法主要針對組織的軸向運動，沒有考慮組織的橫向運動。仔細分析組織橫向運動對合成孔徑序列波束形成方法成像效果的影響，并有針對性地研究橫向運動補償方法，也是下一步研究的方向。

[1] 林軼翰，陳思平. 醫(yī)學超聲影像技術及展望 [J]. 醫(yī)療裝備, 2006, 19(6):13-15.

[2] Fort JR, Landmeier PC, Morgan DJ,etal. Synthetic aperture ultrasound imaging system [P]. U.S. Patent: 5465722. 1995-11-14.

[3] Burckhardt CB, Grandchamp PA, Hoffmann H. An experimental 2 MHz synthetic aperture sonar system intended for medical use [J]. IEEE Transactions on Sonics and Ultrasonics, 1974, 21(1): 1-6.

[4] Nagai K. A new synthetic-aperture focusing method for ultrasonic B-scan imaging by the Fourier transform [J]. IEEE Transactions on Sonics and Ultrasonics, 1985, 32(4): 531-536.

[5] Nikolov SI. Synthetic aperture tissue and flow ultrasound imaging [D]. Lyngby: Technical University of Denmark, 2001.

[6] Jensen JA, Nikolov SI, Gammelmark KL,etal. Synthetic aperture ultrasound imaging [J]. Ultrasonics, 2006, 44: e5-e15.

[7] Karaman M, Li PC, O’Donnell M. Synthetic aperture imaging for small scale systems [J]. IEEE Trans Ultrason Ferroelec Freq Contr, 1995, 42(3): 429-442.

[8] Gammelmark KL, Jensen JA. Multielement synthetic transmit aperture imaging using temporal encoding [J]. IEEE Trans Med Imaging, 2003, 22(4): 552-563.

[9] Gammelmark KL, Jensen JA. Duplex synthetic aperture imaging with tissue motion compensation [C] //Yuhas DE, eds. 2003 IEEE Symposium on Ultrasonics (Volume 2 ). Hawaii: IEEE, 2003: 1569-1573.

[10] Yiu BYS, Tsang IKH, Yu ACH. A modified synthetic aperture imaging approach with axial motion compensation [C] // Waters KR, eds. Proceedings of IUS 2008. Beijing: IEEE, 2008: 1254-1257.

[11] Lockwood GR, Talman JR, Brunke SS. Real-time 3-D ultrasound imaging using sparse synthetic aperture beamforming [J]. IEEE Trans Ultrason Ferroelec Freq Contr, 1998, 45: 980-988.

[12] Hazard CR, Lockwood GR. Theoretical assessment of a synthetic aperture beamformer for real-time 3-D imaging [J]. IEEE Trans Ultrason Ferroelec Freq Contr, 1999, 46: 972-980.

[13] Kortbek J, Jensen JA, Gammelmark KL. Sequential beamforming for synthetic aperture imaging [J]. Ultrasonics, 2013, 53(1): 1-16.

[14] Trahey GE, Nock LP. Synthetic receive aperture imaging with phase correction for motion and for tissue inhomogeneities. II. Effects of and correction for motion [J]. IEEE Trans Ultrason Ferroelec Freq Contr, 1992, 39(4): 496-501.

[15] Angelsen BAJ. A theoretical study of the scattering of ultrasound from blood [J]. IEEE Trans Biomed Eng, 1980, BME-27: 61-67.

[16] Foster SG, Embree PM, O’Brien WD. Flow velocity profile via time-domain correlation: Error analysis and computer simulation [J]. IEEE Trans Ultrason Ferroelec Freq Contr, 1990, 37(3): 164-175.

An Axial Motion Compensation Method for Synthetic Aperture Sequential Beamforming

Meng Xiaohui Li Hua*Wang Junlin

(InstituteofAcoustics,ChineseAcademyofSciences,Beijing100190,China)

Synthetic aperture sequential beamforming (SASB) is a two-stage beamforming procedure, which can be applied to B-mode imaging with an implementation of low complexity compared to the traditional synthetic transmit aperture (STA) imaging. However, like the STA imaging, the SASB is susceptible to motion artifacts due to the summation of a number of RF-data created at different time instances. In this paper, a tissue motion estimation and compensation method was proposed. First, the inter-firing motion was estimated by cross-correlating with extra firings. At last, the second stage beamformer created a set of high resolution image points by combining information from multiple first stage focused scan lines. Field II simulation results showed that the motion compensation method proposed in this paper exhibited enhanced ability to correct location of moving tissue. When the points were moving at 0.1 m/s, the average lateral resolution after motion compensation was only reduced by 2.03% comparing to the static case, and contrast reduction was 2.7 dB. Moreover, when cysts were moving at 0.2 m/s, the contrast resolution of image after motion compensation was only reduced by 8.53% comparing to the static case.

synthetic aperture; synthetic aperture sequential beamforming; motion compensation

10.3969/j.issn.0258-8021. 2015. 05.006

2014-07-14， 錄用日期:2015-07-06

國家自然科學基金(11204346)；中國科學院聲學研究所所長擇優(yōu)基金(Y454221341)；中國科學院聲學研究所創(chuàng)新前瞻項目(Y154211341)

R445.1

A

0258-8021(2015) 05-0558-08

*通信作者(Corresponding author)， E-mail: leehwa@mail.ioa.ac.cn