如何确定vibro_coeff的饱和值与饱和时间？求算法优化及错误排查

一、现有代码的核心错误 #

• 循环变量冲突：外层for i in range(0,701)的i被内部列表推导[vibro_coeff[700-i] for i in range(10)]覆盖，导致索引计算完全混乱，这是最致命的错误。
• 时间计算逻辑错误：saturate_time = i/0.1 -1的推导毫无依据，假设采样间隔是0.1秒，索引idx对应的时间应为idx * 0.1，反向搜索时的时间需要重新对应。
• 稳定判定逻辑颠倒：反向搜索时，当平均斜率大于err，说明这段数据还在变化（未稳定），此时不应直接取这段的平均值作为饱和值，饱和值应该取后续稳定段的均值。

二、修正后的基础实现代码 #

import numpy as np

def saturation_cal(vibro_coeff, err=0.0001, window_size=10, sample_interval=0.1):
    # 转换为numpy数组提升计算效率
    vibro_arr = np.array(vibro_coeff)
    total_points = len(vibro_arr)
    
    # 从后往前滑动窗口遍历
    for start_idx in range(total_points - window_size, -1, -1):
        # 提取当前窗口数据
        window_data = vibro_arr[start_idx:start_idx+window_size]
        # 计算窗口内相邻点的斜率绝对值均值
        slope_abs = np.abs(np.diff(window_data))
        ave_slope = np.mean(slope_abs)
        
        # 当斜率均值大于err，说明当前窗口未稳定，后续窗口为稳定段
        if ave_slope > err:
            # 取第一个稳定窗口的均值作为饱和值
            saturate_window = vibro_arr[start_idx+window_size:start_idx+2*window_size]
            ave_vibro_coeff = np.mean(saturate_window)
            # 饱和时间为稳定窗口的起始时间
            saturate_time = (start_idx + window_size) * sample_interval
            return ave_vibro_coeff, saturate_time
    
    # 若所有窗口都稳定，取最后一段窗口的均值和对应时间
    ave_vibro_coeff = np.mean(vibro_arr[-window_size:])
    saturate_time = (total_points - window_size) * sample_interval
    return ave_vibro_coeff, saturate_time

三、更鲁棒的优化实现 #

针对噪声干扰的情况，增加连续稳定窗口的判定条件，避免单窗口误判：

import numpy as np

def saturation_cal_optimized(vibro_coeff, err=0.0001, window_size=10, sample_interval=0.1, stable_count=3):
    vibro_arr = np.array(vibro_coeff)
    total_points = len(vibro_arr)
    current_stable = 0
    
    # 从后往前遍历窗口，统计连续稳定次数
    for start_idx in range(total_points - window_size, -1, -1):
        window_data = vibro_arr[start_idx:start_idx+window_size]
        slope_abs = np.abs(np.diff(window_data))
        ave_slope = np.mean(slope_abs)
        
        if ave_slope <= err:
            current_stable += 1
            # 达到连续稳定次数，判定为系统稳定
            if current_stable >= stable_count:
                # 取连续稳定窗口的均值作为饱和值
                saturate_window = vibro_arr[start_idx:start_idx+window_size*stable_count]
                ave_vibro_coeff = np.mean(saturate_window)
                saturate_time = start_idx * sample_interval
                return ave_vibro_coeff, saturate_time
        else:
            current_stable = 0
    
    # 若全量数据都稳定，取最后一段连续稳定窗口的结果
    ave_vibro_coeff = np.mean(vibro_arr[-window_size*stable_count:])
    saturate_time = (total_points - window_size*stable_count)*sample_interval
    return ave_vibro_coeff, saturate_time

四、进阶拟合方案 #

如果数据符合指数趋近饱和的规律（如y = A*(1 - e^(-t/τ)) + B），直接拟合该模型能得到更精准的饱和值和时间：

import numpy as np
from scipy.optimize import curve_fit

def saturation_fit(vibro_coeff, sample_interval=0.1):
    # 构造时间轴
    t = np.arange(len(vibro_coeff)) * sample_interval
    # 定义指数饱和模型
    def exp_saturation(t, A, tau, B):
        return A*(1 - np.exp(-t/tau)) + B
    
    # 拟合模型
    popt, _ = curve_fit(exp_saturation, t, vibro_coeff)
    A, tau, B = popt
    # 饱和值为A+B，通常取达到95%饱和值的时间作为饱和时间（t=3τ）
    saturate_value = A + B
    saturate_time = 3 * tau
    return saturate_value, saturate_time

内容的提问来源于stack exchange，提问作者Tanmia