内窥式近红外频域漫射层析成像图像重构算法

秦转萍，赵会娟，杨彦双

(天津大学精密仪器与光电子工程学院，天津 300072)

漫射层析成像(diffuse optical tomography，DOT)是一种无创功能成像技术[1-4]．目前近红外漫射光成像的研究大都应用于一些外部器官，如乳房、头部等[3-4]．采用漫射光进行内窥成像研究开展的较晚，但其研究对于早期宫颈癌和前列腺癌的诊断具有重要的意义．该方法除了具有设备低价、可携带以及对操作者水平依赖性低等优点外，更重要的是具有无创性，可通过观察组织体的新陈代谢状况而实现在体早期癌症诊断．国际上进行内窥式漫射光测量及诊断应用研究较早的当属美国德克萨斯大学的研究组，他们采用近红外连续光照射、连续光检测的方法，但结果表明对原位宫颈癌诊断率较低，对癌前组织几乎没有辨别能力[5]．2006年由美国 MediSpectra Inc公司生产了漫射光宫颈影像系统[6]，该系统虽然可对宫颈上皮内中、高度瘤变(CIN 2、3)进行诊断，但图像空间分辨率低，对癌前病变等级分化不清晰，且仅能对外宫颈进行诊断，无法诊断宫颈管内病变．

近红外漫射光测量系统可分为只提供光强测量量的连续光系统以及可提供多个测量量的时间分辨系统和频域系统．相对于时间分辨系统，频域系统具有测量速度快、系统相对简单等优点．频域系统采用高频信号调制光源，通过测量光经过组织体后光子密度波的幅度衰减和相位延迟实现光学参数重构．2006年 Piao等[7]针对前列腺、直肠等组织首次进行了内窥式近红外频域漫射光成像，并取得了较好的动物实验结果．但受检测设备的限制，其频域系统只获得了直流(DC)数据，用 NIRFAST软件包进行的图像重构结果表明，目标定位精度受目标病灶深度的影响，定位不准确[7-8]．

笔者研究了针对早期宫颈癌诊断的频域漫射光层析成像图像重构方法．首先介绍了所发展的图像重构理论；其次采用所发展的图像重构算法，对吸收系数和约化散射系数分别进行重构，以验证算法的正确性；最后，为与文献[8]仅利用直流数据进行重构的结果进行比较，在无限和有限2种吸收系数对比度情况下，对吸收系数进行了重构，以验证本文所提算法可消除目标深度对目标定位精度的影响．

1 图像重构理论

考虑到实际测量数据的有限性，研究了同时利用频域测量值中幅值和相位信息的重构方法，用于增加测量数据的信息量，减少逆问题的欠定性和病态性．

1.1 正问题

光子在组织中的传播模型采用漫射方程的频域形式，边界条件采用考虑内反射效应的罗宾边界条件[9]．以Φ、0ω、c、0q分别表示光子密度、调制频率、光在组织中的传播速度、迷向光源，Ω表示成像区域，Ω∂表示Ω的边界，s表示一个源点位置，ξ表示一个探测点位置，D为剖分总节点数．

采用有限元法求解频域漫射方程，最终得矩阵方程为

式中：K、B为D×D维矩阵；Φ、Q、P为D×1维矩阵；l、j为节点的整体编号；l, j = 1,2,…,D；μl为第l个节点的形状函数；Φ (ξ , s ,x)为光源在s处激励时在边界ξ处探测到的光子密度；a = ( 1 + Rf) /(1 - Rf)，Rf为扩散传输内反射系数；

通过式(1)可求得光子密度Φ，并通过罗宾边界条件求出探测点处的光子流

频域测量量的输出采用Rytov形式，分别用实部和虚部表示幅值和相位，即Γreal=real(ln(Γ (ξ ,s, x ) ))，Γimag=imag(ln(Γ (ξ ,s, x ) ))．

1.2 图像重构算法

由于组织体的强散射效应、光在组织中的非直线传播、测量数据的有限性、组织体的非均匀性等因素的影响，近红外漫射光图像重构是一个非线性欠定、病态逆问题．采用高斯牛顿法求解该逆问题．正问题算子采用 ()Fx表示．设在宫颈内边界上采用 S个源点照射、M 个点测量，考虑到在每个测量点可获得幅度和相位 2个量，总测量数据量为. 则目标函数为

式中 ΓC和Γ分别为测量的和由正问题计算得到的边界光子流．

对目标函数Ψ ( x)求导，令导数Ψ ′(x ) =0．用泰勒级数对Ψ′(x)在xk附近展开，保留线性项，并采用阻尼最小二乘法(Levenberg-Marquarat，LM)进行正则化．则最终得

式中：δxk+1为迭代更新因子；I为单位阵；λk为第k次迭代的正则化参数；J为雅可比矩阵，J =[Jµa,Jµs′]；Jµa、Jµs′分别为关于μa和μs′的雅可比矩阵．

当解空间维数非常大时，J的计算非常费时．本研究对J的求解采用伴随源法，此方法仅用2次正向求解便可实现对 J的构建，可简化计算，大大减少计算时间．但是 Jµs′的构建涉及到对光子密度求梯度，这不仅耗时，而且在源点附近显示出某些数值上的奇异性，从而导致算法对于组织体深层散射信息不敏感．基于此情况，采用了修正广义脉冲谱技术，此技术利用了强散射媒质 μs′ ＞＞μa的性质、散度定理及频域漫射方程，对进行了近似求解[2,10-11]307-329．同时 J采用 Rytov近似，则最终可获得Jµa、Jµs′的表达式为

式中：r′为区域Ω内任意一点的位置；Ψ*(r ′ ,ξ,x)为源放在探测点ξ处进行激励、由伴随理论求得的r′处的光子密度的共轭值．

在频域内J的结构为

当 J的维数非常大时，不仅求解方程(4)非常费时，而且存储Hessian矩阵(JΤJ)会占用非常大的内存，这一问题对于三维成像尤为重要．为了解决此问题，采用广义最小余量(generalized minimal residual，GMRES)Krylov方法进行求解[1]．

如第k次迭代后满足中止条件，则用 xk来近似真实值 xtrue．

2 光学参数重构精度的验证

2.1 用于图像重构的组织模型

考虑到子宫颈的尺寸和形状，成像区域为圆环模形，其中圆环内半径为 10,mm，外半径为 20,mm，源点和探测点在内环边界上．由于内环尺寸较小，在漫射方程适用范围内[12]107-111，最终选用16个源点和16个探测点，如图 1所示．每次仅在一个源点激励，所有探测点同时获得测量数据，然后换作另一个源点激励并测量，直到所有源点完成激励．

图1 成像区域及源和探测点分布Fig.1 Imaging domain and the distribution of 16 sources and 16 detectors

2.2 吸收系数和约化散射系数的重构精度

为了评价重构出的图像的质量，分别采用保真度和半宽度来评价重构出的光学参数值和重构出的目标尺寸大小，其中保真度定义为 ( max μconstruct−μB)/(μobject−μB)μB，为背景组织的吸收或约化散射系数．

设图1中背景组织体的光学参数为μaB=0.01mm−1、= 1 mm−1，假设在上述成像域内有2个病灶，半径均为 2,mm．病灶的光学参数分别为：μa1=3μaB，μs′1= μs′B，μa2=μaB，μs′2= 2 μs′B．病灶中心坐标分别为(－14，0)和(14，0)，如图 2(a)和(b)所示．

图2 aμ和′sμ的目标图像和重构图像(aμ对比度3∶1，sμ′对比度2∶1)Fig.2 Image ofaμand′sμwith target and recon-structed images from simulated data(contrast levels foraμ and sμ′ are 3∶1 and 2∶1，respectively)

将图 2(c)和(d)的重构图像和目标图像进行对比可见，重构目标病灶方位和尺寸基本同目标图像一致．并且重构目标病灶的位置基本都在真实边界内(实线圆)，证明了重构位置的准确性．分别取图 2(a)和(b)沿过病灶中心的水平方向和垂直方向上的目标值和重构值进行比较，结果如图3所示．可见aμ和sμ′重构的保真度可达81.0%和81.2%．并且aμ图像中，水平方向和垂直方向的重构目标值半高宽分别为3.50,mm和4.67,mm，sμ′2个方向的半高宽分别为3.56,mm和4.83,mm，与真实病灶直径(4,mm)对比可见，重构的目标病灶尺寸基本准确．

图3 图2中水平和垂直方向上目标值和重构值比较Fig.3 Comparison of reconstruction accuracies along the horizontal dashed line and the vertical dashed line in Fig.2

为进一步验证算法的正确性，令病灶与背景组织体的光学参数的对比度发生变化，μa1= 2 μaB，μs′1= μs′B，μa2=μaB，μs′2= 3 μs′B，其他条件与图 2图像重构条件类似，目标图像如图 4(a)和(b)所示．重构图像见4(c)和(d)．

将重构图像和目标图像进行对比可见，重构位置准确．同样，分别取图 4(a) 和(b)沿过病灶中心的水平方向和垂直方向上的目标值和重构值进行比较，结果如图5所示．μa和μs′重构的保真度分别可达82%和66.7%．另外μa图像中水平方向和垂直方向重构目标值的半高宽分别为3.50,mm和4.83,mm，μs′ 2个方向的半高宽分别为 3.81,mm 和 4.67,mm．可见μa和μs′的重构目标值的半高宽基本同真值的直径(4,mm)相同，证明重构目标病灶尺寸基本准确．

图4 aμ和′sμ的目标图像和重构图像(aμ对比度2∶1，sμ′对比度3∶1)Fig.4 Image ofaμand′sμwith the target and the reconstructed images from the simulated data(the contrast levels foraμandsμ′are 2∶1 and 3∶1，respectively)

另外，在对比度与图 2类似的条件下，还分别对背景组织体光学参数的变化、病灶方位的变化、病灶半径的变化对重构精度的影响进行了验证，结果如表1～表3所示．表中R为病灶的半径，(x，y)为病灶的中心坐标，Lx、Ly分别为与x、y轴平行且过病灶中心直线上重构图像的半高宽．表1为R＝2,mm，μa的(x，y)＝(－14，0)，μs′的(x，y)＝(14，0)时，背景组织体光学参数的变化对重构精度的影响；表2为R＝2,mm，μ =0.01mm−1，μ′= 1 mm−1时，病灶方位的aBsB变化对重构精度的影响；表 3为 μ =0.01mm−1，

aBμ′=1mm−1，μ的(x，y)＝(－14，0)，μ′的(x，y)＝sBas(14，0)时，病灶半径R的变化对重构精度的影响．

图5 图4中水平和垂直方向上目标值和重构值比较Fig.5 Comparison of reconstruction accuracies along the horizontal dashed line and the vertical dashed line in Fig.4

表1 背景组织体光学参数不同对aμ和′sμ重构精度影响Tab.1 Independent reconstruction ofaμand′sμfrom simulated data with different background optical properties

表2 病灶方位不同对aμ和′sμ重构精度影响Tab.2 Independent reconstruction ofaμand′sμfrom simulated data with different positions of target

表3 病灶半径不同对aμ和′sμ重构精度影响Tab.3 Independent reconstruction ofaμand′sμfrom simulated data with different sizes of target

由表 1～表 3可见，背景组织体光学参数、病灶方位分别发生变化时，重构的aμ和sμ′保真度的变化都在3%以内，并且水平方向和垂直方向上重构值的半高宽几乎不受上述因素的影响；但病灶半径发生变化时，病灶半径越大，重构精度越高，而病灶半径变小时，重构保真度下降严重且重构目标病灶的尺寸展宽严重．

3 目标深度对目标定位精度的影响

为了和文献[8]进行对比，在以下图像重构中，仅利用频域的DC数据对吸收系数进行图像重构，且设仿组织体的内半径为 10,mm，外半径为 30,mm，源点和探测点均为8个，和图1排列方式相同，背景组织体的光学参数为 μ =0.002mm−1，μ′= 0 .5mm−1，病aBsB灶的直径为7.70,mm．在2种吸收系数对比度下进行图像重构，其中 μa=100μaB，μs′= μs′B代表无限吸收对比度，μ =0.0059mm−1，μ′= 1 .03mm−1代表有限吸收as对比度．

文献[8]中无限吸收对比度的仿真结果显示：当病灶中心在(0，15)时，重构目标病灶的位置基本正确；当中心在(0，20)、(0，25)时，重构目标病灶的位置几乎仍在(0，15)不变，重构出的病灶尺寸变大．在上述条件下，采用本文的算法进行图像重构，结果显示：中心坐标在(0，15)、(0，20)、(0，25)这 3 个位置处，重构目标病灶的方位与真值相同．此结果表明在无限吸收对比度的情况下，本文的算法所重构的结果在目标方位上不受目标深度的影响．

文献[8]中有限吸收对比度的仿真条件为：背景组织体的光学参数不变，病灶中心的坐标分别为(0，15)、(0，17)、(0，19)，如图 6(a)～(c)所示．文献[8]的仿真结果显示：病灶中心在(0，15)时，重构目标病灶的位置基本正确．当中心在(0，17)时，重构目标病灶的位置仍在(0，15)．当中心在(0，19)时，重构目标病灶的中心在(10.6，－10.6)，重构的目标病灶方位产生了很大的偏差．

在上述条件下，采用本文的算法所重构的图像如图 6(d)～(f)所示．可见病灶中心在(0，15)、(0，17)、(0，19)时，重构目标病灶的位置基本与真值位置(实线圆)重合．此结果证明了在有限吸收对比度的情况下，本文算法所重构的结果在目标方位上也不受目标深度的影响．3个位置处aμ重构的保真度分别为 97.4%、69.2%、51.3%，可见随着病灶离内边界距离的增大，aμ重构保真度变差．当病灶中心位于(0，15)时，重构目标病灶的水平和垂直方向的半高宽分别为 7.70,mm 和 5.50,mm；(0，17)时，半高宽分别为10.16,mm 和 8.80,mm；(0，19)时，半高宽分别为13.86,mm和 9.90,mm．与真实直径(7.70,mm)对比可见，随着病灶离内边界距离的增大，重构的病灶尺寸也随之增大．

为了改善保真度和重构目标尺寸的准确性，本文还同时利用频域测量值的幅值和相位信息进行图像重构，结果如图7所示．

图6 病灶处于不同深度时仅利用DC数据的图像重构结果Fig.6 Reconstructed images for targets at different depths using only DC data

可见重构目标病灶的位置完全在真实边界内(实线圆)．3个位置处aμ重构的保真度分别为112.8%、89.7%、66.7%；当病灶中心位于(0，15)时，重构目标病灶的水平和垂直方向的半高宽分别为 7.70,mm和7.43,mm；(0，17)时，半高宽分别为 9.15,mm 和8.21,mm；(0，19)时，半高宽分别为 11.17,mm 和9.75,mm．将图 6和图 7的重构结果进行对比可见，同时利用频域的幅值和相位信息，不仅保真度有所提高，而且重构目标病灶的尺寸更接近于真实病灶．

图7 病灶处于不同深度时同时利用频域幅值和相位信息的图像重构结果Fig.7 Reconstructed images for targets at different depths using both amplitude and phase parts of frequency domain

4 结 语

针对于用于早期宫颈癌诊断的内窥式近红外频域漫射层析成像，本文研究了同时利用频域测量值中幅值和相位信息的重构方法．正问题采用有限元法求解，逆问题采用高斯牛顿理论．雅可比矩阵采用伴随源法进行构建，并采用修正的广义脉冲谱技术．迭代更新因子采用GMRES Krylov方法求解．

图像的重构结果表明：对aμ和sμ′进行单独重构时，重构目标的位置、大小准确度较高．虽然重构值的精度受对比度、背景组织体的光学参数、病灶方位、病灶大小等因素的影响，但重构aμ和sμ′保真度可达 80%，证明了本文算法的正确性；为了与文献[8]进行比较，在无限吸收对比度和有限吸收对比度这 2种情况下仅利用 DC数据进行了图像重构，结果表明：本文所发展的算法可消除目标定位精度受目标深度影响的问题．并且同时利用频域测量值的幅值和相位信息可提高重构的精度．

［1］ Schweiger M，Arridge S R，Nissila I. Gauss-Newton method for image reconstruction in diffuse optical tomography[J]. Physics in Medicine and Biology，2005，50(10)：2365-2386.

［2］ Arridge S R. Optical tomography in medical imaging [J].Inverse Problems，1999，15：R41-R93.

［3］ Brooksby B，Jiang S，Dehghani H，et al. Combining near-infrared tomography and magnetic resonance imaging to study in vivo breast tissue：Implementation of a Laplacian-type regularization to incorporate magnetic resonance structure [J]. Journal of Biomedical Optics，2005，10(5)：1-10.

［4］ Xu H，Dehghani H，Pogue B W，et al. Near-infrared imaging in the small animal brain：Optimization of fiber positions[J]. Journal of Biomedical Optics，2003，8(1)：102-110.

［5］ Mirabal Y N，Chang S K，Atkinson E N，et al. Reflectance spectroscopy for in vivo detection of cervical precancer[J]. Journal of Biomedical Optics，2002，7(4)：587-594.

［6］ Kendrick J E，Huh W K，Alvarez R D，LUMATMcervical imaging system[J]. Expert Review of Medical Devices，2007，4(2)：121-129

［7］ Piao Daqing，Xie Hao，Musgrove C，et al. Nearinfrared optical tomography：Endoscopic imaging approach [J]. Proceedings of SPIE，2007，6431：1-10.

［8］ Musgrove C，Bunting C F，Dehghani H，et al. Computational aspects of endoscopic(Trans-rectal)nearinfrared optical tomography：Initial investigations[J].Proceedings of SPIE，2007，6434：1-10.

［9］ Schweiger M，Arridge S R，Hiraoka M，et al. The finite element method for the propagation of light in scattering media：Boundary and source conditions[J].Medical Physics，1995，22(11)：1779-1792.

［10］ Gao Feng，Zhao Huijuan，Tanikawa Y，et al. Timeresolved diffuse optical tomography using a modified generalized pulse spectrum technique [J]. Ieice Transactions on Information and Systems，2002，E85-D(1)：133-142.

［11］ Sebbah P. Waves and Imaging Through Complex Media[M]. Holland：Klumer，2001.

［12］ Wang Lihong V，Wu Hsin-i. Biomedical Optics[M].New York：John Wiley and Sons Inc，2007.