一种基于离散势能场的危险天气下大量航班改航路径并行化计算方法与流程
本发明属于航班改航路径规划技术，特别是一种基于离散势能场的危险天气下大量航班改航路径并行化计算方法。
背景技术：
危险天气是威胁航空运输安全并导致航班延误的重要原因，当某些空域或航路受危险天气影响而无法正常使用时，对航班实施改航可以有效地解决这一问题。改航就是安排航班临时选择一条不受影响的航线绕过飞行受限区域，航班改航路径规划是改航的一个关键问题。
目前对于航班改航的研究开始于上世纪90年代，经过约20年的研究，国内外学者已经提出了多种相关的算法，主要包括：基于网格的改航路径算法，它是网络改航环境模型中运用路径搜索算法对航路进行规划的方法；基于可视图-权重的改航路径算法，它是将影响飞行的区域划分为不同等级并赋予响应的权重，然后应用可视图确定权重距离最短的航班改航路径；基于已有航路点的a*搜索算法，通过a*算法从航路点集中进行路径搜索。以上这些方法虽可有效进行航班改航路径规划，但大多针对单个航班，计算效率低，未考虑多航班甚至是大量航班时的路径规划问题，同时考虑的限制因素较为单一，不具备限制因素的可扩展性。
技术实现要素：
针对上述问题，本发明提出了一种针对危险天气影响下的大量航班飞行快速生成改航路径的方法，特别是提出来一种基于离散势能场的危险天气下大量航班改航路径并行化计算方法。
本发明的技术方案是：一种基于离散势能场的危险天气下大量航班改航路径并行化计算方法，其操作步骤如下：
步骤1、环境场景栅格化及不适度场构建，即对受天气影响的危险区域，提取危险区域的范围，将危险区域范围划设成多边形区域，并对环境场景进行分割，构成细小单元格并计算危险区域覆盖的单元格数，紧接着根据区域航班密度计算每个单元格的不适度值，当密度高于一定阈值时，不适度值为正无穷，即该单元格无法通行；
步骤2、速度场构建，建立栅格中的航班密度贡献函数，确定每个航班对于单元格的密度贡献值以及每个单元格的密度映射值，紧接着将航班的密度贡献值作为权值，计算单元格内的加权平均速度，并按照上、下、左、右方向分解为向该方向飞行时可达到的最大速度；
步骤3、建立代价函数，即构建与路径长度、耗费时间、不适度三个指标相关的函数，以计算单个航班选择飞向目标点路径时的代价，紧接着对代价函数进行离散化，根据步骤1中的栅格化处理，构建航班从起始单元格进入上、下、左、右的相邻单元格的单位代价函数；
步骤4、建立离散势能场计算函数，即根据步骤3所得的当前位置到目标位置的路径的代价值，建立平面空间内从起始位置s到目标位置g的最优路径的搜索优化问题，通过计算区域内任一单元格的势能值，利用单侧差分方程进行离散化，以负梯度方向作为改航最优路径；
步骤5、并行化势能场求解，即采用分区运算规则，将环境区域的所有单元格划分为预定义大小的区块tile，并行求解每个tile中的单元格的势能值，从而求解单元格中的势能场值，利用gpu的计算效能，快速生成势能场。
其中，步骤1中环境场景栅格化及不适度场构建步骤如下：
步骤1.1、将环境场景分别在x方向和y方向划分单元格，将飞行受限区转化为多边形并映射到单元格中；
步骤1.2、设定单元格内的航班越多该单元格的不适度值越高，当密度高于一定阈值时，不适度值为正无穷；定义如果单元格的不适度值为正无穷，则该单元格代表完全无法通行，飞行受限区所属单元格的不适度值初始化为正无穷；单元格的不适度值的计算由如下公式给出：
其中，步骤2中速度场构建步骤如下：
步骤2.1、建立航班对于单元格的密度贡献，设航班与单元格c的中心点在x轴方向的距离为δx，在y轴方向上的距离为δy；该航班对单元格a、b、c、d的密度贡献可用公式计算：
步骤2.2、计算单元格的密度贡献值，根据步骤2.1将对当前单元格具有密度贡献的航班的密度值相加，即可求得当前单元格的密度映射值，公式如下所示：
ρi,j＝∑a∈aρa，
步骤2.3、计算单元格的加权平均速度，将航班对当前单元格的密度贡献作为权值，计算求得第(i，j)个单元格内的加权平均速度，如以下公式所述：
步骤2.4、计算航班往特定方向飞行可达到的最大速度，其计算与航班密度相关，计算公式为：
步骤3中：航班改航路径的代价的积分运算可以用从单元格进入上、下、左、右的相邻单元格的单位代价值之和进行近似，构建代价函数计算方式如下所述：
步骤4中离散势能场计算函数构建步骤如下：
步骤4.1、离散势能场计算函数构建，根据upwind离散方法进行离散化近似求解，势能场值计算公式如下：
步骤4.2、计算势能场的梯度，根据梯度的定义，梯度方向指向标量场增长最快的方向，采取单侧差分方程进行离散化，针对x方向上的梯度分量，计算的方法为：
步骤5中并行化势能场求解分为算法定义、算法初始化、循环计算三个步骤，具体实现过程如下：
步骤5.1、算法定义，即设定t表示将环境区域进行等尺寸划分的的单个tile，设定运算列表l用于存储需要更新的tile；
步骤5.2、算法初始化，即将所有单元格均匀划分为相同大小的区块tile，将所有包含目标区域单元格的tile移入运算列表l；
步骤5.3、循环计算，即通过多次迭代将各单元格的势能场值计算至收敛以完成更新。
步骤5.3中循环计算具体实现过程如下：
步骤5.3.1、对运算列表l中的每一个区块t进行更新运算：循环n次计算t中的所有单元格的势能值；每计算完一次后比较当前的结果与前一次的结果，如果差值小于阈值，则该单元格标记为收敛；遍历t中的所有单元格，如果所有单元格都被标记为收敛，则t标记为收敛，否则继续进行上一步中的循环；
步骤5.3.2、检查区块t的四向邻tile：对上一步中的已收敛的所有区块t，对其所有相邻区块进行一次势能值计算，检查t的邻区块中是否有单元格势能值发生变化；将上步中发生变化的区块t重新加入运算列表l，势能值没有发生变化的t从l中删除；
步骤5.3.3、将l包含所有未收敛的tile，重复步骤a和b直至l为空。
本发明的有益效果：(1)利用势能场模型概念建立了适用于改航场景下的代价函数，考虑了燃油经济性、飞行受限区等约束，建立了可扩展的代价函数，具备限制条件的扩展性；(2)势能场模型为宏观计算模型，通过建立势能场可为大量航班的改航活动计算改航最优路径；(3)设计了并行化的计算算法，可以充分利用gpu的运算特点与效率，使计算性能提高5倍左右。
附图说明
图1是本发明的流程图；
图2是本发明中离散单元格结构图；
图3是本发明中单元格密度映射规则；
图4是本发明中并行化势能场求解算法流程图；
图5是本发明中tile网格示意图。
具体实施方式
下面结合附图和具体实例对本发明的具体实施方式作进一步说明，但本发明的实施和保护不限于此。
本发明中，一种基于离散势能场的危险天气下大量航班改航路径并行化计算方法，其操作步骤如下：
步骤1、环境场景栅格化及不适度场构建，即对受天气影响的危险区域，提取危险区域的范围，将危险区域范围划设成多边形区域，并对环境场景进行分割，构成细小单元格并计算危险区域覆盖的单元格数，紧接着根据区域航班密度计算每个单元格的不适度值，当密度高于一定阈值时，不适度值为正无穷，即该单元格无法通行；
步骤2、速度场构建，建立栅格中的航班密度贡献函数，确定每个航班对于单元格的密度贡献值以及每个单元格的密度映射值，紧接着将航班的密度贡献值作为权值，计算单元格内的加权平均速度，并按照上、下、左、右方向分解为向该方向飞行时可达到的最大速度；
步骤3、建立代价函数，即构建与路径长度、耗费时间、不适度三个指标相关的函数，以计算单个航班选择飞向目标点路径时的代价，紧接着对代价函数进行离散化，根据步骤1中的栅格化处理，构建航班从起始单元格进入上、下、左、右的相邻单元格的单位代价函数；
步骤4、建立离散势能场计算函数，即根据步骤3所得的当前位置到目标位置的路径的代价值，建立平面空间内从起始位置s到目标位置g的最优路径的搜索优化问题，通过计算区域内任一单元格的势能值，利用单侧差分方程进行离散化，以负梯度方向作为改航最优路径；
步骤5、并行化势能场求解，即采用分区运算规则，将环境区域的所有单元格划分为预定义大小的区块tile，并行求解每个tile中的单元格的势能值，从而求解单元格中的势能场值，利用gpu的计算效能，快速生成势能场。
其中，步骤1中环境场景栅格化及不适度场构建步骤如下：
步骤1.1、将环境场景分别在x方向和y方向划分单元格，将飞行受限区转化为多边形并映射到单元格中；
步骤1.2、设定单元格内的航班越多该单元格的不适度值越高，当密度高于一定阈值时，不适度值为正无穷；定义如果单元格的不适度值为正无穷，则该单元格代表完全无法通行，飞行受限区所属单元格的不适度值初始化为正无穷；单元格的不适度值的计算由如下公式给出：
其中，步骤2中速度场构建步骤如下：
步骤2.1、建立航班对于单元格的密度贡献，设航班与单元格c的中心点在x轴方向的距离为δx，在y轴方向上的距离为δy；该航班对单元格a、b、c、d的密度贡献可用公式计算：
步骤2.2、计算单元格的密度贡献值，根据步骤2.1将对当前单元格具有密度贡献的航班的密度值相加，即可求得当前单元格的密度映射值，公式如下所示：
ρi,j＝∑a∈aρa，
步骤2.3、计算单元格的加权平均速度，将航班对当前单元格的密度贡献作为权值，计算求得第(i，j)个单元格内的加权平均速度，如以下公式所述：
步骤2.4、计算航班往特定方向飞行可达到的最大速度，其计算与航班密度相关，计算公式为：
步骤3中：航班改航路径的代价的积分运算可以用从单元格进入上、下、左、右的相邻单元格的单位代价值之和进行近似，构建代价函数计算方式如下所述：
步骤4中离散势能场计算函数构建步骤如下：
步骤4.1、离散势能场计算函数构建，根据upwind离散方法进行离散化近似求解，势能场值计算公式如下：
步骤4.2、计算势能场的梯度，根据梯度的定义，梯度方向指向标量场增长最快的方向，采取单侧差分方程进行离散化，针对x方向上的梯度分量，计算的方法为：
步骤5中并行化势能场求解分为算法定义、算法初始化、循环计算三个步骤，具体实现过程如下：
步骤5.1、算法定义，即设定t表示将环境区域进行等尺寸划分的的单个tile，设定运算列表l用于存储需要更新的tile；
步骤5.2、算法初始化，即将所有单元格均匀划分为相同大小的区块tile，将所有包含目标区域单元格的tile移入运算列表l；
步骤5.3、循环计算，即通过多次迭代将各单元格的势能场值计算至收敛以完成更新。
步骤5.3中循环计算具体实现过程如下：
步骤5.3.1、对运算列表l中的每一个区块t进行更新运算：循环n次计算t中的所有单元格的势能值；每计算完一次后比较当前的结果与前一次的结果，如果差值小于阈值，则该单元格标记为收敛；遍历t中的所有单元格，如果所有单元格都被标记为收敛，则t标记为收敛，否则继续进行上一步中的循环；
步骤5.3.2、检查区块t的四向邻tile：对上一步中的已收敛的所有区块t，对其所有相邻区块进行一次势能值计算，检查t的邻区块中是否有单元格势能值发生变化；将上步中发生变化的区块t重新加入运算列表l，势能值没有发生变化的t从l中删除；
步骤5.3.3、将l包含所有未收敛的tile，重复步骤a和b直至l为空。
结合图1，本发明基于离散势能场的并行化改航路径计算方法，步骤如下：
第一步，对环境场景进行场景栅格化，并根据飞行受限区建立不适度场；根据环境场景的范围分别在x方向和y方向划分n个单元格，对于每一个单元格可以用坐标(i，j)来标识；每个单元格存储了密度、速度、势能等数据，如图2所示，其中不适度值代表某组具有相同目标的航班对于进入该单元格飞行的选择倾向度，本发明假设单元格的不适度与单元格内的航班密度相关，即单元格内的航班越多该单元格的不适度值越高；其次将所有待计算航班按照改航目标单元格进行分组；其具体实施步骤如下：
(1)、环境场景分割：将环境场景分别在x方向和y方向划分单元格，将飞行受限区转化为多边形并映射到单元格中，显然，单元格的面积越小，计算结果越精细；单元格面积越大，计算速度越快，为了保证计算的实时性，需要选择适当大小的单元格；
(2)、不适度场构建：为简化计算，本发明假设不适度仅仅与航班密度相关，即单元格内的航班越多该单元格的不适度值越高，当密度高于一定阈值时，不适度值为正无穷；如果单元格的不适度值为正无穷，则该单元格代表完全无法通行，飞行受限区所属单元格的不适度值初始化为正无穷，不适度的计算由如下公式给出：
式中：ci,j——单元格(i，j)的不适度值；ni,j——单元格(i，j)内的航班数；ε——单元格的航班数密度阈值。
第二步，如图3所示，建立航班密度与速度的关联关系，即建立一种密度映射规则，并将航班对当前单元格的密度贡献作为权值，计算求得第(i，j)个单元格内航班的加权平均速度，并将速度分解为表示单元格内的所有航班向上、下、左、右四个方向飞行时受到影响能够达到的最大速度，称为单元格(i，j)的上、下、左、右方向的航班限制速度。其具体步骤如下：
(1)、建立航班对于单元格的密度贡献，每一个航班具有单独的密度场，该密度场在该航班位置处取得最大值，而向其四周方向逐渐降低，设航班与单元格c的中心点在x轴方向的距离为δx，在y轴方向上的距离为δy，如图所示；该航班对单元格a、b、c、d的密度贡献可用公式计算：
式中：ρx——航班对单元格x的密度贡献；δx——航班与其所处单元格的中心点在x轴方向的距离；δy——航班与其所处单元格的中心点在y轴方向的距离；λ——密度的衰减速度，λ决定了航班对其所处单元格的密度贡献不小于
(2)、计算单元格的密度贡献值，通过公式3.3计算得到航班对其可以影响到的单元格的密度贡献，从而将对当前单元格具有密度贡献的航班的密度值相加，即可求得当前单元格的密度映射值，公式如下所示：
ρi,j＝∑a∈aρa，
式中：ρi,j——单元格(i，j)的密度映射值；a——单个航班；a——对单元格(i，j)具有密度贡献的所有航班集合；ρa——航班a对单元格(i，j)的密度贡献值；
(3)、计算单元格的加权平均速度，将航班对当前单元格的密度贡献作为权值，计算求得第(i，j)个单元格内的加权平均速度，如公式所述：
式中：
(4)、当对应飞行方向的相邻单元格内的航班密度较大时，则该航班向该方向的飞行速度上限即为该单元格在该方向上的的航班流限制速度；当对应方向的邻单元格的密度较小时，则表明该航班在该方向的飞行可以不受干扰，其可以保持自己的理想速度飞行；若对应方向的相邻单元格的密度处于一个中间值时，则该航班在该方向上的飞行速度由线性插值法获得。综上可得密度影响速度的公式为：
式中：
第三步，建立代价函数，在栅格化的环境区域内，路径可以看成一系列的有序单元格序列，航班从某个单元格沿着某方向飞行可以进入其上、下、左、右的相邻单元格，因此，单条路径的代价的积分运算可以用从单元格进入上、下、左、右的相邻单元格的单位代价值之和进行近似，即构建代价函数计算方式如下所述：
式中：
第四步，建立离散势能场计算函数，构建环境场景的标量场p，将标量场p的梯度下降曲线做为平面上a点到b点的最优路径，其具体步骤如下：
(1)、势能场的构建函数可用如下公式表示：
式中：g——表示目标区域的单元格集合；
根据upwind离散方法，上述公式可进行离散化进行近似求解，公式如下：
式中：pi,j——单元格(i，j)的势能值；px——x轴方向邻格势能值；py——y轴方向邻格势能值；ex——进入与px相同的x方向邻格单位代价值；ey——进入与py相同的y方向邻格单位代价值。其中ex和ey的值由步骤三计算获得；
(2)、计算势能场的梯度，根据梯度的定义，梯度方向指向标量场增长最快的方向，因而计算单元格(i，j)的梯度，需对比该单元格与其上、下、左、右四个邻格的势能值，本发明采取单侧差分方程进行离散化，针对x方向上的梯度分量，计算的方法为：
y方向上的梯度分量计算方法与之类似。
第五步，采用并行化计算算法对势能场进行数值求解，采用分区运算规则，将环境区域的所有单元格划分为预定义大小的区块(tile)，其中每个tile中的单元格的势能值的求解并行进行，具体步骤如下：
(1)、算法定义：
a.设定t表示将环境区域进行等尺寸划分的的单个tile；
b.设定运算列表l，用于存储需要更新的tile；
(2)、算法初始化：
a.将所有单元格均匀划分为相同大小的区块tile；
b.将所有包含目标区域单元格的tile移入运算列表l；
(3)、循环计算：
a.对运算列表l中的每一个区块t进行更新运算：循环n次计算t中的所有单元格的势能值；每计算完一次后比较当前的结果与前一次的结果，如果差值小于阈值，则该单元格标记为收敛。遍历t中的所有单元格，如果所有单元格都被标记为收敛，则t标记为收敛，否则继续进行上一步中的循环。
b.检查区块t的四向邻tile：对上一步中的已收敛的所有区块t，对其所有相邻区块进行一次势能值计算，检查t的邻区块中是否有单元格势能值发生变化；将上步中发生变化的区块t重新加入运算列表l，势能值没有发生变化的t从l中删除；
c.将l包含所有未收敛的tile，重复步骤a和b直至l为空。
本发明的效果在于：(1)利用势能场模型概念建立了适用于改航场景下的代价函数，考虑了燃油经济性、飞行受限区等约束，建立了可扩展的代价函数，具备限制条件的扩展性；(2)势能场模型为宏观计算模型，通过建立势能场可为大量航班的改航活动计算改航最优路径；(3)设计了并行化的计算算法，可以充分利用gpu的运算特点与效率，使计算性能提高5倍左右。本发明首先将环境场景栅格化，提取危险区域的范围，对环境场景进行分割，构建环境的不适度场；构建速度场，建立栅格中的航班密度贡献函数，计算航班沿上、下、左、右方向飞行时可达到的最大速度；建立代价函数，构建与路径长度、耗费时间、不适度三个指标相关的函数，计算单个航班选择前往目标点路径时的代价；建立离散势能场计算函数，利用单侧差分方程对势能场进行离散化计算，以势能场的负梯度方向作为改航最优路径；并行化势能场求解，采用分区运算规则，并行迭代求解单元格的势能值直至收敛。本发明有效提高了计算效率。
- 还没有人留言评论。精彩留言会获得点赞！