核主成分分析法在測井濁積巖巖性識別中的應用

周 游 張廣智 高 剛 趙 威 易院平 魏紅梅

(①中國石油大學(華東)地球科學與技術學院，山東青島266580； ②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266071； ③油氣資源與勘探技術教育部重點實驗室，湖北武漢 430100；④長江大學地球物理與石油資源學院，湖北武漢 430100； ⑤中冶集團武漢勘察研究院有限公司，湖北武漢 430080； ⑥中國石化勝利油田分公司物探技術研究院，山東東營 257022)

0 引言

巖性識別是儲層評價的一個重要環節，是油藏描述、實時鉆井監控及求取儲層參數的基礎。直接對巖心進行實驗測量是識別巖性最準確的方法，但時間與經費消耗巨大，在實際應用中受到限制。測井曲線包含豐富的儲層信息，如何快速、低耗利用測井資料獲取地層巖性信息越來越受到研究人員的重視。Busch等[1]、Doveton[2]提出基于測井曲線交會的地層巖性經典識別技術，并迅速得到推廣應用。然而隨著勘探難度的不斷加大，交會圖法暴露出較大的局限性，許多學者提出了新的解決方案，巖性識別技術得到了長足的發展。劉愛疆等[3]將多元統計學中的主成分分析法(PCA)應用于碳酸鹽巖巖性識別，潘保芝等[4]將因子分析法成功應用于火成巖巖性的識別，但都是針對特定巖性； Samuel等[5]、吳磊等[6]利用神經網絡對巖性進行了有效識別，但由于黑盒效應，無法控制與分析模型的中間過程； 張翔等[7]、吳施楷等[8]利用支持向量機進行巖性分類，但支持向量機需要尋找最優的分界樣本，樣本的選擇不當會影響識別效果； 葛新民等[9]借助核函數方法和時頻分析理論，建立了一套基于核主成分分析(KPCA)與小波能譜分析的復雜儲層油水界面的識別方法； 陳鋼花等[10]利用粒子群算法對支持向量機的核函數進行優化，提高了對復雜砂礫巖體的識別精度，驗證了粒子群算法尋優性能。

東營凹陷董集洼陷的濁積巖類型多樣、縱向巖性變化大、次生成巖作用改造明顯，具有極強的非均質性，測井參數之間非線性關系嚴重，傳統的交會圖法及常規主成分分析法均不能得到精確的識別結果，是該區的研究難點[11-12]。

為解決該難題，綜合前人的研究思路，采用核主成分分析方法，利用基于動態慣性權重系數的粒子群算法構建合理的核函數變換矩陣，通過核函數把原始測井數據從輸入空間投影到高維線性特征空間，在高維線性特征空間中進行主成分分析，選取能夠代表大部分測井信息的核主成分進行交會分析，建立了一種快速有效的濁積巖巖性識別方法。

1 方法原理

1.1 主成分分析法原理

主成分分析是在保證樣本數據信息損失最小的前提下，運用降維的思想將多個指標問題轉換為較少的新的指標問題，并且這些新的指標既互不相關，又能綜合反映原指標所包含信息的一種分析方法[13-14]。PCA運算實際上是一種確定一個坐標系統的正交變換，在這個新的坐標系統下，變換數據點的方差沿新的坐標軸得到了最大化，這些坐標軸經常被稱為是主成分[15]。PCA的具體計算過程如下。

(1)

為避免不同變量的數量級和量綱的影響，需要將原始數據標準化變換為X。變換公式為

(2)

(2)計算經過標準化變換后的樣本矩陣相關系數矩陣

(3)

式中

(4)

(3)利用雅可比變換求出相關系數矩陣R的特征方程非負特征根，并按大小進行排列λi(λ1>λ2>…>λn>0)。通過施密特正交化方法單位正交化其特征向量

(5)

(4)計算樣本的主成分變量，以表示原樣本的線性組合

F=AX

(6)

即

(7)

(5)利用特征值累積貢獻度確定m的值，即當前m個主分量的方差和占全部總方差的比例超過0.85時，選取前m(m

1.2 核主成分分析法原理

利用主成分分析可以較好地處理變量間的線性關系，但處理非線性關系時會導致各主成分的貢獻率過于分散，不能找到能夠有效代表原樣本的綜合變量，處理效果較差[17]。

核方法是一類有效處理非線性問題的方法，通過非線性映射Ф，將低維輸入空間每一個X=(X1，X2，…，Xp)(Xi∈Rn，i=1，2，…，p)中不可分的數據映射到高維特征空間Y，即

(8)

在高維特征空間中進行數據處理，使輸入空間中不可分的數據在高維特征空間中變得可分，再在高維特征空間采用主成分分析進行特征提取，解決原樣本空間中樣本數據線性不可分的問題[18-19]。

圖1為利用核主成分分析解決非線性問題的效果。如圖所示，樣本一和樣本二在低維的輸入空間中線性不可分，而通過簡單的非線性映射函數Ф，可在高維的特征空間中實現樣本分離。

圖1 KPCA實現樣本分離的原理示意圖

利用核主成分分析計算過程如下。

假設待處理的數據已經零均值化，在高維特征空間Y中利用協方差矩陣

Ф(Xj)ФT(Xj)

(9)

構建特征方程尋找非零特征值λ和特征向量v∈Y，即

λv=Cv

(10)

根據Mercer條件，存在如下的核函數矩陣[20]

K(Xi,Xj)=ФT(Xi)Ф(Xj)=Ф(Xi)·Ф(Xj)

(11)

且在空間Y中任一特征向量都可以由該空間中的所有樣本線性表示，則

Ф(Xi)

(12)

式(10)可變換為

(13)

等式兩邊同時乘以ФT(Xk)，并聯合式(9)可得

(14)

即

(15)

式中：K為p×p維的核函數矩陣；α為包含α1、α2、…、αp的列向量。

由于K是滿足Mercer條件的實正定函數，所以式(15)可化簡為

pλα=Kα

(16)

上式的非負特征根可以利用雅可比變換求出，并通過施密特正交化方法單位正交化其特征向量。

對于輸入空間中某一個樣本Xj(j=1，2，…，p)，利用非線性映射得到特征空間Y的像Φ(Xj)，那么它在特征空間Y中的主成分可以表示為在子特征空間向量上vk的投影

(17)

然而，利用核方法的內積運算把樣本數據從低維空間映射到高維特征空間，雖然方便了非線性特征的提取，但是隨著維數的增加，計算量也隨之增加，甚至會出現“維數災難”。所以要利用核函數取代內積運算來避免“維數災難”[21-22]。

目前常用的有多項式核函數、徑向基(高斯)核函數、神經網絡核函數。

(1)多項式核函數為

K(Xi,Xj)=[b(Xi,Xj)+c]d

(18)

式中：b、c為常數系數；d為多項式核的階次。

(2)徑向基(高斯)核函數為

(19)

式中σ是參數，它決定了輸入變量在學習算法中被縮放的程度。

(3)神經網絡核函數為

K(Xi,Xj)=tanh[s(Xi,Xj)+c]

(20)

式中: tanh為雙曲正切函數；s、c是根據先驗知識指定的參數。

將內積用核函數替換,則式(17)變為

(21)

所以借助核函數就可以在高維特征空間中對樣本進行主成分分析[23]。

上述推導過程均是在Ф(X)進行中心化處理后進行的，但在實際應用中并不知道Ф(X)的顯式表達，也就無法對其中心化，所以要使用對角化后的Kc來代替K進行上述求解過程。Kc的表達式為[24]

Kc=K-lpK-Klp+lpKlp

(22)

式中：lp為p×p的矩陣，其每一個元素都是1/p。

1.3 粒子群算法優化原理

粒子群算法的基本思想是通過群體中個體之間的協作和信息共享尋找最優解，是基于對鳥群覓食行為的生物仿真模擬，目前已被廣泛應用于函數優化、模糊系統控制等領域[25]。

在粒子群算法中，慣性權重w是最重要的參數，較大的w有利于提高算法的全局搜索能力，而較小的w會增強算法的局部搜索能力。為了平衡傳統粒子群算法的全局搜索能力和局部尋優能力，本文采用自適應的動態慣性權重系數公式，具體為

(23)

式中：wmax、wmin分別表示w的最大值和最小值；f為粒子當前的目標函數值；favg和fmin分別表示當前所有粒子的平均和最小目標函數值。

基于自適應的動態慣性權重粒子群算法簡單、容易實現并且沒有太多參數的調節，適用于核主成分分析中的核函數參數的優選[26]。

2 實例研究

2.1 地質概況

董集洼陷位于西段東營凹陷北部，在沙三中段沉積時期，董集洼陷接受的南北兩個方向的物源供給，發育多套濁積巖，這些濁積巖被沙三中段的深湖相暗色泥巖包裹，形成自生自儲的巖性油藏。根據研究區的巖心資料，沙三中段沉積時期的巖石類型主要為深灰色的泥巖、灰質泥巖、淺灰色的細砂巖、灰色的灰質砂巖、褐色的油頁巖。研究區的濁積巖主要包括砂巖與泥巖、砂巖與灰質泥巖、灰質砂巖與泥巖、灰質砂巖與灰質泥巖及砂巖與油頁巖等5種巖性組合方式[27]。

2.2 樣本參數選擇

測井曲線包含多個參數，每個參數對巖性都有不同的區分度，但若選取的測井參數過多，往往會增加分析過程的復雜性。基于前人經驗和測井響應原理，選取自然伽馬(GR)、聲波時差(AC)、補償中子孔隙度(CNL)、密度(DEN)、原狀地層電阻率(RT)值作為輸入樣本參數X′=[GR，AC，CNL，DEN，RT]。

選擇研究區W9井的沙三中段4層(Es3z4)的測井數據作為研究樣本，進行單位長度為0.125m的數據重采樣，得到230個樣本，對樣本數據進行異常值處理，去除無效參數，最終確定用于分析的樣本個數為203。不同巖性測井數據的均值和中值統計結果如表1所示。

表1 不同巖性測井數據統計結果

根據測井數據統計可知，砂質泥巖、灰質泥巖的測井曲線具有高GR、高AC、中CNL、中DEN、中RT即“兩高、三中”的響應特征；細砂巖、粉砂巖表現為低GR、低AC、低CNL、高DEN、高RT即“三低、兩高”的特征；泥巖顯示為高GR、高AC、高CNL、低DEN、低RT即“三高、二低”特征。通過對不同巖性測井值的統計，發現五種巖性的測井響應特征存在一定差異，適合采用多元統計學方法進行分析[28]。

2.3 核函數及其參數的選擇

為了選取合適的核函數，對樣本進行試驗。試驗發現，選取神經網絡核函數來構造KPCA核變換矩陣時，對參數的選擇較為苛刻，易造成病態核矩陣，生成負的特征向量從而導致模型構建失敗；而徑向基核函數和多項式核函數的正定性較好，參數的選擇區間較廣，容錯率較高。

對徑向基核函數和多項式核函數運算性能進行進一步的試驗對比(將迭代進程歸一化為0和1，1為開始，0為結束)，對比結果如圖2所示。從圖中可以看出，相對于多項式核函數，徑向基核函數只用了較少的迭代次數和運算時間就完成了迭代進程，運算效率高，且函數需要確定的參數只有一個，模型構建較為容易。通過以上試驗分析，選取徑向基核函數作為核主成分分析的核函數。

圖2 兩種不同核函數的運算效率對比

為了確定徑向基核函數的核參數σ，避免人工選擇的誤差以及提高識別效率，采用基于動態慣性權重的粒子群算法對σ進行優選。粒子群的運動軌跡如圖3所示，可以看出，在迭代初期粒子擴散到全區在較大范圍內搜索最優解，隨著迭代次數的增加，慣性權重也在發生改變，整個粒子群不斷向最優的位置靠近聚集，在迭代后期鎖定了全局最優解的位置，得到核函數參數的最優解大約為0.6(圖中紅圈處)。

為了驗證該方法選取σ的可行性，分別選取不同的σ值對標準井中較難識別的三種巖性進行識別，結果如圖4所示。從圖中可以看出，隨著σ取值的改變，三種巖性的識別準確率也在改變，在σ值為0.6時，三種巖性的識別準確率均達到最高，驗證了該方法的準確性。

圖3 粒子群的運動軌跡

圖4 σ的取值對識別結果的影響分析

2.4 核主成分分析測井巖性識別方法步驟

(1)數據的標準化預處理。對選取的p個測井數據(每個測井數據包含n個參數)進行標準化預處理，克服數據間量綱和數量級的差別，以保證分析結果的客觀性和準確性。

(2)核函數參數優選。基于動態慣性權重的粒子群算法對σ進行優選，避免人工選擇的誤差以及提高識別的準確率。

(3)計算并修正核矩陣。利用確定最優核參數的徑向基核函數構建核矩陣K，并將核矩陣進行中心化修正得到Kc。

(4)高維空間中核矩陣分解。在選定的核函數的作用下，通過雅可比迭代計算修正后的核矩陣Kc的特征值和對應的特征向量，并利用施密特正交化方法單位正交化其特征向量，得到α1、a2、…、αn。

(5)計算特征值的累積貢獻率B1、B2、…、Bm,根據給定的提取效率h，如果Bm≥h，則提取m個特征分量。

(6)高維空間中的主成分提取。計算已修正的核矩陣在提取出的特征向量上的投影，得到m個主成分。

(7)主成分交會識別巖性。根據得到的前兩個主成分進行交會分析，確定每種巖性在主成分坐標系中的分布區間[29]。

2.5 計算過程

在對輸入參數X′進行分析之前，先對樣本數據進行標準化，使每個測井參數變量的均值為0，方差為1，標準化后的樣本為X。

對標準化后的樣本數據分別進行主成分分析和核主成分分析(流程如圖5所示)，并對得到的各主成分的特征值及方差貢獻率進行比較(表2)。

從表2中可以看出，經過主成分分析計算得到的第一主成分F1的方差貢獻率只有40.13%，貢獻率并不理想； 第二主成分F2的方差貢獻率也僅有21.44%，主成分F1和F2方差之和只占到總方差的61.57%，原有變量的信息損失較大。經過核主成分分析計算后僅第一主成分F1的方差貢獻率就達到86.31%，第一主成分F1和第二主成分F2方差貢獻率之和占到總方差93.83%，有效地保留了原有樣本變量的信息。

表2 PCA 和 KPCA 的特征值貢獻率比較

圖5 基于核主成分分析的巖性識別流程

3 應用效果及實例分析

圖6為常規測井交會圖法對研究區的濁積巖識別結果，可以看出，該方法的識別精度較低，很難區分所有巖性。

圖7a為由主成分變換方程計算的F1和F2交會圖，圖7b為由核主成分變換方程計算的F1和F2交會圖。可以看出，主成分分析法對細砂巖和粉砂巖區分度很高，但對泥巖、灰質泥巖、砂質泥巖的區分度并不理想；而核主成分分析可以有效區分泥巖、灰質泥巖、砂質泥巖、粉砂巖和細砂巖，巖性識別精度明顯提高。

分別運用核主成分分析法和主成分分析法對研究區的W9井進行連續測井剖面的巖性解釋，結果如圖8所示。從圖中可以看出，在量綱一致的條件下，KPCA法的F1、F2分量曲線垂向上變化程度更加劇烈，不同巖性之間的差異性更加明顯，對于巖性的區分更加敏感。根據實際錄井資料，W9井2600～2670m深度區間里共出現20個不同的巖性段。PCA法可以有效識別細砂巖和粉砂巖段，但是對泥巖、灰質泥巖、砂質泥巖段的誤判率很高，20個巖性段誤判8個(見圖中紅色圓圈處)，巖性識別率僅為60%； KPCA法解釋結果與實際錄井資料較為吻合，只是將2622～2624m和2627～2629m處的砂質泥巖誤判為灰質泥巖，20個巖性段僅誤判2個(圖中藍色方框處)，巖性識別準確率達到90%，在濁積巖地區進行巖性識別效果良好。

圖6 常規測井交會圖法巖性識別

圖7 兩種方法主成分巖性交會識別對比

圖8 研究區W9井兩種方法測井巖性解釋結果對比

圖9 兩種方法儲層平面預測結果對比

對區內的多口井應用KPCA法和PCA法進行巖性識別后預測了該區的沙三中段4層有利儲層的分布，結果如圖9所示。研究區砂巖的F1值分布區域大致為-20～-12和40～48，從圖中可以看出，KPCA法預測結果在鉆井處與實際測井資料吻合度更高，更精細地刻畫了儲層的分布范圍，符合該區儲層的分布規律。

4 結論

(1)東營凹陷董集洼陷濁積巖巖性復雜多變，基于粒子群算法以及核函數理論，采用核主成分分析法能有效地從多種測井參數中提取突出巖性特征的主成分分量，實現了降維處理。

(2)利用核主成分分析法算出的第一主成分和第二主成分變量進行交會分析，能夠有效地識別泥巖、灰質泥巖、砂質泥巖、細砂巖和粉砂巖，巖性識別精度明顯高于常規測井交會圖法和主成分分析法。

(3)多口井實際資料的識別結果與錄井資料符合程度較高，為后續的儲層評價奠定了基礎。