CN104021311A - 一种基于Hermite函数约束的数据融合计算方法 - Google Patents

Abstract

本发明提供了一种基于Hermite函数约束的数据融合计算方法，包括：用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得所述原始测元的误差方程；将获得的所述原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对所述联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值；将获得的所述待估参数的向量值代入Hermite拟合函数和测元系统误差模型中，获得目标轨迹坐标参数值、速度参数值以及测元系统误差值。本发明提供了一种基于Hermite函数约束的数据融合计算方法，该方法在保证目标轨迹参数估计精度、计算稳定度和计算效率的同时，简化了数学模型构造。

Description

一种基于Hermite函数约束的数据融合计算方法

技术领域

本发明涉及测量技术领域，尤其涉及一种基于Hermite函数约束的数据融合计算方法。

背景技术

运动目标轨迹测量精度直接影响工程试验鉴定结果，也影响对测量系统和被测系统故障分析和改进设计的支持能力。因此，测量精度一直是测量领域孜孜不倦追求的目标。随着被测运动目标机动性能越来越复杂，有针对性建成的综合组网测控体系，涵盖多种类测量系统，包括光学经纬仪、摄影经纬仪、脉冲雷达、连续波雷达、相控阵雷达和空间定位遥测接收站等。由众多测量系统提供的海量原始测量数据，为深入挖掘信息资源，研究更科学的数据处理方法提供了重要保障。

当前，“信息融合”技术方兴未艾，但截至到目前，该技术相关理论和方法还没有形成统一的完备体系。“运动目标数据融合”作为“信息融合”技术的一个重要分支，在航空航天测控和武器试验靶场一直是一项研究的热点和难点。经典的误差模型最佳弹道估计(EMBET)数据融合方法，认为目标运动轨迹在时序上不具有任何关联性，即在每一个采样时刻进行轨迹参数融合解算，这种方法没有充分利用目标轨迹时序关联性，使设计矩阵庞大，耗费运算资源，数据处理周期长，效率低。在此理论基础上，多有学者提出了基于样条函数约束和基于动力学约束等的数据融合方法，目的都是通过压缩目标轨迹待估参数空间维数，增强联合方程约束力和误差探查的敏感性，提高目标轨迹参数估计精度、计算稳定度和计算效率，并有大量工程应用证明，但数学模型构建比较复杂。

发明内容

本发明的目的在于提供一种基于Hermite函数约束的数据融合计算方法，以解决上述技术问题。

为实现上述目的本发明提供了一种基于Hermite函数约束的数据融合计算方法，包括：

一种基于Hermite函数约束的数据融合计算方法，其特征在于，包括：

用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得该原始测元的误差方程；其中，该目标轨迹参数包括：目标轨迹坐标参数及速度参数；该原始测元包括：位置测元、距离测元、方位角测元、俯仰角测元以及径向速度测元；

将获得的该原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值；

将获得的该待估参数的向量值代入Hermite拟合函数和测元系统误差模型中，获得目标轨迹坐标参数值、速度参数值以及测元系统误差值。

进一步，用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得该原始测元的误差方程，具体包括以下步骤：

步骤1，用三次Hermite函数分别表征该目标轨迹坐标参数及该速度参数，如下式：

$\left\{\begin{array}{c}{x}_{t_i}=f_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{f}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{f}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{f}_{x3}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{f}_{x4}\left(T_{\mathrm{xj}}\right)),\\ z_{t_i}=f_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{f}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{f}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{f}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{f}_{z4}\left(T_{\mathrm{zj}}\right)),\\ y_{t_i}=f_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{f}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{f}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{f}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{f}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{\stackrel{\·}{f}}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{\stackrel{\·}{f}}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{\stackrel{\·}{f}}_3\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{\stackrel{\·}{f}}_{x4}\left(T_{\mathrm{xj}}\right)),\\ {\stackrel{\·}{z}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{\stackrel{\·}{f}}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{\stackrel{\·}{f}}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{\stackrel{\·}{f}}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{\stackrel{\·}{f}}_{z4}\left(T_{\mathrm{zj}}\right)),\\ {\stackrel{\·}{y}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{\stackrel{\·}{f}}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{\stackrel{\·}{f}}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{\stackrel{\·}{f}}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{\stackrel{\·}{f}}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(2\right)$

式(1)、式(2)中：

设相邻节点之间的时间已经归一化，则 $\left\{\begin{array}{c}{f}_1\left(T\right)=1-{3T}^2+{2T}^3\\ f_2\left(T\right)=T{(1-T)}^2\\ f_3\left(T\right)={3T}^2-{2T}^3\\ f_4\left(T\right)=-(1-T)T^2\end{array}\right.,$ $\left\{\begin{array}{c}{\stackrel{\·}{f}}_1\left(T\right)=-6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_2\left(T\right)=-2T(1-T)\\ {\stackrel{\·}{f}}_3\left(T\right)=6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_4\left(T\right)=-2T+{3T}^2\end{array}\right.;$ 其中， $T=\frac{t_i-t_k}{t_{k+1}{-t}_k}\∈\left[\mathrm{0,1}\right],1-T=\frac{t_{k+1}{-t}_i}{t_{k+1}-t_k},$ t_i为实际采样时刻，t_k为节点时刻，t_k≤t_i＜t_k+1，t_k∈T∈(T_x,T_z,T_y)；(T_x,T_z,T_y)和(m_x,m_z,m_y)分别为目标轨迹坐标参数的Hermite拟合函数节点分布和内节点个数；

(β_x,β_z,β_y)为Hermite拟合函数系数，即待估参数，其中，β_x表示(β_x(j×2),β_x(j×2)+1,β_x(j×2)+2,β_x(j×2)+3),j＝0～m_x，β_z表示(β_z(j×2),β_z(j×2)+1,β_z(j×2)+2,β_z(j×2)+3),j＝0～m_z，β_y表示(β_y(j×2),β_y(j×2)+1,β_y(j×2)+2,β_y(j×2)+3),j＝0～m_y；

步骤2，以式(1)、式(2)为基函数表征该原始测元的误差方程，公式如下：

位置测元误差方程：

距离测元误差方程：

方位角测元误差方程：

俯仰角测元误差方程：

径向速度测元误差方程：

其中，式(3)至式(7)中，(v_x,v_z,v_y,v_R,v_A,v_E,v_V)为各测元残差值，(x_i,z_i,y_i,R_i,A_i,E_i,V_i)为各测元原始测量值，为各测元真值，(s_x,s_z,s_y,s_R,s_A,s_E,s_V)为各测元系统误差模型，为初值，(x₀,z₀,y₀)为站址已知坐标， $R_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2+{({\tilde{y}}_i-y_0)}^2},R_{\mathrm{i\Δ}}=\frac{{\tilde{x}}_i-x_0}{R_{i0}}{\tilde{x}}_i+\frac{{\tilde{z}}_i-z_0}{R_{i0}}{\tilde{z}}_i+\frac{{\tilde{y}}_i-y_0}{R_{i0}}{\tilde{y}}_i,$ $A_{i0}=\arcsin \left(\frac{{\tilde{z}}_i{-z}_0}{L_{i0}}\right),L_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2},A_{\mathrm{i\Δ}}=\frac{-({\tilde{z}}_i-z_0)}{L_{i0}^2}{\tilde{x}}_i+\frac{({\tilde{x}}_i-x_0)}{L_{i0}^2}{\tilde{z}}_i,E_{i0}=\mathrm{arctg}\left(\frac{{\tilde{y}}_i-y_0}{L_{i0}}\right),$ $E_{\mathrm{i\Δ}}=\frac{-({\tilde{y}}_i-y_0)({\tilde{x}}_i-x_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{x}}_i+\frac{-({\tilde{y}}_i-y_0)({\tilde{z}}_i-z_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{z}}_i+\frac{L_{i0}}{{\left(R_{i0}\right)}^2}{\tilde{y}}_i,$ α_T为判象限角。

进一步，将获得的该原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值，具体包括以下步骤：

步骤1，在式(2)至式(7)中至少选取四种测元方程式，组成联合观测方程式组，并以下式表示：

V＝AX+BC+L    (8)

其中，该位置测元和方位角测元不能单独选取或两两选取；

式(8)中，V为由各测元残差值(v_x,v_z,v_y,v_R,v_A,v_E,v_V)组成的误差向量；A为Hermite拟合函数系数(β_x,β_z,β_y)表征各测元的设计矩阵，B为系统误差模型系数矩阵，L为各测元常数向量；X为由Hermite拟合函数系数(β_x,β_z,β_y)组成的待估参数向量，C为由各测元系统误差模型系数组成的待估参数向量，X和C为需被解算的待估参数向量；

步骤2，依据最小二乘法原理对式(8)进行解算，待估参数向量解如下式：

[X^T,C^T]^T＝-([A,B]^TP[A,B])^-1[A,B]^TPL    (9)

式(9)中，P为权值矩阵；

其中，在对式(8)进行解算过程中，还包括对式(8)进行迭代计算。

进一步，将获得的该原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值之后，该方法还包括：

对融合解算后的各测元残差值进行分析，依据残差统计特征，重新设计调整联合观测方程组中关于测元系统误差模型部分，再次解算联合观测方程组，直到所有测元残差值均值为零。

本发明与现有技术相比，具有以下有益效果：通过用Hermite函数表征目标轨迹参数，并以此为基函数对原始测元进行表征，得到该原始测元的误差方程；然后，将该误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算。该方法在保证目标轨迹参数估计精度、计算稳定度和计算效率的同时，简化了数学模型构造。

附图说明

图1为本发明一种基于Hermite函数约束的数据融合计算方法第一实施例的流程图；

图2为本发明对联合观测方程组进行再次解算的流程图；

图3为设备站址与理论弹道的空间几何关系图。

具体实施方式

下面结合附图所示的各实施方式对本发明进行详细说明，但应当说明的是，这些实施方式并非对本发明的限制，本领域普通技术人员根据这些实施方式所作的功能、方法、或者结构上的等效变换或替代，均属于本发明的保护范围之内。

参图1所示，图1为本发明一种基于Hermite函数约束的数据融合计算方法的流程图。

本实施例提供了图1为本发明一种基于Hermite函数约束的数据融合计算方法的流程图，包括：

步骤S102，用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得该原始测元的误差方程；其中，该目标轨迹参数包括：目标轨迹坐标参数及速度参数；该原始测元包括：位置测元、距离测元、方位角测元、俯仰角测元以及径向速度测元。

步骤S104，将获得的该原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值。

步骤S106，将获得的该待估参数的向量值代入Hermite拟合函数和测元系统误差模型中，获得目标轨迹坐标参数值、速度参数值以及测元系统误差值。

下面以具体实施例详细说明本发明的具体实施步骤：

假设有多种测量系统同时跟踪测量运动目标轨迹，事后提供7种独立原始测元数据，即(x_i,z_i,y_i,R_i,A_i,E_i,V_i)，利用这7种测元采用基于Hermite函数约束的数据融合计算方法的具体实施步骤如下：

1.用Hermite函数对目标轨迹参数进行表征。

Hermite函数定义如下：设函数H(x)在n个节点x₀＜x₁＜…＜x_n-1上的函数值为y₀,y₁,…,y_n-1，一阶导数值为y′₀,y′₁,…,y′_n-1，则H(x)可以用(2n-1)次Hermite函数近似代替，如下式，

$H\left(x\right)=P_{2n-1}\left(x\right)=\underset{k=0}{\overset{n-1}{\mathrm{\Σ}}}[y_k+\left({x-x}_k\right)(y_k^{\′}-{2y}_k{l}_k^{\′}\left(x\right))]l_k^2\left(x\right)$

式中， $l_k\left(x\right)=\underset{\underset{j\≠k}{j=0}}{\overset{n-1}{\mathrm{\Π}}}\left[\frac{\left({x-x}_j\right)}{(x_k-x_j)}\right],l_k^{\′}\left(x\right)=\underset{\underset{j\≠k}{j=0}}{\overset{n-1}{\mathrm{\Σ}}}\left[\frac{1}{x_k-x_j}\right].$

2.则目标轨迹坐标、速度参数用三次Hermite函数表征，如(1)式、(2)式所示，。

$\left\{\begin{array}{c}{x}_{t_i}=f_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{f}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{f}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{f}_{x3}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{f}_{x4}\left(T_{\mathrm{xj}}\right)),\\ z_{t_i}=f_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{f}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{f}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{f}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{f}_{z4}\left(T_{\mathrm{zj}}\right)),\\ y_{t_i}=f_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{f}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{f}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{f}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{f}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{\stackrel{\·}{f}}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{\stackrel{\·}{f}}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{\stackrel{\·}{f}}_3\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{\stackrel{\·}{f}}_{x4}\left(T_{\mathrm{xj}}\right)),\\ {\stackrel{\·}{z}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{\stackrel{\·}{f}}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{\stackrel{\·}{f}}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{\stackrel{\·}{f}}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{\stackrel{\·}{f}}_{z4}\left(T_{\mathrm{zj}}\right)),\\ {\stackrel{\·}{y}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{\stackrel{\·}{f}}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{\stackrel{\·}{f}}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{\stackrel{\·}{f}}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{\stackrel{\·}{f}}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(2\right)$

设相邻节点之间的时间已经归一化，则 $\left\{\begin{array}{c}{f}_1\left(T\right)=1-{3T}^2+{2T}^3\\ f_2\left(T\right)=T{(1-T)}^2\\ f_3\left(T\right)={3T}^2-{2T}^3\\ f_4\left(T\right)=-(1-T)T^2\end{array}\right.,$ $\left\{\begin{array}{c}{\stackrel{\·}{f}}_1\left(T\right)=-6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_2\left(T\right)=-2T(1-T)\\ {\stackrel{\·}{f}}_3\left(T\right)=6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_4\left(T\right)=-2T+{3T}^2\end{array}\right.;$ 其中， $T=\frac{t_i-t_k}{t_{k+1}{-t}_k}\∈\left[\mathrm{0,1}\right],1-T=\frac{t_{k+1}{-t}_i}{t_{k+1}-t_k},$ t_i为实际采样时刻，t_k为节点时刻，t_k≤t_i＜t_k+1，t_k∈T∈(T_x,T_z,T_y)；(T_x,T_z,T_y)和(m_x,m_z,m_y)分别为目标轨迹坐标参数的Hermite拟合函数节点分布和内节点个数；

其中，为获得合适的节点分布，针对目标轨迹加速度变化特点来选取，即在加速度变化剧烈的时段，采用较密的节点分布，在加速度变化平稳的时段，采用稀疏节点分布。

(β_x,β_z,β_y)为Hermite拟合函数系数，即待估参数，其中，β_x表示(β_x(j×2),β_x(j×2)+1,β_x(j×2)+2,β_x(j×2)+3),j＝0～m_x，β_z表示(β_z(j×2),β_z(j×2)+1,β_z(j×2)+2,β_z(j×2)+3),j＝0～m_z，β_y表示(β_y(j×2),β_y(j×2)+1,β_y(j×2)+2,β_y(j×2)+3),j＝0～m_y。

3.则以(1)式、(2)式为基函数表征各测量系统多种类测元的误差方程公式如下所示。

位置测元误差方程：

距离测元误差方程：

方位角测元误差方程：

俯仰角测元误差方程：

径向速度测元误差方程：

以上式中，(v_x,v_z,v_y,v_R,v_A,v_E,v_V)为各测元残差值，(x_i,z_i,y_i,R_i,A_i,E_i,V_i为各测元原始测量值，为各测元真值，(s_x,s_z,s_y,s_R,s_A,s_E,s_V)为各测元系统误差模型，如常值、线性或非线性函数模型等，为初值，(x₀,z₀,y₀)为站址已知坐标， $R_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2+{({\tilde{y}}_i-y_0)}^2},R_{\mathrm{i\Δ}}=\frac{{\tilde{x}}_i-x_0}{R_{i0}}{\tilde{x}}_i+\frac{{\tilde{z}}_i-z_0}{R_{i0}}{\tilde{z}}_i+\frac{{\tilde{y}}_i-y_0}{R_{i0}}{\tilde{y}}_i,$ $A_{i0}=\arcsin \left(\frac{{\tilde{z}}_i{-z}_0}{L_{i0}}\right),L_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2},A_{\mathrm{i\Δ}}=\frac{-({\tilde{z}}_i-z_0)}{L_{i0}^2}{\tilde{x}}_i+\frac{({\tilde{x}}_i-x_0)}{L_{i0}^2}{\tilde{z}}_i,E_{i0}=\mathrm{arctg}\left(\frac{{\tilde{y}}_i-y_0}{L_{i0}}\right),$ $E_{\mathrm{i\Δ}}=\frac{-({\tilde{y}}_i-y_0)({\tilde{x}}_i-x_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{x}}_i+\frac{-({\tilde{y}}_i-y_0)({\tilde{z}}_i-z_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{z}}_i+\frac{L_{i0}}{{\left(R_{i0}\right)}^2}{\tilde{y}}_i,$ α_T为判象限角。

4.由以上7种测元误差方程组成联合观测方程组。(对于以上7种测元，本实施例不限于由全部测元组成联合观测方程组，也可以任选4种以上测元组成联合观测方程组，其中，位置测元和方位角测元不能单独或两两单独选用)，可写成矩阵形式，如(8)式所示，

V＝AX+BC+L    (8)

式中，V为由各测元残差(v_x,v_z,v_y,v_R,v_A,v_E,v_V)组成的误差向量；A为Hermite拟合函数系数(β_x,β_z,β_y)表征各测元的设计矩阵，B为系统误差模型系数矩阵，L为各测元常数向量；X为由Hermite拟合函数系数(β_x,β_z,β_y)组成的待估参数向量，C为由各测元系统误差模型系数组成的待估参数向量，X和C为需被解算的待估参数向量。

5.确定Hermite拟合函数节点分布。为获得合适的节点分布，可针对目标轨迹加速度变化特点来选取，即在加速度变化剧烈的时段，采用较密的节点分布，在加速度变化平稳的时段，采用稀疏节点分布。在已知Hermite拟合函数节点分布确定的前提下，依据最小二乘法原理对(8)式进行解算，如(9)式所示：

[X^T,C^T]^T＝-([A,B]^TP[A,B])^-1[A,B]^TPL    (9)

式中，P＝diag(P_x,P_z,P_y,P_R,P_A,P_E,P_V)为权值矩阵。以7个测元统计随机误差方差比为依据设计为对角矩阵。

6.在实际求解中，由于非线性方程级数展开和初始目标轨迹坐标的近似性共同带来的截断误差，需对(8)式进行迭代计算。

7.为了获得更优的融合解算结果，需要对当前融合解算后的各测元残差值进行分析，依据残差统计特征，重新设计调整联合观测方程组中关于测元系统误差模型部分，再次解算联合观测方程组，直到所有测元残差值均值E(V)为零。(参图2所示)

8.最后，将融合解算获得的参数向量X代入式(1)、式(2)中，即获得目标轨迹坐标参数值和速度参数值；将参数向量C带入测元系统误差模型BC中，即获得各测元系统误差值。

为了验证本发明方法的效果，下面通过具体实验进行说明：

参图3所示，图3为设备站址与理论弹道的空间几何关系图。

当本发明用于运动目标轨迹数据融合时，假设产生一条运动目标理论轨迹，数据采样率为10Hz，以其为依据，设计4台光学经纬仪站址(1站、2站、3站和4站)分别布设在目标飞行轨迹两侧，反推产生方位角和俯仰角共8个真值角度测元；设计1部脉冲雷达布设在目标飞行轨迹前端右侧，反推产生斜距、方位角和俯仰角3个真值测元；设计1部测速雷达布设在目标飞行轨迹前端左侧，反推产生1个真值径向速度测元。

给这12个仿真测元加入相应的随机误差和系统误差，具体数值见表1所示。初值为在三个方向上分别加上300m、500m、200m的固定误差和6m、6m、8m的随机误差。

表1  各测元误差分配值

应用本发明方法对加入随机误差和系统误差的12个仿真测元数据进行目标轨迹参数融合计算，融合解算后获得一组目标坐标数据，与目标理论轨迹真值进行比对，误差统计如果见表2：

表2  融合解算目标轨迹坐标与理论真值比对误差统计

应用本发明方法对参与融合计算的12个测元误差探查结果的统计值见表3。

表3  融合解算获得各测元误差统计结果

从解算结果分析来看，采用本发明方法获得的目标轨迹定位精度在0.5米左右，对原始测元随机误差探查的准确度在94.3％以上，对系统误差探查的准确度在78.3％以上。

本发明提供了一种基于Hermite函数约束的数据融合计算方法，通过用Hermite函数表征目标轨迹参数，并以此为基函数对原始测元进行表征，得到该原始测元的误差方程；然后，将该误差方程组成联合观测方程组，并依据最小二乘准则对该联合观测方程组进行融合解算。该方法在保证目标轨迹参数估计精度、计算稳定度和计算效率的同时，简化了数学模型构造。

上文所列出的一系列的详细说明仅仅是针对本发明的可行性实施方式的具体说明，它们并非用以限制本发明的保护范围，凡未脱离本发明技艺精神所作的等效实施方式或变更均应包含在本发明的保护范围之内。

对于本领域技术人员而言，显然本发明不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

Claims (4)

1.一种基于Hermite函数约束的数据融合计算方法，其特征在于，包括：

用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得所述原始测元的误差方程；其中，所述目标轨迹参数包括：目标轨迹坐标参数及速度参数；所述原始测元包括：位置测元、距离测元、方位角测元、俯仰角测元以及径向速度测元；

将获得的所述原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对所述联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值；

将获得的所述待估参数的向量值代入Hermite拟合函数和测元系统误差模型中，获得目标轨迹坐标参数值、速度参数值以及测元系统误差值。

2.如权利要求1所述的方法，其特征在于，用Hermite函数表征目标轨迹参数，并以此为基函数对测量系统提供的原始测元进行表征，获得所述原始测元的误差方程，具体包括以下步骤：

步骤1，用三次Hermite函数分别表征所述目标轨迹坐标参数及所述速度参数，如下式：

$\left\{\begin{array}{c}{x}_{t_i}=f_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{f}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{f}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{f}_{x3}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{f}_{x4}\left(T_{\mathrm{xj}}\right)),\\ z_{t_i}=f_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{f}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{f}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{f}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{f}_{z4}\left(T_{\mathrm{zj}}\right)),\\ y_{t_i}=f_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{f}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{f}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{f}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{f}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(1\right)$

$\left\{\begin{array}{c}{\stackrel{\·}{x}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_x,T_x)=\underset{j=0}{\overset{m_x}{\mathrm{\Σ}}}({\mathrm{\β}}_{x(j\×2)}{\stackrel{\·}{f}}_{x1}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+1)}{\stackrel{\·}{f}}_{x2}\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+2)}{\stackrel{\·}{f}}_3\left(T_{\mathrm{xj}}\right)+{\mathrm{\β}}_{x(j\×2+3)}{\stackrel{\·}{f}}_{x4}\left(T_{\mathrm{xj}}\right)),\\ {\stackrel{\·}{z}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_z,T_z)=\underset{j=0}{\overset{m_z}{\mathrm{\Σ}}}({\mathrm{\β}}_{z(j\×2)}{\stackrel{\·}{f}}_{z1}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+1)}{\stackrel{\·}{f}}_{z2}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+2)}{\stackrel{\·}{f}}_{z3}\left(T_{\mathrm{zj}}\right)+{\mathrm{\β}}_{z(j\×2+3)}{\stackrel{\·}{f}}_{z4}\left(T_{\mathrm{zj}}\right)),\\ {\stackrel{\·}{y}}_{t_i}={\stackrel{\·}{f}}_3(t_i,{\mathrm{\β}}_y,T_y)=\underset{j=0}{\overset{m_y}{\mathrm{\Σ}}}({\mathrm{\β}}_{y(j\×2)}{\stackrel{\·}{f}}_{y1}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+1)}{\stackrel{\·}{f}}_{y2}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+2)}{\stackrel{\·}{f}}_{y3}\left(T_{\mathrm{yj}}\right)+{\mathrm{\β}}_{y(j\×2+3)}{\stackrel{\·}{f}}_{y4}\left(T_{\mathrm{yj}}\right)),\end{array}\right.---\left(2\right)$

式(1)、式(2)中：

设相邻节点之间的时间已经归一化，则 $\left\{\begin{array}{c}{f}_1\left(T\right)=1-{3T}^2+{2T}^3\\ f_2\left(T\right)=T{(1-T)}^2\\ f_3\left(T\right)={3T}^2-{2T}^3\\ f_4\left(T\right)=-(1-T)T^2\end{array}\right.,$ $\left\{\begin{array}{c}{\stackrel{\·}{f}}_1\left(T\right)=-6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_2\left(T\right)=-2T(1-T)\\ {\stackrel{\·}{f}}_3\left(T\right)=6\left({T-T}^2\right)\\ {\stackrel{\·}{f}}_4\left(T\right)=-2T+{3T}^2\end{array}\right.;$ 其中， $T=\frac{t_i-t_k}{t_{k+1}{-t}_k}\∈\left[\mathrm{0,1}\right],1-T=\frac{t_{k+1}{-t}_i}{t_{k+1}-t_k},$ t_i为实际采样时刻，t_k为节点时刻，t_k≤t_i＜t_k+1，t_k∈T∈(T_x,T_z,T_y)；(T_x,T_z,T_y)和(m_x,m_z,m_y)分别为目标轨迹坐标参数的Hermite拟合函数节点分布和内节点个数；

(β_x,β_z,β_y)为Hermite拟合函数系数，即待估参数，其中，β_x表示(β_x(j×2),β_x(j×2)+1,β_x(j×2)+2,β_x(j×2)+3),j＝0～m_x，β_z表示(β_z(j×2),β_z(j×2)+1,β_z(j×2)+2,β_z(j×2)+3),j＝0～m_z，β_y表示(β_y(j×2),β_y(j×2)+1,β_y(j×2)+2,β_y(j×2)+3),j＝0～m_y；

步骤2，以式(1)、式(2)为基函数表征所述原始测元的误差方程，公式如下：

位置测元误差方程：

距离测元误差方程：

方位角测元误差方程：

俯仰角测元误差方程：

径向速度测元误差方程：

其中，式(3)至式(7)中，(v_x,v_z,v_y,v_R,v_A,v_E,v_V)为各测元残差值，(x_i,z_i,y_i,R_i,A_i,E_i,V_i)为各测元原始测量值，为各测元真值，(s_x,s_z,s_y,s_R,s_A,s_E,s_V)为各测元系统误差模型，为初值，(x₀,z₀,y₀)为站址已知坐标， $R_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2+{({\tilde{y}}_i-y_0)}^2},R_{\mathrm{i\Δ}}=\frac{{\tilde{x}}_i-x_0}{R_{i0}}{\tilde{x}}_i+\frac{{\tilde{z}}_i-z_0}{R_{i0}}{\tilde{z}}_i+\frac{{\tilde{y}}_i-y_0}{R_{i0}}{\tilde{y}}_i,$ $A_{i0}=\arcsin \left(\frac{{\tilde{z}}_i{-z}_0}{L_{i0}}\right),L_{i0}=\sqrt{{({\tilde{x}}_i-x_0)}^2+{({\tilde{z}}_i-z_0)}^2},A_{\mathrm{i\Δ}}=\frac{-({\tilde{z}}_i-z_0)}{L_{i0}^2}{\tilde{x}}_i+\frac{({\tilde{x}}_i-x_0)}{L_{i0}^2}{\tilde{z}}_i,E_{i0}=\mathrm{arctg}\left(\frac{{\tilde{y}}_i-y_0}{L_{i0}}\right),$ $E_{\mathrm{i\Δ}}=\frac{-({\tilde{y}}_i-y_0)({\tilde{x}}_i-x_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{x}}_i+\frac{-({\tilde{y}}_i-y_0)({\tilde{z}}_i-z_0)}{{\left(R_{i0}\right)}^2{L}_{i0}}{\tilde{z}}_i+\frac{L_{i0}}{{\left(R_{i0}\right)}^2}{\tilde{y}}_i,$ α_T为判象限角。

3.如权利要求2所述的方法，其特征在于，将获得的所述原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对所述联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值，具体包括以下步骤：

步骤1，在式(2)至式(7)中至少选取四种测元方程式，组成联合观测方程式组，并以下式表示：

V＝AX+BC+L    (8)

其中，所述位置测元和方位角测元不能单独选取或两两选取；

式(8)中，V为由各测元残差值(v_x,v_z,v_y,v_R,v_A,v_E,v_V)组成的误差向量；A为Hermite拟合函数系数(β_x,β_z,β_y)表征各测元的设计矩阵，B为系统误差模型系数矩阵，L为各测元常数向量；X为由Hermite拟合函数系数(β_x,β_z,β_y)组成的待估参数向量，C为由各测元系统误差模型系数组成的待估参数向量，X和C为需被解算的待估参数向量；

步骤2，依据最小二乘法原理对式(8)进行解算，待估参数向量解如下式：

[X^T,C^T]^T＝-([A,B]^TP[A,B])^-1[A,B]^TPL    (9)

式(9)中，P为权值矩阵；

其中，在对式(8)进行解算过程中，还包括对式(8)进行迭代计算。

4.如权利要求3所述的方法，其特征在于，将获得的所述原始测元的误差方程组成联合观测方程组，并依据最小二乘准则对所述联合观测方程组进行融合解算，获得以Hermite拟合函数系数和测元系统误差模型系数为代估参数的向量值之后，该方法还包括：

对融合解算后的所述各测元残差值进行分析，依据残差统计特征，重新设计调整所述联合观测方程组中关于测元系统误差模型部分，再次解算所述联合观测方程组，直到所有测元残差值均值为零。