一种基于CADET的拦截概率快速计算方法

宋 豹，于剑桥，杨 迪，郭斐然

(北京理工大学宇航学院，北京 100081)

1 引言

拦截概率是武器系统拦截效能的重要评价标准，也是制定多弹协同拦截策略的重要依据。在多智能体攻防对抗场景下，只有准确地得到单拦截器对单目标的拦截概率，才能在此基础上为拦截器合理地选择拦截对象，充分发挥拦截器效能，保证拦截方整体的拦截效果达到最优[1]。在攻防对抗过程中，拦截策略的制定一般是动态的，这对拦截概率的计算速度提出了较高的要求[2]。

对于导弹这类非线性时变系统，一般情况下采用蒙特卡罗法(Monte-Carlo)进行精度分析和概率计算[3，4]。蒙特卡洛法是一种基于随机数的方法，为了提高仿真精度需要进行大量的模拟计算，时间成本非常高，不适合在战场上即时运算。

协方差分析描述函数技术(Covariance Analysis Describing Function Technique，CADET)是美国ASC公司提出的一种快速分析非线性制导系统制导精度的方法。自提出以来，该方法得到了广泛应用。Wang X S等[5]利用CADET对初始扰动下的侧向短周期运动进行了仿真分析。蒋瑞民等[6]提出了一种分析存在内部参数摄动的导弹姿态控制系统的方法。张远等[7]利用CADET进行落点预报，提供落点的误差散布。Zhang T T等[8]基于CADET建立了多种影响因素下的视轴角速度误差传播模型。Guo J G等[9]针对一般的高超声速飞行器姿态跟踪系统进行了控制精度分析。以上研究结果表明，CADET运算速度及准确度较高，但所分析的物理系统状态量较少，结构较为简单，在非线性强、耦合程度高的复杂系统上的应用效果仍需进一步验证。

本文首次提出将CADET应用于动能拦截器拦截概率的计算这一复杂非线性问题之中，在保证计算精确度的前提下提高计算效率。首先介绍CADET的基本原理，基于泰勒级数展开提出高维随机变量函数的统计线性化近似求解方案，据此快速获得动能拦截器的三维落点均值矩阵和协方差矩阵；其次跟据这两个矩阵建立拦截概率计算公式，利用分部积分法降低了积分维数；最后将该方法应用于动能拦截器六自由度制导控制系统。本文提出的方法在保证计算准确度的前提下可以大幅缩短拦截概率的计算时间，为多弹协同拦截策略的制定提供有效依据。

2 拦截概率快速计算方法

2.1 CADET原理

一般连续时域非线性系统可以表示为

x(t)=f(x，t)+G(t)w(t)

(1)

其中，f(x，t)为n维状态向量x(t)的非线性矢量函数，x(t)是系统的n维状态向量，w(t)为p维随机输入向量，G(t)为确定性的函数矩阵。

随机状态向量x(t)包括确定性分量m(t)和随机分量r(t)

(2)

其中p(t)为状态向量的协方差矩阵。

同样，随机输入向量w(t)也由均值向量b(t)和u(t)组成

(3)

其中，u(t)是具有谱密度矩阵Q(t)的高斯白噪声。δ(t-τ)为狄拉克函数。

对于非线性系统应用CADET，在利用描述函数理论将非线性系统统计线性化后，再结合协方差分析的方法进行统计性能分析。

对于非线性矢量函数f(x，t)，其拟线性化形式可以表示为+Nr。利用描述函数理论求取合适的N值使得这个近似表达式与f(x，t)的均方差E[eTSe]达到最小，则产生最好近似。其中e=f--Nr，S为任意半正定矩阵。

若假设x(t)呈联合正态分布，描述函数和准线性增益矩阵N可以用下式计算

(4)

其中，g(x)为x(t)的概率密度函数。

在将非线性系统统计线性化后，便可以得到系统状态均值和协方差矩阵的传播方程

(5)

式中，只要m(0)和p(0)已知便可通过一次积分得到任意时刻的状态均值和协方差矩阵。

2.2 高维随机变量函数的统计线性化近似求解方案

若函数fi(x)(i=1，2，…，n)关于状态向量x(t)高阶可微，将fi(x)在m处进行多元函数的泰勒展开，可以得到

(6)

其中，s代表泰勒展开阶数。

针对不同的fi(x)适当选取泰勒级数展开阶数并忽略皮亚诺余项可以简化描述函数i(x)的计算。一般选取泰勒级数展开阶数为1～2阶即可满足精度要求。

若选取泰勒级数展开阶数为2阶，则

(7)

(8)

根据概率论相关知识，式(8)中

则式(8)可化简为

(9)

类似的，若选取泰勒级数展开阶数为1阶，则

i(x)=fi(m)

(10)

考虑到六自由度动能拦截器制导控制系统维数较高，本文选取1阶泰勒级数展开进行近似求解。

矩阵N为

(11)

2.3 降维拦截概率计算公式

利用式(5)可以快速得到弹目相对位置的均值矩阵和协方差矩阵，但落点散布的均值和协方差矩阵并不能作为制定多弹协同策略的直接依据，为此，需要得到拦截器对对应目标的拦截概率。

本小节根据相对位置的均值矩阵和协方差矩阵给出拦截概率公式，并将其降维以提高计算效率。

假设n维随机变量x服从正态分布N(m，p)，其概率密度函数[10]为

(12)

其中，p为随机向量的协方差矩阵，m为随机向量的均值向量。

协方差矩阵中的非对角元素表征着状态变量之间的相关程度。若p为对角矩阵，则x为独立的随机向量，若p为非对角矩阵，则x为非独立的随机向量。一般地，对于有控导弹系统，相对位置向量x=[Δx，Δy，Δz]T为非独立随机向量。

设动能拦截器与目标的脱靶量在R以内认为动能拦截器成功拦截目标，则拦截概率可以表示为

(13)

其中，Ω为一球型区域，方程为

(Δx-mΔx)2+(Δy-mΔy)2+(Δz-mΔz)2≤R2

(14)

将概率密度函数g(x)沿Δz方向积分，则可以将三维概率密度函数转化为二维概率密度函数，进一步提高计算效率。

设

则概率密度函数g(x)可以降维表示为式(15)。具体推导过程见附录A。

(Φ(h(mΔz+R))-Φ(h(mΔz-R)))×

(15)

式(15)中Φ(h(Δz))表示累积分布函数，可通过读取标准正态分布表获得，且

拦截概率相应地可以表示为

(16)

其中，Ω′为一圆型区域，方程为

(Δx-mΔx)2+(Δy-mΔy)2≤R2

(17)

3 快速计算方法在动能拦截器制导控制模型中的应用

将本文提出的拦截概率快速计算方法应用于动能拦截器制导控制模型中，以验证该方法在复杂系统上的可行性。

3.1 动能拦截器制导控制模型

动能拦截器制导控制模型内部关系可由图1表示。

图1 动能拦截器制导控制模型

在动能拦截器制导控制模型中，主要随机状态变量如表1所示，共计26项，组成高维随机状态向量x(t)。

表1 动能拦截器制导控制模型中主要随机变量

3.2 动能拦截器制导控制系统的统计线性化

为保证统计分析的准确性，需对动能拦截器制导控制系统中的非线性方程均进行统计线性化。为便于计算，本文对于低维函数方程进行完全统计线性化，采用2.2节中提出的近似求解方案对高维函数方程统计线性化。

3.3 动能拦截器的拦截概率

将动能拦截器制导控制系统统计线性化后，运用协方差分析技术可以得到系统响应的状态矢量x(t)的均值和协方差的方程式(18)。

(18)

其中，m(t)为状态向量x(t)的均值

m(t)=[mωy4，mωz4，mamy2，mamz2，mm，mVm，mθm，mψvm，

mXm，mYm，mZm，maTx4，maTy4maTz4，mVTx，mVTy，

mVTz，mxT，myT，mzT，mR，mq1，mq2，mΔq1，mΔq2]T

4 仿真结果及分析

弹目相对位置的均值矩阵和协方差矩阵是拦截弹效能评估和命中概率计算的主要依据。本文分别通过快速计算方法和蒙特卡洛方法对典型拦截场景进行计算机仿真，对比两种方法得到的弹目相对位置的均值和方差，以分析快速计算方法的准确性；同时给出两种方法的仿真时间，以说明快速计算方法的有效性。

4.1 仿真条件设置

初始值设置：动能拦截器初始位置(2360km，0km，0km)，初始速度大小6000m/s，初始速度方向指向目标，目标初始位置(2235.5km，79km，0.5km)，目标初始速度(3795m/s，-1881m/s，20m/s)。

扰动量设置：高低角偏差Δq1均值设置为0，标准差设置为5e-5rad；方位角偏差Δq2均值设置为0，标准差设置为5e-5rad。

步长设置：1ms，仿真末端步长减小至0.01ms。

硬件环境：利用同一台工作站仿真。

4.2 仿真结果对比与分析

图2和图3分别为目标不采取机动策略时弹目相对位置的均值和方差对比曲线。从图2可以看出，采用快速计算方法得到的相对位置均值和蒙特卡洛方法得到的结果是一致的，在制导控制系统的作用下逐渐减小。从图3可以看出采用快速计算方法得到的相对位置方差和蒙特卡洛方法得到的结果也是一致的，在测角误差的影响下相对位置的散布逐渐变大。

将落点均值和协方差代入式(16)中并进行积分运算便可以得到脱靶量特定范围以内的拦截概率。从表2中可以看出，根据快速计算方法得到的命中概率和进行500次蒙特卡洛仿真得到的结果非常接近，且从表2中可以看出，采用本文方法对命中概率的计算的快速性较好，大约是500次蒙特卡洛仿真所用时间的1/90，且10秒的运行时间基本可以满足弹上实时运算的要求。

图2 均值对比曲线

表2 蒙特卡洛方法和快速计算方法在不同脱靶量限制下的拦截概率

表3 蒙特卡洛方法和快速计算方法用时对比

图4 t(aT≠0)=2s时拦截概率对比图

图5 t(aT≠0)=4s时拦截概率对比图

图4和图5分别为在目标机动力持续作用时间为2s和4s的条件下，根据快速计算方法得到的拦截概率对比图。从图4和图5可以看出，对于不同的目标机动力大小、机动力持续时间和突防起始位置，拦截概率有很大的区别。这表明只有在探测到目标的机动方式时，才能获得准确地拦截概率。而在实际战场环境下，对于目标的探测时长是有限的，这也进一步说明了提出拦截概率快速计算方法的必要性。

5 结论

对于动能拦截器拦截概率的统计计算，一般通过蒙特卡洛方法实现，时间成本很高。应用本文所提方法对动能拦截器制导控制系统进行处理可以一次性地得到状态变量的均值矩阵和协方差矩阵，将其代入降维拦截概率计算公式即可解算拦截概率，在保证计算准确度的前提下可以大幅提高计算效率，能够为战场实时决策提供有效信息。文中给出的仿真实例说明了所提方法的可行性和有效性。

本文提出的拦截概率计算方法的创新点主要有：为满足多弹协同策略的制定对于获取拦截概率的实时性要求，首次结合协方差分析描述函数法进行拦截概率的计算；基于泰勒级数展开对传统CADET进行改进，对复杂系统的统计线性化采用近似求解方案，能够快速获得状态向量的均值矩阵和协方差矩阵；采用分部积分法对三维拦截概率计算公式进行降维，能够进一步减少计算时间。