0.   引言

国际民航组织文件Doc 8168-OPS/61根据空域安全间隔专家组的研究成果, 制定了空域规划的方法和原则。国际民航组织文件Doc9426-AN/92在附录中关于ATC扇区容量评估方法, 给出扇区容量的量化方法, 成为国际空中交通管理人员进行空域规划设计的理论依据。世界空管领域的不少学者致力于空域优化问题的研究, 日本电子导航研究所的N Tofukuji教授^[1]根据在交通繁忙的情况下, 管制员工作的完成情况来判断容量极限是否达到, 从而推算出终端区管制扇区的实际最大容量; 法国学者D Delahaye^[2,3]将遗传算法和计算几何相结合, 提出空域的优化方法等。这些研究成果一方面仅限于针对相关国家和地区的空域结构的特点进行分析和研究, 研究的对象是区域导航(RNAV)的空域, 而中国的空域基本上是固定航路、航线空域; 另一方面, 算法的鲁棒性和收敛性不如模拟退火算法强。韩松臣、张明^[4]等人曾经对于空域工作负荷的量化进行了研究, 但在扇区优化方面没有提出较为成熟的随机优化算法。通过分析国际民航组织的有关文件和国际上在空域规划设计方面先进的研究成果, 考虑到空域结构的复杂性、空域交通密度、系统可用性等方面的因素, 本文提出了基于模拟退火算法和网络拓扑关系, 对于实际的管制空域进行扇区的最优划分方法。

1.   空域结构拓扑表达

借助于计算机进行空域的优化扇区工作, 首先必须建立管制空域结构的拓扑数学描述和对应的工作负荷计算模型, 这是利用计算机解决此类问题所必须的数学基础。实际的管制空域包括导航台、机场、位置报告点和连接航路节点的航路组成的网络结构。通过对这些点和线段的描述, 建立拓扑网络数据库, 计算机通过读取数据库数据, 便能自动绘制网络图, 计算空域扇区负荷, 根据算法随机选择航路, 完成扇区优化。

拓扑图的数学形式可通过关联矩阵来描述各线段与航线之间的关系。若设E为线段矩阵, C为航线矩阵, G为m×n阶线段和航线的关联矩阵, 则有

$\begin{array}{l}\mathit{E}=[e{}_{1}e{}_{2}e{}_3\cdots e{}_n]{}^{\text{Τ}}(1)\\ \mathit{C}=[c{}_{1}c{}_{2}c{}_3\cdots c{}_m]{}^{\text{Τ}}(2)\\ \mathit{G}=(g{}_{ij})(i=1,\cdots ,m;j=1,\cdots ,n)(3)\end{array}$

式中: 矩阵G的元素g_ij∈{0, 1}, g_ij=1为计划飞行的路线c_i经过航段e_j, g_ij=0为飞行计划的路线c_i与航段e_j不相关。

线段和航线之间关系为

$\mathit{C}=\mathit{G}\mathit{E}(4)$

若Q设为航线上的流量, 即飞机数, 则有

$\mathit{Q}=[q{}_{c{}_1}q{}_{c{}_2}q{}_{c{}_3}\cdots q{}_{c{}_m}]{}^{\text{Τ}}(5)$

D代表不同线段上的流量矩阵, 各航线流经该线段上的流量之和, 则有

$\begin{array}{l}\mathit{D}=[d{}_{1}d{}_{2}d{}_3\cdots d{}_n]{}^{\text{Τ}}(6)\\ \mathit{Q}=\mathit{G}\mathit{D}(7)\end{array}$

若设F为n×1的工作负荷矩阵, 则有

$\mathit{F}=[f{}_{1}f{}_{2}f{}_3\cdots f{}_n](8)$

F表明在每一航段, 管制人员指挥单架飞机的负荷值。

由式(6)、(8), 可得总工作负荷为

$\begin{array}{l}\mathit{W}=\mathit{F}\mathit{D}=[f{}_{1}f{}_{2}f{}_3\cdots f{}_n]\cdot \\ [d{}_{1}d{}_{2}d{}_3\cdots d{}_n]{}^{\text{Τ}}={\displaystyle \sum f}{}_{i}d{}_j(9)\end{array}$

2.   管制工作负荷的统计方法

根据英国运筹与分析理事会提出的Doratask方法^[4], 将管制员的工作负荷分为看得见部分、看不见部分和恢复时间。看得见部分是指管制员进行常规管制工作、解决冲突的通信和填写进程单等能够被观察员记录和记时的工作; 看不见的工作是指管制员监控雷达屏幕、对照进程单、思考的恢复时间等工作, 这部分工作不便进行记录和记时。此外, Doratask方法认为必须为管制员留有一定的恢复时间, 恢复时间对扇区的安全运行极为重要。对于看不见部分和恢复时间进行统计比较困难, 不过, 这部分的负荷可以在统计的时候根据统计分析, 加入适当的余度表示。根据对厦门管制区的统计, 结合Doratask方法提出平均工作负荷强度小于统计时间的80%, 对厦门管制区的实际观察和统计, 认为看得见部分应当占整个工作负荷的50%左右, 即总工作负荷为看得见工作负荷的2倍左右。

基于以上思想, 根据《中国民用航空无线电通话手册》的标准管制用语, 通过统计雷达、语音记录仪记录的历史数据或对实际管制工作进行记时和统计(表 1为部分民航常用管制用语工作负荷), 确定相应通话种类的经验工作负荷指数的期望值(以时间长短衡量)。第k种工作负荷经验负荷指数为

$\overline{t}{}_k={\displaystyle \sum _{j=1}^{{\rm N}}{\rm T}}{}_{kj}$ N (N > 50) (10)

表  1  部分常用管制用语工作负荷统计

Table  1.  Part workload statistics of common control dictions

管制类型	管制种类	管制用语	工作负荷/s
一般用语	高度指令	在(高度)停止上升/下降	4.5
一般用语	活动通报	没有得到活动报告	3.5
一般用语	气象情报	报告飞行条件	3.5
机场管制	开车程序	开车时间自己掌握	4.0
机场管制	滑行程序	滑行时注意	3.5
机场管制	起飞许可	立即起飞, 否则在跑道外等待	4.7
进近管制	进近指令	保持目视间隔	4.2
进近管制	预计进近时间	没有延误	3.2
区域管制	管制许可	ATC许可(航空器呼号)	6.0
区域管制	紧急下降	紧急下降	3.1

下载: 导出CSV 

| 显示表格

式中: T_kj为第k种管制动作第j次统计的时间值; N为一种工作负荷统计的次数, 实际测量的每种工作负荷应当大于50次。

由于总工作负荷为看得见工作负荷的2倍, 若总的工作负荷为S, 在某时间段包含M种看得见工作负荷, 则有

$S=2{\displaystyle \sum _{k=1}^{{\rm M}}(}{\displaystyle \sum _{j=1}^{{\rm N}}{\rm T}}{}_{kj}/{\rm N})({\rm N}>50)(11)$

3.   扇区划分的数学模型

基于以上管制工作负荷原则, 应当保证空域总的工作负荷最小。这就要求在进行扇区优化的时候, 在满足管制员正常负荷的条件下, 尽量减少扇区的数目, 以减少过多的扇区之间的协调工作负荷。设计的扇区数目满足以下关系

$\begin{array}{l}{\rm N}{}_{\text{s}}=\{\begin{array}{cc}\mathrm{int}({\rm M})+1\hfill & \left({\rm M}\ne \mathrm{int}({\rm M})\right)\hfill \\ {\rm M}\hfill & \left({\rm M}=\mathrm{int}({\rm M})\right)\hfill \end{array}(12)\\ {\rm M}=\frac{W{}_{\text{Τ}}}{S\times 80\%}\end{array}$

式中: N_s为最小扇区数; W_T为空域总工作负荷, 用时间表示。根据Doratask方法, 在设计扇区的时候, 应当取平均工作强度, 实际扇区的工作负荷为管制员工作负荷的80%。当M计算结果不为整数时, 对计算结果取整加1, 得到N_s为最小扇区数; 如果M计算结果为整数, 直接得到N_s为最小扇区数。

在满足扇区数最小化的前提下, 还应当保证使扇区之间的工作负荷较为均衡, 建立如下扇区优化的数学模型, 目标函数为

$J=\min ({\displaystyle \sum _{i=1}^{{\rm N}{}_{\text{s}}}{\displaystyle \sum _{j=1}^{{\rm N}{}_{\text{s}}}}}{\rm Z}{}_i-{\rm Z}{}_j)(13)$

约束条件: 根据飞行计划, 固定航班经历的各个航段, 作为算法的优先搜索航段; i≠j, 根据节点和航段的拓扑关系, 保证随机搜索得到的扇区连续性; 满足式(12)的要求。式(13)中: J为目标值; Z表示扇区实际划分的工作负荷。该式表示划分出的任意两扇区的差值最小化。

4.   扇区最优划分的模拟退火算法

网络拓扑问题是NP问题, 随机优化算法, 特别是模拟退火算法对于解决大型网络问题比较有效。模拟退火算法(Simulated Annealing以下简称SA)是将组合优化问题与统计力学中的热平衡问题相类比, 另辟的求解组合优化问题的新途径。

假设优化问题的一个解i及其目标函数f(i)分别与固体的一个微观状态i及其能量E_i等价, 令随算法进程递减的控制参数t担当固体退火过程中温度T的角色, 则对于t的每个取值, 算法持续进行产生新解-判断-接受/舍弃的迭代过程对应着固体在某一温度下趋于热平衡的过程, 就是执行了一次Metropolis算法。模拟退火算法从某个初始解出发, 经过大量解的变换后, 可以求得给定控制参数值时优化问题的相对最优解; 然后减小t的值, 重复执行Metropolis算法, 就可以在t→0时, 得到优化问题最终的整体最优解。模拟退火算法用Metropolis算法产生优化问题解的序列, 并由与Metropolis准则相对应的转移概率确定是否从当前解i到新解j的转移。以下是概率计算公式

${\rm P}{}_i(i\to j)=\{\begin{array}{cc}1\hfill & (f(i)\le f(j))\hfill \\ \exp \left[\frac{f(i)-f(j)}{t}\right]\hfill & (\text{其}\text{他})\hfill \end{array}(14)$

模拟退火算法依据Metropolis准则接受新解, 因此除接受优化解外, 还在一个限定范围内接受恶化解, 这正是模拟退火算法优于局部搜索算法的本质区别。开始时, t的值较大, 可能接受较差的恶化解, 随着t的减小, 只能接受较好的恶化解, 最后在t趋于零时, 就不再接受任何恶化解了。这就使模拟退火算法可以从局部最优的陷阱中跳出, 最后求出整体最优解。

结合前面介绍的网络拓扑关系及数学模型, 采用模拟退火算法对空域扇区优化, 其优化原则如下。

根据两节点纬度值, 高纬度节点搜索的下一个节点优先选择高纬度节点, 低纬度节点下一个节点优先选择低纬度节点, 然后再按照都优先选择经度较小的下一个节点。如果搜索到的节点的经度相同, 设置随机数, 随机选择下一个节点。如果搜索到网络的末端, 需要返回搜索, 要记录上一个已经历的节点及它的相邻节点域, 下次搜索从上一个节点选择新的节点(当然, 必须排除已搜索过的节点), 将搜索选择的线段存入数组中。每次搜索都进行扇区负荷值的计算。如果新搜索到的线段在另一扇区已经被搜索, 并纳入其中, 将结合目标函数, 重新决定该线段的归属。本次搜索的目标值与上次目标值之差根据式(14), 决定是否进行状态和目标值的迭代, 即根据式(14)确定接受非改善解准则。最后, 当满足终止条件和数学模型中的约束条件, 输出三扇区负荷大小及三扇区包括的航段。终止准则为: y值小于a且T值减小到a。T为退火温度, 每次搜索, T变化为0.95T; y为上一次搜索的目标函数值; a为终止条件, 如果管制扇区的航段数为n, 则a设为总管制工作负荷的n分之一。

5.   实例分析

采用Visual C++6.0对厦门终端区高空航路结构图进行了优化, 优化结果均满足终止条件。表 2是根据5月17日上午10:10到10:15统计的一组航段工作负荷值(包括看得见、看不见部分和恢复时间之和, 工作负荷单位为s)。

表  2  厦门终端区航段工作负荷

Table  2.  Routes workloads of Xiamen control airspace

航段	负荷/s	航段	负荷/s
IKATA-PV	27.1	P176-PV	21.3
IKATA-ZPU	27.4	P163-PV	7.2
ZPU-WF	18.9	P176-DO	0.0
WF-P169	34.1	LJG-DAGMO	10.2
FQG-LJG	9.1	IKATA-GUVIA	31.6
P84-DO	10.1	P93-ZPU	22.6
P48-DO	11.9	WF-P84	10.1
LJG-DO	9.3	P169-FQG	10.3
LJG-RUPOX	7.2	FQG-DO	14.0
DO-MAKIB	31.5	DO-P132	17.8
WYS-P249	22.1	DO-WYS	38.1
P47-DO	11.7

下载: 导出CSV 

| 显示表格

对于这组值进行空域扇区优化, 根据扇区优化结果进行扇区划分。图 1为空域扇区优化图。图中所有的21个节点(如IKATA、ZPU等)均为厦门管制区各个航段的航路节点, 即表 1中各个航段的节点; 点划线为优化后的扇区边界, 将空域划分为南扇区、东扇区和北扇区, 根据优化原则, 该扇区优化是最优划分; 虚线段组成的封闭环为厦门管制区边界线; 实线段组成的网络是厦门管制区航路网, 整个优化过程围绕该网络展开。优化后北扇区的工作负荷为133.1 s, 东扇区工作负荷为133.3 s, n为23, 终止条件a为17.5 s, 南扇区工作负荷为137.2 s, 扇区负荷差为4.1 s, 同时T为6.7, 均满足规定的终止条件17.5 s的要求。图 2扇区优化目标函数呈现收敛曲线, 并且在第21次迭代时, 根据式(14)Metropolis准则接受恶化解, 使得目标函数值继续不断递减, 经过109次迭代得到近似全局最优解。图 3表示南扇区、东扇区和北扇区的工作负荷在优化过程中的变化。通过优化, 发现在终止条件和约束条件下, 扇区之间的工作负荷差值在空域的工作负荷扇区优化达到最优。另外, 对于优化的空域, 建立固定航线所包含的航段的数据库, 作为模拟退火算法迭代时的约束条件。

图  1  厦门管制区空域优化

Figure  1.  Optimum of Xiamen control airspace

下载: 全尺寸图片 幻灯片

图  2  收敛曲线

Figure  2.  Convergence curves

下载: 全尺寸图片 幻灯片

图  3  工作负荷变化曲线

Figure  3.  Change curve of workload

下载: 全尺寸图片 幻灯片

应用结果表明该优化算法具有通用性、有较强的鲁棒性和收敛性, 可以保证得到的解是全局最优解或全局近似最优解。