多閾值BIRCH聚類在相控孔隙度預測中的應用

孫歧峰 段友祥 柳 璠 李洪強

(①中國石油大學(華東)計算機科學與技術學院，山東青島 266580；②中國石化勝利油田工程有限公司鉆井工藝研究院，山東東營 257000)

0 引言

孔隙度是儲層評價、油氣預測的重要參數[1-2]。孔隙度不僅受控于巖石的孔隙結構，還與巖石的顆粒組分密切相關[3]，這為利用巖相作為約束條件準確預測孔隙度提供了可靠依據[4-7]。

許多學者以巖石物理特征為基礎，利用地震反演方法預測儲層參數[8-10]。隨著人工智能和機器學習技術的快速發展，多種分類、回歸、聚類算法也被應用于儲層參數預測。早期學者多利用統計學方法，如樸素貝葉斯分類器[11]、 模糊邏輯[12]、 支持向量機[13]等分析測井相。在地震屬性數據基礎上，張陽等[14]提出模糊C均值地震屬性聚類的沉積相分析方法，龐銳等[15]利用K均值聚類進行地震相識別。在儲層參數預測中，應用較為廣泛的是神經網絡[16-17]算法；其他的回歸算法，如支持向量回歸(SVR)[18]、 K近鄰法[19]、 隨機森林[20-21]等也被應用于儲層預測。

由Zhang等[22]提出的BIRCH(Balanced Iterative Reducing and Clustering Using Hierarchies)方法是一種基于距離的無監督學習層次聚類算法，計算效率較高。邵峰晶等[23]提出了可變多閾值聚類(Multi-threshold BIRCH Clustering，M-BIRCH)方法實現了對任意形狀簇的聚類。本文在考慮計算機系統內存以及聚類質量等因素，通過啟發式設定閾值，并基于Agglomerative算法[24]的全局聚類進一步優化該算法，以井點處孔隙度數據和地震聲波阻抗屬性作為輸入，確定巖相類型；然后在同一相帶內使用嶺回歸(Ridge Regression,RR)方法預測儲層孔隙度，以提高預測精度。

1 多閾值BIRCH聚類分析

1.1 BIRCH聚類原理

BIRCH算法是將數據集中的樣本劃分為多個不相交的類，并且這些類具有類內對象相似、類間對象相異的特點。每個類可表示一種相，稱為“簇”。基于距離度量的算法能夠發現具有相近尺度和密度的球狀簇。

BIRCH算法可通過聚類特征(CF)對地震屬性進行聚類。以波阻抗數據為例，假設在某個簇中存在N個波阻抗樣本{xi}，i=1，2，…，N，將聚類特征定義為三元組

CF=(N， LS， SS)

(1)

聚類特征樹(CF樹)用于存儲層次聚類的簇的特征，是一棵高度平衡樹，即它的子樹的高度差的絕對值不超過1，并且所有子樹都是一棵平衡樹。CF樹中的非葉子節點存儲了其子節點的聚類特征的總和，其中節點的每個聚類特征和它指向子節點的指針組成條目。CF樹具有三個參數：分支因子的非葉節點數B、葉節點數L和簇半徑閾值T。其中，B是特征樹中每個非葉節點最多包含的子節點的數量，L是每個葉節點最多包含的CF向量的數量，T為同一葉子節點中的數據點滿足的閾值。圖1給出B=2、L=3的CF樹結構示例。

1.2 改進的多閾值BIRCH聚類

由于BIRCH算法忽視了簇與簇之間的關系，對葉子節點中的簇設定同一閾值，因而聚類效果不穩定且有一定局限性。多閾值BIRCH聚類方法是對每個聚類特征都設定一個閾值T，即將每個簇的CF三元組改為四元組(N，LS，SS，T)，每個簇可以具有不同的閾值。在CF樹重建過程中，當兩個簇需要合并為一個簇時，新簇的閾值由合并定理求取。

圖1 CF樹結構示例

合并定理為：假定將n個簇Ci(i=1，2，…，n)合并，聚類特征CFi=(Ni，LSi，SSi，Ti)，合并后新簇W的聚類特征為

(2)

式中：Ci.mean、W.mean分別表示簇Ci、新簇W的質心； dist(W.mean，Ci.mean)為簇間距離，計算公式為

(3)

式中： SS1、LS1、N1表示簇C的三元組； SS2、LS2、N2表示簇W的三元組。

根據式(2)，通過簇的聚類特征計算簇間距離，動態調整、增加閾值T。閾值的提升使原始CF樹中的條目進一步合并，進而重建一棵較小的CF樹。調整后的閾值不僅包括原有簇的所有波阻抗樣本，而且大大增加了新的波阻抗樣本插入的可能性，最終產生基于不同閾值的聚類結果。

2 嶺回歸原理

回歸分析[26]是確定兩種或兩種以上變量間相互依賴的定量關系的一種統計分析方法。嶺回歸[27]是一種通過放棄最小二乘法的無偏性，以損失部分信息、降低精度為代價獲得回歸系數的更符合實際、更可靠的回歸分析方法，它對數據擾動敏感的病態矩陣的擬合要優于最小二乘法。

λ‖β‖2}

(4)

式中：λ是調節參數；Rp表示p維實數集。

嶺回歸的代價函數是一個凸函數，同最小二乘法一樣，當代價函數梯度(即一階導數)等于0時可求得全局最優解。由下式最終得到嶺回歸模型的最優系數向量

(5)

式中k為嶺參數。

3 相控孔隙度預測

相控孔隙度預測是以沉積相及其特征等作為約束條件而降低儲層孔隙度預測的多解性的一種技術。本文應用M-BIRCH方法對波阻抗數據聚類分析、劃分巖相帶，在此基礎上應用嶺回歸方法對同一相帶中的孔隙度與波阻抗建立回歸關系，最終進行孔隙度的精確預測。

3.1 關鍵技術

3.1.1 初始閾值的設定

根據M-BIRCH聚類過程，閾值是通過簇間距離進行更新的，即閾值的大小與簇間距離、類內平均距離相關，也就是說閾值的選擇與波阻抗數據分布有關。本文采用啟發式方法設定初始閾值，首先在整個數據集中隨機選取N對波阻抗樣本，將每一個樣本都作為一個簇，根據簇的CF條目，利用式(3)計算簇之間的距離； 然后計算距離(dist)的期望E(dist)和方差D(dist)； 最后設定初始閾值

T=[E(dist)+0.25D(dist)]P

(6)

式中P為預先設定的百分比。

預先定義一個較小的P，通過較小的T對波阻抗數據進行分割。分割后的數據集通過提升閾值重新歸并到較小的CF樹，并將數據裝入內存。在生成CF樹的過程中，若出現內存不夠的情況，可以通過提升P定義較大的初始閾值重建較小的CF樹。計算初始閾值T后，就可以建立一棵初始CF樹。

3.1.2 動態構建CF樹

CF樹的構建是BIRCH聚類算法的關鍵。隨著波阻抗數據的不斷加入，將不斷地重建CF樹。波阻抗數據按照與各個CF 的距離依次插入到各條目中合并為簇，這樣距離相近的波阻抗將處于相同的簇中。不同的巖相帶波阻抗存在差異，重建CF樹能夠不斷地區分波阻抗數據，從而完成相帶劃分的目的。把每個波阻抗數據都看作為獨立的簇，動態建立CF樹的過程也就是將簇的聚類特征CF插入到CF樹的過程，具體如下。

(1)從根節點開始遞歸向下，通過計算CF與插入節點包含的各條目中CF的距離，尋找距離最短的路徑及葉節點。

(2)如果CF與該葉節點各條目中的CF之間的距離小于閾值T，則選擇閾值最小的條目，用合并算法將CF與該條目的CF合并，并自下向上相應地修改各節點的條目。

(3)如果CF無法與該條目的CF合并，則判斷該條目所在葉節點的CF數量是否小于L。若小于L，則將CF新建為一個條目，并將該條目插入到距離最近的條目后面，相應修改CF樹的結構。否則分裂該葉節點，原則是以距離最遠的兩個條目為種子進行分裂，剩余的條目按照距離最近原則合并到這兩個條目中，并更新整個CF樹。

3.1.3 全局聚類

對CF使用聚類算法進行全局聚類，能夠得到較好的聚類效果[22]。本文使用層次聚類算法中的凝聚法，即Agglomerative算法[24]。這是一種自底向上的方法，不需要求取目標函數，即不存在局部最小和對初始點敏感的問題。其算法描述如下。

輸入：阻抗樣本、聚類數目

輸出：聚類結果

將每個阻抗樣本都作為一個獨立的簇；

重復：

計算兩兩簇之間的距離，找出距離最小的兩個簇c1和c2；

將c1和c2合并為一個簇；

直到達到聚類數目或其他設定的條件結束。

3.1.4 優化嶺參數

當自變量之間存在線性關系時，自變量相關矩陣行列式近似為0，導致最小二乘估計失效。在嶺回歸訓練過程中，通過選取合適的嶺參數，可以減少這種多重共線性的影響。當X存在多重共線性時，矩陣XTX的特征值至少有一個接近于0，可用XTX行列式的值表示多重共線性的嚴重程度。因此，令

(7)

式中I為p×p的單位矩陣。由式(7)可知，D(k)隨著k的增加而增大，當00.05時，幾乎不具有共線性。因此，認定在D(k)>0.01的范圍內將最小的k值作為最終的嶺參數[28]。

3.2 相控孔隙度預測流程

在以地震、測井資料為基礎的儲層孔隙度預測中，首先需要對數據進行預處理，即對地震、測井數據反演得到波阻抗數據，通過巖心分析及測井解釋得到井點孔隙度數據；然后將波阻抗和孔隙度數據進行歸一化處理，使其處于同一數量級中，以便于分析。

M-BIRCH使用聚類特征CF表征一個簇，即劃分巖相，并采用平衡樹結構，能夠滿足數據的可伸縮性，處理大規模數據集速度較快。利用M-BIRCH聚類算法劃分巖相可以歸納為三個步驟： 第一步是掃描波阻抗數據，確定初始閾值，構造初始CF樹，在CF樹構造過程中，首先要確定兩個參數，即非葉子節點數B和葉子節點數L，CF樹的大小可以通過調節這兩個參數進行調整； 第二步是隨著波阻抗數據的不斷加入，重建CF樹，并且動態更新閾值，在重建過程中，只通過葉節點重新構建一棵新樹，并不需要重新訪問所有數據[23]； 第三步是用CF代替原有數據集進行聚類。當CF樹構建完成之后，選用Agglomerative算法對葉節點中的簇進行全局聚類。

在相帶的約束下，利用井點處孔隙度與波阻抗數據進行非井點處的孔隙度預測。在每個相帶內都會由嶺回歸計算得到最優參數估計量，在此基礎上，由未標簽波阻抗數據可計算得到孔隙度預測值。

圖2 相控孔隙度預測流程

利用交叉驗證方法將井點處預測的孔隙度值與已知的孔隙度進行誤差分析。本文采用交叉驗證方法，即選擇w-1口井的波阻抗和孔隙度數據作為訓練集進行模型訓練，一口井作為測試集驗證模型及參數。這樣循環w次，使每口井都作為測試集進行檢驗，最后選擇損失函數評估最優的模型及參數。

3.3 算法效率分析

4 應用實例

4.1 合成地震數據孔隙度預測

Stanford VI-E模型[29]包含三層前積型河道，三維地震數據具150×200×200個網格，每個網格的尺寸在水平(x和y)方向上均為25m，在垂直(z)方向上均為1m。本文實驗采用第二層，即網格數為150×200×40。為了簡化實驗數據，將該系統中的點壩和河道砂體合并為砂巖相(圖3c中的黃色部分)，將洪泛平原和邊界合并為頁巖相(圖3c中的藍色部分)。

圖3a為水平切片(網格數為150×200)，即輸入數據。圖3b顯示孔隙度與波阻抗之間存在著負相關關系。但是，其相關系數和偏移量很大程度上依賴于不同巖相的不連續性，并且孔隙度在每個巖相帶內具有低可變性。因此，一旦巖相結構已知，便可通過回歸方法根據波阻抗估計孔隙度。

圖3 本文簡化模型

4.1.1 多閾值BIRCH+RR方法預測

實驗中，將15%的輸入數據作為訓練集，85%數據作為測試集。設定用于計算初始閾值的百分比P為0.1、非葉節點數B和葉節點數L均為30。圖4為圖3a作為輸入的巖相預測。與圖3c相比，M-BIRCH比BIRCH[22]、K-均值(K-means)[15]方法能更準確識別河道。圖4a左上方實際為頁巖相，但BIRCH、K-means方法卻將其劃分為砂巖相。

本文使用調整蘭德系數(ARI)[30]進行量化描述。ARI用來描述數據分布的吻合程度，計算公式為

(8)

(9)

圖4 不同方法的巖相預測

如表1所示，M-BIRCH算法的ARI為0.95，明顯優于另外兩種算法的聚類效果。

表1 不同方法預測相帶的性能比較

圖5為不同方法預測的孔隙度。圖5c為使用SVR[31]回歸方法進行的預測，即考慮只存在一種相帶的情況。對比圖3d，可以看出圖5c中的孔隙度低值(藍色部分)并沒有擬合到真實的結果，因此相控法孔隙度預測結果比SVR方法預測結果更符合實際情況。

表2為三種方法孔隙度預測的性能比較。M-BIRCH+RR(本文方法)預測結果的均方誤差(MSE)以及平均絕對誤差(MAE)最小，說明本文方法比另外三種方法預測精度高。從確定性相關系數(R2)也可得知，本文方法預測結果與實際結果匹配程度更高。另外，使用SVR方法預測結果的MSE達到0.92288，相對于M-BIRCH+RR、K-means+RR方法，SVR方法可行性不高，誤差較大。從方法的運行效率來看，SVR方法雖然運行時間較短，但是最后結果不理想； K-means+RR方法的運行時間最長，但是精度不是最高； 與BIRCH+RR相比，M-BIRCH+RR方法雖然犧牲時間效率，但效果更好(R2值高，MSE值低)。

圖5 不同方法預測的孔隙度

表2 不同方法孔隙度預測性能比較

4.1.2 靈敏度分析

分別向波阻抗數據中輸入0～80%的高斯噪聲，利用M-BIRCH+RR方法預測巖相，最終得到孔隙度預測結果如圖6所示。當噪聲不斷增加時，對巖相劃分結果影響不大，M-BIRCH方法仍能較準確地區分出砂巖相和頁巖相。圖7a是在不同的高斯噪聲下中位絕對誤差(MDAE)的變化，可以看出MDAE隨著高斯噪聲線性增加。圖7b是在不同的高斯噪聲下確定性相關系數的變化，R2隨著高斯噪聲的增加不斷降低，即使高斯噪聲達到80%，M-BIRCH+RR方法預測結果的R2仍能達到75%左右。因此，即使存在一定的噪聲，仍然不影響其性能，說明本文方法具有較高穩定性和抗干擾能力。

4.2 實際地震數據預測

研究區主要為三角洲前緣沉積，分支河道砂體發育。目的層埋深大，壓實作用強，巖性致密。研究區內共有8口井，有效儲層預測難度大。針對上述難點，采用本文方法進行相控孔隙度預測。采用8口探井的巖心及測井資料作為標記樣本，儲層數據網格的尺寸在水平(x和y)方向上均為20m。選取其中一個時間切片(網格數為653×600)的波阻抗樣本，如圖8所示。預測過程中，首先對波阻抗數據進行網格化處理，將井點處的孔隙度數據應用到離井點位置最近的網格點中，并對井點處孔隙度和阻抗數據進行歸一化； 然后采用本文方法進行巖相劃分和孔隙度預測； 最后采用交叉驗證方法對預測結果進行誤差分析，即利用7口井作為訓練數據，剩余1口井作為驗證數據。

圖6 不同噪聲情況下波阻抗切片(左)和本文方法預測的巖相(中)、孔隙度(右)

圖7 本文方法靈敏度分析

圖8 波阻抗切片

預測結果如圖9所示。從圖9a可以看出，M-BIRCH方法根據波阻抗數據將儲層劃分為兩類，黃色部分對應原始數據的高波阻抗值，藍色對應低波阻抗值，較好地刻畫了河道砂體的分布。孔隙度預測結果如圖9b所示。利用稀疏標記訓練預測的孔隙度能夠反映河道砂體橫向展布規律，井點處結果與波阻抗分布、砂巖發育位置均有很好的對應關系，驗證了本文方法的合理性。

由三種方法在研究區的實驗結果性能比較(表3)可知，K-means+RR的性能最差，均方根誤差最高，可靠性較差。相對而言，M-BIRCH+RR方法的MSE較BIRCH+RR方法、K-means+RR方法略小，說明M-BIRCH+RR方法是可靠的。

圖9 實際地震數據預測結果

表3 不同方法的性能比較

5 結束語

本文考慮不同巖相之間波阻抗數據的差異，提出一種多閾值BIRCH聚類的無監督學習方法，并引入Agglomerative算法進行全局聚類，得到準確的巖相預測結果； 然后在巖相約束下以井點處的波阻抗和孔隙度為輸入，采用優化的嶺回歸方法在同一巖相內進行孔隙度預測。模型數據實驗表明，與其他劃分巖相的傳統聚類方法相比，M-BIRCH方法受異常數據的干擾較小、適應性更好、運行效率更高。實際數據應用表明，使用M-BIRCH算法與嶺回歸結合的方法對孔隙度預測能夠得到合理且準確的結果。