1. 海洋平台监测数据质量诊断与清洗

任务对象与数据概况

本例针对某导管架平台同一采集周期内的两路监测通道：加速度通道 ACC-01（采样率 50 Hz，单位 m/s²）与应变通道 SG-03（采样率 50 Hz，单位 με）。原始数据为 10 分钟连续记录，共 30000 个采样点。已知该时段内平台处于正常作业工况，未发生碰撞或极端风暴，因此数据中出现的剧烈突变应视为质量问题而非真实结构响应。

四类典型缺陷的判别依据如下：

• 缺失：时间戳不连续，或数值为 NaN / 哨兵值（如 -9999）。
• 跳变：相邻采样点差分绝对值超过通道阈值的突变，且持续时间极短（1~3 点）。
• 漂移：滑动均值在无荷载变化条件下持续单向偏移，表现为低频趋势项。
• 异常尖峰：幅值超过 6 倍稳健标准差（基于 MAD 估计），且与频谱中结构主频无关的孤立高点。

诊断与清洗处理链

清洗流程按“先诊断、后修复、再验证”的顺序组织。诊断阶段对原始序列做分段统计，标记缺陷位置；修复阶段对缺失点采用线性插值，对跳变与尖峰采用 Hampel 滤波（滑动窗口中值替换），对漂移采用高通滤波去除趋势项；验证阶段对比清洗前后的均值、标准差、偏度、峰度以及功率谱在结构主频处的幅值一致性。

关键实现与参数设置

以下 Python 代码完成从原始数据读入到清洗结果输出的完整流程。数据以 NumPy 数组形式给出，实际使用时替换为读取 CSV 或数据库的代码即可。代码中所有阈值均以本例教学参数形式给出，实际工程中应根据通道噪声水平和结构频率特性重新标定。

import numpy as np
from scipy.signal import butter, filtfilt, welch

# ---------- 1. 模拟原始数据（教学示例） ----------
np.random.seed(42)
fs = 50.0
n = 30000
t = np.arange(n) / fs
# 结构响应：0.3 Hz 主频 + 1.2 Hz 次频 + 噪声
acc = (0.8 * np.sin(2 * np.pi * 0.3 * t)
       + 0.3 * np.sin(2 * np.pi * 1.2 * t)
       + 0.05 * np.random.randn(n))
strain = (120 * np.sin(2 * np.pi * 0.3 * t + 0.2)
          + 30 * np.sin(2 * np.pi * 1.2 * t)
          + 5 * np.random.randn(n))

# 注入缺陷
acc_raw = acc.copy()
strain_raw = strain.copy()
# 缺失：随机 200 点置 NaN
miss_idx = np.random.choice(n, 200, replace=False)
acc_raw[miss_idx] = np.nan
# 跳变：3 处短时突变
for i in [5000, 12000, 22000]:
    acc_raw[i:i+3] += 3.0
# 漂移：叠加缓慢线性趋势
acc_raw += 0.0002 * t
# 尖峰：5 处孤立高点
for i in [8000, 15000, 18000, 25000, 28000]:
    strain_raw[i] += 80.0

# ---------- 2. 质量诊断 ----------
def diagnose(x, jump_thr, mad_k=6.0):
    mask_missing = np.isnan(x)
    x_fill = np.where(mask_missing, np.nanmedian(x), x)
    diff = np.abs(np.diff(x_fill, prepend=x_fill[0]))
    mask_jump = diff > jump_thr
    med = np.median(x_fill)
    mad = np.median(np.abs(x_fill - med)) * 1.4826
    mask_spike = np.abs(x_fill - med) > mad_k * mad
    return mask_missing, mask_jump, mask_spike, mad

miss_a, jump_a, spike_a, mad_a = diagnose(acc_raw, jump_thr=0.5)
miss_s, jump_s, spike_s, mad_s = diagnose(strain_raw, jump_thr=15.0)

# ---------- 3. 修复处理 ----------
def hampel(x, win=7, k=3.0):
    x = x.copy()
    half = win // 2
    for i in range(len(x)):
        lo, hi = max(0, i-half), min(len(x), i+half+1)
        seg = x[lo:hi]
        med = np.median(seg)
        mad = np.median(np.abs(seg - med)) * 1.4826
        if mad > 0 and abs(x[i] - med) > k * mad:
            x[i] = med
    return x

def clean_channel(x, jump_thr, hp_cut=0.05):
    # 缺失：线性插值
    mask = np.isnan(x)
    if mask.any():
        idx = np.arange(len(x))
        x = np.interp(idx, idx[~mask], x[~mask])
    # 跳变 + 尖峰：Hampel
    x = hampel(x, win=7, k=3.0)
    # 漂移：高通滤波
    b, a = butter(2, hp_cut / (fs/2), btype='high')
    x = filtfilt(b, a, x)
    return x

acc_clean = clean_channel(acc_raw, jump_thr=0.5)
strain_clean = clean_channel(strain_raw, jump_thr=15.0)

# ---------- 4. 统计特征与频谱一致性验证 ----------
def stats(x):
    return dict(mean=np.mean(x), std=np.std(x),
                skew=float(np.mean(((x-x.mean())/x.std())**3)),
                kurt=float(np.mean(((x-x.mean())/x.std())**4)))

def main_peak_amp(x):
    f, P = welch(x, fs=fs, nperseg=4096)
    band = (f > 0.2) & (f < 0.5)
    return f[band][np.argmax(P[band])], P[band].max()

for name, raw, cln in [("ACC-01", acc_raw, acc_clean),
                        ("SG-03", strain_raw, strain_clean)]:
    s_raw, s_cln = stats(raw[~np.isnan(raw)]), stats(cln)
    f_raw, a_raw = main_peak_amp(raw[~np.isnan(raw)])
    f_cln, a_cln = main_peak_amp(cln)
    dev = abs(a_cln - a_raw) / a_raw * 100
    print(f"{name}: 主频 {f_raw:.3f}→{f_cln:.3f} Hz, "
          f"幅值偏差 {dev:.2f}%, 标准差 {s_raw['std']:.3f}→{s_cln['std']:.3f}")

结果判读与通过条件

运行上述代码后，应重点观察以下指标：