X射線脈沖星光子序列頻域加權(quán)測相方法

焦榮,許錄平,張華,張嬌

(西安電子科技大學(xué)空間科學(xué)與技術(shù)學(xué)院,710126,西安)

?

X射線脈沖星光子序列頻域加權(quán)測相方法

焦榮,許錄平,張華,張嬌

(西安電子科技大學(xué)空間科學(xué)與技術(shù)學(xué)院,710126,西安)

針對X射線脈沖星信號累積輪廓位相測量誤差大的問題,提出了一種直接位相測量的頻域加權(quán)方法。首先,直接對X射線探測器捕獲和采樣形成的光子序列進(jìn)行頻域變換得到幅度譜,分別提取幅度譜各個成分的相位,再用幅度對相位變量加權(quán)平均,得到位相測量結(jié)果;進(jìn)一步采取分段加權(quán)的策略,解決序列過長引起的內(nèi)存溢出問題;最后對所提出方法的計算復(fù)雜度進(jìn)行了理論分析。仿真結(jié)果表明:在觀測時長和采樣間隔相同的條件下,與平均法相比,加權(quán)法的測相精度提高了10%以上;與最大似然(ML)和非線性最小均方誤差(NLS)法相比,加權(quán)法的比相精度提升了1倍以上,計算速度處于ML法和NLS法之間;分段情況下性能會隨著分段數(shù)增加有所改善。

X射線脈沖星;位相測量;快速傅里葉變換;最大似然

脈沖星是一種高速穩(wěn)定自旋的中子星,其旋轉(zhuǎn)軸和輻射軸不重合,當(dāng)探測器處于其輻射掃描區(qū)內(nèi)時,能夠收到穩(wěn)定的周期性的脈沖信號,利用該信號計時可以實現(xiàn)基于脈沖星的航天器導(dǎo)航。由于利用X射線信號有利于探測器小型化,因此X射線脈沖星導(dǎo)航已成為具有潛力的航天器導(dǎo)航方法,它通過測量脈沖星信號到達(dá)航天器和某一選定慣性參考點之間的相位差來獲取航天器相對于坐標(biāo)中心的位置信息[1-3]。其中,相位測量的誤差直接影響導(dǎo)航定位精度,因此相位測量技術(shù)是影響脈沖星導(dǎo)航性能的關(guān)鍵技術(shù)之一。

X射線脈沖星信號的測量技術(shù)大致可以分為兩類,一類是基于脈沖星信號累積輪廓的位相測量方法[1,4-5],另一類是直接對光子序列的初相進(jìn)行測量[3,6]。基于累積輪廓的位相測量方法又分為時域方法[7-8]、頻域方法和其他變換域方法,基于這幾種方法又延伸出了若干改進(jìn)方法如雙譜變換[9]和基于小波變換的輪廓測相等[3,6,10];光子序列初相測量以最大似然方法為主,根據(jù)模型的不同,也發(fā)展了一些改進(jìn)的方法[3,6]。從實際信號處理過程來看,基于時域的周期估計算法耗時長,頻域方法及其改進(jìn)方法[4-5,7-9]往往能在不顯著增加計算量的前提下獲得更高的精度,但這些方法一般都需要先累積出輪廓,而輪廓的累積過程受到脈沖星周期誤差、演化模型等多種因素的影響,這些因素會直接導(dǎo)致輪廓存在相移[11]和畸變,從而增加了測量誤差。因此,從位相測量的角度,直接對光子序列測相,可以避免輪廓累積的過處理引入的附加誤差,從而精度更高。已有文獻(xiàn)表明,針對光子序列的最大似然方法有能力獲得更高的精度[6]。Zhang等人提出的基于輪廓模型的最大似然方法[3]存在以下不足:①需對每個光子求解代價函數(shù)并搜索,計算量大;②最大似然對信號光子和噪聲光子的處理方法相同,對信號、噪聲以及不同的輪廓成分并不加區(qū)分。因此,本文提出基于快速傅里葉變換(FFT)的X射線脈沖星光子序列位相測量的頻域加權(quán)方法?；舅枷胧?通過頻域變換和加權(quán),使強度不同的成分對測量精度的貢獻(xiàn)不同。此外,由于背景輻射噪聲具有較寬、較平坦的頻譜,在頻譜加權(quán)方法中的權(quán)重較小,其對相位測量精度的影響能被一定程度的抑制。

1 信號模型

脈沖星的脈沖信號具有穩(wěn)定的周期,如果已知t0時刻的初相φ0、頻率f以及頻率的多階導(dǎo)數(shù)f(m),太陽系質(zhì)心(SSB)處的相位可以利用相位演化模型精確預(yù)測[12]

(1)

式中:φ(t)表示t時刻的相位;[·]表示取余數(shù)操作,定義為[φ]=int(φ)+mx,mx∈[0,1),int表示取整,[·]保證相位始終為[0,1)的實數(shù);M為整數(shù);通常,f(m)(m≥2)非常小(10-12～10-19),脈沖星在數(shù)小時甚至幾天內(nèi)可以近似認(rèn)為是等周期的,周期為1/f。由于脈沖星的輪廓結(jié)構(gòu)具有多樣性,因此在頻域可以看到明顯的譜線。本文選擇具有穩(wěn)定的周期特性的X射線脈沖星作為導(dǎo)航的參考基準(zhǔn)。

實際信號強度非常微弱,被認(rèn)為是以光子形態(tài)順序到達(dá)探測器的,再經(jīng)過X射線探測器捕獲和采樣形成光子序列,又由于到達(dá)SSB處脈沖星信號的周期特性,因此常用非齊次循環(huán)平穩(wěn)泊松過程對X射線脈沖星信號建模[13-14]。

2 頻域加權(quán)比相

2.1 信號獲取

X射線脈沖星信號的位相測量是脈沖星導(dǎo)航的基本測量量[15-18],由于導(dǎo)航方式的不同,位相測量的對象也有不同。在絕對形式的導(dǎo)航中,探測器測量到的脈沖星信號與標(biāo)準(zhǔn)信號演化模型進(jìn)行相位比對;在相對形式的脈沖星導(dǎo)航中,同一探測器在不同時段獲得的同一脈沖星信號之間,或者不同探測器在不同軌道位置探測的同一脈沖星信號之間進(jìn)行相位比對。在絕對導(dǎo)航方式中,由于任意時刻的標(biāo)準(zhǔn)信號演化模型是不變的,因此在不同時刻脈沖星信號與標(biāo)準(zhǔn)信號演化模型的比對,可以等效于探測器在不同時刻探測的信號之間的相位比對。所以,本文以相對形式導(dǎo)航為背景,對相位測量方法進(jìn)行闡述。

2.2 頻域加權(quán)比相

加權(quán)FFT相位估計的基本思路是:兩路光子序列來自同一顆脈沖星的信號,對兩列光子流量序列進(jìn)行FFT變換,對其每個頻域點求相位,再對兩列信號對應(yīng)頻點相位做差運算后加權(quán)平均,得到延遲相位估計值。令兩路比相序列分別為G1和G2,加權(quán)FFT比相方法示意圖如圖1所示。

圖1 加權(quán)FFT比相方法示意圖

兩路光子序列來自同一顆脈沖星,經(jīng)過等間隔采樣分別得到采樣后序列x1(n)、x2(n),其中n為整數(shù),n∈[0,N],N為序列長度,且有N=NbNp,Np為觀測時間內(nèi)的脈沖周期個數(shù),Nb代表單周期內(nèi)包含的采樣點數(shù)。設(shè)定x2(n)相比于x1(n)延遲相位τ0,即x2(n)=x1(n-τ0)+γ,其中γ為白噪聲,如果暫時不考慮整周期模糊度問題,有τ0∈[0,Nb],則有兩光子序列的N點FFT分別為

(2)

X2(k)=fft[x2(n)]=fft[x1(n-τ0)+γ]=

(3)

式中:fft[·]表示傅里葉變換;fft[γ]為常數(shù)。分別提取X1(k)、X2(k)對應(yīng)頻點相位值進(jìn)行做差運算,可以得到k頻點相位差為

(4)

式中:φ1(k)、φ2(k)分別為k頻點對應(yīng)的相位。容易得到歸一化的k頻點脈沖星的相位差

(5)

(6)

式中:w(k)為權(quán)值,其計算式為

(7)

(8)

(9)

2.3 分段式頻域加權(quán)比相

當(dāng)信號序列較長時,FFT變換可能導(dǎo)致計算機內(nèi)存溢出,因此本文提出分段式頻域加權(quán)比相方法。

設(shè)長序列{λ(1),λ(2),…,λ(N)}為單位時間片段Tb(周期內(nèi)最小采樣間隔時間長度)內(nèi)的光子數(shù)量。為便于FFT計算,將光子強度序列分為L段,令I(lǐng)=1,2,…,L,表示第I段,每段均包含相同的整數(shù)個周期長度,其中每段含M個樣點,即N=LM。若為了獲得適合FFT算法的長度,可對每段數(shù)據(jù)進(jìn)行補零運算。

相應(yīng)地,兩光子強度序列的相位差為

(10)

式中:φI,1(k)、φI,2(k)分別為序列1及序列2在第I段的相位。兩光子的延遲相位值為

(11)

令sI(k)=ΔφI(k)/k,則有

(12)

(13)

式中:k=1,2,…,M。平均所有數(shù)據(jù)段的相位延遲估計值,即可得到長序列的相位延遲估計值為

(14)

理論上,本文方法性能優(yōu)勢體現(xiàn)在:①從式(8)和式(10)可見,頻域方法不需要采樣時間間隔,其測量相位是連續(xù)的;②X射線背景噪聲流量恒定,通?？烧J(rèn)為服從均勻分布,當(dāng)光子數(shù)較多時,其幅度譜一般可近似為常值[14],信號的周期性使其譜線幅度明顯高于噪聲,因此本文通過引入幅度權(quán)重使噪聲對相位測量的影響得到抑制,相比而言,NLS和ML法對待噪聲和信號是等同的。

2.4 計算復(fù)雜度分析

計算復(fù)雜度是指對用于比相估計的算法所采用的加、減、乘、除的總運算次數(shù)的分析。假定觀測時間Tobs內(nèi)約有Np個脈沖星周期P,則有Tobs≈NpP。對該段觀測時間進(jìn)行等間隔采樣,則每間隔長度為Tb,Nb代表一個周期內(nèi)包含的時間片段數(shù),則有Tb=P/Nb。

采用加權(quán)FFT法進(jìn)行比相估計,需要對兩列光子序列進(jìn)行FFT變換,頻域內(nèi)的相位值作差,累積,最后對所得到的能量加權(quán)。因此,完成加權(quán)FFT比相估計總共需要(NbNp)2+6NbNp步運算過程。

分段加權(quán)FFT運算主要使在FFT變換處的計算代價得到了極大的節(jié)省。假定觀測時間內(nèi)的Np個周期被分成了L段,每段包含K(K?N)個周期,則分段加權(quán)FFT總共需要L(NbK)2+6NbNp步運算。

[6]中NLS方法和參考文獻(xiàn)[3]中ML方法的計算量與本文方法作比較,如表1所示。

表1 NLS、ML、加權(quán)FFT及分段加權(quán)FFT方法 的計算復(fù)雜度比較

注:Ng表示將歸一化相位區(qū)間(0,1)分為Ng段。

可見,加權(quán)FFT法的計算復(fù)雜度介于NLS法與ML法之間。此外,分段加權(quán)FFT法使一次FFT變換的實際數(shù)據(jù)量得到了極大的縮減,從而實現(xiàn)了計算復(fù)雜度的進(jìn)一步改善。

3 仿真結(jié)果與分析

本實驗采用RXTE探測器探測到的蟹狀脈沖星PSR B0531+21數(shù)據(jù)作為分析對象。B0531+21脈沖星的周期為33.4 ms,RXTE探測器的有效面積為6 800 cm2;X射線源發(fā)射的光子序列及探測器的背景噪聲主要由高斯模型和泊松模型擬合生成[14],該方法將脈沖星輪廓用高斯函數(shù)擬合,然后結(jié)合信號光子和噪聲密度計算探測器單位采樣間隔內(nèi)探測到的總流量強度,并作為泊松分布的參數(shù),進(jìn)而模擬生成光子序列。仿真中設(shè)置光子流量密度為0.96 ph/(cm2·s)(軟X射線能段),背景噪聲密度為0.45 ph/(cm2·s)。軟件平臺為Matlab R2012a,計算機配置為:Intel(R)Core(TM)i5-3470 CPU@3.20 GHz,內(nèi)存為4.00 GB(3.47 GB可用),操作系統(tǒng)為Window 7專業(yè)版。

3.1 均方根誤差對比分析

比相估計的均方根誤差主要用于反映比相估計值偏移真值的程度,其計算公式為

φ0)2])1/2

(15)

在測試過程中,設(shè)定觀測時間為0.1～1 000 s,采樣間隔Tb為P/1 024 s(P為脈沖星周期,此處為33.4 ms)。實驗過程中相位是按(0,1)歸一化。采用蒙特卡羅方法進(jìn)行多組隨機仿真,得到NLS、ML及加權(quán)FFT方法估計的比相均方根誤差隨觀測時間的變化曲線如圖2a所示。

由圖2a可以看出:隨著觀測時間的增加,這3種方法估計的比相均方根誤差均呈現(xiàn)逐漸減小的趨勢,其中NLS和ML算法估計的均方根誤差基本在相同的數(shù)量級水平變動,而加權(quán)FFT方法比前兩者估計出的均方根誤差要明顯低,即其比相精度要明顯高于NLS和ML算法。

(a)均方根誤差曲線

(b)CPU時間曲線圖2 NLS、ML及加權(quán)FFT法的比相估計對比

由于頻域內(nèi)兩光子強度序列間的相位差與對應(yīng)頻點的比值在每個頻率值處均對應(yīng)一個比值,這些比值十分接近,因此可反映出時域內(nèi)的相移量。為了準(zhǔn)確獲得時域內(nèi)的相移,需要對這些比值進(jìn)行處理,一般可通過平均法或加權(quán)法來得到一個穩(wěn)定的比值?；谄骄ê图訖?quán)法對這些比值處理后得到的相移誤差如圖3所示。加權(quán)法賦予有效光子信號更大的比值,而平均法則對有用信號和噪聲信號進(jìn)行了等比重處理,由圖3可見,加權(quán)法較平均法得到的測相誤差更小。

圖3 平均法和加權(quán)法的測相誤差對比圖

3.2 計算時間對比分析

NLS、ML及加權(quán)FFT方法比相估計的計算代價主要考察各方法仿真占用CPU時間。其中,NLS法的總CPU時間主要由輪廓累積、完成非線性最小均方估計及完成拋物插值運算所需的時間構(gòu)成;ML法的總CPU時間主要由最大似然函數(shù)模型搭建及插值運算所需的時間構(gòu)成;加權(quán)FFT法的總CPU時間主要是FFT變換和相位累積及加權(quán)所需的時間總和。圖2為比相估計所需的總CPU時間。由圖可知,隨著觀測時間的延長,NLS法使用CPU時間呈現(xiàn)緩慢增長的趨勢,ML和加權(quán)FFT法所用CPU時間均呈現(xiàn)顯著增長的趨勢,但ML法比加權(quán)FFT法使用的CPU時間要長約10倍左右。可見,在計算代價方面,隨著觀測時間的延長,加權(quán)FFT算法的計算代價介于NLS與ML算法之間,且相比于ML算法要明顯省時。

3.3 分段加權(quán)FFT與加權(quán)FFT方法

3.3.1 觀測時間及采樣間隔對加權(quán)FFT比相測量精度的影響 在不同的觀測時間(Tobs)和采樣間隔(Ts)條件下,進(jìn)一步驗證它們對加權(quán)FFT法比相測量精度的影響。實驗采取控制變量的方式及蒙特卡羅法對多組[0,0.5]范圍內(nèi)隨機相位進(jìn)行了測試。通過實驗,得到在不同觀測時間和采樣間隔條件下的加權(quán)FFT比相的均方根誤差分布圖,如圖4所示。從圖4可以看出:在一定的采樣間隔下,當(dāng)觀測時間逐漸延長時,加權(quán)FFT法的比相精度逐漸提高;在一定的觀測時間下,對觀測到的脈沖星序列進(jìn)行不同采樣間隔的采樣,當(dāng)采樣數(shù)逐漸增大時,加權(quán)FFT的比相精度隨之提高。可見,在一定的觀測時間和采樣間隔變化范圍內(nèi),脈沖信號的觀測時間越長,采樣間隔數(shù)越多,加權(quán)FFT算法估計出來的比相值越接近真實值,但當(dāng)觀測時間和采樣間隔足夠大時,加權(quán)FFT法比相精度不再明顯受這兩個變量的影響而逐漸趨于穩(wěn)定。

圖4 加權(quán)FFT法在不同觀測時間和采樣間隔下的均方根誤差二維分布圖

3.3.2 加權(quán)FFT法和分段加權(quán)FFT法比相的RMS、計算代價的對比分析 分段加權(quán)FFT法的比相精度仍然從均方根誤差和計算代價兩方面進(jìn)行測試及驗證。實驗中,分別以相同觀測時間內(nèi)32個周期、64個周期及256個周期為分段區(qū)間長度進(jìn)行測試,結(jié)果如圖5所示。

(a)均方根誤差曲線

(b)CPU時間曲線圖5 加權(quán)FFT法在不同分段長度下的比相估計對比圖

由圖5a可以觀察到:初始時刻,分段加權(quán)FFT法比相的均方根誤差幾乎與該小段時間內(nèi)進(jìn)行加權(quán)FFT法比相的均方根誤差接近,但隨著觀測時間的延長,均方根誤差會緩慢下降,雖然下降的速度明顯要低于加權(quán)FFT法的比相,總體上來說仍然高于ML和NLS法。由圖5b可以看出,隨著觀測時間的延長,分段加權(quán)FFT法的CPU時間仍然保持增長趨勢,但明顯要少于加權(quán)FFT法。此外,加權(quán)FFT法在140 s左右的觀測時間處出現(xiàn)了內(nèi)存溢出的現(xiàn)象,無法繼續(xù)進(jìn)行比相估計,而分段加權(quán)FFT法仍然可以實現(xiàn)比相,可有效解決內(nèi)存溢出問題。

4 結(jié) 論

針對脈沖星導(dǎo)航中的比相問題,本文提出了頻域加權(quán)方法。該方法考慮了X射線信號頻率成分對測相貢獻(xiàn)上的差異,對不同的頻率成分,根據(jù)強度賦予不同的權(quán)值,適用于基于脈沖星的相對導(dǎo)航以及相位增量測量等應(yīng)用。仿真結(jié)果表明,本文的頻域加權(quán)法在測相精度上優(yōu)于平均法、NLS和ML方法,計算復(fù)雜度介于NLS和ML方法之間。分段頻域加權(quán)后,性能有所下降,但可以克服序列過長引起的內(nèi)存溢出問題,從而可以適應(yīng)內(nèi)存受限的計算系統(tǒng)。

參考文獻(xiàn):

[1] SHEIKH S I. The use of variable celestial X-ray sources for spacecraft navigation [D]. Maryland, MD, USA: University of Maryland, College Park, 2005.

[2] LIU Jing, MA Jie, TIAN Jinwen, et al. X-ray pulsar navigation method for spacecraft with pulsar direction error [J]. Advances in Space Research, 2010, 46(11): 1409-1417.

[3] ZHANG Hua, XU Luping, SHEN Yanghe, et al. A new maximum-likelihood phase estimation method for X-ray pulsar signals [J]. Journal of Zhejiang University: Science C, 2014, 15(6): 458-469.

[4] RINAURO S, COLONNE S, SCARANO G. Fast near-maximum likelihood phase estimation of X-ray pulsars [J]. Signal Processing, 2013, 93(1): 326-331.

[5] ZHANG Hua, XU Luping. An improved phase measurement method of integrated pulse profile for pulsar [J]. Science China Technological Sciences, 2011, 54(9): 2263-2270.

[6] EMADZADEH A A, SPEYER J L. On modeling and pulse phase estimation of X-ray pulsars [J]. IEEE Transactions on Signal Processing, 2010, 58(9): 4484-4495.

[7] 李建勛, 柯熙政. 基于循環(huán)平穩(wěn)信號相干統(tǒng)計量的脈沖星周期估計新方法 [J]. 物理學(xué)報, 2010, 59(11): 8304-8310. LI Jianxun, KE Xizheng. Period estimation method for weak pulsars based on coherent statistic of cyclostationary signal [J]. Acta Physic Sinica, 2010, 59(11): 8304-8310.

[8] 李建勛, 柯熙政, 趙寶升, 等. 一種脈沖星周期的時域估計新方法 [J]. 物理學(xué)報, 2012, 61(6): 069701. LI Jianxun, KE Xizheng, ZHAO Baosheng, et al. A new time-domain estimation method for period of pulsars [J]. Acta Physica Sinica, 2012, 61(6): 069701.

[9] 謝振華, 許錄平, 倪廣仁, 等. 基于雙譜的脈沖星累積脈沖輪廓時間延遲測量 [J]. 物理學(xué)報, 2008, 57(10): 6683-6688. XIE Zhenhua, XU Luping, NI Guangren, et al. Time offset measurement algorithm based on bispectrum for pulsar integrated pulse profiles [J]. Acta Physic Sinica, 2008, 57(10): 6683-6688.

[10]閻迪, 許錄平, 謝振華, 等. 脈沖星信號模糊閾值小波降噪算法 [J]. 西安交通大學(xué)學(xué)報, 2007, 10(41): 1193-1196. YAN Di, XU Luping, XIE Zhenhua, et al. Wavelet denoising algorithm based on fuzzy threshold for pulsar signal [J]. Journal of Xi’an Jiaotong University, 2007, 41(10): 1193-1196.

[11]ZHANG Hua, XU Luping, XIE Qiang, et al. Modeling and Doppler measurement of X-ray pulsar [J]. Science China Physics, Mechanics and Astronomy, 2011, 54(6): 1068-1076.

[12]TAYLOR J H. Pulsar timing and relativistic gravity [J]. Philosophical Transactions of the Royal society of London: Series A Physical and Engineering Sciences, 1992, 341(1660): 117-134.

[13]SALA J, URRUELA A, VILLARES X, et al. Feasibility study for a spacecraft navigation system relying on pulsar timing information: ARIADNA Study 03/4202 [R]. Barcelona, Spain: European Space Agency Advanced Concepts Team, 2004.

[14]ZHANG Hua, XU Luping, SONG Shibin, et al. A fast method for X-ray pulsar signal simulation [J]. Acta Astronautica, 2014, 98(5): 189-200.

[15]SHEIKH S I, PINES D J, RAY P S, et al. Spacecraft navigation using X-ray pulsars [J]. Journal of Guidance Control and Dynamics, 2006, 29(1): 49-63.

[16]XIONG Zhi, QIAO Li, LIU Jianye, et al. GEO satellite autonomous navigation using X-ray pulsar navigation and GNSS measurements [J]. International Journal of Innovative Computing, Information and Control, 2012, 8(5): 2965-2977.

[17]BERNHARDT M G, PRINZ T, BECKER W, et al. Timing X-ray pulsars with application to spacecraft navigation [C]∥Proceedings of High Time Resolution Astrophysics IV. Garching, Germany: Institute of Astronautics, 2010: 1-5.

[18]ASHBY N, HOWE D A. Relativity and Timing in X-ray pulsar navigation [C]∥Proceedings of 2006 IEEE International Frequency Control Symposium and Exposition. Piscataway, NJ, USA: IEEE, 2006: 767-770.

(編輯 劉楊)

A Frequency Domain Weighted Measurement Method for X-Ray Pulsar Photon Sequences

JIAO Rong,XU Luping,ZHANG Hua,ZHANG Jiao

(School of Aerospace Science and Technology, Xidian University, Xi’an 710126, China)

A weighted direct phase estimation algorithm in frequency domain for photon sequence is proposed to solve the problem that the error of accumulated phase estimation in the average profile is too large for X-ray pulsar signals. The photon sequence that is captured by the X-ray detector and then sampled is directly transformed to the frequency domain, then the phase of each frequency is respectively extracted and is weighted according the amplitude of the frequency component. Furthermore, the segmented fast Fourier transform strategy is adopted to solve the memory overflow problem caused by long sequences. Finally, the computational complexity of the proposed estimator is investigated in theory. Simulations and comparisons with existing methods under the condition of the same observation time and sampling interval show that the accuracy of the weighted approach increases by 10% over the average method and is two times higher than those of the nonlinear least square (NLS) and maximum likelihood (ML) methods, and that the computational complexity of the algorithm is between NLS and ML methods. Moreover, the performance of the algorithm with segments is improved by increasing the number of segments.

X-ray pulsar; phase estimation; fast Fourier transform; maximum likelihood

2015-12-17。 作者簡介:焦榮(1983—),女,博士生;許錄平(通信作者),男,教授,博士生導(dǎo)師。 基金項目:國家自然科學(xué)基金資助項目(61172138,61401340)。

時間:2016-04-03

10.7652/xjtuxb201606023

P128.4

A

0253-987X(2016)06-0152-07

網(wǎng)絡(luò)出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20160403.1820.004.html