作者：CSDN @ _养乐多_

本文将介绍布伊山德U检验（Buishand U test，BUT）突变点检测代码。以 NDVI 时间序列为例。输入数据可以是csv，一列NDVI值，一列时间。代码可以扩展到遥感时间序列突变检测（突变年份、突变幅度等）中。

结果如下图所示，

一、准备数据

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

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

二、BUT介绍和代码

Buishand U test突变点检测中文名为布伊山德U检验，而其原理是基于正态分布变量的单一变点模型。

2.1 原理和步骤

Buishand U test是用于检测时间序列数据中是否存在一个突变点，即某个时刻数据的特性发生了显著的变化。这种检测方法在气象、水文、生态和经济学等领域具有广泛的应用。

Buishand U test的原理是考虑一个正态随机变量X，并假设存在一个单变点将数据集分为两部分。在变点之前，数据遵循均值为μ和方差σ²的正态分布；在变点之后，数据的均值变为μ+δ。测试的零假设H₀是δ=0，即没有变点；备择假设H₁是δ≠0，即存在变点。为了进行检测，计算调整后的局部和统计量Sₖ，使用样本标准偏差D(x)，然后构造U统计量：

$$U=\frac{1}{n\ast (n+1)}\ast \sum _{k=1}^{n-1}(S[k]-Dx{)}^2$$

其中，n是样本大小，Sₖ是到第k个观测值为止的调整后的局部和，D(x)是样本的标准偏差。变点的位置K通过最大化|Sₖ|来确定。

2.2 核心函数

def Buishand_U_change_point_detection(input_data):
    input_data = np.array(input_data)
    mean_value = np.mean(input_data)
    data_length = input_data.shape[0]
    index_range = range(data_length)
    sum_diff = [np.sum(input_data[0:x+1] - mean_value) for x in index_range]
    standard_deviation = np.sqrt(np.sum((input_data-np.mean(input_data))**2)/(data_length-1))
    U_statistic = np.sum((sum_diff[0:(data_length - 2)]/standard_deviation)**2)/(data_length * (data_length + 1))
    absolute_sum_diff = np.abs(sum_diff)
    max_absolute_sum_diff = np.max(absolute_sum_diff)
    change_point_index = list(absolute_sum_diff).index(max_absolute_sum_diff) + 1
    normalized_sum_diff = (sum_diff/standard_deviation)
    return change_point_index

三、读取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('测试数据\\ee-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 Buishand_U_change_point_detection(input_data):
    input_data = np.array(input_data)
    mean_value = np.mean(input_data)
    data_length = input_data.shape[0]
    index_range = range(data_length)
    sum_diff = [np.sum(input_data[0:x+1] - mean_value) for x in index_range]
    standard_deviation = np.sqrt(np.sum((input_data-np.mean(input_data))**2)/(data_length-1))
    U_statistic = np.sum((sum_diff[0:(data_length - 2)]/standard_deviation)**2)/(data_length * (data_length + 1))
    absolute_sum_diff = np.abs(sum_diff)
    max_absolute_sum_diff = np.max(absolute_sum_diff)
    change_point_index = list(absolute_sum_diff).index(max_absolute_sum_diff) + 1
    normalized_sum_diff = (sum_diff/standard_deviation)
    return change_point_index


# 检测变点
change_point = Buishand_U_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()

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