热传导方程两类参数的统计反演方法

Statistical Inversion Methods for Two Types of Parameters in the Heat Conduction Equation

摘要: 本文以热传导方程为背景，分别研究了不同位置上的常数型未知参数与函数型未知参数，并给出统计反演算法，完成数值反演试验。两类方法均能高效完成反演任务，且在面对高噪声观测信息时，具有很强的鲁棒性。

Abstract: Focusing on the heat conduction equation, this paper investigates constant-type unknown parameters and function-type unknown parameters at different positions, and provides statistical inversion algorithms to complete numerical inversion experiments. Both methods can efficiently accomplish the inversion tasks and exhibit strong robustness in the presence of high-noise observational information.

1. 引言

作为一类常见的描述实验现象与生产细节的数学物理方程，热传导方程被广泛研究并加以应用。其参数反演问题也因逆向研究视角的挑战性和高价值性，被学界广泛关注。本文将借用两类统计学中的数值方法，研究如下的两组方程的参数反演问题。

问题P1：对于一个长为L，宽为H的二维薄板来说，其热传导方程具有如下形式：

$\{\begin{array}{ll}{u}_t-a\left(u_{xx}+u_{yy}\right)=0,\hfill & \left(x,y,t\right)\in \Omega \times \left[0,1\right],\hfill \\ u\left(x,y,0\right)=\sin \left(\pi x\right)\sin \left(\pi y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \\ u\left(0,y,t\right)=u\left(L,y,t\right)=0,\hfill & y\in \left[0,H\right],t\in \left[0,1\right],\hfill \\ u\left(x,0,t\right)=u\left(x,H,t\right)=0,\hfill & x\in \left[0,L\right],t\in \left[0,1\right].\hfill \end{array}$ (1)

对于方程(1)，其反问题是在已知终端时刻，全局温度$g=u\left(x,y,T\right)$ 的附加条件下，估计出合适的常扩散系数a。

问题P2：该问题则聚焦于复杂边界条件下的一维导热杆热传递过程，其具有如下形式的表达：

$\{\begin{array}{ll}{u}_t-u_{xx}={\text{e}}^{-t}\sin \left(\pi x/L\right),\hfill & \left(x,t\right)\in \left[0,L\right]\times \left[0,T\right],\hfill \\ u\left(x,0\right)=x\left(1.75-x\right),\hfill & x\in \left[0,L\right],\hfill \\ u\left(0,t\right)=0,\hfill & t\in \left[0,T\right],\hfill \\ u_x\left(L,t\right)+\sigma \left(t\right)\cdot u\left(L,t\right)=0,\hfill & t\in \left[0,T\right].\hfill \end{array}$ (2)

该反问题将利用x₀处不同时刻的温度记录$h=u\left(x_0,t\right)$ ，反演右边界上的热交换系数。

在反问题的发展过程中，形成了一些较为高效的参数反演方法。第一类是变分泛函为代表的优化解法，诸如Landweber、共轭梯度方法等各种梯度型迭代算法。这类方法通过构建一个结合了未知参数和观测数据的代价泛函，并推导代价泛函的梯度表达式，然后选择适当的迭代算法来寻找泛函的极值，寻找到的局部或全局最优解作为参数的反演结果。例如，研究者Yang和Deng等人在文献[1]中，基于Landweber正则化迭代方案建立了一个最优控制框架，成功重建了热传导方程中的热源参数，不仅证实了该反问题解的全局唯一性和稳定性，也为后续研究提供了理论基础。此外，Wei采用广义Tikhonov正则化的方法，处理了分数阶扩散方程中的时间依赖源项，并通过引入广义交叉验证策略来选择最优正则化参数，从而获取稳定的数值解[2]。

另一类重要方法是基于贝叶斯推断的统计抽样技术，该方法的一般步骤为推导未知参数的后验概率分布，再结合高效的统计抽样算法估计参数的值。Zabaras等人便是利用这一方法，有效识别了一维及二维热传导方程Neumann边界条件下的热流密度q，同时证明了贝叶斯方法与最优控制方法在线性反问题求解中的等价性[3]。

近年来，各类基于数值驱动策略的反演方法也逐渐进入学者们的视野。吴自库利用最小二乘支持向量机(LS-SVM)解决热传导方程中的热源反演问题，将源项$f\left(t\right)$ 的推断转变为一个二次优化问题，进而通过求解线性方程组获得最终答案[4]。Raissi则提出了一种名为物理信息神经网络(Physics Informed Neural Networks, PINN)的新框架，它能够在小规模观测数据的条件下，高效地解决偏微分方程的正向和反向问题[5]。PINNs不仅为数据驱动方法与物理模型之间的融合提供了一条新路径，而且在解决具体科学和工程问题上展现出了巨大潜力。文中分别求解了一维薛定谔方程与2维N-S方程组数值解，并重构了后者的流体密度和动力粘度参数。

本文将采取后两种统计反演方法，分别给出未知参数a和$\sigma \left(t\right)$ 的数值反演结果。并验证不同观测误差下的算法稳定性。

2. 统计反演方法及其算法搭建

2.1. 正问题数值解法

在介绍参数反演问题前，还需构建精度高，计算速度快的数值计算方法，以便在后续的反演算法中频繁调用。本文将采用有限差分算法来完成该任务，其核心原理是将网格剖分节点处的差分算子近似为微分算子，以代替方程及其边界上的导数运算。

对于一维热传导方程，常用显格式、隐格式以及C-N差分三种差分格式，下面给出前两种格式的计算公式[6]：

① 显格式：

$\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\Delta t}-a\left(x_i,t_j\right)\cdot \left[\frac{u\left(x_{i+1},t_j\right)-2u\left(x_i,t_j\right)+u\left(x_{i-1},t_j\right)}{\Delta x^2}\right]=f\left(x_i,t_j\right)$

② 隐格式：

$\frac{u\left(x_i,t_{j+1}\right)-u\left(x_i,t_j\right)}{\Delta t}-a\left(x_i,t_j\right)\cdot \left[\frac{u\left(x_{i+1},t_{j+1}\right)-2u\left(x_i,t_{j+1}\right)+u\left(x_{i-1},t_{j+1}\right)}{\Delta x^2}\right]=f\left(x_i,t_{j+1}\right)$

可以看见前两种差分的主要区别便是空间项导数上对于时间节点选择不同，显格式是正向计算的，隐格式则是逆向求解线性方程组。至于C-N格式则是前两种的线性组合，其拥有更高的计算精度，这里将不再展示。

但对于二维方程组而言，上面两种格式不总是适用，会存在剖分要求高、非对角线方程组求解成本大等问题。因此迫切需要引入一种新的差分格式，来提高二维热传导方程的计算速度与精度。

Peaceman和Racford提出这样一种算法[6] [7]，他们在计算高维方程时引入了一个中间层$u_{i,j}^{\ast }$ ，将原本关于$u_{i,j}^{n+1}$ 和$u_{i.j}^n$ 的方程组，拆分成两步进行计算：首先，在x方向上使用隐格式，在y方向上则使用显格式，并求解$u_{i.j}^n$ 和中间层$u_{i,j}^{\ast }$ 构成的方程组，得到$u_{i,j}^{\ast }$ ；第二步则交替隐格式的顺序，在y方向上使用隐格式，在x方向上使用显格式，之后求解$u_{i.j}^{n+1}$ 关于$u_{i,j}^{\ast }$ 的方程组，得到$u_{i.j}^{n+1}$ 。以时间变量t为顺序，重复上述两个步骤，直至遍历完所有时间剖分节点，得到该高维方程组的数值解。下面给出上述两个步骤的计算过程：

第一步：$\frac{u_{i,j}^{\ast }-u_{i,j}^n}{\Delta t/2}=a_{i,j}\cdot \frac{1}{h^2}\left(u_{i+1,j}^{\ast }-2u_{i,j}^{\ast }+u_{i-1,j}^{\ast }+u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n\right)$ ，此处令$\Delta x=\Delta y=h$ ，经化简有，

$\begin{array}{l}\underset{A}{\left[\begin{array}{ccccc}1+a_{2,j}r& -a_{2,j}r/2& & & \\ -a_{3,j}r/2& 1+a_{3,j}r& -a_{3,j}r/2& & \\ & & \ddots & & \\ & & -a_{Nx-2,j}r/2& 1+a_{Nx-2,j}r& -a_{Nx-2,j}r/2\\ & & & -a_{Nx-1,j}r/2& 1+a_{Nx-1,j}r\end{array}\right]}\times \underset{u^{\ast }}{\left[\begin{array}{c}{u}_{2,j}^{\ast }\\ u_{3,j}^{\ast }\\ ⋮\\ u_{Nx-2,j}^{\ast }\\ u_{Nx-1,j}^{\ast }\end{array}\right]}\\ =\underset{b_n}{\left[\begin{array}{c}{u}_{2,j}^n+\frac{a_{2,j}\cdot r}{2}\left(u_{2,j+1}^n-2u_{2,j}^n+u_{2,j-1}^n\right)\\ u_{3,j}^n+\frac{a_{3,j}\cdot r}{2}\left(u_{3,j+1}^n-2u_{3,j}^n+u_{3,j-1}^n\right)\\ ⋮\\ u_{Nx-2,j}^n\text{+}\frac{a_{Nx-2,j}\cdot r}{2}\left(u_{Nx-2,j+1}^n-2u_{Nx-2,j}^n+u_{Nx-2,j-1}^n\right)\\ u_{Nx-1,j}^n+\frac{a_{Nx-1,j}\cdot r}{2}\left(u_{Nx-1,j+1}^n-2u_{Nx-1,j}^n+u_{Nx-1,j-1}^n\right)\end{array}\right]}\end{array}$ (3)

记等号左边分别为矩阵A和$u^{\ast }$ ，右边为$b_n$ ，通过求解$A{u}^{\ast }=b^n$ ，得到中间虚拟层的$u^{\ast }$ 。

第二步：$\frac{u_{i,j}^{n+1}-u_{i,j}^{\ast }}{\Delta t/2}=a\cdot \frac{1}{h^2}\left(u_{i+1,j}^{\ast }-2u_{i,j}^{\ast }+u_{i-1,j}^{\ast }+u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}\right)$ ，经化简，

$\begin{array}{l}\underset{A}{\left[\begin{array}{ccccc}1+a_{i,2}r& -a_{i,2}r/2& & & \\ -a_{i,3}r/2& 1+a_{i,3}r& -a_{i,3}r/2& & \\ & & \ddots & & \\ & & -a_{i,Nxy-2}r/2& 1+a_{i,Nxy-2}r& -a_{i,Nxy-2}r/2\\ & & & -a_{i,Nxy-1}r/2& 1+a_{i,Nxy-1}r\end{array}\right]}\times \underset{u^{n+1}}{\left[\begin{array}{c}{u}_{i,2}^{n+1}\\ u_{i,3}^{n+1}\\ ⋮\\ u_{i,Ny-2}^{n+1}\\ u_{i,Ny-1}^{n+1}\end{array}\right]}\\ =\underset{b^{\ast }}{\left[\begin{array}{c}{u}_{i,2}^{\ast }+\frac{a_{i,2}\cdot r}{2}\left(u_{i+1,2}^{\ast }-2u_{i,2}^{\ast }+u_{i-1,2}^{\ast }\right)\\ u_{i,3}^{\ast }+\frac{a_{i,3}\cdot r}{2}\left(u_{i+1,3}^{\ast }-2u_{i,3}^{\ast }+u_{i-1,3}^{\ast }\right)\\ ⋮\\ u_{i,Ny-2}^{\ast }\text{+}\frac{a_{i,Ny-2}\cdot r}{2}\left(u_{i+1,Ny-2}^{\ast }-2u_{i,Ny-2}^{\ast }+u_{i-1,Ny-2}^{\ast }\right)\\ u_{i,Ny-1}^{\ast }+\frac{a_{i,Ny-1}\cdot r}{2}\left(u_{i+1,Ny-1}^{\ast }-2u_{i,Ny-1}^{\ast }+u_{i-1,Ny-1}^{\ast }\right)\end{array}\right]}\end{array}$ (4)

2.2. 贝叶斯反演方法

贝叶斯反演的核心原理是在现实观测下，对未知参数的后验概率密度进行估计，并采用一些统计抽样技术，将该后验密度转化为参数的实际取值[8]。其思想可以利用如下的公式概括。对于偏微分方程参数反演模型：

$Y=f\left(X,E\right)$ (5)

其中，X为待估参数，E为观测误差，Y为观测值。这里并没有限制方程和参数的形式。该框架为统计反问题中更一般的，普适性更强的反演框架，能有效面对线性与非线性的参数反演问题。

$p\left(x|y\right)=\frac{p\left(x\right)\cdot p\left(y|x\right)}{p\left(y\right)}\propto p\left(x\right)\cdot p\left(y|x\right)$ (6)

式(6)通过贝叶斯定理，利用参数的先验信息$p\left(x\right)$ 和样本观测$p\left(y|x\right)$ 的似然信息，得到参数关于观测数据的后验分布信息。在得到参数后验概率$p\left(x|y\right)$ 的稳定分布后，可使用下面这两种方式表达反演结果。

① 最大后验估计：$x_{\text{MAP}}=\arg \underset{x\in {\rm X}}{\max }p\left(x|y\right)$ 。

② 后验均值估计：$x_{\text{CM}}=E\left(x|y\right)={\displaystyle {\int }_{\Omega }x\cdot p\left(x|y\right)\text{d}x}$ 。

对于最大后验估计，反问题的求解被同样转换成一个优化问题，可以通过极大似然或者拉格朗日方法求解。此外，一些基于马尔科夫链的抽样估计方法，例如Gibbs和Metropolis-Hastings，在这里同样适用。本文也将采取这种方法进行参数反演的可视化表达。

在贝叶斯反演框架中，先验分布的选择至关重要，它反映着研究者对于真实参数的预期分布。对于预期附近的参数，常赋予较高的概率取值；反之对于偏离预期较远的参数，则赋予较低的概率密度。因此将不同种类的参数先验融入反演模型，其参数反演结果也有所差异，或是收敛速度的快慢，又或者是收敛位置的偏差。

常见的先验分布有高斯分布、均匀分布、Markov随机场、脉冲先验和共轭先验等分布[8]，它们的适用范围各不相同，但均能反映出前文所提及的“远小近大”先验选取理念。高斯分布易于构造、通用性高，因此适合绝大多数参数反演场景；均匀分布相对更适合已知参数大致范围的情形，在后验分布统计抽样中可以减少运算量；后几种先验信息则被广泛应用在图像重建领域。

不论选择何种先验信息，它们都是对已有经验的总结，本文也将在之前各类相关问题的研究基础上，给定先验细节。问题P1中的单参数反演问题，在真参范围不明了的前提下，选择高斯分布是较为稳妥的，在文献[9] [10]中高斯先验的选取也保证了参数反演结果。本文将定义如下的先验信息：

$p_{\text{prior}}\left(a\right)\equiv \frac{1}{\sqrt{2\pi }\cdot {\sigma }_p}\cdot \exp \left\{-\lambda \frac{{\left(a-\overline{a}\right)}^2}{2{\sigma }_p^2}\right\}$(7)

其中$\overline{a}$ 和${\sigma }_p^2$ 分别为算法的先验均值与先验方差。在已有研究工作的基础上，经过多次数值实验后发现，较小的先验均值和方差即可使待估参数迅速收敛到真实位置，且后续似然函数对已有观测的拟合情况也更佳。因此，后续实验中$\overline{a}$ 和${\sigma }_p^2$ 分别取0.5和0.1²。$\lambda$ 为正则化参数，引入它可提升算法稳定性。

在推导a的似然函数之前，还需假定模型噪声$\epsilon$ 为符合正态分布$N\left(0,{\sigma }_{\epsilon }^2\right)$ 的加性噪声，这个噪声是由微分方程数值求解误差与观测误差共同造成的，此时附加的观测条件为$g^{\delta }\left(x,y\right)=u\left(x,y,T;a\right)+\epsilon$ 。则似然密度为，

$p_{\text{like}}\left(g|a\right)=\frac{1}{{\left(2\pi {\sigma }_{\epsilon }\right)}^{\frac{n}{2}}}\cdot \exp \left\{-\frac{{\left({u^a|}_{t=1}-g\right)}^{\text{T}}\left({u^a|}_{t=1}-g\right)}{2{\sigma }_{\epsilon }^2}\right\}$ (8)

式(8)中n的为样本数量，$u^a$ 为数值解$u_N\left(x,y,t;a\right)$ 在$t=1$ 时的温度，之后利用式(6)，有，

$\begin{array}{l}{p}_{\text{post}}\left(a|g\right)=\frac{p_{\text{prior}}\left(a\right)\cdot p_{\text{like}}\left(g|a\right)}{p\left(g\right)}\\ \propto \exp \left\{-\lambda \frac{{\left(a-\overline{a}\right)}^2}{2{\sigma }_p^2}\right\}\cdot \exp \left\{-\frac{{\left({u^a|}_{t=1}-g\right)}^{\text{T}}\left({u^a|}_{t=1}-g\right)}{2{\sigma }_{\epsilon }^2}\right\}\end{array}$(9)

此外，上式中的后验表达采用了正比的形式，这是考虑到相比起较大的系数项${\left(2\pi {\sigma }_{\epsilon }\right)}^{-\frac{n}{2}}$ ，在各类统计抽样技术里，指数函数中的均值与方差是更被关注的。图1绘制了两种误差下的后验密度，其均值集中在0.1附近，且越小的加性噪声对应越紧凑的密度函数。

Figure 1. Posterior probability density function of diffusion coefficient for inversion

图1. 扩散系数的后验概率密度函数

针对问题P2，此时待估参数从常系数变化为了时间依赖的函数型参数，此时仅依靠贝叶斯定理与统计抽样技术，将很难反演出符合要求的热交换系数。于是文章借用径向基函数拟合的思想，将函数$\sigma \left(t\right)$ 的推断转变为权重系数${\omega }_i$ 和基函数${\varphi }_i$ 的推断，即

$\sigma \left(t\right)={\displaystyle {\sum }_{i=1}^d{\omega }_i\cdot {\varphi }_i}$

d为基函数数量。且该拟合的基函数往往是固定的，一般为Gauss核或者多项式核，只需选择合适的系数即可。在上述技巧的支持下，P2的反演目标变成了一组系数$\omega =\left({\omega }_1,{\omega }_2,\cdots ,{\omega }_d\right)$ 的估计。对于多参数的Bayes反演而言，其整体流程仍与P1中的单参数反演一致，只不过先验与后验密度变换为了高维的联合密度。在定义这些密度之前，还需提出系数${\omega }_i$ 的独立性假设，则整体的先验可拆分成各分量${\omega }_i$ 的先验乘积。需要注意的是，这里选择高斯先验将不再合适，因为多个独立分量构造出的d维高斯联合密度函数取值往往较小，是由多个小概率密度连乘取得，在后续的M-H抽样中，会导致接受概率极小，Markov链较难更新，参数反演结果无法收敛。且探索d个不同分量的先验均值与方差的超参数组合，难度太高，很难找到合适的取值。又因为在较小区间上的径向基函数拟合，不同基函数对应的权重系数不会太大，故不妨对称性地假定${\omega }_i$ 上下界，将参数的预期取值限定在范围内。这样可大幅基提升算法反演速率，因此取${\omega }_i$ 的分布为$\left[-r,r\right]$ 上的均匀分布，则联合先验为，

$p_{\text{prior}}\left(\omega \right)={\left(2r\right)}^{-d}\cdot \lambda {\Vert \omega \Vert }_2^2$ (10)

h的似然概率密度仍采用高斯分布，即，

$p_{\text{like}}\left(h|\omega \right)={\left(2\pi {\sigma }_{\epsilon }\right)}^{-\frac{n}{2}}\cdot \exp \left\{-\frac{{\left({u^{\omega }|}_{x=x_0}-h\right)}^{\text{T}}\left({u^{\omega }|}_{x=x_0}-h\right)}{2{\sigma }_{\epsilon }^2}\right\}$ (11)

则其联合后验概率密度函数可被表述为，

$p_{\text{post}}\left(\omega |h\right)=c_{\lambda }\cdot \exp \left\{-\frac{{\left({u^{\omega }|}_{x=x_0}-h\right)}^{\text{T}}\left({u^{\omega }|}_{x=x_0}-h\right)}{2{\sigma }_{\epsilon }^2}\right\}$ (12)

其中$c_{\lambda }$ 为一包含正则算子和先验信息的常数。

最后给出该框架下，两种参数均能使用的Metropolis-Hastings反演算法，详细步骤见表1。

Table 1. M-H inversion algorithm for heat conduction equation

表1. 热传导方程M-H反演算法

Step1. 根据建议分布$q\left(a^{\ast }|a_k\right)$ 生成随机数$a^{\ast }$ ，并计算随机数$a^{\ast }$ 与参数当前状态$a_k$ 的建议分布函数值$q\left(a^{\ast }|a_k\right)$ 和$q\left(a_k|a^{\ast }\right)$ ，以及目标分布函数值$p\left(a^{\ast }\right)$ 和$p\left(a_k\right)$ ，假设$a_0$ 为一固定常数，在算法开始前给定。上述密度函数有如下定义， $\begin{array}{ll}q\left(a^{\ast }|a_k\right)=\frac{1}{\sqrt{2\pi }\cdot {\sigma }_k}\cdot \exp \left\{-\frac{{\left(a^{\ast }-a_k\right)}^2}{2{\sigma }_k^2}\right\},\hfill & {\sigma }_k=0.05a_k\hfill \end{array};$ $\begin{array}{cc}p\left(a^{\ast }\right)=p_{\text{post}}\left(a^{\ast }|g\right),& p\left(a_k\right)=p_{\text{post}}\left(a_k|g\right)\end{array}.$ Step2. 计算接受概率 $A\left(a^{\ast },a_k\right)=\min \left\{1,\frac{p\left(a^{\ast }\right)\cdot q\left(a_k|a^{\ast }\right)}{p\left(a_k\right)\cdot q\left(a^{\ast }|a_k\right)}\right\};$ Step3. 生成[0,1]上的均匀分布随机数r； Step4. 如果有$r\le A\left(a^{\ast },a_k\right)$ ，则令Markov链上的下一个状态$a_{k+1}$ 为$a^{\ast }$ ，反之则保持上一个的状态不变，即$a_{k+1}=a_k$ ； Step5. 令$k=k+1$ ，重复上述步骤，直至收集到足够长的Markov链，此时算法达到最大迭代次数$N_m$ 。

问题P1和P2的目标分布分别为后验分布$p\left(a|g\right)$ 与$p\left(\omega |h\right)$ ，建议分布则与各自先验保持一致，具体参数将在下一章数值实验中给出。

2.3. PINN反演方法

几乎所有的机器学习问题都可以被看做一个优化过程，即定义标签数据与预测数据间合适的误差，并在参数的可行集内使损失达到最小。故本节基于数值驱动的PINN反演方法与最优控制框架和Bayes框架相似，均需给出误差最小化的表达。这里仿照前两种框架，定义如下平方项损失，

${ℒ}_{data}=\frac{1}{M}{\displaystyle \sum _{j=1}^M{\left(u_N\left(x_0,y_0,t_j\right)-u_{obs}\left(x_0,y_0,t_j\right)\right)}^2}$ (13)

其中$u_N$ 为网络计算出的对应位置的数值解，为附加的观测信息。其形式与2.2节中观测信息的似然密度函数，以及最优控制框架中参数在$L^2\left(\Omega \right)$ 空间上的范数定义类似，但不同的是PINN并不需要反复计算正问题，而是将数值解和未知参数看做一个整体，同时给出反演结果。因此还需要正问题中涉及到的定解条件包括控制方程本身，和${ℒ}_{data}$ 一起共同限制PINN网络的损失，并使用自动微分技术，替代方程中的各项微分算子。则网络总损失${ℒ}_{total}$ 为，

${ℒ}_{total}={\omega }_e\cdot {ℒ}_{equation}+{\omega }_i\cdot {ℒ}_{initial}+{\omega }_b\cdot {ℒ}_{bound}+{\omega }_d\cdot {ℒ}_{data}$ (14)

其中前三项损失分别为控制方程损失、初值损失以及边界损失，最后一项为观测损失，这里采用的线性加权汇总总损失，各损失前的系数作为网络的超参数，可在实验中调整。

接下来给出问题P1各项损失的具体定义，

$\begin{array}{c}{\text{MSE}}_{P1}={\omega }_b\cdot {\text{MSE}}_B+{\omega }_i\cdot {\text{MSE}}_I+{\omega }_f\cdot {\text{MSE}}_f+{\omega }_d\cdot {\text{MSE}}_d\\ =\frac{2{\omega }_b}{N_x}{\displaystyle \sum _{n=1}^M{\displaystyle \sum _{j=1}^{N_y}{\left(u_N\left(x,y_j,t_n\right)|{}_{x=0,L}\right)}^2}}+\frac{2{\omega }_b}{N_y}{\displaystyle \sum _{n=1}^M{\displaystyle \sum _{i=1}^{N_x}{\left(u_N\left(x_i,y,t_n\right)|{}_{y=0,H}\right)}^2}}\\ +\frac{{\omega }_i}{N_x\cdot N_y}{\displaystyle \sum _{j=1}^{N_y}{\displaystyle \sum _{i=1}^{N_x}{\left(u_N\left(x_i,y_j,0\right)-\sin \left(\pi x_i\right)\text{sin}\left(\pi y_j\right)\right)}^2}}\\ +\frac{{\omega }_f}{N_f}{\displaystyle \sum _{i=1}^{N_f}{\left(u_t^N\left(x_i,y_i,t_i\right)-a\left(u_{xx}^N\left(x_i,y_i,t_i\right)+u_{yy}^N\left(x_i,y_i,t_i\right)\right)\right)}^2}\\ +\frac{{\omega }_d}{N_x\cdot N_y}{\displaystyle \sum _{j=1}^{N_y}{\displaystyle \sum _{i=1}^{N_x}{\left(u_N\left(x_i,y_j,1\right)-g\left(x_i,y_j\right)\right)}^2}}\end{array}$ (15)

在给定损失之后，仍需给出具体的网络结构作为该反演任务的模型。对于问题P1中较为简单的输入特征(一些随机离散点)而言，在奥卡姆剃刀原理的支持下，使用灵活性高且计算效率更快的前馈型全连接网络具有更大的优势。此外，在以往学者的研究基础上发现[5] [11] [12]，得益于PINN中嵌入的物理信息，仅需少量的隐藏层即可达成反演目标。待估参数a则可作为神经网络的固定参数，同神经元的权值与偏置一起训练。

激活函数在训练中将直接影响网络的输出结果和参数更新梯度，如若选择不当，容易造成梯度消失或者爆炸的现象。不同的激活函数，其值域分布，梯度特性略有区别，需要根据输出需求来具体选择。在偏微分方程数值计算领域中，文献[13] [14]引入了一些自适应激活函数，以提升PINN方法的鲁棒性；文献[12]则从数值实验的角度给出Tanh、Sigmoid、ISRU和Soft Plus四类函数的反演表现，其中Tanh函数表现最佳。其[−1, 1]的值域分布，也不会使后续神经元发生输入偏移或输出为0的情况，进而降低收敛速度与精度，因此本文将选择Tanh函数作为激活函数。

随机梯度下降方法，作为最常见的参数训练算法，其实现简单，计算效率高，但其固定的学习率使得整体训练速度较慢，且会出现陷入局部最优的情况。因此在后续的研究中，学者们陆续提出一些学习率自适应调整算法以及梯度动量加速算法，来避免上述情况。Adam则结合了两种参数优化的优点，其训练时间更短，且训练结果的鲁棒性更强，是一种常见优化方法。本文的PINN反演算法也将采取该优化方式。

对于问题P2，将定义如下损失，

$\begin{array}{c}{\text{MSE}}_{\text{P2}}={\omega }_b\cdot {\text{MSE}}_B+{\omega }_i\cdot {\text{MSE}}_I+{\omega }_f\cdot {\text{MSE}}_f+{\omega }_d\cdot {\text{MSE}}_d\\ =\frac{2{\omega }_b}{N_0}{\displaystyle \sum _{n=1}^{N_0}{\left(u_N\left(0,t_n\right)\right)}^2}+\frac{2{\omega }_b}{N_1}{\displaystyle \sum _{n=1}^{N_1}{\left(u_x^N\left(1,t_n\right)+{\sigma }_N\left(t_n\right)\cdot u_N\left(1,t_n\right)\right)}^2}\\ +\frac{{\omega }_i}{N_i}{\displaystyle \sum _{i=1}^{N_i}{\left(u_N\left(x_i,0\right)-\sin \left(\pi x_i/L\right)\right)}^2}\\ +\frac{{\omega }_f}{N_f}{\displaystyle \sum _{i=1}^{N_f}{\left(u_t^N\left(x_i,t_i\right)-u_{xx}^N\left(x_i,t_i\right)-{\text{e}}^{-t_i}\sin \left(\pi x_i/L\right)\right)}^2}\\ +\frac{{\omega }_d}{N_{x0}}{\displaystyle \sum _{n=1}^{x_0}{\left(u_N\left(x_0,t_n\right)-h\left(t_n\right)\right)}^2}\end{array}$ (16)

注意到，这里面对着2.2节中相同的问题，即如何表达函数型参数$\sigma \left(t\right)$ 。在上一节中，选取了基函数模拟等价替换的思路，这里理论上也可使用该方法，将$\sigma \left(t\right)$ 的反演转换为一组d维系数的反演。除此之外，PINN这类深度网络具其余反演理论不具备的优势，即其高度灵活的模型结构，使它在面对复杂任务或者特殊输出需求时，仅调整网络结构即可。这里仅需在最终的输出层上在添加一维输出作为$\sigma \left(t\right)$ 的反演结果，且该维度的输出仍然与(16)中的损失保持一致，只不过要限制它的输入仅包含时间变量，故该问题对应的网络结构被修改为图2中的形式。

Figure 2. PINN network structure corresponding to problem P2

图2. 问题P2对应的PINN网络结构

问题P2的激活函数与参数优化算法将与P1保持一致，对于学习率、隐藏层层数、神经元个数等超参数的选取，将在文献[11] [12]的基础上，结合多次实验结果给出反演表现相对最佳的取值。

至此已经搭建好两类统计反演算法，下一章将给出这两种方法的数值反演表现。

3. 数值反演实验

针对问题P1，设定扩散系数a的真实取值为0.1，先验均值与方程分别为0.5和0.1，算法初值设定为0.5，在第三章中，已经给出了两种观测误差下的后验均值分布范围，图3给出M-H抽样的可视化进程。

Figure 3. Parameter update process of M-H inversion algorithm for problem P1

图3. 问题P1的M-H反演算法参数更新过程

其反演结果与相对精度见表2。

Table 2. M-H algorithm inversion results for problem P1

表2. 问题P1的M-H算法反演结果

观测误差	反演结果	相对精度	首次平稳步长
${\sigma }_{\epsilon }=0.01$	0.1001	0.10%	85
${\sigma }_{\epsilon }=0.05$	0.1003	0.30%	145

问题P1对应的PINN算法超参数取值情况如下，前馈网络具有5层隐藏层，每层32个神经元，dropout取0.05，训练时每个批次均重新随机抽样，初值与四条边界上各选取2000个采样点，控制方程选取10,000个采样点，为保证算法收敛，需提高观测位置采样点数量，这里取5000个，学习率取10⁻³，激活函数为tanh。其反演结果见表3。

Table 3. PINN algorithm inversion results for problem P1

表3. 问题P1的PINN算法反演结果

观测误差	反演结果	相对精度	收敛所需批次
$\delta =0$	0.1001	0.10%	7042
$\delta =3\%$	0.1010	1.00%	15,269
$\delta =10\%$	0.1007	0.70%	9704

对于问题P2中的热交换系数$\sigma \left(t\right)$ ，这里取真值为$4\sin \left(\pi t\right)/3+\cos \left(\pi t\right)/3$ ，若使用径向基函数模拟，当t被限制在较小区间时，d取3便可满足模拟需求，M-H反演算法中的先验均匀分布为[−10, 10]，算法对初值不敏感，任意较小的初值在经历多次迭代后，皆可使结果收敛在真实值附近。建议分布的抽样变化步长为1，是先验范围的5%。反演结果如图4所示。

Figure 4. M-H algorithm inversion results for problem P2

图4. 问题P2的M-H算法反演结果

其中红色的虚线为真实$\sigma \left(t\right)$ 对应的基函数权值，${\omega }_{true}=\left(-5.4462,8.5574,-2.9136\right)$ 上下波动1%后，绘制出来的。两种误差下的反演结果，均落在预测区间内。

对于问题P2的变体PINN，隐层设定为6层，每层32个神经元，为加快收敛进度，学习率取0.005，dropout则适当调大成0.1。又因该问题的反演目标在右边界上，因此取右边界采样点为5000，左边界与初值设为2000个，观测采样点仍与P1保持一致。其反演效果见图5。

Figure 5. Inversion results of variant PINN algorithm for problem P2

图5. 问题P2的变体PINN算法反演结果

Table 4. Inversion errors of two algorithms for problem P2

表4. 问题P2的两类算法的反演误差

观测误差	反演结果	相对精度
M-H抽样反演
${\sigma }_{\epsilon }=0.01$	0.0030	4.55%
${\sigma }_{\epsilon }=0.05$	0.0034	4.62%
变体PINN反演
$\delta =0$	0.0016	3.42%
$\delta =3\%$	0.0021	3.63%
$\delta =10\%$	0.0019	3.62%

在文章的最后，表4还给出了问题P2中两类方法的均方误差与相对误差表现，整体看来两种方法均能较好地反演出参数实际情况，在高噪声的情况下，反演表现也会优于一般的最优控制框架。横向对比两类方法，后者基于变体PINN在精度上更为出色，但训练时长会略久于M-H抽样反演。但这种情况在高维空间中却有所反转。

4. 总结

本文分别研究了热方程不同位置上两类参数的反演问题，给出两类统计反演算法及其反演表现，相比传统最优控制框架反演，这两种算法不需频繁地推导伴随方程与泛函梯度信息，均能高效地完成未知参数的估计。此外，这两种方法的鲁棒性更强，在高噪情形下的反演结果相对光滑[15]。