呂貴洲,羅賢全
(軍械工程學(xué)院電子與光學(xué)工程系,河北石家莊050003)
基于OI-AGCD的含旋轉(zhuǎn)部件目標(biāo)ISAR成像
呂貴洲,羅賢全
(軍械工程學(xué)院電子與光學(xué)工程系,河北石家莊050003)
旋轉(zhuǎn)部件(如直升機(jī)的螺旋槳)嚴(yán)重影響逆合成孔徑雷達(dá)(inverse synthetic aperture radar,ISAR)的成像特性,基于優(yōu)化初值選擇的自適應(yīng)高斯包絡(luò)線性調(diào)頻基分解(optimized initial-adaptive Gaussian chirplet decomposition,OI-AGCD)設(shè)計(jì)了一種自適應(yīng)旋轉(zhuǎn)部件回波分離算法。建立了含旋轉(zhuǎn)部件目標(biāo)步進(jìn)頻率雷達(dá)ISAR成像模型,分析了旋轉(zhuǎn)部件與目標(biāo)主體回波特性差別?;诰嚯x-- 多普勒成像理論、主體和旋轉(zhuǎn)部件回波在高斯基上投影的差異,設(shè)計(jì)了基于OI-AGCD的旋轉(zhuǎn)部件回波分離算法。仿真數(shù)據(jù)驗(yàn)證結(jié)果表明了該算法去除旋轉(zhuǎn)部件回波的可行性,含旋轉(zhuǎn)部件類目標(biāo)ISAR成像質(zhì)量得到了有效提高。
逆合成孔徑雷達(dá);自適應(yīng)分解;時頻分析;直升機(jī)
逆合成孔徑雷達(dá)(inverse synthetic aperture radar,ISAR)通常用于固定平臺對非合作目標(biāo)進(jìn)行成像,如艦船、飛機(jī)、空間星體等。ISAR成像高的距離分辨力利用發(fā)射寬帶信號獲得,高的方位分辨率利用較大的角度觀測獲得。ISAR成像方法分為傅里葉變換法、現(xiàn)代信號譜估計(jì)法和非平穩(wěn)信號處理方法[12]。無論采用何種方法,對目標(biāo)的要求是小觀測角范圍內(nèi)散射點(diǎn)的位置、強(qiáng)度和特性不變,即目標(biāo)必須是鋼體[3]。然而,非鋼體運(yùn)動情況很多,如帶有螺旋槳的直升飛機(jī)、帶有掃描天線的艦船等都包含旋轉(zhuǎn)部件。對于這些情況,鋼體假設(shè)不再成立,旋轉(zhuǎn)部件回波會給目標(biāo)成像帶來負(fù)面影響,因此對旋轉(zhuǎn)部件回波的處理(識別、分離)是得到清晰像的必備步驟。____
針對帶有旋轉(zhuǎn)部件目標(biāo)的雷達(dá)成像,文獻(xiàn)[4]研究了旋轉(zhuǎn)部件目標(biāo)的包絡(luò)對齊問題,提出了通過一組回波形成穩(wěn)定的參考回波、弱化回波中螺旋槳子回波的影響,實(shí)現(xiàn)包絡(luò)對齊;由于旋轉(zhuǎn)部件的存在,目標(biāo)回波具有明顯的非平穩(wěn)特征,時頻分析方法被廣泛應(yīng)用到子回波分離中。文獻(xiàn)[5]采用Wigner-Ville分解(Wigner Ville decomposition,WVD)分析了直升機(jī)回波信號的調(diào)制特性,對單個旋轉(zhuǎn)點(diǎn)模型數(shù)據(jù)取得較好結(jié)果,但當(dāng)同時存在多個旋轉(zhuǎn)點(diǎn)時,WVD的交叉項(xiàng)問題難以解決。
為解決WVD等二次型時頻分析算法交叉項(xiàng)問題,自適應(yīng)分解方法被應(yīng)用到ISAR成像中。文獻(xiàn)[6]將自適應(yīng)高斯基分解(adaptive Gaussian basis representaion,AGR)算法[7]應(yīng)用到散射和非散射效應(yīng)的分離中,消除了非散射點(diǎn)在圖像中造成的模糊現(xiàn)象。
AGR算法采用的高斯基函數(shù)不含調(diào)頻項(xiàng),對包含調(diào)頻信息的信號缺乏有效的處理能力。針對該問題,文獻(xiàn)[8]將調(diào)頻項(xiàng)引入到高斯基函數(shù)中,形成了高斯包絡(luò)線性調(diào)頻基分解。
文獻(xiàn)[12,14]提出自適應(yīng)高斯包絡(luò)線性調(diào)頻基分解(adaptive Gaussian chirplet decomposition,AGCD)算法,利用最優(yōu)化方法解決多維參數(shù)的搜索問題,通過構(gòu)造和求解超越方程實(shí)現(xiàn)多維優(yōu)化搜索,得到基元Chirplet各參數(shù)解析解,在保證分解精度的同時,運(yùn)算量明顯降低,自適應(yīng)分解的實(shí)用性進(jìn)一步增強(qiáng)。時頻分辨率上AGCD算法表現(xiàn)出色,在可由調(diào)頻類信號建模的信號分析中分辨能力超強(qiáng)[12]。
對于在時頻平面上基函數(shù)深相交的復(fù)雜信號,AGCD由于初值選擇問題使得分解結(jié)果不合實(shí)際,文獻(xiàn)[13]提出了基于優(yōu)化初值選擇的AGCD算法,對初值進(jìn)行優(yōu)化,提高了精度,同時降低了運(yùn)算量。
現(xiàn)有的旋轉(zhuǎn)部件回波分離算法[910]主要針對單個旋轉(zhuǎn)點(diǎn)進(jìn)行處理,在現(xiàn)實(shí)情況下采用單個旋轉(zhuǎn)點(diǎn)對如槳葉,尤其是多個葉片的螺旋槳進(jìn)行建模顯然是不合實(shí)際的。由于多個旋轉(zhuǎn)點(diǎn)的回波非平穩(wěn)特征更加復(fù)雜,處理這種復(fù)雜信號AGCD算法效果較差。本文將基于優(yōu)化初值選擇的自適應(yīng)高斯包絡(luò)線性調(diào)頻基分解(optimized initial-adaptive Gaussian chirplet decomposition,OIAGCD)算法應(yīng)用到多旋轉(zhuǎn)點(diǎn)回波分離中,實(shí)現(xiàn)了多旋轉(zhuǎn)點(diǎn)回波信號分離,大大提高了帶有旋轉(zhuǎn)部件的復(fù)雜目標(biāo)的成像質(zhì)量。
針對自適應(yīng)分解中的運(yùn)算量問題,文獻(xiàn)[12]采用解四維變量方程組的方法獲得分解參量,極大地降低了運(yùn)算量。但是分解過程中的初值選擇對結(jié)果影響重大,對深相交復(fù)雜信號AGCD能量峰值策略并不適用,文獻(xiàn)[13]提出的OI-AGCD算法對其初值選擇方法進(jìn)行了補(bǔ)充,提高了AGCD處理復(fù)雜深相交信號的能力。
AGCD及OI-AGCD算法采用的基函數(shù)為復(fù)高斯包絡(luò)線性調(diào)頻信號
式中,tk為基函數(shù)的時間中心;ωk為基函數(shù)的頻率中心;信號寬度由αk決定;βk表示調(diào)頻斜率。OI-AGCD相對于AGCD具有更優(yōu)性能,采樣率1 MHz,分析信號
式中
該信號表示同一距離單元兩個對稱旋轉(zhuǎn)點(diǎn)的多普勒頻移特征,圖1表示不同時頻分解算法的分析結(jié)果。
圖1 復(fù)雜機(jī)動信號不同時頻分析方法分辨性能
從圖1可以看出,OI-AGCD對復(fù)雜變化規(guī)律多普勒信息提取具有更高的分辨力。WVD交叉項(xiàng)明顯,而AGCD對復(fù)雜信號不具備適應(yīng)性。
2.1 信號模型
在雷達(dá)成像中通常采用散射點(diǎn)模型對未知目標(biāo)進(jìn)行散射模型的建模。將目標(biāo)回波表示為一系列散射點(diǎn)回波的疊加,用M個獨(dú)立散射點(diǎn)對未知目標(biāo)進(jìn)行建模,則其回波表達(dá)式為
式中,信號E是頻率f和時間t的函數(shù);σm為第m個散射點(diǎn)的復(fù)散射系數(shù);Rm(t)為散射點(diǎn)的瞬時距離;f為載波頻率。假設(shè)目標(biāo)的運(yùn)動包含平動(t)和轉(zhuǎn)動θm(t),(xm,ym)表示散射點(diǎn)的位置,根據(jù)雷達(dá)成像中通常采納的目標(biāo)為鋼體的假定[10]可以得到
為便于分析帶有旋轉(zhuǎn)部件的復(fù)雜目標(biāo)回波,將目標(biāo)本身和旋轉(zhuǎn)部件分別假定為鋼體,并假定旋轉(zhuǎn)部件上的獨(dú)立散射點(diǎn)個數(shù)為N,得到如下的回波[10]表達(dá)式
式中,下標(biāo)B和R分別表示目標(biāo)主體和旋轉(zhuǎn)部分。顯然旋轉(zhuǎn)部分和目標(biāo)主體共享同一平動分量,對于主體部分,在雷達(dá)成像期間通常滿足小角度假設(shè),即
假設(shè)進(jìn)行成像處理之前對目標(biāo)回波進(jìn)行了標(biāo)準(zhǔn)的運(yùn)動補(bǔ)償,補(bǔ)償后的信號僅僅包含標(biāo)準(zhǔn)的旋轉(zhuǎn)運(yùn)動,即
這里用0代替目標(biāo)的平動分量,ωB為目標(biāo)的歸一化旋轉(zhuǎn)角速度,而旋轉(zhuǎn)部分的旋轉(zhuǎn)角仍然用θR(t)表示,這是因?yàn)樾D(zhuǎn)部件的運(yùn)動不再滿足小角度假設(shè)。將式(6)、式(7)代入式(5),得到
2.2 步進(jìn)頻率雷達(dá)旋轉(zhuǎn)部件回波仿真及分析
對步進(jìn)頻率雷達(dá)旋轉(zhuǎn)部件回波進(jìn)行仿真,參數(shù)和目標(biāo)模型如表1和圖2所示。
表1 步進(jìn)頻率雷達(dá)參數(shù)_
圖2 旋轉(zhuǎn)部件獨(dú)立散射點(diǎn)模型
由表1的雷達(dá)參數(shù)可以得到雷達(dá)的最遠(yuǎn)作用距離為1.5 km,最大非模糊距離窗為150 m,徑向距離分辨力為0.59 m。圖2中的4個旋轉(zhuǎn)部件上的旋轉(zhuǎn)點(diǎn)(a,b,c,d)圍繞旋轉(zhuǎn)中心o作逆時針旋轉(zhuǎn),旋轉(zhuǎn)角速度為6πrad/s。
考察圖2中旋轉(zhuǎn)點(diǎn)b、d的影響,這類似于偶數(shù)螺旋槳直升機(jī)上對稱槳葉的回波數(shù)據(jù)。以第1和第64個脈沖串的數(shù)據(jù)為例,第1個脈沖串形成的徑向距離像,第64個脈沖串發(fā)射時刻的橫向距離像、a,o,c 3點(diǎn)對應(yīng)距離單元的時頻分布,完整二維像如圖3所示。
從圖3(b)可以看出,由于旋轉(zhuǎn)部件的高速旋轉(zhuǎn)使得橫向距離像展寬到若干個距離單元,正是由于這種橫向距離像的展寬使得旋轉(zhuǎn)部件形成的二維像擴(kuò)展到橢圓狀區(qū)域,如圖3(d)所示。圖3(c)顯示的是兩個對稱旋轉(zhuǎn)點(diǎn)回波的多普勒信息的變化規(guī)律,其調(diào)制特征非常明顯。多個旋轉(zhuǎn)點(diǎn)的存在使得目標(biāo)的回波數(shù)據(jù)更加復(fù)雜,時頻平面上這種回波由多個正弦或余弦規(guī)律變化的多普勒時變曲線相互交叉,其信號回波的時變性更強(qiáng)。采用AGCD和OI-AGCD進(jìn)行100條高斯基分解,得到的時頻分布如圖4所示。
圖3 兩個對稱旋轉(zhuǎn)點(diǎn)回波分析
圖4 對稱旋轉(zhuǎn)點(diǎn)回波的AGCD及OI-AGCD時頻分布
對比圖4(a)和圖4(b),OI-AGCD的分解效果明顯強(qiáng)于AGCD,旋轉(zhuǎn)部件回波的高斯基參數(shù)具有如下特征:
(1)在多普勒信息變化斜率較大的階段,得到的高斯基的調(diào)頻斜率β較大,如圖3(c)中由A點(diǎn)到B點(diǎn)之間的變化曲線;
(2)多普勒信息變化的轉(zhuǎn)折點(diǎn)處,高斯基的中心頻率的模較大,如圖3(c)中A點(diǎn)和B點(diǎn)附近;
(3)旋轉(zhuǎn)部件回波的高斯基長度αk較小。
文獻(xiàn)[9- 10]以上述(1)(2)兩個特征設(shè)定準(zhǔn)則對旋轉(zhuǎn)部件回波信號進(jìn)行分離,本文加入高斯基長度判斷,進(jìn)一步去除旋轉(zhuǎn)部件弱回波信號。
3.1 旋轉(zhuǎn)部件回波分離的必要性
旋轉(zhuǎn)部件回波對目標(biāo)像形成干擾,造成目標(biāo)像模糊,給目標(biāo)成像和識別帶來巨大困難。抑制或消除旋轉(zhuǎn)部件回波是處理帶有旋轉(zhuǎn)部件目標(biāo)回波信號的重要手段,現(xiàn)有的算法對處理旋轉(zhuǎn)部件采用了單個旋轉(zhuǎn)點(diǎn)的假設(shè),雖然取得了較好的效果,但是在實(shí)際應(yīng)用中單旋轉(zhuǎn)點(diǎn)的假設(shè)過于簡單。為此,針對多個旋轉(zhuǎn)點(diǎn)情況提出基于OI-AGCD的旋轉(zhuǎn)部件回波分離算法。
3.2 分離算法
根據(jù)OI-AGCD算法對旋轉(zhuǎn)部件回波分解得到的高斯基參數(shù)的特征,旋轉(zhuǎn)部件回波分離的ISAR成像算法如下:
步驟1將回波數(shù)據(jù)按距離單元存儲構(gòu)成二維矩陣S,行為每一個發(fā)射脈沖按距離單元存儲的回波,列為同一距離單元不同發(fā)射脈沖的回波;
步驟2對S的每一列信號作能量剩余小于1%(作為結(jié)束條件)的高斯基自適應(yīng)分解,記分解次數(shù)為K,調(diào)頻斜率最大值為βmax=MAX(βi),i=1,2,…,K,則旋轉(zhuǎn)部件和主體部分與βmax所在的距離單元相同;
步驟3對D中包含旋轉(zhuǎn)部件距離單元(設(shè)其為第i個距離單元)回波的信號作K個高斯基自適應(yīng)分解,分解參數(shù)集合P={Pk,1≤k≤K},其中Pk=(Ak,αk,tk,ωk,βk);
從分解參數(shù)數(shù)組P中去除基函數(shù),得到新分解參數(shù)集
式中,α′0,ω′0和分別對應(yīng)高斯基信號的長度,中心頻率和調(diào)頻斜率的門限值;
步驟7對D1應(yīng)用傳統(tǒng)的距離多普勒成像算法,得到目標(biāo)二維像。
3.3 仿真分析
仿真中目標(biāo)本身及旋轉(zhuǎn)部件的幾何結(jié)構(gòu)如圖5所示。這里用9個“×”型散射點(diǎn)表示目標(biāo)本身,用2個“●”型散射點(diǎn)表示旋轉(zhuǎn)部件。仿真中采用的雷達(dá)參數(shù)如表1所示,目標(biāo)本身的旋轉(zhuǎn)角速度為0.02 rad/s,旋轉(zhuǎn)部件圍繞旋轉(zhuǎn)中心(散射點(diǎn)5)以5πrad/s的角速度轉(zhuǎn)動。
圖5 含旋轉(zhuǎn)部件的仿真目標(biāo)模型
采用傳統(tǒng)距離 多普勒方法得到的二維像及第170個距離單元(對應(yīng)旋轉(zhuǎn)中心)信號的時頻分布分別如圖6所示。
圖6 兩個對稱旋轉(zhuǎn)點(diǎn)的回波分析
從圖6(a)可以看出:中間距離單元3個獨(dú)立散射點(diǎn)被兩個旋轉(zhuǎn)點(diǎn)形成的多普勒噪聲所覆蓋,這極大地影響了目標(biāo)成像質(zhì)量,給目標(biāo)識別和成像帶來了負(fù)面影響。采用OI-AGCD旋轉(zhuǎn)部件回波分離算法,對第170個距離單元目標(biāo)本身和旋轉(zhuǎn)部件的回波數(shù)據(jù)進(jìn)行100個高斯基自適應(yīng)分解,得到時頻分布如圖7所示。
圖7 目標(biāo)散射點(diǎn)和旋轉(zhuǎn)部件回波OI-AGCD譜
主目標(biāo)(對應(yīng)于圖5中的4,5,6散射點(diǎn))和旋轉(zhuǎn)部件回波時頻分布特征差別,正是基于這種差別提出旋轉(zhuǎn)部件回波分離算法。第170個距離單元的復(fù)合目標(biāo)回波信號的時頻分布如圖8(a)所示,高斯基中心頻率門限值ω0=1.5 rad/s,調(diào)頻斜率β0=0.005 rad/s2,得到第170個距離單元分離出旋轉(zhuǎn)部件回波后剩余信號時頻分布,如圖8(b)所示,去除所有距離單元回波信號中旋轉(zhuǎn)點(diǎn)子回波,目標(biāo)二維像如圖8(c)所示。
圖8 第170距離單元復(fù)合目標(biāo)分離后的時頻分布及重構(gòu)二維像
比較圖2(a)和圖8(b)可以看出,經(jīng)過OI-AGCD分離算法得到的剩余回波和目標(biāo)本身的4、5、6 3個獨(dú)立散射點(diǎn)的回波非常接近,證明OI-AGCD算法分離多旋轉(zhuǎn)點(diǎn)回波數(shù)據(jù)是有效的。對比圖6(a)和圖8(c),將旋轉(zhuǎn)部件回波信號分離后,目標(biāo)本身成像質(zhì)量得到有效提高。
成像目標(biāo)包含旋轉(zhuǎn)部件時,相互干擾的旋轉(zhuǎn)部件和目標(biāo)本身回波極大地影響成像質(zhì)量和目標(biāo)識別效果。AGCD在提取旋轉(zhuǎn)部件回波時只能對單個旋轉(zhuǎn)點(diǎn)進(jìn)行處理,而實(shí)際情況采用一個旋轉(zhuǎn)點(diǎn)近似旋轉(zhuǎn)部件是不夠的。為了解決多旋轉(zhuǎn)點(diǎn)的回波分離問題,本文首先給出了帶有旋轉(zhuǎn)部件的目標(biāo)回波模型,優(yōu)化了旋轉(zhuǎn)部件回波分離準(zhǔn)則,提出了基于OI-AGCD的旋轉(zhuǎn)部件回波分離算法。多個仿真實(shí)例表明該算法在多旋轉(zhuǎn)點(diǎn)回波信號分離中的有效性,通過回波分離,目標(biāo)鋼體部分散射點(diǎn)的聚焦性能得到了提高,得到了高質(zhì)量的ISAR像。
[1]Wehner D R.High resolution radar[M].London:Artech House Publishers,1987.
[2]Walker J L.Range-Doppler imaging of rotating objects[J].IEEE Trans.on Aerospace and Electronic Systems,1980,16(1):23- 52.
[3]LüX G,Hao S Q,Leng J F,et al.Synthetic aperture radar imaging based on time frequency analysis[J].Electronic and Information Warfare Technology,2012,27(5):22- 26.(呂旭光,郝士琦,冷蛟鋒,等.基于時頻分析的合成孔徑激光雷達(dá)成像算法[J].電子與信息對抗技術(shù),2012,27(5):22- 26.)
[4]Lu G Y,Bao Z.Range alignment for targets with moving parts in ISAR imaging[J].Systems Engineering and Electronics,2000,22(6):12- 14.(盧光躍,保錚.ISAR成像中具有游動部件目標(biāo)的包絡(luò)對齊[J].系統(tǒng)工程與電子技術(shù),2000,22(6):12- 14.)
[5]LüJJ,Ding J J.Analysis on signal characters modulated by rotating structures of helicopter for air defense radar based on WVD[J].Radar and Countermeasure,2004,10(4):37- 39.(呂金建,丁建江.基于WVD的防空雷達(dá)直升機(jī)信號特征分析[J].雷達(dá)與對抗,2004,10(4):37- 39.)
[6]Sun Z Z,Chen Z P,Zhuang Z W,et al.A method of understanding and processing ISAR image based on time-frequency analysis[J].Journal of Electronics&Information Technology,2003,25(1):1- 8.(孫真真,陳曾平,莊釗文,等.一種基于時頻分解的ISAR圖像理解與處理方法[J].電子與信息學(xué)報(bào),2003,25(1):1- 8.)
[7]Qian S,Chen D.Signal representation using adaptive normalized Gaussian functions[J].Signal Processing,1994,36(6):1- 11.
[8]Qian S,Chen D,Qinye Y.Adaptive chirplet based signal approximation[C]∥Proc.of the IEEE International Conference on Acoustics,Speech and Signal Processing,1998,1781- 1784.
[9]Li J F,Ling H.ISAR Feature extraction from targets with nonrigid body motion using adaptive chirplet representation[C]∥Proc.of the IEEE Antennas and Propagation Society International Symposium,2002:294- 297.
[10]Li J,Ling H.Application of adaptive chirplet repre-sentation for ISAR features extraction from targets with rotating parts[J].Radar,Sonar Navigation,2003,150(4):284- 291.
[11]Stephane G,Mallat H.Matching pursuits with time-frequency dictionaries[J].IEEE Trans.on Signal Processing,1993,41(12):3397- 3415.
[12]Yin Q Y,Qian S,F(xiàn)eng A G.A fast refinement for adaptive Gaussian chirplet decomposition[J].IEEE Trans.on Signal Processing,2002,50(6):1298- 1306.
[13]LüG Z,He Q Wei Z S.Optimized initial-adaptive Gaussian chirplet signal decomposition[J].Signal Processing,2006,22(4):506- 510.(呂貴洲,何強(qiáng),魏震生.基于優(yōu)化初值選擇的自適應(yīng)高斯包絡(luò)線性調(diào)頻基信號分解[J].信號處理,2006,22(4):506- 510.)
[14]Feng A G,Yin Q Y,LüL.An adaptive fast algorithm for decomposition of signal based on Gaussian chirplet[J].Progress in Natural Science,2002,12(9):982- 988.(馮愛剛,殷勤業(yè),呂利.基于Gauss包絡(luò)Chirplet自適應(yīng)信號分解的快速算法[J].自然科學(xué)進(jìn)展,2002,12(9):982- 988.)
[15]Hagelen M T.ISAR-imaging of helicopters at millimeter-wave frequencies[J].Euro Radar Proceeding,2004,16(1),67- 70.
呂貴洲(197-7- ),男,講師,博士,主要研究方向?yàn)殡娮友b備測試與故障診斷、非平穩(wěn)信號處理、雷達(dá)信號處理及雷達(dá)成像。
E-mail:lgz2691@163.com
羅賢全(1977- ),男,講師,博士,主要研究方向?yàn)樾盘柼幚怼㈦娮友b備測試與故障診斷。
E-mail:luoweed@163.com
ISAR imaging of targets with rotating parts based on OI-AGCD
LüGui-zhou,LUO Xian-quan
(Electronic and Optical Department,Ordnance Engineering College,Shijiazhuang 050003,China)
Echoes produced by rotating parts,such as propellers of helicopters,seriously worsen inverse synthesizer aperture radar(ISAR)imaging quality.An algorithm based on optimized initial adaptive Gaussian chirplet decomposition(OI-AGCD)is proposed.Firstly,the ISAR imaging model which involves the target with rotating parts,based on stepped frequency radar,is established.The different performances between the main body and the rotating parts are analized.Secondly,based on range-Doppler imaging theroy and the different projection properties of main body and roating parts on Gaussian bases,the algorithm of cleaning rotating parts’radar echoes based on AGCD is proposed.Finally,simulation results show the effectiveness of the algorithm to eliminate echoes produced by rotating parts and high quality images can be achieved.
inverse synthesizer aperture radar(ISAR);adaptive decomposition;time-frequency analysis;helicopter
TN 911.7
A
10.3969/j.issn.1001-506X.2015.11.11
1001-506X(2015)11-2492-05
2014- 04- 22;
2015- 05- 28;網(wǎng)絡(luò)優(yōu)先出版日期:2015- 07- 07。
網(wǎng)絡(luò)優(yōu)先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20150707.1451.008.html
國家自然科學(xué)基金(61372039)資助課題