核磁共振橫向弛豫時間譜高斯混合聚類及應(yīng)用

葛新民，薛宗安，周軍，胡法龍，李江濤，張恒榮，王爍龍，牛深園，趙吉兒

（1. 中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，山東青島266580；2. 中國石油天然氣集團(tuán)有限公司測井重點(diǎn)實(shí)驗(yàn)室，西安710077；3. 海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價與探測功能實(shí)驗(yàn)室，山東青島266071；4. 中國地質(zhì)調(diào)查局油氣資源調(diào)查中心，北京100083；5. 中國石油集團(tuán)測井有限公司，西安710077；6. 中國石油勘探開發(fā)研究院，北京100083；7. 中國石油青海油田公司，甘肅敦煌736202；8. 中海石油（中國）有限公司湛江分公司，廣東湛江524057；9. 中國石化石油工程技術(shù)研究院，北京100101）

0 引言

核磁共振測井在復(fù)雜及非常規(guī)油氣儲集層的物性參數(shù)計(jì)算、孔隙結(jié)構(gòu)評價和流體識別中發(fā)揮著重要作用。由于弛豫信號僅與孔隙流體的氫核有關(guān)，通過測量多孔介質(zhì)中氫核的宏觀磁化矢量隨時間的衰減關(guān)系，并將其反演成橫向弛豫時間（T2）譜，結(jié)合相關(guān)模型可計(jì)算出孔隙度、滲透率、束縛水飽和度、孔隙半徑、潤濕性等參數(shù)，T2譜的形態(tài)與分布可定性表征孔隙結(jié)構(gòu)及儲集層級別[1-3]，縱向弛豫時間、擴(kuò)散系數(shù)、內(nèi)部梯度等多維度數(shù)據(jù)可定量表征孔隙流體類型、含量及黏度等[4-8]。與自然伽馬、聲波時差、補(bǔ)償密度等常規(guī)測井相比，核磁共振測井T2譜中蘊(yùn)含著更加豐富的地質(zhì)信息。

受國外技術(shù)封鎖和現(xiàn)場作業(yè)條件限制，國內(nèi)大部分油田的核磁共振測井資料以一維T2譜為主，因此，更有效地利用核磁共振測井T2譜來解決油田實(shí)際問題尤為重要。國內(nèi)外學(xué)者針對核磁共振測井資料的應(yīng)用開展了大量研究，在孔隙度、滲透率、束縛水飽和度、T2截止值等參數(shù)計(jì)算和弛豫時間與孔隙半徑的轉(zhuǎn)換等方面取得豐碩成果[9-12]，并在截止值、平均值、孔徑區(qū)間占比和偽毛管壓力曲線轉(zhuǎn)換等傳統(tǒng)方法基礎(chǔ)上開展了創(chuàng)新研究，提出諸如基于T2譜集中度的孔隙結(jié)構(gòu)評價、基于數(shù)學(xué)形態(tài)學(xué)的T2譜量化分解、基于分形理論的T2譜分維數(shù)提取、基于多元高斯函數(shù)擬合的T2譜參數(shù)提取方法、基于自組織神經(jīng)網(wǎng)絡(luò)與獨(dú)立分量分析和非負(fù)矩陣分解等機(jī)器學(xué)習(xí)的核磁共振測井T2譜定量評價技術(shù)[13-24]，為儲集層孔隙結(jié)構(gòu)評價提供了新的思路。但是，上述方法仍難以描述孔隙結(jié)構(gòu)與核磁共振響應(yīng)的非線性映射關(guān)系，在數(shù)據(jù)處理過程中需要人為調(diào)整參數(shù)，需要實(shí)驗(yàn)數(shù)據(jù)來刻度。

為克服當(dāng)前技術(shù)方法存在的局限性，本文提出基于高斯混合模型的核磁共振測井T2譜無監(jiān)督聚類算法，并利用數(shù)值模擬和核磁共振測井資料對方法的有效性和適用性進(jìn)行驗(yàn)證。

1 高斯混合模型聚類原理

高斯混合模型（GMM）是指有限個獨(dú)立的多元高斯分布模型的線性組合，是一元高斯分布在高維空間的擴(kuò)展。GMM是學(xué)習(xí)速度最快的概率模型，通過擬合輸入數(shù)據(jù)集構(gòu)建合適的混合多維高斯分布模型，達(dá)到無監(jiān)督聚類的目的，提供合理高效的數(shù)據(jù)分布規(guī)律與劃分方案。基于 GMM 的聚類方法在圖像分割、信號處理、生物醫(yī)學(xué)和地球科學(xué)等領(lǐng)域取得重要應(yīng)用[25-31]。根據(jù)高斯混合模型的基本概念[26-29]，測井?dāng)?shù)據(jù)體的概率密度函數(shù)可表示為（1）式—（3）式。

GMM聚類需求解出每個高斯基函數(shù)的均值、標(biāo)準(zhǔn)差和權(quán)重矩陣，一般采用期望最大化（EM）算法。EM算法的基本思想是通過引入隱含變量求解高斯混合模型分布函數(shù)的最大似然估計(jì)，再對隱含變量期望公式和模型分布參數(shù)重估公式進(jìn)行迭代，直至似然函數(shù)值收斂[26,31-33]。為降低計(jì)算復(fù)雜度，對似然函數(shù)取對數(shù)，如（4）式。求解高斯混合模型的參數(shù)可轉(zhuǎn)化為求對數(shù)似然函數(shù)的極大值，如（5）式。

EM算法分兩步，第1步求期望值（稱為E步），第 2步求極大值（稱為 M步），具體迭代步驟為：①初始化模型參數(shù)，設(shè)均值矩陣Mk為隨機(jī)矩陣，設(shè)方差矩陣Dk為單位矩陣，模型的權(quán)重矩陣Wk為每個模型的先驗(yàn)概率；②求出每個深度（或測量點(diǎn)）的測井?dāng)?shù)據(jù)屬于各高斯模型的先驗(yàn)概率，如（6）式；③利用先驗(yàn)概率更新高斯混合模型中的各組參數(shù)，如（7）式—（9）式；④重復(fù)步驟②和步驟③，直至滿足（10）式收斂條件。

確定合理的聚類群集數(shù)是求解高斯混合模型的關(guān)鍵，聚類群集數(shù)太大會造成過擬合，聚類群集數(shù)太小則會使模型缺乏靈活性而不能更好地擬合新數(shù)據(jù)。常用的聚類群集數(shù)確定方法是赤池信息（AIC）準(zhǔn)則。AIC準(zhǔn)則是基于信息熵的概念，選擇信息損失最小的模型[34-37]，如（11）式。

可見，AIC準(zhǔn)則是確定最佳聚類群集數(shù)的基礎(chǔ)。前人在確定最佳聚類群集的過程中發(fā)現(xiàn)，單純依賴AIC可能會出現(xiàn)單調(diào)遞減的現(xiàn)象，因此提出采用AIC變化率來確定最佳聚類群集，在一定程度上可反映AIC對聚類群集數(shù)目增加的敏感程度[38]。

本文采用AIC變化率確定最佳聚類群集數(shù)，如（12）式。當(dāng)聚類群集數(shù)目從p增加到p+1時，如果AIC變化率較大，說明聚類群集數(shù)為p描述原始數(shù)據(jù)集的聚類精度不足；反之則說明p個聚類群集和p+1個聚類群集的聚類結(jié)果相差不大；出于計(jì)算成本考慮，可認(rèn)為p為最佳聚類群集，其對應(yīng)的AIC變化率達(dá)到最小。

2 T2譜的高斯混合模型聚類數(shù)值模擬

鑒于T2譜是孔隙結(jié)構(gòu)和流體組分等巖石物理信息的綜合反映，同一個測井深度的核磁共振T2譜可得到多個弛豫時間的占比信息，符合高斯分布。因此，可以通過GMM聚類算法對T2譜進(jìn)行處理，得到無監(jiān)督聚類結(jié)果，為孔隙結(jié)構(gòu)評價和儲集層劃分提供依據(jù)。

由于T2譜數(shù)據(jù)量較大，在進(jìn)行GMM聚類之前需要通過主成分方法進(jìn)行降維以減弱數(shù)據(jù)間的相關(guān)性，并提高計(jì)算效率。首先通過數(shù)值模擬對GMM聚類方法的適用性和有效性進(jìn)行測試，在此基礎(chǔ)上對實(shí)際的核磁共振測井?dāng)?shù)據(jù)進(jìn)行處理。

2.1 T2譜主成分降維

T2譜由不同的弛豫時間所對應(yīng)幅度構(gòu)成，數(shù)據(jù)維度取決于反演布點(diǎn)數(shù)，不同弛豫時間對應(yīng)的幅度矩陣具有一定相關(guān)性。為降低T2譜數(shù)據(jù)的冗余信息，突出不同弛豫時間對應(yīng)幅值的可分析因素，采用主成分分析方法使多維數(shù)據(jù)信息矩陣映射到低維空間，提取最具代表性的主成分[38-39]。

假設(shè)核磁共振測井T2譜的布點(diǎn)數(shù)目為L，某一深度段內(nèi)共有M組觀測數(shù)據(jù)，則可以將這些T2譜表示為（13）式的觀測矩陣。由于不同弛豫時間所對應(yīng)的幅度組成的列向量存在一定的相關(guān)性，采用主成分分析可以將其變換成少數(shù)幾列線性不相關(guān)且包含最大信息量的新矩陣，提升計(jì)算效率。利用（14）和（15）式對原始觀測矩陣進(jìn)行主成分變換。

主成分分析算法的基本原理和計(jì)算方法可見相關(guān)文獻(xiàn)[38-40]，具體步驟包括：①原始數(shù)據(jù)標(biāo)準(zhǔn)化；②計(jì)算相關(guān)矩陣，得到特征向量和對應(yīng)的特征值；③根據(jù)特征值的分布計(jì)算累計(jì)方差貢獻(xiàn)率；④設(shè)置貢獻(xiàn)率截止值，提取主成分信號。

首先采用3個對數(shù)高斯函數(shù)正演400組核磁共振T2譜（見圖1），主要控制T2譜的峰值為1，10，100 ms，主要控制T2譜分布范圍的方差為0.4，0.3，0.2，弛豫時間為0.1～10 000.0 ms的對數(shù)均勻布點(diǎn)（64個點(diǎn)）。核磁共振T2譜數(shù)據(jù)的幅度為隨機(jī)生成，通過歸一化處理降低了孔隙度對計(jì)算結(jié)果的影響。當(dāng)孔隙流體均為水時，T2譜可代表微孔、小孔和中孔。統(tǒng)計(jì)結(jié)果顯示，微孔的T2譜區(qū)間為0.1～14.0 ms，小孔的T2譜區(qū)間為1.3～86.0 ms，中孔的T2譜區(qū)間為24～372 ms。為便于分析，對每組T2譜幅度進(jìn)行歸一化處理。圖2顯示了正演模擬的核磁共振T2譜數(shù)據(jù)的主成分?jǐn)?shù)與貢獻(xiàn)率、累計(jì)貢獻(xiàn)率的關(guān)系，當(dāng)主成分?jǐn)?shù)為 2時，累計(jì)貢獻(xiàn)率可達(dá)99.95%，說明不同弛豫時間對應(yīng)的幅度存在高度的相關(guān)性，在進(jìn)行高斯混合模型聚類時，僅需對前2個主成分?jǐn)?shù)據(jù)進(jìn)行處理。

圖1 對數(shù)高斯基函數(shù)T2譜

圖2 正演T2譜的主成分貢獻(xiàn)率

2.2 聚類群集數(shù)確定及結(jié)果分析

將降維后的2個主成分作為輸入數(shù)據(jù)，用EM算法進(jìn)行 GMM聚類，迭代閾值為 0.01，最大迭代次數(shù)為10 000，最大聚類群集數(shù)為9。圖3顯示不同聚類群集數(shù)的AIC及其變化率，當(dāng)聚類群集數(shù)為4時AIC的變化率達(dá)到最小值，為0.073%，說明此時擬合結(jié)果達(dá)到最優(yōu)。綜合考慮擬合精度和計(jì)算成本，確定最佳聚類群集數(shù)為4。

圖3 AIC及AIC變化率與聚類群集數(shù)的關(guān)系

計(jì)算4個聚類群集的T2幾何平均值和T2算術(shù)平均值，并將兩者進(jìn)行交會（見圖4）。從第1聚類群集至第4聚類群集，T2幾何平均值與T2算術(shù)平均值均呈增大趨勢，但第2聚類群集至第4聚類群集的T2幾何平均值和T2算術(shù)平均值的區(qū)分度低，采用T2平均值或簡單的交會圖難以有效劃分聚類群集。

圖4 不同聚類群集的T2幾何平均值與算術(shù)平均值交會圖

圖5是不同聚類群集的第1主成分和第2主成分交會圖。與圖4相比，不同聚類群集在主成分交會中的區(qū)分度更加明顯，進(jìn)一步驗(yàn)證了基于高斯混合模型的聚類結(jié)果的正確性。

圖5 不同聚類群集的第1主成分和第2主成分交會圖

聚類群集對應(yīng) 4個高斯模型的加權(quán)系數(shù)分布顯示（見圖6），聚類群集與高斯模型的加權(quán)系數(shù)對應(yīng)較好。采用 GMM 聚類不僅能得到分類結(jié)果，還能依據(jù)群集概率分配群集成員，具有更靈活的類簇形狀，得到樣本點(diǎn)隸屬于各類的概率值，能夠有效避免硬分配，所帶信息量優(yōu)于硬分配聚類[30]。

圖6 不同聚類群集的高斯混合模型加權(quán)系數(shù)分布

將上述算法進(jìn)行整合并對核磁共振測井?dāng)?shù)據(jù)進(jìn)行處理（見圖7），與傳統(tǒng)的核磁共振測井解釋剖面相比，本文方法可得到量化的高斯混合聚類結(jié)果及不同聚類群集的隸屬度值，可更加直觀地解剖T2譜及相關(guān)地質(zhì)信息。

圖7 數(shù)值模擬核磁共振T2譜的GMM聚類解釋成果

3 核磁共振測井T2譜的高斯混合模型聚類效果

柴達(dá)木盆地西北部阿爾金山前東段牛東地區(qū)侏羅系為辮狀河三角洲平原亞相沉積，巖性較粗，以中—粗砂巖為主，成分為長石巖屑砂巖或巖屑長石砂巖，填隙物含量低，主要為高嶺石和方解石；儲集空間主要為次生孔，其次為粒間孔，分布均勻且連通性較好；儲集層物性較好，孔隙度為 5%～20%，平均值為11.3%，滲透率為（0.1～100.0）×10-3μm2，平均值為5.0×10-3μm2（見圖8）。

圖8 柴達(dá)木盆地牛東地區(qū)侏羅系儲集層巖性和物性特征

XX-1井采用哈里伯頓公司的MRIL_P型儀器采集核磁共振信號，采集模式為 D9TW。對長等待時間（12.988 s）、短回波間隔（0.9 ms）的測量數(shù)據(jù)進(jìn)行處理，在主成分降維時貢獻(xiàn)率閾值設(shè)為90%，將64個橫向弛豫時間對應(yīng)的列向量壓縮至 7個主成分。對主成分矩陣進(jìn)行 GMM聚類，AIC變化率如圖9所示，結(jié)合AIC變化率準(zhǔn)則得到4個聚類群集，以0.3，10.0，100.0 ms為T2界限，得到弛豫時間為0.3～10.0 ms的相對幅度占比（P1）、弛豫時間為 10～100 ms的相對幅度占比（P2）、弛豫時間大于 100 ms的相對幅度占比（P3）（見圖10）。

圖9 XX-1井侏羅系核磁共振測井T2譜AIC變化率

圖10 XX-1井侏羅系核磁共振測井T2譜聚類結(jié)果

對核磁共振測井T2譜進(jìn)行高斯混合聚類分析，得到 4個聚類群集的不同弛豫組分占比分布顯示（見圖11），聚類群集與不同弛豫組分占比相關(guān)性較好。第 1聚類群集的P1最大、P3最小，說明小孔隙占比較大。第2聚類群集的P3最大、P1最小，說明大孔隙占比較大。第 3聚類群集的弛豫分布較復(fù)雜，無明顯占優(yōu)。第4聚類群集與第2聚類群集基本一致。

圖11 XX-1井侏羅系4個聚類群集的不同弛豫區(qū)間占比

4個聚類群集的不同弛豫組分占比分布顯示（見圖11），聚類群集與不同弛豫組分占比相關(guān)性較好。第1聚類群集的P1最大、P3最小，說明小孔隙占比較大。第 2聚類群集的P3最大、P1最小，說明大孔隙占比較大。第3聚類群集的弛豫分布較復(fù)雜，無明顯占優(yōu)。第4聚類群集與第2聚類群集基本一致。

XX-1井在侏羅系總共有3個試氣層位。2 074.5～2 082.5 m井段6 mm油嘴壓后自噴日產(chǎn)氣2 422 m3，試氣結(jié)論為低產(chǎn)氣層，GMM聚類結(jié)果顯示為第3聚類群集。2 146.5～2 153.5 m井段10 mm油嘴自噴日產(chǎn)氣13 2942 m3，日產(chǎn)油13.27 t，試氣結(jié)論為高產(chǎn)氣層，聚類結(jié)果顯示為第4聚類群集。2 228.5～2 240.5 m井段5 mm油嘴壓后自噴日產(chǎn)氣22 484 m3，日產(chǎn)油0.42 t，試氣結(jié)論為高產(chǎn)氣層，聚類結(jié)果為第4群集。2 228.5～2 240.5 m井段巖心鑄體薄片照片顯示巖石粒度相對較粗，以（含礫）巨—粗粒長石巖屑砂巖為主，以顆粒支撐為主，次棱磨圓度，點(diǎn)接觸，孔隙式膠結(jié)類型，石英表面干凈，具波狀消光，長石風(fēng)化程度深，巖屑以石英巖為主，填隙物少，孔隙較發(fā)育，儲集空間主要為粒間溶孔、粒內(nèi)溶孔、基質(zhì)微孔隙（見圖12）。2 096～2 103 m井段自然伽馬值較高的泥巖段的聚類結(jié)果為第1聚類群集。

圖12 XX-1井侏羅系2 228.5～2 240.5 m試氣層段巖心鑄體薄片照片

4 結(jié)語

依據(jù)核磁共振測井T2譜特征，提出基于高斯混合模型（GMM）聚類的T2譜處理技術(shù)，實(shí)現(xiàn)T2譜的無監(jiān)督自動聚類。采用主成分分析對T2譜數(shù)據(jù)降維壓縮，降低輸入數(shù)據(jù)之間的相關(guān)性，提高聚類算法的計(jì)算效率。采用AIC及AIC變化率確定高斯混合模型的個數(shù)，得到最佳聚類群集，提升 GMM 聚類結(jié)果的合理性。通過數(shù)值模擬和實(shí)測資料的處理解釋，驗(yàn)證了方法的有效性和可行性。T2譜的聚類結(jié)果與孔隙結(jié)構(gòu)、儲集空間類型和儲集層產(chǎn)能具有較好的對應(yīng)性，定量表征結(jié)果能夠更為直觀地反映儲集層類型和孔隙結(jié)構(gòu)。

值得注意的是，數(shù)值模擬數(shù)據(jù)并未考慮流體性質(zhì)差異對T2譜的影響，而實(shí)際儲集層中的核磁共振測井響應(yīng)往往是由巖性、流體和孔隙結(jié)構(gòu)等多種因素耦合而成。下一步將對算法進(jìn)行完善，實(shí)現(xiàn)復(fù)雜條件下流體信息和孔隙結(jié)構(gòu)的區(qū)分。對于巖性和孔隙結(jié)構(gòu)復(fù)雜的儲集層，如果能在細(xì)分巖性的基礎(chǔ)上開展核磁共振測井T2譜的GMM聚類分析，或者通過巖心實(shí)驗(yàn)、試油、生產(chǎn)動態(tài)數(shù)據(jù)來約束聚類數(shù)目，在一定程度上能提升聚類效果以及解釋的合理性。

致謝：中國石油天然氣集團(tuán)有限公司測井重點(diǎn)實(shí)驗(yàn)室運(yùn)行基金的資助，中國石油青海油田公司勘探開發(fā)研究院提供核磁共振測井資料和基礎(chǔ)地質(zhì)成果，中國石油勘探開發(fā)研究院西北分院劉建宇博士提出寶貴建議，在此一并表示感謝。

符號注釋：

AIC——赤池信息準(zhǔn)則量；d——測井?dāng)?shù)據(jù)的物理維度；Dk——第k個高斯基函數(shù)的方差矩陣；i——測井觀測數(shù)據(jù)的序號；j——迭代次數(shù)；k——高斯基函數(shù)的序號；L——核磁共振測井T2譜的布點(diǎn)數(shù)目；M——測井?dāng)?shù)據(jù)的樣本數(shù)；Mk——第k個高斯基函數(shù)的均值矩陣；N——高斯基函數(shù)的個數(shù)；p——聚類群集個數(shù)；p（X）——總體概率密度函數(shù)；Wk——第k個高斯基函數(shù)的權(quán)重矩陣；X——測井?dāng)?shù)據(jù)矩陣；Y——測井?dāng)?shù)據(jù)對應(yīng)的主成分矩陣；U——測井?dāng)?shù)據(jù)對應(yīng)的基矩陣；θ（X）——對數(shù)似然函數(shù)；θ（X）j——第j次迭代的對數(shù)似然函數(shù)；φ（X，Mk，Dk）——第k個高斯基函數(shù)的概率密度函數(shù)；ε——對數(shù)似然函數(shù)的迭代閾值；Θ（X）——對數(shù)似然函數(shù)的極大值。