王福來龐 晨李永禎王雪松②
①(國防科學(xué)技術(shù)大學(xué)電子信息系統(tǒng)復(fù)雜電磁環(huán)境效應(yīng)國家重點實驗室 長沙 410073)
②(國防科學(xué)技術(shù)大學(xué)理學(xué)院 長沙 410073)
一種同時全極化雷達正交多相編碼波形設(shè)計方法
王福來*①龐 晨①李永禎①王雪松①②
①(國防科學(xué)技術(shù)大學(xué)電子信息系統(tǒng)復(fù)雜電磁環(huán)境效應(yīng)國家重點實驗室 長沙 410073)
②(國防科學(xué)技術(shù)大學(xué)理學(xué)院 長沙 410073)
為了獲得精確的極化散射矩陣測量結(jié)果,同時全極化測量體制雷達對發(fā)射信號的正交性能提出了很高要求。傳統(tǒng)的設(shè)計方法得到的正交多相編碼波形正交性能受碼長的限制,同時對多普勒頻移比較敏感。該文提出一種具有較好多普勒容限的正交多相編碼波形優(yōu)化設(shè)計方法,針對目標勻速和勻加速運動狀態(tài),考慮波形的峰值旁瓣比和波形隔離度等指標,將波形設(shè)計問題轉(zhuǎn)化為非線性優(yōu)化問題,利用遺傳算法進行求解。仿真結(jié)果表明,相比于Deng,Khan等人提出的編碼,該文設(shè)計的正交多相編碼波形具有更好的多普勒容限,同時峰值旁瓣比和正交性能提升約為1.5~2 dB,能夠提高同時全極化測量體制雷達的測量精度。
雷達;同時全極化測量;正交多相編碼;多普勒容限
20世紀60年代以來,極化雷達經(jīng)歷了從非全極化測量體制到全極化測量體制的轉(zhuǎn)變,從分時全極化測量體制到同時全極化測量體制的發(fā)展[1,2]。相比于分時全極化測量體制,同時全極化測量體制雷達能夠在一個脈沖周期內(nèi)實現(xiàn)目標極化散射特性的有效測量[3],其核心思想是同時發(fā)射兩路正交波形,并對兩路正交回波進行同時接收。然而,在工程中滿足完全正交的波形難以找到,因此,設(shè)計具有良好正交性的波形對于同時全極化測量體制雷達具有重要的意義。
針對上述問題,國內(nèi)外學(xué)者對同時全極化測量雷達的發(fā)射波形設(shè)計進行了深入的研究。Giuli等人在1993年提出了利用兩路同時發(fā)射的正交波形進行極化散射矩陣測量的理論方法,并且以峰值旁瓣比(Peak Sidelobe Level,PSL)和波形隔離度(Isolation,I)等指標定量分析了正負斜率線性調(diào)頻波形的極化測量性能[2]。Babur等人針對調(diào)頻連續(xù)波波形進行了深入的研究,設(shè)計了準正交的調(diào)頻連續(xù)波波形[4]。調(diào)頻類波形由于樣式較為簡單,易被敵方截獲,同時不能較好地適應(yīng)各種不同的應(yīng)用場景,在一定程度上限制了其應(yīng)用。隨著數(shù)字信號處理技術(shù)的發(fā)展,編碼波形憑借其較好的抗截獲性能以及較高的波形設(shè)計效率等優(yōu)勢逐漸進入人們的視野。特別是隨著現(xiàn)代優(yōu)化理論與方法的飛速發(fā)展,使得設(shè)計具有較低峰值旁瓣比和波形隔離度的正交編碼波形成為可能。Deng等人在2004年首先利用模擬退火算法設(shè)計出了具有較好自相關(guān)和正交性能的正交多相編碼波形,并進行了性能分析[5]。電子科技大學(xué)的劉波等人利用遺傳算法深入研究了正交多相編碼波形和正交離散頻率編碼波形的設(shè)計問題,分析了碼長和相位編碼數(shù)對波形相關(guān)性能的影響[6]。Gao等人基于距離-多普勒準則分析了正交離散頻率編碼的波形設(shè)計問題[7]。Kansas大學(xué)的Z.Wang等人提出了一種基于重疊正交頻分復(fù)用(Interleaved Orthogonal Frequency Division Multiplexing,Interleaved OFDM)的正交編碼波形設(shè)計方法[8]。然而,在被測目標高速運動以及存在較大加速度的情況下,上述正交編碼波形的性能會發(fā)生退化,如回波脈壓輸出的旁瓣增加以及編碼波形的正交性能下降等。
為了充分考慮由目標運動產(chǎn)生的多普勒效應(yīng)對同時全極化測量波形性能的影響,本文在Deng,Khan等人工作的基礎(chǔ)上,采用靈活的子序列分塊編碼方法,增加編碼波形設(shè)計的自由度,綜合考慮波形的自相關(guān)和互相關(guān)性能設(shè)計聯(lián)合代價函數(shù),利用遺傳算法(Genetic Algorithm,GA)對正交多相編碼波形進行優(yōu)化設(shè)計。
同時全極化測量雷達的信號處理流程如圖1所示[9]:
圖1 同時全極化測量雷達信號處理流程圖Fig.1 Signal processing flow chart of simultaneous full-polarization radar
從圖1可以看出,同時全極化測量體制雷達通過雙極化通道同時發(fā)射兩路分別為H(水平)極化和V(垂直)極化的波形,并在接收端對兩路極化信號進行同時全極化接收,經(jīng)過匹配濾波檢測之后,對目標的極化散射矩陣進行反演。在窄帶條件下,對于點目標,天線的接收信號可以表示為[10]:
其中R和T分別表示天線的接收和發(fā)射方向圖,表征著收發(fā)天線的空域極化特性,上標T表示矩陣的轉(zhuǎn)置運算;eh(t)和ev(t)分別表示H極化通道和V極化通道的發(fā)射信號;S表示目標的極化散射矩陣。接收信號經(jīng)過匹配濾波之后的輸出為:
極化測量正交多相編碼波形可以表示為:
其中:
N表示碼長,表示子脈沖寬度,和分別表示H和V極化波形第n個子脈沖的相位。假設(shè)多相編碼信號的可用相位數(shù)為M,這里我們稱其為相位編碼數(shù),那么和只能從下面的集合中取值:
對于一個包含兩組碼長為N,相位編碼數(shù)為M的多相編碼信號集可以用矩陣表示為:
矩陣中的元素只能從集合式(5)中取值。根據(jù)信號的自相關(guān)和互相關(guān)函數(shù)定義[11],可以得到多相編碼信號的非周期自相關(guān)函數(shù)和互相關(guān)函數(shù)分別為:
需要指出的是,傳統(tǒng)的多相編碼波形設(shè)計方法沒有考慮目標運動產(chǎn)生的多普勒頻移對波形性能的影響,Khan等人在2006年提出了一種對多相編碼信號進行分塊編碼處理的波形設(shè)計方法[12,13],考慮了目標運動產(chǎn)生的多普勒頻移,設(shè)計的波形相比于傳統(tǒng)的多相編碼波形具有更好的多普勒容限,但是其設(shè)計方法在碼長較長時性能改善較小,波形的自相關(guān)和互相關(guān)性能相比傳統(tǒng)的多相編碼波形甚至出現(xiàn)了退化。
本文在Khan等人工作的基礎(chǔ)之上,為了進一步增加波形設(shè)計的自由度,以獲得更好的多普勒容限和正交性能,對多相編碼波形采用了如下更為靈活的分塊方式。將每行序列等間隔的分為G個子序列,則每個子序列的長度為N/G,將每個子序列的編碼相位進一步設(shè)計為:
對于傳統(tǒng)的碼長為N的M相編碼,每一個碼元的相位有M種取值,故其總的搜索空間為Q1=MN。對于本文的分塊多相編碼而言,將N長的多相碼分為G個子序列,每個子序列中的單個碼元共有M種相位取值,故對于每一個子序列而言其搜索空間為MN/G,根據(jù)式(9)可以看出,每個子序列的相位取值空間并不相同,故碼長為N的分塊多相編碼的總搜索空間為Q2=(MN/G)G,所以有如下等式成立:
故本文的分塊編碼波形與傳統(tǒng)的多相編碼波形總的搜索空間相同,即本文的分塊編碼波形設(shè)計方法僅僅增加了波形設(shè)計的自由度,而并未增加波形設(shè)計的復(fù)雜度。在此基礎(chǔ)上可以將多相編碼信號集的矩陣寫成:
綜上,得到了進行分塊編碼處理的多相編碼信號矩陣。
當目標高速運動時,會導(dǎo)致在接收端對信號進行匹配濾波接收時出現(xiàn)失配,嚴重時會出現(xiàn)漏檢的現(xiàn)象。假設(shè)目標相對雷達接收端的徑向速度與加速度為分別為v0和a,則信號不同時刻脈內(nèi)多普勒頻偏可以表示為:
將式(12)帶入式(14)可以得到:
根據(jù)式(15)可以發(fā)現(xiàn),在假設(shè)目標為勻加速運動的情況下,給定發(fā)射端的多相編碼波形的載波頻率、碼長和脈寬,則目標回波的多普勒相移差為等差序列。這也就意味著,當目標為勻速或者勻加速運動時,若能夠測得某些時刻的多普勒相移,則可以根據(jù)上述結(jié)論推出整個碼周期任意時刻的多普勒相位偏移,從而可以得出多相編碼波形所需要的最大多普勒容限,即設(shè)計的發(fā)射波形所需滿足的條件。在考慮目標的運動時,由于回波存在多普勒相移,需要對式(7)和式(8)進行修正,修正后的結(jié)果如下所示:
根據(jù)式(2),同時全極化測量雷達的極化散射矩陣測量精度取決于信號的互相關(guān)特性,即正交性;另外根據(jù)雷達信號處理理論,為了提高雷達的檢測性能,降低旁瓣對弱小檢測目標的影響,發(fā)射波形的自相關(guān)性能也十分重要。所以本文以最小化自相關(guān)函數(shù)峰值旁瓣比(PSL)和波形隔離度(I)[2],以及自相關(guān)函數(shù)旁瓣能量和互相關(guān)函數(shù)能量為優(yōu)化指標,設(shè)計的代價函數(shù)如式(18)所示:
其中:
最小化代價函數(shù)式(18)是一個非線性多變量的NP問題,全局搜索算法的計算量隨著碼長的增加呈指數(shù)增加;傳統(tǒng)的貪婪算法計算量小,但常會陷入局部最優(yōu)解,所以傳統(tǒng)的算法并不適用于此類優(yōu)化問題。遺傳算法是模擬自然界生物進化的“適者生存”原理而設(shè)計的,是一種針對多參數(shù),多目標同時優(yōu)化的快速算法,是解決問題的魯棒算法。針對本文式(18)的優(yōu)化問題,可以將其直接作為遺傳算法中的適應(yīng)度函數(shù),利用遺傳算法優(yōu)化尋找適應(yīng)度函數(shù)最大的編碼波形。算法處理的流程如圖2所示:
利用GA算法優(yōu)化設(shè)計多相編碼波形的步驟[14]如下所示:
(1) 參數(shù)編碼:就本文所優(yōu)化的問題而言,當給定多相編碼的可用相位數(shù)M時,則可以直接使用M進制的序列來進行編碼,然后根據(jù)給定的G,利用式(9)來計算出每個碼元對應(yīng)的相位。
(2) 初始群體的生成以及尺度因子W的賦值:由于遺傳算法的需要,必須產(chǎn)生一些由若干初始解構(gòu)成的初始群體,本文的初始解是隨機產(chǎn)生的。利用初始解估算出適應(yīng)度函數(shù)各個約束條件的數(shù)量級,兩兩相除得到數(shù)量級差異,然后以此來調(diào)整尺度因子W。
(3) 選擇算子:根據(jù)個體適應(yīng)度函數(shù)采用輪盤賭方法進行個體選擇。
(4) 交叉:交叉是產(chǎn)生新個體的重要手段,本文采用的是單點交叉方式。
(5) 變異:變異操作是逐位進行的,目的是充分挖掘個體多樣性。
圖2 GA算法流程圖Fig.2 Flow chart of GA algorithm
根據(jù)上文提出的算法,本文基于Matlab對算法進行了性能分析,其中仿真的輸入?yún)?shù)如表1所示:
表1 編碼波形優(yōu)化設(shè)計主要參數(shù)Tab.1 Main parameters for optimization of poly-phase codes
根據(jù)上述仿真參數(shù)可以得到碼長N=B·T=128,利用上述數(shù)據(jù)進行仿真得到結(jié)果如下所示。圖3所示為目標與接收端相對靜止時接收端匹配濾波的輸出結(jié)果,其中圖3(a)和圖3(b)分別為信號的自相關(guān)函數(shù)和兩路極化信號的互相關(guān)函數(shù)。從圖中可以看出,本文的優(yōu)化編碼相比于Deng等人[5]提出的編碼,自相關(guān)函數(shù)的峰值旁瓣比和互相關(guān)函數(shù)的波形隔離度均有一定性能的提升,性能改善約為2 dB左右。
圖4所示為目標相對于接收端做勻速運動時的仿真結(jié)果。當v0T=0.38時,相比于圖3可以明顯的看出信號的自相關(guān)和正交性能都出現(xiàn)了一定程度的退化,性能退化的程度與v0和T的絕對取值無關(guān),而是與二者的乘積相關(guān),同時本文的優(yōu)化編碼的性能退化程度要小于Deng等人[5]提出的編碼。
圖5所示為目標相對于接收端做勻加速運動時的仿真結(jié)果圖。同樣相比于圖3,當aT2=0.55時,兩種編碼波形的匹配輸出信號的性能出現(xiàn)了退化,退化的程度與a和T2的乘積相關(guān),與二者的絕對取值無關(guān)。同樣可以看出,本文編碼的自相關(guān)峰值旁瓣比和正交性能相比于Deng等人[5]提出的編碼約有2 dB左右的性能提升。
圖3 靜止目標仿真結(jié)果圖Fig.3 Simulation results of stationary target
圖4 勻速運動(v0T=0.38)目標仿真結(jié)果圖Fig.4 Simulation results of uniform motional (v0T=0.38) target
圖5 勻加速運動(aT2=0.55)目標仿真結(jié)果圖Fig.5 Simulation results of uniformly accelerated motional (aT2=0.55) target
圖6所示為接收端匹配濾波輸出信號的峰值旁瓣比和互相關(guān)峰值隨著目標運動速度增加而變化的關(guān)系圖。從圖6中可以看出隨著目標回波多普勒相移的增加,兩種信號的自相關(guān)和正交性能都出現(xiàn)了一定程度的下降,但是本文的編碼信號的性能下降更慢,相同條件下的性能優(yōu)于Deng編碼,即多普勒容限要優(yōu)于Deng編碼[5]。
上文分析了在碼長為N=128時本文算法與Deng等人[5]算法的性能差別,為了充分挖掘算法的性能,進一步改變仿真參數(shù),對比Deng[5]以及Khan[12]等人算法性能結(jié)果如表2,表3所示。其中PSL =(PSLh+ PSLv),GA算法參數(shù)同表1,信號帶寬B=1 MHz:
圖6 相關(guān)函數(shù)隨目標運動速度增加變化仿真結(jié)果圖Fig.6 Simulation results of variation chart of correlation function with increase of speed
表2N=40時算法仿真結(jié)果(dB)Tab.2 Simulation results ofN=40 (dB)
綜合上述不同碼長情況下的仿真結(jié)果可以看出,在碼長較短時,若目標處于靜止狀態(tài),則本文算法與Deng等人的算法性能幾乎沒有差別,但是優(yōu)化波形的正交特性優(yōu)于Khan編碼[12];隨著目標運動速度的增加,本文編碼與Khan編碼的自相關(guān)性能幾乎相同,但是要優(yōu)于Deng編碼,而正交性能方面,本文編碼為三者之中性能最佳的編碼。隨著碼長的增加,波形的自相關(guān)性能和正交性能均有一定程度的提升,文獻[15]指出,編碼波形的自相關(guān)性能和互相關(guān)性能隨碼長的增加呈指數(shù)關(guān)系變化。在長碼長時從仿真結(jié)果可以看出,無論目標處于靜止、勻速運動或勻加速運動狀態(tài),本文優(yōu)化編碼的自相關(guān)和正交性能相比于Deng編碼[5]都得到了改善,性能提升約為1.5~2 dB,另外Khan本人也指出其算法在碼長較長時性能惡化較嚴重[12],這里就不在分析其長碼時的性能。
表3N=512時算法仿真結(jié)果(dB)Tab.3 Simulation results ofN=512 (dB)
針對同時全極化測量雷達,本文提出了一種正交多相編碼波形優(yōu)化設(shè)計方法,深入分析了多普勒頻移對極化測量的影響,并以此為基礎(chǔ),基于分塊編碼的方式,針對不同的目標運動狀態(tài),利用遺傳算法對代價函數(shù)進行優(yōu)化。最終設(shè)計得到了具有更好自相關(guān)與互相關(guān)性能的編碼波形,相比于Deng,Khan等人[5,12]提出的編碼波形性能得到了改善,在長碼長時約有1.5~2 dB左右的性能提升,同時該波形的多普勒容限也有一定的改善。
但是需要指出的是,本文采用的分塊編碼方式為等間隔劃分子塊,如果采用不等間隔劃分可以獲得更大的波形設(shè)計自由度,使得編碼更加靈活。另外,增加波形的碼長也可以使信號的自相關(guān)和正交性能有很大的提升,但是通過文中對算法總搜索空間的分析可以看到,碼長的增加會使搜索空間大小呈指數(shù)規(guī)律上升,因此,搜索空間的快速膨脹給設(shè)計性能更為優(yōu)異的多相編碼波形帶來了新的挑戰(zhàn)。另外本文針對單目標代價函數(shù)對波形進行了優(yōu)化設(shè)計,實際應(yīng)用中除了波形的自相關(guān)以及互相關(guān)性能外,例如不同波形對應(yīng)的目標檢測概率也是我們所關(guān)心的問題[16],而針對多代價函數(shù)的優(yōu)化相比于單目標代價函數(shù)優(yōu)化也更為復(fù)雜,因此如何更為高效的設(shè)計復(fù)雜的編碼波形將是下一步研究工作的重點。
[1]Melvin L S and Banner G P.Radars for the detection and tracking of ballistic missiles,satellites,and planets[J].Lincoln laboratory Journal,2000,12(2): 217–244.
[2]Giuli D,Fossi M,and Facheris L.Radar target scattering matrix measurement through orthogonal signals[J].Radar and Signal Processing,IEE Proceeding F,1993,140(4):233–242.
[3]李永禎,王雪松,曾永虎,等.寬帶雷達目標的瞬態(tài)極化時頻分布表征[J].雷達科學(xué)與技術(shù),2004,2(6): 321–326.Li Yong-zhen,Wang Xue-song,Zeng Yong-hu,et al..Time-Frequency Distribution Characterization of Radar Targets’Instantaneous Polarization Scattering[J].Radar Science and Technology,2004,2(6): 321–326.
[4]Babur G,Krasnov O,Yarovoy A,et al..Nearly Orthogonal Waveforms for MIMO FMCW Radar[J].IEEE Transactions onAerospace and Electronic Systems,2013,49(3):1426–1437.DOI: 10.1109/TAES.2013.6557996.
[5]Deng H.Polyphase code design for orthogonal netted radar systems[J].IEEE Trans on Signal Processing,2004,52(11):3126–3135.DOI: 10.1109/TSP.2004.836530.
[6]Liu B.Orthogonal Discrete Frequency-Coding Waveform Set Design with Minimized Autocorrelation Sidelobes[J].IEEE Transactions on Aerospace and Electronic Systems,2009,45(4): 1650–1657.DOI: 10.1109/TAES.2009.5310326.
[7]Gao C,Teh K C,and Liu A.Orthogonal Frequency Diversity Waveform with Range-Doppler Optimization for MIMO Radar[J].IEEE Signal Processing Letters,2014,21(10): 1201–1205.DOI: 10.1109/LSP.2014.2329944.
[8]Wang Z,Tigrek F,Krasnov O,et al..Interleaved OFDM Radar for Simultaneous Polarimetric Measurement[J].IEEE Transactions on Aerospace and Electronic Systems,2012,48(3): 2085–2099.DOI: 10.1109/TAES.2012.6237580.
[9]Brunkow D A,George J,Bring V N,et al..Recent data system and antenna upgrades to the CSU-CHILL radar[C].AMS 32nd Conference on Radar Meteorology,Albuquerque,New Mexico,USA,2005.
[10]Xiao Jin-jun and Arye Nehorai.Joint Transmitter and Receiver Polarization Optimization for Scattering Estimation in Clutter[J].IEEE Transaction on Signal Processing,2009,57(10): 4142–4147.DOI: 10.1109/TSP.2009.2022887.
[11]Nadav Levanon and Eli Mozeson.Radar Signals[M].Hoboken,New Jersey: John Wiley & Sons,2004: 101–103.
[12]Khan H A and Edwards D J.Doppler problems in orthogonal MIMO radar[C].IEEE Conference on Radar,Verona,NY,USA,April 24–27,2006,Vol.1: 24–27.
[13]Frank R L.Polyphase Complementary-Codes[J].IEEE Transactions on Information Theory,1980,26: 641–647.DOI: 10.1109/TIT.1980.1056272.
[14]Lei Y J,Zhang S W,Li X W,et al..MATLAB GA Algorithm Toolbox and It’s Applications[M].Xi’an: Xidian Press,2005: 66–70.
[15]劉波.MIMO雷達正交波形設(shè)計及信號處理研究[D].[博士論文],電子科技大學(xué),2008.Liu Bo.Study on orthogonal waveform design and signal processing for MIMO radar[D].[Ph.D.dissertation],University of Electronic Science and Technology of China,2008.
[16]王璐璐,王宏強,王滿喜,等.雷達目標檢測的最優(yōu)波形設(shè)計綜述[J].雷達學(xué)報,2016,5(5): 487–498.Wang Lulu,Wang Hongqiang,Wang Manxi,et al..An Overview of Radar Waveform Optimization for Target Detection[J].Journal of Radars,2016,5(5): 487–498.
Orthogonal Polyphase Coded Waveform Design Method for Simultaneous Fully Polarimetric Radar
Wang Fulai①Pang Chen①Li Yongzhen①Wang Xuesong①②
①(State Key Laboratory of Complex Electromagnetic Environment Effects on Electronics and Information System,National University of Defense Technology,Changsha410073,China)
②(School of Science,National University of Defense Technology,Changsha410073,China)
To obtain an accurate polarization scattering matrix,simultaneous full polarization radar systems must transmit two signals.The performance of orthogonal polyphase codes designed by the traditional method is limited by the code length and is sensitive to Doppler frequency.In this paper,we propose a design method for orthogonal polyphase codes that have good Doppler tolerance.We consider the peak sidelobe level and isolation and transform the signal design problem into a nonlinear optimization problem,which we solve using a genetic algorithm.Our simulation results show that our proposed orthogonal polyphase codes have better Doppler tolerance and their peak sidelobe levels and orthogonal performances are 1.5~2 dB better than the codes designed by Deng or Khan.As such,the new design can improve the measurement accuracy of simultaneous full polarization radar systems.
Radar; Simultaneous full polarization; Orthogonal poly-phase codes; Doppler tolerance
s: The National Natural Science Foundation of China (61490690,61490694,61501478)
TN95
A
2095-283X(2017)04-0340-09
10.12000/JR16150
王福來,龐晨,李永禎,等.一種同時全極化雷達正交多相編碼波形設(shè)計方法[J].雷達學(xué)報,2017,6(4):340–348.
10.12000/JR16150.
Reference format:Wang Fulai,Pang Chen,Li Yongzhen,et al..Orthogonal polyphase coded waveform design method for simultaneous fully polarimetric radar[J].Journal of Radars,2017,6(4): 340–348.DOI: 10.12000/JR16150.
王福來(1993–),男,遼寧人,國防科技大學(xué)電子科學(xué)與工程學(xué)院在讀碩士生,主要研究方向為極化信息處理、雷達目標分辨與識別技術(shù)。
E-mail: wangfulai12@nudt.edu.cn
龐 晨(1986–),男,湖北人,博士,國防科技大學(xué)電子科學(xué)與工程學(xué)院講師,主要研究方向為極化信息處理、雷達目標分辨與識別技術(shù)。
E-mail: pangchen@nudt.edu.cn
李永禎(1977–),男,內(nèi)蒙古人,博士,國防科技大學(xué)電子科學(xué)與工程學(xué)院研究員,主要研究方向為雷達極化信息處理、空間電子對抗、目標檢測與識別。
王雪松(1972–),男,內(nèi)蒙古人,博士,教授,國防科技大學(xué)理學(xué)院院長,主要研究方向為極化信息處理、雷達目標識別、新體制雷達技術(shù)。
2016-12-22;改回日期:2017-02-20;網(wǎng)絡(luò)出版:2017-03-31
*通信作者: 王福來 wangfulai12@nudt.edu.cn
國家自然科學(xué)基金(61490690,61490694,61501478)