基于S變換的機載雷達穩健空時自適應算法

高志奇, 陶海紅, 趙繼超

(西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071)

?

基于S變換的機載雷達穩健空時自適應算法

高志奇, 陶海紅, 趙繼超

(西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071)

摘要：非均勻雜波環境中的離群點會導致空時自適應處理(space-time adaptive processing，STAP)性能的下降。針對此問題,該文提出一種新的基于S變換的穩健STAP(S transform-STAP, ST-STAP)算法。該算法主要利用雜波和離群點的S變換在時頻平面上分布特性的不同來實現非均勻雜波的抑制。ST-STAP算法首先將樣本數據對應的快時間序列作S變換得到時頻分布矩陣,舍棄部分高頻分量以去掉離群點的影響,然后用時頻矩陣的剩余部分計算相關矩陣和自適應濾波權。蒙特卡羅實驗證明ST-STAP方法的穩健性和動目標檢測性能均優于傳統的STAP算法。

關鍵詞：機載雷達; 空時自適應; S變換; 離群點

0引言

目前,空時自適應處理(space-time adaptive proces-sing, STAP)是機載雷達雜波抑制和動目標檢測的基本方法[1]。STAP的最優濾波權依賴于雜波協方差矩陣R的精確估計。根據RMB準則,若想通過樣本數據精確估計R,必須滿足以下條件:①所有樣本數據必須滿足獨立同分布(independent and identically distributed, IID)假設;②用于估計R的IID樣本數據的數量至少要達到系統自由度的兩倍才能保證系統的輸出信雜噪比損失不大于3 dB;③樣本數據中不能包含動目標等雜質信息[2]。

STAP方法通常是基于均勻雜波環境,所以R可以通過均勻的訓練樣本數據進行估計。但實際的雜波環境是非均勻的,特別是訓練樣本中包含有動目標信號時,R的估計會出現較大誤差,從而導致STAP性能的極大下降。包含在訓練樣本中的動目標信號稱為離群點。離群點導致訓練樣本不滿足IID假設的情況稱為非均勻雜波[3]。為了克服非均勻雜波對于動目標檢測性能的影響,人們一般使用對角加載(diagonal loading, DL)和非均勻樣本檢測(non-homogeneity detection, NHD)等方法來增強STAP的穩健性[4]。DL方法是通過在樣本數據相關矩陣中加入一定功率的噪聲信號從而抑制離群點的干擾。如何選擇合適的加載功率依然是該方法的一個難點[5],實際應用中一般根據經驗進行選擇,所以該方法由于很難獲得合適的加載功率而不能達到最優濾波效果。NHD一般是通過廣義內積(generalized inner products, GIP)等方法從訓練樣本集中挑選并剔除包含離群點的樣本數據[6],但是當離群點較多時,GIP方法的性能下降很大[7]。

本文提出一種基于S變換(S transform, ST)的穩健STAP算法(ST-STAP)用于抑制離群點的影響。該算法主要利用雜波與離群點在時頻平面分布不同的特點,通過S變換將快時間序列在時頻平面上展開,舍棄離群點頻譜較為集中的頻率段,然后計算自相關矩陣,以消除離群點對于濾波權的影響,提高STAP的穩健性。

1問題模型

1.1信號模型[8-10]

圖1　機載雷達工作原理圖

圖2　雷達采樣數據集

每個距離單元的雜波數據可以看作Nc個雜波反射點的反射信號之和。第i個雜波反射點對應的空間頻率ωs,i和歸一化多普勒頻率ωd,i分別定義為

(1)

(2)

式中,i=1,2,…,Nc;d為陣元間距。第i個雜波反射點的空時導向矢量表示為

(3)

式中

(4)

(5)

分別表示時間導向矢量和空間導向矢量;?表示Kronecker積;(·)T表示矩陣轉置。第k個距離單元的雜波加噪聲快拍數據表示為

(6)

式中,εc,i表示第i個雜波反射點的復系數;N0代表噪聲。

1.2非均勻樣本問題[11-13]

傳統STAP方法的最優濾波權由式(7)決定

(7)

式中,S為目標的導向矢量;R為理想的雜波協方差矩陣。實際上R是未知的,只能進行估計。通過樣本數據估計的雜波協方差矩陣Re表示為

(8)

式中,k為待檢測距離單元;(·)H表示共軛轉置。通常用Re代替式(7)中的R來計算最優濾波權。上述通過樣本數據估計協方差矩陣求逆(samples matrix inverse, SMI)計算最優濾波權的方法記為SMI-STAP。

STAP的訓練樣本數據必須滿足IID假設,即要求被觀測地形的特性保持不變且不包含動目標。但是某些特殊的被觀測地形(如高速公路、城市)在不同的距離單元內會同時出現動目標。如果將包含動目標的采樣數據作為訓練樣本,一旦樣本與被檢測距離單元二者中包含的動目標多普勒頻率接近,濾波器會將被檢測的動目標信號當作雜波濾除而出現動目標信號相消現象。為克服非均勻雜波環境的影響,需要對Re進行改進。DL和NHD是改進Re的常用方法。DL是在Re的對角線上加入一個非負常數,然后對其求逆計算最優濾波權,稱為DL-STAP算法。NHD是通過某種規則將包含離群點的訓練樣本挑選出來并進行剔除的一種方法。廣義內積(generalizedinnerproducts,GIP)檢測法是典型的NHD方法。該方法由樣本數據產生GIP序列,并通過一定的門限選擇和剔除包含離群點的樣本,然后利用剩余樣本估計Re并求逆計算最優濾波權。該方法記為GIP-STAP。

2ST-STAP算法

2.1S變換[14-15]

對任意函數x(t)∈L2(R),其ST定義為

(9)

式中,τ為控制高斯窗口在t軸位置的參數。由定義可以看出,ST相當于信號加了寬度和高度可變的高斯窗函數后進行的傅里葉變換,因此其時間和頻率分辨率都是可變的。

S逆變換的定義為

(10)

由式(9)和式(10)可以看出，ST與傅里葉變換在本質上具有一致性,因此可以借助快速傅里葉變換實現ST的快速計算。實際中使用的都是離散形式的ST,其定義為

(11)

(12)

其中

(13)

D點離散信號x(n)的ST結果是列和行分別對應采樣時間點和頻率點的矩陣,其形式為

(14)

式中,V代表頻率采樣點數；D代表時間采樣點數。此處,S的列按照時間增加的順序排列,而行按照頻率遞減的方式排列。

2.2ST-STAP算法

對于圖2所示的雷達采樣數據集,每個距離單元對應的N×M維數據矩陣按列的順序排列為NM×1維的列向量,其中的元素表示某個陣元對應于特定發射脈沖時的接收信號。將所有列向量按距離單元的順序排列為NM×L維的樣本數據矩陣X,其形式為

(15)

X的每行表示特定的空時陣元隨距離單元(快時間)的變化特性,稱為快時間序列。該序列的大多數元素代表的是雜波加噪聲(其中噪聲是高斯分布的加性白噪聲)。在某個或某幾個快時間處,序列元素代表的是雜波、離群點與噪聲之和,即離群點在快時間序列中是孤立存在的。為便于分析,某個空時陣元對應的快時間序列可以近似表示為

(16)

式中,αli和φli(t)分別為第l個距離單元中的第i個雜波的幅度和相位;βj和γj(t)表示某個距離單元中包含的離群點幅度和相位;tl表示第l個距離單元對應的時間;tj表示包含離群點的第j個距離單元對應的時間;P表示包含離群點的距離單元數目;δ代表沖激函數;N0代表噪聲。式(16)表示每個距離單元都包含Nc個雜波,每個雜波都具有不同的幅度和相位;某些特定距離單元中各自包含一個離群點;噪聲按某種規律分布于所有距離單元中。

忽略噪聲的條件下,快時間序列x(t)的ST結果為

(17)

由式(17)可知,某特定距離單元上的雜波和離群點ST的結果在時間維上的分布沒有明顯變化,各個距離單元的數據在時間維上不會相互影響。對第l個距離單元,如果其中包含離群點,則離群點對應的幅度頻譜可以表示為

(18)

由式(18)可知,第l個距離單元在只包含一個離群點的情況下,其ST頻譜中的相位[γ(tl)-2πftl]不會對幅度頻譜造成影響,因此離群點的幅度頻譜正比于頻率的絕對值。第l個距離單元中雜波對應的幅度頻譜可以表示為

(19)

式(19)中每個雜波的ST頻譜都具有隨機性相位[φ(tli)-2πftli],該相位使得在每個頻率點上的Nc個雜波頻譜都為隨機的復數,因此它們的疊加并不一定會使總的幅度頻譜增大,即

(20)

所以雜波的幅度頻譜會呈現隨機分布的特性。根據上述離群點和雜波的ST頻譜分布特性,將快時間序列ST結果的部分高頻成分舍棄可以減輕甚至去除離群點對于STAP的影響。另一方面,舍棄高頻成分顯然也會對雜波抑制性能造成影響,因此該算法需要在消除離群點的影響和雜波抑制性能二者之間進行折中考慮。

假設NM×L維的樣本數據矩陣第i行的快時間序列xi(t)的ST為Si(f,t),其結構如同式(14)。它是K×L維的矩陣,K為頻率點個數,L為快時間個數,即行代表頻率而列代表時間。由于Si(f,t)的各行按照頻率由高到低的順序排列,所以舍棄高頻分量即將Si(f,t)的前J行舍去,以減弱離群點的影響,則剩余部分構成一個(K-J)×L維的矩陣Fi(f,t)。所有快時間序列進行相同的處理,然后分別提取每個Fi(f,t)對應相同的列,按如下方式構成L個(K-J)×NM維的快拍時頻矩陣為

(21)

計算其時頻相關矩陣為

(22)

通過式(22)的相關運算,各快拍時頻矩陣Ui中包含的頻率信息已經被抵消,所以時頻相關矩陣Rt中只包含時間和空間信息,其維數為NM×NM。文獻[16]已經證明,Rt與直接由式(8)估計得到的樣本數據協方差矩陣具有相同的信號子空間和噪聲子空間,因此可以用Rt代替式(7)中的R來計算STAP的最優濾波權Wopt。由于在Ui中已經舍棄了離群點頻譜集中的高頻分量,所以Rt受離群點的影響很小。因此,由Rt得到的Wopt可以在抑制雜波的同時避免出現信號相消現象。按照上述過程計算Rt和Wopt的算法稱為ST-STAP算法。

顯然,Fi(f,t)中還有部分離群點頻譜分量剩余,但其中雜波的頻譜是主要成分。通過式(22)的相關運算之后,Rt中積累得到的雜波能量遠遠超過離群點的能量,因此可以達到抑制離群點的目的。另外,實際雜波的ST頻譜在頻率維呈現隨機分布,而且影響雜波幅度和相位的因素較多,因此目前尚無法用數學表達式精確描述雜波的ST幅度頻譜分布規律,只能通過仿真實驗進行呈現。同時,上述處理過程需要舍棄的高頻成分數目J很難通過公式推導得出,本文通過多次實驗進行分析和選擇。第3節將對比分析J的不同取值對于ST-STAP算法性能的影響。

3性能分析

本節將通過實驗仿真分析ST-STAP算法的性能。主要仿真參數為:N=10; M=10; λ=0.3m; H=3 000m; va=300m/s; fr=4 000Hz; 雜噪比為60dB;離群點的信噪比為40dB。以下仿真結果均為100次蒙特卡羅實驗的平均結果。

3.1雜波與離群點ST頻譜分析

仿真實驗選取了101～300號距離單元的采樣數據,在第130、180、230和260號距離單元內各包含一個離群點。采樣數據構成形同式(15)的數據矩陣,從中隨機抽取某一行作為仿真用的快時間序列,其ST時頻分布如圖3所示。為了體現離群點、雜波和噪聲的幅度頻譜分布特性的不同,將它們的幅度頻譜分別繪制于圖4～圖6中。

圖3　快時間序列的ST譜

圖4　離群點的ST譜

圖5　雜波的ST譜

圖6　噪聲的ST譜

可以看出,離群點的幅度頻譜只分布在對應的距離單元中,并按照正比例關系沿頻率維分布;而雜波的幅度頻譜在時頻平面上呈現隨機分布,其在頻率維的分布并非集中于高頻段;噪聲的幅度頻譜非常微弱,在后續的分析中可以忽略。為了進一步對比雜波和離群點的頻譜分布,分別提取第130、180、230和260號距離單元對應的雜波頻譜和離群點頻譜,結果如圖7～圖10所示。如果將歸一化頻率的0～0.5看作低頻段,而0.5～1看作高頻段,則分別計算上述4個距離單元中雜波的頻域能量比重,結果如表1所示。可以看出,低頻段的雜波能量比重超出高頻段一倍左右。比較而言,離群點的頻譜在高頻段較為集中,而雜波的頻譜分布較為均勻。因此，實驗的結果與式(18)和式(19)推導出的結果相符。

圖7　130號距離單元的雜波與離群點頻譜

圖8　180號距離單元的雜波與離群點頻譜

圖9　230號距離單元的雜波與離群點頻譜

圖10　260號距離單元的雜波與離群點頻譜

距離單元低頻段高頻段第130號0.65190.3081第180號0.65300.1950第230號0.61320.2990第260號0.68990.2369

3.2舍棄高頻分量數目J的選擇

ST-STAP算法中非常關鍵的一個問題是快時間序列的ST中需要舍棄的高頻分量數目J的選擇。根據上述分析,J很難通過公式推導得出,這里先依據多次實驗選取。首先考慮樣本數據集中包含單個離群點的情況。仿真實驗選取200個距離單元的采樣數據,其中包含一個離群點。傳統STAP算法的歸一化改善因子(normalized improvement factor, NIF)在雜波對應的多普勒頻率處形成零陷,也可能在離群點對應多普勒頻率處形成零陷。根據穩健STAP的要求,NIF在雜波多普勒頻率處的值越小越好,在其他位置的值應該為0 dB。圖11為雜波和離群點對應多普勒頻率位置的NIF值與J之間的對應關系。

圖11　單離群點時J對NIF的影響

由圖11可以看出,在舍棄了約10個高頻點之后,離群點對應多普勒頻率位置不再產生零陷,而雜波多普勒頻率位置的NIF約為-75 dB,此時算法的雜波抑制性能較好。

當樣本數據集中包含多個離群點時,主要考慮兩種情況:一是多個離群點的多普勒頻率相同;二是離群點的多普勒頻率各不相同。對于第一種情況,考慮200個距離單元中包含4個多普勒頻率相同的離群點,雜波和離群點對應多普勒頻率位置的NIF值與J之間的對應關系如圖12所示。

圖12　多個頻率相同離群點時J對NIF的影響

由圖12可以看出,由于離群點數目的增多,其對應多普勒頻率處的最大零陷加深了約20 dB,大約舍棄了30個高頻點之后才能保證將動目標的影響完全消除,而此時雜波多普勒頻率位置的NIF約為-65 dB,說明算法對雜波的抑制性能依然較好。對于第2種情況,考慮200個距離單元中包含4個多普勒頻率各不相同的離群點,雜波和離群點對應多普勒頻率位置的NIF值與J之間的對應關系如圖13所示。

圖13　多個頻率不同離群點時J對NIF的影響

由圖13可以看出,4個不同多普勒頻率的離群點對應產生了4個零陷,舍棄大約30個高頻點之后,離群點的影響被完全消除而雜波的抑制性能可以達到-70 dB。

綜合以上情況,ST-STAP算法中需要舍棄的高頻點數目與離群點的個數和多普勒頻率有關。當離群點的數目較多時,需要舍棄的高頻分量較多,則算法的雜波抑制性能會有一定的下降。J的選擇需要在抑制離群點的影響和保證雜波抑制能力之間進行折中。由上述仿真可知,J分別取10和30時,算法的雜波抑制性能差別不超過10 dB,因此后續所有關于ST-STAP算法性能的分析中,每個快時間序列的ST時頻矩陣均舍棄30個高頻點。

3.3快拍數對算法性能的影響

ST的頻率分辨率與序列長度成正比,而ST-STAP算法需要舍棄快時間序列ST中的若干高頻點,因此快拍數目會直接影響算法的性能。這里分析快拍數目為20～200時系統的性能。圖14為樣本數據包含一個離群點時雜波和離群點多普勒頻率位置的NIF隨快拍數目的變化。圖15為樣本數據包含4個相同多普勒頻率的離群點時相關多普勒頻率位置的NIF隨快拍數目的變化。圖16為樣本數據包含4個不同多普勒頻率的離群點時相關多普勒頻率位置的NIF隨快拍數目的變化。可以看出,快拍數小于40的情況下,ST-STAP算法無法克服離群點的影響。而隨著快拍數的增多,雜波的抑制性能增強,同時在舍棄30個高頻點的情況下,ST-STAP算法可以完全克服多個離群點的影響。

圖14　單離群點時快拍數對NIF的影響

圖15　多個頻率相同離群點時快拍數對NIF的影響

圖16　多個頻率不同離群點時快拍數對NIF的影響

3.4ST-STAP算法的穩健性

仿真首先采用101～300號距離單元的采樣數據作為訓練樣本,其中包含有4個離群點。離群點的歸一化多普勒頻率范圍為0.62～0.72。ST-STAP、SMI-STAP、DL-STAP和GIP-STAP算法分別利用該訓練樣本形成最優濾波權。圖17為4種STAP算法的NIF與多普勒頻率的關系曲線。由于多個強離群點的干擾,SMI-STAP、GIP-STAP和DL-STAP均在離群點多普勒頻率位置產生了較深的零陷,表示其對于動目標信號也具有一定程度的抑制作用;而ST-STAP只在主雜波區形成深度為80 dB的零陷,體現出很好的雜波抑制性能。

圖17　NIF與多普勒頻率的關系

為了檢測算法的動目標檢測性能,301～400號距離單元的采樣數據作為測試樣本。在第350號距離單元處注入一個動目標信號,其信噪比為10 dB,歸一化多普勒頻率為0.65。分別使用上述4種STAP 算法形成的最優濾波權對測試樣本進行濾波,得到的輸出功率如圖18所示。可以看到,由于離群點導致的目標相消,SMI-STAP算法未檢測出動目標信號;ST-STAP算法的動目標輸出功率比DL-STAP和GIP-STAP高約15 dB。因此ST-STAP算法具有更好的動目標檢測性能。

圖18　輸出功率

4結論

針對訓練樣本數據中包含離群點而造成雜波環境非均勻的情況,本文研究了基于S變換的穩健STAP算法。ST-STAP算法利用了時頻分布中雜波的頻譜為隨機分布、而離群點的頻譜集中在高頻段的特點,先舍棄快時間序列ST時頻矩陣的部分高頻分量,然后計算其時頻相關矩陣,則可以有效消除離群點的影響。ST-STAP算法在訓練樣本中包含多個強離群點的情況下,形成的最優濾波權對測試樣本濾波后的動目標輸出功率比其余方法大,動目標檢測的準確率較高。因此ST-STAP算法與其他幾種方法相比穩健性更好,更適合于實際的雜波環境。

參考文獻：

[1] Zhang B H, Zhang Y S, Wang Y L. Clutter suppression method based on estimation of inverse covariance matrixes[J].SystemsEngineeringandElectronics,2011,33(5):1027-1031.(張柏華,張永順,王永良.基于逆協方差矩陣估計的雜波抑制方法[J].系統工程與電子技術,2011,33(5):1027-1031.)

[2] Reed I S, Mallet J D, Brennan L E. Rapid convergence rate in adaptive arrays[J].IEEETrans.onAerospaceandElectronicSystems, 1974, 10(6): 853-863.

[3] Liu C F, Liao G S. Robust generalized inner products algorithm using diagonal loading[J].JournalofXidianUniversity, 2009, 36(4): 583-589. (劉聰鋒,廖桂生.穩健的對角加載廣義內積算法[J].西安電子科技大學學報, 2009, 36(4): 583-589.)

[4] Cheng X J. Direct data domain space time adaptive processing algorithm for non-homogeneous clutter[J].ComputerEngineeringandApplications, 2015, 51(5): 217-221. (程向嬌.非均勻雜波中的直接變換域STAP算法研究[J].計算機工程與應用, 2015, 51(5): 217-221.)

[5] Cao J S, Wang X G. Diagonally loaded SMI algorithm based on inverse matrix recursion[J].SystemsEngineeringandElectronics, 2007, 18(1): 160-163.

[6] Zhou Y, Zhang L R, Liu N, et al. Knowledge aided secondary data selection in space time adaptive processing[J].SystemsEngineeringandElectronics,2010,32(2):405-409.(周宇,張林讓,劉楠,等.空時自適應處理中基于知識的訓練樣本選擇策略[J].系統工程與電子技術,2010,32(2):405-409.)

[7] Tong Y L, Wang T, Dai B Q, et al. Heterogeneous clutter suppression method based on the knowledge of clutter ridge[J].SystemsEngineeringandElectronics, 2015, 37(5): 1035-1041. (同亞龍,王彤,代保全,等.基于雜波脊先驗信息的非均勻雜波抑制方法[J].系統工程與電子技術,2015, 37(5): 1035-1041.)

[8] Yang Z C, Li X, Wang H Q, et al. Adaptive clutter suppression based on iterative adaptive approach for airborne radar[J].SignalProcessing, 2013, 93: 3567-3577.

[9] Wang Y L, Liu W J, Xie W C. Research progress of space-time adaptive detection for airborne radars[J].JournalofRadars, 2014,3(2):201-207.(王永良,劉維建,謝文沖.機載雷達空時自適應檢測方法研究進展[J].雷達學報,2014,3(2):201-207.)

[10] Xiang C, Feng D Z, Cao Y, et al. Three-dimensional spatial-temporal dimension-reduced adaptive processing for airborne radar[J].SystemsEngineeringandElectronics, 2011, 33(3): 548-552.(向聰,馮大政,曹楊,等.機載雷達三維空時降維自適應處理[J].系統工程與電子技術,2011,33(3):548-552.)

[11] Fan H Y, Jiang Y M, Kuang G Y. Target detection in non-stationary clutter background and Riemannian geometry[J].IETRadar,SonarandNavigation, 2014, 8(4): 376-381.

[12] Zhuang J, Manikas A. Interference cancellation beamforming robust to pointing errors[J].IETSignalProcessing, 2013, 7(2): 120-127.

[13] Yang X P, Liu Y X, Long T. Robust non-homogeneity detection algorithm based on prolate spheroidal wave functions for space-time adaptive processing[J].IETRadar,SonarandNavigation,2013,7(1):47-54.

[14] Zhu M Z, Ji H B, Lin L. Instantaneous frequency estimation of multi-component FM signals based on directional S transform[J].SystemsEngineeringandElectronics, 2013, 35(1): 29-33. (朱明哲,姬紅兵,林琳.基于方向性S變換的多分量FM信號瞬時頻率估計[J].系統工程與電子技術,2013,35(1):29-33.)

[15] Sun Z S, Wang J, Bi Y X, et al. Analysis of multi-human micro-Doppler signatures based on generalized S transform[J].SystemsEngineeringandElectronics,2014,36(7):1291-1297.(孫忠勝, 王俊, 畢嚴先, 等.基于廣義S變換的多人微多普勒特征分析[J].系統工程與電子技術,2014,36(7):1291-1297.)

[16] Zhang Y M, Mu W F, Amin M G. Subspace analysis of spatial time-frequency distribution matrices[J].IEEETrans.onSignalProcessing, 2001, 49(4): 747-759.

高志奇(1980-),男,講師,博士研究生,主要研究方向為空時自適應信號處理。

E-mail:gzqnd@163.com

陶海紅(1976-),女,教授,博士,主要研究方向為陣列信號處理。

E-mail:hhtao@xidian.edu.cn

趙繼超(1989-),男,博士研究生,主要研究方向為陣列信號DOA估計。

E-mail:13891894124@139.com

Robust space-time adaptive processing based on S transform for airborne radar

GAO Zhi-qi, TAO Hai-hong, ZHAO Ji-chao

(NationalLabofRadarSignalProcessing,XidianUniversity,Xi’an710071,China)

Abstract:Performance of space-time adaptive processing (STAP) is always corrupted by outliers. Hence, a novel robust STAP algorithm based on the S transform (ST) is proposed. It exploits the characteristic that ST spectrum distribution of the outlier is different from the clutter in the time-frequency domain. Firstly, the fast time sequence is transferred to the time-frequency matrix by ST. Some high frequency points of the ST matrix are abandoned to eliminate outliers. Then the remainder of the ST matrix is used to estimate the clutter covariance matrix and filter weight. Monte Carlo experiments prove that the proposed algorithm is more robust than other conventional STAP algorithms in non-homogeneous clutter environments.

Keywords:airborne radar; space-time adaptive processing (STAP); S transform (ST); outlier

收稿日期：2015-03-25；修回日期：2015-09-23；網絡優先出版日期：2016-02-05。

基金項目：國家自然科學基金(60971108)；國家重點基礎研究發展計劃(973計劃)(2011CB707001)資助課題

中圖分類號：TN 958

文獻標志碼：A

DOI:10.3969/j.issn.1001-506X.2016.06.08

作者簡介：

網絡優先出版地址：http://www.cnki.net/kcms/detail/11.2422.TN.20160205.1429.002.html