校正光通量变化的定量光声内窥成像_孟琪.pdf

0111001-1研究论文第 43 卷 第 1 期/2023 年 1 月/光学学报校正光通量变化的定量光声内窥成像孟琪1，孙正1，2*，侯英飒1，孙美晨11华北电力大学电子与通信工程系，河北 保定 071003；2华北电力大学河北省电力物联网技术重点实验室，河北 保定 071003摘要 在光声内窥成像中，不均匀或不稳定的照明，以及生物组织的复杂光学特性均会导致成像平面内光通量分布不均匀，从而造成重建图像质量和成像深度的下降。提出了一种校正光通量变化的定量光声内窥成像方法，对光吸收能量分布进行稀疏分解，采用贪婪算法重构得到光吸收系数和光通量的稀疏表示，并通过稀疏矩阵分解实现光吸收系数与光通量分布的联合重建。仿真和体模实验结果表明，与一步法和基于模型的定量重建方法相比，采用所提方法估算光吸收系数的均方根误差（RMSE）可降低约 48%，重建图像的归一化平均绝对距离（NMSAD）和结构相似度（SSIM）可分别降低约 25%和提高约 24%。与其他校正光通量的重建算法相比，所提方法估计光吸收系数的 RMSE 可降低约 22%、NMSAD可降低约 20%、SSIM 可提高约 10%。关键词 成像系统；图像重建技术；光声成像；内窥成像；光吸收系数；光通量；定量成像中图分类号 TP391 文献标志码 A DOI：10.3788/AOS2212351引 言光声成像（PAI）是一种基于光声效应的新型复合功能成像方法1。该方法的原理是用短脉冲（ns量级）激光照射组织，组织吸收部分光能量后受热膨胀，激发出宽带（MHz量级）超声波并向组织表面传播，超声换能器接收声压信号后送入计算机，通过求解声学逆问题重建出初始声压分布图或光吸收能量分布图，从而显示出目标的形态结构。其中，常用的方法有反投影2、时间反演（TR）3、傅里叶变换法（FTR）4和反卷积法（DR）5等。由于光吸收能量与光吸收系数和局部光通量的乘积成正比，其中光通量不仅与组织的光学特性有关，还与光强分布有关，因此光吸收能量分布图无法准确反映组织的光学特性。通过求解光学逆问题定量估算组织的光学特性参数（如光吸收系数和散射系数）6、热膨胀系数（Grneisen 系数）和功能特性参数（如发色团质量分数和血氧饱和度）等，可实现定量成像，获得组织的功能成分信息7。定量光声成像（qPAI）是一个非线性病态问题，一般采用基于模型的最优化方案求解，即建立前向成像物理模型，并迭代求解使声压测量值或光吸收能量重建值与前向模型输出的理论值之间误差最小的非线性最小二乘问题，在每次迭代中更新待估计参数，直到前向模型的输出与测量数据相匹配7-11。实现 qPAI的主要困难是光通量未知，为了简化问题，常规方法通常基于理想照明假设和简单光子传输模式，在光源位置固定的情况下，假设入射光在组织表面和内部形成均匀光照。该假设对于较小的成像目标具有一定的合理性，但对于具有复杂结构的大型目标而言，不合理的照明模式会造成光照不均匀，只有靠近和面向入射光的部分组织可被相对均匀地照射。此外，大部分生物组织都是高散射且不透明的混浊介质，组织结构与成分具有高度复杂性和差异性，不同组织成分对光的吸收和散射存在明显差异，也会导致光通量分布不均匀。同时，随着入射光在组织中传播深度的增加，光的穿透能力在减弱，光能量趋于分散，在不同位置处的光通量呈现不同的衰减趋势。上述原因会导致图像质量和成像深度的下降，在低光通量区域中会出现目标模糊和精度降低等问题。目前解决非均匀光照问题的方法主要包括优化照明设备、校正光通量变化和估计光通量三种方案12。优化照明设备的目的是保证成像目标内各个视角上具有足够的光穿透性，如旋转照明13、定制多模光纤照明14、光捕捉器15和光纤扩散器16等，均需改变成像系统。在图像重建过程中对光通量变化进行校正可提高成像精度，如：Rosenthal 等17将图像分解为低空间频率全局分量和高空间频率局域分量，再通过少量基函数的和逼近图像；Bu等18将基于模型的光吸收系数重收稿日期：2022-06-01；修回日期：2022-06-17；录用日期：2022-06-29；网络首发日期：2022-07-09基金项目：国家自然科学基金（62071181）通信作者：*0111001-2研究论文第 43 卷 第 1 期/2023 年 1 月/光学学报建与光通量补偿相结合，减少因光通量变化所致的重建误差；Nyknen 等19采用基于模型的方法利用初始声压和表面光测量数据同时估计光吸收系数、散射系数和 Grneisen系数，但无法确定最小化问题是否具有唯一最优解。对于光声光谱成像，Den-Ben等20利用可逆光开关荧光蛋白的波长开关特性校正散射介质内光通量分布的不均匀性，而 Zhou等21验证了三种基于模型的光通量校正方法的有效性，即基于一维扩散近似模型、基于扩散偶极子模型和基于蒙特卡罗（MC）模拟的方法。此外，将 PAI与其他成像模态相结合，可实现对光通量的估计或测量，如光声-声光（AO）测量22-23、光声-超声成像24、光声-被动超声成像25和光声-扩散光层析成像（DOT）26等。双模态成像方法增加了系统及其操作的复杂度，故其应用范围受限。光声内窥成像将无创的 PAI和内窥检测技术相结合，可用于腔体病变（如消化道病变和冠状动脉粥样硬化等）的早期诊断。提出了一种定量光声内窥成像方法，该方法不依赖于前向光传输模型，而是将 Curvelet和 Haar联合基与傅里叶基作为基函数，对光吸收能量分布进行稀疏表示，通过矩阵分解实现腔体横截面上光吸收系数和光通量分布图的联合重建，从而校正光通量变化对定量重建精度的影响。2原理与方法如图 1所示，所提方法基于两步策略：1）采用常规图像重建方法（如反投影或 TR）根据超声探测器测量的光声信号重建腔体横截面上的光吸收能量密度（AOED）分布图；2）建立光吸收能量的稀疏表达式，得到光吸收系数和光通量的稀疏表示，并通过稀疏矩阵分解实现光吸收系数与光通量分布的联合重建。在 PAI中，光吸收能量、光吸收系数和光通量之间的关系为A(r)=a(r)(r)，（1）式中：A(r)、a(r)和(r)分别为成像平面中位置 r处的光吸收能量、光吸收系数和光通量。成像平面极坐标系中光吸收能量的稀疏表达式为lgA(，l)=lga(，l)+lg(，l)=n=1Ncnn(，l)+m=1Mdmm()，l，（2）式中：(，l)为位置 r 的极坐标，如图 2 所示；n(，l)和m(，l)分别为lga(，l)和lg(，l)稀疏表示的基函数；cn和dm为稀疏分解系数，其中n=1，2，N，m=1，2，M。基函数需满足稀疏条件，即每个基函数只能单一地表示lga(，l)或lg(，l)。由于光吸收系数和光通量分别满足局部高频率和全局低频率特性17，因此本文采用离散 Curvelet 和 Haar 联合基27-28表示lga(，l)，采用离散傅里叶基29表示lg(，l)，构成基函数库。此时，光吸收能量的稀疏分解问题就转化成了寻找满足稀疏表达式的、最少数量系数cn和dm的问题。定义优化问题为Y=minkk0，s.t.Xk-Ckk2 k，（3）图 1联合重建光吸收系数和光通量分布图的流程Fig.1Process for simultaneously reconstructing optical absorption coefficient and optical fluence distribution图 2成像平面极坐标系Fig.2Polar coordinate system of imaging plane0111001-3研究论文第 43 卷 第 1 期/2023 年 1 月/光学学报式中：k=1，2，360为成像平面极坐标系中的第 k个角度；k=c1cNd1dMT为由分解系数cn和dm构成的列向量；k为重构误差阈值；0和 2分 别 为 零 范 数 和 2 范 数；Ck=Cn，kCm，k为L（N+M）维的稀疏矩阵，其各列分量分别对应基函数库n和m，L为一个角度上离散采样点的个数，第 L个采样点对应的极径为L l（单位为 mm），l为一个角度上相邻采样点之间的径向距离；Xk为由第 k 个角度上的光吸收能量对数值组成的列向量。式（3）是非凸优化问题，常用的求解方法有凸松弛算法和贪婪算法30：凸松弛算法是将原目标函数中难以解决的部分用松弛变量代替，使目标函数成为凸函数，如基追踪（BP）算法31；贪婪算法是将求解的问题分成若干子问题，求解子问题的局部最优解并进行合成，得到当前条件下原问题整体最优解的近似解。与凸松弛算法相比，贪婪算法的复杂度低、运算量小、重构效率高且相对稳定，故其在解决稀疏分解问题中的应用最为广泛30。贪婪算法的典型代表是匹配追踪（MP）算法32，它采用迭代的方式根据基函数库重构图像，每次迭代时将基函数库中的一个元素添加到图像的稀疏表示中，当重构精度满足要求或超过预先设定的迭代次数时即停止迭代。由于每个库元素的系数只选择一次，因此无法在当前迭代步骤中修正上一次迭代中的误差。MP 算法的改进算法是正交匹配追踪（OMP）算法30，每次迭代都会重新计算所有基函数的最佳系数，可以补偿之前迭代步骤中的误差。同时，OMP算法的搜索过程具有较好的随机性，可有效避免陷入局部最优，但每次迭代只能选择基函数库中的一个元素，故效率较低。OMP算法的改进算法有多候选正交匹配追踪（MOMP）33和基于 Dice 系数的正交匹配追踪（DOMP）算法34等，这些算法不仅保留了 OMP的优点，还具有更高的重构精度和效率。因此，本文采用 DOMP 算法求解式（3）的优化问题，在第 4.1 节中将结合仿真实验结果讨论其相对于其他非凸优化求解算法的优势。DOMP 算法的原理是：在每一次迭代中，根据Dice系数匹配准则求解 Dice值最大的一个元素，并加入到稀疏矩阵中，更新残差判定重构结果，达到精度要求或允许的迭代次数时停止迭代。求解式（3）的具体步骤如下。第一步，初始化迭代次数 i=1，设置初始稀疏矩阵Ck，0和稀疏分解系数向量k，0均为空矩阵，则Xk的初始值Xk，0为Xk，0=Ck，0k，0=，（4）残差Rk，0为Rk，0=Xk-Xk，0=Xk。（5）第二步，计算第 i 次迭代后的稀疏矩阵。根据Dice系数匹配准则在基函数矩阵库中搜索与残差相似程度最大的库函数作为最匹配的库函数，即 i=arg max|Dice()，Rk，i-1i=arg max|Dice()，Rk，i-1，（6）式中：和 为lga(k，i)和lg(k，i)的库函数；i和i为第 i次迭代后搜索出的最匹配库函数；Rk，i-1为第 i-1次迭代的残差；Dice()为 Dice系数匹配准则函数，其表达式为Dice(x，y)=2i=1P()x2i+y2ii=1Pxiyi，（7）式中：xi和 yi分别为向量 x 和向量 y 中的元素；P 为向量中的元素个数。更新稀疏矩阵为Ck，i=Cn，k，i-1iCm，k，i-1i，（8）式中：Cn，k，i-1和Cm，k，i-1分别为Ck，i-1中对应基函数库n和m的矩阵。第三步，计算第 i次迭代后的稀疏分解系数矩阵k，i=C+k，iXk=(CHk，iCk，i)-1CHk，iXk，（9）式中：C+k，i和CHk，i分别为Ck，i的 Moore-Penrose伪逆矩阵和共轭转置矩阵。第四步，计算第 i次迭代后的残差图 3矩阵分解示意图Fig.3Schematic diagram of matrix decomposition0111001-4研究论文第 43 卷 第 1 期/2023 年 1 月/光学学报Rk，i=Xk-Xk，i=Xk-Ck，ik，i。（10）第五步，判断迭代是否终止。如果第 i 次迭代后的 误 差 满 足Rk，i-Rk，i-12 0.01Xk2，则 终 止 迭代，输出稀疏矩阵Ck，i和稀疏分解系数矩阵k，i。否则，更新迭代次数 ii+1，并转向第二步。按照上述步骤求解出稀疏矩阵和稀疏分解系数矩阵，进而得到光吸收系数和光通量对数的向量表示lg a，k=Cn，kn，klg k=