在环境科学、水处理、海洋化学等领域，溶解性有机质（DOM）扮演着极其关键的角色，它影响着碳循环、污染物迁移转化以及水体生态系统的健康。然而，DOM成分复杂、来源多样，传统化学分析方法往往耗时费力且难以全面表征。你是否也曾为如何快速、有效地解析水样中DOM的来源与组成而困扰？

三维荧光光谱法结合平行因子分析，为我们提供了一把强大的“钥匙”。它不仅能对DOM进行“指纹识别”，还能通过数学方法剥离出混合信号中的独立组分，从而定量解析其来源（如陆源、微生物源、蛋白类等）。本文将围绕这一技术，从EEM图谱的解读到PARAFAC建模的完整流程，为你拆解一套可复现的实战方案。无论你是环境专业的学生、科研人员，还是从事水处理工艺研发的工程师，都能通过本文掌握从原始数据到科学结论的分析全链路。

1. 背景与核心概念：为什么是三维荧光光谱与PARAFAC？

在深入实操之前，我们有必要厘清几个核心概念，理解这项技术为何能成为DOM研究的利器。

1.1 溶解性有机质（DOM）及其复杂性

溶解性有机质是自然水体中一类重要的、能通过0.45微米滤膜的异质有机物混合物。它并非单一物质，而是包含腐殖酸、富里酸、蛋白质、氨基酸、碳水化合物等多种组分的“大杂烩”。这些组分来源广泛：

• 陆源（陆生）：主要由植物残体经微生物分解产生，如腐殖质，具有芳香性结构。
• 微生物源：由水体中藻类、细菌等生命活动产生，如蛋白质、酪氨酸、色氨酸。
• 人为源：来自生活污水、工业废水的排放，成分更为复杂。

传统方法如总有机碳（TOC）测定只能给出总量，无法区分来源；液相色谱等手段则前处理复杂、成本高。因此，我们需要一种能快速、无损、提供“指纹信息”的表征技术。

1.2 三维荧光光谱（EEM）与荧光峰区域

三维荧光光谱法通过扫描不同激发波长（Ex）和发射波长（Em）下的荧光强度，得到一个三维数据矩阵（即EEM图谱）。不同结构的DOM分子在特定Ex/Em波长对下会产生特征荧光峰。

国际上通常采用Coble（1996）提出的经典荧光区域积分（FRI）法，将EEM图谱划分为五个主要区域，每个区域对应一类典型的DOM组分：

1. 区域Ⅰ (Ex: 220-250 nm, Em: 280-330 nm)：酪氨酸类蛋白。
2. 区域Ⅱ (Ex: 220-250 nm, Em: 330-380 nm)：色氨酸类蛋白。
3. 区域Ⅲ (Ex: 220-250 nm, Em: 380-500 nm)：富里酸类物质。
4. 区域Ⅳ (Ex: 250-400 nm, Em: 280-380 nm)：可溶性微生物代谢产物。
5. 区域Ⅴ (Ex: 250-400 nm, Em: 380-500 nm)：腐殖酸类物质。

通过目视观察EEM图谱中峰的分布，可以初步判断DOM的主要来源。例如，样品若在区域Ⅳ和Ⅴ有强峰，可能指示较强的微生物活动或陆源输入。

1.3 平行因子分析（PARAFAC）的核心价值

然而，实际水样中的EEM信号是所有荧光组分叠加的结果，峰与峰之间经常重叠，仅靠目视分区积分误差大、主观性强。平行因子分析（PARAFAC）是一种多维数据分析模型，它能将三维数据阵列分解为一系列独立的荧光组分（因子）及其对应的相对浓度（得分）和激发/发射载荷谱。

简单来说，PARAFAC就像一位高超的“调音师”，能从一首混合交响乐（EEM数据）中，分离出小提琴、大提琴、长笛等每件乐器（独立荧光组分）的单独乐谱（载荷谱）和音量大小（相对浓度）。

它的核心优势在于：

• 数学分解，客观定量：减少主观判断，提供各组分的定量贡献。
• 解析重叠信号：有效分离光谱重叠的组分。
• 揭示潜在来源：每个解析出的组分可与已知的DOM类型（如类腐殖质C1、类色氨酸C2等）对应，从而更精确地溯源。

2. 环境准备与数据处理流程概述

在进行数据分析前，我们需要搭建合适的软件环境并理解完整的数据处理链条。本节将给出一个通用的、基于Python生态的解决方案。

2.1 核心工具与版本说明

本文的实战案例将主要使用Python语言，因其在科学计算和数据分析领域的强大生态。关键库及其作用如下：

• NumPy & Pandas: 用于数值计算和数据表格操作。
• Matplotlib & Seaborn: 用于绘制二维、三维图表。
• Scikit-learn: 用于数据预处理（如标准化）。
• DRUID (或类似工具): 用于PARAFAC建模的核心库。这里我们使用一个功能强大且社区活跃的库scikit-learn的扩展，或专门处理EEM数据的phee或eem包。为了演示通用性，我们将以scikit-learn结合自定义函数为例。在实际研究中，DRUID(DOM Fluorophore Identification by PARAFAC) 是常用工具。
• Jupyter Notebook/Lab: 交互式编程环境，非常适合数据探索和可视化。

版本建议：

Python >= 3.8 numpy >= 1.20 pandas >= 1.3 matplotlib >= 3.5 scikit-learn >= 1.0

请注意，PARAFAC的具体实现库可能需要单独安装，例如通过pip install phee。本文示例将侧重于方法流程，代码具有通用性，你需要根据实际使用的库调整API调用。

2.2 从仪器到模型：全流程概览

一个完整的分析流程通常包含以下步骤，理解它有助于我们后续分步实施：

1. 样品测量与数据导出：使用荧光分光光度计获取原始EEM数据。
2. 数据预处理：包括空白扣除、拉曼散射与瑞利散射校正、内滤效应校正（必要时）、荧光强度标准化（通常以拉曼单位表示）。
3. 构建三维数据阵：将多个样品的EEM数据组织成(样品数, 激发波长数, 发射波长数)的三维数组。
4. PARAFAC模型构建与验证：确定组分数，运行PARAFAC分解，并通过分裂半分析、残差分析等方法验证模型稳健性。
5. 结果解析与可视化：提取各组分的激发/发射载荷谱及样品得分，进行生物学或环境学解释。
6. 统计分析：将组分得分与环境因子（如pH、DOC、营养盐）进行相关性分析等。

3. 数据预处理：清洗与校正的关键步骤

原始EEM数据包含各种干扰信号，必须经过严格的预处理才能用于PARAFAC分析。这是保证结果可靠性的基石。

3.1 散射扣除：移除非目标信号

水分子和溶剂的拉曼散射以及溶剂的瑞利散射会在EEM图谱上形成明显的斜线或带状干扰，必须扣除。常见方法是将其影响区域的数据置为NaN（非数）或通过插值替代。

import numpy as np import pandas as pd def remove_scattering(eem_data, ex_wavelengths, em_wavelengths): """ 一个简单的散射区域掩码函数。 eem_data: 2D数组，单个样品的EEM矩阵 ex_wavelengths: 1D数组，激发波长 em_wavelengths: 1D数组，发射波长 """ # 创建与eem_data同样形状的布尔掩码，初始为False（保留区域） mask = np.zeros_like(eem_data, dtype=bool) # 示例：扣除一阶瑞利散射区域 (Ex = Em) for i, ex in enumerate(ex_wavelengths): for j, em in enumerate(em_wavelengths): if abs(em - ex) <= 10: # 设定一个带宽，例如±10 nm mask[i, j] = True # True表示该点需要被扣除 # 示例：扣除二阶瑞利散射区域 (Em = 2 * Ex) for i, ex in enumerate(ex_wavelengths): for j, em in enumerate(em_wavelengths): if abs(em - 2*ex) <= 10: mask[i, j] = True # 将散射区域的值设置为NaN eem_data_corrected = eem_data.copy() eem_data_corrected[mask] = np.nan # （可选）使用周围值进行插值填充NaN # 这里可以使用pandas或scipy的插值方法，为简化示例，我们仅返回NaN值 return eem_data_corrected # 假设我们有一个样品的数据 # sample_eem 是二维数组 # ex_vec, em_vec 是波长向量 # corrected_eem = remove_scattering(sample_eem, ex_vec, em_vec)

注意：实际应用中，散射扣除算法更复杂，可能涉及插值。许多专业软件（如Origin、MATLAB工具箱）或Python包（如eem）已内置成熟函数，建议优先使用。

3.2 空白扣除与标准化

每个样品都需要扣除超纯水空白样的EEM信号，以消除仪器背景和溶剂拉曼散射的影响。校正后的荧光强度通常用拉曼单位（R.U.）表示，即用样品荧光强度除以空白样在激发波长350 nm处、发射波长371 nm附近的拉曼散射峰面积进行归一化。这一步是为了使不同仪器、不同时间测量的数据具有可比性。

def raman_normalization(sample_eem, blank_eem, ex_index_350): """ 简化的拉曼归一化示例。 sample_eem: 扣除散射后的样品EEM blank_eem: 超纯水空白的EEM（同仪器条件） ex_index_350: 激发波长350 nm对应的索引号 """ # 1. 空白扣除 corrected_eem = sample_eem - blank_eem # 2. 计算空白在Ex=350 nm处的拉曼峰面积（近似） # 假设发射波长371 nm附近有数据点，找到其索引 # 这里需要根据你的实际波长向量计算 # em_index_371 = np.argmin(np.abs(em_wavelengths - 371)) # raman_area = np.trapz(blank_eem[ex_index_350, em_index_371-5:em_index_371+5], # x=em_wavelengths[em_index_371-5:em_index_371+5]) # 为了示例，我们假设已计算出raman_area raman_area = 100.0 # 这是一个示例值，实际需计算 # 3. 归一化到拉曼单位 normalized_eem = corrected_eem / raman_area return normalized_eem

4. PARAFAC建模实战：从数据到组分

预处理完成后，我们进入核心环节——使用PARAFAC模型解析DOM组分。我们将使用scikit-learn的decomposition模块中的PARAFAC实现（需注意，sklearn官方未直接提供，但可使用tensorly库，它是处理张量分解的强大工具）。这里以tensorly为例。

4.1 构建三维数据阵列

首先，将多个预处理后的样品EEM数据堆叠成一个三维数组X。

import tensorly as tl from tensorly.decomposition import parafac import matplotlib.pyplot as plt # 假设我们有3个预处理好的样品数据，每个都是 (num_ex, num_em) 的二维数组 # sample1, sample2, sample3 是已经过校正和归一化的EEM矩阵 sample_list = [sample1, sample2, sample3] # 你的实际数据列表 # 检查所有样品维度是否一致 num_samples = len(sample_list) num_ex, num_em = sample_list[0].shape # 构建三维数据阵: (样品数, 激发波长数, 发射波长数) X = np.zeros((num_samples, num_ex, num_em)) for i, sample in enumerate(sample_list): X[i, :, :] = sample print(f"三维数据阵X的形状: {X.shape}")

4.2 确定组分数与模型拟合

确定最佳组分数（n_components）是PARAFAC分析中最关键且最具挑战性的一步。常用方法包括：

• 分裂半分析（Split-half analysis）：将数据集随机分成两半，分别建模，比较两组分载荷谱的相似性。高度相似说明模型稳健。
• 残差分析：观察模型拟合残差是否随机分布。
• 核心一致性诊断（Core Consistency Diagnostic）：值接近100%表明模型有效。
• 通过解释方差：增加组分数对方差贡献的提升是否显著变小（类似PCA的碎石图）。

由于tensorly的parafac函数不直接提供核心一致性诊断，我们通常需要结合多种方法，或使用专门的DOMFluor工具箱（MATLAB）。以下演示如何用tensorly拟合并评估不同组分数。

# 尝试不同的组分数 possible_components = [2, 3, 4, 5] results = {} for rank in possible_components: # 执行PARAFAC分解 # init='random' 表示随机初始化，多次运行取最优可避免局部最优解 factors = parafac(X, rank=rank, init='random', random_state=42, n_iter_max=5000) # factors 是一个包含三个矩阵的元组: (A, B, C) # A: 样品得分矩阵 (n_samples, rank) -> 组分相对浓度 # B: 激发载荷矩阵 (n_ex, rank) -> 组分激发光谱 # C: 发射载荷矩阵 (n_em, rank) -> 组分发射光谱 # 计算模型重建的数据 X_reconstructed = tl.kruskal_to_tensor(factors) # 计算残差平方和 (SSE) 和解释方差 sse = np.sum((X - X_reconstructed) ** 2) tss = np.sum((X - np.mean(X)) ** 2) var_explained = 1 - (sse / tss) results[rank] = { 'factors': factors, 'sse': sse, 'var_explained': var_explained, 'X_reconstructed': X_reconstructed } print(f"组分数 {rank}: 解释方差 = {var_explained:.4f}")

4.3 模型验证与组分识别

拟合后，需要对模型进行验证。一个简单但重要的检查是可视化载荷谱。一个有效的荧光组分，其激发和发射载荷谱应该是平滑的、具有单一峰或合理形状的，不应出现多个不相关的峰或剧烈震荡。

def plot_component_loadings(factors, ex_wavelengths, em_wavelengths, component_idx): """ 绘制指定组分的激发和发射载荷谱。 factors: PARAFAC分解结果 component_idx: 要绘制的组分索引 (从0开始) """ A, B, C = factors # A:得分， B:激发载荷， C:发射载荷 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4)) # 激发载荷谱 (Excitation Loadings) ax1.plot(ex_wavelengths, B[:, component_idx], 'b-', linewidth=2) ax1.set_xlabel('激发波长 (nm)') ax1.set_ylabel('载荷强度') ax1.set_title(f'组分 {component_idx+1} - 激发光谱') ax1.grid(True, alpha=0.3) # 发射载荷谱 (Emission Loadings) ax2.plot(em_wavelengths, C[:, component_idx], 'r-', linewidth=2) ax2.set_xlabel('发射波长 (nm)') ax2.set_ylabel('载荷强度') ax2.set_title(f'组分 {component_idx+1} - 发射光谱') ax2.grid(True, alpha=0.3) plt.tight_layout() plt.show() # 假设我们选择3组分模型，并查看第一个组分 best_rank = 3 best_factors = results[best_rank]['factors'] plot_component_loadings(best_factors, ex_wavelengths, em_wavelengths, component_idx=0)

组分识别：将绘制出的载荷谱与文献中已报道的标准DOM荧光组分进行比对。例如：

• 组分C1(Ex/Em ~ 240(330)/420 nm): 常被认为是类紫外腐殖质。
• 组分C2(Ex/Em ~ 260(370)/470 nm): 类可见腐殖质。
• 组分C3(Ex/Em ~ 280/340 nm): 类色氨酸（蛋白类）。
• 组分C4(Ex/Em ~ 270/300 nm): 类酪氨酸（蛋白类）。
• 组分C5(Ex/Em ~ 320-340/400-420 nm): 类微生物代谢产物。

通过比对峰位置，可以为每个解析出的因子赋予环境意义。

5. 结果解析与可视化：赋予环境意义

得到验证后的模型和组分后，下一步是将数学结果转化为环境科学解释。

5.1 样品得分分析与来源解析

A矩阵（样品得分矩阵）中的每一列代表一个组分在所有样品中的相对浓度（荧光强度）。我们可以通过分析不同样品在各组分上的得分差异，来推断DOM来源的空间或时间变化。

# 提取得分矩阵 A, B, C = best_factors score_df = pd.DataFrame(A, columns=[f'C{i+1}' for i in range(best_rank)]) score_df['Sample_ID'] = [f'Sample_{i+1}' for i in range(num_samples)] # 假设的样品ID print(score_df) # 可视化得分（例如，用条形图展示某个样品中各组分贡献） sample_idx = 0 # 查看第一个样品 fig, ax = plt.subplots(figsize=(8,5)) components = [f'C{i+1}' for i in range(best_rank)] contributions = A[sample_idx, :] ax.bar(components, contributions, color=['skyblue', 'lightgreen', 'salmon']) ax.set_ylabel('相对荧光强度 (得分)') ax.set_title(f'样品 {score_df.loc[sample_idx, "Sample_ID"]} 的DOM组分贡献') ax.grid(axis='y', alpha=0.3) plt.show()

5.2 绘制组分在EEM等高线图上的位置

将解析出的组分峰位置（激发/发射载荷最大值）叠加在传统的EEM等高线图上，可以直观展示其分布。

def plot_eem_with_components(sample_eem, ex_wavelengths, em_wavelengths, factors, sample_idx=0): """ 绘制指定样品的EEM等高线图，并标记PARAFAC解析出的组分峰位置。 """ A, B, C = factors num_components = B.shape[1] fig, ax = plt.subplots(figsize=(10, 8)) # 绘制EEM等高线图 # 需要将二维数据转置以满足contourf的输入要求 (Em为X轴， Ex为Y轴) XX, YY = np.meshgrid(em_wavelengths, ex_wavelengths) contour = ax.contourf(XX, YY, sample_eem, levels=20, cmap='jet') plt.colorbar(contour, ax=ax, label='荧光强度 (R.U.)') # 标记每个组分的峰位置（激发和发射载荷最大值对应的波长） for comp in range(num_components): ex_max_idx = np.argmax(B[:, comp]) em_max_idx = np.argmax(C[:, comp]) ex_max = ex_wavelengths[ex_max_idx] em_max = em_wavelengths[em_max_idx] # 在图上标记点 ax.scatter(em_max, ex_max, s=150, edgecolors='white', facecolors='none', linewidths=2) # 添加文本标签 ax.text(em_max+5, ex_max+5, f'C{comp+1}', fontsize=12, color='white', fontweight='bold', bbox=dict(boxstyle="round,pad=0.3", facecolor='black', alpha=0.7)) ax.set_xlabel('发射波长 (nm)') ax.set_ylabel('激发波长 (nm)') ax.set_title(f'样品 {sample_idx+1} EEM图谱与PARAFAC解析组分') ax.set_xlim(em_wavelengths.min(), em_wavelengths.max()) ax.set_ylim(ex_wavelengths.min(), ex_wavelengths.max()) ax.grid(True, alpha=0.3) plt.tight_layout() plt.show() # 绘制第一个样品的EEM图并标记组分 plot_eem_with_components(X[0, :, :], ex_wavelengths, em_wavelengths, best_factors, sample_idx=0)

5.3 统计分析：关联环境因子

为了深化理解，可以将PARAFAC得到的组分得分（A矩阵）与同步测量的环境参数（如DOC浓度、pH、硝酸盐浓度等）进行相关性分析（如Pearson相关）或多元统计分析（如冗余分析RDA）。

import pandas as pd import seaborn as sns # 假设我们有一个环境因子DataFrame: env_df # 索引与样品顺序一致，列包括：'DOC', 'pH', 'NO3', 'Temperature'等 env_df = pd.DataFrame({ 'DOC': [5.2, 8.7, 3.1], # 示例数据 'pH': [7.1, 6.8, 7.5], 'NO3': [0.5, 2.1, 0.8] }, index=[f'Sample_{i+1}' for i in range(num_samples)]) # 合并得分和环境因子 analysis_df = score_df.set_index('Sample_ID').join(env_df) print(analysis_df) # 计算相关系数矩阵 corr_matrix = analysis_df.corr() print("\n相关系数矩阵 (PARAFAC得分 vs 环境因子):") print(corr_matrix) # 可视化相关系数热图 plt.figure(figsize=(10, 8)) sns.heatmap(corr_matrix, annot=True, cmap='coolwarm', center=0, square=True, linewidths=0.5, cbar_kws={"shrink": .8}) plt.title('PARAFAC组分得分与环境因子的相关性热图') plt.tight_layout() plt.show()

强正相关或负相关可以揭示DOM组分的潜在来源或转化过程。例如，类腐殖质组分（C1, C2）可能与DOC浓度正相关，而类蛋白组分（C3）可能与微生物活动指标正相关。

6. 常见问题、陷阱与排查思路

在实际操作中，从数据预处理到模型解释，每一步都可能遇到问题。下表汇总了常见问题及其解决方案。

7. 最佳实践与工程化建议

将三维荧光-PARAFAC分析从一次性的科研探索转化为可重复、可比较的常规分析工具，需要建立规范化的流程。

7.1 数据质量控制与标准化流程

1. 建立标准操作程序（SOP）：详细规定样品前处理、仪器参数设置（ slit width, scan speed）、空白和标准样测量频率。
2. 仪器性能监控：定期测量拉曼散射标准（如超纯水）以监控仪器稳定性，确保不同批次数据可比。
3. 统一的预处理脚本：将散射扣除、空白校正、拉曼归一化等步骤封装成可复用的脚本或函数，确保团队内所有成员处理方式一致。
4. 数据备份与版本管理：对原始数据、预处理后数据、模型参数和结果进行妥善备份，并使用版本控制（如Git）管理分析代码。

7.2 PARAFAC建模的稳健性策略

1. 多次初始化与模型选择：由于PARAFAC对初始值敏感，务必使用不同的随机种子运行多次（如50-100次），选择残差平方和（SSE）最小的解作为最终模型。
2. 强制非负约束：荧光强度不应为负。在调用分解函数时，确保启用非负约束（non_negative=True），这能使物理意义更明确且通常提高模型稳定性。
3. 系统性的模型验证：
   • 分裂半分析：将数据集随机分为两半，分别建模。比较两组分载荷谱的相似性（通过计算相关系数）。相关系数>0.95通常认为模型稳健。
   • 残差检查：模型残差应在整个EEM区域内随机、均匀分布，不应有规律性的结构。
   • 留一交叉验证：对于小样本数据集，可以依次剔除一个样品后建模，检验模型预测该样品的能力。

7.3 结果解释与报告的科学性

1. 谨慎进行组分识别：将解析出的组分与已发表的、使用相似仪器和预处理方法的文献中的组分进行比对，而不是简单地套用峰位范围。在文章中应明确说明比对依据。
2. 量化不确定性：报告组分得分时，可以尝试通过自助法（Bootstrapping）来估计其置信区间，从而评估结果的可靠性。
3. 结合多源数据：不要孤立地解释PARAFAC结果。务必与DOC、UV-Vis吸收光谱、稳定同位素等其他DOM表征数据，以及水文、化学、微生物等环境数据相结合，构建更完整的叙事。
4. 清晰可视化：在论文或报告中，应至少提供：
   • 所有样品的EEM等高线图（可选取代表性样品）。
   • 所有解析组分的激发和发射载荷谱图。
   • 组分得分的空间/时间变化图（如折线图、柱状图、空间插值图）。
   • 组分得分与环境因子的相关矩阵热图。

掌握三维荧光光谱结合PARAFAC分析，就相当于为DOM研究装备了高分辨率的“化学显微镜”。从繁琐的预处理代码编写，到谨慎的模型参数调优，再到严谨的环境学解释，每一步都需要耐心和细致。建议初学者从一个公开的小型数据集开始，完整复现整个流程，再应用到自己的研究数据中。随着经验的积累，你将能越来越熟练地运用这把利器，从复杂的水环境信号中，精准地解读出溶解性有机质的来源故事与转化历程。