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

(1.中國石油大學(xué)(北京)地球科學(xué)學(xué)院, 北京 102249;2.中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室, 北京 102249)

0 引 言

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

1 核嶺回歸方法概述

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

1.1 嶺回歸算法

設(shè)有n個樣品,每個樣品有m個自變量和1個因變量,它們的觀測值構(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個自變量的觀測值;yj(j=1,2,…,n)為因變量的第j個觀測值。

最小二乘法線性回歸依據(jù)式(1)數(shù)據(jù),求解線性回歸模型y=b0+∑mi=1bixi中待定常數(shù)與待定系數(shù)b0、b1、…、bm的最佳估計值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

達到最小。

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

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

(3)

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

對于標準差標準化變量,式(3)中a0等于0。因此?？蓪p失函數(shù)寫為矩陣形式

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

(4)

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

A=(XXT+γI)-1XYT

(5)

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

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

1.2 核嶺回歸算法

1.2.1 核嶺回歸算法基本原理

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

j=A·φ(Xj)+a0

(6)

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

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

(7)

對于標準差標準化變量,損失函數(shù)的矩陣形式為

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

(8)

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

λA+XTAX=XTY

(9)

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

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

(10)

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

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

(11)

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

對于新的樣品Xi,y的估計值

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

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

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

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

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

①多項式核函數(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個樣品的相似程度的參數(shù)。

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

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

圖1 λ對擬合誤差的影響

2 實例研究

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

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

2.1.1 核嶺回歸

(1) 樣本曲線。選擇與孔隙度相關(guān)密切的樣本曲線,以減少輸入變量,提高訓(xùn)練速度和精度。該地區(qū)三孔隙度測井響應(yīng)與孔隙度關(guān)系較為密切,其中聲波時差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) 核嶺回歸解釋孔隙度。以三孔隙度測井曲線作為孔隙度解釋模型的輸入。利用核嶺回歸解釋孔隙度,還需要確定模型中的正則化參數(shù)γ和核函數(shù)參數(shù)σ。對于這2個參數(shù)的選取,還沒有統(tǒng)一公認的最優(yōu)方法,目前常用網(wǎng)格搜索法和交叉驗證法進行優(yōu)選(見圖3)。

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

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

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

2.1.2 解釋孔隙度精度比較

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

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

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

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

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

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

3 結(jié) 論

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

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

參考文獻:

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

[2] 張永浩, 杜環(huán)虹, 李新, 等. 蘇里格氣田致密砂巖儲層條件下聲波速度與孔隙度實驗研究 [J]. 測井技術(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] 印興耀, 孔國英, 張廣智. 基于核主成分分析的地震屬性優(yōu)化方法及應(yīng)用 [J]. 石油地球物理勘探, 2008(2): 179-183.

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

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

[7] 李琦, 邵誠. 基于核嶺回歸的非線性系統(tǒng)辨識及其應(yīng)用 [J]. 系統(tǒng)仿真學(xué)報, 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.