中国指挥与控制学会会刊 
随着全球定位技术和移动通信的快速发展,在执行飞行任务时,飞行器能够实时将轨迹数据传输到指挥所的态势系统中,而较高的采集频率使得飞行轨迹数据呈爆炸式增长,给轨迹数据存储、分析挖掘和可视化带来了巨大挑战[1]。据统计,飞行器飞行时1 s会采集4~6个轨迹点,飞行1小时约产生两万个轨迹点,而这些轨迹点通常包含大量的冗余数据,为降低轨迹数据存储和传输压力,提升轨迹数据分析和可视化效果,研究一种有效的飞行轨迹数据压缩算法是十分必要的。
目前,轨迹数据压缩算法有很多,考虑压缩效率和限制条件等因素,本文侧重于研究线段简化压缩算法,该类算法通过比较轨迹点最大误差与设定阈值或特定条件的大小来简化轨迹点[2]。Douglas-Peucker(DP)算法[3]是一种经典的线段简化压缩算法,该算法从全局轨迹点出发,使用垂直欧氏距离作为轨迹误差度量,由于DP算法需要进行递归操作,其时间复杂度为O(N2)。Meratnia等[4]在DP算法基础上考虑时间属性,提出了自顶向下的时间比压缩算法 (TD-TR),虽然有更高的压缩精度,但时间复杂度仍为O(N2)。后续,许多学者对DP算法进行了改进[5],针对DP算法存在阈值不确定问题,赵永清[6]提出了一种自动阈值的DP算法,张志伟等[7]提出了DP算法阈值选取优化方法。近年来,许多学者尝试在轨迹压缩时使用速度、方向等信息来提升压缩质量,Lin等[8]提出了保留轨迹速度特征的压缩算法,但在轨迹划分后仍使用DP进行递归;田智慧等[9]提出了一种同时考虑速度和角度的轨迹简化算法,增加了累积变向角阈值,有效减少了轨迹关键点的丢失。
由于飞行器飞行的特殊性,采集的轨迹点除了位置信息外,还包含高度、横滚角、俯仰角、航向角等信息,如表1所示。然而,只有少数学者考虑如何有效利用飞行轨迹特征进行数据压缩,邬鹏等[10]针对飞行轨迹特点,提出了一种基于飞行参数的分段压缩算法,在每个分段轨迹中,使用DP算法进行轨迹压缩;雷祥等[11]使用俯仰角和航向角的变化率作为阈值条件,提出了一种限定俯仰角和航向角的改进DP算法。虽然两种方法考虑了飞行参数特征,但存在阈值选取不确定、时间复杂度高等问题。针对上述问题,本文在有效利用飞行参数的基础上,选取方向变化率高的点作为轨迹拐点,通过最大限度保留携带信息多的拐点,提升轨迹数据压缩质量和效率。
表1 飞行轨迹点数据格式Tab.1 Flight trajectory point data format |
| 时间 | 经度 | 纬度 | 高度 | 横滚角 | 俯仰角 | 航向角 | 航速 |
|---|---|---|---|---|---|---|---|
| 1697161293 | 121.377 6 | 37.161 2 | 373 | 0.421 3 | 6.652 2 | 54.878 3 | 303.327 8 |
| 1697161293 | 121.377 6 | 37.161 2 | 373 | 0.414 6 | 6.638 4 | 54.452 6 | 304.754 9 |
| 1697161294 | 121.377 8 | 37.161 2 | 371 | -0.327 8 | 6.523 8 | 53.502 1 | 316.876 3 |
| 1697161295 | 121.377 9 | 37.161 2 | 370 | -0.156 1 | 6.286 7 | 53.304 2 | 310.372 3 |
| 1697161296 | 121.377 9 | 37.161 2 | 368 | -0.452 3 | 6.383 2 | 54.012 4 | 313.346 |
1 飞行轨迹数据压缩算法
1.1 符号定义
假设T={tid,fid,P1,P2,…,Pi,…PN}表示一条飞行轨迹,其中,tid表示该轨迹的唯一标识,fid表示产生该轨迹的飞机型号,Pi表示轨迹中第i个轨迹点,每个轨迹点包含时间、经度、纬度、高度、横滚角、俯仰角、航向角和航速8个参数,即Pi=(ti,xi,yi,hi,αi,βi,θi,vi),N为原轨迹点数,T*={tid,fid,P1,P2,…,Pi,…PM}表示压缩后的轨迹,M为压缩后的轨迹点数。
1.2 飞行轨迹数据处理
通过对实际飞行轨迹数据进行分析,我们发现轨迹中存在重复数据、异常数据和某一属性值缺失等情况,在数据压缩之前,应先对飞行轨迹数据进行预处理。
1.2.1 数据清洗转换
从态势数据中获取的飞行轨迹通常是分批次的,需要人工或通过算法进行批次合并,在这个过程中容易出现重复数据,因此,数据清洗首先应删除轨迹中唯一标识、时间相同的轨迹点;其次,在态势数据解码的过程中,会出现乱码或者属性为空值的问题,需删除轨迹点中飞机唯一标识、型号、经度、纬度等属性有明显错误的数据。
在飞行轨迹压缩过程中,可能需要计算两个轨迹点间的距离,而态势数据中是以经纬度来表示飞机的位置,本文采用墨卡托投影方法将其转换到直角坐标系。
1.2.2 异常点检测
一般情况下,飞机的位置由时间、航速等属性决定,由于通信设备异常、信号干扰等因素会导致飞行轨迹中出现跳点,轨迹点位置异常意味着速度异常,因此,跳点检测可以采用如下公式:
$ v_{i}>v_{\max }$
其中,vi表示轨迹中第i个点的速度,vmax表示当前环境下飞机飞行理论上的最大速度。通常,飞机在飞行中高度不会从一个具体数值突然变为0,只有在飞机降落时高度值才为0,而在分析飞行轨迹数据时,人们发现存在很多飞行过程中高度为0的轨迹点,本文也将这些认定为异常轨迹点。
1.2.3 数据插值修复
在检测出轨迹异常点后,最简单的处理办法是将这些点删除。然而,直接删除异常点可能导致部分轨迹间隔较大且不均匀,不利于后续的数据压缩操作。因此,本文采用插值技术对异常点进行修复。
针对位置异常点,本文采用的是结合航向和航速的插值方法[12]。该方法通过异常点的前后两个点分别计算异常点位置的预测值。然后,通过轨迹点之间的时间属性计算权重。最终,采用加权平均的方式得到异常值的新位置。针对高度异常点,本文利用线性插值方式预测新的高度属性值。
1.3 飞行轨迹数据压缩算法
很多研究表明保留轨迹方向的压缩方法优于保留轨迹位置的压缩方法[13],本文在充分利用轨迹方向特征的前提下,提出了一种有效的飞行轨迹数据压缩算法。
1.3.1 轨迹拐点提取
飞行轨迹的方向参数包含横滚角、俯仰角、航向角,为提取具有代表性的方向拐点,本文采用聚类的方法对方向参数进行类别划分。成熟的聚类方法有很多,考虑聚类效率,本文采用经典的k-means聚类方法,聚簇数目k值使用肘部法来确定。
为保证飞行轨迹起止位置准确,我们需提前去除轨迹起点和终点。设n=N-2,轨迹点Pi的方向特征向量为(αi,βi,θi),方向特征集H可表示为
$ \boldsymbol{H}=\left[\begin{array}{ccc}\alpha_{1} & \beta_{1} & \theta_{1} \\\alpha_{2} & \beta_{2} & \theta_{2} \\\vdots & \vdots & \vdots \\\alpha_{n} & \beta_{n} & \theta_{n}\end{array}\right]=\left[\boldsymbol{H}_{i j}\right]_{n \times 3}$
为避免数据数值域大小给聚类带来影响,本文对矩阵H进行标准化处理,计算公式如下:
$ \boldsymbol{H}_{i j}^{*}=\frac{\boldsymbol{H}_{i j}-\overline{\boldsymbol{H}_{j}}}{\sigma_{j}}$
$ \sigma_{j}=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left(\boldsymbol{H}_{i j}-\overline{\boldsymbol{H}_{j}}\right)^{2}}$
$ \overline{\boldsymbol{H}_{j}}=\frac{1}{n} \sum_{i=1}^{n} \boldsymbol{H}_{i j}$
其中,σj、 分别为第j个分量的标准差和平均值。对方向特征向量k-means聚类后,每个轨迹点会获得一个类别标签。将原始轨迹映射为一个标签序列。选择标签变化处的轨迹点为轨迹拐点,如图1所示,一条轨迹含有3个不同类别,红色标记的轨迹点为拐点。
1.3.2 综合得分计算
由于飞机在飞行过程中方向变化频繁,上一节中得到的轨迹拐点数量会非常多,人们需要采用计算轨迹拐点的综合得分来进一步筛选。综合得分通过轨迹拐点的位置、高度和航速等指标来计算。
首先,定义拐点Pi到相邻轨迹点Pi-1和Pi+1位置距离、高度差、航速差的计算公式,具体如下:
$ d_{i}=\frac{\left|\left(y_{i+1}-y_{i-1}\right) x_{i}+y_{i}\left(x_{i+1}-x_{i-1}\right)+y_{i+1} y_{i-1}+x_{i+1} x_{i-1}\right|}{\sqrt{\left(x_{i+1}-x_{i-1}\right)^{2}+\left(y_{i+1}-y_{i-1}\right)^{2}}}$
$ \mu_{i}=\left|h_{i-1}-h_{i}\right|+\left|h_{i}-h_{i+1}\right|$
$ \varphi_{i}=\left|v_{i-1}-v_{i}\right|+\left|v_{i}-v_{i+1}\right|$
图2展示了拐点Pi位置距离误差计算示例。
其次,为避免各指标间的量纲差异,需对指标进行归一化处理,具体公式如下:
$ D_{i j}^{*}=\frac{D_{i j^{-}} \min _{i=1, \ldots, p} D_{i j}}{\max _{i=1, \ldots, p} D_{i j^{-}} \min _{i=1, \ldots, p} D_{i j}}$
式中, 为第i个拐点的第j项指标的归一化数值,而i=1,…,p为拐点的数量,j对应位置、高度和航速三个指标,Dij为第i个拐点的第j项指标值。
接着,计算各指标的熵值hj:
$ h_{j}=-\frac{1}{\ln (p)} \sum_{i=1}^{p} f_{i j} \ln \left(f_{i j}\right)$
$ f_{i j}=\frac{\boldsymbol{D}_{i j}^{*}}{\sum_{i=1}^{p} \boldsymbol{D}_{i j}^{*}}$
其中,fij为第i项指标下第j个数据点占该指标的比重。
最后,轨迹拐点的综合得分si为:
$ s_{i}=\sum_{j=1}^{3} w_{j} f_{i j}$
$ w_{j}=\frac{1-h_{j}}{\sum_{j=1}^{3} h_{j}}$
其中,wj为第j个指标的熵值权重。将各个拐点的综合得分按照降序排序,选取前M-2个拐点,加上轨迹起点、轨迹终点组成压缩后的轨迹点。算法流程如表2所示。
表2 飞行轨迹数据压缩算法Tab.2 Flight trajectory data compression algorithm |
| 输入:原始轨迹T={P1,P2,…,PN},压缩后的轨迹点数M,聚簇数目k |
| 输出:压缩后的轨迹T*={P1,P2,…,PM} |
| 1:构造方向特征矩阵H,按照公式(3)进行标准化处理; |
| 2:随机挑选k个方向特征向量为初始聚类中心; |
| 3:使用k-means方法对方向特征进行聚类,得到各个轨迹点的类别标签; |
| 4:将轨迹点映射为标签序列,标签变化处的前后两个点作为轨迹拐点,组成轨迹拐点集合; |
| 5:遍历轨迹拐点集合,根据公式(6)—(8)计算各个拐点的位置、高度和航速的指标值; |
| 6:根据公式(9)—(13)计算各个拐点的综合得分; |
| 7:取得分最高的M-2个轨迹拐点、轨迹起点和轨迹终点作为最终压缩后的轨迹T*。 |
本算法的时间复杂度取决于轨迹拐点提取和综合得分计算两个步骤,轨迹拐点提取的时间复杂度由k-means算法决定,为O(lkN),其中,l为迭代次数,k为聚簇数,N为轨迹点数;综合得分计算的时间复杂度取决于计算轨迹拐点各个指标的循环操作,均为O(p),其中,p为轨迹拐点数,且p<N。综合上述两个步骤,本算法的时间复杂度为O(lkN)。
2 实验结果与分析
2.1 实验环境及数据来源
实验硬件配置为Intel i7-1160G7处理器,16 GB内存; 实验环境为Windows 10 操作系统,算法开发语言为Python 3.10。
实验数据来自态势系统中飞行器的飞行轨迹,提取一段飞行时间为30分钟的轨迹作为最终的实验数据,包含7392个轨迹点。实验对比算法为DP算法和TD-TR算法。
2.2 性能度量指标
本文选取运行时间、压缩率、平均距离误差和平均方向误差作为度量指标。压缩率计算公式为
$ r=(N-M) / N \times 100 \%$
平均距离误差和平均方向误差分别表示压缩后的轨迹在欧氏距离和方向上的损失,计算公式分别为:
$ e_{p}=\frac{1}{N} \sum_{i=1}^{N}\left(p_{T\left(P_{i}\right)}-p_{T^{*}\left(P_{i}\right)}\right)$
$ e_{d}=\frac{1}{N} \sum_{i=1}^{N}\left(d_{T\left(P_{i}\right)}-d_{T^{*}\left(P_{i}\right)}\right)$
2.3 实验结果及分析
表3给出了在不同压缩率下3种算法的运行时间,从表中可以看出,压缩率为20%时,本文算法的运行时间略优于其他两种算法,随着压缩率不断增大,本文算法的优势逐渐显现,特别是在压缩率为80%时,本文算法的运行时间明显低于其他两种算法。主要原因是DP算法和TD-TR算法都是通过递归实现的,时间复杂度为O(N2),而本文算法时间复杂度为O(lkN)。
表3 在不压缩率下的算法执行时间Tab.3 Algorithm execution time in different compression ratio 单位:s |
| 算法 | 压缩率/% | ||||
|---|---|---|---|---|---|
| 20 | 40 | 50 | 60 | 80 | |
| DP | 0.151 | 0.984 | 1.324 | 2.391 | 5.356 |
| TD-TR | 0.136 | 0.845 | 1.421 | 2.781 | 5.049 |
| 本文算法 | 0.104 | 0.823 | 0.971 | 1.705 | 3.498 |
图3对3种算法在不同压缩率下平均距离误差进行对比,整体来看,3种算法的距离误差随着压缩率的增加而提高,TD-TR算法在距离误差方面略优于DP算法,而使用本文算法进行压缩产生的距离误差在不同压缩率下都低于其他两种算法,说明本文算法提取的轨迹拐点精度更高,能用相当的轨迹点保留更多的位置信息。
图4比较了不同压缩率下3种算法的平均方向误差,从图中可以看出:压缩率低于40%时,本文算法的方向误差较低且变化较小;压缩率为30%时,本文算法的方向误差为TD-TR算法的1/3;压缩率高于70%时,本文算法与其他两个算法相比,方向误差均最小。其主要原因是本文算法在选取轨迹拐点时,充分考虑了飞行方向参数变化因素,而DP算法和TD-TR算法只是简单进行了距离计算。
在实验过程中,DP算法和TD-TR算法要不断调整阈值大小,以确定相应压缩率下的特定阈值,而本文算法只需给出压缩后的轨迹点数量,不需要人为设定阈值,因此在实际应用中更灵活、高效。
3 算法应用
前期,我们基于GIS服务实现了目标与态势三维可视化系统[14],在系统使用过程中,随着飞机飞行次数和时间的增加,产生的态势数据越来越多,严重影响了飞行轨迹可视化的效果。本文算法对飞行轨迹数据进行有效压缩,在态势系统中的应用体现在轨迹回放和轨迹包络可视化两个方面。
3.1 轨迹回放
在大型模拟训练后,通常需要基于态势回放系统进行复盘总结,飞行轨迹回放是将所有轨迹点按照时间顺序进行三维可视化,而飞机在飞行过程中会频繁采集轨迹数据,飞行时间超过1小时可达上万个轨迹点,给态势回放带来了巨大挑战。在实际应用中,如果飞机数量达到20架以上,态势回放功能要加载60 s左右,严重影响用户的使用体验。本文算法通过保留最优轨迹拐点的压缩方式,不仅将飞机航迹压缩至指定长度,也保留原轨迹的方向、位置等重要特征。使用本文算法后,20架飞机态势回放加载时间缩短至15 s左右,并且压缩后的航迹跟原航迹的可视化效果基本一致。
3.2 轨迹包络可视化
飞机执行任务后,根据飞行轨迹和雷达探测范围,利用线缓冲区生成算法[15]在态势展示系统中生成飞行轨迹包络图,能够直观呈现飞机此次任务的侦察区域。然而,一次飞行任务持续时间长,产生轨迹点数量巨大,在轨迹展示和包络范围计算时消耗资源多,直接导致了此功能效率变低,时常出现卡顿,严重影响了用户的使用体验。使用本文算法可在后端服务将飞行轨迹压缩到指定轨迹点数,方便前端服务调用,在保证压缩质量的前提下,有效降低了轨迹点的数量。在实际应用中,可将飞行轨迹包络可视化功能的效率平均提升3倍,且原飞行轨迹点越多,提升效果越明显。
4 结束语
本文针对飞行轨迹数据压缩算法未有效利用飞行参数和阈值设定难等问题,提出了一种有效的飞行轨迹数据压缩算法。算法对每条轨迹的方向特征向量进行聚类,选取方向变化的点作为轨迹拐点,通过引入信息熵理论,按照各个轨迹拐点的综合得分排序,确定最终压缩后的轨迹点。实验结果表明,算法压缩后的轨迹点精度更高且运行时间优于其他算法。最后,本文将算法应用到态势展示系统的飞行轨迹回放和轨迹包络可视化功能中,应用效果显著。