采用状态转移矩阵的拦截中段制导方法*

柴 华，仲 明，梁彦刚

采用状态转移矩阵的拦截中段制导方法*

柴 华1，仲 明2，梁彦刚1

(1.国防科技大学 航天科学与工程学院， 湖南 长沙 410073；

2.中国人民解放军69006部队， 新疆 乌鲁木齐 830001)

基于相对运动理论提出拦截中段制导方法。方法原理为，将拦截器初始轨迹与施加修正后的轨迹视为一主一从两空间目标，以相对运动理论描述从目标相对于主目标的运动规律，于是初始的修正量可由终点处相对状态求解。给出了一般形式的相对运动模型，运用几何法与变分法推导得到了J2摄动影响下相对运动的状态转移矩阵，在此基础上，提出了采用状态转移矩阵的拦截中段制导方法。仿真算例表明，提出的方法能够为工程实际中的拦截中段制导提供有效支持。

中段制导；相对运动；J2摄动；状态转移矩阵

(1.CollegeofAerospaceScienceandEngineering,NationalUniversityofDefenseTechnology,Changsha410073,China;

2.ThePLAUnit69006,Wulumuqi830001,China)

在导弹防御作战中，中段制导是指拦截器在助推段结束之后、末段导引开始之前，为消除环境噪声等因素造成的飞行状态误差，采用脉冲推力等方式实施一次或多次状态修正的过程。中段制导的效果好坏直接关系到拦截器能否捕获目标并在末段导引中占据优势的拦截几何，因而具有十分重要的作用[1-3]。脉冲推力式的中段制导事实上可以抽象为一固定时间拦截问题，即给定初、终位置与飞行时间，求解初始需用速度的问题。在二体假设下，求解固定时间拦截问题有许多经典方法，如Herrick方法[4]、Godal方法[4]、Lambert飞行时间定理[5]等。然而，当引入更为精确的引力场模型时，上述方法不再适用，需要探索更为有效的求解方法。

考虑到对于拦截中段制导这一特定的固定时间拦截问题而言，飞行时间较短(约200s)、速度修正量较小(约100m/s)，本文提出基于相对运动理论的求解方法。相对运动理论广泛地应用于卫星编队飞行、空间交会对接等领域，涌现出了CW方程、TH方程等经典模型。通过引入三大假设：参考轨道偏心率为0，地球引力为中心引力场，两目标距离远小于地心距，Clohessy与Wiltshire[6]将非线性的相对运动模型简化为便于解析求解的线性时变微分方程组，即CW方程(亦称Hill[7]方程)。舍弃圆参考轨道假设，Tschauner与Hempel推导得到了适用于椭圆参考轨道的线性相对运动模型，给出了以参考轨道真近点角/偏近点角表示的解析解，其相对运动模型被称为TH方程。Yamanaka与Ankersen[8]引入变量替换，简化了TH方程的求解过程，推导了一种形式简洁的状态转移矩阵，并探讨了其与CW方程的一致性。尽管CW方程与TH方程在工程实际中得到了广泛应用，但两者均以二体中心引力场为前提，因此难以满足要求。Gim与Alfriend等[9-11]针对J2摄动影响下的相对运动模型展开了研究，避免了直接求解复杂微分方程组，推导给出了相对运动方程的状态转移矩阵，为实现高精度的相对运动建模创造了条件。本文将以此状态转移矩阵为基础，提出一种有效的固定时间拦截问题求解策略，并探讨其在拦截中段制导中的实用性。

1 坐标系定义

所涉及的坐标系为地心惯性坐标系与当地垂直/当地水平坐标系，两者定义如下：

地心惯性(EarthCenteredInertial，ECI)坐标系，记为I。 原点OE位于地心，OExI轴在赤道平面内指向平春分点，OEzI轴垂直于赤道平面指向北极，OEyI轴与OExI，OEzI轴垂直并满足右手定则，如图1所示。由于岁差与章动影响，春分点具有进动性。在建立地心惯性系时需指定春分点基准。本文采用的ECI系是以2000年1月1日12:00:00的平春分点为基准，因此又称为J2000坐标系。

当地垂直/当地水平(LocalVerticalLocalHorizontal，LVLH)坐标系，记为L。原点O位于飞行器质心，Ox轴指向飞行器地心矢径方向，Oy轴在飞行器弹道(轨道)平面内垂直于Ox轴并指向速度方向，Oz轴与Ox，Oy轴垂直并满足右手定则，如图1所示。

图1 ECI坐标系与LVLH坐标系Fig.1 ECI coordinate system and LVLH coordinate system

由ECI系到LVLH系的转换矩阵可由式(1)给出。

MI→L=M3(θ)M1(i)M3(Ω)

(1)

式中，θ=ω+f为飞行器纬度幅角(ω为近地点角距， f为真近点角)，i为轨道倾角，Ω为升交点赤经。

2 相对运动方程

设空间中存在主目标与从目标两飞行器，两者在ECI系下的位置矢量分别为Rc与Rd，由经典轨道力学可知

(2)

式中，μ为地球引力常数，fc与fd分别为两目标所受摄动力产生的加速度。

将从目标相对于主目标的位置矢量记作r=Rd-Rc，那么有

式中，f=fd-fc为两目标摄动加速度之差。

以主目标LVLH系作为基准坐标系，记作LC。式(3)在LC系中描述需要考虑绝对导数与相对导数的相互转换，即

(4)

结合式(3)、式(4)可得

(5)

式(5)即为LC系下一般形式的相对运动方程。对应于不同的引力场模型，该方程具有不同的表现形式，其解析求解是十分困难的。在二体引力场假设下，可将式(5)线性化(即TH方程)，并推导给出形式简洁的状态转移矩阵，相关研究已经较为成熟。当考虑更为精确的引力场模型时，问题变得更加复杂。Alfriend等学者对这一问题进行了研究，推导了J2摄动影响下相对运动的状态转移矩阵。

3 状态转移矩阵

为避免直接求解微分方程带来的困难，Alfriend等学者采用几何法来推导相对运动方程的状态转移矩阵。该方法的主要思路为，将从目标相对于主目标的位置速度转换为轨道根数之差，以轨道根数之差随时间的变化来刻画相对状态的迁移。在推导过程中，需引入差分轨道根数的概念。首先给出两组相互等价的轨道根数：

(6)

式中：eCL为经典轨道根数，a为半长轴，e为偏心率；eNS为非奇异轨道根数，q1=ecosω，q2=esinω。为避免经典轨道根数在描述圆轨道时出现的奇异现象，下文将采用非奇异轨道根数实现状态转移矩阵的推导。

差分轨道根数定义为从目标与主目标轨道根数之差，即

δe=ed-ec

=[δa,δθ,δi,δq1,δq2,δΩ]T

(7)

式中，下标d与c分别表示从目标与主目标。需要指出，在下文推导中，所有不含下标c的轨道根数(或其衍生项)均对应于主目标。例如，式(8)中的R等价于式(2)中的Rc。

3.1 由相对位置速度到差分密切轨道根数的转换

以LC系作为基准坐标系，主目标地心矢径与速度矢量可分别表示为

(8)

从目标地心矢径与速度矢量可分别表示为

(9)

式中：x，y，z分别为从目标相对位置矢量r在LC系三方向的分量；ϖx，ϖy与ϖz分别为主目标角速度矢量在LC系三方向的分量，即

ϖc=[ϖx,ϖy,ϖz]T

(10)

与此同时，从目标满足式(11)所示关系。

(11)

式(1)已经给出了ECI系与LVLH系之间的转换关系，考虑到主目标与从目标的轨道根数足够接近，可将矩阵MLD→I在MLC→I处泰勒展开且仅保留一阶项，结合式(8)、式(9)有

(12)

由轨道力学可知

(13)

(14)

(15)

将式(13)～(15)代入式(12)，可将LC系下从目标的相对状态表示为差分轨道根数的形式。于是有

(16)

3.2 由差分密切轨道根数到差分平均轨道根数的转换

在J2摄动影响下，密切轨道根数随时间的变化较为复杂，难以给出解析描述。因此，在推导轨道根数之差的时间迁移规律时，采用平均轨道根数作为替代对象。本节将推导给出差分密切轨道根数与差分平均轨道根数之间的转换关系。

许多文献对密切轨道根数与平均轨道根数的转换问题进行了研究，此处采用文献[10]的方法，即有

(17)

式中，e为密切轨道根数，em为平均轨道根数，elp为长周期项，esp1与esp2为短周期项。需要指出的是，elp，esp1与esp2均以平均轨道根数表示，也就是说，由密切轨道根数计算平均轨道根数需通过迭代来实现。下文将在4.1节对这一问题进行详细阐述。

由式(17)可知

(18)

式中，I为单位矩阵。

于是有

δe(t)=D(t)δem(t)

(19)

3.3 差分平均轨道根数的时间迁移

与密切轨道根数相比，平均轨道根数的变化规律更为简洁，有

(20)

式中，上标m表示平均根数，下标0表示初始时刻，M为平近点角，λ=M+ω为平纬度幅角。ωm，Mm与Ωm的时间变化率均可由初始时刻的平均轨道根数表示。

对式(20)中最后一式等号两边进行变分操作，保留一阶项，有

(21)

式中，Δt=t-t0。

δem(t)=Φe(t,t0)δem(t0)

(22)

结合式(16)、式(19)、式(22)，可以得到J2假设下相对运动方程的状态转移矩阵，即

(23)

4 拦截中段制导方法

4.1 平均轨道根数的计算

式(17)给出了密切轨道根数与平均轨道根数之间的关系，依据该式计算平均轨道根数需要通过迭代来实现。迭代步骤如下：

(24)

4) 重复上述过程，直到满足收敛准则为止。收敛准则可取为

(25)

式中，ε为容许误差限。

4.2 制导方法流程

前文指出，拦截中段制导问题事实上可以抽象为固定时间拦截问题。考虑到对于中段制导而言，飞行时间较短而速度修正量较小，因而本文提出基于状态转移矩阵的制导方法。

问题描述：

在ECI坐标系下，已知t0时刻拦截器的位置、速度矢量分别为R0，V0，预测命中时刻为tf，预测命中点为Rf。 由于误差因素的影响，拦截器以初始状态难以准确到达预测命中点，需要实施中段导引，确定速度修正量Δv。

方法流程：

1)由初始状态R0、V0计算初始时刻的密切轨道根数es，建立初始的LVLH坐标系；

4) 计算终点LVLH系下预测命中点的相对位置rf，即

(26)

6)初始时刻的速度修正量Δv满足

(27)

由式(27)可知

(28)

图2 制导方法流程

Fig.2Procedureofthemidcourseguidancemethod

5 仿真算例

上文提出了采用状态转移矩阵的拦截中段制导方法。本节将引入仿真算例，验证所提方法的性能，分析其适用条件。

5.1 状态转移矩阵的精度

首先考察本文给出的状态转移矩阵在刻画拦截器相对运动时的精度。设t0=1003.6s时刻拦截器在J2000系下的位置速度矢量分别为

R0=[2 196 308.118, -2 361 659.419, 6 082 688.107]T；V0=[853.466 678，-6 083.398 166，2 767.032 528]T。

在t0时刻向拦截器施加速度冲量，拦截器的飞行轨迹将发生变化。表1给出了LVLH系下3组不同的速度冲量值。以解析的J2轨迹外推结果作为真值，将偏离真值的程度作为衡量状态转移矩阵精度的指标。图3给出了不同速度冲量下状态转移矩阵的误差随时间的变化情况。

表1 LVLH系下的初始速度冲量

图3 状态转移矩阵的误差随时间的变化Fig.3 Error of state transition matrix versus time

由图3可知，在本文给出的仿真条件下(时间不超过200s、速度修正量不大于100m/s)，采用状态转移矩阵刻画拦截器的飞行轨迹具有较高的精度：与真值相比，位置误差不超过300m，速度误差不超过0.4m/s。需要注意的是，在图3中，随着速度冲量的增加，状态转移矩阵的性能出现下降。这是由于当速度冲量增加时，新轨道与参考轨道的偏差在增大，推导状态转移矩阵时引入的小偏差假设受到挑战，忽略非线性项的代价越来越明显地体现出来。

5.2 中段制导方法的性能

设预测命中时刻tf=1186.9s，预测命中点在J2000系下的位置矢量为

Rf=[2 322 135.847, -3 439 448.051, 6 478 344.336]T。

采用两种方法求解初始时刻所需的速度修正量Δv。 一种是图2给出的中段制导方法，一种是文献[5]给出的基于Lambert飞行时间定理的联合方程方法。将两种方法求解得到的速度修正量施加至拦截器，得到预测命中时刻拦截器与Rf点的距离Δr作为衡量制导方法精度的指标。仿真结果如表2所示，表2同时还给出了两种方法的计算耗时Δt。

表2 制导方法性能对比

由表2可知，由于考虑了J2项摄动的影响，本文提出的中段制导方法比联合方程方法更加精确。本文方法的计算耗时约为20ms，尽管比联合方程方法高出一个量级，但是仍处于可接受的范围内。

5.3 中段制导方法的适用性

需要指出，本文采用线性化的相对运动模型来实现拦截中段制导，这决定了制导方法的性能与Δv的大小密切相关(即主从目标是否足够接近从而满足线性化条件)。改变预测命中点Rf的位置，图4给出了需用速度修正量Δv与方法误差Δr的对应关系，以此考察制导方法的适用性。由图可知，随着Δv的增加，方法误差Δr也不断增大。当Δv达到100m/s时，制导方法产生的Δr不超过200m，考虑到拦截器末制导段的修正能力，这一误差量级是足够精确的。

图4 Δr与Δv的对应关系Fig.4 The relations between Δr and Δv

6 结论

对拦截中段制导问题展开了研究。针对高精度引力场模型条件下固定时间拦截问题难以求解的现状，提出了基于相对运动模型的求解策略。在J2引力场假设下，推导得到了相对运动的状态转移矩阵，从而可将固定时间拦截问题转化为相对状态的转移问题。仿真算例表明，提出的方法为快速求解速度修正量提供了一条有效途径。

由于推导过程引入了线性化，因此本文的制导方法仅适用于拦截器实际飞行轨迹与标称飞行轨迹足够接近的情形。但仿真算例分析表明，这一约束并不影响本文方法在拦截中段制导中的应用。

References)

[1]ZarchanP.Tacticalandstrategicmissileguidance[M]. 4thed.USA:AmericanInstituteofAeronauticsandAstronautics, 2002.

[2]ZarchanP.Midcourseguidancestrategiesforexoatmosphericintercept[R].CharlesStarkDraperLaboratory, 1998.

[3]NewmanB.Strategicinterceptmidcourseguidanceusingmodifiedzeroeffortmisssteering[J].JournalofGuidance,Control,andDynamics, 1996, 19(1): 107-112.

[4] 任萱. 人造地球卫星轨道力学 [M]. 长沙: 国防科技大学出版社, 1988.RENXuan.Orbitalmechanicsofartificialearthsatellite[M].Changsha:NationalUniversityofDefenseTechnologyPress, 1988.(inChinese)

[5]BattinRH.Anintroductiontothemathematicsandmethodsofastrodynamics[M].USA:AmericanInstituteofAeronauticsandAstronautics, 1999.

[6]ClohessyWH,WiltshireRS.Terminalguidancesystemforsatelliterendezvous[J].JournaloftheAerospaceScience, 1960, 27(9): 653-658.

[7]HillGW.Researchesinthelunartheory[J].AmericanJournalofMathematics, 1878, 1(3): 245-260.

[8]YamanakaK,AnkersenF.Newstatetransitionmatrixforrelativemotiononanarbitraryellipticalorbit[J].JournalofGuidance,Control,andDynamics, 2002, 25(1): 60-66.

[9]AlfriendKT,SchaubH,GimDW.Gravitationalperturbations,nonlinearityandcircularorbitassumptioneffectsonformationflyingcontrolstrategies[J].AdvancesintheAstronauticalSciences, 2000, 104: 139-158.

[10]GimDW,AlfriendKT.Statetransitionmatrixofrelativemotionfortheperturbednoncircularreferenceorbit[J].JournalofGuidance,Control,andDynamics, 2003, 26(6): 956-971.

[11]AlfriendKT,VadaliSR,GurfilP,etal.Spacecraftformationflying:dynamics,controlandnavigation[M].USA:Butterworth-Heinemann, 2009.

Midcourse guidance of interception using state transition matrix

CHAI Hua1, ZHONG Ming2, LIANG Yangang1

Onthebasisofrelativemotiontheory,amidcourseguidancemethodofexoatmosphericinterceptionwasdeveloped.Themechanismofthismethodisthattreatingtheinitialandmodifiedtrajectoriesaschiefanddeputyspaceobjectsrespectively,modelingthedynamicsofthedeputyrelativetothechiefbytherelativemotiontheoryandsolvingtheinitialcorrectionsthroughthefinalrelativestates.TheuniversalrelativemotionmodelwasprovidedandthenthestatetransitionmatrixofrelativemotionundertheinfluenceofJ2disturbancewasderivedbyemployingthegeometryandthevariationmethod.Lastlythemidcourseguidancemethodusingstatetransitionmatrixwasproposedonthisbasis.Simulationexampleshowsthattheproposedmethodservesasanefficientsupportformidcourseguidanceofinterceptorsinpractice.

midcourseguidance;relativemotion;J2disturbance;statetransitionmatrix

2015-06-10

国家自然科学基金资助项目(11272346)作者简介：柴华(1988—)，男，山西沁水人，博士研究生，E-mail:chaihua@nudt.edu.cn；梁彦刚(通信作者)，男，副教授，博士，E-mail：liangyg@nudt.edu.cn

10.11887/j.cn.201504023

http://journal.nudt.edu.cn

文献标志码： 文章编号：1001-2486(2015)04-137-06