核嶺回歸方法解釋致密砂巖儲(chǔ)層孔隙度

陳文浩, 王志章, 董少群, 侯加根

(1.中國(guó)石油大學(xué)(北京)地球科學(xué)學(xué)院, 北京 102249;2.中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249)

0 引 言

常規(guī)砂巖儲(chǔ)層通常依據(jù)聲波線性解釋儲(chǔ)層的孔隙度,其合理解釋是儲(chǔ)量估算及儲(chǔ)層評(píng)價(jià)工作的重點(diǎn)和難點(diǎn)。對(duì)于致密砂巖儲(chǔ)層,聲波與孔隙度不具線性關(guān)系,已有的孔隙度線性解釋方法不適于解釋致密砂巖儲(chǔ)層的孔隙度。王洪輝等[1]認(rèn)為簡(jiǎn)單的線性解釋致密砂巖儲(chǔ)層孔隙度難以滿足致密砂巖儲(chǔ)層評(píng)價(jià)的要求。張永浩等[2]模擬致密砂巖儲(chǔ)層條件時(shí),對(duì)聲波速度與孔隙度進(jìn)行實(shí)驗(yàn)分析,發(fā)現(xiàn)孔隙度隨壓力增大呈冪函數(shù)減小,而聲波速度隨壓力增大呈冪函數(shù)增大,所以孔隙度與聲波速度的線性關(guān)系隨壓力增大而變差。常規(guī)砂巖儲(chǔ)層孔隙度線性解釋模型不適于致密砂巖儲(chǔ)層的孔隙度解釋。與核函數(shù)融合形成的核嶺回歸算法[3-6]不僅可以有效處理非線性問(wèn)題,而且具有參數(shù)少、建模時(shí)間短、運(yùn)算速度快等優(yōu)點(diǎn)[7],具有廣泛的應(yīng)用空間,但在油氣儲(chǔ)層預(yù)測(cè)中的應(yīng)用研究較少。本文闡述核嶺回歸建立致密砂巖孔隙度解釋模型,并用來(lái)解釋致密砂巖儲(chǔ)層的孔隙度,分析多種方法的解釋結(jié)果,認(rèn)為核嶺回歸是建立致密砂巖儲(chǔ)層孔隙度解釋模型的有效方法。

1 核嶺回歸方法概述

線性回歸中的欠定問(wèn)題,最小二乘回歸的效果較差。為此,眾多學(xué)者提出了最小二乘回歸的改進(jìn)方法,其中比較經(jīng)典的是Hoerl和Kennard[8-9]提出的嶺回歸算法。之后Saunders等嘗試用核技巧改進(jìn)嶺回歸算法,核嶺回歸(Kernel Ridge Regression,KRR)由Cristianini和Shawe-Taylor正式使用[10],是簡(jiǎn)化的支持向量回歸機(jī),相對(duì)于支持向量回歸機(jī),具有引用參數(shù)少、節(jié)省機(jī)時(shí)等優(yōu)點(diǎn)。

1.1 嶺回歸算法

設(shè)有n個(gè)樣品,每個(gè)樣品有m個(gè)自變量和1個(gè)因變量,它們的觀測(cè)值構(gòu)成數(shù)據(jù)矩陣X和Y,即

X=x11x12…x1n

x21x22…x2n

???

xm1xm2…xmn

Y=y1y2…yn

(1)

式中,xij(i=1,2,…,m;j=1,2,…,n)為j樣品第i個(gè)自變量的觀測(cè)值;yj(j=1,2,…,n)為因變量的第j個(gè)觀測(cè)值。

最小二乘法線性回歸依據(jù)式(1)數(shù)據(jù),求解線性回歸模型y=b0+∑mi=1bixi中待定常數(shù)與待定系數(shù)b0、b1、…、bm的最佳估計(jì)值a0、a1、…、am,得線性回歸方程

j=A·Xj+a0(j=1,2,…,n)

(2)

式中,A=a1a2…am為線性回歸方程的系數(shù)向量;Xj為矩陣X的第j元素構(gòu)成的列向量;A·Xj為二者的內(nèi)積;a0為常數(shù)。求解線性回歸方程常數(shù)和系數(shù)的原則是使損失

Q=∑nj=1(yj-j)2

達(dá)到最小。

對(duì)于線性回歸中的欠定問(wèn)題,在損失最小原則下求解A和a0得到的方程會(huì)是病態(tài)方程。為解決該問(wèn)題,Hoerl和Kennard引入了損失函數(shù)

L(A,a0)=γA·A+Q

(3)

式中,A·A為正則化項(xiàng);Q為擬合損失項(xiàng);γ為正則化參數(shù),γ>0,它對(duì)A·A和Q起平衡作用,稱這種有權(quán)重因子的回歸算法為嶺回歸算法。若式(3)中γ=0,則為最小二乘意義下的回歸分析。

對(duì)于標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化變量,式(3)中a0等于0。因此。可將損失函數(shù)寫為矩陣形式

L(A)=γAAT+(Y-AX)(Y-AX)T

(4)

對(duì)式(4)中參數(shù)A求偏導(dǎo),并令其為0,可以得到嶺回歸模型的待定系數(shù)

A=(XXT+γI)-1XYT

(5)

式中,I為m×m的單位矩陣。

線性嶺回歸算法中正則化項(xiàng)的引入使回歸方程具有較好的穩(wěn)定性和泛化能力。但是,這種算法不能有效處理自變量間存在非線性相關(guān)的情況,這是線性回歸算法普遍存在的一個(gè)局限性。實(shí)際應(yīng)用中,用核技巧改進(jìn)后的嶺回歸算法,即核嶺回歸算法可以有效解決上述問(wèn)題。

1.2 核嶺回歸算法

1.2.1 核嶺回歸算法基本原理

將原始空間Rm中的n個(gè)樣品用某一非線性映射函數(shù)φ映射到一個(gè)高維特征空間F中,樣本數(shù)據(jù)變?yōu)棣?Xj),(j=1,2,…,n)。映射后可以在F中進(jìn)行嶺回歸,回歸方程的形式與式(2)類似,即

j=A·φ(Xj)+a0

(6)

相應(yīng)的損失函數(shù)也有類似的形式

L(A,a0)=γA·A+Q

(7)

對(duì)于標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化變量,損失函數(shù)的矩陣形式為

L(A)=γAAT+(Y-AX)(Y-AX)T

(8)

式中,X=φ(X1)φ(X2)…φ(Xn),權(quán)系數(shù)向量為A=a1a2…。對(duì)式(8)中A求偏導(dǎo),并令其為0,得

λA+XTAX=XTY

(9)

因非線性映射φ未知,故X、A未知。在實(shí)際應(yīng)用中一般求不出A的顯式表達(dá)式。將式(9)重新改寫后可得

A=λ-1XT(Y-AX)=XTα

(10)

其中,α=λ-1(Y-AX)。式(10)表明,A可以變換成樣品點(diǎn)的線性組合,α又可改寫為

α=(XXT+λI)-1Y=(K+λI)-1Y

(11)

其中,K稱為核矩陣,kij=φ(Xi)·φ(Xj)為核函數(shù),記為K(Xi,Xj)

對(duì)于新的樣品Xi,y的估計(jì)值

i=α·φ(Xj)=YT(K+λI)-1k

其中,k是n維的列向量,向量的元素為φ(Xi)與φ(Xj)的內(nèi)積。

通過(guò)上述變換,可以避免直接顯式地計(jì)算映射φ(xi),解決了非線性回歸的問(wèn)題。

1.2.2 核嶺回歸方法中的參數(shù)

(1) 核函數(shù)中參數(shù)。應(yīng)用核嶺回歸算法,首先要確定核函數(shù),總體思路是將輸入向量表示為內(nèi)積的形式。目前,應(yīng)用較廣泛的核函數(shù)有3種類型。

①多項(xiàng)式核函數(shù)K(x,y)d,d為常數(shù);

②神經(jīng)網(wǎng)絡(luò)核函數(shù)K(x,y)=tanh(β0(x,y)+β1)

③高斯徑向基核函數(shù)K(x,y)=e-‖x-y‖2σ2,σ為控制著2個(gè)樣品的相似程度的參數(shù)。

高斯徑向基核函數(shù)具有良好的泛化性能[11],本文將使用高斯徑向基核函數(shù),σ過(guò)小時(shí)易出現(xiàn)過(guò)渡擬合,過(guò)大時(shí)又會(huì)使擬合精度下降,因此,σ的取值涉及算法的擬合精度。

(2) KRR中正則化參數(shù)。對(duì)于核嶺回歸,除了核函數(shù),還需要確定正則化參數(shù)γ。從均方根誤差圖(見圖1)可知,正則化參數(shù)γ過(guò)大或過(guò)小都將減低擬合精度。

圖1 λ對(duì)擬合誤差的影響

2 實(shí)例研究

以紅崗油田H90區(qū)塊致密砂巖儲(chǔ)層為例進(jìn)行孔隙度解釋研究,該儲(chǔ)層主要分布于扶余油層(泉四段)。區(qū)內(nèi)發(fā)育辮狀河三角洲平原和前緣亞相,主要砂體為平原亞相的分流河道砂、溢岸薄層砂及前緣亞相的水下分流河道砂。主要巖性為長(zhǎng)石巖屑砂巖,多中砂質(zhì)細(xì)粒。泉四段油層物性較差,孔隙度為3%～13%,最高14.6%,平均8.03%;滲透率多為0.024～1 mD*非法定計(jì)量單位,1 mD=9.87×10-4 μm2; 1 ft=12 in=0.304 8 m; 下同,最高3.6 mD,平均0.294 mD,為低孔隙度超低滲透率儲(chǔ)層。

2.1 不同孔隙度解釋方法的解釋精度

2.1.1 核嶺回歸

(1) 樣本曲線。選擇與孔隙度相關(guān)密切的樣本曲線,以減少輸入變量,提高訓(xùn)練速度和精度。該地區(qū)三孔隙度測(cè)井響應(yīng)與孔隙度關(guān)系較為密切,其中聲波時(shí)差A(yù)C與孔隙度的相關(guān)系數(shù)為0.716 2,中子CNL與孔隙度的相關(guān)系數(shù)為0.664 1,密度DEN與孔隙度的相關(guān)系數(shù)為-0.730 1(見圖2)。

(2) 核嶺回歸解釋孔隙度。以三孔隙度測(cè)井曲線作為孔隙度解釋模型的輸入。利用核嶺回歸解釋孔隙度,還需要確定模型中的正則化參數(shù)γ和核函數(shù)參數(shù)σ。對(duì)于這2個(gè)參數(shù)的選取,還沒(méi)有統(tǒng)一公認(rèn)的最優(yōu)方法,目前常用網(wǎng)格搜索法和交叉驗(yàn)證法進(jìn)行優(yōu)選(見圖3)。

參數(shù)優(yōu)選的實(shí)現(xiàn)過(guò)程是給出γ和σ的取值區(qū)間,對(duì)于取定的γ和σ,把訓(xùn)練集作為原始數(shù)據(jù)集利用K-CV方法得到該組γ和σ下訓(xùn)練集預(yù)測(cè)結(jié)果,根據(jù)預(yù)測(cè)精度選擇γ和σ。網(wǎng)格搜索即在選定的范圍內(nèi),讓?duì)煤挺页手笖?shù)增長(zhǎng),首先選取粗網(wǎng)格進(jìn)行搜索,得到一組最優(yōu)參數(shù);然后逐步減小步長(zhǎng),最終獲得所需參數(shù)。圖3中紅點(diǎn)對(duì)應(yīng)解釋模型的最優(yōu)參數(shù),其中γ=0.064 8,σ=1.393 7。

圖2 三孔隙度參數(shù)與孔隙度交會(huì)圖

圖3 網(wǎng)格搜索法與交叉驗(yàn)證法確定γ和σ

2.1.2 解釋孔隙度精度比較

為了比較不同孔隙度解釋方法的解釋精度,分別采用線性回歸、逐步回歸、BP神經(jīng)網(wǎng)絡(luò)、支持向量回歸機(jī)方法解釋孔隙度,與核嶺回解釋孔隙度精度比較(見表1)。

表1 不同孔隙度解釋方法的解釋精度統(tǒng)計(jì)表

(1) 一元線性回歸解釋孔隙度模型分別為

φ=0.140AC-22.37

φ=0.664CNL+2.077

φ=-17.72DEN+53.94

(2) 多元線性回歸解釋孔隙度解釋模型為

φ=0.0753AC+0.0040CNL-9.8790DEN+12.8709

(3) 逐步回歸解釋孔隙度模型為

φ=0.076AC-10.1711DEN+17.8858

φ=-0.0013AC2-35.1002DEN2+

0.6739AC+0.1996CNL+169.809DEN

(4) BP神經(jīng)網(wǎng)絡(luò)解釋孔隙度

(5) 支持向量回歸機(jī)方法解釋孔隙度

根據(jù)表1各種方法的均方根誤差值可以看出,線性回歸解釋孔隙度方法雖然快捷,但是預(yù)測(cè)精度相對(duì)較差,4種非線性方法預(yù)測(cè)結(jié)果精度提升明顯,但相對(duì)而言,KRR方法的預(yù)測(cè)精度要略高于這些方法。表明針對(duì)致密砂巖儲(chǔ)層孔隙度的測(cè)井解釋,采用非線性算法提取的特征比線性算法提取的特征包含更準(zhǔn)確的信息,證明致密砂巖孔隙度解釋呈現(xiàn)非線性趨勢(shì)大于呈現(xiàn)線性趨勢(shì)。

為了對(duì)比3種非線性模型的執(zhí)行效率,對(duì)訓(xùn)練時(shí)間進(jìn)行了統(tǒng)計(jì)。表1中運(yùn)算時(shí)間為模型訓(xùn)練時(shí)間,統(tǒng)計(jì)結(jié)果在i5 2.0 GHz CPU、1 G RAM、MATLAB R2009a環(huán)境下得到。從表1中運(yùn)算時(shí)間可以看出支持向量回歸算法需要求解二次規(guī)劃問(wèn)題,運(yùn)算時(shí)間大大高于只需進(jìn)行線性方程組求解的核嶺回歸算法;同樣,神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練耗時(shí)1.9 s,核嶺回歸算法僅需0.5 s。核嶺回歸算法具有參數(shù)少、訓(xùn)練時(shí)間短、計(jì)算效率高的特點(diǎn)。

圖4 不同孔隙度解釋方法對(duì)H90井孔隙度解釋效果對(duì)比

由不同孔隙度解釋方法解釋H90井的孔隙度圖(見圖4)可以看出,單參數(shù)線性回歸和多參數(shù)線性回歸解釋的上部砂體的孔隙度均大于真實(shí)的孔隙度,3種非線性方法解釋的孔隙度與樣品孔隙度較接近。

3 結(jié) 論

(1) 核嶺回歸擬合方法能較好地解決非線性問(wèn)題。對(duì)于致密砂巖儲(chǔ)層,采用KRR解釋孔隙度具有更好的預(yù)測(cè)效果。

(2) 基于核函數(shù)的非線性方法方法建立嶺回歸解釋模型能夠盡可能地提取測(cè)井信息與儲(chǔ)層孔隙度之間的非線性映射關(guān)系,對(duì)其他地區(qū)致密砂巖儲(chǔ)層的孔隙度解釋具有一定的參考作用。

參考文獻(xiàn):

[1] 王洪輝, 黎鵬, 段新國(guó). 四川盆地須家河組低孔致密砂巖孔隙度測(cè)井解釋研究 [J]. 成都理工大學(xué)學(xué)報(bào): 自然科學(xué)版, 2009(3): 249-252.

[2] 張永浩, 杜環(huán)虹, 李新, 等. 蘇里格氣田致密砂巖儲(chǔ)層條件下聲波速度與孔隙度實(shí)驗(yàn)研究 [J]. 測(cè)井技術(shù), 2013, 37(3): 229-234.

[3] Vladimir N V. The Nature of Statistical Learning Theory [J]. New York: Springer-Verlag, 1995: 20-30.

[4] 印興耀, 孔國(guó)英, 張廣智. 基于核主成分分析的地震屬性優(yōu)化方法及應(yīng)用 [J]. 石油地球物理勘探, 2008(2): 179-183.

[5] 唐耀華, 張向君, 高靜懷. 基于地震屬性優(yōu)選與支持向量機(jī)的油氣預(yù)測(cè)方法 [J]. 石油地球物理勘探, 2009(1): 75-80.

[6] 鐘儀華, 李榕. 基于主成分分析的最小二乘支持向量機(jī)巖性識(shí)別方法 [J]. 測(cè)井技術(shù), 2009, 33(5): 425-429.

[7] 李琦, 邵誠(chéng). 基于核嶺回歸的非線性系統(tǒng)辨識(shí)及其應(yīng)用 [J]. 系統(tǒng)仿真學(xué)報(bào), 2009(8): 2152-2155.

[8] Arthur E Hoerl, Robert W Kennard. Ridge Regression: Biased Estimation for Nonorthogonal Problems [J]. Technometrics, 1970, 12(1): 55-67.

[9] Arthur E Hoerl, Robert W Kennard. Ridge Regression: Applications to Nonorthogonal Problems [J]. Technometrics, 1970, 12(1): 69-82.

[10] Nello Cristianini, John Shawe-taylor. An Introduction to Support Vector Machines and Other Kernel-based Learning Methods [M]. Cambridge: Cambridge University Press, 2000.

[11] Kristiaan Pelckmans, Jos De Brabanter, J A K Suykens. The Differogram: Non-parametric Noise Variance Estimation and its Use for Model Selection [J]. Neurocomputing (S0925-2312), 2005, 69(3): 100-122.