本發(fā)明涉及一種彈道規(guī)劃
技術(shù)領(lǐng)域：
，尤其涉及一種通過曲線外形計算彈道參數(shù)的草圖交互虛擬域動態(tài)逆彈道規(guī)劃方法。
背景技術(shù)：
：隨著導(dǎo)彈防御系統(tǒng)全面部署的臨近，各軍事強(qiáng)國目前競相研究具有一定機(jī)動能力的導(dǎo)彈以提高導(dǎo)彈突防能力從而保持導(dǎo)彈打擊能力。然而這類導(dǎo)彈多在大氣層內(nèi)飛行，在滿足發(fā)射約束與任務(wù)需求的同時，必須考慮能量、過載、熱流、動壓等約束，從而大大增加了導(dǎo)彈飛行過程復(fù)雜度。如何對飛行彈道進(jìn)行規(guī)劃，在滿足復(fù)雜約束條件下提高作戰(zhàn)效能，是此類導(dǎo)彈作戰(zhàn)應(yīng)用面臨的重大挑戰(zhàn)。常規(guī)彈道規(guī)劃方法應(yīng)用于此類導(dǎo)彈彈道規(guī)劃時存在如下困難：(1)以控制量為彈道規(guī)劃設(shè)計變量，通過調(diào)整控制量改變導(dǎo)彈受力狀態(tài)，從而實現(xiàn)彈道修改，無法根據(jù)任務(wù)需求直接對彈道位置特征進(jìn)行調(diào)整，導(dǎo)致彈道規(guī)劃過程過于復(fù)雜；(2)在導(dǎo)彈受力狀態(tài)改變之后，一般需要通過彈道積分獲得導(dǎo)彈的飛行特性與約束條件的滿足情況，由于彈道規(guī)劃中積分運(yùn)算需要大量重復(fù)迭代，導(dǎo)致彈道規(guī)劃過程效率較低；(3)通過改變控制量間接改變彈道的方式，在實際應(yīng)用中難以獲得與目標(biāo)彈道完全一致的規(guī)劃結(jié)果，導(dǎo)致彈道規(guī)劃過程無法全面反映設(shè)計者意愿。LuPing[1]將動態(tài)逆思想應(yīng)用于制導(dǎo)領(lǐng)域，使制導(dǎo)算法更加簡潔，提高了航天飛機(jī)上升段制導(dǎo)效率。動態(tài)逆方法由于能夠從路徑得到控制量，因而在機(jī)器人路徑規(guī)劃、戰(zhàn)斗機(jī)和無人機(jī)航跡規(guī)劃及巡航導(dǎo)彈彈道規(guī)劃等領(lǐng)域受到關(guān)注，并取得初步的階段成果[2]-[4]。由于動態(tài)逆思想需要將軌跡方程不斷微分以得到速度、加速度等必要信息，所以要求系統(tǒng)方程各階導(dǎo)數(shù)在時域中有特定含義，而這樣做又難以滿足“軌跡可調(diào)節(jié)”的意愿。對此，Yakimenko[5],[6]通過引入虛擬路徑，將軌跡優(yōu)化問題從時域轉(zhuǎn)化到虛擬域，再用動態(tài)逆的方法優(yōu)化軌跡，提出了虛擬域動態(tài)逆(InverseDynamicsintheVirtualDomain，IDVD)方法。該方法通過虛擬速度將虛擬域路徑和時域路徑對應(yīng)起來。虛擬域動態(tài)逆方法首先將虛擬域參數(shù)微分，得到虛擬域速度、加速度和過載等信息，但這些信息沒有時域的時間標(biāo)簽。接下來利用虛擬速度將虛擬域加速度和過載等信息對應(yīng)到時域上，并加上時域的時間標(biāo)簽。這樣做可以得到解析的參數(shù)求解方程，因此計算速度非常快。由于虛擬域軌跡的設(shè)計是通過多項式擬合的，虛擬軌跡并不能保證全局最優(yōu)，進(jìn)而也無法保證時域軌跡的全局最優(yōu)。JohnA.Lukacs[7],IanD.Cowling[8],GeorgeBoyarko等[9]分別與Yakimenko合作，研究了虛擬域動態(tài)逆方法在四旋翼無人機(jī)軌跡規(guī)劃和制導(dǎo)、攔截彈在線制導(dǎo)、飛行器交會等的應(yīng)用。NAkhtar等[10]研究了虛擬域動態(tài)逆方法在滑翔機(jī)實時軌跡規(guī)劃仿真中的應(yīng)用，結(jié)果表明該方法計算速度快，軌跡規(guī)劃方案有利于節(jié)省燃料。閆梁等[3]使用變節(jié)點算法對虛擬域動態(tài)逆方法進(jìn)行了改進(jìn)，在節(jié)點總數(shù)不變的情況下提高了運(yùn)算效率，但算法在大部分情況下，可能需經(jīng)多次迭代積分確定節(jié)點位置，運(yùn)算復(fù)雜度高。ZhangYu等[2]將虛擬域動態(tài)逆方法應(yīng)用在無人戰(zhàn)斗機(jī)的軌跡規(guī)劃仿真中，獲得了較高的計算速度和較好的收斂特性。MarcoCiarcià等[11]基于虛擬域動態(tài)逆方法提出了一種近憂制導(dǎo)策略，并應(yīng)用于航天器交會對接地面試驗中，通過對比實驗證明了方法的有效性。JacopoVentura等[12]研究了快速最優(yōu)姿態(tài)軌跡生成問題，評估了不同情況下虛擬域動態(tài)逆方法的計算性能。技術(shù)實現(xiàn)要素：傳統(tǒng)彈道規(guī)劃方法應(yīng)用于具有變向能力的導(dǎo)彈，例如助推滑翔導(dǎo)彈、有機(jī)動能力的彈道導(dǎo)彈和巡航導(dǎo)彈等存在的不適應(yīng)性。如現(xiàn)有虛擬域動態(tài)逆方法在最優(yōu)控制理論框架下提出，其目的是為了解決制導(dǎo)問題，采用了處理兩點邊值問題的方法解決軌跡規(guī)劃問題，一般能夠?qū)ζ瘘c和終點的位置參數(shù)進(jìn)行調(diào)節(jié)，但對過程參數(shù)如飛行軌跡不能進(jìn)行人機(jī)交互調(diào)整，因此不具備對彈道的交互調(diào)控能力。針對現(xiàn)有技術(shù)存在的以上缺陷，本發(fā)明提出了一種彈道規(guī)劃方法。該方法引入人的智能因素，基于解析曲線直接設(shè)計彈道特征，逆向計算控制量與約束條件，并通過草圖交互技術(shù)將設(shè)計人員意愿與計算機(jī)規(guī)劃過程有機(jī)結(jié)合，顯著改善彈道規(guī)劃過程的直觀性，提高復(fù)雜約束下助推滑翔導(dǎo)彈彈道規(guī)劃的效率。本發(fā)明采用的技術(shù)方案是：一種彈道規(guī)劃方法，包括以下步驟：S1.建立助推滑翔導(dǎo)彈運(yùn)動模型；設(shè)地平面為一平面，導(dǎo)彈飛行側(cè)滑角為零，得到助推滑翔導(dǎo)彈三自由度動力學(xué)計算方程：x·1=Vcosγ·cosψx·2=Vcosγ·sinψx·3=-V·sinψV·=g(nx-sinγ)γ·=gV-1(ny-cosγ)ψ·=gnzV-1·cos-1γ---(1)]]>其中，x1、x2、x3分別表示助推滑翔導(dǎo)彈在北東地坐標(biāo)系下沿三個坐標(biāo)軸的位置，V為導(dǎo)彈飛行速度，γ為彈道傾角，ψ為彈道偏角，g為重力加速度；nx、ny、nz分別為彈道坐標(biāo)系下沿坐標(biāo)軸三個方向的過載：nx=Tcosα-Dmgny=Tsinα+Lmgnz=Zmg---(2)]]>其中，T為導(dǎo)彈發(fā)動機(jī)推力，α為導(dǎo)彈飛行攻角，D為導(dǎo)彈飛行阻力，L為導(dǎo)彈飛行升力，Z為導(dǎo)彈飛行側(cè)向力，m為助推滑翔導(dǎo)彈的質(zhì)量；T、α和m均為時間的函數(shù)，通過插值或擬合后的公式求得：T＝T(t)(3)α＝α(t)(4)m＝m(t)(5)速度系下導(dǎo)彈的阻力、升力、側(cè)向力為D=12ρV2SCD---(6)]]>L=12ρV2SCL---(7)]]>Z=12ρV2SCZ---(8)]]>其中ρ為大氣密度，S為特征面積，CD、CL、CZ分別為阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)，通過插值或擬合后的公式得到：CD＝CD(α,V,-x3)(9)CL＝CL(α,V,-x3)(10)CZ＝CZ(α,V,-x3)(11)S2.約束條件建模(1)駐點熱流密度駐點熱流密度qws采用Kemp-Riddell公式計算：qws=131884.2RN(ρ∞ρo)1/2(v∞vc)3.25(1-hw/hs)≤qwsmax---(12)]]>式中ρο＝1.225(kg/m3)；νc＝7900(m/s)；RN－駐點曲率半徑(m)；qws－駐點熱流密度(kW/m2)；hw－壁面焓值(kJ/kg)；hs－駐點焓值(kJ/kg)；qwsmax－最大駐點熱流密度(kW/m2)；ρ∞－來流密度；ν∞－來流速度；(2)過載nnx≤nxmaxny≤nymaxnz≤nzmax---(13)]]>其中nxmax、nymax、nzmax為助推滑翔導(dǎo)彈在彈道坐標(biāo)系下三個方向可承受的最大過載；(3)動壓qq=12ρV2≤qmax---(14)]]>qmax為助推滑翔導(dǎo)彈可承受的最大動壓；(4)控制量當(dāng)僅考慮助推滑翔導(dǎo)彈縱平面運(yùn)動時，取導(dǎo)彈飛行攻角α為控制量，需要對其幅值進(jìn)行約束，如下：αmin≤α≤αmax(15)αmin為最小攻角，αmax為最大攻角；(5)起始條件和終端約束助推滑翔導(dǎo)彈的初始位置x1=x10,x2=x20,x3=x30---(16)]]>和為給定的發(fā)射點坐標(biāo)；S3.虛擬域草圖快速生成方法；S3.1助推滑翔導(dǎo)彈垂直上升段草圖生成方法取發(fā)射點坐標(biāo)為原點，即O(0,0,0)，那么可通過控制點V0(0,0,V0z)來定義垂直上升段軌跡，其中V0z為控制點V0的z軸坐標(biāo)，O和V0即為垂直上升段彈道的控制點；不失一般性，設(shè)τ∈[a,a+1)，其中a＝0,1,2,…，τ為虛擬域的自變量即虛擬域時間，那么在北東地坐標(biāo)系Oxyz下，虛擬域軌跡為：x(τ)=0y(τ)=0z(τ)=V0z(τ-a)2---(17)]]>其中x(τ)、y(τ)、z(τ)分別為北東地坐標(biāo)系Oxyz中x軸、y軸和z軸的坐標(biāo)；然后得到x(τ)、y(τ)、z(τ)在虛擬域的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)：x′(τ)=0y′(τ)=0z′(τ)=2V0z(τ-a)x′′(τ)=0x′′(τ)=0z′′(τ)=2V0z---(18)]]>S3.2助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段草圖生成方法助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段的草圖可通過C1連續(xù)曲線、C2連續(xù)曲線或者C3連續(xù)曲線生成；S4.草圖交互虛擬域動態(tài)逆彈道參數(shù)求解方法；通過S3獲得一條完整曲線后，可將曲線表示成為只含有一個自變量的函數(shù)r(τ)，以虛擬域時刻τ為自變量，τ的取值范圍是[0,n]，其中n為曲線段數(shù)；假設(shè)有函數(shù)f(·)，τ和t分別為虛擬域時刻和時域時刻，τ和t有數(shù)學(xué)關(guān)系：df(τ)dt=df(τ)dτdτdt---(34)]]>在此定義λ(τ)=dτdt---(35)]]>λ(τ)為虛擬速度，并約定表示f(·)在時域內(nèi)求導(dǎo)，即f′(·)表示f(·)在虛擬域求導(dǎo)，即根據(jù)(34)、(35)兩式有：f·(τ)=f′(τ)λ(τ)---(36)]]>進(jìn)而有：f··(τ)=d(f′(τ)λ(τ))dτdτdt=f′′(τ)λ2(τ)+f′(τ)λ(τ)λ′(τ)---(37)]]>確定虛擬域曲線r(τ)后，在τ的取值范圍取N個節(jié)點將虛擬域曲線r(τ)離散；若τ的取值范圍為[0,τf]，則可通過如下方法獲得節(jié)點之間在虛擬域的時間間隔Δτ=τfN-1---(38)]]>每一個節(jié)點對應(yīng)的虛擬域時刻可表示為τj＝τj-1+Δτ,j＝2,…,N(39)其中N為偶數(shù)；助推滑翔導(dǎo)彈采用固體發(fā)動機(jī)作為動力，時間-質(zhì)量關(guān)系m(t)、時間-推力關(guān)系T(t)視為已知；首先求出當(dāng)前節(jié)點處虛擬域基本狀態(tài)xi、x′i、x″i(i＝1,2,3)，由r(τ)以及其一階、二階導(dǎo)數(shù)得到虛擬域基本狀態(tài)，整理后結(jié)果形式如下：x1,j-1x1,j-1′x1,j-1′′x2,j-1x2,j-1′x2,j-1′′x3,j-1x3,j-1′x3,j-1′′=r(τj-1)r′(τj-1)r′′(τj-1)---(40)]]>當(dāng)τj-1∈[0,1)時r(τj-1)＝r1(τj-1)，當(dāng)τj-1∈[1,2]時r(τj-1)＝r2(τj-1)，其中xi,j-1表示虛擬域基本狀態(tài)xi在節(jié)點(j-1)處的值，即xi(τj-1)；同樣的x′i,(j-1)表示虛擬域基本狀態(tài)x′i在節(jié)點(j-1)處的值，x″i,(j-1)表示虛擬域基本狀態(tài)x″i在節(jié)點(j-1)處的值；其中各節(jié)點處的狀態(tài)量和控制量的計算方法如下：(1)彈道傾角γj-1=arctan-x3,j-1′x1,j-1′2+x2,j-1′2---(41)]]>(2)彈道偏角ψj-1=arctanx2,j-1′x1,j-1′---(42)]]>(3)彈道傾角虛擬域?qū)?shù)γ′j-1和彈道偏角虛擬域?qū)?shù)ψ′j-1γj-1′=x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32cos2γj-1---(43)]]>ψj-1′=x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1---(44)]]>(4)y方向過載ny,j-1=Vj-1g-1λj-1x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32·cos2γj-1+cosγj-1---(45)]]>(5)z方向過載nz,j-1=Vj-1g-1λj-1x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1cosγj-1---(46)]]>(6)大氣密度hj-1＝-x3,j-1(47)ρj-1＝fair(hj-1)(48)其中fair(hj-1)為大氣密度函數(shù)；(7)升力采用攻角的一次函數(shù)表示升力系數(shù)即：CL,j-1＝l1αj-1+l2(49)那么升力Lj-1=12ρj-1Vj-12S(l1αj-1+l2)---(50)]]>其中S為導(dǎo)彈特征面積，假設(shè)這一數(shù)值在飛行過程中不變；其中l(wèi)1和l2為擬合系數(shù)；(8)攻角αj-1=ny,j-1mj-1g-12ρj-1Vj-12Sl2Tj-1+12ρj-1Vj-12Sl1---(51)]]>(9)阻力nx,j-1=Tj-1cosαj-1-12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)mj-1g---(52)]]>其中CD,j-1=d1Vj-12+d2Vj-1+Vj-1---(53)]]>Dj-1=12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)---(54)]]>d1和d2為擬合系數(shù)；(8)速度的虛擬域?qū)?shù)V′(τ)=V·(τ)/λ(τ)---(55)]]>(9)下一節(jié)點處速度Vj＝Vj-1+g(nx,j-1-sinγj-1)/λj-1Δτ(56)(10)節(jié)點(j-1)到節(jié)點j的飛行時間Δtj=Σi=13(xi,j-xi,j-1)2(Vj+Vj-1)/2---(57)]]>(11)下一節(jié)點時域時刻tj＝tj-1+Δtj(58)(12)Δt2的估計值Δt~2=TboostN-1---(59)]]>(13)λ1的估計值λ~1=ΔτΔt2≈ΔτΔt~2=τf(N-1)(N-1)Tboost=τfTboost.---(60)]]>本發(fā)明S3.2中，助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段的草圖通過C1連續(xù)曲線生成的方法如下：首先定義下列3個關(guān)于t的多項式為帶形狀參數(shù)λb和μb的三次調(diào)配函數(shù)：X0(t)=1-λbt+(2λb-3)t2+(2-λb)t3X1(t)=λbt+(μb-2λb)t2+(λb-μb)t3X2(t)=(3-μb)t2+(μb-2)t3,t∈[0,1]---(19)]]>其中0≤λb,μb≤3；然后給定二維或三維空間中3個控制頂點Vi(i＝0,1,2)，稱曲線r(t)為可調(diào)控的三次參數(shù)曲線：r(t)＝X0(t)V0+X1(t)V1+X2(t)V2,t∈[0,1](20)其中Xi(t)(i＝0,1,2)為按式(19)定義的調(diào)配函數(shù)；或者，將上式(20)寫成矩陣形式為：r(t)=1tt2t3100-λbλb02λb-3μb-2λb3-λb2-λbλb-μbμb-2V0V1V2,t∈[0,1].---(21)]]>本發(fā)明S3.2中，助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段的草圖通過C2連續(xù)曲線生成的方法如下：給定一組控制點<V0,V1,…,Vn＞，以其構(gòu)成控制切線多邊形，取重構(gòu)控制點：b2i=Vib2i+1=λciVi+(1-λci)Vi+1,i=0,1,...,n-1---(22)]]>式(22)中λci為切點控制參數(shù)，且0＜λci＜1；那么由重構(gòu)控制點b0,b1,…,b2n構(gòu)成的曲線為：ri(t)=Σj=03Bj(t)bi+j,i=0,1,...,2n-2,0≤t≤αc---(23)]]>其中αc為調(diào)節(jié)參數(shù)，且0≤αc≤π；B0(t),B1(t),B2(t),B3(t)為C-B樣條基函數(shù)，分別定義如下：B0(t)=(αc-t)-sin(αc-t)2αc(1-C)B3(t)=t-sin(t)2αc(1-C)B1(t)=B3(t)-2B0(t)+2(αc-t)(1-C)2αc(1-C)B2(t)=B0(t)-2B3(t)+2t(1-C)2αc(1-C)---(24)]]>其中C＝cosαc，S＝sinαc；或者，將C2連續(xù)的C-B樣條曲線寫成矩陣形式為：ri(t)=Σj=03Bj(t)bi+j=12αc(1-C)sintcostt1]]>C-(1+2C)2+C-1-S2S-S0-11+2C-(1+2C)1αc-2αcCαc0bibi+1bi+2bi+3---(25)]]>其中i＝0,1,…,2n。本發(fā)明S3.2中，助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段的草圖通過C3連續(xù)曲線生成的方法如下：給定一組控制點<V0,V1,…,Vn>，以其構(gòu)成控制多邊形，構(gòu)造邊矢量：ai＝Vi-Vi-1,i＝1,2,…,n(27)則可計算控制多邊形每個頂點處的切矢量T0=t0(-a2)+(1-t0)a1Ti=tiai+(1-ti)ai+1Tn=tnan+(1-tn)(-an-1),i=1,2,...,n-1---(28)]]>令Vi*=Vi+|ai+1×Ti+1||Ti×Ti+1|Ti,i=1,2,...,n-1---(29)]]>及l(fā)0=|V0*-V0|li=min{|Vi*-Vi|,|Vi-Vi-1*|}ln=|Vn-Vn-1*|,i=1,2,...,n-1---(30)]]>取重構(gòu)控制點d4i=Vi-λiliTi|Ti|d4i+1=Vi-12λiliTi|Ti|d4i+2=Vi+12λiliTi|Ti|d4i+3=Vi+λiliTi|Ti|,i=0,1,...,n-1---(31)]]>其中λi為控制調(diào)節(jié)參數(shù)，且0＜λi＜1；那么由重構(gòu)控制點可構(gòu)成四次B樣條曲線ri(t)=Σj=04Ej(t)di+j,i=0,1,...,4n-1,0≤t≤1---(32)]]>其中E0(t)=(t4-4t3+6t2-4t+1)/24E1(t)=(-4t4+12t3-6t2-12t+11)/24E2(t)=(6t4-12t3-6t2+12t+11)/24E3(t)=(-4t4+4t3+6t2+4t+1)/24E4(t)=t4/24.]]>本發(fā)明的有益效果在于：(1)本發(fā)明脫離最優(yōu)控制理論框架，基于完整的解析曲線直接設(shè)計彈道特征，逆向計算控制量與約束條件，不僅能夠?qū)ζ瘘c和終點的位置參數(shù)進(jìn)行調(diào)節(jié)，更能對過程參數(shù)如飛行軌跡進(jìn)行調(diào)整。與現(xiàn)有最好技術(shù)相比，本發(fā)明不局限于解決導(dǎo)彈的制導(dǎo)問題，除此之外還可解決導(dǎo)彈的彈道規(guī)劃、彈道設(shè)計、實現(xiàn)彈道能力評估等問題。(2)本發(fā)明將人類智慧通過草圖彈道的形式引入算法中，與現(xiàn)有技術(shù)相比，實現(xiàn)了人與計算機(jī)的高效交互溝通，極大方便了彈道規(guī)劃者利用計算機(jī)實現(xiàn)本人意愿，降低了因計算機(jī)理解偏差導(dǎo)致復(fù)算的可能性，提高了計算效率。(3)現(xiàn)有技術(shù)是通過調(diào)整導(dǎo)彈飛行的控制量實現(xiàn)對彈道的調(diào)整，調(diào)整過程復(fù)雜、直觀性差，難以進(jìn)行形式復(fù)雜彈道的調(diào)整或規(guī)劃。本發(fā)明實現(xiàn)了規(guī)劃者對彈道的直接調(diào)整，簡化了調(diào)整過程、直觀性強(qiáng)，易于進(jìn)行形式復(fù)雜彈道的調(diào)整或規(guī)劃。(4)本發(fā)明使用解析的曲線不斷微分得到彈道參數(shù)，與目前常用的積分法相比省去了積分環(huán)節(jié)，通過仿真實驗驗證，每條彈道的平均計算時長為0.013s，較積分法提高三個數(shù)量級。附圖說明圖1為本發(fā)明的實現(xiàn)過程示意圖；圖2為本發(fā)明的流程圖；具體實施方式下面結(jié)合附圖和具體實施方式對本發(fā)明作進(jìn)一步說明。一種彈道規(guī)劃方法，包括以下步驟：S1.建立助推滑翔導(dǎo)彈運(yùn)動模型；設(shè)地平面為一平面，導(dǎo)彈飛行側(cè)滑角為零，可得到助推滑翔導(dǎo)彈三自由度動力學(xué)計算方程：x·1=Vcosγ·cosψx·2=Vcosγ·sinψx·3=-V·sinψV·=g(nx-sinγ)γ·=gV-1(ny-cosγ)ψ·=gnzV-1·cos-1γ---(1)]]>其中，x1、x2、x3分別表示助推滑翔導(dǎo)彈在北東地坐標(biāo)系下沿三個坐標(biāo)軸的位置，V為導(dǎo)彈飛行速度，γ為彈道傾角，ψ為彈道偏角，g為重力加速度，文中統(tǒng)一取為9.8m/s2；nx、ny、nz分別為彈道坐標(biāo)系下沿坐標(biāo)軸三個方向的過載：nx=Tcosα-Dmgny=Tsinα+Lmgnz=Zmg---(2)]]>其中，T為導(dǎo)彈發(fā)動機(jī)推力，α為導(dǎo)彈飛行攻角，D為導(dǎo)彈飛行阻力，L為導(dǎo)彈飛行升力，Z為導(dǎo)彈飛行側(cè)向力，m為助推滑翔導(dǎo)彈的質(zhì)量。T、α和m均為時間的函數(shù)，一般通過插值或擬合后的公式求得：T＝T(t)(3)α＝α(t)(4)m＝m(t)(5)速度系(定義見文獻(xiàn)[15])下導(dǎo)彈的阻力、升力、側(cè)向力為D=12ρV2SCD---(6)]]>L=12ρV2SCL---(7)]]>Z=12ρV2SCZ---(8)]]>其中ρ為大氣密度，S為特征面積，CD、CL、CZ分別為阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)，它們與攻角α、飛行速度V和飛行高度h(h＝-x3)有關(guān)，一般在進(jìn)行實驗后制作成插值表，通過插值或擬合后的公式得到：CD＝CD(α,V,-x3)(9)CL＝CL(α,V,-x3)(10)CZ＝CZ(α,V,-x3)(11)S2.約束條件建模(1)駐點熱流密度因為駐點是助推滑翔導(dǎo)彈加熱較為嚴(yán)重的部分，因此通常作為彈道規(guī)劃的重要約束條件。本發(fā)明中采用Kemp-Riddell公式計算：qws=131884.2RN(ρ∞ρo)1/2(v∞vc)3.25(1-hw/hs)≤qwsmax---(12)]]>式中ρο＝1.225(kg/m3)；νc＝7900(m/s)；RN－駐點曲率半徑(m)；qws－駐點熱流密度(kW/m2)；hw－壁面焓值(kJ/kg)；hs－駐點焓值(kJ/kg)；qwsmax－最大駐點熱流密度(kW/m2)；ρ∞－來流密度；ν∞－來流速度；(2)過載n為使導(dǎo)彈結(jié)構(gòu)在飛行過程中不發(fā)生破壞，一般需要對過載進(jìn)行約束。nx≤nxmaxny≤nymaxnz≤nzmax---(13)]]>其中nxmax、nymax、nzmax為助推滑翔導(dǎo)彈在彈道坐標(biāo)系下三個方向可承受的最大過載。(3)動壓q動壓一定程度上反映助推滑翔導(dǎo)彈飛行環(huán)境的優(yōu)劣，一般也作為重要的約束條件給出。q=12ρV2≤qmax---(14)]]>qmax為助推滑翔導(dǎo)彈可承受的最大動壓。(4)控制量當(dāng)僅考慮助推滑翔導(dǎo)彈縱平面運(yùn)動時，一般取導(dǎo)彈飛行攻角α為控制量，這里需要對其幅值進(jìn)行約束。αmin≤α≤αmax(15)αmin為最小攻角，αmax為最大攻角。(5)起始條件和終端約束根據(jù)任務(wù)要求，可以對位置、速度、速度傾角、攻角等進(jìn)行約束。大多數(shù)情況下，彈道規(guī)劃問題中都含有助推滑翔導(dǎo)彈的初始位置這一起始條件。x1=x10,x2=x20,x3=x30---(16)]]>和為給定的發(fā)射點坐標(biāo)。上述(1)至(5)約束條件中，不等式兩端的值，如qwsmax、αmin等均為導(dǎo)彈參數(shù)，由導(dǎo)彈的飛行性能決定，在導(dǎo)彈定型后，這些參數(shù)也隨之確定，可視為已知。S3.虛擬域草圖快速生成方法；(1)助推滑翔導(dǎo)彈垂直上升段草圖生成方法取發(fā)射點坐標(biāo)為原點，即O(0,0,0)，那么可通過控制點V0(0,0,V0z)來定義垂直上升段軌跡，其中V0z為控制點V0的z軸坐標(biāo)，O和V0即為垂直上升段彈道的控制點?？刂泣c為控制一條曲線外形的關(guān)鍵節(jié)點，可以通過改變它們的位置調(diào)整曲線形狀，從而構(gòu)建所需的曲線。不失一般性，設(shè)τ∈[a,a+1)(a＝0,1,2,…)為虛擬域的自變量，即虛擬域時間，那么在北東地坐標(biāo)系Oxyz下，虛擬域軌跡為：x(τ)=0y(τ)=0z(τ)=V0z(τ-a)2---(17)]]>其中x(τ)、y(τ)、z(τ)分別為北東地坐標(biāo)系Oxyz中x軸、y軸和z軸的坐標(biāo)。然后得到x(τ)、y(τ)、z(τ)在虛擬域的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)：x′(τ)=0y′(τ)=0z′(τ)=2V0z(τ-a)x′′(τ)=0x′′(τ)=0z′′(τ)=2V0z---(18)]]>(2)助推滑翔導(dǎo)彈轉(zhuǎn)彎及滑翔段草圖生成方法首先對曲線連續(xù)性作必要說明。設(shè)有m段多項式曲線構(gòu)成一條完整曲線C(τ)，其中第i條多項式曲線為Ci(τ),i＝1,2,…,m.若為Ci(τ)的j階導(dǎo)矢，τi為Ci(τ)的末端點自變量，如果0≤j≤k對每段曲線都成立，則稱曲線C(τ)在τi處是Ck連續(xù)的[23],[24]。轉(zhuǎn)彎及滑翔段的草圖可通過C1連續(xù)曲線、C2連續(xù)曲線或者C3連續(xù)曲線生成，下面分別介紹它們的生成方法。(a)C1連續(xù)曲線草圖生成方法貝塞爾曲線[25]具有操作直觀，連接平滑等特點，適于解決本發(fā)明中的彈道表示問題。這里介紹其中一類可調(diào)控的三次多項式曲線。定義1下列3個關(guān)于t的多項式為帶形狀參數(shù)λb和μb的三次調(diào)配函數(shù)：X0(t)=1-λbt+(2λb-3)t2+(2-λb)t3X1(t)=λbt+(μb-2λb)t2+(λb-μb)t3X2(t)=(3-μb)t2+(μb-2)t3,t∈[0,1]---(19)]]>其中0≤λb,μb≤3。定義2給定二維或三維空間中3個控制頂點Vi(i＝0,1,2)，稱曲線r(t)為可調(diào)控的三次參數(shù)曲線：r(t)＝X0(t)V0+X1(t)V1+X2(t)V2,t∈[0,1](20)其中Xi(t)(i＝0,1,2)為按式(19)定義的調(diào)配函數(shù)。上式(20)改寫成矩陣形式為：r(t)=1tt2t3100-λbλb02λb-3μb-2λb3-λb2-λbλb-μbμb-2V0V1V2,t∈[0,1]---(21)]]>(b)C2連續(xù)曲線草圖生成方法給定一組控制點<V0,V1,…,Vn＞，以其構(gòu)成控制切線多邊形，取重構(gòu)控制點[31][32]b2i=Vib2i+1=λciVi+(1-λci)Vi+1,i=0,1,...,n-1---(22)]]>式(22)中λci為切點控制參數(shù)[33]，且0＜λci＜1。那么由重構(gòu)控制點b0,b1,…,b2n構(gòu)成的曲線為：ri(t)=Σj=03Bj(t)bi+j,i=0,1,...,2n-2,0≤t≤αc---(23)]]>其中αc為調(diào)節(jié)參數(shù)，且0≤αc≤π。B0(t),B1(t),B2(t),B3(t)為C-B樣條基函數(shù)，分別定義如下：B0(t)=(αc-t)-sin(αc-t)2αc(1-C)B3(t)=t-sin(t)2αc(1-C)B1(t)=B3(t)-2B0(t)+2(αc-t)(1-C)2αc(1-C)B2(t)=B0(t)-2B3(t)+2t(1-C)2αc(1-C)---(24)]]>這里C＝cosαc，S＝sinαc。也可將C2連續(xù)的C-B樣條曲線寫成矩陣形式ri(t)=Σj=03Bj(t)bi+j=12αc(1-C)sintcostt1]]>C-(1+2C)2+C-1-S2S-S0-11+2C-(1+2C)1αc-2αcCαc0bibi+1bi+2bi+3---(25)]]>其中i＝0,1,…,2n。若要使曲線通過端點處的控制點可以令初始兩個控制點和結(jié)尾末尾兩個控制點重合[34]。V0=V1Vn-1=Vn,n≥4---(26)]]>(c)C3連續(xù)草圖生成方法給定一組控制點<V0,V1,…,Vn>，以其構(gòu)成控制多邊形[35],[36]。構(gòu)造邊矢量ai＝Vi-Vi-1,i＝1,2,…,n(27)則可以計算控制多邊形每個頂點處的切矢量T0=t0(-a2)+(1-t0)a1Ti=tiai+(1-ti)ai+1Tn=tnan+(1-tn)(-an-1),i=1,2,...,n-1---(28)]]>令Vi*=Vi+|ai+1×Ti+1||Ti×Ti+1|Ti,i=1,2,...,n-1---(29)]]>及l(fā)0=|V0*-V0|li=min{|Vi*-Vi|,|Vi-Vi-1*|}ln=|Vn-Vn-1*|,i=1,2,...,n-1---(30)]]>取重構(gòu)控制點d4i=Vi-λiliTi|Ti|d4i+1=Vi-12λiliTi|Ti|d4i+2=Vi+12λiliTi|Ti|d4i+3=Vi+λiliTi|Ti|,i=0,1,...,n-1---(31)]]>其中λi為控制調(diào)節(jié)參數(shù)，且0＜λi＜1。那么由重構(gòu)控制點可構(gòu)成四次B樣條曲線ri(t)=Σj=04Ej(t)di+j,i=0,1,...,4n-1,0≤t≤1---(32)]]>其中E0(t)=(t4-4t3+6t2-4t+1)/24E1(t)=(-4t4+12t3-6t2-12t+11)/24E2(t)=(6t4-12t3-6t2+12t+11)/24E3(t)=(-4t4+4t3+6t2+4t+1)/24E4(t)=t4/24---(33)]]>S4.草圖交互虛擬域動態(tài)逆彈道參數(shù)求解方法；通過S3獲得一條完整曲線后，可將曲線表示成為只含有一個自變量的函數(shù)r(τ)，以虛擬域時刻τ為自變量，τ的取值范圍是[0,n]，(n為曲線段數(shù))，它不是時域內(nèi)的時刻。r(τ)雖然能夠反映導(dǎo)彈飛行軌跡的形狀，但不能描述軌跡上每一點狀態(tài)與時域內(nèi)時刻的對應(yīng)關(guān)系。下面根據(jù)之前介紹的方法先行建立軌跡上每一點狀態(tài)與時域內(nèi)時刻的對應(yīng)關(guān)系。假設(shè)有函數(shù)f(·)，τ和t分別為虛擬域時刻和時域時刻，τ和t有數(shù)學(xué)關(guān)系：df(τ)dt=df(τ)dτdτdt---(34)]]>這里定義λ(τ)=dτdt---(35)]]>λ(τ)為虛擬速度，由于虛擬域時刻和時域時刻是關(guān)于計算節(jié)點(下文簡稱節(jié)點)的單調(diào)遞增函數(shù)，故λ(τ)＞0，即λ(τ)具有正定性。并約定表示f(·)在時域內(nèi)求導(dǎo)，即f′(·)表示f(·)在虛擬域求導(dǎo)，即下文中其它函數(shù)任意階數(shù)的導(dǎo)數(shù)以此類推。根據(jù)(34)、(35)兩式有：f·(τ)=f′(τ)λ(τ)---(36)]]>更進(jìn)一步有：f··(τ)=d(f′(τ)λ(τ))dτdτdt=f′′(τ)λ2(τ)+f′(τ)λ(τ)λ′(τ)---(37)]]>確定虛擬域曲線r(τ)后，為方便數(shù)值計算，一般的做法是對連續(xù)函數(shù)進(jìn)行離散，這里是在τ的取值范圍取N個節(jié)點將曲線離散。節(jié)點數(shù)目可自行定義，但需注意節(jié)點數(shù)目增大時計算量也隨之增大。本發(fā)明中為方便助推段分段取N為偶數(shù)。若τ的取值范圍為[0,τf]，則可通過如下方法獲得節(jié)點之間在虛擬域的時間間隔Δτ=τfN-1---(38)]]>每一個節(jié)點對應(yīng)的虛擬域時刻可表示為τj＝τj-1+Δτ,j＝2,…,N(39)其中N為偶數(shù)。助推滑翔導(dǎo)彈采用固體發(fā)動機(jī)作為動力，時間-質(zhì)量關(guān)系m(t)、時間-推力關(guān)系T(t)視為已知。首先求出當(dāng)前節(jié)點處虛擬域基本狀態(tài)xi、x′i、x″i(i＝1,2,3)。由r(τ)以及其一階、二階導(dǎo)數(shù)得到虛擬域基本狀態(tài)，整理后結(jié)果形式如下：x1,j-1x1,j-1′x1,j-1′′x2,j-1x2,j-1′x2,j-1′′x3,j-1x3,j-1′x3,j-1′′=r(τj-1)r′(τj-1)r′′(τj-1)---(40)]]>當(dāng)τj-1∈[0,1)時r(τj-1)＝r1(τj-1)，當(dāng)τj-1∈[1,2]時r(τj-1)＝r2(τj-1)。其中x1,j-1表示x1在節(jié)點(j-1)處的值，即x1(τj-1)，其余類似下標(biāo)寫法意義相同，本文在不致混淆的前提下均在用此記法。另外需要說明，當(dāng)前節(jié)點(j-1)對應(yīng)時域時刻tj-1，在計算上一個節(jié)點(j-2)各狀態(tài)量和控制量時(本發(fā)明將這些量統(tǒng)稱節(jié)點變量)給出，下面介紹具體方法。(1)彈道傾角γj-1=arctan-x3,j-1′x1,j-1′2+x2,j-1′2---(41)]]>(2)彈道偏角ψj-1=arctanx2,j-1′x1,j-1′---(42)]]>(3)彈道傾角虛擬域?qū)?shù)γ′j-1和彈道偏角虛擬域?qū)?shù)ψ′j-1由于此步較為繁瑣，為使算式更加清晰，省略γ、xi的下標(biāo)(j-1)。γj-1′=x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32cos2γj-1---(43)]]>ψj-1′=x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1---(44)]]>(4)y方向過載ny,j-1=Vj-1g-1λj-1x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32·cos2γj-1+cosγj-1---(45)]]>(5)z方向過載nz,j-1=Vj-1g-1λj-1x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1cosγj-1---(46)]]>(6)大氣密度使用楊炳尉的大氣參數(shù)模型[38]得到大氣密度ρj-1。hj-1＝-x3,j-1(47)ρj-1＝fair(hj-1)(48)其中fair(hj-1)為大氣密度函數(shù)，已知高度可直接求取大氣密度。(7)升力本文采用攻角的一次函數(shù)表示升力系數(shù)即：CL,j-1＝l1αj-1+l2(49)那么升力Lj-1=12ρj-1Vj-12S(l1αj-1+l2)---(50)]]>其中S為導(dǎo)彈特征面積，假設(shè)這一數(shù)值在飛行過程中不變。其中l(wèi)1和l2為擬合系數(shù)。(8)攻角αj-1=ny,j-1mj-1g-12ρj-1Vj-12Sl2Tj-1+12ρj-1Vj-12Sl1---(51)]]>(9)阻力nx,j-1=Tj-1cosαj-1-12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)mj-1g---(52)]]>其中CD,j-1=d1Vj-12+d2Vj-1+Vj-1---(53)]]>Dj-1=12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)---(54)]]>d1和d2為擬合系數(shù)。(8)速度的虛擬域?qū)?shù)V′(τ)=V·(τ)/λ(τ)---(55)]]>(9)下一節(jié)點處速度Vj＝Vj-1+g(nx,j-1-sinγj-1)/λj-1Δτ(56)其中Δτ已通過(38)式給出。(10)節(jié)點(j-1)到節(jié)點j的飛行時間Δtj=Σi=13(xi,j-xi,j-1)2(Vj+Vj-1)/2---(57)]]>(11)下一節(jié)點時域時刻tj＝tj-1+Δtj(58)(12)Δt2的估計值Δt~2=TboostN-1---(59)]]>(13)λ1的估計值λ~1=ΔτΔt2≈ΔτΔt~2=τf(N-1)(N-1)Tboost=τfTboost.---(60)]]>實際計算時，λ1可在附近取值。在此附上本發(fā)明涉及的參考文獻(xiàn)：[1]LuP.InverseDynamicsApproachtoTrajectoryOptimizationforanAerospacePlane[J].JournalofGuidance,Control,andDynamics,1993,16(4):726-732.[2]ZhangYu,ChenJing,ShenLincheng.Real-timetrajectoryplanningforUCAVair-to-surfaceattackusinginversedynamicsoptimizationmethodandrecedinghorizoncontrol[J].ChineseJournalofAeronautics,2013,26(4):1038–1056.[3]閆梁,李轅,趙繼廣,杜小平.基于變節(jié)點虛擬域動態(tài)逆的軌跡實時優(yōu)化[J].航空學(xué)報,2013,34(12):2794-2803.[4]陳法龍.高超聲速滑翔飛行器彈道快速規(guī)劃研究[D].長沙:國防科學(xué)技術(shù)大學(xué)研究生院,2012.[5]YakimenkoO.DirectMethodforRapidPrototypingofNear-OptimalAircraftTrajectories[J].JournalofGuidanceControlandDynamics,2002,3(5):865-875.[6]YakimenkoO.A.,XuY.J.,GarethB.,ComputingShort-TimeAircraftManeuversUsingDirectMethodsAIAAGuidance,NavigationandControlConferenceandExhibit,Honolulu,Hawaii,18-21August2008,AIAA2008-6632.[7]LukacsJ.A.,YakimenkoO.A.Trajectory-Shape-VaryingMissileGuidanceforInterceptionofBallisticMissilesduringtheBoostPhase[C].AIAAGuidance,NavigationandControlConferenceandExhibit,HiltonHead,SouthCarolina,20-23August2007,AIAA2007-6538.[8]CowlingI.D.,YakimenkoO.A.,WhidborneJ.F.,CookeA.K.DirectMethodBasedControlSystemforanAutonomousQuadrotor[J].IntellRobotSyst,2010,60:285–316.[9]GeorgeB.,YakimenkoO.A.,MarcelloR.Real-Time6DoFGuidanceForofSpacecraftProximityManeuveringandCloseApproachwithaTumblingObject[C].AIAA/AASAstrodynamicsSpecialistConferenceToronto,OntarioCanada,2-5August2010,AIAA2010-7666.[10]AkhtarN.,WhidborneJ.F.,CookeA.K.Real-timeoptimaltechniquesforunmannedairvehiclesfuelsaving[J].PartG:AerospaceEngineering,2011,226:1315-1328.[11]MarcoC.,AlessioG.,MarcelloR.Anear-optimalguidanceforcooperativedockingmaneuvers[J].ActaAstronautica,2014,102:367–377.[12]JacopoV.,MarcelloR.,UlrichW.Performanceevaluationoftheinversedynamicsmethodforoptimalspacecraftreorientation[J].ActaAstronautica,2015,110:266-278.[13]YakimenkoO.DirectMethodforRapidPrototypingofNear-OptimalAircraftTrajectories[J].JournalofGuidanceControlandDynamics,2002,3(5):865-875.[14]錢杏芳,林瑞雄,趙亞男.導(dǎo)彈飛行力學(xué)[M].北京:北京理工大學(xué)出版社,2008:28-77.[15]賈沛然,陳克俊,何力.遠(yuǎn)程火箭彈道學(xué)[M].長沙:國防科技大學(xué)出版社,1993:57-73.[16]謝愈.復(fù)雜約束下高超聲速滑翔飛行器彈道規(guī)劃方法研究[D].長沙:國防科學(xué)技術(shù)大學(xué)研究生院,2012.[17]YLi,TMHosedales,YZSong,SGong.Free-handsketchrecognitionbymulti-kernelfeaturelearning[J].ComputerVision&ImageUnderstanding,2015,137:1-11.[18]吳玲達(dá),鄧維,張友根等.在線草圖識別研究綜述[J].計算機(jī)應(yīng)用研究,2015,32(6):1601-1607.[19]HammondT,PaulsonB.Recognizingsketchedmultistrokeprimitives[J].ACMTransonInteractiveIntelligentSystems,2011,4(1):1-34.[20]DelayeA,AnquetilE.Fuzzyrelativepositioningtemplatesforsymbolrecognition[C].Procofthe11thInternationalConferenceonDocumentAnalysisandRecognition.WashingtonDC:IEEEComputerSociety,2011:1220-122.[21]張友根,吳玲達(dá),宋漢辰等.手繪非規(guī)則軍標(biāo)圖形的結(jié)構(gòu)化識別方法[J].國防科技大學(xué)學(xué)報,2013,35(3):72-77.[22]張友根,吳玲達(dá),鄧維等.融合時空上下文的手繪筆畫圖文分類[J].電子與信息學(xué)報,2013,35(1):113-118.[23]莫蓉,常智勇.計算機(jī)輔助幾何造型技術(shù)(第二版)[M].北京:科學(xué)出版社,2009:27-91.[24]LesPiegl,WayneTiller.非均勻有理B樣條(第二版)[M].北京:清華大學(xué)出版社,2010:32-33.[25]吳曉勤,韓旭里.三次Bézier曲線的擴(kuò)展[J].工程圖學(xué)學(xué)報,2005,26(6):98-102.[26]吳曉勤.帶形狀參數(shù)的Bézier曲線[J].中國圖象圖形學(xué)報,2006,11(2):269-274.[27]李軍成.一類可調(diào)控的三次多項式曲線[J].計算機(jī)工程與科學(xué),2010,32(4):52-61.[28]DelayeA,LiuChenglin.Contextualtext/non-textstrokeclassificationinonlinehandwrittennoteswithconditionalrandomfields[J].PatternRecognition,2014,47(3):959-968.[29]PetersonEJ,StahovichTF,DoiE,etal.Groupingstrokesintoshapesinhand-drawndiagrams[C].Procofthe24thAAAIConferenceonArtificialIntelligence.2010:974-979.[30]胡偉峰,趙江洪.用戶期望意象驅(qū)動的汽車造型基因進(jìn)化[J].機(jī)械工程學(xué)報,2011,47(16):176-181.[31]FangKui,CaiFang,TanJian-Rong.C2andC3ContinuousBéziersplinecurveWithgiventangentpolygon.JournalofComputer-AidedDesign&ComputerGraphics,2000,12(5):330-332.[32]王成偉.帶有給定切線多邊形的C2連續(xù)的C-B樣條曲線[J].計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報,2001,13(12):1-4.[33]趙江洪,譚浩,譚征宇等.汽車造型設(shè)計:理論、方法與應(yīng)用[M].北京理工大學(xué)出版社,2010.[34]韓旭里,劉圣軍.二次Bézier曲線的擴(kuò)展[J].中南工業(yè)大學(xué)學(xué)報(自然科學(xué)版),2003,34(2):214-17.[35]FangKui.Closed(G2-Continuous)BéziercurveWithgiventangentpolygon.ChineseJournalofNumericalMathematics,1991,13(2):34-37.[36]王成偉,姚云.C3連續(xù)的保凸四次B樣條插值曲線[J].北京服裝學(xué)院學(xué)報(自然科學(xué)版),2002,22(1):66-70.[37]LHering.Closed(C2-andC3-Continuous)BézierandB-splinecurveWithgiventangentpolygons.CAD,1983,15(1):3-6.[38]楊炳尉.標(biāo)準(zhǔn)大氣參數(shù)的公式表示[J].宇航學(xué)報,1983,(1).當(dāng)前第1頁1 2 3