本發明涉及一種利用生物腔體組織產生的光聲信號重建組織光吸收系數空間分布的方法，屬于醫學成像技術領域。

背景技術：

光聲層析成像是一種新型的非電離式生物醫學功能成像方法，它結合了超聲成像的高分辨率和光學成像的高對比度的優勢。pat成像方法的逆問題包括聲學和光學兩方面：聲學逆問題是指根據超聲探測器采集到的光聲信號(本質是超聲波)重建組織內部的初始聲壓分布或空間光吸收能量密度；光學逆問題是指根據超聲探測器采集到的光聲信號或者據此重建出的光吸收能量密度，估算組織光學特性參數的空間分布，得到對組織光學特性的定量評價。光吸收能量密度是由局部的光吸收系數和光子數分布共同決定的，并不能反映組織的光學特性，而光吸收系數與組織的化學成分密切相關，正常和病變組織的光學特性參數之間通常有較明顯的差異，因此組織光吸收系數的空間分布圖可為疾病的早期診斷提供更加準確可靠的基礎參考信息。但到目前為止，組織光吸收系數的重建方法仍不成熟，還有必要進一步進行研究。

技術實現要素：

本發明的目的在于針對現有技術之弊端，提供一種用于生物光聲內窺成像的光吸收系數重建方法，以獲取不同組織的光吸收系數之間的差異，為疾病的早期診斷提供機體組織變化的準確可靠的參考信息。

本發明所述問題是以下述技術方案解決的：

一種用于生物光聲內窺成像的光吸收系數重建方法，所述方法首先根據超聲探測器從各個角度采集的生物腔體組織產生的光聲信號重建出光吸收能量分布的測量值；然后通過前向仿真獲得光吸收能量分布的理論值；最后根據光吸收能量分布的測量值和理論值構造目標函數，并通過對目標函數進行優化，重建出腔體橫截面上組織光吸收系數的空間分布。

上述用于生物光聲內窺成像的光吸收系數重建方法，所述重建按以下步驟進行：

a.重建光吸收能量分布的測量值

在直線掃描模式下，根據超聲探測器從各個測量角度采集到的周圍組織產生的光聲信號，得到光吸收能量分布的測量值：

其中，r是θ-l平面極坐標系中的一點，c是光在組織內的傳播速度，hm(r)是位置r處的光吸收能量分布的測量值，z0是位置r處與組織表面θ軸之間的垂直距離，ri是與角度θi(i＝1,2,...,m)相對應的位置，cs是超聲信號在組織中傳播速度，β是組織的等壓膨脹系數，cp是組織的比熱容，pi(r,t)是超聲探測器在時刻t、角度θi(i＝1,2,...,m)、位置r處采集到的組織產生的光聲信號的聲壓；

b.通過前向仿真獲得光吸收能量分布的理論值

首先，用由節點表示的面積元對成像區域進行離散化，并根據成像區域內每個節點處光擴散方程的差分形式，計算出位置r處的光子密度函數φ(r)；

然后根據位置r處的光子密度函數φ(r)得到光吸收能量密度的理論值h(r,μa(r))：

h(r,μa(r))＝μa(r)·h·f·φ(r)

其中，μa(r)是位置r處組織的光吸收系數，h是普朗克常數，f是入射光的頻率；

c.重建光吸收系數的空間分布

假定組織的光散射系數已知，利用下式求得光吸收系數的空間分布：

其中，c是與光聲信號采集系統有關的標定因子，是估算出的、待測區域內位置r處的光吸收系數。

本發明能夠利用生物腔體組織產生的光聲信號完整地重建出腔體橫截面上組織光吸收系數的空間分布，為疾病的早期診斷提供機體組織變化的準確可靠的參考信息。

附圖說明

下面結合附圖對本發明作進一步說明。

圖1是含有脂質斑塊的血管橫截面示意圖；

圖2是θ-l極坐標平面內網格劃分的內點示意圖；

圖3是θ-l極坐標平面內網格劃分的邊界點示意圖；

圖4是腔體成像區域內面向光源的邊界點示意圖。

文中各符號為：θ、極角；l、極徑；θ-l、平面極坐標系；m、血管橫截面被等角度分割的總份數；θi、成像導管的第i個測量角度，其中i＝1,2,...,m；r、θ-l平面極坐標系中的一點；hm(r)、位置r處的光吸收能量分布的測量值；z0、位置r處與組織表面θ軸之間的垂直距離；ri、θ-l平面中與角度θi(i＝1,2,...,m)相對應的位置；cs、超聲信號在生物軟組織中的傳播速度；β、組織的等壓膨脹系數；cp、組織的比熱容；pi(r,t)、超聲探測器在時刻t、角度θi、位置r處采集到的組織產生的光聲信號的聲壓，其中i＝1,2,...,m；哈密頓算子；φ(r)、位置r處的光子密度函數；φ(r,rs,q)、時刻t位置r處的時域光子密度函數的laplace變換；δ(r-rs)、位置r處無向點源的laplace變換；q、laplace變換的復頻率因子，當q＝0時，de是穩態擴散方程；ω、待成像的組織區域；c、光在組織內的傳播速度；組織表面點處的外法向量；qs、光源強度；rs、θ-l平面極坐標系中光源所在位置；rf、擴散傳輸內反射系數；n、組織相對于環境的相對折射率；μa(r)、位置r處組織的光吸收系數；μs(r)、位置r處組織的光散射系數；μ′s(r)、位置r處組織的約化散射系數；g、組織的各向異性因子；d、ω中一個充分光滑的區域；hθ、沿θ軸的正向步長；hl、沿l軸的正向步長；p、區域d中的一點；(ihθ,jhl)、網格劃分后θ-l坐標系中點p的坐標；φi,j、節點(ihθ,jhl)處的光子密度函數值；ki+1,j、點((i+1)hθ,jhl)處的k值；ki,j+1、點(ihθ,(j+1)hl)處的k值；ki-1,j、點((i-1)hθ,jhl)處的k值；ki,j-1、點(ihθ,(j-1)hl)處的k值；|ab|、線段ab的長度；|be|、線段be的長度；|ea|、線段ea的長度；|d|、區域d的面積；h(r,μa(r))、位置r處光吸收能量密度的理論值；h、普朗克常數；f、入射光的頻率；c、與光聲信號采集系統有關的標定因子；待測區域內位置r處估算出的光吸收系數；x、待求的μa(r)；f(x)、待優化的目標函數；位置r處光吸收系數的估計值；λ、正則化參數；x0、p0、光吸收系數和次梯度的初始值；xn+1、pn+1、n次迭代后的光吸收系數和次梯度；xn、pn、n-1次迭代后的光吸收系數和次梯度；<pn,x>、pn與x的內積；m、ne×n維的稀疏矩陣；ne、迭代次數；n、網格個數；|·|、l1準則；α、懲罰參數；ν、第n次迭代產生的殘余量，其初始值v0＝0；shrink(·)、shrink算子；wn、第n次搜索的方向；b0、近似逆hessian矩陣的初始值；i、單位矩陣；bn、第n次迭代的近似逆hessian矩陣；an、步長；g(x)在xn的梯度。

具體實施方式

本發明的處理步驟包括：

1.由超聲探測器接收到的光聲信號重建光吸收能量分布的測量值

如附圖1所示，以血管內光聲成像為例，成像導管(超聲探測器位于成像導管頂端，用于接收周圍組織產生的光聲信號)位于血管橫截面的中心，周圍由內向外依次為管腔、粥樣硬化斑塊、血管壁內膜/中膜和外膜。忽略超聲探測器的孔徑效應，將其看作理想的點換能器，其掃描軌跡為平行于成像平面的圓形軌跡。

血管橫截面所在的坐標系是θ-l極坐標系，其中θ是極角，l是極徑，坐標原點是成像導管中心，水平向右的方向為θ軸正方向。以坐標原點為中心，將血管橫截面等角度分成m個扇區，成像導管的第i個測量角度是θi＝360(i-1)/m，其中i＝1,2,...,m。

在直線掃描模式下，根據超聲探測器采集到的光聲信號，得到光吸收能量分布的測量值：

其中，r是θ-l平面極坐標系中的一點，c是光在組織內的傳播速度，hm(r)是位置r處的光吸收能量分布的測量值，z0是位置r處與組織表面θ軸之間的垂直距離，ri是與角度θi(i＝1,2,...,m)相對應的位置，cs是超聲信號在生物軟組織中的傳播速度，β是組織的等壓膨脹系數，cp是組織的比熱容，pi(r,t)是超聲探測器在時刻t、角度θi(i＝1,2,...,m)、位置r處采集到的組織產生的光聲信號的聲壓。

2.通過前向仿真獲得光吸收能量分布的理論值

首先，用由節點表示的面積元對成像區域進行離散化，并根據成像區域內每個節點處光擴散方程的差分形式，計算出位置r處的光子密度函數φ(r)。具體步驟如下：

采用準直光源模型，將光源視為組織邊界處的內向光子流，進而嵌入到robin邊界條件中，則邊界含有光源項的de復頻域表達式為：

其中，是哈密頓算子；r是θ-l坐標系中的一點；φ(r,rs,q)是時刻t位置r處的時域光子密度函數的laplace變換；δ(r-rs)是位置r處無向點源的laplace變換；q是laplace變換的復頻率因子，當q＝0時，de是穩態擴散方程；ω是待成像的組織區域；c是光在組織內的傳播速度；是組織表面點處的外法向量；qs是光源強度；rs是θ-l平面極坐標系中光源所在位置；rf是擴散傳輸內反射系數：

rf≈-1.4399n^-2+0.7099n^-1+0.6681+0.0636n(3)

式(3)中，n是組織相對于環境的相對折射率；μa(r)和μs(r)分別是位置r處組織的光吸收系數和散射系數；當μs'＞＞μa時

其中，μs′(r)是位置r處組織的約化散射系數：

μ′s(r)＝μs(r)(1-g)(5)

其中g是組織的各向異性因子。

在ω中取一個充分光滑的區域d∈ω，在θ-l平面內，對區域d進行矩形網格劃分，沿θ軸和l軸的正向步長分別為hθ和hl，d中的一點p在θ-l坐標系中的坐標是(ihθ,jhl)，在區域d上，對式(2)進行積分，得到：

如附圖2所示，陰影部分為區域d，若p是內點，則式(6)的差分形式為：

其中，φi,j是節點(ihθ,jhl)處的光子密度函數值；根據式(4)求出ki,j，其中，ki+1,j是在點((i+1)hθ,jhl)處的k值，ki,j+1是在點(ihθ,(j+1)hl)處的k值，ki-1,j是在點((i-1)hθ,jhl)處的k值，ki,j-1是在點(ihθ,(j-1)hl)處的k值。

如附圖3所示，由弧線線段ab和be組成區域d，若p是邊界點，則式(6)的差分形式為：

式中，|ab|、|be|和|ea|分別是線段ab、be和ea的長度；|d|是區域d的面積。

如附圖4所示，光源沿l軸即徑向入射，由弧線線段ab和be組成區域d，若p是面向光源的邊界點，則式(6)的差分形式為：

然后，根據位置r處的光子密度函數φ(r)得到光吸收能量密度的理論值h(r,μa(r))：

h(r,μa(r))＝μa(r)·h·f·φ(r)(10)

其中，h是普朗克常數，f是入射光的頻率。

3.重建光吸收系數的空間分布

假定組織的光散射系數已知，求解光吸收系數分布的問題表述為如下的非線性最小二乘問題：

其中，c是與光聲信號采集系統有關的標定因子，是估算出的、待測區域內位置r處的光吸收系數。為表述方便，用x表示μa(r)，并定義：

f(x)＝||hm(r)-c·h(r,x)||²(12)

該最優化問題是一個病態問題，本發明方法采用tv正則化(gaoh,zhaoh.multilevelbioluminescencetomographybasedonradiativetransferequationpart1:l1regularization.opticsexpress,2010,18(3):1854-1871.)使其良態化，將待優化的目標函數改寫為：

其中，是位置r處光吸收系數的估計值；λ是正則化參數，在迭代過程中設定λ＝1。

本發明采用bregman方法(gaoh,zhaoh,oshers.bregmanmethodsinquantitativephotoacoustictomography.camreport,2010,30(6):3043-3054)迭代求解式(13)的最優解。具體步驟如下：

設定各參數的初始值：光吸收系數x0＝0，次梯度p0＝0。第n次迭代后，光吸收系數xn+1與次梯度pn+1的結果為：

其中，<pn,x>是次梯度pn與x的內積。

在離散網格條件下，式(14)右端的tv項簡化為

||x||tv＝|mx|(16)

其中，m是ne×n維的稀疏矩陣，ne是迭代次數，n是網格個數，|·|是l1準則。

令

dn＝mxn(17)

將式(14)轉化為無約束問題：

其中，α是懲罰參數，取值為1。對式(18)進行n次迭代后的結果為：

νn+1＝νn+dn+1-mxn+1(21)

其中，vn+1是第n次迭代產生的殘余量，其初始值ν0＝0；shrink算子的定義為：

設

求解xn+1＝argmin[g(x)]最小化問題的具體步驟如下：

步驟1：初始化，設x0＝0，近似逆hessian矩陣的初始值b0＝i；

步驟2：計算搜索方向：其中bn是第n次迭代的近似逆hessian矩陣；

步驟3：沿搜索方向找到下一個迭代點：

其中，an是步長；是g(x)在xn的梯度：

的更新方法是：

步驟4：根據判斷迭代是否停止；

步驟5：更新搜索方向wn。