作者：CSDN @ _养乐多_

本文将介绍标准正态同质性检验（Standard Normal Homogeneity Test，SNHT) 突变点检测代码。以 NDVI 时间序列为例。输入数据可以是csv，一列NDVI值，一列时间。代码可以扩展到遥感时间序列突变检测（突变年份、突变幅度等）中。

结果如下图所示，

一、准备数据

测试数据（0积分下载）：https://download.csdn.net/download/qq_35591253/88895803

该数据是GEE上提取的，参考博客《GEE：基于Landsat5/7/8/9数据提取一个点的NDVI时间序列（1986-2024）》

二、Pettitt介绍和代码

2.1 原理和步骤

Pettitt方法是一种非参数检验方法，用于检测时间序列数据中的突变点。这种方法的核心思想是比较序列中每个数据点与其前面所有数据点的值的大小关系，从而判断是否存在突变点。

具体来说，Pettitt方法通过计算一个统计量Ut来检测突变点。这个统计量是基于序列中每个数据点的值与序列平均值的差异累积得到的。如果在某个位置，统计量Ut显著地高于其他位置，那么这个位置就可能是一个突变点。

在实际应用中，通常会结合其他方法（如滑动T检验、Mann-Kendall趋势检验等）来进行综合判断，以提高突变点检测的准确性和可靠性。

2.2 核心函数

def Pettitt_change_point_detection(X):
    inputdata = np.array(X)
    n = inputdata.shape[0]
    U = [];
    s = 0
    U.append(sum(np.sign(X[0] - np.array(X[:]))))
    for t in range(1, len(X)):
        v_t = sum(np.sign(X[t] - np.array(X[:])))
        U.append(U[t - 1] + v_t)
    Kt = max(np.abs(U))
    max_idx = list(np.abs(U)).index(max(np.abs(U)))
    pvalue = 2 * np.exp((-6 * (U ** 2)) / (n ** 3 + n ** 2))
    if pvalue <= 0.05:
        change_point_desc = 0#'显著'
    else:
        change_point_desc = 1#'不显著'
    return max_idx#,change_point_desc

三、读取csv格式时序数据的示例

示例代码以csv格式数据为例子，当然可以扩展到遥感时间序列中，遥感时序数据的分析代码框架参考博客《python：处理遥感时间序列（代码框架），并保存结果》

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 从CSV文件中读取数据
df = pd.read_csv('测试数据\\eendvi-chart.csv')
df = pd.DataFrame(df).dropna()

df["system:time_start"] = pd.to_datetime(df["system:time_start"], format='mixed')
df['time'] = df['system:time_start']
seasonal_ts = df
nSTS = len(seasonal_ts)
# 取序列的前80%
percentile_start = 0 # 74
percentile_end = 100   # 86
percentile_index_start = int(nSTS * percentile_start / 100)
percentile_index_end = int(nSTS * percentile_end / 100)
seasonal_ts1 = seasonal_ts[percentile_index_start : percentile_index_end]
seasonal_ts = seasonal_ts1['NDVI'].values
time = seasonal_ts1['time'].values

def Pettitt_change_point_detection(X):
    inputdata = np.array(X)
    n         = inputdata.shape[0]
    U = [];
    s = 0
    U.append(sum(np.sign(X[0] - np.array(X[:]))))
    for t in range(1, len(X)):
        v_t = sum(np.sign(X[t] - np.array(X[:])))
        U.append(U[t - 1] + v_t)
    Kt = max(np.abs(U))
    max_idx = list(np.abs(U)).index(max(np.abs(U)))
    pvalue = 2 * np.exp((-6 * (U ** 2)) / (n ** 3 + n ** 2))
    if pvalue <= 0.05:
        change_point_desc = 0#'显著'
    else:
        change_point_desc = 1#'不显著'
    return max_idx#,change_point_desc

# 检测变点
change_point = Pettitt_change_point_detection(seasonal_ts)

# 绘制图形
plt.figure()
plt.scatter(time, seasonal_ts, color='black', label='NDVI', s=1.5)
plt.axvline(x = seasonal_ts1['time'].values[change_point], color='r', linestyle='--', label='Change Point')

plt.legend()
plt.show()

声明：
本人作为一名作者，非常重视自己的作品和知识产权。在此声明，本人的所有原创文章均受版权法保护，未经本人授权，任何人不得擅自公开发布。
本人的文章已经在一些知名平台进行了付费发布，希望各位读者能够尊重知识产权，不要进行侵权行为。任何未经本人授权而将付费文章免费或者付费（包含商用）发布在互联网上的行为，都将视为侵犯本人的版权，本人保留追究法律责任的权利。
谢谢各位读者对本人文章的关注和支持！