巖石磁共振多指數(shù)反演算法研究

巖石磁共振多指數(shù)反演算法研究

磁共振石油巖心分析和鉆探技術(shù)可以提供豐富的巖石和流量特征參數(shù)，例如巖石中的空隙大小分布；有效孔隙度、滲透率、可動流和不可動流之間的關(guān)系等。通過差譜法和移譜法，可以定性或定量地提供地層含油氣的指示。通過測量擴散系數(shù)，我們可以獲得各種重要信息，如巖層的潤濕度和流粘度。這些獨特的巖石和材料的物理信息為地質(zhì)勘探和油田開發(fā)工程師以及制定儲層評價、產(chǎn)量預(yù)測和開發(fā)方案提供了重要依據(jù)。磁共振技術(shù)在油田勘探和開發(fā)中的應(yīng)用引起了國內(nèi)外的高度重視。核磁共振石油應(yīng)用技術(shù)是一個正在發(fā)展的新技術(shù),本身還不完善,近年來,我國花費了巨額外匯引進(jìn)多臺美國NUMAR公司核磁共振成像測井儀,但應(yīng)用表明:國外核磁測井解釋軟件由于反演指數(shù)少、弛豫譜精度差、無法反映裂縫及溶洞等復(fù)雜的地質(zhì)情況,應(yīng)用效果不好.要解決這個問題,必須研究更精確的多弛豫反演算法,并形成適合于我國復(fù)雜陸相儲層地質(zhì)條件的巖石核磁測量數(shù)據(jù)的解釋標(biāo)準(zhǔn)、解釋模型.在自然科學(xué)中,很多參數(shù)都是采用間接方法獲得的,也就是說,與觀測數(shù)據(jù)相關(guān)的預(yù)測函數(shù)或預(yù)測矢量是由某種算符確定的,這種算符一般可采用線性積分算符來表示.而在觀測數(shù)據(jù)中通常帶有未知的噪聲,問題關(guān)鍵是怎樣從含有噪聲的數(shù)據(jù)中估計或預(yù)測出我們所期望獲得的函數(shù)或矢量.現(xiàn)代核磁共振石油工業(yè)應(yīng)用測試技術(shù)也是一種間接測量技術(shù),它所采集到的原始數(shù)據(jù)是地下巖石孔隙中所含流體的核磁共振信號,即所謂的自旋回波信號.我們所需的巖石物性參數(shù)以及流體屬性等信息是從原始回波信號中提取出來的,提取這些信息的關(guān)鍵在于對巖石自旋回波信號進(jìn)行多弛豫反演,這個反演實際上是求解第一類Fredholm線性積分方程,這種積分方程通常的解不是惟一的.在對算符求逆時,積分核作用到函數(shù)上,起著光滑的作用,輸入原始數(shù)據(jù)中的小變化(如核磁共振測量過程中的噪聲等)對解的結(jié)果有較大的影響.光滑的結(jié)果常常是丟失了信息,而逆運算又無法找回它,因此,我們必須利用解的一些先驗知識(如解的連續(xù)性、數(shù)據(jù)的信噪比、解的非負(fù)限制等)來最大限度地恢復(fù)我們所需的信息.本文從第一類積分方程的反演求解入手,推導(dǎo)出適合于巖石核磁弛豫多指數(shù)反演的兩種算法:奇異值分解(singularvaluedecomposition,SVD)方法和變換反演算法,分析了兩種算法優(yōu)缺點,給出了最優(yōu)反演模型的選取原則.1孔隙磁化強度信號有兩個參數(shù)可用來描述核(在巖石核磁共振中指氫核)自旋信號的弛豫特征:縱向弛豫時間T1和橫向弛豫時間T2.由于T1的測量時間很長,所以在巖石核磁共振中一般進(jìn)行T2測量.根據(jù)核磁共振理論分析,單個孔隙中的磁化強度信號的衰減滿足單指數(shù)衰減規(guī)律,但由于巖石內(nèi)部是由一系列大小不等的孔隙群體組成的,所以在巖石核磁共振中測得的總磁化強度信號y(t)是一系列單個孔隙磁化強度信號的疊加,同時在實際測量過程中,不可避免地產(chǎn)生隨機噪聲,因此y(t)可描述為y(t)=∑ifi?exp(-tΤ2i)+ε(t),t=n?τ,(1)y(t)=∑ifi?exp(?tT2i)+ε(t),t=n?τ,(1)其中fi為第i類孔隙在總孔隙中所占的份額,T2i為第i類孔隙的T2弛豫時間,τ為回波間隔時間,ε(t)為隨機噪聲序列.巖石核磁共振分析的關(guān)鍵是如何從(1)式中求解出各類孔隙的T2弛豫時間T2i以及各類孔隙在總孔隙中所占的份額fi,我們有時也稱之為T2譜.2不同f2弛豫時間分布函數(shù)的估計結(jié)果求解方程(1),實際上是解第一類Fredholm積分方程,這是一個非適定問題,也就是說,在允許的測量誤差ε條件下,存在不同的f(T2)弛豫時間分布函數(shù)都能相當(dāng)好地擬合原始回波衰減曲線,這種相容性增大了反演過程的難度.2.1svd分解的明確構(gòu)造SVD算法可用來求解大多數(shù)的線性最小二乘法問題.SVD算法基于如下分解定理:對任意的矩陣Am×n,當(dāng)其行數(shù)m大于等于列數(shù)n時,可以分解為正交矩陣Um×n,非負(fù)對角矩陣Wn×n以及正交矩陣Vn×n的轉(zhuǎn)置的乘積,即Am×n=Um×n?[diag(wj)]n×n?VΤn×n,(2)Am×n=Um×n?[diag(wj)]n×n?VTn×n,(2)其中wj≥0(1≤j≤n);U,V為正交矩陣,即滿足m∑i=1uijuik=δjk?n∑i=1vijvik=δjk(1≤j,k≤n),(3)∑i=1muijuik=δjk?∑i=1nvijvik=δjk(1≤j,k≤n),(3)當(dāng)m<n時,SVD也可以執(zhí)行,在這種情況下,奇異值wj=0(m+1≤j≤n),并且U中相應(yīng)的列都是零,這時(3)式僅對j,k≤m時成立.故不管矩陣A是否是奇異,(2)式的分解總可以進(jìn)行,而且這個分解幾乎是惟一的.也就是說,其分解形式惟一到:對矩陣U的列、W的元素和V的列能做相同的置換,或者矩陣U和V的任意列的線性組合,在W中對應(yīng)的元素仍恰好完全相同.SVD分解明確地構(gòu)造了矩陣零空間和值域的正交標(biāo)準(zhǔn)化基.特別地,對U的列,若與其標(biāo)號相同的元素wj為非零元,則其列為值域的一個正交標(biāo)準(zhǔn)化的基礎(chǔ)矢量;對V的列,若與其標(biāo)號相同的wj為零,則其列為零空間的一個正交標(biāo)準(zhǔn)化基.對于如下的多指數(shù)衰減T2模型,有y=Μ?f,(4)y=M?f,(4)其中y=(y1,y2,…,yn)T為測量的自旋回波衰減信號,M=[mij]n×m=[exp(-ti/T2j)]n×m;f=(f1,f2,…,fm)T為弛豫時間T2j對應(yīng)的各點的幅度值,T2j(j=1,2,…,m)為預(yù)先指定的T2時間分布系列,典型的取法為在(T2min,T2max)區(qū)間內(nèi)對數(shù)均勻地選取m個點,我們稱為弛豫時間布點,也可采用2的冪指數(shù)布點、線性均勻布點等方式.矩陣條件數(shù)的定義為矩陣的最大特征值與最小特征值的比值.若矩陣的條件數(shù)為無窮大,則該矩陣奇異;若矩陣的條件數(shù)太大,即其倒數(shù)超出了機器的浮點精度,則稱該矩陣為病態(tài)的矩陣.若直接采用Gauss分解求上式,幾乎是不可能的,原因是矩陣M的條件數(shù)相當(dāng)大,例如:若回波間隔τ=1.2×10-3s,T2在0.1×10-3~10000×10-3s內(nèi)對數(shù)均勻地取50個點,則矩陣M的條件數(shù)可達(dá)1016數(shù)量級,很明顯,矩陣M是高度病態(tài)的.采用SVD分解法來求解上式,系數(shù)矩陣Mn×m=Un×m·[diag(wj)]m×m·VTm×m,這里U,V為正交矩陣,diag(wj)為對角矩陣,其對角元遞減排列,則我們就可以很容易地求得最小二乘意義下的解為?f=V?[diag(1w1,1w2,??SΝRw1,0,??0)]?(UΤ?y).(5)f?=V?[diag(1w1,1w2,??SNRw1,0,??0)]?(UT?y).(5)這里給出了矩陣條件數(shù)小于等于SNR的限制,避免了解的不穩(wěn)定性.其中SNR為從測量數(shù)據(jù)中估計出的信噪比.SNR定義為第1個回波的幅度值除以誤差矢量r(r=y-Μ??f)的標(biāo)準(zhǔn)差σ.2.2im1,,首先給定如下的目標(biāo)函數(shù):χ2=n∑i=1[yi-m∑j=1(fj?mij)]2+λ?m∑j=1f2j=∥y-Μf∥2+λ∥f∥2,(6)這里M=[mij]=[exp(-ti/T2j)],λ為平滑因子.對幅度f=(f1,f2,…,fm)T的第k分量求極值并令其等于0,則有?χ2?fk=-2n∑i=1(yi-m∑j=1fj?mij)?mik+2λ?fk=0.(7)交換求和順序,并移項整理,可得m∑j=1[fj?n∑i=1(mik?mij)]+λfk=n∑i=1(mik?yi).(8)很容易驗證,k=1,2,…,m的m個等式組成的方程組滿足(采用矩陣形式給出)(ΜΤΜ)?f+λΙm×m?f=ΜΤ?y,(9)上式中Im×m為m×m單位矩陣.我們對方程(9)作如下線性變換,令f=ΜΤ?c,(10)注意(10)式中未知變量c=(c1,c2,…,cn)T為n×1維的,而不是m×1維的.采用上述變換我們將m維T2域空間的解變換到n維時域空間來求解.將(10)式代入(9)式,則有ΜΤ?(Μ?ΜΤ+λΙn×n)?c=ΜΤ?y,(11)則原問題的解f可以通過求解出方程(Μ?ΜΤ+λΙn×n)?c=y(12)的解c,再通過線性變換(10)式回代而獲得.選擇合適的λ,以保證矩陣(M·MT+λIn×n)的可逆,則我們就可以很容易地求得方程(1)的最小二乘解?f=ΜΤ?(Μ?ΜΤ+λΙn×n)-1?y.(13)平滑因子λ的選取:增大λ會使解的殘差r=∥y-Μ?f∥遠(yuǎn)離最小的情形,趨向于使‖?f‖最小;反之,減小λ又會遇到矩陣(Μ?ΜΤ+λΙn×n)病態(tài)的問題,使我們獲得的解變得不可靠.λ的選取應(yīng)從主觀和客觀兩個方面加以判斷,客觀上要求λ盡可能地小,主觀上要求解盡可能地光滑連續(xù).根據(jù)先驗信息,若我們已知測量誤差σ,則比較理想的平滑因子為λ=√n?σ∥c∥,(14)還可以選擇不同平滑項,以實現(xiàn)不同的平滑方式.2.3svd算法的改進(jìn)若對上述兩種算法不加任何限制條件,則求解出的T2譜的幅度?f中很可能存在小于零的分量,這在物理意義上來說是不合邏輯的.故在求解模型中,還必須加上非負(fù)性限制條件,即fj≥0(j=1,2,…,m).兩種算法實現(xiàn)非負(fù)性限制的方法是不一樣的.SVD算法采用迭代法逐漸消去負(fù)的組分,即在我們求出的?f中,找出具有最大負(fù)數(shù)的組分,假定為第k項,說明第k項最不合理,應(yīng)該在原始模型中去除該組分,令與fk有關(guān)的所有系數(shù)為零,消去第k項,然后重新進(jìn)行反演計算,直到滿足所有的fj≥0(j=1,2,…,m)為止.這種方法的缺點是會破壞T2分布的連續(xù)性,造成T2譜的畸形.變換反演算法也是采用迭代法逐漸消去負(fù)數(shù)的組分,但該方法通過T2域到時域的變換后,將第k項相鄰點的T2分布信息疊加到負(fù)數(shù)的T2k組分中,重新進(jìn)行反演計算,直到滿足所有的fj≥0(j=1,2,…,m)為止.這種方法的優(yōu)點是可以保持T2分布的連續(xù)性,使反演結(jié)果更符合實際.3t2反演實驗變換反演算法采用變換的方式,將T2域空間的解變換到時域空間進(jìn)行求解,優(yōu)點是算法在實現(xiàn)非負(fù)性限制條件時,可利用相鄰點的T2分布信息來保持T2分布的連續(xù)性,并且其設(shè)計矩陣(MMT+λI)的大小固定不變,總是n×n的,算法穩(wěn)定且容易實現(xiàn).SVD算法是對設(shè)計矩陣直接進(jìn)行SVD分解,采用原始數(shù)據(jù)的信噪比SNR作為矩陣條件數(shù)的限制,且SNR可以在擬合過程中動態(tài)的獲取,避免了變換反演算法中提到的平滑參數(shù)λ的選取問題.但SVD算法在非負(fù)限制條件的實現(xiàn)過程中,有可能造成T2譜的畸形,特別是當(dāng)原始數(shù)據(jù)信噪比較低時,反演結(jié)果的分辨率比較低,解會出現(xiàn)不規(guī)則的跳動,其原因在于低信噪比的數(shù)據(jù)會使解的自由度增大.我們可用如下的實例來說明這個問題.假設(shè)核磁測量的回波間隔τ=1.2×10-3s,T2弛豫時間在0.1×10-3~10000×10-3s內(nèi)對數(shù)均勻地選取50個點,采用SVD算法進(jìn)行反演時,在不同信噪比下所獲得的解的自由度見表1.解的自由度越大,解中所包含的線性無關(guān)的信息含量就越少,也就是說,解的分辨率越低,我們所獲得的信息量就越小.同樣,解的自由度太小,也會造成嚴(yán)重的誤差,原因是解對數(shù)據(jù)的適應(yīng)性差,表明我們所選的模型不合適,如T2布點數(shù)太少等.目前的核磁測井?dāng)?shù)據(jù)的信噪比一般為30左右,采用SVD算法反演出的T2譜中所含有的線性無關(guān)的信息量不超過4,T2譜的分辨率比較低.這就是美國Numar核磁解釋處理軟件所反演出的T2譜大多呈現(xiàn)一個峰的緣故.同時,SVD算法計算處理相當(dāng)費時,尤其是T2布點數(shù)較多時就更為明顯.我們的計算結(jié)果表明當(dāng)原始數(shù)據(jù)的信噪比大于80時,SVD算法可以給出好的結(jié)果.為了便于對比分析,我們構(gòu)造了一個具有雙峰結(jié)構(gòu)的典型石油巖心T2譜(圖1),假定回波間隔時間τ=1.2×10-3s,回波數(shù)NE=1024,可構(gòu)造出無噪聲的回波衰減信號.在4種不同T2對數(shù)均勻布點數(shù)(100,40,15和8)條件下,用SVD法進(jìn)行反演計算,圖1給出了反演結(jié)果的對比.從圖中可以看出,T2布點數(shù)設(shè)為100和40時所得的結(jié)果差異較小,而T2布點數(shù)為15和8的譜分辨率明顯變差.因此為保證T2分布的真實性,在反演過程中應(yīng)保持一定的T2點數(shù)和覆蓋范圍.我們的研究結(jié)果表明,T2布點數(shù)為30~50較為合適,布點區(qū)間應(yīng)跨越3個數(shù)量級.為了對兩種解法進(jìn)行直觀對比分析,在上述構(gòu)造的無噪聲回波衰減信號中,逐點加上隨機Gauss白噪聲,構(gòu)造出不同信噪比的回波串.分別對SNR=∞,100,50,30,20,10的6個回波串(其中SNR=∞為無噪聲的回波串)進(jìn)行反演,反演中均采用對數(shù)均勻布點,布點區(qū)間為0.1×10-3~10000×10-3s,布點數(shù)為50.通過反演共獲得12個不同的T2譜.圖2給出了SVD反演的6個T2譜的對比.信噪比≥50的3條T2譜線具有雙峰結(jié)構(gòu),與真實T2譜線非常接近;SNR為30,20的2條T2譜線只具有單峰結(jié)構(gòu),其單峰仍能較好地反映出真實譜的大峰結(jié)構(gòu);而SNR為10的T2譜線與真實譜線相差較遠(yuǎn),且存在不連續(xù)的現(xiàn)象.從表1中可以發(fā)現(xiàn)上述幾種情況的解所包含的線性無關(guān)信息量是不同的,這就是它們的譜與真實譜差異的原因.圖3給出了變換反演算法反演的結(jié)果對比.SNR≥20的T2譜線均呈現(xiàn)雙峰態(tài)結(jié)構(gòu),與真實譜吻合較好,而SNR為10的T2譜線偏離真實譜線較遠(yuǎn).4精細(xì)解釋及其應(yīng)用作為一個應(yīng)用實例,這里介紹本文變換反演算法在青海柴達(dá)木盆地獅子溝構(gòu)造獅32井中的應(yīng)用情況,獅子溝儲層巖石孔隙結(jié)構(gòu)比較復(fù)雜,既存在基質(zhì)孔,也存在裂縫和微裂縫,還存在少量溶蝕孔,其T2譜大部分呈三峰態(tài).青海油田采用美國NUMAR公司生產(chǎn)的C型核磁共振成像測井儀對該井進(jìn)行了測井,并采用該公司的解釋軟件對核磁回波串進(jìn)行了處理,NUMAR公司所采用的是SVD算法,T2布點方式為2的冪次方,布點數(shù)只有8個,因而很難反映復(fù)雜孔隙結(jié)構(gòu)的細(xì)節(jié),解釋結(jié)果無法用于石油儲層評價和儲量計算.針對上述問題,我們以本文變換反演算法為基礎(chǔ),編寫了一套核磁測井?dāng)?shù)據(jù)精細(xì)解釋軟件,對該井的核磁測