一种结合最小二乘法和梯度下降法的机动目标轨迹实时预测方法

1.本发明涉及航空航天领域，具体来说，涉及一种复杂非线性函数参数优化的方法和飞行轨迹动态预测的方案。


背景技术：

2.对空中机动目标未来时刻轨迹进行高精度预测，将能更好地满足空中对抗决策制定以及对目标进行高效拦截的需求。对空中目标进行轨迹预测，需要目标的历史轨迹数据，并采用一定的方法获取目标未来时刻的位置信息。
3.目前机动目标轨迹的预测研究主要集中于基于数据驱动的方法和行为驱动的方法。基于数据驱动的移动目标轨迹预测方法主要通过海量历史轨迹数据来挖掘数据背后隐藏的移动目标行为特征，再与当前位置数据进行融合匹配，进而对目标的运动趋势进行预判。基于数据驱动的轨迹预测方法，预测效果取决于历史数据的质量和训练效果且计算开销大。
4.基于行为驱动的移动目标轨迹预测方法能够根据移动目标的相关运动特征对未来某时间段内的轨迹展开预测，一般分为动力学模型和意图识别两大类。基于动力学模型的移动目标轨迹预测方法考虑移动目标真实的运动特征，能够较好地预测未来运动情况。由于该模型计算量较大，预测的准确性又依赖于理想的环境和状态假设，在实际的预测任务中难以达到较好的预测精度。当前针对一些特定场景，基于意图识别的移动目标轨迹预测方法在短时预测中表现出了较好的效果，但在目标意图不确定的情况下，预测轨迹存在偏差较大的问题。
5.与上述两种方法相比，基于曲线拟合轨迹预测方法具有所需数据少，计算量小，响应速度快等优点。但是，传统的曲线拟合方法仅采用多项式作为拟合函数，存在拟合效果差，预测精度低等问题。本发明在采用多项式的基础上引入超越函数，能够有效提高拟合效果，进而提高轨迹预测精度。另外，针对利用复杂函数拟合数据序列的参数优化问题，本发明提出了基于最小二乘法和梯度下降法相结合的新型参数优化方法，能解决利用复杂函数对轨迹序列进行实时高精度拟合的难题，从而保证在线预测飞行轨迹。


技术实现要素：

6.本发明方法根据空中机动目标的实际运动轨迹具有连续时变过程的特征，基于飞行轨迹曲线拟合的思想，采用包含多项式函数项和超越函数项的复杂非线性函数对机动目标的飞行轨迹数据序列进行拟合，从而用该函数外推目标未来时刻的飞行轨迹。针对复杂非线性函数的参数优化问题，设计了一种参数优化方法，目的是快速对函数模型中的参数进行优化，根据所设计的参数优化方法，设计动态预测方案，实现对目标飞行轨迹的实时动态预测。
7.基于上述思想和实现对机动目标轨迹的动态预测，本发明的技术方案主要包括以
下步骤:
8.1)获得当前时刻t＝tm的数据和前k-1个飞行轨迹数据序列{p
m-k 1
,p
m-k 2
,...,pm}和对应的时间戳序列{t
m-k 1
,t
m-k 2
,...,tm}；采用包含多项式函数和超越函数的复杂非线性函数拟合目标飞行轨迹数据序列，将轨迹预测问题转换为复杂非线性函数的参数优化问题；
9.2)进行数据预处理：以获取到的第一个历史轨迹点为基准，将{p
m-k 1
,p
m-k 2
,...,pm}与p
m-k 1
做差分，获得差分后的飞行轨迹序列{0,p
m-k 2|m-k 1
,...,p
m|m-k 1
}；
10.3)基于最小二乘方法，求解多项式函数的参数；
11.4)根据已求多项式函数参数，基于梯度下降法进一步求解超越函数项参数，获得优化参数后的复杂非线性函数模型f
x
(t)、fy(t)和fz(t)；
12.5)取未来时间戳序列{t
m 1
,t
m 2
,...,t
m n
},代入优化参数后的函数f
x
(t)、fy(t)和fz(t)，外推输出增量序列{p
m 1|m-k 1
,p
m 2|m-k 1
,...,p
m n|m-k 1
}；
13.6)进行数据后处理：将增量序列与第一个历史轨迹点p
m-k 1
相加，输出预测轨迹点序列
14.7)若有新数据p
m 1
输入，用新数据p
m 1
代替预测值并转到步骤(1)，否则结束预测。
15.作为本发明的优选方案，所述步骤1)中采用包含多项式函数和超越函数的复杂非线性函数拟合目标的飞行轨迹数据(p1,p2,p3,...,pk)，飞行轨迹数据序列对应时间戳序列(t1,t2,t3,...,tk),其中pi＝{xi,yi,zi},i＝1,2,...,k，包含目标位置分量信息；
16.包含多项式函数和超越函数的复杂非线性函数f
x
(t)、fy(t)和fz(t)的具体形式为：
[0017][0018]
其中f
x
(t)、fy(t)和fz(t)分别拟合x、y和z位置分量数据；f
x1
(t)、f
y1
(t)和f
z1
(t)为n阶多项式函数，f
x2
(t)、f
y2
(t)和f
z2
(t)为超越函数；
[0019]
为保证函数能够精确拟合目标飞行轨迹，复杂非线性函数的参数优化问题定义如下：
[0020][0021]
其中(a0,a1,a2,...,a
n q
)、(b0,b1,b2,...,b
n q
)和(c0,c1,c2,...,c
n q
)为待优化参数，t为时间自变量，多项式函数的待优化参数为：(a0,...,an)、(b0,...,bn)和(c0,...,cn)，超越函数待优化参数为(a
n 1
,...,a
n q
)、(bn,...,b
n q
)和(cn,...,c
n q
)。
[0022]
作为本发明的优选方案，所述步骤3)中，基于最小二乘方法，求解多项式参数的具
体步骤如下：
[0023]
用n阶多项式函数对x位置分量飞行轨迹序列{x1,x2,x3,...,xk}进行拟合，n阶多项式模型为：
[0024]fx1
(t)＝a0 a1t a2t2 ... antnꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0025]
为使该函数在ti处的值f(ti)与轨迹数据xi接近，多项式函数在时间节点ti处的残差平方为：
[0026][0027]
根据最小二乘准则以及多元函数求极值的方法，得到如下公式：
[0028][0029]
基于最小二乘方法，结合公式(4),可以获得如下方程：
[0030][0031]
进一步根据所获目标飞行轨迹数据序列，基于最小二乘方法构建方程组：
[0032]
axa＝b
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0033]
其中a、xa和b的具体形式为：
[0034][0035]
根据公式(5)至公式(7),参数矩阵xa的最小二乘解为：
[0036]
xa＝(a
t
a)-1atbꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0037]
根据公式(8)即可求得参数(a0,...,an)，多项式函数f
y1
(t)和f
z1
(t)的参数(b0,...,bn)和(c0,...,cn)的求解与(a0,...,an)求解步骤一致。
[0038]
作为本发明的优选方案，所述4)中，基于梯度下降法求解超越函数项参数的具体步骤如下：
[0039]
根据所述s1中的复杂非线性函数的具体形式，设置需要优化的目标函数为：
[0040][0041]
根据目标函数形式，目标函数l1的梯度为：
[0042][0043]
根据已求多项式函数的参数，设定参数矩阵a＝[a0,...,an,a
n 1
,...,a
n q
]
t
，初始
参数矩阵为：
[0044][0045]
其中为由所述s2根据最小二乘法求解的多项式函数项参数值，o1×q代表1
×
q的零矩阵；
[0046]
根据所述参数矩阵的各参数初始值设定，参数矩阵a的更新的迭代公式为：
[0047][0048]
其中a(i)代表第i次迭代的参数矩阵h＝j
t
j，i为单位矩阵，μ是惩罚因子，采用系数ρ的值来控制μ，系数ρ的求解公式如下：
[0049][0050]
其中h的具体形式为：
[0051]
根据ρ的值对参数矩阵和惩罚因子进行更新；根据公式(9)至公式(12),对函数f
x
(t)的参数(a0,...,an,a
n 1
,...,a
n q
)进行优化求解，fy(t)和fz(t)的参数(b0,...,bn,b
n 1
,...,b
n q
)和(c0,...,cn,c
n 1
,...,c
n q
)的优化求解方法同上述步骤；
[0052]
根据所述参数优化过程可得优化参数后的函数f
x
(t)、fy(t)和fz(t)，进一步将时间戳序列{t
k 1
,t
k 2
,...,t
k n
}分别代入f
x
(t)、fy(t)和fz(t)中，预测目标未来时刻飞行轨迹
[0053]
作为本发明的优选方案，所述的根据ρ的值对参数矩阵和惩罚因子进行更新，具体为：
[0054]
如果ρ＞0,则且惩罚因子更新如下：
[0055][0056]
如果ρ＜0,则a(i)＝a
(i-1)
且惩罚因子μ(i)＝μ
(i-1)
·
v，v＝2。
[0057]
与现有技术相比，本发明所具有的有益效果包括：
[0058]
1)基于曲线拟合思想，针对利用复杂函数拟合飞行轨迹数据序列的问题，在采用多项式的基础上引入超越函数，能够有效提高拟合效果，进而提高轨迹预测精度；
[0059]
2)针对利用复杂函数拟合数据序列的参数优化问题，本发明提出了基于最小二乘法和梯度下降法相结合的新型参数优化方法，能解决利用复杂函数对飞行轨迹数据序列进行实时高精度拟合的难题，从而保证在线预测飞行轨迹。
附图说明
[0060]
图1为本发明所提出的方案原理示意图；
[0061]
图2为是优化参数迭代过程的梯度下降法求解参数的流程图；
[0062]
图3是动态的轨迹预测的流程图；
[0063]
图4a是预测盘旋机动目标的效果图；
[0064]
图4b是预测盘旋机动目标的误差图；
[0065]
图5a是预测l机动的效果图；
[0066]
图5b是预测l机动的误差图；
[0067]
图6是预测s机动的效果图及误差图。
具体实施方式
[0068]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0069]
图1为本发明所提出的方案原理示意图；本实施例提供的面向机动目标运动轨迹的动态预测方法，具体建立步骤为：
[0070]
步骤一.针对机动目标轨迹的预测问题，本实例采用结合二阶多项式函数项和正弦函数项的复杂非线性函数来拟合目标飞行轨迹，将轨迹预测问题转换为复杂非线性函数的参数求解问题；
[0071]
步骤二.根据最小二乘方法，求解二阶多项式项的参数；
[0072]
步骤三.根据以求得的二阶多项式参数，基于梯度下降法求解正弦函数项参数；
[0073]
步骤四.根据建立的机动目标轨迹的动态预测流程，对机动目标的飞行轨迹进行动态预测。
[0074]
其中，步骤一的具体实施过程如下：
[0075]
机动目标的历史飞行轨迹序列(p1,p2,p3,...,pk),pi＝{mi,ni,li},i＝1,2,...,k，包含目标在三维笛卡尔空间坐标系下的位置信息，进行数据预处理：p
i|1
＝p
i-p1,i＝1,2,...,k，获得增量序列{0,p
i|1
,...,p
k|1
}，p
i|1
＝{xi,yi,zi},i＝1,2,...,k；
[0076]
首先确定多项式函数项的阶数，由于受到各种条件的限制，实际的运动目标在当前时刻的时间邻域内加速度可以近似认为不变，这样目标的位置和速度的平方相关，同时考虑到目标存在周期性运动，因此本实例选择二阶多项式函数项和正弦函数项组成的非线性函数来拟合观测到的序列数据。
[0077]
根据上述，确定拟合增量序列所用的复杂非线性函数f
x
(t)、fy(t)和fz(t)具体形式如下：
[0078][0079]
其中(a0,a1,...,a5)、(b0,b1,...,b5)和(c0,c1,...,c5)为待求解参数，t为时间自变量。其中二阶多项式函数项的参数为：(a0,a1,a2)、(b0,b1,b2)和(c0,c1,c2)，正弦函数项的参数为(a3,a4,a5)、(b3,b4,b5)和(c3,c4,c5)。
[0080]
步骤二的具体实施过程如下：
[0081]
二阶多项式函数部分的形式为：
[0082][0083]
为使函数f
x1
(t)在ti处的值f
x1
(ti)与历史数据xi接近，定义函数模型在时间节点ti处的残差为：
[0084][0085]
为使f(ti)能反映给定数据的变化趋势，需要残差δ尽可能小。运用最小二乘准则以及多元函数求极值的方法，即：
[0086][0087]
上式的矩阵形式为：
[0088][0089]
根据上述原理，构建线性方程组：
[0090]
axa＝b
[0091]
其中a、xa和b的具体形式为：
[0092][0093]
根据公式(7)参数矩阵xa的最小二乘解为：
[0094]
xa＝(a
t
a)-1atb[0095]
根据公式(8)即可求得公式(1)中函数模型f
x1
(t)的参数(a0,a1,a2)，同理f
y1
(t)和f
z1
(t)的参数(b0,b1,b2)和(c0,c1,c2)的求解步骤与(a0,a1,a2)求解步骤一致。
[0096]
公式(1)中所包含的正弦项参数(a3,a4,a5)，(b3,b4,b5)，(c3,c4,c5)的值通过构造非线性方程组无法直接求解，转步骤三。
[0097]
步骤三，根据已求的二阶多项式函数参数，基于梯度下降法对正弦函数的参数进行求解，图2为本发明实施例提供的基于梯度下降法对超越函数项参数优化的流程图，具体步骤如下：
[0098]
设置目标函数的形式为：
[0099][0100][0101][0102]
对l1进行梯度值计算，l1的梯度表达式为：
[0103][0104]
根据步骤二计算所得二阶多项式参数，设置初始参数矩阵其中为所述步骤二计算得出的多项式参数，o1×3代表1
×
3的零矩阵。
[0105]
基于梯度下降法，利用参数迭代公式计算(a3,a4,a5)，迭代公式为：
[0106][0107]
h＝j
tj[0108][0109]
其中h为hessian矩阵，μ是惩罚因子，第一次迭代初始初始μ
(0)
＝τ
·
max{h0}，max{h0}代表取h0矩阵中最大的元素值。迭代过程中τ值设定为10-6
，用系数ρ来控制μ：
[0110][0111]
其中
[0112]
根据ρ的值对迭代步长和惩罚因子进行更新，具体策略如下：
[0113]
如果ρ＞0,则且惩罚因子更新如下：
[0114][0115]
如果ρ＜0,则α(i)＝α
(i-1)
且惩罚因子μ＝μ
·
v，v＝2。
[0116]
由以上步骤可求得(a3,a4,a5)的值，函数fy的参数(b3,b4,b5)和函数fz的参数(c3,c4,c5)的优化求解方法同参数(a3,a4,a5)的优化求解步骤一致。
[0117]
获取全部参数后，将{t
k 1
,t
k 2
,...,t
k n
}代入函数f
x
(t)、fy(t)和fz(t)得增量序列{p
k 1|1
,p
k 2|1
,...,p
k n|1
}，数据后处理：j＝k 1,k 2,...,k n，输出预测轨迹序列
[0118]
步骤四：采用建立针对空中机动目标轨迹的动态预测流程，完成飞行轨迹的动态预测，图3为本发明实施例提供的针对空中机动目标的轨迹进行动态预测的流程图，具体步骤如下：
[0119]
(1)获得当前时刻t＝tm的数据和前k-1个历史数据，组成飞行轨迹数据序列{p
m-k 1
,p
m-k 2
,...,pm}和时间戳序列{t
m-k 1
,t
m-k 2
,...,tm}；
[0120]
(2)进行数据预处理，获取到的第一个历史轨迹点为基准，将{p
m-k 1
,p
m-k 2
,...,pm}与p
m-k 1
做差分，获得差分后的序列
[0121]
(3)根据最小二乘法，求解获得二次多项式函数项的参数；
[0122]
(4)根据已求二阶多项式函数项参数，基于梯度下降法进一步求解正弦函数项参数，获得复杂非线性函数模型f
x
(t)、fy(t)和fz(t)；
[0123]
(5)取未来时间戳序列{t
m 1
,t
m 2
,...,t
m n
},代入f
x
(t)、fy(t)和fz(t)，输出增量序列
[0124]
(6)进行数据后处理，将增量序列与第一个历史轨迹点相加，输出预测轨迹点序列
[0125]
(7)若有新的定位数据p
m 1
到达，用测量值代替预测值，并转到步骤(1)，否则结束预测。
[0126]
为衡量本发明方法的预测效果，定义预测的精度指标为：
[0127][0128]ei
表示目标未来时刻预测轨迹点和目标未来时刻实际轨迹点的欧几里得距离。
[0129]
基于amd ryzen 9 5900hs@3.30ghz处理器、64位操作系统、40.0gb内存的计算平台上的matlab r2020b平台，对本发明技术方案进行仿真实验，空中机动目标的轨迹数据序列由飞行控制仿真软件给出。
[0130]
对做s机动，桶滚机动以及l机动和盘旋机动的目标分别进行预测实验，在本实例中目标始终以160m/s的恒定速率飞行。仿真环境设置数据的采样间隔时间t为0.05秒，对s机动、l机动以及盘旋机动进行预测时滑动窗口长度设置为40，时间戳序列为{0.05：0.05：2}，预测时间为10s，未来时间戳序列为{2.05：0.05：12}，预测未来时刻数据个数为100，因桶滚机动属于大过载短时机动动作，因此对桶滚机动进行预测时，滑动窗口截取的数据长度设置为20，时间戳序列为{0.05：0.05：1}，预测时间为10s，未来时间戳序列为{2.05：0.05：11}。本实例中的梯度下降法求解参数的目标函数阈值设置为5，参数更新的最大迭代次数设置为50次。
[0131]
对空中机动目标轨迹的预测效果如图4a至图6。在图4a中空中目标做盘旋机动，预测出的目标未来时刻轨迹紧贴实际轨迹，结合图4b误差图可以看到10s后的误差在150m以内；在图5a中空中目标l机动，结合图5b误差图可以看到10s后的误差在450m以内，这是由于l机动属于大过载机动，其运动量突变较大；在图6中空中目标做s机动，选取s机动段上两转弯点进行预测，可以看到10s后的误差在100m以内；实例展示了本发明提出的方法具有较高的预测精度。
[0132]
对不同机动动作预测所消耗的时间进行统计分析，其结果表如下：
[0133][0134]
对二阶多项式参数和正弦项参数求解所消耗的计算时间进行统计分析，单次预测的计算耗时在0.3s～0.35s之间。本发明提出用最小二乘法求解二阶多项式函数部分的参数，进一步基于梯度下降法对正弦函数项参数进行求解，可以快速求解所提复杂非线性函数模型的全部参数，再通过函数模型外推获得目标未来时刻轨迹，进一步再结合所提的动态预测方案，可动态实时预测目标的飞行轨迹。本发明的方法通过较少的计算开销就可精确预测出空中机动目标未来短时的飞行轨迹，具有良好的实时性，满足实际工程需求。
[0135]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。