CN103336863B - 基于雷达飞行航迹观测数据的飞行意图识别方法 - Google Patents
基于雷达飞行航迹观测数据的飞行意图识别方法 Download PDFInfo
- Publication number
- CN103336863B CN103336863B CN201310251867.XA CN201310251867A CN103336863B CN 103336863 B CN103336863 B CN 103336863B CN 201310251867 A CN201310251867 A CN 201310251867A CN 103336863 B CN103336863 B CN 103336863B
- Authority
- CN
- China
- Prior art keywords
- flight
- intention
- model
- time
- lambda
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 41
- 239000011159 matrix material Substances 0.000 claims description 19
- 230000007704 transition Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 14
- 230000009194 climbing Effects 0.000 claims description 11
- 238000007476 Maximum Likelihood Methods 0.000 claims description 8
- 230000033001 locomotion Effects 0.000 claims description 8
- 238000005096 rolling process Methods 0.000 claims description 7
- 238000010606 normalization Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 2
- 230000001133 acceleration Effects 0.000 claims description 2
- 238000007781 pre-processing Methods 0.000 claims description 2
- 241000287196 Asthenes Species 0.000 claims 1
- 230000019771 cognition Effects 0.000 abstract 3
- 238000001514 detection method Methods 0.000 abstract 1
- 238000005516 engineering process Methods 0.000 description 5
- 238000012544 monitoring process Methods 0.000 description 3
- 230000005540 biological transmission Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Landscapes
- Radar Systems Or Details Thereof (AREA)
- Traffic Control Systems (AREA)
Abstract
本发明公开了一种基于雷达飞行航迹观测数据的飞行意图识别方法,其中包括:根据雷达航迹位置观测数据,计算航向角、飞行速度、爬升率,建立飞行航迹样本库,作为训练意图识别模型的基础数据;根据飞行计划信息、航路点位置数据,建立典型飞行意图模型,并标注飞行航迹样本的意图类别;根据隐马尔科夫模型(Hidden?Markov?Model,HMM)原理,建立几种典型飞行意图识别模型,采用期望最大学习算法,训练识别模型的参数;采用飞行意图识别模型,根据前向算法,计算飞行航迹样本当前时刻的飞行意图,采用滚动时间窗,对一段时间内的局部飞行意图进行加权求和,得到全局飞行意图。
Description
技术领域
本发明涉及空中交通管理技术,尤其涉及一种基于雷达飞行航迹观测数据的飞行意图识别方法。
背景技术
空中交通管理技术的主要目标是保障空中交通能够安全、高效、有序地在空域中运行。在空中交通管理领域,空中交通态势感知能力(AirTrafficSituationalAwareness,ATSA)是保障飞行安全、提高空中交通运行效率的基础技术。根据空中交通态势的时间属性,ATSA可以分为两个主要阶段:(1)当前时刻飞行器的飞行航迹运动状态感知(2)飞行器未来一段时间的飞行意图识别。其中,当前时刻的飞行航迹运动状态可以直接从空中交通监视系统中获取,例如,雷达监视系统、自动相关监视系统、多点定位系统等。飞行航迹运动状态包括经度、纬度、高度等三维位置信息,以及航向角、速度、爬升率等飞行航迹特征。在获取当前时刻飞行航迹运动状态基础上,管制员需要预测未来一段时间内飞行器可能到达的位置,以探测飞行器之间可能存在的飞行冲突,作为安全飞行的决策基础,确保飞行器保持安全飞行间隔。因此,识别未来一段时间飞行器可能到达的位置,即飞行意图识别,是ATSA的关键技术之一。
现有的飞行意图识别方法主要分为(1)基于飞行员离散操作指令的飞行意图识别(2)基于雷达飞行航迹观测数据的飞行意图识别两类。其中,基于飞行员离散操作指令的飞行意图识别方法是指通过收集飞行员的实时操作指令,采用计划识别(PlanRecognition)技术,推测飞行意图。该方法虽然能够较为准确识别出飞行意图,但是由于从机舱内部采集飞行员实时操作指令数据非常困难,加之空地数据传输的困难,该方法在工程应用中受到很大限制。为了避免数据采集和传输等问题,有学者提出了基于雷达飞行航迹观测数据的飞行意图识别方法,即根据雷达系统获取的飞行器飞行航迹的观测数据,与飞行意图模型进行最优匹配计算,将匹配概率最大的飞行意图,作为飞行意图的识别结果。近年来,随着空中交通流量的快速增长,在主干航路和繁忙飞行终端区空域内空中交通态势呈现高密度、高动态等新特点,为意图识别带来了新挑战。现有的识别方法在进行意图模型匹配时,仅仅考虑单个飞行航迹运动状态,直接进行向量点乘运算,在识别精度方面难以满足高动态空域环境需求,特别是针对存在频繁机动操作的空域环境。
发明内容
本发明的目的是为了解决上述问题,提供一种基于雷达飞行航迹观测数据的飞行意图识别方法,综合考虑多种飞行航迹特征,有效提高飞行意图识别精度。
本发明提供了一种基于雷达飞行航迹观测数据的飞行意图识别方法,包括以下几个步骤:
步骤一:根据雷达航迹位置观测数据,计算航向角、飞行速度、爬升率,建立飞行航迹样本库,作为训练意图识别模型参数的基础数据。
步骤二:根据飞行计划信息、航路点位置数据,建立典型飞行意图模型,标注飞行航迹样本的意图类别;
步骤三:根据隐马尔科夫模型原理,建立飞行意图识别模型,根据期望最大学习算法,训练识别模型的参数;
步骤四:采用飞行意图识别模型,根据前向算法,计算飞行航迹样本当前时刻的飞行意图。采用滚动时间窗,对一段时间的飞行意图进行加权,得到最终飞行意图;
本发明的优点在于:
本发明提供的一种基于雷达飞行航迹观测数据的飞行意图识别方法,通过对雷达航迹历史观测数据进行分析,计算关键航迹特征,建立飞行航迹样本库,并根据飞行计划信息、航路点位置数据等辅助信息,提炼出空中交通管理中的几种典型飞行意图模型,并对飞行航迹样本库中的样本进行类别标注,作为训练飞行意图识别方法具体参数的基础数据。根据隐马尔科夫模型原理,对几种典型飞行意图分别建立识别模型,根据期望最大化学习方法,采用飞行航迹样本集训练出识别模型的参数。根据前向算法,判断未知意图类别的飞行航迹样本当前时刻的飞行意图,并采用滚动时间窗,对时间窗内的飞行意图识别结果进行加权,得到最终的全局飞行意图。
附图说明
图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)所示:
在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秒。
然后,本实施例中采用低通滤波器对航迹特征滤波,并进行数据归一化处理,以避免真实雷达航迹观测数据中的数据错误,减小雷达数据噪声。低通滤波器如公式(5)所示:
其中,feat(k)是k时刻滤波前航迹特征(航向角、飞行速度、爬升率),是经过低通滤波器后的航迹特征,θ取0.4。
数据归一化过程通过公式(6),将数据映射到[0,1]区间内:
对空中交通管理系统中采集到的真实雷达航迹观测数据,经过上述坐标转换、特征计算、数据滤波和归一化等处理,可以建立基于关键飞行特征的飞行航迹样本库,包括了上述3种航迹特征,用以训练意图识别模型中关键参数的基础数据。
步骤二,根据飞行计划信息、航路点位置数据,建立典型飞行意图模型,标注飞行航迹样本的意图类别。
民航航班执行飞行任务之前,空中交通管制单位会预先制定飞行计划(FlightPlan),包括了飞行过程中经过的一系列航路点位置坐标、速度限制、飞行高度限制等信息,因此,采用飞行计划信息可以较为准确地建立飞行意图模型。本步骤根据飞行计划信息、航路点位置坐标等数据,构建典型的飞行意图模型。然后,在步骤一基础上,将飞行航迹样本与飞行意图模型进行关联,标注出飞行航迹样本所属的意图类别。
飞行航迹通常包括一系列的直线段和弧线段,线段之间由航迹改变点(TrajectoryChangePoint,TCP)相互连接。因此,本实施例采用飞行计划中的一系列TCP点建立飞行意图模型,一般考虑飞行计划上距离当前位置最近的航迹改变点TCP和之后连续3个航迹改变点,记为TCP+1,TCP+2,TCP+3,这样可以有效地降低计算量。同时,按照空间位置,将意图模型分为水平意图(IHi,i=1,2,3,4,5)和垂直意图(IVi,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点”,记作IH2和“飞往TCP+1点”记作IH3。图3表示在垂直飞行剖面Z上,航班有两个可能到达的意图,分别为“爬升到TCP高度”记作IV2和“下降到TCP+1高度”记作IV3。根据飞行场景的实际情况,可确定上述意图模型的具体参数取值。然后根据实际的雷达飞行航迹样本,判断样本所属的意图类别,对样本集进行类别标注,扩展为带有意图类别的样本集。
步骤三,根据隐马尔科夫模型原理,建立飞行意图识别模型,根据期望最大学习算法,训练识别模型的参数。
为了准确识别飞行航迹当前飞行状态所属的意图类别,本实施例采用隐马尔科夫模型(HMM)建立上述10种飞行意图对应的识别模型。然后,将每种意图模型对应航迹样本集作为训练数据,根据期望最大学习算法,训练得到识别模型参数的最优取值。
HMM是一个状态不可直接观测的马尔科夫模型,记为λ(π,A,e),λ表示HMM模型参数取值的集合,π表示初始时刻k=1状态分布矩阵,A表示k时刻到k+1时刻的状态转移矩阵,e表示状态发出的观测量的分布矩阵。对于包括M个离散状态s的HMM模型,在k=1时刻,第i个状态记为si(1),其概率分布记为P(si(1))=πi,因此,λ(π,A,e)模型的初始状态分布矩阵π={πi,1≤i≤M}。在λ(π,A,e)模型中,当k时刻的i状态si(k)到k+1时刻转移为j状态sj(k+1),状态转移过程表示为条件概率Aij=P(sj(k+1)|si(k)),因此,状态转移矩阵A={Aij,1≤i≤M,1≤j≤M}。在λ(π,A,e)模型中,i状态si发出的观测量集合记为Y,Y的概率分布由ei=P(Y|si)刻画,通常设为高斯分布,均值为μi,方差为因此,e={ei,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):
因此,λ*(π,A,e)为给定观测量集合F的条件下,HMM模型λ(π,A,e)的最优参数。
为了计算方便,采用前向-后向算法,首先定义前向算法:
αi(k)表示HMM模型参数取值为λ条件下,第k时刻状态为i,且时刻1到时刻k的观测量集合为feat(1),...,feat(k)的条件概率。使用初始状态分布π对αi(k)进行初始化,如公式(9):
αi(1)=πiei(feat(1)),i=1,...,M.(9)
为有效地计算αi(k),采用迭代公式求解,如公式(10):
其中,假设第k+1时刻的状态为j,观测量为feat(k+1)。采用第k时刻的αi(k),状态转移概率Aij,以及状态j的发出观测量的分布ej,计算出第k+1时刻的αj(k+1)。
然后,定义后向算法,在时刻k(k≤N)剩余观测序列的概率记为:
βi(k)表示HMM模型参数取值为λ条件下,第k时刻状态为i,且时刻k+1到时刻N的观测量集合为feat(k+1),...,feat(N)的条件概率。
后向算法的迭代计算规则如公式(12)、(13)所示:
βi(N)=1(12)
其中,假设第k时刻的状态为j,观测量为feat(k+1)。采用第k+1时刻的βi(k+1),状态转移概率Aij,以及状态j的发出的观测量的分布ej,计算出第k时刻的βj(k)。而且,第N时刻的βi(N)为1。
使用αi(k)和βi(k),可以进一步计算出γi(k)和ξij(k),分别如(14)和(15)所示:
其中,γi(k)表示HMM模型参数取值为λ,且给定时刻k的观测值feat(k)条件下,状态为i的条件概率。由式(14)可知,γi(k)可通过前向-后向算法定义的αi(k)和βi(k)得到。
ξij(k)表示HMM模型参数取值为λ,且给定观测量集合F条件下,时刻k状态si(k)和时刻k+1状态sj(k+1)的条件概率。由式(15)可知,ξij(k)可通过前向-后向算法定义的αi(k)和βi(k)得到。
结合上述公式,可以得到HMM更新后的模型参数λnew(πnew,Anew,enew),其中,初始分布矩阵πnew、状态转移矩阵Anew和状态发出的观测量分布矩阵enew的具体参数计算如下:
在训练过程中,随着样本的不断增加,重复上述步骤,用λnew(πnew,Anew,enew)替换λ(π,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)所示
其中,表示30个HMM模型分别与观测序列F进行匹配得到的似然值。在实际应用中,根据具体场景,从典型意图中选择相关性较高的典型意图模型作为候选意图,排除相关性低的意图模型,因此,每次计算并不需要遍历全部30个HMM模型,这样可以有效降低计算量。
局部意图Ilocal判别如公式(21)所示
通过比较似然值的大小,将似然值最大的HMM模型对应的意图作为局部意图Ilocal。Ilocal可以较好地表示当前时刻的飞行意图,但是为了更准确的识别航班一段时间的飞行趋势,还需要识别全局飞行意图,具体方法如下。
以当前时刻k=N为起始位置,考虑k=N,N+1,...,N+m时刻共m个观测序列组成滚动时间窗[N,N+m],并设定局部似然值的权重w(k)=ek-(N+m),即随时间增加,权重不断增大,使得接近意图点的航迹特征具有最大的权重。
全局飞行意图Iglobe计算如公式(22)所示:
通过对滚动时间窗[N,N+m]内的意图模型与观测序列F匹配的似然值进行加权求和,然后比较似然值和的大小,将似然值和最大的HMM模型对应的意图作为全局意图Iglobe,可以较好的刻画一段时间内的飞行趋势。
Claims (4)
1.一种基于雷达飞行航迹观测数据的飞行意图识别方法,具体包括如下步骤:
步骤一、根据雷达航迹三维位置观测数据,提取飞行航迹的航向角、飞行速度、爬升率,并对采集到的航迹样本进行预处理,建立飞行航迹样本库,作为飞行意图识别方法的基础数据;
步骤二、根据飞行计划信息、航路点位置坐标,构建典型的飞行意图模型;在步骤一的基础上,将飞行航迹样本与飞行意图模型进行关联,标注出飞行航迹样本所属的飞行意图模型类别;
具体为:采用飞行计划中的航迹改变点,即TCP点,描述飞行意图模型,考虑飞行计划上距离当前位置最近的航迹改变点TCP和之后连续3个航迹改变点,记为TCP+1,TCP+2,TCP+3,根据空间位置,将飞行意图模型分为水平意图IHi,i=1,2,3,4,5和垂直意图IVi,i=1,2,3,4,5两类,如下所示:
根据飞行场景的实际情况,确定上述飞行意图模型的具体参数取值;然后根据实际的雷达飞行航迹样本,判断样本所属的意图类别,对样本集进行类别标注,扩展为带有意图类别的样本集;
步骤三、采用隐马尔科夫模型HMM,建立步骤二中10种飞行意图模型对应的30个飞行意图识别模型;然后,将每种飞行意图模型对应飞行航迹样本集作为训练数据,根据期望最大学习算法,训练得到飞行意图识别模型参数的最优取值;
步骤四、在识别未知飞行意图的过程中,首先采用步骤三训练得到30个最优HMM模型π表示初始时刻k=1状态分布矩阵,A表示k时刻到k+1时刻的状态转移矩阵,e表示状态发出的观测量的分布矩阵;根据前向算法,估计当前时刻N的观测值与之前N-1个观测值组成的序列值F={feat(k)}的似然值,其中k=1,...,N-1,N,feat(k)表示HMM模型给定时刻k的观测值,似然值最大的作为当前时刻的局部飞行意图;然后,以当前时刻为基准,采用滚动时间窗,考虑对未来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)所示:
在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)、飞行速度v(k)、爬升率c(k)三个关键特征参数,其中k=1,...,N;其中,ψ(k)用来描述飞行航迹在水平平面上转弯机动特征,v(k)用来描述飞行航迹在水平平面上加减速运动特征,c(k)用来描述飞行航迹在垂直平面上高度变化特征,计算公式分别如(2)(3)(4)所示,其中Δt是2倍雷达数据更新率;
然后,采用低通滤波器对航迹特征滤波,航迹特征包括航向角、飞行速度、爬升率,并进行数据归一化处理,低通滤波器如公式(5)所示:
其中,feat(k)是k时刻滤波前航迹特征,是经过低通滤波器后的航迹特征,θ取为系数,θ取0.4;
数据归一化过程通过公式(6),将数据映射到[0,1]区间内:
是经过低通滤波器后的航迹特征最大值,是经过低通滤波器后的航迹特征最小值;
对空中交通管理系统中采集到的真实雷达航迹观测数据,经过上述坐标转换、特征计算、数据滤波和归一化处理,建立基于关键飞行特征的飞行航迹样本库,包括了上述航迹特征,用以训练飞行意图识别模型中关键参数的基础数据。
3.根据权利要求1所述的一种基于雷达飞行航迹观测数据的飞行意图识别方法,所述的步骤三中,具体为:
隐马尔科夫模型HMM记为λ(π,A,e),λ表示HMM模型参数取值的集合,π表示初始时刻k=1状态分布矩阵,A表示k时刻到k+1时刻的状态转移矩阵,e表示状态发出的观测量的分布矩阵;对于包括M个离散状态s的HMM模型,在k=1时刻,第i个状态记为si(1),其概率分布记为P(si(1))=πi,因此,λ(π,A,e)模型的初始状态分布矩阵π={πi,1≤i≤M};在λ(π,A,e)模型中,当k时刻的i状态si(k)到k+1时刻转移为j状态sj(k+1),状态转移过程表示为条件概率Aij=P(sj(k+1)|si(k)),因此,状态转移矩阵A={Aij,1≤i≤M,1≤j≤M};在λ(π,A,e)模型中,i状态si发出的观测量集合记为Y,Y的概率分布由ei=P(Y|si)刻画,设为高斯分布,均值为μi,方差为e={ei,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模型,记为然后,采用对应意图类别的样本集作为观测量,应用期望最大学习算法,分别对上述30个HMM模型进行训练,最终得到上述30个模型对应的最优参数;
期望最大学习算法具体为:
期望最大学习算法的为给定观测量集合F={feat(k)},k=1,...,N时,计算HMM模型参数的最大似然估计,如公式(7):
因此,λ*(π,A,e)为给定观测量集合F的条件下,HMM模型λ(π,A,e)的最优参数;
为了计算方便,采用前向-后向算法,首先定义前向算法:
αi(k)表示HMM模型参数取值为λ条件下,第k时刻状态为i,且时刻1到时刻k的观测量集合为feat(1),...,feat(k)的条件概率;使用初始状态分布π对αi(k)进行初始化,如公式(9):
αi(1)=πiei(feat(1)),i=1,...,M.(9)
为有效地计算αi(k),采用迭代公式求解,如公式(10):
其中,假设第k+1时刻的状态为j,观测量为feat(k+1);采用第k时刻的αi(k),状态转移概率Aij,以及状态j的发出观测量的分布ej,计算出第k+1时刻的αj(k+1);
然后,定义后向算法,在时刻k剩余观测序列的概率记为,k≤N:
βi(k)表示HMM模型参数取值为λ条件下,第k时刻状态为i,且时刻k+1到时刻N的观测量集合为feat(k+1),...,feat(N)的条件概率;
后向算法的迭代计算规则如公式(12)、(13)所示:
βi(N)=1(12)
其中,假设第k时刻的状态为j,观测量为feat(k+1);采用第k+1时刻的βi(k+1),状态转移概率Aij,以及状态j的发出的观测量的分布ej,计算出第k时刻的βj(k);而且,第N时刻的βi(N)为1;
使用αi(k)和βi(k),进一步计算出γi(k)和ξij(k),分别如(14)和(15)所示:
其中,γi(k)表示HMM模型参数取值为λ,且给定时刻k的观测值feat(k)条件下,状态为i的条件概率;由式(14)可知,γi(k)通过前向-后向算法定义的αi(k)和βi(k)得到;
ξij(k)表示HMM模型参数取值为λ,且给定观测量集合F条件下,时刻k状态si(k)和时刻k+1状态sj(k+1)的条件概率;由式(15)可知,ξij(k)通过前向-后向算法定义的αi(k)和βi(k)得到;
结合上述公式,得到HMM更新后的模型参数λnew(πnew,Anew,enew),其中,初始分布矩阵πnew、状态转移矩阵Anew和状态发出的观测量分布矩阵enew的具体参数计算如下:
在训练过程中,随着样本的不断增加,重复上述步骤,用λnew(πnew,Anew,enew)替换λ(π,A,e),将收敛到模型参数的局部最大值,随着迭代次数的增加,最终可得到HMM模型的最优参数λ*;
通过上述期望最大学习算法,分别对30个HMM模型进行训练,得到30个模型对应的最优参数,记为将上述最优参数分别带入30个典型意图识别的HMM模型中,得到30个最优的飞行意图识别模型,用于步骤四的飞行意图识别。
4.根据权利要求1所述的一种基于雷达飞行航迹观测数据的飞行意图识别方法,所述的步骤四具体为:
给定观测量序列F={feat(k)},k=1,...,N-1,N.,基于HMM模型,采用前向算法,估计当前时刻N的局部意图似然值,如公式(20)所示
其中,表示30个HMM模型分别与观测序列F进行匹配得到的似然值;
局部意图Ilocal判别如公式(21)所示
通过比较似然值的大小,将似然值最大的HMM模型对应的意图作为局部意图Ilocal;
识别全局飞行意图具体方法如下;
以当前时刻k=N为起始位置,考虑k=N,N+1,...,N+m时刻共m+1个观测序列组成滚动时间窗[N,N+m],并设定局部似然值的权重w(k)=ek-(N+m),即随时间增加,权重不断增大,使得接近意图点的航迹特征具有最大的权重;
全局飞行意图Iglobe计算如公式(22)所示:
通过对滚动时间窗[N,N+m]内的飞行意图模型与观测序列F匹配的似然值进行加权求和,然后比较似然值和的大小,将似然值和最大的HMM模型对应的意图作为全局意图Iglobe。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310251867.XA CN103336863B (zh) | 2013-06-24 | 2013-06-24 | 基于雷达飞行航迹观测数据的飞行意图识别方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310251867.XA CN103336863B (zh) | 2013-06-24 | 2013-06-24 | 基于雷达飞行航迹观测数据的飞行意图识别方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103336863A CN103336863A (zh) | 2013-10-02 |
CN103336863B true CN103336863B (zh) | 2016-06-01 |
Family
ID=49245027
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310251867.XA Active CN103336863B (zh) | 2013-06-24 | 2013-06-24 | 基于雷达飞行航迹观测数据的飞行意图识别方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103336863B (zh) |
Families Citing this family (29)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104504934B (zh) * | 2014-12-30 | 2017-03-08 | 江苏理工学院 | 一种航海交通管制方法 |
CN106228850A (zh) * | 2014-12-30 | 2016-12-14 | 江苏理工学院 | 基于滚动规划策略的船舶轨迹实时预测方法 |
CN106205221A (zh) * | 2015-01-07 | 2016-12-07 | 江苏理工学院 | 基于4d航迹运行的航空器轨迹预测方法 |
CN106297419B (zh) * | 2015-01-07 | 2019-01-25 | 江苏理工学院 | 一种基于4d的航空器轨迹预测方法 |
CN106846924A (zh) * | 2015-01-07 | 2017-06-13 | 江苏理工学院 | 用于冲突预警的空中交通管制系统 |
CN106340209A (zh) * | 2015-01-07 | 2017-01-18 | 江苏理工学院 | 基于4d航迹运行的空中交通管制系统的管制方法 |
CN104794345B (zh) * | 2015-04-21 | 2017-12-01 | 四川大学 | 航迹跟踪中扁平型多假设关联处理方法 |
CN104964690A (zh) * | 2015-05-22 | 2015-10-07 | 西北工业大学 | 一种基于期望最大化的空中机动目标航迹粘连方法 |
CN105185163B (zh) * | 2015-06-02 | 2017-07-25 | 北京航空航天大学 | 飞行路径选择方法和装置、飞行器和空中交通管理系统 |
CN106057195A (zh) * | 2016-05-25 | 2016-10-26 | 东华大学 | 一种基于嵌入式音频识别的无人机探测系统 |
CN109031275A (zh) * | 2018-07-02 | 2018-12-18 | 北京理工大学 | 一种基于高分辨全极化雷达的昆虫水平飞行速度提取方法 |
CN110428830B (zh) * | 2019-07-17 | 2021-09-21 | 上海麦图信息科技有限公司 | 一种基于正则表达式的空管指令意图识别方法 |
CN110349294B (zh) * | 2019-07-22 | 2021-12-24 | 北京润科通用技术有限公司 | 飞行数据的飞行模式识别方法及装置 |
CN110502817B (zh) * | 2019-08-13 | 2022-04-08 | 成都飞机工业(集团)有限责任公司 | 一种三维飞行剖面参数化设计方法 |
CN111142126B (zh) * | 2019-12-11 | 2022-03-01 | 四川九洲空管科技有限责任公司 | 一种基于单通道ads-b地面站的防欺骗综合解决方法 |
CN112541608B (zh) * | 2020-02-19 | 2023-10-20 | 深圳中科保泰空天技术有限公司 | 无人机起飞点预测方法及装置 |
CN111477036B (zh) * | 2020-04-08 | 2021-01-29 | 中国电子科技集团公司第二十八研究所 | 一种空管自动化系统航空器高度异常检测方法 |
CN111638646B (zh) * | 2020-05-29 | 2024-05-28 | 平安科技(深圳)有限公司 | 四足机器人行走控制器训练方法、装置、终端及存储介质 |
CN111782755B (zh) * | 2020-07-20 | 2021-05-25 | 中国人民解放军国防科技大学 | 基于虚拟网格字典的目标行进类意图识别方法和装置 |
CN111783231B (zh) * | 2020-07-20 | 2021-04-27 | 中国人民解放军国防科技大学 | 基于单元分布热力网格的目标任务意图识别方法和装置 |
CN112330982B (zh) * | 2020-10-15 | 2024-06-21 | 中国民用航空中南地区空中交通管理局 | 一种应用于终端区的中期冲突预警方法、设备、存储介质 |
CN112232567A (zh) * | 2020-10-16 | 2021-01-15 | 南京智慧航空研究院有限公司 | 不明意图飞行轨迹对航路影响概率预测方法 |
CN112817017A (zh) * | 2021-02-06 | 2021-05-18 | 北京航空航天大学 | 一种考虑涌浪影响的gnss反射信号建模方法 |
CN113095504B (zh) * | 2021-04-13 | 2022-10-21 | 中国人民解放军空军工程大学 | 一种目标轨迹预测系统及预测方法 |
CN113450599B (zh) * | 2021-05-31 | 2022-12-23 | 北京军懋国兴科技股份有限公司 | 一种飞行动作实时识别方法 |
CN114038242B (zh) * | 2021-11-18 | 2023-12-12 | 中国航空无线电电子研究所 | 一种基于多智能体的大规模航空器运动仿真方法及装置 |
CN114078339B (zh) * | 2022-01-07 | 2022-05-31 | 北京航空航天大学杭州创新研究院 | 基于飞行趋势推理的时空概率分布网格飞行冲突检测方法 |
CN114358211B (zh) * | 2022-01-14 | 2022-08-23 | 中科世通亨奇(北京)科技有限公司 | 基于多模态深度学习的航空器行为意图识别方法 |
CN116580620B (zh) * | 2023-04-14 | 2024-06-04 | 中国人民解放军空军特色医学中心 | 模拟飞行姿态显示的地平经纬网方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102819966A (zh) * | 2012-09-10 | 2012-12-12 | 南京航空航天大学 | 一种基于雷达航迹的偏航航班过点时间预测方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2482269B1 (en) * | 2011-01-28 | 2017-03-22 | The Boeing Company | Providing data for predicting aircraft trajectory |
-
2013
- 2013-06-24 CN CN201310251867.XA patent/CN103336863B/zh active Active
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102819966A (zh) * | 2012-09-10 | 2012-12-12 | 南京航空航天大学 | 一种基于雷达航迹的偏航航班过点时间预测方法 |
Non-Patent Citations (1)
Title |
---|
高空空域飞机高度保持性能分析;蔡开泉等;《航空学报》;20090430;第30卷(第4期);726-731 * |
Also Published As
Publication number | Publication date |
---|---|
CN103336863A (zh) | 2013-10-02 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103336863B (zh) | 基于雷达飞行航迹观测数据的飞行意图识别方法 | |
CN114048889B (zh) | 基于长短期记忆网络的飞行器轨迹预测的方法 | |
CN109686125B (zh) | 一种基于hmm的v2x车联网车辆防撞预警系统 | |
Wang et al. | Short-term 4d trajectory prediction using machine learning methods | |
CN109597864B (zh) | 椭球边界卡尔曼滤波的即时定位与地图构建方法及系统 | |
Lin et al. | An algorithm for trajectory prediction of flight plan based on relative motion between positions | |
CN105489069B (zh) | 一种基于svm的低空空域通航飞机冲突检测方法 | |
CN103345587B (zh) | Ads‑b监视数据与雷达航迹的直觉模糊关联方法、装置 | |
CN110288835B (zh) | 一种基于运动学预测补偿机制的周边车辆行为实时识别方法 | |
Legrand et al. | Robust aircraft optimal trajectory in the presence of wind | |
CN111695299A (zh) | 一种中尺度涡轨迹预测方法 | |
Mahboubi et al. | Learning traffic patterns at small airports from flight tracks | |
CN109613530B (zh) | 一种低小慢空中目标多源信息融合的管控方法 | |
CN113283653B (zh) | 一种基于机器学习和ais数据的船舶轨迹预测方法 | |
CN115064009B (zh) | 一种终端区无人机与有人机冲突风险等级划分方法 | |
CN105701352A (zh) | 空间运动目标的轨迹预测方法 | |
CN117311393B (zh) | 一种无人机自主飞行路径规划方法及系统 | |
CN115147790A (zh) | 一种基于图神经网络的车辆未来轨迹预测方法 | |
Wang et al. | Data-driven conflict detection enhancement in 3d airspace with machine learning | |
Kant et al. | Long short-term memory auto-encoder-based position prediction model for fixed-wing UAV during communication failure | |
CN113470441B (zh) | 一种高机动试飞航空器实时智能防相撞检测方法 | |
KR20160036310A (ko) | 궤적 패턴을 이용한 항공기 도착시간 예측 장치 및 방법 | |
Graas et al. | Quantifying accuracy and uncertainty in data-driven flight trajectory predictions with gaussian process regression | |
Xu et al. | Aircraft trajectory prediction for terminal airspace employing social spatiotemporal graph convolutional network | |
Lowe et al. | Learning and predicting pilot behavior in uncontrolled airspace |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |