譚一廷，荊武興，高長生

(哈爾濱工業(yè)大學(xué)航天工程系，哈爾濱 150001)

0 引 言

近幾年，臨近空間高超聲速飛行器已逐漸成為各國的重點研究對象，相對傳統(tǒng)彈道導(dǎo)彈，其具有機動靈活、高速飛行等優(yōu)勢，且當(dāng)前臨近空間處于各種空天武器的攔截能力之外，成為空天防御的重大現(xiàn)實威脅之一[1]。由于傳統(tǒng)防空反導(dǎo)武器難以應(yīng)對臨近空間高超聲速目標(biāo)，研究針對臨近空間高超聲速目標(biāo)的防御技術(shù)已迫在眉睫。

現(xiàn)有防空反導(dǎo)系統(tǒng)在應(yīng)對臨近空間高超聲速武器作戰(zhàn)時面臨的主要挑戰(zhàn)可歸結(jié)為：1)臨近空間高超聲速武器具有較高的飛行速度，導(dǎo)致防御方攔截窗口較小，可攔截條件轉(zhuǎn)瞬即逝[2]。尤其在末制導(dǎo)段，可供攔截彈調(diào)整的時間有限，中制導(dǎo)需保證攔截彈中末交接班時能處于較好的攔截態(tài)勢，以減小由初、中制導(dǎo)誤差帶來的末制導(dǎo)初始時刻額外增加的機動過載；2)末制導(dǎo)段多采用動能攔截器進行直接碰撞殺傷，動能攔截器末速大，氣動加熱易造成紅外導(dǎo)引頭溫強干擾，因此多把導(dǎo)引頭安裝在攔截器側(cè)邊進行側(cè)窗探測[3]。大氣層高層及以外攔截器，中末交班對視場角的要求可以通過調(diào)整姿態(tài)來滿足，而在大氣層低層以內(nèi)，姿態(tài)的微小調(diào)整便會產(chǎn)生較大的阻尼力矩，對視窗角的要求不能簡單通過姿態(tài)調(diào)整來滿足[4]。末制導(dǎo)導(dǎo)引頭開機時，此時目標(biāo)若在視場沒有較好位置，便要求攔截彈在末制導(dǎo)初始時刻就要產(chǎn)生較大法向過載，從而保證目標(biāo)穩(wěn)定在視場內(nèi)，否則由于高超聲速目標(biāo)的高速強機動特性，會使目標(biāo)脫離攔截彈捕獲區(qū)，導(dǎo)致攔截失??；3)為應(yīng)對臨近空間高超聲速目標(biāo)橫向強機動軌跡特性，中制導(dǎo)律需具有較高的制導(dǎo)精度和計算效率，便于將來在線規(guī)劃應(yīng)用。依賴傳統(tǒng)比例導(dǎo)引或者簡單彈道預(yù)測校正制導(dǎo)[5-7]的方法往往具有較低的制導(dǎo)精度或較大的耗時，已不適用臨近高超目標(biāo)攔截特性。綜合以上分析，攔截彈中制導(dǎo)的主要任務(wù)是保證末制導(dǎo)可靠捕獲目標(biāo)，并希望導(dǎo)彈利用較少的能量沿著較為平滑的彈道向著中末交班點飛行，同時為末制導(dǎo)提供良好的初始條件以保證能精確攔截目標(biāo)同時減小自身機動。

目前綜合考慮以上三點的相關(guān)研究較少，但對其中一些具體問題已有所研究。出于對末制導(dǎo)段零控脫靶量最小考慮，文獻[2]基于縱向平面內(nèi)彈目相對運動，推導(dǎo)分析了交班零控攔截條件，并以此構(gòu)建了中末交班速度約束。文獻[8]分析了臨近防御有效零控交班區(qū)域，在僅考慮縱向平面內(nèi)彈目相對運動情況下，設(shè)計最佳交班窗口實現(xiàn)了交班點導(dǎo)引頭視場范圍最大。文獻[9-10]皆基于預(yù)測零控脫靶量進行了中制導(dǎo)律的設(shè)計，文獻[3-4]從理論上研究了反臨近空間高超聲速飛行器中末交班需用的導(dǎo)引頭關(guān)鍵技術(shù)及視角選擇問題。這些研究為中末交班時攔截彈狀態(tài)約束的建立提供了有力借鑒，但多基于二維平面，缺乏對交班視窗角和末制導(dǎo)段零控脫靶量的綜合考量。

目前最優(yōu)制導(dǎo)，滑模制導(dǎo)等方法[11-13]在滿足含終端約束的攔截制導(dǎo)中均有較成熟的研究，但其通常需要對動力學(xué)模型進行假設(shè)簡化，最終制導(dǎo)精度相對較低。近幾年偽譜法、凸優(yōu)化(Convex optimization，CVX)等軌跡規(guī)劃方法[14-15]在攔截制導(dǎo)中同樣有著較大發(fā)展，它們提高了制導(dǎo)精度以及很好解決了路徑約束問題，但在規(guī)劃效率上存在較大不足。近些年發(fā)展的模型預(yù)測靜態(tài)規(guī)劃(Model predictive static programming，MPSP)方法在綜衡計算效率和制導(dǎo)精度上有著較大的優(yōu)勢，目前在航天器變軌、飛行器軌跡跟蹤及再入制導(dǎo)等方面[16-18]有著較多的研究。為進一步提高MPSP方法在制導(dǎo)應(yīng)用上的計算效率和精度，不少學(xué)者對其進行了改進。文獻[19-20]設(shè)計了一種廣義模型預(yù)測靜態(tài)規(guī)劃方法(Generalized model predictive static programm-ing，G-MPSP)用于導(dǎo)彈再入制導(dǎo)，避免了對動力學(xué)方程的離散，提高了算法求解精度。文獻[21]設(shè)計了一種擬譜模型預(yù)測靜態(tài)規(guī)劃方法(Quasi-spectral model predictive static programming，QS-MPSP)用于攔截末制導(dǎo)，減少了靜態(tài)規(guī)劃求解變量，提高了算法求解速率。文獻[22]提出了一種時間自由的廣義擬譜模型預(yù)測靜態(tài)規(guī)劃(Generalized quasi-spectral model predictive static programming，GS-MPSP)方法用于攔截中制導(dǎo)，結(jié)合了G-MPSP方法和QS-MPSP方法兩者特點，大大提高了計算效率和求解精度，但其中涉及到譜系數(shù)初值難確定的問題，成為實際應(yīng)用的困擾。

考慮到以上問題和需求，本文首先在三維空間分析了攔截彈中末交班零控攔截條件，并利用該條件對中制導(dǎo)交班視窗角約束同末制導(dǎo)段零控脫靶量約束向交班點位置、速度方向約束進行了轉(zhuǎn)化，進而構(gòu)建了中制導(dǎo)多終端約束兩點邊值問題。其次結(jié)合Legendre偽譜法和自適應(yīng)Gauss-Lobatto積分法，推導(dǎo)了一種時間固定下的GS-MPSP方法，避免了拉格朗日乘子難求解及譜系數(shù)初值難確定的問題。將該方法用于上述中制導(dǎo)問題求解，仿真結(jié)果表明了本文方法有著高計算效率和制導(dǎo)精度。

1 中制導(dǎo)問題分析和建立

為便于后文描述，首先進行以下幾點合理假設(shè)：

假設(shè)1.考慮臨近空間滑翔目標(biāo)飛行狀態(tài)及其預(yù)報狀態(tài)完全能由預(yù)警探測系統(tǒng)獲得并給出，即給定中制導(dǎo)段飛行時間，認(rèn)為中末交班時刻目標(biāo)狀態(tài)已知。

假設(shè)2.考慮攔截中制導(dǎo)段飛行過程中攻角、側(cè)滑角保持為零，即認(rèn)為攔截彈速度方向與彈體體軸重合，則側(cè)窗探測視窗角既定為攔截彈速度與彈目視線夾角[8]。

假設(shè)3.考慮到臨近空間滑翔目標(biāo)典型飛行特點，認(rèn)為攔截末制導(dǎo)段目標(biāo)速度大于攔截彈速度。

假設(shè)4.考慮攔截彈導(dǎo)引頭最大探測距離為ρmax，當(dāng)彈目距離到達ρmax時，認(rèn)為攔截彈完成中末交班并直接進入末制導(dǎo)段。

1.1 中末交班零控攔截條件

考慮自上而下的再入攔截方式，攔截末制導(dǎo)段無控飛行，攔截高度較高，且由于末制導(dǎo)段彈目相對速度較大，攔截時間較短，因此忽略末制導(dǎo)段氣動力對攔截彈和目標(biāo)的影響，以及目標(biāo)機動情況。只受引力作用下，基于零引力差模型，攔截末制導(dǎo)段彈目相對運動學(xué)解析解形式如下[9]：

V=V0

(1)

R=V0(t-t′0)+R0

(2)

R(t′f)=V0(t′f-t′0)+R0

(3)

根據(jù)脫靶量的定義[9]，在脫靶時刻有：

R(t′f)·V(t′f)=0

(4)

由此解得t′f和脫靶時刻空間彈目相對位置矢量R(t′f)的預(yù)測值為：

(5)

(6)

要使最終R(t′f)為零，則有：

(7)

圖1所示為中末交班零控攔截示意圖，其中O-XYZ為地面發(fā)射坐標(biāo)系，M0為中制導(dǎo)起始點，VM0為攔截彈中制導(dǎo)起始點空間速度矢量，點M、T分別為交班時刻的攔截彈和目標(biāo)位置點，兩者連線即為彈目視線(Line of sight, LOS)，VM,VT分別為攔截彈和目標(biāo)空間速度矢量，兩者延長線交于點N，并構(gòu)成空間彈目相對速度矢量V0=VT-VM,γc為VM與LOS夾角，即側(cè)窗探測視窗角ξ,ηc則為VT與LOS夾角，q為LOS與水平面夾角，即視線角。在ΔM′NT′中由正弦定理可得：

(8)

圖1 零控攔截條件示意圖Fig.1 Zero effort intercept condition

(9)

(10)

式中：γ,η表示末制導(dǎo)段攔截彈和目標(biāo)速度各自與LOS夾角，R表示末制導(dǎo)段彈目距離。將式(8)代入式(9)至(10)中有：

(11)

(12)

通過以上分析可知，考慮中末交班側(cè)窗探測視窗角約束的零控攔截條件可總結(jié)為：1)交班時刻彈目速度大小約束，即需滿足目標(biāo)速度大于攔截彈速度； 2)交班時刻目標(biāo)和攔截彈速度比δ及兩者速度各自與視線夾角γc、ηc之間應(yīng)滿足的約束關(guān)系，即需滿足式(8)。文獻[2,8]等還分析了零控攔截有效交班區(qū)域，為方便后文分析推導(dǎo)，本文借用其結(jié)論，即有效零控攔截邊界為-π/2≤γc≤π/2，根據(jù)式(8)計算分析，ηc的范圍為-arcsin(1/δ)≤ηc≤arcsin(1/δ)。

1.2 中制導(dǎo)終端約束

本節(jié)目的是將中末交班視窗角及零控脫靶量兩個非等時約束轉(zhuǎn)化為中制導(dǎo)終端狀態(tài)約束，使得中制導(dǎo)問題得以簡化。

中制導(dǎo)需要在中末交班時保證目標(biāo)在側(cè)窗視場范圍內(nèi)，為簡化計算，這里將視場角范圍設(shè)為[ξmin,ξmax]，其中ξmin為導(dǎo)引頭視場錐在體軸與光軸所處平面投影的下邊界與體軸夾角，ξmax為導(dǎo)引頭視場錐在體軸與光軸所處平面投影的上邊界與體軸夾角，如圖2所示。

圖2 交班點理想視窗示意圖Fig.2 Ideal window of handover point

由于末制導(dǎo)段飛行時間較短，認(rèn)為攔截彈被動飛行過程中引力對速度變化影響較小，攔截彈速度方向基本不變。又由于滿足上述零控攔截條件時，末制導(dǎo)段視線保持平行，可近似認(rèn)為末制導(dǎo)段側(cè)窗探測視窗角ξ保持不變，幾乎維持中末交班時刻的初始視窗角度ξc。因此要使末制導(dǎo)段攔截彈對目標(biāo)始終能獲得較好探測視場，這里將視窗角約束強化為在中制導(dǎo)終端，目標(biāo)正好落在導(dǎo)引頭的視場中央，使視線與導(dǎo)引頭光軸重合，即ξc=(ξmin+ξmax)/2。

如圖2所示，記地面發(fā)射坐標(biāo)系為坐標(biāo)系g，為了使交班時刻獲得較好探測視窗角，攔截彈以高拋彈道再入的方式自上而下打擊高超滑翔目標(biāo)最佳[3]，因此認(rèn)為速度矢量VT與VM0并不平行。為減小攔截彈機動，將中制導(dǎo)終端點構(gòu)建在空間速度矢量VT與VM0所組成的平面內(nèi)。

首先以點T為原點建立直角坐標(biāo)系h′：

(13)

設(shè)Chg為發(fā)射系g到坐標(biāo)系h′的坐標(biāo)轉(zhuǎn)換矩陣,則其計算如下：

(14)

由于有效零控攔截條件中存在-arcsin1/δ≤ηc≤arcsin1/δ約束，則|ηc|<π/2，因此理想中末交班彈目相對位置ρ在坐標(biāo)系h′下可表示為：

(15)

當(dāng)視窗角|ξc|<π/2時，滿足有效零控交班約束，則可由零控攔截條件得出：

(16)

將其轉(zhuǎn)化到坐標(biāo)系g下：

(17)

最后得到終端點M理想位置RM為：

(18)

(19)

將其轉(zhuǎn)換到地面發(fā)射系g下：

(20)

最后得到終端點理想速度方向為：

(21)

通過以上分析，由式(18)及式(21)構(gòu)成轉(zhuǎn)化后的中制導(dǎo)終端點狀態(tài)約束。

1.3 中制導(dǎo)研究問題

認(rèn)為攔截彈主推力發(fā)動機不工作，攔截彈僅由直接力軌控發(fā)動機提供法向和側(cè)向機動加速度，在地面發(fā)射坐標(biāo)系下建立動力學(xué)模型：

(22)

(23)

(24)

本文在假設(shè)中末交班時刻目標(biāo)狀態(tài)已知情況下，利用零控攔截條件將交班時視窗角約束及末制導(dǎo)段零控脫靶量約束轉(zhuǎn)化為中制導(dǎo)終端速度方向約束和位置約束，從而將中制導(dǎo)問題進一步簡化為：在已知飛行時間Δt=tf-t0內(nèi)，以能耗最省為指標(biāo)，給出最佳控制序列u=[ay,az]T，使得滿足終端速度方向約束式(21)及位置約束式(18)。

2 基于自適應(yīng)正交配點的GS-MPSP算法

上述中制導(dǎo)問題可歸結(jié)為典型的兩點邊值問題，綜合計算效率和制導(dǎo)精度兩方面考慮，本節(jié)結(jié)合G-MPSP方法和QS-MPSP方法特點推導(dǎo)了固定時間下的GS-MPSP算法，相關(guān)推導(dǎo)及物理量的計算如下。

2.1 GS-MPSP算法

一般多輸入多輸出(MIMO)非線性系統(tǒng)狀態(tài)空間表達式如下：

(25)

Y(t)=h(X,t)

(26)

式中：X∈Rn,U∈Rm,Y∈Rp，分別代表狀態(tài)變量，控制變量及輸出變量。終端時刻tf對應(yīng)的輸出偏差可近似表示為Y(X(tf))的變分如下：

ΔYf=Y(tf)-Y(tf)cf≈δY(X(tf))

(27)

首先對式(25)左右同時乘以W(t)，再左右同時進行t0→tf積分，最后加上Y(X(tf))整理可得：

(28)

式中：W(t)∈Rp×n為權(quán)重矩陣。W(t)將系統(tǒng)動力學(xué)映射到維數(shù)更小的輸出空間，從而減小了狀態(tài)方程維數(shù)，減少了后續(xù)靜態(tài)規(guī)劃時求解變量個數(shù)，以提高計算效率。

對上式最后一項按分部積分得：

(29)

將式(29)代入式(28)得：

Y(X(tf))=Y(X(tf))-W(tf)X(tf)+W(t0)·

(30)

對式(30)左右兩邊同時取變分，并結(jié)合式(27)整理得：

(31)

分析上式，可通過合適選擇W(t)，使得以消除與?X(t)相關(guān)的系數(shù)，從而權(quán)重矩陣W(t)存在具有終端時刻tf相關(guān)邊界條件的微分方程如下：

(32)

(33)

將式(32)和式(33)代入式(31)整理得：

(34)

對于確定的初始狀態(tài)，有δX(t)|t=t0=0，結(jié)合式(34)可得：

(35)

式中：Bs(t)∈Rp×m，稱為靈敏度矩陣，表達式如下：

(36)

式(35)反映了G-MPSP方法中，終端輸出偏差ΔYf與t0到tf飛行過程中連續(xù)控制δU(t)的關(guān)系，其由靈敏度矩陣Bs(t)關(guān)聯(lián)?；谶@種關(guān)系，可通過構(gòu)建靜態(tài)規(guī)劃問題預(yù)測更新每一時刻控制偏差δU(t)，并補償?shù)缴弦徊筋A(yù)測控制U(t)上，從而達到消除終端輸出偏差目的。

假設(shè)每步預(yù)測過程的控制更新為：

Ul+1(t)=Ul(t)+δUl(t)

(37)

式中：Ul(t)、δUl(t)為第l步迭代過程中控制量和控制偏差，Ul+1(t)為第l+1步迭代過程中控制量，從而式(35)可寫成如下形式：

(38)

令，

(39)

式中：Bλ∈Rp，稱為過程矩陣，則:

(40)

為了減少優(yōu)化變量的數(shù)目，采用文獻[21]中提出的控制向量譜表示方法，即將控制向量表示為一些基本譜函數(shù)的加權(quán)和形式：

(41)

式中：μj∈Rm(j=1,2,…,Np)，表示第j個譜函數(shù)的系數(shù)向量，Np代表譜函數(shù)個數(shù)，Pj(t)代表譜函數(shù)基，一般可選擇為Legendre多項式、Chebychev多項式等形式。

假設(shè)每步預(yù)測過程的譜函數(shù)系數(shù)向量的更新為：

(42)

(43)

將式(43)代入式(40)可得：

(44)

令，

j=1,2,…,Np

(45)

式中：Aj∈Rp×m，為了與MPSP、QS-MPSP中靈敏度矩陣加以區(qū)分，稱其為譜靈敏度矩陣。最終終端時刻輸出偏差ΔYf表達式如下：

(46)

下面通過建立靜態(tài)規(guī)劃問題完成對每步預(yù)測過程中控制量Ul+1(t)的更新。

注意到控制量Ul+1(t)的上邊界存在如下不等式：

(47)

由于每一時刻譜函數(shù)Pj(t)大小固定，因此最小化μ=[μ1,μ2, …,μNp]，便可達到最小控制消耗目的。構(gòu)建指標(biāo)函數(shù)[21]如下：

(48)

式中：Rj(t)∈Rm×m(j=1,2,…,Np)，通常選為單位矩陣。引入拉格朗日乘子λ∈Rp，結(jié)合式(46)，增廣指標(biāo)函數(shù)構(gòu)建如下：

(49)

根據(jù)Hamilton邊值問題一階最優(yōu)性條件有：

(50)

式中：

(51)

(52)

將式(52)代入式(46)得：

(53)

式中：Aλ∈Rp×p，表達式如下：

(54)

進而推導(dǎo)出拉格朗日乘子λ為：

(55)

將式(55)代入式(52)得到譜函數(shù)系數(shù)的更新為：

j=1,2,…,Np

(56)

將式(56)代入式(43)得到控制量更新為：

(57)

2.2 Legendre偽譜法計算加權(quán)矩陣W(t)

上一節(jié)提到控制量的更新需要計算譜靈敏度矩陣Aj，而矩陣Aj的計算首先需要確定權(quán)重矩陣W(t)。權(quán)重矩陣W(t)受微分方程式(32)和終端邊界條件式(33)約束，為一初值問題，通?？捎赡嫦蚍e分求解[20]，本文借助偽譜法解微分方程思想，對權(quán)重矩陣W(t)進行解算，以提高計算效率。

假設(shè)時域[t0,tf]等分為Ns個時間區(qū)間，即t0

i=0,1,2,…,Ns-1

(58)

(59)

令，

(60)

則有：

(61)

令，

(62)

式(62)可寫成如下形式：

(63)

(64)

對式(64)求導(dǎo)得：

(65)

其在區(qū)間t∈[ti,ti+1]內(nèi)前N-1個配點上的微分方程如下：

(66)

令，

(67)

式中：D∈R(N-1)×N，微分矩陣D中各元素計算如下：

(68)

則式(66)可寫成如下形式：

(69)

(70)

將式(70)寫成如下形式：

(71)

(72)

2.3 自適應(yīng)Gauss-Lobatto積分計算矩陣Aj, Bλ

在確定加權(quán)矩陣W(t)后，按照式(39)和式(45)，矩陣Aj,Bλ可依靠時域[t0,tf]內(nèi)的模型積分[19-20]而獲得。為使本文所設(shè)計制導(dǎo)方法適應(yīng)于將來在線規(guī)劃應(yīng)用，出于平衡計算精度和計算效率的目的，本文采用自適應(yīng)Gauss-Lobatto積分進行矩陣Aj,Bλ計算。在積分精度不滿足要求時，增加分段個數(shù)，而每個分段區(qū)間內(nèi)配點數(shù)不變來保證積分精度。

首先將時域區(qū)間[t0,tf]等分為Ns個時間區(qū)間，即t0

(73)

則譜靈敏度矩陣Aj和矩陣Bλ由各個時域區(qū)間內(nèi)的Gauss-Lobatto積分值累加得到，即:

hkBs(τk)Pj(τk) (j=1,2,…,Np)

(74)

hkBs(τk)Ul(τk) (j=1,2,…,Np)

(75)

式中：hk為積分權(quán)重，其計算如下：

(76)

自適應(yīng)積分法就是按照積分區(qū)間上被積函數(shù)變化的劇烈程度，動態(tài)調(diào)整積分節(jié)點的分布密度，以達到節(jié)約計算成本的目的。綜合積分精度和計算效率兩者考慮，本文采用以下區(qū)間數(shù)劃分準(zhǔn)則：

記區(qū)間[ti,ti+1]長度為h，區(qū)間[ti,ti+1]上采用N個LGL配點進行Gauss-Lobatto積分的結(jié)果為R1(ti,ti+1)，將區(qū)間[ti,ti+1]對分，在區(qū)間[ti,ti+h/2], [ti+h/2,ti+1]上分別采用同樣N個LGL配點的復(fù)化Gauss-Lobatto積分結(jié)果為R2(ti,ti+1)。設(shè)整個區(qū)間[t0,tf]上要求的積分精度為ε，若:

(77)

則區(qū)間[ti,ti+1]不再對分，否則將區(qū)間[ti,ti+1]對分，并令ti+1=ti+h/2，這時時域[t0,tf]上的區(qū)間個數(shù)Ns增加為Ns+1個，并重復(fù)上述過程。

3 基于GS-MPSP的中制導(dǎo)律設(shè)計

首先為了保證所有狀態(tài)變量具有類似的數(shù)值變化范圍，將攔截彈狀態(tài)變量進行歸一化，并將攔截彈動力學(xué)模型式(22)和式(23)進行歸一化[20]。

利用GS-MPSP算法進行制導(dǎo)律設(shè)計時，需要對初始控制量進行猜測，本文運用比例導(dǎo)引法(PN)給出控制量初始解。文獻[20]給出了PN制導(dǎo)俯仰、偏航方向的指令加速度如下：

(78)

(79)

如圖3所示，本文基于GS-MPSP規(guī)劃的中制導(dǎo)實現(xiàn)流程描述如下：

圖3 基于GS-MPSP規(guī)劃的中制導(dǎo)實現(xiàn)流程Fig.3 Process of mid-guidance using GS-MPSP

1)假設(shè)給定中制導(dǎo)飛行時間tf，并由預(yù)警探測系統(tǒng)給出tf時間后的目標(biāo)預(yù)報狀態(tài)RT,VT，結(jié)合式(13)至式(21)，計算終端點期望狀態(tài)YNf。

2)采用PN制導(dǎo)律式(78)和式(79)計算控制初始解Ul。

3)采用四階龍格庫塔進行攔截彈彈道積分，計算終端點狀態(tài)偏差，若滿足精度則輸出控制序列Ul，否則進入下一步，這一步也對應(yīng)著GS-MPSP方法的模型預(yù)測模塊。

4)結(jié)合式(59)、式(68)、式(70)和式(72)計算矩陣W(t)。

5)結(jié)合式(74)至式(77)計算矩陣Aj,Bλ。

6)結(jié)合式(57)進行控制更新Ul+1(t)，令Ul(t)=Ul+1(t)并返回第3)步，直到滿足精度結(jié)束，這一步也對應(yīng)著GS-MPSP方法的靜態(tài)規(guī)劃模塊。

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

為驗證本文所設(shè)計方法的有效性，開展以下兩種情形進行仿真。設(shè)定仿真中攔截彈導(dǎo)引頭最大探測距離為ρmax=50 km，并假設(shè)中制導(dǎo)過程中目標(biāo)在交接班時刻的狀態(tài)信息依靠遠(yuǎn)程預(yù)警系統(tǒng)或攔截彈的火控雷達預(yù)測得到。仿真平臺相關(guān)參數(shù)為：Intel Core i5 3.00 GHz處理器，window10操作系統(tǒng)，matlab18b作仿真軟件。

4.1 情形一

情形一設(shè)定合理中制導(dǎo)飛行時間為76 s，考慮同一中末交班目標(biāo)與攔截彈速度比δ=1.31以及側(cè)窗探測視場范圍3°～33°，自動駕駛儀一階時間常數(shù)取0.3， GS-MPSP仿真中取譜函數(shù)個數(shù)Np=4，時域[t0,tf]內(nèi)等分區(qū)間初始個數(shù)為Ns=2，各分段時間區(qū)間內(nèi)選擇的LGL配點數(shù)為N=4，其他仿真場景參數(shù)以及利用式(13)至式(21)計算所得的終端期望狀態(tài)結(jié)果如表1所示，仿真結(jié)果如表2和圖4至圖8所示，其中表2中Δrf表示終端點位置偏差，Δθmf表示終端彈道傾角偏差，Δφvmf表示終端彈道偏角偏差，ξf表示交班點視窗角。

情形一采用GS-MPSP、MPSP[17]、CVX軌跡規(guī)劃[14]以及PN制導(dǎo)[20]四種方法進行了對比分析。其中前三種方法對終端位置[xf,yf,zf]的收斂閾值均設(shè)為0.5 m，對終端彈道傾角或彈道偏角偏差收斂閾值均設(shè)為0.5°，并設(shè)置迭代步數(shù)不超過3步，CVX軌跡規(guī)劃同MPSP皆采用歐拉法對動力學(xué)方程進行離散。

表1 仿真場景參數(shù)設(shè)置Table 1 Parameter setting of simulation scene

表2 情形一仿真結(jié)果統(tǒng)計Table 2 Case 1 Statistics of simulation results

如圖5、圖6和表2所示，PN制導(dǎo)雖能成功到達期望終端點，但未能滿足終端速度方向約束。以PN制導(dǎo)結(jié)果作為猜測初值用于前三種方法仿真，本文提出的GS-MPSP方法同MPSP方法和CVX軌跡規(guī)劃法在準(zhǔn)確到達期望終端點的同時，滿足了終端速度方向約束，且中末交班時刻視窗角在側(cè)窗探測視場中央附近，給末制導(dǎo)段提供了較好的攔截態(tài)勢，其中GS-MPSP制導(dǎo)精度相對較高。

從圖4至圖6可以看出，GS-MPSP同MPSP產(chǎn)生的軌跡及指令加速度基本一致，這也間接說明GS-MPSP方法原理上與MPSP方法之間的等效性。但也存在著輕微的不同，這是由于兩者在性能指標(biāo)的構(gòu)建上存在不同，本文設(shè)計的GS-MPSP方法以譜函數(shù)系數(shù)的加權(quán)和最小為指標(biāo)，而MPSP方法以控制量的加權(quán)和最小為指標(biāo)，在取加權(quán)矩陣都是單位陣的情況下，前者在制導(dǎo)中期有著相對較小的控制成本，相比之下，利用CVX軌跡規(guī)劃方法的控制成本明顯增大，且軌跡波動相對較大。

圖4 縱向和側(cè)向指令加速度曲線Fig.4 Longitudinal and lateral command acceleration change

圖5 側(cè)向和縱向平面內(nèi)二維軌跡Fig.5 Two-dimensional trajectory in lateral and longitudinal planes

圖6 彈道傾角和偏角變化曲線Fig.6 Ballistic inclination and ballistic deflection change

如前幾節(jié)所述，譜靈敏度矩陣Aj和過程矩陣Bλ的計算是GS-MPSP方法的重要組成部分。圖7和圖8展示了第一次迭代中采用自適應(yīng)正交配點法計算矩陣Aj,Bλ時，其中各元素的計算結(jié)果，橫坐標(biāo)為元素序號，縱坐標(biāo)為元素值，其中矩陣Aj個數(shù)為4，矩陣Bλ個數(shù)為1，矩陣Aj為5行2列，矩陣Bλ為5行1列，因此矩陣Aj中元素個數(shù)為10，矩陣Bλ中元素個數(shù)為5。

將采用四階龍格庫塔法解算式(32)、式(33)、式(39)和式(45)所得的矩陣Aj,Bλ中各元素值作為標(biāo)準(zhǔn)值，并與本文方法進行對比可看出，兩者相對誤差較小，這表明利用自適應(yīng)正交配點法能保證動力學(xué)方程較高積分精度。由于GS-MPSP方法的收斂性和收斂精度與譜靈敏度矩陣的計算精度相關(guān)，從終端狀態(tài)偏差結(jié)果來看，基于自適應(yīng)正交配點的GS-MPSP方法相比采用歐拉積分法對動力學(xué)方程進行離散的MPSP方法和CVX規(guī)劃方法具有更高制導(dǎo)精度。

圖7 各譜靈敏度矩陣計算值Fig.7 Each spectral sensitivity matrix value

圖8 過程矩陣Bλ計算值Fig.8 Process matrix Bλ value

表3記錄了GS-MPSP、MPSP和CVX三種方法迭代一步所占的計算耗時，其中GS-MPSP方法計算效率最高，大約為MPSP方法的二分之一，同時是CVX規(guī)劃方法的二十分之一，這是由于本文采用自適應(yīng)正交配點法，在保證動力學(xué)方程積分精度的同時較大地減少了離散點個數(shù)，從而較大地提高了計算效率，對于將來攔截彈在線規(guī)劃等方面應(yīng)用具有較大優(yōu)勢。

綜合情形一仿真結(jié)果分析，本文基于自適應(yīng)正交配點的GS-MPSP方法在滿足制導(dǎo)精度的同時控制消耗最小，且計算效率最高。

表3 計算耗時統(tǒng)計Table 3 Calculation time consumption statistics

4.2 情形二

本節(jié)在4.1節(jié)仿真條件的基礎(chǔ)上，考慮攔截彈初始狀態(tài)擾動、氣動阻力系數(shù)攝動及質(zhì)量偏差，偏差項設(shè)置見表4。利用蒙特卡洛打靶仿真200次，統(tǒng)計結(jié)果見表5和圖9，其中表5中各終端參數(shù)物理意義同表2。由于考慮到本文所設(shè)計制導(dǎo)律在實際應(yīng)用時的實時性問題，出于對收斂精度和計算效率的衡量，同樣將GS-MPSP迭代步數(shù)限制在3步以內(nèi)，以滿足實時性要求，進而評估在保證實時性前提下的制導(dǎo)律的魯棒性能。

表4 仿真偏差設(shè)定Table 4 Simulation deviation setting

圖9 終端狀態(tài)偏差分布圖Fig.9 Distribution of terminal state deviation

表5 打靶結(jié)果統(tǒng)計Table 5 Shooting result statistics

可以看出，在考慮初始狀態(tài)擾動和參數(shù)攝動的情況下，到達期望終端約束的偏差在合理范圍內(nèi)，并保證中末交班時視窗角在視場中央，表明本文設(shè)計的GS-MPSP制導(dǎo)方法具有一定魯棒性。

在考慮其他分系統(tǒng)理想工作、攔截末制導(dǎo)段攔截時間較短及彈目速度變化較小情況下，情形一和情形二驗證了本文考慮交班視窗角約束的中制導(dǎo)算法的有效性。在后續(xù)研究中，還需綜合考慮目標(biāo)強機動、目標(biāo)預(yù)報誤差及攔截彈控制系統(tǒng)動態(tài)響應(yīng)特性等對制導(dǎo)律進行進一步設(shè)計。

5 結(jié) 論

本文為解決臨近空間防御中考慮零控交班及交班視窗角約束的中制導(dǎo)問題提供了一種制導(dǎo)精度和計算效率較高的中制導(dǎo)方法，具體總結(jié)如下：

1)推導(dǎo)分析了攔截末制導(dǎo)段零控攔截條件，此條件由中末交班時刻目標(biāo)和攔截彈速度比及兩者速度分別與視線夾角確定，并利用此條件提出了一種將中末交班側(cè)窗視窗角約束同零控脫靶量約束轉(zhuǎn)化為中制導(dǎo)終端點位置、速度方向約束的方法，為滿足多約束的中制導(dǎo)設(shè)計提供了更簡便的思路。

2)提出了一種固定時間條件下的廣義擬譜模型預(yù)測靜態(tài)規(guī)劃算法。通過引入低維權(quán)重矩陣和控制量的譜表達式以降低靜態(tài)規(guī)劃問題的維數(shù)，并推導(dǎo)得到了控制更新顯示表達式，該表達式直接由依靠起始到終端時刻連續(xù)積分的譜靈敏度矩陣、過程矩陣和終端偏差確定，而拉格朗日乘子和譜系數(shù)僅作為迭代過程中的中間量，解決了拉格朗日乘子難求解，以及譜系數(shù)初值難確定的問題。相比PN制導(dǎo)、MPSP方法及CVX軌跡規(guī)劃方法，GS-MPSP方法控制成本較小且彈道更平滑，制導(dǎo)精度更高。

3)提出了基于自適應(yīng)正交配點法的譜靈敏度矩陣和過程矩陣的計算方法。利用自適應(yīng)Gauss-Lobatto積分法實現(xiàn)時間分區(qū)間個數(shù)自適應(yīng)分配來計算譜靈敏度矩陣和過程矩陣，從而保證計算精度；用Legendre偽譜法實現(xiàn)了權(quán)重矩陣的解算。相比CVX軌跡規(guī)劃方法及MPSP方法，有著更高的求解效率，對將來實時在線規(guī)劃應(yīng)用有著較強適用性。

4)在考慮攔截彈初始狀態(tài)擾動和參數(shù)攝動的條件下，本文方法在保證實時性要求情況下，依然能夠滿足中末交班制導(dǎo)要求，具有較強魯棒性。