从稀疏数据中学习控制方程的物理启发

Physics-informed learning of governing equations from scarce data

摘要：通过数据发现描述复杂物理系统的控制方程或规律，可以极大地促进各个科学和工程领域对这些系统的建模､仿真和理解｡本文提出了一种新颖的方法，称为稀疏回归的物理启发神经网络，用于从稀疏和噪声数据中发现非线性时空系统的控制偏微分方程｡特别地，该发现方法无缝集成了深度神经网络在丰富表示学习､物理嵌入､自动微分和稀疏回归方面的优势，以逼近系统变量的解，计算必要的导数，以及识别形成方程结构和显式表达式的关键导数项和参数｡该方法的有效性和鲁棒性在多种偏微分方程系统的数值和实验上得到了验证，这些系统考虑了不同的数据稀缺性和噪声水平，以及不同的初始 / 边界条件｡所得的计算框架显示了在大规模准确数据难以获取的实际应用中发现封闭形式模型的潜力｡

复杂动态系统建模的当前实践

目前对复杂动态系统建模的实践主要基于描述系统行为的常微分方程 (ODE) 和偏微分方程 (PDE)｡这些控制方程通常是通过严格的第一性原理 (如守恒定律) 或基于知识的现象学推导而获得的｡然而，仍然存在许多未被充分探索的现实复杂系统，其解析描述尚未发现，并且控制方程的简约闭合形式仍不清楚或部分未知｡幸运的是，观察数据集变得越来越丰富，提供了一种从数据中提取潜在方程的替代方案｡利用数据揭示控制规律或方程，可以显著推进和转变我们对各种科学和工程领域复杂物理系统的建模､仿真和理解｡例如，从观察数据 (如卫星遥感图像) 获取控制海冰演变的数学方程，为更好地理解和预测北极冰盖的生长､融化和移动提供了明显的好处｡从现场传感数据 (如多普勒雷达记录) 提炼出明确的公式将加速对天气和气候模式的更准确预测｡最近，机器学习理论､计算能力和数据可用性方面的进展激发了对基于数据的物理法则和控制方程发现的显著热情和努力｡

Bongard 和 Lipson 的开创性贡献利用分层符号回归和遗传编程成功地从数据中提取了控制非线性系统动态的潜在微分方程｡然而，这种优雅的方法在系统的维度上无法很好地扩展，计算成本高，并且可能遭遇过拟合问题｡最近，Brunton 等人的一项令人印象深刻的突破，提出了一种称为非线性动力学稀疏识别 (SINDy) 的创新性稀疏促进方法，通过稀疏回归从高维非线性函数空间中选择主导候选函数，以揭示简约的控制方程，特别是 ODE｡稀疏性是通过顺序阈值岭回归 (STRidge) 算法实现的，该算法递归确定受硬阈值约束的稀疏解｡这样的做法能够在所识别模型的复杂性和准确性之间取得平衡，从而实现简约｡近年来，SINDy 引起了极大的关注，导致了一系列变体算法的出现，用于识别以一阶 ODE 形式的投影低维代理模型，或通过线性嵌入用于广泛范围的非线性动态系统，如流体流动､结构系统､生物和化学系统､活性物质､非线性动态的预测控制､多时间尺度系统､捕食 - 猎物系统和随机过程等｡还有许多其他 SINDy 的扩展，通过发现隐式动态､纳入物理约束以及嵌入随机抽样来提高对稀疏发现高维动态的鲁棒性｡

SINDy 框架在偏微分方程的数据驱动发现方面的一个关键瓶颈在于对测量数据的质量和数量的强依赖，因为需要计算导数以构建控制方程｡特别地，使用有限差分或过滤来计算导数导致了一个关键挑战，降低了算法的鲁棒性｡这特别限制了 SINDy 在当前形式下应用于数据高度不完整､稀疏和嘈杂的场景｡值得注意的是，变分系统识别表明，基于同质分析计算导数的鲁棒性在发现 PDE 的弱形式时表现令人满意｡然而，这种方法在可用数据的保真度方面不易扩展｡另一个研究表明，弱形式能够显著提高对噪声的发现鲁棒性，但需要精心设计测试函数，这在高维时空系统中难以实现｡

自动微分被证明是解决上述问题的有效工具，这在物理启发神经网络 (PINN) 中成功应用于非线性 PDE 的正向和反向分析｡特别地，深度神经网络 (DNN) 用于逼近由 PDE 和少量可用数据约束的解｡PINN 在解决各种科学问题 (如流体流动､涡旋诱导振动､心血管系统等) 时引起了越来越多的关注，特别是在 PDE 的显式形式已知的情况下｡最近，Raissi 的一项重要工作介绍了一种深度隐藏物理模型，基于稀疏数据进行时空动态的数据驱动建模，其中未知的潜在物理通过可能的 PDE 项以弱约束的方式施加并由辅助神经网络隐式学习｡然而，所得到的模型仍然是一个 “黑箱”, 缺乏足够的可解释性，因为无法揭示闭合形式的控制方程｡最新研究显示，在预定义的 PDE 项库的约束搜索空间中，使用 DNN 和自动微分从嘈杂数据中获取闭合形式的 PDE 的潜力；然而，由于使用较不严格的稀疏回归和 DNN 训练而发生的假阳性识别也时有发生｡事实上，同时优化 DNN 参数和稀疏 PDE 系数，同时准确地施加稀疏性，是一个复杂的任务，并且在闭合形式 PDE 的发现中仍然是一项重要挑战｡

为此，我们利用这些进展，提出了一种新颖的 PINN-SR 方法 (即，带有稀疏回归的 PINN), 具有解释性和泛化性的显著特征，从稀疏和嘈杂的数据中发现非线性时空系统的控制 PDE｡我们的方法整合了 DNN 在丰富表示学习方面的优势､自动微分在准确导数计算方面的能力，以及 ℓ0 稀疏回归以应对现有方法在数据噪声和稀缺性方面的根本限制｡具体而言，本文涉及两个方法论贡献:

1.一个 “根 - 分支” 网络，由统一的潜在物理约束，能够处理来自不同初始 / 边界条件的少量多数据集｡

2.一种简单但有效的交替方向训练策略，用于异质参数 (即 DNN 可训练参数和稀疏 PDE 系数) 的优化｡

我们的研究在各种 PDE 系统上验证了该方法的有效性和鲁棒性，基于数值和实验数据集｡

结果

带稀疏回归的 PINN 用于 PDE 发现

我们考虑一个多维时空系统，其控制方程可以用以下形式的非线性耦合参数化 PDE 描述:

其中，u=u(x,t)∈R^1×n 是多维潜在解 维度 , 而 u_t 是一阶时间导数项；t∈[0,T] 表示时间，x∈Ω 指定空间；F[⋯] 是 u 及其空间导数的复杂非线性函数，由 λ 参数化；∇ 是相对于 x 的梯度算子；p=p(x,t) 是源项 (注意，在许多常见情况下，p=0 表示系统没有源输入)｡如果已知，PDE 也受初始和边界条件 (I/BCs) 的约束，表示为 I[x∈Ω,t=0;u,u_t]=0 和B[x∈∂Ω;u,▽_xu]=0。

我们的目标是从假设为不完整､稀疏和嘈杂的可用时空测量中找到 F[⋯] 的闭合形式｡我们假设物理定律仅由少数重要项控制，这些项可以从大型候选函数库中选择，其中可以应用稀疏回归｡该假设固有的将方程重构为以下形式 (为简化起见假设零或未知源):

u_t=ϕΛ

这里，ϕ=ϕ(u)∈R^1×s 是由许多候选项构成的广泛符号函数库，例如关于每个空间维度的常数､多项式和三角函数项，按行向量组合为 ϕ=1,u,u²,∇_xu,∇_yu,...,u³∇_xy,sin(u),... , 其中 ⊙ 表示逐元素的哈达玛乘积；表示库中的候选项总数；在 x,y,z 的上下文中，x,y,z 的下标表示导数；Λ∈R^s×n 是稀疏系数矩阵 (ϕ 中仅有活动候选项的值为非零), 例如 Λ=[λ_u,λ_v,λ_w]∈R^s×3 对于 u=u,v,w ｡如果有未知的源输入，源项的候选函数 p 也可以纳入 ϕ 中以进行发现 (见补充说明)｡因此，发现问题可以表述为：给定时空测量数据 Du , 寻找稀疏的Λ 使得方程 (2) 成立｡

我们提出一种新的 PINN-SR 范式，同时对系统响应建模并识别控制 PDE 的简约闭合形式｡该方法的创新算法架构如图 1 所示，其中考虑了来自两个不同 I/BC 场景的采样数据集：(1) 来自单一 I/BC 的一个数据集；(2) 来自多个 I/BC 的 r≥2 个独立数据集｡对于单一数据集的情况，我们通过 DNN (记为 N) 对潜在解 u 进行解释，即 u_θ=u(x,t;θ) , 其中 θ 表示包括权重和偏差的 DNN 可训练参数，如图 1a 所示｡当有多个独立数据集可用时，设计了一种 “根 - 分支” DNN, 如图 1b 所示，用于逼近与不同 I/BC 相关的潜在解 u_i(i=1,...,r) ), 即 u_θi=u(x,t;θ⁽⁰⁾,θ⁽ⁱ⁾) , 其中 θ⁽⁰⁾ 和 θ⁽ⁱ⁾ 分别表示根层 N⁽⁰⁾ 和分支层 N⁽ⁱ⁾ 的可训练参数｡值得注意的是，I/BC 不必事先已知或测量，因为测量数据已经反映了特定的 I/BC (例如，I/BC 和 PDE 解之间存在一一映射)｡DNN 本质上作为非线性函数，用于逼近潜在解及数据损失函数 L_d(θ;D_u)

图. 1: PINN-SR 数据驱动 PDE 发现框架的示意架构｡a. 来自单一 I/BC 的一个数据集的网络，b. 来自多个 I/BC 的 r≥2 个独立数据集的 “根 - 分支” 网络，c. 基于交替方向优化的网络训练示意图｡网络由两个组成部分构成：一个由可训练参数 θ 控制的 DNN, 将时空坐标 x,t 映射到潜在解 u=u,v,w , 以及由一组非线性 PDE 描述的物理定律，这些 PDE 由以未知稀疏系数 Λ 参数化的导数候选函数 ϕ 形成｡注意，对于多个独立数据集的情况，库 ϕ(i) 被串联以构建用于构建统一控制 PDE 的 ϕc ｡总损失函数 L (θ, Λ, Du, Dc) 由数据损失 Ld (θ, Du)､物理损失 αLp (θ, Λ, Dc) 和促进稀疏性的 ℓ0 正则化项 β||Λ||0 组成｡这里，α a 和 β 表示损失函数的相对权重，而 Du 和 Dc 分别代表测量数据和聚合样本｡注意，物理损失以残差形式仅在时空聚合样本上进行评估｡右侧稀疏系数矩阵 (或向量) 中的有色点表示非零值｡未知参数 {θ, Λ} 的同时优化导致了用于推断数据驱动全场解的训练 DNN 和发现的简约闭合 PDE｡

通过自动微分，在机器精度下评估 u 的导数，可以从 DNN 计算候选函数库 ϕ_θ ｡在多个独立数据集的情况下，来自分支网络的库 ϕ⁽ⁱ⁾ 被串联，以构建统一的 来构建统一的控制 PDE｡因此，重构 PDE 的稀疏表示可以用残差形式表示，即 Rθ:=u_θ^t−ϕ_θΛ=0 , 其中 R_θ∈R^1×n 表示 PDE 残差｡基本概念是同时调整 DNN 可训练参数 θ 和 PDE 系数 Λ , 使得神经网络可以拟合测量数据，同时满足由基础 PDE 定义的约束｡PDE 残差将在大量在时空空间中随机抽样的聚合点 Dc={xi,ti}_i=1^Nc 上评估，形成残差物理损失函数 L_p(θ,Λ,D_c) ｡当考虑多个 I/BC 时，计算数据损失和物理损失时，测量数据和聚合点将叠加 (基于统一的物理残差公式 R_θ→0 )。

PINN-SR 网络的总损失函数因此由数据损失 Ld ､残差物理损失 Lp 和正则化项组成，表示为:

其中，α 是残差物理损失函数的相对权重；β 是正则化参数；表示 ℓ0 范数｡优化总损失函数可以产生一个 DNN, 它不仅能够预测基于数据的全场系统响应，还能够揭示简约的闭合 PDE, 即 f _θ˙∗,Λ^∗:=argmin_θ,ΛL(θ,Λ,Du,Dc) , 其中 θ^∗,Λ^∗ 表示最优参数集合｡值得注意的是，总损失函数具有隐式复杂形式，因此直接求解优化问题是高度不可行的，因为 ℓ0 正则化使这个问题变成 NP 难题｡为了解决这个挑战，我们提出了一种交替方向优化 (ADO) 算法，将整体优化问题分解为一组可处理的子问题，以便在几次交替迭代中顺序优化可训练参数，如图 1c 所示｡在运行 ADO 算法进行发现之前，PINN-SR 的预训练通过简单地用 ∥Λ∥1 替换 Eq. (3) 中的 ∥Λ∥0 来进行，其中对 θ 和 的强制梯度优化变得适用｡ℓ1 正则化的预训练可以通过提供一个可接受的 “初始猜测” 来加速 ADO 的收敛｡更详细的公式和算法描述见方法部分和补充说明 1｡

DNN 和稀疏回归的协同作用导致以下结果：DNN 提供潜在解､其导数和候选函数项的准确建模，作为构建控制 PDE 的基础，而稀疏表示的 PDE 反过来约束了 DNN 的建模并投影出正确的候选函数，最终将测量系统转变为闭合形式的 PDE｡

单一数据集的基准 PDE 发现

我们观察到我们的方法在一组经典 PDE 上的有效性和鲁棒性，这些 PDE 用于表示具有非线性､周期性和 / 或混沌行为的广泛物理系统｡特别地，我们从由多个传感器在固定位置记录的稀疏和嘈杂时间序列测量中发现了 Burgers 方程､Kuramoto–Sivashinsky (KS) 方程､非线性薛定谔方程､Navier–Stokes (NS) 方程和 λ−(1) 反应 - 扩散 (RD) 方程的闭合形式 (数据受高斯白噪声污染), 这些数据来自单一 I/BC｡结果在表 1､图 2 和补充说明 2.1 中展示，表明该方法在数据稀缺性和噪声方面表现出良好的发现精度和鲁棒性｡

表 1: PINN-SR 发现结果的总结，包括多种基准模型的准确性｡

图. 2: 针对稀疏采样测量数据 (10% 噪声) 发现的选定基准 PDE｡a. 发现的 Burgers’ 方程：用于构建 PDE 的 16 个候选函数 ϕ∈R16 的稀疏系数 Λ 的演变，其中颜色代表系数值｡b. 发现的 KS 方程：36 个候选函数 ϕ∈R36 的稀疏系数 Λ 的演变｡c. 发现的非线性薛定谔方程：候选函数 ϕ∈R40 的稀疏系数 Λ 的演变｡d. 发现的 NS 方程：60 个候选函数 ϕ∈R60 的稀疏系数 Λ 的演变｡e. 发现的 RD 方程：用于重建 u 方程和 v 方程的 110 个候选函数 ϕ∈R110 的稀疏系数 λu∈R110 和 λv∈R110 的演变｡

Burgers 方程：我们首先考虑一个由 1D 粘性 Burgers 方程表示的耗散系统，其动态由

ut=−uux+νuxx 控制，其中 ν (等于 0.1) 表示扩散系数｡该方程描述了在有限时间后系统的衰减稳态粘性冲击，通常出现在简化流体力学､非线性声学和气体动力学中｡我们在 Burgers 方程的解中测试 PINN-SR 方法，该解在高斯初始条件下记录的行进冲击波中得到｡在特定情况下，10 个传感器随机放置在 256 个空间网格中的固定位置，并记录波形 101 个时间步长，导致用于参考文献 6 的数据集的 3.19%｡数据集的完整描述､候选函数库的设计 (16 项) 和模型训练见补充说明 2.1.1｡图 2a 显示了在 10% 噪声数据集上的 Burgers 方程发现｡系数 Λ∈R16×1 的演变展示了对真实值的鲁棒收敛 (误差约 0.88%), 导致精确发现｡训练后的 PINN-SR 能够正确重现来自嘈杂测量的动态响应 (例如，完整场的 ℓ2 预测误差为 1.32%), 如补充图 1 所示｡ADO 算法在第一次交替迭代后收敛，并显示出恢复 PDE 的正确稀疏模式的能力｡我们还在具有 10% 噪声的稀疏 u 测量中发现了 Burgers 方程，具有未知 / 未测量源，在发现潜在控制方程时，源项应该同时考虑和重建｡在这种情况下，我们将 14 个源候选函数 (由 和它们的组合组成) 纳入上述库，总计 30 个候选项，以同时发现 PDE 和重建未知源｡相应的发现结果总结在补充图 12 中，包括发现的方程和源函数､稀疏系数 Λ∈R30×1 的演变和预测的全场响应｡结果表明 PDE 和源项及其系数均被良好识别｡

Kuramoto–Sivashinsky (KS) 方程：考虑另一个具有固有不稳定性的耗散系统，由 1D KuramotoSivashinsky (KS) 方程控制 ut=−uux−uxx−uxxxx , 其中反扩散项 −uxx 导致破坏性行为，而四阶导数 uxxxx 产生混沌模式，使其成为理想的方程发现测试问题｡KS 方程广泛用于建模层流火焰前沿和耗散捕获离子模式等不稳定性｡我们随机选择 320 个点作为固定传感器，并记录波响应 101 个时间步长，导致参考文献 6 使用数据集的 12.6%｡共使用 36 个候选函数来构建潜在 PDE｡该示例的详细描述见补充说明 2.1.2｡值得注意的是，混沌行为在从测量不良的数据中逼近全场时空导数方面提出了显著挑战，尤其是高阶导数 uxxxx ｡现有方法 (例如 SINDy 方法系列) 在这种情况下由于测量非常粗糙和嘈杂而最终失败｡然而，PINN-SR 成功从具有 10% 噪声的子采样稀疏数据中提炼出 KS 方程的闭合形式，如图 2b 所示｡系数 Λ∈R³⁶×1 的演变表明，候选项及其相应系数均被正确识别 (接近原始参数；误差约为 0.94%) ADO 迭代内完成｡训练后的 PINN-SR 预测的全场波形也与精确解相吻合，相对 ℓ2 误差为 2.14%(补充图 2)｡

非线性薛定谔方程：在第三个示例中，我们发现非线性薛定谔方程 iut=−0.5u_xx−∣u∣²u , 其中 u 是一个复场变量｡这个著名的方程广泛用于建模光在非线性光纤中的传播､玻色 - 爱因斯坦凝聚､热等离子体中的朗缪尔波等｡我们从数据集中随机选择 37.5% 的子样本 (例如，从空间网格中随机选择), 如表 1 所示，用 40 个候选函数构建 PDE｡由于该函数是复值的，我们分别对 DNN 输出中的解的实部 (uR) 和虚部 (uI) 进行建模，将它们组合以获得复解 u=uR+iuI , 并构建用于 PDE 发现的复值候选函数｡为避免优化中的复杂梯度问题，我们在残差物理损失 Lp 中使用模值 而不是在公式 (5) 中显示的 ℓ2 范数 (有关更多细节，见补充说明 2.1.3)｡图 2c 显示了在 10% 噪声情况下发现的薛定谔方程｡稀疏系数 Λ∈R40×1 的演变历史清楚地表明 (图 2c; 误差约为 0.08%), 导致 PDE 的闭合形式被准确识别，同时重建的全场响应在实部和虚部上都很好地与精确解匹配，相对 ℓ2 误差为 0.26%(补充图 3)｡

Navier-Stokes (NS) 方程：我们考虑一个通过圆柱体的 2D 流体流动，其局部旋转动态由著名的 Navier-Stokes 旋涡方程控制 wt=−(u⋅∇)w+ν∇2w , 其中 w 是时空变化的涡度，u=u,v 表示流体速度，ν 是运动粘度 (在雷诺数 100 时 ν=0.01 )｡我们利用开放模拟数据和子采样数据集在指示区域内随机选择 500 个空间位置的流动响应 u,v,w , 记录 60 个时间步长｡所得数据集仅占参考文献 6 中数据集的 10%｡该示例的全面讨论见补充说明 2.1.4｡图 2d 总结了在 10% 噪声数据集上发现的 NS 方程的结果｡令人鼓舞的是，即使在 10% 噪声污染下，发现的 PDE 表达式也几乎与真实值相同，导数项及其系数均正确｡对应于 60 个候选函数 ϕ∈R1×60 的系数 Λ∈R60×1 在第一轮 ADO 迭代后迅速收敛到正确值 (图 2d)｡涡度模式和幅度也得到很好的预测，如补充图 5 所示 (所有快照的全场 ℓ2 误差约为 2.58%)｡这个例子为所提出的 PINN-SR 方法提供了一个令人信服的测试案例，能够在稀疏和嘈杂的数据中发现闭合形式的 NS 方程｡

反应 - 扩散 (RD) 方程：上述示例大多是低维模型，复杂性有限｡我们这里考虑一个 λ−(1) 反应 - 扩散 (RD) 系统，在 2D 域中由两个耦合 PDE 控制: ut=0.1∇²u+λ(g)u−ω(g)v 和 vt=0.1∇²v+ω(g)u+λ(g)v , 其中 u 和 是两个场变量，g=u²+v² , ω=−g² , λ=1−g² RD 方程表现出广泛的行为，包括波动现象和化学与生物系统中发现的自组织模式｡考虑的特定 RD 方程显示了在周期性边界条件下的螺旋波｡关于数据集､候选函数选择和 PINN-SR 模型的超参数设置的详细信息见补充说明 2.1.5｡图 2e 显示了在 10% 噪声情况下 110 个候选函数 ϕ∈R^1×110 的稀疏系数 λ_u,λ_v∈R^110×1 的演变｡稀疏项和相关系数都被准确识别，以形成闭合方程 (如图 2e 所示)｡由于 PDE 的复杂性和高维性，在 ADO 中需要稍微更多的训练周期以保持可靠的收敛｡训练后的 PINN-SR 在补充图 6 中的预测响应快照 (例如，在 t=2.95 时) 与真实值接近｡这个例子特别展示了我们方法在高维系统中从高度嘈杂数据中发现控制 PDE 的能力和鲁棒性｡

多个独立数据集的 PDE 发现

为了展示图 1b 中提出的 “根 - 分支” 网络在基于多个独立数据集发现 PDE 时的有效性，我们考虑 (1) 具有冲击行为的 1D Burgers 方程，以及 (2) 描述外部刺激激发的激活 - 抑制神经活动的 2D Fitzhugh– Nagumo (FN) 型反应 - 扩散系统｡测量数据在三种不同 I/BC 下以 10% 的噪声稀疏采样 (例如，时间序列或快照)｡注意，由于测量已经反映了特定的 I/BC, 因此 I/BC 不必事先已知或测量，这与系统响应之间存在唯一的一对一映射｡发现结果如下所述｡

带冲击行为的 Burgers 方程：在此示例中，我们测试了之前讨论的 Burgers 方程，其扩散 / 粘度参数较小 ( (ν=0.01/π≈0.0032) ), 基于通过施加三种不同 I/BC 生成的数据集｡如此小的系数会在具有陡峭梯度的紧凑区域内产生冲击形成 (见图 3c), 这可能会挑战 DNN 的逼近能力，从而影响发现｡三组初始和 Dirichlet 边界条件包括:

其中 G 表示高斯函数｡尽管不同 I/BC 的测量数据集表现出完全不同的系统响应，但它们遵循相同的潜在 PDE, 即 ut=−uux+0.0032u_xx ｡对于所有 I/BC, 我们假设在空间 ( (x∈[−1,1]) ) 中随机部署 30 个传感器测量波传播 (例如 u) 500 个时间瞬时 ( (t∈[0,1]) )｡在此，需要与前一个 Burgers 示例相比更密集的传感器网格，以捕获冲击行为｡图 3a 显示了在 10% 噪声下六个典型传感器记录的一些测量｡采用图 1b 中的三分支网络 ( (r=3) ) 进行发现｡数据集､候选函数库 (16 项) 和模型训练的完整描述见补充说明 2.3.1｡图 3b 描述了候选函数的系数演变 ( (Λ∈R^16×1) ), 其中库中正确的项 ( uu_x 和u_xx ) 被成功提取，而其他冗余项被 ADO 逐渐消除 (例如，几乎阈值到零)｡活动项的系数也被准确识别 (特别是导致冲击形成的小粘度参数，例如 0.0039)｡发现的 PDE 为 ut=−1.002uu_x+0.0032u_xx 。图 3c 和 3d 显示了三种 I/BC 情况的预测响应和误差，堆叠的全场 ℓ2 误差为 0.65%｡

图. 3: 基于在 10% 噪声下采样的三个 I/BC 数据集发现的 Burgers’ 方程 (小粘度)｡a. 三个数据集的嘈杂测量可视化｡注意，虽然有 30 个传感器，但仅在此图中显示了一部分｡b. 用于构建 PDE 的 16 个候选函数 ϕ∈R16 的稀疏系数 Λ∈R¹⁶ 的演变，其中颜色表示系数值｡正确的项 (uu_x 和 u_xx) 及其系数成功识别，而其他冗余项通过 ADO 被消除｡c, d. 三个 I/BC 情况下的预测响应和误差｡由于可视化几乎与预测无法区分，因此未列出真实值 (见补充图 7)｡堆叠预测的全场 ℓ2 相对误差为 0.65%｡

FitzHugh–Nagumo (FN) 反应 - 扩散系统：我们考虑一个 FN 型反应 - 扩散系统，在具有周期性边界条件的 2D 域 Ω=[0,150]×[0,150] 中，其控制方程由两个耦合 PDE 表达: ut=γuΔu+u−u3−v+α 和 vt=γvΔv+β(u−v) ｡这里，u 和 v 代表两个相互作用的组分 (例如，生物组分), γ_u=1 和 γ_v=100 是扩散系数，α=0.01 和 β=0.25 是反应项的系数，Δ 是拉普拉斯算子｡FN 方程通常用于描述由外部刺激 () 激发的生物神经活动，表现出激活 - 抑制系统，因为一个方程促进两个组分的生成，而另一个方程则消散其新增长｡我们采取三种随机场作为初始条件，以生成用于发现的三个独立数据集，每个数据集由 31 个低分辨率快照组成 (投影到 31 × 31 网格), 在 10% 噪声条件下从高保真模拟中降采样 (见补充图 8)｡我们假设 PDE 中的扩散项 ( Δu 和 Δv ) 是已知的，但系数 ( (γ_u 和 γ_v ) 仍需识别｡设计了一个包含 72 个候选函数的库 ( (ϕ∈R^1×72) ) 用于发现耦合 PDE, 尤其是非线性反应项｡与前一个示例类似，图 1b 中的根 - 分支网络用于发现｡有关数据生成､特定候选函数和模型训练的更多描述见补充说明 2.3.2｡图 4a 和 4b 描述了 72 个候选函数的稀疏系数演变 ( λu,λv∈R^72×1 )｡预训练步骤提供了系统在 72 个候选项上的冗余投影；然而，在第一次 ADO 迭代之后，次要候选项被剔除｡其余的 ADO 迭代继续细化所有可训练参数，包括 θ , λu 和 λv ｡最终发现的 PDE 如图 4 所示，与真实值比较｡可以看出，PDE 的形式被精确揭示，所有正确的活动项 (包括第一个方程中的未知外部刺激) 被包含在内｡相应的识别系数一般接近真实值，除了 v 的扩散系数 (即 γ_v ), 根据我们的测试，这似乎是一个不太敏感的参数｡需要注意的是，在此示例中，由于测量数据极其稀疏且嘈杂，根 - 分支 DNN 面临着准确建模具有尖锐传播前沿的解的挑战 (见图 4c)｡DNN 的不太准确的解近似影响了发现的精度｡这个问题可以通过增加时空测量的分辨率自然缓解 (即使在相当大的噪声污染下，例如 10%)｡尽管如此，在这一具有挑战性的示例中，PDE 的确切形式成功被发现，这被认为是更重要的，因为在额外数据到达时，系数可以进一步调优 / 校准｡图 4c 显示了预测的 u 和 v 组分的典型快照､真实值参考和某一未测量时间点 ( (t=18.72) ) 的误差分布｡堆叠的全场 ℓ2 误差为 5.02%｡

图. 4: 基于在 10% 噪声下采样的三种初始条件 (ICs) 发现的 Fitzhugh–Nagumo 方程｡a. 用于构建第一个 PDE (u 方程) 的 72 个候选函数的稀疏系数 λu∈R72 的演变，其中颜色表示系数值｡b. 第二个 PDE (v 方程) 的稀疏系数 λv∈R72 的演变｡为了可视化目的，我们将 u 方程中的常数刺激项 “1” 的识别系数乘以 100, 并将 V 方程中的扩散项 Δv 的识别系数除以 50。c. 在未测量的时间实例 (t=18.72) 下，对于所有三个 IC 的预测响应､真实值和误差分布的快照｡预测的全场响应 (堆叠的 u 和 V v) 的相对 ℓ2 误差为 5.02%｡

实验发现细胞迁移和增殖的 PDE

最后一个示例展示了提出的方法用于发现描述细胞迁移和增殖的控制 PDE, 基于从体外细胞迁移 (划痕) 实验中收集的稀疏和嘈杂的实验数据｡不同时间点 (0 小时､12 小时､24 小时､36 小时､48 小时) 的 1D 细胞密度分布通过高分辨率成像提取，采用图像分割和细胞计数的方法进行处理｡根据不同的初始细胞密度 (例如，总细胞数从 10,000 到 20,000) 进行了多次实验，按照在补充图 10 中显示的测试孔的指定初始分布进行设置｡实验设置和数据集的更详细描述见参考文献 46｡我们在这里的目标是揭示用于建模细胞密度动态的简约 PDE｡这里我们考虑了细胞初始数量范围从 14,000､16,000､18,000 到 20,000 的四种场景｡我们从三次相同准备的实验重复中取平均值，以训练我们在图 1a 中展示的模型用于 PDE 发现｡鉴于我们先前知道细胞动态可以通过扩散 (迁移) 和反应 (增殖) 过程来描述，我们假设 PDE 具有以下形式:

其中 γ 是未知的扩散系数，F 表示潜在的非线性反应函数｡我们使用 8 个附加候选项 (例如，1,ρ,ρ²,ρ³,ρ_x,ρρ_x,ρ²ρ_x,ρ³ρ_x) ) 重建 F , 其系数是稀疏的｡因此，所需的可训练系数总数仍然是 9 (例如，Λ∈R^9×1 )｡我们认为，结合我们的领域特定先验知识是合理的，并且在可解释的模型发现中应得到鼓励，这可能有助于在可用数据非常稀疏和嘈杂时提高我们的解的信心 (例如，在此示例中)｡

图 5a 显示了在 18,000 细胞示例中的 9 个系数的演变，其中冗余候选项通过将相应系数硬阈值化为零而在第一次 ADO 迭代后被剔除｡随后的 ADO 迭代经过后续调整，精细化了活动项的系数，以最终重构 PDE｡图 5b 描述了不同细胞数量的识别活动项系数和相应 PDE, 其中均统一为 ρ_t=γρ_xx+λ₁ρ+λ₂ρ² , 与著名的 Fisher-Kolmogorov 模型完全匹配｡迁移 (扩散) 和增殖 (反应) 速率通常随着细胞数量的增加而增加，如图 5b 中识别的系数所示｡通过发现的 PDE, 我们模拟 / 预测不同时间点 (12 小时､24 小时､36 小时和 48 小时) 的细胞密度演变，如图 5c–f 所示，其中在 0 小时的测量用作初始条件，而 ρx(x=0,t)=ρx(x=1900,t)=0 被用作 Neumann 边界条件｡预测和测量之间的满意一致性为我们发现的 PDE 提供了清晰的验证｡值得注意的是，极其稀疏和嘈杂的实验数据集不幸地给任何现有方法 (如 SINDy) 带来了难以处理的挑战，无法产生合理的发现｡这个实验示例进一步展示了所提出方法在处理高度数据稀缺和噪声方面的能力和强大｡

图. 5: 细胞迁移和增殖的发现结果｡a. 对于 18,000 细胞的例子，9 个候选函数的稀疏系数 Λ∈R9 的演变，冗余候选项通过将相应的系数硬阈值化为零而在第一次 ADO 迭代后被剔除｡后续的 ADO 迭代随后的后调优精细化了活跃项的系数以最终重构 PDE｡b. 不同细胞数量的识别活跃项系数和相应 PDE, 其中统一形式为 ρt=γρxx+λ1ρ+λ2ρ2 , 完全匹配著名的 Fisher-Kolmogorov 模型。c,f. 基于发现的 PDEs 对 14,000､16,000､18,000 和 20,000 细胞的不同时间点的细胞密度进行模拟 / 预测，其中 0 小时的测量作为初始条件，而 ρx(x=0,t)=ρx(x=1900,t)=0 被用作 Neumann 边界条件｡模拟结果以实线表示，标记表示测量数据｡

讨论

总之，我们提出了一种新颖的深度学习方法，用于从稀疏和嘈杂的数据中发现物理规律，特别是简约闭合形式的 PDE (在科学研究和现实应用中普遍存在)｡这种方法结合了 DNN 在非线性函数丰富表示学习方面的优势､自动微分在准确导数计算中的应用以及 ℓ0 稀疏回归，以应对现有稀疏促进方法在数据噪声和稀缺性方面的根本限制｡使用聚合点 (与测量数据无关) 可以使所提出的框架对稀疏和嘈杂的测量具有良好的容忍性，使 DNN 的 PDE 解近似具有良好的泛化能力｡特殊的网络架构设计能够处理在不同初始 / 边界条件下采样的多个独立数据集｡提出了一种交替方向优化策略，以同时训练 DNN 并确定所选候选项的最佳稀疏系数，从而重构 PDE｡DNN 和稀疏 PDE 表示的协同作用导致以下结果：DNN 提供了解的准确建模及其导数作为构建控制方程的基础，而稀疏表示的 PDE 反过来告知和约束 DNN, 这使得其具有良好的泛化能力并进一步增强了发现能力｡总体方法根植于自下而上 (数据驱动) 和自上而下 (物理启发) 过程的综合整合，融合了物理启发的深度学习､稀疏回归和优化｡我们在多个表现出非线性时空行为 (例如混沌､冲击､传播前沿等) 的动态系统上演示了该方法，这些系统由多维 PDE 控制，基于单个或多个数据集，数值或实验数据｡结果突出显示，该方法能够在信息贫乏的空间中准确发现控制方程的确切形式，即使在多维测量稀缺且嘈杂的情况下｡所提出的方法还保持对不同类型噪声 (高斯和非高斯) 的满意鲁棒性，用于 PDE 的发现｡

然而，目前的 PINN-SR 框架在物理规律发现中仍然存在一些潜在限制｡首先，我们必须承认与最先进的 SINDy 方法相比，PINN-SR 的计算成本高得多，主要是由于 DNN 训练的耗时｡然而，SINDy 的关键瓶颈在于其对大规模高质量 (干净) 结构化测量数据的需求，因其使用数值微分，这在数据稀疏且嘈杂的实际应用中 (例如，细胞迁移和增殖示例中的实验数据) 造成了关键限制｡显然，计算效率与对高质量数据需求之间存在权衡｡另一个限制是，尽管本研究中使用的全连接 DNN 在通过自动微分对 PDE 导数进行解析逼近方面具有优势，但直接将其应用于建模高维系统的解 (例如 3D 域中的长期 / 短期响应演变) 会导致计算瓶颈和优化挑战，例如，为了保持满意的准确性，需要大量的聚合点｡具有时空离散化的离散 DNN (例如卷积长短期记忆网络 ConvLSTM 或类似网络) 有潜力帮助解决这个挑战，这将在我们的未来工作中得到展示｡此外，当可用大量独立数据集时，根 - 分支方案可能会面临可扩展性问题，导致 PDE 解近似网络的许多分支｡可训练变量的数量､保持解准确性所需的聚合点以及因此计算内存，通常将随着独立数据集数量的增加而线性增长 (例如，O(r)) )｡尽管如此，这一问题可以通过多 GPU 并行化得到有效解决｡理想情况下，如果 I/BC 是已知的，并且可以在大量多样化数据集可用的条件下进行参数化，则可以在所提出的 PINN-SR 中开发参数化 DNN 学习方案，或神经算子学习以考虑不同的 I/BC｡然而，我们强调，假设有大量数据集超出了我们当前研究的范围，因为在方程发现相关应用中通常难以满足这一要求｡

当前版本的 PINN-SR 不适用于 PDE 系数变化的情况 (例如，时间和 / 或空间依赖)｡然而，鉴于 PINN 在识别 PDE 的变化系数方面的能力，PINN-SR 可以自然扩展以发现闭合形式的 PDE, 其中变化系数分别建模和识别｡此外，PINN 在建模具有混沌行为或尖锐传播波前的系统方面不太理想，主要由于其解场近似的全局基础方式｡当缺少标记数据时 (例如，给定 I/BC 解决 PDE) 或模型形式未知时 (例如，受隐藏物理约束的数据驱动建模), 这种限制尤其明显｡然而，当标记数据相对丰富且明确给出 PDE 模型时 (例如，基于库的模型), 这种限制可以显著缓解｡然而，如果训练数据稀疏和嘈杂，学习的全场响应仍然存在传播波前的误差｡尽管这些误差对发现的 PDE 结构影响不大，但导致 PDE 系数的识别不够准确｡具有局部基础支持的网络可能有助于解决此问题｡最后，虽然 PINN-SR 依赖于预定义的候选项库，但设计一个先验包含但不必过大的库仍然是一项困难的任务 (见方法部分中的更多细节)｡结合表达树或符号神经网络与 PINN 和自动微分有潜力突破库基方法在稀疏和嘈杂数据条件下的限制｡

方法

本工作的创新建立在深度神经网络的优势之上，以进行丰富的表示学习､物理嵌入､自动微分和稀疏回归，来 (1) 近似系统变量的解，(2) 计算必要的导数，以及 (3) 识别形成 PDE (偏微分方程) 结构和显式表达式的关键导数项和参数｡主要方法组成部分如下｡

网络架构

所提出的 PINN-SR 网络架构如图 1a､1b 所示，分别处理单一初始 / 边界条件 (I/BC) 数据集和多个 I/BC (r) 独立数据集。潜在解 u 通过一个密集 (完全连接的) DNN 表示，如图 1a 所示，即 u_θ=u(x,t;θ) , 用于单一数据集的情况，而 “根 - 分支” 密集 DNN 如图 1b 所示，旨在近似对应不同 I/BC 的潜在解 ui (i=1,...,r) ), 即 u_θi=u(x,t;θ⁽⁰⁾,θ⁽ⁱ⁾) , 用于多个独立数据集｡在这里，表示 DNN 的可训练参数｡DNN 接受时空域坐标 x,t 作为输入，并经过多个完全连接的前馈隐藏层 (每层有几十个节点)｡我们使用双曲正切 (tanh) 或正弦 (sin) 作为通用激活函数，因为它们在高阶微分和对正负值的无偏估计方面具有优势｡当系统响应表现出周期性模式时，使用正弦函数｡输出则基于线性激活进行通用幅度映射｡

当有多个数据集可用时，例如从不同 I/BC 采样，域坐标输入到 “根” 网络 (共享隐藏层), 然后经过个 “分支” 网络 (各自的隐藏层), 预测与每个 I/BC 或数据集对应的系统响应｡“根” 学习所有数据集的共同模式 (例如，解的均匀部分), 而 “分支” 学习由每个 I/BC 确定的特定细节 (例如，特定 I/BC 导致的因果关系)｡结果的 “根 - 分支” 网络在统一的基础物理限制下，能够考虑不同的 I/BC｡这种架构整合了来自不同测量的数据，代价是更大的计算工作，并产生满足统一物理 (例如，控制 PDE) 的解近似，这本质上增强了 PINN 进行多源数据驱动建模的能力｡DNN 在此基本上起到非线性泛函的作用，用以近似潜在解｡

DNN 与物理法则 (PDE 重建) 通过自动微分器连接，在机器精度下评估的导数｡候选函数库 ϕ0 可以从 DNN 计算得出｡对于多个独立数据集的情况，来自 “分支” 网络的库 ϕ(i) 被串联以构建统一的 ϕ_θ 如果存在未知的源输入，则可以将候选函数 p 也纳入库中进行发现。重构的 PDE 的稀疏表示随后以残差形式表示:

其中 Rθ∈R1×n 表示 PDE 残差，S 表示稀疏性约束集，n 是系统变量的维度 (例如，u∈R^1×n )｡因此，整体网络架构由异质可训练变量组成，即 DNN 参数 θ∈R^n_θ×1 和 PDE 系数 Λ∈S⊆R^s×n , 其中 nθ 表示 DNN 可训练参数的数量，且 n_θ≫sn

物理约束稀疏正则化损失函数

物理约束稀疏正则化损失函数，如公式 (3) 所示，由三部分组成：数据损失 Ld ､残差物理损失 Lp 和施加在 上的稀疏性正则化项｡数据损失函数表示为:

（4）

其中 um 为测量数据，uθ 为相应的 DNN 近似解，Nm 为数据点总数，∥⋅∥2 表示 Frobenius 范数｡当多个数据集可用时，例如 um=um1,...,umr 和 uθ=uθ1,...,uθr , r≥2 , 如图 1b 所示｡PDE 残差 Rθ 在大量随机采样的协作点 Dc 上评估，并用于形成残差物理损失函数，给出如下:

（5）

其中 和 分别表示一阶时间导数项的离散化和在协作点上评估的候选函数库；Nc 为时空协作点的总数｡对于多个数据集的情况，U˙ 和 在不同 I/BC 的索引上进行串联，以确保施加相同的基础物理法则 (特别是控制 PDE)｡注意，Ld 确保 DNN 通过拟合数据准确解释 PDE 的潜在解，而 Lp 则通过重构 PDE 的闭合形式为 DNN 提供约束和推广｡公式 (3) 中的 ℓ0 正则化项促进了系数 的稀疏性，以稀疏表示 PDE｡

交替方向优化

通过解决公式 (3) 中定义的优化问题进行网络的暴力训练是高度不可行的，因为 ℓ0 正则化使得这个问题为 NP 困难｡虽然通过较不严格的 ℓ1 正则化放松 ℓ0 项，提高了良好解的可行性，并使优化在连续空间中得以进行，但会出现假阳性识别，导致无法实现 PDE 系数的准确稀疏性｡为了解决这个挑战，我们提出了一种交替方向优化 (ADO) 算法，将整体优化问题分解为一组可处理的子问题，以顺序优化和，如公式 (6) 所示:

（6a、b）

ADO 算法的基本概念类似于 (或可以看作是) 乘子交替方向方法｡每次交替迭代 中，通过 STRidge (一个顺序阈值回归过程，作为 ℓ0 正则化的代理) 更新公式 (6a) 中的稀疏 PDE 系数，基于上一次迭代的 DNN 参数 (例如，θk∗ )｡STRidge 的收敛分析可以在文献中找到｡当前迭代中的 DNN 参数随后通过标准的 DNN 训练算法更新 (特别是结合 Adam + L-BFGS 优化器), 将 Λ_k+1^∗ 视为已知｡注意，在训练网络时应使用足够的迭代次数，以实现 θ_k+1^∗ 的满意解精度｡子优化解之间的交替将导致满足全局收敛的高质量优化解。ADO 序列以 q 线性收敛 (见定理 1), 其中 q 表示 "商"。关于广义交替优化的详细理论分析可在相关文献中找到｡值得注意的是，Adam 优化器在全局搜索中发挥作用，而 L-BFGS 优化器则负责在局部解区域进行微调｡测试示例中的 Adam 学习率范围从 10−5 到 10−3 。

定理 1

设 Θ∗=θ∗,Λ∗ 为总损失函数 L(θ,Λ;Du,Dc) 的局部最小值，且 在邻域 N(Θ∗;δ) 内严格凸，其中 表示可训练参数的数量｡我们选择 0<ϵ≤δ , 使得 L 在 N(Θ∗;ϵ) 上严格凸｡如果 y=f(θ;Λ∗)∈N(Θ∗;ϵ) , 且 θ∗ 在 L(θ;Λ∗;Du;Dc) 中局部最小，则 为唯一的全局最小值｡对 Λ∗ 也同样适用｡对于任意可接受的初始解 Θ0∈N(Θ∗;ϵ) , 相应的 ADO 迭代序列理论上以 q 线性收敛。实际收敛速率取决于每次 ADO 迭代中的误差传播｡

PINN-SR 的预训练在运行 ADO 算法进行发现之前进行，通过简单地用 ∥Λ∥0 替换公式 (3) 中的 ∥Λ∥1 , 使得对和 的暴力梯度基优化成为可能:

其中 γ 表示 ℓ1 正则化参数｡ ℓ1 正则化预训练可以通过提供可接受的 “初始猜测” 来加速 ADO 的收敛｡在预训练过程中，DNN 学习稀疏和噪声数据下的物理模式，弱约束于控制 PDE 的回归形式｡后训练 (或后调优) 也是适用的，可以在发现 PDE 的闭合形式 (结构) 之后进行｡可以通过训练 DNN 同时识别发现的非零系数实现后训练，即:

（8）

其中未知参数 θ,Λ 的初始化可以继承自 ADO 结果｡后训练步骤是完全可选的，因为 ADO 方法已经可以提供高质量的解，如测试示例所示｡尽管如此，后训练通过微调可以进一步提高发现的准确性｡

值得一提的是，多步训练的基本直觉在深度学习社区中得到了广泛应用并被证明有效，特别是在 DNN 压缩 (例如，网络预训练､权重修剪和后训练) 方面｡所提出的训练策略类似于这一常用过程｡提出的三步训练策略的启发性理由是：预训练阶段学习良好的 PDE 解近似，ADO 揭示稀疏的 PDE 结构，而后训练阶段对发现的 PDE 结构的系数进行微调｡

超参数选择

适当选择超参数 (例如，α ,B,) 及补充算法 1 和算法 2 所需的超参数) 确保所提出方法在 PDE 发现中的成功｡在本研究中，超参数的选择遵循以下启发性一致标准｡

α: 此超参数平衡网络训练中数据和物理正则化的损失贡献，一般可以根据测量响应 um 和其时间导数 umt (通过有限差分估计 / 近似) 之间的规模比率来估计｡具体而言，的大小设置为与 um 和 umt 之间的偏差比率相似，即 α~rσ=σ(umt)2/σ(um) ｡注意，为了突出测量数据的 PDE 解近似，我们通常在预训练阶段将 α 的值降低数倍 (例如，), 以放宽物理约束｡在 ADO 和后训练中，\alpha 的值增加 (例如，\alpha \sim r\sigma), 以增强 PDE 结构的发现和 PDE 系数的微调｡然而，我们也发现一些例外，例如 λ−ω方程，其中α在预训练和 ADO 阶段都应设置大于缩放比｡由于螺旋模式中 u 和 v 之间的高度相似性，如果测量中包含来自不同 I/BC 的数据集，这种情况可以得到缓解｡尽管如此，我们必须提到，如何选择此超参数是 PINN 社区中的一个普遍且关键的开放性问题｡

β: 此超参数是正则化项的系数，用于促进系数 的稀疏性，以便于稀疏表示 PDE｡它的选择可以根据具体问题的需要进行调整，以确保在学习过程中保持一定的模型复杂性，同时避免过拟合｡

γ : 这是 ℓ1 正则化的系数，用于在预训练阶段增强解的稀疏性｡通过合适的选择，可以帮助模型在面对噪声数据时更好地保持稳定性｡

在测试中，所有超参数的选择均经过多次试验和调优，以找到最优的配置｡我们采用了不同的超参数组合，以测试其对模型性能的影响，并记录了每次实验的结果，以确保在实际应用中能获得最佳效果｡

选择候选函数

候选函数库是 PINN-SR 中的重要组成部分｡设计一个包含先验知识但不必要庞大的候选库是一项困难的任务｡一方面，我们希望使候选库尽可能多样化；另一方面，平衡日益增加的理论和计算复杂性对应用至关重要｡我们相信，基于领域特定知识和统计经验的专门库能够限制搜索空间并减少 PDE 发现的复杂性｡

尽管库的维度越高，从数据中发现精确项的可能性就越大，但一个规模庞大的库 (例如，组件数量在 103 量级) 在本质上会由 DNN 近似，因此很可能会出现秩缺陷和条件差的问题，此外还伴随着日益增加的理论复杂性和计算负担｡平衡这些问题并基于领域特定知识找到建立有效候选库的数学原理仍然是一个开放问题｡此外，未能包括重要的候选函数将导致对简约闭合 PDE 的假阳性发现，尽管可以找到一个 “最佳拟合” 形式 (见补充说明 3.2)｡由于大多数知名的一阶 PDE 可以通过几个活动的线性 / 非线性项的线性组合表示，因此我们在本研究中尝试尽可能多地包含常见的多项式基函数项｡

---

文章转载自微信公众号：智能流体计算

原文链接：https://mp.weixin.qq.com/s/wEPN8OaSOWI7tGIDD9apbQ

从稀疏数据中学习控制方程的物理启发

Physics-informed learning of governing equations from scarce data

摘要：通过数据发现描述复杂物理系统的控制方程或规律，可以极大地促进各个科学和工程领域对这些系统的建模､仿真和理解｡本文提出了一种新颖的方法，称为稀疏回归的物理启发神经网络，用于从稀疏和噪声数据中发现非线性时空系统的控制偏微分方程｡特别地，该发现方法无缝集成了深度神经网络在丰富表示学习､物理嵌入､自动微分和稀疏回归方面的优势，以逼近系统变量的解，计算必要的导数，以及识别形成方程结构和显式表达式的关键导数项和参数｡该方法的有效性和鲁棒性在多种偏微分方程系统的数值和实验上得到了验证，这些系统考虑了不同的数据稀缺性和噪声水平，以及不同的初始 / 边界条件｡所得的计算框架显示了在大规模准确数据难以获取的实际应用中发现封闭形式模型的潜力｡

复杂动态系统建模的当前实践

目前对复杂动态系统建模的实践主要基于描述系统行为的常微分方程 (ODE) 和偏微分方程 (PDE)｡这些控制方程通常是通过严格的第一性原理 (如守恒定律) 或基于知识的现象学推导而获得的｡然而，仍然存在许多未被充分探索的现实复杂系统，其解析描述尚未发现，并且控制方程的简约闭合形式仍不清楚或部分未知｡幸运的是，观察数据集变得越来越丰富，提供了一种从数据中提取潜在方程的替代方案｡利用数据揭示控制规律或方程，可以显著推进和转变我们对各种科学和工程领域复杂物理系统的建模､仿真和理解｡例如，从观察数据 (如卫星遥感图像) 获取控制海冰演变的数学方程，为更好地理解和预测北极冰盖的生长､融化和移动提供了明显的好处｡从现场传感数据 (如多普勒雷达记录) 提炼出明确的公式将加速对天气和气候模式的更准确预测｡最近，机器学习理论､计算能力和数据可用性方面的进展激发了对基于数据的物理法则和控制方程发现的显著热情和努力｡

Bongard 和 Lipson 的开创性贡献利用分层符号回归和遗传编程成功地从数据中提取了控制非线性系统动态的潜在微分方程｡然而，这种优雅的方法在系统的维度上无法很好地扩展，计算成本高，并且可能遭遇过拟合问题｡最近，Brunton 等人的一项令人印象深刻的突破，提出了一种称为非线性动力学稀疏识别 (SINDy) 的创新性稀疏促进方法，通过稀疏回归从高维非线性函数空间中选择主导候选函数，以揭示简约的控制方程，特别是 ODE｡稀疏性是通过顺序阈值岭回归 (STRidge) 算法实现的，该算法递归确定受硬阈值约束的稀疏解｡这样的做法能够在所识别模型的复杂性和准确性之间取得平衡，从而实现简约｡近年来，SINDy 引起了极大的关注，导致了一系列变体算法的出现，用于识别以一阶 ODE 形式的投影低维代理模型，或通过线性嵌入用于广泛范围的非线性动态系统，如流体流动､结构系统､生物和化学系统､活性物质､非线性动态的预测控制､多时间尺度系统､捕食 - 猎物系统和随机过程等｡还有许多其他 SINDy 的扩展，通过发现隐式动态､纳入物理约束以及嵌入随机抽样来提高对稀疏发现高维动态的鲁棒性｡

SINDy 框架在偏微分方程的数据驱动发现方面的一个关键瓶颈在于对测量数据的质量和数量的强依赖，因为需要计算导数以构建控制方程｡特别地，使用有限差分或过滤来计算导数导致了一个关键挑战，降低了算法的鲁棒性｡这特别限制了 SINDy 在当前形式下应用于数据高度不完整､稀疏和嘈杂的场景｡值得注意的是，变分系统识别表明，基于同质分析计算导数的鲁棒性在发现 PDE 的弱形式时表现令人满意｡然而，这种方法在可用数据的保真度方面不易扩展｡另一个研究表明，弱形式能够显著提高对噪声的发现鲁棒性，但需要精心设计测试函数，这在高维时空系统中难以实现｡

自动微分被证明是解决上述问题的有效工具，这在物理启发神经网络 (PINN) 中成功应用于非线性 PDE 的正向和反向分析｡特别地，深度神经网络 (DNN) 用于逼近由 PDE 和少量可用数据约束的解｡PINN 在解决各种科学问题 (如流体流动､涡旋诱导振动､心血管系统等) 时引起了越来越多的关注，特别是在 PDE 的显式形式已知的情况下｡最近，Raissi 的一项重要工作介绍了一种深度隐藏物理模型，基于稀疏数据进行时空动态的数据驱动建模，其中未知的潜在物理通过可能的 PDE 项以弱约束的方式施加并由辅助神经网络隐式学习｡然而，所得到的模型仍然是一个 “黑箱”, 缺乏足够的可解释性，因为无法揭示闭合形式的控制方程｡最新研究显示，在预定义的 PDE 项库的约束搜索空间中，使用 DNN 和自动微分从嘈杂数据中获取闭合形式的 PDE 的潜力；然而，由于使用较不严格的稀疏回归和 DNN 训练而发生的假阳性识别也时有发生｡事实上，同时优化 DNN 参数和稀疏 PDE 系数，同时准确地施加稀疏性，是一个复杂的任务，并且在闭合形式 PDE 的发现中仍然是一项重要挑战｡

为此，我们利用这些进展，提出了一种新颖的 PINN-SR 方法 (即，带有稀疏回归的 PINN), 具有解释性和泛化性的显著特征，从稀疏和嘈杂的数据中发现非线性时空系统的控制 PDE｡我们的方法整合了 DNN 在丰富表示学习方面的优势､自动微分在准确导数计算方面的能力，以及 ℓ0 稀疏回归以应对现有方法在数据噪声和稀缺性方面的根本限制｡具体而言，本文涉及两个方法论贡献:

1.一个 “根 - 分支” 网络，由统一的潜在物理约束，能够处理来自不同初始 / 边界条件的少量多数据集｡

2.一种简单但有效的交替方向训练策略，用于异质参数 (即 DNN 可训练参数和稀疏 PDE 系数) 的优化｡

我们的研究在各种 PDE 系统上验证了该方法的有效性和鲁棒性，基于数值和实验数据集｡

结果

带稀疏回归的 PINN 用于 PDE 发现

我们考虑一个多维时空系统，其控制方程可以用以下形式的非线性耦合参数化 PDE 描述:

其中，u=u(x,t)∈R^1×n 是多维潜在解 维度 , 而 u_t 是一阶时间导数项；t∈[0,T] 表示时间，x∈Ω 指定空间；F[⋯] 是 u 及其空间导数的复杂非线性函数，由 λ 参数化；∇ 是相对于 x 的梯度算子；p=p(x,t) 是源项 (注意，在许多常见情况下，p=0 表示系统没有源输入)｡如果已知，PDE 也受初始和边界条件 (I/BCs) 的约束，表示为 I[x∈Ω,t=0;u,u_t]=0 和B[x∈∂Ω;u,▽_xu]=0。

我们的目标是从假设为不完整､稀疏和嘈杂的可用时空测量中找到 F[⋯] 的闭合形式｡我们假设物理定律仅由少数重要项控制，这些项可以从大型候选函数库中选择，其中可以应用稀疏回归｡该假设固有的将方程重构为以下形式 (为简化起见假设零或未知源):

u_t=ϕΛ

这里，ϕ=ϕ(u)∈R^1×s 是由许多候选项构成的广泛符号函数库，例如关于每个空间维度的常数､多项式和三角函数项，按行向量组合为 ϕ=1,u,u²,∇_xu,∇_yu,...,u³∇_xy,sin(u),... , 其中 ⊙ 表示逐元素的哈达玛乘积；表示库中的候选项总数；在 x,y,z 的上下文中，x,y,z 的下标表示导数；Λ∈R^s×n 是稀疏系数矩阵 (ϕ 中仅有活动候选项的值为非零), 例如 Λ=[λ_u,λ_v,λ_w]∈R^s×3 对于 u=u,v,w ｡如果有未知的源输入，源项的候选函数 p 也可以纳入 ϕ 中以进行发现 (见补充说明)｡因此，发现问题可以表述为：给定时空测量数据 Du , 寻找稀疏的Λ 使得方程 (2) 成立｡

我们提出一种新的 PINN-SR 范式，同时对系统响应建模并识别控制 PDE 的简约闭合形式｡该方法的创新算法架构如图 1 所示，其中考虑了来自两个不同 I/BC 场景的采样数据集：(1) 来自单一 I/BC 的一个数据集；(2) 来自多个 I/BC 的 r≥2 个独立数据集｡对于单一数据集的情况，我们通过 DNN (记为 N) 对潜在解 u 进行解释，即 u_θ=u(x,t;θ) , 其中 θ 表示包括权重和偏差的 DNN 可训练参数，如图 1a 所示｡当有多个独立数据集可用时，设计了一种 “根 - 分支” DNN, 如图 1b 所示，用于逼近与不同 I/BC 相关的潜在解 u_i(i=1,...,r) ), 即 u_θi=u(x,t;θ⁽⁰⁾,θ⁽ⁱ⁾) , 其中 θ⁽⁰⁾ 和 θ⁽ⁱ⁾ 分别表示根层 N⁽⁰⁾ 和分支层 N⁽ⁱ⁾ 的可训练参数｡值得注意的是，I/BC 不必事先已知或测量，因为测量数据已经反映了特定的 I/BC (例如，I/BC 和 PDE 解之间存在一一映射)｡DNN 本质上作为非线性函数，用于逼近潜在解及数据损失函数 L_d(θ;D_u)

图. 1: PINN-SR 数据驱动 PDE 发现框架的示意架构｡a. 来自单一 I/BC 的一个数据集的网络，b. 来自多个 I/BC 的 r≥2 个独立数据集的 “根 - 分支” 网络，c. 基于交替方向优化的网络训练示意图｡网络由两个组成部分构成：一个由可训练参数 θ 控制的 DNN, 将时空坐标 x,t 映射到潜在解 u=u,v,w , 以及由一组非线性 PDE 描述的物理定律，这些 PDE 由以未知稀疏系数 Λ 参数化的导数候选函数 ϕ 形成｡注意，对于多个独立数据集的情况，库 ϕ(i) 被串联以构建用于构建统一控制 PDE 的 ϕc ｡总损失函数 L (θ, Λ, Du, Dc) 由数据损失 Ld (θ, Du)､物理损失 αLp (θ, Λ, Dc) 和促进稀疏性的 ℓ0 正则化项 β||Λ||0 组成｡这里，α a 和 β 表示损失函数的相对权重，而 Du 和 Dc 分别代表测量数据和聚合样本｡注意，物理损失以残差形式仅在时空聚合样本上进行评估｡右侧稀疏系数矩阵 (或向量) 中的有色点表示非零值｡未知参数 {θ, Λ} 的同时优化导致了用于推断数据驱动全场解的训练 DNN 和发现的简约闭合 PDE｡

通过自动微分，在机器精度下评估 u 的导数，可以从 DNN 计算候选函数库 ϕ_θ ｡在多个独立数据集的情况下，来自分支网络的库 ϕ⁽ⁱ⁾ 被串联，以构建统一的 来构建统一的控制 PDE｡因此，重构 PDE 的稀疏表示可以用残差形式表示，即 Rθ:=u_θ^t−ϕ_θΛ=0 , 其中 R_θ∈R^1×n 表示 PDE 残差｡基本概念是同时调整 DNN 可训练参数 θ 和 PDE 系数 Λ , 使得神经网络可以拟合测量数据，同时满足由基础 PDE 定义的约束｡PDE 残差将在大量在时空空间中随机抽样的聚合点 Dc={xi,ti}_i=1^Nc 上评估，形成残差物理损失函数 L_p(θ,Λ,D_c) ｡当考虑多个 I/BC 时，计算数据损失和物理损失时，测量数据和聚合点将叠加 (基于统一的物理残差公式 R_θ→0 )。

PINN-SR 网络的总损失函数因此由数据损失 Ld ､残差物理损失 Lp 和正则化项组成，表示为:

其中，α 是残差物理损失函数的相对权重；β 是正则化参数；表示 ℓ0 范数｡优化总损失函数可以产生一个 DNN, 它不仅能够预测基于数据的全场系统响应，还能够揭示简约的闭合 PDE, 即 f _θ˙∗,Λ^∗:=argmin_θ,ΛL(θ,Λ,Du,Dc) , 其中 θ^∗,Λ^∗ 表示最优参数集合｡值得注意的是，总损失函数具有隐式复杂形式，因此直接求解优化问题是高度不可行的，因为 ℓ0 正则化使这个问题变成 NP 难题｡为了解决这个挑战，我们提出了一种交替方向优化 (ADO) 算法，将整体优化问题分解为一组可处理的子问题，以便在几次交替迭代中顺序优化可训练参数，如图 1c 所示｡在运行 ADO 算法进行发现之前，PINN-SR 的预训练通过简单地用 ∥Λ∥1 替换 Eq. (3) 中的 ∥Λ∥0 来进行，其中对 θ 和 的强制梯度优化变得适用｡ℓ1 正则化的预训练可以通过提供一个可接受的 “初始猜测” 来加速 ADO 的收敛｡更详细的公式和算法描述见方法部分和补充说明 1｡

DNN 和稀疏回归的协同作用导致以下结果：DNN 提供潜在解､其导数和候选函数项的准确建模，作为构建控制 PDE 的基础，而稀疏表示的 PDE 反过来约束了 DNN 的建模并投影出正确的候选函数，最终将测量系统转变为闭合形式的 PDE｡

单一数据集的基准 PDE 发现

我们观察到我们的方法在一组经典 PDE 上的有效性和鲁棒性，这些 PDE 用于表示具有非线性､周期性和 / 或混沌行为的广泛物理系统｡特别地，我们从由多个传感器在固定位置记录的稀疏和嘈杂时间序列测量中发现了 Burgers 方程､Kuramoto–Sivashinsky (KS) 方程､非线性薛定谔方程､Navier–Stokes (NS) 方程和 λ−(1) 反应 - 扩散 (RD) 方程的闭合形式 (数据受高斯白噪声污染), 这些数据来自单一 I/BC｡结果在表 1､图 2 和补充说明 2.1 中展示，表明该方法在数据稀缺性和噪声方面表现出良好的发现精度和鲁棒性｡

表 1: PINN-SR 发现结果的总结，包括多种基准模型的准确性｡

图. 2: 针对稀疏采样测量数据 (10% 噪声) 发现的选定基准 PDE｡a. 发现的 Burgers’ 方程：用于构建 PDE 的 16 个候选函数 ϕ∈R16 的稀疏系数 Λ 的演变，其中颜色代表系数值｡b. 发现的 KS 方程：36 个候选函数 ϕ∈R36 的稀疏系数 Λ 的演变｡c. 发现的非线性薛定谔方程：候选函数 ϕ∈R40 的稀疏系数 Λ 的演变｡d. 发现的 NS 方程：60 个候选函数 ϕ∈R60 的稀疏系数 Λ 的演变｡e. 发现的 RD 方程：用于重建 u 方程和 v 方程的 110 个候选函数 ϕ∈R110 的稀疏系数 λu∈R110 和 λv∈R110 的演变｡

Burgers 方程：我们首先考虑一个由 1D 粘性 Burgers 方程表示的耗散系统，其动态由

ut=−uux+νuxx 控制，其中 ν (等于 0.1) 表示扩散系数｡该方程描述了在有限时间后系统的衰减稳态粘性冲击，通常出现在简化流体力学､非线性声学和气体动力学中｡我们在 Burgers 方程的解中测试 PINN-SR 方法，该解在高斯初始条件下记录的行进冲击波中得到｡在特定情况下，10 个传感器随机放置在 256 个空间网格中的固定位置，并记录波形 101 个时间步长，导致用于参考文献 6 的数据集的 3.19%｡数据集的完整描述､候选函数库的设计 (16 项) 和模型训练见补充说明 2.1.1｡图 2a 显示了在 10% 噪声数据集上的 Burgers 方程发现｡系数 Λ∈R16×1 的演变展示了对真实值的鲁棒收敛 (误差约 0.88%), 导致精确发现｡训练后的 PINN-SR 能够正确重现来自嘈杂测量的动态响应 (例如，完整场的 ℓ2 预测误差为 1.32%), 如补充图 1 所示｡ADO 算法在第一次交替迭代后收敛，并显示出恢复 PDE 的正确稀疏模式的能力｡我们还在具有 10% 噪声的稀疏 u 测量中发现了 Burgers 方程，具有未知 / 未测量源，在发现潜在控制方程时，源项应该同时考虑和重建｡在这种情况下，我们将 14 个源候选函数 (由 和它们的组合组成) 纳入上述库，总计 30 个候选项，以同时发现 PDE 和重建未知源｡相应的发现结果总结在补充图 12 中，包括发现的方程和源函数､稀疏系数 Λ∈R30×1 的演变和预测的全场响应｡结果表明 PDE 和源项及其系数均被良好识别｡

Kuramoto–Sivashinsky (KS) 方程：考虑另一个具有固有不稳定性的耗散系统，由 1D KuramotoSivashinsky (KS) 方程控制 ut=−uux−uxx−uxxxx , 其中反扩散项 −uxx 导致破坏性行为，而四阶导数 uxxxx 产生混沌模式，使其成为理想的方程发现测试问题｡KS 方程广泛用于建模层流火焰前沿和耗散捕获离子模式等不稳定性｡我们随机选择 320 个点作为固定传感器，并记录波响应 101 个时间步长，导致参考文献 6 使用数据集的 12.6%｡共使用 36 个候选函数来构建潜在 PDE｡该示例的详细描述见补充说明 2.1.2｡值得注意的是，混沌行为在从测量不良的数据中逼近全场时空导数方面提出了显著挑战，尤其是高阶导数 uxxxx ｡现有方法 (例如 SINDy 方法系列) 在这种情况下由于测量非常粗糙和嘈杂而最终失败｡然而，PINN-SR 成功从具有 10% 噪声的子采样稀疏数据中提炼出 KS 方程的闭合形式，如图 2b 所示｡系数 Λ∈R³⁶×1 的演变表明，候选项及其相应系数均被正确识别 (接近原始参数；误差约为 0.94%) ADO 迭代内完成｡训练后的 PINN-SR 预测的全场波形也与精确解相吻合，相对 ℓ2 误差为 2.14%(补充图 2)｡

非线性薛定谔方程：在第三个示例中，我们发现非线性薛定谔方程 iut=−0.5u_xx−∣u∣²u , 其中 u 是一个复场变量｡这个著名的方程广泛用于建模光在非线性光纤中的传播､玻色 - 爱因斯坦凝聚､热等离子体中的朗缪尔波等｡我们从数据集中随机选择 37.5% 的子样本 (例如，从空间网格中随机选择), 如表 1 所示，用 40 个候选函数构建 PDE｡由于该函数是复值的，我们分别对 DNN 输出中的解的实部 (uR) 和虚部 (uI) 进行建模，将它们组合以获得复解 u=uR+iuI , 并构建用于 PDE 发现的复值候选函数｡为避免优化中的复杂梯度问题，我们在残差物理损失 Lp 中使用模值 而不是在公式 (5) 中显示的 ℓ2 范数 (有关更多细节，见补充说明 2.1.3)｡图 2c 显示了在 10% 噪声情况下发现的薛定谔方程｡稀疏系数 Λ∈R40×1 的演变历史清楚地表明 (图 2c; 误差约为 0.08%), 导致 PDE 的闭合形式被准确识别，同时重建的全场响应在实部和虚部上都很好地与精确解匹配，相对 ℓ2 误差为 0.26%(补充图 3)｡

Navier-Stokes (NS) 方程：我们考虑一个通过圆柱体的 2D 流体流动，其局部旋转动态由著名的 Navier-Stokes 旋涡方程控制 wt=−(u⋅∇)w+ν∇2w , 其中 w 是时空变化的涡度，u=u,v 表示流体速度，ν 是运动粘度 (在雷诺数 100 时 ν=0.01 )｡我们利用开放模拟数据和子采样数据集在指示区域内随机选择 500 个空间位置的流动响应 u,v,w , 记录 60 个时间步长｡所得数据集仅占参考文献 6 中数据集的 10%｡该示例的全面讨论见补充说明 2.1.4｡图 2d 总结了在 10% 噪声数据集上发现的 NS 方程的结果｡令人鼓舞的是，即使在 10% 噪声污染下，发现的 PDE 表达式也几乎与真实值相同，导数项及其系数均正确｡对应于 60 个候选函数 ϕ∈R1×60 的系数 Λ∈R60×1 在第一轮 ADO 迭代后迅速收敛到正确值 (图 2d)｡涡度模式和幅度也得到很好的预测，如补充图 5 所示 (所有快照的全场 ℓ2 误差约为 2.58%)｡这个例子为所提出的 PINN-SR 方法提供了一个令人信服的测试案例，能够在稀疏和嘈杂的数据中发现闭合形式的 NS 方程｡

反应 - 扩散 (RD) 方程：上述示例大多是低维模型，复杂性有限｡我们这里考虑一个 λ−(1) 反应 - 扩散 (RD) 系统，在 2D 域中由两个耦合 PDE 控制: ut=0.1∇²u+λ(g)u−ω(g)v 和 vt=0.1∇²v+ω(g)u+λ(g)v , 其中 u 和 是两个场变量，g=u²+v² , ω=−g² , λ=1−g² RD 方程表现出广泛的行为，包括波动现象和化学与生物系统中发现的自组织模式｡考虑的特定 RD 方程显示了在周期性边界条件下的螺旋波｡关于数据集､候选函数选择和 PINN-SR 模型的超参数设置的详细信息见补充说明 2.1.5｡图 2e 显示了在 10% 噪声情况下 110 个候选函数 ϕ∈R^1×110 的稀疏系数 λ_u,λ_v∈R^110×1 的演变｡稀疏项和相关系数都被准确识别，以形成闭合方程 (如图 2e 所示)｡由于 PDE 的复杂性和高维性，在 ADO 中需要稍微更多的训练周期以保持可靠的收敛｡训练后的 PINN-SR 在补充图 6 中的预测响应快照 (例如，在 t=2.95 时) 与真实值接近｡这个例子特别展示了我们方法在高维系统中从高度嘈杂数据中发现控制 PDE 的能力和鲁棒性｡

多个独立数据集的 PDE 发现

为了展示图 1b 中提出的 “根 - 分支” 网络在基于多个独立数据集发现 PDE 时的有效性，我们考虑 (1) 具有冲击行为的 1D Burgers 方程，以及 (2) 描述外部刺激激发的激活 - 抑制神经活动的 2D Fitzhugh– Nagumo (FN) 型反应 - 扩散系统｡测量数据在三种不同 I/BC 下以 10% 的噪声稀疏采样 (例如，时间序列或快照)｡注意，由于测量已经反映了特定的 I/BC, 因此 I/BC 不必事先已知或测量，这与系统响应之间存在唯一的一对一映射｡发现结果如下所述｡

带冲击行为的 Burgers 方程：在此示例中，我们测试了之前讨论的 Burgers 方程，其扩散 / 粘度参数较小 ( (ν=0.01/π≈0.0032) ), 基于通过施加三种不同 I/BC 生成的数据集｡如此小的系数会在具有陡峭梯度的紧凑区域内产生冲击形成 (见图 3c), 这可能会挑战 DNN 的逼近能力，从而影响发现｡三组初始和 Dirichlet 边界条件包括:

其中 G 表示高斯函数｡尽管不同 I/BC 的测量数据集表现出完全不同的系统响应，但它们遵循相同的潜在 PDE, 即 ut=−uux+0.0032u_xx ｡对于所有 I/BC, 我们假设在空间 ( (x∈[−1,1]) ) 中随机部署 30 个传感器测量波传播 (例如 u) 500 个时间瞬时 ( (t∈[0,1]) )｡在此，需要与前一个 Burgers 示例相比更密集的传感器网格，以捕获冲击行为｡图 3a 显示了在 10% 噪声下六个典型传感器记录的一些测量｡采用图 1b 中的三分支网络 ( (r=3) ) 进行发现｡数据集､候选函数库 (16 项) 和模型训练的完整描述见补充说明 2.3.1｡图 3b 描述了候选函数的系数演变 ( (Λ∈R^16×1) ), 其中库中正确的项 ( uu_x 和u_xx ) 被成功提取，而其他冗余项被 ADO 逐渐消除 (例如，几乎阈值到零)｡活动项的系数也被准确识别 (特别是导致冲击形成的小粘度参数，例如 0.0039)｡发现的 PDE 为 ut=−1.002uu_x+0.0032u_xx 。图 3c 和 3d 显示了三种 I/BC 情况的预测响应和误差，堆叠的全场 ℓ2 误差为 0.65%｡

图. 3: 基于在 10% 噪声下采样的三个 I/BC 数据集发现的 Burgers’ 方程 (小粘度)｡a. 三个数据集的嘈杂测量可视化｡注意，虽然有 30 个传感器，但仅在此图中显示了一部分｡b. 用于构建 PDE 的 16 个候选函数 ϕ∈R16 的稀疏系数 Λ∈R¹⁶ 的演变，其中颜色表示系数值｡正确的项 (uu_x 和 u_xx) 及其系数成功识别，而其他冗余项通过 ADO 被消除｡c, d. 三个 I/BC 情况下的预测响应和误差｡由于可视化几乎与预测无法区分，因此未列出真实值 (见补充图 7)｡堆叠预测的全场 ℓ2 相对误差为 0.65%｡

FitzHugh–Nagumo (FN) 反应 - 扩散系统：我们考虑一个 FN 型反应 - 扩散系统，在具有周期性边界条件的 2D 域 Ω=[0,150]×[0,150] 中，其控制方程由两个耦合 PDE 表达: ut=γuΔu+u−u3−v+α 和 vt=γvΔv+β(u−v) ｡这里，u 和 v 代表两个相互作用的组分 (例如，生物组分), γ_u=1 和 γ_v=100 是扩散系数，α=0.01 和 β=0.25 是反应项的系数，Δ 是拉普拉斯算子｡FN 方程通常用于描述由外部刺激 () 激发的生物神经活动，表现出激活 - 抑制系统，因为一个方程促进两个组分的生成，而另一个方程则消散其新增长｡我们采取三种随机场作为初始条件，以生成用于发现的三个独立数据集，每个数据集由 31 个低分辨率快照组成 (投影到 31 × 31 网格), 在 10% 噪声条件下从高保真模拟中降采样 (见补充图 8)｡我们假设 PDE 中的扩散项 ( Δu 和 Δv ) 是已知的，但系数 ( (γ_u 和 γ_v ) 仍需识别｡设计了一个包含 72 个候选函数的库 ( (ϕ∈R^1×72) ) 用于发现耦合 PDE, 尤其是非线性反应项｡与前一个示例类似，图 1b 中的根 - 分支网络用于发现｡有关数据生成､特定候选函数和模型训练的更多描述见补充说明 2.3.2｡图 4a 和 4b 描述了 72 个候选函数的稀疏系数演变 ( λu,λv∈R^72×1 )｡预训练步骤提供了系统在 72 个候选项上的冗余投影；然而，在第一次 ADO 迭代之后，次要候选项被剔除｡其余的 ADO 迭代继续细化所有可训练参数，包括 θ , λu 和 λv ｡最终发现的 PDE 如图 4 所示，与真实值比较｡可以看出，PDE 的形式被精确揭示，所有正确的活动项 (包括第一个方程中的未知外部刺激) 被包含在内｡相应的识别系数一般接近真实值，除了 v 的扩散系数 (即 γ_v ), 根据我们的测试，这似乎是一个不太敏感的参数｡需要注意的是，在此示例中，由于测量数据极其稀疏且嘈杂，根 - 分支 DNN 面临着准确建模具有尖锐传播前沿的解的挑战 (见图 4c)｡DNN 的不太准确的解近似影响了发现的精度｡这个问题可以通过增加时空测量的分辨率自然缓解 (即使在相当大的噪声污染下，例如 10%)｡尽管如此，在这一具有挑战性的示例中，PDE 的确切形式成功被发现，这被认为是更重要的，因为在额外数据到达时，系数可以进一步调优 / 校准｡图 4c 显示了预测的 u 和 v 组分的典型快照､真实值参考和某一未测量时间点 ( (t=18.72) ) 的误差分布｡堆叠的全场 ℓ2 误差为 5.02%｡

图. 4: 基于在 10% 噪声下采样的三种初始条件 (ICs) 发现的 Fitzhugh–Nagumo 方程｡a. 用于构建第一个 PDE (u 方程) 的 72 个候选函数的稀疏系数 λu∈R72 的演变，其中颜色表示系数值｡b. 第二个 PDE (v 方程) 的稀疏系数 λv∈R72 的演变｡为了可视化目的，我们将 u 方程中的常数刺激项 “1” 的识别系数乘以 100, 并将 V 方程中的扩散项 Δv 的识别系数除以 50。c. 在未测量的时间实例 (t=18.72) 下，对于所有三个 IC 的预测响应､真实值和误差分布的快照｡预测的全场响应 (堆叠的 u 和 V v) 的相对 ℓ2 误差为 5.02%｡

实验发现细胞迁移和增殖的 PDE

最后一个示例展示了提出的方法用于发现描述细胞迁移和增殖的控制 PDE, 基于从体外细胞迁移 (划痕) 实验中收集的稀疏和嘈杂的实验数据｡不同时间点 (0 小时､12 小时､24 小时､36 小时､48 小时) 的 1D 细胞密度分布通过高分辨率成像提取，采用图像分割和细胞计数的方法进行处理｡根据不同的初始细胞密度 (例如，总细胞数从 10,000 到 20,000) 进行了多次实验，按照在补充图 10 中显示的测试孔的指定初始分布进行设置｡实验设置和数据集的更详细描述见参考文献 46｡我们在这里的目标是揭示用于建模细胞密度动态的简约 PDE｡这里我们考虑了细胞初始数量范围从 14,000､16,000､18,000 到 20,000 的四种场景｡我们从三次相同准备的实验重复中取平均值，以训练我们在图 1a 中展示的模型用于 PDE 发现｡鉴于我们先前知道细胞动态可以通过扩散 (迁移) 和反应 (增殖) 过程来描述，我们假设 PDE 具有以下形式:

其中 γ 是未知的扩散系数，F 表示潜在的非线性反应函数｡我们使用 8 个附加候选项 (例如，1,ρ,ρ²,ρ³,ρ_x,ρρ_x,ρ²ρ_x,ρ³ρ_x) ) 重建 F , 其系数是稀疏的｡因此，所需的可训练系数总数仍然是 9 (例如，Λ∈R^9×1 )｡我们认为，结合我们的领域特定先验知识是合理的，并且在可解释的模型发现中应得到鼓励，这可能有助于在可用数据非常稀疏和嘈杂时提高我们的解的信心 (例如，在此示例中)｡

图 5a 显示了在 18,000 细胞示例中的 9 个系数的演变，其中冗余候选项通过将相应系数硬阈值化为零而在第一次 ADO 迭代后被剔除｡随后的 ADO 迭代经过后续调整，精细化了活动项的系数，以最终重构 PDE｡图 5b 描述了不同细胞数量的识别活动项系数和相应 PDE, 其中均统一为 ρ_t=γρ_xx+λ₁ρ+λ₂ρ² , 与著名的 Fisher-Kolmogorov 模型完全匹配｡迁移 (扩散) 和增殖 (反应) 速率通常随着细胞数量的增加而增加，如图 5b 中识别的系数所示｡通过发现的 PDE, 我们模拟 / 预测不同时间点 (12 小时､24 小时､36 小时和 48 小时) 的细胞密度演变，如图 5c–f 所示，其中在 0 小时的测量用作初始条件，而 ρx(x=0,t)=ρx(x=1900,t)=0 被用作 Neumann 边界条件｡预测和测量之间的满意一致性为我们发现的 PDE 提供了清晰的验证｡值得注意的是，极其稀疏和嘈杂的实验数据集不幸地给任何现有方法 (如 SINDy) 带来了难以处理的挑战，无法产生合理的发现｡这个实验示例进一步展示了所提出方法在处理高度数据稀缺和噪声方面的能力和强大｡

图. 5: 细胞迁移和增殖的发现结果｡a. 对于 18,000 细胞的例子，9 个候选函数的稀疏系数 Λ∈R9 的演变，冗余候选项通过将相应的系数硬阈值化为零而在第一次 ADO 迭代后被剔除｡后续的 ADO 迭代随后的后调优精细化了活跃项的系数以最终重构 PDE｡b. 不同细胞数量的识别活跃项系数和相应 PDE, 其中统一形式为 ρt=γρxx+λ1ρ+λ2ρ2 , 完全匹配著名的 Fisher-Kolmogorov 模型。c,f. 基于发现的 PDEs 对 14,000､16,000､18,000 和 20,000 细胞的不同时间点的细胞密度进行模拟 / 预测，其中 0 小时的测量作为初始条件，而 ρx(x=0,t)=ρx(x=1900,t)=0 被用作 Neumann 边界条件｡模拟结果以实线表示，标记表示测量数据｡

讨论

总之，我们提出了一种新颖的深度学习方法，用于从稀疏和嘈杂的数据中发现物理规律，特别是简约闭合形式的 PDE (在科学研究和现实应用中普遍存在)｡这种方法结合了 DNN 在非线性函数丰富表示学习方面的优势､自动微分在准确导数计算中的应用以及 ℓ0 稀疏回归，以应对现有稀疏促进方法在数据噪声和稀缺性方面的根本限制｡使用聚合点 (与测量数据无关) 可以使所提出的框架对稀疏和嘈杂的测量具有良好的容忍性，使 DNN 的 PDE 解近似具有良好的泛化能力｡特殊的网络架构设计能够处理在不同初始 / 边界条件下采样的多个独立数据集｡提出了一种交替方向优化策略，以同时训练 DNN 并确定所选候选项的最佳稀疏系数，从而重构 PDE｡DNN 和稀疏 PDE 表示的协同作用导致以下结果：DNN 提供了解的准确建模及其导数作为构建控制方程的基础，而稀疏表示的 PDE 反过来告知和约束 DNN, 这使得其具有良好的泛化能力并进一步增强了发现能力｡总体方法根植于自下而上 (数据驱动) 和自上而下 (物理启发) 过程的综合整合，融合了物理启发的深度学习､稀疏回归和优化｡我们在多个表现出非线性时空行为 (例如混沌､冲击､传播前沿等) 的动态系统上演示了该方法，这些系统由多维 PDE 控制，基于单个或多个数据集，数值或实验数据｡结果突出显示，该方法能够在信息贫乏的空间中准确发现控制方程的确切形式，即使在多维测量稀缺且嘈杂的情况下｡所提出的方法还保持对不同类型噪声 (高斯和非高斯) 的满意鲁棒性，用于 PDE 的发现｡

然而，目前的 PINN-SR 框架在物理规律发现中仍然存在一些潜在限制｡首先，我们必须承认与最先进的 SINDy 方法相比，PINN-SR 的计算成本高得多，主要是由于 DNN 训练的耗时｡然而，SINDy 的关键瓶颈在于其对大规模高质量 (干净) 结构化测量数据的需求，因其使用数值微分，这在数据稀疏且嘈杂的实际应用中 (例如，细胞迁移和增殖示例中的实验数据) 造成了关键限制｡显然，计算效率与对高质量数据需求之间存在权衡｡另一个限制是，尽管本研究中使用的全连接 DNN 在通过自动微分对 PDE 导数进行解析逼近方面具有优势，但直接将其应用于建模高维系统的解 (例如 3D 域中的长期 / 短期响应演变) 会导致计算瓶颈和优化挑战，例如，为了保持满意的准确性，需要大量的聚合点｡具有时空离散化的离散 DNN (例如卷积长短期记忆网络 ConvLSTM 或类似网络) 有潜力帮助解决这个挑战，这将在我们的未来工作中得到展示｡此外，当可用大量独立数据集时，根 - 分支方案可能会面临可扩展性问题，导致 PDE 解近似网络的许多分支｡可训练变量的数量､保持解准确性所需的聚合点以及因此计算内存，通常将随着独立数据集数量的增加而线性增长 (例如，O(r)) )｡尽管如此，这一问题可以通过多 GPU 并行化得到有效解决｡理想情况下，如果 I/BC 是已知的，并且可以在大量多样化数据集可用的条件下进行参数化，则可以在所提出的 PINN-SR 中开发参数化 DNN 学习方案，或神经算子学习以考虑不同的 I/BC｡然而，我们强调，假设有大量数据集超出了我们当前研究的范围，因为在方程发现相关应用中通常难以满足这一要求｡

当前版本的 PINN-SR 不适用于 PDE 系数变化的情况 (例如，时间和 / 或空间依赖)｡然而，鉴于 PINN 在识别 PDE 的变化系数方面的能力，PINN-SR 可以自然扩展以发现闭合形式的 PDE, 其中变化系数分别建模和识别｡此外，PINN 在建模具有混沌行为或尖锐传播波前的系统方面不太理想，主要由于其解场近似的全局基础方式｡当缺少标记数据时 (例如，给定 I/BC 解决 PDE) 或模型形式未知时 (例如，受隐藏物理约束的数据驱动建模), 这种限制尤其明显｡然而，当标记数据相对丰富且明确给出 PDE 模型时 (例如，基于库的模型), 这种限制可以显著缓解｡然而，如果训练数据稀疏和嘈杂，学习的全场响应仍然存在传播波前的误差｡尽管这些误差对发现的 PDE 结构影响不大，但导致 PDE 系数的识别不够准确｡具有局部基础支持的网络可能有助于解决此问题｡最后，虽然 PINN-SR 依赖于预定义的候选项库，但设计一个先验包含但不必过大的库仍然是一项困难的任务 (见方法部分中的更多细节)｡结合表达树或符号神经网络与 PINN 和自动微分有潜力突破库基方法在稀疏和嘈杂数据条件下的限制｡

方法

本工作的创新建立在深度神经网络的优势之上，以进行丰富的表示学习､物理嵌入､自动微分和稀疏回归，来 (1) 近似系统变量的解，(2) 计算必要的导数，以及 (3) 识别形成 PDE (偏微分方程) 结构和显式表达式的关键导数项和参数｡主要方法组成部分如下｡

网络架构

所提出的 PINN-SR 网络架构如图 1a､1b 所示，分别处理单一初始 / 边界条件 (I/BC) 数据集和多个 I/BC (r) 独立数据集。潜在解 u 通过一个密集 (完全连接的) DNN 表示，如图 1a 所示，即 u_θ=u(x,t;θ) , 用于单一数据集的情况，而 “根 - 分支” 密集 DNN 如图 1b 所示，旨在近似对应不同 I/BC 的潜在解 ui (i=1,...,r) ), 即 u_θi=u(x,t;θ⁽⁰⁾,θ⁽ⁱ⁾) , 用于多个独立数据集｡在这里，表示 DNN 的可训练参数｡DNN 接受时空域坐标 x,t 作为输入，并经过多个完全连接的前馈隐藏层 (每层有几十个节点)｡我们使用双曲正切 (tanh) 或正弦 (sin) 作为通用激活函数，因为它们在高阶微分和对正负值的无偏估计方面具有优势｡当系统响应表现出周期性模式时，使用正弦函数｡输出则基于线性激活进行通用幅度映射｡

当有多个数据集可用时，例如从不同 I/BC 采样，域坐标输入到 “根” 网络 (共享隐藏层), 然后经过个 “分支” 网络 (各自的隐藏层), 预测与每个 I/BC 或数据集对应的系统响应｡“根” 学习所有数据集的共同模式 (例如，解的均匀部分), 而 “分支” 学习由每个 I/BC 确定的特定细节 (例如，特定 I/BC 导致的因果关系)｡结果的 “根 - 分支” 网络在统一的基础物理限制下，能够考虑不同的 I/BC｡这种架构整合了来自不同测量的数据，代价是更大的计算工作，并产生满足统一物理 (例如，控制 PDE) 的解近似，这本质上增强了 PINN 进行多源数据驱动建模的能力｡DNN 在此基本上起到非线性泛函的作用，用以近似潜在解｡

DNN 与物理法则 (PDE 重建) 通过自动微分器连接，在机器精度下评估的导数｡候选函数库 ϕ0 可以从 DNN 计算得出｡对于多个独立数据集的情况，来自 “分支” 网络的库 ϕ(i) 被串联以构建统一的 ϕ_θ 如果存在未知的源输入，则可以将候选函数 p 也纳入库中进行发现。重构的 PDE 的稀疏表示随后以残差形式表示:

其中 Rθ∈R1×n 表示 PDE 残差，S 表示稀疏性约束集，n 是系统变量的维度 (例如，u∈R^1×n )｡因此，整体网络架构由异质可训练变量组成，即 DNN 参数 θ∈R^n_θ×1 和 PDE 系数 Λ∈S⊆R^s×n , 其中 nθ 表示 DNN 可训练参数的数量，且 n_θ≫sn

物理约束稀疏正则化损失函数

物理约束稀疏正则化损失函数，如公式 (3) 所示，由三部分组成：数据损失 Ld ､残差物理损失 Lp 和施加在 上的稀疏性正则化项｡数据损失函数表示为:

（4）

其中 um 为测量数据，uθ 为相应的 DNN 近似解，Nm 为数据点总数，∥⋅∥2 表示 Frobenius 范数｡当多个数据集可用时，例如 um=um1,...,umr 和 uθ=uθ1,...,uθr , r≥2 , 如图 1b 所示｡PDE 残差 Rθ 在大量随机采样的协作点 Dc 上评估，并用于形成残差物理损失函数，给出如下:

（5）

其中 和 分别表示一阶时间导数项的离散化和在协作点上评估的候选函数库；Nc 为时空协作点的总数｡对于多个数据集的情况，U˙ 和 在不同 I/BC 的索引上进行串联，以确保施加相同的基础物理法则 (特别是控制 PDE)｡注意，Ld 确保 DNN 通过拟合数据准确解释 PDE 的潜在解，而 Lp 则通过重构 PDE 的闭合形式为 DNN 提供约束和推广｡公式 (3) 中的 ℓ0 正则化项促进了系数 的稀疏性，以稀疏表示 PDE｡

交替方向优化

通过解决公式 (3) 中定义的优化问题进行网络的暴力训练是高度不可行的，因为 ℓ0 正则化使得这个问题为 NP 困难｡虽然通过较不严格的 ℓ1 正则化放松 ℓ0 项，提高了良好解的可行性，并使优化在连续空间中得以进行，但会出现假阳性识别，导致无法实现 PDE 系数的准确稀疏性｡为了解决这个挑战，我们提出了一种交替方向优化 (ADO) 算法，将整体优化问题分解为一组可处理的子问题，以顺序优化和，如公式 (6) 所示:

（6a、b）

ADO 算法的基本概念类似于 (或可以看作是) 乘子交替方向方法｡每次交替迭代 中，通过 STRidge (一个顺序阈值回归过程，作为 ℓ0 正则化的代理) 更新公式 (6a) 中的稀疏 PDE 系数，基于上一次迭代的 DNN 参数 (例如，θk∗ )｡STRidge 的收敛分析可以在文献中找到｡当前迭代中的 DNN 参数随后通过标准的 DNN 训练算法更新 (特别是结合 Adam + L-BFGS 优化器), 将 Λ_k+1^∗ 视为已知｡注意，在训练网络时应使用足够的迭代次数，以实现 θ_k+1^∗ 的满意解精度｡子优化解之间的交替将导致满足全局收敛的高质量优化解。ADO 序列以 q 线性收敛 (见定理 1), 其中 q 表示 "商"。关于广义交替优化的详细理论分析可在相关文献中找到｡值得注意的是，Adam 优化器在全局搜索中发挥作用，而 L-BFGS 优化器则负责在局部解区域进行微调｡测试示例中的 Adam 学习率范围从 10−5 到 10−3 。

定理 1

设 Θ∗=θ∗,Λ∗ 为总损失函数 L(θ,Λ;Du,Dc) 的局部最小值，且 在邻域 N(Θ∗;δ) 内严格凸，其中 表示可训练参数的数量｡我们选择 0<ϵ≤δ , 使得 L 在 N(Θ∗;ϵ) 上严格凸｡如果 y=f(θ;Λ∗)∈N(Θ∗;ϵ) , 且 θ∗ 在 L(θ;Λ∗;Du;Dc) 中局部最小，则 为唯一的全局最小值｡对 Λ∗ 也同样适用｡对于任意可接受的初始解 Θ0∈N(Θ∗;ϵ) , 相应的 ADO 迭代序列理论上以 q 线性收敛。实际收敛速率取决于每次 ADO 迭代中的误差传播｡

PINN-SR 的预训练在运行 ADO 算法进行发现之前进行，通过简单地用 ∥Λ∥0 替换公式 (3) 中的 ∥Λ∥1 , 使得对和 的暴力梯度基优化成为可能:

其中 γ 表示 ℓ1 正则化参数｡ ℓ1 正则化预训练可以通过提供可接受的 “初始猜测” 来加速 ADO 的收敛｡在预训练过程中，DNN 学习稀疏和噪声数据下的物理模式，弱约束于控制 PDE 的回归形式｡后训练 (或后调优) 也是适用的，可以在发现 PDE 的闭合形式 (结构) 之后进行｡可以通过训练 DNN 同时识别发现的非零系数实现后训练，即:

（8）

其中未知参数 θ,Λ 的初始化可以继承自 ADO 结果｡后训练步骤是完全可选的，因为 ADO 方法已经可以提供高质量的解，如测试示例所示｡尽管如此，后训练通过微调可以进一步提高发现的准确性｡

值得一提的是，多步训练的基本直觉在深度学习社区中得到了广泛应用并被证明有效，特别是在 DNN 压缩 (例如，网络预训练､权重修剪和后训练) 方面｡所提出的训练策略类似于这一常用过程｡提出的三步训练策略的启发性理由是：预训练阶段学习良好的 PDE 解近似，ADO 揭示稀疏的 PDE 结构，而后训练阶段对发现的 PDE 结构的系数进行微调｡

超参数选择

适当选择超参数 (例如，α ,B,) 及补充算法 1 和算法 2 所需的超参数) 确保所提出方法在 PDE 发现中的成功｡在本研究中，超参数的选择遵循以下启发性一致标准｡

α: 此超参数平衡网络训练中数据和物理正则化的损失贡献，一般可以根据测量响应 um 和其时间导数 umt (通过有限差分估计 / 近似) 之间的规模比率来估计｡具体而言，的大小设置为与 um 和 umt 之间的偏差比率相似，即 α~rσ=σ(umt)2/σ(um) ｡注意，为了突出测量数据的 PDE 解近似，我们通常在预训练阶段将 α 的值降低数倍 (例如，), 以放宽物理约束｡在 ADO 和后训练中，\alpha 的值增加 (例如，\alpha \sim r\sigma), 以增强 PDE 结构的发现和 PDE 系数的微调｡然而，我们也发现一些例外，例如 λ−ω方程，其中α在预训练和 ADO 阶段都应设置大于缩放比｡由于螺旋模式中 u 和 v 之间的高度相似性，如果测量中包含来自不同 I/BC 的数据集，这种情况可以得到缓解｡尽管如此，我们必须提到，如何选择此超参数是 PINN 社区中的一个普遍且关键的开放性问题｡

β: 此超参数是正则化项的系数，用于促进系数 的稀疏性，以便于稀疏表示 PDE｡它的选择可以根据具体问题的需要进行调整，以确保在学习过程中保持一定的模型复杂性，同时避免过拟合｡

γ : 这是 ℓ1 正则化的系数，用于在预训练阶段增强解的稀疏性｡通过合适的选择，可以帮助模型在面对噪声数据时更好地保持稳定性｡

在测试中，所有超参数的选择均经过多次试验和调优，以找到最优的配置｡我们采用了不同的超参数组合，以测试其对模型性能的影响，并记录了每次实验的结果，以确保在实际应用中能获得最佳效果｡

选择候选函数

候选函数库是 PINN-SR 中的重要组成部分｡设计一个包含先验知识但不必要庞大的候选库是一项困难的任务｡一方面，我们希望使候选库尽可能多样化；另一方面，平衡日益增加的理论和计算复杂性对应用至关重要｡我们相信，基于领域特定知识和统计经验的专门库能够限制搜索空间并减少 PDE 发现的复杂性｡

尽管库的维度越高，从数据中发现精确项的可能性就越大，但一个规模庞大的库 (例如，组件数量在 103 量级) 在本质上会由 DNN 近似，因此很可能会出现秩缺陷和条件差的问题，此外还伴随着日益增加的理论复杂性和计算负担｡平衡这些问题并基于领域特定知识找到建立有效候选库的数学原理仍然是一个开放问题｡此外，未能包括重要的候选函数将导致对简约闭合 PDE 的假阳性发现，尽管可以找到一个 “最佳拟合” 形式 (见补充说明 3.2)｡由于大多数知名的一阶 PDE 可以通过几个活动的线性 / 非线性项的线性组合表示，因此我们在本研究中尝试尽可能多地包含常见的多项式基函数项｡

---

文章转载自微信公众号：智能流体计算

原文链接：https://mp.weixin.qq.com/s/wEPN8OaSOWI7tGIDD9apbQ