数值微分在设备组合弹道精度评估中的应用

何京江，杨涵，郭茂耘，柴毅; HE Jingjiang; YANG Han; GUO Maoyun; CHAI Yi

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

数值微分在设备组合弹道精度评估中的应用 PDF

- ORCID：
何京江 ¹
✉
- ORCID：
杨涵 ²
- ORCID：
郭茂耘 ²
- ORCID：
柴毅 ²

1. 西昌卫星发射中心，四川西昌 615000； 2. 重庆大学自动化学院，重庆 400044

中图分类号： TM571.6

最近更新：2025-02-27

DOI：10.11835/j.issn.1000-582X.2025.01.001

摘要

设备组合弹道精度评估是测量设备布站设计和设备组合弹道选优方案制定的重要基础。设备组合弹道精度评估通常采用基于误差传播原理的弹道精度评估算法,求解弹道分量关于测元的雅克比矩阵是该算法的核心。然而，对于复杂的弹道方程很难求得雅克比矩阵的解析式。为了解决复杂弹道方程的雅克比矩阵求解难题，提出基于样条数值微分的设备组合弹道精度评估算法，通过构造弹道分量数值队列，进行基于样条插值的数值微分计算出弹道分量关于测元的雅克比矩阵，计算得到设备组合弹道理论精度。最后，通过与单台设备定位方程的理论雅克比矩阵和弹道精度进行对比,验证基于数值微分的设备组合弹道精度评估算法的有效性和实用性。

关键词

弹道精度评估模型; 雅克比矩阵; 样条插值; 数值微分

组合测量设备精度分析评估对航天发射具有重要意义。快速准确地完成测量设备精度分析评估，直接影响到设备测量站配置设计、弹道选优策略制定、安全控制裕量确定和入轨精度评估等航天发射指挥决策活动的开展。近年来，随着中国商业航天、空间站和大型星座建设等任务进入密集实施阶段，对测控数据分析的精度提出了更高要求。针对原有弹道精度评估模型中，求取复杂弹道方程的雅克比矩阵困难，采用微商近似计算，精度不够理想的情况，笔者提出了一种基于样条插值数值微分的雅克比矩阵求取方法，高效准确地满足航天发射任务对弹道精度处理分析要求，为设备组合选优方案论证应用和指挥决策智能优化提供技术支持。

目前，在航天发射场设备组合弹道理论精度评估算法的弹道精度分析、评估和修正，以及数值计算应用等领域，国内外研究学者开展了如下研究：1）在精度分析处理总体研究方面，张金槐，刘利生，唐雪梅等^[

1⁃3]系统论述了外测设备数据处理的理论、方法和技术，对外弹道测量设备的测量精度影响因素、精度评估和鉴定等方法开展了相应研究；2）在精度误差模型分析评估方面，常晓华等^{[参考文献 4

百度学术}4]结合弹道飞行器质点动力学模型，分析了垂线偏差对弹道参数等的影响机理，建立了数学模型，并给出弹道落点偏差计算方法。薛斌等^{[参考文献 5

百度学术}5]基于误差源对弹道精度影响权重，提出了一维弹道精度评估方法。雷晓云^{[参考文献 6

百度学术}6]利用蒙特卡洛法分析研究影响射击精度的因素，为提高修正精度提供数据分析支持。冯燕来^{[参考文献 7

百度学术}7]分析了测量设备误差、地球引力势能等影响弹道导弹落点预测精度的因素，较准确得到弹道预报的落点偏差。郭军海^{[参考文献 8

百度学术}8]就弹道测量误差模型、弹道模型、测量数据实时和事后融合处理方法，以及组合弹道处理精度评估方法等开展研究；3）在系统误差与精度修正方面，宫志华^{[参考文献 9

百度学术}9]研究了非合作目标高精度弹道测量系统误差，提出4项修正内容，为雷达组网和数据处理提供了参考。王子鉴^{[参考文献 10

百度学术}10]提出一种等权归一化最小二乘法的弹道误差系数辨识方法，用于全程弹道落点偏差折合精度评估，获得了更高的辨识精度和更好的多源观测信息拟合结果。王若晗等^{[参考文献 11

百度学术}11]在对光电经纬仪误差分析的基础上，对测量数据进行分段最小二乘和3次样条平滑处理，减小随机误差对弹道解算结果的影响，提高弹道处理精度。Ji等^{[参考文献 12

百度学术}12]针对助推段误差分析，利用高斯过程模型，对助推段弹道开展预测研究，仿真表明该方法提高了预测精度；4）在样条插值等数值计算方法用于求解方程方面，张哲^{[参考文献 13

百度学术}13]针对经典数值积分算法精度严重依赖小积分步长，难以快速计算的问题，提出一种局部配点反馈迭代算法，提高了轨道动力学方程解算速度。Saumya等^{[参考文献 14

百度学术}14]在薛定谔方程的数值求解中应用B样条有限元分析法，该方法不仅具有良好的守恒性，还降低了求解过程中非对角矩阵计算的复杂性。Hosseinian等^{[参考文献 15

百度学术}15]基于薄板样条函数，给出了一种求解第二类非线性Volterra积分方程的数值法，该方法具有稳定收敛且对计算资源要求低的特点。

由以上研究可看出，弹道理论精度分析评估具有重要的工程应用价值，其分析评估工作通常基于误差传播原理的弹道精度评估模型开展，该模型的核心是求解弹道分量关于测量元素（简称测元）的雅克比矩阵^[

2]。然而，在数据处理软件研制和应用实践中，对于复杂的弹道方程很难求得其雅克比矩阵解析式，难以编程实现自动计算设备组合弹道精度。为避免求解复杂弹道方程的雅克比矩阵解析式，并能易于编程实现，提出一种基于数值微分的设备组合弹道理论精度评估算法。算法的核心是通过构造测量元素和弹道分量数值队列，对弹道分量数值队列进行基于样条插值的数值微分计算，得到弹道分量关于测元的雅克比矩阵和设备组合弹道理论精度。

1 设备组合弹道精度评估模型

设备组合弹道精度是指由测量设备的基本误差传播到测量数据，再经过组合测量和数学方法计算传播到弹道的总误差统计量。因此，通常采用误差传播模型估算设备组合弹道理论精度。

1.1 弹道精度评估内容

设备组合弹道精度受许多与测量设备自身有关和无关的因素影响：飞行弹道、测量体制、测量设备数量、精度和误差指标、跟踪性能、设备测量站配置和弹道估计方法等^[

2]，如图1所示。

图1 弹道精度影响因素示意

Fig. 1 The illustration of factors affecting trajectory accuracy

常用的弹道精度评估有误差传播估算和精度检验2种途径。其中，误差传播估算包括仿真分析和理论分析2大类，主要用于设备选优组合方案设计和论证；精度检验方法可分为“内符合”和“外符合”2大类，分别用于实时处理弹道和事后处理弹道的精度评估^[

2]。

本文研究的弹道精度评估属于理论分析类别。在航天发射任务测量方案的论证和设计中，通常采用基于测量误差传播估计的理论分析方法进行弹道理论精度评估。

1.2 测量误差传播分析

设备组合测量的误差传播是由多个测量数据（元素）经数学方法计算得到所求量的过程。设待求量弹道为 $X = {[X_{1}, X_{2}, \dots, X_{m}]}^{T}$ ，设备测元为 $U = {[U_{1}, U_{2}, \dots, U_{n}]}^{T}$ 。设测元 $U$ 的误差为 $Δ U$ ，引起的弹道 $X$ 的误差为 $Δ X$ ， $A$ 是弹道 $X$ 关于测元 $U$ 的雅克比矩阵。通常测元 $Δ U$ 较小，则弹道误差可由下式计算

Δ X = A \times Δ U

，

经推导，得到测量误差传播关系为

P_{X} = A P_{U} A^{T}

，

其中： $P_{X}$ 为弹道 $X$ 的协方差矩阵； $P_{U}$ 为测元 $U$ 的协方差矩阵。

通常，测元 $U$ 的误差由多个误差源引起。如弹道测量设备的误差源有设备误差、电波折射误差、时统误差、跟踪部位误差等。需要将这些误差源引起的误差带入到误差传播公式中。设误差源向量 $a = {[a_{1}, a_{2}, \dots, a_{k}]}^{T}$ ，由误差源传播到测元 $U$ 的关系式为

U = g （ a ）

，

经推导，可得测量误差源 $a$ 传播到测元 $U$ 的误差传播矩阵形式为

P_{U} = \frac{\partial U}{\partial a} P_{a} {(\frac{\partial U}{\partial a})}^{T} = B P_{a} B^{T}

，

式中： $B = \frac{\partial U}{\partial a}$ 为测元 $U$ 关于误差源 $a$ 的雅克比矩阵。进一步可得引入测量误差源后的误差传播关系式，即

P_{X} = \frac{\partial X}{\partial U} P_{U} {(\frac{\partial X}{\partial U})}^{T} = \frac{\partial X}{\partial U} \frac{\partial U}{\partial a} P_{a} {(\frac{\partial U}{\partial a})}^{T} {(\frac{\partial X}{\partial U})}^{T} = A B P_{a} B^{T} A^{T}

，

其中， $P_{a}$ 为测量误差源的协方差矩阵。

1.3 弹道精度理论分析

在航天发射任务中，设备组合弹道理论精度分析是设备测量站配置几何优化和设备组合弹道方案论证设计的基础。在飞行器弹道测量中，设备测元通常有以下几类

$A$ ：目标相对测站的方位角；

$E$ ：目标相对测站的俯仰角；

$R$ ：目标至测站（主站）的径向距离；

$\dot{R}$ ：目标至测速设备测站（主站）的径向距离变化率；

$\dot{r}$ ：目标至测速设备主、副站间的径向距离差变化率；

$\dot{S}$ ：目标至测速设备主、副站间的径向距离和变化率。

设测量向量为 $U = {[u_{1}, u_{2}, \dots, u_{N}]}^{T}$ ，测量元素 $u_{j}$ 为上述的某类测元， $j = 1,2, \dots, N$ ，其中 $N$ 表示参与弹道计算的测元数量；设弹道向量为 $X = {[x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6}]}^{T}$ ， $x_{i}$ 为弹道参数， $i = 1,2, \dots, 6$ 。设由测量元素解算弹道的关系式为

X = f (U)

，

（1）

由式（1）可得测量向量 $U$ 的测量误差传播到弹道向量 $X$ 的理论精度估算公式为

P_{X} = \frac{\partial X}{\partial U} P_{U} {(\frac{\partial X}{\partial U})}^{T}

，

（2）

其中， $P_{X}$ 为弹道参数向量 $X$ 误差协方差矩阵， $P_{U}$ 为测量向量 $U$ 的测量误差协方差矩阵，即

P_{X} = [\begin{matrix} σ_{x_{1}}^{2} & σ_{x_{1} x_{2}} & \dots & σ_{x_{1} x_{6}} \\ σ_{x_{2} x_{1}} & σ_{x_{2}}^{2} & \dots & σ_{x_{2} x_{6}} \\ \dots & \dots & \dots & \dots \\ σ_{x_{6} x_{1}} & σ_{x_{6} x_{2}} & \dots & σ_{x_{6}}^{2} \end{matrix}]

，

（3）

P_{U} = [\begin{matrix} σ_{u_{1}}^{2} & σ_{u_{1} u_{2}} & \dots & σ_{u_{1} u_{N}} \\ σ_{u_{2} u_{1}} & σ_{u_{2}}^{2} & \dots & σ_{u_{1} u_{N}} \\ \dots & \dots & \dots & \dots \\ σ_{u_{N} u_{1}} & σ_{u_{N} u_{2}} & \dots & σ_{u_{N}}^{2} \end{matrix}]

。

（4）

$\frac{\partial X}{\partial U}$ 为弹道参数向量 $X$ 关于测量向量 $U$ 的雅克比矩阵，即

\frac{\partial X}{\partial U} = \frac{\partial (x_{1}, x_{2}, x_{3}, x_{4}, x_{5}, x_{6})}{\partial (u_{1}, u_{2}, \dots, u_{N})} = [\begin{matrix} \frac{\partial x_{1}}{\partial u_{1}} & \frac{\partial x_{1}}{\partial u_{2}} & \dots & \frac{\partial x_{1}}{\partial u_{N}} \\ \frac{\partial x_{2}}{\partial u_{1}} & \frac{\partial x_{2}}{\partial u_{2}} & \dots & \frac{\partial x_{2}}{\partial u_{N}} \\ \dots & \dots & \dots & \dots \\ \frac{\partial x_{6}}{\partial u_{1}} & \frac{\partial x_{6}}{\partial u_{2}} & \dots & \frac{\partial x_{6}}{\partial u_{N}} \end{matrix}]

。

（5）

由参考文献[

2]可知，影响弹道理论精度的关键因素为观测方程、测量站配置几何、目标飞行轨迹和测量设备精度等。研究是在观测方程、测量站配置几何、目标飞行轨迹确定的情况下，分析测量设备精度对设备组合弹道精度的影响。影响测量设备精度的误差源有设备误差、电波折射误差、时统误差、跟踪部位误差等。在方案设计阶段时，只考虑设备系统误差和随机误差对弹道精度的影响，即

P_{U}

为测量误差之间的协方差矩阵。

通常，测量方程 $X = f (U)$ 是多元、非常规、复杂方程，式（5）所示的雅克比矩阵求解复杂，偏微分方程不易采用解析方法求得。并且，不同测量体制的测量方程不一样，不利于计算机编程实现，有必要研究简便有效、易于编程实现的通用算法来求解雅克比矩阵，以实现通用的弹道精度评估算法。

2 基于样条插值的数值微分方法

为解决复杂弹道方程的雅克比矩阵求解难题，利用基于样条插值的数值微分方法，完成雅克比矩阵中偏导数的数值计算，满足工程需求。其中，样条插值是为数值微分做准备的。

2.1 样条插值

插值就是利用已知数据求出给定点的函数值。可选用的插值函数有代数多项式（多项式）、三角函数（傅里叶级数）和有理函数等。由于代数多项式结构不复杂，数值计算和理论分析简单，因此，研究选用代数多项式作为插值函数。样条插值函数作为目前被认为的一种有效的代数多项式插值工具^[

16]，能够克服拉格朗日插值的不收敛性和提高分段线性插值在节点处的不光滑性。采用3次样条插值

S_{3} (u)

，定义如下

对于给定的函数值 $x_{i} = f (u_{i}) (i = 0,1, \dots, n)$ ，

a = u_{0} < u_{1} < \dots < u_{n} = b

。

若存在函数 $S (u)$ 满足：

1) $S (u)$ 在每个子区间 $[u_{i}, u_{i + 1}] (i = 0,1, \dots, n - 1)$ 上具有不高于3次的多项式，

2) $S (u)$ 满足插值条件 $S (u_{i}) = x_{i} (i = 0,1, \dots, n)$ ，即过插值点，

3) $S (u)$ 在每个内节点 $u_{i} (i = 1,2, \dots, n - 1)$ 上具有二阶连续导数，则称 $S (u)$ 为函数 $f (u)$ 关于节点 $u_{i} (i = 0,1, \dots, n)$ 的3次样条函数。

通过推导， $S (u)$ 在子区间 $[u_{i - 1}, u_{i}] (i = 1,2, \dots, n)$ 上的表达式为

S (u) = \frac{(u_{i} {- u)}^{3}}{6 h_{i}} M_{i - 1} + \frac{(u - u_{i - 1})^{3}}{6 h_{i}} M_{i} + (y_{i - 1} - \frac{h_{i}^{2}}{6} M_{i - 1}) \frac{u_{i} - u}{h_{i}} + (y_{i} - \frac{h_{i}^{2}}{6} M_{i}) \frac{u - u_{i - 1}}{h_{i}}

，

（6）

其中， $u_{i - 1} \leq u \leq u_{i}$ ， $h_{i} = u_{i} - u_{i - 1}$ 。只要确定 $M_{i} (i = 0,1, \dots, n)$ 这 $n + 1$ 个值，便可确定，3次样条插值函数。按边界条件，给定2个断点处的一阶导数值，即

\dot{S} (u_{0}) = {\dot{x}}_{0}

，

\dot{S} (u_{n}) = {\dot{x}}_{n}

，

（7）

可推导出求解 $M_{i} (i = 0,1, \dots, n)$ 的方程组，即

[\begin{matrix} 2 & 1 \\ μ_{1} & 2 & λ_{1} \\ μ_{2} & 2 & λ_{2} \\ ⋱ & ⋱ & ⋱ \\ μ_{n - 1} & 2 & λ_{n - 1} \\ 1 & 2 \end{matrix}] [\begin{matrix} M_{0} \\ M_{1} \\ M_{2} \\ ⋮ \\ M_{n - 1} \\ M_{n} \end{matrix}] = [\begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ ⋮ \\ d_{n - 1} \\ d_{n} \end{matrix}]

，

（8）

其中

μ_{i} = \frac{h_{i}}{h_{i} + h_{i + 1}} ， \begin{matrix} λ_{i} = 1 - μ_{i} ， \end{matrix}

d_{i} = \frac{6}{h_{i} + h_{i + 1}} (\frac{y_{i + 1} - y_{i}}{h_{i + 1}} - \frac{y_{i} - y_{i - 1}}{h_{i}}) ， \begin{matrix} i = 1,2, \dots, n - 1 \end{matrix}

，

d_{0} = \frac{6}{h_{1}} (\frac{y_{1} - y_{0}}{h_{1}} - {\dot{y}}_{0}) ， \begin{matrix} d_{n} = \frac{6}{h_{n}} ({\dot{y}}_{n} - \frac{y_{n} - y_{n - 1}}{h_{n}}) \end{matrix}

。

2.2 数值微分

由式（6）推导出，3次样条插值 $S (u)$ 的导数 $\dot{S} (u)$ 在子区间 $[u_{i - 1}, u_{i}] (i = 1,2, \dots, n)$ 上的表达式为

\dot{S} (u) = \frac{(u - u_{i - 1})^{2}}{2 h_{i}} M_{i} - \frac{(u_{i} {- u)}^{2}}{2 h_{i}} M_{i - 1} + \frac{h_{i}^{}}{6} (M_{i - 1} - M_{i}) + \frac{y_{i} - y_{i - 1}}{h_{i}}

，

u_{i - 1} \leq u \leq u_{i}

，

特别的， $\dot{S} (u_{i})$ 即为所求的 $u_{i} (i = 1,2, \dots, n)$ 点处的数值微分值，即

\dot{S} (u_{i}) = \frac{h_{i}}{6} (M_{i - 1} + 2 M_{i}) + \frac{y_{i} - y_{i - 1}}{h_{i}}

，

i = 1, \dots, n

。

（9）

3 基于数值微分的设备组合弹道理论精度评估算法设计

3.1 问题描述

在设备组合弹道理论精度评估时，已知各设备测量元素的理论测量精度、各测量设备在时刻 $t_{k}$ 的测元 $U (t_{k})$ ，根据测量方程 $X = f (U)$ 求得设备组合弹道值 $X (t_{k})$ 。因此，本算法要解决的问题可描述如下：

已知：测元 $U (t_{k})$ ，式（1）所示的测量方程 $X = f (U)$ ，式（4）所示的测量误差协方差矩阵 $P_{U}$

求：式（3）所示的表征弹道 $X$ 精度的误差协方差矩阵 $P_{X} (t_{k})$ 。

其中

P_{X} (t_{k}) = \frac{\partial X (t_{k})}{\partial U (t_{k})} P_{U} {(\frac{\partial X (t_{k})}{\partial U (t_{k})})}^{T}

，

可见，求出上式中各时刻雅克比矩阵 $\frac{\partial X (t_{k})}{\partial U (t_{k})}$ ，是解决问题的关键。算法中涉及的测元、弹道值、弹道精度、雅克比矩阵都是飞行时 $t_{k}$ 的函数，为了表述简洁，下面公式描述中将略去变量 $t_{k}$ 。

3.2 算法步骤

为求解矩阵 $\frac{\partial X}{\partial U}$ 中弹道分量对测元的偏导数 $\frac{\partial x_{i}}{\partial u_{j}}$ ( $i = 1,2, \dots, 6$ ， $j = 1,2, \dots, N$ , $N$ 表示参与弹道计算的测元数量),基于数值微分的设备组合弹道理论精度评估算法步骤如下。

1）构造测元变量数值队列

基于测量向量 $U = {[u_{1}, u_{2}, \dots, u_{N}]}^{T}$ 构造变量队列 $\{U_{j, w}^{'}\}$ ， $w = - W - 1, \dots, 0, \dots, W + 1$ ， $W$ 为正整数， $2 W + 1$ 为队列中变量个数。定义向量 $G = {[g_{1}, g_{2}, \dots, g_{N}]}^{T}$ ， $g_{j}$ 表示测元分量 $u_{j}$ 的单位增量。定义向量 $Q_{j} = {[q_{1}, \dots, q_{k}, \dots, q_{N}]}^{T}$ ， $Q_{j}$ 的第 $j$ 位元素 $q_{j}$ 为1，其余为0，即

q_{k} = \{\begin{matrix} \begin{matrix} 1, & k = j \end{matrix} ， \\ \begin{matrix} 0, & k \neq j \end{matrix} 。 \end{matrix}

按如下公式构造变量 $U_{j, w}^{'}$

U_{j, w}^{'} = U + w Q_{j}^{T} Q_{j}^{} G = [\begin{matrix} u_{1} \\ \begin{matrix} ⋮ \\ u_{j - 1} \end{matrix} \\ u_{j} \\ \begin{matrix} u_{j + 1} \\ ⋮ \end{matrix} \\ u_{N} \end{matrix}] + w [\begin{matrix} 0 \\ \begin{matrix} ⋮ \\ 0 \end{matrix} \\ g_{j} \\ \begin{matrix} 0 \\ ⋮ \end{matrix} \\ 0 \end{matrix}]

[\begin{matrix} u_{1} \\ \begin{matrix} ⋮ \\ u_{j - 1} \end{matrix} \\ u_{j, w}^{'} \\ \begin{matrix} u_{j + 1} \\ ⋮ \end{matrix} \\ u_{N} \end{matrix}]

，

（10）

当 $w = 0$ 时， $U_{j, 0}^{'} = U$ ， $u_{j, w}^{'} = u_{j} + w g_{j}$ 。

2）构造弹道分量数值队列

将变量 $U_{j, w}^{'}$ 代入测量方程 $X = f (U)$ ，得到弹道向量 $X_{j, w}^{'}$

X_{j, w}^{'} = f (U_{j, w}^{'})

，

其中， $X_{j, w}^{'} = {[x_{1, j, w}^{'}, x_{2, j, w}^{'}, \dots, x_{6, j, w}^{'}]}^{T}$ ， $x_{i, j, w}^{'}$ 为弹道向量 $X_{j, w}^{'}$ 的第 $i$ 个弹道分量， $i = 1, \dots, 6$ 。

进一步由数值队列 $\{U_{j, w}^{'}\}$ 和 $\{X_{j, w}^{'}\}$ ，建立起测量变量和第 $i$ 个弹道分量组成的数值队列 $\{(U_{j, w}^{'}, x_{i, j, w}^{'})\}$ ， $i = 1, \dots, 6$ ， $j = 1, \dots, N$ ， $w = - W - 1, \dots, 0, \dots, W + 1$ 。

3）求数值队列的样条插值函数

要为 $\{(U_{j, w}^{'}, x_{i, j, w}^{'})\}$ 构造如式（6）所示的样条插值函数 $S_{i} (u)$ ，需要式（7）所示的样条曲线边界条件。根据弹道分量平滑变化特点，提出如下边界条件计算方法

{\dot{S}}_{i} (u_{- W}) = (\frac{x_{i, j, - W}^{'} - x_{i, j, - W - 1}^{'}}{u_{j, - W}^{'} - u_{j, - W - 1}^{'}} + \frac{x_{i, j, - W + 1}^{'} - x_{i, j, - W}^{'}}{u_{j, - W + 1}^{'} - u_{j, - W}^{'}}) / 2 = \frac{x_{i, j, - W + 1}^{'} - x_{i, j, - W - 1}^{'}}{2 g_{j}}

，

{\dot{S}}_{i} (u_{W}) = (\frac{x_{i, j, W + 1}^{'} - x_{i, j, W}^{'}}{u_{j, W + 1}^{'} - u_{j, W}^{'}} + \frac{x_{i, j, W}^{'} - x_{i, j, W - 1}^{'}}{u_{j, W}^{'} - u_{j, W - 1}^{'}}) / 2 = \frac{x_{i, j, W + 1}^{'} - x_{i, j, W - 1}^{'}}{2 g_{j}}

。

在此边界条件计算方法之下，方程式（8）表述为

[\begin{matrix} 2 & 1 \\ μ_{- W + 1} & 2 & λ_{- W + 1} \\ μ_{- W + 2} & 2 & λ_{- W + 2} \\ ⋱ & ⋱ & ⋱ \\ μ_{W - 1} & 2 & λ_{W - 1} \\ 1 & 2 \end{matrix}] [\begin{matrix} M_{- W} \\ M_{- W + 1} \\ M_{- W + 2} \\ ⋮ \\ M_{W - 1} \\ M_{W} \end{matrix}] = [\begin{matrix} d_{- W} \\ d_{- W + 1} \\ d_{- W + 2} \\ ⋮ \\ d_{W - 1} \\ d_{W} \end{matrix}]

，

（11）

进一步计算得到计算方程式（11）中所需的参数

μ_{w} = = \frac{1}{2} ， \begin{matrix} λ_{w} = 1 - μ_{w} = \frac{1}{2} \begin{matrix} ， & w = - W, \dots, W \end{matrix} \end{matrix}

，

d_{w} = \frac{3}{g_{j}^{2}} (x_{i, j, w + 1}^{'} - 2 x_{i, j, w}^{'} + x_{i, j, w - 1}^{'}) ， \begin{matrix} w = - W + 1, \dots, 0, \dots, W - 1 \end{matrix}

，

d_{- W} = \frac{6}{g_{j}^{2}} (x_{i, j, - W + 1}^{'} - 2 x_{i, j, - W}^{'} + x_{i, j, - W - 1}^{'})

，

d_{W} = \frac{6}{g_{j}^{2}} (x_{i, j, W + 1}^{'} - 2 x_{i, j, W}^{'} + x_{i, j, W - 1}^{'})

。

求解方程式（11）即可得到 ${[M_{- W}, M_{- W + 1}, \dots, 0, \dots, M_{W - 1}, M_{W}]}^{T}$ ，从而得到 $S (u)$ 在子区间 $[u_{j, w - 1}^{'}, u_{j, w}^{'}]$ （ $w = - W + 1, \dots, W$ ）上的表达式为

S (u) = \frac{(u_{j, w}^{'} {- u)}^{3}}{6 h_{w}} M_{w - 1} + \frac{(u - u_{j, w - 1}^{'})^{3}}{6 h_{w}} M_{w} + (x_{i, j, w - 1}^{'} - \frac{h_{w}^{2}}{6} M_{w - 1}) \frac{u_{j, w}^{'} - u}{h_{w}} + (x_{i, j, w}^{'} - \frac{h_{w}^{2}}{6} M_{w}) \frac{u - u_{j, w - 1}^{'}}{h_{w}}

，

其中， $u_{j, w - 1}^{'} < u \leq u_{j, w}^{'}$ 。 $S (u)$ 即为数值队列 $\{(U_{j, w}^{'}, x_{i, j, w}^{'})\}$ 的样条函数。

4）计算雅克比矩阵

3次样条插值 $S (u)$ 的导数 $\dot{S} (u)$ 在子区间 $[u_{j, w - 1}^{'}, u_{j, w}^{'}]$ （ $w = - W + 1, \dots, W$ ）上的表达式为

\dot{S} (u) = \frac{(u - u_{j, w - 1}^{'})^{2}}{2 h_{w}} M_{w} - \frac{(u_{j, w}^{'} {- u)}^{2}}{2 h_{w}} M_{w - 1} + \frac{h_{w}}{6} (M_{w - 1} - M_{w}) + \frac{x_{i, j, w}^{'} - x_{i, j, w - 1}^{'}}{h_{w}}

，

$\dot{S} (t)$ 在 $u_{j, w}^{'} (w = - W + 1, \dots, 0, \dots, W - 1)$ 点处的表达式为

\dot{S} (u_{j, w}^{'}) = \frac{h_{w}}{6} (M_{w - 1} + 2 M_{w}) + \frac{x_{i, j, w}^{'} - x_{i, j, w - 1}^{'}}{h_{w}} ， \begin{matrix} w = - W + 1, \dots, 0, \dots, W - 1 \end{matrix}

，

得到雅克比矩阵元素 $\frac{\partial x_{i}}{\partial u_{j}}$ 的值，即

\frac{\partial x_{i}}{\partial u_{j}} = \dot{S} (U_{j, 0}^{'}) = \frac{h_{0}}{6} (M_{- 1} + 2 M_{0}) + \frac{x_{i, j, 0}^{'} - x_{i, j, - 1}^{'}}{h_{0}}

，

（12）

由式（12）计算得到的 $\frac{\partial x_{i}}{\partial u_{j}} (i = 1, \dots, 6 ， j = 1, \dots, N)$ ，就组成了式（5）所示的弹道关于测元的雅克比矩阵 $\frac{\partial X}{\partial U}$ 。

5）评估弹道精度

由式（2）可得测量向量 $U$ 的测量误差传播到弹道参数向量 $X$ 的精度 $P_{X}$ ，即

P_{X} = \frac{\partial X}{\partial U} P_{U} {(\frac{\partial X}{\partial U})}^{T}

，

（13）

在进行设备组合弹道理论精度计算时，测量向量 $U$ 的误差协方差矩阵 $P_{U}$ 里的测量误差应考虑测量系统误差和随机误差。

4 算法验证

如果利用设计的算法计算得到的雅克比矩阵和弹道精度评估值，与理论上的雅克比矩阵和弹道精度值相差无几，则可以得出算法是有效可行的结论。

4.1 验证方法

为了便于得到理论上的雅克比矩阵和弹道精度值，采用简单的单台雷达或光学定位方程。由单台设备定位方程直接解析求导，计算得到弹道序列 $X_{k}$ 相对测元序列 $U_{k}$ 的理论雅克比矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{0}$ 和弹道理论精度 $P_{X}^{0}$ ；利用提出的方法计算得到弹道序列 $X_{k}$ 相对测元序列 $U_{k}$ 的雅克比矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{S}$ 和弹道精度 $P_{X}^{S}$ ； ${(\frac{\partial X_{k}}{\partial U_{k}})}^{S}$ 和 $P_{X}^{S}$ 的各分量值分别与 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{0}$ 和 $P_{X}^{0}$ 的各分量值比较，如果误差很小满足任务要求，则提出的算法是有效可行的。

4.2 验证过程

为突出重点和简化计算，采用设备测量坐标系，假设目标初始点位于测量系原点。

1）由理论弹道反算测元

选取某次任务的 $0 \sim 280.0$ s的理论弹道序列 $[t_{k} ， x_{k} ， y_{k} ， z_{k}]$ $^{T}$ ，反算测元序列 $[t_{k} ， R_{k} ， A_{k} ， E_{k}]$ $^{T}$ ， $T_{k}$ 为弹道点 $k$ 的弹道时间， $k = 1,2, \dots, K$ ， $K = 288$ ，为弹道点数。设备测元计算如式（14）所示

\{\begin{array}{l} R_{k} = \sqrt[]{x_{k}^{2} + y_{k}^{2} + z_{k}^{2}} ， \\ A_{k} = a r c s i n (z_{k} / \sqrt[]{x_{k}^{2} + z_{k}^{2}}) ， \\ E_{k} = a r c t a n (y_{k} / \sqrt[]{x_{k}^{2} + z_{k}^{2}}) 。 \end{array}

（14）

2）计算理论雅克比矩阵

单台设备定位方程为

\{\begin{array}{l} x = R \times c o s E \times c o s A ， \\ y = R \times s i n E ， \\ z = R \times c o s E \times s i n A 。 \end{array}

（15）

由式（15）可推导出弹道位置相对于测元的偏导数解析表达式，得到 $T_{k}$ 时刻理论雅克比矩阵为

{(\frac{\partial X_{k}}{\partial U_{k}})}^{0} = (\begin{matrix} c o s E_{k} c o s A_{k} & - R_{k} c o s E_{k} s i n A_{k} & - R_{k} s i n E_{k} c o s A_{k} \\ s i n E & 0 & R_{k} c o s E_{k} \\ c o s E_{k} s i n A_{k} & R_{k} c o s E_{k} c o s A_{k} & - R_{k} s i n E_{k} s i n A_{k} \end{matrix})

，

其中， $R_{k}$ 、 $A_{k}$ 、 $E_{k}$ 由式（14）计算得到。

3）利用提出的数值微分方法计算雅克比矩阵

在各个弹道点时刻 $T_{k}$ ，按照算法步骤，计算弹道位置分量 $X_{k}$ 相对于测元 $U_{k}$ 的偏导数矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{S}$ 。式（10）中的变量取值如下 $W = 2$ ， $g_{R} = 10$ m， $g_{A} = 0.000 15$ rad， $g_{E} = 0.000 15$ rad。

4）比较雅克比矩阵精度

雅克比矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{S}$ 计算精度用绝对误差 $σ_{M}$ 来衡量，计算如下式所示

σ_{M}^{} = \frac{1}{K} \sum_{k = 1}^{K} σ_{M_{k}}^{} ， σ_{M k}^{2} = \frac{1}{9} \sum_{i = 1}^{3} \sum_{j = 1}^{3} {({(\frac{\partial X_{k}}{\partial U_{k}})}_{i, j}^{S} - {(\frac{\partial X_{k}}{\partial U_{k}})}_{i, j}^{0})}^{2}

，

计算得到的绝对误差值为： $σ_{M}^{} = 0.007 9$ 。其中： ${(\frac{\partial X_{k}}{\partial U_{k}})}_{i, j}^{S}$ 表示矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{S}$ 的第 $i$ 行第 $j$ 列的元素； ${(\frac{\partial X_{k}}{\partial U_{k}})}_{i, j}^{0}$ 表示矩阵 ${(\frac{\partial X_{k}}{\partial U_{k}})}^{0}$ 的第 $i$ 行第 $j$ 列的元素。

5）计算弹道精度

按照式（13）分别计算出弹道位置分量评估精度 $P_{X}^{0}$ 和本算法的弹道评估精度 $P_{X}^{S}$

P_{X}^{0} = {(\frac{\partial X_{k}}{\partial U_{k}})}^{0} P_{U} {({(\frac{\partial X_{k}}{\partial U_{k}})}^{0})}^{T}

，

P_{X}^{S} = {(\frac{\partial X_{k}}{\partial U_{k}})}^{S} P_{U} {({(\frac{\partial X_{k}}{\partial U_{k}})}^{S})}^{T}

。

假设测元误差之间相互独立，测元误差协方差矩阵 $P_{U} = d i a g (σ_{R}^{2}, σ_{A}^{2}, σ_{E}^{2})$ ，按根据设备测量精度要求，测量误差取值为： $σ_{R} = 46$ m， $σ_{A} = 0.000 45$ rad， $σ_{E} = 0.000 49$ rad。

6）比较弹道精度

协方差矩阵 $P_{X}^{S}$ 计算精度用绝对误差 $σ_{X}$ 来衡量，计算方法如下式所示

\begin{matrix} σ_{X}^{} = \frac{1}{K} \sum_{k = 1}^{K} σ_{X_{k}}^{} ， & σ_{X_{k}}^{2} = \sum_{i = 1}^{3} {({(P_{X}^{S})}_{i, i} - {(P_{X}^{0})}_{i, i})}^{2} \end{matrix}

，

其中： ${(P_{X}^{S})}_{i, i}$ 表示矩阵 $P_{X}^{S}$ 的第 $i$ 行第 $i$ 列的元素； ${(P_{X}^{0})}_{i, i}$ 表示矩阵 $P_{X}^{0}$ 的第 $i$ 行第 $i$ 列的元素。

统计的位置精度与理论位置精度的误差均值为： $σ_{X}^{} = 2.09 \times 10^{- 5}$ m。任务方案中对弹道精度要求的描述是精确到m。因此，本算法的计算精度满足工程要求。从上面验证过程和结果可知，采用研究提出的基于数值微分的弹道理论精度评估算法计算精度较高，算法有效可行。

5 结束语

基于误差传播原理的精度评估方法广泛应用于精度评估领域。但算法中的雅克比矩阵的求解比较复杂，复杂组合弹道方程的偏导数不易采用解析方法求得，难以编码实现计算机自动处理。研究提出基于数值微分实现简便有效、易于编程实现的通用方法来求解雅克比矩阵，设计了基于数值微分的设备组合弹道理论精度评估算法，验证算法的有效性和可行性。提出的算法可为设备组合弹道精度评估系统，为测控设备测量站配置方案，设备组合选优方案的制定提供算法基础。

参考文献

张金槐. 远程火箭精度分析与评估[M]. 长沙: 国防科技大学出版社, 1995. [百度学术]

Zhang J H. Accuracy analysis and evaluation of long range rocket[M]. Changsha: National University of Defense Science and Technology Press, 1995. (in Chinese) [百度学术]

刘利生. 外弹道测量数据处理[M]. 北京: 国防工业出版社, 2002. [百度学术]

Liu L S. Data processing of exterior ballistic measurement[M]. Beijing: National Defense Industry Press, 2002. (in Chinese) [百度学术]

唐雪梅, 蔡洪, 杨华波, 等. 导弹武器精度分析与评估[M]. 北京: 国防工业出版社, 2015. [百度学术]

Tang X M, Chai H, Yang H B, et al. Accuracy analysis and evaluation of missile weapon[M]. Beijing: National Defense Industry Press, 2015. (in Chinese) [百度学术]

常晓华, 蒋鲁佳, 杨锐, 等. 垂线偏差对弹道落点精度的影响分析[J]. 国防科技大学学报, 2017, 39(4): 1-5. [百度学术]

Chang X H，Jiang L J，Yang Ｒ, et al. Analysis on effects of vertical deflection for trajectory impact accuracy[J]. Journal of National University of Defense Technology. 2017, 39(4):1-5. (in Chinese) [百度学术]

薛冰, 霍鹏飞, 杨小会. 基于误差合成的射程修正系统精度评估[J]. 探测与控制学报, 2016, 38(4): 58-61. [百度学术]

Xue B, Huo P F, Yang X H. Accuracy evaluation of range correction system based on error synthesis[J]. Journal of Detection & Control, 2016, 38(4): 58-61. (in Chinese) [百度学术]

雷晓云, 张志安. 基于蒙特卡罗法的一维弹道修正弹落点精度分析[J]. 系统仿真学报, 2016, 28(7): 1685-1691. [百度学术]

Lei X Y, Zhang Z A. Research of impact point accuracy of one-dimensional trajectory correction projectile based on Monte Carlo[J]. Journal of System Simulation, 2016, 28(7): 1685-1691. (in Chinese) [百度学术]

冯燕来, 王红杰, 李旭等. 基于外弹道修正理论的导弹落点预测精度评估方法[J]. 指挥信息系统与技术, 2017, 8(4): 48-52. [百度学术]

Feng Y L, Wang H J, Li X, et al. Precision evaluation method for ballistic missile impact-point prediction based on exterior ballistic correction theory[J]. Command Information System and Technology, 2017, 8(4): 48-52. (in Chinese) [百度学术]

郭军海. 弹道测量数据融合技术[M]. 北京: 国防工业出版社, 2012. [百度学术]

Guo J H. Ballistic measurement data fusion technology[M]. Beijing: National Defense Industry Press, 2012. (in Chinese) [百度学术]

宫志华, 刘洋, 陈春江. 分布式雷达对非合作目标弹道测量精度分析[J]. 弹道学报, 2017, 29(3): 43-48. [百度学术]

Gong Z H, Liu Y, Chen C J. Analysis on precision of distributed radar measuring[J]. Ballistic Journal, 2017, 29(3): 43-48. (in Chinese) [百度学术]

王子鉴. 多弹头惯性/星光复合制导精度评估与弹道折合方法研究[D]. 长沙: 国防科技大学, 2018. [百度学术]

Wang Z J. Study on accuracy evaluation and trajectory conversion of inertial-stellar integrated guidance system based on multiple warheads[D]. Changsha: National University of Defense Science and Technology, 2018. (in Chinese) [百度学术]

王若晗. 光电经纬仪系统误差分析及修正方法研究[D]. 西安: 西安理工大学, 2023. [百度学术]

Wang R H. Research on system error analysis and correction method of photoelectric theodolite[D]. Xi’an: Xi’an University of Technology, 2023. (in Chinese) [百度学术]

Ji R P, Liang Y, Xu L F, et al. Trajectory prediction of ballistic missiles using Gaussian process error model[J]. Chinese Journal of Aeronautics, 2022, 35(1): 458-469. [百度学术]

张哲, 代洪华, 冯浩阳, 等. 初值约束与两点边值约束轨道动力学方程的快速数值计算方法[J]. 力学学报, 2022, 54(2):503-516. [百度学术]

Zhang Z，Dai H H, Feng H Y, et al. Efficient numerical method for orbit dynamic functions with initial value and two-point boundary-value constraints[J]. Journal of Theoretical and Applied Mechanics, 2022, 54(2): 503-516. (in Chinese) [百度学术]

Saumya R J, Senapati A, Stability, convergence and error analysis of B-spline collocation with Crank-Nicolson method and finite element methods for numerical solution of Schrodinger equation arises in quantum mechanics[J]. Physica Scripta, 2023, 98(11):115232. [百度学术]

Hosseinian A, Assari P, Dehghan M. The numerical solution of nonlinear delay Volterra integral equations using the thin plate spline collocation method with error analysis[J].Computational and Applied Mathematics,Computational and Applied Mathematics 2023, 42(2): 83. [百度学术]

吴开腾. 数值计算方法及其程序实现[M]. 北京: 科学出版社, 2015. [百度学术]

Wu K T. Numerical calculation method and its program implementation[M]. Beijing: Science Press, 2015. (in Chinese) [百度学术]

期刊社主页

编辑部首页

期刊介绍

编委会

数据库收录

过刊浏览

联系我们