谭一廷，荆武兴，高长生

(哈尔滨工业大学航天工程系，哈尔滨 150001)

0 引 言

近几年，临近空间高超声速飞行器已逐渐成为各国的重点研究对象，相对传统弹道导弹，其具有机动灵活、高速飞行等优势，且当前临近空间处于各种空天武器的拦截能力之外，成为空天防御的重大现实威胁之一[1]。由于传统防空反导武器难以应对临近空间高超声速目标，研究针对临近空间高超声速目标的防御技术已迫在眉睫。

现有防空反导系统在应对临近空间高超声速武器作战时面临的主要挑战可归结为：1)临近空间高超声速武器具有较高的飞行速度，导致防御方拦截窗口较小，可拦截条件转瞬即逝[2]。尤其在末制导段，可供拦截弹调整的时间有限，中制导需保证拦截弹中末交接班时能处于较好的拦截态势，以减小由初、中制导误差带来的末制导初始时刻额外增加的机动过载；2)末制导段多采用动能拦截器进行直接碰撞杀伤，动能拦截器末速大，气动加热易造成红外导引头温强干扰，因此多把导引头安装在拦截器侧边进行侧窗探测[3]。大气层高层及以外拦截器，中末交班对视场角的要求可以通过调整姿态来满足，而在大气层低层以内，姿态的微小调整便会产生较大的阻尼力矩，对视窗角的要求不能简单通过姿态调整来满足[4]。末制导导引头开机时，此时目标若在视场没有较好位置，便要求拦截弹在末制导初始时刻就要产生较大法向过载，从而保证目标稳定在视场内，否则由于高超声速目标的高速强机动特性，会使目标脱离拦截弹捕获区，导致拦截失败；3)为应对临近空间高超声速目标横向强机动轨迹特性，中制导律需具有较高的制导精度和计算效率，便于将来在线规划应用。依赖传统比例导引或者简单弹道预测校正制导[5-7]的方法往往具有较低的制导精度或较大的耗时，已不适用临近高超目标拦截特性。综合以上分析，拦截弹中制导的主要任务是保证末制导可靠捕获目标，并希望导弹利用较少的能量沿着较为平滑的弹道向着中末交班点飞行，同时为末制导提供良好的初始条件以保证能精确拦截目标同时减小自身机动。

目前综合考虑以上三点的相关研究较少，但对其中一些具体问题已有所研究。出于对末制导段零控脱靶量最小考虑，文献[2]基于纵向平面内弹目相对运动，推导分析了交班零控拦截条件，并以此构建了中末交班速度约束。文献[8]分析了临近防御有效零控交班区域，在仅考虑纵向平面内弹目相对运动情况下，设计最佳交班窗口实现了交班点导引头视场范围最大。文献[9-10]皆基于预测零控脱靶量进行了中制导律的设计，文献[3-4]从理论上研究了反临近空间高超声速飞行器中末交班需用的导引头关键技术及视角选择问题。这些研究为中末交班时拦截弹状态约束的建立提供了有力借鉴，但多基于二维平面，缺乏对交班视窗角和末制导段零控脱靶量的综合考量。

目前最优制导，滑模制导等方法[11-13]在满足含终端约束的拦截制导中均有较成熟的研究，但其通常需要对动力学模型进行假设简化，最终制导精度相对较低。近几年伪谱法、凸优化(Convex optimization，CVX)等轨迹规划方法[14-15]在拦截制导中同样有着较大发展，它们提高了制导精度以及很好解决了路径约束问题，但在规划效率上存在较大不足。近些年发展的模型预测静态规划(Model predictive static programming，MPSP)方法在综衡计算效率和制导精度上有着较大的优势，目前在航天器变轨、飞行器轨迹跟踪及再入制导等方面[16-18]有着较多的研究。为进一步提高MPSP方法在制导应用上的计算效率和精度，不少学者对其进行了改进。文献[19-20]设计了一种广义模型预测静态规划方法(Generalized model predictive static programm-ing，G-MPSP)用于导弹再入制导，避免了对动力学方程的离散，提高了算法求解精度。文献[21]设计了一种拟谱模型预测静态规划方法(Quasi-spectral model predictive static programming，QS-MPSP)用于拦截末制导，减少了静态规划求解变量，提高了算法求解速率。文献[22]提出了一种时间自由的广义拟谱模型预测静态规划(Generalized quasi-spectral model predictive static programming，GS-MPSP)方法用于拦截中制导，结合了G-MPSP方法和QS-MPSP方法两者特点，大大提高了计算效率和求解精度，但其中涉及到谱系数初值难确定的问题，成为实际应用的困扰。

考虑到以上问题和需求，本文首先在三维空间分析了拦截弹中末交班零控拦截条件，并利用该条件对中制导交班视窗角约束同末制导段零控脱靶量约束向交班点位置、速度方向约束进行了转化，进而构建了中制导多终端约束两点边值问题。其次结合Legendre伪谱法和自适应Gauss-Lobatto积分法，推导了一种时间固定下的GS-MPSP方法，避免了拉格朗日乘子难求解及谱系数初值难确定的问题。将该方法用于上述中制导问题求解，仿真结果表明了本文方法有着高计算效率和制导精度。

1 中制导问题分析和建立

为便于后文描述，首先进行以下几点合理假设：

假设1.考虑临近空间滑翔目标飞行状态及其预报状态完全能由预警探测系统获得并给出，即给定中制导段飞行时间，认为中末交班时刻目标状态已知。

假设2.考虑拦截中制导段飞行过程中攻角、侧滑角保持为零，即认为拦截弹速度方向与弹体体轴重合，则侧窗探测视窗角既定为拦截弹速度与弹目视线夹角[8]。

假设3.考虑到临近空间滑翔目标典型飞行特点，认为拦截末制导段目标速度大于拦截弹速度。

假设4.考虑拦截弹导引头最大探测距离为ρmax，当弹目距离到达ρmax时，认为拦截弹完成中末交班并直接进入末制导段。

1.1 中末交班零控拦截条件

考虑自上而下的再入拦截方式，拦截末制导段无控飞行，拦截高度较高，且由于末制导段弹目相对速度较大，拦截时间较短，因此忽略末制导段气动力对拦截弹和目标的影响，以及目标机动情况。只受引力作用下，基于零引力差模型，拦截末制导段弹目相对运动学解析解形式如下[9]：

V=V0

(1)

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

(2)

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

(3)

根据脱靶量的定义[9]，在脱靶时刻有：

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

(4)

由此解得t′f和脱靶时刻空间弹目相对位置矢量R(t′f)的预测值为：

(5)

(6)

要使最终R(t′f)为零，则有：

(7)

图1所示为中末交班零控拦截示意图，其中O-XYZ为地面发射坐标系，M0为中制导起始点，VM0为拦截弹中制导起始点空间速度矢量，点M、T分别为交班时刻的拦截弹和目标位置点，两者连线即为弹目视线(Line of sight, LOS)，VM,VT分别为拦截弹和目标空间速度矢量，两者延长线交于点N，并构成空间弹目相对速度矢量V0=VT-VM,γc为VM与LOS夹角，即侧窗探测视窗角ξ,ηc则为VT与LOS夹角，q为LOS与水平面夹角，即视线角。在ΔM′NT′中由正弦定理可得：

(8)

图1 零控拦截条件示意图Fig.1 Zero effort intercept condition

(9)

(10)

式中：γ,η表示末制导段拦截弹和目标速度各自与LOS夹角，R表示末制导段弹目距离。将式(8)代入式(9)至(10)中有：

(11)

(12)

通过以上分析可知，考虑中末交班侧窗探测视窗角约束的零控拦截条件可总结为：1)交班时刻弹目速度大小约束，即需满足目标速度大于拦截弹速度； 2)交班时刻目标和拦截弹速度比δ及两者速度各自与视线夹角γc、ηc之间应满足的约束关系，即需满足式(8)。文献[2,8]等还分析了零控拦截有效交班区域，为方便后文分析推导，本文借用其结论，即有效零控拦截边界为-π/2≤γc≤π/2，根据式(8)计算分析，ηc的范围为-arcsin(1/δ)≤ηc≤arcsin(1/δ)。

1.2 中制导终端约束

本节目的是将中末交班视窗角及零控脱靶量两个非等时约束转化为中制导终端状态约束，使得中制导问题得以简化。

中制导需要在中末交班时保证目标在侧窗视场范围内，为简化计算，这里将视场角范围设为[ξmin,ξmax]，其中ξmin为导引头视场锥在体轴与光轴所处平面投影的下边界与体轴夹角，ξmax为导引头视场锥在体轴与光轴所处平面投影的上边界与体轴夹角，如图2所示。

图2 交班点理想视窗示意图Fig.2 Ideal window of handover point

由于末制导段飞行时间较短，认为拦截弹被动飞行过程中引力对速度变化影响较小，拦截弹速度方向基本不变。又由于满足上述零控拦截条件时，末制导段视线保持平行，可近似认为末制导段侧窗探测视窗角ξ保持不变，几乎维持中末交班时刻的初始视窗角度ξc。因此要使末制导段拦截弹对目标始终能获得较好探测视场，这里将视窗角约束强化为在中制导终端，目标正好落在导引头的视场中央，使视线与导引头光轴重合，即ξc=(ξmin+ξmax)/2。

如图2所示，记地面发射坐标系为坐标系g，为了使交班时刻获得较好探测视窗角，拦截弹以高抛弹道再入的方式自上而下打击高超滑翔目标最佳[3]，因此认为速度矢量VT与VM0并不平行。为减小拦截弹机动，将中制导终端点构建在空间速度矢量VT与VM0所组成的平面内。

首先以点T为原点建立直角坐标系h′：

(13)

设Chg为发射系g到坐标系h′的坐标转换矩阵,则其计算如下：

(14)

由于有效零控拦截条件中存在-arcsin1/δ≤ηc≤arcsin1/δ约束，则|ηc|<π/2，因此理想中末交班弹目相对位置ρ在坐标系h′下可表示为：

(15)

当视窗角|ξc|<π/2时，满足有效零控交班约束，则可由零控拦截条件得出：

(16)

将其转化到坐标系g下：

(17)

最后得到终端点M理想位置RM为：

(18)

(19)

将其转换到地面发射系g下：

(20)

最后得到终端点理想速度方向为：

(21)

通过以上分析，由式(18)及式(21)构成转化后的中制导终端点状态约束。

1.3 中制导研究问题

认为拦截弹主推力发动机不工作，拦截弹仅由直接力轨控发动机提供法向和侧向机动加速度，在地面发射坐标系下建立动力学模型：

(22)

(23)

(24)

本文在假设中末交班时刻目标状态已知情况下，利用零控拦截条件将交班时视窗角约束及末制导段零控脱靶量约束转化为中制导终端速度方向约束和位置约束，从而将中制导问题进一步简化为：在已知飞行时间Δt=tf-t0内，以能耗最省为指标，给出最佳控制序列u=[ay,az]T，使得满足终端速度方向约束式(21)及位置约束式(18)。

2 基于自适应正交配点的GS-MPSP算法

上述中制导问题可归结为典型的两点边值问题，综合计算效率和制导精度两方面考虑，本节结合G-MPSP方法和QS-MPSP方法特点推导了固定时间下的GS-MPSP算法，相关推导及物理量的计算如下。

2.1 GS-MPSP算法

一般多输入多输出(MIMO)非线性系统状态空间表达式如下：

(25)

Y(t)=h(X,t)

(26)

式中：X∈Rn,U∈Rm,Y∈Rp，分别代表状态变量，控制变量及输出变量。终端时刻tf对应的输出偏差可近似表示为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为权重矩阵。W(t)将系统动力学映射到维数更小的输出空间，从而减小了状态方程维数，减少了后续静态规划时求解变量个数，以提高计算效率。

对上式最后一项按分部积分得：

(29)

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

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

(30)

对式(30)左右两边同时取变分，并结合式(27)整理得：

(31)

分析上式，可通过合适选择W(t)，使得以消除与∂X(t)相关的系数，从而权重矩阵W(t)存在具有终端时刻tf相关边界条件的微分方程如下：

(32)

(33)

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

(34)

对于确定的初始状态，有δX(t)|t=t0=0，结合式(34)可得：

(35)

式中：Bs(t)∈Rp×m，称为灵敏度矩阵，表达式如下：

(36)

式(35)反映了G-MPSP方法中，终端输出偏差ΔYf与t0到tf飞行过程中连续控制δU(t)的关系，其由灵敏度矩阵Bs(t)关联。基于这种关系，可通过构建静态规划问题预测更新每一时刻控制偏差δU(t)，并补偿到上一步预测控制U(t)上，从而达到消除终端输出偏差目的。

假设每步预测过程的控制更新为：

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

(37)

式中：Ul(t)、δUl(t)为第l步迭代过程中控制量和控制偏差，Ul+1(t)为第l+1步迭代过程中控制量，从而式(35)可写成如下形式：

(38)

令，

(39)

式中：Bλ∈Rp，称为过程矩阵，则:

(40)

为了减少优化变量的数目，采用文献[21]中提出的控制向量谱表示方法，即将控制向量表示为一些基本谱函数的加权和形式：

(41)

式中：μj∈Rm(j=1,2,…,Np)，表示第j个谱函数的系数向量，Np代表谱函数个数，Pj(t)代表谱函数基，一般可选择为Legendre多项式、Chebychev多项式等形式。

假设每步预测过程的谱函数系数向量的更新为：

(42)

(43)

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

(44)

令，

j=1,2,…,Np

(45)

式中：Aj∈Rp×m，为了与MPSP、QS-MPSP中灵敏度矩阵加以区分，称其为谱灵敏度矩阵。最终终端时刻输出偏差ΔYf表达式如下：

(46)

下面通过建立静态规划问题完成对每步预测过程中控制量Ul+1(t)的更新。

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

(47)

由于每一时刻谱函数Pj(t)大小固定，因此最小化μ=[μ1,μ2, …,μNp]，便可达到最小控制消耗目的。构建指标函数[21]如下：

(48)

式中：Rj(t)∈Rm×m(j=1,2,…,Np)，通常选为单位矩阵。引入拉格朗日乘子λ∈Rp，结合式(46)，增广指标函数构建如下：

(49)

根据Hamilton边值问题一阶最优性条件有：

(50)

式中：

(51)

(52)

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

(53)

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

(54)

进而推导出拉格朗日乘子λ为：

(55)

将式(55)代入式(52)得到谱函数系数的更新为：

j=1,2,…,Np

(56)

将式(56)代入式(43)得到控制量更新为：

(57)

2.2 Legendre伪谱法计算加权矩阵W(t)

上一节提到控制量的更新需要计算谱灵敏度矩阵Aj，而矩阵Aj的计算首先需要确定权重矩阵W(t)。权重矩阵W(t)受微分方程式(32)和终端边界条件式(33)约束，为一初值问题，通常可由逆向积分求解[20]，本文借助伪谱法解微分方程思想，对权重矩阵W(t)进行解算，以提高计算效率。

假设时域[t0,tf]等分为Ns个时间区间，即t0

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

(58)

(59)

令，

(60)

则有：

(61)

令，

(62)

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

(63)

(64)

对式(64)求导得：

(65)

其在区间t∈[ti,ti+1]内前N-1个配点上的微分方程如下：

(66)

令，

(67)

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

(68)

则式(66)可写成如下形式：

(69)

(70)

将式(70)写成如下形式：

(71)

(72)

2.3 自适应Gauss-Lobatto积分计算矩阵Aj, Bλ

在确定加权矩阵W(t)后，按照式(39)和式(45)，矩阵Aj,Bλ可依靠时域[t0,tf]内的模型积分[19-20]而获得。为使本文所设计制导方法适应于将来在线规划应用，出于平衡计算精度和计算效率的目的，本文采用自适应Gauss-Lobatto积分进行矩阵Aj,Bλ计算。在积分精度不满足要求时，增加分段个数，而每个分段区间内配点数不变来保证积分精度。

首先将时域区间[t0,tf]等分为Ns个时间区间，即t0

(73)

则谱灵敏度矩阵Aj和矩阵Bλ由各个时域区间内的Gauss-Lobatto积分值累加得到，即:

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

(74)

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

(75)

式中：hk为积分权重，其计算如下：

(76)

自适应积分法就是按照积分区间上被积函数变化的剧烈程度，动态调整积分节点的分布密度，以达到节约计算成本的目的。综合积分精度和计算效率两者考虑，本文采用以下区间数划分准则：

记区间[ti,ti+1]长度为h，区间[ti,ti+1]上采用N个LGL配点进行Gauss-Lobatto积分的结果为R1(ti,ti+1)，将区间[ti,ti+1]对分，在区间[ti,ti+h/2], [ti+h/2,ti+1]上分别采用同样N个LGL配点的复化Gauss-Lobatto积分结果为R2(ti,ti+1)。设整个区间[t0,tf]上要求的积分精度为ε，若:

(77)

则区间[ti,ti+1]不再对分，否则将区间[ti,ti+1]对分，并令ti+1=ti+h/2，这时时域[t0,tf]上的区间个数Ns增加为Ns+1个，并重复上述过程。

3 基于GS-MPSP的中制导律设计

首先为了保证所有状态变量具有类似的数值变化范围，将拦截弹状态变量进行归一化，并将拦截弹动力学模型式(22)和式(23)进行归一化[20]。

利用GS-MPSP算法进行制导律设计时，需要对初始控制量进行猜测，本文运用比例导引法(PN)给出控制量初始解。文献[20]给出了PN制导俯仰、偏航方向的指令加速度如下：

(78)

(79)

如图3所示，本文基于GS-MPSP规划的中制导实现流程描述如下：

图3 基于GS-MPSP规划的中制导实现流程Fig.3 Process of mid-guidance using GS-MPSP

1)假设给定中制导飞行时间tf，并由预警探测系统给出tf时间后的目标预报状态RT,VT，结合式(13)至式(21)，计算终端点期望状态YNf。

2)采用PN制导律式(78)和式(79)计算控制初始解Ul。

3)采用四阶龙格库塔进行拦截弹弹道积分，计算终端点状态偏差，若满足精度则输出控制序列Ul，否则进入下一步，这一步也对应着GS-MPSP方法的模型预测模块。

4)结合式(59)、式(68)、式(70)和式(72)计算矩阵W(t)。

5)结合式(74)至式(77)计算矩阵Aj,Bλ。

6)结合式(57)进行控制更新Ul+1(t)，令Ul(t)=Ul+1(t)并返回第3)步，直到满足精度结束，这一步也对应着GS-MPSP方法的静态规划模块。

4 仿真结果与分析

为验证本文所设计方法的有效性，开展以下两种情形进行仿真。设定仿真中拦截弹导引头最大探测距离为ρmax=50 km，并假设中制导过程中目标在交接班时刻的状态信息依靠远程预警系统或拦截弹的火控雷达预测得到。仿真平台相关参数为：Intel Core i5 3.00 GHz处理器，window10操作系统，matlab18b作仿真软件。

4.1 情形一

情形一设定合理中制导飞行时间为76 s，考虑同一中末交班目标与拦截弹速度比δ=1.31以及侧窗探测视场范围3°～33°，自动驾驶仪一阶时间常数取0.3， GS-MPSP仿真中取谱函数个数Np=4，时域[t0,tf]内等分区间初始个数为Ns=2，各分段时间区间内选择的LGL配点数为N=4，其他仿真场景参数以及利用式(13)至式(21)计算所得的终端期望状态结果如表1所示，仿真结果如表2和图4至图8所示，其中表2中Δrf表示终端点位置偏差，Δθmf表示终端弹道倾角偏差，Δφvmf表示终端弹道偏角偏差，ξf表示交班点视窗角。

情形一采用GS-MPSP、MPSP[17]、CVX轨迹规划[14]以及PN制导[20]四种方法进行了对比分析。其中前三种方法对终端位置[xf,yf,zf]的收敛阈值均设为0.5 m，对终端弹道倾角或弹道偏角偏差收敛阈值均设为0.5°，并设置迭代步数不超过3步，CVX轨迹规划同MPSP皆采用欧拉法对动力学方程进行离散。

表1 仿真场景参数设置Table 1 Parameter setting of simulation scene

表2 情形一仿真结果统计Table 2 Case 1 Statistics of simulation results

如图5、图6和表2所示，PN制导虽能成功到达期望终端点，但未能满足终端速度方向约束。以PN制导结果作为猜测初值用于前三种方法仿真，本文提出的GS-MPSP方法同MPSP方法和CVX轨迹规划法在准确到达期望终端点的同时，满足了终端速度方向约束，且中末交班时刻视窗角在侧窗探测视场中央附近，给末制导段提供了较好的拦截态势，其中GS-MPSP制导精度相对较高。

从图4至图6可以看出，GS-MPSP同MPSP产生的轨迹及指令加速度基本一致，这也间接说明GS-MPSP方法原理上与MPSP方法之间的等效性。但也存在着轻微的不同，这是由于两者在性能指标的构建上存在不同，本文设计的GS-MPSP方法以谱函数系数的加权和最小为指标，而MPSP方法以控制量的加权和最小为指标，在取加权矩阵都是单位阵的情况下，前者在制导中期有着相对较小的控制成本，相比之下，利用CVX轨迹规划方法的控制成本明显增大，且轨迹波动相对较大。

图4 纵向和侧向指令加速度曲线Fig.4 Longitudinal and lateral command acceleration change

图5 侧向和纵向平面内二维轨迹Fig.5 Two-dimensional trajectory in lateral and longitudinal planes

图6 弹道倾角和偏角变化曲线Fig.6 Ballistic inclination and ballistic deflection change

如前几节所述，谱灵敏度矩阵Aj和过程矩阵Bλ的计算是GS-MPSP方法的重要组成部分。图7和图8展示了第一次迭代中采用自适应正交配点法计算矩阵Aj,Bλ时，其中各元素的计算结果，横坐标为元素序号，纵坐标为元素值，其中矩阵Aj个数为4，矩阵Bλ个数为1，矩阵Aj为5行2列，矩阵Bλ为5行1列，因此矩阵Aj中元素个数为10，矩阵Bλ中元素个数为5。

将采用四阶龙格库塔法解算式(32)、式(33)、式(39)和式(45)所得的矩阵Aj,Bλ中各元素值作为标准值，并与本文方法进行对比可看出，两者相对误差较小，这表明利用自适应正交配点法能保证动力学方程较高积分精度。由于GS-MPSP方法的收敛性和收敛精度与谱灵敏度矩阵的计算精度相关，从终端状态偏差结果来看，基于自适应正交配点的GS-MPSP方法相比采用欧拉积分法对动力学方程进行离散的MPSP方法和CVX规划方法具有更高制导精度。

图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规划方法的二十分之一，这是由于本文采用自适应正交配点法，在保证动力学方程积分精度的同时较大地减少了离散点个数，从而较大地提高了计算效率，对于将来拦截弹在线规划等方面应用具有较大优势。

综合情形一仿真结果分析，本文基于自适应正交配点的GS-MPSP方法在满足制导精度的同时控制消耗最小，且计算效率最高。

表3 计算耗时统计Table 3 Calculation time consumption statistics

4.2 情形二

本节在4.1节仿真条件的基础上，考虑拦截弹初始状态扰动、气动阻力系数摄动及质量偏差，偏差项设置见表4。利用蒙特卡洛打靶仿真200次，统计结果见表5和图9，其中表5中各终端参数物理意义同表2。由于考虑到本文所设计制导律在实际应用时的实时性问题，出于对收敛精度和计算效率的衡量，同样将GS-MPSP迭代步数限制在3步以内，以满足实时性要求，进而评估在保证实时性前提下的制导律的鲁棒性能。

表4 仿真偏差设定Table 4 Simulation deviation setting

图9 终端状态偏差分布图Fig.9 Distribution of terminal state deviation

表5 打靶结果统计Table 5 Shooting result statistics

可以看出，在考虑初始状态扰动和参数摄动的情况下，到达期望终端约束的偏差在合理范围内，并保证中末交班时视窗角在视场中央，表明本文设计的GS-MPSP制导方法具有一定鲁棒性。

在考虑其他分系统理想工作、拦截末制导段拦截时间较短及弹目速度变化较小情况下，情形一和情形二验证了本文考虑交班视窗角约束的中制导算法的有效性。在后续研究中，还需综合考虑目标强机动、目标预报误差及拦截弹控制系统动态响应特性等对制导律进行进一步设计。

5 结 论

本文为解决临近空间防御中考虑零控交班及交班视窗角约束的中制导问题提供了一种制导精度和计算效率较高的中制导方法，具体总结如下：

1)推导分析了拦截末制导段零控拦截条件，此条件由中末交班时刻目标和拦截弹速度比及两者速度分别与视线夹角确定，并利用此条件提出了一种将中末交班侧窗视窗角约束同零控脱靶量约束转化为中制导终端点位置、速度方向约束的方法，为满足多约束的中制导设计提供了更简便的思路。

2)提出了一种固定时间条件下的广义拟谱模型预测静态规划算法。通过引入低维权重矩阵和控制量的谱表达式以降低静态规划问题的维数，并推导得到了控制更新显示表达式，该表达式直接由依靠起始到终端时刻连续积分的谱灵敏度矩阵、过程矩阵和终端偏差确定，而拉格朗日乘子和谱系数仅作为迭代过程中的中间量，解决了拉格朗日乘子难求解，以及谱系数初值难确定的问题。相比PN制导、MPSP方法及CVX轨迹规划方法，GS-MPSP方法控制成本较小且弹道更平滑，制导精度更高。

3)提出了基于自适应正交配点法的谱灵敏度矩阵和过程矩阵的计算方法。利用自适应Gauss-Lobatto积分法实现时间分区间个数自适应分配来计算谱灵敏度矩阵和过程矩阵，从而保证计算精度；用Legendre伪谱法实现了权重矩阵的解算。相比CVX轨迹规划方法及MPSP方法，有着更高的求解效率，对将来实时在线规划应用有着较强适用性。

4)在考虑拦截弹初始状态扰动和参数摄动的条件下，本文方法在保证实时性要求情况下，依然能够满足中末交班制导要求，具有较强鲁棒性。