CN110060513A - 基于历史轨迹数据的空中交通管制员工作负荷评估方法 - Google Patents

Abstract

本发明公开了一种基于历史轨迹数据的空中交通管制员工作负荷评估方法，包括以下步骤：历史飞行轨迹的获取与处理；遍历轨迹集合中各个航空器的历史飞行轨迹，对轨迹进行线性拟合；推算可能产生的管制指令、协调信息和管制监控航空器架次；计算总的工作时间即工作负荷，完成空中交通管制员工作负荷评估。本发明从航空器历史轨迹入手，通过分析历史航迹中航空器航行诸元变化的趋势，推算管制指令类型，并根据管制指令次数、发布时机、实际飞行流量分布及变化，量化管制员管制指挥工作业务量，评估单位时间内管制员工作负荷。

Description

基于历史轨迹数据的空中交通管制员工作负荷评估方法

技术领域

本发明涉及一种空中交通管制员工作负荷评估方法，尤其涉及一种基于历 史轨迹数据的空中交通管制员工作负荷评估方法。

背景技术

空中交通管制员是空中交通管理直接参与者，通过通信、导航及监视设备 掌握和预测管制空域内的交通态势，以对空中交通运行实施实时、科学的管制 指挥，确保空域内交通安全、有序和高效运行。随着飞行流量递增，管制员工 作负荷也随之上涨甚至使管制员处于超负荷工作状态，导致工作效率骤降，威 胁航空安全。根据统计数据表明，近年以来，由空中交通管制原因所造成的空 中交通事故征候呈逐年上升趋势。2016年10月11日，由于监控不到位，管制 员错发指令导致虹桥机场发生A类跑道侵入事件，为空中交通运行安全敲响警 钟。因此面向管制员技能水平、管制员工作强度量化的相关研究成为近年来民航工作的重点。

目前，管制员工作负荷相关数据方式的差异性可将评定方法分为三大类： 主观评定法、客观评定法、相关指标评定法。

主观评定法主要运用心理学以设计、制作问卷，对空中交通管制员工作负 荷的主观评定进行引导，以问卷调查的方式结合使用李特式量表(常用5点式 或7点式)予以实施，主要包括：空中交通负荷输入技术(Air traffic workload input technique,ATWIT)；主观工作负荷评估技术(Subjective Workload Analysis Technique,SWAT)；NASA的工作负荷指数(Task Load Index,TLX)；库柏-哈柏 修正法(Modified Cooper-Harper Ratings,MCH)。

客观评定法主要在于使用一些可以量化的客观参数表征空中交通管制员一 段时间内的工作负荷，最为主要的可量化客观参数即为时耗，通过对管制员执 行某项管制指挥任务(如：填写进程单、发布管制许可、协调等)的时耗进行 独立重复统计，使用时耗大小反映空中交通管制员工作负荷，该类方法使用的 是客观数据，客观性显著，现存较为成熟的方法包括：英国运筹与分析理事会 的预测管制员工作负荷法(The method developed byCAA’sdirectorate of operational research and analysis,DORATASK)；波_布分析法(Messerschmitt MBB)。

(3)相关指标评定法，使用特定指标测评管制员工作负荷的方法，如工作表 现指标评估法以及生理、心理指标评估法等。

由于管制员对空指挥工作和协调工作等是通过无线电陆空通话向航空器驾 驶员发布管制指令的方式完成。因此目前第二类客观评定法被广泛使用。即将 单位时间内管制员的陆空无线电通话时间总量作为管制员负荷的表征指标。因 此我国空管系统将管制工作负荷定义为单位时间管制工作量，通常以单位小时 内的管制工作时间进行衡量。由此中国国内相关研究多以DORATASK技术为基 础，将管制工作分类，采集并分析管制通话录音，计时量化通话中的指令及协 调信息发布时间，并累加统计得到管制工作符合。此外，也有部分学者采用利 用飞行架次、进出港航班量等飞行流量统计信息，结合空域结构，构造回归模 型预测或推算管制员工作负荷。

管制工作负荷评估是航空人因领域的热点问题，其评估结构能用于评价管 制员工作强度，还可用于空域容量评估、航路布局优化以及管制员培训和管理 的强化等方面。

主观评定法需要管制员从自身疲劳度、工作强度和心理压力等方面评价负 荷等级，有可能由于监控不到位意识不到飞行冲突，导致自我评价与实际情况 不符。同时，该方法实施时会在管制员工作时打断工作流程，实时获取评价数 据，某种程度上会干扰管制员监控与思考，增加了工作难度，导致管制负荷增 加。

相关指标评定法需要基于丰富的运行经验制定工作指标评定表，而生理、 心理指标评估法要求装设完备、先进的生理、心理测试仪器，因此，该类方法 推广难度系数较高。

客观法相关研究及应用较多，但目前研究成果存在以下缺点：

要量化管制指令时间，需要现场采集的管制员陆空通话，对管制话音计时 统计。目前这部分工作主要通过人工收听通话，区分飞行员发话与管制员指令， 并对管制指令计时，该过程需要耗费大量人力与时间。由于管制通话专业性强， 一个管制员对多个飞行员发布指令、飞行情报，因此只有了解管制业务的人员 才能分析通话内容，区分指令、飞行情报与飞行员请求等不同信息。由于信号 干扰等原因，大部分通话录音存在噪声干扰，因此人工识别有一定困难，效率 不高。

根据DORATASK方法，管制员工作分为看得见与看不见两种，目前方法及 方案多对看得见的，例如指令发布、进程单操作等看得见工作进行统计，而管 制员监控态势、冲突预判等看不见的工作难以统计工作量。

管制员会与航空器驾驶员一次通话中包含多个不同指令信息。目前工作负 荷计算方法计算发指令的总时间，并未详细区分不同类型指令对管制员工作负 荷影响。

由于大流量空域内雷达引导较多，航空器按照管制员临时安排的路线飞行， 严格按照程序飞行的比例较小，由于忽略了实际飞行轨迹与计划航迹的差异， 利用飞行流量和空域结构推算管制工作负荷的方法会导致较大误差。

发明内容

本发明的目的就在于为了解决上述问题而提供一种基于历史轨迹数据的空 中交通管制员工作负荷评估方法。

本发明通过以下技术方案来实现上述目的：

一种基于历史轨迹数据的空中交通管制员工作负荷评估方法，包括以下步 骤：

步骤1：历史飞行轨迹的获取与处理：选择关注的管制空域及时间范围，获 取该时间段内该空域的运行数据，获取航空器的历史飞行轨迹；将历史飞行轨 迹中关键字段进行对比，按照时间顺序产生不同航空器的历史飞行轨迹，形成 轨迹集合；按照各历史飞行轨迹的运行方向和高度、速度变化，判断该航班为 起飞或者降落航班；同时对每个航班在空域内的飞行距离、飞行时间进行统计；

步骤2：遍历轨迹集合中各个航空器的历史飞行轨迹，对轨迹进行线性拟合， 分析轨迹中的速度、高度、航向随时间变化的趋势，得到高度、速度、航向突 然产生变化的时刻及位置，同时对航空器首次出现时刻位置、最后出现时刻位 置、进港航空器到达起始进近定位点并到达起始进近进近高度的时刻进行记录；

步骤3：利用步骤2的记录数据，推算可能产生的管制指令、协调信息和管 制监控航空器架次；

步骤4：根据步骤3得到的推算结果计算总的工作时间即工作负荷，完成空 中交通管制员工作负荷评估。

作为优选，所述步骤1中，所述运行数据是来自空管自动化设备或空管记 录仪中的监视数据，包括但不限于二次航管雷达数据、ADS-B数据即广播式自 动相关监视数据、融合后的系统航迹数据；所述形成轨迹集合的方法是：按行 读取运行数据的每一条记录，查询容器中是否存在以该呼号为唯一标识的航空 器对象，如果没有则新建航空器对象并加入容器，否则读取该条记录中与航空 器起降相关的时间、位置、速度、使用跑道信息到航迹点对象中，并将航迹点 数据加入航空器对象的航迹链表，最后在容器中得到以时间排序并且以呼号为 标识的航空器对象序列，形成航空器航迹集合；所述判断该航班的起降类型包括以下步骤：(1)同时满足下述三个条件被判断为离港航班即起飞航班：航空 器第一个航迹点在机场跑道中心点2km范围内，且高度低于1000m；航空器最 后一个航迹点高度、速度大于第一个航迹点高度、速度；航空器航迹点个数大 于50个；(2)同时满足下述三个条件被判断为进港航班即降落航班：航空器最 后一个航迹点在机场跑道中心点2km范围内，且高度低于1000m；航空器最后 一个航迹点高度、速度小于第一个航迹点高度、速度；航空器航迹点个数大于 50个；(3)若航空器航迹不满足上述两种条件，该航迹太短，或者数据异常无法判断，则从航迹集合中去除该航空器轨迹。

所述步骤2中，所述对轨迹进行线性拟合包括以下步骤：

步骤(2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点按 照p_i(t)升序排序；

步骤(2.2)：从p_i的第一个p_i,0元素开始，初始化第一条轨迹线段lf_i,0，令 lf_i,0(n_s)＝0，求两轨迹点连线的夹角angle_i,j，令lf_i,0(a)＝angle_i,0；

步骤(2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ＝0；

步骤(2.4)：设lf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(2.5)：求从k至k+1的轨迹点连线的夹角angle_i,k，若|angle_i,k-lf_i,λ(a)|＞θ轨迹出现拐点，若λ＝0，转至步骤(2.7)，若λ>0，转至步骤(2.6)，若|angle_i,k-lf_i,λ(a)|≤θ，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角 转至步骤(2.9)；

步骤(2.6)：判断当前路径与上一条路径的夹角，若转至步骤(2.7)；若转至步骤(2.8)；

步骤(2.7)：当前路径结束，p_i,k为线段lf_i,λ的结束点lf_i,λ(n_e)＝k，求lf_i,λ的其 他元素值，令λ＝λ+1增加一条新的轨迹线段，lf_i,λ(n_s)＝k，以p_i,k为开始点记录 新的线段，转至步骤(2.9)；

步骤(2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点 f_i,λ-1(n_e)＝k，更新lf_i,λ-1的其他元素值，令lf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(2.9)；

步骤(2.9)：令k＝k+1，重复步骤(2.4)，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令lf_i,λ(n_e)＝n，lf_i,λ(n_s)＝n，求lf_i,λ的其他元素值；

步骤(2.10)：将lf_i,0,lf_i,1,...,lf_i,λ写入集合LF_i；

步骤(2.11)：重复步骤(2.2)-(2.7)，直至遍历完F中所有航班，将 LF₀,LF₁,...,LF_n写入集合LF。

所述步骤3中，所述管制指令包括雷达引导指令、高度指令、速度指令；

所述雷达引导指令的推算方法包括以下步骤：

步骤(3.1.1)：标准飞行程序转弯点拟合，包括以下步骤：

步骤(3.1.1.1)：取F中航班f_i的轨迹线段集合LF_i，λ为轨迹线段序号，令 λ＝0；

步骤(3.1.1.2)：若λ≤num(LF_i),设lf_i,λ的起点为n_s，若λ＞num(LF_i)，设lf_i,λ的 终点为n_s，遍历S中的各个元素，求p_i,ns与S中的各个元素的距离，若转至步骤(3.1.1.3)，否则转至步骤(3.1.1.4)；

步骤(3.1.1.3)：线段起点与标准飞行程序转弯点距离小于门限值的一半， 用标准飞行程序转弯点修正历史飞行轨迹拐点，若令并 更新s_j的经纬度λ＝λ+1，转至步骤(3.1.1.2)；

步骤(3.1.1.4)：创建一个标准飞行程序中路径点迹集合S，令s.long＝p_{i, ns}.long； s.lat＝p_i,ns.lat；F′＝{f_i}，令s∈S。λ＝λ+1，若λ＞num(LF_i)，转至步骤(3.1.1.5)， 否则转至步骤(3.1.1.2)；

步骤(3.1.1.5)：依次遍历轨迹线段集合LF中所有元素，直至结束；

步骤(3.1.1.6)：依次遍历S中元素，计算标准程序点s_j的利用率δ_j，若 δ_j＞μ，保留该关键点，最终得到标准飞行程序转弯点，否则去除该点，令 S＝S-{s_j}，重复本步骤；

步骤(3.1.2)：将标准飞行程序转弯点与历史航空器飞行轨迹比较，将历史 航空器飞行轨迹拐点集合中与标准飞行程序转弯点集合相交的点去除后重新更 新航空器轨迹拐点集合，集合中的元素为航空器飞行过程中机动飞行产生的航 向变化点，机动拐点个数即为管制员为该航空器发布雷达引导指令的个数；

所述高度指令的推算方法包括以下步骤：

步骤(3.2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点 并按照p_i(t)升序排序；

步骤(3.2.2)：从p_i的第一个元素p_i,0开始，初始化高度剖面第一条轨迹线段 hf_i,0及轨迹线段起止点的高度值p_i,k(h)和经度坐标值p_i,k(long)；

步骤(3.2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号，令k＝1， λ＝0；

步骤(3.2.4)：设hf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(3.2.5)：求从k至k+1的轨迹点连线的高度差hf_i,k(Δh)，若hf_i,k(Δh)≥60轨迹出现拐点，若λ＝0，转至步骤(3.2.7)，若λ>0，转至步骤(3.2.6)，若 hf_i,k(Δh)＜60，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角 转至步骤(3.2.9)；

步骤(3.2.6)：判断当前路径与上一条路径的夹角，若 |hf_i,λ-1(Δh)-hf_i,λ(Δh)|＞30，转至步骤(3.2.7)；若|hf_i,λ-1(Δh)-hf_i,λ(Δh)|≤30，转至 步骤(3.2.8)；

步骤(3.2.7)：当前路径结束，p_i,k为线段hf_i,λ的结束点hf_i,λ(n_e)＝k，令λ＝ λ+1增加一条新的高度剖面轨迹线段，hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线 段，转至步骤(3.2.9)；

步骤(3.2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点 f_i,λ-1(n_e)＝k，令hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转 至步骤(3.2.9)；

步骤(3.2.9)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令hf_i,λ(n_e)＝n,hf_i,λ(n_e)＝k；

步骤(3.2.10)：将hf_i,0,hf_i,1,...,hf_i,λ写入集合HF_i；

步骤(3.2.11)：重复步骤(3.2.2)-(3.2.7)，直至遍历完F中所有航班，将 HF₀,HF₁,...,HF_n写入集合HF；

步骤(3.2.11)：集合HF_i中元素个数即为航空器在空域内飞行时高度指令个 数。

所述速度指令的推算方法包括以下步骤：

步骤(3.3.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点 并按照p_i(t)升序排序；

步骤(3.3.2)：从p_i的第一个元素p_i,0开始，初始化速度剖面第一条轨迹线段 sf_i,0及轨迹线段起止点的高度值p_i,k(s)和经度坐标值p_i,k(long)；

步骤(3.3.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ ＝0；

步骤(3.3.4)：设sf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(3.3.5)：求从k至k+1的轨迹点连线的斜率sangle_i,j，若sangle_i,j≠0轨迹 出现拐点，若λ＝0，转至步骤(3.3.7)，若λ>0，转至步骤(3.3.6)，若sangle_i,j＝0， 轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角转至步骤 (3.3.9)；

步骤(3.3.6)：判断当前路径与上一条路径的斜率，若 |sf_i,λ-1(slope)-sf_i,λ(slope)|＝0，转至步骤(3.3.7)；若转至步骤(3.3.8)；

步骤(3.3.7)：当前路径结束，p_i,k为线段sf_i,λ的结束点sf_i,λ(n_e)＝k，令λ＝ λ+1增加一条新的速度剖面轨迹线段，sf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线 段，转至步骤(3.3.8)；

步骤(3.3.8)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令sf_i,λ(n_e)＝n，sf_i,λ(n_e)＝k，转至步骤(3.3.9)；

步骤(3.3.9)：将sf_i,0,sf_i,1,...,sf_i,λ写入集合SF_i；

步骤(3.3.10)：重复步骤(3.3.2)-(3.3.7)，直至遍历完F中所有航班，将 SF₀,SF₁,...,SF_n写入集合SF；

步骤(3.3.11)：SF_i集合中元素个数即为航空器在空域内飞行时速度指令个 数。

所述步骤3中，所述管制指令还包括但不限于以下其他指令或飞行情报： 雷达识别、进离场程序通报、气象条件通报、进近许可、脱波，这些指令在管 制过程中仅为航空器提供一次指令，故推算一架航空器以上类型的指令各发布 一次。

所述步骤4中，按以下公式对中交通管制员的工作负荷进行评估：

其中，为单位时间t内的管制员工作负荷量，单位为秒，σ_i为单位时 间t内管制服务i的次数，t_i表示管制员为提供管制服务i的工作时间。

本发明的有益效果在于：

本发明从航空器历史轨迹入手，通过分析历史航迹中航空器航行诸元变化 的趋势，推算管制指令类型，并根据管制指令次数、发布时机、实际飞行流量 分布及变化，量化管制员管制指挥工作业务量，评估单位时间内管制员工作负 荷；具体优势如下：

1、利用计算机进行轨迹海量处理与挖掘，省略了人工收听通话录音及人工 区分并计时的过程，可节省人力与时间成本；

2、分析历史轨迹可得到航空器航行诸元变化的准确时刻，避免由于通话噪 音或干扰产生的指令统计误差；

3、按照不同的类型解析航迹、推算指令，因此可针对不同管制指令进行分 类统计和评估，比传统的进统计总的指令发布时间更细致，为后期分析不同指 令的影响等工作打下基础；

4、除了管制指令之外，通过飞行轨迹处理与挖掘还能统计飞行流量等信息， 用于量化管制员监控工作量和协调工作量，使工作负荷评估结果更加准确。

附图说明

图1是本发明实施例某日成都双流机场进近空域飞行流量统计显示界面示 意图；

图2是本发明实施例单一航迹点集合选取及4D信息显示示意图；

图3是本发明实施例航空器历史飞行线性轨迹拟合示意图；

图4是本发明实施例根据航迹拟合结果推算管制雷达引导服务位置及次数 示意图。

图5是本发明实施例根据航迹拟合结果推算管制指令高度变更次数及位置 示意图。

图6是本发明实施例根据航迹拟合结果推算管制指令速度变更次数及位置 示意图。

具体实施方式

下面结合实施例和附图对本发明作进一步说明：

实施例：

一种基于历史轨迹数据的空中交通管制员工作负荷评估方法，总体包括以 下四个步骤：

步骤1：历史飞行轨迹的获取与处理：选择关注的管制空域及时间范围，获 取该时间段内该空域的运行数据，获取航空器的历史飞行轨迹；将历史飞行轨 迹中关键字段进行对比，按照时间顺序产生不同航空器的历史飞行轨迹，形成 轨迹集合；按照各历史飞行轨迹的运行方向和高度、速度变化，判断该航班为 起飞或者降落航班；同时对每个航班在空域内的飞行距离、飞行时间进行统计；

步骤2：遍历轨迹集合中各个航空器的历史飞行轨迹，对轨迹进行线性拟合， 分析轨迹中的速度、高度、航向随时间变化的趋势，得到高度、速度、航向突 然产生变化的时刻及位置，同时对航空器首次出现时刻位置、最后出现时刻位 置、进港航空器到达起始进近定位点并到达起始进近进近高度的时刻进行记录；

步骤3：利用步骤2的记录数据，推算可能产生的管制指令、协调信息和管 制监控航空器架次；

步骤4：根据步骤3得到的推算结果计算总的工作时间即工作负荷，完成空 中交通管制员工作负荷评估。

下面对上述各步骤的优选方法进行具体说明。

所述步骤1中，所述运行数据是来自空管自动化设备或空管记录仪中的监 视数据，包括但不限于二次航管雷达数据、ADS-B数据即广播式自动相关监视 数据、融合后的系统航迹数据；所述形成轨迹集合的方法是：按行读取运行数 据的每一条记录，查询容器中是否存在以该呼号为唯一标识的航空器对象，如 果没有则新建航空器对象并加入容器，否则读取该条记录中与航空器起降相关 的时间、位置、速度、使用跑道信息到航迹点对象中，并将航迹点数据加入航 空器对象的航迹链表，最后在容器中得到以时间排序并且以呼号为标识的航空 器对象序列，形成航空器航迹集合；所述判断该航班的起降类型包括以下步骤： (1)同时满足下述三个条件被判断为离港航班即起飞航班：航空器第一个航迹 点在机场跑道中心点2km范围内，且高度低于1000m；航空器最后一个航迹点 高度、速度大于第一个航迹点高度、速度；航空器航迹点个数大于50个；(2) 同时满足下述三个条件被判断为进港航班即降落航班：航空器最后一个航迹点 在机场跑道中心点2km范围内，且高度低于1000m；航空器最后一个航迹点高 度、速度小于第一个航迹点高度、速度；航空器航迹点个数大于50个；(3)若 航空器航迹不满足上述两种条件，该航迹太短，或者数据异常无法判断，则从 航迹集合中去除该航空器轨迹。

举例来说，图1是某日成都双流机场进近空域飞行流量统计显示功能界面， 利用此界面上的信息可以将其进离场航空器飞行航迹信息从历史监视数据记录 文件里提取出来，对其航行信息进行统计，如获取航空器呼号、二次雷达应答 机编码、机型、尾流类型、经度、纬度、起降时间、进近时间、飞行距离和降 落次序等信息。本文提取航空器历史雷达数据之后将用于轨迹显示、飞行流量 统计、识别管制业务等。

图2是单一航空器对象航迹点集合示意图，此航迹点集合的获取是将监视 系统中进场、离港、飞越航班的航迹单独分离出来而输出形成的航迹点集合， 从中可以清楚的掌握每一个航班的过点时间，其关键点和拐点的经度、纬度、 高度和航向的信息进而可以对其做飞行轨迹拟合并分析掌握其飞行性质。

所述步骤2中，所述对轨迹进行线性拟合包括以下步骤：

步骤(2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点按 照p_i(t)升序排序；

步骤(2.2)：从p_i的第一个p_i,0元素开始，初始化第一条轨迹线段lf_i,0，令 lf_i,0(n_s)＝0，求两轨迹点连线的夹角angle_i,j，令lf_i,0(a)＝angle_i,0；

步骤(2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ＝0；

步骤(2.4)：设lf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(2.5)：求从k至k+1的轨迹点连线的夹角angle_i,k，若|angle_i,k-lf_i,λ(a)|＞θ轨迹出现拐点，若λ＝0，转至步骤(2.7)，若λ>0，转至步骤(2.6)，若 |angle_i,k-lf_i,λ(a)|≤θ，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角 转至步骤(2.9)；

步骤(2.6)：判断当前路径与上一条路径的夹角，若转至步骤(2.7)；若转至步骤(2.8)；

步骤(2.7)：当前路径结束，p_i,k为线段lf_i,λ的结束点lf_i,λ(n_e)＝k，求lf_i,λ的其 他元素值，令λ＝λ+1增加一条新的轨迹线段，lf_i,λ(n_s)＝k，以p_i,k为开始点记录 新的线段，转至步骤(2.9)；

步骤(2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点 f_i,λ-1(n_e)＝k，更新lf_i,λ-1的其他元素值，令lf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(2.9)；

步骤(2.9)：令k＝k+1，重复步骤(2.4)，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令lf_i,λ(n_e)＝n，lf_i,λ(n_s)＝n，求lf_i,λ的其他元素值；

步骤(2.10)：将lf_i,0,lf_i,1,...,lf_i,λ写入集合LF_i；

步骤(2.11)：重复步骤(2.2)-(2.7)，直至遍历完F中所有航班，将 LF₀,LF₁,...,LF_n写入集合LF。

上述步骤2的基本思想是对航空器进近飞行航迹点集合P中的各个元素进 行遍历，依次计算出相邻两航班轨迹点形成的线段。再比较相邻线段的角度差， 如果角度差大于门限值，则认为轨迹产生拐点，计算上一条路径的终点，和下 一条路径的节点，如此往复，直到遍历完所有航路点，可将有限的离散航迹点 转化成为线性轨迹线段。相关变量定义描述如表1所示：

表1历史轨迹点线性拟合算法变量定义

航空器历史轨迹线性拟合相关算法如下：

将相邻的两航班轨迹前后两点的经纬度坐标转化为直角坐标，计算可得两 轨迹点p_i,j和p_i,j+1组成的线段的斜率并用角度angle_i,j来表示。由数学知识知，斜 率angle_i,j大小范围在[0,π]之间，因此将angle_i,j转换到[0,π]中。

0≤angle_i,j≤π

若angle_i,j＝π-angle_i,j

两轨迹点p_i,j和p_i,j+1之间的距离fdis(p_i,j,p_i,j+1)计算过程如下：

设p1,p2为地理坐标上的两个点，均包含经纬度两个元素{long,lat}，

令p1.lat＝p_i,j.lat；p1.long＝p_i,j.long

p2.lat＝p_i,j+1.lat；p2.long＝p_i,j+1.long

其中，r＝6378.137，单位公里。

通过斜率角度判断，可将轨迹中的点分段，各轨迹线段由轨迹点拟合产生。 采用最小二乘法将航迹点数据拟合成飞行轨迹。

设拟合线段直角坐标系的直线方程为L＝ax+by+c，即L＝{a,b,c}

定义拟合点为集合{(x₁,y₁),(x₂,y₂),…,(x_n,y_n)}，各线性拟合因子由以下公式计算 得出：

将轨迹中从第m到第n轨迹点的经纬度转化为直角坐标，对这n-m+1个点 进行最小二乘法得到拟合线段：

上式为线段斜率计算公式，并用下式做角度转换。

langle_i,m,n＝arctan(a),0≤langle_i,m,n≤π

若：a≤0，则langle_i,m,n＝π-langle_i,m,n

以下公式为确定了起止点后，lf_i，k中各元素的计算方法。

起止点飞行时间：

起止点高度差：

起止点距离：

起止点路径的角度：

lf_i,k轨迹序号的合理化约束：

n_e-n_s≥1。

所述步骤3中，所述管制指令包括雷达引导指令、高度指令、速度指令；

图4是本发明实施例根据航迹拟合结果推算管制雷达引导服务位置及次数 示意图，结合图4，所述雷达引导指令的推算方法包括以下步骤：

步骤(3.1.1)：标准飞行程序转弯点拟合，包括以下步骤：

步骤(3.1.1.1)：取F中航班f_i的轨迹线段集合LF_i，λ为轨迹线段序号，令 λ＝0；

步骤(3.1.1.2)：若λ≤num(LF_i),设lf_i,λ的起点为n_s，若λ＞num(LF_i)，设lf_i,λ的 终点为n_s，遍历S中的各个元素，求p_i,ns与S中的各个元素的距离，若转至步骤(3.1.1.3)，否则转至步骤(3.1.1.4)；

步骤(3.1.1.3)：线段起点与标准飞行程序转弯点距离小于门限值的一半， 用标准飞行程序转弯点修正历史飞行轨迹拐点，若令并 更新s_j的经纬度λ＝λ+1，转至步骤(3.1.1.2)；

步骤(3.1.1.4)：创建一个标准飞行程序中路径点迹集合S，令s.long＝p_{i, ns}.long； s.lat＝p_i,ns.lat；F′＝{f_i}，令s∈S。λ＝λ+1，若λ＞num(LF_i)，转至步骤(3.1.1.5)， 否则转至步骤(3.1.1.2)；

步骤(3.1.1.5)：依次遍历轨迹线段集合LF中所有元素，直至结束；

步骤(3.1.1.6)：依次遍历S中元素，计算标准程序点s_j的利用率δ_j，若 δ_j＞μ，保留该关键点，最终得到标准飞行程序转弯点，否则去除该点，令 S＝S-{s_j}，重复本步骤；

步骤(3.1.2)：将标准飞行程序转弯点与历史航空器飞行轨迹比较，将历史 航空器飞行轨迹拐点集合中与标准飞行程序转弯点集合相交的点去除后重新更 新航空器轨迹拐点集合，集合中的元素为航空器飞行过程中机动飞行产生的航 向变化点，机动拐点个数即为管制员为该航空器发布雷达引导指令的个数。

上述雷达引导指令的推算思想是：由于飞行轨迹的拐点处航空器发生了航 向的改变，所以每个航班飞行拐点个数的可以视为其飞行过程中至少改变飞行 方向的次数。图3所示为进近航空器历史飞行轨迹线性拟合图,分析知经过拐点 BHS点的航空器都有航向的改变,但BHS点是标准飞行程序的转弯点,航空器通 过标准飞行程序自主改变航向,故此时不能认为是管制员发布雷达引导指令指挥 航空器改变航向。综上分析拟合过后的飞行轨迹拐点既包括了标准进离场程序 中的转弯点又包括了管制员雷达引导的转弯点。因此雷达引导拐点判断的方法 是首先进行标准程序转弯点进行拟合，将与历史飞行轨迹对比后去除标准程序 转弯点后剩下的飞行拐点用于雷达引导指令识别。

标准飞行程序转弯点拟合算法首先要遍历历史飞行轨迹中各轨迹线段，依 次计算航迹中拐点p_i,j的距离差，对不超过距离门限的拐点进行数据拟合。再依 次遍历集合S中的标准飞行程序转弯点s_j，若s_j与航迹拐点p_i,j距离小于距离门 限值，则认为p_i,j是该雷达引导点。若s_j与航迹拐点p_i,j距离大于距离门限值则认 为p_i,j是非雷达引导点，将拐点p_i,j加入集合S，作为新的标准飞行程序转弯点。 在遍历完所有历史飞行轨迹中各轨迹线段和集合S中的元素后，再计算各标准 飞行程序转弯点s_j的利用率，不满足利用率门限的点舍去，最后生成满足一般降 落轨迹特征的标准飞行程序转弯点集合。标准飞行程序转弯点拟合算法核心公 式及约束条件如下：

设F′和F″为历史航班数据中实际和计划经过标准飞行转弯点的航班集合， num(s_j.F′)为历史数据中经过标准飞行转弯点的航班个数，num(s_j.F″)为飞行计划 中经过s_j的航班个数。μ为利用率门限，其他变量定义见表1。

若s_j与p_i,j合并。s_j的经纬度修正公式：

各标准飞行程序转弯点s_j的利用率δ_j为:

图5是本发明实施例根据航迹拟合结果推算管制指令高度变更次数及位置 示意图，结合图5，所述高度指令的推算方法包括以下步骤：

步骤(3.2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点 并按照p_i(t)升序排序；

步骤(3.2.2)：从p_i的第一个元素p_i,0开始，初始化高度剖面第一条轨迹线段 hf_i,0及轨迹线段起止点的高度值p_i,k(h)和经度坐标值p_i,k(long)；

步骤(3.2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号，令k＝1， λ＝0；

步骤(3.2.4)：设hf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(3.2.5)：求从k至k+1的轨迹点连线的高度差hf_i,k(Δh)，若hf_i,k(Δh)≥60轨迹出现拐点，若λ＝0，转至步骤(3.2.7)，若λ>0，转至步骤(3.2.6)，若 hf_i,k(Δh)＜60，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角 转至步骤(3.2.9)；

步骤(3.2.6)：判断当前路径与上一条路径的夹角，若 |hf_i,λ-1(Δh)-hf_i,λ(Δh)|＞30，转至步骤(3.2.7)；若|hf_i,λ-1(Δh)-hf_i,λ(Δh)|≤30，转至 步骤(3.2.8)；

步骤(3.2.7)：当前路径结束，p_i,k为线段hf_i,λ的结束点hf_i,λ(n_e)＝k，令λ＝ λ+1增加一条新的高度剖面轨迹线段，hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线 段，转至步骤(3.2.9)；

步骤(3.2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点 f_i,λ-1(n_e)＝k，令hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转 至步骤(3.2.9)；

步骤(3.2.9)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令hf_i,λ(n_e)＝n,hf_i,λ(n_e)＝k；

步骤(3.2.10)：将hf_i,0,hf_i,1,...,hf_i,λ写入集合HF_i；

步骤(3.2.11)：重复步骤(3.2.2)-(3.2.7)，直至遍历完F中所有航班，将 HF₀,HF₁,...,HF_n写入集合HF；

步骤(3.2.11)：集合HF_i中元素个数即为航空器在空域内飞行时高度指令个 数。

上述高度指令的推算思想是：首先对高度剖面线进行拟合：

从雷达监视系统中提取相邻的两航班轨迹前后两点的高度和经度坐标数据, 并将其转化为直角坐标。

定义HF＝{HF₁,HF₂,...,HF_m}为F中各航班在高度剖面内进近的轨迹线段集 合。

HF_i＝{hf_i,1,hf_i,2,...,hf_i,m}为高度剖面内航班f_i的进近轨迹线段。

hf_i,m＝{n_s,n_e,Δt,Δh,Δd,a}为高度剖面内每个轨迹线段hf_i,m均包含6个元素，n_s,n_e分别为该线段起点终点在p_i中的序号，Δt为该路径的飞行时间，Δh为该路 径的高度差，Δd为起止点距离，a为该路径的角度。

确定了起止点后，计算高度剖面内hf_i，k中起止点高度差。

根据规定，中国二次航管雷达其C模式下，管制员为航空器之间飞行分配 的高度层间隔为300m，并且，只要航空器的高度变化值在60m以内，就认为该 航空器在此高度层上飞行，故进近飞行高度变化状态解析约束条件为：

若认为航空器处于平飞状态。

若认为航空器处于上升状态。

若认为航空器处于下降状态。

通过高度差判断，可将高度剖面中的轨迹点分段，高度剖面线由轨迹点拟 合产生。采用最小二乘法将高度剖面点迹数据拟合高度剖面线。

设高度剖面线拟合线段直角坐标系的直线方程为L′＝a′x+b′y+c′，即 L′＝{a′,b′,c′}。

定义拟合点为集合{(x₁,y₁),(x₂,y₂),...,(x_n,y_n)}，各线性拟合因子由公式(3.5)计算得出。

将高度剖面轨迹中从第m到第n轨迹点的高度和经度转化为直角坐标，对 这n-m+1个点进行最小二乘法得到拟合线段：

得到

hangle_i,m,n＝arctan(a′),0≤hangle_i,m,n≤π

若：a′≤0，则hangle_i,m,n＝π-hangle_i,m,n。

图6是本发明实施例根据航迹拟合结果推算管制指令速度变更次数及位置 示意图，结合图6，所述速度指令的推算方法包括以下步骤：

步骤(3.3.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点 并按照p_i(t)升序排序；

步骤(3.3.2)：从p_i的第一个元素p_i,0开始，初始化速度剖面第一条轨迹线段 sf_i,0及轨迹线段起止点的高度值p_i,k(s)和经度坐标值p_i,k(long)；

步骤(3.3.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ ＝0；

步骤(3.3.4)：设sf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步 骤；

步骤(3.3.5)：求从k至k+1的轨迹点连线的斜率sangle_i,j，若sangle_i,j≠0轨迹 出现拐点，若λ＝0，转至步骤(3.3.7)，若λ>0，转至步骤(3.3.6)，若sangle_i,j＝0， 轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角转至步骤 (3.3.9)；

步骤(3.3.6)：判断当前路径与上一条路径的斜率，若 |sf_i,λ-1(slope)-sf_i,λ(slope)|＝0，转至步骤(3.3.7)；若|sf_i,λ-1(slope)-sf_i,λ(slope)|≠0， 转至步骤(3.3.8)；

步骤(3.3.7)：当前路径结束，p_i,k为线段sf_i,λ的结束点sf_i,λ(n_e)＝k，令λ＝ λ+1增加一条新的速度剖面轨迹线段，sf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线 段，转至步骤(3.3.8)；

步骤(3.3.8)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完， 记录最后一条路径，令sf_i,λ(n_e)＝n，sf_i,λ(n_e)＝k，转至步骤(3.3.9)；

步骤(3.3.9)：将sf_i,0,sf_i,1,...,sf_i,λ写入集合SF_i；

步骤(3.3.10)：重复步骤(3.3.2)-(3.3.7)，直至遍历完F中所有航班，将 SF₀,SF₁,...,SF_n写入集合SF；

步骤(3.3.11)：SF_i集合中元素个数即为航空器在空域内飞行时速度指令个 数。

上述速度指令的推算思想是：速度剖面线拟合算法：

从雷达监视系统中提取相邻的两航班轨迹前后两点的速度和经度坐标数据, 并将其转化为直角坐标。

定义SF＝{SF₁,SF₂,...,SF_m}为F中各航班在速度剖面内进近的轨迹线段集合。

SF_i＝{sf_i,1,sf_i,2,…,sf_i,m}为速度剖面内航班f_i的进近轨迹线段。

sf_i,m＝{n_s,n_e,Δt,Δs,Δd,slope}为速度剖面内每个轨迹线段sf_i,m均包含的6个元 素，n_s,n_e分别为该线段起点终点在p_i中的序号，Δt为该路径的飞行时间，Δs为 该路径的速度差，Δd为起止点距离，slope为该路径的斜率。

确定了起止点后，在速度剖面内计算得到两轨迹点p_i,j和p_i,j-1组成线段的斜率 用slope_i,j来表示。

由数学知识，当slope_i,j＝0时轨迹线段平行于坐标轴，航空器此时可以认为 是等速运行，也即管制员未进行调速，当slope_i,j＞0时航空器处于增速状态，即 管制员给了航空器一个增速指令，当slope_i,j＜0时航空器处于减速状态，即管制 员给了航空器一个减速指令。当在一个轨迹线段内slope_i,j值不断的处于正负交替 的变化时，即此时航空器在不断的进行机动飞行，此时拟合的速度曲线会出现 间断性的振荡。

通过斜角判断，可将速度剖面中的轨迹点分段，速度剖面线由轨迹点拟合 产生。采用最小二乘法将速度剖面的点迹数据拟合成速度剖面线。设速度剖面 线拟合线段直角坐标系的直线方程为L″＝a″x+b″y+c″，即 L″＝{a″,b″,c″}。

定义拟合点为集合{(x₁,y₁),(x₂,y₂),…,(x_n,y_n)}。将速度剖面轨迹中从 第m到第n轨迹点的速度和经度转化为直角坐标，对这n-m+1个点进行最小二 乘法得到拟合线段：

得到

sangle_i,m,n＝arctan(a″),0≤sangle_i,m,n≤π

若a″≤0，则sangle_i,m,n＝π-sangle_i,m,n。

所述步骤3中，所述管制指令还包括但不限于以下其他指令或飞行情报： 雷达识别、进离场程序通报、气象条件通报、进近许可、脱波，这些指令在管 制过程中仅为航空器提供一次指令，故推算一架航空器以上类型的指令各发布 一次。其他协调信息包括：和相邻管制单位的进行接收协调和移交协调，通常 为一个航空器进入、离开本空域时进行协调，因此可推算一架航空器以上协调 过程各进行一次。

所述步骤4中，按以下公式对中交通管制员的工作负荷进行评估：

其中，为单位时间t内的管制员工作负荷量，单位为秒，σ_i为单位时 间t内管制服务i的次数，t_i表示管制员为提供管制服务i的工作时间。

中交通管制员的工作负荷评估思想是：按照管制指令解析的方法步骤对历 史的航空器飞行轨迹做位置、高度和速度等解析后识别得到一系列的管制指令， 根据表2提供的部分常用的管制用语的工作负荷统计数据表将雷达管制过程中 的管制服务进行分类，以此来推算航空器进近飞行状态下的单位时间优化的和 实际的管制工作负荷。将已经统计得到的进近管制员的单位时间实际管制工作 负荷和优化的管制工作负荷进行数学运算，进而可以精确的量化雷达管制效率。

下表2示出了常用进近雷达管制用语的工作负荷统计：

表2常用进近雷达管制用语的工作负荷统计表

上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制， 只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视 为落入本发明专利的权利保护范围内。

Claims (6)

1.一种基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：包括以下步骤：

步骤1：历史飞行轨迹的获取与处理：选择关注的管制空域及时间范围，获取该时间段内该空域的运行数据，获取航空器的历史飞行轨迹；将历史飞行轨迹中关键字段进行对比，按照时间顺序产生不同航空器的历史飞行轨迹，形成轨迹集合；按照各历史飞行轨迹的运行方向和高度、速度变化，判断该航班为起飞或者降落航班；同时对每个航班在空域内的飞行距离、飞行时间进行统计；

步骤2：遍历轨迹集合中各个航空器的历史飞行轨迹，对轨迹进行线性拟合，分析轨迹中的速度、高度、航向随时间变化的趋势，得到高度、速度、航向突然产生变化的时刻及位置，同时对航空器首次出现时刻位置、最后出现时刻位置、进港航空器到达起始进近定位点并到达起始进近进近高度的时刻进行记录；

步骤3：利用步骤2的记录数据，推算可能产生的管制指令、协调信息和管制监控航空器架次；

步骤4：根据步骤3得到的推算结果计算总的工作时间即工作负荷，完成空中交通管制员工作负荷评估。

2.根据权利要求1所述的基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：所述步骤1中，所述运行数据是来自空管自动化设备或空管记录仪中的监视数据，包括但不限于二次航管雷达数据、ADS-B数据即广播式自动相关监视数据、融合后的系统航迹数据；所述形成轨迹集合的方法是：按行读取运行数据的每一条记录，查询容器中是否存在以该呼号为唯一标识的航空器对象，如果没有则新建航空器对象并加入容器，否则读取该条记录中与航空器起降相关的时间、位置、速度、使用跑道信息到航迹点对象中，并将航迹点数据加入航空器对象的航迹链表，最后在容器中得到以时间排序并且以呼号为标识的航空器对象序列，形成航空器航迹集合；所述判断该航班的起降类型包括以下步骤：(1)同时满足下述三个条件被判断为离港航班即起飞航班：航空器第一个航迹点在机场跑道中心点2km范围内，且高度低于1000m；航空器最后一个航迹点高度、速度大于第一个航迹点高度、速度；航空器航迹点个数大于50个；(2)同时满足下述三个条件被判断为进港航班即降落航班：航空器最后一个航迹点在机场跑道中心点2km范围内，且高度低于1000m；航空器最后一个航迹点高度、速度小于第一个航迹点高度、速度；航空器航迹点个数大于50个；(3)若航空器航迹不满足上述两种条件，该航迹太短，或者数据异常无法判断，则从航迹集合中去除该航空器轨迹。

3.根据权利要求1所述的基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：所述步骤2中，所述对轨迹进行线性拟合包括以下步骤：

步骤(2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点按照p_i(t)升序排序；

步骤(2.2)：从p_i的第一个p_i,0元素开始，初始化第一条轨迹线段lf_i,0，令lf_i,0(n_s)＝0，求两轨迹点连线的夹角angle_i,j，令lf_i,0(a)＝angle_i,0；

步骤(2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ＝0；

步骤(2.4)：设lf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步骤；

步骤(2.5)：求从k至k+1的轨迹点连线的夹角angle_i,k，若|angle_i,k-lf_i,λ(a)|＞θ轨迹出现拐点，若λ＝0，转至步骤(2.7)，若λ>0，转至步骤(2.6)，若|angle_i,k-lf_i,λ(a)|≤θ，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角转至步骤(2.9)；

步骤(2.6)：判断当前路径与上一条路径的夹角，若转至步骤(2.7)；若转至步骤(2.8)；

步骤(2.7)：当前路径结束，p_i,k为线段lf_i,λ的结束点lf_i,λ(n_e)＝k，求lf_i,λ的其他元素值，令λ＝λ+1增加一条新的轨迹线段，lf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(2.9)；

步骤(2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点f_i,λ-1(n_e)＝k，更新lf_i,λ-1的其他元素值，令lf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(2.9)；

步骤(2.9)：令k＝k+1，重复步骤(2.4)，直至k＝n，p_i中点迹已经遍历完，记录最后一条路径，令lf_i,λ(n_e)＝n，lf_i,λ(n_s)＝n，求lf_i,λ的其他元素值；

步骤(2.10)：将lf_i,0,lf_i,1,...,lf_i,λ写入集合LF_i；

步骤(2.11)：重复步骤(2.2)-(2.7)，直至遍历完F中所有航班，将LF₀,LF₁,...,LF_n写入集合LF。

4.根据权利要求1所述的基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：所述步骤3中，所述管制指令包括雷达引导指令、高度指令、速度指令；

所述雷达引导指令的推算方法包括以下步骤：

步骤(3.1.1)：标准飞行程序转弯点拟合，包括以下步骤：

步骤(3.1.1.1)：取F中航班f_i的轨迹线段集合LF_i，λ为轨迹线段序号，令λ＝0；

步骤(3.1.1.2)：若λ≤num(LF_i),设lf_i,λ的起点为n_s，若λ＞num(LF_i)，设lf_i,λ的终点为n_s，遍历S中的各个元素，求p_i,ns与S中的各个元素的距离，若转至步骤(3.1.1.3)，否则转至步骤(3.1.1.4)；

步骤(3.1.1.3)：线段起点与标准飞行程序转弯点距离小于门限值的一半，用标准飞行程序转弯点修正历史飞行轨迹拐点，若令并更新s_j的经纬度λ＝λ+1，转至步骤(3.1.1.2)；

步骤(3.1.1.4)：创建一个标准飞行程序中路径点迹集合S，令s.long＝p_i,ns.long；s.lat＝p_i,ns.lat；F′＝{f_i}，令s∈S。λ＝λ+1，若λ＞num(LF_i)，转至步骤(3.1.1.5)，否则转至步骤(3.1.1.2)；

步骤(3.1.1.5)：依次遍历轨迹线段集合LF中所有元素，直至结束；

步骤(3.1.1.6)：依次遍历S中元素，计算标准程序点s_j的利用率δ_j，若δ_j＞μ，μ为利用率门限，保留该关键点，最终得到标准飞行程序转弯点，否则去除该点，令S＝S-{s_j}，重复本步骤；

步骤(3.1.2)：将标准飞行程序转弯点与历史航空器飞行轨迹比较，将历史航空器飞行轨迹拐点集合中与标准飞行程序转弯点集合相交的点去除后重新更新航空器轨迹拐点集合，集合中的元素为航空器飞行过程中机动飞行产生的航向变化点，机动拐点个数即为管制员为该航空器发布雷达引导指令的个数；

所述高度指令的推算方法包括以下步骤：

步骤(3.2.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点并按照p_i(t)升序排序；

步骤(3.2.2)：从p_i的第一个元素p_i,0开始，初始化高度剖面第一条轨迹线段hf_i,0及轨迹线段起止点的高度值p_i,k(h)和经度坐标值p_i,k(long)；

步骤(3.2.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号，令k＝1，λ＝0；

步骤(3.2.4)：设hf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步骤；

步骤(3.2.5)：求从k至k+1的轨迹点连线的高度差hf_i,k(Δh)，若hf_i,k(Δh)≥60轨迹出现拐点，若λ＝0，转至步骤(3.2.7)，若λ>0，转至步骤(3.2.6)，若hf_i,k(Δh)＜60，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角转至步骤(3.2.9)；

步骤(3.2.6)：判断当前路径与上一条路径的夹角，若|hf_i,λ-1(Δh)-hf_i,λ(Δh)|＞30，转至步骤(3.2.7)；若|hf_i,λ-1(Δh)-hf_i,λ(Δh)|≤30，转至步骤(3.2.8)；

步骤(3.2.7)：当前路径结束，p_i,k为线段hf_i,λ的结束点hf_i,λ(n_e)＝k，令λ＝λ+1增加一条新的高度剖面轨迹线段，hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(3.2.9)；

步骤(3.2.8)：当前路径与上一条路径可以合并，修改上一条路径的结束点f_i,λ-1(n_e)＝k，令hf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(3.2.9)；

步骤(3.2.9)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完，记录最后一条路径，令hf_i,λ(n_e)＝n,hf_i,λ(n_e)＝k；

步骤(3.2.10)：将hf_i,0,hf_i,1,…,hf_i,λ写入集合HF_i；

步骤(3.2.11)：重复步骤(3.2.2)-(3.2.7)，直至遍历完F中所有航班，将HF₀,HF₁,...,HF_n写入集合HF；

步骤(3.2.11)：集合HF_i中元素个数即为航空器在空域内飞行时高度指令个数。

所述速度指令的推算方法包括以下步骤：

步骤(3.3.1)：遍历所有进近航班中的轨迹集合p_i，将轨迹集合中各轨迹点并按照p_i(t)升序排序；

步骤(3.3.2)：从p_i的第一个元素p_i,0开始，初始化速度剖面第一条轨迹线段sf_i,0及轨迹线段起止点的高度值p_i,k(s)和经度坐标值p_i,k(long)；

步骤(3.3.3)：设k为轨迹点遍历次数序号，λ为轨迹线段序号,令k＝1，λ＝0；

步骤(3.3.4)：设sf_i,λ的起点为n_s，求点与p_i,k的距离若转至下一步骤，若令k＝k+1，重复本步骤；

步骤(3.3.5)：求从k至k+1的轨迹点连线的斜率sangle_i,j，若sangle_i,j≠0轨迹出现拐点，若λ＝0，转至步骤(3.3.7)，若λ>0，转至步骤(3.3.6)，若sangle_i,j＝0，轨迹继续延伸，将p_i,k加入轨迹中，求线段的斜角转至步骤(3.3.9)；

步骤(3.3.6)：判断当前路径与上一条路径的斜率，若|sf_i,λ-1(slope)-sf_i,λ(slope)|＝0，转至步骤(3.3.7)；若|sf_i,λ-1(slope)-sf_i,λ(slope)|≠0，转至步骤(3.3.8)；

步骤(3.3.7)：当前路径结束，p_i,k为线段sf_i,λ的结束点sf_i,λ(n_e)＝k，令λ＝λ+1增加一条新的速度剖面轨迹线段，sf_i,λ(n_s)＝k，以p_i,k为开始点记录新的线段，转至步骤(3.3.8)；

步骤(3.3.8)：令k＝k+1，重复步骤4，直至k＝n，p_i中点迹已经遍历完，记录最后一条路径，令sf_i,λ(n_e)＝n，sf_i,λ(n_e)＝k，转至步骤(3.3.9)；

步骤(3.3.9)：将sf_i,0,sf_i,1,...,sf_i,λ写入集合SF_i；

步骤(3.3.10)：重复步骤(3.3.2)-(3.3.7)，直至遍历完F中所有航班，将SF₀,SF₁,...,SF_n写入集合SF；

步骤(3.3.11)：SF_i集合中元素个数即为航空器在空域内飞行时速度指令个数。

5.根据权利要求4所述的基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：所述步骤3中，所述管制指令还包括但不限于以下其他指令或飞行情报：雷达识别、进离场程序通报、气象条件通报、进近许可、脱波，这些指令在管制过程中仅为航空器提供一次指令，故推算一架航空器以上类型的指令各发布一次。

6.根据权利要求1所述的基于历史轨迹数据的空中交通管制员工作负荷评估方法，其特征在于：所述步骤4中，按以下公式对中交通管制员的工作负荷进行评估：

其中，为单位时间t内的管制员工作负荷量，单位为秒，σ_i为单位时间t内管制服务i的次数，t_i表示管制员为提供管制服务i的工作时间。