徐逸帆,曹樂,黃經(jīng)緯,闞秀
上海工程技術(shù)大學(xué)電子電氣工程學(xué)院,上海201620
生物阻抗頻譜(BIS)是一種基于電路的生物組織阻抗測量和監(jiān)測技術(shù),通過分析BIS信號可以知道生物組織的生理和病理狀態(tài)[1]。近年來,BIS 技術(shù)在在早期疾病診斷、生物組織液監(jiān)測、生理狀態(tài)評估等方面應(yīng)用廣泛[2-4]。2003年,Cox-Reijven[5]利用BIS技術(shù)對胃腸道疾病患者實際無脂肪質(zhì)量進行評價,解決了傳統(tǒng)方法評價困難的問題。2009年,Jaffein等[6]通過分析50 kHz下腕踝的BIS特性,對細胞外液進行水分測定,為人體成分分析奠定基礎(chǔ)。2013年,Hornero 等[7]設(shè)計了一種用于評價肌肉及心血管活動的生物阻抗監(jiān)控系統(tǒng),通過對單個細胞進行BIS 分析,研究其他細胞的生理參數(shù)。2016年,何想等[8]基于所采集不同人手掌的BIS數(shù)據(jù),通過提取特征參數(shù)對手掌所屬身份進行識別。
Cole-Cole 生物電阻抗特征方程的建立,為生物阻抗參數(shù)提取和分析提供了有效理論工具,越來越多的學(xué)者關(guān)注Cole 參數(shù)提取方法。2013年,Yang[9]對經(jīng)典最小二乘法(LS)和最小一乘法(LAD)的Cole參數(shù)提取效果進行實驗對比,結(jié)果證明LAD 在處理受噪聲擾動數(shù)據(jù)時具備更好的魯棒性。2015年,Nejadghol 等[10]采用奇異值分解法對BIS 數(shù)據(jù)進行處理,降低了Cole 參數(shù)的提取誤差和標(biāo)準(zhǔn)差,并通過手臂實驗進行了證明。2016年,陳恒等[11]利用模擬退火(SA)算法和阻尼(LS)算法進行聯(lián)合反演求解Cole參數(shù),解決了經(jīng)典LS 算法對于初始值依賴的問題。2018年,郭玥[12]將LAD 算法和SA 算 法 結(jié)合提取Cole參數(shù),規(guī)避了LAD 迭代求導(dǎo)的問題,并對初值選擇誤差和易陷局部最優(yōu)問題進行了分析。雖然現(xiàn)有方法在一定程度上解決了Cole 參數(shù)的提取問題,但是普遍存在對奇異點和隨機噪聲敏感度高、魯棒性差和擬合時間長等缺點。本文采用精英保留策略和LAD 算法對傳統(tǒng)的遺傳算法(GA)進行改進,給出一種魯棒性強、收斂速度快的Cole 參數(shù)提取新方法,并將仿真結(jié)果與現(xiàn)有方法進行對比,驗證所得結(jié)果的實用性與有效性。
研究發(fā)現(xiàn)低頻電流(低于1 MHz)只能流經(jīng)細胞外液,而高頻電流可以穿透細胞膜流經(jīng)細胞內(nèi)外液[13]。根據(jù)生物細胞組織的這種阻抗特性,通常采用如圖1所示的等效電路分析生物電性質(zhì)及其變化規(guī)律,Re是細胞外液電阻,Ri是細胞內(nèi)液電阻,C是細胞膜電容。
圖1 生物體等效電路模型Fig.1 Biological equivalent circuit model
1940~1941年,美國研究學(xué)者Cole等[14]基于生物組織等效電路提出Cole-Cole理論,即生物電阻抗實際測量值的軌跡是復(fù)平面第4象限中的一段圓弧且圓心在第1象限內(nèi)(圖2)。
圖2 生物組織Cole圓弧曲線Fig.2 Biological tissue Cole arc curve
建立如下生物電阻抗特征方程:
其中,Z(ω)為生物復(fù)阻抗,ω為角頻率,R∞為頻率無窮大時的電阻,R0為直流(頻率為0)時的電阻,α為散射系數(shù)(0<a<1),τ為弛豫時間,τ對應(yīng)特征頻率fc:
為深入分析生物電阻抗性質(zhì)及生物組織變化規(guī)律,需要基于測量數(shù)據(jù)求得生物電阻抗等效電路的模型參數(shù)。然而,由于受到測量儀器精度和內(nèi)外界噪聲等因素的影響,實際測量到的BIS數(shù)據(jù)大多分布在標(biāo)準(zhǔn)圓弧軌跡的兩側(cè)[15],這種情況下需要對采集到的BIS 數(shù)據(jù)進行擬合,求取Cole-Cole 曲線對應(yīng)的圓心(x0,y0)和半徑r0,再由幾何關(guān)系進一步確定生物電阻抗模型的特征參數(shù),如下所示:
從式(6)可以看出,測量誤差直接導(dǎo)致弛豫時間τ不確定。實際中通常采用如下式表示的多數(shù)據(jù)點對應(yīng)求均值的方式確定弛豫時間τ:
綜合以上分析,Cole-Cole 模型特征參數(shù)提取的關(guān)鍵是對實測BIS 數(shù)據(jù)的有效擬合?,F(xiàn)有的擬合結(jié)果大多通過LS、LAD 和SA 等方法得到,這些方法存在迭代過程中陷入局部最優(yōu)、前期設(shè)置參數(shù)較多和人機交互依賴性強等缺點。本文提出精英保留策略遺傳算法(EGA)-LAD 彌補現(xiàn)有擬合算法的不足,提取更為精準(zhǔn)的Cole-Cole模型特征參數(shù)。
GA模擬了自然界中生物“優(yōu)勝劣汰”的進化過程[16],利用適應(yīng)度函數(shù)對種群中個體進行評價,淘汰適應(yīng)度低的個體,經(jīng)過隨機選擇、交叉和變異等遺傳操作,實現(xiàn)生物種群的優(yōu)化。但這一過程中,種群中適應(yīng)度最高的個體可能會遭到破壞,這在一定程度上降低了GA全局搜索能力。本文采用精英保留策略對GA進行優(yōu)化,將當(dāng)前代適應(yīng)度最高的個體作為精英個體保留,其余個體進行隨機選擇、交叉和變異等遺傳操作,將精英個體插入到經(jīng)過遺傳操作后的個體集合形成新的種群。EGA-LAD提高了傳統(tǒng)GA的全局收斂性和魯棒性,該算法的步驟包括:種群初始化,設(shè)置適應(yīng)度函數(shù),選擇、交叉和變異,精英保留策略。
2.1.1 種群初始化初始化種群通常以某一分布概率密度函數(shù)隨機生成,為提高算法的反應(yīng)速度和工作效率,本文使用經(jīng)典LS[9]擬合值(x0,y0,r0)作為基礎(chǔ)值,根據(jù)多次實驗數(shù)據(jù)統(tǒng)計計算出合適閾值(hx,hy,hr),確定初始化種群的取值范圍,在取值范圍內(nèi)隨機生成種群。LS 提取Cole 參數(shù)的基本原理是:使得各離散點與擬合圓上對應(yīng)點的徑向誤差的平方和達到最小,對應(yīng)Cole參數(shù)擬合的LS表達式為:
其中,M為實際測量數(shù)據(jù)數(shù)目,xi為激勵頻率fi下所測生物電阻抗的實部,yi為激勵頻率fi下所測生物電阻抗的虛部,(x0,y0)為擬合圓的圓心,r0表示擬合圓的半徑,ei表示徑向誤差即實際測量數(shù)據(jù)點在徑向方向與擬合圓心之間距離。根據(jù)迭代法計算出的LS擬合Cole-Cole 圓的圓心坐標(biāo)和半徑(x0,y0,r0),經(jīng)過多次實驗統(tǒng)計出合適的閾值(hx,hy,hr),確定初始化種群取值范圍hr,y0+hr],再利用均勻分布概率密度函數(shù)在所確定取值范圍內(nèi)隨機生成初始化種群,具體操作如下:
2.1.2 適應(yīng)度函數(shù)在GA 中,通過適應(yīng)度函數(shù)對種群中個體進行評估,選取適應(yīng)度高的優(yōu)秀個體遺傳給下一代,適應(yīng)度函數(shù)的選擇直接影響算法的收斂性能和輸出結(jié)果。下面引入引理1,以便理解下文中適應(yīng)度函數(shù)的定義。
引理1[17]:若存在θ=θ*∈Rn,使目標(biāo)函數(shù)
成立,稱這n個零偏點是“0”點,反之若上式構(gòu)成的方程組不成立,則Q的極小化也不成立。
定義如下適應(yīng)度函數(shù)f(x0,y0,r0):
由引理1,參數(shù)(x0,y0,r0)可由最小化F(x0,y0,r0)求得。
2.1.3 選擇算子本文采用錦標(biāo)賽選擇法對遺傳基因進行選擇。這種選擇策略在一定程度上可以避免過早收斂和停滯現(xiàn)象的產(chǎn)生[15]。錦標(biāo)賽選擇法的具體策略為:從種群中隨機選擇k個個體(每個個體被選擇的概率相同),選擇適應(yīng)度最好的個體作為生成下一代的父體,重復(fù)操作直到新的種群規(guī)模達到原來的種群規(guī)模。
2.1.4 交叉算子本文采用兩點交叉對選出的兩個個體進行交叉操作。兩點交叉是指在兩個配對個體基因串中隨機設(shè)置兩個交叉位,然后再進行交叉位之間基因塊的交換。交叉發(fā)生概率記為Pc,一般Pc范圍設(shè)置為0.7~0.9。對于染色體長度為l的兩個配對個體x1=λ11λ12…λ1l和x2=λ21λ22…λ2l,設(shè)隨機設(shè)置的兩個交叉位為c1、c2,且1 ≤c1<c2≤l- 1,滿足下列公式:
則,生成的新個體為x1' =λ11'λ12'…λ1l' 和x2' =λ21'λ22'…λ2l'。
2.1.5 變異算子生物在遺傳和進化過程中會出現(xiàn)變異,即父代種群中某些基因按一定概率進行改變生成新的個體[18]。這里的一定概率指的是變異概率Pm,由于實際中變異的可能性非常小,一般將Pm范圍設(shè)置為0.01~0.20。本文選用二進制編碼符號串表示個體,利用基本位變異法對個體進行變異操作,即按照概率Pm指定變異位并將該位對應(yīng)的二進制字符取反。對于個體x=λ1λ2…λl,變異運算公式如下:
其中,Δx∈[0, 1]為任意隨機數(shù)。則,變異后生成新的個體為x' =λ1'λ2'…λl'。
2.1.6 精英保留策略將每代種群進化中搜索到的最優(yōu)適應(yīng)度值的個體best保存為精英個體,再對剩下的N-1個個體進行遺傳操作,以此避免目前種群中最好的基因受到丟失和破壞。采用EGA擬合Cole-Cole模型參數(shù)具有如下優(yōu)點:(1)保留傳統(tǒng)GA的可并行性,算法效率高,節(jié)省阻抗參數(shù)提取時間。(2)精英保留策略的引入,提高了傳統(tǒng)GA的搜索速度,全局收斂性更快。
EGA的具體步驟如圖3所示。
圖3 精英保留算法流程圖Fig.3 Elitist preservation algorithm flow chart
Step1,設(shè)置最大迭代次數(shù)L,生成包含N個個體的初始種群P(t)。
Step2,設(shè)置適應(yīng)度函數(shù),計算每一個個體的適應(yīng)度,將適應(yīng)度最高的個體best保存為精英個體。
Step3,判斷迭代次數(shù)t是否達到L次,若滿足執(zhí)行Step7,否則執(zhí)行Step4。
Step4,使用錦標(biāo)賽法對剩下的N-1 個個體進行選擇操作,重新生成規(guī)模為N-1的種群Q(t)。
Step5,對Q(t)進行交叉、變異操作,得到規(guī)模為N-1的種群Q(t)'。
Step6,將精英個體best隨機插入到Q(t)'中,形成得到新一代規(guī)模為N的種群P(t)',令t=t+1,并執(zhí)行Step2。
Step7,輸出計算結(jié)果。
1995年,Rigaud和Hamzaoui等[19]通過實驗得到的一組肌肉阻抗特征參數(shù)R0=150 Ω,R∞=50 Ω,α=0.8,τ=3.0×10-6,從頻率1 kHz~1 MHz 中選取等相角間隔的32 個頻率點( )f1,f2, …,f32,32 個頻率對應(yīng)的即為標(biāo)準(zhǔn)BIS 數(shù)據(jù)點,且標(biāo)準(zhǔn)BIS 數(shù)據(jù)點呈均勻的對數(shù)分布,標(biāo)準(zhǔn)數(shù)據(jù)點如表1 所示。由公式(2)得到特征頻率fc=53.051 6 kHz,選取的BIS數(shù)據(jù)點均勻分布在圓心為(100,-16.246 0)、半徑為r0=52.573 1的Cole圓上。
表1 32個頻率點對應(yīng)的BIS數(shù)據(jù)點Tab.1 BIS data points corresponding to 32 frequency points
在標(biāo)準(zhǔn)BIS數(shù)據(jù)集的基礎(chǔ)上,通過在標(biāo)準(zhǔn)數(shù)據(jù)集中添加奇異值和不同的噪聲生成3組仿真數(shù)據(jù)集。
(1)存在奇異點的仿真數(shù)據(jù)集D1。
在BIS 標(biāo)準(zhǔn)數(shù)據(jù)集頻率f6和f22對應(yīng)的數(shù)據(jù)點處按以下公式增加30%的徑向誤差,其余30 個頻率點的數(shù)據(jù)保持不變:
其中,Psi(xsi,ysi)表示當(dāng)頻率為fi時的標(biāo)準(zhǔn)數(shù)據(jù),表示當(dāng)頻率為fi時的添加噪聲和奇異點的仿真數(shù)據(jù),θi表示頻率fi下橫軸與徑向間夾角,計算公式如下:
(2)存在隨機噪聲的仿真數(shù)據(jù)集D2。
將標(biāo)準(zhǔn)BIS 數(shù)據(jù)中奇數(shù)頻率點處數(shù)據(jù)按下式隨機添加0 至±10%滿足正態(tài)分布的徑向誤差,其余16個偶數(shù)頻率點處數(shù)據(jù)保持不變:
(3)奇異點和隨機噪聲同時存在的仿真數(shù)據(jù)集D3。
將標(biāo)準(zhǔn)數(shù)據(jù)的兩個頻率點f6和f22增加30%的徑向誤差,同時在所有的奇數(shù)頻率點處數(shù)據(jù)隨機添加0至±10%滿足正態(tài)分布的徑向誤差,構(gòu)成仿真數(shù)據(jù)集D3。
為驗證本文方法EGA-LAD 對Cole 參數(shù)提取的有效性和精確性,實驗中選用經(jīng)典LS、SA-LS、SALAD 和EGA-LS 4 種算法分別對D1、D2 和D3 中的BIS 數(shù)據(jù)進行擬合,并與EGA-LAD 的對應(yīng)擬合結(jié)果進行比較。在實驗結(jié)果中,實線代表擬合圓,虛線代表標(biāo)準(zhǔn)圓,eR0、eR∞、eα和efc分別表示Cole參數(shù)R0、R∞、α和fc的擬合誤差。
3.2.1 基于數(shù)據(jù)集D1 的仿真實驗結(jié)果對于仿真數(shù)據(jù)集D1,將兩個偶數(shù)頻率點處標(biāo)準(zhǔn)BIS 數(shù)據(jù)增加30%的徑向誤差。LS、SA-LS、SA-LAD 和EGA-LAD的擬合實驗結(jié)果如圖4 所示,表2 給出BIS 標(biāo)準(zhǔn)數(shù)據(jù)集和4種算法提取的Cole參數(shù)和擬合誤差的對比值。從實驗結(jié)果可以看出,LS算法的擬合效果最差,使用SA-LS 算法的擬合效果明顯改善,但由于LS 算法對數(shù)據(jù)中奇異點敏感度高,使得SA-LS 算法擬合精度低于SA-LAD 算法。本文EGA-LAD 算法對Cole 圓擬合效果最好,Cole 參數(shù)提取精度最高,4 個特征參數(shù)R0、R∞、α和fc擬合誤差分別為0.024 5%、0.118 6%、0.106 3%、0.036 9%。
3.2.2 基于數(shù)據(jù)集D2 的仿真實驗結(jié)果對于仿真數(shù)據(jù)集D2,將標(biāo)準(zhǔn)BIS 奇數(shù)頻率點處隨機添加0至±10%滿足正態(tài)分布的徑向誤差。由于添加徑向誤差的隨機性導(dǎo)致每一次的擬合結(jié)果會有波動,為了更好地評價4種算法對隨機噪聲的魯棒性,本文進行50 次的擬合實驗,最終在表3 中分別列出50 次實驗結(jié)果的平均值進行對比,LS、SA-LS、SA-LAD 和EGA-LAD的擬合實驗結(jié)果如圖5所示。從實驗結(jié)果可以看出,LS算法的擬合效果最差,本文EGA-LAD 算法比兩種改進的SA-LS 和SA-LAD 算法擬合效果好,Cole 參數(shù)提取精度最高,4 個特征參數(shù)R0、R∞、α和fc擬合誤差分別為0.000 7%、0.001 8%、0、0,其中擬合的α和fc參數(shù)與標(biāo)準(zhǔn)值一致,擬合曲線與標(biāo)準(zhǔn)曲線基本重合,實驗結(jié)果說明EGA-LAD 對隨機噪聲的魯棒性最好。
圖4 數(shù)據(jù)集D1實驗結(jié)果Fig.4 Data set D1 experiment results
表2 D1數(shù)據(jù)集下實驗結(jié)果對比Tab.2 Comparison of experimental results on D1 data set
表3 D2數(shù)據(jù)集下實驗結(jié)果對比Tab.3 Comparison of experimental results on D2 data set
圖5 數(shù)據(jù)集D2實驗結(jié)果Fig.5 Data set D2 experiment results
3.2.3 基于數(shù)據(jù)集D3的仿真實驗結(jié)果對于仿真數(shù)據(jù)集D3,將兩個偶數(shù)頻率點處標(biāo)準(zhǔn)BIS數(shù)據(jù)增加30%的徑向誤差,并將奇數(shù)頻率點處隨機添加0至±10%滿足正態(tài)分布的徑向誤差。同樣進行50次的擬合實驗,最終在表4中分別列出50次實驗結(jié)果的平均值進行對比,LS、SA-LS、SA-LAD和EGA-LAD的擬合實驗結(jié)果如圖6所示。從實驗結(jié)果可以看出,相比于與仿真數(shù)據(jù)集D2的擬合結(jié)果,4種算法的擬合精度均有下降,說明奇異點對擬合結(jié)果的影響很大,本文EGA-LAD算法對Cole圓擬合效果依然最好,4個特征參數(shù)R0、R∞、α和fc擬合誤差分別為0.001 6%、0.002 2%、0、0.002 1%,但在數(shù)據(jù)點出現(xiàn)大量偏差時EGA-LAD算法的魯棒性會降低。
表4 D3數(shù)據(jù)集下實驗結(jié)果對比Tab.4 Comparison of experimental results on D3 data set
圖6 數(shù)據(jù)集D3實驗結(jié)果Fig.6 Data set D3 experiment results
表5 給出4種算法對于3組數(shù)據(jù)集的平均擬合時間。綜合對應(yīng)仿真實驗結(jié)論,可知作為經(jīng)典擬合方法的LS 算法用時最短速度最快,但由于其對奇異點的敏感度較高,在奇異點存在的條件下擬合精度會受到嚴重影響,而在BIS 應(yīng)用中,對提取參數(shù)的精度要求優(yōu)先于時間復(fù)雜度要求。改進的SA 改善了魯棒性差的問題,但由于其在每一溫度下均需充分搜索最優(yōu)解,使得擬合時間變長,并在一定程度上影響擬合精度。本文EGA-LAD 算法具有隱式并行性,在保證擬合精度的情況下,縮短了擬合時間,提高了計算效率。
表5 不同算法擬合平均仿真時間Tab.5 Average simulation time of different fitting algorithms
本文綜合考慮精英保留策略和LAD,改進傳統(tǒng)的GA 算法,在生物電阻抗Cole-Cole 特征模型基礎(chǔ)上給出Cole 參數(shù)提取的新方法。通過添加奇異點和隨機噪聲,在標(biāo)準(zhǔn)BIS 數(shù)據(jù)集的基礎(chǔ)上生成3 組仿真實驗數(shù)據(jù)集,利用LS、SA-LS、SA-LAD 和EGA-LAD 4 種算法分別對3 組仿真數(shù)據(jù)集進行仿真擬合,比較了4 種算法的擬合精度和時間復(fù)雜度。實驗結(jié)果表明,本文EGA-LAD 算法具有擬合效果好、對奇異點和隨機噪聲的魯棒性強、Cole 參數(shù)提取精度高的優(yōu)點,驗證了EGA-LAD 算法對生物電阻抗Cole參數(shù)提取的有效性和實用性。