基于红外检测乳腺肿瘤的多参数反演研究*

黄鉴，卢玫，李博汉，张涛

(上海理工大学，能源与动力工程学院，上海 200093)

1 引 言

乳腺癌为女性中最常见的恶性肿瘤，其发病率自上世纪70年代末起一直呈上升趋势，已成为当前社会的重大公共卫生问题。早期发现、早期诊断是提高疗效、降低乳腺癌死亡率的关键。目前诊断的方法主要有X射线、CT检测、MRI检测等方法，这些方法或者昂贵，或者会对人体有一定的伤害。随着红外热像仪测量精度的不断提高，使得热像仪越来越多的被用到乳腺肿瘤诊断中。

通过红外热像仪检测获得的红外热像图来诊断乳腺肿瘤[1]的方法主要是比对病变乳腺体表温度与正常乳腺体表温度，这种诊断方法仅是一种定性分析。其实，乳腺体表温度不仅与肿瘤的性质有关，也与肿瘤的大小、位置、个体乳腺组织的导热系数等诸多因素有关。国内外学者在红外检测乳腺肿瘤方面已做了一些研究，Gonz′alez[2]研究了不同热像仪测量精度条件下，各种大小肿瘤能被检测到的最大深度。Mital等[3]在二维计算模型下，根据单一的目标函数对乳腺肿瘤的位置、大小及代谢产热进行反演计算，并获得较高精度的解。Bezerra[4]等对乳腺肿瘤的导热系数及内部的血液灌注率进行反演计算，并实验研究了热像仪与人体之间距离对反演精度的影响。研究者在求解此类反问题时用到的优化方法主要有序列二次规划法[4]，遗传算法[5]，神经网络法[6]，模式搜索法[7]等。在求解反问题时，不论采用哪种优化算法，都需利用测点温度信息建立一个目标函数，而单一目标函数在反演多个参数时，对各个参数的优化调整存在一定的盲目性。

本研究以红外热像仪检测乳腺肿瘤为应用背景，通过红外热像仪测得的体表温度，建立与反演参数相关的多目标函数，并对所采用的粒子群优化算法中的更新方式进行了相应的改进，以提高多参数反演的准确度和计算速度。

2 物理模型和数学描述

2.1 物理模型

本研究所讨论的物理问题原型是根据红外检测乳腺肿瘤所得到的体表温度定量分析肿瘤的位置、肿瘤代谢热强度和和乳腺组织的导热系数。这一问题可以抽象为一个含有内热源的多参数导热反问题。首先将检测对象简化为如图1所示的半球体，直径d=0.144 m；肿瘤及肿瘤代谢热抽象为半径r=0.01 m的球体热源；底面Γb为腔体内平面，半球面Γt为红外检测乳腺体表面。人体表面的散热包括体表与周围环境的对流换热、辐射换热和汗液的蒸发，检测时人体静坐，因此忽略汗液蒸发。根据测量所得热像图，选用若干观测点温度，反演热源位置、热源强度及导热系数。

图1 肿瘤检测抽象模型Fig 1 Abstract model for tumor detection

2.2 数学描述

为便于建立数学模型和计算，坐标系设置见图1。半球底面置于坐标系xoy平面，红外热像图中表面温度最高点置于坐标系yoz平面。由于半球体内部物性参数均匀，半球表面热边界条件相同，所以热源中心也在yoz平面，且半球体温度场对称于yoz平面。描述上述物理问题的微分方程采用如式(1)所示的生物体传热方程—PENNES[8]方程，对应的定解条件如式(2)～(4)所示。

λ2T+wbiρbθb(Tb-T)+Qmi=0

(1)

∈Ω=Ω1+Ω2
T=T0∈Γb

(2)

(3)

T=Tm，,jj=1，2，3……n

(4)

式中，T为温度，℃；Tb为血液温度，℃；λ为导热系数，W/m·℃；wbi为血液灌注率，s-1；ρb为血液密度，kg/m3；cb为血液比热容，J/kg·℃；Qmi为组织的代谢产热，W/m3；下标i=1，2，分别对应正常乳腺组织区域Ω1及肿瘤区域Ω2；T0为半球底面温度，℃；Tf为环境温度，℃；h为考虑对流换热和辐射换热的总换热系数，W/m2·℃；Tm,j为边界观测点测量温度，℃；n为测点数。

3 多参数导热反问题的求解

导热反问题求解的一般流程为：首先对反演参数随机预测一个值，根据预测值求解式(1)～(3)所描述的导热正问题，得到求解区域的温度场，然后跟据观测点的计算温度值Tc,j和式(4)给出的测量温度值计算由式(5)定义的目标函数值：

(5)

若目标函数F值满足收敛条件，则反演结束；若不满足，则通过优化算法调整反演参数值，重新代入式(1)～(3)计算，直至目标函数F满足收敛条件。

3.1 辅助目标函数

在多参数反演过程中，仅以式(5)所示单一目标函数作为多个参数调整的依据存在盲目性。不同的物理参数对边界温度存在不同的影响，鉴于这一特点，考虑在式(5)所示目标函数基础上，增加与被反演参数相应的辅助目标函数，用于多参数反问题求解过程中参数的调整、优化。针对本研究的导热反问题，则分别研究热源位置、导热系数等参数对表面温度影响的规律，然后在此基础上建立相应的热源位置目标函数和导热系数目标函数。

(1) 热源位置目标函数

李博汉[9]等在反演热源位置的导热反问题中采用了相关度作为目标函数，其表达式为：

(6)

在同时反演多个参数的情况下，如果相关度受热源强度、导热系数影响不大，则仍可沿用此相关度。为此计算了热源位置、导热系数和代谢热强度等反演参数分别偏离真实值时的相关度，结果见表1。可以看出，相关度受肿瘤位置影响较大：预测位置与真实位置偏差越大，相关度就越小；而相关度受导热系数及热源强度的影响较小。因此，式(6)定义的相关度也适合应用在多参数导热反问题中，作为搜索热源位置的辅助目标函数。

表1 各参数对相关度的影响

(2)导热系数目标函数

导热系数反应了物体热传导的能力，物体置于同样的环境中，较强的热传导能力使内部的热量容易向外传递，边界散热的热流密度增大，边界温度增加，反之边界温度将降低。为分析导热系数与热源强度、热源位置对边界温度的相对影响，计算了表2所示工况的温度场，并将边界上观测点温度值绘制成图2所示曲线。

表2 热源参数及导热系数取值

图2不同导热系数对应边界测点温度示意图

Fig2Testpointstemperaturefordifferentthermalconductivity

图2曲线1、2、3分别对应表2计算工况1、2、3。可以看出，随着导热系数增大，边界上各观测点温度值均相应增大。图2曲线3、4、5分别对应表2计算工况3、4、5，即导热系数不变、热源位置和热源强度变化的3个计算工况。由曲线3、4、5可以看出，热源强度、热源位置的变化对靠近热源位置的边界温度影响较大，而在距离热源位置较远的边界区域，如测点10～20，三个工况对应的温度曲线基本重合，这说明该区域温度受热源位置、热源强度影响很小。因此，当计算得到的导热系数趋于真值时，该区域观测点的计算温度应趋于测量温度值。根据这一特点可以建立式(7)所示的众数为反演导热系数的辅助目标函数：

M=mode(Tm,1-Tc,1，Tm,2-Tc,2，…Tm,n-Tc,n)

(7)

所有边界测点温度计算值与测量值之差的绝对值中出现频率最高的数即为众数，众数越接近0，则说明反演参数中的导热系数越接近真值。

3.2 粒子群算法更新方式改进

粒子群算法PSO(particle swarm optimization，PSO)[10]是Kennedy和Eberhart提出的一种基于群体智能的全局随机搜索算法。粒子群算法中，每个粒子在求解空间中的位置代表了优化过程中的一组预测值，并以各自的速度在求解空间内搜索，其基本位置和速度更新方程为:

xi(k+1)=xi(k)+vi(k+1)

(8)

vi(k+1)=ωvi(k)+c1r1[xpbesti(k)-xi(k)]+c2r2[xgbest(k)-xi(k)]

(9)

式中，k为迭代步计数；xi为粒子i在求解空间中位置矢量，x=(x1,x2,…xN)，N为反演参数数目，即求解空间维数；i=1,2,3,…,s，s为粒子总数；vi为粒子i飞行速度矢量，v=(v1,v2,…vN)；ω为惯性权重；c1为认知系数；c2为社会系数；r1、r2为[0,1]均匀分布的随机数；xpbesti为粒子i所经过的所有位置中目标函数F值最小的位置；xgbest为全部粒子所经过的所有位置中目标函数F值最小的位置。

由式(8)～(9)所描述的更新方程可知，仅有一个目标函数时，粒子搜索下一个空间位置的三个飞行速度分量，由目标函数F控制。而本研究提出了与反演参数相关的多目标函数后，则可以根据辅助目标函数对飞行速度分量进行调整。因此，需要对上述更新方程做出相应的改进(改进粒子群算法—modified particle swarm optimization，MPSO)。

针对热源位置和导热系数反演参数，计算过程中分别记录下所有粒子经过的相似度最高的位置xgpbest1、xgpbest2和众数最小的位置xgpbest3，这些位置中分别具有相对最优的热源位置和导热系数。飞行速度分量更新方程为：

vij(k+1)=ωvij(k)+c1r1[xpbestij(k)-xij(k)]+c2r2[xgbest j(k)-xij(k)]+c3r3[xgpbest j(k)-xij(k)]

(10)

式(10)中，j=(1,2,3)；c3为多目标社会系数；r3为均匀分布随机数。与式(9)相比，式(10)对速度的调整增加了由辅助目标函数确定的第四项。

对于热源位置参数，定义xgpbestj，j=1,2的更新方程为：

(11)

对于导热系数，定义xgpbestj,j=3的更新方程为：

(12)

4 数值计算结果及分析

计算时，设半球底面温度Tc=37℃，半球体表面总对流换热系数h=5W/(m2·℃)，环境温度Tf=22℃。 PENNES方程中各参数设置见表3。

表3 PENNES方程物性参数

表4粒子群算法系数取值及参数搜索范围

Table4ValuesofcoefficientsinMPSO,PSOandestimatedparameters’searchsope

粒子数s10最大迭代次数kmax300惯性权重ω0.8-kkmax×0.8 认知系数c12.0社会系数c2kkmax×0.2 多目标社会系数c30.2-kkmax×0.2 热源强度搜索范围700≤Qm≤70000导热系数搜索范围0.46≤λ≤0.500.22≤y2+z2≤0.58热源位置搜索范围θ=arctan(zy)根据热像图上温度最高测点及相邻两个测点确定

分别采用MPSO和PSO两种方法对不同工况进行反演计算，其计算结果列于表5。

表5所列结果显示，三个反演参数中导热系数及热源位置的相对误差较小, 均在1%以内，而热源强度的相对误差较大。图3所示为各工况中反演得到的热源强度的误差绝对值，图4为对应工况的测点温度分布曲线。对照图3、图4可以发现，体表温度变化越大，计算结果误差就越小。相同病变情况下，MPSO反演得到的结果，其相对误差普遍小于PSO，而且在体表温度变化微弱的情况下，更能显示出MPSO的优势。通过四组不同算式的计算，MPSO达到收敛条件的迭代步数约为100，相较PSO约为150的平均迭代步数，反演所需时间缩短33%。

图3 Qm误差绝对值

工况参数真实值(Qm, λ, y, z)反演结果(Qm, λ, y, z)相对误差εrel/%迭代步数1MPSO(29000 , 0.480 ,0.039 , 0.039)(28944 , 0.480 , 0.03901 , 0.03900)(0.19 , 0.0 , -0.03 , 0.0)72PSO(28976 , 0.480 , 0.03889 , 0.03889)(0.08 , 0.0 , 0.28 , 0.28)2342MPSO(20000 , 0.470 ,0.025 , 0.0433)(19529 , 0.470 , 0.02505 , 0.04335)(2.36 , 0.0 , -0.2 , -0.12)87PSO(19455 , 0.470 , 0.02505 , 0.04336)(2.73 , 0.0 , -0.2 , -0.14)1003MPSO(45000 , 0.0475 ,0.010 , 0.040)(44188 , 0.475 , 0.01000 , 0.04015)(1.8 , 0.0 , 0.0 , -0.38)113PSO(43253 , 0.476 , 0.01008 , 0.04021)(3.88 , -0.21 , -0.8 , -0.53)1014MPSO(1800 , 0.485 ,0.020 , 0.050)(1634 , 0.485 , 0.02002 , 0.05002)(9.22 , 0.0 , -0.1 , -0.04)154PSO(700 , 0.486 , 0.02011 , 0.05025)(61.1 , -0.21 , -0.55 , -0.5)194

图5所示为分别采用PSO和MPSO求解算例4时各反演参数时最优解(xgbest)的收敛曲线，在图5(a)、5(b)中一个较为明显的特征是， PSO搜索得到最优热源位置收敛曲线变化缓和，下一个迭代步获得的最优解近似于上一步。而MPSO有较多的锯齿型波动，下一步获得的最优解与上一步相比更似一个突变。这是由于PSO在位置更新时仅受单一目标函数控制，而MPSO在搜索热源位置的过程中，位置更新受到两个目标函数的影响，并且这两个目标函数所确定的热源位置最优解通常是不一致的，这势必会增大粒子的搜索范围，使其能不拘泥于原来的最优解，因此，MPSO中粒子有着更强的搜索能力，并最终得到更接近真值(y=0.02,z=0.05)的解；图5(c)显示的最优导热系数收敛曲线所反映的情况与热源位置收敛曲线近似，MPSO在搜索导热系数时同样受辅助目标函数影响，最终的反演结果也是MPSO优于PSO。而从图5(d)曲线中可以看出PSO在搜索热源强度时陷入了搜索范围的下界限，并在之后的搜索过程中一直困于局部最优解。MPSO搜索结果也一度到达下界限，但随后便跳出局部最优，并最终获得与真值接近的解，这说明辅助目标函数虽没有直接影响粒子在搜索热源强度时的更新方式，但间接的对搜索过程带来了积极的影响。图6为目标函数F收敛曲线，其结果表明PSO在经过30次左右的迭代步后，目标函数值下降非常缓慢，说明求解已经陷入了局部最优。而MPSO在30次迭代步之后的搜索过程中，目标函数值依然不断下降。

图4各工况测点温度分布图

Fig4Thedistributionoftestpointstemperatureinfourcases

图5各反演参数收敛曲线
(a)热源位置y坐标收敛曲线； (b)热源位置z坐标收敛曲线；
(c)导热系数λ收敛曲线； (d)热源强度Qm收敛曲线；
Fig5Convergenceplotofmultivariables

(a) convergence plot of y; (b) convergence plot of z; (c) convergence plot of λ; (d) convergence plot of Qm

图6 目标函数F收敛曲线Fig 6 Convergence plot of fitness fuction F

粒子群算法是一种基于概率的算法，反演求解过程中具有随机性，引入辅助目标函数，增大了粒子在位置更新过程中遇到全局最优解的概率。相似度目标函数及众数目标函数所确定的最优热源位置，及最优导热系数以xgpbest项的形式参与到粒子的更新过程中，使得粒子群反演参数中的热源位置项及导热系数项往xgpbest靠拢，从而增大得到全局最优解的概率。

5 结论

本研究采用粒子群算法，根据红外热像图上观测点的温度，可准确反演出内部肿瘤未知参数。同时根据热源位置、导热系数对边界温度的不同影响，建立相关度及众数作辅助目标函数的多目标粒子群算法，通过对不同的工况的数值模拟结果表明：与PSO相比，MPSO计算收敛速度快，获得的解更精确；在搜索过程中，最优解更新频次高，不易陷入局部最优；在观测点温度差异小的工况下，MPSO的优势更为明显，提高了反演的准确性。