一种基于方位和频率的单站水下机动目标运动分析方法与流程

1.本发明属于机动目标运动分析技术领域，具体涉及一种基于方位和频率的单站水下机动目标运动分析方法。


背景技术：

2.目标运动分析(tma)是对未知目标的运动参数(航向、航速、航迹等)进行解算，以达到对目标定位跟踪的目的。水下目标运动分析主要分为主动与被动两种方法，相较于主动目标运动分析，被动目标运动分析的观测站无需向被探测的目标发射信号，只需要通过目标辐射的信号或者噪声信息来预测和估计目标运动参数，大大提高了系统的实用性与隐蔽性，更能满足现代作战防御需求。
3.在水下被动目标运动分析中，纯方位(bot)问题受到了广泛的研究。目前纯方位tma主要面临以下难点：一，由于在海洋环境中存在着大量的噪声以及干扰，水下目标的状态参数和定位观测参数之间基本上会变成非线性关系，无法继续利用传统的线性知识去求解，存在较大的研究困难；二，研究表明，纯方位tma系统是一个不完全可观测系统，无法唯一得出确定的目标运动轨迹(nardone s c,aidala v j.observability criteria for bearings-only target motion analysis[j].ieee trans on electron syst,1981,17(2):162～166.)。若要形成完全可观测系统，需观测站做特定的机动。但现实中观测站可能无法完成机动(如无机动能力的小平台，浮标等)，或者观测站很难在短时间内完成特定的机动使tma系统完全可观测。为了解决这一问题，学者们将频率信息引入，与方位联合观测形成方位频率tma系统，经研究验证，该系统是完全可观测系统(胡友峰,孙进才,朱维杰.基于方位频率测量的水下被动目标运动分析及其特性研究[j].西北工业大学学报,2001,(04):537-42.)。
[0004]
方位频率tma系统通过解形式大致可分为以下三种：
[0005]
1.解析解：解析解的代表算法为传统最小二乘算法。文献(rosenqvist,p.a.passive doppler-bearing tracking using a pseudo-linear estimator[j].ieee journal of oceanic engineering,1995,20(2):114-8.以及胡青,宫先仪.方位/频率目标运动分析实验研究[j].声学学报,2005,(02):120-4.)利用了伪线性最小二乘算法对匀速直线运动的目标进行了目标运动分析。结果表明，该算法最终可以估计出目标的运动轨迹。但当量测信息逐渐增多时，量测矩阵的可逆性变差，最小二乘算法运行速度变慢，跟踪精度降低。
[0006]
2.迭代解：迭代解的常用准则为最大似然准则。文献(ho k c,chan y t.an unbiased estimator for bearings-only tracking and doppler-bearing tracking[j].ieee transactions on pattern analysis&machine intelligence,2003,32(9):1721-8.)利用了最大似然准则下的高斯牛顿迭代算法对目标进行了目标运动分析。虽然目标跟踪轨迹可以收敛到目标真实轨迹附近，但该算法需要预先输入迭代初值，大多数情况下无法单独使用。常与最小二乘算法结合使用以获得初值，但两算法结合易导致计算量变
大，运算时间增长，难以实时对目标进行跟踪。
[0007]
3.递归解：递归解的代表算法为贝叶斯类滤波算法。文献(杜选民,姚蓝.基于方位-频率及多阵方位的无源目标跟踪性能研究[j].声学学报,2001,(02):127-34.以及m.shasha,s.ning.an iterative extended kalman filter algorithm applying doppler and bearing measurements for underwater passive target tracking[c]//ieee international conference on information technology,big data and artificial intelligence(iciba),.ieee,2020.)利用了扩展卡尔曼滤波(ekf)算法对方位频率tma模型进行了目标运动分析，取得了较好的跟踪结果。但贝叶斯类滤波算法也需预先给定初值，若初值选取不好，算法极易发散，对解的影响较大，跟踪精度降低。
[0008]
但是，以往的相关研究与文献，大多对基于方位频率tma系统下静止、匀速直线运动的目标进行了目标运动分析，不适用于机动目标，无法实现机动目标运动分析，因而具有一定的局限性；而且需要事先给出中心频率f0，或利用量测频率的均值来替代中心频率f0。由于实际中很难确知f0的先验信息，取均值造成的误差较大，严重影响了算法的稳定性。


技术实现要素：

[0009]
本发明的目的是为解决现有目标运动分析方法无法实现机动目标运动分析，以及无法获取到中心频率的先验信息时运动分析方法的稳定性差的问题，而提出的一种基于方位和频率的单站水下机动目标运动分析方法。
[0010]
本发明为解决上述技术问题所采取的技术方案是：一种基于方位和频率的单站水下机动目标运动分析方法，所述方法具体包括以下步骤：
[0011]
步骤一、构建由状态向量与状态方程组成的目标状态模型；
[0012]
步骤二、基于测量方位与频率构建方位频率tma模型；
[0013]
步骤三、根据构建的目标状态模型和方位频率tma模型，利用imm-ekf算法对机动目标进行目标运动分析，输出机动目标运动分析结果。
[0014]
进一步地，所述步骤一的具体过程为：
[0015]
k时刻的状态向量表示为：
[0016]
x(k|k)＝[x
tk-x
ok y
tk-y
ok v
txk-v
oxk v
tyk-v
oyk 1/f
0k
]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0017]
其中，x(k|k)代表k时刻的状态向量，x
tk
代表k时刻目标在x方向的位置坐标，y
tk
代表k时刻目标在y方向的位置坐标，x
ok
代表k时刻观测站在x方向的位置坐标，y
ok
代表k时刻观测站在y方向的位置坐标，v
txk
代表k时刻目标在x方向的速度分量，v
tyk
代表k时刻目标在y方向的速度分量，v
oxk
代表k时刻观测站在x方向的速度分量，v
oyk
代表k时刻观测站在y方向的速度分量，f
0k
代表k时刻目标辐射中心频率；
[0018]
状态方程由式(2)表示：
[0019]
x(k|k)＝f(k)x(k-1|k-1) γv(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0020]
其中，f(k)为k时刻的状态转移矩阵，γ为过程噪声分布矩阵，v(k)为k时刻的过程噪声向量，x(k-1|k-1)代表k-1时刻的状态向量。
[0021]
进一步地，所述方位频率tma模型为：
[0022]
z(k|k)＝h(x(k|k)) r(k)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0023]
其中，z(k|k)为k时刻的量测，r(k)是k时刻的量测噪声协方差矩阵，z(k|k)、h(
·
)
和r(k)满足：
[0024][0025]
其中，为k时刻的量测方位，c代表水下声速，βk为中间变量，为k时刻的量测频率，ε
βk
为k时刻的量测方位误差，ε
fk
为k时刻的量测频率误差。
[0026]
进一步地，所述目标辐射中心频率f
0k
满足：
[0027]
(a
t
a)-1at
b＝[f
0k
(v
txk-v
oxk
)f
0k
(v
tyk-v
oyk
)f
0k
]
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0028]
其中，中间变量中间变量t＝1,2,
…
,k，为t时刻的量测方位，为t时刻的量测频率。
[0029]
进一步地，所述步骤三的具体过程为：
[0030]
步骤三一、假设目标共有n个运动模型，运动模型i到运动模型j的转移概率为p
ij
，其中i,j＝1,2,
…
,n，令k时刻运动模型i的初始目标状态估计为k时刻运动模型i的初始协方差估计为运动模型i在k时刻的概率为μi(k)，利用公式(6)对n个运动模型进行交互，输出交互后运动模型i在k时刻的目标状态估计和协方差估计；
[0031][0032]
其中：为交互后运动模型i在k时刻的目标状态估计，为交互后运动模型i在k时刻的协方差估计，中间变量μ
ij
(k|k)为：
[0033][0034]
步骤三二、利用ekf算法对机动目标进行跟踪；
[0035][0036]
其中，代表交互后运动模型i的目标状态估计的一步预测；
[0037][0038]
其中，代表交互后运动模型i的协方差估计的一步预测，q(k)为k时刻的状态噪声协方差矩阵，上角标t代表转置；
[0039][0040]
其中，zi(k 1|k)为交互后运动模型i量测的一步预测；
[0041]vi
(k 1)＝zi(k 1|k)-z(k 1|k 1)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0042]
其中，z(k 1|k 1)为k 1时刻的量测，vi(k 1)为交互后运动模型i在k 1时刻的滤波残差；
[0043][0044]
其中，si(k 1)为交互后运动模型i在k 1时刻的新息协方差矩阵，hi(k 1)为交互后运动模型i在k 1时刻的雅可比矩阵，r(k 1)是k 1时刻的量测噪声协方差矩阵；
[0045][0046]
其中，上角标-1代表矩阵的逆，ki(k 1)为交互后运动模型i在k 1时刻的增益；
[0047][0048]
其中，为交互后运动模型i在k 1时刻的目标状态估计；
[0049][0050]
其中，为交互后运动模型i在k 1时刻的协方差估计；
[0051]
步骤三三、利用式(11)计算出的vi(k 1)以及式(12)计算出的si(k 1)计算运动模型i在k 1时刻的概率μi(k 1)：
[0052][0053]
其中，的定义与的定义相同，中间变量c和为：
[0054]
[0055]
步骤三四、将式(14)计算出的式(15)计算出的和式(16)计算出的μi(k 1)代入式(18)：
[0056][0057]
其中，为最终输出的在k 1时刻机动目标的目标状态估计结果，为最终输出的在k 1时刻机动目标的协方差估计结果；
[0058]
将式(14)计算出的作为式(6)中的将式(15)计算出的作为式(6)中的并利用μi(k 1)返回步骤三一；
[0059]
步骤三五、重复步骤三一至步骤三四的过程，完成对机动目标的目标运动分析。
[0060]
更进一步地，所述雅可比矩阵hi(k 1)的计算方式为：
[0061][0062]
其中，代表一步预测结果中目标在y方向的位置坐标，y
o(k 1)
代表k 1时刻观测站在y方向的位置坐标，代表一步预测结果中目标在x方向的位置坐标，x
o(k 1)
代表k 1时刻观测站在x方向的位置坐标，为利用状态的一步预测值计算的k 1时刻方位值，为k 1时刻的量测频率。
[0063]
本发明的有益效果是：
[0064]
本发明构建了新的方位频率tma模型，减弱了中心频率对跟踪精度的影响，不需要已知中心频率f0就能对机动目标进行很好的目标运动分析，解决了传统算法需要预先已知中心频率的弊端。同时本发明方法适用于机动目标运动分析，给方位频率下的机动目标运动分析提供了一种新的方法。
附图说明
[0065]
图1为本发明仿真实验涉及的目标位置轨迹示意图；
[0066]
图2a为本发明仿真实验涉及的方位量测示意图；
[0067]
图2b为本发明仿真实验涉及的频率量测示意图；
[0068]
图3a为本发明仿真实验涉及的位置均方根误差曲线图；
[0069]
图3b为本发明仿真实验涉及的速度均方根误差曲线图；
[0070]
图4为本发明仿真实验涉及的算法与传统算法性能的比较结果图；
[0071]
图5为本发明仿真实验涉及的目标跟踪结果图。
具体实施方式
[0072]
本实施方式提供了一种基于方位和频率的单站水下机动目标运动分析方法，在不需要中心频率f0的情况下，结合方位与频率信息就可以对机动目标进行目标运动分析。即针对imm机动模型，本发明利用最小二乘算法估计中心频率f0并结合ekf算法，对机动目标进行目标运动分析。本发明方法具体包括以下步骤：
[0073]
步骤一：构建目标状态模型，输出状态向量与状态方程。
[0074]
由于较难确知f0的先验信息，故将f0作为一个未知的参数来进行估计，k时刻的状态量表示为：
[0075]
xk＝[x
tk-x
ok y
tk-y
ok v
txk-v
oxk v
tyk-v
oyk 1/f
0k
]
t
[0076]
状态方程下式表示：
[0077]
xk＝fx
k-1
γvk[0078]
其中，f为状态转移矩阵，γ为过程噪声分布矩阵，v为过程噪声向量。x
t
和y
t
代表目标在x,y方向的位置坐标，xo和yo代表观测站在x,y方向的位置坐标。v
tx
和v
ty
代表目标在x,y方向的速度分量，v
ox
和v
oy
代表观测站在x,y方向的速度分量。f0代表目标辐射中心频率。
[0079]
步骤二：构建基于测量方位与频率的目标量测模型，输出方位频率tma模型。
[0080]
构建的量测方程可以表示为：
[0081]
zk＝h(xk) rk[0082]
其中：
[0083][0084]f0k
满足于(a
t
a)-1at
b，其中：
[0085][0086]
其中，ε
β
,εf分别为量测方位、频率误差。c代表水下声速。
[0087]
步骤三：构建方位频率tma模型，利用imm-ekf算法对机动目标进行目标运动分析，
输出机动目标运动分析结果。
[0088]
模型初始化，输入交互：
[0089]
假设目标共有n个运动模型，模型i到模型j的转移概率为p
ij
，其中i,j＝1,2,
…
,n。令k时刻模型i的初始目标状态估计与协方差估计分别为模型i在k时刻的概率为μi(k)，其中，i＝1,2,
…
,n。利用下式对n个模型进行交互，输出交互后模型i在k时刻的目标状态和协方差估计。
[0090][0091]
其中：
[0092][0093]
模型滤波：
[0094]
将模型i在k时刻交互后得到的目标状态和协方差估计输入给ekf算法。
[0095]
首先计算出基于方位与频率信息的雅可比矩阵：
[0096][0097]
其次利用ekf算法对目标进行跟踪，输出模型i在k 1时刻的目标状态估计和协方差估计
[0098]
x(k 1|k)＝f(k)x(k|k)
[0099]
p(k 1|k)＝f(k)p(k|k)f
t
(k) q(k)
[0100]
z(k 1|k)＝h[x(k 1|k)]
[0101]
v(k 1)＝z(k 1|k)-z(k 1)
[0102]
s(k 1)＝h(k 1)p(k 1|k)h
t
(k 1) r(k 1)
[0103]
k(k 1)＝p(k 1|k)h
t
(k 1)s-1
(k 1)
[0104]
x(k 1|k 1)＝x(k 1|k) k(k 1)v(k 1)
[0105]
p(k 1|k 1)＝p(k 1|k)-k(k 1)s(k 1)k(k 1)
t
[0106]
模型概率更新：
[0107]
由公式计算得出输入模型i在k 1时刻的滤波残差vi(k 1)与新息协方差矩阵si(k 1)。若量测误差服从高斯分布，则输出模型i的更新概率：
[0108][0109]
其中，
[0110][0111]
模型输出综合：
[0112]
输入模型i在k 1时刻的更新概率μi(k 1)、目标状态估计和协方差估计利用下式计算出k 1时刻目标的状态估计以及协方差估计x(k 1|k 1)，p(k 1|k 1)。
[0113][0114]
输出k 1时刻目标的状态估计后，将模型i的概率μi(k 1)、目标状态估计和协方差估计输入回步骤三中的模型初始化，形成循环。以此逻辑即可完成对机动目标的目标运动分析。
[0115]
仿真实验
[0116]
仿真条件：假设有1个做匀速直线运动(cv)与匀速转弯运动(ct)的机动目标，其初始位置为(-5000,5000)m，初始速度为(10,5)m/s。匀速直线运动、匀速转弯运动过程的状态转移矩阵以及过程噪声分布矩阵分别为：
[0117][0118]
过程噪声向量v为[0.00010.0001]
t
。目标总共运动2000帧，每帧间隔1s。其中，目标在前1000帧进行匀速直线运动，在后1000帧进行匀速转弯运动，目标位置轨迹的示意图如图1所示，角速度ω大小为-0.1π/180。目标辐射中心频率为1000hz。静止观测站坐标为(0,0)m。观测站方位量测误差标准差为1
°
，频率量测误差标准差为0.1hz，观测站测得的方
位量测如图2a所示，观测站测得的频率量测如图2b所示，水下声速为1500m/s。设置模型跳变概率统一为p
ij
＝[0.9，0.1；0.1，0.9]，模型的初始概率均为μ＝[0.5，0.5]。蒙特卡洛仿真100次。在仿真过程中，滤波初始位置为在真实初始位置的基础上赋予一定的随机误差。
[0119]
本发明采用均方根误差来衡量滤波器的性能：
[0120][0121]
其中，mc为蒙特卡洛仿真次数，x(k)为目标真实状态，为目标估计状态，n为总采样帧数。
[0122]
仿真总结：
[0123]
本发明对基于方位频率信息的机动目标进行仿真并给出目标运动分析结果。由图3a和图3b可以看出，利用本发明算法可以在未知中心频率f0的情况下对机动目标进行跟踪，最终均方根误差曲线可以达到收敛状态，位置均方根误差收敛到大约80m。在与传统的ekf算法对比中，传统算法分别被赋予了预估的中心频率f0和准确的中心频率f0，本发明算法仅赋予了预估的中心频率f0。其中，预估的中心频率f0由公式(5)计算得出。图4表明传统算法只有在已知确切的中心频率时滤波器才能够收敛，最终位置均方根误差收敛到大约200m。一旦中心频率不准确，滤波器将很难达到收敛状态，无法对机动目标进行跟踪。通过本发明方法进行跟踪获得的位置轨迹示意图如图5所示，仿真实验验证了本发明的可行性，具有很高的实用价值。
[0124]
本发明的上述算例仅为详细地说明本发明的计算模型和计算流程，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。