CN103336863A - 基于雷达飞行航迹观测数据的飞行意图识别方法 - Google Patents

Abstract

本发明公开了一种基于雷达飞行航迹观测数据的飞行意图识别方法，其中包括：根据雷达航迹位置观测数据，计算航向角、飞行速度、爬升率，建立飞行航迹样本库，作为训练意图识别模型的基础数据；根据飞行计划信息、航路点位置数据，建立典型飞行意图模型，并标注飞行航迹样本的意图类别；根据隐马尔科夫模型（Hidden Markov Model，HMM）原理，建立几种典型飞行意图识别模型，采用期望最大学习算法，训练识别模型的参数；采用飞行意图识别模型，根据前向算法，计算飞行航迹样本当前时刻的飞行意图，采用滚动时间窗，对一段时间内的局部飞行意图进行加权求和，得到全局飞行意图。

Description

基于雷达飞行航迹观测数据的飞行意图识别方法

技术领域

本发明涉及空中交通管理技术，尤其涉及一种基于雷达飞行航迹观测数据的飞行意图识别方法。

背景技术

空中交通管理技术的主要目标是保障空中交通能够安全、高效、有序地在空域中运行。在空中交通管理领域，空中交通态势感知能力（Air Traffic Situational Awareness，ATSA）是保障飞行安全、提高空中交通运行效率的基础技术。根据空中交通态势的时间属性，ATSA可以分为两个主要阶段：（1）当前时刻飞行器的飞行航迹运动状态感知（2）飞行器未来一段时间的飞行意图识别。其中，当前时刻的飞行航迹运动状态可以直接从空中交通监视系统中获取，例如，雷达监视系统、自动相关监视系统、多点定位系统等。飞行航迹运动状态包括经度、纬度、高度等三维位置信息，以及航向角、速度、爬升率等飞行航迹特征。在获取当前时刻飞行航迹运动状态基础上，管制员需要预测未来一段时间内飞行器可能到达的位置，以探测飞行器之间可能存在的飞行冲突，作为安全飞行的决策基础，确保飞行器保持安全飞行间隔。因此，识别未来一段时间飞行器可能到达的位置，即飞行意图识别，是ATSA的关键技术之一。

现有的飞行意图识别方法主要分为（1）基于飞行员离散操作指令的飞行意图识别（2）基于雷达飞行航迹观测数据的飞行意图识别两类。其中，基于飞行员离散操作指令的飞行意图识别方法是指通过收集飞行员的实时操作指令，采用计划识别（Plan Recognition）技术，推测飞行意图。该方法虽然能够较为准确识别出飞行意图，但是由于从机舱内部采集飞行员实时操作指令数据非常困难，加之空地数据传输的困难，该方法在工程应用中受到很大限制。为了避免数据采集和传输等问题，有学者提出了基于雷达飞行航迹观测数据的飞行意图识别方法，即根据雷达系统获取的飞行器飞行航迹的观测数据，与飞行意图模型进行最优匹配计算，将匹配概率最大的飞行意图，作为飞行意图的识别结果。近年来，随着空中交通流量的快速增长，在主干航路和繁忙飞行终端区空域内空中交通态势呈现高密度、高动态等新特点，为意图识别带来了新挑战。现有的识别方法在进行意图模型匹配时，仅仅考虑单个飞行航迹运动状态，直接进行向量点乘运算，在识别精度方面难以满足高动态空域环境需求，特别是针对存在频繁机动操作的空域环境。

发明内容

本发明的目的是为了解决上述问题，提供一种基于雷达飞行航迹观测数据的飞行意图识别方法，综合考虑多种飞行航迹特征，有效提高飞行意图识别精度。

本发明提供了一种基于雷达飞行航迹观测数据的飞行意图识别方法，包括以下几个步骤：

步骤一：根据雷达航迹位置观测数据，计算航向角、飞行速度、爬升率，建立飞行航迹样本库，作为训练意图识别模型参数的基础数据。

步骤二：根据飞行计划信息、航路点位置数据，建立典型飞行意图模型，标注飞行航迹样本的意图类别；

步骤三：根据隐马尔科夫模型原理，建立飞行意图识别模型，根据期望最大学习算法，训练识别模型的参数；

步骤四：采用飞行意图识别模型，根据前向算法，计算飞行航迹样本当前时刻的飞行意图。采用滚动时间窗，对一段时间的飞行意图进行加权，得到最终飞行意图；

本发明的优点在于：

本发明提供的一种基于雷达飞行航迹观测数据的飞行意图识别方法，通过对雷达航迹历史观测数据进行分析，计算关键航迹特征，建立飞行航迹样本库，并根据飞行计划信息、航路点位置数据等辅助信息，提炼出空中交通管理中的几种典型飞行意图模型，并对飞行航迹样本库中的样本进行类别标注，作为训练飞行意图识别方法具体参数的基础数据。根据隐马尔科夫模型原理，对几种典型飞行意图分别建立识别模型，根据期望最大化学习方法，采用飞行航迹样本集训练出识别模型的参数。根据前向算法，判断未知意图类别的飞行航迹样本当前时刻的飞行意图，并采用滚动时间窗，对时间窗内的飞行意图识别结果进行加权，得到最终的全局飞行意图。

附图说明

图1为本发明基于雷达飞行航迹观测数据的飞行意图识别方法实施例的流程图；

图2水平意图识别场景；

图3垂直意图识别场景；

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于雷达飞行航迹观测数据的飞行意图识别方法，流程如图1所示，具体可以包括如下步骤：

步骤一，根据雷达航迹位置观测数据，计算航向角、飞行速度、爬升率，建立飞行航迹样本库。

本步骤根据雷达航迹三维位置观测数据，提取飞行航迹的航向角、飞行速度、爬升率，并对采集到的航迹样本进行预处理，建立飞行航迹样本库，作为飞行意图识别方法的基础数据。

在本实施例中，飞行航迹数据输入采用WGS-84坐标系下的三维位置观测数据，记为Traj(k)＝{lon(k),lat(k),alt(k)},k＝1,...,N。由于雷达航迹数据的更新率为5秒/次，因此航迹数据可表示为离散形式，其中，lon(k)是k时刻的经度，lat(k)是k时刻的纬度，alt(k)是k时刻的高度。

为了方便提取航迹数据中的关键飞行特征，首先将航迹数据从WGS-84坐标系转换到地心地固直角坐标系（Earth-Centered,Earth-Fixed,ECEF），具体坐标转换如公式（1）所示：

$\left\{\begin{array}{c}x\left(k\right)=(N_e+\mathrm{alt}\left(k\right))\cos \left(\mathrm{lat}\left(k\right)\right)\cos \left(\mathrm{lon}\left(k\right)\right)\\ y\left(k\right)=(N_e+\mathrm{alt}\left(k\right))\cos \left(\mathrm{lat}\left(k\right)\right)\sin \left(\mathrm{lon}\left(k\right)\right)\\ z\left(k\right)=[N_e(1-e^2)+\mathrm{alt}\left(k\right)]\sin \left(\mathrm{lat}\left(k\right)\right)\end{array}\right.---\left(1\right)$

在ECEF坐标系下,原点为地球质心，航迹表示为Traj(k)＝{x(k),y(k),z(k)},k＝1,...,N，其中，x(k)为由原点指向经纬度(0,0)位置的X轴上的坐标点，y(k)为由原点指向90°经线的Y轴上的坐标点，z(k)为由原点向北沿地球自转方向的Z轴上的坐标点。式（1）中，lon(k)是WGS-84坐标系下k时刻的经度，lat(k)是WGS-84坐标系下k时刻的纬度，alt(k)是WGS-84坐标系下k时刻的高度。其他参数的具体含义如下：

是主垂直面的曲率半径；a是地球椭球的长半轴，即地球赤道半径，取6378137.0米；

是地球椭球偏心率，且f＝1-b/a是地球椭球扁平率，其中，b是地球椭球的短半轴，即地球极半径，取6356752.3米。

在ECEF坐标系下，为了更好地描述飞行航迹的运动特征，本实施例根据Traj(k)＝{x(k),y(k),z(k)},k＝1,...,N计算航向角ψ(k)、飞行速度v(k)、爬升率c(k)三个关键特征参数。其中，ψ(k)用来描述飞行航迹在水平平面上转弯机动特征，v(k)用来描述飞行航迹在水平平面上加减速运动特征，c(k)用来描述飞行航迹在垂直平面上高度变化特征，计算公式分别如（2）（3）（4）所示，其中Δt是2倍雷达数据更新率，取10秒。

$\mathrm{\ψ}\left(k\right)=\arctan \frac{y(k+1)-y(k-1)}{x(k+1)-x(k-1)},k=2,...,N-1.---\left(2\right)$

$v\left(k\right)=\frac{y(k+1)-y(k-1)}{x(k+1)-x(k-1)},k=2,...,N-1.---\left(3\right)$

$c\left(k\right)=\frac{z(k+1)-z(k-1)}{\mathrm{\Δt}},k=2,...,N-1.---\left(4\right)$

然后，本实施例中采用低通滤波器对航迹特征滤波，并进行数据归一化处理，以避免真实雷达航迹观测数据中的数据错误，减小雷达数据噪声。低通滤波器如公式（5）所示：

$\tilde{f}\mathrm{eat}\left(1\right)=\mathrm{feat}\left(1\right)$                                 （5）

$\tilde{f}\mathrm{eat}\left(k\right)=\mathrm{\θfeat}\left(k\right)+(1-\mathrm{\θ})\tilde{f}\mathrm{eat}(k-1),k=2,...,N-1.$

其中，feat(k)是k时刻滤波前航迹特征（航向角、飞行速度、爬升率），

是经过低通滤波器后的航迹特征，θ取0.4。

数据归一化过程通过公式（6），将数据映射到[0,1]区间内：

$\hat{f}\mathrm{eat}\left(k\right)=\frac{\tilde{f}\mathrm{eat}\left(k\right)-\tilde{f}\mathrm{ea}{t}_{\min }}{\hat{f}\mathrm{ea}{t}_{\max }-\hat{f}\mathrm{ea}{t}_{\min }}---\left(6\right)$

对空中交通管理系统中采集到的真实雷达航迹观测数据，经过上述坐标转换、特征计算、数据滤波和归一化等处理，可以建立基于关键飞行特征的飞行航迹样本库，包括了上述3种航迹特征，用以训练意图识别模型中关键参数的基础数据。

步骤二，根据飞行计划信息、航路点位置数据，建立典型飞行意图模型，标注飞行航迹样本的意图类别。

民航航班执行飞行任务之前，空中交通管制单位会预先制定飞行计划（Flight Plan），包括了飞行过程中经过的一系列航路点位置坐标、速度限制、飞行高度限制等信息，因此，采用飞行计划信息可以较为准确地建立飞行意图模型。本步骤根据飞行计划信息、航路点位置坐标等数据，构建典型的飞行意图模型。然后，在步骤一基础上，将飞行航迹样本与飞行意图模型进行关联，标注出飞行航迹样本所属的意图类别。

飞行航迹通常包括一系列的直线段和弧线段，线段之间由航迹改变点（TrajectoryChange Point，TCP）相互连接。因此，本实施例采用飞行计划中的一系列TCP点建立飞行意图模型，一般考虑飞行计划上距离当前位置最近的航迹改变点TCP和之后连续3个航迹改变点，记为TCP+1,TCP+2,TCP+3，这样可以有效地降低计算量。同时，按照空间位置，将意图模型分为水平意图（I_Hi,i＝1,2,3,4,5）和垂直意图（I_Vi,i＝1,2,3,4,5）两类，如表1所示。

表1意图模型

序号	意图模型描述
		IH1	在水平飞行剖面上，保持当前航向，沿直线飞行
IH2	在水平飞行剖面上，飞往TCP
		IH3	在水平飞行剖面上，飞往TCP+1
IH4	在水平飞行剖面上，飞往TCP+2
		IH5	在水平飞行剖面上，飞往TCP+3
IV1	在垂直飞行剖面上，保持TCP高度
		IV2	在垂直飞行剖面上，爬升/下降到TCP高度
IV3	在垂直飞行剖面上，爬升/下降到TCP+1高度
		IV4	在垂直飞行剖面上，爬升/下降到TCP+2高度
IV5	在垂直飞行剖面上，爬升/下降到TCP+3高度

图2和3表示水平和垂直意图识别场景，其中，图2表示在水平飞行剖面X-Y上，航班有两个可能到达的意图，分别为“飞往TCP点”，记作I_H2和“飞往TCP+1点”记作I_H3。图3表示在垂直飞行剖面Z上，航班有两个可能到达的意图，分别为“爬升到TCP高度”记作I_V2和“下降到TCP+1高度”记作I_V3。根据飞行场景的实际情况，可确定上述意图模型的具体参数取值。然后根据实际的雷达飞行航迹样本，判断样本所属的意图类别，对样本集进行类别标注，扩展为带有意图类别的样本集。

步骤三，根据隐马尔科夫模型原理，建立飞行意图识别模型，根据期望最大学习算法，训练识别模型的参数。

为了准确识别飞行航迹当前飞行状态所属的意图类别，本实施例采用隐马尔科夫模型（HMM）建立上述10种飞行意图对应的识别模型。然后，将每种意图模型对应航迹样本集作为训练数据，根据期望最大学习算法，训练得到识别模型参数的最优取值。

HMM是一个状态不可直接观测的马尔科夫模型，记为λ(π,A,e)，λ表示HMM模型参数取值的集合，π表示初始时刻k=1状态分布矩阵，A表示k时刻到k+1时刻的状态转移矩阵，e表示状态发出的观测量的分布矩阵。对于包括M个离散状态s的HMM模型，在k＝1时刻，第i个状态记为s_i(1)，其概率分布记为P(s_i(1))＝π_i，因此，λ(π,A,e)模型的初始状态分布矩阵π＝{π_i，1≤i≤M}。在λ(π,A,e)模型中，当k时刻的i状态s_i(k)到k+1时刻转移为j状态s_j(k+1)，状态转移过程表示为条件概率A_ij＝P(s_j(k+1)|s_i(k))，因此，状态转移矩阵A＝{Ai_j,1≤i≤M,1≤j≤M}。在λ(π,A,e)模型中，i状态s_i发出的观测量集合记为Y，Y的概率分布由e_i＝P(Y|s_i)刻画，通常设为高斯分布，均值为μ_i，方差为因此，e＝{e_i,1≤i≤M}。

对应步骤二中的10个典型飞行意图，选择对应意图类别的飞行航迹样本库进行模型训练。每个样本包含3个关键特征序列（即航向角、飞行速度、爬升率三个序列）和1个确定的意图类别，其中，每个关键特征序列可以表示为关键特征feat(k)从时刻k=1到时刻k＝N的时间序列，记为F＝{feat(k)},k＝1,...,N。在训练HMM的过程中，关键特征序列即为状态对应的观测量集合，即Y=F。因此，为了识别10种典型飞行意图，需建立30个不同参数的HMM模型，记为

I＝1,2,...,10,j＝1,2,3.。然后，采用对应意图类别的样本集作为观测量，应用期望最大学习算法，分别对上述30个HMM模型进行训练，最终得到上述30个模型对应的最优参数。期望最大学习算法具体为：

期望最大学习算法的基本思想是给定观测量集合F＝{feat(k)},k＝1,...,N时，计算HMM模型参数的最大似然估计，如公式（7）：

${\mathrm{\λ}}^{*}(\mathrm{\π},A,e)=\underset{\mathrm{\λ}}{\arg \max }P\left(F\right|\mathrm{\λ}(\mathrm{\π},A,e))---\left(7\right)$

因此，λ^*(π,A,e)为给定观测量集合F的条件下，HMM模型λ(π,A,e)的最优参数。

为了计算方便，采用前向-后向算法，首先定义前向算法：

${\mathrm{\α}}_i\left(k\right)\stackrel{\mathrm{def}}{=}P(\mathrm{feat}\left(1\right),...,\mathrm{feat}\left(k\right),s_i\left(k\right)|\mathrm{\λ})---\left(8\right)$

α_i(k)表示HMM模型参数取值为λ条件下，第k时刻状态为i，且时刻1到时刻k的观测量集合为feat(1),...,feat(k)的条件概率。使用初始状态分布π对α_i(k)进行初始化，如公式（9）：

α_i(1)＝π_ie_i(feat(1)),i＝1,...,M.    （9）

为有效地计算α_i(k)，采用迭代公式求解，如公式（10）：

${\mathrm{\α}}_j(k+1)=\left[\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}\right]e_j\left(\mathrm{feat}(k+1)\right),i=1,...,M.---\left(10\right)$

其中，假设第k+1时刻的状态为j，观测量为feat(k+1)。采用第k时刻的α_i(k)，状态转移概率A_ij，以及状态j的发出观测量的分布e_j，计算出第k+1时刻的α_j(k+1)。

然后，定义后向算法，在时刻k（k≤N）剩余观测序列的概率记为：

${\mathrm{\β}}_i\left(k\right)\stackrel{\mathrm{def}}{=}P(\mathrm{feat}(k+1),...,\mathrm{feat}\left(N\right),s_i\left(k\right)|\mathrm{\λ})---\left(11\right)$

β_i(k)表示HMM模型参数取值为λ条件下，第k时刻状态为i，且时刻k+1到时刻N的观测量集合为feat(k+1),...,feat(N)的条件概率。

后向算法的迭代计算规则如公式（12）、（13）所示：

β_i(N)＝1    （12）

${\mathrm{\β}}_j\left(k\right)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_i(k+1)---\left(13\right)$

其中，假设第k时刻的状态为j，观测量为feat(k+1)。采用第k+1时刻的β_i(k+1)，状态转移概率A_ij，以及状态j的发出的观测量的分布e_j，计算出第k时刻的β_j(k)。而且，第N时刻的β_i(N)为1。

使用α_i(k)和β_i(k)，可以进一步计算出γ_i(k)和ξ_ij(k)，分别如（14）和（15）所示：

${\mathrm{\γ}}_i\left(k\right)=P\left(s_i\left(k\right)\right|\mathrm{feat}\left(k\right),\mathrm{\λ})=\frac{{\mathrm{\α}}_i\left(k\right){\mathrm{\β}}_i\left(k\right)}{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\α}}_i\left(k\right){\mathrm{\β}}_i\left(k\right)}---\left(14\right)$

${\mathrm{\ξ}}_{\mathrm{ij}}\left(k\right)=P(s_i\left(k\right),s_j(k+1)|F,\mathrm{\λ})$

$=\frac{{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_j(k+1)}{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\Σ}}_{j=1}^M{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_j(k+1)}---\left(15\right)$

其中，γ_i(k)表示HMM模型参数取值为λ，且给定时刻k的观测值feat(k)条件下，状态为i的条件概率。由式（14）可知，γ_i(k)可通过前向-后向算法定义的α_i(k)和β_i(k)得到。

ξ_ij(k)表示HMM模型参数取值为λ，且给定观测量集合F条件下，时刻k状态s_i(k)和时刻k+1状态s_j(k+1)的条件概率。由式（15）可知，ξ_ij(k)可通过前向-后向算法定义的α_i(k)和β_i(k)得到。

结合上述公式，可以得到HMM更新后的模型参数λ^new(π^new,A^new,e^new)，其中，初始分布矩阵π^new、状态转移矩阵A^new和状态发出的观测量分布矩阵e^new的具体参数计算如下：

${\mathrm{\π}}_i^{\mathrm{new}}={\mathrm{\γ}}_i\left(1\right)---\left(16\right)$

$A_{\mathrm{ij}}^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^{N-1}{\mathrm{\ξ}}_{\mathrm{ij}}\left(k\right)}{{\mathrm{\Σ}}_{k=1}^{N-1}{\mathrm{\γ}}_{\mathrm{ij}}\left(k\right)}---\left(17\right)$

${\mathrm{\μ}}_i^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)\mathrm{feat}\left(k\right)}{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)}---\left(18\right)$

${\mathrm{\σ}}_i^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right){(\mathrm{feat}\left(k\right)-{\mathrm{\μ}}_i)}^2}{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)}---\left(19\right)$

在训练过程中，随着样本的不断增加，重复上述步骤，用λ^new(π^new,A^new,e^new)替换λ(π,A,e)，将收敛到模型参数的局部最大值，随着迭代次数的增加，最终可得到HMM模型的最优参数λ^*。

通过上述期望最大学习算法，分别对30个HMM模型进行训练，得到30个模型对应的最优参数，记为I=1,2,...,10,j＝1,2,3.。将上述最优参数分别带入30个典型意图识别的HMM模型中，得到30个最优的意图识别模型，用于步骤四的飞行意图识别。

步骤四，采用飞行意图识别模型，根据前向算法，计算飞行航迹样本当前时刻的局部飞行意图。采用滚动时间窗，进行加权求和，得到最终飞行意图。

在识别未知飞行意图的过程中，首先采用步骤三训练得到30个最优HMM模型

(π,A,e),I＝1,2,...,10,j＝1,2,3.，根据步骤三中定义的前向算法，估计当前时刻N的观测值与之前N-1个观测值组成的序列值F＝{feat(k)},k＝1,...,N-1,N的似然值，似然值最大的作为当前时刻的局部飞行意图。然后，以当前时刻为基准，采用滚动时间窗，考虑对未来m个观测序列的局部意图似然值，进行加权求和，似然值最大的意图判定为全局飞行意图。

给定观测量序列F＝{feat(k)},k＝1,...,N-1,N.，基于HMM模型，采用前向算法，估计当前时刻N的局部意图似然值，如公式（20）所示

$P_{I_j,N}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(N\right),I=\mathrm{1,2},...,10,j=\mathrm{1,2,3}.---\left(20\right)$

其中，

表示30个HMM模型分别与观测序列F进行匹配得到的似然值。在实际应用中，根据具体场景，从典型意图中选择相关性较高的典型意图模型作为候选意图，排除相关性低的意图模型，因此，每次计算并不需要遍历全部30个HMM模型，这样可以有效降低计算量。

局部意图I_local判别如公式（21）所示

$I_{\mathrm{local}}={\widetilde{\mathrm{\λ}}}_{I_j}^{*}=\underset{{\mathrm{\λ}}_{I_j}^{*}}{\arg \max }{P}_{I_j,N}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})---\left(21\right)$

通过比较似然值的大小，将似然值最大的HMM模型

对应的意图作为局部意图I_local。I_local可以较好地表示当前时刻的飞行意图，但是为了更准确的识别航班一段时间的飞行趋势，还需要识别全局飞行意图，具体方法如下。

以当前时刻k＝N为起始位置，考虑k＝N,N+1,...,N+m时刻共m个观测序列组成滚动时间窗[N,N+m]，并设定局部似然值的权重w(k)＝e^k-(N+m)，即随时间增加，权重不断增大，使得接近意图点的航迹特征具有最大的权重。

全局飞行意图I_globe计算如公式（22）所示：

$I_{\mathrm{globe}}={\widetilde{\mathrm{\λ}}}_{I_j}^{*}=\underset{{\mathrm{\λ}}_{I_j}^{*}}{\arg \max }\underset{k=n}{\overset{N+m}{\mathrm{\Σ}}}w\left(k\right)\·P_{I_j,k}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})---\left(22\right)$

通过对滚动时间窗[N,N+m]内的意图模型与观测序列F匹配的似然值进行加权求和，然后比较似然值和的大小，将似然值和最大的HMM模型对应的意图作为全局意图I_globe，可以较好的刻画一段时间内的飞行趋势。

Claims (4)

1.一种基于雷达飞行航迹观测数据的飞行意图识别方法，具体包括如下步骤：

步骤一、根据雷达航迹位置观测数据，计算航向角、飞行速度、爬升率，建立飞行航迹样本库；

根据雷达航迹三维位置观测数据，提取飞行航迹的航向角、飞行速度、爬升率，并对采集到的航迹样本进行预处理，建立飞行航迹样本库，作为飞行意图识别方法的基础数据；

步骤二、根据飞行计划信息、航路点位置数据，建立典型飞行意图模型，标注飞行航迹样本的意图类别；

根据飞行计划信息、航路点位置坐标，构建典型的飞行意图模型；在步骤一的基础上，将飞行航迹样本与飞行意图模型进行关联，标注出飞行航迹样本所属的飞行意图模型类别；

具体为：采用飞行计划中的航迹改变点，即TCP点，描述飞行意图模型，考虑飞行计划上距离当前位置最近的航迹改变点TCP和之后连续3个航迹改变点，记为TCP+1,TCP+2,TCP+3，根据空间位置，将意图模型分为水平意图I_Hi,i＝1,2,3,4,5和垂直意图I_Vi,i＝1,2,3,4,5两类，如表1所示：

表1意图模型

序号 意图模型描述 I_H1 在水平飞行剖面上，保持当前航向，沿直线飞行 I_H2 在水平飞行剖面上，飞往TCP I_H3 在水平飞行剖面上，飞往TCP+1 I_H4 在水平飞行剖面上，飞往TCP+2 I_H5 在水平飞行剖面上，飞往TCP+3 I_V1 在垂直飞行剖面上，保持TCP高度 I_V2 在垂直飞行剖面上，爬升/下降到TCP高度 I_V3 在垂直飞行剖面上，爬升/下降到TCP+1高度

I_V4 在垂直飞行剖面上，爬升/下降到TCP+2高度 I_V5 在垂直飞行剖面上，爬升/下降到TCP+3高度

根据飞行场景的实际情况，确定上述意图模型的具体参数取值；然后根据实际的雷达飞行航迹样本，判断样本所属的意图类别，对样本集进行类别标注，扩展为带有意图类别的样本集；

步骤三、根据隐马尔科夫模型原理，建立飞行意图识别模型，根据期望最大学习算法，训练识别模型的参数；

采用隐马尔科夫模型HMM，建立步骤二中10种飞行意图模型对应的30个飞行意图识别模型；然后，将每种意图模型对应飞行航迹样本集作为训练数据，根据期望最大学习算法，训练得到识别模型参数的最优取值；

步骤四、采用飞行意图识别模型，根据前向算法，计算飞行航迹样本当前时刻的局部飞行意图，采用滚动时间窗，进行加权求和，得到最终飞行意图；

在识别未知飞行意图的过程中，首先采用步骤三训练得到30个最优HMM模型

I＝1,2,...,10,j＝1,2,3.，根据前向算法，估计当前时刻N的观测值与之前N-1个观测值组成的序列值F＝{feat(k)},k＝1,...,N-1,N的似然值，似然值最大的作为当前时刻的局部飞行意图；然后，以当前时刻为基准，采用滚动时间窗，考虑对未来m个观测序列的局部意图似然值，进行加权求和，似然值最大的意图判定为全局飞行意图。

2.根据权利要求1所述的一种基于雷达飞行航迹观测数据的飞行意图识别方法，所述的步骤一中，当雷达航迹三维位置观测数据采用WGS-84坐标系下的三维位置观测数据时，具体为：

雷达航迹三维位置观测数据采用WGS-84坐标系下的三维位置观测数据，记为Traj(k)＝{lon(k),lat(k),alt(k)},k＝1,...,N，其中，lon(k)是k时刻的经度，lat(k)是k时刻的纬度，alt(k)是k时刻的高度；

首先将航迹数据从WGS-84坐标系转换到地心地固直角坐标系ECEF，具体坐标转换如公式（1）所示：

$\left\{\begin{array}{c}x\left(k\right)=(N_e+\mathrm{alt}\left(k\right))\cos \left(\mathrm{lat}\left(k\right)\right)\cos \left(\mathrm{lon}\left(k\right)\right)\\ y\left(k\right)=(N_e+\mathrm{alt}\left(k\right))\cos \left(\mathrm{lat}\left(k\right)\right)\sin \left(\mathrm{lon}\left(k\right)\right)\\ z\left(k\right)=[N_e(1-e^2)+\mathrm{alt}\left(k\right)]\sin \left(\mathrm{lat}\left(k\right)\right)\end{array}\right.---\left(1\right)$

在ECEF坐标系下，原点为地球质心，航迹表示为Traj(k)＝{x(k),y(k),z(k)},k＝1,...,N，其中，x(k)为由原点指向经纬度(0,0)位置的X轴上的坐标点，y(k)为由原点指向90°经线的Y轴上的坐标点，z(k)为由原点向北沿地球自转方向的Z轴上的坐标点；式（1）中，lon(k)是WGS-84坐标系下k时刻的经度，lat(k)是WGS-84坐标系下k时刻的纬度，alt(k)是WGS-84坐标系下k时刻的高度；

是主垂直面的曲率半径；a是地球椭球的长半轴，即地球赤道半径，取6378137.0米；

是地球椭球偏心率，且f＝1-b/a是地球椭球扁平率，其中，b是地球椭球的短半轴，即地球极半径，取6356752.3米；

在ECEF坐标系下，根据Traj(k)＝{x(k),y(k),z(k)},k＝1,...,N计算航向角ψ(k)、飞行速度v(k)、爬升率c(k)三个关键特征参数；其中，ψ(k)用来描述飞行航迹在水平平面上转弯机动特征，v(k)用来描述飞行航迹在水平平面上加减速运动特征，c(k)用来描述飞行航迹在垂直平面上高度变化特征，计算公式分别如（2）（3）（4）所示，其中Δt是2倍雷达数据更新率；

$\mathrm{\ψ}\left(k\right)=\arctan \frac{y(k+1)-y(k-1)}{x(k+1)-x(k-1)},k=2,...,N-1.---\left(2\right)$

$v\left(k\right)=\frac{y(k+1)-y(k-1)}{x(k+1)-x(k-1)},k=2,...,N-1.---\left(3\right)$

$c\left(k\right)=\frac{z(k+1)-z(k-1)}{\mathrm{\Δt}},k=2,...,N-1.---\left(4\right)$

然后，采用低通滤波器对航迹特征滤波，航迹特征包括航向角、飞行速度、爬升率，并进行数据归一化处理，低通滤波器如公式（5）所示：

$\tilde{f}\mathrm{eat}\left(1\right)=\mathrm{feat}\left(1\right)$     （5）

$\tilde{f}\mathrm{eat}\left(k\right)=\mathrm{\θfeat}\left(k\right)+(1-\mathrm{\θ})\tilde{f}\mathrm{eat}(k-1),k=2,...,N-1.$

其中，feat(k)是k时刻滤波前航迹特征，

是经过低通滤波器后的航迹特征，θ取为系数；

数据归一化过程通过公式（6），将数据映射到[0,1]区间内：

$\hat{f}\mathrm{eat}\left(k\right)=\frac{\tilde{f}\mathrm{eat}\left(k\right)-\tilde{f}\mathrm{ea}{t}_{\min }}{\hat{f}\mathrm{ea}{t}_{\max }-\hat{f}\mathrm{ea}{t}_{\min }}---\left(6\right)$

对空中交通管理系统中采集到的真实雷达航迹观测数据，经过上述坐标转换、特征计算、数据滤波和归一化等处理，建立基于关键飞行特征的飞行航迹样本库，包括了上述3种航迹特征，用以训练意图识别模型中关键参数的基础数据。

3.根据权利要求1所述的一种基于雷达飞行航迹观测数据的飞行意图识别方法，所述的步骤三中，具体为：

隐马尔科夫模型HMM记为λ(π,A,e)，λ表示HMM模型参数取值的集合，π表示初始时刻k=1状态分布矩阵，A表示k时刻到k+1时刻的状态转移矩阵，e表示状态发出的观测量的分布矩阵；对于包括M个离散状态s的HMM模型，在k＝1时刻，第i个状态记为s_i(1)，其概率分布记为P(s_i(1))＝π_i，因此，λ(π,A,e)模型的初始状态分布矩阵π＝{π_i，1≤i≤M}；在λ(π,A,e)模型中，当k时刻的i状态s_i(k)到k+1时刻转移为j状态s_j(k+1)，状态转移过程表示为条件概率A_ij＝P(s_j(k+1)|s_i(k))，因此，状态转移矩阵A＝{A_ij,1≤i≤M,1≤j≤M}；在λ(π,A,e)模型中，i状态s_i发出的观测量集合记为Y，Y的概率分布由e_i＝P(Y|s_i)刻画，设为高斯分布，均值为μ_i，方差为

e＝{e_i,1≤i≤M}；

对应步骤二中的10个典型飞行意图模型，选择对应意图类别的飞行航迹样本库对飞行意图识别模型进行模型训练；每个样本包含3个关键特征序列和1个确定的意图类别，3个关键特征序列即航向角、飞行速度、爬升率三个序列，其中，每个关键特征序列表示为关键特征feat(k)从时刻k=1到时刻k＝N的时间序列，记为F＝{feat(k)},k＝1,...,N；在训练HMM的过程中，关键特征序列即为状态对应的观测量集合，即Y=F；因此，为了识别10种典型飞行意图，需建立30个不同参数的HMM模型，记为

I＝1,2,...,10,j＝1,2,3.；然后，采用对应意图类别的样本集作为观测量，应用期望最大学习算法，分别对上述30个HMM模型进行训练，最终得到上述30个模型对应的最优参数；

期望最大学习算法具体为：

期望最大学习算法的为给定观测量集合F＝{feat(k)},k＝1,...,N时，计算HMM模型参数的最大似然估计，如公式（7）：

${\mathrm{\λ}}^{*}(\mathrm{\π},A,e)=\underset{\mathrm{\λ}}{\arg \max }P\left(F\right|\mathrm{\λ}(\mathrm{\π},A,e))---\left(7\right)$

因此，λ^*(π,A,e)为给定观测量集合F的条件下，HMM模型λ(π,A,e)的最优参数；

为了计算方便，采用前向-后向算法，首先定义前向算法：

${\mathrm{\α}}_i\left(k\right)\stackrel{\mathrm{def}}{=}P(\mathrm{feat}\left(1\right),...,\mathrm{feat}\left(k\right),s_i\left(k\right)|\mathrm{\λ})---\left(8\right)$

α_i(k)表示HMM模型参数取值为λ条件下，第k时刻状态为i，且时刻1到时刻k的观测量集合为feat(1),...,feat(k)的条件概率；使用初始状态分布π对α_i(k)进行初始化，如公式（9）：

α_i(1)＝π_ie_i(feat(1)),i＝1,...,M.    （9）

为有效地计算α_i(k)，采用迭代公式求解，如公式（10）：

${\mathrm{\α}}_j(k+1)=\left[\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}\right]e_j\left(\mathrm{feat}(k+1)\right),i=1,...,M.---\left(10\right)$

其中，假设第k+1时刻的状态为j，观测量为feat(k+1)；采用第k时刻的α_i(k)，状态转移概率A_ij，以及状态j的发出观测量的分布e_j，计算出第k+1时刻的α_j(k+1)；

然后，定义后向算法，在时刻k剩余观测序列的概率记为，k≤N：

${\mathrm{\β}}_i\left(k\right)\stackrel{\mathrm{def}}{=}P(\mathrm{feat}(k+1),...,\mathrm{feat}\left(N\right),s_i\left(k\right)|\mathrm{\λ})---\left(11\right)$

β_i(k)表示HMM模型参数取值为λ条件下，第k时刻状态为i，且时刻k+1到时刻N的观测量集合为feat(k+1),...,feat(N)的条件概率；

后向算法的迭代计算规则如公式（12）、（13）所示：

β_i(N)＝1    （12）

${\mathrm{\β}}_j\left(k\right)=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{A}_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_i(k+1)---\left(13\right)$

其中，假设第k时刻的状态为j，观测量为feat(k+1)；采用第k+1时刻的β_i(k+1)，状态转移概率A_ij，以及状态j的发出的观测量的分布e_j，计算出第k时刻的β_j(k)；而且，第N时刻的β_i(N)为1；

使用α_i(k)和β_i(k)，进一步计算出γ_i(k)和ξ_ij(k)，分别如（14）和（15）所示：

${\mathrm{\γ}}_i\left(k\right)=P\left(s_i\left(k\right)\right|\mathrm{feat}\left(k\right),\mathrm{\λ})=\frac{{\mathrm{\α}}_i\left(k\right){\mathrm{\β}}_i\left(k\right)}{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\α}}_i\left(k\right){\mathrm{\β}}_i\left(k\right)}---\left(14\right)$

${\mathrm{\ξ}}_{\mathrm{ij}}\left(k\right)=P(s_i\left(k\right),s_j(k+1)|F,\mathrm{\λ})$

$=\frac{{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_j(k+1)}{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\Σ}}_{j=1}^M{\mathrm{\α}}_i\left(k\right)A_{\mathrm{ij}}{e}_j\left(\mathrm{feat}(k+1)\right){\mathrm{\β}}_j(k+1)}---\left(15\right)$

其中，γ_i(k)表示HMM模型参数取值为λ，且给定时刻k的观测值feat(k)条件下，状态为i的条件概率；由式（14）可知，γ_i(k)通过前向-后向算法定义的α_i(k)和β_i(k)得到；

ξ_ij(k)表示HMM模型参数取值为λ，且给定观测量集合F条件下，时刻k状态s_i(k)和时刻k+1状态s_j(k+1)的条件概率；由式（15）可知，ξ_ij(k)通过前向-后向算法定义的α_i(k)和β_i(k)得到；

结合上述公式，得到HMM更新后的模型参数λ^new(π^new,A^new,e^new)，其中，初始分布矩阵π^new、状态转移矩阵A^new和状态发出的观测量分布矩阵e^new的具体参数计算如下：

${\mathrm{\π}}_i^{\mathrm{new}}={\mathrm{\γ}}_i\left(1\right)---\left(16\right)$

$A_{\mathrm{ij}}^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^{N-1}{\mathrm{\ξ}}_{\mathrm{ij}}\left(k\right)}{{\mathrm{\Σ}}_{k=1}^{N-1}{\mathrm{\γ}}_{\mathrm{ij}}\left(k\right)}---\left(17\right)$

${\mathrm{\μ}}_i^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)\mathrm{feat}\left(k\right)}{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)}---\left(18\right)$

${\mathrm{\σ}}_i^{\mathrm{new}}=\frac{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right){(\mathrm{feat}\left(k\right)-{\mathrm{\μ}}_i)}^2}{{\mathrm{\Σ}}_{k=1}^N{\mathrm{\γ}}_i\left(k\right)}---\left(19\right)$

在训练过程中，随着样本的不断增加，重复上述步骤，用λ^new(π^new,A^new,e^new)替换λ(π,A,e)，将收敛到模型参数的局部最大值，随着迭代次数的增加，最终可得到HMM模型的最优参数λ^*；

通过上述期望最大学习算法，分别对30个HMM模型进行训练，得到30个模型对应的最优参数，记为

I＝1,2,...,10,j＝1,2,3.；将上述最优参数分别带入30个典型意图识别的HMM模型中，得到30个最优的意图识别模型，用于步骤四的飞行意图识别。

4.根据权利要求1所述的一种基于雷达飞行航迹观测数据的飞行意图识别方法，所述的步骤四具体为：

给定观测量序列F＝{feat(k)},k＝1,...,N-1,N.，基于HMM模型，采用前向算法，估计当前时刻N的局部意图似然值，如公式（20）所示

$P_{I_j,N}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_i\left(N\right),I=\mathrm{1,2},...,10,j=\mathrm{1,2,3}.---\left(20\right)$

其中，

表示30个HMM模型分别与观测序列F进行匹配得到的似然值；

局部意图I_local判别如公式（21）所示

$I_{\mathrm{local}}={\widetilde{\mathrm{\λ}}}_{I_j}^{*}=\underset{{\mathrm{\λ}}_{I_j}^{*}}{\arg \max }{P}_{I_j,N}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})---\left(21\right)$

通过比较似然值的大小，将似然值最大的HMM模型

对应的意图作为局部意图I_local；

识别全局飞行意图具体方法如下；

以当前时刻k＝N为起始位置，考虑k＝N,N+1,...,N+m时刻共m个观测序列组成滚动时间窗[N,N+m]，并设定局部似然值的权重w(k)＝e^k-(N+m)，即随时间增加，权重不断增大，使得接近意图点的航迹特征具有最大的权重；

全局飞行意图I_globe计算如公式（22）所示：

$I_{\mathrm{globe}}={\widetilde{\mathrm{\λ}}}_{I_j}^{*}=\underset{{\mathrm{\λ}}_{I_j}^{*}}{\arg \max }\underset{k=n}{\overset{N+m}{\mathrm{\Σ}}}w\left(k\right)\·P_{I_j,k}\left(F\right|{\mathrm{\λ}}_{I_j}^{*})---\left(22\right)$

通过对滚动时间窗[N,N+m]内的意图模型与观测序列F匹配的似然值进行加权求和，然后比较似然值和

的大小，将似然值和最大的HMM模型

对应的意图作为全局意图I_globe。

Images