• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于電阻抗層析成像的高強度聚焦超聲溫度監(jiān)測技術(shù)?

    2017-09-07 20:56:32郭各樸宿慧丹丁鶴平馬青玉
    物理學報 2017年16期
    關(guān)鍵詞:電導(dǎo)率平面電極

    郭各樸 宿慧丹 丁鶴平 馬青玉

    (南京師范大學物理科學與技術(shù)學院,南京 210023)

    基于電阻抗層析成像的高強度聚焦超聲溫度監(jiān)測技術(shù)?

    郭各樸 宿慧丹 丁鶴平 馬青玉?

    (南京師范大學物理科學與技術(shù)學院,南京 210023)

    (2017年4月8日收到;2017年6月7日收到修改稿)

    作為一種對正常組織無損傷且不易引起癌細胞轉(zhuǎn)移的非入侵腫瘤治療手段,高強度聚焦超聲(HIFU)治療過程中焦域的溫度監(jiān)測是實現(xiàn)劑量精準控制的關(guān)鍵.本文基于生物組織的溫度-電阻抗的關(guān)系,將電阻抗層析成像(EIT)和HIFU治療相結(jié)合,提出了一種利用組織焦平面的表面電壓實現(xiàn)電阻抗重構(gòu)的檢測技術(shù).建立了HIFU治療和EIT綜合系統(tǒng)模型,在考慮組織的聲吸收條件下,對三維Helmholtz方程在柱坐標下的聲場計算進行了二維簡化,并引入Pennes生物熱傳導(dǎo)方程來計算HIFU焦域的聲壓和溫升分布特性;引入生物組織的溫度-電阻抗關(guān)系,基于麥克斯韋電磁場理論,建立了具有溫度分布HIFU焦域的電流和電壓計算模型,利用恒流注入的邊界條件實現(xiàn)電場計算,獲得焦平面的表面電壓分布.在數(shù)值計算中,利用實驗聚焦換能器參數(shù),模擬了在固定聲功率下組織焦域的聲場和溫度場分布,以及中心和偏心聚焦條件下不同治療時刻的電導(dǎo)率分布;然后通過對稱電極的循環(huán)電流注入,計算了組織模型焦平面內(nèi)的電流密度和電勢分布,獲得了焦平面圓周分布的表面電極電壓;進一步采用修正的牛頓-拉夫遜算法,利用32×32的表面電極電壓實現(xiàn)了焦平面內(nèi)電導(dǎo)率分布的重建.結(jié)果表明,基于溫度-電阻抗關(guān)系的EIT電導(dǎo)率重建技術(shù)不但能準確定位HIFU焦域中心,還能恢復(fù)HIFU治療中焦域的溫度分布,證明了EIT用于HIFU治療中溫度監(jiān)測的可行性,為其療效評估和劑量控制提供了一種無創(chuàng)電阻抗測量和成像新方法.

    高強度聚焦超聲,電阻抗層析成像,溫度監(jiān)測,表面電極電壓

    1 引 言

    高強度聚焦超聲[1?4](high intensity focused ultrasound,HIFU)是一種具有廣闊應(yīng)用前景的無創(chuàng)腫瘤治療技術(shù),它利用超聲波在組織中的穿透性和易聚焦性,將換能器發(fā)射的超聲波匯聚到腫瘤靶區(qū),利用組織的聲熱效應(yīng)產(chǎn)生65?C以上的高溫,實現(xiàn)腫瘤組織在短時間內(nèi)凝固性壞死從而達到治療腫瘤的目的.在HIFU治療過程中,既要殺滅腫瘤細胞,又不損傷周圍的正常組織,準確的溫度控制是關(guān)鍵,其中實時溫度和療效監(jiān)測對HIFU的臨床應(yīng)用具有重要的意義.

    為了避免在HIFU治療中插入測溫探針,減少探針對HIFU聲場的影響并降低癌細胞轉(zhuǎn)移概率[5],國內(nèi)外學者提出了多種無創(chuàng)測溫技術(shù).微波測溫技術(shù)利用組織溫度和熱輻射的關(guān)系,通過體外測量體內(nèi)的熱輻射來推測體內(nèi)溫度,但滲透深度有限,測量精度較差.磁共振成像測溫技術(shù)[6]通過溫度相關(guān)的擴散系數(shù)、質(zhì)子共振頻率或弛豫時間的測量實現(xiàn)組織溫度圖像的重建,具有無創(chuàng)傷、無電離輻射、高溫度分辨率的優(yōu)點,但其時間分辨率不高.通過將超聲測溫技術(shù)[7?11]和HIFU系統(tǒng)進行融合,利用聲速和回波時移以及非線性等參數(shù)實現(xiàn)溫度監(jiān)測,但它們的溫度變化較小,測量精度較低.近年來B超[12]被用來進行HIFU定位引導(dǎo),監(jiān)測治療前后組織的供血變化,但是還不能實現(xiàn)高精度的溫度監(jiān)控和實時療效評價.

    研究表明,生物組織也是一種導(dǎo)電體,在低頻信號(<1 MHz)激勵下正常組織的電導(dǎo)率為0—0.5 S/m[13],其導(dǎo)電能力和溫度有著明顯的對應(yīng)關(guān)系.國內(nèi)外研究發(fā)現(xiàn)生物組織的電阻抗與溫度變化存在近似線性關(guān)系[14?16],其電導(dǎo)率隨著溫度的升高而增大,其溫度-電阻抗系數(shù)(teMperature iMpedance variation factor,TIVF)約為?2%/?C,在70?C熱凝固變性時突變到?14.7%/?C.37?C和70?C時生物組織的電導(dǎo)率分別為0.41 S/m和0.79 S/m,其變化幾乎達到100%.和聲阻抗相比,生物組織的電阻抗范圍和變化率更高,這為無創(chuàng)測溫提供了物理基礎(chǔ).

    在本實驗室近期的研究中,利用組織模型的電阻抗相對變化[17](relative iMpedance variation,RIV)進行了HIFU焦域的電阻抗測量,結(jié)果表明當HIFU聲功率一定時,組織模型的RIV和治療時間呈線性關(guān)系;在達到HIFU治療療效時(焦域徑向±0.4 mm區(qū)域達到70?C),組織模型的RIV和HIFU聲功率呈現(xiàn)反比關(guān)系,而所需治療時間和聲功率的平方成反比,證明組織模型的RIV可以用來實現(xiàn)HIFU治療過程中組織焦域的溫度監(jiān)測和療效評估.但是由于HIFU焦域的尺寸很小,其電阻抗變化對模型RIV的影響較小,因此對電阻抗測量系統(tǒng)的靈敏度提出了更高的要求.同時,由于RIV測量的是HIFU治療過程中組織模型的整體電阻抗變化,很難實現(xiàn)高精度的焦點定位和焦域溫度反演,進一步將HIFU治療過程中的溫度-電阻抗關(guān)系和焦域的準確定位相結(jié)合,進行電阻抗/溫度圖像的重建在HIFU治療的溫度監(jiān)測中具有重大的研究價值和應(yīng)用前景.

    近年來,電阻抗成像技術(shù)(electrical iMpedance tomography,EIT)[18?20]得到了迅速的發(fā)展,并在生物醫(yī)學領(lǐng)域得到了廣泛的應(yīng)用.EIT根據(jù)生物組織內(nèi)部電阻抗分布的不同,通過對物體表面電流和電壓的測量來重建物體內(nèi)部電導(dǎo)率的分布.由于EIT技術(shù)不使用射線,激勵電流在安全范圍以內(nèi),對人體無害,成本低廉,且可以實現(xiàn)HIFU治療聲學系統(tǒng)和電阻抗測量的電學系統(tǒng)的完全分離,減少相互干擾,提高測量的可靠性,因此用EIT來測量HIFU焦域的電阻抗分布是一種新方法.

    本文基于組織的溫度-電阻抗關(guān)系,建立了HIFU治療和EIT測量系統(tǒng)模型,研究了HIFU的聲場和溫度場分布特性,分析了HIFU焦域的電阻抗分布,提出了一種基于EIT的HIFU焦域溫度監(jiān)測技術(shù).首先,利用有限元建立了系統(tǒng)的仿真模型,計算了HIFU過程中組織模型的聲場和溫度場分布;然后基于組織的溫度-電阻抗關(guān)系,將聲學系統(tǒng)擴展到電學系統(tǒng),利用焦平面表面對稱電極的循環(huán)電流注入,計算HIFU焦域的電流密度和電勢分布,得到焦平面上的表面電極電壓;最后采用修正的牛頓-拉夫遜算法,利用模擬的32×32表面電極電壓進行焦平面的電導(dǎo)率分布重建.中心和偏心聚焦條件下不同治療時間的仿真結(jié)果表明,重建的電導(dǎo)率分布能較好地反映出HIFU的焦點位置和焦域的電導(dǎo)率變化,間接反映了HIFU治療過程中焦域的溫升情況.研究結(jié)果證明了EIT用于HIFU焦域電阻抗分布重建的可行性,為HIFU治療中焦域的精確定位和溫度監(jiān)測以及療效評估提供了一種電阻抗測量和成像的新方法.

    2 原理與方法

    Crum等[21]對超聲非線性對HIFU焦域溫升的研究表明,由于熱傳導(dǎo)的影響,超聲非線性對HIFU焦域的組織溫升影響較小;另外Myers和Soneson[22?24]證明n階諧波所產(chǎn)生的溫升和其諧波熱源的比值為log(n)/n,在HIFU治療中非線性效應(yīng)產(chǎn)生明顯焦點溫升誤差(達到10%)的聲功率臨界值為115W.本研究所用的聲功率為15.68W時,非線性溫升誤差遠小于2%,可以采用線性模型進行聲場計算和溫升估計,為進一步HIFU焦域的電壓測量和EIT重建提供基礎(chǔ).

    考慮到媒質(zhì)的黏滯性對聲波的吸收,聚焦超聲換能器所激發(fā)的聲波在生物組織內(nèi)部傳播時會產(chǎn)生能量衰減和溫度升高,其聲壓p滿足一維波動方程[25?27]:

    圖1 HIFU治療和EIT測量系統(tǒng)原理圖Fig.1.Sketch Map of the coMp rehensive systeMof H IFU therapy and EIT Measu reMent.

    HIFU治療系統(tǒng)的原理如圖1(a)所示,聚焦換能器表面任意超聲陣元的波動方程可以簡化為Helmholtz[26]方程:

    由于換能器聲場的軸對稱性,(3)式可以在二維軸對稱圓柱坐標下簡化為

    其中,er和ez分別是沿r和z方向的單位矢量.將cc= ω/k′代入(4)式,將公式左右兩邊同乘(?r),得到黏滯介質(zhì)中的簡化方程:

    眾所周知,超聲的聲熱效應(yīng)[26?28]是超聲熱療的基本原理.一定強度的超聲在組織媒介中傳播時,部分聲能被組織吸收轉(zhuǎn)化為熱能,這種被吸收的聲能就是HIFU治療的熱源.假設(shè)聲波在沿z軸傳播過程中通過一個面積為d S、長度為d z的體積元媒質(zhì),如z處的聲強為I1,z+d z處的聲強為I2,考慮到媒質(zhì)對聲波能量的吸收,可知I1>I2.在一維聲傳播條件下,聲場中單位時間內(nèi)單位體積的媒質(zhì)所吸收的聲能量可以由聲強表示:

    在三維空間的聲傳播中,熱源則可以由聲強的空間梯度[28]計算,

    在平面波近似下,聲波沿z方向傳播,若入射超聲強度I0,傳播距離z處的聲強可以表示為I=I0exp(?2αaz),其中αa為吸收系數(shù),代入(7)式得到[28]:

    在生物組織中,入射聲能量的損失一般由聲衰減系數(shù)α表征,其為吸收系數(shù)αa和散射系數(shù)αs的總和.實際計算中很難區(qū)分組織對聲能量的吸收和散射分量,一般直接將衰減系數(shù)等效為吸收系數(shù)[28],即α≈αa.在忽略介質(zhì)中的熱耗散和熱傳導(dǎo),傳播距離z處的超聲熱源表示為[28]

    在聲場計算中,HIFU焦點處的聲強可以由聲壓各次諧波幅值表示[28]:

    其中,ρt和ct分別為組織密度和聲速,Cn是n次諧波的聲壓幅值.在本研究中所用超聲功率不大,忽略非線性影響,(10)式可以簡化為用基波聲壓計算.

    在不考慮血液流動影響的前提下,將組織所吸收的熱量應(yīng)用到Pennes生物熱傳導(dǎo)方程[28,29]中,

    其中Kt是組織熱導(dǎo)率,T為組織溫度,T0為初始溫度.由于HIFU的熱效應(yīng),在焦點處形成中心溫度最高、周圍溫度較低、具有明顯的溫度梯度分布的橢球狀(mm)的溫升焦域.為了定量分析組織電阻抗隨溫度的變化,將TIVF[14?16]應(yīng)用到電導(dǎo)率中,得到組織的溫度-電導(dǎo)率分段函數(shù)為

    可見,隨著HIFU治療時間的增長,生物組織焦域的溫度逐漸升高,電導(dǎo)率隨溫度的升高而增大,在70?C組織凝固時產(chǎn)生快速的電導(dǎo)率提升,在焦域內(nèi)形成電導(dǎo)率的梯度分布,這為EIT在HIFU中的溫度監(jiān)測提供了物理基礎(chǔ).

    基于EIT的HIFU焦域溫度監(jiān)測技術(shù)的正問題研究包括聲場、溫度和電導(dǎo)率分布,焦平面內(nèi)的電勢分布以及邊界電極的電壓分布.圖1(b)顯示了HIFU治療中組織模型(半徑為R)焦平面的電導(dǎo)率分布,在邊界圓周上均勻設(shè)置N個電極,用對稱電極1/17電流注入,計算焦平面內(nèi)的電流密度和電勢分布,得到32個邊界電極的電壓值;然后利用對稱電極的循環(huán)電流注入,重復(fù)以上過程,得到邊界電極電壓Vij(i,j=1,2,...,N),其中i和j表示激勵電極和測量電極的序號.

    EIT正問題求解是在已知HIFU焦域的內(nèi)部電阻抗分布的前提下,且具有特殊邊界條件的電場計算[30,31],組織模型可以等效為導(dǎo)體,其電場分布滿足麥克斯韋方程組[32],在10—100 kHz的低頻電流激勵下,忽略介電常數(shù)的影響,得到場域的表達式[33]:

    其中,σ(T)為場域內(nèi)電導(dǎo)率分布,T為節(jié)點的溫度;Γ1為第一類邊值條件,表示已知位函數(shù)在場域邊界上各點的值,φ為場域內(nèi)電勢分布函數(shù),φ0為給定的邊界電位,初始時認為給定點邊界電位為0;Γ2為第二類邊值問題,表示已知位函數(shù)在場域邊界上各點的法向?qū)?shù)值;Jn是給定邊界注入的電流密度.對(13)式采用變分法求解[33],建立拉普拉斯方程的泛函

    根據(jù)格林定理,考慮模型的邊界條件,并使整個場域的總能量最小,得到

    在HIFU治療中,焦域的電導(dǎo)率分布隨著溫升而改變,結(jié)合激勵電極的邊界條件,利用有限元算法可以計算出組織模型內(nèi)各剖分單元和節(jié)點的電壓,進一步獲得邊界電極測量電壓Vij.

    基于EIT的HIFU焦域溫度監(jiān)測的逆問題是通過已經(jīng)獲得的模型邊界電極電壓和電流激勵模式,利用Vij重建出模型內(nèi)的電阻抗分布.修正的牛頓-拉夫遜(MNR)算法[18,34,35]是通過不斷迭代來改變電阻率分布,進而使目標函數(shù)(重建的邊界電極電壓和測量電壓之間誤差范數(shù)的平方)最小.如組織的電導(dǎo)率為σ(T),則相應(yīng)的電阻率分布為T.在HIFU治療的正問題中,通過對稱電極的循環(huán)電流注入得到32×32邊界電極測量電壓Vij;在MNR電阻抗重建的逆問題中,假設(shè)相應(yīng)對稱電極電流注入時計算得到邊界電極電壓為Uij(ρ),則目標函數(shù)為

    其中N=32.通過不斷迭代使目標函數(shù)最小,得到焦平面的穩(wěn)定電阻率分布.為使目標函數(shù)最小,令f′(ρ)=U′(ρ)(U(ρ)?V)=0,其中U′(ρ)稱為Jacobian矩陣[36].對f′(ρ)泰勒級數(shù)展開,并保留線性項得到f′≈ f′(ρk)+f′′(ρk)?ρk,其中?ρk=ρk+1? ρk,f′′(ρk)=[J(ρk)]TJ(ρk)為Hessian矩陣,?ρk= ?[[J(ρk)]TJ(ρk)]?1[J(ρk)]T[U(ρk)?V],得到MNR的迭代公式為

    其中k是迭代次數(shù). 為了避免對病態(tài)[J(ρk)]TJ(ρk)求逆,在重建中引入Tikhonov正則化[33]修正,并將其補償項表示為ρ的平方函數(shù)形式,得到重構(gòu)電阻率分布的迭代公式為

    將MNR重建方法應(yīng)用到HIFU焦域的電阻率重建中,先基于EIT正問題得到的邊界電極測量電壓Vij,假設(shè)模型焦平面內(nèi)電阻率分布均勻,利用有限元模型計算得到邊界電極電壓Uij(ρ),建立目標函數(shù)f(ρ),通過求解不同電阻率分布下的雅克比矩陣和海森矩陣,獲得新的迭代電阻率分布,直到目標函數(shù)小于預(yù)設(shè)值,此時的電阻率分布ρk即為重建出的組織模型焦平面的電阻率分布,進一步利用σ(T)=1/ρk可以計算出電導(dǎo)率分布.

    3 數(shù)值計算

    如圖1所示,由于換能器和聲傳播的對稱性,在有限元[37]計算中,HIFU聲場和溫度場以及組織電導(dǎo)率分布采用二維軸對稱的柱坐標模型,其中軸向z是聲傳播方向,r是半徑方向,仿真區(qū)域為超聲換能器、水域環(huán)境以及組織模型.通過調(diào)整換能器表面振速控制輸出聲功率,計算組織焦域的聲場、溫度場以及電導(dǎo)率分布.為了在保證計算精度的前提下提高計算速度,水域環(huán)境區(qū)域剖分網(wǎng)格尺寸為λWater/4,組織區(qū)域和焦域的剖分網(wǎng)格尺寸分別λTissue/4和λTissue/8.設(shè)置聚焦超聲換能器的直徑和焦距均為10 cm,中心頻率1.13 MHz.直徑和高度分別為32mm和35 mm的圓柱形組織模型中心放在超聲換能器的焦域處.在模型焦平面的表面均勻設(shè)置32個電極測量表面電極的電壓.在HIFU治療中,隨著超聲的作用,模型內(nèi)焦域的溫度升高,其電導(dǎo)率隨之提高,形成電導(dǎo)率的梯度分布.模型采用仿組織透明凝膠[38],其物理參數(shù)和人體組織較為接近.計算中水和凝膠組織以及人體組織的相關(guān)參數(shù)列于表1.結(jié)合HIFU系統(tǒng)的實驗參數(shù),將換能器表面振幅設(shè)定為10.2 nm,通過有限元計算得到焦域處的聲壓以及聲強分布,定義焦平面內(nèi)聲壓衰減6 dB的面積為有效截面面積,通過I(x,y)d A計算得到聲功率為15.68 W[17],并以此作為聲源參數(shù)進行相關(guān)計算.

    表1 溫度為293 K(20?C)時的仿真參數(shù)Tab le 1.ParaMeters used in siMu lation at the temperature of 293 K(20?C).

    通過有限元仿真得到固定聲功率15.68 W時HIFU焦域的軸向剖面和徑向焦平面的聲壓分布,結(jié)果如圖2(a)和圖2(b)所示,可見HIFU焦域呈現(xiàn)橢球狀,焦點處的聲壓最大,隨著軸向和徑向范圍的擴大,聲壓逐漸降低.結(jié)合(11)式,計算不同治療時間(?t)組織焦域的溫升,得到如圖3所示的焦平面溫度分布.在?t=1 s時,焦域半徑小于0.5 mm,焦域中心的溫度較低,未達到70?C.隨著治療時間的延長,焦域的能量逐漸積累,焦點及周圍組織的溫度不斷升高,同時由于組織的熱擴散,焦域面積不斷擴大,在焦平面上形成圓形的焦斑.在?t=2 s時,焦域中心0.2 mm半徑范圍內(nèi)的溫度超過70?C,達到了治療效果.在?t=3 s時,焦域半徑超過1.3 mm,其中超過在半徑0.5 mm的范圍內(nèi)溫度超過70?C.進一步將組織的溫度-電阻抗關(guān)系應(yīng)用到圖3(a)中,得到如圖3(b)所示的不同治療時間的焦平面電導(dǎo)率分布.可見隨著的延長,HIFU焦域的溫度逐漸升高,電導(dǎo)率相應(yīng)提高,形成焦域中心導(dǎo)電能力強,周圍導(dǎo)電性能弱的梯度分布.

    圖2 (網(wǎng)刊彩色)固定聲功率為15.68W時,HIFU的(a)軸向剖面和(b)焦平面的聲壓分布Fig.2.(color on line)Pressu re distributions of(a)axial p rofi le and(b)focal p lane for HIFU at a fi xed acoustic power of 15.68 W.

    圖3 (網(wǎng)刊彩色)不同治療時刻HIFU焦平面的(a)溫度和(b)電導(dǎo)率分布及(c)偏心聚焦的電導(dǎo)率分布Fig.3.(color on line)Focal d istribu tions of(a)teMperatu re and(b)conductivity for centric HIFU therapy,as well as(c)the conductivity distributions for eccentric HIFU therapy at diff erent treatMent tiMes.

    為了證明研究的普適性,將焦點沿著徑向移動8 mm形成HIFU偏心聚焦,得到了如圖3(c)所示的不同治療時間HIFU焦平面的電導(dǎo)率分布.由于組織媒質(zhì)相同,模型內(nèi)不同位置的聲學特性和電阻抗特性相同,因此除了焦域位置的移動,焦平面上焦點處的電導(dǎo)率分布和圖3(b)基本一致.

    在HIFU治療的同時,通過對稱電極的電流注入,在模型內(nèi)部形成相應(yīng)的電場分布,其變化能反映焦域的溫度變化.圖4(a)和圖4(b)分別顯示了HIFU中心和偏心聚焦時,在電極1/17電流注入,?t=1,2,3 s時焦平面內(nèi)的電流密度和電勢分布.如圖4(a)的虛線環(huán)所示,隨著治療時間的延長,焦域的溫度和電導(dǎo)率升高,原來較為均勻分布的電流向焦域的中心匯聚,形成向上彎曲的電勢等位線.對于如圖4(b)所示的偏心聚焦,在焦域內(nèi)也形成匯聚的電流和彎曲的電勢等位線,其彎曲程度由焦域的電導(dǎo)率分布決定.

    圖4 (網(wǎng)刊彩色)對稱電極1/17電流注入,HIFU(a)中心和(b)偏心聚焦時焦平面的電流密度和電勢分布Fig.4.(color on line)D istributions of cu rrent density and electrical potential of the focal p lane With(a)centric and(b)eccentric focusing for the current in jection froMthe electrodes 1/17.

    為了進行基于EIT的電阻抗分布重建,采用如圖1所示的三角形網(wǎng)格剖分[39],其中剖分節(jié)點數(shù)為1225,剖分單元數(shù)為2320,并將有限元計算得到的焦平面電導(dǎo)率分布和和電勢分布導(dǎo)入到所建立的剖分網(wǎng)格中,將組織模型的焦平面圓周32等分,獲得循環(huán)對稱電極電流注入時的邊界電極測量電壓Vij.圖5(a)和圖5(b)分別顯示了對稱電極1/17和9/25電流注入時32個邊界電極上的測量電壓分布.可見隨著治療時間的延長,焦域溫度升高,電阻率降低,相同電流注入時邊界電極電壓降低.在改變電流注入電極后,邊界電極的電壓分布發(fā)生相應(yīng)的變化.然而,由于相對于模型尺寸,焦域面積很小,電阻抗差異較小,因此電極電壓分布隨著治療時間和溫度以及位置的變化不明顯.圖5(a)和圖5(b)的內(nèi)插局部放大圖分別顯示了兩種情況下測量得到電極7/8和17/18的電壓分布,可以看出相應(yīng)電極的電壓產(chǎn)生了幅度較小的變化(mV量級).

    圖5 (網(wǎng)刊彩色)對稱電極(a)1/17和(b)9/25電流注入時,模型焦平面表面電極的電壓分布Fig.5.(color on line)E lectrical voltage distributions of the surface electrodes in the focal p lane With the current in jections froMthe symMetric electrodes of 1/17 and 9/25.

    圖6 (網(wǎng)刊彩色)H IFU中心聚焦時,(a)模型焦平面的網(wǎng)格化電導(dǎo)率分布和(b)重建的電導(dǎo)率分布Fig.6.(color on line)(a)G rid d istributions and(b)reconstructed distributions of conductivity in the focal p lane for centric focusing.

    在EIT的電阻抗重建中,根據(jù)HIFU治療中組織焦域的電導(dǎo)率改變,設(shè)定初始均勻電導(dǎo)率為0.73 S/m,利用Matlab編程計算焦平面上的電極電位分布Uij(ρ),建立目標函數(shù),并通過求解每次迭代電阻率分布下的雅克比矩陣和海森矩陣獲得新的迭代電阻率分布,當?shù)螖?shù)達到30次時,獲得了較穩(wěn)定的重建結(jié)果.圖6(a)和圖6(b)分別顯示了HIFU中心聚焦時模型焦平面的網(wǎng)格化電導(dǎo)率分布和基于邊界電極電壓的重建電導(dǎo)率分布,為了進行重建效果比較,HIFU偏心聚焦時的相應(yīng)結(jié)果在圖7(a)和圖7(b)中給出.可見重建結(jié)果能夠較為準確地反映出治療區(qū)域的位置和形狀以及電導(dǎo)率的變化趨勢,能夠?qū)崿F(xiàn)準確的焦點定位,重建焦域的尺寸和電導(dǎo)率隨治療時間的增長而變大,和模型計算結(jié)果的變化趨勢一致,能夠反映出不同治療時間HIFU焦域的電導(dǎo)率差異.但是由于中心聚焦的焦域具有完全對稱性,其邊界電極電壓的偏差較小,因此重建焦域的精度不如偏心聚焦效果好.

    由于在電導(dǎo)率的網(wǎng)格化處理中使用了三角形內(nèi)部的數(shù)據(jù)平均,所得到的最大電導(dǎo)率小于有限元的仿真結(jié)果.同時由于重建算法的病態(tài)性,所得到的電導(dǎo)率小于實際值.為了比較重建效果,將重建結(jié)果用對應(yīng)時刻模型電導(dǎo)率的最大值進行歸一化處理,得到如圖8所示的不同治療時間模型焦平面的網(wǎng)格化電導(dǎo)率徑向分布和重建電導(dǎo)率的徑向分布.雖然重建結(jié)果比實際電導(dǎo)率范圍寬,但是仍然可以精確定位焦點位置,反映焦平面內(nèi)焦域溫度和電導(dǎo)率的分布;隨著治療時間的增加,焦域的電導(dǎo)率變化增大,重建的電導(dǎo)率分布和模型電導(dǎo)率分布趨勢一致,能夠較好地反映組織電導(dǎo)率隨溫度的變化特性.

    圖7 (網(wǎng)刊彩色)HIFU偏心聚焦時,(a)模型焦平面的網(wǎng)格化電導(dǎo)率分布和(b)重建的電導(dǎo)率分布Fig.7.(color on line)(a)G rid d istributions and(b)reconstructed distributions of conductivity in the focal p lane for eccentric focusing.

    圖8 (網(wǎng)刊彩色)H IFU(a)中心和(b)偏心聚焦時,不同治療時間組織模型焦平面的網(wǎng)格化電導(dǎo)率歸一化徑向分布和重建電導(dǎo)率的歸一化徑向分布Fig.8.(color on line)NorMalized rad ial grid d istributions and reconstructed d istribu tions of conductivity in the focal p lane for(a)centric and(b)eccentric focusing.

    4 討 論

    本研究中由于HIFU的聲功率不太高,超聲非線性所產(chǎn)生溫升對線性聲場所產(chǎn)生的誤差較小,為了簡化溫度場的計算,選用黏滯媒質(zhì)中的線性方程求解,取得較好的計算效果.但實際中,HIFU聲場一般為非線性,要獲得更為準確的聲場和溫度場,甚至考慮聲空化的影響,在進一步的研究中采用Westervelt[40],KZK[41,42]和SBE[43]方程求解,會得到更為精確HIFU焦域的溫度分布.

    理論和模擬結(jié)果證明,雖然測量的邊界電極的電壓變化較小,但利用EIT可以精準定位HIFU焦域的位置和焦域內(nèi)電阻抗的變化趨勢,重建圖像中電導(dǎo)率的差異可以反映焦域的溫度分布特性,在HIFU的療效監(jiān)測中有著廣泛的應(yīng)用前景.另外,當前EIT逆問題求解算法主要包括動態(tài)重構(gòu)算法(例如等位線反投影法)和靜態(tài)重構(gòu)算法(如牛頓-拉夫遜迭代算法)等.鑒于動態(tài)重構(gòu)算法計算速度快和效率高的優(yōu)點,本研究首先采用等位線投影法對32電極的測量數(shù)據(jù)進行了仿真計算,然而由于HIFU焦域的面積小,電導(dǎo)率變化范圍大,而邊界電極的電壓變化小,因此重建算法對電阻抗的變化不敏感,重建圖像的精度差.雖然修正的牛頓-拉夫遜算法是一種靜態(tài)的重建算法,但具有良好的收斂性,重建結(jié)果較動態(tài)算法更好,但由于該算法重建問題的病態(tài)性會隨剖分網(wǎng)格的增加而增加,電阻抗重建的網(wǎng)格剖分不能過密,限制了重建精度的提高,因此需要進一步研究EIT重建算法,在較少的網(wǎng)格剖分和表面電極條件下,獲得更高精度的重建圖像.

    本研究使用循環(huán)對稱電極電流注入,其電流密度和電勢分布對位置較為敏感,當HIFU焦域在模型中心時,不同位置激勵下的表面電極電壓的差異較小,因此電阻抗重建的精度和分辨率較低.為了獲得較大的表面電極電壓差異,需要改進電極注入方式,實現(xiàn)了多種激勵模式下的電阻抗分布重建.當前EIT理論較為成熟,算法簡便,設(shè)備廉價,對生物組織傷害小,重建速度較快,本研究把EIT和HIFU相結(jié)合,實現(xiàn)電學系統(tǒng)和聲學系統(tǒng)完全分離,對HIFU治療系統(tǒng)干擾小,可以同步實現(xiàn)HIFU治療和電阻抗測量,利用組織溫度-電導(dǎo)率的關(guān)系,可以實現(xiàn)HIFU治療過程中焦域溫度的反演.進一步研究中,將模型對稱電極的RIV測量和邊界電極電壓監(jiān)測相結(jié)合,可以有效提高焦點定位和焦域溫度監(jiān)測的準確性,為HIFU療效評估和劑量控制提供實時電阻抗測量和成像新方法.

    5 結(jié) 論

    本文針對HIFU治療中的無創(chuàng)溫度監(jiān)測難題,基于組織的溫度-電阻抗關(guān)系,提出了一種基于EIT的HIFU焦域的溫度和療效監(jiān)測技術(shù).本研究將HIFU治療的聲學系統(tǒng)和電阻抗測量的電學系統(tǒng)有機結(jié)合,利用組織模型表面的電阻抗測量來實現(xiàn)HIFU焦域的電阻抗分布的重建,實現(xiàn)溫度和療效估計.建立了HIFU治療和EIT測量系統(tǒng)模型,模擬了固定聲功率(15.68W)時組織焦域的聲場和溫度場分布,以及中心和偏心聚焦的電導(dǎo)率分布.通過對稱電極的旋轉(zhuǎn)電流注入,計算了圓柱組織模型焦平面內(nèi)的電流密度和電勢分布,獲得了焦平面的表面電極電壓分布.在電阻抗重建中,采用修正的牛頓-拉夫遜算法,通過迭代計算實現(xiàn)了焦平面的電導(dǎo)率分布重建.結(jié)果表明,重建的電導(dǎo)率分布能準確反映HIFU焦域的位置和范圍,體現(xiàn)隨治療時間增長的電導(dǎo)率變化,證明EIT在溫度監(jiān)測中應(yīng)用的可行性.本研究為HIFU治療中焦點的精確定位,焦域的實時溫度監(jiān)測和療效評估以及劑量控制提供了一種無創(chuàng)電阻抗測量和成像新方法.

    [1]Hutchinson L 2011 Nat.Rev.C lin.Oncol.8 385

    [2]K ennedy J E 2005 Nat.Rev.Canc.5 321

    [3]Q ian S Y,Wang H Z 2001 Acta Phys.Sin.50 501(in Chinese)[錢盛友,王鴻樟 2001物理學報 50 501]

    [4]Gavrilov L R 2013 J.Acoust.Soc.AMer.133 4348

    [5]Jiang L X,Hu B 2006 Tech.Acoust.25 43(in Chinese)[姜立新,胡兵2006聲學技術(shù) 25 43]

    [6]Shen J,Shen J L,Zou J Z 2007 J.U ltrasound C lin.Med.9 486(in Chinese)[沈潔,申俊玲,鄒建中2007臨床超聲醫(yī)學雜志9 486]

    [7]Ye G,SMith P P,Noble J A 2010 U ltrasound Med.Biol.36 234

    [8]Daniels MJ,Varghese T,Madsen E L,Zagzebski J A 2007 Phys.Med.Biol.52 4827

    [9]Fan T B,Zhang D,Zhang Z,Ma Y,Gong X F 2008 Chin.Phys.B 17 3372

    [10]Anand A,KaczkoWski P J 2004 J.Acoust.Soc.AMer.115 2490

    [11]Ma Y,Zhang D,Gong X F,Liu X Z,Ma Q Y,Q iu Y Y 2007 Chin.Phys.16 2745

    [12]Fan L Z,Luo F,Yu D Y,Liu X T,Zhang J,X ie MX 2005 C lin.J.Med.Instru.29 115(in Chinese)[范良志,羅飛,喻道遠,劉夏天,張靜,謝明星 2005中國醫(yī)療器械雜志29 115]

    [13]Gabriel C,PeyMan A,G rant E H 2009 Phys.Med.Biol.54 4863

    [14]Zurbuchen U,HolMer C,LehMann K S,Stein T,Roggan A,Seifarth C,Buh r H J,Ritz J P 2010 In t.J.HypertherMia 26 26

    [15]G riffi ths H,AhMed A 1987 C lin.Phys.Physiol.Meas.8 147

    [16]Cai H,You F S,Shi X T,Fu F,Liu R G,Tang C,Dong X Z 2010 Chin.Med.Equip.J.31 8(in Chinese)[蔡華,尤富生,史學濤,付峰,劉銳崗,湯池,董秀珍2010醫(yī)療衛(wèi)生裝備31 8]

    [17]Su H D,Guo G P,Ma Q Y,Tu J,Zhang D 2017 Chin.Phys.B 26 054302

    [18]Li K Q 2015 M.S.Dissertation(Nan jing:Nan jing University of Posts and TelecomMunications)(in Chinese)[李凱強2015碩士學位論文 (南京:南京郵電大學)]

    [19]Xu G X 2004 Ph.D.D issertation (Chongqing:Chongqing University)(in Chinese)[徐管鑫 2004博士學位論文(重慶:重慶大學)]

    [20]Zhang L 2013 Ph.D.D issertation(Nan jing:Nan jing University of Science and Technology)(in Chinese)[張麗2013博士學位論文(南京:南京理工大學)]

    [21]Curra F P,Mourad P D,Khokh lova V A,CruML A 2000 IEEE Trans.U ltrason.Ferroelect.Freq.Con trol 47 1077

    [22]Soneson J E,Myers MR 2010 IEEE Trans.U ltrason.Ferroelect.Freq.Con trol 57 2450

    [23]Myers MR,Soneson J E 2009 J.Acoust.Soc.AMer.126 425

    [24]Soneson J E,Myers MR 2007 J.Acoust.Soc.AMer.122 2526

    [25]Du G H,Zhu Z M,Gong X F 2012 FundaMentals of Acoustics(Nan jing:Nan jing University Press)pp292–305(in Chinese)[杜功煥,朱哲民,龔秀芬 2012聲學基礎(chǔ)(南京:南京大學出版社)第292—305頁]

    [26]Cheng J C 2013 Principles of Acoustics(Beijing:Science Press)pp571–576(in Chinese)[程建春 2013聲學原理 (北京:科學出版社)第571—576頁]

    [27]Wan MX,Zong J Y,Wang S P 2010 BioMedical U ltrasound(Beijing:Science Press)pp649–669(in Chinese)[萬明習,宗瑜瑾,王素品2010生物醫(yī)學超聲學(北京:科學出版社)第649—669頁]

    [28]Zhang D,Guo X S,Ma Q Y,Tu J 2014 FundaMen ta ls of Medical U ltrasound(Beijing:Science Press)pp415–418(in Chinese)[章東,郭霞生,馬青玉,屠娟 2014醫(yī)學超聲基礎(chǔ)(北京:科學出版社)第415—418頁]

    [29]Pennes H H 1948 J.Appl.Physiol.1 93

    [30]Chen Z J,Zhou G L 2014 Inform.ComMun.4 36(in Chinese)[陳姝君,周廣麗 2014信息通信 4 36]

    [31]Chen Z J 2008 Ph.D.D issertation(Nan jing:Nan jing University of Science and Technology)(in Chinese)[陳姝君2008博士學位論文(南京:南京理工大學)]

    [32]Li G Y 2011 M.S.D issertation(Beijing:Beijing Jiaotong University)(in Chinese)[黎光宇2011碩士學位論文(北京:北京交通大學)]

    [33]Xu G Z,Li Y,Yang S,Wu H L,Zhang S,Zhang J J 2010 BioMedical E lectrical IMpedance ToMograph y(Beijing:Machinery Industry Press)pp33–34(in Chinese)[徐桂芝,李穎,楊碩,吳煥麗,張帥,張劍軍2010生物醫(yī)學電阻抗成像技術(shù)(北京:機械工業(yè)出版社)第33—34頁]

    [34]Lin MF 2014M.S.Dissertation(Nan jing:Nanjing University of Posts and TelecomMunications)(in Chinese)[林明鋒2014碩士學位論文(南京:南京郵電大學)]

    [35]Zhang L 2014 Electr.Design Eng.22 184(in Chinese)[張麗2014電子設(shè)計工程22 184]

    [36]Q in S L 2000 CoMputat.Phys.17 314(in Chinese)[秦世倫2000計算物理17 314]

    [37]Ma H,Wang G 2009 COMSOL Mu ltiphysics Basic Operation Guide and FAQ(Beijing:China ComMunications Press)(in Chinese)[馬慧,王剛 2009 COMSOL Mu ltiphysics基本操作指南和常見問題解答 (北京:人民交通出版社)]

    [38]Hu B,Jiang L X,Huang Y 2006 Tech.Acoust.25 613(in English)[胡兵,姜立新,黃瑛2006聲學技術(shù) 25 613]

    [39]Li G 2013 M.S.D issertation(T ian jing:T ian jing University of Science and Technology)(in Chinese)[黎鴿2013碩士學位論文(天津:天津科技大學)]

    [40]Bessonova O,Wilkens V 2013 J.Acoust.Soc.AMer.134 4213

    [41]Sun JM,Yu J,Guo X S,Zhang D 2013 Acta Phys.Sin.62 054301(in Chinese)[孫健明,于潔,郭霞生,章東2013物理學報62 054301]

    [42]Chen T,Fan T,Zhang W,Q iu Y Y,Tu J,Guo X S,Zhang D 2014 J.Appl.Phys.115 114902

    [43]Chen T 2014 Ph.D.D issertation(Nan jing:Nan jing University)(in Chinese)[陳濤 2014博士學位論文 (南京:南京大學)]

    PACS:43.80.Ev,43.35.YbDOI:10.7498/aps.66.164301

    *Pro ject supported by the National Natural Science Foundation of China(G rant Nos.11474166,11604156),the Natu ral Science Foundation of Jiangsu Province,China(G rant No.BK 20161013),the Postdoctoral Science Foundation of China(G rant No.2016M591874)and the Priority AcadeMic PrograMDevelopMent of Jiangsu H igher Education Institu tions,China.

    ?Corresponding author.E-Mail:Maqingyu@njnu.edu.cn

    N on invasive teMperatu re Mon itoring for high intensity focused u ltrasound therapy based on electrical iMpedance toMography?

    Guo Ge-Pu Su Hui-Dan Ding He-Ping Ma Qing-Yu?
    (School of Physics and Technology,Nanjing NorMal University,Nanjing 210023,China)

    8 Ap ril 2017;revised Manuscrip t

    7 June 2017)

    As a neWtreatment modality With little thermal damage and feWcellmetastases to surrounding normal tissues,high intensity focused ultrasound(HIFU)therapy is considered to be one of theMost proMising technologies for tuMor therapy in the 21stcentury.However,noninvasive teMperatureMonitoring for the focal region exhibits great significance of p recise thermal dosage control in HIFU treatment.By combining electrical iMpedancemeasurement and HIFU,an electrical iMpedance toMography(EIT)based teMperature Monitoring Method using surface voltages is proposed to reconstruct the distribution of electrical conductivity inside the focal p lane on the basis of the teMperature dependent electrical iMpedance of tissues.In theoretical study,a coMp rehensive systeMof EIT measurement during HIFU therapy is established.With the consideration of acoustic absorp tion in viscous tissues,three-diMensional HelMholtz equation for HIFU is siMp lified into two-diMensional axisymMetric cylindrical coordinates,and the characteristics of teMperature rising in the focal region are derived using Pennes bio-heat transfer equation.Then,by introducing the teMperatureconductivity relation into tissues,the processing Methods for electrical field and surface voltage in the focal region are constructed With constant current injection froMtwo symmetrical electrodes.In simu lation study,by app lying the experimental paraMeters of the focused transducer,the distributions of acoustic p ressure and teMperature are simu lated at a fixed acoustic power,and then the corresponding distributions of conductivity in the focal p lane are achieved at diff erent treatment times for centric and eccentric focusing.Furthermore,With the simu lations of current density and electricalpotentialgenerated by the rotating current injection froM16 pairsof symMetricalelectrodes,32×32 voltages are detected by the 32 surface electrodes p laced around the focal p lane of theModel.In conductivity iMage reconstruction,themodified NeWton-Raphson(MNR)algorithMis eMp loyed to conduct iterative calcu lation.It shows that With the increase of HIFU treatMent tiMe,the electrical conductivity in the focal region increases accordingly and reaches a MaximuMvalue in the center due to the highest acoustic p ressure and theMost energy accumulation.It is p roved that not only the position of the focal center,but also the conductivity distribution inside the focal region can be restored accurately by the p roposed EIT based reconstruction algorithm.The favorable results deMonstrate the feasibility of teMperaturemonitoring during HIFU therapy,and also provide a neWmethod of evaluating the noninvasive effi cacy and controlling the dose based on electrical iMpedanceMeasureMents.

    high intensity focused ultrasound,electrical impedance tomography,temperaturemonitoring,surface electrode voltage

    10.7498/aps.66.164301

    ?國家自然科學基金(批準號:11474166,11604156)、江蘇省自然科學基金(批準號:BK 20161013)、國家博士后基金(批準號:2016M591874)和江蘇高校優(yōu)勢學科資助的課題.

    ?通信作者.E-Mail:Maqingyu@n jnu.edu.cn

    ?2017中國物理學會C h inese P hysica l Society

    http://Wu lixb.iphy.ac.cn

    猜你喜歡
    電導(dǎo)率平面電極
    基于比較測量法的冷卻循環(huán)水系統(tǒng)電導(dǎo)率檢測儀研究
    低溫脅迫葡萄新梢電導(dǎo)率和LT50值的研究
    參考答案
    關(guān)于有限域上的平面映射
    三維電極體系在廢水處理中的應(yīng)用
    三維鎳@聚苯胺復(fù)合電極的制備及其在超級電容器中的應(yīng)用
    Ti/SnO2+Sb2O4+GF/MnOx電極的制備及性能研究
    參考答案
    高電導(dǎo)率改性聚苯胺的合成新工藝
    電導(dǎo)率法快速測定榨菜鹽分含量
    食品科學(2013年24期)2013-03-11 18:30:38
    精品人妻一区二区三区麻豆 | 嫩草影院新地址| 99久久精品国产国产毛片| 最新在线观看一区二区三区| 麻豆一二三区av精品| 在线a可以看的网站| 亚洲欧美日韩无卡精品| 亚洲无线在线观看| 亚洲四区av| 蜜桃久久精品国产亚洲av| 亚洲成av人片在线播放无| 成年女人看的毛片在线观看| 国产精品人妻久久久影院| 亚洲婷婷狠狠爱综合网| 高清毛片免费看| 免费看光身美女| 午夜福利在线在线| 欧美性感艳星| 欧美丝袜亚洲另类| 日韩欧美在线乱码| 日产精品乱码卡一卡2卡三| 国产成年人精品一区二区| 久久午夜福利片| 特大巨黑吊av在线直播| 亚洲国产精品合色在线| 精品国产三级普通话版| 美女内射精品一级片tv| 又黄又爽又刺激的免费视频.| 婷婷精品国产亚洲av| 3wmmmm亚洲av在线观看| 69av精品久久久久久| 老司机午夜福利在线观看视频| 搡老熟女国产l中国老女人| 嫩草影视91久久| 麻豆国产97在线/欧美| 99在线视频只有这里精品首页| 国产高清三级在线| 99久久久亚洲精品蜜臀av| 五月玫瑰六月丁香| 国产真实乱freesex| 又黄又爽又刺激的免费视频.| 亚洲人成网站在线播放欧美日韩| 国产一区亚洲一区在线观看| 国产精华一区二区三区| 国产精品野战在线观看| 自拍偷自拍亚洲精品老妇| 国产成人福利小说| 亚洲av第一区精品v没综合| 高清午夜精品一区二区三区 | 国产精品一区二区免费欧美| 蜜桃久久精品国产亚洲av| 欧洲精品卡2卡3卡4卡5卡区| 国产综合懂色| 亚洲av五月六月丁香网| 国产女主播在线喷水免费视频网站 | 精品一区二区三区av网在线观看| 99久久精品一区二区三区| 久久久欧美国产精品| 国产黄片美女视频| 99久久久亚洲精品蜜臀av| 97人妻精品一区二区三区麻豆| 久久精品久久久久久噜噜老黄 | 国产精品乱码一区二三区的特点| 国产精品1区2区在线观看.| 日本免费一区二区三区高清不卡| 久久国产乱子免费精品| 亚洲精品色激情综合| 精品人妻熟女av久视频| 午夜视频国产福利| 国产三级在线视频| 国产成人a∨麻豆精品| 别揉我奶头 嗯啊视频| 网址你懂的国产日韩在线| 日韩精品有码人妻一区| 日本一二三区视频观看| 国产探花在线观看一区二区| 1000部很黄的大片| 一级毛片久久久久久久久女| 如何舔出高潮| 久久精品91蜜桃| 99热只有精品国产| 中文字幕精品亚洲无线码一区| 看片在线看免费视频| 午夜福利高清视频| 欧美一区二区亚洲| 九九爱精品视频在线观看| 网址你懂的国产日韩在线| 九九久久精品国产亚洲av麻豆| 久久韩国三级中文字幕| 老女人水多毛片| 中文亚洲av片在线观看爽| 国产精品,欧美在线| 亚洲va在线va天堂va国产| 国产又黄又爽又无遮挡在线| 国产精品综合久久久久久久免费| 精品久久久久久久久av| 最后的刺客免费高清国语| 亚洲美女搞黄在线观看 | 亚洲内射少妇av| 可以在线观看的亚洲视频| 午夜视频国产福利| 久久精品国产自在天天线| 毛片一级片免费看久久久久| 淫妇啪啪啪对白视频| 不卡一级毛片| 成人特级av手机在线观看| 韩国av在线不卡| 欧美激情国产日韩精品一区| 欧美bdsm另类| 午夜福利在线观看吧| 国产又黄又爽又无遮挡在线| 99久久久亚洲精品蜜臀av| 成人特级av手机在线观看| 麻豆精品久久久久久蜜桃| 欧美激情国产日韩精品一区| 久久中文看片网| 久久午夜亚洲精品久久| 欧美人与善性xxx| 色噜噜av男人的天堂激情| 久久天躁狠狠躁夜夜2o2o| 国产又黄又爽又无遮挡在线| 神马国产精品三级电影在线观看| 免费大片18禁| 午夜久久久久精精品| 永久网站在线| 日韩欧美免费精品| 最近最新中文字幕大全电影3| 精品熟女少妇av免费看| 嫩草影院精品99| 久久草成人影院| 黄色日韩在线| a级毛片免费高清观看在线播放| 干丝袜人妻中文字幕| 看十八女毛片水多多多| 十八禁网站免费在线| 最近2019中文字幕mv第一页| 久久精品国产亚洲av涩爱 | 亚洲av第一区精品v没综合| 别揉我奶头 嗯啊视频| 午夜福利在线观看免费完整高清在 | 国产精品爽爽va在线观看网站| 亚洲精品456在线播放app| .国产精品久久| 国产高清视频在线观看网站| 夜夜夜夜夜久久久久| 99久久九九国产精品国产免费| 日韩欧美一区二区三区在线观看| 亚洲欧美中文字幕日韩二区| 欧美激情国产日韩精品一区| 久久精品国产清高在天天线| 一进一出抽搐动态| 三级男女做爰猛烈吃奶摸视频| 五月玫瑰六月丁香| 欧美色欧美亚洲另类二区| 国产精品国产高清国产av| 搡女人真爽免费视频火全软件 | 国产高清不卡午夜福利| 久久久精品欧美日韩精品| 色av中文字幕| 国产精品综合久久久久久久免费| 毛片一级片免费看久久久久| 亚洲精华国产精华液的使用体验 | 久久天躁狠狠躁夜夜2o2o| 免费人成在线观看视频色| 国产真实伦视频高清在线观看| 六月丁香七月| 在线播放无遮挡| 少妇丰满av| 精品欧美国产一区二区三| 免费av观看视频| 在线观看一区二区三区| 最近的中文字幕免费完整| 久久鲁丝午夜福利片| 91麻豆精品激情在线观看国产| 欧美xxxx性猛交bbbb| 国内精品一区二区在线观看| 男女边吃奶边做爰视频| 亚洲av不卡在线观看| 你懂的网址亚洲精品在线观看 | АⅤ资源中文在线天堂| 亚洲成a人片在线一区二区| 黄色日韩在线| 少妇裸体淫交视频免费看高清| 国内久久婷婷六月综合欲色啪| 国产精品综合久久久久久久免费| 欧美xxxx黑人xx丫x性爽| 麻豆一二三区av精品| 欧美bdsm另类| 精品一区二区三区视频在线| 亚洲经典国产精华液单| 18禁在线播放成人免费| 久久人妻av系列| 嫩草影视91久久| 嫩草影院入口| 欧美人与善性xxx| 免费不卡的大黄色大毛片视频在线观看 | 全区人妻精品视频| 岛国在线免费视频观看| 午夜免费男女啪啪视频观看 | 99久久精品一区二区三区| 男人舔女人下体高潮全视频| 在线看三级毛片| 波野结衣二区三区在线| 免费看美女性在线毛片视频| 人妻久久中文字幕网| 日本黄大片高清| 我要搜黄色片| 天堂av国产一区二区熟女人妻| 日韩精品有码人妻一区| 美女高潮的动态| 亚洲第一区二区三区不卡| 乱人视频在线观看| av专区在线播放| 少妇丰满av| 一进一出抽搐动态| 小说图片视频综合网站| 一夜夜www| 天堂网av新在线| 美女大奶头视频| 老司机影院成人| 成人高潮视频无遮挡免费网站| 亚洲综合色惰| 男人的好看免费观看在线视频| 精品一区二区三区av网在线观看| 欧美一区二区精品小视频在线| 我要看日韩黄色一级片| 人妻久久中文字幕网| 亚洲中文字幕日韩| 成年女人永久免费观看视频| 欧美成人a在线观看| 久久精品国产亚洲网站| 亚洲中文日韩欧美视频| 大型黄色视频在线免费观看| a级毛片a级免费在线| 最近2019中文字幕mv第一页| 午夜日韩欧美国产| 人人妻,人人澡人人爽秒播| 日本免费a在线| 成人亚洲精品av一区二区| 大香蕉久久网| 色综合亚洲欧美另类图片| 99riav亚洲国产免费| 精品一区二区三区av网在线观看| 精品人妻偷拍中文字幕| 乱人视频在线观看| 欧美极品一区二区三区四区| 国产视频一区二区在线看| 天堂√8在线中文| 热99re8久久精品国产| 中国美白少妇内射xxxbb| 少妇被粗大猛烈的视频| 男人的好看免费观看在线视频| 国产又黄又爽又无遮挡在线| 夜夜看夜夜爽夜夜摸| 亚洲熟妇中文字幕五十中出| 一进一出抽搐动态| 中文资源天堂在线| 美女大奶头视频| 天天躁日日操中文字幕| 尤物成人国产欧美一区二区三区| 九九在线视频观看精品| 亚洲av免费在线观看| 国产人妻一区二区三区在| 18+在线观看网站| 成人三级黄色视频| 国产精品久久久久久av不卡| 午夜久久久久精精品| 一进一出抽搐gif免费好疼| 久久综合国产亚洲精品| 欧美国产日韩亚洲一区| 久久精品91蜜桃| 淫秽高清视频在线观看| 精品99又大又爽又粗少妇毛片| 国产亚洲精品久久久久久毛片| 亚洲国产欧洲综合997久久,| 国产片特级美女逼逼视频| 免费人成在线观看视频色| 男人舔奶头视频| 欧美一区二区亚洲| 成人亚洲精品av一区二区| 精品熟女少妇av免费看| 国产一区二区激情短视频| 久久99热6这里只有精品| 欧美激情久久久久久爽电影| 久久久成人免费电影| 春色校园在线视频观看| 联通29元200g的流量卡| 国产av在哪里看| av卡一久久| 国产欧美日韩精品一区二区| 午夜影院日韩av| 高清毛片免费观看视频网站| 亚洲性久久影院| 麻豆国产97在线/欧美| 最近最新中文字幕大全电影3| 色噜噜av男人的天堂激情| 久久99热6这里只有精品| 久久久久久九九精品二区国产| 一本久久中文字幕| 国产精品亚洲一级av第二区| 一边摸一边抽搐一进一小说| 看免费成人av毛片| av视频在线观看入口| 亚洲国产日韩欧美精品在线观看| 国产一区二区在线av高清观看| 国产色爽女视频免费观看| 久久这里只有精品中国| 国产精品伦人一区二区| 少妇人妻一区二区三区视频| 国产单亲对白刺激| 少妇人妻一区二区三区视频| 男人狂女人下面高潮的视频| 亚洲欧美日韩无卡精品| 成年女人毛片免费观看观看9| 久久久午夜欧美精品| 精品一区二区三区视频在线| 一区二区三区高清视频在线| 国产成人影院久久av| 婷婷亚洲欧美| 毛片一级片免费看久久久久| 婷婷精品国产亚洲av| 成人午夜高清在线视频| 欧美不卡视频在线免费观看| 岛国在线免费视频观看| 婷婷亚洲欧美| 性插视频无遮挡在线免费观看| 特级一级黄色大片| 天美传媒精品一区二区| 日韩,欧美,国产一区二区三区 | 日韩欧美三级三区| 国产在线精品亚洲第一网站| 国产一区亚洲一区在线观看| 成人综合一区亚洲| 少妇的逼好多水| 国内精品久久久久精免费| 一进一出好大好爽视频| 久久人妻av系列| 久久精品国产清高在天天线| 真实男女啪啪啪动态图| 美女黄网站色视频| 国产成人a区在线观看| 高清日韩中文字幕在线| 国产 一区 欧美 日韩| 男女做爰动态图高潮gif福利片| 蜜桃亚洲精品一区二区三区| eeuss影院久久| 久久精品国产99精品国产亚洲性色| 午夜福利在线在线| 亚洲一区高清亚洲精品| 搡老岳熟女国产| 成人二区视频| 黄片wwwwww| 极品教师在线视频| 精品国产三级普通话版| 精品人妻一区二区三区麻豆 | 特大巨黑吊av在线直播| 一区福利在线观看| 精品久久久久久久久久久久久| 国产伦精品一区二区三区视频9| 特大巨黑吊av在线直播| 国产精品av视频在线免费观看| 欧美成人精品欧美一级黄| 亚洲欧美日韩无卡精品| 一卡2卡三卡四卡精品乱码亚洲| 国产 一区精品| 久久精品国产99精品国产亚洲性色| 国产午夜精品久久久久久一区二区三区 | 美女黄网站色视频| 成人美女网站在线观看视频| 国产黄a三级三级三级人| 免费人成在线观看视频色| 久久人人爽人人爽人人片va| 看片在线看免费视频| 欧美一区二区国产精品久久精品| 91在线观看av| 日本a在线网址| 一个人看视频在线观看www免费| 一个人免费在线观看电影| av在线亚洲专区| 免费观看人在逋| 久久精品国产亚洲av涩爱 | 精品久久久久久成人av| 菩萨蛮人人尽说江南好唐韦庄 | 国产亚洲精品久久久com| 精品一区二区三区视频在线观看免费| 国产片特级美女逼逼视频| 一进一出抽搐gif免费好疼| 国产精品美女特级片免费视频播放器| 国产av在哪里看| 热99re8久久精品国产| 日韩精品青青久久久久久| 麻豆乱淫一区二区| 精华霜和精华液先用哪个| 午夜福利成人在线免费观看| 男女做爰动态图高潮gif福利片| 久久精品国产鲁丝片午夜精品| 香蕉av资源在线| 国产 一区 欧美 日韩| 国产成年人精品一区二区| 精品国内亚洲2022精品成人| 噜噜噜噜噜久久久久久91| 中文资源天堂在线| 丝袜美腿在线中文| 插逼视频在线观看| 国产淫片久久久久久久久| 日韩人妻高清精品专区| 免费看av在线观看网站| 日韩精品有码人妻一区| 久久人妻av系列| 久久人人爽人人爽人人片va| 97超碰精品成人国产| 日本免费a在线| 少妇高潮的动态图| 欧美成人一区二区免费高清观看| 国产成人精品久久久久久| 午夜爱爱视频在线播放| 久久人人精品亚洲av| 特大巨黑吊av在线直播| 国产成人91sexporn| 又黄又爽又刺激的免费视频.| 欧美精品国产亚洲| 亚洲国产日韩欧美精品在线观看| 女人十人毛片免费观看3o分钟| 免费av不卡在线播放| 久久久久九九精品影院| 国产乱人视频| 91精品国产九色| 精品国产三级普通话版| 久久久久久久亚洲中文字幕| 丰满乱子伦码专区| 国产探花极品一区二区| 性插视频无遮挡在线免费观看| 夜夜看夜夜爽夜夜摸| 亚洲国产精品成人久久小说 | 免费大片18禁| 欧美最新免费一区二区三区| 最近在线观看免费完整版| 又粗又爽又猛毛片免费看| 欧美色视频一区免费| 国产亚洲精品久久久久久毛片| 好男人在线观看高清免费视频| 97碰自拍视频| 男人狂女人下面高潮的视频| 国产成人91sexporn| 欧美一级a爱片免费观看看| 九九热线精品视视频播放| 特级一级黄色大片| av卡一久久| 精品一区二区三区视频在线观看免费| 99热6这里只有精品| 午夜福利18| 高清毛片免费看| 亚洲18禁久久av| 午夜免费男女啪啪视频观看 | 永久网站在线| 97超级碰碰碰精品色视频在线观看| 欧美日韩精品成人综合77777| 久久久久精品国产欧美久久久| 亚洲熟妇中文字幕五十中出| av在线亚洲专区| 毛片一级片免费看久久久久| 亚洲一区二区三区色噜噜| 少妇熟女欧美另类| eeuss影院久久| 丝袜喷水一区| 欧美性猛交黑人性爽| 天天一区二区日本电影三级| 亚洲精品国产av成人精品 | 成年av动漫网址| 深夜精品福利| av在线亚洲专区| 六月丁香七月| 久久热精品热| 欧美极品一区二区三区四区| 波多野结衣高清作品| 18+在线观看网站| 欧美在线一区亚洲| 精品一区二区三区人妻视频| 啦啦啦啦在线视频资源| 国产精品一及| 亚洲成av人片在线播放无| 久久韩国三级中文字幕| 91久久精品电影网| 成年女人毛片免费观看观看9| 99热网站在线观看| 18+在线观看网站| 少妇猛男粗大的猛烈进出视频 | 亚洲国产欧美人成| 国产毛片a区久久久久| 午夜精品在线福利| 国产真实伦视频高清在线观看| 精品人妻一区二区三区麻豆 | 天堂网av新在线| 欧美性感艳星| 在线天堂最新版资源| 精品无人区乱码1区二区| 国产一区二区在线观看日韩| 午夜精品一区二区三区免费看| 18禁在线播放成人免费| 插逼视频在线观看| 乱系列少妇在线播放| www日本黄色视频网| 国产成人a区在线观看| 欧美性感艳星| 欧美日韩一区二区视频在线观看视频在线 | 成人无遮挡网站| 日韩欧美 国产精品| 韩国av在线不卡| 久久99热这里只有精品18| 亚洲欧美日韩高清在线视频| 欧美一区二区亚洲| 你懂的网址亚洲精品在线观看 | 能在线免费观看的黄片| 久久精品国产自在天天线| 国内精品久久久久精免费| 午夜亚洲福利在线播放| 国产极品精品免费视频能看的| 麻豆av噜噜一区二区三区| 久久婷婷人人爽人人干人人爱| 免费一级毛片在线播放高清视频| 日韩av在线大香蕉| 中文字幕av成人在线电影| 99热这里只有精品一区| 国产精品日韩av在线免费观看| 久久久久国产精品人妻aⅴ院| 亚洲欧美日韩高清专用| 精品久久久久久久久亚洲| 国产色婷婷99| 欧美zozozo另类| 99久久久亚洲精品蜜臀av| 特大巨黑吊av在线直播| 国产黄a三级三级三级人| 国产毛片a区久久久久| 亚洲国产色片| 成人国产麻豆网| 国产精品野战在线观看| 嫩草影院入口| av天堂在线播放| 午夜久久久久精精品| av在线天堂中文字幕| 亚洲成人av在线免费| 99热全是精品| 国产精品免费一区二区三区在线| 国产精品亚洲美女久久久| 久久九九热精品免费| 日本色播在线视频| 99精品在免费线老司机午夜| 日韩av在线大香蕉| 欧美+亚洲+日韩+国产| 亚洲图色成人| 听说在线观看完整版免费高清| 在线天堂最新版资源| 性欧美人与动物交配| 亚洲自偷自拍三级| 久久久久国产网址| 能在线免费观看的黄片| 国产蜜桃级精品一区二区三区| 97超碰精品成人国产| 亚洲欧美清纯卡通| 亚洲国产精品sss在线观看| 日韩大尺度精品在线看网址| 亚洲国产精品成人综合色| 国产精品亚洲美女久久久| 午夜精品在线福利| 在线观看一区二区三区| 在线国产一区二区在线| 男人的好看免费观看在线视频| 深夜a级毛片| 久久精品国产亚洲av天美| 久久久久久国产a免费观看| 国产高清不卡午夜福利| 最新在线观看一区二区三区| 国产亚洲91精品色在线| 国产白丝娇喘喷水9色精品| 秋霞在线观看毛片| 色综合亚洲欧美另类图片| 久久久久免费精品人妻一区二区| 日本欧美国产在线视频| www日本黄色视频网| 精品久久久久久久人妻蜜臀av| 亚洲av熟女| 日本免费一区二区三区高清不卡| 99久久成人亚洲精品观看| 亚洲av一区综合| 99久久久亚洲精品蜜臀av| 色尼玛亚洲综合影院| 日韩欧美免费精品| 成人毛片a级毛片在线播放| 国产精品永久免费网站| 欧美性感艳星| 成人av一区二区三区在线看| 国内久久婷婷六月综合欲色啪| 欧美一区二区亚洲| 国产高清视频在线播放一区| 在线播放国产精品三级| 欧美一区二区亚洲| 一进一出抽搐gif免费好疼| 国产乱人偷精品视频| 日韩成人av中文字幕在线观看 | 欧美极品一区二区三区四区| 国产精品爽爽va在线观看网站| 成年女人看的毛片在线观看| 精品久久久久久成人av| 99九九线精品视频在线观看视频| 国产免费一级a男人的天堂| 亚洲国产精品国产精品| 国产黄色小视频在线观看| 大型黄色视频在线免费观看| 真实男女啪啪啪动态图| 欧美在线一区亚洲| 2021天堂中文幕一二区在线观| 99精品在免费线老司机午夜| 最后的刺客免费高清国语| av国产免费在线观看| 欧美高清成人免费视频www|