南京地湖所最新研究揭示内蒙古黄旗海近1600年来生态系统演化及驱动机制

STARS 的核心逻辑是设定一个截断长度（Cut-off length, L），通过计算当前点与后续 L 个点的 T 检验，计算出情势转换指数（RSI，Regime Shift Index）。以下是精简版 Python 实现代码：

import numpy as npimport matplotlib.pyplot as pltfrom scipy import stats# ==========================================# 1. STARS 算法核心函数 (稳态转换检测)# ==========================================def calculate_stars(series, L=20, p_value=0.05):    """    STARS (Sequential T-test Algorithm for Regime Shifts) 的Python实现    :param series: 一维时间序列数据 (如降水、温度等序列)    :param L: 样本窗口长度 (论文中设定为 19.4，此处建议取整数，如 20)    :param p_value: 显著性水平 (论文中设定为 0.05)    :return: 包含突变点索引和对应 RSI 值的字典    """    series = np.array(series)    n = len(series)    rsi_dict = {}    # 检查数据长度是否足够    if n < 2 * L:        print("警告：数据长度不足以支撑设定的窗口长度 L")        return rsi_dict    # 根据自由度(2L-2)和p_value计算t检验的临界值    t_crit = stats.t.ppf(1 - p_value/2, 2*L - 2)    # 初始化第一个稳态的均值和方差    current_mean = np.mean(series[:L])    current_var = np.var(series[:L], ddof=1)    for i in range(1, n - L + 1):        # 提取接下来的L个点作为潜在的新稳态窗口        new_window = series[i : i+L]        new_mean = np.mean(new_window)        new_var = np.var(new_window, ddof=1)        # 计算合并方差 (Pooled Variance)        pooled_var = ((L-1)*current_var + (L-1)*new_var) / (2*L - 2)        if pooled_var == 0 or current_var == 0:            continue        # 计算 T 统计量        t_stat = abs(new_mean - current_mean) / np.sqrt(pooled_var * (2/L))        # 如果T统计量大于临界值，说明可能发生了稳态转换 (突变)        if t_stat > t_crit:            # 计算临界均值差异            diff_crit = t_crit * np.sqrt(pooled_var * (2/L))            rsi = 0            # 计算 RSI (稳态转换指数 Regime Shift Index)            for j in range(L):                if new_mean > current_mean:                    rsi += (series[i+j] - (current_mean + diff_crit)) / np.sqrt(current_var)                else:                    rsi += ((current_mean - diff_crit) - series[i+j]) / np.sqrt(current_var)            # 如果 RSI > 0，确认稳态转换发生            if rsi > 0:                # 记录突变点的索引和 RSI 值 (加入简单去重逻辑，避免同一突变期连续报告)                if len(rsi_dict) == 0 or (i - list(rsi_dict.keys())[-1]) > L/2:                    rsi_dict[i] = rsi                    # 突变发生后，将新窗口确立为当前的稳态                    current_mean = new_mean                    current_var = new_var            else:                # 若RSI未过阈值，平滑滚动更新当前稳态均值                current_mean = np.mean(series[i:i+L])                current_var = np.var(series[i:i+L], ddof=1)    return rsi_dict# ==========================================# 2. 数据输入区 (这里替换为实际序列即可)# ==========================================# 这是一组模拟数据：前30个数据在1.5左右波动，后30个数据突变跃迁到3.0左右波动my_data = [    1.52, 1.48, 1.55, 1.41, 1.60, 1.45, 1.51, 1.58, 1.49, 1.42,     1.53, 1.57, 1.46, 1.50, 1.44, 1.59, 1.43, 1.54, 1.47, 1.56,     1.40, 1.61, 1.48, 1.52, 1.45, 1.55, 1.49, 1.51, 1.42, 1.58,     3.01, 2.95, 3.08, 2.85, 3.10, 2.92, 3.05, 3.12, 2.98, 2.88,     3.02, 3.09, 2.91, 3.00, 2.94, 3.15, 2.89, 3.04, 2.96, 3.11,     2.86, 3.14, 2.97, 3.03, 2.90, 3.13, 2.99, 3.06, 2.87, 3.07]# ==========================================# 3. 运行突变检测# ==========================================# 注意：窗口长度 L 的设置需要小于数据总长度的一半。# 这里演示数据总长60，所以 L 设为 15。# (如果实际数据有上百个点，建议按照论文设置为 L=20)shift_points = calculate_stars(my_data, L=15, p_value=0.05)print("检测到的稳态转换(突变)点及对应的 RSI 值：")for idx, rsi in shift_points.items():    print(f"数据索引位置: {idx}, 突变强度(RSI)值: {rsi:.2f}")# ==========================================# 4. 结果可视化# ==========================================plt.figure(figsize=(10, 5))plt.plot(my_data, label='My Data Sequence', color='#2ca02c', marker='o', markersize=4, alpha=0.8)# 标注检测到的突变点for idx, rsi in shift_points.items():    plt.axvline(x=idx, color='red', linestyle='--', linewidth=2,                 label='Regime Shift (稳态突变点)')    # 在线上方标注 RSI 值    plt.text(idx + 1, max(my_data)*0.9, f'RSI: {rsi:.2f}', color='red', fontweight='bold')plt.title("STARS Regime Shift Detection (稳态转换检测)")plt.xlabel("Sequence Index (时间/深度序列索引)")plt.ylabel("Value (如: Chl.a 浓度)")plt.legend()plt.grid(True, linestyle=':', alpha=0.6)plt.show()

三、研究亮点