国际民航组织（ICAO）有几个出版物是数学家写的，其中DOC8168是讲如何用数学防止撞山撞障碍物，DOC9689则是讲如何用数学防止航空器空中相撞和危险接近。飞行离不开数学，从达芬奇开始。而当今时代，数学离不开计算机，上述出版物里面的所有公式和算法，都可以经由计算机打包成函数库，从而极大的简化阅读，投入实战。
中国民航的AC-91-27标题是《飞行程序》，来自DOC8168第一卷，基本目标应该是飞行员使用、签派员熟悉、管制员了解。AC-97-5《飞行程序设计规范》则来自DOC8168第二卷，供程序设计人员使用。
本篇结合计算机编程讨论上述两个AC都提到的风螺旋线，以使飞行员进一步了解数学如何帮助我们避障，减少航图中的疑惑，并鼓励飞行员广泛参与到程序设计中来。
传统的程序设计人员使用AUTO CAD来分析飞行程序并画图，底层逻辑不可控，今天我们改用Python徒手画，并利用Folium直接呈现在地形图上，地理数据采用天地图，且Python、Folium、天地图均可合法免费获取。
什么是风螺旋线
先提个问题：如上图，飞机飞向某个点，然后在这个点做一个盘旋，盘旋改出的时候会不会有地形问题？
根据速度和坡度，可以计算出盘旋的半径：
r = TAS / (20 * np.pi * R) # 无风半径 公式I-2-3-3
从而确定基本航迹，看上去没问题，但是风却是不确定因素，最差的情况是刮东风，导致飞机改出位置在红三角处，很可能受地形影响。
因此，能不能在这里盘旋，需要考虑最差情况，即此处的东风极值，风速×盘旋时间，得到出圈距离，然后比较障碍物或地形因素。
如果只转半圈，则需要考虑西风极值，以及盘旋半圈所需时间，得到出圈距离；如果只转90度，考虑南风，以及四分之一圈所需时间。
所有方向，考虑最不利风，最终出圈距离连成线，就是风螺旋线，如下图（蓝色、黄色三角分别是180度、90度对应的出圈位置）：
再来看8168描述的风螺旋线图：
实现代码：
from geopy.distance import geodesic
import geopy
import numpy as np
def getLatLon(f1, f2, center, hdg): # 根据位置点对中心点的偏移位置计算经纬度
d1 = np.arctan2(f1, f2) / np.pi * 180
d = hdg + d1 # 实际方位
e = (f1 ** 2 + f2 ** 2) ** 0.5 / 1000 # 距离
p = geopy.distance.distance(kilometers=e).destination(center, bearing=d) # 求坐标
return [p.latitude, p.longitude]
def spiral(r,RWhdg,p,R,w,Edir,delta,lr): # r值，跑道真航向，圆心，R值，风速，E偏流增减方向，E偏流角，左右旋
Spiral=[]
start=0
for bearing in np.arange(start,start+lr*(Turn+135), lr*1): # 设置角度范围
pos = geopy.distance.distance(kilometers=r).destination(p, bearing=bearing + RWhdg -lr*180 ) # 求盘旋轨迹坐标
x0 = pos.latitude
y0 = pos.longitude
E = abs((bearing-start) / R * w / 3600) # 计算Etheta 来自公式I-2-3-4
pos1 = geopy.distance.distance(kilometers=E).destination((x0, y0),
bearing=bearing + RWhdg -lr*180 - lr*Edir * delta) # 求风螺旋轨迹坐标
x1 = pos1.latitude
y1 = pos1.longitude
Spiral.append([x1, y1]) # 风螺旋线
return Spiral
if __name__ == '__main__':
# 基本参数
lat=39.3382
lon=112.3105
center = (lat, lon) # 绘图中心
waypoint=center # 航路点，也是螺旋中心点
RWhdg =97 # 跑道方向，或下一段航迹方向
magvar=6.8 # 磁差
RWhdg = RWhdg - magvar # 消除磁差
IAS = 380 # km/h 最大表速 来自表I-4-1-1，查航图速度限制
Alt = 2400 # 关键点高（米） 来自航图
ISAd = 15 # 与ISA的温度差
bank = 25 # 坡度，来自表I-2-3-1
w = 30*1.852 # 风速，来自表I-1-2-3-1
Turn = 90 # 转弯角度
lr=1 # 螺旋方向 左旋-1，右旋 1
# 初步计算
bank = bank / 180 * np.pi # 转为弧度
TAS = IAS * 171233 * ((288 + ISAd) - 0.006496 * Alt) ** 0.5 / (288 - 0.006496 * Alt) ** 2.628 # 公式I-2-1 附录1
R = 6355 * np.tan(bank) / (np.pi * TAS) # 转弯率 公式I-2-3-1
if R > 3:
R = 3
r = TAS / (20 * np.pi * R) # 无风半径 公式I-2-3-3
C = (TAS + 56) * 6 / 3600 # 飞行员延迟距离 来自表I-2-3-2
v = TAS
delta = np.arcsin(w / v) / np.pi * 180 # w/v对应的角度 来自图I-2-3-4
E0 = Turn / R * w / 3600 # 计算E0 对应表I-2-3-2的E值
print('IAS=', IAS, 'TAS=', TAS)
print('R=', R)
print('r=', r)
print('C=', C)
print('E0=', E0)
Edir=lr # E偏流方向
TT=180-Turn #
po = getLatLon(-r * 1000, 0, center, RWhdg + TT) # 确定转弯圆心，假设为过点转弯
# 风螺旋线坐标数据
spiraline=spiral(r,RWhdg+TT-90,po,R,w,Edir,delta,lr)
print(spiraline)
这段代码直接用地球经纬度计算，结果可以直接输出到各类地图；打印参数用于与8168表格核对，以验证计算准确性；两个函数则可以直接打包引用。
小结：风螺旋线，就是盘旋或转弯过程中，在最不利风的影响下，飞机的出圈位置，转的越久，距离越远。用这条线测算障碍物的影响，比用盘旋半径做的轨迹靠谱的多。
风速打哪儿来
一年365天，天天风不同，风向可以取0-360°，那么风速怎么办？
8168给出了几种方案：
指定风：相对机场高度较低的地点，使用指定风，由8168的表I-2-3-1直接指定，40kt、30kt、25kt，随高度降低。不考虑风向，均按照最不利风向。本案例为指定风30kt。
小结：懒人可以直接用标准风和指定风；想获取实测风数据的，可以找航科院，那里有10年以上的QAR数据，覆盖所有的航路点、转弯点，挖掘即得。
转弯的膝盖
假设要建一座新机场，这是它的进场程序，几个航路点依次是IAF、IF、FAF、MAPt。我们知道飞机速度很大，不可能沿灰色线直角转弯；如果到了IF点（实心标记）转弯，就需要转半个圆才能回到下一段航路上，这称为过点转弯。通常我们不这么飞。
如果提前一个盘旋半径的距离转弯，转90°即可回到下一个航段，最经济最安全，这称为旁切转弯。这是最常见的转弯方式。
但是研究真实飞行轨迹的话，还要考虑两个因素：定位容差和飞行员误差。
定位容差意思是任何时候飞机的定位没那么准，有时候是卫星信号问题，你拿着手机不动，获得的坐标可能时不时小范围跳变；有时候是飞机设备问题，也有可能是坐标本身没那么精确，所以8168给出的方案是：
TSE=1.852 # RNP-1，1NM
ATT=0.8*TSE # 见2.2.1
ATT表示定位容差，也就是转弯起点有可能前后差0.8海里。
飞行员误差是说飞行员有最多6秒延迟转弯的可能：
C = (TAS + 56) * 6 / 3600 # 飞行员延迟距离 来自表I-2-3-2
两个因素叠加，最早的转弯起点比上图起点早ATT那么远（蓝线），最晚的转弯起点则比上图起点晚ATT+C（红线）。
上图红线表示最晚转弯的情况，蓝线表示最早转弯的情况。现在可以把风的因素考虑进去，得出最晚转弯螺旋线：
这时蓝线和红线之间的部分，是飞机有可能飞过的区域，且呈高斯分布状态，中央多，边缘少；之外的区域则属于飞机基本不可能去的地方。
小结：飞机因为速度大，不可能转直角弯；旁切转弯是最佳方案；由于定位容差和飞行员误差的原因，实际转弯轨迹会在某个区间，这个区间可以被计算出来，其形状有点像人的膝盖。
如何用风螺旋线避障
先回顾一下超障的概念。关于MOC的概念，在“低温修正批判”有介绍，记不清的可以去翻翻。上图则是截面图，可以看到整个航路宽度中，中间的50%是主区，障碍物高度和飞行高度间隔是MOC，例如IAF段MOC=300米；副区则是两侧各25%，间隔高度呈线性缩减直至0。
本篇的案例，假定是RNP1进场，RNP0.3进近，那么起始段总宽是5nm，主区宽度是2.5nm，如果用灰色表示副区，整体画出来是上面这个样子，夹在灰区中间的部分是主区。
显然，主区副区都严重受地形影响，需要换地方。
地方政府给划了新的地块建设机场，现在看主区副区都可以接受。
现在处理膝盖问题。为了便于理解，先说过点转弯：前面讲过，转弯有最早和最晚的情况，上图黄色区域的上下边界即表示定位容差；转弯还有左右误差的可能性，所以要研究沿主区边界转弯的极端情况，上图黄区的左右边界即表示主区宽度；那么从黄区的四个角分别画转弯轨迹，可以得到扩展的主区范围。此处省略，很少有人这么飞。
常用的旁切转弯，考虑到定位容差、飞行员误差，根据前面的计算，黄格转为粉红格区域，从该区域四个角分别为起点画转弯轨迹，可以得到扩展的主区范围：
四条轨迹中，左下角的转弯轨迹完全在原主区，因此忽略；其余三条（两红一蓝）转弯轨迹分别侵入了主区之外，因此有必要重新调整主区，确保对正五边之前三条轨迹不会侵入灰区；主区调整之后，再调整相应副区：
基于风螺旋线的转弯保护区就此成型。现在，由于地形的影响肉眼可见，这个程序也不可用，需要继续找新地块。
跑道位置和方向调整之后，地形对主区副区来讲都满足要求。最后，还要分析人工建筑物：
总结：飞机进场到着陆期间的程序，必须严格评估障碍物，8168提供了整套数学方法，既保证飞行中不会碰撞障碍物，也不会过度限制空域，导致经济建设受阻。其中风螺旋线是比较难啃的一块，搞明白之后其余数学公式都可以轻松理解运用。
后记
上面构建的函数和图形，仍然需要人工去判别障碍物是否超限，有没有更好的办法，让计算机自动判别呢？那必须。
地形数据，介绍下SRTM：本世纪初，奋进号航天飞机对地球表面做了全面扫描，获取了南北纬60度之间全部地形的雷达影像数据，后经处理成为DEM，即全球数字高程模型，面向公众免费下载。国内可在国家基础学科公共科学数据中心网站下载。其高程误差最大也就20米，完全满足地形分析的需求。
人工建筑物信息，可以找当地政府获取。
障碍物信息收集齐全后，可以用以下函数判定是否超限：
def fence(point, fence): # 检测位置点是否出圈，不支持异构
c = point # 障碍物坐标
f = fence.copy() # 围界坐标列表
eage = len(f)
f.append(fence[0]) # 闭合圈
yes = True
for k in range(eage): # 遍历每个边
a = f[k] # 边线起点
b = f[k + 1] # 边线终点
hdg1 = geoDegree(a, b, 0) # 计算每个边线的航迹角
d = [(a[0] + b[0]) / 2, (a[1] + b[1]) / 2] # 边线中点坐标
hdg2 = geoDegree(c, d, 0) # 计算障碍物和边线中点的航迹角
diff = hdg2 - hdg1 # 航迹差
if diff < 0:
diff += 360
if 180 >= diff >= 0: # 任何一边突破，认为在区域外，结束判断
yes = False
break
return yes
深入了解数学防撞，请进入下列链接：
