拉曼光譜顏料分析入門

第一篇　觀念：光譜在說什麼（不寫程式）

第 1 章　拉曼光譜到底在量什麼

拿一道雷射光照在顏料上，絕大部分的光會原封不動彈回來，但有極少數（大約一億分之一）的光子， 會把一點點能量交給分子，讓分子裡的原子「振動」起來，於是彈回來的光子能量變小了一點點。 這個能量差，正好等於那個振動模式的能量。

不同的化學鍵、不同的晶體結構，振動的頻率就不同。所以量出這些能量差， 就等於直接讀到分子的「指紋」。這就是拉曼光譜。

第 2 章　一張光譜長什麼樣子

打開你手上任何一個光譜檔，裡面就是兩欄數字：

77.3545 0
80.4019 924
83.4475 967
...

左邊那欄叫拉曼位移（Raman shift），單位是「每公分幾個波」，寫成 cm⁻¹，唸作「倒數公分」。 它就是前面說的那個能量差。右邊那欄是強度，單位是偵測器數到的光子數（counts）， 這個絕對值沒有意義，重要的是相對高低。

光譜裡一根一根往上凸的東西叫峰（peak 或 band，中文也叫「帶」）。 一個峰要記三個數字：

第 3 章　為什麼一定要扣基線

真實的光譜長這樣：峰坐在一大坨圓滾滾的隆起上面。那坨隆起叫螢光背景， 是樣品裡的有機物被雷射激發後發出的螢光，強度往往是拉曼訊號的好幾倍甚至幾十倍。

它會造成兩個麻煩：第一，峰的「高度」要從哪裡開始量？第二， 兩張光譜的背景不一樣高，就沒辦法互相比較。所以第一步一定是把它扣掉。

我們用的方法叫 ALS（Asymmetric Least Squares，非對稱最小平方）。名字嚇人， 概念其實一句話：找一條又平滑又儘量貼著「谷底」走的線。

「平滑」跟「貼著谷底」是兩個互相拉扯的要求，所以有兩個旋鈕：

• λ（lambda）控制平滑度。λ 越大，線越硬越直；越小越軟越會扭來扭去。
• p 控制不對稱程度。演算法會問每個點「你在線的上面還是下面？」 在上面（可能是峰）的點權重壓到 p，在下面的點權重給 1-p。 p 越小，線越死命往谷底鑽。

第 4 章　訊噪比：怎麼知道這是峰還是雜訊

扣完基線之後，光譜看起來是平的，但仔細看還是有細細的抖動——那是雜訊。 問題來了：一個小凸起，到底是真的峰，還是雜訊剛好連續往上幾點？

判斷方法是訊噪比（signal-to-noise ratio，簡稱 S/N）：

S/N ＝ 峰的突出高度 ÷ 雜訊的標準差 σ

雜訊 σ 怎麼估？最省事又可靠的做法是看相鄰兩點的差：

σ ≈ standard_deviation(diff(y)) / √2

因為真實的峰是平滑變化的，相鄰點差很小；雜訊卻是隨機跳動，相鄰點差就是雜訊本身。 除以 √2 是因為「兩個獨立雜訊相減」的變異數是單個的兩倍。

第 5 章　鑑定的邏輯：三道門檻

比對譜庫聽起來很簡單：把實測的峰位跟文獻的帶位對一對，對上了就是。 但實際上這樣做非常容易出錯。我們用三道門檻來把關。

門檻一：關鍵帶必須全部出現

每個顏料的文獻帶位有一長串，但其中只有一兩個是「非有不可」的。 例如硫酸鋇最強的帶在 988 cm⁻¹，若一張光譜沒有 988，那無論其他帶怎麼湊，都不可能是硫酸鋇。 我們把這種帶叫關鍵帶（key band）。

門檻二：那個位置必須真的是一個「峰」

這是最容易被忽略、也最容易出錯的一點。看下面這張圖：

解法：一個文獻帶要算「出現」，必須同時滿足強度夠高而且該位置本身被偵測為一個峰。 斜坡上的點不是峰，這樣就擋掉了。

門檻三：主帶強度要跟「含量」相稱

如果某個顏料是樣品的主成分，那它文獻上最強的那一帶， 在你這張光譜裡通常也會是數一數二高的峰。我們要求主帶強度 ≥ 全譜最強峰的 25%， 過了才叫「主成分」；沒過但前兩關過了，就標成「次要成分」。

第二篇　Python 路線

第 6 章　把環境準備好

只需要三個套件：numpy（數值運算）、scipy（訊號處理與最佳化）、 matplotlib（畫圖）。

pip install numpy scipy matplotlib

如果你裝的是 Anaconda，這三個本來就有，可以直接跳過。要確認裝好了沒：

import numpy, scipy, matplotlib
print(numpy.__version__, scipy.__version__, matplotlib.__version__)

import matplotlib.pyplot as plt
plt.rcParams["font.family"] = ["Microsoft JhengHei"]
plt.rcParams["axes.unicode_minus"] = False   # 不然負號也會變方框

第 7 章　讀檔與畫出第一張圖

最單純的寫法：

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt("cinnabar_1.csv", delimiter=",", skiprows=1)
x, y = data[:, 0], data[:, 1]      # 第一欄波數，第二欄強度

plt.figure(figsize=(10, 4))
plt.plot(x, y, lw=1)
plt.xlabel("拉曼位移 (cm$^{-1}$)")
plt.ylabel("強度")
plt.xlim(150, 1400)                # 低於 150 是濾光片邊緣，不是訊號
plt.show()

但真實世界的檔案格式五花八門，所以我們寫一個會自己判斷的版本。 關鍵是：把「逗號 / Tab / 分號 / 空白」四種都試一遍，看哪一種能切出兩個數字。

def load_spectrum(path):
    with open(path, encoding="utf-8", errors="ignore") as f:
        head = [f.readline() for _ in range(5)]

    def try_split(line, d):
        parts = line.strip().split(d) if d is not None else line.split()
        if len(parts) < 2:
            return False
        try:
            float(parts[0]); float(parts[1]); return True
        except ValueError:
            return False

    delim, skip, found = None, 0, False
    for s in (0, 1):                      # s=0 沒表頭，s=1 有一行表頭
        for d in (",", "\t", ";", None):  # None = 以任意空白切
            if try_split(head[s], d):
                delim, skip, found = d, s, True
                break
        if found:
            break
    if not found:
        raise ValueError("看不懂這個檔案的格式")

    data = np.loadtxt(path, delimiter=delim, skiprows=skip, usecols=(0, 1))
    x, y = data[:, 0], data[:, 1]
    order = np.argsort(x)      # 有的儀器由高波數往低存，統一排序
    return x[order], y[order]

第 8 章　ALS 扣基線

把第 3 章的觀念寫成程式。核心是解一個線性方程組：

from scipy.sparse import diags
from scipy.sparse.linalg import spsolve

def als_baseline(y, lam=1e5, p=0.01, n_iter=12):
    L = len(y)
    # D 是二階差分矩陣。lam * D'D 就是「不准彎太兇」的懲罰項
    D = diags([1, -2, 1], [0, -1, -2], shape=(L, L - 2), dtype=float)
    DtD = lam * D.dot(D.T)

    w = np.ones(L)                 # 每個點的權重，一開始都是 1
    W = diags(w, 0)
    for _ in range(n_iter):
        W.setdiag(w)
        z = spsolve((W + DtD).tocsc(), w * y)      # 解出這一輪的基線 z
        w = p * (y > z) + (1 - p) * (y < z)        # 在 z 上面的點（峰）權重壓到 p
    return z

那一行 w = p * (y > z) + (1 - p) * (y < z) 是整個演算法的靈魂。 (y > z) 是布林陣列，在 numpy 裡跟數字相乘會自動變成 0/1， 所以這一行的意思是「在基線上面的點給權重 p=0.01，下面的點給 0.99」。 下一輪解方程時，基線就會被拉向那些權重高的點——也就是谷底。

接著平滑一下，並估雜訊：

from scipy.signal import savgol_filter

def preprocess(x, y, lam=1e5, p=0.01):
    base = als_baseline(y, lam, p)
    yc = savgol_filter(y - base, 9, 3)   # 視窗 9 點、三次多項式
    noise = np.std(np.diff(yc)) / np.sqrt(2)
    return yc, base, noise

第 9 章　找峰

scipy.signal.find_peaks 幫我們做完了。重點是選對判斷條件—— 要用 prominence（突出高度）而不是 height（絕對高度）。

from scipy.signal import find_peaks

def detect_peaks(x, yc, noise, snr=6.0, xmin=150.0, xmax=1800.0, min_sep=4):
    m = (x >= xmin) & (x <= xmax)
    xm, ym = x[m], yc[m]
    idx, props = find_peaks(ym, prominence=snr * noise, distance=min_sep)
    return [{"position": float(xm[j]),
             "height":   float(ym[j]),
             "snr":      float(props["prominences"][i] / noise)}
            for i, j in enumerate(idx)]

想要更精確的峰位和半高寬，用 Lorentzian 擬合：

from scipy.optimize import curve_fit

def lorentzian(x, a, x0, g, c):
    return a * g**2 / ((x - x0)**2 + g**2) + c

def refine_peak(x, yc, center, halfwidth=20.0):
    m = (x > center - halfwidth) & (x < center + halfwidth)
    xi, yi = x[m], yc[m]
    p0 = [yi.max() - np.median(yi), xi[np.argmax(yi)], 8.0, np.median(yi)]
    popt, _ = curve_fit(lorentzian, xi, yi, p0=p0, maxfev=40000)
    fwhm = abs(2 * popt[2])
    # 擬合失控要擋掉：峰心跑掉，或寬度寬到超過擬合視窗
    if abs(popt[1] - center) > halfwidth or fwhm > 4 * halfwidth:
        return None
    return {"center": popt[1], "fwhm": fwhm}

第 10 章　比對譜庫

譜庫就是一個字典。每個顏料記四件事：所有文獻帶位、關鍵帶、主帶、說明。

LIBRARY = {
    "辰砂 cinnabar (HgS)": {
        "bands": [253, 284, 343], "key": [253, 343], "main": 253,
        "note": "硃砂／銀硃。古代至今通用的紅色顏料。"},
    "金紅石 rutile (TiO2)": {
        "bands": [144, 232, 447, 609], "key": [447, 609], "main": 447,
        "note": "鈦白 PW6 的一種晶型。顏料級 1938 年起商業生產。"},
    "普魯士藍 PB27": {
        "bands": [276, 538, 950, 2091, 2154], "key": [2154], "main": 2154,
        "note": "亞鐵氰化鐵。1704 年發明。2154 落在拉曼靜默區，極好認。"},
    # ...腳本裡共 15 種
}

比對時，把第 5 章的三道門檻寫成程式：

def band_height(x, yc, center, tol=6.0):
    m = (x > center - tol) & (x < center + tol)
    return float(yc[m].max())

def check_band(x, yc, noise, peak_positions, b, tol=6.0, strong=8.0):
    # 門檻一：強度要夠
    tall  = band_height(x, yc, b, tol) > strong * noise
    # 門檻二：那裡真的要有一個峰（這行就是擋掉「肩部假陽性」的關鍵）
    is_pk = np.min(np.abs(peak_positions - b)) <= tol
    return tall and is_pk

第 11 章　混合物解混：NNLS

真實的顏料層幾乎都是混合物。假設未知樣的光譜可以寫成幾個已知標準品的加權和：

未知樣 ≈ c₁ × 參考譜₁ + c₂ × 參考譜₂ + ... + 常數項

求那些係數 c 就是最小平方問題。但有個額外要求：係數不能是負的 （「含有負 30% 的鈦白」沒有意義）。這就是 NNLS（非負最小平方）。

from scipy.optimize import nnls

def unmix(x, yc, refs, lo=350.0, hi=1750.0):
    g = (x >= lo) & (x <= hi)
    X, Y = x[g], yc[g]

    cols = []
    for name, xr, yr in refs:
        r = np.interp(X, xr, yr)      # 內插到同一個 x 格點（很重要！）
        cols.append(r / r.max())      # 各自正規化，係數才好比較
    cols.append(np.ones_like(X))      # 常數項，吸收殘餘背景

    A = np.vstack(cols).T
    coef, _ = nnls(A, Y)
    fit = A @ coef
    return {"coef": coef, "fit": fit, "residual": Y - fit,
            "r": np.corrcoef(fit, Y)[0, 1]}

第 12 章　直接用打包好的腳本

以上全部整理成一支 raman_pipeline.py，附在教材資料夾裡。

# 分析單一檔案
python raman_pipeline.py spectra/UNK2_blue.csv

# 分析整個資料夾
python raman_pipeline.py spectra/ -o 結果

# 加上參考譜做 NNLS 解混
python raman_pipeline.py spectra/UNK3_blue.csv \
       --refs spectra/PB15_phthalo_blue.csv spectra/rutile_titanium_white.csv

# 看看內建譜庫有哪些顏料
python raman_pipeline.py --list-library

每個樣品會產生三個檔案：_report.txt（文字報告）、 _peaks.csv（峰位表，可以直接丟進 Excel）、_report.png（三步驟圖）。 文字報告長這樣：

樣品：UNK2_blue
波數範圍：170.0 – 3200.0 cm⁻¹　資料點：3031
雜訊 σ ≈ 45.1　訊背比＝22.44

── 偵測到的峰（S/N ≥ 6，150–1800 cm⁻¹）──
      279.0 cm⁻¹   高度=     1954   S/N=  44.5　FWHM=38.7
      532.0 cm⁻¹   高度=     1627   S/N=  37.3　FWHM=35.8
     1001.0 cm⁻¹   高度=      951   S/N=  22.1　FWHM=13.5
     ...

── 內建譜庫比對 ──
   ★ 主成分  　普魯士藍 PB27　命中 5✔+0△/5（100%）　主帶 2154 強度佔比 100%
       276:1954✔  538:1615✔  950:373✔  2091:2180✔  2154:8976✔
   ◆ 次要成分　聚苯乙烯系樹脂　命中 3✔+6△/9（67%）　主帶 1001 強度佔比 11%

── 初步判定 ──
   主成分　：普魯士藍 PB27
   次要成分：聚苯乙烯系樹脂

你也可以把它當成模組來用，在自己的程式裡呼叫：

import raman_pipeline as rp

x, y = rp.load_spectrum("spectra/cinnabar_1.csv")
yc, base, noise = rp.preprocess(x, y)
peaks = rp.detect_peaks(x, yc, noise, snr=6)
print(f"最強的峰在 {peaks[0]['position']:.1f} cm-1，S/N = {peaks[0]['snr']:.0f}")

第三篇　R 路線

第 13 章　R 環境與讀檔

好消息：這條路線一個套件都不用裝。我們只用 base R 加上 Matrix， 而 Matrix 是 R 官方隨附的 recommended 套件，裝好 R 就有了。

library(Matrix)

df <- read.table("cinnabar_1.csv", sep = ",", skip = 1)
x <- df[[1]]; y <- df[[2]]

plot(x, y, type = "l", xlim = c(150, 1400),
     xlab = "拉曼位移 (cm-1)", ylab = "強度")

自動判斷格式的版本，邏輯和 Python 完全一樣：

load_spectrum <- function(path) {
  head5 <- readLines(path, n = 5, warn = FALSE)
  ok <- function(line, sep) {
    p <- if (is.null(sep)) strsplit(trimws(line), "[[:space:]]+")[[1]]
         else strsplit(trimws(line), sep, fixed = TRUE)[[1]]
    if (length(p) < 2) return(FALSE)
    all(!is.na(suppressWarnings(as.numeric(p[1:2]))))
  }
  seps <- list(",", "\t", ";", NULL)
  delim <- NA; skip <- 0; found <- FALSE
  for (s in 0:1) {
    for (d in seps) {
      if (length(head5) > s && ok(head5[s + 1], d)) {
        delim <- d; skip <- s; found <- TRUE; break
      }
    }
    if (found) break
  }
  if (!found) stop("看不懂這個檔案的格式")
  df <- if (is.null(delim)) read.table(path, skip = skip)
        else read.table(path, sep = delim, skip = skip)
  o <- order(df[[1]])
  list(x = df[[1]][o], y = df[[2]][o])
}

第 14 章　用 Matrix 做 ALS 基線

ALS 需要解一個很大的線性方程組（有幾千個未知數），但矩陣是「帶狀」的—— 只有主對角線附近有非零元素。Matrix 套件的稀疏矩陣就是為此而生， 不然幾千乘幾千的稠密矩陣會把記憶體吃光。

als_baseline <- function(y, lambda = 1e5, p = 0.01, n_iter = 12) {
  L <- length(y)
  # 二階差分矩陣：每一列是 (1, -2, 1) 往右移一格
  D <- bandSparse(L - 2, L, k = c(0, 1, 2),
                  diagonals = list(rep(1, L-2), rep(-2, L-2), rep(1, L-2)))
  DtD <- lambda * crossprod(D)          # crossprod(D) 就是 t(D) %*% D

  w <- rep(1, L); z <- y
  for (i in seq_len(n_iter)) {
    W <- Diagonal(x = w)
    z <- as.numeric(solve(W + DtD, w * y))
    w <- p * (y > z) + (1 - p) * (y < z)
  }
  z
}

那一行 w <- p * (y > z) + (1 - p) * (y < z) 跟 Python 版一字不差—— R 的邏輯向量跟數字相乘時也會自動變成 0/1。

第 15 章　自己寫 Savitzky-Golay 與找峰

R 的 base 沒有現成的 SG 平滑，但它的原理很好寫：在每個點的鄰域套一條多項式， 取中心值。這等於用一組固定係數做卷積，而那組係數就是最小平方解的第一列。

sg_filter <- function(y, window = 9, poly = 3) {
  if (window %% 2 == 0) window <- window + 1
  half <- (window - 1) / 2
  t <- -half:half
  A <- outer(t, 0:poly, "^")               # Vandermonde 矩陣
  coef <- solve(crossprod(A), t(A))[1, ]   # 取第一列＝多項式在中心的值

  n <- length(y)
  ypad <- c(rep(y[1], half), y, rep(y[n], half))   # 兩端補值，避免邊界縮短
  as.numeric(stats::filter(ypad, rev(coef), sides = 2))[(half+1):(half+n)]
}

找峰也要自己寫。先找出「比左右鄰居都高」的候選點，再算 prominence：

peak_prominence <- function(y, i) {
  n <- length(y); h <- y[i]
  j <- i; lmin <- h
  while (j > 1) { j <- j - 1; if (y[j] > h) break; lmin <- min(lmin, y[j]) }
  j <- i; rmin <- h
  while (j < n) { j <- j + 1; if (y[j] > h) break; rmin <- min(rmin, y[j]) }
  h - max(lmin, rmin)          # 兩邊谷底取「較高」的那個
}

detect_peaks <- function(x, yc, noise, snr = 6, xmin = 150, xmax = 1800) {
  m <- which(x >= xmin & x <= xmax)
  xs <- x[m]; ys <- yc[m]; n <- length(ys)
  cand <- which(ys[2:(n-1)] > ys[1:(n-2)] & ys[2:(n-1)] >= ys[3:n]) + 1
  prom <- vapply(cand, function(i) peak_prominence(ys, i), numeric(1))
  keep <- prom >= snr * noise
  data.frame(position = xs[cand[keep]], height = ys[cand[keep]],
             prominence = prom[keep], snr = prom[keep] / noise)
}

第 16 章　自己實作 NNLS

base R 沒有 NNLS。經典演算法叫 Lawson–Hanson 主動集法，想法是：

1. 一開始假設所有係數都是 0；
2. 找出「最想變成正數」的那一個變數，把它加進「活躍集合」；
3. 只對活躍集合解一般的最小平方；
4. 如果解出來有負值，就沿著方向縮回去，把變負的踢出集合；
5. 重複，直到沒有變數想再變正。

nnls_fit <- function(A, b, max_iter = 300, tol = 1e-10) {
  n <- ncol(A); P <- logical(n); xx <- rep(0, n)
  w <- as.numeric(crossprod(A, b - A %*% xx))     # 梯度：誰最想變正
  it <- 0
  while (any(!P) && max(w[!P]) > tol && it < max_iter) {
    it <- it + 1
    j <- which(!P)[which.max(w[!P])]
    P[j] <- TRUE                                  # 加進活躍集合
    s <- rep(0, n); s[P] <- qr.solve(A[, P, drop = FALSE], b)
    while (min(s[P]) <= 0) {                      # 出現負值 → 往回縮
      neg <- P & (s <= 0)
      alpha <- min(xx[neg] / (xx[neg] - s[neg]))
      xx <- xx + alpha * (s - xx)
      P[P & abs(xx) < tol] <- FALSE
      s <- rep(0, n); if (!any(P)) break
      s[P] <- qr.solve(A[, P, drop = FALSE], b)
    }
    xx <- s
    w <- as.numeric(crossprod(A, b - A %*% xx))
  }
  pmax(xx, 0)
}

第 17 章　直接用打包好的 R 腳本

Rscript raman_pipeline.R spectra/UNK2_blue.csv -o 結果
Rscript raman_pipeline.R spectra/ -o 結果
Rscript raman_pipeline.R spectra/UNK3_blue.csv \
        --refs spectra/PB15_phthalo_blue.csv spectra/rutile_titanium_white.csv
Rscript raman_pipeline.R --list-library

或在 RStudio 裡當函式庫用：

source("raman_pipeline.R")

sp <- load_spectrum("spectra/cinnabar_1.csv")
pr <- preprocess(sp$x, sp$y)
pk <- detect_peaks(sp$x, pr$yc, pr$noise, snr = 6)
head(pk[order(-pk$prominence), ], 5)

第四篇　Orange Data Mining 路線（不寫程式）

第 18 章　Orange 與 orange-spectra 附加元件

Orange 是一套用「拉方塊、連線」來做資料分析的工具。 每個方塊叫一個 widget（元件），把它們串起來就是一條分析流程，完全不用寫程式。

orange-spectra 是專門處理光譜的附加元件，提供 11 個 widget。安裝方式二選一：

選單 Options → Add-ons → Add more...
輸入 orange-spectra → 勾選 → OK → 重新啟動 Orange

pip install orange-spectra
python -m Orange.canvas

需求：Python ≥ 3.9、Orange3。

第 19 章　載入光譜、疊圖看一眼

1. 拖出 Load Spectra Files，雙擊打開，選擇教材的 spectra/ 資料夾。 它會一次把裡面 10 個 CSV 全部讀進來，每一列是一張光譜。
2. 接上 Merge Spectra。因為這批檔案的波數格點不一樣 （標準品非等距、未知樣 1 或 2 cm⁻¹ 等距），Merge Spectra 會自動內插到共同範圍。 在 Normalization 選 max，讓每張譜的最高點都變成 1，形狀才好互相比較。
3. 接上 Orange 內建的 Data Table 看看資料長什麼樣，確認 10 列都在。

第 20 章　Peak Finder：自動找峰

Peak Finder 會輸出每個峰的位置、高度、FWHM、prominence 和面積—— 就是第 9 章那段 Python 程式做的事，只是不用寫。

要調的參數就是門檻。門檻設太低會抓到一堆雜訊，太高會漏掉弱帶。 把門檻從低往高慢慢拉，看峰的數量什麼時候穩定下來，那附近就是好設定。

第 21 章　相似度比對與建立譜庫

做法一：兩兩相似度

Spectra Similarity 提供四種指標。初學就先用 Pearson 相關係數： 數值 1 代表形狀完全一樣，0 代表無關。

做法二：建立譜庫再查詢

Spectral Library 有兩個輸入孔：Reference（參考譜）與 Query（未知樣）。 把標準品接到 Reference，未知樣接到 Query，它會輸出排名清單（Hits）與最佳匹配（Best Match）， 並可存成 .speclib 檔重複使用。

參考譜不夠？用 Import Spectrum URL 從 IRUG（Infrared and Raman Users Group） 或 SOPRANO 線上譜庫抓——這兩個都是文物保存領域公開的標準譜庫。

第 22 章　Mixture Analysis：不寫程式的 NNLS

這個 widget 做的事跟第 11 章的 Python 程式一模一樣：把未知樣拆成 mixture ≈ Σ cᵢ · refᵢ，強制係數非負，輸出組成百分比與 R² 擬合度， 以及 Fit data（擬合曲線）。

照著做一次：把 UNK3_blue 當 mixture， PB15_phthalo_blue 與 rutile_titanium_white 當 references。你會看到 R² 只有 0.4 左右。

第五篇　三個真實案例，以及程式會犯的錯

這三件都是現場採樣的藍色顏料，來自不同時間、不同儀器的歸檔資料。 同樣叫「藍色」，材料卻完全不同。

第 23 章　案例一 未知樣 A：酞菁藍 + 鈦白

資料：逗號分隔、300–3200 cm⁻¹、每 2 cm⁻¹ 一點，共 1451 點。

觀察到的峰：592 / 682 / 748 / 1004 / 1042 / 1126 / 1142 / 1340 / 1452 / 1528（最強） / 1598， 另外還有 444 與 610，以及 1718、2870、2928、3072。

推理：

1. 1528 那組十一個帶，與酞菁藍 PB15 的文獻帶位全部對上 → 藍色來源確定。
2. 444 與 610 是金紅石型 TiO₂。怎麼確定不是銳鈦礦？因為銳鈦礦的 396 與 516 完全沒有出現。 （銳鈦礦最強的 143 帶因為檔案從 300 開始所以看不到，但光憑 396/516 缺席就夠了。）
3. 1718 是酯類的 C=O，2870/2928 是 C–H，3072 是芳香族 C–H → 有機黏合劑。

第 24 章　案例二 未知樣 B：普魯士藍

決定性證據只有一個數字：2154 cm⁻¹，S/N 高達 202，是全譜最強的峰， 旁邊還有一個 2091。

黏合劑方面：1001 是又尖又強的單取代苯環呼吸振動，配上 1032 / 1449 / 1597 / 2937 / 3066， 聚苯乙烯的十二個特徵帶命中十個 → 苯乙烯系合成樹脂。而且這件完全沒有 TiO₂（447 缺席）。

第 25 章　案例三 未知樣 C：兩種酞菁的混合

這件最複雜，也最能示範「殘差分析」的威力。

第一個線索：主帶落在 1535。 酞菁藍 PB15 的主帶是 1527， 酞菁綠 PG7（氯化銅酞菁）是 1538。1535 卡在正中間——這是兩者混合的典型表現。

第二個線索：殘差。 用第 11 章的方法，拿「孔雀藍(PB15) + 鈦白(rutile)」去解混， 扣掉之後殘差裡還剩 687 / 777 / 813 / 1086 / 1218 / 1293 / 1539—— 正好是 PG7 的帶位，而且七個帶彼此自洽。

第三個線索：強度比。 776 與 747 兩帶的高度比，在這件是 0.52， 而純 PB15 標準品只有 0.26。多出來的 776 就是 PG7 貢獻的。

第 26 章　程式一定會犯的三種錯

本教材的腳本已經加了三道門檻和資料品質檢查，但它仍然會出錯。 你必須知道它會怎麼錯，才有資格用它的輸出寫報告。

錯誤型一：肩部假陽性

強帶的斜坡上，強度隨便都超過門檻。這是最常見的一種。 我們用「該位置必須本身是個峰」來擋，但相鄰很近的帶還是可能誤判。

怎麼查：把那個區間的原始數值印出來逐點看，確認峰心到底在哪裡。

錯誤型二：對著雜訊硬湊答案

拿 PB15_3_1（量測失敗那張）跑比對，早期版本的程式給出「赤鐵礦，主成分」。 那完全是幻覺——那張光譜根本沒有任何 S/N ≥ 20 的峰。

怎麼防：腳本現在會先檢查資料品質，沒有夠強的峰就印出警告並把所有判定降級。 你自己看報告時，第一件事就是看有沒有那個警告。

錯誤型三：庫外成分被最接近的那個吃掉

譜庫只有 15 種顏料。真實世界的現代顏料有好幾百種。 當樣品的成分不在庫裡，程式不會說「不知道」，它會給你最像的那個。

教材裡有兩個現成的例子：carmine 被判成石膏（次要成分）， 未知樣 C 被判出一個根本不存在的石膏——因為樹脂的苯環帶 1001 跟石膏的 1008 只差 7 cm⁻¹，落在容忍範圍內。

怎麼查：問自己「這個顏料的其他帶呢？」真石膏會有 414 / 493 / 619 / 1135 一起出現， 而且它的 1008 是又尖又窄的帶。只有一個帶對上，就不要相信。

延伸練習

1. 把 cinnabar_1.csv 和 cinnabar_2.csv 都跑一次， 比較兩次量測的峰位差幾個 cm⁻¹。這個差就是你的量測重複性。
2. 把 ALS 的 λ 從 10³ 掃到 10⁸，記錄硃砂 257 帶的高度怎麼變。畫成一張圖。
3. 教材的硃砂實測峰位一致比文獻高約 4 cm⁻¹，但同一台儀器的鈦白卻跟文獻吻合。 想想有哪些可能的原因，並設計一個實驗來區分它們。
   （提示：矽晶片在 520.7 cm⁻¹ 有一個非常穩定的帶，是拉曼的標準校正物。）
4. 用 Orange 的 PLS-DA，看看哪些波段最能把「含酞菁」和「不含酞菁」的樣品分開。 把它找出的重要波段，跟第 10 章譜庫裡 PB15 的帶位對照。

資料與文獻來源
顏料帶位：Burgio, L. & Clark, R. J. H. (2001). Library of FT-Raman spectra of pigments, minerals, pigment media and varnishes, and supplement to existing library of Raman spectra of pigments with visible excitation. Spectrochimica Acta Part A, 57, 1491–1521.
顏料年代：MFA CAMEO 材料資料庫（Phthalocyanine blue、Titanium dioxide 條目）。
Orange 附加元件：orange-spectra 0.5.0（PyPI）。
光譜資料：2026-08-10 標準品量測，以及 2022–2023 現場採樣歸檔資料。 教材中所有數字均為實際執行結果。