1. 项目概述：从数据噪声中挖掘物理方程

在化学工程、环境科学乃至生物制药领域，我们常常面对一个经典难题：手头有一堆从实验或模拟中得来的、带着噪声的数据，比如一个吸附柱出口的浓度随时间变化的“突破曲线”，我们如何从中反推出背后那个控制着物质吸附速率的精确数学模型？传统做法是，我们先根据物理直觉猜一个模型形式，比如朗缪尔等温线配合线性驱动力模型，然后把数据套进去拟合参数。但问题来了，如果我们的“直觉”不准，或者系统本身非常复杂（比如存在非线性、多机制耦合），预设的模型可能根本无法捕捉真实动态，导致预测外推时一塌糊涂。

这正是我最近深度折腾的一个方向： 利用稀疏回归和符号回归，直接从数据中“挖掘”出可解释的数学模型 。这不仅仅是另一个机器学习拟合任务，它的终极目标是 模型发现 ——让数据自己告诉我们它遵循什么样的数学规律。想象一下，你不需要告诉算法“请用一个指数函数加一个线性项来拟合”，而是对它说：“这里有一堆关于吸附量 q 、平衡吸附量 q* 和吸附速率 ∂q/∂t 的数据，你试试看，能不能找出它们之间最简洁、最准确的数学关系式？”

我手头的案例，正是围绕一个非线性的对流-扩散-吸附偏微分方程系统展开的。通过数值模拟，我生成了包含噪声的突破曲线数据，然后用人工神经网络去学习这个复杂系统的动态。但神经网络是个“黑箱”，预测虽准，却无法告诉我们内在的物理机制。这时，稀疏回归和符号回归就登场了，它们像两把手术刀，剖开神经网络的预测结果，试图提取出里面隐藏的、人类可读的数学表达式。这个过程，对于从单点、有噪声的实验数据中逆向识别吸附动力学模型，比如区分是经典的LDF模型还是更复杂的Vermeulen模型，具有极高的实用价值。

2. 核心思路与技术选型：为什么是稀疏与符号回归？

面对“从数据找方程”这个问题，可选路径其实不少。为什么最终聚焦在稀疏回归和符号回归上？这背后是一系列针对工程实际问题的考量。

2.1 问题本质与挑战

我们的核心数据是“突破曲线”，即吸附柱出口流体浓度随时间的变化。这个曲线背后，是一个描述床层内浓度分布的偏微分方程，其中最关键也最难以先验确定的，是描述吸附速率的源项 g(q, q*) 。传统方法要求我们预先指定 g 的函数形式（如 k(q* - q) ），但实际系统可能涉及表面扩散、孔内扩散、非线性平衡等多种机制的耦合，真实函数形式可能非常复杂，甚至包含奇异项（如 1/q ）或指数项。

更棘手的是，实验数据总是有噪声的，且通常我们只能获得有限位置（如出口）的数据。这就好比只通过听一首交响乐在音乐厅某个角落的回声，去反推整个乐队的乐谱和每个乐器的演奏法则，信息是高度压缩且模糊的。因此，我们需要的工具必须满足：

1. 抗噪声能力强 ：不能因为数据有轻微扰动就拟合出完全不同的方程。
2. 偏好简洁性 ：根据奥卡姆剃刀原理，在解释力相近的情况下，更简单的模型通常更接近真理，也更具泛化能力。
3. 无需强先验 ：最好能自动探索多种数学结构，而不是被限制在几种预设形式里。
4. 结果可解释 ：最终输出的应该是像 0.22*(q* - q) 这样的数学表达式，而不是一堆无法解读的神经网络权重。

2.2 稀疏回归：用“惩罚”寻找关键特征

稀疏回归，特别是结合了L1正则化的LASSO方法，其核心思想是 特征选择 。我们首先构造一个巨大的“特征库”。这个库包含所有可能相关的数学项，比如对于变量 q 和 q* ，我们可以构造： 1 , q , q* , q^2 , q*^2 , q*q* , q^3 , sin(q) , exp(q) 等等，理论上可以到很高的阶次和复杂的非线性组合。

然后，我们用线性模型 ∂q/∂t ≈ β_0 + β_1*q + β_2*q* + β_3*q^2 + ... 去拟合数据。L1正则化的妙处在于，它会在损失函数中加入一项惩罚项 λ * Σ|β_i| 。这个 λ 参数就像一把“剪刀”，它会将那些对预测贡献不大的特征对应的系数 β_i 直接压缩为零 。通过调整 λ ，我们可以控制模型的稀疏程度（即非零系数的个数）。

为什么选它？

• 自动化的模型简化 ：它自动完成了科学家手动尝试的“哪些项可以忽略”的过程。
• 缓解过拟合 ：在特征远多于样本数的情况下，强制稀疏性可以有效防止模型记住噪声。
• 可解释的输出 ：最终模型只保留少数几个项，如 -0.535 - 0.225*q + 0.234*q* ，物理意义清晰（常数项、与 q 负相关、与 q* 正相关）。

在实际操作中，我们需要通过交叉验证或信息准则（如贝叶斯信息准则BIC）来选取最优的 λ ，从而确定最佳的非零项数量。这个过程虽然需要调参，但逻辑清晰，计算相对高效。

2.3 符号回归：让算法自己“写公式”

如果说稀疏回归是在一个庞大的、预先定义好的“零件库”里挑选几个来组装，那么符号回归就是让算法自己从基本的数学“原子”（加、减、乘、除、乘方、指数、对数等）开始，从头发明组装方式。

它通常基于遗传编程等进化算法。算法随机生成一大批数学表达式树（例如，一个树根是加法，左子树是 0.22*q* ，右子树是乘法，乘法的子树是 0.22 和 q ，这样就代表了 0.22*q* - 0.22*q ），然后像自然选择一样：

1. 评估 ：计算每个表达式对数据的拟合误差（如均方误差）和复杂度。
2. 选择 ：保留那些拟合好且相对简单的表达式。
3. 交叉与变异 ：让这些“优秀”的表达式相互交换部分子树（交叉），或者随机改变某个节点（如把加法变乘法，把常数0.22变成0.23），产生新一代表达式。
4. 迭代 ：重复这个过程成千上万次，最终逼近最优解。

为什么选它？

• 结构发现 ：这是它最大的优势。它可能发现我们从未想到过的函数形式，比如 (q*^2 - q^2)/(2.0*q) 这样的结构，而这正是Vermeulen动力学模型的精确形式。
• 无需预设特征库 ：避免了人工构造特征库时可能存在的偏见或遗漏。
• 一步到位 ：直接得到完整的数学表达式，包括结构和参数。

它的缺点是计算成本通常比稀疏回归高，而且搜索空间巨大，可能陷入局部最优。但在本案例中，由于潜在的真实模型相对明确（是几种已知动力学模型的变体），符号回归展示出了强大的发现能力。

2.4 混合策略：神经网络作为“数据净化器”

你可能会问，为什么不直接对原始噪声数据做稀疏/符号回归？因为噪声会严重干扰方程结构的识别。这里，我们采用了一个巧妙的 混合策略 ：

1. 神经网络学习系统动态 ：首先，用一个物理信息神经网络或通用微分方程框架，去学习并逼近整个复杂的对流-扩散-吸附PDE系统。神经网络擅长在噪声中学习平滑的映射关系，其预测结果可以看作是对真实系统动态的一个“去噪”版本。
2. 对网络预测进行回归 ：然后，我们不再对原始噪声数据，而是对神经网络在大量空间-时间点上生成的、平滑的预测值 (q, q*, ∂q/∂t) 进行稀疏回归或符号回归。

这样做的好处是，我们用于方程发现的数据质量大大提高了，显著提升了回归方法找到正确数学结构的能力。神经网络充当了一个强大的“数据模拟器”或“降噪滤波器”。

3. 实操全流程：从数据生成到方程提取

下面，我将以朗缪尔等温线结合改进型LDF动力学这一具体案例，拆解整个操作流程。你可以把这个流程看作一个模板，应用于你自己的吸附或类似动力学系统。

3.1 第一步：生成高质量“在硅”数据集

一切始于数据。为了可控地验证方法，我们通常先进行数值模拟，生成“在硅”数据。

3.1.1 构建机理模型 我们使用一个标准的固定床吸附柱模型，核心PDE如下：

∂c/∂t = - ((1-ε)/ε) * g(q, q*) * τ_s - ∂c/∂x + (1/Pe) * ∂²c/∂x²
∂q/∂t = g(q, q*) * τ_s

其中， c 是液相浓度， q 是固相吸附量， q* 是由等温线方程（如朗缪尔： q* = Q_max * b*c / (1 + b*c) ）计算出的平衡吸附量。 g(q, q*) 就是我们想要发现的 动力学速率方程 。 ε 是床层孔隙率， τ_s 是空时， Pe 是佩克莱特数。

3.1.2 数值求解与加噪

• 空间离散 ：采用 正交配置有限元法 。将吸附柱长度划分为60个均匀的有限元，在每个单元内用三次Hermite多项式近似浓度分布。这种方法精度高，且能方便地处理边界条件。
• 时间积分 ：对离散后得到的大型刚性ODE系统，采用 变阶变步长的BDF方法 进行积分。我常用的绝对和相对容差设置为 1e-8 ，以确保数值解的精确性。
• 模拟场景 ：设定入口浓度阶梯变化，模拟吸附和脱附全过程，得到床层内每个位置每个时刻的 c, q, q* 。
• 提取与加噪 ：我们只“观察”柱出口（x=L）的浓度 c_out(t) ，这就是“突破曲线”。为了模拟真实实验，对此曲线加入一定比例（例如5%）的高斯白噪声。 这就是我们假设唯一能从实验获得的数据 。

实操心得 ：噪声水平的选择很重要。太小的噪声不真实，太大的噪声会给后续的神经网络学习和方程发现带来巨大困难。通常根据实际仪器测量精度来设定。生成数据时，务必保存下所有内部状态（全空间的 q, q*, ∂q/∂t ），作为验证回归结果的“金标准”。

3.2 第二步：训练神经网络学习系统动态

现在，我们只有一条带噪声的出口突破曲线 c_out(t) ，目标是反推出内部的动力学模型 g(q, q*) 。这里使用 物理信息神经网络 或 通用微分方程 框架。

3.2.1 网络结构与损失函数

• 输入 ：空间坐标 x 和时间 t 。
• 输出 ：浓度 c(x,t) 和吸附量 q(x,t) 的预测值。
• 物理约束 ：损失函数包含两部分：
  1. 数据损失 ：网络在出口处预测的 c(L,t) 与带噪声的观测数据之间的均方误差。
  2. 物理损失 ：将网络预测的 c, q 代入前述的PDE中，计算残差（即PDE左边减右边）。要求这个残差尽可能小。这就将物理定律作为软约束引入了学习过程。
• 训练 ：使用自适应矩估计优化器，并可能结合伴随敏感性分析来高效计算梯度。

3.2.2 获取学习后的动力学数据 训练完成后，神经网络已经学会了系统动态。我们在整个计算域（所有x和t）上查询训练好的网络，得到它预测的 q_pred(x,t) 和 c_pred(x,t) 。

• 根据 c_pred 和已知的等温线参数，计算 q*_pred 。
• 利用自动微分，计算时间导数 ∂q_pred/∂t 。 至此，我们获得了回归分析所需的高质量数据集： (q_pred, q*_pred, ∂q_pred/∂t) 。这个数据集平滑、连续，且隐含了真实的物理规律。

3.3 第三步：应用稀疏回归挖掘方程

现在进入核心环节。我们以改进型LDF动力学（真实模型为 g = 0.22*(q* + 0.2789*q* * exp(-q/(2*q* - q))) ）为例。

3.3.1 构建多项式特征库 我们假设动力学方程 g 可以用 q 和 q* 的多项式来近似。构建一个至三阶的完整多项式特征库：

特征项 = [1, q, q*, q², q*q*, q*², q³, q²*q*, q*q*², q*³]

这里一共10个特征。注意，我们没有加入 1/q 或 exp(q) 这类项，因为我们想测试多项式能否逼近一个本身包含指数项的复杂函数。

3.3.2 执行LASSO与超参数调优

1. 组织数据 ：将成千上万个数据点 (q, q*, ∂q/∂t) 整理成矩阵形式。 ∂q/∂t 是目标变量，特征矩阵的每一列对应一个特征项。
2. 标准化 ：对每个特征列进行标准化（减去均值，除以标准差），以防止量纲不同导致系数偏差。
3. 路径计算 ：使用LARS或坐标下降法，计算整个正则化参数 λ 路径上的解。随着 λ 增大，越来越多的特征系数会变为零。
4. 选择最优模型 ：使用 贝叶斯信息准则 作为选择标准。BIC同时考虑了模型的拟合优度（误差）和复杂度（非零参数个数），其值越小越好。我们遍历不同的 λ （即不同的非零项数量），计算每个模型的BIC。
   • 在案例中，对于改进型LDF，BIC在非零项数为4时最小。
5. 系数精炼 ：LASSO得到的系数可能会有轻微偏差。确定最优的4个特征项后，我们去掉L1正则化，仅对这些选中的项用普通最小二乘法或拟牛顿法（如BFGS）重新拟合一次，得到更精确的系数。

3.3.3 结果解读 最终，稀疏回归给出的方程为：

g_learned = -0.549 - 0.221*q + 0.278*q* - 0.000212*q*q*

对比真实模型，我们发现：

• 它抓住了核心驱动项： q* 的正系数和 q 的负系数，这与 (q* - q) 的直觉一致。
• 它包含了一个很小的交叉项 -0.000212*q*q* ，这可能是为了捕捉真实指数项带来的微小非线性修正。
• 虽然形式不同，但将这个多项式代入PDE进行预测，得到的突破曲线与原始数据（包括训练集和测试集）吻合得非常好（如图13所示）。这说明 在数据覆盖的范围内，这个简单的多项式是复杂指数函数的一个高度精确的代理模型 。

3.4 第四步：应用符号回归发现方程结构

并行地，我们将同一组神经网络输出的 (q, q*, ∂q/∂t) 数据喂给符号回归算法。

3.4.1 算法配置 使用基于遗传编程的 SymbolicRegression.jl 库。

• 运算符集 ： + , - , * , / , ^ （幂）。
• 函数集 ：这里我们尝试不提供指数、对数等函数，仅用基本算术和幂运算，看能否发现结构。
• 超参数 ：种群大小设为30，运行30代。算法会维护一个帕累托前沿，平衡表达式复杂度和拟合误差。

3.4.2 结果与对比 符号回归最终发现的表达式为：

g_learned_symbolic = -0.554 - 0.234*q + 0.281*q*

这个结果比稀疏回归的结果更简洁！它只有三个项，并且与经典LDF模型 k*(q* - q) 的形式惊人地相似（可以写成 0.281*q* - 0.234*q - 0.554 ）。虽然多了一个常数项，且系数与真实LDF的0.22有差异，但其 数学结构 被准确地发现了：吸附速率是 q* 和 q 的线性组合。

注意事项 ：符号回归的结果可能每次运行都有细微差别，因为进化算法具有随机性。通常需要多次运行，选择在帕累托前沿上最优且稳定的表达式。此外，结果的简洁性很大程度上依赖于所给的运算符集。如果加入 exp 函数，它有可能直接发现原始的指数形式。

4. 结果深度分析与工程启示

通过多个案例（朗缪尔/Sips等温线 × LDF/改进LDF/Vermeulen动力学）的测试，我们得到了一系列从数据中挖掘出的多项式方程。将它们汇总并与真实模型对比，能给我们带来更深的洞察。

4.1 案例结果汇总与对比

等温线模型	真实动力学模型	稀疏回归学习结果	符号回归学习结果	关键发现
朗缪尔	LDF: 0.22(q* - q)	-0.535 -0.225q +0.234q*	-0.554 -0.234q +0.281q*	结构高度还原 。两者都得到了 q* 和 q 的线性组合，与真实模型核心结构一致。常数项可视为对线性关系的微小调整。
朗缪尔	改进LDF: 0.22[q*+0.2789q*exp(...)]	-0.549 -0.221q +0.278q* -0.000212q q*	-0.554 -0.234q +0.281q*	稀疏回归发现了交叉项 q q* 以捕捉非线性；符号回归则用一个更简单的线性式逼近。两者预测效果俱佳。
朗缪尔	Vermeulen: 0.22(q*²-q²)/(2q)	-1.253 -0.429q +0.0034q² +0.506q* -0.0023q*² -0.0021q q*	-0.6098 +0.0122q +0.263q* -0.00526q q*	真实模型含 1/q 奇异项。稀疏回归用到了 q², q*², q q* 等高阶项逼近；符号回归用了 q q* 项。 多项式成功逼近了奇异函数 。
Sips	LDF: 0.22(q* - q)	-0.275 -0.210q +0.214q*	0.198q* - 0.200q	符号回归取得了完美胜利 ！它几乎精确地还原了真实模型的形式和参数，证明了其在简单线性结构发现上的强大能力。
Sips	改进LDF	-0.183 -0.215q -0.00029q² +0.272q*	0.277q* - 0.241q	与朗缪尔案例类似，符号回归找到了更简洁的线性近似。
Sips	Vermeulen	-0.309 -0.449q +0.00037q² +0.438q* -6.326e-5q*³ +6.214e-5q q*²	-0.003557q*² -0.216q +0.395q*	两者都用多项式逼近了复杂分式。符号回归发现了 q*² 项，这是一个有趣的结构性发现。

4.2 泰勒展开：为什么多项式能行？

一个核心疑问是：为什么一个简单的多项式（如 a + b*q + c*q* ）能够很好地逼近像 exp(-q/(2q*-q)) 或 (q*²-q²)/q 这样复杂的函数？ 泰勒展开 给出了完美的数学解释。

以Vermeulen模型 g = 0.22*(q*² - q²) / (2.0*q) 为例。它在某个工作点 (q*, q) = (50.13, 48.11) 附近进行二阶泰勒展开：

g ≈ 0.454 + 0.229q* - 0.229q + 0.002286q*² - 0.004764q*q + 0.002482q² + 高阶小量

看，展开后的形式正是一个关于 q 和 q* 的二次多项式！如果系统的 (q, q*) 在整个动态过程中变化范围不是特别大（通常吸附过程确实如此），那么这个多项式展开式在整个区间内都能对原函数给出极好的近似。图22直观展示了这一点：即使在远离展开点的地方，泰勒多项式（虚线）也与真实吸附速率曲线（实线）几乎重合。

工程启示 ：对于许多在合理操作范围内变化平滑的物理/化学过程，其本构关系（如动力学方程）完全可以用一个中低阶的多项式来高精度地代理。这为基于数据的模型简化提供了坚实的数学基础。

4.3 方法优势与适用场景总结

1. 从单点数据反推全局模型 ：这是本方法最大的价值。传统上，要识别动力学模型，可能需要设计多个不同条件的实验（如不同初始浓度、流速）。而本方法仅利用 单次实验 得到的单点突破曲线（且允许有噪声），结合物理约束的神经网络，就能反演出整个床层内部的动力学规律。
2. 超越黑箱，获得可解释模型 ：最终得到的是简洁的数学公式，工程师和科学家可以直观理解其物理意义（如“速率与推动力 (q*-q) 成正比”），并可用于传统的机理模型分析和流程模拟器中。
3. 对噪声鲁棒 ：通过神经网络进行初步学习，有效过滤了噪声，使后续的方程发现更加稳定可靠。
4. 两种回归方法互补 ：
   • 稀疏回归 ：更稳定，计算快，通过超参数调优可以控制模型复杂度。适合快速获取一个可靠的代理模型。
   • 符号回归 ：具有强大的结构发现能力，可能找到意想不到的、更简洁或更贴近原始机理的表达式。适合进行深入的机理探索和模型发现。

5. 避坑指南与实战技巧

在实际操作中，我踩过不少坑，也积累了一些确保成功的关键技巧。

5.1 数据准备与神经网络训练

• 数据量要足够 ：神经网络需要足够的数据来学习PDE系统。时间-空间网格不能太稀疏。在我的案例中，通常需要数千到上万个 (x, t) 数据点。
• 物理损失权重是关键 ：在训练PINN时，数据损失和物理损失之间的权重需要仔细调整。权重过大，网络可能难以拟合噪声数据；权重过小，物理约束可能不起作用。建议从1:1开始，根据验证集上的预测误差进行调整。
• 验证泛化能力 ：务必使用“测试集”——即模拟时预留的一部分时间区间或参数范围的数据，来检验训练好的神经网络和最终挖掘出的方程的预测能力。防止过拟合。

5.2 稀疏回归实操要点

• 特征库的构建要有物理直觉 ：虽然可以暴力生成高阶项，但结合对过程的物理理解来构建特征库，效率更高。例如，在吸附中， (q* - q) 作为推动力很可能是一个重要特征，可以显式加入库中。
• 标准化必不可少 ： q 和 q* 可能数值差异很大，务必进行标准化，否则L1正则化惩罚会不公平地偏向数值大的特征。
• 用BIC/AIC而不是只看拟合误差 ：随着特征项增加，拟合误差总会下降，但模型会变复杂。BIC/AIC能更好地平衡拟合与复杂度，避免过拟合。
• 检查系数稳定性 ：可以用Bootstrap抽样（从数据中有放回地多次采样并重复回归）来观察所选特征和其系数的稳定性。如果每次结果波动很大，说明模型可能不可靠。

5.3 符号回归调参心得

• 种群大小和迭代次数 ：不宜过小。对于中等复杂度问题，种群大小30-50，迭代50-100代是一个不错的起点。如果发现结果不稳定或未收敛，应增大这两个参数。
• 帕累托前沿分析 ：不要只看误差最小的那个表达式。分析帕累托前沿上的一系列解，选择一个在误差和复杂度之间取得良好平衡的模型。一个误差稍大但极其简洁的模型，往往比一个误差极小但极其复杂的模型更有用、更可靠。
• 运算符集的先验知识 ：如果你怀疑过程中存在指数衰减（如某些吸附动力学），可以把 exp 加入运算符集。如果怀疑有饱和效应，可以加入 /(1+...) 结构。合理的先验能极大缩小搜索空间，提高发现真实模型的概率。

5.4 结果验证与工程转化

• 必须进行数值验证 ：将挖掘出的方程代回原始的PDE求解器中，重新模拟突破曲线，与原始观测数据对比。这是检验模型有效性的 金标准 。如图13、15、17等所示，预测曲线与观测点的匹配度是最终的评判依据。
• 物理合理性检查 ：检查方程的物理意义。例如，吸附速率 ∂q/∂t 是否在合理范围内恒为正（吸附过程）？当 q 接近 q* 时，速率是否趋近于零？不符合物理常识的方程，即使拟合再好，也可能只是数据巧合。
• 用于机理洞察 ：比较从不同条件（如不同等温线）数据中挖掘出的方程。如果它们都具有 a*q* + b*q 的形式，只是系数不同，那就强有力地暗示了该系统的动力学可能普遍遵循线性驱动力模式，系数 a, b 可能与等温线参数有关。这本身就是一种深刻的机理发现。

最后想说的是，稀疏回归和符号回归并不是要完全取代传统的机理建模，而是成为了连接“数据海洋”与“物理定律”之间一座强大的新桥梁。当系统过于复杂、机理不甚明朗时，这套“数据驱动模型发现”的组合拳，能为工程师提供一个强有力的假设生成器。它给出的不是一个黑箱预测，而是一个白箱公式，你可以分析它、理解它、批判它，并在此基础上设计新的实验去验证它。这个过程，本身就是科学探索的乐趣所在。