基于COMSOL的乳腺癌光声检测仿真研究

尚栋琳 赵欣洒 韩建宁

摘  要：在COMSOL Multiphysics中建立光声检测的有限元模型，模拟在激光脉冲作用下乳腺组织中光—热—声转换的过程以及该过程中光能量、温度、位移和声压数值的变化。结果表明，乳腺组织在接受脉冲激光辐照时，组织内部肿瘤处的光能量值大约是其周围正常组织的81%，温度与正常组织相比略高，最高约为0.28 K，由瞬态热膨胀产生的超声压力波最高约为5 kPa，中心频率为10 MHz。进而说明所建立的光声转换过程切实有效，有助于区分正常组织和病变肿瘤组织，实现乳腺癌的早期诊断与治疗。

关键词：乳腺癌；光声检测；COMSOL

中图分类号：TP391.4  文献标识码：A  文章编号：2096-4706（2023）15-0039-07

Simulation Study of Photoacoustic Detection of Breast Cancer Based on COMSOL

SHANG Donglin， ZHAO Xinsa， HAN Jianning

（School of Information and Communication Engineering， North University of China， Taiyuan  030051， China）

Abstract： Establish a finite element model for photoacoustic detection in COMSOL Multiphysics to simulate the process of light-heat-sound conversion in breast tissue under the action of laser pulses， as well as the changes in light energy， temperature， displacement， and sound pressure values during this process. The results show that when breast tissue is irradiated by pulsed laser， the light energy value at the tumor site inside the tissue is about 81% of that of the surrounding normal tissue. The temperature is slightly higher than that of normal tissue， with a maximum of about 0.28 K. The ultrasonic pressure wave generated by transient thermal expansion is about 5 kPa， with a central frequency of 10 MHz. This further shows that the established photoacoustic conversion process is effective and helpful to distinguish between normal tissue and diseased tumor tissue， achieving early diagnosis and treatment of breast cancer.

Keywords： breast cancer; photoacoustic detection; COMSOL

0  引  言

我國乳腺癌的发病率增速非常快，是全球增速的2倍，城市地区更为明显[1]，由于发现肿瘤的时期普遍较晚，女性在患乳腺癌之后的死亡率非常高[2]，1990—2019年中国女性乳腺癌标化发病率从17.07万/10万升至35.61万/10万，标化死亡率从7.22万/10万升至13.40万/10万，而早期的检测与诊断能够有效地发现和治疗乳腺癌，大大提升患者的生存率。

当前针对乳腺癌的非入侵式检测手段主要是乳腺超声检查和乳腺钼靶检查[3]，超声检查难以识别小于10 mm的肿块，且探头频率与检查者的操作技术和思维分析对诊断结果影响较大，钼靶检查对高密度乳房中癌症肿块与非癌症区域区分度低；细胞学穿刺检查作为一种入侵式的检测手段，有使癌细胞转移的风险，这些方法或是会产生误诊漏诊，或是会对人体产生危害，作为一种非电离无损伤，同时兼顾高对比度与检测深度的新型成像方式，光声成像（Photo acoustic Imaging， PAI）近些年在全世界备受研究人员的关注。

光声成像以光声效应[4]为理论基础，生物组织被脉冲激光照射后，吸收光子转化为热能引起局部温度上升，产生瞬态热膨胀，最后产生宽带（兆赫兹级）的超声波，超声换能器可以接收到此超声波[5]。传统的超声诊断需要利用造影剂将散射的回声信号加强，使诊断结果的分辨力和特异性得到提升[6]，而在癌症早期，PAI则不需要注射造影剂，而是直接利用人体内不同组份的光吸收能力不同来改变接收到的超声信号强度[7]，血液中的血红蛋白是一种很强的光吸收体，肿瘤在生成、发展和转移的过程中边缘区域新生血管密度增大，血红蛋白比例增加，使得光声成像更加适合成为对乳腺癌的一种检测方式[8]，因此，许多研究小组对乳腺癌的光声成像进行大量研究，Li [9]课题组开发出一套360°全环光声成像系统—人体乳腺光声断层成像系统，Siphanto [10]团队利用光声检测成功将乳腺肿瘤微小血管可视化，中科院Song[11]等学者研发了光声显微镜和光声内窥成像，在成像方面取得了不俗的成绩，这些研究成果说明了在肿瘤检测方面光声成像具有很大的潜力和可行性。

然而，医学成像实验成本巨大，且对于实验人体风险高，为了降低成本与安全风险，缩短实验周期，尽可能最优的调整参数，本文利用有限元法在COMSOL中的数学、生物传热、固体力学和瞬态压力声学模块，模拟了以光能量、温度、位移和声压为量化参数的光声肿瘤检测，实现了光热-热膨胀-超声波的转化过程，最后得到光声效应所产生的光声信号，成功将正常组织与病变肿瘤组织区分开来，为乳腺癌的光声早期检测提供仿真基础。

1  理论分析与模型建立

1.1  光传输

路德维希·玻尔兹曼于1872年提出玻尔兹曼传输方程：

该方程有效描述了在人体内光的传输过程[12]，然而此方程难以解出，对传输方程进行漫射近似可以得到形式简单而又容易解出的漫方程，含时间的漫射方程可以写作：

其中D表示漫射系数，由光吸收系数μa，散射系数μs，各向异性因子g决定：

COMSOL数学模块中有偏微分方程：

对比偏微分方程和漫射近似方程，将其中相应的项进行替换，α、β、γ和ea分别设置为0，da设置为n/c，n表示折射率，c表示真空中的光速，a设置为吸收系数μa，扩散系数c则替换为式（3）中的D，完成相应的替换后即可在COMSOL中描述出光在组织中的传输，将其中的源f设置为高斯函数：

完成点源发射出激光束的设置，Wp表示激光能量4 mJ/cm2，τp表示脉冲宽度5 ns，τc表示峰值时间30 ns，其余部分光学参数如表1所示。

1.2  热传输

生物组织被激光辐照时，会吸收大量光子能量，光子能量会转化为热能，这一过程称之为生物组织的光热效应[13]。接受激光辐照时温度的差异是区分正常乳腺组织和肿瘤的重要标志，在COMSOL中有Pennes“生物传热方程”[14]：

其中ρ表示血液密度，c表示血液比热容，ωb表示血液灌注率，Tart表示血液温度，Qr表示代谢热源，因为激光辐照时间非常短，所以在模型设置中忽略体温调节机制。在COMSOL生物传热模块中对Pennes方程进行调整和设置，完成模型的热传输设置。

1.3  热膨胀

在前两个模块的基础上添加固体力学模块，将生物传热与固体力学两个模块通过求解下列方程进行耦合，以模拟生物组织吸收光子后产生的瞬态热膨胀：

εth = α（T ）（T - Tref）                     （7）

其中，εth表示热应变，α表示热膨胀系数，T表示实际温度，Tref表示参考温度，设置Tref为310.5 K，耦合接口设置为生物传热与固体力学，所需的部分热学与力学部分参数如表2所示。

1.4  光声信号

固体力学处的设置使得生物组织发生热膨胀效应，热弹性膨胀压力不能及时传递出去，最终产生以超声波为载体的光声信号。COMSOL中的瞬态压力声学模块可以很好地模拟出超声波的产生和传播，声场的产生与传播满足以下方程：

其中，c表示声音在人体组织中的传播速度，设置为1 450 m/s，ρ表示组织密度，pt表示背景压力场，qd和Qm分别表示单极和偶極域源，在瞬态压力声学模块中分别将这三个参数设置为0，即可完成对声场的设置。

1.5  模型的几何参数

为了减小计算量，此仿真将乳腺组织设置为半径为10 mm×5 mm的矩形居中放置，在乳腺组织模型内部居中设置了半径为0.3 mm的圆形肿瘤，乳腺组织下方设置10 mm×3 mm的水层，使得激光可以均匀辐照在乳腺组织表面，在水层下方（0，-5.5）处设置点光源，整个模型的几何形状如图1所示。

本文中所建立的模型利用自由三角形进行网格剖分，单元大小设置为常规。构建完成后图形如图2所示。

2  仿真结果及讨论

2.1  光能量分布

组织中光能量的变化如图3所示，截取在四个不同时刻乳腺组织内部光分布图。

16 ns时，激光辐照之初，乳腺组织和肿瘤中光能量较低，并且激光辐照深度增加光能量值减小，肿瘤及其周围的光能量值差距不大，由于模型底部有一层水，使得整个乳腺组织底部能够被激光脉冲均匀辐照，所以该时刻乳腺组织的二维图像中光能量颜色呈现横向条状图形；30 ns时，激光能量达到峰值，组织中光能量明显增加且为全周期最高，此时组织内部光能量值变大，但肿瘤及其周围颜色差距变大；37 ns时，激光能量峰值已过，相比于30 ns，此时整个组织内部光能量值降低，肿瘤的光能量值明显小于周围组织的光能量值，并且肿瘤周围正常组织的光能量分布很不均匀，这是因为肿瘤的光能量吸收系数远大于乳腺组织，散射系数远小于乳腺组织，导致肿瘤可以吸收周围组织大量的光子，对周围的组织产生了影响；41 ns时，随着激光能量降低，乳腺组织中光能量值逐渐下降，此时肿瘤周围正常组织的光能量分布依旧不均匀。组织的光能量分布是光声成像对肿瘤检测的开始，肿瘤和乳腺组织吸收系数和散射系数的差别把二者的光能量分布变得可视化，为癌症检测提供了一个可靠的量化指标。

以肿瘤的中心为中点，选取贯穿肿瘤中心的径向和轴向截线，在此截线上对乳腺组织和肿瘤的光能量进行定量分析，因激光在30 ns时达到顶峰，所以时间选取t = 30 ns，如图4（a）所示，可以得到径向截线上的光能量分布，弧长值最低时距离激光源最近，光能量值最高，随着弧长长度的增加，截线上光能量分布数值呈指数下降，在弧长约为2.2 mm到2.8 mm之间的光能量值发生突变，快速下降，此处为肿瘤的位置，2.8 mm之后上升并逐渐稳定；图4（b）为轴向截线上的光能量分布情况，在弧长约为4.7 mm和5.3 mm处，即肿瘤与正常组织的分界处，光能量值有轻微突变现象，其周围正常组织中的光能量分布不均匀，肿瘤处的光能量分布最低约为周围正常组织光能量分布的81%，这与之前二维光能量图像的结果相符合。

2.2  温度分布

光能量的变化会引起组织中温度的变化，本节研究了乳腺组织和肿瘤随时间变化的温度分布情况，选取时间为0 ns、28 ns、30 ns和35 ns时组织的乳腺温度分布，结果如图5所示，t = 0 ns时，乳腺组织和肿瘤温度均为初始温度310.5 K，28 ns时随着激光辐照，乳腺组织和肿瘤的温度开始增加，直到35 ns时，激光辐照结束，温度达到最高值，肿瘤内部的温度达到310.63 K以上，在接受激光辐照之后肿瘤处的温度相比于正常组织高，这是因为肿瘤吸收了大量光子，光能转换为热能，且由于肿瘤边界处介质发生突变，导致肿瘤自身温度的分布也并不均匀，其内部温度相比于边界较低。

图6中含有径向截线上t = 20 ns到t = 38 ns之间共10条温度曲线，可大致描述在激光脉冲的作用下温度随时间的变化过程，在t = 28 ns之前，由于激光脉冲能量较小，乳腺组织和肿瘤内部温度几乎不变，维持在310.5 K左右，t = 28 ns时，随着激光脉冲能量增加，组织内部温度有了变化，此时肿瘤内部温度比正常组织温度高0.02 K，t = 30 ns时激光脉冲达到峰值，由于光能量的不断增加，组织内部温度继续升高，在t = 36 ns时达到最高值，此时肿瘤内部温度为310.76 K，正常组织温度310.52 K，二者相差0.24 K，弧长为2.2 mm和2.8 mm处温度达到了310.8 K与310.78 K，这两处是肿瘤和乳腺组织的分界点，肿瘤边界与内部由于介质不同而引起温度突变形成凹坑状，且靠近激光源一侧的温度比远离激光源一侧的温度略高，与正常组织温度相差0.28 K，温度分布图符合预期结果。

2.3  肿瘤表面位移

肿瘤在吸收大量光子后温度升高，产生瞬态热膨胀，本节在原有物理场的基础上加入固体力学模块，在COMSOL中模拟出随着温度的升高肿瘤在乳腺组织中的热膨胀现象，选取时间为0 ns、30 ns、50 ns、100 ns的仿真结果如图7所示。

由图7可以看出，在0 ns时，没有激光辐照，肿瘤内部和边界的颜色基本一致，肿瘤边界没有产生位移，30 ns时，组织受到激光辐照，肿瘤表面开始产生位移形变，此时位移量很小，50 ns时，肿瘤表面位移相较于30 ns时比较明显，最大是30 ns时的50倍左右，此时肿瘤内部大部分区域均有位移，由边界到中心位移量逐渐减小，直到肿瘤的圆心处，可以看出此处无位移，100 ns时，肿瘤的位移量和位移区域相较于50 ns时均减小，此时肿瘤处于收缩状态。

图8为径向截线上100 ns时的位移值，可以看出，弧长为2.2 mm与2.8 mm处产生明显位移，这两处为肿瘤边界，2.2 mm处靠近激光源，此处位移最大，越靠近2.5 mm处位移越小，2.5 mm处位移值约为0，符合100 ns时位移二维图像的结果。

为了更好地观察肿瘤表面位移大小随时间的变化，我们选取肿瘤表面上下边界两点进行观察，将其位移分为X分量与Y分量，位移与时间的关系如图9所示。

从图9可以看出，前30 ns由于温度变化量较小，肿瘤边界几乎无位移，30 ns时肿瘤开始产生热膨胀，肿瘤边界开始发生位移，上下两边界点处位移Y分量逐渐增加，随着激光辐照结束，温度到达顶峰，在70 ns时肿瘤的热膨胀到達顶峰，此时位移值最大，随后开始减小，在100 ns时达到最小值，之后略有波动直至稳定，波动的原因是肿瘤与乳腺组织都设置为线弹性材料，在位移回复过程中存在阻尼振动，稳定后的位移值相较于开始也有很大的变化，原因是稳定后的温度相较于初始温度较高，所以肿瘤边界仍会保留一定的位移量。本节对于肿瘤位移模拟的结果与上节温度变化结果相符合，同时是光声信号产生的基础。

2.4  光声信号的产生

肿瘤在经过热膨胀之后产生超声波，即光声信号，本节在前文三个模块的基础上加入瞬态压力声学模块，并将声与固体力学进行耦合，模拟光声信号的产生与传播，如图10为时间在60 ns、120 ns、300 ns、600 ns时组织表面光声信号传播情况。

由图10可以看出，光声信号从肿瘤表面产生，由于乳腺组织材料为各向同性，光声信号在其中传播速度相同，所以随着时间的推移，光声信号以环状的形式向四周扩散；由颜色图例可知，t = 60 ns时，光声信号产生之初，信号强度最大，随着信号距离肿瘤表面越来越远，信号强度逐渐变低，这是因为声波在生物组织中传播存在吸收，散射和衰减的现象。在肿瘤下方（0，-0.5）坐标处放置一个探针，检测此处的光声信号强度，其时域与频谱图如图11所示。

可以看到，该点处光声信号的初始值为0，在大约180 ns时开始出现光声信号并逐渐变大，在250 ns时光声信号达到最大值，大约为5 kPa，250 ns到310 ns之间光声信号值急剧下降，310 ns时到达最低值-3 kPa左右，之后再次上升并由于阻尼振动产生部分波动，最后信号逐渐回归到0，在这个过程中增速越来越慢，从信号的频谱图中可以得到声信号的中心频率为10 MHz，本节得到的时域和频域两方面的声压信号图在其他文献中得以验证[15]。

3  结  论

本文在COMSOL中建立了一个完整的乳腺癌光声可视化模型，在二维图像中观察到超声信号的扩散，成功模拟出从激光脉冲辐照导致光能量发生变化、温度升高、产生瞬态热膨胀、最后产生声压信号的全过程，有利于通过光声信号的参数变化区分正常组织与病变组织，对乳腺癌的早期检测具有一定的参考意义。

参考文献：

[1] FAN L，ZHENG Y，YU K D，et al. Breast cancer in a transitional society over 18 years：trends and present status in Shanghai，China [J].Breast Cancer Research and Treatment，2009，117：409-416.

[2] MCDANIEL S M，RUMER K K，BIROC S L，et al. Remodeling of the Mammary Microenvironment after Lactation Promotes Breast Tumor Cell Metastasis [J].The American journal of Pathology，2006，168（2）：608-620.

[3] 孟淑萍，张正平，王霈，等.CT、超声、X线钼靶在乳腺癌诊断中的应用价值研究[J].中国CT和MRI杂志，2014，12（7）：33-35.

[4] GRINVALD A，LIEKE E，FROSTIG R D，et al. Functional architecture of cortex revealed by optical imaging of intrinsic signals [J].Nature，1986，324（6095）：361-364.

[5] KU G，FORNAGE B D，JIN X，et al. Thermoacoustic and Photoacoustic Tomography of Thick Biological Tissues toward Breast Imaging [J].Technology in Cancer Research &Treatment，2005，4（5）：559-565.

[6] BUSHBERG J T，SEIBERT J A，LEIDHOLDT E M，et al. The Essential Physics of Medical Imaging [M].Lippincott Williams & Wilkins，2011.

[7] 张建英，谢文明，曾志平，等.光声成像技术的最新进展 [J].中国光学，2011，4（2）：111-117.

[8] MALLIDI S，LUKE G P，EMELIANOV S. Photoacoustic imaging in cancer detection，diagnosis，and treatment guidance [J].Trends in Biotechnology，2011，29（5）：213-221.

[9] LI L，WANG L H，MORALES C S. Integration of photoacoustic computed tomography with multitargeted polymer-based nanoparticles visualizes breast cancer intratumor heterogeneity [C]//Photons Plus Ultrasound：Imaging and Sensing.San Francisco：SPIE，2022：1196006. https：//doi.org/10.1117/12.2610584.

[10] SIPHANTO R I，KOLKMAN R G M，HUISJES A，et al. Imaging of small vessels using photoacoustics：an in vivo study [J].Lasers in surgery and medicine，2004，35（5）：354-362.

[11] SONG L，MASLOV K，WANG L H. Multifocal optical-resolution photoacoustic microscopy in vivo [J].Optics letters，2011，36（7）：1236-1238.

[12] FLOCK S T，PATTERSON M S，WILSON B C，et al. Monte Carlo modeling of light propagation in highly scattering tissues. I. Model predictions and comparison with diffusion theory [J].IEEE Transactions on Biomedical Engineering，1989，36（12）：1162-1168.

[13] 谷懷民，杨思华，向良忠.光声成像及其在生物医学中的应用 [J].生物化学与生物物理进展，2006（5）：431-437.

[14] 李艾伦，傅卓佳，李柏纬，等.含肿瘤皮肤组织传热分析的广义有限差分法[J].力学学报，2018，50（5）：1198-1205.

[15] SOWMIYA C，THITTAI A K. Simulation of photoacoustic tomography （PAT） system in COMSOL and comparison of two popular reconstruction techniques [C]//Medical Imaging 2017：Biomedical Applications in Molecular，Structural，and Functional Imaging.Orlando，2017，10137.[2023-02-08].https：//www.spiedigitallibrary.org/conference-proceedings-of-spie/10137/1/Simulation-of-photoacoustic-tomography-PAT-system-in-COMSOLR-and-comparison/10.1117/12.2254450.short.