2.4 ALS 自动基线实战

2.4 ALS 自动基线实战

本节摘要：不对称最小二乘（ALS）把基线估计写成一个带两项诉求的目标函数——贴着谱的底部走，但又要保持平滑。两个超参数（平滑度因子与非对称权重）各管一头，十来轮迭代自动收敛。它是批量处理数百条谱时唯一现实的手法，也是本章从手工走向规模化的拐点。

2.3 的迭代多项式已经让算法接手了一半判断，但锚点窗口的选取仍要人盯。设想一个真实场景：三百家橄榄油作坊各送一份近红外谱，你明天要交鉴定结论——逐条挑点的方案当场破产。ALS（Asymmetric Least Squares，不对称最小二乘）把整条谱一次交给算法，人只调两个参数。

一、两个诉求写进同一个目标函数

ALS 的想法朴素：基线应该"贴底"且"平滑"。把两个诉求写成目标函数，就是最小化下面这个量：

加权残差平方和 + 平滑度因子 × 二阶差分平方和

第一项要求基线贴近数据（贴底），但加权不对称——高于当前基线的点（可能是峰）给小权重，低于它的点给大权重，这样基线被"压"在谱的底部而不是穿峰而过；第二项惩罚基线的弯曲（平滑），二阶差分衡量曲率，因子越大基线越懒得起伏。两个诉求天然打架：太贴底会爬进峰脚，太平滑会垫到峰下。讨价还价的结果由两个超参数裁定——平滑度因子 lam 管第二项的嗓门，非对称权重 p 管第一项的偏心程度。

实现只有十几行，稀疏矩阵让它在几千个点的谱上瞬间跑完：

# ALS 基线校正（Eilers 与 Boelens 的经典实现）
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve

def als(y, lam=1e6, p=0.01, niter=15):
    L = len(y)
    D = sparse.diags([1, -2, 1], [0, -1, -2], shape=(L, L-2))  # 二阶差分算子
    DTD = lam * (D @ D.T)                                      # 平滑惩罚项
    w = np.ones(L)                                             # 初始等权
    for _ in range(niter):
        W = sparse.diags(w)
        z = spsolve((W + DTD).tocsc(), w * y)                  # 解加权最小二乘
        w = p * (y > z) + (1 - p) * (y < z)                    # 不对称重新加权
    return z

# 沿用 2.3 的模拟档案与真值
rng = np.random.default_rng(21)
x = np.linspace(0, 100, 400)
true_bl = 0.30 + 0.006*x - 6.0e-5*x**2 + 0.25*np.exp(-0.5*((x-88)/25)**2)
peaks = (1.0*np.exp(-0.5*((x-30)/2.5)**2) + 0.6*np.exp(-0.5*((x-52)/3.0)**2)
         + 0.35*np.exp(-0.5*((x-70)/2.0)**2))
y = true_bl + peaks + rng.normal(0, 0.01, x.size)

z = als(y, lam=1e6)
print(f"lam=1e6: bl_RMSE={np.sqrt(np.mean((z-true_bl)**2)):.4f}")
​

输出：

lam=1e6: bl_RMSE=0.0270
​

不挑一个锚点、不给阶数，基线误差 0.027——与 2.3 手工三次多项式的 0.0164 同一量级，而人工介入为零。

二、两个旋钮各自拧过头会怎样

ALS 只有 lam 与 p 两个旋钮，但拧过头的样子必须亲眼见过才认得。先扫 lam（三个数量级），再扫 p（两个数量级）：

# lam 与 p 的敏感度实验
m_pk1 = (x >= 25) & (x <= 35)
m_val = (x >= 40) & (x <= 47)
def rmse(a, b): return float(np.sqrt(np.mean((a-b)**2)))

for lam in [1e3, 1e6, 1e9]:
    z = als(y, lam=lam)
    print(f"lam={lam:.0e}: bl_RMSE={rmse(z, true_bl):.4f}  "
          f"rec_peak1={(y-z)[m_pk1].max():.3f}  valley={(y-z)[m_val].min():.4f}")
for p in [0.001, 0.01, 0.1]:
    z = als(y, lam=1e6, p=p)
    print(f"p={p}: bl_RMSE={rmse(z, true_bl):.4f}  "
          f"rec_peak1={(y-z)[m_pk1].max():.3f}  valley={(y-z)[m_val].min():.4f}")
​

输出：

lam=1e3: bl_RMSE=0.0912  rec_peak1=0.741  valley=-0.0180
lam=1e6: bl_RMSE=0.0270  rec_peak1=0.987  valley=-0.0242
lam=1e9: bl_RMSE=0.0760  rec_peak1=1.044  valley=0.0539
p=0.001: bl_RMSE=0.0629  rec_peak1=1.019  valley=0.0200
p=0.01: bl_RMSE=0.0270  rec_peak1=0.987  valley=-0.0242
p=0.1: bl_RMSE=0.0458  rec_peak1=0.929  valley=-0.0819
​

lam=1e3，基线太灵活：它爬进了峰里，峰高只剩 0.741——四分之一的峰被当成背景拆走了，这是"咬峰"的极端形态；lam=1e9，基线太僵硬：一条几乎不弯的线被峰顶起来，垫在峰下，表观峰高虚高到 1.044、峰谷还残留正背景 0.054——峰没坏，但定量基准错了。p 的方向恰好相反：p 越大，峰上点的权重越大，基线越被往上拉（p=0.1 时峰高被压到 0.929）；p 越小基线越沉底（p=0.001 时谷底残留正偏差 0.020）。三档 lam、三档 p，把两个旋钮的失败模式各演示了一遍——**调参不是玄学，是认出过拟合与欠拟合在基线问题上的长相**。

实用起点由此可以给出：拉曼荧光背景取 lam 在 1e5 到 1e7、p 在 0.001 到 0.01；近红外散射背景形状更平，lam 可以更大。定了起点后，用峰高恢复率与峰谷残差两个指标微调——2.3 的"标定试金石"方法在真实数据上依然适用：挑两三条有代表性的谱，人工校一次真值，其余全靠它对账。

三、批处理与它的边界

ALS 真正的价值在批处理。同一套参数跑三百条谱，一致性本身就是质检指标：把每条谱估出的基线在某个无峰波数处的取值排成一列，若有谱的基线明显离群，先怀疑装样或仪器异常，而不是急着改参数——基线是采集状态的记录仪，第 6 章会诊样品不均匀问题时还要用它倒查现场。

边界也要说清。其一，ALS 默认背景缓变，遇到极窄的强峰（半高宽只有几个点）可能把峰脚当背景；其二，两边对称的双肩峰之间如果凹陷很深，基线可能"抄近路"横穿谷底；其三，它与 2.2 的去噪有顺序依赖——噪声太大会干扰不对称加权的判断，所以工序表上它排在去噪之后不是偶然。

💡 关键直觉：lam 管基线"多懒"，p 管基线"多偏心"。记不住方向时想一句话——基线是惰性的：它能不动就不动，除非你允许它爬。

本节要点回顾

• 目标函数两项诉求：贴底（不对称加权）与平滑（二阶差分惩罚），由 lam 与 p 裁定平衡；
• lam 过小爬峰（实验里峰高掉到 0.741），lam 过大垫峰（虚高到 1.044）；
• p 过大基线上抬压峰（0.929），p 过小沉底留残差；
• 批处理才是它的主场：基线一致性可当采集质检指标用；
• 窄峰与深谷是它的边界，工序上排在去噪之后。

背景拆完，档案只剩"形状"这一个问题：不同批次、不同装样测出的谱，整体高低不一。下一节把五把标尺摆上桌，讲清哪把尺子量哪种任务。