基于改進(jìn)A*算法的滑翔飛行器軌跡規(guī)劃

盧 青, 周 軍, 呼衛(wèi)軍

(西北工業(yè)大學(xué)精確制導(dǎo)與控制研究所, 陜西 西安 710072)

?

基于改進(jìn)A*算法的滑翔飛行器軌跡規(guī)劃

盧 青, 周 軍, 呼衛(wèi)軍

(西北工業(yè)大學(xué)精確制導(dǎo)與控制研究所, 陜西 西安 710072)

為了研究滑翔飛行器再入過(guò)程中的側(cè)向機(jī)動(dòng)問(wèn)題,A*算法的相關(guān)理論被引入到軌跡規(guī)劃中。飛行距離較遠(yuǎn)時(shí),地球的曲率必須加以考慮,且再入過(guò)程中飛行器會(huì)受到各種過(guò)程約束。針對(duì)現(xiàn)有規(guī)劃方法中在平面上進(jìn)行A*算法推演的局限性,本文提出了一種基于橢球面A*算法的軌跡規(guī)劃方法。通過(guò)推導(dǎo)將過(guò)程約束引入A*算法搜索過(guò)程中,從而規(guī)劃出合理的軌跡。仿真結(jié)果表明,本文所提出的方法能夠快速地對(duì)再入段的軌跡進(jìn)行有效規(guī)劃,具有很強(qiáng)的實(shí)用性。

A*算法; 滑翔飛行器; 軌跡規(guī)劃; 橢球面

0 引 言

再入段是航天器利用地球大氣層這一天然資源,使再入器減速下降,并消耗它具有的巨大能量。再入飛行器在再入飛行的過(guò)程中還需要同時(shí)考慮路徑點(diǎn)和禁飛區(qū)等實(shí)際約束的突防優(yōu)化問(wèn)題。

隨著高超聲速飛行器發(fā)展的深入進(jìn)行,飛行器橫向的軌跡規(guī)劃正在越來(lái)越多地引起國(guó)內(nèi)外學(xué)者的廣泛關(guān)注。2007年Jorris T R在文獻(xiàn)[1]中首次考慮高超聲速滑翔飛行器在再入過(guò)程中通過(guò)航路點(diǎn)和回避禁飛區(qū)的側(cè)向軌跡優(yōu)化策略。在這之后的兩年Jorris T R等人又先后提出了二維[2]的以及三維[3]的滿足航路點(diǎn)和禁飛區(qū)約束的再入飛行器軌跡規(guī)劃方法。這些方法并沒(méi)有對(duì)飛行過(guò)程中的約束給出明確處理,同時(shí)沒(méi)有考慮多個(gè)禁飛區(qū)的情形。

在計(jì)算機(jī)技術(shù)不斷發(fā)展的基礎(chǔ)上,以A*算法為代表的啟發(fā)式搜索方法越來(lái)越多地被運(yùn)用到軌跡規(guī)劃的過(guò)程中。文獻(xiàn)[4]采用了線性的實(shí)時(shí)A*算法來(lái)搜尋可行的路徑樹。最終得到滿足目標(biāo)要求的最小最短路徑。Richard N D則通過(guò)對(duì)標(biāo)準(zhǔn)的A*算法進(jìn)行次優(yōu)化修改可以在少量增加路徑代價(jià)的基礎(chǔ)上獲得可行的次優(yōu)軌跡,這就大大提高了軌跡規(guī)劃的效率[5],A*算法的運(yùn)用縮短了軌跡規(guī)劃的時(shí)間,但是并沒(méi)有將飛行器在飛行過(guò)程中受到的各種過(guò)程約束考慮到規(guī)劃的過(guò)程中。文獻(xiàn)[6]通過(guò)使用變步長(zhǎng)的A*算法對(duì)飛行器飛行過(guò)程中的各種約束條件進(jìn)行了考慮。這些利用A*算法進(jìn)行軌跡規(guī)劃過(guò)程都是以平面規(guī)劃為基準(zhǔn),當(dāng)飛行器的飛行距離較遠(yuǎn)時(shí),顯然不能繼續(xù)使用平面的相關(guān)計(jì)算方式。為此,本文考慮各種過(guò)程約束條件,將飛行器可達(dá)區(qū)域的推演引入到橢球面A*算法的搜索過(guò)程中,從而完成完整的軌跡規(guī)劃過(guò)程。

1 問(wèn)題描述

1.1 滑翔飛行器數(shù)學(xué)模型

忽略發(fā)動(dòng)機(jī)推力,得到飛行器在再入過(guò)程中半速度坐標(biāo)系下較為完整的數(shù)學(xué)模型[7-8]為

(1)

式中,ωe為地球自轉(zhuǎn)角速度;r為地心矢;θT為彈道傾角;σT為彈道偏角;γv為傾側(cè)角;gr為地心矢方向上的引力加速度;f1,f2,f3,f4,f5,f6為右函數(shù);V為速度;λ為經(jīng)度;φ為緯度。Y=CLSρV2/2,X=CxSρV2/2,Z=CxSρV2tanγv/2,分別是升力、阻力、側(cè)向力,ρ為大氣密度,S為參考面積。

1.2 軌跡約束

利用再入飛行器在受到的各種約束與高度的關(guān)系得到高度-速度剖面,建立再入走廊[9]:

(2)

2 算法流程

2.1 任務(wù)與假設(shè)

再入飛行器的主要任務(wù)是在上述各種過(guò)程約束的限制下完成整個(gè)軌跡規(guī)劃過(guò)程。為保證推演的嚴(yán)謹(jǐn)性作出如下假設(shè):

假設(shè) 1 禁飛區(qū)為具有特定半徑的圓域,表示為(L,B,R),L，B分別表示威脅區(qū)中心的經(jīng)度、緯度,R表示威脅區(qū)的半徑。

保證了算法推演的可行性,其他形狀可以被近似為圓形,符合實(shí)際。

假設(shè) 2 僅考慮地心矢徑方向的地球引力。

對(duì)數(shù)學(xué)模型進(jìn)行了簡(jiǎn)化,方便后續(xù)推導(dǎo)。

假設(shè) 3 再入過(guò)程中速度變化較小。

符合再入過(guò)程的特點(diǎn),可將速度表示為航程的線性函數(shù),使后續(xù)算法推演中狀態(tài)量可解。

2.2 代價(jià)函數(shù)的引入

引入代價(jià)函數(shù)[10]:

f=w1S+w2L

(3)

式中,S表示當(dāng)前待擴(kuò)展點(diǎn)到目標(biāo)點(diǎn)大圓弧長(zhǎng)度;L表示航跡點(diǎn)與威脅區(qū)最近距離。依次計(jì)算待擴(kuò)展點(diǎn)與威脅區(qū)中心大地線長(zhǎng)度,取最小者。

權(quán)重值取為w2=R/L,w1=1-w2,這樣可以保證在不同的階段,自動(dòng)調(diào)整權(quán)重值。

由于各個(gè)代價(jià)分量的單位各不相同,這樣不便于權(quán)重值的選取。因此,對(duì)上述變量進(jìn)行歸一化處理:

(4)

(5)

式中,m表示待擴(kuò)展的節(jié)點(diǎn)總數(shù)。歸一化處理后代價(jià)函數(shù)可以寫成:

(6)

2.3 路徑生成

設(shè)起點(diǎn)經(jīng)緯度為(L0,B0),目標(biāo)點(diǎn)經(jīng)緯度為(Lf,Bf),當(dāng)前節(jié)點(diǎn)經(jīng)緯度為(L,B)。航跡點(diǎn)與目標(biāo)點(diǎn)距離小于ε,搜索過(guò)程結(jié)束。

(1) 初始判斷

若初始點(diǎn)或目標(biāo)點(diǎn)位于任何一個(gè)威脅區(qū)內(nèi),那么搜索過(guò)程直接終止,無(wú)可行軌跡。

(2) 航跡點(diǎn)生成

具體航跡點(diǎn)生成過(guò)程如下[11-14]:

①當(dāng)前節(jié)點(diǎn)為(L,B),L為經(jīng)度,B為緯度。全方向的8個(gè)航跡點(diǎn)為(L+ΔL,B+ΔB),(L+ΔL,B),(L,B-ΔB),(L+ΔL,B-ΔB),(L,B+ΔB),(L-ΔL,B+ΔB),(L-ΔL,B),(L-ΔL,B-ΔB),將這8點(diǎn)放入OPEN表中,成為待擴(kuò)展航跡點(diǎn)。再添加一個(gè)由當(dāng)前點(diǎn)指向目標(biāo)的大地圓弧上的參考點(diǎn)(L+ΔL,StdB),其中,StdB表示指向目標(biāo)圓弧上點(diǎn)的緯度值。如圖1所示,將這個(gè)點(diǎn)加入OPEN表中。

圖1 航跡點(diǎn)擴(kuò)展示意圖Fig.1 Diagram of extended flight points

②利用9個(gè)航跡點(diǎn)對(duì)CLOSE表進(jìn)行擴(kuò)展。待擴(kuò)展點(diǎn)位于威脅區(qū)域內(nèi)時(shí),刪除該點(diǎn)。

③可達(dá)性判定:判定OPEN表中的6個(gè)航跡點(diǎn)是否位于當(dāng)前點(diǎn)可達(dá)范圍內(nèi)。判斷方法推導(dǎo)如下。

取狀態(tài)向量為

x=(v,θT,σT,φ,λ,r)T

取定哈米爾頓函數(shù)為

H=λ1f1+λ2f2+λ3f3+λ4f4+λ5f5+λ6f6

(7)

λ=[λ1,λ2,λ3,λ4,λ5,λ6]為共軛向量且有

(8)

式中，?fi/?xj(i,j=1,2,…,6)是右函數(shù)偏導(dǎo)數(shù)。

飛行器縱向航程主要由攻角指令決定,橫向最大航程主要由傾側(cè)角γv決定,γv受到再入走廊的約束。需要找到傾斜角和高度之間的關(guān)系,從而將再入走廊的限制引入到橫程最大彈道的計(jì)算過(guò)程中。由質(zhì)心動(dòng)力學(xué)方程中航跡傾角的式子,并將gr=-(fM/r2)[1+J(ae/r)2(1-5sin2φ)]代入,其中ae為地球長(zhǎng)半軸,fM=μ為地球引力系數(shù),J為常值,可得

(9)

將式(9)簡(jiǎn)記(用Matlab程序方便求解)為

(10)

根據(jù)式(1)的第6式可得

(11)

將式(11)代入式(9)即可得到傾斜角γvopt和矢徑r之間的關(guān)系。下面對(duì)地心矢徑和飛行器高度之間的關(guān)系進(jìn)行求解。這里結(jié)合大地橢球的相關(guān)知識(shí),直接給出結(jié)果[15]：

(12)

式中，re為星下點(diǎn)地心矩,計(jì)算式為

(13)

式中,ae為赤道半徑;f為地球扁率;e為地球第一偏心率。

式(12)中φ為地心緯度,φ為地理緯度:

φ=arctan(tanφ/(1-e2))

(14)

將式(13)和式(14)代入式(12)即可得到地心矢和高度之間的關(guān)系,結(jié)合式(9)、式(11)、式(12)推導(dǎo)出傾斜角γvopt和高度之間的關(guān)系,由此將再入走廊的限制引入橫程解算中。

對(duì)于橫向最大位置彈道,飛行時(shí)間tf不作要求,狀態(tài)變量終端約束除要求r(tf)=rf外,其余5個(gè)狀態(tài)變量vf,θTf,σTf,φf(shuō),λf均無(wú)要求,僅要求落點(diǎn)的橫向位移最大[7]。

(15)

從而

λ1(tf)=0,λ2(tf)=0,λ3(tf)=0,

λ4(tf)=k1,λ5(tf)=k2,λ6(tf)=v1

(16)

通過(guò)上述推導(dǎo)過(guò)程,只要確定當(dāng)前航跡點(diǎn)的6個(gè)狀態(tài)量,就可以用橫向最大位置的推導(dǎo)方法并結(jié)合最大升阻比攻角推導(dǎo)最大縱程，從而得到飛行器在當(dāng)前節(jié)點(diǎn)的可達(dá)區(qū)域。

當(dāng)前狀態(tài)量的選取方法如下:

在進(jìn)行航跡規(guī)劃時(shí)認(rèn)為飛行器的規(guī)劃速度隨著飛行器的航程線性變化,則速度為

V=V0-(V0-Vf)·Sgo/(Sgo+Stogo)

(17)

式中,V0為再入初始時(shí)刻飛行器的速度;Vf為再入段結(jié)束時(shí)飛行器的速度;Sgo,Stogo分別為已飛航程、待飛航程。

飛行器在再入飛行過(guò)程當(dāng)前的彈道傾角取為θ<ζ,ζ為設(shè)定的一個(gè)小角度。

航跡點(diǎn)的航跡偏角取為航跡點(diǎn)與初始點(diǎn)連接的圓弧方向與初始航向之間的夾角,即

σ=|Adpointn-Ad0|sign(η)

(18)

式中,Adpointn為航跡點(diǎn)與初始點(diǎn)圓弧的方向;Ad0為初始方向角;sign(η)表示當(dāng)前的偏角符號(hào)由與航跡點(diǎn)位置相關(guān)的變量η決定。

當(dāng)前航跡點(diǎn)的經(jīng)緯度及地心矢都已知。在此基礎(chǔ)上就可以得到可達(dá)區(qū)域,刪除OPEN表中位于可達(dá)區(qū)域之外的點(diǎn)。同時(shí)對(duì)OPEN表中的待擴(kuò)展航跡點(diǎn)也進(jìn)行可達(dá)區(qū)域解算,若目標(biāo)點(diǎn)位于待擴(kuò)展航跡點(diǎn)所形成的可達(dá)區(qū)域之外，則刪除OPEN表中的這個(gè)點(diǎn)。

對(duì)剩下的航跡點(diǎn)選取代價(jià)函數(shù)較小的點(diǎn),將該點(diǎn)添加至CLOSE表中。

(3) 到達(dá)條件判斷

Stogo<ε

(19)

式中,ε是一個(gè)大于一部搜索距離ΔS的常量。當(dāng)獲取的航跡點(diǎn)與目標(biāo)點(diǎn)之間的距離小于ε時(shí)搜索過(guò)程結(jié)束。

(4)收斂性證明

航跡的收斂性是指所規(guī)劃航跡對(duì)于任意初始點(diǎn)出發(fā)最終都能到達(dá)目標(biāo)點(diǎn)。待飛航程表示為

Stogo=Stogolast+ξΔS

(20)

式中,ξ為符號(hào)函數(shù)；Stogolast為上一步的待飛航程。由于指向目標(biāo)的航跡點(diǎn)的存在,能夠保證ξ≤0。又因?yàn)棣為變化量,且ΔS<ε,所以只要經(jīng)過(guò)若干步迭代之后就一定能夠使航跡點(diǎn)滿足結(jié)束條件,到達(dá)目標(biāo)點(diǎn)。

2.4 航路光順化

為了使規(guī)劃航跡具有更強(qiáng)的可行性,需要對(duì)航跡進(jìn)行光順化處理,處理過(guò)程如下。

由以上規(guī)劃方法形成的軌跡由一系列航跡點(diǎn)構(gòu)成,可以表示成如下形式:

ROAD=(Point0,Point1,Point2,…,Pointn)

(21)

從這些航跡點(diǎn)中任意選取兩個(gè)航跡點(diǎn),判斷連線是否經(jīng)過(guò)威脅區(qū),若不經(jīng)過(guò)則把兩個(gè)航跡點(diǎn)之間的航跡點(diǎn)刪除。從初始節(jié)點(diǎn)開始依次進(jìn)行，如圖2所示。

圖2 軌跡規(guī)劃流程圖Fig.2 Flowchart of trajectory planning

可進(jìn)行如下判斷:

上述全部規(guī)劃過(guò)程可用圖2表示。

圖3 橢球面三角形示意圖Fig.3 Diagram of elliptical triangle

3 仿真結(jié)果與分析

進(jìn)行仿真研究,氣動(dòng)參數(shù)通過(guò)氣動(dòng)數(shù)據(jù)表插值獲得。初始狀態(tài)和目標(biāo)點(diǎn)的狀態(tài)如表1和表2所示。

表1 軌跡規(guī)劃狀態(tài)量

仿真過(guò)程中取4個(gè)半徑為500 km的圓域建立威脅區(qū)的模型。

表2 威脅區(qū)的信息

分別采用傳統(tǒng)的A*算法和本文中的A*算法進(jìn)行比較,仿真結(jié)果如圖4所示。

圖4 與傳統(tǒng)算法對(duì)比的軌跡規(guī)劃示意圖Fig.4 Trajectory planning diagram comparing with traditional algorithm

從圖4中可以看出傳統(tǒng)的A*算法由于缺少光順化的過(guò)程,規(guī)劃出的軌跡彎折嚴(yán)重,不適用于實(shí)際的飛行狀況，本文的算法則能夠使軌跡更加平滑。與此同時(shí),通過(guò)記錄仿真時(shí)間發(fā)現(xiàn),由于本文算法中指向目標(biāo)的航跡點(diǎn)的引入,算法可以在1 s內(nèi)完成規(guī)劃任務(wù),而傳統(tǒng)的A*算法花費(fèi)的時(shí)間在2～3 s之間。

為驗(yàn)證算法的適應(yīng)性,普通的計(jì)算機(jī)上,在起點(diǎn)不變的情況下取多個(gè)目標(biāo)點(diǎn)進(jìn)行仿真分析,同時(shí)記錄算法的收斂時(shí)間,如表3和圖5所示。

表3 目標(biāo)點(diǎn)的信息

圖5 多目標(biāo)情況下的軌跡規(guī)劃示意圖Fig.5 Trajectory planning diagram with the situation of multiple targets

通過(guò)仿真結(jié)果可以看出算法的收斂時(shí)間一般在0.6 s左右,針對(duì)不同的目標(biāo)點(diǎn)都能夠準(zhǔn)確快速地規(guī)劃出合理的軌跡。

為進(jìn)一步驗(yàn)證算法的可靠性,對(duì)威脅區(qū)的位置進(jìn)行了調(diào)整,采用初次仿真的初始點(diǎn)和目標(biāo)點(diǎn)坐標(biāo)進(jìn)行驗(yàn)證,仿真結(jié)果如圖6所示。

圖6 變威脅區(qū)情況下的軌跡規(guī)劃示意圖Fig.6 Trajectory planning diagram with the situation of variable threatened zones

從圖6中可以看出本算法在威脅區(qū)改變時(shí)仍然具有較好的適應(yīng)性。在上述仿真的基礎(chǔ)上,針對(duì)本文的方法又進(jìn)行了蒙特卡羅打靶的大樣本仿真,軌跡規(guī)劃成功率接近100%,表明本文的方法具有很強(qiáng)的實(shí)用性。

4 結(jié) 論

本文針對(duì)滑翔飛行器這一研究對(duì)象,對(duì)傳統(tǒng)的A*算法進(jìn)行了改進(jìn)。成功地將滑翔飛行器飛行過(guò)程中的各種過(guò)程約束引入到算法的推演過(guò)程中。同時(shí)結(jié)合實(shí)際應(yīng)用特點(diǎn)進(jìn)行橢球面上相關(guān)量的推演。通過(guò)仿真分析,證明了本文算法收斂速度較快,針對(duì)不同目標(biāo),不同威脅區(qū)都具有較好的適應(yīng)性。

[1] Jorris T R. Common aero vehicle automous reentry trajectory optimization satisfying waypoint and no-fly zone constrains[D]. Ohio: Air Force Instituite of Technology, 2007.

[2] Jorris T R, Cobb R G. Multiple method 2-D trajectory optimization satifying waypoints and no-fly zone contraints[J].JournalofGuidance,ControlandDynamic, 2008, 31(3):543-553.

[3] Jorris T R, Cobb R G. Three-dimensional trajectory optimization satisfying waypoint and no-fly zone contraints[J].JournalofGuidance,ControlandDynamic, 2009, 32(2):551-572.

[4] Howlett J, Goodrich M, Mclain T. Learning real-time A*path planner for sensing closely-spaced targets from an aircraft[C]∥Proc.oftheAIAAGuidance,NavigationandControlConferenceandExihibit, 2003:AIAA2003-5338.

[5] Richards N D, Sharma M, Ward D G. A hybrid A*/automaton approach to on-line path planning with obstacle avoidance[C]∥Proc.oftheAIAA1stIntelligentSystemsTechnicalConference, 2004: AIAA2004-6229.

[6] Huang X, Huang P F, Yan J, et al. Route planning for hypersonic vehicle using A*algorithm[J].ComputerSimulation, 2009, 26(9): 62-127. (黃雄, 黃攀峰, 閆杰, 等. 基于A*算法的高超聲速飛行器航跡規(guī)劃方法[J].計(jì)算機(jī)仿真, 2009, 26(9): 62-65，127.)

[7] Zhao H Y.Reentryguidanceanddynamicsofaircraft[M]. Changsha: National University of Defense Technology Press, 1997. (趙漢元. 飛行器再入動(dòng)力學(xué)和制導(dǎo)[M]. 長(zhǎng)沙: 國(guó)防科技大學(xué)出版社, 1997.)

[8] Barton J D. Fundamentals of small unmanned aircraft flight[J].JohnsHopkinsAppliedPhysicsLaboratoryTechnicalDigest, 2012,31(2): 132-149.

[9] Grant M J, Clark I G, Braun R D. Rapid entry corridor trajectory optimization for conceptual design[C]∥Proc.oftheAIAAAtmosphericFlightMechanicsConference, 2010:AIAA2010-7810.

[10] Li M, Wang D B, Sheng S Z, et al. Multiple route planning based on particle swarm optimization and weighted K-means clustering[J].SystemsEngineeringandElectronics, 2012, 34(3): 512-516. (李猛, 王道波, 盛守照, 等. 基于加權(quán)K-均值聚類與粒子群優(yōu)化的多航跡規(guī)劃[J].系統(tǒng)工程與電子技術(shù), 2012, 34(3): 512-516.)

[11] Li X Y, Zhou D Y, Feng Q. Multiple routes planning for A*algorithm based on hierarchical planning[J].SystemsEngineeringandElectronics, 2015, 37(2): 318-322. (李梟揚(yáng), 周德云, 馮琦. 基于分級(jí)規(guī)劃策略的A*算法多航跡規(guī)劃[J]. 系統(tǒng)工程與電子技術(shù), 2015, 37(2): 318-322.)

[12] Pehlivanoglu Y V. A new vibrational genetic algorithm enhanced with a voronoi diagram for path planning of autonomous UAV[J].AerospaceScienceandTechnology, 2012, 16(1): 47-55.

[13] Foo J L, Knutzon J, Kalivarapu V, et al. Path planning of unmanned aerial vehicles using B-splines and particle swarm optimization[J].JournalofAerospaceComputing,Information,andCommunication, 2009, 6(4): 271-290.

[14] Poli R. Analysis of the publications on the applications of particle swarm optimization[J].JournalofArtificialEvolutionandApplications, 2008, 2008(1): 1-10.

[15] Lv Z P, Qiao S B.Foundationofgeodesy[M]. Wuhan: Wuhan University Press, 2010. (呂志平, 喬書波. 大地測(cè)量學(xué)基礎(chǔ)[M]. 武漢: 武漢大學(xué)出版社, 2010.)

盧 青(1990-), 男, 博士研究生, 主要研究方向?yàn)轱w行器航跡規(guī)劃及制導(dǎo)控制理論。

E-mail: lamaxiya1990@163.com

周 軍(1966-), 男, 教授, 博士, 主要研究方向?yàn)轱w行器制導(dǎo)控制與先進(jìn)控制理論。

E-mail: zhoujun@nwpu.edu.cn

呼衛(wèi)軍(1979-), 男, 副教授, 博士, 主要研究方向?yàn)轱w行器數(shù)字仿真技術(shù)。

E-mail: huyanwj@126.com

Trajectory planning for gliding aircraft based on improved A*algorithm

LU Qing, ZHOU Jun, HU Wei-jun

(InstituteofPrecisionGuidanceandControl,NorthwesternPolytechnicalUniversity,Xi’an710072,China)

In order to solve the lateral maneuver problem of gliding aircraft during reentry process, the theory of A*algorithm is introduced into trajectory planning. Earth curvature must be considered when flying distance is far, besides, there are varies constrains during reentry process. According to the boundedness of existing trajectory planning using A*algorithm based on plane, an A*algorithm for trajectory planning based on ellipsoid is put forward. Process constrains are also introduced into the searching procedure of A*algorithm through deduction. Finally, a reasonable trajectory can be found. The results of simulation verify that the method of this paper can find out a trajectory rapidly and efficiently, and this method has strong commodity.

A*Algorithm; gliding aircraft; trajectory planning; ellipsoid

2016-02-22;

2016-10-18;網(wǎng)絡(luò)優(yōu)先出版日期:2016-10-27。

國(guó)家自然科學(xué)基金(61473226)資助課題

TP 273

A

10.3969/j.issn.1001-506X.2016.12.12

網(wǎng)絡(luò)優(yōu)先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20161027.1612.024.html