彭 科,胡 凡,張為華,武澤平
(國(guó)防科技大學(xué) 航天科學(xué)與工程學(xué)院, 湖南 長(zhǎng)沙 410073)
?
序列近似優(yōu)化方法及其在火箭外形快速設(shè)計(jì)中的應(yīng)用*
彭科,胡凡,張為華,武澤平
(國(guó)防科技大學(xué) 航天科學(xué)與工程學(xué)院, 湖南 長(zhǎng)沙410073)
摘要:針對(duì)序列近似優(yōu)化方法在代理模型構(gòu)造與采樣策略等方面的不足,基于采樣點(diǎn)局部密度,引入與局部密度成反比的樣本點(diǎn)影響體積概念,通過優(yōu)化總影響體積確定徑向基函數(shù)最優(yōu)核寬度,滿足序列近似優(yōu)化過程不同規(guī)模、非均勻樣本條件下的函數(shù)近似需要;建立潛在可行域最大距離加點(diǎn)準(zhǔn)則,并與潛在最優(yōu)加點(diǎn)準(zhǔn)則結(jié)合平衡算法的探索能力與開發(fā)能力;建立三步收斂判定準(zhǔn)則;構(gòu)建改進(jìn)序列近似優(yōu)化的算法流程。對(duì)于Golinski減速器的優(yōu)化設(shè)計(jì)問題,算法在目標(biāo)函數(shù)調(diào)用42次后便搜索到全局最優(yōu)解,體現(xiàn)了其良好的全局尋優(yōu)能力與搜索效率。以“天航二號(hào)”火箭為例,建立其外形優(yōu)化問題數(shù)學(xué)模型,所提優(yōu)化方法在調(diào)用原始計(jì)算模型165次之后便搜索到全局最優(yōu)解,大大提高了設(shè)計(jì)效率,同時(shí)飛行試驗(yàn)也表明設(shè)計(jì)結(jié)果滿足要求。
關(guān)鍵詞:序列近似優(yōu)化;代理模型;采樣策略;局部密度;收斂判定準(zhǔn)則;“天航二號(hào)”火箭;外形優(yōu)化
序列近似優(yōu)化(Sequential Approximate Optimization,SAO)方法基于少量初始采樣點(diǎn)構(gòu)造初始代理模型,采取一定的加點(diǎn)策略更新采樣點(diǎn),逐步提高代理模型對(duì)最優(yōu)解的近似精度,根據(jù)收斂準(zhǔn)則終止算法并輸出最優(yōu)解[1-4]。相比智能優(yōu)化算法,序列近似優(yōu)化方法可在保證全局最優(yōu)的前提下大幅降低原模型調(diào)用次數(shù),顯著提高優(yōu)化效率[1]。代理模型構(gòu)造與采樣策略是序列近似優(yōu)化方法的兩項(xiàng)關(guān)鍵技術(shù),國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了大量研究[3-7]。
徑向基插值模型在精度和魯棒性方面皆較為可靠[8],是序列近似優(yōu)化方法中被廣泛使用的代理模型,其核寬度的合理確定對(duì)近似精度有決定性影響。文獻(xiàn)[6]解決了采樣點(diǎn)較為均勻條件下的徑向基函數(shù)核寬度確定問題;文獻(xiàn)[4]提出了一種采樣點(diǎn)個(gè)數(shù)趨于無窮情況下的非均勻核寬度確定方法;文獻(xiàn)[1]提出了基于采樣點(diǎn)局部密度的核寬度確定方法,但算法中總影響體積這一關(guān)鍵參數(shù)依賴經(jīng)驗(yàn)確定,難以實(shí)現(xiàn)最優(yōu)的近似精度。序列近似優(yōu)化過程中樣本點(diǎn)數(shù)量由少逐漸增多,樣本分布較不均勻,建立不同規(guī)模、非均勻樣本的徑向基函數(shù)(Radial Basis Function, RBF)核寬度確定方法,對(duì)提高序列近似優(yōu)化過程代理模型精度、提高優(yōu)化效率具有重要意義。
合理的加點(diǎn)策略是平衡算法探索能力與開發(fā)能力以保證全局最優(yōu)前提下提高算法收斂速率的關(guān)鍵因素,目前常用的再采樣策略主要有:EI準(zhǔn)則[9]、潛在最優(yōu)加點(diǎn)準(zhǔn)則、最大曲率準(zhǔn)則、最小密度加點(diǎn)準(zhǔn)則等,現(xiàn)有準(zhǔn)則各有特點(diǎn),但皆難以同時(shí)保證優(yōu)化結(jié)果的全局性與收斂速率。
火箭外形設(shè)計(jì)問題是總體設(shè)計(jì)過程的重要問題,不斷提高外形設(shè)計(jì)效率是設(shè)計(jì)人員的不懈追求。本研究基于樣本點(diǎn)局部密度,引入與局部密度成反比的樣本點(diǎn)影響體積概念,通過優(yōu)化總影響體積確定徑向基函數(shù)最優(yōu)核寬度,并驗(yàn)證近似效果;建立潛在可行域最大距離加點(diǎn)準(zhǔn)則并與潛在最優(yōu)加點(diǎn)準(zhǔn)則結(jié)合,平衡序列近似優(yōu)化過程的收斂速度與全局搜索精度;建立收斂判定準(zhǔn)則;構(gòu)建改進(jìn)序列近似優(yōu)化的算法流程并通過算例分析算法性能;以“天航二號(hào)”(TianHang-Ⅱ, TH-Ⅱ)火箭為例,建立其外形優(yōu)化問題數(shù)學(xué)模型,采用序列近似優(yōu)化方法求解并分析結(jié)果。
1徑向基函數(shù)最優(yōu)核寬度確定方法
基本徑向基函數(shù)的數(shù)學(xué)模型為:
(1)
取基函數(shù)為Gauss函數(shù),即:
(2)
式中,σi為基函數(shù)的核寬度。
基本徑向基函數(shù)模型是一種非線性代理模型,對(duì)線性模型預(yù)測(cè)誤差較大,故將線性項(xiàng)引入徑向基代理模型,可以增加其線性泛化能力。其數(shù)學(xué)模型為:
(3)
式中,n為設(shè)計(jì)空間維數(shù),λj(j=0,1,2,…,n)為線性回歸系數(shù)。
式(3)中前兩項(xiàng)為樣本點(diǎn)的線性回歸模型:
(4)
(5)
此時(shí)將N個(gè)訓(xùn)練樣本代入式(5),求解線性方程組即可求得基函數(shù)的權(quán)系數(shù)wi(i=1,2,…,N)。
下面選取線性模型
y=x,-3 (6) 與非線性模型 (7) 對(duì)基本徑向基函數(shù)與含線性項(xiàng)徑向基函數(shù)的近似能力進(jìn)行比較,結(jié)果如圖1、圖2所示,后者對(duì)線性函數(shù)的近似能力顯著高于前者,對(duì)非線性函數(shù)的近似能力也略優(yōu)于前者,體現(xiàn)了其很好的線性和非線性泛化能力。 (a) 線性函數(shù)(a) Linear function (b) 非線性函數(shù)(b) Nonlinear function圖1 基本徑向基函數(shù)對(duì)測(cè)試函數(shù)的近似Fig.1 Elementary RBF for testing function approximation 式(2)中基函數(shù)的核寬度σi對(duì)代理模型預(yù)測(cè)精度有決定性影響,文獻(xiàn)[1]提出了基于樣本點(diǎn)局部密度的核寬度σi計(jì)算方法,局部密度定義如式(8)所示。 (a) 線性函數(shù)(a) Linear function (b) 非線性函數(shù)(b) Nonlinear function圖2 含線性項(xiàng)徑向基函數(shù)對(duì)測(cè)試函數(shù)的近似Fig.2 RBF containing linear terms for testing function approximation (8) (9) (10) 在樣本點(diǎn)分布較密的區(qū)域應(yīng)使每個(gè)樣本點(diǎn)的影響適當(dāng)減小,即核寬度σi應(yīng)相對(duì)較小,而在樣本點(diǎn)分布較為稀疏區(qū)域應(yīng)取相對(duì)較大的核寬度σi。 算法步驟如下[1]: 1)根據(jù)式(8)計(jì)算每個(gè)樣本點(diǎn)的密度ρ(xi); 2)按Vi/Vj=ρj/ρi計(jì)算每個(gè)樣本點(diǎn)的影響體積之比; Vt為所有樣本點(diǎn)的影響體積總和,其取值直接決定各樣本點(diǎn)對(duì)應(yīng)核寬度,從而對(duì)徑向基代理模型精度產(chǎn)生影響,文獻(xiàn)[1]中Vt取值根據(jù)經(jīng)驗(yàn)給出。圖3以函數(shù)f(x)=xsin(1.5x)為例,給出了10個(gè)隨機(jī)采樣點(diǎn)情況下不同Vt值對(duì)應(yīng)代理模型與原函數(shù)的對(duì)比情況;代理模型均方根誤差(Root Mean Square Error, RMSE)隨Vt變化規(guī)律如圖4所示。由圖可知,當(dāng)Vt過小或過大皆不能實(shí)現(xiàn)較好的近似效果,確定最優(yōu)Vt即可得到徑向基函數(shù)最優(yōu)核寬度,對(duì)提高代理模型精度具有關(guān)鍵意義。 圖3 不同Vt值對(duì)應(yīng)代理模型與原函數(shù)的對(duì)比Fig.3 Comparison between surrogate models applying different Vtvalues and original model 圖4 不同Vt值對(duì)應(yīng)代理模型均方根誤差Fig.4 RMSE of different Vt 最優(yōu)Vt值與樣本在設(shè)計(jì)空間的分布有關(guān),對(duì)樣本數(shù)量不敏感[1]。將樣本點(diǎn)均分為分區(qū)情況基本一致的兩份分別作為近似建模樣本與驗(yàn)證樣本,以驗(yàn)證樣本點(diǎn)代理模型預(yù)測(cè)值均方根誤差最小為目標(biāo),采用黃金分割法[10]迭代搜索最優(yōu)Vt。 近似建模樣本各點(diǎn)局部密度ρ(xi)確定后,黃金分割法搜索最優(yōu)Vt的單步計(jì)算過程如下: 1)根據(jù)Vt值計(jì)算近似建模樣本點(diǎn)的核寬度ri; 2)求解式(3)中λj與wi得代理模型s(x); 3)根據(jù)代理模型s(x)計(jì)算驗(yàn)證樣本點(diǎn)預(yù)測(cè)值; 4)驗(yàn)證樣本點(diǎn)計(jì)算預(yù)測(cè)值均方根誤差。 采用以上方法得圖3對(duì)應(yīng)樣本點(diǎn)最優(yōu)Vt值為1.17,對(duì)應(yīng)核寬度如圖5所示,代理模型與原函數(shù)的對(duì)比如圖6所示,二者均方根誤差為0.25,吻合較好,表明給出的徑向基函數(shù)最優(yōu)核寬度確定方法有效。 圖5 最優(yōu)樣本點(diǎn)核寬度Fig.5 Optimized kernel width of each sampling points 圖6 最優(yōu)Vt值對(duì)應(yīng)代理模型與原函數(shù)的對(duì)比Fig.6 Comparison between surrogate model applying optimized Vt and original model 2改進(jìn)序列近似優(yōu)化方法 2.1序列近似優(yōu)化方法基本流程 文獻(xiàn)[1]建立了如圖7所示的序列近似優(yōu)化方法基本流程,主要步驟為:①將設(shè)計(jì)空間線性映射到n維單位立方體內(nèi),采用拉丁超立方方法得到初始采樣點(diǎn),并計(jì)算采樣點(diǎn)目標(biāo)函數(shù)與約束值;②基于徑向基函數(shù)構(gòu)建目標(biāo)函數(shù)與各約束值的代理模型;③基于代理模型,求解樣本點(diǎn)更新的優(yōu)化問題;④收斂判斷。 圖7 序列近似優(yōu)化方法基本流程Fig.7 Flow chart of sequential approximate optimization 2.2采樣點(diǎn)更新準(zhǔn)則 一般優(yōu)化問題可表述為如式(11)所示形式: (11) 式中,i=0表示無不等式約束,j=0表示無等式約束。 最大距離加點(diǎn)準(zhǔn)則是序列近似優(yōu)化過程常用的加點(diǎn)準(zhǔn)則,可最大程度保證算法的探索能力。 令求解域中點(diǎn)x到所有樣本點(diǎn)的最小距離dmin(x)為: (12) (13) 為平衡算法的探索能力與開發(fā)能力,兼顧序列近似優(yōu)化過程的收斂速度與全局搜索精度,將潛在最優(yōu)加點(diǎn)準(zhǔn)則與潛在可行域最大距離加點(diǎn)準(zhǔn)則同時(shí)使用。潛在最優(yōu)加點(diǎn)將式(14)的最優(yōu)解作為新的樣本點(diǎn)加入到模型中。 (14) 采樣點(diǎn)更新涉及式(13)、式(14)定義的兩個(gè)優(yōu)化問題,皆采用粒子群優(yōu)化算法求解。 2.3收斂判定準(zhǔn)則 按以下3個(gè)步驟判斷算法是否收斂: (15) 式中,i=N-Sxl+1,N-Sxl+2,…,N-1。 (16) 以上三步皆滿足時(shí)則近似序列優(yōu)化收斂。 改進(jìn)后的序列近似優(yōu)化方法流程圖如圖8所示。 圖8 改進(jìn)序列近似優(yōu)化算法流程圖Fig.8 Flow chart of enhanced sequential approximate optimization 2.4算例分析 Golinski減速器優(yōu)化設(shè)計(jì)問題是美國(guó)國(guó)家航空航天局提供的一個(gè)多學(xué)科設(shè)計(jì)優(yōu)化測(cè)試問題[11],其目標(biāo)函數(shù)為重量最小,共有7個(gè)設(shè)計(jì)變量:齒輪面寬度x1、牙膜x2、小齒輪牙齒數(shù)目x3、軸1長(zhǎng)度x4、軸2長(zhǎng)度x5以及軸1與軸2的直徑x6和x7。其數(shù)學(xué)模型表述如式(17)所示。 (17) 式中: 各設(shè)計(jì)變量取值方位如表1所示。 表1 Golinski減速器設(shè)計(jì)變量取值范圍 注:圓點(diǎn)為可行解;三角形為不可行解。圖9 Golinski減速器優(yōu)化問題目標(biāo)函數(shù)與誤差收斂曲線Fig.9 Convergence curve of objective function and relative error for Golinski retarder optimization 所用方法在迭代至第4步,目標(biāo)函數(shù)調(diào)用15次之后便搜索到全局最優(yōu)解附近;迭代至17步,目標(biāo)函數(shù)調(diào)用42次后搜索到全局最優(yōu)解:x=[3.5,0.7,17,7.3,7.715 3,3.350 2,5.286 7]T;代理模型函數(shù)值與真實(shí)模型函數(shù)值的相對(duì)誤差僅為1.02×10-9。文獻(xiàn)[12]提出的并行模擬退火單純形算法是目前搜索效率相對(duì)很高的智能優(yōu)化算法,該方法搜索到Golinski減速器優(yōu)化問題最優(yōu)解附近約需調(diào)用1000次原始模型。本文方法能在保證全局搜索精度的前提下顯著減少模型調(diào)用次數(shù),可大幅提高火箭外形設(shè)計(jì)優(yōu)化問題求解效率。 3“天航二號(hào)”火箭外形快速設(shè)計(jì)優(yōu)化 3.1優(yōu)化問題 根據(jù)總體要求與任務(wù)特點(diǎn),“天航二號(hào)”火箭氣動(dòng)布局為“—·×”布局,總長(zhǎng)為2000 mm,頭部長(zhǎng)度為450 mm,發(fā)動(dòng)機(jī)長(zhǎng)度取為320 mm,彈翼、空氣舵皆為矩形,空氣舵緊貼儀器艙后端面安裝,如圖10所示。“天航二號(hào)”滿載質(zhì)量為60 kg,對(duì)應(yīng)質(zhì)心系數(shù)0.56,耗盡質(zhì)量為56 kg,對(duì)應(yīng)質(zhì)心系數(shù)為0.53。 圖10 “天航二號(hào)”火箭外形布局Fig.10 TH-Ⅱ rocket configuration 1)設(shè)計(jì)變量為:翼距頭部距離、翼弦長(zhǎng)、翼展長(zhǎng)、舵弦長(zhǎng)、舵展長(zhǎng),取值范圍依次為[1000, 1250]、[150, 190]、[350, 500]、[100, 150]、[150, 200](單位mm)。 綜上,“天航二號(hào)”火箭外形優(yōu)化問題數(shù)學(xué)模型為: (18) 式中,g4(x)為幾何約束。 3.2優(yōu)化結(jié)果與分析 優(yōu)化過程氣動(dòng)特性計(jì)算采用工程方法[13]完成,取初始樣本點(diǎn)為100個(gè),迭代至第33步,原始計(jì)算模型調(diào)用165次之后便搜索到全局最優(yōu)解,所用機(jī)時(shí)僅為190.6 s(AMD Athlon 四核處理器,主頻2.8 GHz),大大提高了設(shè)計(jì)效率。優(yōu)化過程目標(biāo)函數(shù)收斂曲線如圖11所示,相對(duì)誤差收斂曲線如圖12所示,外形優(yōu)化結(jié)果如圖13所示。需要說明的是,式(18)中約束條件與目標(biāo)函數(shù)計(jì)算分三步,依次為:先g1(x)和g2(x),然后g3(x)和g4(x),最后f(x),圖11中三角形對(duì)應(yīng)樣本點(diǎn)因其g2(x)值不滿足約束,無須計(jì)算目標(biāo)函數(shù)f(x)值,f(x)為人為給定值。 注:圓點(diǎn)為可行解;三角形與方行為不可行解。圖11 火箭外形優(yōu)化目標(biāo)函數(shù)收斂曲線Fig.11 Convergence curve of objective function for rocket shape optimization 圖12 火箭外形優(yōu)化相對(duì)誤差收斂曲線Fig.12 Convergence curve of relative error for rocket shape optimization 圖13 “天航二號(hào)”火箭外形優(yōu)化結(jié)果Fig.13 TH-II rocket optimized shape 圖14為采用本軟件系統(tǒng)設(shè)計(jì)優(yōu)化氣動(dòng)外形得到的“天航二號(hào)”火箭飛行試驗(yàn)概況,飛行試驗(yàn)取得圓滿成功,表明本文設(shè)計(jì)優(yōu)化得到的“天航二號(hào)”試驗(yàn)火箭靜穩(wěn)定性、操縱性、配平升阻力特性等性能指標(biāo)符合設(shè)計(jì)要求,滿足工程需要。 4結(jié)論 研究改進(jìn)了序列近似優(yōu)化方法,并將其應(yīng)用于火箭快速設(shè)計(jì),主要研究工作和結(jié)論如下: (a) 總裝后外形(a) Shape after assembling (b) 點(diǎn)火出架(b) Firing (c) 開傘回收(c) Recovery using parachute圖14 “天航二號(hào)”火箭飛行試驗(yàn)概況Fig.14 TH-Ⅱ rocket flight test 1)建立了含線性項(xiàng)的徑向基函數(shù),提高了其對(duì)線性函數(shù)的近似能力;基于樣本點(diǎn)的局部密度,采用黃金分割法得到了最優(yōu)總影響體積數(shù)值,進(jìn)而計(jì)算最優(yōu)核寬度,仿真結(jié)果表明了該方法的有效性。 2)建立了潛在可行域最大距離加點(diǎn)準(zhǔn)則并與潛在最優(yōu)加點(diǎn)準(zhǔn)則結(jié)合,兼顧了序列近似優(yōu)化算法的探索能力與開發(fā)能力;提出了序列近似優(yōu)化方法三步收斂判定準(zhǔn)則;構(gòu)建了改進(jìn)序列近似優(yōu)化的算法流程;通過Golinski減速器優(yōu)化設(shè)計(jì)問題分析了算法性能,目標(biāo)函數(shù)調(diào)用42次后便搜索到全局最優(yōu)解,體現(xiàn)了良好的全局尋優(yōu)能力與搜索效率。 3)以“天航二號(hào)”火箭為例,建立了其外形優(yōu)化問題數(shù)學(xué)模型,采用改進(jìn)的序列近似優(yōu)化方法迭代至第33步,原始計(jì)算模型調(diào)用165次之后便搜索到全局最優(yōu)解,所用機(jī)時(shí)僅為190.6 s,大大提高了設(shè)計(jì)效率。“天航二號(hào)”火箭飛行試驗(yàn)取得圓滿成功,表明設(shè)計(jì)優(yōu)化得到的“天航二號(hào)”試驗(yàn)火箭各項(xiàng)氣動(dòng)性能指標(biāo)符合設(shè)計(jì)要求,滿足工程需要。 參考文獻(xiàn)(References) [1]Wang D H, Wu Z P, Fei Y, et al. Structural design employing a sequential approximation optimization approach[J]. Computers and Structures, 2014, 134:75-87. [2]Kitayama S, Arakawa M, Yamazaki K. Sequential approximate optimization for discrete design variable problems using radial basis function network[J]. Applied Mathematics and Computation, 2012, 219(8): 4143-4156. [3]Nakayama H, Arakawa M, Sasaki R. Simulation-based optimization using computational intelligence[J]. Optimization and Engineering, 2002, 3(2): 201-214. [4]Kitayama S, Arakawa M, Yamazaki K. Sequential approximate optimization using radial basis function network for engineering optimization[J]. Optimization and Engineering, 2011, 12(4): 535-557. [5]Deng Y M, Lam I C, Tor S B, et al. A CAD-CAE integrated injection molding design system [J].Engineering with Computers, 2002, 18: 80-92. [6]Kitayama S, Yamazaki K. Simple estimate of the width in Gaussian kernel with adaptive scaling technique[J]. Applied Soft Computing, 2011, 11(8): 4726-4737. [7]Luo C T, Zhang S L, Wang C, et al. A metamodel-assisted evolutionary algorithm for expensive optimization[J]. Journal of Computational and Applied Mathematics, 2011, 236(5): 759-764. [8]Jin R, Chen W, Simpson T W. Comparative studies of metamodeling techniques under multiple modeling criteria[J]. Journal of Structural and Multidisciplinary Optimization, 2001, 23(1): 1-13. [9]Jones D R. A taxonomy of global optimization methods based on response surfaces[J]. Journal of Global Optimization, 2001, 21: 345-383. [10]謝政, 李建平, 陳摯. 非線性最優(yōu)化理論與方法[M]. 北京: 高等教育出版社, 2010. XIE Zheng, LI Jianping, CHEN Zhi. Nonlinear optimization theory and methods[M]. Beijing: Higher Education Press, 2010. (in Chinese) [11]Padula S L, Alexandrov N, Green L L. MDO test suite at NASA langley research center[C]//Proceedings of 6th AIAA/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, Bellevue, 1996. [12]Luo Y Z,Tang G J. Parallel simulated annealing using simplex method [J].AIAA Journal, 2006, 44(12): 3143-3146. [13]嚴(yán)恒元. 飛行器氣動(dòng)特性分析與工程計(jì)算[M]. 西安: 西北工業(yè)大學(xué)出版社, 1990. YAN Hengyuan. Aircraft aerodynamic characteristics analysis and engineering calculation[M]. Xi′an: Northwestern Polytechnic University Press, 1990. (in Chinese) Sequential approximate optimization method and its application in rapid design of rocket shape PENGKe,HUFan,ZHANGWeihua,WUZeping (College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China) Abstract:Sequential approximate optimization method has shortcomings in several respects, such as surrogate model establishing and infill strategy at present. Basing on local density of sampling points, the influence volume concept which is inversely proportional to local density was introduced and then the optimal kernel width of radial basis function was obtained by means of total influence volume optimization, thus, the function approximation needs in sequential approximate optimization process under the conditions of different scales and heterogeneous samples were satisfied. Potential feasible region infill strategy was proposed and potential optimal strategy was applied together, both exploration and exploitation capacity of the algorithm were satisfied. Three-step convergence criterion was set up. The algorithm flow process of sequential approximate optimization was constructed. For Golinski reducer optimization problem,the global optimal solution was solved after calculating original model 42 times, which embodied the good global optimization capacity and searching efficiency of the algorithm. Shape optimization mathematical model was established for TH-II rocket, global optimal shape was gained after 165 times of original model calling using the proposed method. The design efficiency was increased greatly and TH-II rocket aerodynamic shape was proved reliable by flight testing. Key words:sequential approximate optimization; surrogate model; infill strategy; local density; convergence criterion; TH-II rocket; shape optimization 中圖分類號(hào):V421.1;TP391.7 文獻(xiàn)標(biāo)志碼:A 文章編號(hào):1001-2486(2016)01-129-08 作者簡(jiǎn)介:彭科(1989—),男,四川鄰水人,博士研究生,E-mail:pengke_pk@163.com;張為華(通信作者),男,教授,博士,博士生導(dǎo)師,E-mail:zhangweihua@nudt.edu.cn 基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(51105368) *收稿日期:2015-02-03 doi:10.11887/j.cn.201601021 http://journal.nudt.edu.cn