Xử lý phổ Raman và mô hình NNLS bằng Python: Theo dõi quá trình kết tinh zeolite

Phổ Raman hiếm khi đi thẳng từ máy quang phổ vào mô hình. Trước khi trích xuất được thông tin định lượng, ta thường phải trải qua một loạt bước tiền xử lý, bao gồm cắt vùng phổ hữu ích, loại bỏ đỉnh nhiễu (spike), trừ nền và chuẩn hóa dữ liệu.

Trong hướng dẫn này, chúng ta sẽ đi qua toàn bộ quy trình đó bằng dữ liệu phổ Raman thực tế, lấy từ một nghiên cứu gần đây về theo dõi quá trình kết tinh zeolite. Điểm cốt lõi là phương pháp Bình phương tối thiểu không âm (NNLS) — một công cụ đơn giản nhưng mạnh mẽ để "tách" hỗn hợp phổ thành các thành phần riêng lẻ.

Bối cảnh: Quá trình chuyển pha giữa các zeolite

Trong tổng hợp zeolite, hiện tượng chuyển pha liên zeolite (interzeolite transformation) xảy ra khi một pha tinh thể này biến đổi thành một pha khác. Việc theo dõi tỷ lệ đóng góp của từng pha theo thời gian là rất hữu ích.

Phổ Raman có thể ghi nhận quá trình này, nhưng các phổ đo được là sự chồng chập của nhiều tín hiệu. Đây chính là lúc NNLS phát huy tác dụng: phổ đo được có thể được mô tả như một tổ hợp không âm của các phổ tham chiếu từ từng pha riêng lẻ.

Bài viết tập trung vào hai phần: xử lý dữ liệu phổ Raman và áp dụng NNLS, bắt đầu từ phổ thô và khôi phục lại diễn biến tương đối của hai pha faujasite (FAU) và gismondine (GIS) trong quá trình kết tinh.

Bộ dữ liệu thí nghiệm đầy đủ được công bố mở dưới định dạng NeXus/HDF5 trên Zenodo. Toàn bộ hướng dẫn cũng có sẵn dưới dạng Jupyter notebook trên GitHub của tác giả.

Bước 1: Nạp phổ Raman từ tệp NeXus

Định dạng NeXus rất thuận tiện cho loại thí nghiệm này vì phổ và siêu dữ liệu thí nghiệm (bao gồm cả thời gian kết tinh) được lưu cùng nhau. Sau khi tải tệp raman_crystallization_solid.nxs về máy, ta nạp thư viện và dữ liệu:

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"]

Tiếp theo, vòng lặp duyệt qua các nhóm dữ liệu và lưu mỗi phổ cùng thời điểm đo và trục bước sóng:

spectra = []
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"])
times_all = np.array([s["time"] for s in spectra])

Các phổ đo ở 0 giờ và 16 giờ không được dùng trong phân tích tiếp theo vì đặc trưng Raman của các pha zeolite tinh thể chưa đủ rõ. Ta bắt đầu phân tích từ mốc 18 giờ trở đi.

Minh họa quy trình xử lý dữ liệu phổ

Bước 2: Cắt phổ và nội suy về trục chung

Bước tiền xử lý đầu tiên là loại bỏ các vùng phổ không dùng đến. Với bộ dữ liệu này, ta giữ vùng 250–1250 cm⁻¹.

Vì các phổ không được ghi trên cùng một trục bước sóng tuyệt đối, sau khi cắt, mỗi phổ cần được nội suy về một trục chung. Điều này là bắt buộc trước khi gộp các phổ thành một ma trận dữ liệu duy nhất cho các bước tiếp theo.

Điểm quan trọng: cả việc cắt và nội suy đều được thực hiện trước khi loại đỉnh nhiễu và trừ nền, để các thuật toán đó chỉ thao tác trên phần tín hiệu thực sự liên quan và đã được đồng nhất.

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]

Bước 3: Loại bỏ đỉnh nhiễu (spike removal)

Phổ Raman thực tế thường xuất hiện các đỉnh nhiễu hẹp do tia vũ trụ hoặc nhiễu cảm biến. Thuật toán ở đây dựa trên độ nổi bật (prominence) và độ rộng (width) của đỉnh để phát hiện các đỉnh bất thường, sau đó thay thế vùng nhiễu bằng phép nội suy từ dữ liệu sạch hai bên.

Hai tham số cần tinh chỉnh trực quan là WIDTH_THRESHOLD và PROMINENCE_THRESHOLD, phụ thuộc vào độ phân giải phổ và độ rộng của các dải Raman thật.

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
])

Thuật toán sẽ bỏ qua các vùng nhiễu nằm sát biên hoặc những nơi cần ngoại suy, đồng thời đưa ra cảnh báo để người dùng kiểm tra lại.

Bước 4: Trừ nền bằng AsLS

Tiếp theo là xử lý nền phổ biến thiên chậm — một dạng nhiễu phổ biến. Ở đây ta dùng bình phương tối thiểu bất đối xứng (AsLS), được cài đặt trong thư viện pybaselines.

Ba tham số chính được giữ gần với phân tích gốc:

• lam = 3e3: độ mượt của đường nền
• p = 0.001: tham số bất đối xứng
• max_iter = 300: số vòng lặp tối đa

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)

Y_corrected = np.array(Y_corrected)

Phổ sau khi trừ nền

Bước 5: Tập trung vào vùng chẩn đoán FAU/GIS

Với phân tích NNLS, ta không cần toàn bộ phổ. Diễn biến hữu ích nhất nằm trong vùng hẹp chứa các dải khung đặc trưng của FAU và GIS. Do đó ta:

• Cắt về vùng 460–540 cm⁻¹
• Chuẩn hóa mỗi phổ bằng tổng cường độ của nó

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)

Trong vùng này, hai vạch đánh dấu khoảng 490 cm⁻¹ (đặc trưng FAU) và 509,5 cm⁻¹ (đặc trưng GIS) giúp quan sát rõ sự chuyển dịch tín hiệu.

Bước 6: Mô hình Bình phương tối thiểu không âm (NNLS)

Một phổ Raman chứa hai pha có thể được xấp xỉ như một tổ hợp tuyến tính của hai phổ tham chiếu:

y ≈ C_FAU · S_FAU + C_GIS · S_GIS

trong đó C_FAU và C_GIS là các hệ số tương ứng. NNLS ràng buộc các hệ số phải không âm, giải bài toán:

min với c ≥ 0 của ‖SC − y‖²₂

Điều này khiến NNLS rất phù hợp cho việc tách phổ đơn giản: mô hình nhanh, minh bạch và dễ diễn giải.

Chọn phổ tham chiếu

Ở đây, phổ đo tại 18 giờ được dùng làm tham chiếu FAU và phổ tại 144 giờ làm tham chiếu GIS:

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]

Kết quả

Đồ thị diễn biến cho thấy rõ sự phát triển của hai pha. Vào đầu quá trình kết tinh, tín hiệu Raman chủ yếu là FAU, nhưng tỷ lệ đóng góp của nó giảm dần khi GIS hình thành.

Diễn biến tỷ lệ hai pha FAU và GIS theo thời gian

Kết quả NNLS nắm bắt được quá trình chuyển pha FAU → GIS, chứng minh rằng phổ Raman có thể được dùng để theo dõi trực tiếp diễn biến pha ngay trong lúc tổng hợp.

Tổng kết quy trình

Toàn bộ luồng xử lý gồm các bước:

Phổ Raman thô → cắt phổ → loại đỉnh nhiễu → trừ nền → cắt vùng chẩn đoán → chuẩn hóa → NNLS

Điều thú vị là mô hình cuối cùng cực kỳ đơn giản. Một khi phổ đã được tiền xử lý đúng cách, NNLS chỉ cần các dấu hiệu phổ đại diện của từng pha để khôi phục diễn biến tương đối của chúng trong suốt quá trình.

Quan trọng hơn, quy trình này không giới hạn ở phản ứng tổng hợp cụ thể trên. Cùng một cách tiếp cận có thể áp dụng như một công nghệ phân tích quá trình (PAT) cho nhiều quy trình khác, nơi phổ đo được biến đổi như một tổ hợp của các thành phần đã biết — ví dụ chuyển pha, quá trình kết tinh, phản ứng hóa học hay thay đổi thành phần vật liệu.

Khi có sẵn phổ tham chiếu phù hợp, việc kết hợp phổ Raman tại chỗ với các phương pháp xử lý dữ liệu như NNLS mang lại một cách tương đối đơn giản để biến diễn biến phổ phức tạp thành thông tin có ý nghĩa về mặt hóa học.

Tài liệu tham khảo

• N. Coca-Lopez và cộng sự, "Process monitoring of zeolite crystallization: Bridging in situ Raman spectroscopy, ex situ XRD, and predictive machine learning models", Chemical Engineering Journal 543 (2026) 177865.
• Bộ dữ liệu mở trên Zenodo: https://doi.org/10.5281/zenodo.18972297
• N. Coca-Lopez, "An intuitive approach for spike removal in Raman spectra based on peaks' prominence and width", Analytica Chimica Acta 1295 (2024) 342312.
• N. Coca-Lopez và cộng sự, "Artificial Intelligence-Powered Raman Spectroscopy through Open Science and FAIR Principles", ACS Nano 19 (44): 38189–38218, 2025.