基于BG理論確定Phillips方法的正則化參數并評價反演結果

高 陽,肖立志,謝慶明

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

基于BG理論確定Phillips方法的正則化參數并評價反演結果

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

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

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

核磁共振測井;Phillips-Twomey方法;分辨率矩陣;協方差矩陣;折中準則

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

1 Phillips-Twomey方法及其折中解

1.1 Phillips-Twomey方法

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

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

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

1.2 正則參數對反演結果的影響

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

1.3 正則參數的選取

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

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

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

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

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

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

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

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

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

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

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

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

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

如圖7所示,點A即為ε取極小的點即最佳折中點,但曲線在A點附近極為平滑,因此雖然B、C兩點距離A點較遠,但是它們的ε值與A點相差甚微。考慮到r值選取應盡量的小,結合折中準則和r值最小準則,定義

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

圖7 ε-r關系圖Fig.7 ε-r curve

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

2 反演結果評價

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

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

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

3 NMR巖心資料處理

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

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

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

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

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

4 結 論

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

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

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

[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].武漢:中國地質大學出版社,1998.

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

[5] 肖立志,張恒榮,廖廣志,等.基于Backus-Gilbert理論的孔隙介質核磁共振弛豫反演[J].地球物理學報, 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方法對核磁共振弛豫測量結果的反演及其效果分析[J].西安石油大學學報:自然科學版,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] 肖立志,于慧俊,劉化冰,等.新型核磁共振孔隙介質分析儀的研制[J].中國石油大學學報:自然科學版,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譜與電阻率指數的對應性研究[J].中國石油大學學報:自然科學版,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] 王忠東,肖立志,劉堂宴.核磁共振弛豫信號多指數反演新方法及其應用[J].中國科學: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] 廖廣志,肖立志,謝然紅,等.孔隙介質核磁共振弛豫測量多指數反演影響因素研究[J].地球物理學報,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].中國石油大學學報:自然科學版, 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] 肖立志.核磁共振成像測井與巖石核磁共振及其應用[M].北京:科學出版社,1998.

[18] 鄧克俊.核磁共振測井理論及應用[M].東營:中國石油大學出版社,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

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

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