基于自適應(yīng)核密度的貝葉斯概率模型巖性識別方法研究

蔡澤園，魯寶亮，熊盛青，王萬銀

(1.長安大學(xué) 重磁方法技術(shù)研究所,陜西 西安 710054；2.長安大學(xué) 地質(zhì)工程與測繪學(xué)院，陜西 西安 710054；3.長安大學(xué) 西部礦產(chǎn)資源地質(zhì)工程教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054；4.中國自然資源航空物探遙感中心, 北京 100083)

0 引言

巖性識別是地質(zhì)研究過程中非常重要的基礎(chǔ)工作，尤其是在近地表以及深部無法直接采樣區(qū)的地質(zhì)研究中，準(zhǔn)確地刻畫深部巖石類型及其結(jié)構(gòu)關(guān)系，可以為能源礦產(chǎn)勘探、深部結(jié)構(gòu)與構(gòu)造等研究提供重要的地質(zhì)信息。因此采用什么數(shù)據(jù)、什么方法來進(jìn)行巖性識別是一項(xiàng)極具價(jià)值的研究工作。傳統(tǒng)(地表)地質(zhì)填圖已遠(yuǎn)不能滿足深部探測的需要。隨著多源地學(xué)數(shù)據(jù)(地質(zhì)、地球物理、地球化學(xué)、遙感以及鉆井?dāng)?shù)據(jù)等)的獲取、地球物理三維反演技術(shù)和三維地質(zhì)建模的發(fā)展，深部巖性識別已成為現(xiàn)實(shí)。但目前對多源數(shù)據(jù)的利用并不充分，尚未有效地利用多種物性參數(shù)進(jìn)行巖性識別。因此，如何利用多源不同類型、屬性地學(xué)數(shù)據(jù)所反映的巖性信息，從高維數(shù)據(jù)空間準(zhǔn)確地進(jìn)行巖性識別是亟待解決的難題，對深部探測和地學(xué)數(shù)據(jù)融合也具有重要的科學(xué)意義。

目前巖性識別主要有直接手段(巖芯、手標(biāo)本以及薄片鑒定等)和間接手段(重磁、地震、電磁、測井、遙感、地球化學(xué)等)。從深部探測來講，鉆井巖芯可能是唯一的巖石標(biāo)本，但只是點(diǎn)數(shù)據(jù)，而遙感(面數(shù)據(jù))只能探測地表情況，且受地表覆蓋影響大；因此地球物理方法在地下深部巖性識別研究中將發(fā)揮重要作用。而物性(密度、磁化率、電導(dǎo)率、縱橫波速度等)是巖性和地球物理場之間的紐帶，因此，通過地球物理反演物性，再聯(lián)合其他地質(zhì)數(shù)據(jù)進(jìn)行巖性的識別具有可行性。地震反演是進(jìn)行巖性識別的有效方法，可以利用地震不同的彈性參數(shù)，如利用縱橫波速度、密度[1-2]、波阻抗、振幅、頻率、相位[3-5]對目標(biāo)巖體(流體)進(jìn)行識別，但是當(dāng)巖相間的地震響應(yīng)差別不明顯時(shí)，依靠波形的分類結(jié)果無法準(zhǔn)確刻畫巖性且地震方法成本高。在測井識別巖性方面，可以利用交會圖版[6-9]或者基于統(tǒng)計(jì)、聚類、支持向量機(jī)以及神經(jīng)網(wǎng)絡(luò)等技術(shù)[10-14]進(jìn)行巖性識別，測井方法在精度、算法以及技術(shù)上都有明顯的優(yōu)勢，但是只能識別鉆孔附近小范圍內(nèi)的巖性，難以進(jìn)行大面積或是沒有鉆孔地區(qū)的巖性識別。利用重磁技術(shù)識別巖性主要依靠密度與磁化率比值，例如采用交會圖版結(jié)合邏輯運(yùn)算識別巖性[15-16]，重磁數(shù)據(jù)覆蓋面積廣，采樣密度高，容易獲得大規(guī)模的巖性識別，但是在無約束情況下，垂向分辨較差，多解性強(qiáng)。

單一的地球物理方法往往很難獲得理想識別結(jié)果，因此聯(lián)合多個(gè)地球物理方法的巖性識別成為主流方式。在多源數(shù)據(jù)識別巖性方面，對于存在巖性與物性對應(yīng)關(guān)系的模糊區(qū)域，以地震反射特征為約束，利用重磁電資料識別了具有密度、磁性和電阻率等特征差異的火成巖巖性[17]。近年來隨著機(jī)器學(xué)習(xí)(machine learning)的興起，該算法已被廣泛應(yīng)用于巖性識別。如利用多種地球物理數(shù)據(jù)基于模糊聚類分析[18]將巖石進(jìn)行分類識別。利用地球化學(xué)數(shù)據(jù)和地球物理數(shù)據(jù)對巖性進(jìn)行多元回歸分析[19-20]或者多準(zhǔn)則決策方法識別。利用無監(jiān)督模糊分區(qū)聚類對航空伽馬射線數(shù)據(jù)與陸地衛(wèi)星波段數(shù)據(jù)聯(lián)合進(jìn)行巖性識別[21]。應(yīng)用受限玻爾茲曼機(jī)(restricted boltzmann machine)和隨機(jī)森林模型到區(qū)域尺度多參數(shù)地球科學(xué)數(shù)據(jù)集，從而預(yù)測斑巖銅金礦床的遠(yuǎn)景區(qū)域。還有采用隨機(jī)森林法(random forests)和自組織映射技術(shù)(SOM,self-organising maps)來識別連續(xù)的火山單元子類，由于算法只針對部分地球科學(xué)或地質(zhì)參數(shù)進(jìn)行巖性識別，因此無法針對不同類型的參數(shù)進(jìn)行統(tǒng)一處理。雖然各類算法都有改進(jìn)，但是巖性識別結(jié)果唯一，對模糊區(qū)間的多種可能性無法準(zhǔn)確表述。

理論上講，通過不同地質(zhì)、地球物理技術(shù)可以獲得地下物性結(jié)構(gòu)，如密度結(jié)構(gòu)、速度結(jié)構(gòu)、電阻率結(jié)構(gòu)、磁性結(jié)構(gòu)等，那么如何將這些物性結(jié)構(gòu)轉(zhuǎn)換為巖性是值得研究的問題，實(shí)質(zhì)上這是一個(gè)模式識別問題。通常來講，巖性與物性的對應(yīng)關(guān)系并不總是明確的，在交會圖上存在較大的重疊區(qū)域，從而使得基于規(guī)則的分類方法難以解決該問題。在前人的研究基礎(chǔ)之上，考慮到貝葉斯方法是非規(guī)則分類，該方法依據(jù)類的概率、概率密度，并按照某種規(guī)則使得分類結(jié)果從統(tǒng)計(jì)上達(dá)到最佳?；诖?，筆者提出了基于自適應(yīng)核密度估計(jì)的貝葉斯概率巖性識別方法，完成了從物性到巖性的轉(zhuǎn)換。該方法具有較強(qiáng)的泛化能力，預(yù)測的巖性分類結(jié)果帶有概率參數(shù)，可以存在模糊區(qū)間，提供多種巖性分類結(jié)果。該方法具有較強(qiáng)可擴(kuò)展性，可以有任意數(shù)量類型的輸入?yún)?shù)(允許存在缺省參數(shù))以及任意數(shù)量的巖性分類輸出。通過實(shí)驗(yàn)對比了傳統(tǒng)的高斯密度、固定帶寬核密度以及自適應(yīng)帶寬核密度對巖石物性數(shù)據(jù)判別的效果，說明了該方法具有良好巖性識別效果。

1 基本原理

1.1 貝葉斯分類

貝葉斯算法是基于統(tǒng)計(jì)學(xué)的基本算法之一，假設(shè)各個(gè)條件之間相互獨(dú)立，可以得到樸素貝葉斯算法。如果我們將巖石的類型表示為c事件，將巖石屬性表示為x事件，巖石類型c和對應(yīng)屬性x是發(fā)生在同一空間的兩個(gè)事件，假設(shè)某研究區(qū)的完整巖石類型c是由兩種巖石類型c1、c2、…、cn構(gòu)成，c1、c2、…、cn中一個(gè)巖石種類出現(xiàn)必然伴隨著某一屬性x的發(fā)生，即若x發(fā)生，則c必然有一個(gè)會發(fā)生，根據(jù)概率可以得到樣本集的已知各類別ci的先驗(yàn)概率以及各類條件概率P(x/ci)。對于未知樣本，貝葉斯公式可以計(jì)算出待測樣本分屬各類的概率，稱為后驗(yàn)概率。

使用貝葉斯定理來預(yù)測后驗(yàn)概率最大的類，主要是估計(jì)每一類的概率密度函數(shù)，通過多元正態(tài)分布來建模。樸素貝葉斯分類器基于條件獨(dú)立性假設(shè)，是概率分類器中最簡單的分類器，在很多情況下具有相當(dāng)高的分類準(zhǔn)確率，因此以高效率和良好的泛化能力而著稱。對于某個(gè)測區(qū)，已知該地區(qū)的物性分布特征(如密度、磁化率、電阻率等)以及部分的巖石樣本，假設(shè)屬性之間相互獨(dú)立，可以使用貝葉斯分類器來對整個(gè)測區(qū)的巖石類型進(jìn)行預(yù)測。概率分類器的優(yōu)點(diǎn)在于在得到分類結(jié)果的同時(shí)，會對每一種類別進(jìn)行概率計(jì)算，對于巖性的識別而言，可以通過概率或相對概率來人為判定分類結(jié)果的可信度而不僅僅依靠算法本身的置信度來決定，小概率區(qū)間會提供一個(gè)模糊帶，模糊帶的類別區(qū)分結(jié)果可能不唯一。

貝葉斯分類器通過對每個(gè)未知樣點(diǎn)x和每個(gè)類ci來估計(jì)其后驗(yàn)概率P(ci/x)，即計(jì)算未知樣本x屬于ci類巖石的概率，從而選擇最大概率的類作為未知樣本x最終的預(yù)測類型。利用貝葉斯定理，后驗(yàn)概率P(ci/x)可以表示為：

(1)

對于給定的一點(diǎn)x，P(x)是確定的。用fi(x)表示未知樣本x屬于ci類的概率密度，似然用概率密度可以表示為P(x/ci)=2afi(x)，其中a表示鄰近x的一個(gè)極小區(qū)間，因此可以得到后驗(yàn)概率的計(jì)算公式：

∝fi(x)P(ci)，

(2)

由式(2)可知，概率密度函數(shù)是影響分類結(jié)果的一個(gè)重要因素，由于大部分巖石物性參數(shù)是基本遵循正態(tài)分布的規(guī)律，這里考慮了傳統(tǒng)的高斯公式作為概率密度函數(shù)，然而樣本本身并不完全遵循正態(tài)分布，為了極大地避免由于選定高斯函數(shù)的影響，同時(shí)考慮使用核密度估計(jì)的方法，核密度估計(jì)的優(yōu)點(diǎn)在于核函數(shù)的選取對于最終分類結(jié)果影響不大，更適用于類正態(tài)數(shù)據(jù)樣本，這會在之后的模型測試中加以驗(yàn)證。

1.2 自適應(yīng)核密度估計(jì)

核密度估計(jì)[22-23]是統(tǒng)計(jì)學(xué)中用來估計(jì)隨機(jī)變量的概率密度函數(shù)的方法，屬于非參數(shù)檢驗(yàn)方法的一種。核密度估計(jì)的基本表達(dá)式為：

(3)

核密度帶寬的選擇是得到最佳估計(jì)結(jié)果的關(guān)鍵。固定帶寬是比較常見的方法[25-26]。不同的帶寬選取，會對概率密度函數(shù)產(chǎn)生不一樣的影響，由于概率密度函數(shù)存在一個(gè)必然性，即概率密度之和為1，因此，不同的帶寬在影響曲線光滑程度的同時(shí)，也會對峰值產(chǎn)生影響，曲線越平滑，峰值就會越低，曲線越抖，峰值就會越高。帶寬選擇的基本思想是基于最小平方差，根據(jù)積分均方誤差最小，求出最優(yōu)帶寬。

圖1展示了對于1 200組完全正態(tài)分布的樣本，選取不同帶寬得到的概率密度函數(shù)曲線圖，可以得到，對于該組數(shù)據(jù)而言，下列所選帶寬的結(jié)果擬合程度最高的是h=5時(shí)，當(dāng)所選帶寬h過大，會造成核密度估計(jì)曲線過于平滑，從而失去應(yīng)有的特征細(xì)節(jié)。帶寬h過小，會導(dǎo)致核密度估計(jì)曲線光滑性差，過于粗糙，會產(chǎn)生過擬合的問題。對于高斯核函數(shù)，最優(yōu)帶寬為：

圖1 不同帶寬下概率密度函數(shù)曲線Fig.1 Probability density function curve under different bandwidth

hopt=1.059σn-0.2。

(4)

巖石的物性參數(shù)往往存在較大的差異，而每種參數(shù)的誤差范圍也不同，因此實(shí)際操作中固定帶寬的方法可能并不適用于巖性的識別，加上在實(shí)際采樣中，受限于各種地理和人為因素的影響，無法保證樣本的稀疏程度一致，對于不同質(zhì)量的樣本，無法用同一帶寬來進(jìn)行密度估計(jì)，需要對不同稀疏度的樣本分別討論帶寬，因此相比于固定帶寬法，滑動變帶寬更能滿足巖性預(yù)測的要求。令帶寬值隨著數(shù)據(jù)的密度變化自動進(jìn)行適當(dāng)?shù)恼{(diào)整，實(shí)現(xiàn)滑動變帶寬的方法，主要是基于積分均方誤差，通過計(jì)算每個(gè)點(diǎn)的最優(yōu)帶寬值，得到變帶寬函數(shù)h(x)。針對滑動窗口的實(shí)現(xiàn)過程，于傳強(qiáng)等[27]給出了詳細(xì)的推導(dǎo)過程，關(guān)于滑動變帶寬的算法流程如下：

第一，根據(jù)固定帶寬的經(jīng)驗(yàn)公式，選擇樣本組的最優(yōu)固定帶寬hopt以及對應(yīng)的核密度估計(jì)函數(shù)fopt(x)，

(5)

第二， 根據(jù)給定的估計(jì)點(diǎn)，求取優(yōu)化后的帶寬

(6)

第三，優(yōu)化后的核密度估計(jì)函數(shù)記為

(7)

上述計(jì)算的帶寬在數(shù)據(jù)質(zhì)量相對較好的情況下往往會出現(xiàn)過擬合現(xiàn)象，反而結(jié)果不盡如人意，因此在實(shí)際的分類過程中，當(dāng)計(jì)算得到固定帶寬以及優(yōu)化后的帶寬后，需要對帶寬的差值進(jìn)行判定，對于二者相差不大的情況(差值需要根據(jù)樣本的數(shù)量和質(zhì)量人為給定)，采用優(yōu)化前的帶寬。將優(yōu)化后的核密度估計(jì)函數(shù)作為貝葉斯模型的概率密度函數(shù)，針對某一未知樣本，可以得到該樣本歸屬于每一類的概率密度，通過比較得到最大概率類別作為預(yù)測的類。

2 模型測試

本次實(shí)驗(yàn)選取密度、磁化率和電阻率作為分類依據(jù)，測區(qū)內(nèi)巖石類型為板巖、片麻巖和花崗巖，模型具體參數(shù)見表1，從表中可以看出，整個(gè)模型的密度變化較小，但是相比于其他兩種巖石的密度，花崗巖的密度范圍變化更小，主要集中在 2 600 kg/m3左右，可以作為劃分花崗巖的物性依據(jù)。圖2、3分別是模型區(qū)內(nèi)的巖石物性分布情況以及該模型下的巖石分布示意圖，為了對比分類器模型的準(zhǔn)確率、穩(wěn)定性以及計(jì)算復(fù)雜度，從所有的樣本中選取不同位置、不同數(shù)據(jù)量的樣本作為訓(xùn)練集和測試集進(jìn)行分類，判定分類結(jié)果。

表1 模型參數(shù)統(tǒng)計(jì)

a—模型密度；b—模型磁化率；c—模型電阻率a— density of the model；b—magnetic susceptibility of the model；c—resistivity of the model

圖3 巖石分布Fig.3 rock distribution of the model

該模型共有8 690個(gè)樣本點(diǎn)，分別選取250、435、870個(gè)點(diǎn)作為訓(xùn)練樣本，其他點(diǎn)作為測試樣本，圖4～9分別為不同訓(xùn)練樣本的物性參數(shù)交互圖以及分類結(jié)果，表2是關(guān)于不同訓(xùn)練樣本錯誤率統(tǒng)計(jì)表。

表2 模型錯誤率統(tǒng)計(jì)

圖4 250點(diǎn)訓(xùn)練樣本物性參數(shù)交互圖Fig.4 Interaction diagram of physical parameters of 250 training samples

對照紅色散點(diǎn)和灰色的巖性邊界可以看出，識別錯誤的樣本集中分布在邊界的低概率區(qū)，而巖性邊界可以看作整個(gè)區(qū)域的模糊區(qū)間，因此概率密度圖在低概率區(qū)刻畫了整個(gè)區(qū)域的模糊區(qū)，對比圖5、圖7、圖9的概率圖，針對同一訓(xùn)練樣本，核密度估計(jì)算法對于模糊區(qū)的刻畫整體優(yōu)于傳統(tǒng)的高斯算法估計(jì)模型;對于同一算法下的不同訓(xùn)練樣本，傳統(tǒng)高斯算法的低密度帶并沒有表現(xiàn)出由于訓(xùn)練樣本增加其結(jié)果得到改善，而隨著訓(xùn)練樣本的增加，自適應(yīng)帶寬核密度估計(jì)算法的概率圖上對于分類結(jié)果正確區(qū)域的低概率帶有了明顯改善。

a～f分別表示對于250個(gè)訓(xùn)練樣本點(diǎn)的傳統(tǒng)高斯分類、固定帶寬的核密度估計(jì)、自適應(yīng)帶寬的核密度估計(jì)的貝葉斯分類結(jié)果以及其對應(yīng)的概率分布Figures a～f respectively represent the Bayesian classification results， corresponding probability distribution map of the traditional Gaussian classification, fixed bandwidth kernel density estimation， adaptive bandwidth kernel density estimation for 250 training sample points

圖6 435點(diǎn)訓(xùn)練樣本物性參數(shù)交互圖Fig.6 Interaction diagram of physical parameters of 435 training samples

a～f分別表示對于435個(gè)訓(xùn)練樣本點(diǎn)的傳統(tǒng)高斯分類、固定帶寬的核密度估計(jì)、自適應(yīng)帶寬的核密度估計(jì)的貝葉斯分類結(jié)果以及其對應(yīng)的概率分布Figures a～f respectively represent the Bayesian classification results， corresponding probability distribution map of the traditional Gaussian classification, fixed bandwidth kernel density estimation， adaptive bandwidth kernel density estimation for 435 training sample points

圖8 869點(diǎn)訓(xùn)練樣本物性參數(shù)交互圖Fig.8 Interaction diagram of physical parameters of 869 training samples

從表2的錯誤率上而言，在同樣的訓(xùn)練樣本下，基于自適應(yīng)核密度估計(jì)的貝葉斯概率模型明顯優(yōu)于其他兩類模型。由此可知，貝葉斯概率密度模型對于深部巖性識別的方法是有效可行的，而基于自適應(yīng)核密度估計(jì)的貝葉斯概率模型相比于其他概率密度模型而言效果相對較好。

通過測試不同訓(xùn)練樣本對于分類結(jié)果的影響，最終可以得到針對該模型，訓(xùn)練樣本每類都少于40個(gè)時(shí)便無法進(jìn)行合理的巖性識別，當(dāng)然并非訓(xùn)練樣本的數(shù)量越多越好，樣本質(zhì)量對于識別結(jié)果也有著決定性作用。

通過對比3種概率密度方法得到的分類結(jié)果可知，基于自適應(yīng)帶寬的貝葉斯概率模型下的分類錯誤率是相對較低的且在模糊區(qū)有一個(gè)明顯的低概率帶，因此，模型二將針對模型一中435點(diǎn)訓(xùn)練樣本在自適應(yīng)帶寬下的測試樣本進(jìn)行概率差下的統(tǒng)計(jì)分析。

a～f分別表示對于869個(gè)訓(xùn)練樣本點(diǎn)的傳統(tǒng)高斯分類、固定帶寬的核密度估計(jì)、自適應(yīng)帶寬的核密度估計(jì)的貝葉斯分類結(jié)果以及其對應(yīng)的概率分布Figures a～frespectively represent the Bayesian classification results， corresponding probability distribution map of the traditional Gaussian classification, fixed bandwidth kernel density estimation， adaptive bandwidth kernel density estimation for 869 training sample points

若兩種類型巖石的預(yù)測概率差在擬定的差值之內(nèi)，就認(rèn)為該未知樣本屬于這兩類巖性的概率相同，對于3種巖性也同樣適用，對于該模型而言,若3種巖性的概率差都在擬定的差值內(nèi)，則認(rèn)為無法對該樣本進(jìn)行巖性識別，選擇更大的概率差，意味著模型的容錯率降低，同時(shí)，識別類型唯一的樣本可信度也隨之變大，針對不同的概率差，預(yù)測結(jié)果不同。

圖10為概率差在10%、20%、30%、40%、50%、60%概率差下的巖性預(yù)測結(jié)果，白色區(qū)域表示3種巖性被認(rèn)為是等概率的情況，即無法識別的區(qū)域，同時(shí)，淺色區(qū)域代表可能為兩種巖性。在無法進(jìn)行準(zhǔn)確的巖性識別的區(qū)域進(jìn)行人工識別或者二次識別，避免了機(jī)器學(xué)習(xí)在巖性識別過程中由于算法本身的局限性，導(dǎo)致預(yù)測結(jié)果唯一且無法進(jìn)行人工推斷可信度的問題，在機(jī)器學(xué)習(xí)和人工識別中找尋一個(gè)平衡，發(fā)揮二者優(yōu)勢，在提高巖性識別效率和準(zhǔn)確度的同時(shí)使得預(yù)測結(jié)果明朗可操作。

a～f依次為概率差在10%、20%、30%、40%、50%、60%時(shí)的分類結(jié)果a～f are the classification results when the probability difference is 10%, 20%, 30%, 40%, 50% and 60% respectively

3 結(jié)論與建議

筆者將基于自適應(yīng)核密度估計(jì)的貝葉斯概率模型應(yīng)用到巖性識別中，該方法的優(yōu)勢在于，對于有一定重合區(qū)域的物性參數(shù)，可以有效地進(jìn)行未知區(qū)域的巖性識別，提高計(jì)算的精度和效率，將多種物性參數(shù)進(jìn)行合理的處理解釋，提高數(shù)據(jù)利用率，通過各個(gè)數(shù)據(jù)之間的相互作用，最終獲得巖性識別結(jié)果。

從本文的模型可以看出，基于自適應(yīng)核密度估計(jì)的貝葉斯概率模型能夠較好地刻畫未知地區(qū)巖性的類別和輪廓，對于模糊區(qū)也可以較好地進(jìn)行判定。事實(shí)上，該方法并不局限于以上幾種參數(shù),對于其他的物性參數(shù)也可以進(jìn)行同樣的處理，識別結(jié)果的精度也依靠于更多類型的物性參數(shù)和更好質(zhì)量的訓(xùn)練樣本。