首次提出了基于多島遺傳算法進行生物組織三維溫度場無損重構的重要思想, 把復雜的生物熱傳導反問題的求解轉換為正問題的求解過程, 并通過實驗對該重構思想的可行性和可靠性進行了驗證。以生物組織內點熱源的位置P(x, y, z)和溫度t為優化變量, 把同一表面各個對應點的實驗溫度值和計算溫度值相減并取絕對值之和, 以此為目標函數逐次迭代。目標值越小, 則當前變量, 即熱源位置和溫度值最優。多島遺傳算法可以很好地應用于生物組織三維溫度場的重構, 無需提取生物組織全部的表面溫度數據, 為熱傳導反問題的研究提供了一條具有借鑒性的方法與思路。
引用本文： 葉福麗, 李凱揚, 史貴連. 基于多島遺傳算法的生物組織三維溫度場重構研究. 生物醫學工程學雜志, 2016, 33(4): 666-673. doi: 10.7507/1001-5515.20160110 復制
0 引言
由生物體體表溫度推導體內的溫度場分布及熱源分布是一個典型的熱傳導反問題[1]。熱傳導反問題是反問題的重要組成部分，求解這類反問題的方程組往往是高度病態的，也就是說熱傳導反問題往往是一類高度不適定問題[2]。不適定問題并不意味著它們是不可解的，只是傳統的解析方法很難甚至無法直接應用于不適定性問題的求解或分析，即便硬性地采用了解析方法，所求得的解或許是毫無實際意義的，因為這些解與問題的真實解之間可能存在很大的差距[3]。為使所求的解逼近實際情況，需要對不適定性問題進行正則化處理[4]。目前，在生物熱傳導反問題的建模仿真中，通常采用的方法有邊界元法、有限元法和有限差分法等，這些方法要么對仿真對象的幾何外形要求十分苛刻，要么計算量很大，尤其是在分析求解之前，必須獲取全部表面的溫度數據，而對于生物體某些部位或器官而言，獲取其全部的表面溫度數據往往十分困難[5]。盡管如此，這些數值方法在生物熱傳導反問題領域有著重要的指導和借鑒價值。
生物傳熱領域的一個積極的研究成果是利用計算機來進行人體三維溫度場的實時重構和模擬，這項技術的發展和進步必將很大程度上提升高溫熱療的效果，為廣大的癌癥患者帶來低成本而高效的治療途徑[6]。對于溫度場反問題的求解，目前通常采取解析和數值求解兩種方法。解析方法的優勢在于系統的行為和特征可以通過方程及其解的形式來解釋和揭示，主要缺點是只有簡單問題才能得到極其有限的解，比如線性或常參數系統等[7]；數值方法可以解決復雜的問題，但這類方法通常面臨繁瑣的邊界條件處理問題，計算量也十分龐大，這也是這類方法應用于臨床的主要難題。
不同的熱源(包括位置、形狀和能量)對應著不同的溫度場和表面溫度分布，它們之間存在著唯一映射關系[8]。基于該理論，本文提出了基于多島遺傳算法(multi-island genetic algorithm, MIGA)進行生物體三維溫度場無損重構的思想，并對該重構思想的可行性和可靠性進行了驗證。由物體表面溫度推導求解內部的溫度場分布和熱源分布是一個典型的熱傳導反問題，本文以模型內熱源的溫度和位置為優化變量，將這一反問題轉換成正問題的求解過程，從而避免了傳統數值方法中繁瑣的邊界條件處理，使得求解溫度場分布和熱源分布的過程得到了很大的簡化。本文的主旨是一種生物組織三維溫度場重構思想的提出及其可行性驗證，因此本文在實驗中采用與人體脂肪組織熱物性參數相近的高純度聚丙烯材料作為熱傳導分析對象，避開了其它因素對驗證結果的影響。
1 基于多島遺傳算法的生物體三維溫度場重構思想
傳統的遺傳算法是單種群的，很容易讓優化過早地陷入成熟。作為一種偽并行算法，多島遺傳算法不僅可以避免早熟，而且極大地提升了收斂速度，同時在優化域中還可以更高效地尋找全局最優解[9]。多島遺傳算法增加了樣本的多樣性, 而計算量僅有少許增加。在該算法中，種群的數量并沒有增加，僅僅是把種群分解成為幾個島, 在各個島上分別進行傳統的遺傳計算并以一定概率在島間進行基因交叉, 其主要意圖是增加算法的多峰搜索能力，避免陷入局部最優[10]。
1.1 多島遺傳算法的基本理論與操作步驟
多島遺傳算法首先將整個進化群體分解為若干個子群體(也可稱之為“島”)，然后進行相應的遺傳操作，它是在傳統遺傳算法基礎上衍生而來的。對于分解后所得到的每個島上的子群體，分別進行獨立的標準遺傳操作，包括選擇、交叉和變異等等。多島遺傳算法的示意如圖 1所示。算法在選擇機制上一般仍然采用輪盤賭的方法，采用優秀個體保留策略，將上一代中的最佳個體復制到下一代，即子代中。為了保證群體的多樣性，進而克制算法的早熟現象，多島遺傳算法會選擇某些個體進行“島間遷移”的操作，也就是將個體從一個島嶼遷移到另一個島嶼上，這種遷移是定期隨機的，而島間遷移系數則決定了需隔代遷移的代數。由此可見，相對于傳統遺傳算法，多島遺傳算法具有更優良的優化效率和全局搜索能力[11]。
圖1
多島遺傳算法中的子種群(即“島”生成原理)
Figure1.
Creating principle of island in MIGA
多島遺傳算法通過反復地循環使用遺傳算法的算子和選擇原則，逐代繁殖并進化，從而使種群更加適應環境。一個多島遺傳算法的操作步驟如下[12]：①初始化群體；②分解成多島；③計算群體上每個個體的適應度值；④以個體適應度值所決定的某種規則，來選擇將進入下一代的個體；⑤按概率Pc進行交叉操作；⑥按概率Pm進行變異操作；⑦按概率Pmig進行島間遷移；⑧如果沒有滿足某種停止條件，那么轉入第②步，否則進入第⑨步；⑨找出目標值最優的個體并輸出，以之作為問題的最優解。
多島遺傳算法是一種偽并行的優化算法，因為個體存在一定概率的島間遷移，增加了基因復制的多樣性，所以與傳統遺傳算法比較起來，有著計算效率高、避免過早收斂等優勢。
1.2 基于多島遺傳算法的生物體三維溫度場重構思想
設定生物體內點熱源的位置P(x, y, z)和溫度t為優化變量，把同一表面各個對應點的實驗溫度值和計算溫度值相減并取絕對值之和，以此為目標函數逐次迭代。目標值越小，則當前變量，即熱源位置和溫度值最優。目標函數如式(1)所示：
| $ {f_{\min }}\left( {x,y,z,t} \right) = \sum {|{T_n} - {T_{En}}|} $ |
式中，Tn和TEn分別為同一表面對應點的仿真溫度和實驗溫度，n為所取的對應點的個數。當點熱源的最優位置和溫度值確定下來之后，生物體三維溫度場的分布信息也隨之確定下來。由物體表面溫度推導求解內部的溫度場分布和熱源分布是一個典型的熱傳導反問題，本文提出的優化模型將這一反問題轉換成正問題的求解過程，從而避免了傳統數值方法中繁瑣的邊界條件處理過程。如果以表面溫度作為邊界條件，傳統的數值方法則必須首先得到物體全部表面的溫度數據，對于生物體而言，這個過程是十分繁瑣的甚至是無法獲得的，而本文提出的優化模型只需要有限的表面溫度數據即可。
事實上，在建立優化模型時，本文提出了一種新穎的目標函數的建立方法，即在實驗模型的紅外熱圖中，選取同一表面的兩條垂直相交直線上的若干個等距點(點的數量根據模型的實際大小自由合理設置)，提取這些點的溫度值，然后減去每次迭代仿真的相同表面對應點的溫度值，并取絕對值之和。這種方法簡單易行，避免了繁瑣的目標函數的建立過程，極大地提高了智能優化速度，在分析生物體熱源分布和溫度場分布等方面，具有很大的理論和應用價值。所以式(1)可修正為：
| $ {f_{\min }}\left( {x,y,z,t} \right) = \sum {|{T_1n} - {T_{E1n}}|} + \sum {|{T_{2n}} - {T_{E2n}}|} $ |
式中，T1n和T2n分別表示位于模型某表面的兩條直線上的點的仿真溫度值，TE1n和TE2n則分別表示實驗模型同一表面對應點的實驗溫度值。
1.3 集成優化平臺ISIGHT與大型有限元分析軟件ANSYS
ANSYS是一款大型有限元軟件，ANSYS熱分析涵蓋了三種典型的熱傳遞方式，能夠進行各種復雜環境下的熱模擬、仿真和分析。然而，僅靠ANSYS本身的GUI功能難以滿足目前各種各樣以及日益復雜的工業需求，為此產生了ANSYS參數化設計語言APDL(ANSYS Parametric Design Language)。作為一種有力的方法和工具，它為廣大用戶提供和展現了一個自動實施有限元分析的平臺，指定的函數或變量可以設定為程序的輸入。應用APDL語言，不僅可以通過參數化變量的方式建立新的分析模型，還可以完成ANSYS中常規的分析操作，根據選用的分析類型，用戶可以自由地做出選擇。此外，APDL是ANSYS實現優化設計和自適應網格劃分的一種不可或缺的工具[13]。
ISIGHT是一個集合多學科設計優化的框架。它把設計探索技術、推理技術和數字技術有機結合并進行了綜合性開發，其功能得到了前所未有的拓展，就好像一個軟件機器人在替代工程設計人員進行重復性的數字處理和設計工作，也就是由軟件來自動化處理大量本需要人工來完成的工作。ISIGHT軟件可以集成多學科仿真代碼并提供相應的智能支持，實現對多個設計備選方案進行研究和評估，從而極大地減少了產品設計的周期，工作效率也得到明顯的提升。
2 實驗及數據處理
2.1 實驗材料與過程
采用與人體脂肪組織密度和導熱率近似的高純度聚丙烯材料(密度：920 kg/m3；導熱率：0.236 W/(m·K)[14])材料進行熱傳導實驗，其尺寸為8.5 cm×6.2 cm×6.6 cm(坐標設置如圖 2所示)，電源采用5 V穩壓電源，內置熱源分別采用15.5 Ω和30.5 Ω的電阻，細導線4根。在長方體聚丙烯模塊的表面選擇兩點，位置為x=3.3 cm, y=6.2 cm, z=3.6 cm和x=5.8 cm, y=6.2 cm, z=3.2 cm，分別打兩個深度為2.0 cm和2.3 cm的孔，為熱源的放置做好準備。模型表面的測溫工具為非制冷式焦平面紅外熱像儀(溫度分辨率0.05 ℃，空間分辨率1.4 mrad，像素640×480，聚焦范圍0.5 m～∞)。熱源的溫度探測采用JM222型高精度數字點溫計，其探針可以穿過孔內的聚丙烯粉末到達電阻中間處的孔內壁(植入電阻后，以聚丙烯粉末填實)。電阻表面溫度也可以用點溫計探測，但由于有少部分熱量由孔壁與電阻之間的聚丙烯粉末流出，所以具有分析價值的是處于電阻中間位置的孔內壁溫度。
圖2
聚丙烯模塊及其坐標設置
Figure2.
Polypropylene module and its setting of coordinate
在其中一個孔內插入15.5 Ω的電阻熱源，另一個孔內插入30.5 Ω的電阻熱源，此時兩個電阻的中心所在位置分別為x=3.3 cm, y=4.2 cm, z=3.6 cm和x=5.8 cm, y=3.9 cm, z=3.2 cm，填入聚丙烯粉末并壓實，盡量避免電阻所產生的熱通過孔的縫隙傳遞出來，從而影響熱分析與仿真結果。通過四根導線將2個電阻與5 V的恒壓源連接。同樣，為了保證聚丙烯的6個面是自由面，實驗中采用4根木牙簽對聚丙烯模塊進行了支撐。開通電源后，電阻給聚丙烯加熱，待熱能從聚丙烯內部傳導至表面并達到熱平衡后(此過程持續40 min)，采用紅外探測器采集聚丙烯表面溫度。當得到滿意的熱圖后，存儲并進行溫度數據的提取，為優化算法中的目標函數做好準備。打開點溫計，通過孔內的聚丙烯粉末把點溫計的探針插入到電阻的中心位置，依次獲得兩個熱源的溫度值。根據流體力學原理，當氣流速度僅為V≤0.15 m/s，環境溫度為18~20 ℃時，空氣中的物體表面熱傳導歸屬自然層流對流。實驗過程中，保持環境溫度為19.4 ℃，保持門窗關門，并盡量減少實驗人員的走動，在此情況下，表面對流傳熱系數取3.368 W/(m2·K)。
2.2 數據處理
在該聚丙烯模型內，通過改變2個熱源的溫度值(熱源的位置不變)，我們進行了5次實驗并相應地做了仿真，本文給出了其中一例的數據采集、處理方法及其相應的仿真過程。采用紅外探測器獲取的聚丙烯表面(z=6.6 cm)，并通過探測器專用溫度數據提取軟件，提取z=6.6 cm表面上的分布在兩條直線上28個點的溫度(如圖 3所示)，這兩條直線分別為:①y=4.2 cm, z=6.6 cm和②x=4.25 cm, z=6.6 cm。在兩條直線上分別等距取14個點，28個點坐標及其對應的溫度值如表 1所示，為目標函數的列寫做好準備。這里選取的兩條直線也盡可能地通過了紅外熱圖的高溫區，使其盡可能反映聚丙烯表面溫度的分布特征。通過點溫計探測到相應于兩個電阻中間位置的孔內壁的溫度分別為84.2 ℃(15.5 Ω)和60.9 ℃(30.5 Ω)。
圖3
聚丙烯表面z=6.6 cm上所取的兩條直線及其等距點(2個熱源)
Figure3.
Two straight lines acquired in polypropylene module surface z=6.6 cm and their equally spaced points (2 heat sources)
表1
所取28個點的坐標及其對應的溫度值(2個熱源)
Table1.
Coordinate value and the temperature value of 28 equally spaces points (2 heat sources)
3 仿真與結果
3.1 仿真過程
基于ISIGHT多島遺傳算法進行溫度場重構的仿真流程如圖 4所示。
圖4
基于多島遺傳算法溫度場重構的ISIGHT流程圖(2個熱源)
Figure4.
ISIGHT flow diagram of reconstruction of temperature field based on MIGA (2 heat sources)
圖 4中，三個數據處理模塊之間的數據傳遞是單向交互的。開始時模塊“MIGA Optimization”將初始化變量輸出給模塊“ANSYS 2-heat point”，模塊“ANSYS 2-heat point”將接收到的變量代入模型進行有限元計算，接著將計算結果傳遞給模塊“Objective Calculator”，模塊“Objective Calculator”存儲有自編的適應度函數計算程序，可以自動計算每個個體的適應度值，之后將適應度值傳遞給優化模塊“MIGA Optimization”，從而對每個個體進行適應度評價，模塊“MIGA Optimization”根據優化設置(如迭代次數、適應度達到某個最小值等)決定是否更新變量并傳遞給下一個模塊。三個模塊的變量(數據)交換關系具體如圖 5所示。
圖5
三個模塊的變量(數據)交換關系圖(2個熱源)
Figure5.
Commutation relation schema of variables (data) of three modules (2 heat sources)
在“Objective Calculator模塊”中，我們需要設置適應度函數并進行計算, 結合實驗中(2個熱源)采集到的兩條直線上的28個等距點的溫度數據，在“Objective Calculator模塊”中設置的適應度函數具體如下：
| $ \begin{gathered} {f_{\min }} = |{T_{11}} - 25.32| + |{T_{12}} - 26.28| + |{T_{13}} - 27.48| + |{T_{14}} - 28.21| + |{T_{15}} - 28.74| + |{T_{16}} - 30.27| + \\ |{T_{17}} - 30.05| + |{T_{18}} - 30.22| + |{T_{19}} - 29.69| + |{T_{110}} - 28.43| + |{T_{111}} - 27.23| + |{T_{112}} - 26.21| + \\ |{T_{113}} - 25.17| + |{T_{114}} - 24.61| + |{T_{21}} - 25.78| + |{T_{22}} - 26.33| + |{T_{23}} - 27.41| + |{T_{24}} - 27.96| + \\ |{T_{25}} - 28.01| + |{T_{26}} - 28.86| + |{T_{27}} - 29.74| + |{T_{28}} - 29.91| + |{T_{29}} - 29.93| + |{T_{210}} - 30.24| + \\ |{T_{211}} - 30.15| + |{T_{212}} - 29.97| + |{T_{213}} - 29.46| + |{T_{214}} - 28.96| \\ \end{gathered} $ |
式(3)中，T1n和T2n(1≤n≤14)為28個點的溫度計算值，來自“ANSYS 2-heat point模塊”，而適應度值的計算結果fmin將傳遞給“MIGA Optimization模塊”。
為了避免ISIGHT在調用ANSYS的過程中反復出現重復性操作，首先需要編寫批處理文件并調入到模塊“ANSYS 2-heat point”中，批處理文件如下：“C:\Program Files\ANSYS Inc\v130\ansys\bin\intel\ansys130”-b-i MIGA bioheat2. txt-o output.txt。
執行多島遺傳算法所采用的參數配置如表 2所示。
表2
多島遺傳算法采用的參數配置
Table2.
Parameter configuration adopted in MIGA
3.2 仿真結果
經過250代10 000次迭代后，2個熱源的仿真位置分別為(3.12 cm, 4.06 cm, 3.55 cm)和(5.54 cm, 3.82 cm, 3.03 cm)，熱源溫度分別為98.94 ℃和66.72 ℃。而實驗中，所設置的熱源位置分別為(3.3 cm, 4.2 cm, 3.6 cm)和(5.8 cm, 3.9 cm, 3.2 cm)，相應于電阻中間位置的孔內壁的溫度為84.2 ℃和60.9 ℃。從以上結果可知，熱源的仿真位置與實驗位置十分接近，但仿真溫度與實驗測量溫度相差較大，原因有如下兩點：①實驗中，所測量的熱源溫度其實是相應于電阻中心位置的孔內壁溫度，這個測量點其實與電阻的中心有一定距離；②孔內植入電阻后，盡管以聚丙烯粉末填入并壓實，但還是有一部分熱量從粉末縫隙間流出。基于以上兩點，并結合聚丙烯材料的導熱性能，熱源的仿真溫度與實驗溫度之間的誤差是完全合理的。優化過程中，我們設置目標函數的方式是以聚丙烯實體z=6.6 cm表面28個等距點的溫度值減去仿真模型z=6.6 cm表面相應28個等距點的溫度值，并取絕對值之和，其值(目標值)越小，結果越優。
當優化仿真的熱源位置和溫度值確定下來之后，模型內部的三維溫度場也隨之確定下來。圖 6顯示了任意截面的溫度分布。這里需要說明的是，由于模型內部的溫度要高于其表面溫度，在所取的截面溫度分布圖中，高出表面最高溫度的區域會顯示成灰色，因此我們需要在ANSYS中，重新設置色帶顯示的溫度范圍，即重新調整ANSYS中Min contour value和Max coutour value的值，從而避免在截面圖中出現不含溫度意義的灰色區域，使圖 6中色帶所表示的溫度與與截面溫度分布一致。
圖6
模型內部截面溫度分布圖(2個熱源)
Figure6.
Cross section of temperature distribution in internal module (2 heat sources)
3.3 討論與分析
盡管我們設置了總的迭代次數為10 000次，而事實上當迭代次數達到3 766時，解開始收斂。通過多島遺傳算法優化計算出來的熱源位置及其對應溫度值與實驗設置的熱源位置及其測量溫度值之間存在一定誤差。分析其產生的原因有以下幾點：①實驗中，我們分別采用15.5 Ω和30.5 Ω的電阻來模擬兩個點熱源，事實上，這兩個電阻的體積相對聚丙烯模型來說是比較大的；而仿真過程中，我們設定熱源為理想點熱源。②實驗過程中，為了設置熱源，我們首先在聚丙烯的同一表面上打孔，置入電阻后，采用聚丙烯粉末填入孔與導線之間的縫隙并壓緊，但事實上，我們發現還是有少部分熱量從填壓的聚丙烯粉末的微小縫隙間流出。此外，孔內壁與電阻中心還是有一定的距離，因此我們仿真的熱源中心溫度比實際測量的相應于電阻中間位置的孔內壁溫度要高出許多，基于仿真聚丙烯材料的導熱性，這個差距是合理的。③連接電阻的導線也會產生熱，這就相當于置于聚丙烯內的部分導線成為了一個近似的線熱源。與電阻所產生的熱量相比，這部分熱所占的比例很小，但也是誤差產生的一個重要因素。④聚丙烯模型中，熱源實際位置的測量存在一定誤差。孔的直徑大于電阻直徑，電阻在孔內所處的不同姿態勢必會影響其位置的測量結果。⑤盡管我們已經盡量使用了高純度的聚丙烯材料，但事實上，其內部結構也不是完全各向同性的；而在仿真過程中，我們設置了聚丙烯模塊是理想的各向同性。⑥多島遺傳算法本身也屬于一種數值算法，解的空間本身是離散的，這是數值方法與解析方法之間的一種重要差別。在同一個聚丙烯模型內，保持2個熱源的位置不變，通過改變熱源的溫度值，我們進行了5次實驗并相應地做了仿真，誤差分析結果表明，2個熱源的仿真位置和實際位置距離均不超過0.4 cm，而仿真溫度和實際溫度之間的誤差較大。
目前，我們的研究主要集中在基于導熱機制的生物組織溫度場實時仿真，后續的研究將進一步接近人體的真實結構，如在實驗模型中構建肌肉、脂肪、皮膚等結構，引入代謝熱、血液灌注率以及熱傳導的各向異性等熱物性參數。
4 結論
本文提出了一種基于多島遺傳算法進行生物體三維溫度場重構的思想，將復雜的生物熱傳導反問題轉化為正問題的求解過程。在優化仿真過程中，我們又提出了一種新穎的目標函數建立方法，使得優化速度和效率得到了很大的提高。采用與生物脂肪組織的密度和導熱率十分近似的聚丙烯材料作為實驗和仿真對象，進行了2個熱源的熱傳導實驗，采集了相關點的溫度數據，為多島遺傳算法中目標函數的列寫做好準備。仿真過程中，我們采用了自己編寫的APDL程序進行了建模、網格劃分、熱源設置及有限元計算，并基于ISIGHT優化集成平臺對模型的熱源位置和溫度進行了優化。優化結果表明，熱源的仿真位置與實際位置比較接近，但仿真溫度與實際溫度之間的誤差較大。考慮到實驗過程中的環境因素和誤差因素，以及仿真過程中參數設置的理想因素和仿真溫度與實際溫度之間存在較大的差別，這樣的誤差是在合理范圍之內的。當熱源的最優位置及其溫度確定下來以后，其對應的三維溫度場也隨之確定下來。本文提出的生物體三維溫度場重構的思想和方法方便、快捷，蘊藏著廣闊的應用前景，為熱傳導反問題的研究提供了一條具有借鑒性的方法與思路。
0 引言
由生物體體表溫度推導體內的溫度場分布及熱源分布是一個典型的熱傳導反問題[1]。熱傳導反問題是反問題的重要組成部分，求解這類反問題的方程組往往是高度病態的，也就是說熱傳導反問題往往是一類高度不適定問題[2]。不適定問題并不意味著它們是不可解的，只是傳統的解析方法很難甚至無法直接應用于不適定性問題的求解或分析，即便硬性地采用了解析方法，所求得的解或許是毫無實際意義的，因為這些解與問題的真實解之間可能存在很大的差距[3]。為使所求的解逼近實際情況，需要對不適定性問題進行正則化處理[4]。目前，在生物熱傳導反問題的建模仿真中，通常采用的方法有邊界元法、有限元法和有限差分法等，這些方法要么對仿真對象的幾何外形要求十分苛刻，要么計算量很大，尤其是在分析求解之前，必須獲取全部表面的溫度數據，而對于生物體某些部位或器官而言，獲取其全部的表面溫度數據往往十分困難[5]。盡管如此，這些數值方法在生物熱傳導反問題領域有著重要的指導和借鑒價值。
生物傳熱領域的一個積極的研究成果是利用計算機來進行人體三維溫度場的實時重構和模擬，這項技術的發展和進步必將很大程度上提升高溫熱療的效果，為廣大的癌癥患者帶來低成本而高效的治療途徑[6]。對于溫度場反問題的求解，目前通常采取解析和數值求解兩種方法。解析方法的優勢在于系統的行為和特征可以通過方程及其解的形式來解釋和揭示，主要缺點是只有簡單問題才能得到極其有限的解，比如線性或常參數系統等[7]；數值方法可以解決復雜的問題，但這類方法通常面臨繁瑣的邊界條件處理問題，計算量也十分龐大，這也是這類方法應用于臨床的主要難題。
不同的熱源(包括位置、形狀和能量)對應著不同的溫度場和表面溫度分布，它們之間存在著唯一映射關系[8]。基于該理論，本文提出了基于多島遺傳算法(multi-island genetic algorithm, MIGA)進行生物體三維溫度場無損重構的思想，并對該重構思想的可行性和可靠性進行了驗證。由物體表面溫度推導求解內部的溫度場分布和熱源分布是一個典型的熱傳導反問題，本文以模型內熱源的溫度和位置為優化變量，將這一反問題轉換成正問題的求解過程，從而避免了傳統數值方法中繁瑣的邊界條件處理，使得求解溫度場分布和熱源分布的過程得到了很大的簡化。本文的主旨是一種生物組織三維溫度場重構思想的提出及其可行性驗證，因此本文在實驗中采用與人體脂肪組織熱物性參數相近的高純度聚丙烯材料作為熱傳導分析對象，避開了其它因素對驗證結果的影響。
1 基于多島遺傳算法的生物體三維溫度場重構思想
傳統的遺傳算法是單種群的，很容易讓優化過早地陷入成熟。作為一種偽并行算法，多島遺傳算法不僅可以避免早熟，而且極大地提升了收斂速度，同時在優化域中還可以更高效地尋找全局最優解[9]。多島遺傳算法增加了樣本的多樣性, 而計算量僅有少許增加。在該算法中，種群的數量并沒有增加，僅僅是把種群分解成為幾個島, 在各個島上分別進行傳統的遺傳計算并以一定概率在島間進行基因交叉, 其主要意圖是增加算法的多峰搜索能力，避免陷入局部最優[10]。
1.1 多島遺傳算法的基本理論與操作步驟
多島遺傳算法首先將整個進化群體分解為若干個子群體(也可稱之為“島”)，然后進行相應的遺傳操作，它是在傳統遺傳算法基礎上衍生而來的。對于分解后所得到的每個島上的子群體，分別進行獨立的標準遺傳操作，包括選擇、交叉和變異等等。多島遺傳算法的示意如圖 1所示。算法在選擇機制上一般仍然采用輪盤賭的方法，采用優秀個體保留策略，將上一代中的最佳個體復制到下一代，即子代中。為了保證群體的多樣性，進而克制算法的早熟現象，多島遺傳算法會選擇某些個體進行“島間遷移”的操作，也就是將個體從一個島嶼遷移到另一個島嶼上，這種遷移是定期隨機的，而島間遷移系數則決定了需隔代遷移的代數。由此可見，相對于傳統遺傳算法，多島遺傳算法具有更優良的優化效率和全局搜索能力[11]。
圖1
多島遺傳算法中的子種群(即“島”生成原理)
Figure1.
Creating principle of island in MIGA
多島遺傳算法通過反復地循環使用遺傳算法的算子和選擇原則，逐代繁殖并進化，從而使種群更加適應環境。一個多島遺傳算法的操作步驟如下[12]：①初始化群體；②分解成多島；③計算群體上每個個體的適應度值；④以個體適應度值所決定的某種規則，來選擇將進入下一代的個體；⑤按概率Pc進行交叉操作；⑥按概率Pm進行變異操作；⑦按概率Pmig進行島間遷移；⑧如果沒有滿足某種停止條件，那么轉入第②步，否則進入第⑨步；⑨找出目標值最優的個體并輸出，以之作為問題的最優解。
多島遺傳算法是一種偽并行的優化算法，因為個體存在一定概率的島間遷移，增加了基因復制的多樣性，所以與傳統遺傳算法比較起來，有著計算效率高、避免過早收斂等優勢。
1.2 基于多島遺傳算法的生物體三維溫度場重構思想
設定生物體內點熱源的位置P(x, y, z)和溫度t為優化變量，把同一表面各個對應點的實驗溫度值和計算溫度值相減并取絕對值之和，以此為目標函數逐次迭代。目標值越小，則當前變量，即熱源位置和溫度值最優。目標函數如式(1)所示：
| $ {f_{\min }}\left( {x,y,z,t} \right) = \sum {|{T_n} - {T_{En}}|} $ |
式中，Tn和TEn分別為同一表面對應點的仿真溫度和實驗溫度，n為所取的對應點的個數。當點熱源的最優位置和溫度值確定下來之后，生物體三維溫度場的分布信息也隨之確定下來。由物體表面溫度推導求解內部的溫度場分布和熱源分布是一個典型的熱傳導反問題，本文提出的優化模型將這一反問題轉換成正問題的求解過程，從而避免了傳統數值方法中繁瑣的邊界條件處理過程。如果以表面溫度作為邊界條件，傳統的數值方法則必須首先得到物體全部表面的溫度數據，對于生物體而言，這個過程是十分繁瑣的甚至是無法獲得的，而本文提出的優化模型只需要有限的表面溫度數據即可。
事實上，在建立優化模型時，本文提出了一種新穎的目標函數的建立方法，即在實驗模型的紅外熱圖中，選取同一表面的兩條垂直相交直線上的若干個等距點(點的數量根據模型的實際大小自由合理設置)，提取這些點的溫度值，然后減去每次迭代仿真的相同表面對應點的溫度值，并取絕對值之和。這種方法簡單易行，避免了繁瑣的目標函數的建立過程，極大地提高了智能優化速度，在分析生物體熱源分布和溫度場分布等方面，具有很大的理論和應用價值。所以式(1)可修正為：
| $ {f_{\min }}\left( {x,y,z,t} \right) = \sum {|{T_1n} - {T_{E1n}}|} + \sum {|{T_{2n}} - {T_{E2n}}|} $ |
式中，T1n和T2n分別表示位于模型某表面的兩條直線上的點的仿真溫度值，TE1n和TE2n則分別表示實驗模型同一表面對應點的實驗溫度值。
1.3 集成優化平臺ISIGHT與大型有限元分析軟件ANSYS
ANSYS是一款大型有限元軟件，ANSYS熱分析涵蓋了三種典型的熱傳遞方式，能夠進行各種復雜環境下的熱模擬、仿真和分析。然而，僅靠ANSYS本身的GUI功能難以滿足目前各種各樣以及日益復雜的工業需求，為此產生了ANSYS參數化設計語言APDL(ANSYS Parametric Design Language)。作為一種有力的方法和工具，它為廣大用戶提供和展現了一個自動實施有限元分析的平臺，指定的函數或變量可以設定為程序的輸入。應用APDL語言，不僅可以通過參數化變量的方式建立新的分析模型，還可以完成ANSYS中常規的分析操作，根據選用的分析類型，用戶可以自由地做出選擇。此外，APDL是ANSYS實現優化設計和自適應網格劃分的一種不可或缺的工具[13]。
ISIGHT是一個集合多學科設計優化的框架。它把設計探索技術、推理技術和數字技術有機結合并進行了綜合性開發，其功能得到了前所未有的拓展，就好像一個軟件機器人在替代工程設計人員進行重復性的數字處理和設計工作，也就是由軟件來自動化處理大量本需要人工來完成的工作。ISIGHT軟件可以集成多學科仿真代碼并提供相應的智能支持，實現對多個設計備選方案進行研究和評估，從而極大地減少了產品設計的周期，工作效率也得到明顯的提升。
2 實驗及數據處理
2.1 實驗材料與過程
采用與人體脂肪組織密度和導熱率近似的高純度聚丙烯材料(密度：920 kg/m3；導熱率：0.236 W/(m·K)[14])材料進行熱傳導實驗，其尺寸為8.5 cm×6.2 cm×6.6 cm(坐標設置如圖 2所示)，電源采用5 V穩壓電源，內置熱源分別采用15.5 Ω和30.5 Ω的電阻，細導線4根。在長方體聚丙烯模塊的表面選擇兩點，位置為x=3.3 cm, y=6.2 cm, z=3.6 cm和x=5.8 cm, y=6.2 cm, z=3.2 cm，分別打兩個深度為2.0 cm和2.3 cm的孔，為熱源的放置做好準備。模型表面的測溫工具為非制冷式焦平面紅外熱像儀(溫度分辨率0.05 ℃，空間分辨率1.4 mrad，像素640×480，聚焦范圍0.5 m～∞)。熱源的溫度探測采用JM222型高精度數字點溫計，其探針可以穿過孔內的聚丙烯粉末到達電阻中間處的孔內壁(植入電阻后，以聚丙烯粉末填實)。電阻表面溫度也可以用點溫計探測，但由于有少部分熱量由孔壁與電阻之間的聚丙烯粉末流出，所以具有分析價值的是處于電阻中間位置的孔內壁溫度。
圖2
聚丙烯模塊及其坐標設置
Figure2.
Polypropylene module and its setting of coordinate
在其中一個孔內插入15.5 Ω的電阻熱源，另一個孔內插入30.5 Ω的電阻熱源，此時兩個電阻的中心所在位置分別為x=3.3 cm, y=4.2 cm, z=3.6 cm和x=5.8 cm, y=3.9 cm, z=3.2 cm，填入聚丙烯粉末并壓實，盡量避免電阻所產生的熱通過孔的縫隙傳遞出來，從而影響熱分析與仿真結果。通過四根導線將2個電阻與5 V的恒壓源連接。同樣，為了保證聚丙烯的6個面是自由面，實驗中采用4根木牙簽對聚丙烯模塊進行了支撐。開通電源后，電阻給聚丙烯加熱，待熱能從聚丙烯內部傳導至表面并達到熱平衡后(此過程持續40 min)，采用紅外探測器采集聚丙烯表面溫度。當得到滿意的熱圖后，存儲并進行溫度數據的提取，為優化算法中的目標函數做好準備。打開點溫計，通過孔內的聚丙烯粉末把點溫計的探針插入到電阻的中心位置，依次獲得兩個熱源的溫度值。根據流體力學原理，當氣流速度僅為V≤0.15 m/s，環境溫度為18~20 ℃時，空氣中的物體表面熱傳導歸屬自然層流對流。實驗過程中，保持環境溫度為19.4 ℃，保持門窗關門，并盡量減少實驗人員的走動，在此情況下，表面對流傳熱系數取3.368 W/(m2·K)。
2.2 數據處理
在該聚丙烯模型內，通過改變2個熱源的溫度值(熱源的位置不變)，我們進行了5次實驗并相應地做了仿真，本文給出了其中一例的數據采集、處理方法及其相應的仿真過程。采用紅外探測器獲取的聚丙烯表面(z=6.6 cm)，并通過探測器專用溫度數據提取軟件，提取z=6.6 cm表面上的分布在兩條直線上28個點的溫度(如圖 3所示)，這兩條直線分別為:①y=4.2 cm, z=6.6 cm和②x=4.25 cm, z=6.6 cm。在兩條直線上分別等距取14個點，28個點坐標及其對應的溫度值如表 1所示，為目標函數的列寫做好準備。這里選取的兩條直線也盡可能地通過了紅外熱圖的高溫區，使其盡可能反映聚丙烯表面溫度的分布特征。通過點溫計探測到相應于兩個電阻中間位置的孔內壁的溫度分別為84.2 ℃(15.5 Ω)和60.9 ℃(30.5 Ω)。
圖3
聚丙烯表面z=6.6 cm上所取的兩條直線及其等距點(2個熱源)
Figure3.
Two straight lines acquired in polypropylene module surface z=6.6 cm and their equally spaced points (2 heat sources)
表1
所取28個點的坐標及其對應的溫度值(2個熱源)
Table1.
Coordinate value and the temperature value of 28 equally spaces points (2 heat sources)
3 仿真與結果
3.1 仿真過程
基于ISIGHT多島遺傳算法進行溫度場重構的仿真流程如圖 4所示。
圖4
基于多島遺傳算法溫度場重構的ISIGHT流程圖(2個熱源)
Figure4.
ISIGHT flow diagram of reconstruction of temperature field based on MIGA (2 heat sources)
圖 4中，三個數據處理模塊之間的數據傳遞是單向交互的。開始時模塊“MIGA Optimization”將初始化變量輸出給模塊“ANSYS 2-heat point”，模塊“ANSYS 2-heat point”將接收到的變量代入模型進行有限元計算，接著將計算結果傳遞給模塊“Objective Calculator”，模塊“Objective Calculator”存儲有自編的適應度函數計算程序，可以自動計算每個個體的適應度值，之后將適應度值傳遞給優化模塊“MIGA Optimization”，從而對每個個體進行適應度評價，模塊“MIGA Optimization”根據優化設置(如迭代次數、適應度達到某個最小值等)決定是否更新變量并傳遞給下一個模塊。三個模塊的變量(數據)交換關系具體如圖 5所示。
圖5
三個模塊的變量(數據)交換關系圖(2個熱源)
Figure5.
Commutation relation schema of variables (data) of three modules (2 heat sources)
在“Objective Calculator模塊”中，我們需要設置適應度函數并進行計算, 結合實驗中(2個熱源)采集到的兩條直線上的28個等距點的溫度數據，在“Objective Calculator模塊”中設置的適應度函數具體如下：
| $ \begin{gathered} {f_{\min }} = |{T_{11}} - 25.32| + |{T_{12}} - 26.28| + |{T_{13}} - 27.48| + |{T_{14}} - 28.21| + |{T_{15}} - 28.74| + |{T_{16}} - 30.27| + \\ |{T_{17}} - 30.05| + |{T_{18}} - 30.22| + |{T_{19}} - 29.69| + |{T_{110}} - 28.43| + |{T_{111}} - 27.23| + |{T_{112}} - 26.21| + \\ |{T_{113}} - 25.17| + |{T_{114}} - 24.61| + |{T_{21}} - 25.78| + |{T_{22}} - 26.33| + |{T_{23}} - 27.41| + |{T_{24}} - 27.96| + \\ |{T_{25}} - 28.01| + |{T_{26}} - 28.86| + |{T_{27}} - 29.74| + |{T_{28}} - 29.91| + |{T_{29}} - 29.93| + |{T_{210}} - 30.24| + \\ |{T_{211}} - 30.15| + |{T_{212}} - 29.97| + |{T_{213}} - 29.46| + |{T_{214}} - 28.96| \\ \end{gathered} $ |
式(3)中，T1n和T2n(1≤n≤14)為28個點的溫度計算值，來自“ANSYS 2-heat point模塊”，而適應度值的計算結果fmin將傳遞給“MIGA Optimization模塊”。
為了避免ISIGHT在調用ANSYS的過程中反復出現重復性操作，首先需要編寫批處理文件并調入到模塊“ANSYS 2-heat point”中，批處理文件如下：“C:\Program Files\ANSYS Inc\v130\ansys\bin\intel\ansys130”-b-i MIGA bioheat2. txt-o output.txt。
執行多島遺傳算法所采用的參數配置如表 2所示。
表2
多島遺傳算法采用的參數配置
Table2.
Parameter configuration adopted in MIGA
3.2 仿真結果
經過250代10 000次迭代后，2個熱源的仿真位置分別為(3.12 cm, 4.06 cm, 3.55 cm)和(5.54 cm, 3.82 cm, 3.03 cm)，熱源溫度分別為98.94 ℃和66.72 ℃。而實驗中，所設置的熱源位置分別為(3.3 cm, 4.2 cm, 3.6 cm)和(5.8 cm, 3.9 cm, 3.2 cm)，相應于電阻中間位置的孔內壁的溫度為84.2 ℃和60.9 ℃。從以上結果可知，熱源的仿真位置與實驗位置十分接近，但仿真溫度與實驗測量溫度相差較大，原因有如下兩點：①實驗中，所測量的熱源溫度其實是相應于電阻中心位置的孔內壁溫度，這個測量點其實與電阻的中心有一定距離；②孔內植入電阻后，盡管以聚丙烯粉末填入并壓實，但還是有一部分熱量從粉末縫隙間流出。基于以上兩點，并結合聚丙烯材料的導熱性能，熱源的仿真溫度與實驗溫度之間的誤差是完全合理的。優化過程中，我們設置目標函數的方式是以聚丙烯實體z=6.6 cm表面28個等距點的溫度值減去仿真模型z=6.6 cm表面相應28個等距點的溫度值，并取絕對值之和，其值(目標值)越小，結果越優。
當優化仿真的熱源位置和溫度值確定下來之后，模型內部的三維溫度場也隨之確定下來。圖 6顯示了任意截面的溫度分布。這里需要說明的是，由于模型內部的溫度要高于其表面溫度，在所取的截面溫度分布圖中，高出表面最高溫度的區域會顯示成灰色，因此我們需要在ANSYS中，重新設置色帶顯示的溫度范圍，即重新調整ANSYS中Min contour value和Max coutour value的值，從而避免在截面圖中出現不含溫度意義的灰色區域，使圖 6中色帶所表示的溫度與與截面溫度分布一致。
圖6
模型內部截面溫度分布圖(2個熱源)
Figure6.
Cross section of temperature distribution in internal module (2 heat sources)
3.3 討論與分析
盡管我們設置了總的迭代次數為10 000次，而事實上當迭代次數達到3 766時，解開始收斂。通過多島遺傳算法優化計算出來的熱源位置及其對應溫度值與實驗設置的熱源位置及其測量溫度值之間存在一定誤差。分析其產生的原因有以下幾點：①實驗中，我們分別采用15.5 Ω和30.5 Ω的電阻來模擬兩個點熱源，事實上，這兩個電阻的體積相對聚丙烯模型來說是比較大的；而仿真過程中，我們設定熱源為理想點熱源。②實驗過程中，為了設置熱源，我們首先在聚丙烯的同一表面上打孔，置入電阻后，采用聚丙烯粉末填入孔與導線之間的縫隙并壓緊，但事實上，我們發現還是有少部分熱量從填壓的聚丙烯粉末的微小縫隙間流出。此外，孔內壁與電阻中心還是有一定的距離，因此我們仿真的熱源中心溫度比實際測量的相應于電阻中間位置的孔內壁溫度要高出許多，基于仿真聚丙烯材料的導熱性，這個差距是合理的。③連接電阻的導線也會產生熱，這就相當于置于聚丙烯內的部分導線成為了一個近似的線熱源。與電阻所產生的熱量相比，這部分熱所占的比例很小，但也是誤差產生的一個重要因素。④聚丙烯模型中，熱源實際位置的測量存在一定誤差。孔的直徑大于電阻直徑，電阻在孔內所處的不同姿態勢必會影響其位置的測量結果。⑤盡管我們已經盡量使用了高純度的聚丙烯材料，但事實上，其內部結構也不是完全各向同性的；而在仿真過程中，我們設置了聚丙烯模塊是理想的各向同性。⑥多島遺傳算法本身也屬于一種數值算法，解的空間本身是離散的，這是數值方法與解析方法之間的一種重要差別。在同一個聚丙烯模型內，保持2個熱源的位置不變，通過改變熱源的溫度值，我們進行了5次實驗并相應地做了仿真，誤差分析結果表明，2個熱源的仿真位置和實際位置距離均不超過0.4 cm，而仿真溫度和實際溫度之間的誤差較大。
目前，我們的研究主要集中在基于導熱機制的生物組織溫度場實時仿真，后續的研究將進一步接近人體的真實結構，如在實驗模型中構建肌肉、脂肪、皮膚等結構，引入代謝熱、血液灌注率以及熱傳導的各向異性等熱物性參數。
4 結論
本文提出了一種基于多島遺傳算法進行生物體三維溫度場重構的思想，將復雜的生物熱傳導反問題轉化為正問題的求解過程。在優化仿真過程中，我們又提出了一種新穎的目標函數建立方法，使得優化速度和效率得到了很大的提高。采用與生物脂肪組織的密度和導熱率十分近似的聚丙烯材料作為實驗和仿真對象，進行了2個熱源的熱傳導實驗，采集了相關點的溫度數據，為多島遺傳算法中目標函數的列寫做好準備。仿真過程中，我們采用了自己編寫的APDL程序進行了建模、網格劃分、熱源設置及有限元計算，并基于ISIGHT優化集成平臺對模型的熱源位置和溫度進行了優化。優化結果表明，熱源的仿真位置與實際位置比較接近，但仿真溫度與實際溫度之間的誤差較大。考慮到實驗過程中的環境因素和誤差因素，以及仿真過程中參數設置的理想因素和仿真溫度與實際溫度之間存在較大的差別，這樣的誤差是在合理范圍之內的。當熱源的最優位置及其溫度確定下來以后，其對應的三維溫度場也隨之確定下來。本文提出的生物體三維溫度場重構的思想和方法方便、快捷，蘊藏著廣闊的應用前景，為熱傳導反問題的研究提供了一條具有借鑒性的方法與思路。
