CN105550497A - 一种高精度的弹道修正方法 - Google Patents

Abstract

本发明公开了一种高精度的弹道修正方法，其特征在于，包括：利用噪声估计器对状态噪声和量测噪声进行估计和校正的方法、利用噪声估计器完成一步预测的过程和基于噪声估计器和渐消因子完成弹道修正的过程。本发明提供的一种高精度的弹道修正方法，设计的弹道修正方法是在拟线性最优滤波过程中，利用量测数据一方面对状态值进行预测，另一方面通过过程噪声估计器和量测噪声估计器自适应的估计和校正噪声的特性参数以调节滤波增益矩阵，降低了噪声参数估计不精确对弹道修正结果的影响；同时引入渐消因子对陈旧的数据进行合理的渐消，增加新的量测数据的权重，达到减小滤波估计误差、抑制由于系统模型不精确以及状态突变导致的滤波发散的目的。

Description

一种高精度的弹道修正方法

技术领域

本发明涉及一种高精度的弹道修正方法，具体地说是一种涉及弹道修正过程中抑制发散并提高精度的方法，尤其考虑弹道模型的非线性、时变性和噪声干扰等因素，旨在抑制由于弹道模型的建模误差及噪声统计特性的估计误差在弹道数据修正过程中引起的发散，提高弹道形成的精度。本发明能有效保证弹道修正过程的稳定性以及弹道轨迹形成的精确性，适用于炮位雷达的弹道轨迹形成系统。

背景技术

炮位雷达利用获取的弹丸的飞行参数进行弹道修正，形成弹道轨迹以便进行目标跟踪。在大空域中飞行的弹丸受到各种因素的影响，如弹丸飞行的高度、马赫数、气体动力参数的测量值和计算值不精确等会导致弹丸的模型参数变化很大。另外，弹道修正过程中设置的经验估计值和实际弹丸飞行的参数之间也会存在一定的差异。在使用拟线性最优平滑滤波算法进行弹道修正时，弹道模型离散化和线性化所带来的误差、弹道模型描述不精确、所采用的过程噪声和量测噪声的统计特性与实际情况的偏差大等因素都可能导致状态估计值与弹丸实际飞行数据偏离，致使滤波精度下降，无法满足实际工程的要求。实际环境变化等随机因素还会引起系统过程噪声和量测噪声异常变化，在滤波过程中会产生增益矩阵和预测误差方差的错误估计，这种情况下增益矩阵将失去合适的加权作用，新的量测数据对状态估计值的调节修正作用将逐渐降低，从而出现滤波发散现象，导致弹道修正失败，无法继续跟踪目标。

发明内容

本发明所要解决的技术问题是，提供一种利用噪声估计器对状态噪声和量测噪声进行估计和校正的高精度的弹道修正方法；更进一步地，本发明提供一种利用噪声估计器完成一步预测，降低拟线性最优平滑滤波过程中不确定的噪声统计特性和其他随机干扰等因素对弹道修正精度的影响的高精度的弹道修正方法；更进一步地，本发明提供一种基于噪声估计器和渐消因子完成弹道修正，渐消因子对状态估计的一步预测误差的协方差和滤波增益矩阵进行自适应的调整，抑制由于弹道模型不精确以及滤波过程中状态突变导致的发散现象的高精度的弹道修正方法。

为解决上述技术问题，本发明采用的技术方案为：

一种高精度的弹道修正方法，其特征在于，包括：利用噪声估计器对状态噪声和量测噪声进行估计和校正的方法，该法包括以下步骤：

S01，弹道模型：

弹道模型采用的状态方程为：

$\frac{d}{dt}X\left(t\right)=f\left(X\right(t\left)\right)+\mathrm{\Γ}\left(t\right)W\left(t\right)---\left(1\right)$

式(1)中，X(t)表示弹丸的状态，W(t)是过程噪声，Γ(t)是噪声驱动阵；

弹道模型采用的量测方程为：

Z(t)＝h(X(t))+V(t)(2)

式(2)中Z(t)表示量测变量，V(t)为量测噪声；

对式(1)、式(2)进行离散化后得到：

X_k+1＝G(X_k)+Γ(X_k)W_k(3)

Z_k+1＝h(X_k+1)+V_k+1(4)

式(3)、式(4)中， $G\left(X_k\right)=X_k+f\left(X_k\right)T+A\left(X_k\right)f\left(X_k\right)\frac{T^2}{2},$ X_k＝X(t)|_t＝kT，T＝t_k+1-t_k为采样间隔，t_k为第k个采样时间；W_k＝W(t)|_t＝kT，V_k+1＝V(t)|_t＝(k+1)T；

S02，弹道修正的初始值计算：

对弹道轨迹进行修正时，定义状态矢量为：

X＝[xyzv_xv_yv_zca_c](5)

式(5)中x，y和z分别表示弹丸的射向距离、高度和侧偏；v_x，v_y和v_z分别表示弹丸速度的三个分量；c和a_c分别表示弹道系数和弹道修正系数；

设雷达测得炮弹的最初的两点采样值为(x₁,y₁,z₁)和(x₂,y₂,z₂)，由此得到状态矢量X的初始估值为：

${\hat{X}}_0=\[\begin{array}{cccccccc}{x}_2& y_2& z_2& \frac{x_2-x_1}{T}& \frac{y_2-y_1}{T}& \frac{z_2-z_1}{T}& c& a_c\end{array}\]---\left(6\right)$

初始预测误差方差阵为：

式(7)中，为系统随机误差的方差，是弹道系数和弹道修正系数均方差的经验估算值。

本发明还包括利用噪声估计器完成一步预测的过程，该过程包括以下步骤：

S03，弹道修正：

S03-1，计算一步预测平滑值：

S03-1-1，利用初始值和式(7)的预测结果计算相关参数，

令 ${\hat{X}}_k^{*}={\hat{X}}_k,(k=0,1,2,...)---\left(8\right)$

计算状态函数：

$f\left({\hat{X}}_k^{*}\right)=f\left(X\right(t\left)\right){|}_{X\left(t\right)={\hat{X}}_k^{*}}=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mathrm{cN}}_1{x}_{4r}+a_c{v}_x/v_r\\ -{\mathrm{cN}}_1{v}_y-g+a_c{v}_y/v_r\\ -{\mathrm{cN}}_1{x}_{6r}+B_Z+a_c{v}_z/v_r\\ 0\\ 0\end{array}\right]{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(9\right)$

式(9)中x_4r＝v_x-w_x和x_6r＝v_z-w_z分别表示弹丸在x和z方向受风速影响后的实际速度，w_x和w_z为风速在x和z方向的分量，ρ为空气密度，S为弹丸最大截断面积即d为弹径，m为弹体质量；C_x(Ma)为空气阻力系数，是马赫数Ma的函数，g为海拔y处的重力加速度，B_Z为侧向升力加速度；

计算状态函数的偏导数

$A\left({\hat{X}}_k^{*}\right)=\frac{\∂f\left(X\right(t\left)\right)}{\∂X\left(t\right)}{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(10\right)$

计算一步状态转移矩阵

$\mathrm{\φ}\left({\hat{X}}_k^{*}\right)=\frac{\∂G\left(X\right)}{\∂X}{|}_{X={\hat{X}}_k^{*}}---\left(11\right)$

计算过程噪声的方差矩阵和量测噪声的方差阵：

${\hat{Q}}_k=E\left\{W_k{W}_k^H\right\}---\left(12\right)$

${\hat{R}}_k=E\left\{V_k{V}_k^H\right\}---\left(13\right)$

式(12)、(13)中，H表示共轭转置；

S03-1-2，进行预测估计：

计算一步预测估计值：

$\begin{array}{c}{\hat{X}}_{k+1/k}=G\left({\hat{X}}_k^{*}\right)+\mathrm{\φ}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\\ ={\hat{X}}_k^{*}+f\left({\hat{X}}_k^{*}\right)T+A\left({\hat{X}}_k^{*}\right)f\left({\hat{X}}_k^{*}\right)\frac{T^2}{2}+\mathrm{\φ}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\end{array}---\left(14\right)$

计算一步预测误差协方差：

${\hat{P}}_{k+1/k}=\mathrm{\φ}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\φ}}^H\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(15\right)$

计算量测转移矩阵：

$B\left({\hat{X}}_{k+1/k}\right)=\frac{\∂h\left(X\right)}{\∂X}{|}_{X={\hat{X}}_{k+1/k}}---\left(16\right)$

计算增益矩阵：

$K_k={\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right){\[B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right)+{\hat{R}}_k^{*}\]}^{-1}---\left(17\right)$

S03-1-3，利用噪声估计器修正噪声：

对过程噪声的方差矩阵进行修正：

${\hat{Q}}_{k+1}=(1-d_{k+1}){\hat{Q}}_k+d_{k+1}(K_k{\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T{K}_k^T+P_{k+1/k}-{\mathrm{\Φ}}_k{P}_k{\mathrm{\Φ}}_k^T)---\left(18\right)$

对量测噪声的方差阵进行修正：

${\hat{R}}_{k+1}=(1-d_{k+1}){\hat{R}}_k+d_{k+1}({\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T-B_{k+1/k}{\hat{P}}_{k+1/k}{B}_{k+1/k}^T)---\left(19\right)$

式(18)、(19)中，残差新息的计算为：

${\tilde{Z}}_{k+1}=Z_{k+1}-H\left(X_{k+1/k}\right)---\left(20\right)$

式(18)、(19)中，d_k＝(1-α)/(1-α^k)，0＜α＜1为遗忘因子；

S03-1-4，进行自适应平滑滤波：

计算滤波平滑值：

$\begin{array}{c}{\hat{X}}_{k/k+1}={\hat{X}}_k+P_k{\mathrm{\φ}}_k\left({\hat{X}}_k^{*}\right)B^T\left({\hat{X}}_{k+1/k}\right)\\ \·{\[B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\]}^{-1}\\ \·\[Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\]\end{array}---\left(21\right)$

S03-1-5，迭代计算

令重复S03-1-1到S03-1-5计算过程至少三次，如此进行至少三次迭代后即可求得的近似值，完成一步预测平滑值的计算。

本发明还包括基于噪声估计器和渐消因子完成弹道修正的过程，该过程包括：

S03-2，利用渐消因子修正滤波结果：

S03-2-1，计算渐消因子，令渐消因子

${\mathrm{\&lambda;}}_k=\frac{Tr\left(N_k\right)}{Tr\left(M_k\right)}=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_k,& {\mathrm{\&lambda;}}_k\&GreaterEqual;1\\ 1,& {\mathrm{\&lambda;}}_k<1\end{array}\right.---\left(22\right)$

式(22)中，Tr表示求矩阵的迹； $N_k=b_k-H\left({\hat{X}}_{k/k+1}\right)Q_k^{*}{H}^T\left({\hat{X}}_{k/k+1}\right)-l_k{R}_k;$ $M_k=H\left({\hat{X}}_{k/k+1}\right)\mathrm{\&phi;}\left({\hat{X}}_{k/k+1}\right){\hat{P}}_{k+1/k}{\mathrm{\&phi;}}^T\left({\hat{X}}_{k/k+1}\right)H^T\left({\hat{X}}_{k/k+1}\right);$ l_k＝1-d_k；

S03-2-2，进行一步预测结果修正，

计算带渐消因子的一步预测的误差协方差矩阵

$P_{k+1/k}={\mathrm{\&lambda;}}_k\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\&phi;}}^T\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(23\right)$

计算增益矩阵

$K_{k+1}=P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right){\&lsqb;B\left({\hat{X}}_{k+1/k}\right)P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\&rsqb;}^{-1}---\left(24\right)$

计算滤波估计值

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}\&lsqb;Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\&rsqb;---\left(25\right)$

计算滤波误差方差阵

$P_{k+1}=\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;P_{k+1/k}{\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;}^H+K_{k+1}{R}_{k+1}{K}_{k+1}^H---\left(26\right)$

S03-3，令k+1→k，重复S03-1到S03-3步骤，完成弹道修正过程。

所述迭代计算的次数为三次。

本发明提供的一种高精度的弹道修正方法，本发明设计的弹道修正方法是在拟线性最优滤波过程中，利用量测数据一方面对状态值进行预测，另一方面通过过程噪声估计器和量测噪声估计器自适应的估计和校正噪声的特性参数以调节滤波增益矩阵，降低了噪声参数估计不精确对弹道修正结果的影响；同时引入渐消因子对陈旧的数据进行合理的渐消，增加新的量测数据的权重，通过对预测误差的协方差矩阵和滤波增益矩阵进行实时在线调整，达到减小滤波估计误差、抑制由于系统模型不精确以及状态突变导致的滤波发散的目的。

说明书附图

图1为本发明基于自适应渐消因子的弹道修正流程图；

图2为现有技术基于拟线性最优平滑滤波算法弹道修正的发散现象图；

图3本发明中高精度弹道修正方法的仿真结果图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1～图3所示，为了验证本发明的基于噪声估计器和渐消因子的弹道修正方法有效性，在炮弹初速度选为1100m/s，射向角选为25o条件下模拟生成雷达量测数据，分别采用本发明的弹道轨迹修正算法以及拟线性最优平滑滤波算法，对雷达测量得到的炮弹飞行数据进行滤波分析，形成的弹道轨迹的三维仿真结果如图2、图3所示，其中“仿真数据”为模拟的雷达量测数据，“滤波数据”为滤波估计的弹道轨迹。如图2所示，图2中出现了滤波发散现象，而图3中使用基于噪声估计器和渐消因子的弹道修正方法有效的抑制了发散。

一种高精度的弹道修正方法，其特征在于，包括：利用噪声估计器对状态噪声和量测噪声进行估计和校正的方法，该法包括以下步骤：

S01，弹道模型：

弹道模型采用的状态方程为：

$\frac{d}{dt}X\left(t\right)=f\left(X\right(t\left)\right)+\mathrm{\&Gamma;}\left(t\right)W\left(t\right)---\left(1\right)$

式(1)中，X(t)表示弹丸的状态，W(t)是过程噪声，Γ(t)是噪声驱动阵；

弹道模型采用的量测方程为：

Z(t)＝h(X(t))+V(t)(2)

式(2)中Z(t)表示量测变量，V(t)为量测噪声；

对式(1)、式(2)进行离散化后得到：

X_k+1＝G(X_k)+Γ(X_k)W_k(3)

Z_k+1＝h(X_k+1)+V_k+1(4)

式(3)、式(4)中， $G\left(X_k\right)=X_k+f\left(X_k\right)T+A\left(X_k\right)f\left(X_k\right)\frac{T^2}{2},$ X_k＝X(t)|_t＝kT，T＝t_k+1-t_k为采样间隔，t_k为第k个采样时间；W_k＝W(t)|_t＝kT，V_k+1＝V(t)|_t＝(k+1)T；

S02，弹道修正的初始值计算：

对弹道轨迹进行修正时，定义状态矢量为：

X＝[xyzv_xv_yv_zca_c](5)

式(5)中x，y和z分别表示弹丸的射向距离、高度和侧偏；v_x，v_y和v_z分别表示弹丸速度的三个分量；c和a_c分别表示弹道系数和弹道修正系数；

设雷达测得炮弹的最初的两点采样值为(x₁,y₁,z₁)和(x₂,y₂,z₂)，由此得到状态矢量X的初始估值为：

${\hat{X}}_0=\&lsqb;\begin{array}{cccccccc}{x}_2& y_2& z_2& \frac{x_2-x_1}{T}& \frac{y_2-y_1}{T}& \frac{z_2-z_1}{T}& c& a_c\end{array}\&rsqb;---\left(6\right)$

初始预测误差方差阵为：

式(7)中，为系统随机误差的方差，是弹道系数和弹道修正系数均方差的经验估算值。

本发明还包括利用噪声估计器完成一步预测的过程，该过程包括以下步骤：

S03，弹道修正：

S03-1，计算一步预测平滑值：

S03-1-1，利用初始值和式(7)的预测结果计算相关参数，

令 ${\hat{X}}_k^{*}={\hat{X}}_k,(k=0,1,2,...)---\left(8\right)$

计算状态函数：

$f\left({\hat{X}}_k^{*}\right)=f\left(X\right(t\left)\right){|}_{X\left(t\right)={\hat{X}}_k^{*}}=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mathrm{cN}}_1{x}_{4r}+a_c{v}_x/v_r\\ -{\mathrm{cN}}_1{v}_y-g+a_c{v}_y/v_r\\ -{\mathrm{cN}}_1{x}_{6r}+B_Z+a_c{v}_z/v_r\\ 0\\ 0\end{array}\right]{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(9\right)$

式(9)中x_4r＝v_x-w_x和x_6r＝v_z-w_z分别表示弹丸在x和z方向受风速影响后的实际速度，w_x和w_z为风速在x和z方向的分量，ρ为空气密度，S为弹丸最大截断面积即d为弹径，m为弹体质量；C_x(Ma)为空气阻力系数，是马赫数Ma的函数，g为海拔y处的重力加速度，B_Z为侧向升力加速度；

计算状态函数的偏导数

$A\left({\hat{X}}_k^{*}\right)=\frac{\&part;f\left(X\right(t\left)\right)}{\&part;X\left(t\right)}{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(10\right)$

计算一步状态转移矩阵

$\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)=\frac{\&part;G\left(X\right)}{\&part;X}{|}_{X={\hat{X}}_k^{*}}---\left(11\right)$

计算过程噪声的方差矩阵和量测噪声的方差阵：

${\hat{Q}}_k=E\left\{W_k{W}_k^H\right\}---\left(12\right)$

${\hat{R}}_k=E\left\{V_k{V}_k^H\right\}---\left(13\right)$

式(12)、(13)中，H表示共轭转置；

S03-1-2，进行预测估计：

计算一步预测估计值：

$\begin{array}{c}{\hat{X}}_{k+1/k}=G\left({\hat{X}}_k^{*}\right)+\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\\ ={\hat{X}}_k^{*}+f\left({\hat{X}}_k^{*}\right)T+A\left({\hat{X}}_k^{*}\right)f\left({\hat{X}}_k^{*}\right)\frac{T^2}{2}+\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\end{array}---\left(14\right)$

计算一步预测误差协方差：

${\hat{P}}_{k+1/k}=\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\&phi;}}^H\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(15\right)$

计算量测转移矩阵：

$B\left({\hat{X}}_{k+1/k}\right)=\frac{\&part;h\left(X\right)}{\&part;X}{|}_{X={\hat{X}}_{k+1/k}}---\left(16\right)$

计算增益矩阵：

$K_k={\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right){\&lsqb;B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right)+{\hat{R}}_k^{*}\&rsqb;}^{-1}---\left(17\right)$

S03-1-3，利用噪声估计器修正噪声：

对过程噪声的方差矩阵进行修正：

${\hat{Q}}_{k+1}=(1-d_{k+1}){\hat{Q}}_k+d_{k+1}(K_k{\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T{K}_k^T+P_{k+1/k}-{\mathrm{\&Phi;}}_k{P}_k{\mathrm{\&Phi;}}_k^T)---\left(18\right)$

对量测噪声的方差阵进行修正：

${\hat{R}}_{k+1}=(1-d_{k+1}){\hat{R}}_k+d_{k+1}({\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T-B_{k+1/k}{\hat{P}}_{k+1/k}{B}_{k+1/k}^T)---\left(19\right)$

式(18)、(19)中，残差新息的计算为：

${\tilde{Z}}_{k+1}=Z_{k+1}-H\left(X_{k+1/k}\right)---\left(20\right)$

式(18)、(19)中，d_k＝(1-α)/(1-α^k)，0＜α＜1为遗忘因子；

S03-1-4，进行自适应平滑滤波：

计算滤波平滑值：

$\begin{array}{c}{\hat{X}}_{k/k+1}={\hat{X}}_k+P_k{\mathrm{\&phi;}}_k\left({\hat{X}}_k^{*}\right)B^T\left({\hat{X}}_{k+1/k}\right)\\ \&CenterDot;{\&lsqb;B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\&rsqb;}^{-1}\\ \&CenterDot;\&lsqb;Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\&rsqb;\end{array}---\left(21\right)$

S03-1-5，迭代计算

令重复S03-1-1到S03-1-5计算过程至少三次，如此进行至少三次迭代后即可求得的近似值，完成一步预测平滑值的计算。

本发明还包括基于噪声估计器和渐消因子完成弹道修正的过程，该过程包括：

S03-2，利用渐消因子修正滤波结果：

S03-2-1，计算渐消因子，令渐消因子

${\mathrm{\&lambda;}}_k=\frac{Tr\left(N_k\right)}{Tr\left(M_k\right)}=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_k,& {\mathrm{\&lambda;}}_k\&GreaterEqual;1\\ 1,& {\mathrm{\&lambda;}}_k<1\end{array}\right.---\left(22\right)$

式(22)中，Tr表示求矩阵的迹； $N_k=b_k-H\left({\hat{X}}_{k/k+1}\right)Q_k^{*}{H}^T\left({\hat{X}}_{k/k+1}\right)-l_k{R}_k;$ $M_k=H\left({\hat{X}}_{k/k+1}\right)\mathrm{\&phi;}\left({\hat{X}}_{k/k+1}\right){\hat{P}}_{k+1/k}{\mathrm{\&phi;}}^T\left({\hat{X}}_{k/k+1}\right)H^T\left({\hat{X}}_{k/k+1}\right);$ l_k＝1-d_k；

S03-2-2，进行一步预测结果修正，

计算带渐消因子的一步预测的误差协方差矩阵

$P_{k+1/k}={\mathrm{\&lambda;}}_k\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\&phi;}}^T\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(23\right)$

计算增益矩阵

$K_{k+1}=P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right){\&lsqb;B\left({\hat{X}}_{k+1/k}\right)P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\&rsqb;}^{-1}---\left(24\right)$

计算滤波估计值

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}\&lsqb;Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\&rsqb;---\left(25\right)$

计算滤波误差方差阵

$P_{k+1}=\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;P_{k+1/k}{\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;}^H+K_{k+1}{R}_{k+1}{K}_{k+1}^H---\left(26\right)$

S03-3，令k+1→k，重复S03-1到S03-3步骤，完成弹道修正过程。

所述迭代计算的次数为三次。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种高精度的弹道修正方法，其特征在于，包括：利用噪声估计器对状态噪声和量测噪声进行估计和校正的方法，该法包括以下步骤：

S01，弹道模型：

弹道模型采用的状态方程为：

$\frac{d}{dt}X\left(t\right)=f\left(X\right(t\left)\right)+\mathrm{\&Gamma;}\left(t\right)W\left(t\right)---\left(1\right)$

式(1)中，X(t)表示弹丸的状态，W(t)是过程噪声，Γ(t)是噪声驱动阵；

弹道模型采用的量测方程为：

Z(t)＝h(X(t))+V(t)(2)

式(2)中Z(t)表示量测变量，V(t)为量测噪声；

对式(1)、式(2)进行离散化后得到：

X_k+1＝G(X_k)+Γ(X_k)W_k(3)

Z_k+1＝h(X_k+1)+V_k+1(4)

式(3)、式(4)中， $G\left(X_k\right)=X_k+f\left(X_k\right)T+A\left(X_k\right)f\left(X_k\right)\frac{T^2}{2},$ X_k＝X(t)|_t＝kT，T＝t_k+1-t_k为采样间隔，t_k为第k个采样时间；V_k+1＝V(t)|_t＝(k+1)T；

S02，弹道修正的初始值计算：

对弹道轨迹进行修正时，定义状态矢量为：

X＝[xyzv_xv_yv_zca_c](5)

式(5)中x，y和z分别表示弹丸的射向距离、高度和侧偏；v_x，v_y和v_z分别表示弹丸速度的三个分量；c和a_c分别表示弹道系数和弹道修正系数；

设雷达测得炮弹的最初的两点采样值为(x₁,y₁,z₁)和(x₂,y₂,z₂)，由此得到状态矢量X的初始估值为：

${\hat{X}}_0=\&lsqb;\begin{array}{cccccccc}{x}_2& y_2& z_2& \frac{x_2-x_1}{T}& \frac{y_2-y_1}{T}& \frac{z_2-z_1}{T}& c& a_c\end{array}\&rsqb;---\left(6\right)$

初始预测误差方差阵为：

$P_0=\left[\begin{array}{cccccccc}{\mathrm{\&sigma;}}_x^2& & & & & & & \\ & {\mathrm{\&sigma;}}_y^2& & & & 0& & \\ & & {\mathrm{\&sigma;}}_z^2& & & & & \\ & & & 2{\mathrm{\&sigma;}}_x^2/T^2& & & & \\ & & & & 2{\mathrm{\&sigma;}}_y^2/T^2& & & \\ & & 0& & & 2{\mathrm{\&sigma;}}_z^2/T^2& & \\ & & & & & & {\mathrm{\&sigma;}}_c^2& \\ & & & & & & & {\mathrm{\&sigma;}}_{ac}^2\end{array}\right]---\left(7\right)$

式(7)中， 为系统随机误差的方差， 是弹道系数和弹道修正系数均方差的经验估算值。

2.根据权利要求1所述的一种高精度的弹道修正方法，其特征在于：还包括利用噪声估计器完成一步预测的过程，该过程包括以下步骤：

S03，弹道修正：

S03-1，计算一步预测平滑值：

S03-1-1，利用初始值和式(7)的预测结果计算相关参数，

令 ${\hat{X}}_k^{*}={\hat{X}}_k(k=0,1,2,...)---\left(8\right)$

计算状态函数：

$f\left({\hat{X}}_k^{*}\right)=f\left(X\right(t\left)\right){|}_{X\left(t\right)={\hat{X}}_k^{*}}=\left[\begin{array}{c}{v}_x\\ v_y\\ v_z\\ -{\mathrm{cN}}_1{x}_{4r}+a_c{v}_x/v_r\\ -{\mathrm{cN}}_1{v}_y-g+a_c{v}_y/v_r\\ -{\mathrm{cN}}_1{x}_{6r}+B_Z+a_c{v}_z/v_r\\ 0\\ 0\end{array}\right]{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(9\right)$

式(9)中x_4r＝v_x-w_x和x_6r＝v_z-w_z分别表示弹丸在x和z方向受风速影响后的实际速度，w_x和w_z为风速在x和z方向的分量， ρ为空气密度，S为弹丸最大截断面积即d为弹径，m为弹体质量；C_x(Ma)为空气阻力系数，是马赫数Ma的函数，g为海拔y处的重力加速度，B_Z为侧向升力加速度；

计算状态函数的偏导数

$A\left({\hat{X}}_k^{*}\right)=\frac{\&part;f\left(X\right(t\left)\right)}{\&part;X\left(t\right)}{|}_{X\left(t\right)={\hat{X}}_k^{*}}---\left(10\right)$

计算一步状态转移矩阵

$\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)=\frac{\&part;G\left(X\right)}{\&part;X}{|}_{X={\hat{X}}_k^{*}}---\left(11\right)$

计算过程噪声的方差矩阵和量测噪声的方差阵：

${\hat{Q}}_k=E\left\{W_k{W}_k^H\right\}---\left(12\right)$

${\hat{R}}_k=E\left\{V_k{V}_k^H\right\}---\left(13\right)$

式(12)、(13)中，H表示共轭转置；

S03-1-2，进行预测估计：

计算一步预测估计值：

$\begin{array}{c}{\hat{X}}_{k+1/k}=G\left({\hat{X}}_k^{*}\right)+\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\\ ={\hat{X}}_k^{*}+f\left({\hat{X}}_k^{*}\right)T+A\left({\hat{X}}_k^{*}\right)f\left({\hat{X}}_k^{*}\right)\frac{T^2}{2}+\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)({\hat{X}}_k-{\hat{X}}_k^{*})\end{array}---\left(14\right)$

计算一步预测误差协方差：

${\hat{P}}_{k+1/k}=\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\&phi;}}^H\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(15\right)$

计算量测转移矩阵：

$B\left({\hat{X}}_{k+1/k}\right)=\frac{\&part;h\left(X\right)}{\&part;X}{|}_{X={\hat{X}}_{k+1/k}}---\left(16\right)$

计算增益矩阵：

$K_k={\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right){\&lsqb;B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^H\left({\hat{X}}_{k+1/k}\right)+{\hat{R}}_k^{*}\&rsqb;}^{-1}---\left(17\right)$

S03-1-3，利用噪声估计器修正噪声：

对过程噪声的方差矩阵进行修正：

${\hat{Q}}_{k+1}=(1-d_{k+1}){\hat{Q}}_k+d_{k+1}(K_k{\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T{K}_k^T+P_{k+1/k}-{\mathrm{\&Phi;}}_k{P}_k{\mathrm{\&Phi;}}_k^T)---\left(18\right)$

对量测噪声的方差阵进行修正：

${\hat{R}}_{k+1}=(1-d_{k+1}){\hat{R}}_k+d_{k+1}({\tilde{Z}}_{k+1}{\tilde{Z}}_{k+1}^T-B_{k+1/k}{\hat{P}}_{k+1/k}{B}_{k+1/k}^T)---\left(19\right)$

式(18)、(19)中，残差新息的计算为：

${\tilde{Z}}_{k+1}=Z_{k+1}-H\left(X_{k+1/k}\right)---\left(20\right)$

式(18)、(19)中，d_k＝(1-α)/(1-α^k)，0＜α＜1为遗忘因子；

S03-1-4，进行自适应平滑滤波：

计算滤波平滑值：

$\begin{array}{c}{\hat{X}}_{k/k+1}={\hat{X}}_k+P_k{\mathrm{\&phi;}}_k\left({\hat{X}}_k^{*}\right)B^T\left({\hat{X}}_{k+1/k}\right)\\ \&CenterDot;{\&lsqb;B\left({\hat{X}}_{k+1/k}\right){\hat{P}}_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\&rsqb;}^{-1}\\ \&CenterDot;\&lsqb;Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\&rsqb;\end{array}---\left(21\right)$

S03-1-5，迭代计算

令重复S03-1-1到S03-1-5计算过程至少三次，如此进行至少三次迭代后即可求得的近似值，完成一步预测平滑值的计算。

3.根据权利要求2所述的一种高精度的弹道修正方法，其特征在于：还包括基于噪声估计器和渐消因子完成弹道修正的过程，该过程包括：

S03-2，利用渐消因子修正滤波结果：

S03-2-1，计算渐消因子，令渐消因子

${\mathrm{\&lambda;}}_k=\frac{Tr\left(N_k\right)}{Tr\left(M_k\right)}=\left\{\begin{array}{cc}{\mathrm{\&lambda;}}_k,& {\mathrm{\&lambda;}}_k\&GreaterEqual;1\\ 1,& {\mathrm{\&lambda;}}_k<1\end{array}\right.---\left(22\right)$

式(22)中，Tr表示求矩阵的迹； $N_k=b_k-H\left({\hat{X}}_{k/k+1}\right)Q_k^{*}{H}^T\left({\hat{X}}_{k/k+1}\right)-l_k{R}_k;$

$M_k=H\left({\hat{X}}_{k/k+1}\right)\mathrm{\&phi;}\left({\hat{X}}_{k/k+1}\right){\hat{P}}_{k+1/k}{\mathrm{\&phi;}}^T\left({\hat{X}}_{k/k+1}\right)H^T\left({\hat{X}}_{k/k+1}\right);l_k=1-d_k;$

S03-2-2，进行一步预测结果修正，

计算带渐消因子的一步预测的误差协方差矩阵

$P_{k+1/k}={\mathrm{\&lambda;}}_k\mathrm{\&phi;}\left({\hat{X}}_k^{*}\right)P_k{\mathrm{\&phi;}}^T\left({\hat{X}}_k^{*}\right)+Q_k^{*}---\left(23\right)$

计算增益矩阵

$K_{k+1}=P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right){\&lsqb;B\left({\hat{X}}_{k+1/k}\right)P_{k+1/k}{B}^T\left({\hat{X}}_{k+1/k}\right)+R_{k+1}\&rsqb;}^{-1}---\left(24\right)$

计算滤波估计值

${\hat{X}}_{k+1}={\hat{X}}_{k+1/k}+K_{k+1}\&lsqb;Z_{k+1}-H\left({\hat{X}}_{k+1/k}\right)\&rsqb;---\left(25\right)$

计算滤波误差方差阵

$P_{k+1}=\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;P_{k+1/k}{\&lsqb;I-K_{k+1}B\left({\hat{X}}_{k+1/k}\right)\&rsqb;}^H+K_{k+1}{R}_{k+1}{K}_{k+1}^H---\left(26\right)$

S03-3，令k+1→k，重复S03-1到S03-3步骤，完成弹道修正过程。

4.根据权利要求2所述的一种高精度的弹道修正方法，其特征在于：所述迭代计算的次数为三次。