3. 主机缸套水温基线建模与残差预警

为什么缸套水温不能直接用固定阈值

主机缸套冷却水出口温度受负荷、海水温度、三通阀开度共同影响。同一台主机在 50% 负荷与 85% 负荷下的正常水温可能相差 6~10 ℃，若用单一固定上限报警，低负荷时漏报、高负荷时误报。本例采用工况分层基线 + 残差 + 动态阈值的思路：先按负荷与海水温度把运行工况分箱，在每个箱内用清洗后的历史数据拟合正常水温基线，实时计算残差 r = T_obs − T_base，再用残差的滑动统计量构造随工况变化的动态阈值，最后用一段人工注入的缓慢漂移数据检验预警提前量。

工况分箱与基线拟合

分箱维度选负荷率与海水温度，因为二者对缸套水温的影响最直接且可实时获取。负荷按 25% 间隔分 4 档（0–25、25–50、50–75、75–100），海水温度按 5 ℃ 间隔分 3 档（<20、20–25、>25），共 12 个箱。每箱样本不足 200 点时合并相邻箱，避免小样本导致回归系数不稳定。箱内基线用线性回归：

T_base = a * load + b * T_sea + c

其中 load 为负荷率（0–1），T_sea 为海水温度（℃）。回归后同时记录该箱残差标准差 σ_r，作为动态阈值的基础。以下为可直接运行的 Python 实现（示例数据为教学构造，非实测）：

import numpy as np

# 教学示例：构造 12 个工况箱的历史数据（非实测）
rng = np.random.default_rng(42)
load_bins = [(0.0, 0.25), (0.25, 0.5), (0.5, 0.75), (0.75, 1.0)]
sea_bins = [(15, 20), (20, 25), (25, 30)]

baselines = {}
for li, (l0, l1) in enumerate(load_bins):
    for si, (s0, s1) in enumerate(sea_bins):
        n = 300
        load = rng.uniform(l0, l1, n)
        sea = rng.uniform(s0, s1, n)
        # 真实关系：负荷每升 10% 水温升约 1.2℃，海温每升 1℃ 水温升约 0.6℃
        t_out = 32 + 12 * load + 0.6 * sea + rng.normal(0, 0.35, n)
        X = np.column_stack([load, sea, np.ones(n)])
        coef, *_ = np.linalg.lstsq(X, t_out, rcond=None)
        resid = t_out - X @ coef
        baselines[(li, si)] = {
            "coef": coef,
            "sigma": resid.std(ddof=1),
            "n": n,
        }

def pick_bin(load, sea):
    li = min(int(load / 0.25), 3)
    si = 0 if sea < 20 else (1 if sea < 25 else 2)
    return (li, si)

def residual(load, sea, t_obs):
    li, si = pick_bin(load, sea)
    b = baselines[(li, si)]
    t_base = b["coef"][0] * load + b["coef"][1] * sea + b["coef"][2]
    return t_obs - t_base, b["sigma"]

动态阈值与漂移验证

动态阈值取 UCL = μ_r + k·σ_r，μ_r 用残差滑动均值（窗口 30 点）估计，k=3。为抑制单点噪声，判定条件设为连续 3 点超过 UCL 才触发预警。验证时在 50% 负荷、海温 22 ℃ 的稳定工况下注入 0.05 ℃/min 的缓慢温升，模拟冷却器结垢或流量下降导致的渐进漂移：

# 漂移验证：50% 负荷、海温 22℃，注入 0.05 ℃/min 温升
dt_min = 1
duration = 120
drift_rate = 0.05
load, sea = 0.5, 22.0

li, si = pick_bin(load, sea)
b = baselines[(li, si)]
t_base = b["coef"][0] * load + b["coef"][1] * sea + b["coef"][2]
sigma = b["sigma"]

resid_hist = []
alert_time = None
for t in range(duration):
    t_obs = t_base + drift_rate * t + rng.normal(0, sigma)
    r, _ = residual(load, sea, t_obs)
    resid_hist.append(r)
    if len(resid_hist) >= 30:
        mu = np.mean(resid_hist[-30:])
        ucl = mu + 3 * sigma
        if len(resid_hist) >= 3 and all(x > ucl for x in resid_hist[-3:]):
            alert_time = t
            break

print(f"基线温度 {t_base:.2f} ℃，箱内 σ_r={sigma:.3f} ℃")
print(f"首次预警时刻 t_alert = {alert_time} min")
print(f"相对漂移起点提前量 = {alert_time} min")

示例运行中，基线温度约 41.2 ℃，σ_r 约 0.35 ℃，预警在漂移开始后约 95 min 触发，相对 120 min 的漂移终点提前约 25 min。这一提前量取决于漂移速率与 σ_r 的比值：漂移越慢、噪声越大，越限越晚。实际应用中应把 k 与连续点数作为可调参数，在误报率与提前量之间权衡。