6.1 谱峰漂移与样品不均匀
6.1 谱峰漂移与样品不均匀
本节摘要：两类"位置病"的会诊记录。峰位漂移用导数域互相关测量、抛物线插值把精度压到亚采样间隔（实测误差 0.1 波数内），再决定对齐还是回仪器端；样品不均匀用方差分解拆开"样品间"与"装样间"两个方差，组内相关系数一个数回答取样代表性够不够。
复诊登记簿的第一页贴着两张谱图：同一份标样，周一测一次、周五测一次，峰位肉眼看着没动，谱库比对分数却掉了四成——典型的漂移投诉。第二页是一瓶粉末装样三次的三条谱，最高与最低差了一成——不均匀投诉。两张单子开错了药方都是白忙：漂移病拿平滑治（治标不治本），不均匀病拿基线校正治（方向全错）。先诊断，后开药。
一、漂移：先量出来，再谈治
漂移的来源在 1.3 节的五段链路里点过名：温度变化让光栅与干涉仪的几何微变、仪器搬运后的重新定标、甚至样品池换批——半天的温差能让峰位漂零点几个波数。定性比对首看峰位（1.1 节的"峰位即身份"），漂移一两个波数就足以让相关系数明显掉分。
测量漂移的标准手法是互相关：把两条谱做相关，相关峰的位置就是相对位移。两个工程细节决定成败。在导数谱上做——原始谱的互相关会被基线与整体幅度干扰，一阶导数把这些压掉，只留峰形信息；做亚采样插值——互相关的位移分辨率卡在采样间隔上（本案 0.5 波数），用相关峰及其两侧两点拟合抛物线，顶点位置就是细化后的位移：
# 漂移测量：导数域互相关 + 抛物线插值 import numpy as np from scipy.signal import correlate, savgol_filter rng = np.random.default_rng(13) x = np.linspace(900, 1200, 600) dx = x[1] - x[0] def ref_spec(): return (0.8*np.exp(-0.5*((x-1000)/8)**2) + 0.5*np.exp(-0.5*((x-1080)/12)**2) + 0.35*np.exp(-0.5*((x-1140)/10)**2)) ref = ref_spec() def shifted(shift_wn, noise=0.012): y = np.interp(x, x + shift_wn, ref_spec()) # 沿波数轴平移 return y + rng.normal(0, noise, x.size) def xcorr_shift(y): d_ref = savgol_filter(ref, 11, 3, deriv=1) # 窄窗导数压基线 d_y = savgol_filter(y, 11, 3, deriv=1) cc = correlate(d_y - d_y.mean(), d_ref - d_ref.mean(), mode='same') k = np.argmax(cc) ym, y0, yp = cc[k-1], cc[k], cc[k+1] # 相关峰及两侧 delta = 0.5*(ym - yp)/(ym - 2*y0 + yp) # 抛物线顶点偏移 return (k - len(x)//2 + delta) * dx for s in [0.0, 1.5, -2.0, 4.0, -3.5]: print(f"true shift={s:+.1f} xcorr estimate={xcorr_shift(shifted(s)):+.2f}")
输出：
true shift=+0.0 xcorr estimate=-0.02 true shift=+1.5 xcorr estimate=+1.48 true shift=-2.0 xcorr estimate=-1.99 true shift=+4.0 xcorr estimate=+4.02 true shift=-3.5 xcorr estimate=-3.42
五档漂移全部测准到 0.1 波数以内，包括零漂移那档没有把噪声当漂移（-0.02）。插值那一行值得多看一眼：采样间隔 0.5 波数本来是分辨率下限，抛物线拟合把精度又往下压了一个量级——免费的精度，只差三行代码。
测出来之后怎么治，看两个数的关系。漂移小于最窄峰半峰宽的一半：软件对齐即可——用测得的位移把谱插值回标准坐标，谱库比对照常进行；漂移接近或超过半峰宽：软件对齐会开始扭曲峰形，病根大概率在仪器端（温控失效、需重新定标），回 1.3 节的链路逐段排查。本案的五档漂移（最多 4 波数）对半峰宽近 20 的峰属于轻症，对齐即愈；但若委托的是高分辨气体谱（半峰宽 0.05），同样的漂移就是绝症级——同样的病，峰宽不同，处方完全不同。
对齐本身一行插值：y_aligned = np.interp(x, x + shift, y_raw)。注意对齐后的谱在边缘各让出漂移量宽度，别把外推段喂进比对。
二、不均匀：方差分解出诊断书
第二张投诉单的病不一样。同一瓶粉末三次装样谱高低差一成——这不是噪声（噪声是逐点随机），是装样差异（整条谱的系统偏移）。要回答的问题是：样品与样品之间的真实差异，够不够盖过装样带来的抖动？如果盖不过，任何建模都是在给装样噪声建模型。
量化工具是方差分解。安排八份样品、每份装样测三次，主峰高度记成八乘三的矩阵。把总波动拆成两层：样品间方差（真实化学差异）与装样间方差（操作抖动），两者之比构成组内相关系数 ICC——它是"你测到的差异里有多少属于样品本身"的占比：
# 样品不均匀：方差分解与 ICC n_samples, n_reps = 8, 3 peak = np.zeros((n_samples, n_reps)) for i in range(n_samples): amp = 1.0 + rng.normal(0, 0.05) # 样品间真实差异 5% for r in range(n_reps): peak[i, r] = amp*(1 + rng.normal(0, 0.04)) # 装样间抖动 4% overall = peak.mean() between = n_reps*((peak.mean(1) - overall)**2).sum()/(n_samples - 1) # 样品间 within = peak.var(1, ddof=1).mean() # 装样间 icc = between/(between + within) print(f"sample-level std={peak.mean(1).std():.4f} rep-level std={peak.std(1).mean():.4f}") print(f"ICC={icc:.3f}")
输出：
sample-level std=0.0466 rep-level std=0.0205 ICC=0.896
ICC 等于 0.896：测得差异的九成属于样品本身、一成来自装样。判定线的惯例是 0.75 以上取样代表性可接受、0.5 到 0.75 边缘、0.5 以下说明装样抖动反客为主。本案 0.896 属健康——注意它背后是"样品差 5% 对装样差 4%"的组合，样品间差异越大 ICC 越好看，反之若一批样品本来就相近（掺假比例都低于 2%），同样的装样手法 ICC 会直接掉线。所以 ICC 要按任务分段评估：想分辨 2% 的差异，先算清这个浓度段的 ICC 够不够。
ICC 低于 0.75 时的处方分两级。抖动是乘性为主（谱整体缩放）：SNV 或 MSC 能压掉大半（2.5 节），模型端有救；取样差为主（每次舀到的粉末成分真的不同）：任何后处理都无济于事，只能扩取样——混匀后分多点取样合并、或每份样品测三到五次取平均，用根号 N 法则把装样方差摊薄。会诊经验：ICC 低时先补三次重复测量，比换任何预处理方法都便宜有效。
⚠️ 常见坑：用"三次测量的平均谱"进训练集、却用单次测量做预测。平均把装样方差压掉了，模型从没见过单次测量的抖动，遇到预测样本直接水土不服——训练与预测的测量协议必须一致。
本节要点回顾
- 漂移先量后治：导数域互相关加抛物线插值，0.1 波数内测准，含零漂移档不误报；
- 处方看峰宽：小于半峰宽之半软件对齐，超过则回仪器端，同样的病峰宽不同处方不同；
- ICC 量化取样代表性：样品间方差占比，0.75 可接受、0.5 以下抖动反客为主；
- ICC 低时分两路：乘性抖动上 SNV 或 MSC，取样差异扩重复测量，后者更便宜；
- 训练与预测的测量协议必须一致，平均谱训练单次预测是大坑。
位置病诊完，最后一间诊室治"选择病"：几十种预处理组合摆在面前，怎么选才不算掷骰子。