孙兴龙，马克茂，姜宇，侯振乾
（1.哈尔滨工业大学航天学院，黑龙江 哈尔滨 150001；
2.上海机电工程研究所，上海 201109）
近年来，由于临近空间飞行器等高速武器的迅猛发展，相应拦截技术和防御手段的研究已迫在眉睫。针对此类高空、高速机动目标，拦截弹大多采用复合制导策略，即“中制导+末制导”的远程、梯次的拦截方式。在复合制导的过程中，对中、末制导交班时刻提出了一定条件，首先是导引头交班，即导引头可靠的截获目标；
其次是弹道交班，即弹道应平滑过渡。以上分别对中制导结束时刻的弹-目相对距离和弹道前置角进行了约束。
目前常用的制导策略包括比例制导、滑模制导和最优制导等。其中，最优制导的求解包含了对约束条件的考量，因此适用于具有终端约束的中制导律设计。文献［1］根据中制导末端的视场角约束，提出一种新的虚拟目标的设置方法，并运用最优控制理论设计了满足视场角约束的中制导律。文献［2-3］对中制导问题进行了描述，在满足相应约束条件下，根据不同的优化指标采用GPM对最优中制导问题进行了离散求解。其中，前者寻求飞行时间最短，后者寻求对制导角度的约束。周觐等［4］采用GPM对含有相关约束条件的中制导律进行了求解并将此作为基准最优弹道，通过改变终端约束条件，并根据邻域最优控制理论对基准最优弹道进行调整，最终得到中制导最优弹道簇，有效减少了弹道优化时间。在最优中制导律求解的相关成果中，多数旨在提高弹道优化的求解效率。然而，在中制导阶段，所需制导信息多数依赖于载机或地面基站来提供，因此，可能存在由于通讯故障等原因导致的制导指令无法更新的风险。
滑模控制因其对外部干扰具有强鲁棒性且算法简单的优点，也被广泛应用于制导律的设计［5-7］。对于高速目标而言，为保证制导回路性能不受目标机动的影响，特别是末制导阶段，需在设计制导律时加入对目标机动的补偿项［8］。文献［9］利用扩张状态观测器对目标机动等干扰项进行估计并补偿，设计了带角度约束的滑模制导律。然而，末制导阶段对导弹的响应速度及拦截精度有较高要求，因此，形式简单及制导精度较高的比例导引常常被用于末制导律的设计。文献［10］根据精确的虚拟目标比例导引模型设计了比例导引律，在保证导弹以给定落角击中目标的同时，避免了导引末端控制量发散的问题。文献［11］提出了一种带落角约束的偏置比例导引律，并设计了一种盲区控制方案以减少终点处的法向过载。
基于以上分析，本文基于GPM离散化方法及滚动时域的方式，对含有终端约束的中制导律进行优化求解，以此实现拦截弹道的在线修正。采用高增益观测器对目标机动等信息进行估计，并引入末制导律的设计中，实现对临近空间飞行器等机动目标的有效拦截。
在三维空间中，弹-目相对运动几何关系的描述如图1所示。图中：M和T分别表示导弹和目标的质心；
R表示弹-目之间的相对距离；
qε和qβ分别表示视线高低角和偏转角；
θm和ψm分别表示导弹弹道倾角和弹道偏角；
Vm和Vt分别表示导弹和目标的速度大小。坐标系Mxyz平行于参考惯性坐标系，Mx4y4z4为视线坐标系，Mx4轴以弹-目视线方向为正，My4轴位于包含视线的铅锤面内，与Mx4轴垂直且向上为正，Mz4轴与Mx4轴、My4轴满足右手定则。
图1 三维空间导弹-目标相对运动几何关系Fig.1 Geometric relationship of relative motion between missile-target in 3-D space
在以上空间关系中，弹-目相对运动方程可描述为
式中：(atr，atε，atβ)、(am4r，am4ε，am4β)分别表示目标加速度和导弹加速度在视线坐标系的3个轴上的分量。
在对临近空间高速目标的拦截过程中，中制导律的设计通常以能量消耗最小或飞行时间最短为前提，目标是将导弹以较好姿态导引至相应空域并为末制导阶段的成功拦截创造有利条件。因此，首先对中制导律求解过程的相关约束进行描述。
2.1 约束条件
设中制导过程的初始时刻为0，终端时刻为tf。初始时刻制导系统的状态变量为
1）交班距离约束
依据导引头的探测距离，设定中、末制导交班时刻弹-目相对距离为Rc，则可设定
2）弹道前置角约束
为实现中、末制导的顺利交班，在末制导终端时刻要保证拦截弹前置角小于30°。由图1所示，拦截弹前置角在俯仰平面为弹道倾角θm与视线高低角qε的差值，在偏航平面为弹道偏角ψm与视线偏转角qβ的差值。在此将该约束转换为视线角的终端约束，设qε和qβ的期望值分别为qεd和qβd，则有
式 中：θm-30° 3）视线角速率约束 在制导过程中令视线角速率趋于零，可保证导弹以准平行接近的方式靠近目标。于是令 4）拦截过载约束 由于导弹物理条件所限，拦截过载满足 式中：um为拦截过载的上限值。 以拦截所需能量为性能指标 式中：tgo=tf-t表示剩余飞行时间； 最后，最优中制导律的设计问题归结为：确定uε、uβ和终端时间tf，在满足式（1）及以上约束条件的前提下，使得式（6）最小。 针对类似以上优化问题，传统求解过程通常要构造相应Hamiltonian函数，这样势必引入协态变量并增加了运算量，因而不利于高速目标的拦截。现采用GPM对以上优化问题进行离散化处理，并转化为NLP的形式进行求解。 在进行中制导律求解之前，先对GPM离散化求解过程进行简单描述。 2.3.1 非线性最优反馈控制 考虑一般形式的非线性系统模型 式 中 ：x(t)∈X⊂RNx，X为 状 态 空 间 ； 针对以上非线性系统，最优反馈控制求解的实质为下述Bolza类型代价函数的优化问题： 式中：E(⋅)表示终端性能指标； 2.3.2 基于GPM的连续Bolza问题的离散求解 GPM通过离散化处理，将连续Bolza问题转化为NLP，通过求解NLP来确定原始最优问题的解。由于GPM的配点都分布在区间[-1，1]上，因此伪谱法的第一步都是将上一节中最优控制问题的时间区间由t∈[t0，tf]转换至τ∈[-1，1]，对时间变量t作如下变换： 上节所述最优控制问题转化为 GPM选取K个Legendre-Gauss（LG）点以及τ0=-1为节点，构成K+1个Lagrange插值多项式Li(τ)(i=0，…，K)，并以此为基函数近似状态变量： 式中：Lagrange插值基函数为 使得节点上的近似状态与实际状态相等。 同 理，采 用Lagrange插 值 多 项 式Lˉi(τ)(i=1，…，k)作为基函数近似控制变量： 式中：τi(i=1，…，K)为LG点。 对近似状态变量求导可得状态变量导数，从而将动力学微分方程约束转换为代数约束，即 式中：微分矩阵D∈Rk×(k+1)可离散确定，即 式中：τk(k=1，…，K)为集合κ中的点； 这样，就将最优控制问题的动力学微分方程约束转换为代数约束 GPM中的节点包括K个配点、初始点τ0=-1以及终点τf=1。上文中近似状态变量表达式未定义终端状态Xf，终端状态也应满足动力学方程约束： 将终端状态约束条件离散并用Gauss积分近似，得 将性能指标函数中的积分项用Gauss积分来近似，得 边界条件 路径约束 至此，以上微分方程代数约束、离散性能指标、离散边界条件和离散路径约束就定义了一个NLP。即在满足式（14）、式（17）以及式（20）~（21）的情况下，通过求解式（11）、式（13），保证式（19）最小。该NLP的解即为连续Bolza问题的解。 2.3.3 实时最优反馈控制算法 序列二次规划算法多被用来求解约束优化问题，采用序列二次规划算法对以上NLP进行求解。考虑到在中制导过程中，可能会出现载机与拦截弹通讯故障等情况，因制导信息缺失不能及时更新制导指令。为应对此类突发状况，基于所求得的NLP开环最优解序列，并考虑实际应用过程中可能存在的信息中断问题，设计了一种实时最优闭环反馈控制算法。具体步骤如下： 步骤1选定节点个数N，采样周期为T，即区间[t0，T1]的间隔为定值T，初始状态为x0，离线计算得出开环最优控制量(x(t0)，t)。 步骤2在区间[t0，T1+Δt1]内，将最优控制量(x(t0)，t)应用于系统上，记T1时刻的状态测量值为x(T1)，令i=1； 步骤3在区间[T1+Δt1，T2+Δt2]内，于t1时刻发生通讯中断导致数据丢失，且通讯中断时常未知，此时将以状态x(T1)求得的开环最优控制量(x(T1)，t)进行t∈[T1，tf]时间域内插值，并将插值所得控制量取出应用于系统上，直至t2时刻通讯恢复并重新对系统状态进行采样。 步骤4在区间[Ti+Δti，Ti+1+Δti+1]内，将开环最优控制量(x(Ti)，t)应用于系统上，记Ti+1时刻的状态测量值为x(Ti+1)； 步骤5令i=i+1，返回步骤4。 中制导律求解过程如图2所示，其中t1为数据丢失初始时刻，t2为数据恢复时刻。 图2 信息缺失情况下的制导律实现Fig.2 Realization of guidance law in the absence of information 注1：该算法将状态采样周期设为T，可将此作为控制周期，根据系统实际需要进行调节。同时，能够很好地应对制导过程中出现的信息中断问题。 注2：如在算法步骤2和步骤4中的描述所示，基于前一采样时刻所得的系统状态计算出开环最优解序列，将该最优解序列的第一值应用于制导系统，直至基于下一采样时刻的系统状态计算出新的开环最优解序列。这种不断滚动的动态优化过程，可达到基于开环最优解来实现闭环控制的目的。在图3所示的制导系统流程中也有体现，这种处理思想称为滚动时域控制。 图3 闭环制导系统流程Fig.3 Process of closed loop guidance system 至此，可保证拦截弹在中制导结束时满足交班条件，并顺利进入末制导阶段。为降低末制导阶段目标机动对拦截所需过载的要求，本文基于高增益观测器及比例导引策略，设计了带有目标机动补偿的末制导律。 在 末 制 导 阶 段，令x=[qεβ]T，y=[qεqβ]T，u=[amεamβ]T，d=[atεatβ]T， 则相对运动方程可描述为 式中： 考虑到实际物理条件限制，假设目标机动连续且幅值有限。同时，导弹也存在最大机动过载aM，即满足 假设所有变量均可获得且目标机动已知，针对式（22）所示末制导模型，将系统存在的非线性项进行抵消，并对目标机动进行补偿，可得如下末制导律： 式中：k1>0、k2>0为设计参数。将式（24）代入式（22）可得 通过适当选取设计参数k1、k2，可以保证视线转率的快速收敛，且不会受到目标机动的影响。在式（24）所示制导律中，导引头可提供变量包括R、Ṙ、qε和qβ，式中其他不可直接测量的变量则采用下节中描述的高增益观测器进行估计。 为对目标加速度等制导信息进行估计，针对式（22），建立高增益扩展状态观测器： 式中：̂观测器状态变量为观测器扩展状态变量。观测器矩阵增益H(ε)和F(ε)分别取为 式中：参数aij(i=1，2，3； 目标机动信息不可直接测量，需由观测器的相应状态变量进行替换。由于3.1节中观测器进行设计时，扩展观测变量σ̂是对a(y，t)d的估计，因此末制导律中的atε和atβ用-Rσ̂1和σ̂2Rcosqε代替，得到制导律的实现如下： 观测器设计中，由于引入了高增益反馈，在观测误差向量中会产生尖峰现象，利用饱和函数对控制量进行限幅： 式中：γ≥0为设计中预留的过载裕量； 在末制导律的实现过程中，通过引入饱和函数保证了控制的有界性。由于尖峰现象持续时间较短，并不影响系统模型的变化，因此考虑时间区间[T∗(ε)，T]内，假设控制量处于线性区内，将控制量u1、u2代入相对运动模型中，可得闭环系统如下： 式（31）可看作对式（25）所示标称系统的摄动，其中摄动项满足 式中：Lε>0和Lβ>0为常数。对于观测误差收敛性，文献［11］中已经给出，在此不再累述。此处，考虑制导回路的闭环稳定性时，只考虑式（31）所示系统。 定义Lyapunov函数 经计算可知，V()沿式（31）所示闭环系统对时间的导数满足 由此可得，式（31）所示闭环系统的状态轨迹将趋向于集合Ωq。 以上中、末制导律的完整闭环制导过程如图3所示。 以再入机动目标为例，拦截弹与目标的初始条件见表1。假设弹-目初始相对距离R(0)=350 km，目标纵向、侧向机动加速度为atε=3gsin(πt/20)、atε=5gcos(πt/20)，中、末交班时刻弹-目相对距离Rc=30 km，中制导控制周期为T=1 s，制导指令上限um=7g，性能指标权重系数n=1，节点个数N=25。同时假设中制导开始5 s后，出现持续时长为5 s的通讯中断，以模拟制导信息缺失情况下制导算法的有效性，进行拦截仿真验证。本文中制导指令序列插值算法采用三次样条插值，仿真结果如图4~5所示。 表1 拦截弹与目标初始条件Tab.1 Initial conditions of interceptor and target 图4 拦截弹中制导指令Fig.4 Interceptor midcourse guidance command 图4为拦截弹中制导指令，图5为弹道前置角变化。由图4可以看出，本文所求拦截弹中制导指令为锯齿状。这是因为中制导弹道优化周期采用可调定值T，即每隔一个时间T，对系统状态进行一次采样，并以采样时刻状态为初值进行一次优化求解。此举避免了弹载制导系统的计算负荷，并保证了中制导阶段的弹道约束。 图5 拦截弹弹道前置角曲线Fig.5 Heading angle curve of interceptor trajectory 由图5可以看出整个中制导阶段的弹道前置角均维持在20°以内，满足中、末交班弹道前置角约束。从图4~5中可以看出，拦截开始后的第5 s时刻，本文所设计的中制导优化算法可以很好地应对因通讯中断等原因造成的制导信息缺失状况。而通常情况下，当制导信息无法获得时，制导系统会以故障发生前的指令来导引拦截弹，此时会令弹道产生偏差甚至影响拦截效果，如图5中蓝色曲线所示。 为了验证所设计的最优中制导律在满足相应约束的前提下，可以保证以所需过载的积分为性能指标的最小化，现以比例导引为对照组进行中制导阶段的仿真验证。选取式（6）中n=0，比例导引的俯仰通道比例系数为4、偏航通道比例系数为3，其他初始条件同上文，验证结果如图6~7所示。图6为拦截弹所需控制量的对比曲线，图7为中制导阶段性能指标的对比曲线。 由图6可以看出，在拦截过程的前期阶段，比例导引所需控制量远高于最优中制导，这是由比例导引的特性决定的。反观最优中制导，在前期阶段所需控制量不算很大，有效降低了不必要的能量浪费。图7所示的性能指标变化曲线更好地反映了这一点。由图7可以看出，与比例导引相比，最优中制导的性能指标更小。在中制导结束时刻，比例导引的性能指标大小为73 022 m2/s3，而最优中制导的性能指标大小仅为40 692 m2/s3，该结果显示最优中制导具有式（6）所示指标的最优性。 图6 最优中制导与比例制导控制量对比Fig.6 Comparison of control quantities between optimal midcourse guidance and proportional guidance 图7 最优中制导与比例制导性能指标对比Fig.7 Comparison of performance indexes between optimal midcourse guidance and proportional guidance 中制导结束后，进入末制导阶段，此时高增益观测器开始对目标机动等信息进行估计。选取观测器参数为α1j=3、α2j=3、α3j=1，其中j=1、2，参数ε=0.01，观测器状态初值为0，控制器参数k1=4、k2=3。 图8为观测器对目标机动加速度的观测曲线。从图中可以看出，在观测初期有较为明显的振荡，但很快能跟踪上真实值，在临近拦截结束时，目标加速度观测误差开始变大。这是由于弹-目相对距离R趋近于零，且视线角速率急剧增大所致。 图8 目标加速度观测曲线Fig.8 Target acceleration observation curve 图9为拦截弹全程制导指令曲线，从图中可以看出，过载指令在66.1 s处出现尖峰，此为观测器中高增益反馈所致。而本文进行末制导律设计时，通过饱和函数对该尖峰现象进行了抑制。图10为三维拦截的轨迹，最终脱靶量为1.05 m，可以看出本文所设计制导律可有效拦截机动目标。 图9 拦截弹制导指令Fig.9 Interceptor guidance command 图10 三维拦截轨迹Fig.10 3-D interception trajectory 本文对最优中制导律和基于目标机动补偿的末制导律进行了设计。首先将中制导律的设计描述为多约束条件的最优求解过程，采用GPM将连续最优控制问题转化为NLP形式，在满足相应约束的前提下实现了过载指令的在线求解，并基于滚动时域的方式完成了最优中制导律的闭环实现，同时保证了在制导信息短时间缺失等工况下的适用性。考虑末制导阶段目标机动对拦截性能的影响，采用高增益观测器对视线转率和目标机动等信息进行估计，并将估计值作为反馈补偿项引入末制导律，以此降低了高机动目标拦截时对过载的要求。
2.2 性能指标
1tgno(n≥0)为时变权重系数，随tgo的减小而不断增大，n越大，权重系数的影响也越大。2.3 GPM离散求解最优控制问题
u(t，x)∈U⊂RNu，U为控制空间；
p0为系统标称参数；
x为状态向量；
u为控制向量；
f⊂RNx，关于其变量是Lipschitz连续的。
F(⋅)表示动态性能指标；
ẋ(t)=f(⋅)表示动态约束；
Φ=0表示边界约束；
C≤0表示路径约束。
τi(i=0，…，K)属于集合κ0={τ0，τ1，…，τK}。
同时，根据时刻T1的状态测量值x(T1)，实 时 计 算 开 环 最 优 控 制 量(x(T1)，t)，t∈[T1，tf]，设Δt1为优化计算时间。
同时，根据时刻Ti+1的状态测量值x(ti)，实时计算开环最优控制量(x(Ti+1)，t)，t∈[Ti+1，tf]，设Δti为优化计算时间。
3.1 制导信息估计
j=1，2)均选择为实数，且使得多项式s3+a1js2+a2js+a3j(j=1，2)均为Hurwitz多项式；
ε>0为小的设计参数。3.2 基于目标机动补偿的导引律实现
sat(⋅)为饱和函数，定义为
3.3 末制导回路稳定性分析
