郭付陽 張子敬 楊林森

(西安電子科技大學(xué) 雷達(dá)信號(hào)處理國家重點(diǎn)實(shí)驗(yàn)室,西安 710071)

基于尺度變換的寬帶線性調(diào)頻信號(hào)時(shí)差/尺度差估計(jì)算法

郭付陽 張子敬 楊林森

(西安電子科技大學(xué) 雷達(dá)信號(hào)處理國家重點(diǎn)實(shí)驗(yàn)室,西安 710071)

提出了一種基于尺度變換的寬帶線性調(diào)頻信號(hào)時(shí)差和尺度差的快速算法. 根據(jù)兩路接收到的線性調(diào)頻信號(hào)間調(diào)頻率之比為尺度差的平方的特點(diǎn),利用分?jǐn)?shù)階傅里葉變換分別估計(jì)出兩路信號(hào)的調(diào)頻率,即可獲得尺度差的估計(jì). 將估計(jì)的尺度差對(duì)一路信號(hào)進(jìn)行伸縮,并計(jì)算伸縮后信號(hào)與另一接收信號(hào)的時(shí)域相關(guān),根據(jù)相關(guān)峰的位置估計(jì)出時(shí)差. 相比于傳統(tǒng)基于寬帶互模糊函數(shù)的方法,該方法避免了二維搜索寬帶互模糊函數(shù)的峰值,只需若干次快速傅里葉變換即可實(shí)現(xiàn),能夠顯著降低運(yùn)算量. 仿真結(jié)果顯示該方法在高信噪比下逐漸接近克拉美-羅下界.

分?jǐn)?shù)階傅里葉變換;寬帶互模糊函數(shù);時(shí)/尺度差估計(jì);寬帶線性調(diào)頻信號(hào); 克拉美-羅下界

引 言

在無源定位中,一種常用的定位方法是利用多個(gè)位于不同位置的接收機(jī)來接收輻射源信號(hào),通過估計(jì)各接收信號(hào)間的時(shí)差來實(shí)現(xiàn)定位. 當(dāng)接收機(jī)和輻射源存在相對(duì)運(yùn)動(dòng)時(shí),接收機(jī)所接收到的信號(hào)波形相比于發(fā)射信號(hào)表現(xiàn)出被拉伸或者壓縮的特性. 對(duì)于窄帶信號(hào),該影響可近似認(rèn)為接收信號(hào)相比發(fā)射信號(hào)的載頻產(chǎn)生了多普勒頻移. 因此,對(duì)于窄帶信號(hào),在相對(duì)運(yùn)動(dòng)存在時(shí),常通過估計(jì)兩路信號(hào)的時(shí)差和多普勒頻差來定位未知的輻射源[1-4]. 然而,對(duì)于寬帶信號(hào),繼續(xù)沿用多普勒頻移的設(shè)定會(huì)導(dǎo)致定位出現(xiàn)較大誤差. 為了對(duì)輻射源進(jìn)行精確定位,需要估計(jì)的參數(shù)為接收信號(hào)間的時(shí)差和尺度差.

寬帶互模糊函數(shù)(Wideband Cross Ambiguity Function, WBCAF)是估計(jì)兩路信號(hào)間時(shí)差和尺度差的常用工具之一[5-7]. WBCAF將兩路信號(hào)在時(shí)域和尺度域進(jìn)行二維相關(guān),其峰值的位置對(duì)應(yīng)真實(shí)的時(shí)差和尺度差. 由于WBCAF與連續(xù)小波變換具有相同的表達(dá)式,因此WBCAF的計(jì)算可通過小波變換來實(shí)現(xiàn). 然而,在無源定位中,接收信號(hào)為未知輻射源輻射的信號(hào),其解析式往往未知,直接應(yīng)用小波變換需要對(duì)接收信號(hào)進(jìn)行變采樣率處理,運(yùn)算量過大,尤其當(dāng)尺度差接近1時(shí),往往難以實(shí)現(xiàn)[8]. 文獻(xiàn)[9]利用互小波變換的性質(zhì),分別計(jì)算出兩路接收信號(hào)相對(duì)于某一給定的小波的小波變換,再由小波變換的結(jié)果獲得兩路接收信號(hào)的WBCAF. 在該方法中,伸縮變換的對(duì)象為解析式已知的小波,無需對(duì)接收信號(hào)進(jìn)行變采樣率處理,因而降低了部分運(yùn)算量. 然而,該方法需要聯(lián)合兩路接收信號(hào)小波變換的結(jié)果來計(jì)算WBCAF,運(yùn)算量仍然較大. 文獻(xiàn)[10]提出了一種對(duì)離散信號(hào)進(jìn)行尺度變換的快速計(jì)算方法,并將之用于WBCAF的計(jì)算,能夠有效降低運(yùn)算量. 不過,該方法仍需二維搜索WBCAF峰值,在尺度差需要搜索的范圍較大時(shí),運(yùn)算量偏大,不利于實(shí)時(shí)處理.

本文提出了一種基于尺度變換的時(shí)差/尺度差的快速估計(jì)方法. 由于線性調(diào)頻(Linear Frequency Modulation, LFM)信號(hào)在雷達(dá)、聲吶等領(lǐng)域的廣泛應(yīng)用[11-13],本文所針對(duì)的信號(hào)為寬帶LFM信號(hào). 對(duì)于寬帶LFM信號(hào),接收機(jī)所接收到的同一輻射源信號(hào)間的尺度差的平方等于兩路接收信號(hào)調(diào)頻率之比. 通過分?jǐn)?shù)階傅里葉變換(Fractional Fourier Transform, FrFT)分別估計(jì)出兩路接收信號(hào)的調(diào)頻率,即可獲得尺度差的估計(jì). 在搜索最優(yōu)角度的過程中,為了提高搜索效率,采用了分級(jí)搜索的結(jié)構(gòu)．利用估計(jì)到的尺度差對(duì)第一路接收信號(hào)進(jìn)行伸縮,并將伸縮后的信號(hào)與另一路接收信號(hào)做時(shí)域相關(guān),由時(shí)域相關(guān)峰的位置估計(jì)出時(shí)差. 與傳統(tǒng)基于WBCAF的方法相比,該方法只需一維搜索出接收的LFM信號(hào)所對(duì)應(yīng)的最優(yōu)角度即可估計(jì)得到尺度差,再一維搜索一次時(shí)域相關(guān)峰的峰值獲得時(shí)差的估計(jì). 由于避免了二維搜索模糊函數(shù)峰值,并且只需對(duì)接收信號(hào)做一次尺度變換,該方法能夠顯著降低估計(jì)所需的運(yùn)算量. 仿真實(shí)驗(yàn)表明隨著信噪比的提高,該方法所估計(jì)的時(shí)差和尺度差的均方誤差逐漸接近克拉美-羅下界.

1 信號(hào)模型

兩路接收信號(hào)可以表示為[14]

(1)

式中:

(2)

常用估計(jì)時(shí)差/尺度差的方法是二維搜索WBCAF的峰值. 兩路接收信號(hào)r1(t)和r2(t)的WBCAF可以表示為

(3)

式中,(*)表示復(fù)共軛. 不考慮噪聲時(shí),式(3)表示為

(4)

針對(duì)上述問題,文獻(xiàn)[10]給出了一種對(duì)離散信號(hào)進(jìn)行快速尺度變換的算法,并將該方法用于估計(jì)寬帶信號(hào)的時(shí)差和尺度差. 尺度變換的思路是對(duì)接收到的離散信號(hào)用Sinc函數(shù)重構(gòu)其對(duì)應(yīng)的連續(xù)信號(hào),并對(duì)重構(gòu)后的信號(hào)進(jìn)行伸縮,最后再對(duì)伸縮后的信號(hào)進(jìn)行采樣. 尺度變換的快速實(shí)現(xiàn)如下所示,假定接收到的離散信號(hào)為

x(k)=s(k)+n(k),-N′≤k≤N′,

(5)

則希望得到的尺度變換后的信號(hào)為[10]

(6)

式(4)中的WBCAF可寫為

(7)

雖然文獻(xiàn)[10]提出的尺度變換方法能夠用于估計(jì)時(shí)差和尺度差,然而該方法需要對(duì)接收信號(hào)r1(t)在不同尺度下進(jìn)行變換,當(dāng)需要變換的尺度數(shù)量變大時(shí),其運(yùn)算量急劇增大,不利于實(shí)時(shí)處理. 針對(duì)該問題,本文提出了基于尺度變換的時(shí)差和尺度差估計(jì)方法. 該方法首先利用FrFT估計(jì)出兩路接收信號(hào)的尺度差,然后用該尺度差對(duì)第一路接收信號(hào)做伸縮,并將伸縮后的信號(hào)與第二路接收信號(hào)做時(shí)域相關(guān),由相關(guān)峰的位置估計(jì)出兩路接收信號(hào)的時(shí)差. 該方法避免了二維搜索模糊函數(shù)的峰值,只需對(duì)信號(hào)進(jìn)行一次尺度變換,并且可利用快速傅里葉變換(Fast fourier Transform,FFT)實(shí)現(xiàn),能夠顯著降低估計(jì)所需的運(yùn)算量.

2 基于尺度變換的時(shí)差/尺度差估計(jì)

2.1尺度差的估計(jì)

根據(jù)式(2),對(duì)s1(t)做尺度為1/σ0的伸縮變換并將其延時(shí)τ0即可得到

(8)

LFM信號(hào)調(diào)頻率的估計(jì)常用FrFT來實(shí)現(xiàn). 某個(gè)信號(hào)x(t)的FrFT為

(9)

式中,

(10)

cotαopt=-m．

(11)

因此,估計(jì)某個(gè)LFM信號(hào)x(t)的最優(yōu)角度可通過搜索x(t)在不同角度上的FrFT的峰值來實(shí)現(xiàn),FrFT取到最大值所對(duì)應(yīng)的角度即為最優(yōu)的角度,有

(12)

當(dāng)LFM信號(hào)x(t)湮沒在白噪聲中時(shí),式(12)仍然可以直接用于估計(jì)最優(yōu)角度,這是因?yàn)榘自肼曉谌我饨嵌鹊腇rFT不會(huì)形成峰值. 假設(shè)s1(t)所對(duì)應(yīng)的最優(yōu)角度為α1,s2(t)所對(duì)應(yīng)的最優(yōu)角度為α2,則兩路信號(hào)最優(yōu)角度的估計(jì)為

(13)

式中,R1,α(u)和R2,α(u)分別為接收信號(hào)r1(t)和r2(t)在角度α的FrFT. 結(jié)合式(11)和式(13),獲得尺度差的估計(jì)為

(14)

在搜索LFM信號(hào)對(duì)應(yīng)的最優(yōu)角度時(shí),為了降低運(yùn)算量,可以采用分級(jí)搜索來減少搜索角度的個(gè)數(shù),一般只需三級(jí)分級(jí)搜索即可搜索到最優(yōu)角度．FrFT需要搜索的初始角度范圍為[-π/2,π/2],假設(shè)初始搜索步長為Δα,且有0<Δα<π,采用三級(jí)分級(jí)搜索的步驟為:

分級(jí)搜索的優(yōu)勢(shì)在于可以用較低的運(yùn)算量獲得同樣的角度精度．在采用三級(jí)分級(jí)搜索時(shí),可以看到,當(dāng)搜索角度的精度為Δα/100時(shí),步驟1需要搜索的角度數(shù)為π/Δα,步驟2和步驟3分別需要搜索20個(gè)角度,因此采用三級(jí)搜索只需搜索π/Δα+40個(gè)角度,而如果直接以Δα/100搜索整個(gè)角度范圍需要搜索100π/Δα個(gè)角度.由于0<Δα<π,顯然在搜索最優(yōu)角度時(shí)采用分級(jí)搜索能夠顯著減少運(yùn)算量.

2.2時(shí)差的估計(jì)

(15)

(16)

(17)

3 運(yùn)算量分析

本節(jié)首先分析了本文所提方法的運(yùn)算量,接著與文獻(xiàn)[10]中尺度變換方法的運(yùn)算量進(jìn)行了對(duì)比. 假設(shè)處理的信號(hào)長度為N,需要搜索的尺度個(gè)數(shù)為L.

估計(jì)尺度差主要的運(yùn)算量為計(jì)算FrFT,根據(jù)文獻(xiàn)[18]中FrFT的快速算法,一次FrFT所需的復(fù)乘次數(shù)為3Nlog2N+3N, 由于需要對(duì)兩路信號(hào)都做FrFT,因此一個(gè)角度下的FrFT需要進(jìn)行6Nlog2N+6N次復(fù)乘. 雖然需要搜索的角度個(gè)數(shù)為L,然而通過分級(jí)搜索可將搜索的角度減少為M(M一般為幾十),因此估計(jì)尺度差需要的復(fù)乘次數(shù)約為6MNlog2N+6MN. 時(shí)差估計(jì)的運(yùn)算量分為兩部分:文獻(xiàn)[10]中給出了信號(hào)進(jìn)行一次尺度變換所需的復(fù)乘次數(shù)為4Nlog2N+4N;時(shí)域相關(guān)利用FFT來實(shí)現(xiàn),包括三次FFT和一次信號(hào)點(diǎn)乘,一次時(shí)域相關(guān)需要計(jì)算3Nlog2N+N次復(fù)乘. 因此,運(yùn)用本文方法所需的復(fù)乘次數(shù)約為6MNlog2N+6MN+7Nlog2N+5N.

文獻(xiàn)[10]中的尺度變換方法在一個(gè)尺度下所需的運(yùn)算量與本文方法估計(jì)時(shí)差時(shí)候的運(yùn)算量一致,需要的復(fù)乘次數(shù)為7Nlog2N+5N,因此L個(gè)尺度下需要計(jì)算7LNlog2N+5LN次復(fù)乘.

表1為不同信號(hào)點(diǎn)數(shù)N下的本文方法與尺度變換方法的運(yùn)算量對(duì)比,假設(shè)需要搜索的尺度個(gè)數(shù)L=1 000,則采用三級(jí)分級(jí)搜索后需要搜索的尺度個(gè)數(shù)M=50. 從表1可以看到,運(yùn)用本文方法估計(jì)時(shí)差和尺度差所需的運(yùn)算量要遠(yuǎn)少于尺度變換方法.

表1 本文方法與尺度變換方法運(yùn)算量對(duì)比

4 仿真實(shí)驗(yàn)

本節(jié)首先給出了本文方法在低信噪比下所估計(jì)到的時(shí)差和尺度差的仿真;接著,在不同信噪比下,對(duì)比了本文方法和尺度變換方法的估計(jì)性能,并與克拉美-羅下界進(jìn)行了比較.

4.1尺度差和時(shí)差的估計(jì)

圖1 第一路接收信號(hào)的分?jǐn)?shù)階傅里葉變換

圖2 第二路接收信號(hào)的分?jǐn)?shù)階傅里葉變換

圖3 時(shí)域相關(guān)結(jié)果

圖4 第一路接收信號(hào)的分?jǐn)?shù)階傅里葉變換

圖5 第二路接收信號(hào)的分?jǐn)?shù)階傅里葉變換

圖6 時(shí)域相關(guān)結(jié)果

4.2均方根誤差以及估計(jì)時(shí)間的對(duì)比

在不同信噪比下,仿真了利用本文方法和尺度變換方法所估計(jì)到的時(shí)差和尺度差的均方根誤差,并與克拉美-羅下界進(jìn)行了對(duì)比. 克拉美-羅下界的仿真曲線根據(jù)文獻(xiàn)[9]中給出的寬帶信號(hào)時(shí)差/尺度差估計(jì)的克拉美-羅下界公式所獲得. 圖7為時(shí)差估計(jì)的均方根誤差仿真圖,圖8為尺度差估計(jì)的均方根誤差仿真圖. 從圖7和圖8可以看到:本文方法與尺度變換方法所估計(jì)的均方根誤差曲線基本一致;隨著信噪比的提高,兩種方法所估計(jì)到的時(shí)差和尺度差的均方根誤差均明顯降低,并逐漸接近克拉美-羅下界.

圖7 時(shí)差估計(jì)的均方誤差曲線

圖8 尺度差估計(jì)的均方誤差曲線

圖9為本文方法與尺度變換方法估計(jì)時(shí)差和尺度差所需CPU時(shí)間的對(duì)比圖.仿真中采樣率fs=500 MHz,所需搜索的尺度個(gè)數(shù)L=1 000,運(yùn)用本文方法估計(jì)時(shí),通過采用三級(jí)分級(jí)搜索處理結(jié)構(gòu),所需搜索的角度個(gè)數(shù)減少為M=50. 由圖9可看出:本文方法所需的估計(jì)時(shí)間要明顯少于尺度變換的方法;尤其在處理長信號(hào)時(shí),運(yùn)用本文方法能夠更好滿足實(shí)時(shí)處理的要求.

圖9 本文方法與尺度變換方法估計(jì)時(shí)間對(duì)比

5 結(jié) 論

針對(duì)傳統(tǒng)的時(shí)差/尺度差估計(jì)方法存在運(yùn)算量大導(dǎo)致無法實(shí)時(shí)估計(jì)的問題,本文提出了一種寬帶線性調(diào)頻信號(hào)的時(shí)差和尺度差估計(jì)的快速算法．該方法通過估計(jì)信號(hào)的調(diào)頻率來獲得尺度差的估計(jì),進(jìn)而計(jì)算WBCAF在尺度差處的切片來獲得時(shí)差的估計(jì)．與傳統(tǒng)的尺度變換方法相比,該方法僅需對(duì)信號(hào)進(jìn)行一次伸縮,能夠有效提高估計(jì)效率．仿真實(shí)驗(yàn)證明,本文提出的方法能夠在不損失估計(jì)性能的同時(shí)顯著提高估計(jì)的效率,能夠?qū)崟r(shí)估計(jì)出時(shí)差和尺度差．

[1] 孫正波, 葉尚福. 一種時(shí)差/頻率差快速聯(lián)合估計(jì)方法[J]. 電波科學(xué)學(xué)報(bào), 2006, 21(5): 641-646.

SUN Z B, YE S F. Fast algorithm for joint estimation of DTO and DFO[J]. Chinese journal of radio science, 2006, 21(5): 641-646. (in Chinese)

[2] KIM Y H, KIM D G, KIM H N. Two-step estimator for moving-emitter geolocation using time difference of arrival/frequency-difference of arrival measurements[J]. IET radar, sonar & navigation, 2015, 9(7): 881-887.

[3] HO K C, LU X, KOVAVISARUCH L. Source localization using TDOA and FDOA measurements in the presence of receiver location errors: Analysis and solution[J]. IEEE transactions on signal processing, 2007, 55(2): 684-696.

[4] AMAR A, LEUS G, FRIEDLANDER B. Emitter localization given time delay and frequency shift measurements[J]. IEEE transactions on aerospace and electronic systems, 2012, 48(2): 1826-1837.

[5] STEIN S. Algorithms for ambiguity function processing[J]. IEEE transactions on aerospace and signal processing, 1981, 29(3): 588-599.

[6] SHARIF M R, ABEYSEKERA S S. Efficient wideband signal parameter estimation using a radon-ambiguity transform slice[J]. IEEE transactions on AES, 2007, 43(2): 673-688.

[7] 孫正波, 葉尚福. 利用互模糊函數(shù)實(shí)現(xiàn)衛(wèi)星干擾源定位[J]. 電波科學(xué)學(xué)報(bào), 2004, 19(5): 525-529.

SUN Z B, YE S F. Satellite interference location using cross ambiguity function[J]. Chinese journal of radio science, 2004, 19(5): 525-529. (in Chinese)

[8] WEISS L G. Wavelets and wideband correlation processing[J]. IEEE signal processing magazine, 1994, 11(1): 13-32.

[9] NIU X X, CHING P C. Wavelet based approach for joint time delay and Doppler stretch measurements[J]. IEEE transactions on aerospace and electronic systems, 1999, 35(3): 1111-1119.

[10] HO K C, CHAN Y T. Optimum discrete wavelet scaling and its application to delay and Doppler estimation[J]. IEEE transactions on signal processing, 1998, 46(9): 2285-2290.

[11] CHEN X L, GUAN J, LIU N B, et al. Maneuvering target detection via radon-fractional Fourier transform-based long-time coherent integration[J]. IEEE transactions on signal processing, 2014, 62(4): 939-953.

[12] 朱文濤, 蘇濤, 楊濤, 等. 線性調(diào)頻連續(xù)波信號(hào)檢測(cè)與參數(shù)估計(jì)算法[J]. 電子與信息學(xué)報(bào), 2014, 36(3): 552-558.

ZHU W T, SU T, YANG T, et al. Detection and parameter estimation of linear frequency modulation continuous wave signal[J]. Journal of electronics & information technology, 2014, 36(3): 552-558. (in Chinese)

[13] 朱文濤, 蘇濤, 楊濤, 等. 低信噪比下線性調(diào)頻連續(xù)波信號(hào)的參數(shù)估計(jì)[J]. 電波科學(xué)學(xué)報(bào), 2013, 28(6): 1158-1164.

ZHU W T, SU T, YANG T, et al. Parameter estimation of linear frequency modulated continuous wave signal in low SNR[J]. Chinese journal of radio science, 2013, 28(6): 1158-1164. (in Chinese)

[14] FOWLER M L, HU X. Signal models for TDOA/FDOA estimation[J]. IEEE transactions on aerospace and electronic systems, 2008, 44 (4): 1543-1549.

[15] 杜雨洺, 楊建宇. 基于FRFT的LFMCW雷達(dá)加速動(dòng)目標(biāo)檢測(cè)與參數(shù)估計(jì)[J]. 電波科學(xué)學(xué)報(bào), 2005, 20(6): 815-819.

DU Y M, YANG J Y. Novel method of moving target detection and parameter estimation for LFMCW radar[J]. Chinese journal of radio science, 2005, 20(6): 815-819. (in Chinese)

[16] 李會(huì)勇, 許丁文, 胡進(jìn)峰, 等. 基于FRFT的天波雷達(dá)多機(jī)動(dòng)目標(biāo)檢測(cè)[J]. 系統(tǒng)工程與電子技術(shù), 2014, 36(9): 1725-1730.

LI H Y, XU D W, HU J F, et al. FRFT based algorithm for maneuvering target detection with OTH radar[J]. Systems engineering and electronics, 2014, 36(9): 1725-1730. (in Chinese)

[17] OONINCX P J. Joint time-frequency offset detection using the fractional Fourier transform[J]. Signal processing, 2008, 88: 2936-2942.

[18] OZAKTAS H M, ARIKAN O, KUTAY M A, et al. Digtal computation of the fractional Fourier transform[J]. IEEE transactions on signal processing, 1996, 44(9): 2141-2150.

郭付陽(1991—),男,江西人,博士研究生,研究方向?yàn)闊o源定位、雷達(dá)信號(hào)處理.

張子敬(1967—),男,陜西人,教授,博士生導(dǎo)師,研究方向?yàn)闊o源定位、雷達(dá)信號(hào)處理.

楊林森(1988—),男,陜西人,博士研究生,研究方向?yàn)闊o源定位、雷達(dá)信號(hào)處理.

Scaling-basedTDOA/SDOAestimationalgorithmforwidebandchirpsignals

GUOFuyangZHANGZijingYANGLinsen

(NationalLabofRadarSignalProcessing,XidianUniversity,Xi’anShaanxi710071,China)

A method based on scaling is proposed to estimate the time difference of arrival (TDOA) and scale difference of arrival (SDOA) between two

wideband chirp signals. Using the relation that the square of SDOA equals to the ratio of chirp-rates of two received chirp signals, each chirp-rate of two chirp signals is first estimated using the fractional Fourier transform, and then the SDOA can be evaluated. By scaling one received chirp signal with the estimated SDOA and evaluating the correlation of the scaled chirp signal and the other received chirp signal, the TDOA is finally estimated. Since the 2-D searching of the peak position of the wideband cross ambiguity function is avoided and it can be evaluated using only few fast Fourier transforms, the computational cost is significantly reduced. Simulation results show that the root mean square errors of the estimated TDOA and SDOA using proposed method closely meet to the Cramer-Rao lower bound under high signal-to-noise ratios.

fractional Fourier transform(FrFT); wideband cross ambiguity ambiguty(WBCAF); TDOA/SDOA estimation; wideband chirp signal;Cramer-Rao lower bound

郭付陽, 張子敬, 楊林森. 基于尺度變換的寬帶線性調(diào)頻信號(hào)時(shí)差/尺度差估計(jì)算法[J]. 電波科學(xué)學(xué)報(bào),2017,32(4):441-448.

10.13443/j.cjors.2017060201

GUO F Y, ZHANG Z J, YANG L S. Scaling-based TDOA/SDOA estimation algorithm for wideband chirp signals[J]. Chinese journal of radio science,2017,32(4):441-448. (in Chinese). DOI: 10.13443/j.cjors.2017060201

TN911.7

A

1005-0388(2017)04-0441-08

DOI10.13443/j.cjors.2017060201

2017-06-02

國家自然科學(xué)基金(No.61571349)

聯(lián)系人： 張子敬 E-mail: zjzhang@xidian.edu.cn