基于差分進(jìn)化算法的核磁共振T_2譜多指數(shù)反演

1、第57卷第9期2008年9月100023290200857(09)5956206物理學(xué)報(bào)ACTAPHYSICASINICAVol.57,No.9,September,20082008Chin.Phys.Soc.基于差分進(jìn)化算法的核磁共振T2譜多指數(shù)反演潘克家1)󰂍3陳華譚永基2)1)1)(復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院,上海200433)2)(中國石油大學(xué)(華東),東營257061)(2007年11月15日收到;2008年1月19日收到修改稿)提出了一種基于差分進(jìn)化(DE)算法的核磁共振弛豫信號(hào)多指數(shù)反演新方法.將核磁共振T2譜反演問題轉(zhuǎn)化為帶非負(fù)約束的非線性優(yōu)化問題,不需要預(yù)先給定橫向

2、弛豫時(shí)間T2分布,直接利用差分進(jìn)化算法進(jìn)行反演計(jì)算.在測量信號(hào)低信噪比的情況下,計(jì)算機(jī)模擬和實(shí)驗(yàn)數(shù)據(jù)反演都表明了該方法在分析處理NMR弛豫信號(hào)中的有效性.關(guān)鍵詞:核磁共振,多指數(shù)反演,差分進(jìn)化,巖心分析PACC:7660,受控參數(shù)少,無需設(shè)置初值和進(jìn)11引自旋回波串言1,2行導(dǎo)數(shù)計(jì)算,且易于理解和實(shí)現(xiàn),是一種比較有發(fā)展前途的算法.核磁共振弛豫信號(hào)多指數(shù)反演問題,可以轉(zhuǎn)化為一個(gè)帶非負(fù)約束的非線性優(yōu)化問題,然后利用差分進(jìn)化算法求解.在核磁共振測井中,一般采用CMPG方法測量,縱向弛豫時(shí)間T1和橫向弛豫時(shí)間T2都是用來描述自旋回波信號(hào)的弛豫特征.由于和橫向弛豫時(shí)間T2相對應(yīng)的T2譜能夠提供許多巖石物

3、性和流體特性的信息,越來越受到人們的關(guān)37注.而T2譜的求解本質(zhì)上是一個(gè)多指數(shù)反演的問題,對于多指數(shù)反演,近年來先后出現(xiàn)了罰函數(shù)法(BRD)、奇異值分解法(SVD)、非負(fù)最小二乘(NNLS)以及聯(lián)合迭代重建算法(SIRT),并且在罰函數(shù)法和奇異值分解法的基礎(chǔ)上又發(fā)展了一些改進(jìn)算法,這些算法從不同角度來解決多指數(shù)反演問10題.但這些算法除了受到橫向弛豫時(shí)間布點(diǎn)數(shù)、原始回波采集個(gè)數(shù)等因素影響外,還易受信噪比的影響,在不做任何參數(shù)校正的情況下,這些算法都要求信噪比SNR>401對于低信噪比(SNR<30)情況,上述算法會(huì)出現(xiàn)基線偏離的情況,導(dǎo)致孔隙度11估計(jì)不準(zhǔn),另外,在處理T2譜的非負(fù)

4、約束條件時(shí),它們均采用迭代法修正,大大增加了計(jì)算時(shí)間.12,13差分進(jìn)化(differentialevolution,DE)算法是由Store和Price于1996年共同提出的一種采用浮點(diǎn)矢量編碼在連續(xù)空間中進(jìn)行隨機(jī)搜索的優(yōu)化算法,3國家自然科學(xué)基金(批準(zhǔn)號(hào):10431030)資助的課題.󰂍E2mail:kjpan8,921核磁共振T2弛豫信號(hào)數(shù)學(xué)描述根據(jù)核磁共振理論,核磁共振信號(hào)與橫向弛豫時(shí)間譜滿足第一類Fredholm積分方程:y(t)=Tmax2Tmin2f(T2)e-T2dT2+(t),(1)其中橫向弛豫時(shí)間譜f(T2)為橫向弛豫時(shí)間T2對應(yīng)的幅度值,從物理意義上講,f

5、(T2)為非負(fù)值.T2min與T2max是測量的衰減信號(hào)所能分辨的最短和最長弛豫時(shí)間,一般T2min=011ms,T2max=10000ms,(t)為隨機(jī)噪聲.把(1)式離散化后,所測的核磁共振信號(hào)y(ti)可以認(rèn)為是由一系列滿足單指數(shù)衰減規(guī)律的回波信號(hào)疊加,即mty(ti)=j=1fje-T2j+(ti),1in,(2)其中m,n分別為測量到的回波個(gè)數(shù)和橫向弛豫時(shí)間T2分布個(gè)數(shù),ti為采集時(shí)刻(通常為回波間隔TE的整數(shù)倍),f(ti)為ti時(shí)刻測得的回波幅度.9期潘克家等:基于差分進(jìn)化算法的核磁共振T2譜多指數(shù)反演5957(6)反演T2譜.31多指數(shù)反演優(yōu)化模型41差分進(jìn)化算法橫向弛豫時(shí)間T

6、2分布跨度大,通常范圍從011ms到10000ms(因儀器不同而有所差異),導(dǎo)致方程(2)的求解是一個(gè)典型的不適定問題.測量數(shù)據(jù)很小的誤差也會(huì)導(dǎo)致最后的計(jì)算結(jié)果偏差很大,這也增大了反演過程的難度.現(xiàn)有的各種T2譜反演方法1416,T2分布都是事先采取某種對數(shù)均勻布點(diǎn)得到,然后由方程(2)得到一個(gè)線性方程組,再利用SVD,SIRT,NNLS等方法求解.而對f(T2)非負(fù)約束的處理,都是采用某種迭代修正的思想.由于T2分布事先給定,導(dǎo)致這些方法都存在一個(gè)普遍的問題:就是在測量T2組分離散且分布較寬時(shí),反演得到的T2譜分辨率較低.并且,理非負(fù)約束的方法,求解使用了Levenberg2Marquard

7、t方法,對初值的依賴性比較強(qiáng),可能陷入局部極值.為了避免這些問題,方程(2)的求解轉(zhuǎn)化為一個(gè)帶非負(fù)約束的優(yōu)化問題,同時(shí)反演出T2和f(T2),得到橫向弛豫時(shí)間T2和f(T2)的關(guān)系圖,即T2譜圖.考慮到原問題的病態(tài)性質(zhì)和測量的誤差,通過18,19如下兩種正則化方法,將原問題轉(zhuǎn)化為如下帶非負(fù)約束的優(yōu)化問題.范數(shù)平滑nmt差分進(jìn)化算法是一種基于實(shí)數(shù)編碼的具有保優(yōu)思想的全局優(yōu)化進(jìn)化算法.在1996年舉行的第一屆國際IEEE進(jìn)化優(yōu)化競賽上,對提出的各種方法進(jìn)行了現(xiàn)場驗(yàn)證,DE被證明是最快的進(jìn)化算法.與傳統(tǒng)的進(jìn)化算法如遺傳算法等相較,除了具有收斂速度快,控制變量少,更易于理解和編程實(shí)現(xiàn)等特點(diǎn)之外,其主要

8、區(qū)別在于:1)傳統(tǒng)的GA算法采取二進(jìn)制編碼,而DE算法采取實(shí)數(shù)編碼;2)在GA中,)傳統(tǒng)GA產(chǎn)生的子個(gè)體通過某種概率選擇直接取代父個(gè)體,而DE中新產(chǎn)生的子個(gè)體,只有當(dāng)它比相應(yīng)的父代個(gè)體優(yōu)良時(shí)才替換父個(gè)體.4111基本差分進(jìn)化算法為求非線性函數(shù)f(x1,x2,xD)的最小值,進(jìn)化過程中的第G代利用NP個(gè)D維參數(shù)向量構(gòu)成種群xi,G,i=1,2,NP,其中種群大小NP在整個(gè)進(jìn)化過程中保持不變.411111種群初始化在D維空間,隨機(jī)產(chǎn)生滿足約束條件的NP個(gè)向量xji,0=randj(bj,U-bj,L)+bj,L構(gòu)成初始種群.其中,randj為區(qū)間(0,1)上的隨機(jī)數(shù),bj,U和bj,L分別為第j個(gè)

9、變量xj的上界和下界.411121變異mini=1yi-j=1fje-T2jm+j=1fj,(3)(4)s.t.fj0,1jm.曲率平滑nmimini=1ym-1j=2-j=1fje-tT2j對每一個(gè)目標(biāo)向量xi,G,在1,2,NP范圍內(nèi)隨機(jī)選擇3個(gè)互異整數(shù)r1,r2,r3,且使得rji(1(5)(6)+fj+1-2fj+fj-1,s.t.fj0,1jm.j3),得到變異向量vi,G+1=xr1,G+F(xr2,G-xr3,G),其中,變異因子F(0,2)為一常數(shù),控制偏差向量的放大程度.411131交叉交叉是為了增加種群的多樣性,通過如下方式得到試驗(yàn)向量:vji,G+1,randjCR或j=

10、rnbri,uji,G+1=xji,G+1,randj>CR和jrnbri,其中,為正則化參數(shù),它的選取與原始數(shù)據(jù)的誤差水平相關(guān).對上述非線性優(yōu)化問題,采取差分進(jìn)化算法求解.計(jì)算結(jié)果表明,基于曲率光滑的正則化優(yōu)化算法能更好地抑制f(T2)的振蕩,對正則化參數(shù)的選取不如范數(shù)平滑時(shí)敏感,對此類多指數(shù)反演問題更加有效.和文獻(xiàn)17采取范數(shù)平滑不同,最后的反演實(shí)例都采用曲率平滑,求解優(yōu)化問題(5),5958物理學(xué)報(bào)57卷j=1,2,D.2)變異因子.在理論上也沒有一個(gè)最佳的選式中randj為(0,1)之間的隨機(jī)數(shù),rnbrj1,2,D為隨機(jī)選擇的序列,用它來確保ui,G+1從vi,G+1至少獲取一

11、個(gè)分量.交叉因子CR0,1.411141選擇擇準(zhǔn)則,比較好的容許區(qū)間為014,1,不落在此區(qū)間的變異因子只在極少的情況下有效.通常取初始值F=015,若種群過早收斂,則增大F或NP.3)交叉因子.通常DE算法如果收斂,則交叉因子CR越大,收斂速度越快.故為得到較快的收斂速度,一般先選取CR=110或者CR=019.4)加權(quán)因子.加權(quán)因子(0,1),通常取作018.差分進(jìn)化的選擇方案基于局部競爭機(jī)理,如果ui,G+1的目標(biāo)函數(shù)值小于xi,G+1的目標(biāo)函數(shù)值,則令xi,G+1等于ui,G+1;反之,令xi,G+1等于xi,G.411151邊界條件處理在有邊界約束的問題中,必須確保產(chǎn)生新個(gè)體的參數(shù)值

12、落在問題的可行域中,一個(gè)簡單方法是將不符合邊界約束的新個(gè)體用在可行域中隨機(jī)產(chǎn)生的參數(shù)向量代替.4121差分進(jìn)化算法其他形式51反演數(shù)值結(jié)果這一節(jié),采用差分進(jìn)化算法反演核磁共振T2譜分布.對反演結(jié)DE反演算,收斂速度比較快.第二小節(jié)給出了運(yùn)用差分進(jìn)化算法反演T2譜的一個(gè)實(shí)例.5111理論模型試算結(jié)果與分析上面介紹的是最基本的,用中還發(fā)展了DE,xy加以區(qū)分,其中x“隨機(jī)的”或“最佳的”;y是所利用的差向量的個(gè)數(shù).利用這個(gè)表示方法,前面敘述的基本DE策略可描述為DErand11其他多種形式,主要體現(xiàn)在變異操作的不同,如:1)DEbest1vi,G+1=xbest,G+F(xr1,G-xr2,G).

13、在這一節(jié)的計(jì)算實(shí)例中,先給定一個(gè)具有雙峰特性的T2譜分布,由(2)式計(jì)算出不同時(shí)刻ti的信號(hào)強(qiáng)度值yi,再對y加上一定大小SNR的高斯白噪音得到󰁥y,然后由數(shù)據(jù)(ti,󰁥yi)用DE算法反演T2分布.信噪比SNR按如下方式定義:SNR=20log102)DEbest2vi,G+1=xbest,G+F(xr1,G-xr2,G+xr3,G-xr4,G).3)DErand2vi,G+1=xr1,Gy.y-󰁥y2+F(xr2,G-xr3,G+xr4,G-xr5,G).當(dāng)T2組分m=10時(shí),構(gòu)造了512個(gè)測量數(shù)據(jù)進(jìn)行反演;當(dāng)T2組分m=20時(shí),我們構(gòu)造了1

14、024個(gè)測量數(shù)據(jù)進(jìn)行反演.從表1可以看出,對于反演20個(gè)未知參數(shù)的情形,不加噪聲時(shí)反演結(jié)果非常精確,相對誤差l2范數(shù)為0103%;即使加入信噪比為80的噪聲后,反演結(jié)果也還很精確,相對誤差僅為1185%.這說明了DE反演算法對于比較小的m,高信噪比數(shù)據(jù)反演非常有效,精確.下面研究m=20(反演40個(gè)參數(shù))時(shí),正則化參數(shù)對反演結(jié)果的影響.反演結(jié)果如圖1所示:當(dāng)取10-8時(shí)計(jì)算結(jié)果相當(dāng)好;但不加正則化項(xiàng)(m取作0)時(shí),計(jì)算出來的f值振蕩得比較厲害;而正則化參數(shù)過大取10時(shí)又會(huì)把解磨得太光,甚至失去函數(shù)本該有的雙峰特性.所以,選擇合適的正則化參數(shù)對此類問題的求解起著至關(guān)重要的作用.目前還沒有非常好的

15、正則化參數(shù)選取方法,下一步-64)DErand2to2best1(xbest,G-xi,G)+F(xr1,G-xr2,G).vi,G+1=xi,G+實(shí)際反演計(jì)算表明,DErand2to2best1算法既保持了種群的多樣性,又具有較快的收斂速度,對此類NMR多指數(shù)反演問題最為有效.4131DE算法控制參數(shù)選擇控制參數(shù)對于一個(gè)全局進(jìn)化優(yōu)化算法的影響是比較大的,DE算法的控制參數(shù)選擇有如下經(jīng)驗(yàn)規(guī)則12,13:1)種群數(shù)量.根據(jù)問題規(guī)模,NP的合理選擇在5D10D之間.過大的NP導(dǎo)致算法收斂很慢,過小的NP可能導(dǎo)致DE算法過早收斂,得不到目標(biāo)函數(shù)的全局極小值.9期潘克家等:基于差分進(jìn)化算法的核磁共振T

16、2譜多指數(shù)反演5959希望能研究出一種自適應(yīng)的正則化參數(shù)選取算法,來更好地解決此類實(shí)際問題.表1m=10,SNR分別為+和80時(shí)反演結(jié)果T21T22T23T24T25T26T27T28T29T210Exact11000×10-311995×10-331981×10-371943×10-311585×10-231162×10-261310×10-211259×10-121512×10-151012×10-1SNR=11000×10-311995×10-331981×1

17、0-371943×10-311585×10-231163×10-261311×10-211259×10-121513×10-151017×10-1SNR=10091987×10-421013×10-341083×10-381232×10-311659×10-231356×10-261535×10-211279×10-121535×10-141998×10-1f1f2f3f4f5f6f7f8f9f10Exact51819

18、5;10-271660×10-111576115761164221370213712142521407214072142721038210392105231347313483131951028510294194621374213722130421383×10-121373×10-121362×10-1SNR=51821×10-271660×10-1SNR=10051888×10-271922×10-1感.在文獻(xiàn)10中給出了幾種傳統(tǒng)算法的信噪比穩(wěn)定區(qū)間,其中SVDS30,40<SNR,并且,發(fā)展了相應(yīng)的校正方法

19、提高反演質(zhì)量.文獻(xiàn)17中提到NNLS方法得到的T2譜可信度較高,為了與DE反演算法進(jìn)行比較,下面采用NNLS方法和改進(jìn)的SVD方法15反演m=20時(shí)的理論模型.并且,還假設(shè)T2弛豫時(shí)間布點(diǎn)已知,在程序中布點(diǎn)數(shù)為30,且其中的20個(gè)點(diǎn)正好與理論模型中給定的數(shù)據(jù)相圖1正則化參數(shù)對反演結(jié)果的影響同,然后利用NNLS方法和改進(jìn)SVD方法反演T2譜,計(jì)算結(jié)果分別如圖3和圖4所示.從圖中可以看出,NNLS算法和SVD算法在沒有誤差,或者信噪比較高的情況下計(jì)算非常準(zhǔn)確,但在有測量誤差比較大的情況下,反演得到的結(jié)果振蕩比較厲害.當(dāng)信噪比為60時(shí),NNLS方法已經(jīng)得不到滿意的反演結(jié)果,當(dāng)信噪比為30時(shí),改進(jìn)的S

20、VD算法得到的T2譜振蕩也比較厲害;而基于DE的優(yōu)化算法由于加入了曲率平滑項(xiàng),能夠很好地抑制T2譜的振蕩,即使在信噪比為10的情況下,也能得到如圖2所示的比較可信的T2譜.前面所有反演計(jì)算,都是采用DErand2to2best1求解優(yōu)化問題.DE控制參數(shù)選取如下:變異因子F取為018,交叉因子CR取為110.加權(quán)因子取為018,NP取為未知參數(shù)個(gè)數(shù)D的10倍.圖2中在各圖2信噪比對反演結(jié)果的影響圖2為在不同信噪比下選擇合適的正則化參數(shù)進(jìn)行反演,得到的一組反演結(jié)果.從圖中可以看出,DE算法對測量數(shù)據(jù)的信噪比不敏感,即使在噪聲比較大,信噪比僅為10時(shí),仍能反演出比較好的結(jié)果.而傳統(tǒng)的反演算法,如傳

21、統(tǒng)的SVD,NNLS等,由于沒有引進(jìn)正則化項(xiàng),對測量數(shù)據(jù)的誤差比較敏種信噪比情況下的計(jì)算時(shí)間和目標(biāo)函數(shù)值如表2所示.表中計(jì)算時(shí)間為用1024個(gè)數(shù)據(jù)反演40個(gè)參數(shù)值,Fortran程序在Pentium(R)D312G電腦上的運(yùn)行時(shí)間.從表2中可以看出,尤其在信噪比比較小的情況下,DE算法的收斂得非?？?當(dāng)SNR=105960物理學(xué)報(bào)57卷圖3信噪比對NNLS方法反演結(jié)果的影響圖5實(shí)際數(shù)據(jù)反演結(jié)果比較,對長慶某井巖心進(jìn)行T2譜測量.試驗(yàn)參數(shù)TE=017ms,回波采樣個(gè)數(shù)1024個(gè),經(jīng)數(shù)據(jù)壓縮后,實(shí)際回波采樣個(gè)數(shù)為128-10個(gè).DE反演算法正則化參數(shù)取為710×10,T2布點(diǎn)數(shù)m取為20

22、;而BRD和SVD算法T2均采取對數(shù)均勻布點(diǎn),布點(diǎn)數(shù)為301計(jì)算結(jié)果如圖5中所示,DE算法在不預(yù)先給定T2的情況下,得到的結(jié)果與已有算法取得的結(jié)果非常符合.61結(jié)論圖4信噪比對改進(jìn)SVD方法反演結(jié)果的影響時(shí),僅用了10s,就得到了如圖2中比較好的反演結(jié)果.由此也說明了DE反演算法非常適合低信噪比核磁共振弛豫信號(hào)反演.表2DE算法收斂速度SNR10102021406980104100295通過計(jì)算機(jī)模擬和實(shí)際算例,得到如下結(jié)論:11和傳統(tǒng)算法如SVD相比,基于DE的優(yōu)化算法不依賴于初值選取,計(jì)算穩(wěn)定;更容易實(shí)現(xiàn)非負(fù)約束,計(jì)算速度快.21優(yōu)化算法能同時(shí)反演出T2分布,得到分辨率比較高的T2譜.31

23、優(yōu)化算法引入正則化項(xiàng),抗噪能力比較強(qiáng),計(jì)算時(shí)間s目標(biāo)函數(shù)值11271150×10-11140×10-33170×10-54140×10-7可用于低信噪比T2譜反演.41優(yōu)化算法非常適合T2組分離散且分布較寬(m<40),信噪比較低的NRM弛豫信號(hào)反演;當(dāng)T25121應(yīng)用實(shí)例使用英國ResonanceInstruments公司的MARAN2組分較多甚至連續(xù)分布時(shí),可采用SVD,SIRT,NNLS等方法求解.1GoelmanG,PrammerMG1995J.Magn.Res.113112ZhangGQ,HirasakiGJ2003J.Magn.Res.

24、1635579期345678910潘克家等:基于差分進(jìn)化算法的核磁共振T2譜多指數(shù)反演MoffatBA,PopeJM2002Magn.Res.Imag.283XuF,HuangYR2002ActaPhys.Sin.51415(inChinese)許1213145961莫修文、王學(xué)秋2003吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)33232StornR,PriceK1997J.Glob.Optim.11341KaeloP,AliMM2006Eur.J.Oper.Res.1691176WangZD,XiaoLZ,LiuTY2003ScienceinChina(Ser.G)33323(inChinese)王忠東、

25、肖立志、劉堂宴2003中國科學(xué)G3332315JiangRZ,YaoYP,MiaoS,ZhangCS2005ActaPetro.Sin.2657(inChinese)姜瑞忠、姚彥平、苗盛、張春生2005石油學(xué)峰、黃永仁2002物理學(xué)報(bào)51415MengQA,HuCM1997ActaPhys.Sin.461961(inChinese)孟慶安、胡傳民1997物理學(xué)報(bào)461961ZhengSK,ChenZ,ChenZW,ZhongJH2001Chin.Phys.10558LiGY,XuXC1995ActaPhys.Sin.441847(inChinese)李鯁穎、徐學(xué)誠1995物理學(xué)報(bào)441847B

26、orgiaGC,BrownRJS,FantazziniP1998J.Magn.Res.13265BorgiaGC,BrownRJS,FantazziniP2000J.Magn.Res.147273LiaoGZ,XiaoLZ,XieRH,FuJJ2007ChineseJ.Geophys.50932(inChinese)廖廣志、肖立志、謝然紅、付娟娟2007地球181716報(bào)2657WangWM,LiP,YeZH2001ScienceinChina(Ser.A)31730(inChinese)王為民、李培、葉朝輝2001中國科學(xué)A31730WangH,LiGY2005ActaPhys.Sin.541431(inChinese)王鶴、李鯁穎2005物理學(xué)報(bào)541431ButlerJP,ReedsJA,DawsonSV1981SIAMJ.Numer.Anal.18381DunnKJ,JL,BergmanDJ1994SPE物理學(xué)報(bào)5093211WengAH,LiZB,MOXW,WangXQ2003J.JilinUniv.(EarthSci.Ed.)33232(inChinese)翁愛華、李舟波、Multi2inversionofT2spectrumi