彭 寬, 高新波, 趙 恒, 王曉蕊,屈曉超, 陳多芳, 陳雪利, 梁繼民

1.西安電子科技大學電子工程學院，西安710071；

2.西安電子科技大學生命科學與技術學院，西安710071；

3.西安電子科技大學技術物理學院，西安7100071

基于點加權最小二乘無網格法的光學成像光傳輸模型求解

彭 寬1,2, 高新波1, 趙 恒2, 王曉蕊3,屈曉超2, 陳多芳2, 陳雪利1,2, 梁繼民2

1.西安電子科技大學電子工程學院，西安710071；

2.西安電子科技大學生命科學與技術學院，西安710071；

3.西安電子科技大學技術物理學院，西安7100071

漫射方程是目前光學成像中應用最廣泛的光傳輸模型，通常采用有限元方法進行求解。但是，有限元方法依賴于對整個求解域的網格剖分，而對于類似生物組織的非規(guī)則形狀求解域的網格剖分是非常困難的，甚至昂貴的商業(yè)軟件也難以很好地處理。本文將點加權最小二乘無網格法應用于漫射方程的求解。這種方法只需要在求解域內布置一系列規(guī)則分布的配點即可進行數值求解，從而可以完全避免有限元方法中困難的剖分工作。此外，這種方法通過最小化控制方程和邊界條件在每個配點上產生的殘量加權平方和，建立光源功率和光流率密度之間的關系，不需要進行任何數值積分運算，非常適合應用于非規(guī)則求解域的求解?；跀底质竽Ｐ偷臄抵捣抡鎸嶒烌炞C了該方法的準確性和有效性。

無網格法；漫射方程；光傳輸模型；光學成像

引 言

作為一種新興的成像技術，分子影像在近10年來獲得了飛速發(fā)展，并已經在疾病診斷、腫瘤檢測和藥物研發(fā)中獲得了廣泛應用[1]。它可以從分子和細胞水平上實現(xiàn)對生物有機體生理變化的實時、無創(chuàng)、動態(tài)的成像，為特定基因功能、生物體生長發(fā)育、疾病發(fā)生發(fā)展，以及藥物作用效果和動力學變化等研究提供了獲取信息的有效手段[1]。其中光學分子成像技術由于其在成像靈敏度、性價比和安全性等方面的優(yōu)勢，已經成為分子影像領域新的研究熱點[2～5]。

光學成像的核心問題之一是如何對光在生物組織中的傳輸情況進行準確建模和快速求解。輻射傳輸方程 (radiation transfer equation，RTE)可以準確地描述光在生物組織中的傳輸過程，但它的求解難度和求解計算量非常大，因而目前無法應用于實際問題[6,7]。RTE的一種低階近似——漫射方程(diffusion equation，DE)的求解計算量遠遠小于RTE，在組織光學特性參數滿足漫射條件，即約化散射系數遠大于吸收系數時，具有可接受的計算精度。由于很多生物組織的光學特性參數滿足漫射條件，因此DE成為了目前光學成像中應用最為廣泛的一種光傳輸模型[6,7]。

考慮到光學成像的實際需要，對DE的求解大多采用數值方法。其中有限元方法 (finite element method，F(xiàn)EM)是最為常用的數值求解方法[6～8]。有限元方法嚴重依賴于對整個求解域的網格剖分，但對于光學成像中常見的非規(guī)則形狀求解域，網格剖分是一個極為復雜和困難的過程[6]，通常需要借助專業(yè)的商業(yè)軟件，在某些極端情況下，商業(yè)軟件甚至都會出現(xiàn)剖分失敗的現(xiàn)象。本文提出了一種基于點加權最小二乘無網格法 (point weighted least-squares，PWLS)[9]的DE數值求解方法。通過移動最小二乘方法為每個配點構造近似函數，用在求解域內規(guī)則布置的配點代替有限方法中的網格單元，從而徹底避免了對求解域內的網格剖分問題。此外，本文提出的方法通過最小化控制方程和邊界條件在所有配點上產生的殘量加權平方和，構造生物組織中光源功率和光流率密度之間的關系，不需要進行任何數值積分，因而能夠很容易地處理非規(guī)則形狀的求解域。

理論與方法

漫射方程及其邊界條件

如果生物組織的光學特性參數滿足漫射條件，則穩(wěn)定光源發(fā)出的光在其中的傳輸情況可以用如下穩(wěn)態(tài)漫射方程表示[10]：

其中，μa為吸收系數，D=1/(3μa＋3μs′)為漫射系數，μs′為約化散射系數，Φ 為光流率密度，S為光源功率密度。

如果光學成像在全黑的環(huán)境中進行，即沒有環(huán)境光從外部進入生物組織，如下Robin邊界條件可以用于描述光在生物組織邊界的邊界效應[11]：

其中，An為生物組織與空氣之間的折射率失配系數，V為生物組織邊界單位外法向量。

移動最小二乘近似函數的構造

移動最小二乘近似是一種常見的無網格近似函數構造方法，它根據配點周圍采樣點上的二次誤差和最小的原則，為每個配點構造近似函數[12]。在這種方法中，首先通過如下多項式基函數的線性組合來表示配點i上光流率密度Φi的近似函數：

其中，Pig為基函數向量在第i個配點的第g個采樣點上的值；G為一個配點周圍采樣點的個數；wig是一個非負的權值函數，由配點-采樣點間距與采樣點的影響域半徑的比值r決定。本文采用了如下的權值函數：

在本文中，離某個配點最近的27個配點構成了它的采樣點，采樣點的影響域半徑定義為離采樣點最近的27個配點中最遠的采樣點-配點間距的1.8倍。通過?J/?=0最小化J，可以獲得如下的線性方程組：

其中，為配點周圍采樣點上的光流率密度向量，ai為待求的系數列向量。通過式(6)可求解出待求系數向量ai=(Ai)－1BiΦi，將其帶入式(3)，即可獲得某個配點上的近似函數：

其中，Ni=P(Ai)－1Bi為移動最小二乘近似的形函數。

點加權最小二乘無網格法

將求解域中配點i和邊界上配點j的移動最小二乘近似函數代入DE和Robin邊界條件，即可獲得其對應的控制方程殘量和邊界條件殘量：

其中，列向量Si為配點i周圍采樣點上的光源功率密度。然后對求解域中所有配點的控制方程殘量與邊界上所有配點的邊界條件殘量進行加權平方求和[9]：

其中，Φ和S分別為所有區(qū)域配點上的光流率密度和光源功率密度。在光學成像的光傳輸模型求解中，S被認為是已知量，通過式(12)即可求解出Φ。

實驗和結果

我們分別用基于數字鼠模型的勻質和非勻質數值仿體測試了本文方法的性能，并通過與FEM方法求解結果的對比來驗證方法的有效性。此外，由于蒙特卡洛 (Monte Carlo，MC)方法在描述生物組織光傳輸方面具有很好的精度[6,13～15]，它的求解結果被作為評價求解精度的基準。均方根誤差e被用來衡量兩組歸一化數據d1和d2之間的誤差：

圖1 基于小鼠胃部幾何形狀構建的勻質仿體及其離散情況 (A)小鼠胃的幾何結構；(B)PWLS方法中基于配點的仿體離散；(C)有限元方法中基于網格的仿體離散Fig.1 The homogeneous phantom based on the mouse stomach and its discretization (A)Geometrical structure of mouse stomach;(B)The phantom discretization based on collocation points for PWLS;(C)The phantom discretization based on mesh for FEM

勻質仿體的求解

在本實驗中，文獻[16]中提供的數字鼠胃部CT數據被用來構造勻質仿體，其形狀如圖1A所示。首先提取區(qū)域的邊界，并通過三角形進行離散。基于這個離散化的邊界，對整個仿體進行四面體網格剖分以生成FEM方法所需的網格。本實驗中，PWLS方法所需的配點通過自己的代碼布設，其中邊界配點直接使用了邊界離散三角形的頂點，內部的配點則以0.5 mm的分辨率均勻布設，對整個仿體的離散一共使用了3283個配點，如圖1B所示。FEM方法所用的網格中包含了4110個節(jié)點和19572個四面體，如圖1C所示。MC方法則不需要對求解域的內部進行離散，直接使用前面提到的三角形離散表面即可，MC方法中使用的光源光子隨機樣本數為108。一個半徑為0.5 mm的球形光源被放置于(25,55,10)mm處。根據文獻[17]，整個仿體被設置為小鼠胃壁的光學參數，如表1最后一行所示。由于MC方法難以獲取仿體內部的光流率密度分布，我們只比較三種方法獲取的仿體邊界光流率密度。PWLS、FEM與MC三種方法獲得的歸一化的邊界光流率密度如圖2所示，可見三種方法獲得了非常相似的結果。為了進一步比較三種方法獲得的結果，我們取出三種方法在位于z=(10±0.3)mm和z=(13±0.3)mm兩個區(qū)域內的邊界離散三角形頂點上所獲結果進行曲線比較，如圖3所示。在圖3A中，PWLS方法與MC方法在兩條曲線上的均方根誤差為1.08%；在圖3B中，F(xiàn)EM方法與MC方法在兩條曲線上的均方根誤差為1.63%。由于PWLS方法中所用的配點個數少于有限元網格中的節(jié)點數，可見PWLS方法在離散精細度較低的情況下反而獲得了比FEM方法稍好的結果。

表1 數字鼠各組織的光線參數Table 1 Optical parameters of the mouse tissue

圖2 PWLS、FEM與MC三種方法在勻質仿體表面獲得的歸一化的光流率密度 (A)PWLS方法的結果；(B)FEM方法的結果，(C)MC方法的結果Fig.2 Normalized photon flux density on the surface of homogeneous phantom (A)Result of PWLS;(B)Result of FEM;(C)Result of MC

圖3 對PWLS、FEM與MC三種方法在勻質仿體表面獲得的歸一化的光流率密度的曲線比較 (A)PWLS方法與MC方法的曲線結果；(B)FEM方法與MC方法的曲線結果Fig.3 Curve comparisonsof the normalized photon flux densityon homogeneous phantom surface between PWLS,FEM and MC results (A)Comparison between PWLS and MC results;(B)Comparison between FEM and MC results

非勻質仿體的求解

在這個實驗中，我們使用自主研制的Micro-CT系統(tǒng) (X-ray tube：series 5000 apogee，OXFORD instruments；Flat plane X-ray sensor：C7921CA-02，HAMAMATSU) 所獲取的小鼠軀干CT數據來構造非勻質仿體。這個CT數據切片數為256，分辨率為512×512，其中一共包含了脂肪、心臟、肺、肝和腎五種不同的生物組織，如圖4所示。由于各種組織

圖4 基于小鼠軀干幾何形狀構建的非勻質仿體及其離散情況 (A)小鼠軀干的幾何結構；(B)PWLS方法中基于配點的仿體離散；(C)有限元方法中基于網格的仿體離散Fig.4 The heterogeneous phantom based on the structure of mouse trunk (A) Tissuegeometrical information of the mouse trunk; (B) The phantom discretization with collocation points;(C) The phantom discretization with mesh

之間CT數據對比度不高，各種組織的分割需要人工干預。根據文獻[17]，設置各個組織光學參數，如表1中1～5行所示。與上個實驗類似，由三角形離散的各個組織的表面首先被提取出來，然后再進行四面體剖分或配點布置。其中，F(xiàn)EM網格中包含6357個節(jié)點和34791個四面體；而PWLS方法則使用了分辨率為2 mm的4673個配點。兩種方法的仿體離散結果如圖4所示。MC方法需要使用各個組織經過三角形離散的表面，所用的光源光子隨機樣本數仍為108。本實驗中所使用的光源為橢球形，三軸半徑分別為1.5、1.5和5 mm，其中心位于(29,23,26)mm。PWLS、FEM與MC三種方法獲得的歸一化的邊界光流率密度如圖5所示。從圖中可見，雖然三種方法獲得了非常相似的數據分布趨勢，但在光斑的中心區(qū)域出現(xiàn)了一定的差別。這是由于小鼠軀干中所包含的肝組織的吸收系數較大，明顯不滿足約化散射系數遠大于吸收系數的漫射條件，從而導致DE方程的準確性下降造成的。與上個實驗類似，我們同樣取出三種方法在位于z=(27±1)mm和z=(35±1)mm兩個區(qū)域內的邊界離散三角形頂點上所獲結果進行曲線比較，如圖6所示。從圖中可以看出，由于前面提到的原因，PWLS和FEM的計算精度相比上一個實驗都有所下降，兩個方法在兩條曲線上與MC方法的均方根誤差分別達到了1.83%和2.26%，但PWLS方法仍然用較低的離散精細度達到了稍好于FEM的計算精度。

圖5 PWLS、FEM與MC三種方法在非勻質仿體表面獲得的歸一化的光流率密度(A)PWLS方法的結果；(B)FEM方法的結果；(C)MC方法的結果Fig.5 Normalized photon flux density on the surface of heterogeneous phantom(A)Result of PWLS;(B)Result of FEM;(C)Result of MC

圖6 對PWLS、FEM與MC三種方法在非勻質仿體表面獲得的歸一化的光流率密度的曲線比較 (A)PWLS方法與MC方法的曲線結果；(B)FEM方法與MC方法的曲線結果Fig.6 Comparisons of the normalized photon flux density on heterogeneous phantom surface between PWLS,FEM and MC results (A)Comparison between PWLS and MC results;(B)Comparison between FEM and MC results

討 論

有限元方法是目前最為常用的光學成像光傳輸模型——DE的求解方法。這種方法依賴于對整個求解域的網格剖分，但對于生物組織這種形狀非規(guī)則區(qū)域，進行網格剖分通常非常困難，即使對商業(yè)軟件而言也是如此。在本文實驗中，在提取了離散的仿體邊界之后，對整個求解域進行網格剖分時常常會出現(xiàn)剖分失敗的現(xiàn)象，需要反復調整仿體邊界離散三角形的數量，才能使網格剖分成功進行，可見對形狀非規(guī)則區(qū)域進行網格剖分的難度。本文提出的基于加權最小二乘無網格法的求解方法，只需布置一系列規(guī)則排列的配點來對求解域進行離散，難度遠遠小于對求解域的網格剖分。此外，本文的方法通過最小化控制方程和邊界條件在每個配點上產生的殘量加權平方和，建立光源功率和光流率密度之間的關系，不需要進行任何形式的數值積分，大大地方便了其在生物組織這種形狀非規(guī)則區(qū)域中的應用。

在實驗中可以發(fā)現(xiàn)，本文的方法可以用較低的求解域離散度，實現(xiàn)與有限元方法相當的求解精度。這是因為無網格法中配點近似函數的構造考慮了周圍27個配點，而有限元網格中每個節(jié)點的近似函數構造只受與其共單元的節(jié)點影響，其數量通常不到20個。因此，每個配點近似函數中可以包含更多的信息量，從而降低離散整個求解域所需的配點數量。由于最后獲得的系統(tǒng)方程組的維數與用于離散求解域的節(jié)點或配點個數相等，在保證計算精度的前提下使用了較少的配點數量，顯然可以降低求解的內存消耗與計算量。

最后，在第二個實驗中，我們發(fā)現(xiàn)DE作為小鼠全身成像光傳輸模型時，也存在一定的局限性：由于小鼠中一部分組織的光學特性參數并不滿足漫射條件，應用DE有可能導致求解精度的下降。這就要求在小鼠成像中使用不受漫射條件限制的RTE的高階近似形式，如簡化球諧波、球諧波和離散坐標等[6,18]，但這些高階近似方程的計算量都較大，如何提高它們的求解速度將是一個非常值得研究的問題。

1. Cherry S.In vivomolecular and genomic imaging:New challenges for imaging physics.Phys Med Biol,2004,49:13～48

2. Weissleder R,NtziachristosV.Shedding light onto live molecular targets.Nat Med,2003,9:124～128

3. Bhaumik S, GambhirSS. Optical imaging of Renilla luciferase reporter gene expression in living mice. Proc Natl Acad Sci USA,2002,99:377～382

4. Massoud TF,Gambhir SS.Molecular imaging in living subjects:Seeing fundamental biological processes in a new light.Genes Dev,2003,17:545～580

5. RiceW,CableMD,NelsonMB.Invivoimagingof light-emitting probes.J Biomed Opt,2001,6:432～440

6. Gibson AP,Hebden JC,Arridge SR.Recent advances in diffuse optical imaging.Phys Med Biol,2005,50:1～43

7.Wang G,Cong W,Li Y,Han W,Kumar D,Qian X,Shen H,Jiang M,Zhou T,Cheng J,Tian J,Lv Y,Li H,Luo J.Recent development in bioluminescence tomography.Curr Imaging Rev,2006,2:453～457

8. Arridge SR,Dehghani H,Schweiger M,Okada E.The finite element modelfor the propagation of light in scattering media: A directmethod fordomains with nonscattering regions.Med Phys,2000,27:252～264

9. Wang Q, LiH, Lam K. Developmentofa new meshless-point weighted least-squares(PWLS)method for computationalmechanics. ComputMech, 2005, 35:170～181

10. Ntziachristos V, Tung C, BremerC, WeisslederR.Fluorescence molecular tomography resolves protease activityin vivo.Nat Med,2002,8:757～760

11.Schweiger M,Arridge SR,Hiraoka M,Delpy DT.The finite element method for the propagation of light in scattering media:Boundary and source conditions.Med Phys,1995,22:1779～1792

12.Belytschko T,Krongauz Y,Organ D.Meshless methods:An overview andrecent developments. Comput Method Appl Mech Eng,1996,139:3～47

13.WangL,JacquesSL,ZhengL.MCML-MonteCarlo modeling ofphoton transportin multi-layered tissues.Comput Meth Prog Biomed,1995,47:131～146

14.Boas D,Culver J,Stott J,Dunn A.Three dimensional MonteCarlo code for photon migration through complex heterogeneous media including the adult human head.Opt Express,2002,10:159～169

15.Ren NN,Liang JM,Qu XC,Li JF,Lu BJ,Tian J.GPU-based Monte Carlo simulation for light propagation in complex heterogeneous tissues.Opt Express,2010,18:6811～6823

16. Dogdas B, StoutD, Chatziioannou A, Leahy RM.Digimouse:A 3D whole body mouse atlas from CT and cryosection data.Phys Med Biol,2007,52:577～587

17. Alexandrakis G, Rannou FR, Chatziioannou AF.Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: A computer simulation feasibility study. Phys Med Biol, 2005, 50:4225～4241

18.Klose AD,Larsen EW.Light transport in biological tissue based on the simplified spherical harmonics equations.J Comput Phys,2006,220:441～470

Point Weighted Least-Squares Meshless Method for Photon Transport in Complex Biological Tissues

PENG Kuan1,2,GAO Xinbo1,ZHAO Heng2,WANG Xiaorui3,QU Xiaochao2,CHEN Duofang2,CHEN Xueli1,2,LIANG Jimin2
1.School of Electronic Engineering,Xidian University,Xi'an,710071,China;
2.Life Sciences Research Center,School of Life Sciences and Technology,Xidian University,Xi'an 710071,China;
3.School of Technology Physics,Xidian University,Xi'an 710071,China

This work was supported by grants from the"973"Program(2011CB707702),the National Natural Science Foundation of China(81090272,81000632,30900334,60771068),the Shaanxi Provincial Natural Science Foundation Research Project(2009JQ8018),and the Fundamental Research Funds for the Central Universities

Nov 1,2010 Accepted:Jan 18,2011

LIANG Jimin,Tel:+86(29)81891060,E-mail:jimleung@mail.xidian.edu.cn

Diffusion equation is a widely used model for photon transport in biological tissues.Commonly,it can be solved by the finite element method,which depends on the mesh discretization of interesting area.However,this kind of discretization is quite difficult for complex geometries such as biological tissues,even using commercial software.In this paper,a point weighted least-squares(PWLS)meshless method for the photon transport in diffusive biological tissueswith complex and irregular geometriesis presented. The proposed method doesthe calculation with some collocation pointsdistributed regularly in the interesting area,so that complicated mesh generation required by finite element method can be avoided.Moreover,our method establishes the relationship between the source distribution and the photon flux density by minimizing the weighted residuals quadratic sum of all collocation points. Thus, it doesn't need the complicated integration calculation,which provides more convenience for the application in irregular shaped region like biological tissues.Numerical simulations with the phantoms constructed based on mouse tissues demonstrate the accuracy and effectiveness of the proposed method.

Meshless method;Diffusion equation;Light transport model;Optical imaging

2010-11-01；接受日期：2011-01-18

“973”計劃項目(2011CB707702)，中國科學院百人計劃，國家自然科學基金項目(81090272，81000632，30900334，60771068)，陜西省自然科學基礎研究計劃項目(2009JQ8018)，中央高?；究蒲袠I(yè)務費專項資金資助

梁繼民，電話：(029)81891060，E-mail：jimleung@mail.xidian.edu.cn

TP-391

10.3724/SP.J.1260.2011.00373