本發(fā)明涉及一種利用生物腔體組織產(chǎn)生的光聲信號(hào)重建組織光吸收系數(shù)空間分布的方法，屬于醫(yī)學(xué)成像技術(shù)領(lǐng)域。

背景技術(shù)：

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

技術(shù)實(shí)現(xiàn)要素：

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

本發(fā)明所述問題是以下述技術(shù)方案解決的：

一種用于生物光聲內(nèi)窺成像的光吸收系數(shù)重建方法，所述方法首先根據(jù)超聲探測(cè)器從各個(gè)角度采集的生物腔體組織產(chǎn)生的光聲信號(hào)重建出光吸收能量分布的測(cè)量值；然后通過前向仿真獲得光吸收能量分布的理論值；最后根據(jù)光吸收能量分布的測(cè)量值和理論值構(gòu)造目標(biāo)函數(shù)，并通過對(duì)目標(biāo)函數(shù)進(jìn)行優(yōu)化，重建出腔體橫截面上組織光吸收系數(shù)的空間分布。

上述用于生物光聲內(nèi)窺成像的光吸收系數(shù)重建方法，所述重建按以下步驟進(jìn)行：

a.重建光吸收能量分布的測(cè)量值

在直線掃描模式下，根據(jù)超聲探測(cè)器從各個(gè)測(cè)量角度采集到的周圍組織產(chǎn)生的光聲信號(hào)，得到光吸收能量分布的測(cè)量值：

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

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

首先，用由節(jié)點(diǎn)表示的面積元對(duì)成像區(qū)域進(jìn)行離散化，并根據(jù)成像區(qū)域內(nèi)每個(gè)節(jié)點(diǎn)處光擴(kuò)散方程的差分形式，計(jì)算出位置r處的光子密度函數(shù)φ(r)；

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

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

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

c.重建光吸收系數(shù)的空間分布

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

其中，c是與光聲信號(hào)采集系統(tǒng)有關(guān)的標(biāo)定因子，是估算出的、待測(cè)區(qū)域內(nèi)位置r處的光吸收系數(shù)。

本發(fā)明能夠利用生物腔體組織產(chǎn)生的光聲信號(hào)完整地重建出腔體橫截面上組織光吸收系數(shù)的空間分布，為疾病的早期診斷提供機(jī)體組織變化的準(zhǔn)確可靠的參考信息。

附圖說明

下面結(jié)合附圖對(duì)本發(fā)明作進(jìn)一步說明。

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

圖2是θ-l極坐標(biāo)平面內(nèi)網(wǎng)格劃分的內(nèi)點(diǎn)示意圖；

圖3是θ-l極坐標(biāo)平面內(nèi)網(wǎng)格劃分的邊界點(diǎn)示意圖；

圖4是腔體成像區(qū)域內(nèi)面向光源的邊界點(diǎn)示意圖。

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

具體實(shí)施方式

本發(fā)明的處理步驟包括：

1.由超聲探測(cè)器接收到的光聲信號(hào)重建光吸收能量分布的測(cè)量值

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

血管橫截面所在的坐標(biāo)系是θ-l極坐標(biāo)系，其中θ是極角，l是極徑，坐標(biāo)原點(diǎn)是成像導(dǎo)管中心，水平向右的方向?yàn)棣容S正方向。以坐標(biāo)原點(diǎn)為中心，將血管橫截面等角度分成m個(gè)扇區(qū)，成像導(dǎo)管的第i個(gè)測(cè)量角度是θi＝360(i-1)/m，其中i＝1,2,...,m。

在直線掃描模式下，根據(jù)超聲探測(cè)器采集到的光聲信號(hào)，得到光吸收能量分布的測(cè)量值：

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

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

首先，用由節(jié)點(diǎn)表示的面積元對(duì)成像區(qū)域進(jìn)行離散化，并根據(jù)成像區(qū)域內(nèi)每個(gè)節(jié)點(diǎn)處光擴(kuò)散方程的差分形式，計(jì)算出位置r處的光子密度函數(shù)φ(r)。具體步驟如下：

采用準(zhǔn)直光源模型，將光源視為組織邊界處的內(nèi)向光子流，進(jìn)而嵌入到robin邊界條件中，則邊界含有光源項(xiàng)的de復(fù)頻域表達(dá)式為：

其中，是哈密頓算子；r是θ-l坐標(biāo)系中的一點(diǎn)；φ(r,rs,q)是時(shí)刻t位置r處的時(shí)域光子密度函數(shù)的laplace變換；δ(r-rs)是位置r處無向點(diǎn)源的laplace變換；q是laplace變換的復(fù)頻率因子，當(dāng)q＝0時(shí)，de是穩(wěn)態(tài)擴(kuò)散方程；ω是待成像的組織區(qū)域；c是光在組織內(nèi)的傳播速度；是組織表面點(diǎn)處的外法向量；qs是光源強(qiáng)度；rs是θ-l平面極坐標(biāo)系中光源所在位置；rf是擴(kuò)散傳輸內(nèi)反射系數(shù)：

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

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

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

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

其中g(shù)是組織的各向異性因子。

在ω中取一個(gè)充分光滑的區(qū)域d∈ω，在θ-l平面內(nèi)，對(duì)區(qū)域d進(jìn)行矩形網(wǎng)格劃分，沿θ軸和l軸的正向步長(zhǎng)分別為hθ和hl，d中的一點(diǎn)p在θ-l坐標(biāo)系中的坐標(biāo)是(ihθ,jhl)，在區(qū)域d上，對(duì)式(2)進(jìn)行積分，得到：

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

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

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

式中，|ab|、|be|和|ea|分別是線段ab、be和ea的長(zhǎng)度；|d|是區(qū)域d的面積。

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

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

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

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

3.重建光吸收系數(shù)的空間分布

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

其中，c是與光聲信號(hào)采集系統(tǒng)有關(guān)的標(biāo)定因子，是估算出的、待測(cè)區(qū)域內(nèi)位置r處的光吸收系數(shù)。為表述方便，用x表示μa(r)，并定義：

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

該最優(yōu)化問題是一個(gè)病態(tài)問題，本發(fā)明方法采用tv正則化(gaoh,zhaoh.multilevelbioluminescencetomographybasedonradiativetransferequationpart1:l1regularization.opticsexpress,2010,18(3):1854-1871.)使其良態(tài)化，將待優(yōu)化的目標(biāo)函數(shù)改寫為：

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

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

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

其中，<pn,x>是次梯度pn與x的內(nèi)積。

在離散網(wǎng)格條件下，式(14)右端的tv項(xiàng)簡(jiǎn)化為

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

其中，m是ne×n維的稀疏矩陣，ne是迭代次數(shù)，n是網(wǎng)格個(gè)數(shù)，|·|是l1準(zhǔn)則。

令

dn＝mxn(17)

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

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

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

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

設(shè)

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

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

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

步驟3：沿搜索方向找到下一個(gè)迭代點(diǎn)：

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

的更新方法是：

步驟4：根據(jù)判斷迭代是否停止；

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