Raman Spectra Preprocessing and Non-Negative Least Squares Modelling in Python: Monitoring Zeolite Crystallization

# 如果需要：
# !pip install nexusformat scipy pybaselines matplotlib
import numpy as np
import matplotlib.pyplot as plt
from nexusformat.nexus import nxload
filename = "raman_crystallization_solid.nxs"
root = nxload(filename)
entry = root["entry"]

spectra = []
# 遍历 NeXus 数据组，将每条拉曼光谱与其采集时间和光谱轴一起存储：
for name in entry.keys():
    if name.startswith("data_"):
        d = entry[name]
        spectra.append({
            "name": name,
            "time": float(np.asarray(d["time_hours"])) if "time_hours" in d else np.nan,
            "x": np.asarray(d["wavenumber"], dtype=float),
            "y": np.asarray(d["intensity"], dtype=float),
        })
# 按结晶时间对光谱进行排序：
spectra_all = sorted(spectra, key=lambda s: s["time"])
# 将结晶时间存储为 NumPy 数组：
times_all = np.array([s["time"] for s in spectra])
print(f"光谱数量: {len(spectra_all)}")
print("时间 (h):", times_all)

plt.figure(figsize=(10, 5))
for s in spectra:
    plt.plot(s["x"], s["y"], label=f"{s['time']:g} h")
    plt.xlabel(r"拉曼位移 (cm$^{-1}$)")
    plt.ylabel("强度")
    plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
    plt.tight_layout()
    plt.savefig("figure.png", dpi=600, bbox_inches="tight")
    plt.show()

spectra = [s for s in spectra_all if s["time"] >= 18]
times = np.array([s["time"] for s in spectra])

# ----------------------------------------------------------
# 将所有光谱裁剪并插值到统一的拉曼位移轴上
# ----------------------------------------------------------
low_end = 250
high_end = 1250
# 找出所有光谱实际共有的光谱范围
x_min = max(low_end, max(np.min(s["x"]) for s in spectra))
x_max = min(high_end, min(np.max(s["x"]) for s in spectra))
# 使用第一个光谱作为参考轴
x0 = spectra[0]["x"]
mask_ref = (x0 >= x_min) & (x0 <= x_max)
x_crop = x0[mask_ref]
# 确保统一的轴从低到高拉曼位移排列
x_crop = np.sort(x_crop)
Y_raw = []
for s in spectra:
x = s["x"]
y = s["y"]
# 单独裁剪每个光谱
mask = (x >= x_min) & (x <= x_max)
x_i = x[mask]
y_i = y[mask]
# np.interp 需要递增的x
order = np.argsort(x_i)
x_i = x_i[order]
y_i = y_i[order]
# 插值到统一轴上
y_interp = np.interp(x_crop, x_i, y_i)
Y_raw.append(y_interp)
Y_raw = np.array(Y_raw)
# 绘图
plt.figure(figsize=(10, 5))
for t, y in zip(times, Y_raw):
plt.plot(x_crop, y, label=f"{t:g} h")
plt.xlabel(r"Raman shift (cm$^{-1}$)")
plt.ylabel("Intensity (a.u.)")
plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.savefig("figure_2.png", dpi=600, bbox_inches="tight")
plt.show()
print("共同拉曼范围:", x_crop.min(), "-", x_crop.max(), "cm-1")
print("共同点数:", len(x_crop))
print("Y_raw 形状:", Y_raw.shape)

from scipy.signal import find_peaks, peak_widths
from scipy import interpolate
import warnings
def spike_removal(
y,
width_threshold,
prominence_threshold,
interpolation_window=10,
width_param_rel=0.8,
interp_kind="linear",
remove_positive=True,
remove_negative=True,
):
"""使用峰显著度和宽度从拉曼光谱中去除尖峰。尖峰检测基于Coca-Lopez在《Analytica Chimica Acta 1295 (2024) 342312》中描述的显著度/FWHM方法。https://doi.org/10.1016/j.aca.2024.342312
显著度高于'prominence_threshold'且宽度低于'width_threshold'的峰被归类为尖峰，并通过插值替换为邻近未损坏点的值。
该实现通过将相同过程应用于-y来扩展原始方法，可选择检测负尖峰。它还使用保守检查，以避免在检测到的尖峰两侧没有干净数据时进行外推或插值。"""
y = np.asarray(y, dtype=float)
def detect_spikes(signal):
spike_mask = np.zeros(len(signal), dtype=bool)
peaks, _ = find_peaks(signal, prominence=prominence_threshold)
if len(peaks) == 0:
return spike_mask
widths = peak_widths(signal, peaks, rel_height=0.5)[0]
width_results = peak_widths(signal, peaks, rel_height=width_param_rel)
for width, left, right in zip(widths, width_results[2], width_results[3]):
if width < width_threshold:
start = max(int(np.floor(left)) - 1, 0)
end = min(int(np.ceil(right)) + 2, len(signal))
spike_mask[start:end] = True
return spike_mask
spikes = np.zeros(len(y), dtype=bool)
if remove_positive:
spikes |= detect_spikes(y)
if remove_negative:
spikes |= detect_spikes(-y)
if not np.any(spikes):
return y.copy()
y_out = y.copy()
padded = np.r_[False, spikes, False]
changes = np.diff(padded.astype(int))
starts = np.where(changes == 1)[0]
ends = np.where(changes == -1)[0] - 1
required = {"linear": 2, "quadratic": 3, "cubic": 4}
if interp_kind not in required:
raise ValueError("interp_kind必须为'linear'、'quadratic'或'cubic'。")
skipped = []
for start, end in zip(starts, ends):
window = np.arange(
max(start - interpolation_window, 0),
min(end + interpolation_window + 1, len(y)),
)
good = window[~spikes[window]]
bad = np.arange(start, end + 1)
if len(good) < required[interp_kind]:
skipped.append((start, end, "清洁点不足"))
continue
if not np.any(good < start) or not np.any(good > end):
skipped.append((start, end, "两侧均无清洁数据"))
continue
if bad[0] < good.min() or bad[-1] > good.max():
skipped.append((start, end, "需要外推"))
continue
try:
f = interpolate.interp1d(
good,
y[good],
kind=interp_kind,
bounds_error=True,
)
y_out[bad] = f(bad)
except ValueError:
skipped.append((start, end, "插值失败"))
if skipped:
warnings.warn(
f"{len(skipped)}个尖峰区域未被修正，因为插值被认为不安全: {skipped[:10]}"
)
return y_out

WIDTH_THRESHOLD = 2
PROMINENCE_THRESHOLD = 100
WINDOW = 5
Y_despiked = np.array([
spike_removal(
y,
width_threshold=WIDTH_THRESHOLD,
prominence_threshold=PROMINENCE_THRESHOLD,
interpolation_window=WINDOW,
)
for y in Y_raw
])
# 找出经过尖峰去除处理的光谱
affected = np.where(
np.any(~np.isclose(Y_raw, Y_despiked), axis=1)
)[0]
print(f"包含已修正尖峰的光谱数量: {len(affected)}")
print("时间:", times[affected])
# 绘制所有受影响的光谱
for i in affected:
plt.figure(figsize=(10, 4))
plt.plot(
x_crop,
Y_raw[i],
label=f"原始数据 — {times[i]:g} h",
alpha=0.7
)
plt.plot(
x_crop,
Y_despiked[i],
label="去尖峰处理",
lw=2
)
plt.xlabel(r"拉曼位移 (cm$^{-1}$)")
plt.ylabel("强度 (任意单位)")
plt.legend()
plt.tight_layout()
plt.savefig("figure_3_spike.png", dpi=600, bbox_inches="tight")
plt.show()

from pybaselines import Baseline
baseline_fitter = Baseline(x_data=x_crop)
baselines = []
Y_corrected = []
for y in Y_despiked:
baseline, _ = baseline_fitter.asls(
y,
lam=3e3,
p=0.001,
max_iter=300
)
baselines.append(baseline)
Y_corrected.append(y - baseline)
baselines = np.array(baselines)
Y_corrected = np.array(Y_corrected)
# 绘制所有光谱及其估计基线
for i in range(len(Y_despiked)):
plt.figure(figsize=(10, 4))
plt.plot(x_crop, Y_despiked[i], label=f"光谱 — {times[i]:g} h")
plt.plot(x_crop, baselines[i], label="估计基线", lw=2)
plt.xlabel(r"拉曼位移 (cm$^{-1}$)")
plt.ylabel("强度 (任意单位)")
plt.legend()
plt.tight_layout()
plt.savefig("figure_4_baseline.png", dpi=600, bbox_inches="tight")
plt.show()

plt.figure(figsize=(10, 5))
for t, y in zip(times, Y_corrected):
plt.plot(x_crop, y, label=f"{t:g} h")
plt.xlabel(r"拉曼位移 (cm$^{-1}$)")
plt.ylabel("基线校正后的强度 (任意单位)")
plt.tight_layout()
plt.savefig("figure_5_baselined.png", dpi=600, bbox_inches="tight")
plt.show()

mask_nnls = (x_crop >= 460) & (x_crop <= 540)
x_nnls = x_crop[mask_nnls]
Y_nnls = Y_corrected[:, mask_nnls].copy()
Y_nnls /= Y_nnls.sum(axis=1, keepdims=True)
cmap = plt.cm.plasma
plt.figure(figsize=(9, 5))
for i, (t, y) in enumerate(zip(times, Y_nnls)):
plt.plot(x_nnls, y, color=cmap(i / (len(times) - 1)), label=f"{t:g} h")
plt.axvline(490, color="magenta", alpha=0.35, ls="--")
plt.axvline(509.5, color="green", alpha=0.35, ls="--")
plt.xlabel(r"拉曼位移 (cm$^{-1}$)")
plt.ylabel("归一化强度")
plt.legend(bbox_to_anchor=(1.02, 1), loc="upper left")
plt.tight_layout()
plt.savefig("figure_spectra_evolution.png", dpi=600, bbox_inches="tight")
plt.show()

from scipy.optimize import nnls
FAU_TIME = 18
GIS_TIME = 144
i_fau = np.argmin(np.abs(times - FAU_TIME))
i_gis = np.argmin(np.abs(times - GIS_TIME))
R = np.column_stack([Y_nnls[i_fau], Y_nnls[i_gis]])
coef = np.array([nnls(R, y)[0] for y in Y_nnls])
# 从系数到加权分数：
fractions = 100 * coef / coef.sum(axis=1, keepdims=True)
fau_fraction = fractions[:, 0]
gis_fraction = fractions[:, 1]

plt.figure(figsize=(9, 4.5))
plt.plot(times, fau_fraction, "o-", color="green", label="FAU", lw=2, ms=8)
plt.plot(times, gis_fraction, "^-", color="purple", label="GIS", lw=2, ms=8)
plt.xlabel("Crystallization time (h)")
plt.ylabel("Relative Raman phase fraction (%)")
plt.legend()
plt.show()