探空指数产品-python实现

最近有任务需要计算探空指数产品，在经历了漫长扒文献、找公式实现公式的过程后，终于是凑齐了所需的指数产品。总结一下这一段的历程，同时提供给可能需要的人，提高效率减少大海中找针的工作量。

分类

本篇文章中涉及到的指数产品分为三类：一、强对流天气指数；二、高空特征层；三、综合观测探空。

强对流天气指数

强对流天气指数共有10种，分别是：

1、对流有效位能（CAPE）

2、A指数

3、K指数

4、沙氏指数（SI）

5、强天气威胁指数（SWEAT）

6、风暴强度指数（SSI）

7、抬升指数（LI）

8、总指数（TT）

9、垂直风切边指数（VWS）

10、山崎指数（KYI）

对流有效位能CAPE

环境正浮力对绝热上升到自由对流高度后，处于自由对流高度和平衡高度之间的气块累计做的功，反映了不稳定能量以及对流发展的强弱程度。它是一个能定量反映大气环境中是否可能发生深厚对流的热力变量，反映了从低层到高层的能量变化，表示从自由对流高度到平衡高度对流层的累积浮力能。

表示干空气比气体常数，，、分别表示气块虚温和环境虚温，具体实现可以参考这篇文中的公式关于虚温订正CAPE求算方法的讨论，、分别表示平衡高度气压和自由对流高度气压。

介绍一个很强大的用于处理气象数据的工具metpy，这个工具封装了大量气象业务用到的公式，可以避免重复造轮子，提高开发效率。

metpy中内置了cape计算公式，通过metpy.calc模块中的cape_cin()方法计算得到。需要注意的是通过metpy计算得到的物理量是一个units对象，通过magnitude方法可以得到对应的数值。

A指数

Ａ指数不仅包括各层之间的温度递减率而且包括低层、中低层和中层的温湿度条件。Ａ指数愈大，大气愈不稳定也愈潮湿。可反映整层大气层结不稳定程度。

不知道为什么公式编辑器崩溃了！！！

A = T850 + T500 - (Td850+Td700+Td500)

A指数等于850hpa和500hpa处的温度减去850hpa，700hpa和500hpa处的露点温度。

K指数

K指数计算式中第一项表示温度直减率，第二项表示低层水汽条件，第三项表示中层饱和程度，所以K指数能够反映大气的层结稳定情况，K指数越大，层结越不稳定。

沙氏指数 

沙氏指数指850hPa 等压面上的湿空气块沿干绝热上升到抬升凝结高度（LCL），再沿湿绝热线上升到500hPa。在500hPa 上的大气实际温度（ T500）与该上升气块达到500hPa 时的温度（ T′）之差。

SI=T500-T′

def showalter(self,t8,td8,t5):
        """
        calculating Shoalter Index
        Input:
            t8: temperature at 850hpa [degree C]
            td8: dew point temperature at 850hpa [degree C]
            t5: temperature at 500hpa [degree C]
        Output:
            showalter index [deg C]
        """
        P5 = 500
        P8 = 850
        Cpd = 0.2403
        Cpv = 0.445
        Rd = 6.85578*0.01
        Rw = 11.017874*0.01
        T0 = 273.15
        L0 = 597.40
        C = 1.002
        C1 = 0.57
        step = 10
        T8 = T0 + t8
        Etd8 = self.E_water(t8)
        w = (Rd/Rw)*Etd8/(P8-Etd8)
        ml = (Cpd*(1+Cpv*w/Cpd))/(Rd*(1+w/(Rd/Rw)))
        Ta = self.Tc(P8,t8,td8) + T0
        Pa = P8*(pow(Ta/T8,ml))
        m2 = (Cpd/Rd)*(1+C*w/Cpd)
        ETa = self.E_water(Ta-T0)
        Qc = math.log((Pa-ETa)/pow(Ta,m2))-(0.622/Rd)*(((L0+C1*(T0-Ta))/Ta)*(ETa/(Pa-ETa)))
        out = T8
        Eout = self.E_water(out-T0)
        Q5 = math.log((P5-Eout)/pow(out, m2))-(0.622/Rd) * (((L0+C1*(T0-out))/out)*(Eout/(P5-Eout)))
        while (abs(Qc-Q5) > 0.0001):
            if (Qc > Q5):
                out -= step 
            else:
                out += step
                step /= 5
                out -= step 

            Eout = self.E_water(out-T0)
            Q5 = math.log((P5-Eout)/pow(out, m2))-(0.622/Rd) * (((L0+C1*(T0-out))/out)*(Eout/(P5-Eout)))

        return t5-(out-T0)

强天气威胁指数

强天气威胁指数（SWEAT）是根据20世纪70年代Miller和Maddox引入的，是根据328次龙卷风资料和日常预报经验得出的一个预报指数。

其中I代表SWEAT；为850hpa露点温度，如果是负数，此项为0；TT =(T+Td)-2*T500表示总指数，如果TT小于49，则20*(TT-49)项等于0；f850为850hpa风速；S=sin(a500-a850),a500和a850分别代表500hpa风向和850hpa风向；最后一项125*(S+0.2)在下列4个条件中任何一条件不具备时为零；850hpa风向在135度～250度之间；500hpa风向在210度～310度之间；500hpa风向减850hpa风向为正；850hpa及500hpa风速至少等于7.5m/s。

风暴强度指数

 一种可用于区分强雷暴与非强雷暴的经验物理参数，表示浮力能与垂直风切变之间的线性关系。

其中SHR为平均垂直风切变，CAPE为对流有效位能，当CAPE越强时，SSI越大，产生强风暴天气的可能性以及风暴强度就越大（大气热力不稳定条件贡献更大）；当SHR越大时，对CAPE数值的削减作用就越小，使得SSI变大（强垂直风切变的动力作用贡献更大）。

抬升指数

指一个气块从自由对流高度出发沿湿绝热线上升到500hPa处所示的温度（T）与500hPa等压面上的环境温度（T′）之差。

总指数
表示垂直温度直减率，TT越大越容易发生对流天气。

<44，无雷暴活动； 
45 ~ 50 ，有雷暴活动；
51~55，强雷暴可能发生
>55 ，强雷暴极有可能发生，或可能出现龙卷；

垂直风切变指数（VWS）

环境风的垂直切变对风暴单体的组织、结构及演变特征影响极大，是对流系统发展的转换条件。

U200,U850,V200,V850分别表示200hpa和850hpa处的东向风和纬向风。

山崎指数

该指数根据对流发生的三个条件：大气稳定度、低层水汽和上升运动归纳而成，用于预报强降水的发生，预报时效为12小时。

TA <=SI; KYI = 0

TA > SI,

其中TA为850hpa-500hap的温度平流，SI为沙氏指数，T-Td为850hpa的温度露点差。

下篇文章将会介绍高空特征层的提取方式