楊宇祥 白世展 林海軍 李建閩 張甫

(湖南師范大學(xué)工程與設(shè)計(jì)學(xué)院,長(zhǎng)沙 410081)

本文從周期信號(hào)的整周期采樣無(wú)頻譜泄露這一原理出發(fā),提出基于multisine 信號(hào)的整周期采樣理論,從理論上推導(dǎo)出滿足multisine 整周期采樣的采樣率設(shè)置條件,構(gòu)建了基于FPGA+數(shù)模轉(zhuǎn)換器+模數(shù)轉(zhuǎn)換器的整周期采樣實(shí)現(xiàn)方法,研制了一種基于multisine 激勵(lì)和整周期采樣的新型多頻電阻抗成像(mfEIT)系統(tǒng);設(shè)計(jì)了胡蘿卜棒+黃瓜棒的雙目標(biāo)成像模型,并進(jìn)行了多頻時(shí)差成像和頻差成像實(shí)驗(yàn).實(shí)驗(yàn)表明,本mfEIT 系統(tǒng)能夠在一個(gè)基波周期(1 ms)內(nèi)實(shí)現(xiàn)20 個(gè)頻率點(diǎn)(2—997 kHz)多目標(biāo)組織邊界的全頻阻抗測(cè)量,成像結(jié)果可區(qū)分具有不同電特性生物組織的結(jié)構(gòu)與位置.本文提出的基于multisine 信號(hào)的整周期采樣理論及其實(shí)現(xiàn)方法,只需一個(gè)multisine 基波周期即可完成一次全頻阻抗測(cè)量,為研制高速mfEIT 系統(tǒng)奠定了理論和技術(shù)基礎(chǔ).

1 引言

電阻抗成像(electrical impedance tomography,EIT)是一種通過(guò)生物組織邊界電阻抗重建其內(nèi)部電導(dǎo)率分布的可視化圖像方法[1].相比于其他成熟的、非侵入性成像技術(shù),如磁共振成像(MRI)[2]和X 射線計(jì)算機(jī)斷層掃描(CT)[3],EIT 具有無(wú)輻射、低成本、小型化等優(yōu)勢(shì),尤其適合連續(xù)和實(shí)時(shí)監(jiān)測(cè)場(chǎng)合,在乳腺癌篩查[4]、肺部急性損傷檢測(cè)[5]、腦部快速電活動(dòng)檢測(cè)[6]及體外血栓檢測(cè)成像[7]等醫(yī)學(xué)臨床檢測(cè)中得到了日益廣泛的應(yīng)用.傳統(tǒng)的EIT系統(tǒng)通?；趩晤l阻抗測(cè)量和時(shí)差成像算法(timedifference EIT,td-EIT)[8],單頻EIT 的一個(gè)基本問(wèn)題是生物組織的絕對(duì)阻抗很難確定,只能觀測(cè)生理或病理過(guò)程引起的相對(duì)阻抗變化[9];td-EIT 可反映心肺引起的體內(nèi)阻抗的前后變化,但在急性卒中治療等特殊的臨床應(yīng)用場(chǎng)合,往往因缺乏中風(fēng)前的參照阻抗而使得td-EIT 成像結(jié)果無(wú)法區(qū)分缺血性卒中和出血性卒中,從而耽誤治療[10].

隨著生物電阻抗譜(bioimpedance spectroscopy,BIS)測(cè)量技術(shù)的進(jìn)步,基于頻率依賴特性的生物組織表征已經(jīng)成為可能[11],多頻電阻抗成像(multi-frequency EIT,mfEIT)技術(shù)和頻差成像算法(frequency-difference EIT,fd-EIT)應(yīng)運(yùn)而生[12].mfEIT 是EIT 與BIS 結(jié)合的技術(shù),根據(jù)阻抗隨頻率的變化特性來(lái)成像,可滿足不同生物組織對(duì)測(cè)量敏感頻率的要求,并實(shí)現(xiàn)對(duì)生物組織電學(xué)特性的精確表征和二維/三維圖像的實(shí)時(shí)重建[13].與td-EIT算法相比,fd-EIT 算法利用同一時(shí)間內(nèi)測(cè)量得到的多頻阻抗進(jìn)行圖像重構(gòu),不需要過(guò)去的參照阻抗,解決了實(shí)際臨床環(huán)境中參考量不可獲得的問(wèn)題[14].此外,新的td-EIT 算法也將多頻阻抗信息融合到時(shí)差成像中,降低了逆問(wèn)題的自由度和病態(tài)性,可獲得增強(qiáng)的重構(gòu)圖像[15].因此,mfEIT 成為EIT 研究的一個(gè)重要轉(zhuǎn)變,世界各地的研究小組已經(jīng)在開發(fā)mfEIT 系統(tǒng)方面做了大量的努力[16].

mfEIT 技術(shù)的基礎(chǔ)是對(duì)生物體多頻阻抗的快速準(zhǔn)確測(cè)量,其測(cè)量方法主要分為掃頻測(cè)量法和多頻同步測(cè)量法[8].掃頻測(cè)量法利用分時(shí)單頻正弦激勵(lì)信號(hào)步進(jìn)掃描獲取不同頻率下的生物電阻抗信息[17].掃頻測(cè)量法實(shí)現(xiàn)簡(jiǎn)單,目前仍是線性時(shí)不變(LTI)假設(shè)條件下測(cè)量靜態(tài)阻抗的主要方法,但完成一次掃頻測(cè)量的時(shí)間相對(duì)較長(zhǎng);而對(duì)于生命體連續(xù)變化的時(shí)變系統(tǒng)(如心血管系統(tǒng)),掃頻測(cè)量法很難準(zhǔn)確獲取這種非平穩(wěn)條件下的動(dòng)態(tài)阻抗,從而丟失重要的診斷信息[18].多頻同步測(cè)量法利用寬頻激勵(lì)信號(hào)一次獲取多個(gè)頻率點(diǎn)的生物電阻抗信息,測(cè)量速度是掃頻測(cè)量法的2.76 倍[19],可準(zhǔn)確記錄生命時(shí)變系統(tǒng)某時(shí)刻的瞬時(shí)阻抗譜信息,因此多頻同步測(cè)量法是mfEIT 的發(fā)展趨勢(shì)[20].

多頻同步激勵(lì)信號(hào)的選擇是實(shí)現(xiàn)mfEIT 系統(tǒng)中多頻阻抗快速測(cè)量的關(guān)鍵,需同時(shí)滿足LTI 系統(tǒng)的線性前提假設(shè)和活體組織的安全標(biāo)準(zhǔn)[21].本文作者前期研究證明,具有稀疏頻譜分布的寬帶信號(hào)(信號(hào)能量較均勻地分布在頻率間隔較大的有限個(gè)頻點(diǎn)上)是多頻同步激勵(lì)信號(hào)的理想選擇[22].Chirp 信號(hào)(線性調(diào)頻脈沖)是當(dāng)前mfEIT 系統(tǒng)經(jīng)常采用的一類多頻同步激勵(lì)信號(hào),Kusche 等[23]開發(fā)的寬頻EIT 系統(tǒng)和天津大學(xué)譚超團(tuán)隊(duì)[24]研發(fā)的敏感帶寬SWEIT 系統(tǒng)均采用Chirp 作為多頻阻抗測(cè)量的激勵(lì)信號(hào).Chirp 信號(hào)的優(yōu)點(diǎn)是具有平坦寬泛的頻譜范圍,使用者可彼此獨(dú)立地選擇激勵(lì)脈沖的頻率范圍和持續(xù)時(shí)間;缺點(diǎn)是頻譜過(guò)于密集且存在多余的高次諧波,每條譜線上的激發(fā)能量往往太低,對(duì)測(cè)量的信噪比(SNR)影響較大[25].

最新研究證明,多頻正弦(multisine)信號(hào)沒(méi)有多余的高次諧波,因此具有更高的能效,對(duì)信號(hào)放大器件的壓擺率要求較低(約0.3 V/μs)[19].西班牙學(xué)者Sanchez 等[26]早前已證明,相比于Chirp激勵(lì)信號(hào),采用multisine 激勵(lì)信號(hào)的生物電阻抗測(cè)量系統(tǒng)具有高出20—30 dB 的SNR.近年來(lái),國(guó)內(nèi)外多個(gè)EIT 團(tuán)隊(duì)陸續(xù)開發(fā)了基于multisine 激勵(lì)的mfEIT 系統(tǒng),如RWTH Aachen 大學(xué)Aguiar等[8]通過(guò)9 個(gè)不同頻率正弦波疊加生成多正弦同步信號(hào),愛(ài)丁堡大學(xué)Yunjie 和Jiabin[13]通過(guò)數(shù)字配置正弦波頻率與相位實(shí)現(xiàn)兩個(gè)頻率離散波形疊加,天津大學(xué)Yan 等[27]使用5 個(gè)單頻正弦求和而合成多諧波波形.然而,上述mfEIT 系統(tǒng)的multisine激勵(lì)信號(hào)均通過(guò)多個(gè)直接數(shù)字合成器(direct digital synthesizer,DDS)產(chǎn)生不同的單頻正弦波再經(jīng)累加器疊加而成,這樣合成的multisine 信號(hào)所能包含的頻率點(diǎn)個(gè)數(shù)嚴(yán)格受到系統(tǒng)DDS 硬件資源的限制,同時(shí)因參與疊加的多個(gè)正弦波相位未經(jīng)優(yōu)化而可能使得合成的multisine 信號(hào)具有過(guò)高的波峰因數(shù)(crest factor,CF),從而對(duì)LTI 系統(tǒng)的線性假設(shè)和活體組織的安全標(biāo)準(zhǔn)造成威脅[19].

本文作者早前提出基于相位迭代優(yōu)化的multisine 合成算法[28],該算法可合成具有業(yè)界最小CF 值的multisine 信號(hào).在此基礎(chǔ)上,本文從周期信號(hào)的整周期采樣無(wú)頻譜泄露這一原理出發(fā),提出基于multisine 信號(hào)的整周期采樣理論,構(gòu)建基于現(xiàn)場(chǎng)可編程門陣列(FPGA)+數(shù)模轉(zhuǎn)換器(DAC)+模數(shù)轉(zhuǎn)換器(ADC)的整周期采樣實(shí)現(xiàn)方法,研制了一種基于multisine 激勵(lì)和整周期采樣技術(shù)的新型mfEIT 系統(tǒng),該系統(tǒng)只需一個(gè)multisine 信號(hào)基波周期的整周期采樣即可完成一次全頻阻抗測(cè)量,可大幅提升mfEIT 的成像速度.論文設(shè)計(jì)了胡蘿卜棒與黃瓜棒的多目標(biāo)阻抗測(cè)量、時(shí)差與頻差成像實(shí)驗(yàn),以驗(yàn)證該mfEIT 系統(tǒng)的有效性.

2 Multisine 整周期采樣理論與多頻阻抗快速測(cè)量原理

2.1 Multisine 整周期采樣理論

設(shè)時(shí)域multisine 信號(hào)x(t)包含M個(gè)諧波分量,可表示為

式中Am,fm,φm分別表示multisine 信號(hào)第m次諧波的幅值、模擬頻率和初相位,m為正整數(shù).用采樣率fs對(duì)x(t)進(jìn)行采樣,得到離散的multisine信號(hào)x(n):

式中N為采樣點(diǎn)數(shù).對(duì)信號(hào)x(n)進(jìn)行N點(diǎn)離散傅里葉變換(DFT)變換得

將(2)式代入(3)式,并利用歐拉公式ejθ=cosθ+j sinθ得

定義基波頻率f0、采樣率fs和采樣點(diǎn)數(shù)N滿足關(guān)系式

由于multisine 各頻率分量滿足諧波分布,即

式中,qm為正整數(shù),表示multisine 信號(hào)所包含的各次諧波分量相對(duì)于基波頻率f0的倍數(shù).將(6)式代入(4)式得

由(7)式可知,當(dāng)采樣點(diǎn)數(shù)N滿足關(guān)系式(5)時(shí),multisine 信號(hào)的頻譜X(k)只在譜線k=qm(m=1,···,M)處才有非零值,頻率分辨率 Δf=f0,即DFT 運(yùn)算后只在qmf0(m=1,···,M) 處具有非零譜線,頻譜無(wú)泄漏.一般地,若采樣率fs滿足:

式中fmax表示multisine 信號(hào)中的最大諧波頻率,則只需一個(gè)基波周期T0=1/f0即可獲得N點(diǎn)的整周期采樣,(8)式就是整周期采樣的條件.為了方便快速傅里葉變換(FFT)運(yùn)算,N一般取1024,2048 等2 的指數(shù)倍數(shù)值.

2.2 基于multisine 整周期采樣的多頻阻抗快速測(cè)量原理

基于整周期采樣的多頻阻抗快速測(cè)量原理如圖1 所示.

圖1 中,FPGA 將存儲(chǔ)在其內(nèi)部ROM 中的一個(gè)基波周期的N點(diǎn)離散multisine 信號(hào)順序循環(huán)讀出,經(jīng)數(shù)模轉(zhuǎn)換器(DAC)后轉(zhuǎn)換成連續(xù)的multisine 信號(hào),并經(jīng)恒流源驅(qū)動(dòng)后變?yōu)殡娏餍盘?hào)i,注入到被測(cè)阻抗ZX中,得到響應(yīng)電壓信號(hào)v.對(duì)i和v進(jìn)行N點(diǎn)整周期采樣,分別得到i(n)和v(n).根據(jù)信號(hào)與系統(tǒng)理論,對(duì)系統(tǒng)輸入信號(hào)i(n)和輸出信號(hào)v(n)進(jìn)行傅里葉變換,即可得系統(tǒng)的頻率響應(yīng)H(ω):

圖1 基于multisine 整周期采樣的多頻阻抗測(cè)量原理圖Fig.1.Schematic diagram of multi-frequency impedance measurement based on integer-period sampling.

式中,I(ω)表示電流激勵(lì)i的傅里葉變換,V(ω)表示電壓響應(yīng)v的傅里葉變換,因此,頻率響應(yīng)H(ω)的物理意義就是被測(cè)阻抗的阻抗譜ZX(ω).

圖1 中,為了實(shí)現(xiàn)整周期采樣,激勵(lì)電流i的生成和信號(hào)的采樣均受FPGA 發(fā)出的整周期采樣同步時(shí)鐘信號(hào)CLK 統(tǒng)一控制,采樣率fs被精確控制為基波頻率f0的N倍,因此在一個(gè)基波周期內(nèi)采樣后得到的離散激勵(lì)電流信號(hào)i(n)和響應(yīng)電壓信號(hào)v(n)都包含N個(gè)點(diǎn),這時(shí)在頻域?qū)(n)和v(n)進(jìn)行N點(diǎn)DFT 運(yùn)算時(shí),就不會(huì)引起頻譜泄漏.令I(lǐng)k,Φk(k=0,1,···,N— 1)分別代表激勵(lì)電流i(n)進(jìn)行N點(diǎn)DFT 運(yùn)算后的幅值譜和相位譜,Vk,Ψk(k=0,1,···,N—1)分別代表響應(yīng)電壓v(n)N點(diǎn)DFT 運(yùn)算后的幅值譜和相位譜.由(7)式可知,基于整周期采樣的DFT 運(yùn)算只在qmf0(m=1,···,M)處具有非零譜線,則被測(cè)阻抗ZX落在multisine激勵(lì)信號(hào)(1)式所包含的M個(gè)諧波頻率點(diǎn)上的多頻阻抗可由下式計(jì)算得到:

式中,Zk,θk分別表示被測(cè)阻抗ZX的幅值和相位.

3 多頻電阻抗成像(mfEIT)系統(tǒng)設(shè)計(jì)

3.1 基于FPGA 的mfEIT 硬件系統(tǒng)設(shè)計(jì)

本文構(gòu)建了基于FPGA 的mfEIT 系統(tǒng),其硬件結(jié)構(gòu)原理如圖2 所示.系統(tǒng)主要包括現(xiàn)場(chǎng)可編程門陣列(FPGA)模塊、DAC 與ADC 模塊以及模擬前端與電極陣列(含恒流源、差分放大電路、電極切換電路等).

圖2 基于FPGA 的mfEIT 系統(tǒng)結(jié)構(gòu)原理圖Fig.2.System structure diagram of the mfEIT system based on FPGA.

圖2 中,一個(gè)基波周期的multisine 信號(hào)被離散化成4096 個(gè)點(diǎn)預(yù)先存儲(chǔ)在FPGA 的ROM 中,DAC 在鎖相環(huán)(PLL)的控制下順序讀取ROM 中的波形值生成模擬multisine 信號(hào),經(jīng)低通濾波器(LPF)濾波后送入恒流源轉(zhuǎn)變?yōu)閙ultisine 電流信號(hào)i,通過(guò)模擬開關(guān)選擇一對(duì)相鄰電極作為電流激勵(lì)電極,將i注入到被測(cè)生物邊界阻抗ZX中,并通過(guò)參考串聯(lián)電阻Rf返回恒流源;同時(shí),通過(guò)模擬開關(guān)選擇除電流激勵(lì)電極以外的一對(duì)電極作為電壓采集電極,測(cè)量經(jīng)差分放大后得到的響應(yīng)電壓信號(hào)v,并同步測(cè)量Rf上的壓降經(jīng)差分放大后得到的電流信號(hào)i;PLL 同步控制兩路ADC 分別對(duì)i和v進(jìn)行同步整周期采樣,分別得到離散序列i(n)和v(n)并緩存在FPGA 片內(nèi)資源配置的兩路先入先出(FIFO)隊(duì)列中;在FPGA 上開辟FFT 運(yùn)算單元,分別對(duì)i(n)和v(n)進(jìn)行4096 點(diǎn)FFT運(yùn)算,分別得到對(duì)應(yīng)的傅里葉系數(shù)(即激勵(lì)信號(hào)的幅值譜Ik、相位譜φk,響應(yīng)信號(hào)的幅值譜Vk、相位譜Ψk),并根據(jù)(10)式計(jì)算得到本次電極配置下生物邊界阻抗ZX的幅值譜與相位譜信息,即完成一次全頻阻抗測(cè)量;依次切換電壓采集電極,對(duì)于一個(gè)N電極的EIT 系統(tǒng)來(lái)說(shuō),一種激勵(lì)模式下共需進(jìn)行(N— 3)次全頻阻抗測(cè)量.隨后,切換激勵(lì)電極,重復(fù)上述電壓采集過(guò)程;N個(gè)電極可切換產(chǎn)生N種激勵(lì)模式,共需進(jìn)行 (N— 3) ×N個(gè)測(cè)量通道的全頻阻抗測(cè)量.在獲得各種電極配置下的待測(cè)場(chǎng)邊界阻抗譜數(shù)據(jù)后,最終根據(jù)成像算法計(jì)算待測(cè)場(chǎng)內(nèi)部各單元電導(dǎo)率分布,并重構(gòu)測(cè)量對(duì)象時(shí)差與頻差圖像.

mfEIT 系統(tǒng)實(shí)物如圖3 所示,其中FPGA 平臺(tái)選用火龍果(Red Pitaya)125-14Starter Kit FPGA開發(fā)套件,該套件搭載了Xilinx 公司FPGA Zynq-7010,并集成了14 位125Msps 的雙路同步ADC及DAC;設(shè)計(jì)了一個(gè)包含16 電極的圓柱體形水槽作為成像模型,該成像模型通過(guò)電極陣列切換共可構(gòu)成 (16 — 3) × 16=208 種基于4 電極法的阻抗測(cè)量通道配置,FPGA 平臺(tái)控制電極陣列依次完成208 個(gè)測(cè)量通道在各個(gè)頻率下的邊界阻抗,進(jìn)而對(duì)水槽進(jìn)行網(wǎng)格剖析,計(jì)算各單元電導(dǎo)率的變化并重構(gòu)圖像.

圖3 mfEIT 系統(tǒng)實(shí)物圖Fig.3.Photo of the mfEIT system.

3.2 Multisine 信號(hào)的合成

如(1)式所示,multisine 是一種通過(guò)有限個(gè)具有不同幅值、頻率、相位的正弦波疊加合成的周期信號(hào),若各正弦波分量的相位隨機(jī)設(shè)置,則合成的multisine 信號(hào)往往具有較高的波峰因數(shù)(CF,峰值/有效值),這對(duì)于要求盡量壓低激勵(lì)信號(hào)峰值從而使被測(cè)體保持線性的生物電阻抗測(cè)量很不利.在峰值一定的情況下,擁有較高CF 的激勵(lì)信號(hào)意味著較低的能量注入被測(cè)體;反之,激勵(lì)信號(hào)擁有較低的CF 意味著可以提供更多的能量從而提高測(cè)量精度[29].因此,為了使multisine 激勵(lì)信號(hào)在測(cè)量生物電阻抗時(shí)擁有最大的測(cè)量精度,通過(guò)合適的相位組合優(yōu)化以最小化其CF 值是必然選擇.關(guān)于multisine 的相位優(yōu)化算法,國(guó)內(nèi)外學(xué)者已經(jīng)研究很多年,先后發(fā)展出解析法[30]和迭代法[31]兩大類別,但各有其優(yōu)缺點(diǎn).

本文作者綜合解析法和迭代法的優(yōu)點(diǎn),提出一種改進(jìn)相位迭代優(yōu)化的multisine 合成算法[28],該算法可以按用戶所需合成任意頻譜分布的multisine信號(hào),并具有業(yè)界最小的CF 值,可顯著提高阻抗測(cè)量系統(tǒng)的精度.本文利用此算法合成了一種包含20 個(gè)等幅值質(zhì)數(shù)偽對(duì)數(shù)頻譜分布的multisine信號(hào),其各個(gè)諧波分量對(duì)應(yīng)的頻率、初相位如表1所列,其中各個(gè)頻率分量的歸一化幅值均為0.3162,基波周期f0=1 kHz.multisine 信號(hào)一個(gè)完整周期的時(shí)域波形和頻譜分布分別如圖4(a)和圖4(b)所示.此multisine 信號(hào)通過(guò)偽對(duì)數(shù)頻譜分布保證寬頻測(cè)量范圍,而質(zhì)數(shù)頻點(diǎn)則保證各頻點(diǎn)之間無(wú)諧波關(guān)系,從而減小非線性諧波畸變對(duì)阻抗譜測(cè)量的影響.

圖4 Multisine 信號(hào) (a)時(shí)域波形;(b)頻譜分布Fig.4.Multisine signal:(a) Time domain waveform;(b) spectrum distribution.

表1 合成的包含20 個(gè)等幅值偽對(duì)數(shù)頻譜分布的multisine 信號(hào)頻率、相位Table 1.Frequencies and phases of the synthesized multisine signal with equivalent amplitude and pseudo-logarithmic spectral distribution.

3.3 系統(tǒng)信噪比分析

基于本文合成的multisine 激勵(lì),對(duì)所研制的mfEIT 系統(tǒng)在不同頻率下進(jìn)行了信噪比(SNR)評(píng)估.測(cè)量選用上述水槽均質(zhì)場(chǎng)模型,采用電導(dǎo)率為0.07 S/m 的鹽水作為傳感對(duì)象.采用相鄰電極法,對(duì)生成一幀圖像所需的208 個(gè)通道進(jìn)行阻抗譜測(cè)量,通過(guò)下式評(píng)估各頻率點(diǎn)的SNR:

式中,Zk為通道在頻率fk時(shí)對(duì)應(yīng)的阻抗測(cè)量平均幅值;為頻率fk時(shí)對(duì)應(yīng)的阻抗測(cè)量模值方差,分別定義為

式中,P為測(cè)量次數(shù),每個(gè)通道進(jìn)行50 次阻抗譜重復(fù)測(cè)量;分別對(duì)應(yīng)于輸入電流激勵(lì)和輸出電壓響應(yīng)在頻率fk時(shí)對(duì)應(yīng)的測(cè)量信號(hào)幅值平均值與方差.

各頻率點(diǎn)SNR 的估計(jì)值如表2 所列,結(jié)果表明,各頻率點(diǎn)的SNR 較為均衡,平均SNR 為55.3 dB,平均標(biāo)準(zhǔn)差為 ± 6.2 dB.由于鹽水的電導(dǎo)率隨測(cè)量頻率理論上近似不變,因此高頻與低頻段的SNR 相差不大,而2 kHz 有明顯變小的主要原因可能是低頻電極接觸阻抗影響較大,導(dǎo)致測(cè)量阻抗的準(zhǔn)確性下降.

表2 Multisine 信號(hào)20 個(gè)頻率點(diǎn)的通道信噪比平均值及標(biāo)準(zhǔn)差Table 2.Average and standard deviation of the channel SNR at 20 frequency points of the multisine signal.

4 mfEIT 成像實(shí)驗(yàn)與驗(yàn)證

4.1 mfEIT 成像方法

待測(cè)場(chǎng)生物組織可視為離子導(dǎo)電體,內(nèi)部單元電導(dǎo)率σ(x,y)和可測(cè)參數(shù)邊界電極間電位φ的函數(shù)關(guān)系為

式中?Ω為場(chǎng)域邊界,n表示場(chǎng)域Ω的外法向單位向量,j表示流入場(chǎng)域Ω的激勵(lì)電流密度.

為了計(jì)算邊界阻抗ΔZ與場(chǎng)域Ω內(nèi)部各單元Δσ(x,y)分布之間的關(guān)系,利用聯(lián)合仿真軟件(COMSOL Multiphysics 5.3a with MATLAB)仿真正演模型.為更接近于真實(shí)實(shí)驗(yàn)效果,構(gòu)建實(shí)驗(yàn)水槽二維仿真模型(直徑130 mm),16 個(gè)弧形電極傳感器寬度與相互之間間隔呈1∶1.5 比例均勻圍繞模型一圈,同時(shí)設(shè)置仿真模型內(nèi)溶液電導(dǎo)率(液體的電導(dǎo)率參數(shù)設(shè)置為0.01 S/m).在研究穩(wěn)態(tài)物理電流場(chǎng)中利用有限元網(wǎng)格剖分分析模型場(chǎng)域內(nèi)電位分布情況,計(jì)算得到場(chǎng)域內(nèi)各節(jié)點(diǎn)電位與電導(dǎo)率之間關(guān)系,用靈敏度矩陣J近似表示為

1)時(shí)差成像:待測(cè)場(chǎng)生物組織阻抗Zf可用mfEIT 系統(tǒng)測(cè)量得到:

式中fm表示第m個(gè)諧波頻率分量.當(dāng)邊界電壓產(chǎn)生波動(dòng)時(shí),時(shí)差成像中用起始時(shí)刻t0從均質(zhì)場(chǎng)獲取邊界阻抗信息作為參考數(shù)據(jù),ti時(shí)刻同一頻率阻抗變化ΔZf(t)用于時(shí)差圖像重建,ΔZf(t)計(jì)算公式為

2)頻差成像:與時(shí)差成像方法相比,頻差成像不需要過(guò)去時(shí)刻的阻抗參考值,可解決臨床環(huán)境中無(wú)法獲取參考時(shí)刻邊界阻抗的問(wèn)題.以任意時(shí)刻ti的不同頻率下阻抗數(shù)據(jù)變化ΔZt(f)進(jìn)行圖像重建,ΔZt(f)計(jì)算公式為

3)圖像重建:利用邊界阻抗變化ΔZ[ΔZf(t),ΔZt(f)]計(jì)算內(nèi)部各單元電導(dǎo)率變化Δσ(x,y),重構(gòu)出目標(biāo)組織結(jié)構(gòu).利用Tikhonov-Noser 組合正則化算法優(yōu)化計(jì)算目標(biāo)場(chǎng)內(nèi)部各單元電導(dǎo)率的變化:

式中,εT是Tikhonov 正則化參數(shù),εN是Noser 正則化參數(shù),E是與JTJ維度相同的單位矩陣.

EIT 逆問(wèn)題的求解是一個(gè)非線性病態(tài)問(wèn)題,解存在嚴(yán)重的不穩(wěn)定性.本文選擇Tiknonov-Noser組合正則化算法作為圖像重構(gòu)的算法,目標(biāo)為改善其非線性的病態(tài)問(wèn)題,變成最小化目標(biāo)函數(shù)的求解,這樣可以得到存在且唯一的穩(wěn)定解.選用Tiknonov-Noser 正則化算法對(duì)EIT 逆問(wèn)題的求解,既可以提供正確的成像目標(biāo)位置又可以很好地去除噪聲干擾[32,33].正則化參數(shù)是一個(gè)需要人為設(shè)置的經(jīng)驗(yàn)常數(shù),對(duì)重構(gòu)圖像的質(zhì)量至關(guān)重要,經(jīng)實(shí)驗(yàn)驗(yàn)證正則化參數(shù)εT=1×10-6,εN=100 是比較合適的值.

4)圖像質(zhì)量評(píng)價(jià):為了定量評(píng)估重建圖像的準(zhǔn)確性,采用歸一化的實(shí)驗(yàn)?zāi)Ｐ碗妼?dǎo)率變化與仿真模型電導(dǎo)率變化之間的相關(guān)系數(shù)(CC)作為圖像質(zhì)量的評(píng)價(jià)指標(biāo),定義如下:

式中,Δσi,ΔσTrue分別表示第i個(gè)元素的重建電導(dǎo)率和仿真電導(dǎo)率的變化值,分別表示重建電導(dǎo)率和仿真電導(dǎo)率的變化平均值,n為根據(jù)COMSOL 有限元剖分方法得到的網(wǎng)格節(jié)點(diǎn)數(shù),本文中n=3360.

CC 與電導(dǎo)率的具體值無(wú)關(guān),僅與電導(dǎo)率的空間分布有關(guān).因被測(cè)目標(biāo)的電導(dǎo)率隨頻率變化,在不同基礎(chǔ)參考頻率下選擇不同頻率間隔可以獲得不同的圖像,頻差成像的質(zhì)量難以評(píng)價(jià)[34],故本文僅對(duì)時(shí)差圖像進(jìn)行質(zhì)量評(píng)價(jià).

4.2 mfEIT 系統(tǒng)成像實(shí)驗(yàn)

4.2.1 時(shí)差成像

在成像實(shí)驗(yàn)開始前,先測(cè)量圓柱體形水槽內(nèi)只有鹽水時(shí)的均質(zhì)場(chǎng)邊界阻抗譜數(shù)據(jù),作為時(shí)差成像中的參考數(shù)據(jù)Zf(t0).均質(zhì)場(chǎng)邊界阻抗譜如圖5 所示.由圖5 可見(jiàn),不同頻率下的均質(zhì)場(chǎng)邊界阻抗譜變化規(guī)律基本一致,只是隨著頻率增加,溶液的邊界阻抗譜降低,邊界阻抗從低頻(2 kHz)時(shí)的150 Ω降至高頻(997 kHz)時(shí)的不足50 Ω.

圖5 均質(zhì)場(chǎng)邊界阻抗譜Fig.5.Boundary impedance spectroscopy of homogeneous field.

在測(cè)得均質(zhì)場(chǎng)邊界阻抗譜數(shù)據(jù)基礎(chǔ)上,將胡蘿卜棒和黃瓜棒放入水槽(如圖3 所示),利用(17)式和(19)式得到多頻阻抗時(shí)差圖像,圖6 中依次顯示了胡蘿卜棒和黃瓜棒在20 個(gè)頻率激勵(lì)下(2—997 kHz)的時(shí)差圖像,圖像下方的小數(shù)值分別代表各個(gè)頻率下圖像重建結(jié)果與仿真模型相關(guān)度(即(20)式表示的圖像評(píng)價(jià)指標(biāo)CC),CC 越接近于1 則表示成像結(jié)果相關(guān)度越高,成像質(zhì)量越好.圖6 中,圖像質(zhì)量最好的頻率范圍是73—373 kHz(CC ＞ 0.700).由圖6 可見(jiàn),胡蘿卜棒和黃瓜棒的電導(dǎo)率變化值由低于溶液電導(dǎo)率的變化逐漸升高,成像顏色由藍(lán)色漸變到紅色的變化.其中胡蘿卜棒在13—53 kHz,黃瓜棒在513 kHz,電導(dǎo)率變化與溶液電導(dǎo)率變化相近,故其在此頻率下目標(biāo)圖像并不明顯甚至“消失”.而在低頻和高頻區(qū),由于蘿卜棒和黃瓜棒與溶液電導(dǎo)率差值的影響,靠近邊界的靈敏度高于中心區(qū)域靈敏度,重構(gòu)圖像的兩個(gè)目標(biāo)之間產(chǎn)生了陰影.

圖6 mfEIT 系統(tǒng)的時(shí)差成像Fig.6.Time difference images of the mfEIT system.

4.2.2 頻差成像

頻差成像選擇同時(shí)刻2 kHz 激勵(lì)下的邊界阻抗數(shù)據(jù)作為基礎(chǔ)參考數(shù)據(jù),利用(18)式和(19)式得到同時(shí)刻的多頻頻差阻抗圖像,圖7 是胡蘿卜棒和黃瓜棒在multisine 激勵(lì)(2 —997 kHz)下的多目標(biāo)阻抗頻差圖像.圖7 中,在激勵(lì)頻率小于269 kHz 時(shí),胡蘿卜棒和黃瓜棒的電導(dǎo)率變化明顯均高于溶液,且黃瓜棒的電導(dǎo)率變化高于胡蘿卜棒,因此二者的成像顏色不同.

圖7 mfEIT 系統(tǒng)的頻差成像Fig.7.Frequency difference images of the mfEIT system.

從上述時(shí)差和頻差成像實(shí)驗(yàn)的結(jié)果來(lái)看,本文所設(shè)計(jì)的mfEIT 系統(tǒng)能夠?qū)Χ嗄繕?biāo)生物組織進(jìn)行檢測(cè)與成像,并可通過(guò)對(duì)目標(biāo)阻抗譜與多頻圖像重建結(jié)果進(jìn)行分析,以識(shí)別和區(qū)分不同的生物組織,突出不同目標(biāo)生物組織在寬頻率范圍下的阻抗變化特點(diǎn)和成像變化趨勢(shì).

5 結(jié)論

本文從周期信號(hào)的整周期采樣無(wú)頻譜泄露這一原理出發(fā),提出基于multisine 信號(hào)的整周期采樣理論,首次從理論上推導(dǎo)出滿足multisine 整周期采樣的采樣率設(shè)置條件;構(gòu)建了基于FPGA+DAC+ADC 的整周期采樣實(shí)現(xiàn)方法,研制了一種基于multisine 激勵(lì)和整周期采樣的新型mfEIT系統(tǒng);設(shè)計(jì)了胡蘿卜棒+黃瓜棒的雙目標(biāo)成像模型,并進(jìn)行了多頻時(shí)差成像和頻差成像實(shí)驗(yàn),實(shí)驗(yàn)表明,本mfEIT 系統(tǒng)能夠在一個(gè)基波周期(1 ms)內(nèi)完成一次包含20 個(gè)頻率點(diǎn)(2—997 kHz)的全頻阻抗測(cè)量,成像結(jié)果可區(qū)分具有不同導(dǎo)電特性生物組織的結(jié)構(gòu)與位置.本文提出的基于multisine 信號(hào)的整周期采樣理論及其實(shí)現(xiàn)方法,只需一個(gè)multisine 基波周期即可完成一次全頻阻抗測(cè)量,為研制高速多頻EIT 系統(tǒng)奠定了理論和技術(shù)基礎(chǔ).下一步將通過(guò)提高基波頻率大幅提升單次全頻阻抗的測(cè)量速度,并通過(guò)設(shè)計(jì)并行阻抗測(cè)量方式大幅縮減EIT 多通道阻抗測(cè)量的時(shí)間,實(shí)現(xiàn)基于FPGA 的高速mfEIT 系統(tǒng),并將其應(yīng)用于生命體的動(dòng)態(tài)實(shí)時(shí)成像,如肺通氣監(jiān)測(cè)、心搏血量檢測(cè)等.