基于數(shù)據(jù)驅(qū)動(dòng)模型求解熱傳導(dǎo)反問(wèn)題

陳豪龍, 柳占立

(清華大學(xué) 航天航空學(xué)院 應(yīng)用力學(xué)教育部實(shí)驗(yàn)室，北京 100084)

1 引 言

熱傳導(dǎo)反問(wèn)題是采用物體內(nèi)部或表面的溫度信息來(lái)反求物體的某些未知信息，如熱物性參數(shù)、幾何形狀、邊界條件以及熱源項(xiàng)等。熱傳導(dǎo)反問(wèn)題在化工、材料、航天、無(wú)損探傷、動(dòng)力工程、冶金工程和生物工程等領(lǐng)域有著廣泛的應(yīng)用。

熱傳導(dǎo)反問(wèn)題屬于典型的不適定性問(wèn)題[1,2]。通常求解熱傳導(dǎo)反問(wèn)題的方法有Tikhonov正則化[3]、梯度類算法[4,5]和元啟發(fā)式優(yōu)化算法[6,7]等。正則化方法用一組與原不適定問(wèn)題相鄰近的適定問(wèn)題的解去逼近，從而降低反問(wèn)題的不適定性；而迭代類算法是通過(guò)搜索待演變量，找到使計(jì)算溫度和測(cè)量溫度分布最接近的變量作為待求的解。但是迭代類算法在求解過(guò)程中需要不斷更新搜索變量，特別是對(duì)于邊界幾何形狀識(shí)別問(wèn)題，需要對(duì)域內(nèi)網(wǎng)格重剖分。文獻(xiàn)[8]提出了一種基于參數(shù)轉(zhuǎn)化的思想，將邊界幾何形狀識(shí)別問(wèn)題轉(zhuǎn)化為有效導(dǎo)熱系數(shù)的識(shí)別問(wèn)題，從而避免了域內(nèi)網(wǎng)格的重剖分，提高了計(jì)算效率。但是，這種方法結(jié)合迭代類算法在求解反問(wèn)題過(guò)程中，并沒(méi)有使用歷史數(shù)據(jù)信息，對(duì)于不同問(wèn)題，需要重復(fù)計(jì)算。

近年來(lái)，隨著計(jì)算機(jī)科學(xué)的發(fā)展，機(jī)器學(xué)習(xí)和數(shù)據(jù)科學(xué)取得了長(zhǎng)足進(jìn)步[9]。借助數(shù)據(jù)驅(qū)動(dòng)模型，能夠降低與本構(gòu)模型相關(guān)的人力和時(shí)間成本，極大地拓展工程師分析力學(xué)的能力和進(jìn)行數(shù)值計(jì)算的效率，特別是針對(duì)特定的力學(xué)問(wèn)題[10,11]。

本文采用數(shù)據(jù)驅(qū)動(dòng)模型求解熱傳導(dǎo)反問(wèn)題。首先，使用有限元法求解熱傳導(dǎo)正問(wèn)題，并利用隨機(jī)模型產(chǎn)生訓(xùn)練數(shù)據(jù)。在此基礎(chǔ)上，運(yùn)用有效導(dǎo)熱系數(shù)轉(zhuǎn)化的思想，建立深度學(xué)習(xí)模型，求解了測(cè)點(diǎn)溫度與有效導(dǎo)熱系數(shù)之間的抽象映射關(guān)系，進(jìn)而識(shí)別了管道內(nèi)壁的幾何形狀。并進(jìn)一步擴(kuò)展深度學(xué)習(xí)模型，識(shí)別了皮膚腫瘤的生長(zhǎng)參數(shù)，同時(shí)討論了不同測(cè)量誤差對(duì)計(jì)算結(jié)果的影響。

2 數(shù)據(jù)驅(qū)動(dòng)基本理論

2.1 人工神經(jīng)網(wǎng)絡(luò)ANN

ANN(Artificial Neural Network)是一種受自然啟發(fā)的模擬生物神經(jīng)網(wǎng)絡(luò)模型的方法，因其具有優(yōu)良的擬合性能，逐漸成為機(jī)器學(xué)習(xí)的主流技術(shù)之一。首先簡(jiǎn)要介紹ANN的基本原理。

ANN的核心思想是構(gòu)建一個(gè)含有大量待定參數(shù)的網(wǎng)絡(luò)系統(tǒng)，并通過(guò)梯度下降法尋求最佳的網(wǎng)絡(luò)參數(shù)，從而使誤差函數(shù)最小。ANN中包含許多神經(jīng)元。一組神經(jīng)元中的每個(gè)神經(jīng)元都與另一組神經(jīng)元相連。圖1是ANN中單個(gè)神經(jīng)元的數(shù)學(xué)模型。xi為輸入元素，ωi j為權(quán)重，bj為偏差，則輸入與權(quán)重的乘積X可以表示為

圖1 神經(jīng)元的數(shù)學(xué)模型

(1)

然后通過(guò)激活函數(shù)f產(chǎn)生輸出

(2)

常用的非線性激活函數(shù)包括ReLU和Tansig。

如圖2所示，ANN的基本結(jié)構(gòu)由三種類型的層組成，分別為輸入層、隱藏層和輸出層。ANN可以不需要任何額外的信息，尋找輸入層和輸出層之間的抽象映射關(guān)系，這個(gè)過(guò)程叫做訓(xùn)練。在訓(xùn)練過(guò)程中，根據(jù)相應(yīng)的輸入和輸出(訓(xùn)練集)來(lái)確定網(wǎng)絡(luò)中不同層次神經(jīng)元的權(quán)重ωi j和偏差bj等網(wǎng)絡(luò)參數(shù)。重復(fù)訓(xùn)練過(guò)程，直到損失函數(shù)值下降到一個(gè)可接受的范圍，從而確定網(wǎng)絡(luò)參數(shù)。損失函數(shù)定義為

圖2 ANN的結(jié)構(gòu)

(3)

梯度下降法是一種求解損失函數(shù)最小值的有效方法。假設(shè)L是向量v的函數(shù)，則L的變化可以用v的微小變化來(lái)近似,即

ΔL=L·Δv

(4)

(5)

(6)

式中上標(biāo)k為訓(xùn)練次數(shù)。一般來(lái)說(shuō)，訓(xùn)練ANN的方法有隨機(jī)梯度下降法SGD、均方根傳播算法RMPprop以及自適應(yīng)梯度算法Adam等。本文采用Adam訓(xùn)練ANN。

雖然訓(xùn)練ANN過(guò)程耗時(shí)較長(zhǎng)，但在激活階段能夠根據(jù)輸入快速計(jì)算輸出，不需要任何迭代過(guò)程。需要指出的是，在求解熱傳導(dǎo)反問(wèn)題過(guò)程中，網(wǎng)絡(luò)的輸入為溫度數(shù)據(jù)，輸出則為未知參數(shù)或函數(shù)。

2.2 卷積神經(jīng)網(wǎng)絡(luò)CNN

CNN(Convolutional Neural Network)是ANN的一個(gè)重要分支，主要具有以下特點(diǎn)，同一個(gè)卷積層的所有神經(jīng)元共享一組相同的權(quán)重和偏差；每一層只有部分神經(jīng)元與下一層的神經(jīng)元連接；在卷積層之后引入池化層。CNN的結(jié)構(gòu)如圖3所示?？梢钥闯觯珻NN的結(jié)構(gòu)與傳統(tǒng)的ANN非常相似，不同的是，神經(jīng)元之間的連接是空間稀疏的。

圖3 CNN的結(jié)構(gòu)

在卷積運(yùn)算中，卷積核作為一個(gè)濾波器，應(yīng)用于輸入層的所有空間，以提取每個(gè)空間位置的特征。此外，前一層的特征可以由卷積核獲得，更復(fù)雜的特征由后一層的卷積核得到，然后進(jìn)行池化操作。經(jīng)卷積和池化后的每個(gè)特征都與全連接層相連，這些層之間的連接與ANN相同。

2.3 自動(dòng)編碼器AE

自20世紀(jì)80年代以來(lái)，機(jī)器學(xué)習(xí)領(lǐng)域中AE(Auto Encoder)的概念就一直廣受關(guān)注。AE通常由編碼器和解碼器兩個(gè)部分組成，有如下特點(diǎn)，(1) 有損壓縮和解壓數(shù)據(jù)；(2) 提取數(shù)據(jù)特征；(3) 對(duì)數(shù)據(jù)本身進(jìn)行學(xué)習(xí)。AE的結(jié)構(gòu)如圖4所示。

圖4 AE的結(jié)構(gòu)

AE的基本思想是構(gòu)建一個(gè)與輸入具有相同維度的輸出層。編碼器函數(shù)可以表示為

h=f(x)

(7)

式中h為密碼。將密碼h傳遞給解碼器，可以重構(gòu)得到數(shù)據(jù)r。解碼器函數(shù)可以表示為

r=g(h)=g(f(x))

(8)

給定AE的損失函數(shù)如下，

(9)

式中S為輸出層參數(shù)的個(gè)數(shù)。本文編碼器的作用是提取數(shù)據(jù)特征，解碼器的作用是從特征中復(fù)現(xiàn)樣本數(shù)據(jù)。

3 數(shù)據(jù)驅(qū)動(dòng)求解熱傳導(dǎo)反問(wèn)題

3.1 管道內(nèi)壁識(shí)別問(wèn)題[11]

管道內(nèi)壁識(shí)別問(wèn)題的模型如圖5所示。假定管道外壁的半徑為定值r1，管道內(nèi)壁的半徑r2是極角θ的函數(shù)，導(dǎo)熱系數(shù)為常數(shù)k0。則管道熱傳導(dǎo)問(wèn)題的控制方程表示為

圖5 管道內(nèi)壁識(shí)別幾何模型

(10)

給定邊界條件如下,

T(r,θ)=Tatr=r2

(11)

-k0[?T(r,θ)/?r]=qatr=r1

(12)

在幾何形狀識(shí)別問(wèn)題中，除管道內(nèi)邊界形狀r2=p(θ)未知外，所有參數(shù)和邊界條件均已知。管道內(nèi)壁幾何形狀識(shí)別問(wèn)題的目標(biāo)是通過(guò)測(cè)量表面或區(qū)域內(nèi)的溫度來(lái)反求管道內(nèi)壁的幾何形狀。

本文采用一種轉(zhuǎn)換方法如圖6所示，在求解過(guò)程中，將未知的邊界幾何形狀轉(zhuǎn)換為規(guī)則形狀，但導(dǎo)熱系數(shù)是邊界坐標(biāo)的函數(shù)。然后,將幾何形狀識(shí)別問(wèn)題轉(zhuǎn)化為有效導(dǎo)熱系數(shù)的識(shí)別問(wèn)題，從而避免了求解過(guò)程中網(wǎng)格的重構(gòu)。有效導(dǎo)熱系數(shù)可計(jì)算為

圖6 轉(zhuǎn)換后的管道模型

(13)

然后,將原問(wèn)題轉(zhuǎn)化為求解功能梯度材料的熱傳導(dǎo)問(wèn)題，控制方程為

(14)

式中k為式(13)的有效導(dǎo)熱系數(shù)，同樣給定邊界條件為式(11,12)。

得到有效導(dǎo)熱系數(shù)后，經(jīng)轉(zhuǎn)換，可得未知的邊界幾何形狀為

(15)

管道內(nèi)壁識(shí)別的流程如圖7所示。所有算例都在3.4 GHz CPU(Interl Xeon E5-1680 v4 128.00 GB RAM)的個(gè)人計(jì)算機(jī)上進(jìn)行。

圖7 管道內(nèi)壁識(shí)別的計(jì)算流程

3.1.1 測(cè)量誤差對(duì)計(jì)算結(jié)果的影響

(16)

討論不同隱藏層數(shù)和訓(xùn)練次數(shù)對(duì)計(jì)算結(jié)果的影響，隱藏層分別取2,3和4，訓(xùn)練次數(shù)分別設(shè)置為500次、1000次和2000次，使用隨機(jī)產(chǎn)生的5000個(gè)模型作為訓(xùn)練樣本。計(jì)算結(jié)果的損失函數(shù)值列入表1。

表1 不同隱藏層及訓(xùn)練次數(shù)的損失函數(shù)值

由表1可知，隨著訓(xùn)練次數(shù)的增加，損失函數(shù)值逐漸減小。當(dāng)隱藏層數(shù)為3時(shí)，經(jīng)過(guò)2000次訓(xùn)練得到的損失函數(shù)值最小。所以，在后續(xù)的計(jì)算中設(shè)置隱藏層數(shù)為3，訓(xùn)練的次數(shù)為2000次。

討論不同測(cè)量誤差對(duì)計(jì)算結(jié)果的影響。測(cè)量誤差分別取為1%和5%。表2和圖8給出了不同測(cè)量誤差的反演結(jié)果，圖9展示了反演結(jié)果的相對(duì)誤差。

表2 不同測(cè)量誤差時(shí)參數(shù)的反演結(jié)果

表2表明，在不考慮測(cè)量誤差的情況下，反演結(jié)果與精確解非常接近。隨著測(cè)量誤差的增加，反演結(jié)果偏離精確解。從圖8和圖9可以看出，隨著測(cè)量誤差的增加，反演結(jié)果愈加偏離精確的邊界形狀，相對(duì)誤差變大。在測(cè)量誤差為5%時(shí)，仍然可以得到較為準(zhǔn)確的反演結(jié)果。

圖8 不同測(cè)量誤差的反演結(jié)果

圖9 不同測(cè)量誤差時(shí)，反演結(jié)果的相對(duì)誤差

3.1.2 測(cè)點(diǎn)數(shù)量對(duì)計(jì)算結(jié)果的影響

考慮一個(gè)橢圓管道內(nèi)壁，橢圓的長(zhǎng)軸和短軸分別為 1.5 m 和 1.0 m。真實(shí)的管道內(nèi)邊界幾何形狀可以表示為

式中α1=1.5,α2=1.0,α3=0.25。討論不同測(cè)點(diǎn)數(shù)量對(duì)計(jì)算結(jié)果的影響。如圖10所示，分別采用10個(gè)、20個(gè)和40個(gè)測(cè)點(diǎn)。轉(zhuǎn)化后的有效導(dǎo)熱系數(shù)可以表示為

圖10 不同測(cè)點(diǎn)數(shù)量的有限元模型

(17)

反演結(jié)果列入表3，反演得到的管道內(nèi)邊界形狀如圖11所示。圖12展示了反演結(jié)果的相對(duì)誤差。表3表明，計(jì)算結(jié)果與精確解較為吻合。從 圖11 和圖12可以看出，隨著測(cè)點(diǎn)數(shù)量的增加，反演結(jié)果更接近于準(zhǔn)確的幾何形狀，反演結(jié)果的相對(duì)誤差變小。

表3 不同測(cè)點(diǎn)數(shù)量的反演結(jié)果

圖11 不同測(cè)點(diǎn)數(shù)量的反演結(jié)果

圖12 不同測(cè)點(diǎn)數(shù)量時(shí)反演結(jié)果的相對(duì)誤差

3.2 皮膚腫瘤生長(zhǎng)參數(shù)識(shí)別問(wèn)題

如圖13(a)所示，在笛卡爾坐標(biāo)系下建立三維皮膚模型，考慮了五種主要的皮膚組織，分別為表皮、乳頭狀真皮、網(wǎng)狀真皮、脂肪和肌肉。假設(shè)每一層組織都是有限厚度dm的均勻介質(zhì)。皮膚組織的有限元模型如圖13(b)所示。皮膚組織的生物熱傳導(dǎo)控制方程可以表示為

圖13 皮膚組織模型

(18)

表4 皮膚組織的熱物性參數(shù)[5,12]Tab.4 Properties of the skin tissue[5,12]

皮膚腫瘤生長(zhǎng)參數(shù)識(shí)別的深度學(xué)習(xí)模型如 圖14 所示。首先，將皮膚腫瘤的生長(zhǎng)參數(shù)形成矩陣；然后利用AE的編碼器部分提取矩陣的抽象特征，利用AE的解碼器部分重構(gòu)參數(shù)矩陣。

圖14 腫瘤生長(zhǎng)參數(shù)識(shí)別流程

在此基礎(chǔ)上，利用隨機(jī)模型，對(duì)帶有不同生長(zhǎng)參數(shù)的有限元模型進(jìn)行分析，得到了皮膚表面隨時(shí)間變化的溫度數(shù)據(jù)。同時(shí)，利用數(shù)據(jù)訓(xùn)練CNN，建立皮膚表面溫度與AE特征之間的抽象映射關(guān)系。之后，利用真實(shí)的皮膚表面溫度，使用已訓(xùn)練的CNN獲得生長(zhǎng)參數(shù)的抽象特征。最后，由AE的解碼器部分得到腫瘤的生長(zhǎng)參數(shù)信息。

本文采用18900個(gè)隨機(jī)模型作為樣本，同時(shí)訓(xùn)練AE和CNN，分別識(shí)別了Clark II和Clark IV腫瘤的生長(zhǎng)參數(shù)。同時(shí)，考慮不同測(cè)量誤差對(duì)計(jì)算結(jié)果的影響。表5和表6分別給出了Clark II和Clark IV腫瘤在1%和2%測(cè)量誤差下，深度學(xué)習(xí)模型的識(shí)別結(jié)果。

表5 不同測(cè)量誤差下的Clark II反演結(jié)果

表6 不同測(cè)量誤差下的Clark IV反演結(jié)果

由表5和表6可知，深度學(xué)習(xí)模型能夠準(zhǔn)確估計(jì)腫瘤的生長(zhǎng)參數(shù)。數(shù)值結(jié)果表明，皮膚表面的溫度信息包含了腫瘤生長(zhǎng)參數(shù)的特征。深度學(xué)習(xí)模型可以提取皮膚表面溫度的抽象特征，建立起生長(zhǎng)參數(shù)與特征之間的抽象關(guān)系。隨著測(cè)量誤差的減小，反演結(jié)果更加精確。在2%的測(cè)量誤差時(shí)，仍可以得到較為準(zhǔn)確的反演結(jié)果。值得注意的是，本文提出的方法能夠準(zhǔn)確識(shí)別腫瘤的生熱率和血液灌注率，為腫瘤的醫(yī)療診斷提供參考。

4 結(jié) 論

本文提出了一種基于數(shù)據(jù)驅(qū)動(dòng)模型求解熱傳導(dǎo)反問(wèn)題的方法，實(shí)現(xiàn)了對(duì)導(dǎo)熱管道內(nèi)壁幾何形狀和皮膚腫瘤生長(zhǎng)參數(shù)的識(shí)別。通過(guò)數(shù)值算例，對(duì)提出的兩個(gè)模型進(jìn)行了驗(yàn)證，同時(shí)討論了不同測(cè)量誤差對(duì)反演結(jié)果的影響。計(jì)算結(jié)果表明，基于數(shù)據(jù)驅(qū)動(dòng)模型，能夠有效求解熱傳導(dǎo)反問(wèn)題。隨著測(cè)量誤差的減小，反演結(jié)果更加精確。

此外，本文提出的求解模型可以進(jìn)一步擴(kuò)展，用于求解其他類型的熱傳導(dǎo)反問(wèn)題，如導(dǎo)熱系數(shù)識(shí)別、邊界條件重構(gòu)以及熱源識(shí)別等問(wèn)題。