一种航空器运行可视化方法与流程

1.本发明属于空中交通管理技术领域，特别涉及一种基于网格的航空器运行可视化方法。


背景技术：

2.随着国民经济的快速发展，近些年中国民航业也迎来了飞速发展的时期，截至2020年，中国民航旅客运输量已连续15年稳居世界第二。民航事业的快速发展，使得空中交通流量大大增加。为了保证民航飞行的安全，需要飞行管制员在显示屏上时刻关注空中航空器的运行状态。
3.目前表征航空器运行的方法为：在航线规划过程中预先确定好起始点、转弯点和终止点，其中起始点到转弯点、转弯点到转弯点、转弯点到终止点之间的连线称为航段，将航段首尾连接即为航线。飞行管制员通过获取的雷达信号在飞行管理系统的显示屏上监控航空器的运行状态。由于航空器在空中的空间位置随着时间的推移是不断变化的，因此航空器的运行态势需要时间维度和空间维度共同表征，目前的飞行管理系统中表征航空器运行的方法是在二维(水平面)或三维(水平面 高度)表征航空器的运行状态，飞行管制员很难通过显示屏实时直观地看到航空器的四维运行状态以及航空器运行到各个点的时刻。


技术实现要素：

4.为了克服现有技术的不足，实现在三维空间中显示四维信息，需要通过一定的方法显示时间信息。本发明提供一种航线网格化的方法，网格为四维时空网格，通过四维时空网格表征出航空器的运行状态，包括：
5.步骤1：确定地球表面二维平面网格的最小边长，即根据航空器的巡航速度确定网格的最小边长；
6.步骤2：建立面向航空器运行的地球表面剖分网格体系，即将地球表面剖分为多层级、多尺度的网格，然后对各层级的网格进行编码；
7.步骤3：对网格高度进行剖分并编码，即覆盖航空器所有的飞行高度进行剖分编码；
8.步骤4：对航空器运行航线进行二维平面网格化表征，即将连续的航线离散化，用无重叠、无缝隙、首尾相接的一串网格表征出来，每个网格代表航空器不同时刻所在的位置；
9.步骤5：构建面向航空器运行的四维时空网格，即将航段经过二维网格化后，构建三维立体空间网格，进而构建四维时空网格；
10.步骤6：对四维时空网格着色，即将不同时间信息的时空网格着上不同的颜色。
11.进一步的，
12.步骤1中航空器的巡航速度范围为v1～v2千米/小时，v1<v2，则地球表面二维平面网格的最小边长的范围为a～b，其中：
[0013][0014]
步骤2中地球表面剖分网格体系剖分层级设定为7级；
[0015]
步骤3中选取0～40000米的高度范围进行剖分编码，将高度层级设定为7级；
[0016]
步骤4包括：
[0017]
步骤4
‑
1：选取航段；
[0018]
步骤4
‑
2：计算各航段的外包矩形；
[0019]
步骤4
‑
3：计算各航段的最小外包网格及编码，得到各航段的最小外包网格层级数level_n；
[0020]
步骤4
‑
4：在各航段的最小外包网格内建立直角坐标系；
[0021]
步骤4
‑
5：计算每个航段所在第7层级网格的行列坐标集合；
[0022]
步骤4
‑
6：计算出表征各航段的第7层级二维平面网格编码的集合；
[0023]
步骤4
‑
7：航段网格化；
[0024]
步骤5包括：
[0025]
步骤5
‑
1：构建三维立体空间网格；
[0026]
步骤5
‑
2：构建四维时空网格；
[0027]
步骤6中将航空器运行时刻与色光的r值、g值和b值对应。
[0028]
进一步的，
[0029]
步骤1中v1＝780千米/小时、v2＝1000千米/小时，地球表面二维平面网格的最小边长a＝210米、b＝270米；
[0030]
步骤2包括：
[0031]
步骤2
‑
1：进行第1层级网格剖分与编码；
[0032]
步骤2
‑
2：进行第2～6层级网格剖分与编码；
[0033]
步骤2
‑
3：进行第7层级剖分与编码；
[0034]
步骤2
‑
4：将航空器所在经纬度坐标转化为网格编码；
[0035]
步骤3包括：
[0036]
步骤3
‑
1：第1层级高度范围编码；
[0037]
步骤3
‑
2：第2层级高度范围编码；
[0038]
步骤3
‑
3：第3层级高度范围编码；
[0039]
步骤3
‑
4：进行第4到第7层级高度范围编码；
[0040]
步骤3
‑
5：确定出航空器的飞行高度所在的第7层级高度范围。
[0041]
进一步的，
[0042]
步骤2
‑
4中将航空器所在经纬度坐标转化为网格编码包括：
[0043]
步骤2
‑4‑
1：计算第1层级网格编码；
[0044]
步骤2
‑4‑
2：计算第2层级网格编码；
[0045]
步骤2
‑4‑
3：计算第3层级网格编码；
[0046]
步骤2
‑4‑
4：计算出第4层级网格的编码；
[0047]
步骤2
‑4‑
5：计算出第5层级网格的编码；
[0048]
步骤2
‑4‑
6：计算出第6层级网格的编码；
[0049]
步骤2
‑4‑
7：计算出第7层级网格的编码。
[0050]
进一步的，
[0051]
步骤4
‑
6中设某个第7层级网格的行列坐标为(x，y)，然后根据行列坐标(x，y)计算得到网格第level_n 1～7级编码计算出表征各航段的第7层级二维平面网格编码的集合包括：
[0052]
步骤4
‑6‑
1：计算坐标(x，y)对应的网格第7级编码；
[0053]
步骤4
‑6‑
2：计算坐标(x，y)对应的网格第6级编码；
[0054]
步骤4
‑6‑
3：计算坐标(x，y)对应的网格第5级编码；
[0055]
步骤4
‑6‑
4：计算坐标(x，y)对应的网格第4级编码；
[0056]
步骤4
‑6‑
5：计算坐标(x，y)对应的网格第3级编码；
[0057]
步骤4
‑6‑
6：计算坐标(x，y)对应的网格第2级编码。
[0058]
进一步的，
[0059]
步骤2
‑
1中选择本初子午线与赤道交点为剖分原点，用4
°×4°
的第1层级网格将纬度范围为[
‑
88
°
,88
°
]、经度范围为[
‑
180
°
,180
°
]的地球表面剖分为44
×
90份；
[0060]
步骤2
‑
2中层级数小的为上级，层级数大的为下级，上一层级网格为高层级网格，作为父网格，下一层级的网格为低层级网格，作为子网格，每1个上一层级父网格平均分成16份，得到下一层级16个子网格，对16个子网格按照皮亚诺空间填充曲线进行16进制编码，直到剖分至第6层级；
[0061]
步骤2
‑
3中将第6层级网格平均分成4份，得到第7层级网格。
[0062]
进一步的，
[0063]
横坐标为数字自西向东从00～89，间隔为1，共90个数字；纵坐标中北纬从低到高为字母a～v，共22个大写英文字母；南纬从低到高为字母a～v，共22个小写英文字母；第1层级网格序号在纬度方向上用a～v和a～v共44个英文字母从赤道按照纬度由低到高编码，其中北纬大写、南纬小写；第1层级网格序号在经度方向上从本初子午线开始由西向东用00～89共90个数字编码。
[0064]
进一步的，
[0065]
步骤2
‑
4中设航空器所在的经纬度坐标为(l，b)，其中l为纬度，b为经度，l和b的表示形式均为度分秒，记作l＝l
d
°
l
m
′
l
s
″
,b＝b
d
°
b
m
′
b
s
″
，其中l
d
为纬度值的度数，l
m
为纬度值的分数，l
s
为纬度值的秒数；b
d
为经度值的度数，b
m
为经度值的分数，b
s
为经度值的秒数；
[0066]
步骤2
‑4‑
1中计算第1层级网格编码，计算公式为：
[0067][0068]
l
d/
/4
°
的商为m1、余数为n1，m1 1对应的英文字母即为该点所在的第1层级网格的纬度方向编码，北纬大写，南纬小写；b
d/
/4
°
商为m1、余数为n1，m1即为该点所在的第1层级网格的经度方向编码；
[0069]
步骤2
‑4‑
2中计算第2层级网格编码，计算公式为：
[0070][0071]
n1为步骤2
‑4‑
1的余数，m2即为该点所在的第2层级网格的纬度方向编号；n1为步骤2
‑4‑
1的余数，m2即为该点所在的第2层级网格的经度方向编号；
[0072]
步骤2
‑4‑
3中计算第3层级网格编码，计算公式为：
[0073][0074]
l
m
/15
′
的商为m3、余数为n3，m3即为该点所在的第3层级网格的纬度方向编号；b
m
/15
′
商为m3、余数为n3，m3即为该点所在的第3层级网格的经度方向编号；
[0075]
步骤2
‑4‑
4，计算出第4层级网格的编码，计算公式为：
[0076][0077]
n3为步骤2
‑4‑
3的余数，故n3/3.75
′
的商为m4、余数为n4，m4即为该点所在的第4层级网格的纬度方向编号；n3为步骤2
‑4‑
3的余数，n3/3.75
′
商为m4、余数为n4，m4即为该点所在的第4层级网格的经度方向编号；
[0078]
步骤2
‑4‑
5中计算出第5层级网格的编码，计算公式为：
[0079][0080]
步骤2
‑4‑
4的余数为n4，故(n4×
60 l
s
)/56
″
的商为m5，余数为n5，m5即为该点所在的第5层级网格的纬度方向编号；(n4×
60 b
s
)/56
″
商为m5、余数为n5，m5即为该点所在的第5层级网格的经度方向编号；
[0081]
步骤2
‑4‑
6中计算出第6层级网格的编码，计算公式为：
[0082][0083]
n5/14
″
的商为m6、余数为n6，m6即为该点所在的第6层级网格的纬度方向编号；n5/14
″
商为m6、余数为n6，m6即为该点所在的第6层级网格的经度方向编号；
[0084]
步骤2
‑4‑
7中计算出第7层级网格的编码，计算公式为：
[0085][0086]
n6/7
″
的商为m7，余数为n7，m7即为该点所在的第7层级网格的纬度方向编号；n6/7
″
商为m7，余数为n7，m7即为该点所在的第7层级网格的经度方向编号。
[0087]
进一步的，
[0088]
步骤3
‑
1中将整个0～40000米的高度范围进行二等分，高度在[0米，20000米)之内编码为0，高度在[20000米，40000米]之内编码为1；
[0089]
步骤3
‑
2对第1层级上下2个高度范围编码均进行二等分，第2层级的高度差为10000米；将[0米，20000米)的高度范围进行二等分，高度在[0米，10000米)之内编码为00，高度在[10000米，20000米)之内编码为01；同理，将[20000米，40000米]的高度范围进行二等分，高度在[20000米，30000米)之内编码为10，高度在[30000米，40000米]之内编码为11；
[0090]
步骤3
‑
4中第4到第7层级高度差分别为2500米、1250米、625米、312.5米；层级数小的为上级，层级数大的为下级，下一层级在上一层级各个高度范围内均进行二等分，其中，二等分之后高度范围较低的在上一层级高度范围编码后面加0，高度范围较高的在上一层级高度范围编码后面加1；上下级高度范围编码之间具有继承性，编码的长度代表该高度范围所在的层级；
[0091]
步骤3
‑
5中由高度值h计算出该高度值所在的第7层级高度范围编码的方法为：h除以312.5，若能整除，所得结果为t；若不能整除，所得结果的整数部分为t，小数部分为t；将十进制数值t转化为7位二进制编码，即为高度值h所在的第7层级高度范围编码。
[0092]
进一步的，
[0093]
步骤6中r值对应小时，g值对应分钟，b值对应秒；r值按照步长为8从0开始取30个值，最小为0，最大为232；g值和b值按照步长为4取60个值，最小为0，最大为236。
[0094]
本发明方法通过将航空器的运行时刻与四维时空网格的颜色rgb值一一对应，为航空器不同运行时刻所在的网格着上不同的颜色，在显示屏上直观地一体化显示出四维时空信息，实现航空器飞行过程中四维状态可视化，方便管制员及时了解航空器的位置和运行时刻，能够为管制员提供直观、动态、全面的航空器运行状态显示，使得管制员对航空器未来一段时间的飞行动态有一个全局的认识，更好地保证飞行安全。
附图说明
[0095]
图1为本发明方法流程框图；
[0096]
图2为地球表面剖分网格体系第1层级剖分和编码方案示意图；
[0097]
图3为地球表面剖分网格体系第2～7层级剖分和编码方案示意图；
[0098]
图4为地球表面剖分网格的完整编码结构图；
[0099]
图5为网格高度剖分与编码方案示意图；
[0100]
图6为航空器运行航线的二维平面网格化表征示意图；
[0101]
图7为确定各航段的外包矩形示意图；
[0102]
图8为确定各航段的最小外包网格示意图；
[0103]
图9为在各航段的最小外包网格内建立直角坐标系示意图；
[0104]
图10为基于四维时空网格着色的航空器运行可视化表征示意图；
[0105]
图11为本发明的一个具体实施例中航段的二维平面网格化表征示意图；
[0106]
图12为本发明的一个具体实施例中对某航段进行四维时空网格表征的场景图。
具体实施方式
[0107]
为使本发明的目的、技术方案以及优势更加明晰，下面结合附图和实施例对本发明的内容进行进一步详细说明。
[0108]
本发明的基本原理是通过构造面向航空器运行的四维时空网格体系，将连续的航线进行离散化，从而能够显示出航空器每一秒的运行位置，并将航空器到达各网格的预计时间的时分秒三个参数与色光的rgb三个参数一一对应，对航空器运行过程中不同时刻所在的四维时空网格进行着色，直观地显示出航空器到达不同位置的预计到达时刻。
[0109]
如图1所示，本发明方法包括以下步骤：
[0110]
步骤1：确定地球表面二维平面网格的最小边长。
[0111]
地球表面二维平面网格的最小边长为航空器巡航阶段每秒的飞行距离。由于航空器的速度会根据安全、经济等因素综合决定，形成一种相对固定的常态速度来飞行，即巡航状态。巡航状态对应的速度为巡航速度。航空器在飞行过程中，单位距离上以巡航速度飞行所消耗的燃油量最少，而且巡航状态在整个飞行中占据时间最长，所以保持巡航速度飞行才最划算。设航空器的巡航速度范围为v1～v2千米/小时，v1<v2，则网格的最小边长(单位为米)的范围为a～b，其中：
[0112][0113]
民航航班的巡航速度一般为780～1000千米/小时，如空客a320系列巡航速度828km/小时，a340
‑
200巡航速度为896km/小时，而波音系列的巡航速度为900公里/小时左右，比如747
‑
800正常巡航速度是917km/小时。若航班以780～1000千米/小时的巡航速度飞行，每分钟的飞行距离为13千米～16.7千米，每秒的飞行距离为210米～270米。为了更加精确地表征出航空器的运行时刻，本发明将网格所表征的时间信息精确到秒级，因此二维平面网格的最小边长为210米～270米，二维平面网格为边长210米～270米的正方形。目前国际上主流的民航航空器的尺寸如表1所示。二维平面网格的最小尺寸能够完全包含表1中的主流航空器。
[0114]
表1
[0115]
型号长宽(翼展)高空客a38072.75米79.75米24.09米空客a32037.57米34.10米11.76米波音747
‑
876.4米68.5米19.4米波音73728.6米28.3米11.3米c91937.57米34.10米11.76米
[0116]
步骤2：建立面向航空器运行的地球表面剖分网格体系。即将地球表面剖分为多层级、多尺度的网格，各层级网格有对应的面积，然后按照编码规则对各层级的网格进行编码。地球表面的任一网格都有唯一的编码，根据编码能够索引出网格，即一个编码确定出唯一的网格所在的位置，且网格编码与经纬度坐标之间能够进行快速转换。步骤2具体为：
[0117]
步骤2
‑
1：进行第1层级网格剖分与编码。如图2所示，选择本初子午线与赤道交点为剖分原点，用4
°×4°
的第1层级网格将纬度范围为[
‑
88
°
,88
°
]、经度范围为[
‑
180
°
,180
°
]的地球表面剖分为44
×
90份。
[0118]
在本发明的一个具体实施例中，横坐标为数字自西向东从00～89，间隔为1，共90个数字；纵坐标中北纬从低到高为字母a～v，共22个大写英文字母；南纬从低到高为字母a～v，共22个小写英文字母。第1层级网格序号在纬度方向上用a～v和a～v共44个英文字母从赤道按照纬度由低到高编码，其中北纬大写、南纬小写；第1层级网格序号在经度方向上从本初子午线开始由西向东用00～89共90个数字编码。因为地球的平均半径为约为6371公里，地球的平均周长为2π
×
6371≈40030千米，因此1
°
网格的边长为2π
×
6371
÷
360≈111千米，故第1层级网格(4
°×4°
网格)的边长约为444千米。
[0119]
步骤2
‑
2：进行第2～6层级网格剖分与编码。层级数小的为上级，层级数大的为下级，上一层级网格为高层级网格，作为父网格，下一层级的网格为低层级网格，作为子网格。每1个上一层级父网格平均分成16份，得到下一层级16个子网格，对16个子网格按照皮亚诺空间填充曲线进行16进制编码(0123456789abcdef)，直到剖分至第6层级，如图3所示。所以子网格编码继承上一层级父网格的编码，并比上一层级父网格编码多1位。第2～6级网格的大小分别为1
°×1°
、15
′×
15
′
、3.75
′×
3.75
′
、56
″×
56
″
、14
″×
14
″
，网格边长分别为111千米、28千米、7千米、1.75千米、440米。
[0120]
在本发明的一个具体实施例中，第3层级某15
′×
15
′
网格的编码为j29ab，将其平均分成16个3.75
′×
3.75
′
的第4层级网格，对这16个子网格按照图3中的编码方式进行编码，编码分别为{j29ab0，j29ab1，j29ab2，j29ab3，j29ab4，j29ab5，j29ab6，j29ab7，j29ab8，j29ab9，j29aba，j29abb，j29abc，j29abd，j29abe，j29abf}。
[0121]
步骤2
‑
3：进行第7层级剖分与编码。因为第6层级网格的边长为440米，如果依然平均分成16份，则第7层级网格的变成为110米，不在步骤1确定的空间网格最小尺寸在210米～270米的范围内，因此第7层级网格的剖分方法不能按照第2～6层级网格的剖分方式。此时，选择将第6层级网格平均分成4份，得到第7层级网格的边长为220米，该尺度在210米～270米范围内，兼顾了航空器每秒的飞行距离和航空器的尺寸，同时第7层级网格的剖分方法减少了计算量，提高了计算机显示的速度，故选定第7层级网格的尺寸为220米
×
220米。地球表面剖分网格体系将剖分层级最终设定为7级。网格的层级数越小，代表网格的层级越靠上、层级越高，整个7层级剖分网格呈金字塔形。网格的层级数level_n与该层级网格编码长度l_code的对应关系为：
[0122]
l_code＝level_n 2
[0123]
地球表面剖分网格的完整编码结构如图4所示。由于地球表面剖分网格编码的继承关系，第level_n层级地球表面剖分网格的编码结构为编码第1级到第level_n级。
[0124]
经过步骤2
‑
1、2
‑
2和2
‑
3，将地球剖分成无缝隙、无重叠的网格，每个网格都有唯一的编码以及相对应的面积。各层级网格大小和边长见表2。
[0125]
表2
[0126][0127]
步骤2
‑
4：将经纬度坐标转化为网格编码。设航空器所在的经纬度坐标为(l，b)，其中l为纬度，b为经度，l和b的表示形式均为度分秒，记作l＝l
d
°
l
m
′
l
s
″
,b＝b
d
°
b
m
′
b
s
″
，其中l
d
为纬度值的度数，l
m
为纬度值的分数，l
s
为纬度值的秒数；b
d
为经度值的度数，b
m
为经度值的分数，b
s
为经度值的秒数。如图3所示，第1层级网格中纵坐标代表纬度方向编码，大写字母为北纬，小写字母为南纬，横坐标代表经度方向编码00～89，图3中只显示了部分第1层级网格及编码。第2～6层级网格中纵坐标代表纬度方向，编号为0、1、2、3；横坐标代表经度方向，编号为0、1、2、3。坐标顺序为先纬向、再经向。第2～6层级坐标与编码的对应关系如表3所示：
[0128]
表3
[0129]
坐标对应的编码坐标对应的编码(0，0)0(2，0)8(0，1)1(2，1)9(1，0)2(3，0)a(1，1)3(3，1)b(0，2)4(2，2)c(0，3)5(2，3)d(1，2)6(3，2)e(1，3)7(3，3)f
[0130]
第7层级网格中纵坐标代表纬度方向，编号为0、1；横坐标代表经度方向，编号为0、1，坐标顺序同样为先纬向、再经向。第7层级坐标与编码的对应关系如表4所示：
[0131]
表4
[0132]
坐标对应的编码坐标对应的编码(0，0)0(1，0)2(0，1)1(1，1)3
[0133]
经纬度坐标转化为网格编码的详细步骤为：
[0134]
步骤2
‑4‑
1：计算第1层级网格编码。计算公式为：
[0135][0136]
由于第1层级网格为4
°×4°
网格，故计算出l
d
/4
°
的商为m1、余数为n1，m1 1对应的英文字母即为该点所在的第1层级网格的纬度方向编码，北纬大写，南纬小写；b
d
/4
°
商为m1、余数为n1，m1即为该点所在的第1层级网格的经度方向编码。两个编码按照先纬向、后经向的顺序组合成该点所在的第1层级网格编码。例如北京天坛祈年殿的经纬度坐标为(39
°
53
′
24
″
n，116
°
25
′
09
″
e)，39
°
/4
°
＝9余3
°
，116
°
/4
°
＝29余0
°
，位于第1层级中北纬向第10个、经向第30个4
°×4°
网格，北纬第10个编码为大写英文字母j，所以北京天坛祈年殿所在的第1层级网格的编码为j29。
[0137]
步骤2
‑4‑
2：计算第2层级网格编码。计算公式为：
[0138][0139]
由于第2层级网格为1
°×1°
网格，故步骤2
‑4‑
1的余数n1除以1
°
得到数值m2，m2即为该点所在的第2层级网格的纬度方向编号；n1除以1
°
得到数值m2，m2即为该点所在的第2层级网格的经度方向编号，根据两个编号形成的坐标(m2，m2)按照表3对应出第2级编码，加在第1层级编码后面得到第2层级网格编码。上例中，3
°
/1
°
＝3，0
°
/1
°
＝0，坐标(3，0)对应的编码为a，所以北京天坛祈年殿所在的第2层级网格的编码为j29a。
[0140]
步骤2
‑4‑
3：计算第3层级网格编码。计算公式为：
[0141][0142]
第3层级网格为15
′×
15
′
网格，故l
m
/15
′
的商为m3、余数为n3，m3即为该点所在的第3层级网格的纬度方向编号；b
m
/15
′
商为m3、余数为n3，m3即为该点所在的第3层级网格的经度方向编号。根据两个编号形成的坐标(m3，m3)按照表3对应出第3级编码，加在第2层级编码后面得到第3层级网格编码。例如53
′
/15
′
＝3余8
′
，25
′
/15
′
＝1余10
′
，坐标(3，1)对应的编码为b，所以北京天坛祈年殿所在的第3层级网格的编码为j29ab。
[0143]
步骤2
‑4‑
4：计算出第4层级网格的编码。计算公式为：
[0144][0145]
第4层级网格为3.75
′×
3.75
′
网格，故n3/3.75
′
的商为m4、余数为n4，m4即为该点所在的第4层级网格的纬度方向编号；n3/3.75
′
商为m4、余数为n4，m4即为该点所在的第4层级网格的经度方向编号。根据两个编号形成的坐标(m4，m4)按照表3对应出第4级编码，加在第3层级编码后面得到第4层级网格编码。例如8
′
/3.75
′
＝2余0.5
′
，10
′
/3.75
′
＝2余2.5
′
，坐标(2，2)对应的编码为c，所以北京天坛祈年殿所在的第4层级网格的编码为j29abc。
[0146]
步骤2
‑4‑
5：计算出第5层级网格的编码。计算公式为：
[0147][0148]
第5层级网格为56
″×
56
″
网格，故(n4×
60 l
s
)/56
″
的商为m5，余数为n5，m5即为该点所在的第5层级网格的纬度方向编号；(n4×
60 b
s
)/56
″
商为m5、余数为n5，m5即为该点所在的第5层级网格的经度方向编号。根据两个编号形成的坐标(m5，m5)按照表3对应出第5级编码，加在第4层级编码后面得到第5层级网格编码。例如(0.5
×
60 24)
″
/56
″
＝0余54
″
，(2.5
×
60 9)
″
/56
″
＝2余47
″
，坐标(0，2)对应的编码为4，所以北京天坛祈年殿所在的第5层级网格的编码为j29abc4。
[0149]
步骤2
‑4‑
6：计算出第6层级网格的编码。计算公式为：
[0150][0151]
第6层级网格为14
″×
14
″
网格，故n5/14
″
的商为m6、余数为n6，m6即为该点所在的第6层级网格的纬度方向编号；n5/14
″
商为m6、余数为n6，m6即为该点所在的第6层级网格的经度方向编号。根据两个编号形成的坐标(m6，m6)按照表3对应出第6级编码，加在第5层级编码后面得到第6层级网格编码。例如54
″
/14
″
＝3余12
″
，47
″
/14
″
＝3余5
″
，坐标(3，3)对应的编码为f，所以北京天坛祈年殿所在的第6层级网格的编码为j29abc4f。
[0152]
步骤2
‑4‑
7：计算出第7层级网格的编码。计算公式为：
[0153][0154]
第7层级网格为7
″×7″
网格，故n6/7
″
的商为m7，余数为n7，m7即为该点所在的第7层级网格的纬度方向编号；n6/7
″
商为m7，余数为n7，m7即为该点所在的第7层级网格的经度方向编号。根据两个编号形成的坐标(m7，m7)按照表4对应出第7级编码，加在第6层级编码后面得到第7层级网格编码。由于剖分到第7层级为止，不再往下进行剖分，故余数n7和n7不再进行继续运算。例如12
″
/7
″
＝1余5
″
，5
″
/7
″
＝0余7
″
，坐标(1，0)对应的编码为2，所以北京天坛祈年殿所在的第7层级网格的编码为j29abc4f2。
[0155]
步骤3：对网格高度进行剖分并编码。航空器运行是有一定的高度范围的，航空器的飞行高度大多在20000米以下，少数类型的航空器飞行高度能超过20000米，本发明为了保证网格的高度范围能够覆盖航空器所有的飞行高度，选取0～40000米的高度范围进行剖分编码。网格高度剖分编码体系如图5所示，详细过程如下：
[0156]
步骤3
‑
1：第1层级高度范围编码。第1层级的高度差为20000米。即将整个0～40000米的高度范围进行二等分，高度在[0米，20000米)之内编码为0，高度在[20000米，40000米]之内编码为1。在本说明书中，符号“[”和“]”代表包含区间的端点，符号“)”代表不包含区间的端点。
[0157]
步骤3
‑
2：第2层级高度范围编码。对第1层级上下2个高度范围编码均进行二等分，第2层级的高度差为10000米。将[0米，20000米)的高度范围进行二等分，高度在[0米，10000米)之内编码为00，高度在[10000米，20000米)之内编码为01。同理，将[20000米，40000米]的高度范围进行二等分，高度在[20000米，30000米)之内编码为10，高度在[30000米，40000米]之内编码为11。
[0158]
步骤3
‑
3：第3层级高度范围编码。对第2层级的4个高度范围编码均进行二等分，第3层级的高度差为5000米。将[0米，10000米)的高度范围进行二等分，高度在[0米，5000米)之内编码为000，高度在[5000米，10000米)之内编码为001；将[10000米，20000米)的高度范围进行二等分，高度在[10000米，15000米)之内编码为010，高度在[15000米，20000米)之内编码为011；将[20000米，30000米)的高度范围进行二等分，高度在[20000米，25000米)之内编码为100，高度在[25000米，30000米)之内编码为101；将[30000米，40000米]的高度范围进行二等分，高度在[30000米，35000米)之内编码为110，高度在[35000米，40000米]之内编码为111。
[0159]
步骤3
‑
4：进行第4到第7层级高度范围编码。第4到第7层级高度差分别为2500米、1250米、625米、312.5米。层级数小的为上级，层级数大的为下级，下一层级在上一层级各个高度范围内均进行二等分，其中，二等分之后高度范围较低的在上一层级高度范围编码后面加0，高度范围较高的在上一层级高度范围编码后面加1。上下级高度范围编码之间具有继承性，编码的长度代表该高度范围所在的层级。
[0160]
步骤3
‑
5：确定出航空器的飞行高度所在的第7层级高度范围。经过步骤3
‑
1到3
‑
4，每个高度范围都对应唯一的编码，编码与高度范围一一对应。航空器的飞行高度一定在某个第7层级高度范围内，计算出航空器飞行高度所在的第7层级高度范围编码，就能根据该编码对应出航空器飞行高度所在的第7层级高度范围。由高度值h计算出该高度值所在的第7层级高度范围编码的方法为：h除以312.5，若能整除，所得结果为t；若不能整除，所得结果的整数部分为t，小数部分为t。将十进制数值t转化为7位二进制编码，即为高度值h所在的第7层级高度范围编码，由该编码对应出第7层级高度范围。举例说明高度剖分编码，某航空器的巡航高度为9000米，9000/312.5＝28.8，将十进制数值28转化7位二进制编码为0011100，故9000米所在的第7层级高度编码为0011100，编码0011100对应的第7层级高度范围为[8750米，9062.5米)。
[0161]
步骤4：对航空器运行航线进行二维平面网格化表征。即将连续的航线离散化，用无重叠、无缝隙、首尾相接的一串网格表征出来，每个网格代表航空器不同时刻所在的位置。将航线表征为一串首尾相连的网格的集合，航线的二维平面网格化表征如图6所示，详细步骤包括：
[0162]
步骤4
‑
1：选取航段。航空器在飞行过程中，除了在转弯点处转弯外，各转弯点之间的航线为直线，故航线可看作由不同段的直线段连接而成，即航段，各航段的连接点为航线的转弯点。
[0163]
步骤4
‑
2：计算各航段的外包矩形。将航空器的预计飞行航线在不同转弯点处分割为不同的航段。如图7所示，各航段的起始点为p1(l1，b1)，终止点为p2(l2，b2)，其中l1为起始点p1的纬度，b1为起始点p1的经度；l2为终止点p2的纬度，b2为终止点p2的经度。方向规定为“上北下南，左西右东”。若l1≠l2且b1≠b2，则过航段的起始点和终止点分别作两条南北向
的直线，同时，过各航段的起始点和终止点分别作两条东西向的直线，4条直线确定的矩形即为该航段的外包矩形。若某航段起始点和终止点的经纬度坐标关系为l1＝l2或者b1＝b2，说明该航段的方向为东西向或者南北向，则无需计算该航段的外包矩形，省略此步。
[0164]
步骤4
‑
3：计算各航段的最小外包网格及编码。对于l1≠l2且b1≠b2的航段，计算出航段的外包矩形的面积s＝d1×
d2，其中d1为外包矩形的长，d2为外包矩形的宽。由于第1～7层级网格都有固定的面积，比较出网格面积大于s且与s最接近的网格层级level_n。对于l1＝l2或者b1＝b2的航段，比较出第1～7层级网格中的边长大于该航段的长度，并且与该航段长度最接近的网格层级level_n。在各航段中都任意取一点p3，按照步骤2
‑
4的方法计算出该点所在的第level_n层级网格编码，该编码对应的网格即为该航段的最小外包网格。以图8中的航段为例，说明最小外包网格的确定方法。图8中航段的外包矩形的宽d1＝7千米，长d2＝5千米，则外包矩形面积为35平方千米。因为第4层级网格边长为7千米、面积为49平方千米，大于35平方千米且与35平方千米最接近，所以该航段的最小外包网格对应的层级为4，如图8所示。假设在该航段中任意选择一点为p3，将p3的经纬度坐标按步骤2
‑
4的方法转化为第4层级的网格编码为k06e75，该编码对应的网格即为该航段的最小外包网格。采用上述方法，计算得到所有航段的最小外包网格。
[0165]
步骤4
‑
4：在各航段的最小外包网格内建立直角坐标系。各航段的最小外包网格的层级数为level_n，将各航段的最小外包网格剖分至第7层级，则各航段的最小外包网格内剖分得到的第7层级网格数量为166‑
level_n
×
4。以最小外包网格的左下顶点所在的第7层级网格为坐标原点，建立直角坐标系，向东，向北为正向，网格横纵坐标轴依次用0，1，2，3，4，5，6，
…
，编号，坐标轴上的坐标代表不同的第7层级网格所在的位置，例如左下角网格的坐标为(0,0)。以图9为例，粗黑线为某航段，该航段的最小外包网格层级数level_n为5，将该航段的最小外包网格剖分至第7层级，得到16
×
4共64个第7层级网格，横纵坐标均为0，1，2，3，4，5，6，7。
[0166]
步骤4
‑
5：计算每个航段所在第7层级网格的行列坐标集合。确定每个航段起始点和终止点的行列坐标后，根据各航段的航向，计算出各航段的斜率。飞机纵轴前方的延长线叫航向线，飞机所在位置的经线北端顺时针测量至航向线的夹角叫做航向，航向范围为[0,360
°
)。由航向c计算得到航段斜率k的公式为：
[0167][0168]
当航向为0
°
时，代表航向为正北，航段的斜率为∞，设步骤4
‑
4得到起始点坐标为(x1，y1)，终止点坐标为(x1，y2)，则此时航段起始点和终止点之间的网格坐标的横坐标均为x1，纵坐标从y1到y2按照步长为1增加，得到航段起始点和终止点之间的网格坐标集合。当航向为180
°
时，代表航向为正南，航段的斜率为∞，设航段的起始点坐标为(x1，y1)，终止点坐标为(x1，y2)，则此时航段起始点和终止点之间的网格坐标的横坐标均为x1，纵坐标从y1到y2按照步长为
‑
1减小，得到航段起始点和终止点之间的网格坐标集合。当航向不为0
°
或者
180
°
，计算出航段斜率k，则各航段起始点和终止点之间的网格坐标集合可以可由直线插值计算得到。对于航向为0
°
或者180
°
的航段，该航段上某时刻航空器的运行位置点正好位于东西两个网格的分界线上，则统一规定用东边网格的坐标代表所在的第7层级网格的行列坐标；同理对于航向为90
°
或者270
°
的航段，该航段上某时刻航空器的运行位置点正好位于南北两个网格的分界线上，则统一规定用北边网格的坐标代表所在的第7层级网格的行列坐标。
[0169]
步骤4
‑
6：计算出表征各航段的第7层级二维平面网格编码的集合。在步骤4
‑
5得到表征航段的第7层级网格的行列坐标集合的基础上，利用网格的行列坐标计算出网格的编码，求得表征各航段的第7层级二维平面网格编码的集合。在之前步骤4
‑
3中计算得到了各航段的最小外包网格及编码，各航段的最小外包网格层级数level_n，其中1≤level_n≤7。设某个第7层级网格的行列坐标为(x，y)，然后根据行列坐标(x，y)计算得到网格第level_n 1～7级编码，采用“由低到高，逐层计算；由高到低，反向添加”的方法，即先计算出第7级编码，然后计算第6，5，
……
，level_n 1级的编码，然后根据各层级网格编码间的继承关系，依次将第level_n 1，level_n 2，
……
，7级的编码添加到航段的最小外包网格的编码后面。详细的计算方法为：
[0170]
步骤4
‑6‑
1：计算坐标(x，y)对应的网格第7级编码。计算公式为：
[0171][0172]
其中等式右边为计算结果，余数组合(b_7，b_7)即表4中的坐标，得到第7级编码。
[0173]
步骤4
‑6‑
2：计算坐标(x，y)对应的网格第6级编码。计算公式为：
[0174][0175]
其中等式右边为计算结果，余数组合(b_6，b_6)即表3中的坐标，得到第6级编码。
[0176]
步骤4
‑6‑
3：计算坐标(x，y)对应的网格第5级编码。计算公式为：
[0177][0178]
其中等式右边为计算结果，余数组合(b_5，b_5)即表3中的坐标，得到第5级编码。
[0179]
步骤4
‑6‑
4：计算坐标(x，y)对应的网格第4级编码。计算公式为：
[0180][0181]
其中等式右边为计算结果，余数组合(b_4，b_4)即表3中的坐标，得到第4级编码。
[0182]
步骤4
‑6‑
5：计算坐标(x，y)对应的网格第3级编码。计算公式为：
[0183][0184]
其中等式右边为计算结果，余数组合(b_3，b_3)即表3中的坐标，得到第3级编码。
[0185]
步骤4
‑6‑
6：计算坐标(x，y)对应的网格第2级编码。计算公式为：
[0186][0187]
其中等式右边为计算结果，余数组合(b_2，b_2)即表3中的坐标，得到第2级编码。
[0188]
根据航段的最小外包网格层级数level_n的数值，采用步骤4
‑6‑
1到步骤4
‑6‑
6的方法只需要计算出第7到level_n 1级的编码，然后依次将第level_n 1到第7级的编码添加在航段的最小外包网格编码后面，即得到坐标(x，y)对应的第7层级网格编码。按照上述方法，将各航段中所有第7层级网格的行列坐标逐个转化为第7层级网格的编码，最终实现各航段第7层级网格的行列坐标集合转化为第7层级网格编码集合。
[0189]
步骤4
‑
7：航段网格化。对各航段的网格编码集合求并集，得到整个航线的网格编码集合。由于网格编码与网格位置的一一对应关系，根据网格编码索引出相应的网格，将航段表征为网格的集合。如图9所示，图9中深色网格为某航段经过二维网格化表征的结果，最终将该航段表征为深色网格的集合。
[0190]
步骤5：构建面向航空器运行的四维时空网格。该步骤包括：
[0191]
步骤5
‑
1：构建三维立体空间网格。经过步骤4之后，航空器飞行航线可通过平面网格来表征，将航线中不同的航段的飞行高度，运用步骤3
‑
5中的方法，计算出飞行高度所在的第7层级高度范围，给各航段所在的平面网格加上相应的高度范围，构建三维立体空间网格。结合步骤1和步骤2确定的第7层级二维平面网格的尺寸220米
×
220米，步骤3确定的第7层级高度差为312.5米，因此三维立体网格的长宽高尺寸分别为220米
×
220米
×
312.5米。三维立体空间网格的编码结构为在二维平面网格编码后面加上高度编码。
[0192]
步骤5
‑
2：构建四维时空网格。由步骤1可知，第7层级网格的尺寸与航空器巡航阶段每秒钟的飞行距离相当，因此，用第7层级网格表征航空器的航线之后，航空器每移动一个网格，运行时间相应增加一秒钟，故每一个第7层级网格都对应唯一的航空器运行时刻。因此，步骤5
‑
1已经构建出的三维立体空间网格，实际上也蕴含航空器的时间信息，每一个小网格都能表征出航空器运行的空间和时间信息，即在三维网格中能表征出四维时空信息，在三维空间网格的基础上构建出四维时空网格。
[0193]
步骤6：对四维时空网格着色。将航空器运行时刻与rgb值对应，为包含不同时间信息的时空网格着上不同的颜色。
[0194]
因为任何一种色光都可以由红、绿、蓝按照一定的比例混合而成，因此，红、绿、蓝被称为光学三原色。这种由光学三原色混合后，组成显示屏显示颜色的模式称为三原色光模式(英语：rgb color model)，又称rgb颜色模型，是一种加色模型，通过调整色光的rgb
值，能够调成我们日常所能看到的各式各样的色光。将色光的r值、g值和b值分别与时间的时、分、秒三个属性对应起来。r值对应小时，g值对应分钟，b值对应秒。r值、g值和b值的可调范围均为0～255，调整的最小幅度为1，因此，r值、g值和b值能够分别对应256个不同的值。能够完全涵盖每小时对应的60分钟，以及每分钟对应的60秒。通过rgb值与时间的对应关系能够表征出航空器256个小时的全部运行时间。而实际上，考虑到航空器的最大航程以及机组人员的疲劳程度，世界上最长直飞航线——纽约到悉尼航线总长1.6万公里，共需飞行20小时，远远小于r值所能表征的256小时。故用rgb值表征时间完全能够表征出航空器完整的预计运行时间。
[0195]
本发明为了能够完整的表示出航空器在巡航过程中完整的预计运行过程，在全球最长直飞航线需飞行20小时的基础上，再加10小时，共表征出航空器自进入巡航状态后的30个小时的运行。由于r值的取值范围为0～255，为使不同时间对应的色光更具辨识度，因此按照步长为8从0开始取30个值，最小为0，最大为232，代表时间从第0小时到第29小时。同理因为1小时等于60分钟，1分钟等于60秒，所以在g值和b值的取值范围(0～255)内按照步长为4取60个值，最小为0，最大为236，分别对应0～59共60分钟以及0～59共60秒。rgb值与时、分、秒的对应关系如表5所示。
[0196]
表5
[0197]
[0198]
[0199][0200]
如图10所示，假设飞机于2021年2月1日早上8点进入巡航状态，航线上某点a所在的时空网格颜色的rgb值分别为(24,92,100)。查表可知，24代表航空器巡航阶段开始后的第3小时，92代表第23分钟，100代表第25秒，综上，该网格对应的时刻是航空器进入巡航状态的时刻加上3小时23分25秒，即航空器运行到该网格的预计到达时间为11时23分25秒。航线上另外一点b所在的时空网格颜色的rgb值分别为(48,120,112)。查表可知，48代表航空器巡航阶段开始后的第6小时，120代表第30分钟，112代表第28秒，综上，该网格对应的时刻是航空器进入巡航状态的时刻加上6小时30分28秒，即航空器运行到该网格的预计到达时间为14时30分28秒。
[0201]
将时间值与rgb值一一对应之后，即可为航空器运行所在的四维时空网格着色，用不同的颜色表征出航空器不同的运行时间，用不同网格的位置表征出航空器的运行位置，航空器的预计到达时间通过不同的色彩显示出来，能使飞行管制员能对航空器在未来一段时间的运行状态有更清楚的全局认识。
[0202]
如图11所示，在本发明的一个具体实施例中，航线上一小段航段为东西向，长为1100米，起始点的经纬度坐标为(45
°
21
′
20
″
n，110
°
37
′
10
″
e)，终止点的经纬度坐标为(45
°
21
′
20
″
n，110
°
37
′
45
″
e)，飞机航向为90
°
，飞行高度为10000米，航空器进入巡航状态1小时10分钟0秒后达到图11中航段的起始点。图11的横纵坐标轴编号按照步骤4
‑
4建立直角坐标系的方法确定，代表第7层级网格所在的行列编号，由横纵编号确定第7层级网格的行列坐标。
[0203]
在经过步骤1到步骤3之后，建立起三维网格体系。然后按照步骤4的流程将该航段进行二维平面网格化表征。由于该行段起始点和终止点的纬度相同，省略步骤4
‑
2，直接进入步骤4
‑
3，计算该航段的最小外包网格及编码。比较出第1～7层级网格中的边长大于该航段的长度，并且与该航段长度最接近的网格层级数。由于该航段的长度为1.1千米，第5层级网格的边长为1.75千米，大于该航段长度且与该行段的长度最接近，所以该航段的最小外包网格的层级数为5。在该航段上任意取一点，本实施例取航线起始点，按照步骤2
‑
4的方法，将起始点的经纬度坐标转化为第5层级网格编码，步骤为：
[0204]
步骤2
‑4‑
1：计算第1层级网格编码
[0205]
航段起始点的经纬度度坐标为(45
°
21
′
20
″
n，110
°
37
′
10
″
e)，45
°
/4
°
＝11余1
°
，110
°
/4
°
＝28余2
°
，位于第1层级中北纬向第12个、经向第29个4
°×4°
网格，北纬第12个编码为大写英文字母l，所以第1层级网格编码为l28。
[0206]
步骤2
‑4‑
2：计算第2层级网格编码
[0207]1°
/1
°
＝1，2
°
/1
°
＝2，坐标(1，2)按照表3对应的编码为6，所以第2层级网格的编码为l286。
[0208]
步骤2
‑4‑
3：计算第3层级网格编码
[0209]
21
′
/15
′
＝1余6
′
，37
′
/15
′
＝2余7
′
，坐标(1，2)按照表3对应的编码为6，所以第3层级网格的编码为l2866。
[0210]
步骤2
‑4‑
4：计算第4层级网格编码
[0211]6′
/3.75
′
＝1余2.25
′
，7
′
/3.75
′
＝1余2.25
′
，坐标(1，1)按照表3对应的编码为3，所以第4层级网格的编码为l28663。
[0212]
步骤2
‑4‑
5：计算第5层级网格编码
[0213]
(2.25
×
60 20)
″
/56
″
＝2余43
″
，(2.25
×
60 10)
″
/56
″
＝2余33
″
，坐标(2，2)按照表3对应的编码为c，所以第5层级网格的编码为l28663c。即该航段的最小外包网格编码为l28663c。
[0214]
按照步骤4
‑
4的方法在该航段的最小外包网格内建立直角坐标系。将该航段的最小外包网格剖分至第7层级，则最小外包网格内剖分得到的第7层级网格数量为16
×
4＝64个。在图11中，横纵坐标轴均依次为0，1，2，3，4，5，6，7。按照步骤4
‑
5的方法计算该航段所在第7层级网格的行列坐标集合为{(1,4)，(2,4)，(3,4)，(4,4)，(5,4)}，在图11中用深色网格表示。按照步骤4
‑
6的方法将第7层级网格的行列坐标集合转化为编码集合，以坐标(1,4)为例，转化方法为：
[0215]
步骤4
‑6‑
1：计算坐标(1,4)对应的网格第7级编码。
[0216]
4/2＝2余0，1/2＝0余1，(0,1)按照表4中的对应关系，得到第7级编码为1。
[0217]
步骤4
‑6‑
2：计算坐标(1,4)对应的网格第6级编码。
[0218]
2/4＝0余2，0/4＝0与0，(2,0)按照表3中的对应关系，得到第6级编码为8。
[0219]
将第6级、第7级编码依次添加在该航段的最小外包网格编码为l28663c后面，得到坐标(1,4)对应的第7层级网格编码为l28663c81。按照步骤4
‑6‑
1和步骤4
‑6‑
2的方法，将坐标集合中的坐标逐个转化为第7层级网格编码，最终得到表征该航段的第7层级二维平面网格编码的集合为{l28663c81，l28663c90，l28663c91，l28663cc0，l28663cc1}。根据编码集合索引出网格集合，得到该航段的二维平面网格化表征。
[0220]
进入步骤5，构建四维时空网格。首先，按照步骤5
‑
1构建三维立体空间网格。因为该航段的飞行高度为10000米，按照步骤3
‑
5计算该高度值所在的第7层级高度范围。10000/312.5＝32，将十进制数值32转化为7位二进制编码为0100000，故高度值10000米对应的第7层级高度范围编码为0100000，该编码对应的第7层级高度范围为[10000米，10312.5米)，在步骤4表征的二维平面网格上加入高度范围，构建出三维立体空间网格。表征该航段的三维立体空间网格的编码集合为{l28663c81
‑
0100000，l28663c90
‑
0100000，l28663c91
‑
0100000，l28663cc0
‑
0100000，l28663cc1
‑
0100000}。然后，按照步骤5
‑
2，根据网格所蕴含
的航空器运行时刻信息，构建四维时空网格。按照步骤6的方法给表征图11中航段的5个时空网格着色。根据表5的对应关系，该航段从起始点到终止点对应的5个时空网格的颜色的rgb值分别为(8，40，0)，(8，40，4)，(8，40，8)，(8，40，12)，(8，40，16)。
[0221]
如图12所示是按照步骤1～步骤6的原理，在功能验证软件中将某航段可视化表征为着色的四维时网格的示意图。该功能验证软件用javascript语言编写，运行环境：操作系统windows10，64位；cpu：intel core i7
‑
105102.3ghz；ram：16g。
[0222]
传统的航空器运行软件只能表征出航空器的二维或三维运行状态。从图12可以看出，本发明提供的方法通过将时间与颜色相对应，能够直观地表征出航空器运行的空间和时间信息，实现在三维网格内表征出四维时空信息的目的，不需要再额外的增加显示界面即可将时间和空间信息结合起来，能够更加直观的表征出航空器在空中的四维时空信息。将时间信息通过不同的色彩表现出来，能够使管制员集中注意力于航空器在空中的位置，有助于管制员及时发现航空器的实际运行状态和预计运行状态之间的差异，及时下达调整指令，保证飞行安全。

技术特征：
1.一种航空器运行可视化方法，包括：步骤1：确定地球表面二维平面网格的最小边长，即根据航空器的巡航速度确定网格的最小边长；步骤2：建立面向航空器运行的地球表面剖分网格体系，即将地球表面剖分为多层级、多尺度的网格，然后对各层级的网格进行编码；步骤3：对网格高度进行剖分并编码，即覆盖航空器所有的飞行高度进行剖分编码；步骤4：对航空器运行航线进行二维平面网格化表征，即将连续的航线离散化，用无重叠、无缝隙、首尾相接的一串网格表征出来，每个网格代表航空器不同时刻所在的位置；步骤5：构建面向航空器运行的四维时空网格，即将航段经过二维网格化后，构建三维立体空间网格，进而构建四维时空网格；步骤6：对四维时空网格着色，即将不同时间信息的时空网格着上不同的颜色。2.如权利要求1所述的一种航空器运行可视化方法，其特征在于：步骤1中航空器的巡航速度范围为v1～v2千米/小时，v1<v2，则地球表面二维平面网格的最小边长的范围为a～b，其中：步骤2中地球表面剖分网格体系剖分层级设定为7级；步骤3中选取0～40000米的高度范围进行剖分编码，将高度层级设定为7级；步骤4包括：步骤4
‑
1：选取航段；步骤4
‑
2：计算各航段的外包矩形；步骤4
‑
3：计算各航段的最小外包网格及编码，得到各航段的最小外包网格层级数level_n；步骤4
‑
4：在各航段的最小外包网格内建立直角坐标系；步骤4
‑
5：计算每个航段所在第7层级网格的行列坐标集合；步骤4
‑
6：计算出表征各航段的第7层级二维平面网格编码的集合；步骤4
‑
7：航段网格化；步骤5包括：步骤5
‑
1：构建三维立体空间网格；步骤5
‑
2：构建四维时空网格；步骤6中将航空器运行时刻与色光的r值、g值和b值对应。3.如权利要求2所述的一种航空器运行可视化方法，其特征在于：步骤1中v1＝780千米/小时、v2＝1000千米/小时，地球表面二维平面网格的最小边长a＝210米、b＝270米；步骤2包括：步骤2
‑
1：进行第1层级网格剖分与编码；步骤2
‑
2：进行第2～6层级网格剖分与编码；步骤2
‑
3：进行第7层级剖分与编码；步骤2
‑
4：将航空器所在经纬度坐标转化为网格编码；
步骤3包括：步骤3
‑
1：第1层级高度范围编码；步骤3
‑
2：第2层级高度范围编码；步骤3
‑
3：第3层级高度范围编码；步骤3
‑
4：进行第4到第7层级高度范围编码；步骤3
‑
5：确定出航空器的飞行高度所在的第7层级高度范围。4.如权利要求3所述的一种航空器运行可视化方法，其特征在于：步骤2
‑
4中将航空器所在经纬度坐标转化为网格编码包括：步骤2
‑4‑
1：计算第1层级网格编码；步骤2
‑4‑
2：计算第2层级网格编码；步骤2
‑4‑
3：计算第3层级网格编码；步骤2
‑4‑
4：计算出第4层级网格的编码；步骤2
‑4‑
5：计算出第5层级网格的编码；步骤2
‑4‑
6：计算出第6层级网格的编码；步骤2
‑4‑
7：计算出第7层级网格的编码。5.如权利要求3所述的一种航空器运行可视化方法，其特征在于：步骤4
‑
6中设某个第7层级网格的行列坐标为(x，y)，然后根据行列坐标(x，y)计算得到网格第level_n 1～7级编码计算出表征各航段的第7层级二维平面网格编码的集合包括：步骤4
‑6‑
1：计算坐标(x，y)对应的网格第7级编码；步骤4
‑6‑
2：计算坐标(x，y)对应的网格第6级编码；步骤4
‑6‑
3：计算坐标(x，y)对应的网格第5级编码；步骤4
‑6‑
4：计算坐标(x，y)对应的网格第4级编码；步骤4
‑6‑
5：计算坐标(x，y)对应的网格第3级编码；步骤4
‑6‑
6：计算坐标(x，y)对应的网格第2级编码。6.如权利要求3所述的一种航空器运行可视化方法，其特征在于：步骤2
‑
1中选择本初子午线与赤道交点为剖分原点，用4
°×4°
的第1层级网格将纬度范围为[
‑
88
°
,88
°
]、经度范围为[
‑
180
°
,180
°
]的地球表面剖分为44
×
90份；步骤2
‑
2中层级数小的为上级，层级数大的为下级，上一层级网格为高层级网格，作为父网格，下一层级的网格为低层级网格，作为子网格，每1个上一层级父网格平均分成16份，得到下一层级16个子网格，对16个子网格按照皮亚诺空间填充曲线进行16进制编码，直到剖分至第6层级；步骤2
‑
3中将第6层级网格平均分成4份，得到第7层级网格。7.如权利要求6所述的一种航空器运行可视化方法，其特征在于：横坐标为数字自西向东从00～89，间隔为1，共90个数字；纵坐标中北纬从低到高为字母a～v，共22个大写英文字母；南纬从低到高为字母a～v，共22个小写英文字母；第1层级网格序号在纬度方向上用a～v和a～v共44个英文字母从赤道按照纬度由低到高编码，其中北纬大写、南纬小写；第1层级网格序号在经度方向上从本初子午线开始由西向东用00～89共90个数字编码。8.如权利要求7所述的一种航空器运行可视化方法，其特征在于：
步骤2
‑
4中设航空器所在的经纬度坐标为(l，b)，其中l为纬度，b为经度，l和b的表示形式均为度分秒，记作l＝l
d
°
l
m
′
l
s
″
,b＝b
d
°
b
m
′
b
s
″
，其中l
d
为纬度值的度数，l
m
为纬度值的分数，l
s
为纬度值的秒数；b
d
为经度值的度数，b
m
为经度值的分数，b
s
为经度值的秒数；步骤2
‑4‑
1中计算第1层级网格编码，计算公式为：l
d
/4
°
的商为m1、余数为n1，m1 1对应的英文字母即为该点所在的第1层级网格的纬度方向编码，北纬大写，南纬小写；b
d
/4
°
商为m1、余数为n1，m1即为该点所在的第1层级网格的经度方向编码；步骤2
‑4‑
2中计算第2层级网格编码，计算公式为：n1为步骤2
‑4‑
1的余数，m2即为该点所在的第2层级网格的纬度方向编号；n1为步骤2
‑4‑
1的余数，m2即为该点所在的第2层级网格的经度方向编号；步骤2
‑4‑
3中计算第3层级网格编码，计算公式为：l
m
/15
′
的商为m3、余数为n3，m3即为该点所在的第3层级网格的纬度方向编号；b
m
/15
′
商为m3、余数为n3，m3即为该点所在的第3层级网格的经度方向编号；步骤2
‑4‑
4，计算出第4层级网格的编码，计算公式为：n3为步骤2
‑4‑
3的余数，故n3/3.75
′
的商为m4、余数为n4，m4即为该点所在的第4层级网格的纬度方向编号；n3为步骤2
‑4‑
3的余数，n3/3.75
′
商为m4、余数为n4，m4即为该点所在的第4层级网格的经度方向编号；步骤2
‑4‑
5中计算出第5层级网格的编码，计算公式为：步骤2
‑4‑
4的余数为n4，故(n4×
60 l
s
)/56
″
的商为m5，余数为n5，m5即为该点所在的第5层级网格的纬度方向编号；(n4×
60 b
s
)/56
″
商为m5、余数为n5，m5即为该点所在的第5层级网格的经度方向编号；
步骤2
‑4‑
6中计算出第6层级网格的编码，计算公式为：n5/14
″
的商为m6、余数为n6，m6即为该点所在的第6层级网格的纬度方向编号；n5/14
″
商为m6、余数为n6，m6即为该点所在的第6层级网格的经度方向编号；步骤2
‑4‑
7中计算出第7层级网格的编码，计算公式为：n6/7
″
的商为m7，余数为n7，m7即为该点所在的第7层级网格的纬度方向编号；n6/7
″
商为m7，余数为n7，m7即为该点所在的第7层级网格的经度方向编号。9.如权利要求3所述的一种航空器运行可视化方法，其特征在于：步骤3
‑
1中将整个0～40000米的高度范围进行二等分，高度在[0米，20000米)之内编码为0，高度在[20000米，40000米]之内编码为1；步骤3
‑
2对第1层级上下2个高度范围编码均进行二等分，第2层级的高度差为10000米；将[0米，20000米)的高度范围进行二等分，高度在[0米，10000米)之内编码为00，高度在[10000米，20000米)之内编码为01；同理，将[20000米，40000米]的高度范围进行二等分，高度在[20000米，30000米)之内编码为10，高度在[30000米，40000米]之内编码为11；步骤3
‑
4中第4到第7层级高度差分别为2500米、1250米、625米、312.5米；层级数小的为上级，层级数大的为下级，下一层级在上一层级各个高度范围内均进行二等分，其中，二等分之后高度范围较低的在上一层级高度范围编码后面加0，高度范围较高的在上一层级高度范围编码后面加1；上下级高度范围编码之间具有继承性，编码的长度代表该高度范围所在的层级；步骤3
‑
5中由高度值h计算出该高度值所在的第7层级高度范围编码的方法为：h除以312.5，若能整除，所得结果为t；若不能整除，所得结果的整数部分为t，小数部分为t；将十进制数值t转化为7位二进制编码，即为高度值h所在的第7层级高度范围编码。10.如权利要求2所述的一种航空器运行可视化方法，其特征在于：步骤6中r值对应小时，g值对应分钟，b值对应秒；r值按照步长为8从0开始取30个值，最小为0，最大为232；g值和b值按照步长为4取60个值，最小为0，最大为236。
技术总结
本发明提供一种基于四维时空网格的航空器运行可视化方法，包括：确定地球表面二维平面网格的最小边长，建立面向航空器运行的地球表面剖分网格体系，对网格高度进行剖分并编码，对航空器运行航线进行二维平面网格化表征，构建面向航空器运行的四维时空网格，对四维时空网格着色；本发明方法可以在显示屏上直观化、一体化的显示航空器四维时空信息，实现航空器飞行过程中四维状态可视化，方便管制员及时了解航空器的位置和运行时刻，能够为管制员提供直观、动态、全面的航空器运行状态显示，使得管制员对航空器未来一段时间的飞行动态有一个全局的认识，更好地保证飞行安全。更好地保证飞行安全。更好地保证飞行安全。


技术研发人员：万路军 徐鑫宇 黄阿倩 蔡明 高志周 戴江斌 沈堤 陈致远 余付平 霍丹
受保护的技术使用者：中国人民解放军空军工程大学
技术研发日：2021.04.02
技术公布日：2021/6/29