本發(fā)明涉及一種電阻率圖像重建方法，特別涉及一種注入電流式電阻率圖像重建方法。

背景技術(shù)：

由于激勵(lì)頻率的限制，傳統(tǒng)電阻抗成像技術(shù)的靈敏度和空間分辨率不高。單一場(chǎng)都有其物理局限性，多物理場(chǎng)成像由一種物理場(chǎng)提供分辨率，另一種物理場(chǎng)提供對(duì)比度，實(shí)現(xiàn)對(duì)比度和分辨率的同時(shí)提高。電磁場(chǎng)和超聲相結(jié)合的多物理場(chǎng)成像技術(shù)正是考慮到電磁場(chǎng)對(duì)人體組織電阻率的高對(duì)比度和超聲波探測(cè)的高分辨率特性，成為人們研究的熱點(diǎn)，磁熱聲成像正是一種新興的多物理場(chǎng)成像技術(shù)。

磁熱聲成像是由新加坡南洋理工大學(xué)的feng在2013年首次提出的新型的電阻抗成像方法，其原理為：通過(guò)對(duì)導(dǎo)電成像體施加mhz量級(jí)的交變磁場(chǎng)，在目標(biāo)體內(nèi)部產(chǎn)生感應(yīng)電場(chǎng)，進(jìn)而產(chǎn)生焦耳熱，激發(fā)熱彈性的超聲信號(hào)，檢測(cè)超聲信號(hào)進(jìn)行成像。與微波熱聲成像相比，允許更低的功率進(jìn)行高效的成像，并且具有便攜式成像的潛力，同時(shí)，激勵(lì)源的頻率降低，使得磁場(chǎng)穿透組織更深，也避免的輻射。

作為一種新型的多物理場(chǎng)成像方法，2013年feng利用金屬銅仿體，檢測(cè)到磁熱聲信號(hào)，并得到銅仿體的熱聲圖像，并未進(jìn)行電阻率圖像重建，且生物組織不同于金屬銅仿體，磁場(chǎng)和感應(yīng)電流作用產(chǎn)生的洛倫茲力較弱。專(zhuān)利“一種磁熱聲成像的電阻率重建方法”(201410771496.2)公布了一種基于磁熱聲效應(yīng)的電阻率重建方法，其在熱函數(shù)基礎(chǔ)上提出了電阻率重建方法，依然為采用線圈激勵(lì)方式。同時(shí)，外加激勵(lì)線圈產(chǎn)生時(shí)變磁場(chǎng)在目標(biāo)體中產(chǎn)生二次磁場(chǎng)和感應(yīng)電流，磁場(chǎng)和感應(yīng)電流作用同時(shí)產(chǎn)生洛倫茲力和焦耳熱，即磁聲效應(yīng)和磁熱聲效應(yīng)共存，如何區(qū)分磁聲效應(yīng)和磁熱聲效應(yīng)是仍需解決的問(wèn)題。

基于此，采用注入電流式熱聲成像方法，可以避開(kāi)磁聲效應(yīng)和磁熱聲效應(yīng)共存問(wèn)題，注入電流式熱聲成像與磁熱聲成像在激勵(lì)方式、從熱函數(shù)到電阻率的重建方法上均不相同。

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

本發(fā)明的目的是克服現(xiàn)有磁熱聲成像存在的問(wèn)題，提出一種基于最優(yōu)化迭代方法的注入電流式熱聲電阻率圖像重建方法。本發(fā)明避免了熱聲成像過(guò)程中磁聲效應(yīng)的干擾，同時(shí)采用注入電流式激勵(lì)，可增強(qiáng)熱聲效應(yīng)，也實(shí)現(xiàn)了目標(biāo)體電阻率的圖像重建。

注入電流式熱聲成像原理為：通過(guò)注入電極向成像目標(biāo)體注入電流，在成像目標(biāo)體中產(chǎn)生焦耳熱，引起熱膨脹，產(chǎn)生超聲信號(hào)，采用超聲換能器進(jìn)行檢測(cè)，根據(jù)檢測(cè)的超聲信號(hào)，重建熱聲源分布和電阻率。

本發(fā)明基于最優(yōu)化迭代方法的注入電流式熱聲電阻率圖像重建方法具體步驟如下：

第一步為獲取注入電流式熱聲信號(hào)：利用超聲換能器獲取目標(biāo)體各斷層熱聲信號(hào)，每個(gè)斷層掃描一周進(jìn)行熱聲信號(hào)的檢測(cè)，即檢測(cè)超聲信號(hào)；第二步根據(jù)注入電流式熱聲聲壓波動(dòng)方程，獲取目標(biāo)體每個(gè)斷層上的熱聲源分布，利用插值法獲取目標(biāo)體整體的熱聲源分布；第三步對(duì)導(dǎo)電物體的電阻率進(jìn)行空間離散，設(shè)定電阻率的初值，根據(jù)歐姆定律，利用線性有限元方程重建的到目標(biāo)體的矢量電位；第四步根據(jù)熱聲源和電流密度關(guān)系，結(jié)合電流連續(xù)性定理，獲取電阻率所滿(mǎn)足的關(guān)于熱聲源與矢量電位的方程，定義目標(biāo)函數(shù)；第五步將求解的矢量電位代入目標(biāo)函數(shù)，利用目標(biāo)函數(shù)最小化原則，求解導(dǎo)電物體的電阻率。具體分述如下：

第一步：獲取導(dǎo)電物體熱聲信號(hào)

激勵(lì)源通過(guò)注入電極a和注入電極b，對(duì)目標(biāo)體注入脈沖電流，目標(biāo)體在電流作用下產(chǎn)生焦耳熱，進(jìn)而產(chǎn)生熱膨脹，激發(fā)超聲信號(hào)，超聲信號(hào)通過(guò)耦合劑耦合到超聲換能器，超聲換能器接收到信號(hào)后，通過(guò)檢測(cè)系統(tǒng)進(jìn)行信號(hào)的放大、濾波、采集和存儲(chǔ)，超聲換能器在控制器控制下對(duì)目標(biāo)體進(jìn)行掃描檢測(cè)；

第二步：求解目標(biāo)體熱聲源分布

已知熱聲成像的聲壓波動(dòng)方程：

其中r為超聲換能器位置，p(r,t)是聲壓，cs為介質(zhì)中的聲速，cp為目標(biāo)體(3)的比熱容，β為目標(biāo)體(3)的熱膨脹系數(shù)，δ(t)是狄拉克函數(shù)，s(r)為熱聲源分布，t為時(shí)間，為拉普拉斯算符；

熱聲源分布的時(shí)間反演法重建公式為：

其中r為標(biāo)量，r＝|r′-r|，r為矢量，er為單位矢量，r'是超聲探頭的位置，r為熱聲源位置，sd是超聲探頭所在的平面，p′是聲壓對(duì)時(shí)間的一階導(dǎo)數(shù)，n是r′位置sd的法線矢量，β為目標(biāo)體的熱膨脹系數(shù)；

選取目標(biāo)體的某一斷層面z＝z1，超聲換能器在此斷層面上進(jìn)行圓周掃描，采集超聲信號(hào)，利用方程(1)、(2)求解z＝z1斷層面上的熱聲源分布s(x,y,z1)。移動(dòng)超聲換能器在不同斷層面上進(jìn)行掃描檢測(cè)，求解不同斷層面上的熱聲源分布，目標(biāo)體(3)上整體熱聲源分布可由分層檢測(cè)計(jì)算或z方向上插值得到；

第三步：獲取矢量電位

熱聲源s同時(shí)是電阻率和電流密度的函數(shù)，熱聲源可以表示為：

s＝ρj²＝ρj·j(3)

其中，ρ為目標(biāo)體(3)的電阻率，j為目標(biāo)體(3)內(nèi)電流密度分布；

根據(jù)電流連續(xù)性定理，電流的散度為零，引入矢量電位，則有：

其中，t為矢量電位，為矢量電位的旋度，為哈密頓算符，j為目標(biāo)體(3)內(nèi)電流密度分布；

利用歐姆定律，j＝σe＝e/ρ，則有：

滿(mǎn)足的邊界條件為：

其中，γa,b為注入電極位置，γg為除去電極之外的目標(biāo)體邊界，a0為電極與目標(biāo)體接觸面積，i為注入電流，n表示區(qū)域的外法向單位矢量；

對(duì)導(dǎo)電物體進(jìn)行空間離散，并給出電阻率的初值[ρ]⁰，將電阻率的初值[ρ]⁰代入公式(5)中，考慮在成像體的邊界施加邊界條件(6)，采用有限元方法求解矢量電位[t]¹；

第四步：定義目標(biāo)函數(shù)

由式(3)和式(4)可得：

根據(jù)公式(7)定義目標(biāo)函數(shù)f(ρ)：

其中，s(r)為通過(guò)測(cè)量聲壓由時(shí)間反演法計(jì)算得到的熱聲源分布，s為可以通過(guò)公式(5)、(6)和(7)，第三步和第四步，由電阻率和矢量電位計(jì)算得到的熱聲源。

第五步：目標(biāo)函數(shù)最優(yōu)化迭代求解

將矢量電位[t]¹代入目標(biāo)函數(shù)，通過(guò)最優(yōu)化算法更新電阻率進(jìn)行迭代，直到誤差滿(mǎn)足最優(yōu)終止條件ε為所要求的精度，最終求解目標(biāo)體電阻率。

第五步的迭代方法為：

(1)若誤差滿(mǎn)足最優(yōu)終止條件則計(jì)算過(guò)程結(jié)束，得到目標(biāo)體電阻率；

(2)若不滿(mǎn)足，采用高斯牛頓法或非線性阻尼最小二乘法，計(jì)算雅克比矩陣，求解電阻率參數(shù)增量△ρ；

(3)將電阻率參數(shù)增量△ρ代入上一個(gè)電阻率值[ρ]ⁱ，得到新的電阻率[ρ]ⁱ⁺¹；

(4)計(jì)算新的矢量電位，反復(fù)迭代下去，直至誤差滿(mǎn)足最優(yōu)終止條件

附圖說(shuō)明

圖1本發(fā)明重建方法圖像重建過(guò)程示意圖；

圖2本發(fā)明重建方法熱聲信號(hào)檢測(cè)示意圖；

圖中：1激勵(lì)源、2注入電極a、3目標(biāo)體、4注入電極b、5超聲換能器。

具體實(shí)施方式

以下結(jié)合附圖和具體實(shí)施方式對(duì)本發(fā)明做進(jìn)一步說(shuō)明。

如圖1所示，本發(fā)明基于最優(yōu)化迭代方法的注入電流式熱聲電阻率圖像重建方法步驟如下：

第一步為獲取注入電流式熱聲信號(hào)：利用超聲換能器獲取目標(biāo)體各斷層熱聲信號(hào)，每個(gè)斷層掃描一周進(jìn)行熱聲信號(hào)的檢測(cè)，即檢測(cè)超聲信號(hào)；第二步根據(jù)注入電流式熱聲聲壓波動(dòng)方程，獲取目標(biāo)體每個(gè)斷層上的熱聲源分布，利用插值法獲取目標(biāo)體整體的熱聲源分布；第三步對(duì)導(dǎo)電物體的電阻率進(jìn)行空間離散，設(shè)定電阻率的初值，根據(jù)歐姆定律，利用線性有限元方程重建的到目標(biāo)體的矢量電位；第四步根據(jù)熱聲源和電流密度關(guān)系，結(jié)合電流連續(xù)性定理，獲取電阻率所滿(mǎn)足的關(guān)于熱聲源與矢量電位的方程，定義目標(biāo)函數(shù)；第五步將求解的矢量電位代入目標(biāo)函數(shù)，利用目標(biāo)函數(shù)最小化原則，求解導(dǎo)電物體的電阻率。具體分述如下：

第一步：獲取導(dǎo)電物體熱聲信號(hào)

如圖2所示，激勵(lì)源1通過(guò)注入電極a2、注入電極b4對(duì)目標(biāo)體3注入脈沖電流，目標(biāo)體3在電流作用下產(chǎn)生焦耳熱，進(jìn)而產(chǎn)生熱膨脹，激發(fā)超聲信號(hào)，超聲信號(hào)通過(guò)耦合劑耦合到超聲換能器5，超聲換能器5接收到信號(hào)后通過(guò)檢測(cè)系統(tǒng)進(jìn)行信號(hào)的放大、濾波、采集和存儲(chǔ)，超聲換能器5在控制器控制下對(duì)目標(biāo)體進(jìn)行掃描檢測(cè)；

第二步：求解目標(biāo)體熱聲源分布

已知熱聲成像的聲壓波動(dòng)方程：

其中r為超聲換能器位置，p(r,t)是聲壓，cs為介質(zhì)中的聲速，cp為目標(biāo)體3的比熱容，β為目標(biāo)體3的熱膨脹系數(shù)，δ(t)是狄拉克函數(shù)，s(r)為熱聲源分布，t為時(shí)間，為拉普拉斯算符；

熱聲源分布的時(shí)間反演法重建公式為：

其中r為標(biāo)量，r＝|r′-r|，r為矢量，er為單位矢量，r'是超聲探頭的位置，r為熱聲源位置，sd是超聲探頭所在的平面，p′是聲壓對(duì)時(shí)間的一階導(dǎo)數(shù)，n是r′位置sd的法線矢量，β為目標(biāo)體的熱膨脹系數(shù)；

選取目標(biāo)體的某一斷層面z＝z1，超聲換能器在此斷層面上進(jìn)行圓周掃描，采集超聲信號(hào)，利用方程(1)、(2)，求解z＝z1斷層面上的熱聲源分布s(x,y,z1)。移動(dòng)超聲換能器在不同斷層面上進(jìn)行掃描檢測(cè)，求解不同斷層面上的熱聲源分布，目標(biāo)體3上整體熱聲源分布可由分層檢測(cè)計(jì)算或z方向上插值得到；

第三步：獲取矢量電位

熱聲源s同時(shí)是電阻率和電流密度的函數(shù)，可以表示為：

s＝ρj²＝ρj·j(3)

其中，ρ為目標(biāo)體3的電阻率，j為目標(biāo)體3內(nèi)電流密度分布；

根據(jù)電流連續(xù)性定理，電流的散度為零，引入矢量電位，則有：

其中，t為矢量電位，為矢量電位的旋度，為哈密頓算符；

利用歐姆定律，j＝σe＝e/ρ，則有：

滿(mǎn)足的邊界條件為：

其中，γa,b為注入電極位置，γg為除去電極之外的目標(biāo)體邊界，a0為電極與目標(biāo)體接觸面積，i為注入電流，n表示區(qū)域的外法向單位矢量；

對(duì)導(dǎo)電物體進(jìn)行空間離散，并給出電阻率的初值[ρ]⁰，將電阻率的初值[ρ]⁰代入公式(5)中，考慮在成像體的邊界施加邊界條件(6)，采用有限元方法求解矢量電位[t]¹。

第四步：定義目標(biāo)函數(shù)

由式(3)和式(4)可得：

根據(jù)公式(7)定義目標(biāo)函數(shù)f(ρ)：

其中，s(r)為通過(guò)測(cè)量聲壓由時(shí)間反演法計(jì)算得到的熱聲源分布，s為可以通過(guò)公式(5)、(6)和(7)，以及第三步和第四步，由電阻率和矢量電位計(jì)算得到的熱聲源。

第五步：目標(biāo)函數(shù)最優(yōu)化迭代求解

將矢量電位[t]¹代入目標(biāo)函數(shù)，通過(guò)最優(yōu)化算法更新電阻率進(jìn)行迭代，直到誤差滿(mǎn)足最優(yōu)終止條件ε為所要求的精度，最終求解目標(biāo)體電阻率。

第五步的迭代方法為：

(1)若誤差滿(mǎn)足最優(yōu)終止條件則計(jì)算過(guò)程結(jié)束，得到目標(biāo)體電阻率；ε為所要求的精度；

(2)若不滿(mǎn)足，采用高斯牛頓法或非線性阻尼最小二乘法，計(jì)算雅克比矩陣，求解電阻率參數(shù)增量△ρ；

(3)將電阻率參數(shù)增量代入上一個(gè)電阻率值[ρ]ⁱ，得到新的電阻率[ρ]ⁱ⁺¹；

(4)計(jì)算新的矢量電位，反復(fù)迭代，直至誤差滿(mǎn)足最優(yōu)終止條件