第４４卷第９期 天津大学学报 Ｖｏｌ ４４ ＮＯ．９ ２０１１年９月 Ｊｏｕｒｎａｌ ｏｆ Ｔｉａｎｊｉｎ Ｕｎｉｖｅｒｓｉｔｙ Ｓｅｐ．２０１１ 基于自然边界条件求解辐射传输方程的方法 金 蒙 ，高 峰１，２，马文娟 ，杨 芳 ，赵会娟。， （１．天津大学精密仪器与光电子＿Ｌ程学院，天津３０００７２： ２．天津市生物医学检测技术与仪器重点实验室，天津３０００７２） 摘要：基于零内向流边界条件的辐射传输方程是目前被广泛采用的描述有限组织体中光子传播过程的模型 由于 生物组织体及其周围介质之间折射率的差异，零内向流边界条件无法准确模拟光子到达边界时的传播过程，从而导 致了实际应用中对光学参数的重建存在明显误差．为解决这一问题，建立了考虑边界内反射效应的自然边界条件下 的辐射传输方程模型，发展了一种结合立体角元离散法和有限差分法的自然边界条件稳态辐射传输方程求解算法， 并通过与蒙特卡洛模拟结果的对比验证了该算法的正确性 关键词：辐射传输方程；有限差分方法；边界条件；扩散光学层析成像 中图分类号：Ｑ６３ 文献标志码：Ａ 文章编号：０４９３—２１３７（２０１ １）０９一Ｏ８１６ ０７ Ａ Ｆｉｎｉｔｅ Ｄｉｆｆｅｒｅｎｃｅ Ｍｅｔｈｏｄ ｔｏ Ｓｏｌｖｅ Ｒａｄｉａｔｉｖｅ Ｔｒａｎｓｆｅｒ Ｅｑｕａｔｉｏｎ ｗｉｔｈ Ｎａｔｕｒａｌ Ｂｏｕｎｄａｒｙ Ｃｏｎｄｉｔｉｏｎ ＪＩＮ Ｍｅｎｇ ，ＧＡＯ Ｆｅｎｇ ２，ＭＡＷｅｎ－ｊｕａｎ，ＹＡＮＧ Ｆａｎｇ 。ＺＨＡＯ Ｈｕｉ－ｊｕａｎ ＇。 （１．Ｓｃｈｏｏｌ ｏｆ Ｐｒｅｃｉｓｉｏｎ Ｉｎｓｔｒｕｍｅｎｔ ａｎｄ Ｏｐｔｏ—Ｅｌｅｃｔｒｏｎｉｃｓ Ｅｎｇｉｎｅｅｒｉｎｇ，Ｔｉａｎｊｉｎ Ｕｎｉｖｅｒｓｉｔｙ，Ｔｉａｎｊｉｎ ３０００７２，Ｃｈｉｎａ ２．Ｔｉａｎｊｉｎ Ｋｅｙ Ｌａｂｏｒａｔｏｒｙ ｏｆＢｉｏｍｅｄｉｃａｌ Ｄｅｔｅｃｔｉｎｇ Ｔｅｃｈｎｉｑｕｅｓ ａｎｄ Ｉｎｓｔｒｕｍｅｎｔｓ，Ｔｉａｎｊｉｎ ３０００７２，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｔｈｅ ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ ｅｑｕａｔｉｏｎ（ＲＴＥ）ｗｉｔｈ ｔｈｅ ｉｎｆｌｏｗ．ｚｅｒｏ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ ｉｓ ｃｕｒｒｅｎｔｌｙ ｉｎｃｒｅａｓｉｎｇｌｖ ａｄｏｐｔｅｄ ａｓ ａ ｕｎｉｖｅｒｓａｌｌｙ ａｐｐｌｉｃａｂｌｅ ｐｈｏｔｏｎ—ｍｉｇｒａｔｉｏｎ ｍｏｄｅｌ ｉｎ ｔｈｅ ｔｉｓｓｕｅｓＨｏｗｅｖｅｒ，ｓｉｎｃｅ ｂｉｏｌｏｇｉｃａｌ ｔｉｓｓｕｅｓ ｈａｖｅ ．ｎｏｒｍａｌｌｙ ｄｉｆｆｅｒｅｎｔ ｒｅｆｒａｃｔｉｖｅ ｉｎｄｅｘｅｓ ｆｒｏｍ ｔｈｅ ｓｕｒｒｏｕｎｄｉｎｇｓ ｓｕｃｈ ａｓ ａｉｒｔｈｉｓ ｉｎｆｌｏｗ—ｚｅｒｏ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ ｃａｎ ｎ０ｔ ，ａｃｃｕｒａｔｅｌｙ ｍｏｄｅｌ ｔｈｅ ｂｅｈａｖｉｏｒ ｏｆ ｐｈｏｔｏｎｓ ｔｈａｔ ｈｉｔ ｔｈｅ ｂｏｕｎｄａｒｙ，ｌｅａｄｉｎｇ ｔｏ ｌａｒｇｅ ｅｒｒｏｒｓ ｉｎ ｒｅｔｒｉｅｖａｌ ｏｒ ｉｍａｇｅ ｒｅｃｏｎ— ｓｔｒｕｃｔｉｏｎ ｏｆ ｔｈｅ ｏｐｔｉｃａｌ ｐｒｏｐｅｒｔｉｅｓ．Ｔｏ ｏｖｅｒｃｏｍｅ ｔｈｅ ｌｉｍｉｔａｔｉｏｎｓ，ｔｈｅ ｎａｔｕｒａｌ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ ｔｈａｔ ｔａｋｅｓ ｉｎｔｏ ａｃｃｏｕｎｔ ｔｈｅ ｉｎｔｅｒｎａｌ ｒｅｆｌｅｃｔｉｏｎ ｅｆｆｅｃｔ ｓｈｏｕｌｄ ｂｅ ｃｏｎｓｉｄｅｒｅｄ．Ａｊｏｉｎｔ ｄｉｓｃｒｅｔｅ—ｏｒｄｉｎａｔｅ ａｎｄ ｉｆｎｉｔｅ—ｄｉｆｆｅｒｅｎｃｅ ｓｏｌｕｔｉｏｎ ｔｏ ｔｈｅ ｓｔｅａｄｖ— ｓｔａｔｅ ＲＴＥ ｗｉｔｈ ｔｈｅ ｎａｔｕｒａｌ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ ｉｓ ｐｒｏｐｏｓｅｄ ａｎｄ ｖａｌｉｄａｔｅｄ ｂｙ ｃｏｍｐａｒｉｓｏｎ ｗｉｔｈ ｔｈｅ Ｍｏｎｔｅ—Ｃａｒｌｏ ｓｉｍｕｌａ— ｔｉｏｎｓ． Ｋｅｙｗｏｒｄｓ：ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ ｅｑｕａｔｉｏｎ；ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅ ｍｅｔｈｏｄ；ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ；ｄｉｆｆｕｓｅ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ 扩散光学层析成像（ｄｉｆｆｕｓｅ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ． 性强等多重优点，因此在新生儿的脑血氧监测、乳腺 肿瘤的早期检测和基于小动物疾病模型的分子成像 ＤＯＴ）是一种新兴的无创功能成像模态ｌＪ 】，其基本原 理是：使用波长范围在６５０～９００ ｒｉｍ的近红外光“扫 捕”组织体，通过测量组织体表面上的光强分布来重 ０刈等领域具有广泛的重要应用前景．ＤＯＴ研究广泛 采用基于模型的图像重建方法（ｍｏｄｅｌ—ｂａｓｅｄ ｉｍａｇｅ 建组织体内部的光学参数图像．ＤＯＴ的成像对象是 ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ ｓｃｈｅｍｅ）［８１，其中扩散方程（ｄｉｆｆｕｓｉｏｎ ｅｑｕａｔｉｏｎ，ＤＥ）和辐射传输方程（ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ ｅｑｕａｔｉｏｎ，ＲＴＥ）为目前应用的２种主要光子输运正向 反映组织体病生理功能信息的光学参数（吸收和散射 系数），具有灵敏度和特异性高、安全廉价以及实时 收稿日期： ２０ｌｏ一１ｏ－】２；修回日期：２０１０—１２．０６， 基金项目： 国家自然科学基金资助项目（３０８７０６５７，３０９７０７７５）；国家重点基础研究发展计划（９７３计划）资助项目（２００６ＣＢ７０５７００）；ｌ目 家高技术研究发展计划（８６３计划）资助项目（２００９ＡＡ０２Ｚ４１３）；天津市自然科学基金资助项目（０９ＪＣＺＤＪＣＩ ８２００， １ＯＪＣＺＤＪＣｌ７３００）． 作者简介： 金通讯作者： 高蒙（１９８０一 ），男，博 研究生，ｋｉｎｇｍａｎｌ９８０＠ｈｏｔｍａｉｌ．ＣＯＩＩ１． 峰，ｇａｏｆｅｎｇ＠ｔｊｕ．ｅｄｕ．ｃｎ． ２０１１年９月 金蒙等：基于自然边界条件求解辐射传输方程的方法 模型．扩散方程是辐射传输方程的球谐函数一阶近 界某一位置上垂直于边界向内的方向上，存在一准直 光源，用 函数可以表示为 似，其适用范围有相当的局限性：仅适用于高散射、 低吸收组织体的远源区光子密度场的描述【１ Ｊ．因此， 方面，对于小尺寸组织体（如小动物）的成像应用， 由于光源和探测器的分布较密，扩散方程的适用性较 差，影响近源场内测量信息的利用和光学参数的重 一ｑ（ｒ， ）＝ （，一ｒｏ）８（ｓ—Ｓ０） ｒｏ∈ａ ，ｓ０：一，ｌ （３） 式中：ｒｎ为光源所在的边界上的一点；ｎ为位于边界 建；另一方面，在生物组织内的某些区域，如高散射、 上源位置处的外法向方向单位向量． 零内向流边界条件为 高吸收的骨骼组织、低散射、低吸收的脑脊液层以及 腹部大量的空腔域等，扩散方程无法实现光子传播行 为的有效建模．为此，国内外学者一直对基于辐射传 输方程的ＤＯＴ成像技术有着广泛的兴趣．目前，基 于辐射传输方程的光子传播模型普遍基于零内向流 边界条件（ｉｎｆｌｏｗ—ｚｅｒｏ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ，ＩＺＢＣ），即 假设边界处向组织体内部的辐射率为零【９ ”．应用该 边界条件可以为求解辐射传输方程提供极大便利，但 由于ＩＺＢＣ未考虑组织体及其周围介质之间折射率 的差异，因此，忽略了光在实际组织体边界上的折射 和反射，从而造成了正向模型和光的实际传播物理过 程之问存在较大误差．在基于扩散方程的扩散层析 成像中，已经发现了组织体和周围介质的折射率差异 对组织体内部光学参数的重建有很大的影响ｌｌ引． 为了准确构建光在组织体内传播过程的正向模 型，必须考虑采用能够有效反映边界反射效应的自然 边界条件（ｎａｔｕｒａｌ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎ，ＮＢＣ）．笔者提 出了一种采用离散立体角元和有限差分方法求解基 于自然边界条件的辐射传输方程，并将数值求解结果 和蒙特卡洛模拟结果进行了比对，验证了所提方法的 正确性． １辐射传输方程 在组织体区域中，精确描述光子输运过程的稳态 辐射传输方程 为 （ ・Ｖ＋ （，）＋ （ｒ）） （ ）＝ ∥ （，．）ｆｐ（Ｓ，Ｓ’） （ ’）出’ （１） ４ｎ 式中： （ｒ，ｓ）表不空ｌ司位置ｒ上、方向ｓ上的辐射率， Ｗ／（ｍｍ ．ｓｒ）；Ｕｓ和 分别为散射系数和吸收系数， ｌｎｉｎ一；ｐ（ｓ， ’１为散射相位函数．本文采用广泛使用 的Ｈｅｎｙｅｙ—Ｇｒｅｅｎｓｔｅｉｎ函数㈣，即 １） 式中ｇ为各向异性因子．式（１）中不含源项，考虑边 （，，Ｓ）＝０ Ｓ・，ｌ＜０，ｒＥ ａ （４） 该边界条件的物理意义为：当光子抵达组织区域 边界时，将从组织体表面彻底溢出，不再返回．显然， 该边界条件仅仅是一种理想化的边界条件，其有效性 建立在组织体与其周围媒质光折射率相等的基础 上．由于实际中组织体与其周围介质之间折射率差 异明显，在组织体边界处必然会发生光折射和反 射．只有发生折射的光子能溢出边界，而反射的光子 则返回组织体内继续传播，从而影响光子密度分 布．显然，应用ＩＺＢＣ的辐射传输方程作为ＤＯＴ成 像的正向模型势必引起图像重建结果的显著误 差．为了正确模拟组织体内光子输运行为，发展基于 ＮＢＣ的辐射传输方程建模方法十分必要． 光在到达边界处的入射角 和折射角 遵循斯 涅耳折射定律（见图１），表示为 ｓｉｎ ７／Ｔ ｓｉｎ缉 一Ｉ ｓｉｎ ＜ （５） ，ｚＩ 式中： 为入射角； 为折射角； 为反射角；ｒ／ 为 模拟组织体的折射率； 为组织体包围介质的折射 率，本文中所讨论的情况仅仅涉及 ≥ 的情况． ａ （ｒ，ｓ） 图１斯涅耳折射定律 Ｆｉｇ．１ Ｓｎｅｌｌ’Ｓ ｌａｗ ｏｆ ｒｅｆｒａｃｔｉｏｎ 根据菲涅尔公式 】，可以得到反射光的辐射率 （，， ）和入射光辐射率 （，， ）的公式为 （，．，Ｓ ）＝ （ ） （，，Ｓ ） ｒＥ ａ （６） 式中 天 津 大 学 学 报 第４４卷第９期 』Ｒ（ ）＝Ｒ（ ） 厂＝ｌ （，．， ） 一Ｉ （，， ） ＝ ５ ０ ｓｎ’＜Ｏ ２ｌ｝ ｓｉｎ ＋（ ＋ ）。 ｔａｎ 。（ ＋珥） Ｉｓ…ｉｎ ＜ 、 等 ∑（１～ ）哝（ｒ）』ｄｓ 趣 （１ １） １ ｓｉｎ ≥ 称为菲涅尔反射系数．式（６）即为自然边界条件．根 据能量守恒定律，可以得到自然边界条件下的表面外 向光子流量率，即探测器在组织体表面测得的光子流 率，其公式为 厂＝Ｉ。 （，， ）出一Ｉ （，， ） （７） ｎ＞０ ０ 当组织与周边媒质光折射率相等时，Ｒ＝０，于 是式（７）第２项为０，从而退化为零内向流边界条件的 情形． ２数值求解 在笔者以前的工作中，采用了联合的立体角元离 散法和有限差分法在零内向流边界条件下对辐射传 输方程进行了求解【】 ．对于自然边界条件下的辐射 传输方程进行数值求解需要在此基础上进行改 进．首先，需要对自然边界条件进行角度离散，本文 对自然边界条件采用了离散立体角元的方法进行离 散．根据文献［１６］中提出的立体角离散方法，对于入 射光所在的某一立体角元 肯定可以找到与入射面 对称的反射光所对应的立体角元 ，且满足关系 ．ｆ ＝』ｄｓ （８） ｎ 假设在某一离散立体角元 内，入射光辐射率 为一恒定值，对式（６）在离散立体角元 内做立体角 积分，可以得到 ．ｆ （ｒ，Ｓ）出＝』 （ ） （ ）出 ，∈ａ （９） ａ＇２Ｒ ．ｃ２１ 对于任意一个入射光方向所在的离散立体角元 ，可以得到离散的反射系数为 ：譬 采用文献［１６】中离散立体角元方法可以求得对 应于不同离散立体角元的 ．同样，对应于自然边界 条件的表面外向光子流率可以离散为 在文献［１６］中提到的采用有限差分方法求解辐 射传输方程的迭代过程中，根据解得的表面外方向的 辐射率以及式（１ １），可以求解出表面内方向上的辐射 率，表面内方向上的辐射率将影响组织体内部光子密 度的分布，具体过程如下． 步骤１ 根据源的位置确定迭代初始值 。． 步骤２第，＋１次迭代，迭代过程为 ：：兰 ：兰 。 ＋等州ｒ）＋州］ ． ．竿 ＋华 ｈ ＋ ＋ ，）＋ （，）Ｉ Ｊｈ ］ ｆ “ ： ｆ ． ．Ｓ ｎ 州 一 Ｓ ｎ 西ｋ＜０，ｚ｜＜０ ｈ Ｓ ～ ｈ ＋ ，）＋ （，）］ ’ ｎ ｄ ｄ （，）∑ ． ＋ 一 ■膏 『Ｉ ｆ。 ｆＣ 《ｌ 式中： 妒 ∞ ｓｉｎ Ｏｄ０ｄＯ； ：Ｉ Ｉ ｓｉｎ ０ｄ０；０表示纬角；妒表示 经角； 、 、 、 分别表示第足个离散立体角 元对应的纬角和经角的积分上下限；ｈ为有限空间 格步长．该迭代的推导过程可以参考文献［１６］，此处 不再赘述． 步骤３ 根据第，＋１次迭代，得到的 “ 及式 （１０），更新 对应于介质表面内方向上的辐射率 “ ＝ 川Ｈ 这里 ， 表示表面位置上向外 方向上的立体角元内的辐射率， ¨ ，表示 … 与介质平面对称的立体角元内的辐射率． ２０１１年９月 金蒙等：基于自然边界条件求解辐射传输方程的方法 步骤４根据终止判据判断是否终止迭代，如果 终止，最后一次迭代的 即为最终数值求解结果，根 据式（１１）可以获得介质表面的表面光子流率．如果 未完成迭代则返回步骤２． 根据上述方法，可以求得有限空间内基于自然边 界条件的辐射传输方程的数值解． ３模拟验证 针对小动物成像问题，在大小为２ ｃｍ×２ Ｃ１ＴＩ的 二维模拟组织体内，对稳态辐射传输方程在２种边 界条件下进行了数值求解．模拟组织体的光学参数 为 ＝１０ ｍｍ～， ＝０．０１ ｍｍ，ｇ＝０．９，分别取３种 折射率ｎ．．＝１．５６，７／ ，＝１．４５…７／＝１．３８，假设模拟组织体 被空气包围，即 ：１．００．光源位于模拟组织体一边 的中心，定义在光源对面的面为透射面，在光源侧面 为侧面，如图２所示．为了验证求解辐射传输方程算 法的正确性，对Ｗａｎｇ ＂Ｊ提出的无限半空间内的光 子传播蒙特卡洛模拟程序进行了扩展，加入了二维边 界，当光子溢出边界时，记录光子的溢出位置，并将 相同位置上的光子权重进行累加，获得该位置上的表 面光子流率．通过蒙特卡洛模拟得到的侧面和透射 面上的表面光子流率，如图３所示．同时，根据辐射 传输方程数值求解的结果以及式（１１）获得模拟组织 体的表面光子流率，然后将结果和蒙特卡洛模拟结果 进行了比对，比对结果如图４所示．从图３中可以看 到，当模拟组织体的折射率从１．３８增加到１．５６的过 程中，在侧面上和透射面上得到的表面光子流率随着 折射率的增加而增加．当模拟组织体折射率为１．５６ 时，在侧面上的表面光子流率比模拟组织体折射率为 １．４５时约高６．２４％，比模拟组织体折射率为１．３８时约 高１１．４８％，比应用零内向流边界条件时约高 ６２．５５％．在透射面上的表面光子流率也有类似的分 布，当模拟组织体折射率为１．５６时，透射面上的表面 侧面 准直光源 射面 图２模拟仿体形状和光源位置 Ｆｉｇ．２ Ｓｈａｐｅ ｏｆ ｓｉｍｕｌａｔｅｄ ｐｈａｎｔｏｍ ａｎｄ ｓｏｕｒｃｅ ｐｏｓｉｔｉｏｎ ４－＂ 吕 旨 ● 熏 薄 旧 槲 （ａ）侧面 －＾ ｇ 吕 ● ０ 料 避 恒 （ｂ）透射面 ——ＭＣ ｗｉｔｈ ＮＢＣ（ｎ『＝１ ５６） ………… ＭＣ ｗｉｔｈ ＮＢＣ（ｎ＿＝１ ３８） 一一一一ＭＣ ｗｉｔｈＮＢＣ（ｎ；１．４５、 一…一ＭＣｗｉｔｈＩＺＢＣ 图３根据蒙特卡洛模拟得到的表面光子流率 Ｆｉｇ．３ Ｏｕｔｇｏｉｎｇ ｐｈｏｔｏｎ ｄｅｎｓｉｔｙ ａｃｑｕｉｒｅｄ ｆｒｏｍ ＭＣ ｓｉｍｕｌａ— ｔｉｏｎｓ 光子流率比模拟组织体折射率为１．４５时约高７％，比 模拟组织体折射率为１．３８时约高１２．３７％，比应用零 内向流边界条件约高６２．５５％．值得一提的是，在图 ３（ａ）中，应用零内向流边界条件的蒙特卡洛模拟得到 的侧面的表面光子流率的峰值相对于应用自然边界 条件的蒙特卡洛模拟得到的侧面的表面光子流率的 峰值在Ｙ轴方向上相对于源的位置远大约１ ｍｍ． 图４为数值求解辐射传输方程所得之透射面和反 射面外向光子流率和蒙特卡洛模拟结果的比对结果， 从图４中可以看到，在考虑组织内外折射率差异的自 然边界条件及零内向流边界条件下数值求解结果分 别与其相应蒙特卡洛模拟结果相吻合，从而证明了本 文所提出的应用离散立体角元和有限差分方法求解 自然边界条件下的辐射传输方程的算法的正确性． 图５所示为在零内向流边界条件和自然边界条 件（ ＝１．５６，／ＴＴ＝１．００）下求解辐射传输方程所得的模 拟组织体内部光子密度分布．可以看出，自然边界条 件下所得到的组织体内部的光子密度比零内向流边 界条件下的值平均高出１７８．２６％，因此，不同的边界 条件对组织体内部的光子密度分布影响极大． ・８２０・ 天 津 大 学 学 报 第４４卷第９期 爵妮 ０ Ｏ ８ Ｏ ６ 恒 ０ ４ ０ ２ ０ 诲避 Ｏ ８ Ｏ ６ 恒 ０ ４ Ｏ ２ ０ 醑 ０ ８ 恒 ０ ６ 静蛙 Ｏ ２ ０ 茸 ０ ６ ０ ４ ０ ８ ０ ４ ０ ２ ０ Ｏ （ａ）侧面＂１＝１，５６ （ｂ）侧面Ｆ／ｉ＝１．４５ ０ Ｏ （Ｃ）侧面ｎ１＝１．３８ （ｄ）侧面ＩＺＢＣ 好堰 ０ Ｏ ８ Ｏ ６ 恒 Ｏ ４ 霉 Ｏ ２ １ ０ 恒 ０ ８ 醑 Ｏ ２ 恒 ０ ６ Ｏ ４ Ｏ ２ Ｏ 料避 Ｏ ８ 旧 ０ ６ ０ ６ Ｏ ４ ０ Ｏ ８ ０ ４ Ｏ ２ Ｏ ｘ／ｍｍ （ｅ）透射面月１＝１．５６ （ｆ）透射面／／＇１＝１．４５ ０ ｘ／ｍｍ ０ ／ｉｉｌｍ （ｇ）透射而 ＝１．３８ （ｈ）透射面ＩＺＢＣ 图４蒙特卡洛模拟和辐射传输方程数值解得到的归一化表面光子流率比对 Ｆｉｇ．４ Ｃｏｍｐａｒｉｓｏｎｓ ｏｆ ｎｏｒｍａｌｉｚｅｄ ｏｕｔｇｏｉｎｇ ｐｈｏｔｏｎ ｄｅｎｓｉｔｙ ａｃｑｕｉｒｅｄ ｆｒｏｍ ＭＣ ｓｉｍｕｌａｔｉｏｎｓ ａｎｄ ＦＤＭ ｓｏｌｕｔｉｏｎ ｔｏ ＲＴＥ ｗｉｔｈ ＮＢＣ ａｎｄＩＺＢＣ ２０ １５ 量 １０ ５ Ｏ ５ ｌ０ ｙ／ｍｍ １５ ２Ｏ ｙ／ｍｍ （ａ）零内向流边界条件 （ｂ）自然边界条件 图５不同边界条件下点光源激励下均匀组织体内的光子密度分布 Ｆｉｇ．５ Ｃｏｍｐａｒｉｓｏｎｓ ｏｆ ｐｈｏｔｏｎ ｄｅｎｓｉｔｙ ｉｎｓｉｄｅ ｔｉｓｓｕｅ ａｃｑｕｉｒｅｄ ｆｒｏｍ ＦＤＭ ｓｏｌｕｔｉｏｎ ｔｏ ＲＴＥ ｗｉｔｈ ＩＺＢＣ ２０１１年９月 金蒙等：基于自然边界条件求解辐射传输方程的方法 ・８２１・ 本文中模拟了２个非均匀模拟组织体在如图２ 所示的准直光源激励的情况下，对辐射传输方程在零 内向流边界条件下和自然边界条件下进行了数值求 解．非均匀模拟组织体内部的吸收系数和散射系数 分布分别如图６（ａ）～６（ｄ）所示．零内向流条件下获得 的模拟组织体内部的光子密度分布分别如图６（ｆ）、 （ｈ）所示，自然边界条件下得到的模拟组织体内部的 光子密度分布如图６（ｅ）、（ｇ）所示．在２种不同边界 条件下得到的２个模拟组织体表面上的表面光子流 率分布如图７所示． 一 ｇ５． 一、槲斌 恒 一 ｇ． ０Ｈ一、斟煺 米恒懈 ｙ／ｍｍ ｙ／ｒａｍ （ｇ）组织体Ｂ内光子密度分布（ＮＢＣ）（ｈ）组织体Ｂ内光子密度分布（ＩＺＢＣ） 图６非均匀组织体内部光学参数分布以及点光源激励下内 部光子密度分布 Ｆｉｇ．６ Ｏｐｔｉｃａｌ ｃｏｅｆｉｆｃｉｅｎｔｓ ｄｉｓｔｒｉｂｕｔｉｏｎ ｏｆ Ｔｉｓｓｕｅ Ａ ａｎｄ Ｂ， ａｎｄ ｐｈｏｔｏｎ ｄｅｎｓｉｔｙ ｉｎｓｉｄｅ Ｔｉｓｓｕｅ Ａ ａｎｄ Ｂ ｗｊｔｈ ｄｉｆ－ ｆｅｒｅｎｔ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎｓ （ａ）组织体Ａ侧面 一 ． ． 一、褂 旧 （ｂ）组织体Ａ透射面 （ｃ）组织体Ｂ侧面 吕 鲁 ● ＇ ！ 一 懈 斌 喧 （ｄ）组织体Ｂ透射面 图７ 不同边界条件下非均匀组织体Ａ和Ｂ在点光源激励下 得到的表面光子流率 Ｆｉｇ．７ Ｏｕｔｇｏｉｎｇ ｐｈｏｔｏｎ ｄｅｎｓｉｔｙ ｏｎ ｓｉｄｅ ｓｕｒｆａｃｅ ａｎｄ ｔｒａｎｓ— ｍｉｓｓｉｏｎ ｓｕｒｆａｃｅ ｏｆ Ｔｉｓｓｕｅ Ａ ａｎｄ Ｂ ｗｉｔｈ ｄｉｆｅｒｅｎｔ ｂｏｕｎｄａｒｙ ｃｏｎｄｉｔｉｏｎｓ ４结语 提出了一种在自然边界条件下求解辐射传输方 程的离散立体角元一有限差分方法．通过对小动物尺 ・８２２・ 天 津 大 学 学 报 第４４卷第９期 寸二维组织域进行求解，并与相应的蒙特卡洛模拟结 果进行比较，证明了该算法的可行性和有效性．在此 基础上，进一步探讨了不同边界条件对基于辐射传输 方程的ＤＯＴ成像正向模型所产生的影响．研究表 ７２（５）：６９１－７１３． ［９］ Ｋｉｍ Ｈｙｕｎ Ｋｅｏｌ，Ｃｈａｒｅｔｔｅ Ａ．Ｆｒｅｑｕｅｎｃｙ ｄｏｍａｉｎ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ ｕｓｉｎｇ ａ ｃｏｎｊｕｇａｔｅ ｇｒａｄｉｅｎｔ ｍｅｔｈｏｄ ｗｉｔｈｏｕｔ ｌｉｎｅ ｓｅａｒｃｈ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｑｕａｎｔｉｔａｔｉｖｅ Ｓｐｅｃｔｒｏｓｃｏｐｙ ａｎｄ Ｒａｄｌ。ａｔｉｖｅ Ｔｒａｎｓｆｅｒ，２００７，１０４（２）：２４８—２５６． １ １Ｏ ｊ Ｂａｌｉｍａ Ｏ，Ｐｉｅｒｒｅ Ｔ，Ｃｈａｒｅｔｔｅ Ａ，ｅｔ ａ１．Ａ ｌｅａｓｔ ｓｑｕａｒｅ 明，不同边界条件下的正向模型数值解差异明显，因 此，在今后的基于辐射传输方程的ＤＯＴ逆问题研究 中应采用更符合实际的自然边界条件． 参考文献 Ａｒｒｉｄｇｅ Ｓ Ｒ．Ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ ｉｎ ｍｅｄｉｃａｌ ｉｍａｇｉｎｇ ［Ｊ］．ＩｎｖＰｒｏｂｌ，１９９９。１５（２）：４１—９３． ［２］ Ｙｏｄｈ Ａ，Ｃｈａｎｃｅ Ｂ．Ｓｐｅｃｔｒｏｓｃｏｐｙ ａｎｄ ｉｍａｇｉｎｇ ｗｉｔｈ ｄｉｆｆｕｓｉｎｇ ｌｉｇｈｔ［Ｊ］．尸 Ｔｏｄａｙ，１９９５，４８：３４—４０． ［３］ Ａ￣ｉｄｇｅ Ｓ Ｒ，Ｈｅｂｄｅｎ Ｊ Ｃ．Ｏｐｔｉｃａｌ ｉｍａｇｉｎｇ ｉｎ ｍｅｄｉ— ｃｉｎｅ（Ｐａｒｔ ２）：Ｍｏｄｅｌｌｉｎｇ ａｎｄ ｒｅｃｏｎｓｔｒｕｃｔｉｏｎ［Ｊ］．ＪＰ ＭｅｄＢｉｏｌ，ｌ９９７，４２（５）：８４１．８５３． ［４］ Ｇｉｂｓｏｎ Ａ，Ｈｅｂｄｅｎ Ｊ Ｃ，Ａｒｒｉｄｇｅ Ｓ Ｒ．Ｒｅｃｅｎｔ ａｄｖａｎｃｅｓ ｉｎ ｄｉｆｆｕｓｅ ｏｐｔｉｃａｌ ｉｍａｇｉｎｇｌＪ］．Ｐ Ｍｅｄ Ｂｉｏｌ，２００５， ５０（４）： １—４３． ［５］ Ｚｈａｏ Ｈｕｉｊｕａｎ，Ｇａｏ Ｆｅｎｇ，Ｙｕｋａｒｉ Ｔａｎｉｋａｗａ，ｅｔ ａ１． Ｔｉｍｅ—ｒｅｓｏｌｖｅｄ ｄｉｆｆｕｓｅ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｉｃ ｉｍａｇｉｎｇ ｆｏｒ ｔｈｅ ｐｒｏｖｉｓｉｏｎ ｏｆ ｂｏｔｈ ａｎａｔｏｍｉｃａｌ ａｎｄ ｆｕｎｃｔｉｏｎａｌ ｉｎｆｏｒｍａ— ｔｉｏｎ ａｂｏｕｔ ｂｉｏｌｏｇｉｃａｌ ｔｉｓｓｕｅ［Ｊ］．Ａｐｐｌ Ｏｐｔ，２００５， ４４（１０）：１９０５—１９１６． ［６］ Ｇａｏ Ｆｅｎｇ，Ｚｈａｏ ＨｕＯｕａｎ，Ｔａｎｉｋａｗａ Ｙｕｋａｒｉ，ｅｔ ａ１． Ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｉｃ ｍａｐｐｉｎｇ ｏｆ ｃｅｒｅｂｒａｌ ｈａｅｍｏｄｙｎａｍ— ｉｃｓ ｂｙ ｍｅａｎｓ ｏｆ ｔｉｍｅ—ｄｏｍａｉｎ ｄｅｔｅｃｔｉｏｎ：Ｍｅｔｈｏｄｏｌｏｇｙ ａｎｄ ｐｈａｎｔｏｍ ｖａｌｉｄａｔｉｏｎ［Ｊ］．Ｐ Ｍｅｄ Ｂｉｏｌ，２００４， ４９（６）：１０５５—１０７８． ［７］ Ｚｈａｏ Ｈ，Ｇａｏ Ｆ，Ｔａｎｉｋａｗａ Ｙ，ｅｔ ａ１．Ｔｉｍｅ—ｒｅｓｏｌｖｅｄ Ｄｉｆｆｕｓｅ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ ａｎｄ ｉｔｓ ａｐｐｌｉｃａｔｉｏｎ ｔｏ ｉｎ ｖｉｔｒｏ ａｎｄ ｖｉｖｏ ｉｍａｇｉｎｇ［Ｊ］．Ｊｏｕｒｎａｌ ｏＪ Ｂｉｏｍｅｄｉｃａｌ Ｏｐｔｉｃｓ， ２００７，１２（６）：Ｎｏ．０６２１０７． ｌ ８ ｊ Ｋｌｏｓｅ Ａ Ｄ，Ｕｗｅ Ｎ，Ｊｔｉｒｇｅｎ Ｂ，ｅｔ ａ１．Ｏｐｔｉｃａｌ ｔｏｍｏｇｒａ— ｐｈｙ ｕｓｉｎｇ ｔｈｅ ｔｉｍｅ—ｉｎｄｅｐｅｎｄｅｎｔ ｅｑｕａｔｉｏｎ ｏｆ ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ（Ｐａｒｔ １）：Ｆｏｒｗａｒｄ ｍｏｄｅｌ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｑｕａｎｔｉ— ｔａｔｉｖｅ Ｓｐｅｃｔｒｏｓｃｏｐｙ ａｎｄ Ｒａｄｉａｔｉｖｅ Ｔｒａｎｓｆｅｒ，２００２， ｉｆｎｉｔｅ ｅｌｅｍｅｎｔ ｆｏｒｍｕｌａｔｉｏｎ ｏｆ ｔｈｅ ｃｏｌｌｉｍａｔｅｄ ｉｒｒａｄｉａｔｉｏｎ ｉｎ ｆｒｅｑｕｅｎｃｙ ｄｏｍａｉｎ ｆｏｒ ｏｐｔｉｃａｌ ｔｏｍｏｇｒａｐｈｙ ａｐｐｌｉｃａｔｉｏｎｓ ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｑｕａｎｔｉｔａｔｉｖｅ Ｓｐｅｃｔｒｏｓｃｏｐｙ ａｎｄ Ｒａｄｉａｔｉｖｅ Ｔｒａｎｓｆｅｒ，２０１０，１１１（２）：２８０—２８６． ［１ １］Ｔａｒｖａｉｎｅｎ Ｔａｎｊａ，Ｖａｕｈｋｏｎｅｎ Ｍａｒｋｏ，Ｋｏｌｅｈｍａｉｎｅｎ Ｖｉｌｌｅ， ｅｔ ａ１．Ｃｏｕｐｌｅｄ ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ ｅｑｕａｔｉｏｎ ａｎｄ ｄｉｆｆｕｓｉｏｎ ａｐｐｒｏｘｉｍａｔｉｏｎ ｍｏｄｅｌ ｆｏｒ ｐｈｏｔｏｎ ｍｉｇｒａｔｉｏｎ ｉｎ ｔｕｒｂｉｄ ｍｅ— ｄｉｕｍ ｗｉｔｈ ｌｏｗ－－ｓｃａｔｔｅｒｉｎｇ ａｎｄ ｎｏｎ－－ｓｃａｔｔｅｒｉｎｇ ｒｅｇｉｏｎｓ ［Ｊ＿．Ｐ ＭｅｄＢｉｏｌ，２００５，５０（２０）：４９１３．４９３０． ［１２］Ｃｈｅｎ Ｂｉｎｇｑｕａｎ，Ｓｔａｍｎｅｓ Ｋ，Ｊａｋｏｂ Ｓ Ｊ．Ｖａｌｉｄｉｔｙ ｏｆ ｔｈｅ ｄｉｆｆｕｓｉｏｎ ａｐｐｒｏｘｉｍａｔｉｏｎ ｉｎ ｂｉｏ—ｏｐｔｉｃａｌ ｉｍａｇｉｎｇ［Ｊ］．Ａｐｐ ｐＯｔ．２００１，４０（３４）：６３５６—６３６６． １ １３］Ｄｕｄｅｒｓｔａｄｔ Ｊ Ｊ，Ｍａｒｔｉｎ Ｗ Ｒ．Ｔｒａｎｓｐｏｒｔ Ｔｈｅｏｒｙ［Ｍ Ｊ． Ｎｅｗ Ｙｏｒｋ：Ｗｉｌｅｙ，１９７９． 【１４ Ｊ Ｈｅｎｙｅｙ Ｌ Ｇ，Ｇｒｅｅｎｓｔｅｉｎ Ｊ Ｌ．Ｄｉｆｆｕｓｅ ｒａｄｉａｔｉｏｎ ｉｎ ｔｈｅ ｇａｌａｘｙ［Ｊ］．Ａｎｎａｌｅｓ ｄ＇Ａｓｔｒｏｐｈｙｓｉｑｕｅ，１９４０，３：１１７— １３７． １ １５］Ｈｅｃｈｔ Ｅ．Ｏｐｔｉｃｓ［Ｍ］．２ｎｄ ｅｄ．ＵＳＡ：Ａｄｄｉｓｏｎ ａｎｄ Ｗｅｓｌｅｙ Ｐｕｂｌｉｓｈｉｎｇ Ｃｏｍｐａｎｙ，Ｉｎｃ，１ ９８７． ［１６］金蒙，高峰，李娇，等．二维稳态辐射传输 ‘ 程的有限差分求解法［Ｊ］．光子学报，２０１０，３９（９）： １５９４。ｌ６０１． Ｊｉｎ Ｍｅｎｇ，Ｇａｏ Ｆｅｎｇ，Ｌｉ Ｊｉａｏ，ｅｔ ａ１．Ａ ｆｉｎｉｔｅ— ｄｉｆｆｅｒｅｎｃｅ—－ｍｅｔｈｏｄ ｓｏｌｕｔｉｏｎ ｔｏ ｔｗｏ・－ｄｉｍｅｎｓｉｏｎａｌ ｓｔｅａｄｙ—－ ｓｔａｔｅ ｒａｄｉａｔｉｖｅ ｔｒａｎｓｆｅｒ ｅｑｕａｔｉｏｎ［ＪＪ．Ａｃｔａ Ｐｈｏｔｏｎｉｃａ Ｓｉｎｉｃａ，２０１０，３９（９）：１５９４．１６０１（ｉｎ Ｃｈｉｎｅｓｅ）． １ １７ ｊ Ｗａｎｇ Ｌｉｈｏｎｇ，Ｊａｃｑｕｅｓ Ｓ Ｌ，Ｚｈｅｎｇ Ｌｉｑｉｏｎｇ．ＭＬＭＣ— Ｍｏｎｔｅ Ｃａｒｌｏ ｍｏｄｅｌｉｎｇ ｏｆ ｌｉｇｈｔ ｔｒａｎｓｐｏｒｔ ｉｎ ｍｕｌｔｉ—ｌａｙｅｒｅｄ ｔｉｓｓｕｅｓ［Ｊ］．Ｃｏｍｐｕｔｅｒ Ｍｅｔｈｏｄｓ＂ａｎｄ Ｐｒｏｇｒａｍｓ ｉｎ Ｂｉｏｍｅｄｉ— ｃｉｎｅ，１９９５．４７（２）：１３１ １４６．
