超高斯光束抽運調Q固體激光器仿真模型研究?

程秋虎 王石語 過振 蔡德芳 李兵斌

(西安電子科技大學物理與光電工程學院,西安 710071)

超高斯光束抽運調Q固體激光器仿真模型研究?

程秋虎?王石語 過振 蔡德芳 李兵斌

(西安電子科技大學物理與光電工程學院,西安 710071)

數值仿真,二極管抽運激光器,光束傳輸,Q開關

1 引 言

對于激光器的各種不同應用目的以及研究方向[1?7],激光光束的參數一直是需要考慮的重要因素.例如,輸出功率和總體效率是衡量激光系統的關鍵性指標.光束M2因子、束腰半徑、遠場發散角同樣是必須考慮的指標.通常,獲取激光系統參數的手段是通過實驗測量.但是如果能夠通過仿真對這些參數直接進行數值計算,則激光系統的設計和分析成本將大幅減小.

實際的激光系統是復雜的非線性開放系統,影響其穩定性的因素包括:抽運功率空間上的分布和時間上的起伏;光學元件的瑕疵和缺陷;晶體的熱效應和應力形變;多模振蕩、模式競爭等.這些因素可以通過Maxwell-Bloch方程進行描述.針對不同類型激光器的模擬仿真,計算結果已經充分驗證了該半經典理論的普遍性和精確性[8?12].但由于Maxwell-Bloch方程的復雜性,仿真計算需占用巨大的計算資源.而在特定的情況下,某些因素對激光系統的影響可以忽略.這時,使用比較簡易的模型對激光系統進行描述可以有效地縮短仿真計算的時間.

基于上述原因,本文提出了一種調Q固體激光器的仿真模型.假定任意的激光光束可由諧振腔的本征模式疊加構成.通過引入諧振腔的速率方程,將這一系列本征模式與激光器晶體的速率方程進行耦合.在耦合方程的基礎上,加入Q開關的時域調制.基于此模型,研究了抽運功率和抽運光場分布對主動調Q多模激光諧振腔的影響.在仿真計算中,抽運光束采用超高斯分布.仿真計算結果與對照實驗測量數據相符合,說明物理模型的有效性.為了全面分析調Q過程中的模式競爭,仿真計算了調Q脈沖產生過程,弛豫振蕩過程和連續振蕩過程,得出了脈沖功率和光束質量因子隨時間變化的規律.此規律可為激光器的設計和光學諧振腔的優化提供參考.

2 物理模型及研究方法

2.1 理論基礎

本文研究Nd:YAG固體激光器.其能級如圖1所示.

圖1 Nd:YAG晶體的能級圖Fig.1.Nd:YAG crystal’s energy level diagram.

粒子初始全部處在能級0,即基態.當Nd:YAG晶體被808 nm的抽運光照射后,一部分粒子吸收抽運光子躍遷到能級3.處在能級3的粒子非常不穩定,快速地躍遷到能級2.同樣,處在能級1的粒子也會以非輻射的形式快速地躍遷到基態.相對于能級2的壽命τf,能級3和能級1的粒子壽命非常短,可以近似認為能級3和能級1的粒子數為零.假設增益介質是均勻摻雜,濃度為Ntot,能級2的粒子濃度為N,則N同時是總的反轉粒子數濃度.基于上述近似,基態的粒子濃度可以表示為Ntot?N.

調Q技術是將激光器從連續運作模式轉換成脈沖模式,產生時域上脈寬短、功率高的激光脈沖.調Q通常分為兩個階段,Q開關關閉和Q開關打開.如圖2所示,Q開關關閉對應圖中的上能級粒子積累階段,Q開關打開對應圖中的脈沖產生階段.

圖2 調Q的時域圖Fig.2.Time scheme for Q-switching.

在上能級粒子積累階段,諧振腔處于高損耗狀態,增益介質被激光二極管持續抽運,反轉粒子數持續增長.脈沖產生階段,Q開關打開,諧振腔處于低損耗狀態,振蕩光迅速起振,反轉粒子數被消耗,產生持續時間短峰值功率高的激光脈沖.當Q開關再次關閉,激光器進入下一個調Q周期.

基于上述調Q原理,本文引入了Q開關的損耗,δq(t).因為δq(t)是隨時間變化的函數,所以通過設定δq(t)函數形式,則可以模擬Q開關的打開和關閉.理想的Q開關可以表示為

(1)式可以理解為:當Q開關關閉時,光束無法穿過Q開關,諧振腔被完全關閉;當Q開關打開時,Q開關完全透明,對穿過的光場無任何影響.

本文研究的激光器如圖3所示.

圖3 激光諧振腔的原理圖Fig.3.Illustration of the laser cavity.

在圖3中,激光諧振腔的腔長定義為Lc,Nd:YAG晶體的長度為Lg,Q開關的長度為Lq.為了方便表示三維空間積分,定義激光諧振腔所占的空間三維區域為?c,Nd:YAG晶體的空間三維區域為?g,Q開關的空間三維區域為?q.ng為增益晶體的折射率,nq為Q開關的折射率.則振蕩光在腔內往返一次所需的時間T,可以近似表示為

即往返一次總的光程差除以速度,其中c為真空中的光速.圖3中的固體激光器的鏡面和增益晶體的橫截面都為圓形,滿足圓對稱關系.為了描述圓對稱諧振腔的激光器,本文采用拉蓋爾-高斯光束表示諧振腔的光場模式.拉蓋爾-高斯光束的數學表達式為

Um,n(r,θ,z)是柱坐標下的拉蓋爾-高斯光束分布,m和n分別是角向參數和徑向參數;ω(z)為光束的光斑半徑,束腰ω0=ω(z=0)設定在輸出鏡面的位置;是拉蓋爾多項式;R(z)是拉蓋爾-高斯光束的等相位面曲率半徑;φ(z)是Gouy相移.

2.2 仿真模型

本文研究的是多模調Q激光器,近似假設每一個模式都是獨立地在諧振腔內進行振蕩.每一個模式的光子總數為m,n(t),則諧振腔內總的光子數為

因為(4)式表示整個諧振腔內的光子數,而本文需要建立諧振腔內部每個位置的速率方程,所以定義每一個模式下的歸一化光子數密度為|Um,n|2,其滿足 ∫

根據(4)和(5)式,空間上某一點的光子數密度可以表示成各個模式的光子數密度的疊加,即從(5)式可以得出,激光腔內所有光學元件的空間位置確定后,歸一化光子數密度是不隨時間變化的.所以總體光子數密度空間分布隨時間的變化,是由于每一個模式的光子總數m,n(t)隨時間變化引起的,即模式競爭.

圖3中,左邊藍色的為輸入鏡,對振蕩光全反射,反射率r1=1;右邊藍色的為輸出鏡,對振蕩光部分反射,反射率滿足0

其中Q開關的損耗δq(t)是時間的函數,所以光子壽命τi也是隨時間變化的函數.

基于上述定義,主動調Q多模固體激光器的速率方程可以表示為

式中,σ為增益晶體的受激發射截面,Rp是抽運速率.(7)式右邊第一項表示受激輻射對單個模式增加的光子數,第二項表示單個模式由于諧振腔內部損耗引起的光子數衰減.(8)式右邊第一項表示所有模式的受激輻射對上能級粒子數的消耗,第二項表示由自發輻射引起的上能級粒子數衰減,最后一項表示抽運過程.

本文研究的抽運光在增益晶體中的空間分布為超高斯分布,其滿足下面的通式:

其中η是總體的抽運效率,它包含了抽運光的耦合效率、增益介質的吸收效率、量子效率;Pp是抽運功率;υp是抽運光的頻率,υp=808 nm;h是普朗克常數;α是吸收系數,滿足α=σNtot;S是增益晶體的橫截面面積;ωp為高斯抽運光在晶體入射面上的光斑半徑;j是超高斯分布的階數,同時j必須是大于等于2的偶數.

兩個常見的超高斯分布分別為均勻分布的抽運光(高斯平頂光)和高斯分布的抽運光.均勻分布的抽運光對應j=∞,則(9)式就簡化成

高斯分布的抽運光對應j=2,則(9)式簡化成

(10)和(11)式都采用柱坐標的表達形式,其原點設在增益晶體入射面的中心.因為(10)式是均勻分布,所以Rp在增益晶體內部的任意橫截面上處處相等,Rp只隨z變化.而(11)式是高斯分布,Rp不僅隨z變化,也有徑向r上的變化.在圖4中,四條曲線分別代表當j=2,4,10,100時,超高斯分布的抽運光在晶體入射面上的歸一化強度分布.因為超高斯分布是以半徑旋轉對稱的,所以圖4只需給出徑向截面的強度分布.可以看出隨著j數值的增大,相對應的超高斯分布抽運光越發接近平頂高斯光.當j=100,此時超高斯分布可以在數值計算中近似的當作平頂高斯光處理.

圖4 歸一化超高斯分布Fig.4.Normalized super-Gaussian distribution.

3 數值方法及對照實驗

3.1 數值方法

模擬激光器的主動調Q過程,需要對(7)式和(8)式組合成的耦合方程進行數值仿真,所以需同時對空間和時間進行數值劃分.因為調Q激光器的研究重心在激光脈沖產生的階段,并且脈沖產生階段本身時間短,所以在數值模擬中,圖2中的脈沖階段需要短的時間間隔.而抽運階段和弛豫階段,由于沒有輸出功率,這時增大時間間隔,可節省計算時間.

因為仿真結果需和實驗數據進行對照,所以實驗測量的變量,如脈沖波形、輸出功率、光束質量因子M2都需要被已知參數表示出來.根據(6)式,激光諧振腔的輸出功率可解釋為輸出鏡對腔內光束的損耗,則單一光束模式的輸出功率可寫成

所以總體的輸出功率,就是每個模式輸出功率的疊加,即

對于單一模式的光束,其光束質量因子M2m,n滿足關系

最后,因為Pm,n(t),Pout(t)和都是和時間有關的變量,所以它們對應時間上的平均值可以表示為

式中的t1和t2是測量的起始時間和終止時間.

3.2 對照實驗

基于圖3,本文的對照實驗采用的是二極管抽運聲光調Q固體激光器.激光器的總體腔長為500 mm.抽運形式為單向端面抽運,抽運光的中心波長為808 nm.抽運光通過光纖耦合輸出進入晶體,在晶體內抽運光半徑ωp約為0.5 mm,抽運光分布近似為j=4的超高斯分布.增益介質使用的摻雜濃度為0.3%的Nd:YAG晶體,激光波長為λ=1064 nm.Nd:YAG晶體的長度為20 mm,橫截面直徑為d=3 mm.輸入鏡對波長1064 nm的激光全反,輸出鏡的反射率為95%.聲光調Q開關的重復頻率為1 kHz,占空比為1%.

表1 不同諧振腔模式的衍射損耗δm,nTable 1.Diffraction loss of different resonant modes δm,n.

針對上述具體實驗參數,可以通過公式N=(d/2)2/λLc計算出諧振腔的菲涅耳數約為4.23.在文獻[14]中,McCumber計算了擁有不同菲涅耳數的圓對稱激光諧振腔的各個模式的衍射損耗.結果表明,在相同菲涅耳數的諧振腔內,階數高的模式在腔內振蕩一次的損耗高,并且所有模式的損耗隨菲涅耳數的增加而減小.在本文第2節中定義了諧振腔對每個振蕩光模式的固有衍射損耗δm,n,根據文獻[14],可以計算出在菲涅耳數為4.23的諧振腔內各個模式的δm,n,如表1所列.

4 計算結果與討論

本節將仿真模擬的結果和實驗測量的結果進行數據對比.圖5中研究抽運功率對調Q過程的影響,給出了平均輸出功率ˉPout的仿真結果和實驗測量結果.隨著抽運功率的增加,仿真結果和測量結果都呈遞增的趨勢,但是遞增的斜率趨于平緩.這是因為在小功率抽運條件下,抽運功率值沒有超過大部分高階模式的振蕩閾值,使得激光器處在低階模式的工作狀態.而低階模式的衍射損耗δm,n小,所以在小功率抽運的條件下,激光器整體能量利用率較高,斜效率較高;當在大功率抽運條件下,抽運功率值超過大部分高階模式的振蕩閾值,使得激光器處在多模式的工作狀態,模式競爭導致高階模吸收利用了一部分低階模的抽運能量.根據表1,相對于低階模,高階模的衍射損耗δm,n大,所以高階模對抽運能量的利用率比低階模低.最終,激光器整體的斜效率被振蕩的高級模式拉低,并且此影響隨著抽運功率的增加而加重.

圖5 抽運功率和平均輸出功率的關系Fig.5.Relation between pump power and average output power.

圖6給出了抽運功率對脈沖寬度(FWHM)的影響.可以看出隨著抽運功率的增加,脈沖寬度呈遞減的趨勢.在抽運功率達到最大值28 W時,實驗測得的脈沖寬度為150 ns,仿真結果為175 ns.抽運功率越高,則在Q開關打開前,儲備的最大反轉粒子數濃度越高.當Q開關打開后,在具有高濃度的反轉粒子數的諧振腔內,光子受激放大過程所需的腔內往返振蕩次數較少,則脈沖的產生過程就會縮短,輸出的脈沖時域寬度也會縮短.

圖6 抽運功率和脈沖寬度的關系Fig.6.Relation between pump power and duration of pulse.

接下來研究抽運光分布對調Q過程的影響.假設抽運功率不變,改變抽運光的超高斯分布,即只改變超高斯分布階數j,仿真計算了不同超高斯分布抽運光對脈沖寬度和平均輸出功率的影響.如圖7所示,當抽運光的超高斯分布階數j逐漸增大時,平均輸出功率逐漸減小.當階數j超過6后,j的取值對平均輸出功率基本沒有影響.圖7說明了具有大階數超高斯分布的抽運光與高階模式振蕩光在晶體內部的重疊區域大,從而使更多的抽運能量被高階模式振蕩光吸收利用.因為高階模的衍射損耗δm,n大于低階模,所以高階模對抽運能量的能量利用率比低階模低,導致平均輸出功率降低.

圖7 j和平均輸出功率的關系Fig.7.Relation between j and average output power.

圖8給出超高斯分布階數j和脈沖寬度的關系.同樣是因為具有大階數的超高斯分布的抽運光與高階模式振蕩光在晶體內部的重疊區域大,所以隨著抽運光的超高斯分布階數j的逐漸增大,使更多的抽運能量被高階模式振蕩光吸收利用,導致脈沖產生過程中的模式競爭加強.強烈的模式競爭在一定程度一下抑制了低階模式的受激放大,而高階模的衍射損耗δm,n大于低階模,所以為了形成脈沖輸出,高階模需要較多次的腔內往返振蕩放大,最終結果是調Q脈沖的整體時間寬度被拉長.

圖8 j和脈沖寬度的關系Fig.8.Relation between j and duration of pulse.

為了更詳細地研究模式競爭對調Q過程的影響,計算了Q開關的一個完整周期,當Q開關打開后,使它一直保持打開的狀態.因此,一個仿真過程既包含了調Q脈沖產生過程,也包含了弛豫振蕩過程和連續振蕩過程.圖9給出了瞬時輸出功率隨時間的變化曲線.圖中t1時間點對應的就是調Q脈沖尖峰,t2到t5這段時間屬于弛豫振蕩過程,t5以后激光器穩定成連續工作模式.可以看出弛豫振蕩過程的輸出功率由開始脈沖間斷的形式,在t3時刻變成了連續的起伏形式,最終在t5時刻穩定下來.在設計連續工作模式的激光器時,一般是根據連續輸出功率值選取具有合適損傷閾值的光學元件.圖9中,t5時刻以后的連續輸出功率約為7.7 W,但t1時刻的瞬時輸出功率約為230 W,功率相差為30倍.基于上述原因,當分析連續工作模式激光器的光學元件損傷問題時,必須考慮弛豫振蕩過程,以最大瞬時輸出功率作為損傷閾值的參考值.

圖10中計算了輸出光場的瞬時光束質量因子M2和時間的關系.可以看出調Q脈沖光場的光束質量因子基本接近1.隨著激光器進入弛豫振蕩過程,光束質量因子也開始增大,并且其數值隨時間的推移進行著無規則的變化.最后當進入到連續工作模式時,光束質量因子M2才穩定下來.圖10說明,由于基模的衍射損耗最小,所以基模能夠在初始脈沖形成過程中抑制其他高階模起振,因此初始調Q脈沖的光場成分基本為激光器的基模;弛豫振蕩過程中,由于與基模光場重疊的晶體增益分布被初始脈沖大量消耗,而與基模光場非重疊的晶體增益分布還沒有被充分利用,因此能夠利用剩余增益分布的高階模開始起振,即諧振腔內開始模式競爭.所以此時段的光場成分包含了基模和高階模式,且模式之間的比例進行著無規則的變化;連續振蕩時,模式競爭達到了動態平衡,即光場的成分同時包含基模和高階模式,但是模式之間的比例是不隨時間變化的.

圖9 輸出功率隨時間的變化Fig.9.Change of output power with time.

圖10 M2隨時間的變化Fig.10.Change of M2with time.

圖11中,單獨計算了在t1,t2,t3,t4和t5這五個時刻的輸出光場強度的歸一化分布.從圖11可以看出,在t1時刻和t2時刻,光場強度分布的差異很小,此時段光場主要有諧振腔的基模構成.但當到了t3時刻,光場強度分布中心出現了空洞,說明此時由于存在模式競爭,光場成分中高階模的比例增加,而基模的比例下降.t4時刻,光場強度分布變為兩個花瓣狀.最終,光場強度分布在t5時刻穩定下來,此時的光場即為連續輸出光場,由基模和高階模共同構成.

圖11 輸出光束強度分布圖 (a)t1時刻;(b)t2時刻;(c)t3時刻;(d)t4時刻;(e)t5時刻Fig.11. Intensity distribution of output beam:(a)Time of t1;(b)time of t2;(c)time of t3;(d)time of t4;(e)time of t5.

5 結 論

本文提出了一種調Q固體激光諧振腔的仿真模型.假定任意的激光光束都可以由諧振腔的本征模式疊加構成,且諧振腔采用圓形鏡面,因此本文的諧振腔本征模式采用拉蓋爾-高斯光束.通過引入Q開關的時域調制函數,將這一系列本征模式與激光器晶體的速率方程耦合.基于此物理模型,研究了抽運功率和抽運光場分布對主動調Q多模激光諧振腔的影響.在仿真計算中,抽運光場在晶體內部的分布滿足超高斯分布.計算結果與對照實驗測量數據相符合,說明此物理模型的有效性.例如,隨著抽運功率的增加,激光器輸出功率增加,但斜效率減小.為了更加全面地分析調Q過程中的模式競爭問題,仿真計算了包含調Q脈沖產生過程,弛豫振蕩過程和連續振蕩過程的總體過程,得出了脈沖功率和光束質量因子隨時間變化的規律,此規律在分析激光器光學元件損傷問題時有一定參考價值.

[1]Yao Y H,Lu C H,Xu S W,Ding J X,Jia T Q,Zhang S A,Sun Z R 2014Acta Phys.Sin.63 184201(in Chinese)[姚云華,盧晨暉,徐淑武,丁晶新,賈天卿,張詩按,孫真榮2014物理學報63 184201]

[2]Wang X F,Wu Z M,Xia G Q 2016Acta Phys.Sin.65 024204(in Chinese)[王小發,吳正茂,夏光瓊2016物理學報65 024204]

[3]Mao Y F,Zhang H L,Xu L,Deng B,Sang S H,He J L,Xing J C,Xin J G,Jiang Y 2015Acta Phys.Sin.64 014203(in Chinese)[毛葉飛,張恒利,徐瀏,鄧波,桑思晗,何京良,邢冀川,辛建國,江毅2015物理學報64 014203]

[4]Zhu S S,Zhang S L,Liu W X,Niu H S 2014Acta Phys.Sin.63 064201(in Chinese)[朱守深,張書練,劉維新,牛海莎2014物理學報63 064201]

[5]Hou L,Han H N,Zhang L,Zhang J W,Li D H,Wei Z Y 2015Acta Phys.Sin.64 134205(in Chinese)[侯磊,韓海年,張龍,張金偉,李德華,魏志義 2015物理學報 64 134205]

[6]Sun Q,Yang Y,Deng Y Q,Meng F,Zhao K 2016Acta Phys.Sin.65 150601(in Chinese)[孫青,楊奕,鄧玉強,孟飛,趙昆2016物理學報65 150601]

[7]Dou Z Y,Tian J R,Li K X,Yu Z H,Hu M T,Huo M C,Song Y R 2015Acta Phys.Sin.64 064206(in Chinese)[竇志遠,田金榮,李克軒,于振華,胡夢婷,霍明超,宋晏蓉2015物理學報64 064206]

[8]Ge L,Chong Y,Stone A D 2010Phys.Rev.A82 063824

[9]Cerjan A,Chong Y D,Stone A D 2015Opt.Express23 6455

[10]Pick A,Cerjan A,Liu D,Rodriguez A W,Stone A D,Chong Y D,Johnson S G 2015Phys.Rev.A91 063806

[11]Cerjan A,Chong Y,Ge L,Stone A D 2012Opt.Express20 474

[12]Türeci H E,Ge L,Rotter S,Stone A D 2008Science320 643

[13]Wohlmuth M,P fl aum C,Altmann K,Paster M,Hahn C 2009Opt.Express17 17303

[14]McCumber D E 1965Bell Sys.Tech.J.44 333

Simulation model of super Gaussian beam pumpedQ-switched solid-state laser?

Cheng Qiu-Hu?Wang Shi-Yu Guo Zhen Cai De-Fang Li Bing-Bin

(School of Physics and Optoelectronic Engineering,Xidian University,Xi’an 710071,China)

28 February 2017;revised manuscript

2 June 2017)

Computer simulation is always an important means for studying laser,while laser theory is the basis of simulation.Although the semi-classical laser theory can accurately describe the generation process of laser,its complexity leads to a need of huge resources and time for computation.However,in particular cases,the in fluence of some factors on the laser system can be neglected.If a simpler model is employed to describe the laser system,the time of simulation can be shortened significantly.

In order to simulate the laser system more efficiently,a simulation model ofQ-switched solid-state laser is proposed in this paper.In this model,the time-domain function ofQswitch is introduced,which represents the modulation ofQswitch loss over time.Because the cross section of the Nd:YAG rod is circularly shaped,the resonator eigenmodes are assumed to be a Laguerre-Gaussian beam for simplicity.Then,any other laser beam can be formed by superposition of the eigenmodes of the resonator.These series of resonator eigenmodes are coupled with the rate equations of laser crystals.Finally,the distribution of pump light field inside the laser crystal is approximated as super Gaussian distribution.

Based on this physical model,the in fluence of pump power and pump light field distribution on the output beam of multimodeQ-switched solid-state laser is investigated.The simulation results are in good agreement with the experimental data,which explains the validity of the proposed model.For instance,with the increase of pump power,the output power of the laser increases,but the overall slope efficiency decreases.This is because the diffraction lossδm,nof the lower order mode is less than the diffraction loss of higher order mode.When the pumping power increases,the higher order mode that starts to oscillate has lower utilization efficiency of pump energy.Therefore,the overall slope efficiency of the laser is reduced.In order to analyze the mode competition in the multimodeQ-switched solid-state laser more comprehensively,the processes of laser pulse generation,relaxation oscillation and continuous oscillation are calculated as one full cycle.The laws of pulse power and beam quality factor versus time are obtained.For example,the maximum instantaneous output power of the relaxation oscillation is about 30 times the steady continuous output power.This law has a certain reference value when analyzing the damage threshold of laser optical element.In the pulse generation stage,the beam quality factor is close to 1,which explains the fact that the pulse field composition is nearly the fundamental mode of the laser.In the relaxation oscillation,the value of the beam quality factor changes irregularly with time,because mode competition is in a non-equilibrium state at this time.When stable continuous oscillation occurs,the mode competition achieves dynamic equilibrium,which means that the proportion of each mode is no longer changed in the output light field.

numerical simulation,diode-pumped lasers,wave propagation,Q-switching

PACS:02.60.Cb,42.55.Xi,42.25.Bs,42.60.GdDOI:10.7498/aps.66.180204

*Project supported by the National Defense Pre-Research Foundation of China(Grant No.9140A020105).

?Corresponding author.E-mail:chengqiouhu@126.com

(2017年2月28日收到;2017年6月2日收到修改稿)

仿真計算一直是研究激光器的重要手段,而激光理論是仿真計算的基礎.雖然半經典激光理論能夠精確地描述激光的產生過程,但是其復雜性導致仿真需要龐大的計算資源和計算時間.為了能夠更加高效地對激光器進行仿真,提出了一種調Q固體激光器的仿真模型.基于此模型,研究了影響主動調Q激光產生過程的因素.這些因素包括抽運功率、抽運光分布和模式競爭.仿真計算結果與對照實驗測量數據相符合,說明了模型的有效性.

10.7498/aps.66.180204

?國防預研究基金(批準號:9140A020105)資助的課題.

?通信作者.E-mail:chengqiouhu@126.com