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

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

0 引言

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

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

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

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

1 基本原理

1.1 貝葉斯分類

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

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

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

(1)

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

∝fi(x)P(ci)，

(2)

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

1.2 自適應核密度估計

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

(3)

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

圖1展示了對于1 200組完全正態(tài)分布的樣本，選取不同帶寬得到的概率密度函數(shù)曲線圖，可以得到，對于該組數(shù)據(jù)而言，下列所選帶寬的結(jié)果擬合程度最高的是h=5時，當所選帶寬h過大，會造成核密度估計曲線過于平滑，從而失去應有的特征細節(jié)。帶寬h過小，會導致核密度估計曲線光滑性差，過于粗糙，會產(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ù)的誤差范圍也不同，因此實際操作中固定帶寬的方法可能并不適用于巖性的識別，加上在實際采樣中，受限于各種地理和人為因素的影響，無法保證樣本的稀疏程度一致，對于不同質(zhì)量的樣本，無法用同一帶寬來進行密度估計，需要對不同稀疏度的樣本分別討論帶寬，因此相比于固定帶寬法，滑動變帶寬更能滿足巖性預測的要求。令帶寬值隨著數(shù)據(jù)的密度變化自動進行適當?shù)恼{(diào)整，實現(xiàn)滑動變帶寬的方法，主要是基于積分均方誤差，通過計算每個點的最優(yōu)帶寬值，得到變帶寬函數(shù)h(x)。針對滑動窗口的實現(xiàn)過程，于傳強等[27]給出了詳細的推導過程，關(guān)于滑動變帶寬的算法流程如下：

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

(5)

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

(6)

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

(7)

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

2 模型測試

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

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

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個樣本點，分別選取250、435、870個點作為訓練樣本，其他點作為測試樣本，圖4～9分別為不同訓練樣本的物性參數(shù)交互圖以及分類結(jié)果，表2是關(guān)于不同訓練樣本錯誤率統(tǒng)計表。

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

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

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

a～f分別表示對于250個訓練樣本點的傳統(tǒng)高斯分類、固定帶寬的核密度估計、自適應帶寬的核密度估計的貝葉斯分類結(jié)果以及其對應的概率分布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點訓練樣本物性參數(shù)交互圖Fig.6 Interaction diagram of physical parameters of 435 training samples

a～f分別表示對于435個訓練樣本點的傳統(tǒng)高斯分類、固定帶寬的核密度估計、自適應帶寬的核密度估計的貝葉斯分類結(jié)果以及其對應的概率分布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點訓練樣本物性參數(shù)交互圖Fig.8 Interaction diagram of physical parameters of 869 training samples

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

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

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

a～f分別表示對于869個訓練樣本點的傳統(tǒng)高斯分類、固定帶寬的核密度估計、自適應帶寬的核密度估計的貝葉斯分類結(jié)果以及其對應的概率分布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

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

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

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

3 結(jié)論與建議

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

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