Python气象数据处理与绘图(5)：气候突变检验(年代际突变检验)

长时间尺度的气候研究中通常涉及到年际变化和年代际变化。文献中也经常出现所谓年代际突变这类的描述，这次便介绍两个用来检验年代际变化的方法，其一是滑动t检验(sliding t-test)，另一个则是曼肯德尔检验(Man-Kendall test)。

滑动t检验 (sliding t-test)

滑动t检验是考察两组样本平均值的差异是否显著来检验突变。
滑动t检验的基本思想是：把一气候序列中两段子序列均值有无显著差异看为来自两个总体均值有无显著差异的问题来检验。如果两段子序列的均值差异超过了一定的显著性水平，可以认为均值发生了质变，有突变发生。

参考魏凤英老师现代气候统计诊断与预测技术一书

原文也提到了该方法的局限性，需要人为设置滑动步长，具有一定主观性，需反复设置不同步长最终确定合适的突变点。

滑动t检验计算步骤

根据以上步骤，可通过如下代码块实现：

def slidet(inputdata,step):
    inputdata = np.array(inputdata)
    n = inputdata.shape[0]
    t = np.zeros(n)
    t1 = np.empty(n)
    n1 = step    #n1, n2为滑动补偿，需调整
    n2 = step
    n11 = 1 / n1
    n22 = 1 / n2
    m = np.sqrt(n11 + n22)
    for i in range (step, n-step-1):
        x1_mean = np.mean(inputdata[i-step : i])   
        x2_mean = np.mean(inputdata[i : i+step])
        s1 = np.var(inputdata[i-step : i])          
        s2 = np.var(inputdata[i : i+step])
        s = np.sqrt((n1 * s1 + n2 * s2) / (n1 + n2 - 2))
        t[i-step] = (x2_mean - x1_mean) / (s * m)
        t1 = np.roll(t , step-1)
        t1[:step]=np.nan  
        t1[n-step+1:]=np.nan             
    return  t1

自由度n1+n2-2，根据置信度检验表查找对应显著性阈值。

曼肯德尔检验(Man-Kendall test)

曼肯检验是一种非参检验方法，避免了滑动t检验的局限性，在年代际变化研究中受到广泛应用。

参考魏凤英老师现代气候统计诊断与预测技术一书

以下是计算步骤：

参考魏凤英老师现代气候统计诊断与预测技术一书

根据以上步骤，可通过如下代码块实现：

def mktest(inputdata):
    inputdata = np.array(inputdata)
    n=inputdata.shape[0]
    Sk             = [0]
    UFk            = [0]
    s              =  0
    Exp_value      = [0]
    Var_value      = [0]
    for i in range(1,n):
        for j in range(i):
            if inputdata[i] > inputdata[j]:
                s = s+1
            else:
                s = s+0
        Sk.append(s)
        Exp_value.append((i+1)*(i+2)/4 )                    
        Var_value.append((i+1)*i*(2*(i+1)+5)/72 )            
        UFk.append((Sk[i]-Exp_value[i])/np.sqrt(Var_value[i]))
    Sk2             = [0]
    UBk             = [0]
    UBk2            = [0]
    s2              =  0
    Exp_value2      = [0]
    Var_value2      = [0]
    inputdataT = list(reversed(inputdata))
    for i in range(1,n):
        for j in range(i):
            if inputdataT[i] > inputdataT[j]:
                s2 = s2+1
            else:
                s2 = s2+0
        Sk2.append(s2)
        Exp_value2.append((i+1)*(i+2)/4 )                    
        Var_value2.append((i+1)*i*(2*(i+1)+5)/72 )            
        UBk.append((Sk2[i]-Exp_value2[i])/np.sqrt(Var_value2[i]))
        UBk2.append(-UBk[i])
    UBkT = list(reversed(UBk2))
    return UFk, UBkT

示例

对于同一组数据：

a = [15.4,14.6,15.8,14.8,15.0,15.1,15.1,15.0,15.2,15.4,
      14.8,15.0,15.1,14.7,16.0,15.7,15.4,14.5,15.1,15.3,
      15.5,15.1,15.6,15.1,15.1,14.9,15.5,15.3,15.3,15.4,
      15.7,15.2,15.5,15.5,15.6,15.1,15.1,16.0,16.0,16.8,
      16.2,16.2,16.0,15.6,15.9,16.2,16.7,15.8,16.2,15.9,
      15.8,15.5,15.9,16.8,15.5,15.8,15.0,14.9,15.3,16.0,
      16.1,16.5,15.5,15.6,16.1,15.6,16.0,15.4,15.5,15.2,
      15.4,15.6,15.1,15.8,15.5,16.0,15.2,15.8,16.2,16.2,
      15.2,15.7,16.0,16.0,15.7,15.9,15.7,16.7,15.3,16.1]
plt.figure(figsize=(10,5))
plt.plot(a,'r')
plt.show()

分别使用两种方法进行检验，首先是滑动t检验：

t = slidet(a,5)
plt.figure(figsize=(10,5))
plt.plot(t,'r')
plt.axhline(1.8595)
plt.axhline(-1.8595)
plt.show()

sliding t-test

存在多个突变点，这时便需要调整滑动补偿，选取合适的步长。
而利用MK检验：

uf,ub = mktest(a)
plt.figure(figsize=(10,5))
plt.plot(uf,'r',label='UFk')
plt.plot(ub,'b',label='UBk')
plt.legend()
plt.axhline(1.96)
plt.axhline(-1.96)
plt.show()

M-K test

对于该组数据，相比之下，MK检验的效果要优于滑动t检验。
两种检验的分析方法如下：
滑动t检验:

滑动t检验

M-K检验：

MK检验

最后编辑于：2020.02.05 17:58:43