基于BG理論確定Phillips方法的正則化參數(shù)并評價反演結(jié)果

高 陽,肖立志,謝慶明

(1.中國石油大學(xué)理學(xué)院,北京102249;2.中國石油大學(xué)地球物理與信息工程學(xué)院北京102249; 3.重慶地質(zhì)礦產(chǎn)研究院國土資源部頁巖氣資源勘查重點實驗室,重慶400042)

基于BG理論確定Phillips方法的正則化參數(shù)并評價反演結(jié)果

高 陽1,肖立志2,謝慶明3

(1.中國石油大學(xué)理學(xué)院,北京102249;2.中國石油大學(xué)地球物理與信息工程學(xué)院北京102249; 3.重慶地質(zhì)礦產(chǎn)研究院國土資源部頁巖氣資源勘查重點實驗室,重慶400042)

正則參數(shù)的選取在Phillips-Twomey方法中起著關(guān)鍵性的作用,但一方面人們很難確定該參數(shù)的選取是否恰當(dāng),另一方面在參數(shù)選取時存在著很大的主觀因素。研究利用Phillips-Twomey方法反演核磁共振橫向弛豫時間譜時正則參數(shù)的確定,其確定方法是基于BG理論中的折中準(zhǔn)則和正則參數(shù)最小準(zhǔn)則;同時,還利用分辨率矩陣和協(xié)方差矩陣來評價反演所得橫向弛豫時間譜。數(shù)值模擬和巖心分析均表明,基于BG理論的折中準(zhǔn)則和正則參數(shù)最小準(zhǔn)則方法適用于確定Phillips-Twomey方法反演核磁共振橫向弛豫時間譜的正則參數(shù)r的選取。巖心分析結(jié)果表明,該方法反演效果比常規(guī)的SVD方法好,并且該方法避免了反演者的主觀判斷,因此反演結(jié)果沒有反演者個人偏好影響。

核磁共振測井;Phillips-Twomey方法;分辨率矩陣;協(xié)方差矩陣;折中準(zhǔn)則

Backus-Gilbert線性評價理論(即BG理論)[1-2]的核心是展布準(zhǔn)則和折中準(zhǔn)則[3-4]。目前已經(jīng)有學(xué)者將BG理論中的分辨率矩陣和協(xié)方差矩陣應(yīng)用到反演當(dāng)中[3]。Phillips-Twomey方法在地球物理反演和核磁共振反演中都有著非常重要的作用[5-9],而正則參數(shù)r的選取往往是根據(jù)反演者的經(jīng)驗,使得反演結(jié)果帶有反演者的主觀判斷?；诖?筆者引入BG理論中的折中準(zhǔn)則。對正則參數(shù)r的不同取值分別計算出分辨率矩陣和協(xié)方差矩陣,再算出展布和協(xié)方差矩陣的模,利用折中準(zhǔn)則和r值最小準(zhǔn)則確定正則參數(shù)r。同時嘗試?yán)肂G理論評價核磁共振測井反演結(jié)果,確定反演結(jié)果的可信度。

1 Phillips-Twomey方法及其折中解

1.1 Phillips-Twomey方法

核磁共振測井多指數(shù)反演結(jié)果的好壞直接影響到有關(guān)巖石儲集物性及流體特性參數(shù)的計算[10-15],該反演實際上是解超定線性方程組[16-17]Gm=d,其中G為M×N矩陣,m為N維列向量,d為M維列向量。假設(shè)d有誤差δd。由于系數(shù)矩陣G的條件數(shù)很高——在核磁共振測井中G的條件數(shù)可能達(dá)到1015數(shù)量級,造成d的微小誤差在m中都會被放大得很大,使得該問題嚴(yán)重不適定。Phillips-Twomey方法認(rèn)為解應(yīng)該是最光滑的,即解的二階導(dǎo)數(shù)的模長應(yīng)該最小。因此Phillips-Twomey方法定義光滑度量

并求s(m)在約束條件δd=Gm-d下的極小值點。由于誤差項δd未知,只能要求‖δd‖=‖Gm-d‖盡量地小。 因此,定義E(m)=rs(m)+‖Gmd‖2,并將E的極小值點m^作為線性方程組Gm=d的解。其中r為正則化參數(shù)或權(quán)重參數(shù),決定兩部分在E中所占比重。 令有

因此解m^=(GTG+rLTL)-1GTd。 此時,G的廣義逆G-g=(GTG+rLTL)-1GT。

1.2 正則參數(shù)對反演結(jié)果的影響

在Phillips-Twomey方法中,正則參數(shù)r的選取是十分關(guān)鍵的。 隨著r的增大,s(m)所占的比重增大,所得解就越光滑,但由于此時誤差項δd所占的比重相應(yīng)減小,因此原方程組在中所占的比重也減小;相應(yīng)地,若r減小,δd所占的比重增大,原方程組在中所占的比重也增大,但這時s(m)所占的比重減小,因此所得解可能不夠光滑。 如果r值選擇不當(dāng),不管是過大還是過小,都會導(dǎo)致反演結(jié)果失真。圖1說明了正則參數(shù)r的選取對反演結(jié)果的影響。從圖中可以看出,恰當(dāng)?shù)膔值為3,過大(如圖中r= 10)或者過小(如圖中r=0.5)都會導(dǎo)致反演結(jié)果失真。 考慮到反演的任務(wù)是解方程組,原方程組應(yīng)起主要作用。 因此在保證解足夠光滑的情況下,r值選取應(yīng)盡量小,并稱這一準(zhǔn)則為r值最小準(zhǔn)則。

1.3 正則參數(shù)的選取

按照BG理論,如果線性方程組Gm=d的廣義逆為G-g,則分辨率矩陣R=G-gG,協(xié)方差矩陣cov[m]=σ2G-g(G-g)T[1-4],這里假設(shè)d的誤差在統(tǒng)計上是獨立的且它們有相同的方差σ2。 核磁共振測井儀器的設(shè)計可以保證誤差獨立。核磁共振測井所采集的回波串,經(jīng)過相敏檢波器可以分離為兩道,一道是含噪聲的回波串?dāng)?shù)據(jù),一道是純噪聲,這樣就可以得到噪聲的方差[18]。 在Phillips-Twomey方法中,給定一個r,就可以求得一個廣義逆G-g= (GTG+rLTL)-1GT,也就可以得到分辨率矩陣和協(xié)方差矩陣。 圖2所示為r=0,數(shù)據(jù)方差σ2=1時的分辨率矩陣和協(xié)方差矩陣。

由圖2可以看出,分辨率矩陣比較接近單位矩陣,分辨率最佳。這是因為=G-gd=G-g(Gm)= (G-gG)m=Rm。當(dāng)R=I時,=m,即就是模型m。當(dāng)R≠I時,的某個分量,它不僅與mj有關(guān),而且還與其他的mi有關(guān),實際上,它是模型m的一個加權(quán)平均。 因此,可以用矩陣R與單位矩陣I的相似程度來衡量的好壞。R與I越接近,與m越接近,就越好;反之,就越壞。在實際應(yīng)用中,一般用矩陣R-I的L2模來衡量分辨率矩陣,并稱它為矩陣R的非對角元素的展布。

圖2 r=0,σ2=1時的分辨率矩陣(左)和協(xié)方差矩陣(右)Fig.2 Resolution matrix and covariance matrix when r=0 and σ2=1

圖3 r=0,σ2=1時的反演結(jié)果與模型對比Fig.3 Comparison of inversion result and model when r=0 and σ2=1

由于r=0,m^=(GTG)-1GTd,實際上就是最小方差解[2],因此是沒有正則化的解。 沒有正則化的解與真實解必然相距甚遠(yuǎn),圖3也正說明了這一點。

通過將r增大來觀察模型分辨率矩陣和協(xié)方差矩陣。 圖4是r=0.01,數(shù)據(jù)方差σ2=1時的模型分辨率矩陣和協(xié)方差矩陣。此時的模型分辨率矩陣與單位矩陣差別比較明顯,說明分辨率減小了。同時,協(xié)方差矩陣的數(shù)量級已經(jīng)從1012降低到了100,說明方差明顯改善。但是這是以犧牲分辨率為代價的。從分辨率矩陣和協(xié)方差矩陣的變化看,這種犧牲還是值得的,因為這是以較小的分辨率為代價換取較大方差的改善。圖5是r=0.01,數(shù)據(jù)方差σ2=1時,模型m與反演結(jié)果的對比。與圖3相比,反演結(jié)果與模型接近了許多,也說明了這種以犧牲分辨率為代價換取方差的改善是值得的。當(dāng)然,從圖中也可以看到,小弛豫分量反演結(jié)果與模型相差還是很大,因此這還不是理想結(jié)果。BG折中理論中的最佳折中解是以犧牲最小分辨率為代價換取最大方差減小的解。

實際上,要使分辨率高,就要使模型分辨率矩陣與單位矩陣差別小,即小;要使解的誤差小,就要使小。BG理論認(rèn)為分辨率和協(xié)方差是一對矛盾,二者不可兼得,只能取折中:要么犧牲一定的分辨率來換取協(xié)方差的減小,要么以較大的協(xié)方差來獲得較高的分辨率。這就是BG折中理論的核心思想。

圖6所示為Phillips-Twomey方法中隨r值的增加,spread(R)和的變化。

隨著r的增加,展布越來越大,即分辨率越來越低;但是方差卻越來越小。 這再一次說明要使spread(R)和同時達(dá)到極小是不可能的,因此BG折中理論定義

并對ε取極小。 此時的r值即為最佳折中解所對應(yīng)的r值。 為了確定哪一個r值使得ε取極小,可以做ε-r關(guān)系圖(圖7)。

圖4 r=0.01,σ2=1時的模型分辨率矩陣(左)和協(xié)方差矩陣(右)Fig.4 Resolution matrix and covariance matrix when r=0.01 and σ2=1

圖5 r=0.01,σ2=1時的反演結(jié)果與模型對比Fig.5 Comparison of inversion result and model when r=0.01 and σ2=1

圖6 spread(R)和‖cov[m]‖與r的變化關(guān)系Fig.6 Variaion of spread(R)and‖cov[m]‖with r

如圖7所示,點A即為ε取極小的點即最佳折中點,但曲線在A點附近極為平滑,因此雖然B、C兩點距離A點較遠(yuǎn),但是它們的ε值與A點相差甚微?？紤]到r值選取應(yīng)盡量的小,結(jié)合折中準(zhǔn)則和r值最小準(zhǔn)則,定義

并對其求最小,將此時的r值確定為Phillips-Twomey方法所需要的r值。圖8是利用最佳折中準(zhǔn)則和r值最小準(zhǔn)則確定Phillips-Twomey方法的正則參數(shù)r,反演不同信噪比(rSN)回波串得到的T2譜。

圖7 ε-r關(guān)系圖Fig.7 ε-r curve

圖8 不同信噪比回波串反演結(jié)果Fig.8 Inversion results of echoes with different rSN

2 反演結(jié)果評價

由圖8可以看出,對于不同的信噪比來說,反演結(jié)果與模型相當(dāng)符合,說明反演結(jié)果非常好。但在核磁共振測井反演時,并不知道真實的T2譜是怎樣的。此時對反演結(jié)果的評價可以通過分辨率矩陣和協(xié)方差矩陣來實現(xiàn)。

圖9是信噪比rSN=100,r=1.7時的分辨率矩陣和協(xié)方差矩陣。從整體上看,分辨率矩陣是對角占優(yōu)矩陣,雖然對角線上的值小于0.3,但對角線及其前后兩個的值加起來大部分都大于0.6(理論上R的行和應(yīng)該為1,由于在做非負(fù)約束時磨平了系數(shù)矩陣,使得R有些行和不為1,這是導(dǎo)致不是所有的值都大于0.6的主要因素)。以R的第18行為例,最大值在第18列為0.2836,這說明m18是18的主要貢獻(xiàn)者;而且第18行的第17列和第19列分別為0.2384和0.245 2,這三者加起來為0.769 2。因為T2譜是連續(xù)的,所以m17、m19與m18非常接近。從這個角度考慮,可以認(rèn)為m18的貢獻(xiàn)為0.7692。因此從分辨率矩陣的角度看,此次反演的可信度是比較高的。再從協(xié)方差矩陣看,絕對值最大值為0.0223,最小值幾乎為0,而信號噪聲的標(biāo)準(zhǔn)差為0.2947,這說明該方法明顯地壓制了信號的噪聲。因此,此次反演的誤差也是相當(dāng)小的。

圖9 rSN=100,r=1.7時的模型分辨率矩陣(左)和協(xié)方差矩陣(右)Fig.9 Resolution matrix and covariance matrix when r=1.7 and rSN=100

3 NMR巖心資料處理

圖10是根據(jù)折中準(zhǔn)則和r值最小準(zhǔn)則反演某地區(qū)3塊巖心所得T2譜(圖中標(biāo)有PH-T字樣的曲線)和用SVD方法得到的T2譜(圖中標(biāo)有SVD字樣的曲線)對比。這3塊巖心都是碳酸鹽巖,孔徑比較大,巖石孔隙度也相對較大,這一點從反演得到的T2譜也可以看出。根據(jù)資料,確定束縛水的T2截止值為92 ms,泥質(zhì)束縛水的T2截止值為4 ms[11-12]。

圖10 巖心反演T2譜Fig.10 Inversed T2spectrums of cores

從圖中可以看出,兩種反演方法在長弛豫時間上基本一致。在短弛豫方面,用SVD方法得到的第一塊巖心的T2譜只能反映出毛管束縛水,第二塊巖心的T2譜根本不能反映出束縛水信息;而根據(jù)折中準(zhǔn)則和r值最小準(zhǔn)則,得到的這兩塊巖心的T2譜都可以清晰地分辨出泥質(zhì)束縛水和毛管束縛水。由此可知,在短弛豫方面根據(jù)折中準(zhǔn)則和r值最小準(zhǔn)則得到的T2譜具有一定的優(yōu)勢。對于第三塊巖心,兩種方法反演的結(jié)果差別不大,說明本方法也不是每次都比SVD方法好。

圖11是第一塊巖心在上述r值下的分辨率矩陣和協(xié)方差矩陣。分辨率矩陣的峰值基本上都在對角線上,且峰值較高;同時協(xié)方差矩陣整體上來說值是比較小的,最大值沒有超過0.04。因此,這次反演的可信度還是比較高的。

圖11 第一塊巖心的分辨率矩陣(左)和協(xié)方差矩陣(右)Fig.11 Resolution matrix and covariance matrix of the first core

4 結(jié) 論

(1)數(shù)值模擬和巖心分析均表明基于BG理論的折中準(zhǔn)則和r值最小準(zhǔn)則方法適用于確定Phillips-Twomey方法反演核磁共振橫向弛豫時間譜的正則參數(shù)r的選取,且一般來說其效果比常規(guī)的SVD方法好。正則參數(shù)r的選取方法如下:首先對正則參數(shù)r的不同取值分別計算出分辨率矩陣和協(xié)方差矩陣,然后計算出分辨率矩陣和協(xié)方差矩陣的模,最后利用折中準(zhǔn)則和r值最小準(zhǔn)則確定正則參數(shù)r,即對公式(2)求極小。這樣在選取r時,避免了反演者的主觀判斷,因此反演結(jié)果沒有反演者個人偏好的影響。

(2)數(shù)值模擬所設(shè)定的地層巖性是砂巖,巖心分析中的巖心為碳酸鹽巖,所以本文所述方法適用于砂巖和碳酸鹽巖地層,對于其他常規(guī)儲層,理論上是能夠正確反演的。對于頁巖等非常規(guī)儲層,尚不清楚是否適用。

(3)利用BG評價理論評價核磁共振測井反演結(jié)果,確定反演結(jié)果可信度差別。其方法是看分辨率矩陣主峰是否都在對角線上,其他值是否接近0,協(xié)方差矩陣是否接近0矩陣。如果是,則反演結(jié)果比較好,反之則較差。如果想評價某一個分量,就看它所在的行是否滿足上述條件。

[1] BACKUS G E,GILBERT J F.The resolving power of gross earth data[J].Geophysical Journal,1968,16(2): 169-205.

[2] BACKUS G E,GILBERT J F.Uniqueness in the inversion of inaccurate gross earth data[J].Philosophical Transactions,1970,266(1173):123-192.

[3] 王家映.地球物理反演理論[M].武漢:中國地質(zhì)大學(xué)出版社,1998.

[4] 楊文采.地球物理反演的理論和方法[M].北京:地質(zhì)出版社,1996.

[5] 肖立志,張恒榮,廖廣志,等.基于Backus-Gilbert理論的孔隙介質(zhì)核磁共振弛豫反演[J].地球物理學(xué)報, 2012,55(11):3821-3828. XIAO Li-zhi,ZHANG Heng-rong,LIAO Guang-zhi,et al.Inversion of NMR relaxation in porous media based on Backus-Gilbert theory[J].Chinese J Geophys,2012,55 (11):3821-3828.

[6] PHILLIPS D L.A technique for the numerical solution of certain integral equations of the first kind[J].J ACM, 1962,9:84-87.

[7] TWOMEY S.On the numerical solution of Fredholm integral equations of the first kind by the inversion of linear systems profuctuced by quadrature[J].J ACM,1963, 10:97-101.

[8] 高陽,肖立志,謝慶明.利用Phillips-Twomey方法對核磁共振弛豫測量結(jié)果的反演及其效果分析[J].西安石油大學(xué)學(xué)報:自然科學(xué)版,2012,27(5):32-38. GAO Yang,XIAO Li-zhi,XIE Qing-ming.Inversion of NMR transverse relaxation time using Phillips-Twomey method[J].Journal of Xi?an Shiyou University(Natural Science Edition),2012,27(5):32-38.

[9] 肖立志,于慧俊,劉化冰,等.新型核磁共振孔隙介質(zhì)分析儀的研制[J].中國石油大學(xué)學(xué)報:自然科學(xué)版,2013,37(3):68-73.XIAO Li-zhi,YU Hui-jun,LIU Hua-bing,et al.A novel low field nuclear magnetic resonance analyzer for porous media[J].Journal of China University of Petroleum(E-dition of Natural Science),2013,37(3):68-73.

[10] MILLER A,CHEN S A.New method for estimating T2distributions from NMR measurements[J].Mag Re Imaging,1998,16:617.

[11] BUTLER J P,REEDS J A,DAWSON S V.Estimating solutions of first kind integral equations with nonnegative constraints and optimal smoothing[J].SIAMJ Numer Anal,1981,18(3):381-397.

[12] 葛新民,范宜仁,吳飛,等.巖石核磁共振T2譜與電阻率指數(shù)的對應(yīng)性研究[J].中國石油大學(xué)學(xué)報:自然科學(xué)版,2012,36(6):53-56. GE Xin-min,FAN Yi-ren,WU Fei,et al.Correspondence of core nuclear magnetic resonance T2spectrum and resistivity index[J].Journal of China University of Petroleum(Edition of Natural Science),2012,36(6): 53-56.

[13] 王忠東,肖立志,劉堂宴.核磁共振弛豫信號多指數(shù)反演新方法及其應(yīng)用[J].中國科學(xué):G輯,2003,33 (4):323-332. WANG Zhong-dong,XIAO Li-zhi,LIU Tang-yan.New inversion method with multi-exponent and its application [J].Science in China(G),2003,33(4):323-332.

[14] 廖廣志,肖立志,謝然紅,等.孔隙介質(zhì)核磁共振弛豫測量多指數(shù)反演影響因素研究[J].地球物理學(xué)報,2007,50(3):932-938. LIAO Guang-zhi,XIAO Li-zhi,XIE Ran-hong,et al. Influence factors of multi-exponential inversion of NMR relaxation measurement in porous media[J].Chinese J Geophys,2007,50(3):932-938.

[15] 謝然紅,肖立志,張建民,等.低滲透儲層特征及測井評價方法[J].中國石油大學(xué)學(xué)報:自然科學(xué)版, 2006,30(1):47-51. XIE Ran-hong,XIAO Li-zhi,ZHANG Jian-min,et al. Low permeability reservoir characteristics and log evaluation method[J].Journal of China University of Petroleum(Edition of Natural Science),2006,30(1):47-51.

[16] COATES G R,XIAO L Z,PRAMMER M G.NMR logging principles and applications[M].Texas:Gulf Publishing Company,1999.

[17] 肖立志.核磁共振成像測井與巖石核磁共振及其應(yīng)用[M].北京:科學(xué)出版社,1998.

[18] 鄧克俊.核磁共振測井理論及應(yīng)用[M].東營:中國石油大學(xué)出版社,2010.

(編輯 修榮榮)

Evaluating regular parameter of Phillips method and assessing inversions based on BG theory

GAO Yang1,XIAO Li-Zhi2,XIE Qing-Ming3
(1.College of Science in China University of Petroleum,Beijing 102249,China; 2.College of Geophysics and Information Engineering in China University of Petroleum,Beijing 102249,China; 3.Key Laboratory of Shale Gas Exploration,Ministry of Land and Resources,Chongqing Institute of Geology and Mineral Resources,Chongqing 400042,China)

In the Phillips-Twomey method,it is essential to find a suitable regular parameter r.The selection of the parameter r remains practically difficult,and is often biased by researcher?s subjectivity.Finding a suitable regular parameter r in using Phillips-Twomey method to inverse nuclear magnetic resonance(NMR)logging transverse relaxation time was focused on.The way to finding the suitable r was based on the compromise criterion in BG theory and the r minimal criterion.The inversed transverse relaxation time was also assessed through the resolution matrix and the covariance matrix.Digital simulation and core analysis show that the approach in determining the regular parameter r of the Phillips-Twomey inverse method works very well for the inverse of NMR transverse relaxation time.Core analysis also shows that this method is better than the classical SVD method without subjective bias.

nuclear magnetic resonance logging(NMR logging);Phillips-Twomey method;resolution matrix;covariance matrix;compromise criterion

O 175.5

:A

1673-5005(2014)03-0050-07

10.3969/j.issn.1673-5005.2014.03.008

2013-10-11

中國石油大學(xué)(北京)科研基金資助(KYJJ2012-06-02);中國石油科技創(chuàng)新基金項目(2011D-5006-0307)

高陽(1979-),男,講師,博士,主要從事核磁共振測井反演方法研究。E-mail:gaoyang1203@163.com。