基于高斯混合模型的粒度數據分析新方法及沉積相應用

雷曉虹, 尹太舉*, 趙紅靜, 王楊君, 李昂達, 李佳洋, 石放, 侯新宇

(1.長江大學地球科學學院, 武漢 430100; 2.長江大學資源與環境學院, 武漢 430100;3.塔里木油田勘探開發研究院, 庫爾勒 841000)

沉積物粒度是衡量沉積介質能量和沉積環境能量的重要指標[1]。而粒度的垂向變化,可以反映地質歷史時期沉積營力的演化特征[2]。

由于母巖風化后的產物在搬運和沉積過程中多個“動力組分”疊置[3],造成某些粒級物質優先富集[4],屬于連續有序的隨機過程的動力學總體,故沉積物粒度分布為一個復雜的多變量函數。 實驗表明[5],自然界的粒度由2個及以上次總體混合分布且符合高斯分布[6],基本呈現為雙峰或多峰(多眾數)的駝峰狀[7]。若想通過粒度反演沉積過程,對粒度分布的敏感組分進行分離和提取[8]才是有效的解決辦法。

目前最常用的粒度分析圖件為粒度概率累計曲線。Visher[9]用概率累計曲線圖上的線段數解釋搬運作用類型,但該方法存在以下問題,①根據統計學理論,當對數坐標系下的樣本累積分布如果近似于n條直線,則不能說明該樣本來自n個正態母體[10];②有很大的主觀性無法分離、提取與沉積過程相關的粒度分布次總體[11]。

為解決當前圖件存在的問題,本文引用高斯混合模型應用至粒度分析中,構建粒度次總體擬合圖。目前,高斯混合模型已被應用至水文學[12]、地球化學[13]等領域,并逐漸演變為分析不同種類異質分布的工具[14]。而該模型尚未應用至粒度與沉積學關系的研究中。

以渤海灣盆地館陶組為例,結合巖心剖面、沉積微相等資料,探尋無監督聚類的高斯混合模型能否精確識別各個粒度分布中的次總體、探尋參數的地質意義并驗證方法的有效性和適用性。

1 研究區概況

研究區所在的渤海灣盆地斷裂體系十分發育,伸展和走滑斷裂系統上下疊置、交錯,控制了盆地“東西分帶、南北分塊”的基本構造格局。其處于渤海最大的富生烴洼陷--渤中洼陷東南部渤南凸起高部位[15],經歷了印支運動、燕山運動及喜馬拉雅運動的演化和改造[16](圖1)。

圖1 PL19-3油田區域構造位置圖Fig.1 Regional structural location map of PL19-3 oilfield

館陶組沉積時期整體屬于亞熱帶季風氣候,溫暖濕潤[17]。早期主要受東部膠遼隆起控制,沉積類型為辮狀河;館陶組沉積中晚期至明下段時期主控物源轉變為遼東隆起,沉積類型逐漸由近緣的辨狀河沉積演變為遠源曲流河沉積,屬于辮狀河與曲流河的過渡相類型[18]。館陶組下部為厚層狀灰色砂礫巖、含礫砂巖、砂巖夾紫紅色、灰綠色和暗灰色泥巖[19],其中砂巖以細-中粒為主,單層厚度2～15 m,孔隙度20%～35%,粒度概率累積曲線為一段式,粒度中值4.0～700.0 μm(圖2),分選中等到差[20];上部為灰色砂礫巖、砂巖與暗紅色泥巖互層,其中粒度為粉-細粒,砂巖單層厚度1～5 m,孔隙度28%～35%,粒度概率累積曲線為兩段式,粒度中值10.0～500.0 μm,如圖2所示,分選中等到較差[20]。整體粒度成下部粗、上部細的特點,與下伏古近系不整合接觸[19]。

φ為粒徑圖2 PL19-3-8ST井館陶組下、中上部粒度概率累計曲線Fig.2 Cumulative curve of grain size probability in the lower, middle and upper part of Guantao Formation in Well PL19-3-8ST

2 方法

2.1 高斯混合模型

高斯混合模型(Gaussian mixture model,GMM) 是將一個事物分解為若干個高斯概率密度函數加權和的模型[21-22],依據群集概率分配群集成員,構建多元高斯分布模型[23-24](圖3),從而找尋合理的數據分布規律和劃分方案。其數學表達式為

φ為粒徑圖3 GMM曲線及參數示意圖Fig.3 Schematic of GMM curves and parameters

(1)

式(1)中:fk(xi|θk)為多元正態密度函數;xi為第k類的密度函數;θk為相應的參數。θk由均值μk和協方差矩陣σk組成,每個向量xi都是P維的矢量。所以該分布P由G個高斯密度函數加權平均后概率密度函數表達式為

(2)

式(2)中:wk為各分量的權重;G為成分數。

Dempster等[25]提出EM算法,可在給定粒度數據集合X的情況下,尋找合適的參數θ,使似然函數最大。研究區粒度數據集合X={x1,x2,…,xn},xi從集合X中抽取,集合中元素相互獨立且符合高斯混合模型。

其中,最優聚類群集數k(即聚類簇數)是求解高斯混合模型拆分次總體數目的關鍵,k常被用作評價聚類結果的性能指標[26-28]。為選擇最佳聚類群集數k,采用均方誤差(mean squared error,MSE)[29-30]及確定系數(R2)2種誤差評價標準,這些標準可定量識別最優聚類群集數k并權衡參數似然估計擬合度好壞。因此MSE越小,模型精確度越高,誤差值越小,值為0即為完美的模型[31]。因R2取值范圍在[0,1],取值越接近1表示模型逼近效果越好[32]。

2.2 粒度數據的升采樣

為解決現有數據僅為激光粒度儀分析報告而導致數據量匱乏的問題,通過數據升采樣模擬一組符合概率分布的粒度分布仿真數據集的方法,實現概率值到精確值的轉換,其中升采樣過程中產生的空值可以通過線性插值等方式進行填充[33]。

本文研究對每組原始數據進行升采樣模擬,生成 2 000個且符合概率分布的仿真數據點,保證次總體擬合圖誤差足夠小。

3 結果分析

3.1 處理數據

以渤海地區館陶組P井1 250.4 m為例。

(1)將激光粒度儀分析報告整理為概率分布直方圖,如圖4(a)所示。

φ為粒徑; λ為面積; μ為均值; σ為方差; M為峰值; C1、C2為次總體1、2; P1/2為交點; S1為重疊面積圖4 1 250.4 m樣本擬合過程和擬合結果Fig.4 Fitting process and fitting results of 1 250.4 m samples

(2)將原始樣本升采樣,模擬符合概率分布的粒度分布仿真數據集,如圖4(b)所示。

(3)計算聚類數目從1～10時,MSE及R2數值大小并比較其大小。當k=2時, MSE最小,R2最大,即k=2為最佳聚類群集數。

(4)將k=2代入EM算法內,求解均值μi(μi為樣本均值,表示一組樣本密度概率最大時對應的粒徑值,μi<0時,為礫巖;μi>0時,μi越大粒徑越小)與方差σi(σi為樣點標準差,表示隨機變量取值的分散程度。σi值越小,峰值Mi越高;μi對應的密度概率值為峰值Mi,Mi是反映其尾數變化判斷沉積環境的有用粒度參數[7])。計算得μ1=1.27,μ2=6.51,σ1=0.65,σ2=2.42。

(5)依據μi與σi,繪制次總體擬合圖,如圖4(c)所示,運用數形結合思想,獲取交點及面積,如圖4(d)所示。

(6)結合巖心剖面,判斷分析沉積微相,探討不同組分的沉積機理和指示意義。

3.2 沉積相識別

截取該研究區館陶組PL19-3-8ST井典型深度段1 414～1 426 m、1 490～1 506 m并繪制巖性剖面和次總體擬合圖(圖5和圖6),1 414～1 426 m下部棕黃色粒質砂巖、細砂巖,發育板狀交錯層理、槽狀交錯層理,氧化及強氧化環境。中部為灰黑色泥巖、灰褐色泥質粉砂巖,還原環境。上部棕黃色含礫砂巖、細砂巖,發育板狀交錯層理,局部可見灰黑色粉砂質泥巖,主要為氧化及強氧化環境。由下至上共有三套正旋回,水動力由強至弱,為辨狀河沉積體系。1 490～1 506 m下部為灰黑色泥巖變為深褐色粉砂巖、灰褐色泥質粉砂巖,還原環境變為氧化及強氧化環境。中部為深褐色細砂巖、灰色泥質粉砂巖,發育板狀交錯層理,氧化及強氧化環境。上部為灰褐色泥質粉砂巖、深褐色粉砂巖,氧化及強氧化條件,局部見灰褐色泥巖及粉砂質泥巖。由下至上發育七套正旋回和兩套反旋回,水動力由強至弱再變強,為辨狀河沉積體系。

φ為粒徑; λ為面積; μ為均值; σ為方差; M為峰值; C1、C2為次總體1、2; P1/2為交點; S1為重疊面積圖5 PL19-3-8ST井1 412～1 426 m基于粒度次總體擬合圖的沉積相識別Fig.5 Identification of sedimentary facies based on grain size subtotal fitting diagram from 1 412 m to 1 426 m in Well PL19-3-8ST

圖6 PL19-3-8ST井1 490～1 506 m基于粒度次總體擬合圖的沉積相識別Fig.6 Identification of sedimentary facies based on grain size subtotal fitting diagram from 1 490 m to 1 506 m in Well PL19-3-8ST

4 討論

分析結果顯示,A～K的次總體擬合圖多為雙峰,偶見三峰,粒徑分布范圍-0.75φ～11.25φ,主峰出現的位置在0.75φ～7.90φ。其中C1、C2(次總體Ck:次總體曲線按照由細到粗的順序分別命名為C1、C2,…,Ck)的μ、σ、λ、M范圍分別為0.75φ～3.94φ,0.39～1.97,23%～67%,0.20～1.04;4.34φ～7.90φ,0.65～2.68,31%～77%,0.15～0.61。C1主要為中砂-粗粉砂,C2主要為粗粉砂-泥。三峰型M1=0.90,M2=0.43,M3=0.19;σ1=0.44,λ1=0.41,μ1=0.89φ,中砂;σ2=0.94,λ2=0.31,μ2=2.26φ,細砂,σ3=2.11,λ3=0.28,μ3=6.68φ,極細粉砂。溢岸砂S1(Sm/n為重疊面積,Sm/n越大,疊合度越高,即粒徑差異性越小水動力越近似;反之,水動力能量越復雜。而當μ1與μ2接近,旁瓣和Sm/n較小時,水動力能量反而近似。)整體較高12.04%～15.07%,粒徑差異性越小,疊合度較高,水動力環境越近似。辮狀河道S1集中11.64%～12.75%,混合度較穩定。心灘S1分布范圍較廣7.15%～13.50%,即水動力環境可能存在兩種,水動力條件近似的下和水動力能力復雜的上游。河漫灘S1為10.57～24.58,J的水動力條件近似,而C的水動力較J復雜。

對比不同沉積環境中的次總體擬合圖的特征,發現各個微相的次總體擬合圖各參數相互關系與巖性剖面、沉積微相均有很好的匹配性且M1差異明顯。整體而言,溢岸砂(F、H、K)中M1集中在0.8～0.9,M2均小于0.20,μ1集中在0.92φ～1.51φ,中-細砂,μ2集中在4.34φ～5.09φ,粗-中粉砂;辮狀河道(A、G)中M1分布兩個區間0.50～0.60或1.00～1.10,M2均小于0.20,μ1集中在0.75φ～1.34φ,中砂-細砂,μ2集中在2.68φ～5.89φ,細砂-中粉砂;心灘(B、D、E、I)雙峰型及三峰型,M1集中在0.60～0.80,M2均在0.20附近,μ1集中在1.42φ～3.63φ,中砂-粗粉砂,μ2集中在5.47φ～7.30φ,中粉砂-極細粉砂;三峰型M1=0.90,M2=0.43,M3=0.19,μ1=0.89φ,中砂,μ2=2.26φ,細砂,μ3=6.68φ,極細粉砂;而河漫灘(C、J)中M1集中在0.2～0.3,M2大于M1,μ1集中在1.57φ～1.97φ,中砂-細砂,μ2集中在7.65φ～7.90φ,極細粉砂-泥。

A、G雖均為辮狀河道,但G的M1大于A的M1,其原因可能為G階段發生一次小規模的洪水事件導致物源供給量短時間的增大,水動力增強,而A的M1均小于心灘(B、D、E、I)的M1,其分選性比心灘差。B、D、E、I雖均為心灘,但I的μ1均小于B、D樣本,大于A樣本且M1遠遠大于河漫灘(C、J),粗碎屑顆粒比河漫灘多,I為心灘中細粒成分。E為三峰型,C1、C2、C3各組分混雜堆積,符合心灘粒度變化范圍寬,成分復雜等地質認識,E可能為非活動性河道砂壩。C、J 為河漫灘,其M1最小,分選性差,整體粒徑較細,水動力弱,主要對應粉砂質泥巖及泥巖。F、H、K整體為細砂～粉砂,巖性剖面觀察到頂底巖性突變面,其次總體擬合圖形態幾乎保持一致,M1集中在0.8～0.9,M2均小于0.20。

從以上分析可以得出,C1為粗碎屑顆粒,C2為較細碎屑顆粒,M1及μ為劃分心灘、辮狀河道、溢岸砂及河漫灘關鍵組分,而σ越小,次總體擬合圖越窄,M1越高、μ1越小,粗碎屑顆粒值越大,則水動力越強。C2、C3為次要組分,主要營力為懸浮搬運,細碎屑顆粒填充在粗砂孔隙中,粗碎屑顆粒孔隙內水動力減弱,當流水動力不足以克服碎屑的重力時,懸浮的粉砂或黏土在流速減小到一定程度時,便開始沉降并滯留在原處,當M2越高、μ1越大,C2、C3的組分相對較高,則水動力越弱。

5 結論

提出了一種快速有效的 GMM粒度數據拆分新方法,在一定程度上解決了粒度曲線模擬過程中無法分離、提取與沉積過程相關的粒度分布次總體的問題,實現粒度資料的無監督自動聚類,準確識別正態母體的個數,精確獲得數學參數并定義其地質含義。

將GMM應用至渤海灣盆地館陶組,該地區共有2種次總體組合樣式,K=2的雙峰型及K=3的三峰型,其中M1及μ為劃分心灘、辮狀河道、溢岸砂及河漫灘關鍵組分,σi越小,Mi越高,μi越小,即次總體的圖像越窄,μi對應的碎屑含量越多,分選性越好,水動力越強。C2、C3為次要組分,M2越高、μ1越大,C2、C3的組分相對較高,水動力越弱。

次總體擬合圖的各參數相互關系反映的水動力條件與巖性剖面、沉積微相均有很好的匹配性,故該圖件在實際應用中具有良好的適用性并提高水動力條件識別效率和精度。