三维稳态热传导方程参数反演问题的PINNs解法

PINNs Method for Parameter Inversion Problem of Three-Dimensional Steady-State Heat Conduction Equation

摘要: 热传导反问题广泛应用于地球科学、医学等领域，本文基于内嵌物理机理神经网络(PINNs)研究了一类三维有界域上热传导方程参数反演问题的数值解法。构建三维参数辨识模型，提出基于PINNs求解立方体区域三维稳态热传导方程参数反演问题的算法。通过同步反演热传导系数与整体温度值，实现热物性参数的高精度识别。在含噪声数据条件下，热传导系数反演的误差精度控制在10⁻⁴数量级，证实了算法对病态问题具有良好的抗噪性。

Abstract: Inverse heat conduction problem is widely used in earth science, medicine and other fields. In this paper, based on embedded physical mechanism neural network (PINNs), a numerical method for parameter inversion of heat conduction equation in three-dimensional bounded domain is studied. A three-dimensional parameter identification model is established, and an algorithm based on PINNs is proposed to solve the parameter inversion problem of three-dimensional steady-state heat conduction equation in cube region. Through simultaneous inversion of heat transfer coefficient and global temperature value, high precision identification of thermal physical property parameters is realized. Under the condition of noisy data, the error accuracy of heat transfer coefficient inversion is controlled at the order of 10⁻⁴, which proves that the algorithm has good anti-noise performance for ill-conditioned problems.

1. 引言

热传导反问题因其广泛的工程应用背景及固有的不适定性，长期以来受到研究者们的高度关注。近年来，传统数值方法、软计算方法以及深度学习方法的不断发展，为热传导反问题的有效求解提供了丰富的途径和方法。

传统数值方法主要包括有限元法、有限差分法及边界元法等方法。这类方法理论成熟，但在处理不适定问题时需引入正则化技术以增强稳定性。薛齐文等人[1]采用Tikhonov正则化方法，结合加权Bregman距离函数，有效解决了非线性热传导反问题中多个参数的辨识问题，显著提高了计算的稳定性与抗噪能力。周焕林等人[2]利用边界元法反演二维瞬态热传导问题的温度相关导热系数，展示了传统数值方法处理非线性问题的有效性。Mohebbi等人[3]则针对空间与温度相关的导热系数估计问题，提出了有效的数值解决方案。

随着计算机技术的不断进步，人工神经网络(ANN)、支持向量机以及模糊神经网络等软计算方法逐渐应用于热传导领域。Yang等人[4]利用ANN精确建模了超临界乙烷的导热系数，验证了ANN在热物性预测中的高效性。吴自库等人[5]则采用最小二乘支持向量机(LS-SVM)高效求解了一维热传导方程中的热源反问题，获得了较好的稳定性与计算精度。Haznedar等人[6]通过模拟退火算法优化ANFIS结构，实现了动态系统的参数辨识。此外，Psichogios等人[7]的神经网络与物理原理融合的早期研究和Lagaris等人[8]提出的神经网络求解偏微分方程的工作，为软计算方法的应用奠定了基础。

深度学习的迅速兴起，尤其是物理信息神经网络(PINNs)的提出，为热传导反问题的研究注入了新的活力。Raissi等人[9]首次提出PINNs框架，通过引入物理约束成功解决了多个复杂非线性偏微分方程，开创了深度学习与物理知识相结合的新研究领域。随后，Cai等人[10]进一步将PINNs应用于热传导问题，处理了相变Stefan问题等多个复杂的实际工程问题。Rasht-Behesht等人[11]则利用PINNs成功求解了波传播与地震波反演问题，进一步验证了该方法的适用性。陈豪龙等人[12]基于PINNs深入研究了非线性瞬态热传导问题的参数辨识，讨论了网络结构及数据配置对结果的影响，显示出PINNs的良好泛化性能。Lu等人[13]开发的DeepONet模型，通过非线性算子学习，有效提高了偏微分方程求解的泛化性能与误差收敛性。此外，Musgrave等人[14]开发的傅里叶域PINNs方法，显著提升了稀疏数据下的求解精度。Borate等人[15]通过PINNs从实验数据中反演了非线性摩擦参数，证实了该方法在实际数据分析中的可靠性。Zhang等人[16]则成功应用深度神经网络解决了热传导方程的逆源项问题。Chanda等人[17]采用非侵入式方法，结合神经网络高效识别了热接触电导。Liu等人[18]则针对二维时间分数阶热传导方程提出了有效的PINNs逆问题求解方法，拓展了PINNs的适用范围。

传统数值方法虽然理论成熟，但在抗噪性和计算复杂性方面存在局限。以PINNs为代表的深度学习方法，通过物理信息融合和非线性拟合能力，展现了更为广阔的应用前景。

2. PINNs方法求解稳态热传导方程参数反演问题

2.1. 三维稳态热传导方程参数反演问题模型

稳态热传导方程在多个工程领域有广泛应用，主要用于分析和优化与热传导相关的稳态问题，本文研究有界区域$\Omega \left(\Omega \subset {ℝ}^3\right)$ 上的稳态热传导方程参数反演问题，数学模型简述如下：

$\{\begin{array}{l}-\nabla \cdot \left(\lambda \nabla u\right)\left(x,y,z\right)=f\left(x,y,z\right),\left(x,y,z\right)\in \Omega ,\\ {u\left(x,y,z\right)|}_S=g\left(x,y,z\right),\left(x,y,z\right)\in S.\end{array}$ (1)

其中，$\nabla$ 为空间梯度算子，$\nabla u$ 为$u$ 的梯度；$u\left(x,y,z\right)$ 表示区域内温度，$g\left(x,y,z\right)$ 为边界$S=\partial \Omega$ 上的Dirichlet边界条件，$\lambda$ 表示热传导系数，可为常数或与位置相关的函数。

已知$f\left(x,y,z\right),g\left(x,y,z\right)$ ，识别未知参数$\lambda$ ，并预测区域$\Omega$ 上的温度值。由于反问题的不适定性，需要额外附加观测数据，本文考虑附加局部区域${\Omega }_D\left({\Omega }_D\subset \Omega \right)$ 上的温度测量值，对$u_D\left(x,y,z\right)=u\left(x,y,z\right)\left(\left(x,y,z\right)\in \Omega \right)$ 进行数值求解。

2.2. 三维稳态热传导方程参数反演问题PINNs方法求解

基于神经网络方法，可以把附加了内部温度测量条件的问题(1)转化为以下的最优控制问题[19]：

$\{\begin{array}{l}\underset{w,b}{\min }\text{ C}\left(u\left(w,b\right);w,b\right),\\ C\left(u\right)={\Vert \Delta u-f\Vert }_{L_2(\Omega )}^2+{\Vert u-g_1\Vert }_{L_2(S_1)}^2+{\Vert u-u_D\Vert }_{L_2(D)，}^2\\ u=\text{NET}\left(x,y,z;w,b\right).\end{array}$ (2)

其中，$g_1$ 为已知边界条件，$u_D$ 为内部局部区域测量值，$w,b$ 为神经网络中的参数，通过神经网络方法将其转化为参数优化问题。

具体来说，利用一个物理信息神经网络函数$\widehat{u}\left(x,y,z;\theta \right)$ 逼近未知函数$u\left(x,y,z\right)$ ，从而得出未知热传导系数$\lambda$ ，其中$\theta$ 为神经网络参数，其中参数$\theta$ 包含了权重向量$w$ 和偏置$b$ ，采用L-BFGS算法[20]更新网络参数，使损失函数达到最小，求解过程见图1。

Figure 1. Technical route of solving PDE based on PINNs

图1. 基于PINNs求解PDE的技术路线

建立神经网络模型，输入层3个，输出层1个，隐藏层$L-1$ 层，每个隐藏层有$N$ 个神经元，如图2所示。

Figure 2. Neural network structure diagram of inverse problem of heat conduction equation

图2. 热传导方程反问题神经网络结构图

定义偏微分方程残差${ℒ}_f\left(\theta \right)$ 为：

${ℒ}_f\left(\theta \right)=\frac{1}{N_f}{\Vert -\nabla \cdot \left(\lambda \nabla \widehat{u}\right)\left(x,y,z;\theta \right)-f\Vert }_{L^2(\Omega )}^2.$ (3)

边界条件残差${ℒ}_{u_b}\left(\theta \right)$ 为：

${ℒ}_{u_b}\left(\theta \right)=\frac{1}{N_b}{\Vert \widehat{u}\left(x,y,z;\theta \right)-g_i\Vert }_{L^2(S_i)}^2.$ (4)

观测数据残差${ℒ}_{{}_{data}}^D\left(\theta \right)$ 为：

${ℒ}_{{}_{data}}^D\left(\theta \right)=\frac{1}{N_{data}}{\Vert \widehat{u}\left(x,y,z;\theta \right)-u_D\Vert }_{L^2(D)}^2.$ (5)

其中，$N_f,N_b,N_{data}$ 分别为方程配置点、边界点，以及观测数据样本集大小，偏导数由自动微分得到。

定义损失函数为：

$ℒ{\left(\theta \right)}_D={\omega }_1{ℒ}_{u_b}\left(\theta \right)+{\omega }_2{ℒ}_f\left(\theta \right)+{\omega }_3{ℒ}_{{}_{data}}^D\left(\theta \right).$ (6)

其中，${\omega }_i\left(i=1,2,3\right)$ 为权重且${\omega }_i>0$ 。损失函数的训练集分为三个部分，分别为边界条件的采样点、偏微分方程残差进行的配置点和观测条件的采样点。其中配置点的选取方式有两种，分别为类网格点选取方式以及伪随机离散取点方式，如图3所示，本文数值实验采用的是伪随机离散点选取方式。

其中损失函数通常都是非线性或者非凸的，利用优化算法(L-BFGS算法等)来求解最小化损失函数，最终得到近似解。

参数反演问题的神经网络求解算法如表1所示。

(a) 类网格取点 (b) 伪随机离散取点

Figure 3. Configuration point selection method

图3. 配置点选取方式

Table 1. Algorithm for solving parameter inversion problems using PINNs

表1. PINNs解参数反演问题算法

算法1：PINNs解参数反演问题算法

Step 1：取边界采样点$N_b$ ，方程配置点$N_f$ 和观测采样点$N_{data}$ 作为训练集；

Step 2：构造输出层有1个神经元的网络结构；

Step 3：初始化参数$\theta$ ，即神经网络的权重$w$ 和偏差$b$ ；

Step 4：利用自动微分技术得到物理约束方程${ℒ}_f\left(\theta \right)$ 、边界残差${ℒ}_{u_b}\left(\theta \right)$ ；

Step 5：再与边界观测残差${ℒ}_{{}_{data}}^D\left(\theta \right)$ 构成损失函数$ℒ{\left(\theta \right)}_D$ ；

Step 6：利用优化算法对网络参数$\theta$ 更新, 最小化损失函数$ℒ{\left(\theta \right)}_D$ ；

Step 7：直到损失函数$ℒ{\left(\theta \right)}_D$ 达到一定精度停止得到最终的近似解$\widehat{u}\left(x,y,z;\theta \right)$ 。

数值算例：求解区域为立方体区域。采用添加随机误差的观测数据模拟计算，表达式为：

$u^{\delta }=u\times \left(1+\delta \times r\right).$ (7)

其中，$\delta$ 为噪声水平，$r$ 为$\left[-1,1\right]$ 区间上的随机数。

数值精度由相对误差来判断，其计算公式如下：

$E_r\left(u\right)=\frac{{\Vert u-\widehat{u}\Vert }_2}{{\Vert u\Vert }_2}.$ (8)

本文所有数值算例的代码运行的版本为TensorFlow 1.15.0，具体实现主要依赖NumPy和TensorFlow两个库。

3. 数值实验

在一个立方体区域$\overline{\Omega }=\left\{\left(x,y,z\right)|\left[0,1\right]\times \left[0,1\right]\times \left[0,1\right]\right\}$ 内，考虑以下热传导方程分段常参数反演问题数值算例：

$\{\begin{array}{l}-\nabla \cdot \left(\lambda \nabla u\left(x,y,z\right)\right)=-\lambda \left({\pi }^2{\text{e}}^{\pi x}+4\right),\left(x,y,z\right)\in \Omega ,\\ {u|}_{S_1}=g_1=u\left(0,y,z\right),{u|}_{S_2}=g_2=u\left(1,y,z\right),y\in \left[0,1\right],z\in \left[0,1\right],\\ {u|}_{S_3}=g_3=u\left(x,0,z\right),{u|}_{S_4}=g_4=u\left(x,1,z\right),x\in \left[0,1\right],z\in \left[0,1\right],\\ {u|}_{S_5}=g_5=u\left(x,y,0\right),{u|}_{S_6}=g_6=u\left(x,y,1\right),x\in \left[0,1\right],y\in \left[0,1\right],\\ {u|}_D=u\left(\frac{1}{2},y,z\right)={\text{e}}^{\frac{1}{2}\pi }+1+y^2+z^2，y\in \left[0,1\right],z\in \left[0,1\right].\end{array}$ (9)

其中，$\lambda \{\begin{array}{l}{\lambda }_1=1,x\in \left(0,0.5\right]\\ {\lambda }_2=2,x\in \left(0.5,1\right)\end{array}$ ，温度真解为$u\left(x,y,z\right)=e^{\pi x}+2x+y^2+z^2$ ，其中边界条件$g_i\left(i=1,2,3,4,5,6\right)$ 已知，$S_i\left(i=1,2,3,4,5,6\right)$ 是指立方体区域的六个面，立方体区域示意图见图4，${u|}_D$ 为附加内部观测条件。

Figure 4. Schematic diagram of a 3D cube

图4. 三维立方体示意图

建立神经网络模型，以反演出未知参数${\lambda }_1,{\lambda }_2$ 。

网络的隐藏层数分别使用2、3、4、6进行测试，各隐藏层网络的神经元个数均为30，激活函数分别使用ReLU函数、Logistic函数和Tanh函数。得到的整体温度近似解误差$e_u$ 和参数近似解误差$e_{{\lambda }_1}$ 和$e_{{\lambda }_2}$ 如表2所示。

Table 2. Experimental results of different activation functions at different layer numbers

表2. 不同激活函数在不同层数下的试验结果

隐藏网络层数	激活函数	$e_u$	$e_{{\lambda }_1}$	$e_{{\lambda }_2}$
2	ReLU	7.29E−3	5.02E−1	5.12E−1
2	Logistic	1.79E−1	9.65E−1	5.33E−1
2	Tanh	2.63E−3	8.73E−4	2.82E−3
3	ReLU	1.13E−2	5.06E−1	5.01E−1
3	Logistic	1.78E−1	8.71E−1	5.43E−1
3	Tanh	1.08E−4	2.39E−4	2.45E−4
4	ReLU	1.18E−2	5.36E−1	5.21E−1
4	Logistic	5.28E−1	3.58E−1	6.17E−1
4	Tanh	1.12E−4	2.35E−4	2.21E−4
6	ReLU	1.03E−2	5.24E−1	5.27E−1
6	Logistic	6.21E−1	4.95E−1	5.16E−1
6	Tanh	1.44E−4	3.56E−4	5.39E−4

由此可以推测，拟合效果虽然随着层数的增加，有一定的提高，但并不会随着层数的增加而无限提高，由表2可以看出，激活函数ReLU函数和Logistic函数在此问题中的效果远远不如Tanh函数；综合考虑，本文的激活函数选用Tanh函数。

确定了激活函数的选取，接下来探索改变网络层数和神经元的数量对PINNs性能的影响，表3展示了在不同网络层数和不同神经元个数下得到的整体温度近似解的误差$e_u$ 。

Table 3. Prediction errors of PINNs with different network layers and neurons for $e_u$

表3. 不同网络层数和神经元的PINNs预测误差$e_u$

网络层数	神经元个数
	10	20	40	50	60
2	1.31E−2	1.07E−3	1.48E−1	1.51E−1	1.49E−1
3	2.69E−1	2.48E−1	1.90E−4	1.47E−4	1.56E−1
4	4.11E−2	1.99E−2	2.27E−4	9.67E−5	1.46E−4
5	8.36E−1	1.72E−1	2.23E−4	1.88E−1	6.97E−1

从表中数据可知，当神经元数为40或50，网络层数为3层或4层时，观察到其对误差值的精度影响不大；由表可以观察到4层50个神经元得到的结果效果最好，出于计算精度和时间成本考虑，本数值实验选用4层50个神经元的网络结构来做数值算例。

先建立包含输入、输出和包含50个神经元的4个隐藏层的神经网络，优化算法采用L-BFGS算法，并且取Tanh函数作为激活函数。共选取2800个点的训练数据：$N_b=400\times 6,N_f=200,N_{data}=200$ 。损失函数为公式(6)，所选取的权重分别为${\omega }_1=7,{\omega }_2=2,{\omega }_3=1$ ，损失函数权重的选取也是一个值得考虑的问题，其中多任务似然损失平衡方法具有较好的数学基础以及理论依据[21]，但是会出现数值不稳定。因此，用经验选取损失函数的权重。而边界条件残差${ℒ}_{u_b}\left(\theta \right)$ 、偏微分方程残差${ℒ}_f\left(\theta \right)$ 以及观测数据残差${ℒ}_{{}_{data}}^D\left(\theta \right)$ 分别为：

$\begin{array}{l}{ℒ}_{u_b}\left(\theta \right)=\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(0,y_i,z_i;\theta \right)-u_1\left(0,y_i,z_i\right)\right|}^2}+\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(1,y_i,z_i;\theta \right)-u\left(1,y_i,z_i\right)\right|}^2}\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}+\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(x_i,0,z_i;\theta \right)-u\left(x_i,0,z_i\right)\right|}^2}+\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(x_i,1,z_i;\theta \right)-u\left(x_i,1,z_i\right)\right|}^2}\\ \begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}+\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(x_i,y_i,0;\theta \right)-u\left(x_i,y_i,0\right)\right|}^2}+\frac{1}{N_b}{\displaystyle \sum _{i=1}^{N_b}{\left|\widehat{u}\left(x_i,y_i,1;\theta \right)-u\left(x_i,y_i,1\right)\right|}^2}.\end{array}$ (10)

${ℒ}_f\left(\theta \right)=\frac{1}{N_f}{{\displaystyle \sum _{i=1}^{N_f}\left|-\left[\frac{\partial }{\partial x_i}\left(\lambda \frac{\partial \widehat{u}\left(p_i;\theta \right)}{\partial x_i}\right)+\frac{\partial }{\partial y_i}\left(\lambda \frac{\partial \widehat{u}\left(p_i;\theta \right)}{\partial y_i}\right)+\frac{\partial }{\partial z_i}\left(\lambda \frac{\partial \widehat{u}\left(p_i;\theta \right)}{\partial z_i}\right)\right]-f\left(p_i\right)\right|}}^2.$ (11)

${ℒ}_{data}^D\left(\theta \right)=\frac{1}{N_{data}}{\displaystyle \sum _{i=1}^{N_{data}}{\left|\widehat{u}\left(\frac{1}{2},y_i,z_i;\theta \right)-u\left(\frac{1}{2},y_i,z_i;\theta \right)\right|}^2}.$ (12)

其中，$p_i=\left(x_i,y_i,z_i\right)$ 表示在区域内部随机选取的点。为了体现该方法的抗噪性，将带有不同噪声水平的附加条件放入训练，得到损失函数的变化如图5所示。

Figure 5. Loss function iteration diagram of piecewise constant parameter inversion problem

图5. 分段常参数反演问题损失函数迭代图

图6为立方体区域利用不同噪声的内部观测数据得到的近似解$\widehat{u}$ 和精确解$u$ 的绝对误差对比图，(a)为$\delta =0$ 时近似解的误差，(b)为$\delta =1\%$ 时近似解的误差，(c)为$\delta =5\%$ 时近似解的误差，(d)为$\delta =10\%$ 时近似解的误差。

(a) 无噪声时绝对误差 (b) 1%噪声时绝对误差

(c) 5%噪声时绝对误差 (d) 10%噪声时绝对误差

Figure 6. Absolute error of approximate solution of piecewise constant parameter inversion problem

图6. 分段常参数反演问题近似解的绝对误差

Table 4. Relative error of the results for the piecewise constant parameter inversion problem

表4. 分段常参数反演问题结果的相对误差

	$\widehat{u}$ 与$u$ 的相对误差	${\widehat{\lambda }}_1$ 与$\lambda$ 相对误差	${\widehat{\lambda }}_2$ 与$\lambda$ 相对误差
无噪声	9.670638E−5	1.01E−4	9.43E−6
1%噪声	9.394715E−5	5.47E−5	1.16E−5
5%噪声	1.971324E−4	8.37E−4	8.18E−4
10%噪声	1.601566E−4	6.86E−4	6.94E−4

如表4所示，对于分段常参数反演问题，反演得到的近似解$\widehat{u}$ 在无噪声时相对误差精度为${10}^{-5}$ ，而有噪声时的相对误差精度为${10}^{-4}$ 。预测的参数的相对误差精度基本稳定在${10}^{-5}$ 左右。不管是近似解$\widehat{u}$ 还是参数$\widehat{\lambda }$ ，利用PINNs方法求解立方体域三维稳态热传导参数反演问题精度较好。

4. 结语

利用PINNs方法求解三维稳态热传导方程参数反演问题，将此问题转化为优化问题，并利用观测数据、边界条件以及偏微分方程约束条件构造损失函数，再通过优化算法将其极小化。针对参数反演的数值算例的结果显示，在给定噪声水平的情况下，所反演得到的参数近似解精度可达到${10}^{-5}$ 到${10}^{-4}$ 数量级，表明该方法对噪声具有较低的敏感性，并在数值实验中展现出良好的稳定性。