基于獨(dú)立分量分析法的稀疏陣列穿墻成像雷達(dá)直達(dá)波干擾抑制

張 馳 李?lèi)傷?周智敏

(國(guó)防科技大學(xué)電子科學(xué)與工程學(xué)院 長(zhǎng)沙 410073)

基于獨(dú)立分量分析法的稀疏陣列穿墻成像雷達(dá)直達(dá)波干擾抑制

張 馳*李?lèi)傷?周智敏

(國(guó)防科技大學(xué)電子科學(xué)與工程學(xué)院 長(zhǎng)沙 410073)

強(qiáng)的墻體直達(dá)波嚴(yán)重干擾了穿墻成像雷達(dá)的目標(biāo)信號(hào)，而采用稀疏陣列天線的穿墻雷達(dá)，由于通道數(shù)少使得基于單次快拍數(shù)據(jù)分離目標(biāo)和干擾信號(hào)難度更大。該文根據(jù)實(shí)際通道每次收發(fā)獲取的距離像存在微弱起伏的現(xiàn)象，提出一種采用單通道多次快拍數(shù)據(jù)建立觀察矩陣，并使用獨(dú)立分量分析分離直達(dá)波與目標(biāo)信號(hào)的方法。仿真及實(shí)測(cè)數(shù)據(jù)處理結(jié)果表明，在低信雜比(-30 dB)的條件下，該方法仍可有效分離目標(biāo)和干擾成分，達(dá)到較好的直達(dá)波抑制效果。

穿墻成像雷達(dá)(TWIR)；直達(dá)波抑制；獨(dú)立分量分析(ICA)；稀疏陣列

1 引言

穿墻成像雷達(dá)(Through-the-Wall Imaging Radar, TWIR)可對(duì)墻后態(tài)勢(shì)進(jìn)行感知，因而在軍事、反恐、救災(zāi)等領(lǐng)域有重大的應(yīng)用前景，近年來(lái)受到人們的廣泛重視[1,2]。如何消除墻體直達(dá)波對(duì)目標(biāo)信號(hào)造成的嚴(yán)重干擾，是TWIR需要解決的關(guān)鍵問(wèn)題之一。對(duì)于靜止目標(biāo)，一種直接的方法是事先錄取背景信號(hào)，然后采用背景相消法去除墻體直達(dá)波在內(nèi)的靜止干擾[3]，但該方法的局限性在于需要獲得成像區(qū)域的先驗(yàn)知識(shí)。以奇異值分解(Singular Value Decomposition, SVD)、主元分析(Principal Component Analysis, PCA)為代表的子空間法，無(wú)需場(chǎng)景的先驗(yàn)知識(shí)，按照奇異值(或特征值)的大小，將原始信號(hào)投影到不同的子空間中，在一定程度上可以達(dá)到抑制干擾的目的[4,5]。然而在信雜比較低的情況下，弱目標(biāo)的子空間特征往往無(wú)法突顯。子空間法的本質(zhì)是把原始信號(hào)分解為互不相關(guān)的分量，相比之下，獨(dú)立分量分析法(Independent Component Analysis, ICA)假設(shè)原始信號(hào)由若干統(tǒng)計(jì)獨(dú)立的源信號(hào)線性混合而成，這在許多應(yīng)用中更符合信號(hào)的本質(zhì)結(jié)構(gòu)[6]，因而能更有效地提取目標(biāo)信號(hào)。

單收發(fā)系統(tǒng)的TWIR常采用合成孔徑模式，通過(guò)掃描獲取 B-scan(距離向-方位向 2維觀察數(shù)據(jù))矩陣進(jìn)行雜波抑制[4－7]。而采用陣列天線的TWIR，則多采用單發(fā)多收或多發(fā)多收模式，通過(guò)快拍方式(指雷達(dá)系統(tǒng)完成一次收發(fā)通道的遍歷)由各收發(fā)通道對(duì)場(chǎng)景進(jìn)行觀察[8,9]。其陣元的稀疏分布將導(dǎo)致觀察信號(hào)的欠缺，這使得單次快拍難以獲取足夠的觀察信號(hào)用于分離其中目標(biāo)和干擾分量。

針對(duì)上述問(wèn)題，本文首先推導(dǎo)了1維距離像的形成過(guò)程及其數(shù)學(xué)表達(dá)，引入了ICA數(shù)據(jù)模型。之后分析了實(shí)際通道每次收發(fā)獲取的信號(hào)存在微弱起伏的現(xiàn)象，并研究了 ICA對(duì)確定性源信號(hào)的適用性，在此基礎(chǔ)上提出一種采用單通道多次快拍數(shù)據(jù)建立觀察矩陣，使用ICA分離直達(dá)波與目標(biāo)信號(hào)的方法。仿真及實(shí)測(cè)數(shù)據(jù)表明，在信雜比低至-30 dB的情況下，該方法仍可分離目標(biāo)和直達(dá)波干擾，解決了稀疏陣列 TWIR因單次快拍觀察信號(hào)不足造成的無(wú)法有效提取墻后目標(biāo)的問(wèn)題。

2 1維距離像形成

假設(shè)稀疏陣列 TWIR采用步進(jìn)頻連續(xù)波(Stepped-Frequency Continuous Wave, SFCW)體制，其頻點(diǎn)數(shù)為N，則陣列天線第l個(gè)收發(fā)通道，在第k個(gè)頻點(diǎn)發(fā)射的子脈沖為[10]：

其中 T,f0和Δf分別為子脈沖持續(xù)時(shí)間、初始頻率和頻率步進(jìn)增量；?l為通道l初始相位(不失一般性，設(shè) ?l=0)；矩形函數(shù)為：

經(jīng)理想解調(diào)，對(duì)輸出的復(fù)信號(hào)在其中心處進(jìn)行采樣的結(jié)果為[11,12]：

若將Ul(k)視為頻域樣點(diǎn)，則對(duì)所有頻點(diǎn)構(gòu)成的序列做N點(diǎn)逆離散傅里葉變換(Inverse Discrete Fourier Transform, IDFT)，可得到理想情況下通道l的1維距離像[12]：

假設(shè)墻體是勻質(zhì)的，則其前表面也可等效為理想散射點(diǎn)[11]，那么Hl(n)中直達(dá)波分量可表示為(n),d∈{1 , 2, … , P }。由于實(shí)際中直達(dá)波能量遠(yuǎn)大于目標(biāo)散射點(diǎn)回波能量，且墻體會(huì)對(duì)后者造成較大衰減[13]，這使得直達(dá)波分量往往將目標(biāo)分量“淹沒(méi)”。為了有效提取目標(biāo)信號(hào)，抑制直達(dá)波干擾，首先需要對(duì)Hl(n)中各分量進(jìn)行有效分離。

3 ICA模型

ICA作為一種盲源分離方法[14－16]，近年來(lái)被引入到TWIR領(lǐng)域中。ICA假設(shè)觀察變量是若干統(tǒng)計(jì)獨(dú)立的源變量線性混合而成，并依據(jù)非高斯性最大化準(zhǔn)則，從觀察變量中估計(jì)出各源變量。其基本模型如下[14]：

由于實(shí)際中只能獲取各隨機(jī)變量的采樣，若采樣數(shù)為N，則式(6)，式(7)中X和S分別擴(kuò)展為M×N維和P×N維采樣矩陣，其每一行,分別表示N點(diǎn)源信號(hào)和觀察信號(hào)，并且：

其中aji表示矩陣A的元素。由式(6)和式(7)可知，上述模型成立的前提是M,N和P存在以下關(guān)系：

其中rank表示矩陣的秩。由于通常信號(hào)的采樣數(shù)較大，于是有 min(M , N)=M ≥ P 。以上關(guān)系說(shuō)明，ICA只能求解不多于觀察信號(hào)數(shù)量的源信號(hào)。此外，由于S,A均未知，使得S中各源信號(hào)si的幅度及次序具有模糊性[15]。

4 稀疏陣列的ICA方法

4.1 距離像起伏分析

在理想情況下，只要收發(fā)通道位置確定且目標(biāo)靜止，每次快拍獲取的Hl(n)是相同的，并不能增加有效的觀察信號(hào)。但是，我們注意到在實(shí)際雷達(dá)系統(tǒng)中，即使是單一收發(fā)通道的距離像也會(huì)受系統(tǒng)相位漂移、雜波環(huán)境變化、多路徑等多種因素的影響[17]，使其在每次收發(fā)中存在起伏。圖1為對(duì)實(shí)際TWIR系統(tǒng)某通道多次收發(fā)獲取的距離像在固定距離單元進(jìn)行取樣并去除均值的結(jié)果(圖中一條曲線對(duì)應(yīng)一個(gè)距離單元)。其中橫軸為收發(fā)次數(shù)，縱軸表示幅度，n為距離單元編號(hào)。可以看出在每一次收發(fā)中，不同距離單元處Hl(n)均存在微弱起伏(一般

在均值的0.1%左右)。這樣可將通道l第j次收發(fā)獲取的距離像寫(xiě)作Hlj(n)。并且通過(guò)比較圖 1(a)與圖1(b)可知相鄰距離單元的變化趨勢(shì)類(lèi)似，而相距較遠(yuǎn)的則差異較大，由于Hlj(n)的能量主要集中在各的峰值處，于是僅考慮各峰值處的變化，可將式(5)改寫(xiě)為：

4.2 ICA分離確定性信號(hào)原理

圖1 距離像起伏Fig. 1 Fluctuation of range profile

在E{|wHX|2}=1的約束下，ICA算法通過(guò)迭代獲取JG的極值點(diǎn)wopt，以求得源信號(hào)wHoptX，其迭代步驟如下[16]：

步驟1 任意選擇初始權(quán)向量w0；

步驟2

如第3節(jié)所述，實(shí)際X為觀察信號(hào)的采樣矩陣，根據(jù)各態(tài)歷經(jīng)性假設(shè)，以上迭代過(guò)程中各統(tǒng)計(jì)量均由時(shí)間平均進(jìn)行計(jì)算[19]，即E(?)≈，其中算子表示時(shí)間平均。進(jìn)一步由式(6)可知 wHX=wHAS= zHS，可見(jiàn)，即使S為確定性信號(hào)，在時(shí)間平均的意義下，上述求解過(guò)程依然成立，但條件是(表示由時(shí)間平均計(jì)算的峰度)不能為 0[16]，否則極值點(diǎn)無(wú)法求取。另外，當(dāng)==…=時(shí)，算法可能得到錯(cuò)誤的結(jié)果[19]。

4.3 仿真分析

表1 源信號(hào)及混合矩陣參數(shù)Tab. 1 Parameters of source signals and mixture matrix

圖2(a)，圖2(d)，圖2(g)為s1,s2峰值位置不同時(shí)的源信號(hào)圖(實(shí)際計(jì)算時(shí)采用的 IDFT點(diǎn)數(shù)大于頻點(diǎn)數(shù)，使得源信號(hào)峰值幅度小于1)。圖2(b)，圖2(e)，圖2(h)為對(duì)s1,s2乘以混合系數(shù)并疊加后得到的觀察值x1,x2，注意到s2成分幾乎被s1淹沒(méi)。圖2(c)，圖2(f)顯示經(jīng)過(guò)ICA處理后，相比混合之前，除幅度發(fā)生變化外，s1,s2均被有效分離，這說(shuō)明所提方法具有可行性。但圖 2(i)結(jié)果顯示當(dāng)目標(biāo)位置較近時(shí)，分離后的s1,s2產(chǎn)生畸變，峰值出現(xiàn)重疊成分，這是由于當(dāng)Δn趨于0時(shí)，,趨同造成的。

4.4 基于ICA的墻體直達(dá)波抑制流程

實(shí)際應(yīng)用時(shí)可通過(guò)M(M＞P)次快拍，使得每個(gè)通道l( l=1,2,… ,Q(Q -1))均獲取足夠觀察信號(hào)構(gòu)成觀察矩陣 Xl，并通過(guò)ICA求得各 Xl中的源信號(hào)Sl及混合矩陣 Al。依據(jù)式(6)三者滿足以下關(guān)系：

圖2 si峰值位置不同時(shí)ICA分離結(jié)果Fig. 2 Results of ICA when locations of peak of si differ

圖3 基于ICA的稀疏陣列直達(dá)波抑制方法流程圖Fig. 3 Flow chart of wall clutter mitigation for sparse array antenna based on ICA

5 實(shí)測(cè)數(shù)據(jù)分析

實(shí)驗(yàn)使用4塊Vivaldi平板天線間隔均勻設(shè)置為 3 m長(zhǎng)的稀疏天線陣列。一次快拍采取順序單發(fā)多收工作模式，即從陣元1到陣元4，任意時(shí)刻只有1個(gè)天線發(fā)射，余下3個(gè)天線順序接收。因此，雷達(dá)系統(tǒng)共有 12個(gè)收發(fā)通道。發(fā)射信號(hào)采用SFCW體制，掃頻范圍為317～1815 MHz，頻率步進(jìn)增量為2 MHz，共750個(gè)頻點(diǎn)。實(shí)驗(yàn)系統(tǒng)參數(shù)如表2。

表2 雷達(dá)系統(tǒng)參數(shù)Tab. 2 Parameters of radar system

實(shí)驗(yàn)場(chǎng)景以及坐標(biāo)系設(shè)置如圖4所示，其中陰影部分均為混凝土墻。實(shí)驗(yàn)中天線陣列放置在左側(cè)空房間中，電磁波需穿透墻A，對(duì)墻后3.2 m×7 m的房間進(jìn)行成像。陣列與墻A平行并與其前表面相距6.25 m，墻體厚度為0.35 m。目標(biāo)為一30 cm三面角，放置在墻后房間中軸線上，距離墻A后表面2 m處，目標(biāo)距離地面1 m，與天線陣列處于同一水平面。

首先在不放置目標(biāo)的情況下對(duì)成像區(qū)域進(jìn)行一次快拍，以獲取各通道的背景干擾。之后將目標(biāo)放置在指定位置進(jìn)行50次快拍。在有目標(biāo)的情況下，某次快拍通道1(陣元1發(fā)射陣元2接收)獲取的成像區(qū)域距離像如圖5(a)所示。在直達(dá)波以及其它雜波的嚴(yán)重干擾下，距離像中目標(biāo)分量幾乎不可見(jiàn)。從中減去該通道的背景干擾后，可觀察到目標(biāo)距離像，如圖5(b)所示，其峰值功率與圖5(a)中直達(dá)波峰值功率之比約為-30 dB。

之后使用本文提出的方法，對(duì)50次快拍中各通道獲取觀察矩陣 Xl進(jìn)行ICA處理，得到 Al及Sl。同樣以通道1為例，經(jīng)分解后的 X1中第20行 x1,20中部分信號(hào)分量如圖6所示?？梢钥吹綀D6(c)中分量3為目標(biāo)分量。而由于直達(dá)波峰值位于第6個(gè)距離單元附近(圖5(a)所示)，由此可判定圖6(a)中分量1是直達(dá)波分量，可將其去除。此外圖6中還存在著其它信號(hào)分量，這是由多路徑信號(hào)以及周?chē)h(huán)境雜波等干擾造成的，這需要進(jìn)一步的檢測(cè)方法將其與目標(biāo)分量予以區(qū)分。以上結(jié)果說(shuō)明ICA可有效分離各通道觀察信號(hào)中的目標(biāo)和干擾分量。

在去除干擾后，取第20次快拍對(duì)應(yīng)的2維矩陣B20，利用后向投影(BackProjection, BP)算法進(jìn)行成像的結(jié)果如圖7所示。其中圖7(a)為對(duì)原始數(shù)據(jù)進(jìn)行成像的結(jié)果，在嚴(yán)重的直達(dá)波干擾下幾乎無(wú)法從中獲取目標(biāo)的信息。而采用本文方法處理后，成像結(jié)果中出現(xiàn)了目標(biāo)成分，如圖7(b)所示。圖7(c)，圖 7(d)分別為對(duì)單次快拍數(shù)據(jù)進(jìn)行 ICA以及采用SVD進(jìn)行直達(dá)波抑制后的圖像(具體方法可參考文獻(xiàn)[6,11])，可以看到經(jīng)過(guò)以上方法處理后，從圖像中仍然無(wú)法分辨出目標(biāo)。通過(guò)計(jì)算圖像信雜比進(jìn)一步評(píng)價(jià)各方法的效果，其計(jì)算式為[20]：

式中，A1為圖像中目標(biāo)信號(hào)的分布區(qū)，N1為區(qū)域A1中的數(shù)據(jù)點(diǎn)數(shù)；A2為包含目標(biāo)在內(nèi)的雜波分布區(qū)，N2為區(qū)域A2中的數(shù)據(jù)點(diǎn)數(shù)；為相應(yīng)像素點(diǎn)的強(qiáng)度。這里選取以目標(biāo)為中心0.5 m×0.3 m的矩形區(qū)域作為A1，(如圖7(b)中紅色方框所示)，整個(gè)圖像作為A2，各圖像信雜比計(jì)算的結(jié)果如表3所示。

圖4 實(shí)驗(yàn)場(chǎng)景Fig. 4 Experimental scene

圖5 通道1及目標(biāo)距離像Fig. 5 Range profile of channel 1 and the target

圖6 部分信號(hào)分量Fig. 6 Parts of signal components

圖7 BP成像結(jié)果Fig. 7 Results of BP imaging

表3 不同圖像信雜比Tab. 3 Signal-to-clutter ratio of different images

相比不做任何處理的情況，使用 SVD算法降低了圖像的信雜比，這是由于本實(shí)驗(yàn)的陣列配置和收發(fā)模式，使得墻體子空間和目標(biāo)子空間存在較大交疊[11]，導(dǎo)致算法在去除墻體干擾的同時(shí)不可避免地?fù)p失較多目標(biāo)能量。而直接使用單次快拍數(shù)據(jù)進(jìn)行ICA，解得的分量混合了各通道的直達(dá)波和目標(biāo)信號(hào)，基于類(lèi)似的原因該方法同樣降低了圖像的信雜比。相比之下，本文方法有效地分離了每個(gè)通道中的干擾和目標(biāo)分量，較好地抑制了直達(dá)波，提高了圖像的信雜比。

6 結(jié)論

針對(duì)稀疏陣列單次快拍觀察信號(hào)不足的問(wèn)題，本文提出采用單通道多次快拍數(shù)據(jù)建立觀察矩陣，使用ICA分離直達(dá)波與目標(biāo)信號(hào)的方法。仿真和實(shí)驗(yàn)結(jié)果證明，該方法可有效提取目標(biāo)分量，達(dá)到較好的直達(dá)波抑制效果。對(duì)于實(shí)測(cè)數(shù)據(jù)中存在的殘余干擾問(wèn)題，需要研究進(jìn)一步的檢測(cè)手段將其與目標(biāo)成分進(jìn)行區(qū)分，以去除由多路徑信號(hào)以及周?chē)h(huán)境雜波產(chǎn)生的虛假目標(biāo)。

[1] Baranoski E J. Through-wall imaging: historical perspective and future directions[J]. Journal of the Franklin Institute, 2008, 345(6): 556-569.

[2] Amin M G and Ahmad F. Through-the-Wall Radar Imaging: Theory and Applications[M]. Oxford: E-Reference Signal Processing, Elsevier, 2013.

[3] Debes C, Amin M G, and Zoubir A M. Target detection in single-and multiple-view through-the-wall radar imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(5): 1349-1361.

[4]Tivive F H C, Bouzerdoum A, and Amin M G. An SVD-based approach for mitigating wall reflections in through-the-wall radar imaging[C]. IEEE Radar Conference (RADAR), Kansas City, USA, 2011: 519-524.

[5] Kabourek V, ?erny P, and Mazánek M. Clutter reduction based on principal component analysis technique for hidden objects detection[J]. Radio Engineering, 2012, 21(1): 464-470. [6] Gaikwad A N, Singh D, and Nigam M J. Application of clutter reduction techniques for detection of metallic and low dielectric target behind the brick wall by stepped frequency continuous wave radar in ultra-wideband range[J]. IET Radar, Sonar & Navigation, 2011, 5(4): 416-425.

[7] Karlsen B, Sorensen H B, Larsen J, et al.. Independent component analysis for clutter reduction in ground penetrating radar data[C]. International Society for Optics and Photonics, Orlando, USA, 2002: 378-389.

[8] Browne K E, Burkholder R J, and Volakis J L. Through-wall radar imaging system utilizing a light-weight low-profile printed array[C]. Proceedings of the Fourth European Conference on Antennas and Propagation (EuCAP), Barcelona, Spain, 2010: 1-5.

[9] Gu F, Chi L, Zhang Q, et al.. Single snapshot imaging method in multiple-input multiple-output radar with sparse antenna array[J]. IET Radar, Sonar & Navigation, 2013, 7(5): 535-543.

[10]Cui G, Kong L, and Yang X. Reconstruction filter design for stepped-frequency continuous wave[J]. IEEE Transactions on Signal Processing, 2012, 60(8): 4421-4426.

[11]Tivive F H C, Amin M G, and Bouzerdoum A. Wall clutter mitigation based on eigen-analysis in through- the-wall radar imaging[C]. IEEE 17th International Conference on Digital Signal Processing (DSP), Corfu, 2011: 1-8.

[12]Wehner D R. High Resolution Radar[M]. Norwood, MA, Artech House, 1987: 197-234.

[13]Jia Y, Kong L, Yang X, et al.. Multi-channel throughwall-radar imaging based on image fusion[C]. IEEE Radar Conference (RADAR), Kansas City, USA, 2011: 103-105.

[14]Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE Transactions on Neural Networks, 1999, 10(3): 626-634.

[15]Hyv?rinen A and Oja E. Independent component analysis: algorithms and applications[J]. Neural Networks, 2000, 13(4): 411-430.

[16]Bingham E and Hyv?rinen A. A fast fixed-point algorithm for independent component analysis of complex valued signals[J]. International Journal of Neural Systems, 2000, 10(1): 1-8.

[17]Du L, Liu H, Bao Z, et al.. A two-distribution compounded statistical model for radar HRRP target recognition[J]. IEEE Transactions on Signal Processing, 2006, 54(6): 2226-2238.

[18]Forootan E and Kusche J. Separation of deterministic signals using independent component analysis (ICA)[J]. Studia Geophysica et Geodaetica, 2013, 57(1): 17-26.

[19]Kirimoto T, Danial K, Amishima T, et al.. Performance analysis of independent component analysis to separate mixtures of complex sinusoidal signals[C]. IEEE SICE Annual Conference, Tokyo, 2008: 3424-3430.

[20]Xu Xiao-yin and Miller E L. Entropy optimized contrast stretch to enhance remote sensing imagery[C]. IEEE 16th International Conference on Pattern Recognition, Quebec, Canada, 2002: 915-918.

張 馳(1990–)，男，四川成都人，現(xiàn)為國(guó)防科技大學(xué)電子科學(xué)與工程學(xué)院碩士研究生，主要研究方向?yàn)榇Τ上窭走_(dá)信號(hào)處理。

E-mail: finderzc@aliyun.com

李?lèi)傷?1973–)，女，湖南瀏陽(yáng)人，博士，現(xiàn)為國(guó)防科技大學(xué)電子科學(xué)與工程學(xué)院副教授，獲國(guó)家科技進(jìn)步二等獎(jiǎng) 1項(xiàng)，部委科技進(jìn)步一、二等獎(jiǎng)各1項(xiàng)，主要研究方向?yàn)楹铣煽讖嚼走_(dá)成像以及實(shí)時(shí)信號(hào)處理系統(tǒng)實(shí)現(xiàn)。

周智敏(1957–)，男，山東濰坊人，現(xiàn)為國(guó)防科技大學(xué)電子科學(xué)與工程學(xué)院教授，博士生導(dǎo)師，主要研究方向?yàn)樾麦w制雷達(dá)系統(tǒng)與技術(shù)。

Wall Clutter Mitigation in Through-the-Wall Imaging Radar with Sparse Array Antenna Based on Independent Component Analysis

Zhang Chi Li Yue-li Zhou Zhi-min
(College of Electronic Science and Engineering, National University of Defense Technology, Changsha 410073, China)

For Through-the-Wall Imaging Radar (TWIR), wall clutter is critical for detecting target signals behind a wall. For a system with a sparse antenna array, the lack of observation channels makes it more difficult to separate the target signals and wall clutter. On the basis of fluctuation of the range profile in real transmit/receive channels, this paper proposes to use Independent Component Analysis (ICA) on multiple down-range observations of each transmit/receive channel to remove the wall clutter. The simulation and experimental results show that the proposed method effectively separate target and clutter components, even though the signal-to-clutter ratio is only -30 dB.

Through-the-Wall Imaging Radar (TWIR); Wall clutter mitigation; Independent Component Analysis (ICA); Sparse array antenna

TN951

A

2095-283X(2014)05-0524-09

10.3724/SP.J.1300.2014.14066

2014-04-04收到，2014-06-05改回；2014-09-01網(wǎng)絡(luò)優(yōu)先出版國(guó)家自然科學(xué)基金(61302146)資助課題

*通信作者： 張馳 finderzc@aliyun.com