陳一暢 張 群
?
一種基于SAR稀疏采樣數(shù)據(jù)的動目標(biāo)運(yùn)動參數(shù)估計方法
陳一暢*①②張 群①③
①(空軍工程大學(xué)信息與導(dǎo)航學(xué)院 西安 710077)②(清華大學(xué)電子工程系 北京 100084)③(復(fù)旦大學(xué)電磁波信息科學(xué)教育部重點(diǎn)實(shí)驗室 上海 200433)
該文針對地面動目標(biāo)運(yùn)動參數(shù)估計問題進(jìn)行研究,提出一種利用單天線合成孔徑雷達(dá)(SAR)稀疏采樣數(shù)據(jù)的動目標(biāo)2維速度估計方法。首先以目標(biāo)2維速度為參數(shù)構(gòu)建一個等效參數(shù)化模型將動目標(biāo)回波數(shù)據(jù)轉(zhuǎn)化為小斜視回波數(shù)據(jù),然后利用改進(jìn)的迭代閾值算法實(shí)現(xiàn)不同參數(shù)條件下的動目標(biāo)2維成像,最后以成像結(jié)果的圖像熵值為優(yōu)化準(zhǔn)則對初始模型參數(shù)進(jìn)行搜索,從而獲得準(zhǔn)確的動目標(biāo)運(yùn)動參數(shù)。該方法以稀疏采樣數(shù)據(jù)為輸入,可以減少所需數(shù)據(jù)量,并且能夠有效避免多普勒模糊問題,在較低信雜比條件下仍然能夠準(zhǔn)確估計出目標(biāo)運(yùn)動參數(shù)。仿真實(shí)驗結(jié)果驗證了所提動目標(biāo)參數(shù)估計方法的有效性。
合成孔徑雷達(dá);參數(shù)化模型;運(yùn)動目標(biāo);參數(shù)估計
合成孔徑雷達(dá)(SAR)是一種裝載在航空航天平臺,具備對地全天候、全天時成像能力的微波遙感設(shè)備[1]。其中對地面動目標(biāo)運(yùn)動參數(shù)估計是SAR的一項重要應(yīng)用,準(zhǔn)確的運(yùn)動參數(shù)估計對于動目標(biāo)檢測、成像、識別都有著重要作用。傳統(tǒng)的動目標(biāo)運(yùn)動參數(shù)估計方法大多是基于Doppler相位分析,通過估計出回波信號的Doppler頻偏和Doppler調(diào)頻率,從而計算出動目標(biāo)在距離向和方位向的2維速度參數(shù)[2,3]。但是這一類方法存在固有缺陷,當(dāng)SAR以一定的脈沖重復(fù)頻率(PRF)在方位向采樣時,只可對區(qū)間內(nèi)的Doppler頻譜進(jìn)行分析,而當(dāng)目標(biāo)運(yùn)動速度較快時,可能使得Doppler頻偏超出此區(qū)間,即Doppler模糊,最終使得目標(biāo)速度估計存在較大誤差。解決Doppler模糊問題最簡單直接的辦法是增大PRF,但是同時又會引入縮短模糊距離、增加數(shù)據(jù)量等新問題。許多文獻(xiàn)基于多通道SAR系統(tǒng)提出了有效的動目標(biāo)速度估計方法,以及動目標(biāo)成像算法[4,5],然而,需要指出的是多通道解決Doppler模糊問題的同時也增加了硬件復(fù)雜度。也有文獻(xiàn)針對單天線SAR數(shù)據(jù),提出了動目標(biāo)參數(shù)估計方法,如文獻(xiàn)[6],將動目標(biāo)信號做距離徙動矯正(Range Cell Migration Correction, RCMC)和距離壓縮操作,然后以距離多普勒域的圖像對比度最大化為準(zhǔn)則,對目標(biāo)的2維速度進(jìn)行搜索。因為RCMC和距離壓縮操作與目標(biāo)速度密切相關(guān),所以只有根據(jù)準(zhǔn)確的目標(biāo)運(yùn)動參數(shù)對信號進(jìn)行RCMC和距離壓縮才能使距離多普勒圖像的對比度最大化。文獻(xiàn)[6]的方法可以有效克服Doppler模糊,但對信雜比有較高的要求。文獻(xiàn)[7]利用Radon變換估計目標(biāo)徑向速度,避免了Doppler模糊問題,然后采用一種“雙向分段成像”方法估計目標(biāo)方位向速度。所謂的“雙向分段成像”方法是利用前后兩段部分孔徑數(shù)據(jù)分別對動目標(biāo)進(jìn)行粗成像,然后根據(jù)兩幅粗像的位置差異計算方位向速度,獲得的估計精度有限。
需要指出的是上述動目標(biāo)運(yùn)動參數(shù)估計方法大都是基于SAR全采樣數(shù)據(jù)的估計方法。近年來,壓縮感知理論被廣泛地引入到SAR應(yīng)用中,一系列基于稀疏降采樣數(shù)據(jù)的方法被提出,主要用于SAR數(shù)據(jù)壓縮[8],提高SAR圖像分辨率[9,10],SAR平臺運(yùn)動補(bǔ)償[11]等方面。稀疏采樣方式是未來SAR實(shí)際應(yīng)用中一個重要發(fā)展方向,因此有必要研究基于稀疏采樣數(shù)據(jù)的動目標(biāo)運(yùn)動參數(shù)估計方法。
參數(shù)化稀疏表征是字典學(xué)習(xí)的一個特殊分支,能夠?qū)崿F(xiàn)雷達(dá)探測過程中未知參數(shù)的動態(tài)學(xué)習(xí)和雷達(dá)信號的最優(yōu)稀疏表征[12],其概念最早在文獻(xiàn)[13]中提到,被用于逆合成孔徑雷達(dá)(Inverse Synthetic Aperture Radar, ISAR)成像。文獻(xiàn)[12]綜述了參數(shù)化稀疏表征在雷達(dá)探測中的應(yīng)用,重點(diǎn)介紹了參數(shù)化稀疏表征在ISAR成像、SAR自聚焦和目標(biāo)識別中的應(yīng)用。本文結(jié)合壓縮感知理論,將目標(biāo)運(yùn)動參數(shù)作為影響觀測模型的唯一變量,通過優(yōu)化動目標(biāo)成像結(jié)果的稀疏性來估計目標(biāo)的運(yùn)動參數(shù)。文章首先分析了動目標(biāo)回波信號稀疏采樣數(shù)據(jù)的特點(diǎn),然后構(gòu)建一個轉(zhuǎn)換模型,將動目標(biāo)回波信號轉(zhuǎn)換為等效的斜視SAR回波信號??梢宰C明,等效后的平臺速度、斜視角與原動目標(biāo)模型中目標(biāo)2維速度成一一對應(yīng)關(guān)系。利用作者前期提出的斜視SAR壓縮感知成像算法[14],實(shí)現(xiàn)動目標(biāo)的粗成像。在目標(biāo)2維速度速度域進(jìn)行搜索,所選取的2維速度越接近真值,所獲得的成像結(jié)果聚焦性能將越好。最終本文以圖像熵最小化為準(zhǔn)則,在2維速度域估計出動目標(biāo)徑向和方位向2維速度。本文方法可以克服Doppler模糊問題,相比于傳統(tǒng)方法,可以采用稀疏采樣數(shù)據(jù)完成對動目標(biāo)的運(yùn)動參數(shù)精確估計,減少了所需數(shù)據(jù)量的同時獲得較高的估計精度。文章采用仿真數(shù)據(jù)進(jìn)行了實(shí)驗,說明了本文算法的有效性,并且討論了信雜比對本文估計方法的影響,分析了本文估計方法對目標(biāo)徑向速度和方位向速度的估計精度。
本文考慮傳統(tǒng)單天線SAR系統(tǒng),工作模式為正側(cè)視條帶式,SAR平臺與觀測場景2維示意圖如圖1(a)所示。假定SAR裝載在飛機(jī)平臺上,并沿著軸以速度做勻速直線飛行。假設(shè)觀測場景中有一動目標(biāo)P,徑向速度分量記為(以靠近SAR平臺飛行軌跡方向為正方向),方位向速度分量記為(以SAR平臺速度方向為正方向),即我們需要估計的參數(shù)。為了便于說明方法推導(dǎo),我們首先采用全采樣數(shù)據(jù)進(jìn)行模型分析。SAR平臺以“走停”模式發(fā)射并接受信號,脈沖重復(fù)頻率記為PRF。圖1中標(biāo)出了各虛擬合成孔徑陣元位置。不失一般性,我們假設(shè)目標(biāo)處于SAR平臺正側(cè)方時坐標(biāo)為,其中0表示方位向坐標(biāo),表示距離向坐標(biāo),并以當(dāng)前時刻為慢時間零點(diǎn)。則動目標(biāo)與SAR瞬時距離為
本文采用線性調(diào)頻(Linear Frequency Modulation, LFM)信號作為SAR發(fā)射信號,運(yùn)動參數(shù)為的運(yùn)動目標(biāo)回波可以表示為
(2)
(4)
(6)
在全采樣數(shù)據(jù)條件下,可以采用多種經(jīng)典算法對斜視模型信號進(jìn)行成像處理,如非線性頻調(diào)變標(biāo)算法(Non-linear Chirp Scaling, NCS),距離徙動算法(Range Migration Algorithm, RMA)等。本文所建立的模型在距離向上對回波信號采用隨機(jī)稀疏采樣方式,其采樣結(jié)果可以看作是全采樣數(shù)據(jù)在距離向上的隨機(jī)降維觀測。對于稀疏采樣數(shù)據(jù),傳統(tǒng)算法將不再適用,因此本文利用前期工作中提出的基于NCS算子的斜視SAR壓縮感知成像方法對動目標(biāo)回波稀疏采樣數(shù)據(jù)進(jìn)行成像[14],然后根據(jù)成像結(jié)果估計動目標(biāo)參數(shù),具體操作步驟在第3節(jié)論述。
(8)
(10)
上述方法是在合速度和斜視角確定情況下的稀疏采樣數(shù)據(jù)目標(biāo)成像方法,針對運(yùn)動參數(shù)未知的運(yùn)動目標(biāo),其等效得出的合速度和斜視角也是未知的。采用不同的運(yùn)動參數(shù)進(jìn)行上述過程求解,可以獲得不同的成像結(jié)果,因此目標(biāo)像可以看作是目標(biāo)速度參數(shù)的隱函數(shù)。當(dāng)與實(shí)際參數(shù)匹配時,才能得到清晰的目標(biāo)2維像。圖像熵是常用的衡量復(fù)圖像清晰度的指標(biāo)[19],因此本文采用熵值最小化為優(yōu)化目標(biāo),在2維速度空間內(nèi)估計目標(biāo)運(yùn)動參數(shù)。復(fù)圖像的計算公式為
(12)
步驟1 在整幅SAR圖像中提取感興趣的動目標(biāo)區(qū)域,并將其轉(zhuǎn)換到原始數(shù)據(jù)域;
步驟2 對2維速度空間離散化,并根據(jù)式(3)將離散點(diǎn)從2維速度空間映射到合速度-斜視角空間,利用ITA算法求解式(10)所述稀疏重構(gòu)問題;
步驟3 根據(jù)式(11)計算重構(gòu)結(jié)果的圖像熵;
關(guān)于上述速度估計方法,需要說明兩點(diǎn):第一,動目標(biāo)在整個SAR觀測場景中往往是稀疏分布的,具有空間稀疏性,且步驟1可以去除部分背景雜波,提高動目標(biāo)局部信雜比;第二,在搜索圖像最小熵時,可以采取變步長搜索方式,首先采用大步長搜索減少運(yùn)算時間,然后在極值附近采用小步長搜索提高估計精度。
本節(jié)采用仿真數(shù)據(jù)驗證本文方法的有效性,并根據(jù)實(shí)驗結(jié)果分析本文方法性能。實(shí)驗中用到的SAR主要參數(shù)如下:SAR平臺速度,載波波長,發(fā)射信號帶寬,脈寬,天線實(shí)際孔徑長度,脈沖重復(fù)頻率,場景中心與SAR平臺飛行軌跡之間的距離為,距離向全采樣點(diǎn)數(shù)為。觀測場景如圖2所示,尺寸為,其中包含若干靜止散射點(diǎn),A為勻速運(yùn)動目標(biāo),其方位向速度分量為,距離向速度分量為。根據(jù)動目標(biāo)Doppler中心計算公式,可以得到,已經(jīng)超出,因此會出現(xiàn)Doppler模糊。
首先,我們通過實(shí)驗驗證本文所構(gòu)建的轉(zhuǎn)換模型的正確性。對距離向數(shù)據(jù)進(jìn)行隨機(jī)稀疏采樣,僅錄取個采樣點(diǎn),定義降采樣率為。圖3(a)是采用模型參數(shù)獲得的成像結(jié)果,其中靜止目標(biāo)可以獲得準(zhǔn)確聚焦結(jié)果,但是運(yùn)動目標(biāo)無法正常聚焦,并且其散焦效應(yīng)使得靜止散射點(diǎn)成像受到干擾。圖3(b)采用本文模型對動目標(biāo)進(jìn)行成像,最終獲得整個場景的準(zhǔn)確重構(gòu)結(jié)果,實(shí)驗結(jié)果說明了采用本文所提出的等效斜視模型,在運(yùn)動參數(shù)確定的情況下可以實(shí)現(xiàn)對運(yùn)動目標(biāo)的準(zhǔn)確成像。
在運(yùn)動目標(biāo)運(yùn)動參數(shù)未知情況下,利用本文方法對運(yùn)動參數(shù)進(jìn)行估計,仿真實(shí)驗中設(shè)定信雜比為。以SAR平臺速度的20%為運(yùn)動參數(shù)估計區(qū)間,即在仿真實(shí)驗中,與搜索范圍均為。計算各搜索點(diǎn)的圖像熵值,結(jié)果如圖4。因為圖中縱軸采用負(fù)熵值坐標(biāo),因此峰值位置對應(yīng)了最小熵位置,從圖4中可以看出峰值位置對應(yīng)的運(yùn)動參數(shù)估計結(jié)果為,這與目標(biāo)運(yùn)動參數(shù)真值完全一致。本文方法提供了目標(biāo)在方位向和距離向2個維度上的速度分量,圖5給出了熵值圖在兩個維度上剖面圖,通過對比方位向和距離向速度-熵值剖面圖,可以看出圖5(a)中方位向熵值圖的峰值寬度明顯小于圖5(b)中距離向熵值圖的峰值寬度。這與實(shí)際情況相符,方位向的速度誤差對圖像散焦影響更大。為了進(jìn)一步分析本文所提估計方法的性能,在不同信雜比條件下,分別考察方位向和距離向速度參數(shù)的估計精度,同時分析降采樣率對估計精度的影響。圖6給出了不同條件下,估計結(jié)果的相對誤差結(jié)果,其中最小搜索步長為。從圖中可以看出隨著信雜比的提高,相對誤差逐漸減小。降采樣率為條件下,當(dāng)信雜比達(dá)到6 dB以上時,方位向速度分量和距離向速度分量估計結(jié)果的相對誤差均可控制在2.5%以內(nèi)。此外,本文算法處理過程中用到了稀疏重構(gòu)理論,降采樣率隊估計精度有一定影響,隨著估計降采樣率增大,相對誤差也有所增大。
圖2 觀測場景???????????圖3 稀疏采樣回波重構(gòu)結(jié)果
圖4 2維速度空間動目標(biāo)粗像熵值分布圖
圖5 目標(biāo)粗像熵值1維剖面圖
圖6 動目標(biāo)速度參數(shù)相對誤差與信雜比(SCR)關(guān)系曲線
當(dāng)場景中有多個不用運(yùn)動參數(shù)的運(yùn)動目標(biāo)時,可以通過圖像分割的方法分別估計不同目標(biāo)的運(yùn)動參數(shù)。若不同目標(biāo)散焦圖像之間有重疊,無法分離,也可采用本文方法同時估計2個目標(biāo)的運(yùn)動參數(shù)。以前一節(jié)仿真實(shí)驗為例,將圖2中B目標(biāo)更改為運(yùn)動目標(biāo),其2維運(yùn)動參數(shù)為。圖7給出了包含估計結(jié)果信息的熵值圖,圖中出現(xiàn)兩個峰值,分別對應(yīng)兩個運(yùn)動目標(biāo)的運(yùn)動參數(shù)。
圖7 雙目標(biāo)2維速度空間動目標(biāo)粗像熵值分布圖
本文根據(jù)速度合成原理將SAR動目標(biāo)觀測模型轉(zhuǎn)換為斜視SAR觀測模型,然后利用改進(jìn)的ITA算法對動目標(biāo)回波稀疏采用數(shù)據(jù)進(jìn)行處理,以動目標(biāo)粗像熵值最小化為準(zhǔn)則,在2維速度空間估計動目標(biāo)的運(yùn)動參數(shù)。利用本文提出的動目標(biāo)運(yùn)動參數(shù)估計方法,可以準(zhǔn)確估計出稀疏采樣數(shù)據(jù)條件下的動目標(biāo)運(yùn)動參數(shù),在降采樣率條件下,當(dāng)信雜比達(dá)到6 dB時,估計值的相對誤差可以控制在2.5%以內(nèi)。此外,本文估計方法還可以克服Doppler模糊,并且能夠同時實(shí)現(xiàn)多個目標(biāo)運(yùn)動參數(shù)的估計,本文采用仿真數(shù)據(jù)對模型正確性和算法性能進(jìn)行了驗證分析。下一步工作重點(diǎn)是設(shè)計相應(yīng)的迭代算法,實(shí)現(xiàn)熵值最小化。
[1] SAEEDI J and FAEZ K. Synthetic aperture radar imaging using nonlinear frequency modulation signal[J]., 2016, 52(1): 99-110. doi: 10.1109/TAES.2015.140310.
[2] WERNESS S A S, CARRARA W G, JOYCE L S,. Moving target imaging algorithm for SAR data[J]., 1990, 26(1): 57-67. doi: 10.1109/7.53413.
[3] PERRY R P, DIPIETRO R C, and FANTE R L. SAR imaging of moving targets[J]., 1999, 35(1): 188-200. doi: 10.1109/7.745691.
[4] ZHANG Shuangxi, XING Mengdao, XIA Xianggen,. Robust clutter suppression and moving target imaging approach for multichannel in azimuth high-resolution and wide-swath synthetic aperture radar[J]., 2015, 53(2): 687-709. doi: 10.1109/TGRS.2014.2327031.
[5] GUO Bin, VU Duc, XU Luzhou,. Ground moving target indication via multichannel airborne SAR[J]., 2011, 49(10): 3753-3764. doi: 10.1109/TGRS.2011.2143420.
[6] LI Gang, XIA Xianggen, XU Jia,. A velocity estimation algorithm of moving targets using single antenna SAR[J]., 2009, 45(3): 1052-1062. doi: 10.1109/TAES.2009.5259182.
[7] HOU Lili, SONG Hongjun, ZHENG Mingjie,. Fast moving target imaging and motion parameters estimation based on Radon transform and bi-directional approach[J].,&, 2016, 10(6): 1013-1023. doi: 10.1049/iet-rsn.2014.0567.
[8] 陳一暢, 張群, 朱麗莉, 等. 基于壓縮感知和矢量量化的 SAR 數(shù)據(jù)級聯(lián)壓縮方法[J]. 現(xiàn)代雷達(dá), 2013, 35(10): 36-40. doi: 10.16592/J.CNKI.1004-7859.2013.10.013.
CHEN Yichang, ZHANG Qun, ZHU Lili,. A cascade compress method of SAR data based on compressed sensing and vector quantization[J]., 2013, 35(10): 36-40. doi: 10.16592/J.CNKI.1004-7859.2013.10.013.
[9] BU Hongxia, TAO Ran, BAI Xia,. A novel SAR imaging algorithm based on compressed sensing[J]., 2015, 12(5): 1003-1007. doi: 10.1109/LGRS.2014.2372319.
[10] 陳一暢, 張群, 陳校平, 等. 多重測量矢量模型下的稀疏步進(jìn)頻率SAR成像算法[J]. 電子與信息學(xué)報, 2014, 36(12): 2986-2993. doi: 10.3724/SP.J.1146.2013.01831.
CHEN Yichang, ZHANG Qun, CHEN Xiaoping,. An imaging algorithm of sparse stepped frequency SAR based on multiple measurement vectors model[J].&, 2014, 36(12): 2986-2993. doi: 10.3724/SP.J.1146.2013.01831.
[11] GU Fufei, ZHANG Qun, CHI Long,. A novel motion compensating method for MIMO-SAR imaging based on compressed sensing[J]., 2015, 15(4): 2157-2165. doi: 10.1109/JSEN.2014.2371451.
[12] 李剛, 夏香根. 參數(shù)化稀疏表征在雷達(dá)探測中的應(yīng)用[J]. 雷達(dá)學(xué)報, 2016, 5(1): 1-7. doi: 10.12000/JR15126.
LIGang and XIA Xianggen. Parametric sparse representation and its applications in radar sensing[J]., 2016, 5(1): 1-7. doi: 10.12000/JR15126.
[13] RAO W, LI G, Wang X,. Parametric sparse representation method for ISAR imaging of rotating targets[J]., 2014, 50(2): 910-919. doi: 10.1109/TAES.2014. 120535.
[14] 顧福飛, 張群, 楊秋, 等. 基于NCS算子的大斜視SAR壓縮感知成像方法[J]. 雷達(dá)學(xué)報, 2016, 5(1): 16-24. doi: 10.12000 /JR15035.
GU Fufei, ZHANG Qun, YANG Qiu,. Compressed sensing imaging algorithm of high-squint SAR based on NCS operator[J]., 2016, 5(1): 16-24. doi: 10. 12000/JR15035.
[15] 李震宇, 梁毅, 邢孟道, 等. 彈載合成孔徑雷達(dá)大斜視子孔徑頻率相位濾波成像算法[J]. 電子與信息學(xué)報, 2015, 37(4): 954-959. doi: 10.11999/JEIT140618.
LI Zhenyu, LIANG Yi, XING Mengdao,. A frequency phase filtering imaging algorithm for highly squint missile-borne synthetic aperture radar with subaperture[J].&, 2015, 37(4): 954-959. doi: 10.11999/JEIT140618.
[16] 江淮, 趙惠昌, 漢敏, 等. 基于DFT濾波器組的大斜視SAR成像算法[J]. 電子與信息學(xué)報, 2016, 38(1): 104-110. doi: 10.11999/JEIT150381.
JIANG Huai, ZHAO Huichang, HAN Min,. Highly squint SAR imaging algorithm based on DFT filter banks[J].&, 2016,38(1): 104-110. doi: 10.11999/JEIT150381.
[17] 楊軍, 李震宇, 孫光才, 等. 一種新的大斜視TOPS SAR全孔徑成像方法[J]. 西安電子科技大學(xué)學(xué)報(自然科學(xué)版), 2015, 42(1): 42-48. doi: 10.3969/j.issn.1001-2400.2015.01.08.
YANG Jun, LI Zhenyu, SUN Guangcai,. Novel full aperture imaging algorithm for highly squinted TOPS SAR[J].(), 2015, 42(1): 42-48. doi: 10.3969 /j.issn.1001-2400.2015.01.08.
[18] FANG Jian, XU Zongben, ZHANG Bingchen,. Fast compressed sensing SAR imaging based on approximated observation[J]., 2014, 7(1): 352-363. doi: 10.1109/JSTARS.2013.2263309.
[19] KRAGH T J and KHARBOUCH A A. Monotonic iterative algorithm for minimum-entropy autofocus[C]. IEEE International Conference on Image Processing, Atlanta, 2006: 645-648.
陳一暢: 男,1988年生,博士生,研究方向為稀疏微波成像.
張 群: 男,1964年生,教授,博士生導(dǎo)師,主要研究方向包括雷達(dá)信號處理、電子對抗.
Parameter Estimation Method of Moving Targets with SAR Sparse Sampling Data
CHEN Yichang①②ZHANG Qun①③
①(,,710077,)②(,,100084,)③(,,200433,)
To solve the problem of motion parameter estimation of ground moving target, a parameter estimation method of moving targets with sparse sampling data of single SAR sensor is proposed. First, based on the 2 dimensional velocity of moving targets, an equivalent parametric model is constructed to transform the moving target echo into squint SAR echo. Then, with different parameters the modified iterative thresholding algorithm is applied to achieving imagery of moving target. Finally, the motion parameters of targets are obtained by minimizing the image entropy. It is shown that, using the proposed method, the required echo sampling can be reduced, the Doppler ambiguity problem can be avoided and accurate velocity estimation can be obtained even in low signal-to-clutter ration scenarios. Simulation results verify the effectiveness of the proposed method.
SAR; Parametric model; Moving target; Parameter estimation
TN957.51
A
1009-5896(2016)12-3049-07
10.11999/JEIT160922
2016-09-12;改回日期:2016-11-11;
2016-12-13
陳一暢 cyc_2007@163.com
國家自然科學(xué)基金(61631019, 61471386)
The National Natural Science Foundation of China (61631019, 61471386)