冀全偉王文磊劉治博祝茂強袁長江
1.自然資源部古地磁與古構造重建重點實驗室,北京 100081;2.中國地質科學院地質力學研究所,北京 100081;3.中國地質大學(北京),北京 100083;4.中國地質科學院礦產資源研究所,北京 100037
信息化時代,社會經濟發展與生態環境治理對地質調查工作提出了新的要求,地質調查工作面臨新的機遇與挑戰。例如,在特殊地質地貌區開展區域地質調查工作將有助于特殊地質景觀區基礎地質問題的研究,服務于多門類自然資源與生態環境問題的解決(胡健民和陳虹,2019)。隨著地質調查工作的持續開展,基礎地質研究程度不斷提高,成果數據資料保持快速積累與更新。如何系統整合已有地質、地球化學、地球物理、遙感等多元、多尺度地質調查數據資料,發展能夠提高工作質量與效率的方法,深度挖掘有用信息,進而優化提升基礎地質、礦產地質、水文地質、災害地質等調查評價技術(楊星辰等,2020;張鑫剛等,2020),被認為是地質調查工作手段升級,提高社會經濟服務能力的突破口之一。亟需學習吸收并引進數學、信息學等學科先進的數據與信息挖掘技術,創新發展地質調查評價思路與方法。
地質填圖作為區域地質調查工作最基本的核心工作內容之一,其效率和精度將直接影響后續研究工作的開展。傳統地質填圖工作主要包括前期資料收集整理、工作方案編制、野外實地勘查、樣品測試分析及數據處理、成圖及報告編寫等階段。其中,前期資料收集整理工作多停留在基本資料了解階段,基礎資料及數據的應用程度不高;而野外工作依靠地質工作者的主觀判斷來確定填圖單元,受限于填圖技術人員的業務水平不同,填圖質量受到一定影響。因此,為保證填圖成果質量,野外實地勘查工作需投入較高的人力、財力和物力成本來完成大量路線調查及剖面測量等實物工作量。此外,在偏遠山區、無人區、高原地區開展野外工作還存在一定風險性。
隨著機器學習方法的快速發展,基于機器學習的巖性填圖方法的提出,取得了較好的研究成果與進展。相較傳統地質填圖技術,機器學習方法中的分類模型或組合算法在巖性分類識別方面具有高效、智能化的特點,可作為具有巨大潛在優勢的輔助手段來提高傳統地質填圖技術方法體系的工作效率與能力。已有基于機器學習方法的巖性填圖研究(吳俊等,2016;陳松等,2017),通過系統整合多源遙感、地震、物探、化探、航磁等數據,建立巖性分類的基礎數據集,利用度量學習、支持向量機(SVM)、人工神經網絡(ANN)、隨機森林(RF)等機器學習分類算法,開展了巖性識別、巖性單元填圖等相關分類問題的研究 (Cracknell and Reading, 2014;Harris and Grunsky, 2015;鄭陽,2017;Othman and Gloaguen, 2017;Kuhn et al., 2018;張艷等,2019;段友祥等,2020;朱明永等,2020;Wang et al., 2020a, 2020b;Wu et al., 2021)。已有研究表明,這一巖性填圖思路在特定地質條件下具有特殊優勢(嚴昊偉等,2017)。
文章主要通過野外基礎地質調查和機器學習分類算法的有機融合,在填圖空白地區或工作程度較低地區開展基于勘查數據分析預測的巖性單元填圖方法探索性研究。選取西藏多龍礦集區開展模型試驗主要是考慮到兩方面原因。首先,多龍礦集區是中國重要成礦區帶班公湖-怒江成礦帶內已發現最大的斑巖型Cu-Au礦產地,具有巨大資源潛力。區內近年來已完成了1∶5萬區域與礦產地質調查工作,對巖性單元劃分具有較為清晰的認識,有利于預測結果的驗證與應用效果評價。其次,多龍礦集區積累了大量基礎圖件和勘查數據資料,可供研究通過選取不同基礎預測數據組合,構建不同工作基礎條件下的模型方法試驗。同時,文中提出的數據填圖方法需要開展多批次小范圍野外填圖支撐巖性單元預測的迭代算法。在模型試驗過程中,已有地質圖件能夠代替野外填圖直接為預測模型提供原始數據和現有知識補充。換而言之,通過從已有地質圖中提取迭代算法所需的小范圍巖性單元分布來實現數據集與知識庫的更新,為模型試驗節省了實際野外填圖的時間成本。因此,研究以多龍礦集區為模型試驗區,選擇1∶5萬勘查地球化學數據為基礎預測數據,以1∶5萬區域地質圖為參考,進行基于梯度提升決策樹算法的巖性預測填圖模型試驗。
多龍礦集區位于西藏阿里地區改則縣境內,所處的大地構造位置為班公湖-怒江成礦帶西段,班公湖-怒江縫合帶北側、羌塘-三江復合板片南緣(郭娜等,2018;李興奎等,2018;任紀瞬等,2019)。地層分區屬于羌南-保山地層區多瑪地層分區,區內地層(圖1)以中生界為主,主要有中侏羅統曲色組 (J2q)和色哇組 (J2s)濁積巖建造、下白堊統美日切錯組(K1m)火山碎屑巖建造以及新生界新近系康托組(N1k)陸源碎屑巖建造和第四系殘坡積物(Q4)(江少卿等,2014;陳紅旗等,2015)。其中,J2q組巖性為粉砂質板巖夾變長石石英砂巖(李云強等,2020),J2s組的巖石主要由砂巖、砂礫巖和變長石石英砂巖等組成(符家駿等,2014),同時兩組地層也是含礦巖體的主要圍巖(王勤等,2018)。K1m組的巖石主要為安山巖、英安巖、玄武巖、火山角礫巖和碎屑巖等。N1k組以礫巖、含礫砂巖、紅色泥巖為主要巖性(韋少港等,2017)。多龍礦集區巖漿活動極為發育,總體上以噴發、噴溢和淺成、超淺成侵入為主,具多期活動特征,形成時代為燕山中—晚期(江少卿等,2014;李紅梅,2017)。噴出巖主要由玄武巖、安山巖和流紋巖組成, (孫嘉等,2019)。侵入巖主要為基性、中酸性侵入巖,基性巖主要為輝長巖和輝綠巖,中酸性淺成巖主要為閃長巖、英安巖、花崗閃長斑巖,侵入時代以早白堊為主(陳紅旗等,2015;王勤等,2018)。區內接觸變質巖變質程度不高,巖體周邊廣泛發育熱液蝕變及少量石英脈(王繼斌,2018)。

圖1 多龍礦集區巖性分布圖Fig.1 Spatial distribution of the lithologic units in the Duolong mineral district, Tibet, China
自20世紀70年代以來,先后有多家地勘和研究單位對多龍礦集區開展了1∶100萬、1∶25萬和1∶5萬圖幅的區域地質調查工作。該區的區域物探、化探、遙感、礦床勘查工作以及相關巖石地球化學(韋少港等,2019)、年代學(王勤等,2015)、控礦構造識別(劉治博等,2017)、遙感異常信息提取(代晶晶等,2013;別小娟等,2013)、蝕變礦物學(趙子歐等,2020)等方面研究取得了較好的成果進展。通過近些年多方面研究,對多龍礦集區的地質背景、成礦規律、礦床模型等有了新的認識(楊歡歡等,2019;王勤等,2019;石洪召等,2019;孫嘉等,2020),目前正根據已有資料開展進一步綜合研究。
基于機器學習方法的巖性填圖對研究區的基礎地質數據積累與研究程度具有較高要求,大多針對特定的數據資料類型且依賴高質量數據集,在空白區或數據資料不充分地區開展工作,將會面臨缺乏基礎地質支撐的困難。文中通過野外地質調查與機器學習方法的有機融合,提出了一種基于梯度提升決策樹 (Gradient boosting decision tree, GBDT)算法的巖性單元填圖方法(圖2):①選擇研究區內小范圍已填圖區作為假想野外填圖區,建立原始數據集并初步構建巖性單元與預測數據(遙感、化探、物探)對應關系;②利用機器學習方法對預測數據進行多分類任務,進而開展目標填圖區預測填圖工作;③通過概率選區選定概率較低目標區,開展進一步的小范圍野外地質調查假想填圖,對原始數據和現有知識進行補充;④迭代循環以上流程,直至預測填圖達到要求。

圖2 基于機器學習的巖性填圖思路Fig.2 Flowchart of machine learning-based lithologic mapping
數據預處理時,若將研究區整體定義為單一柵格作為目標選區基本單元,代表性較弱,同時易受模型分類過程中分類準確率的影響。因此,需要通過對研究區進行網格化劃分(圖1),將基本單元由單一柵格分解為w×h個網格單元,并在此基礎上進行概率均值的統計,以此作為迭代填圖目標選區的評判基礎。文中將多龍礦集區內填圖范圍劃分成90個網格單元,網格單元面積為3.5 km×3.5 km。
針對研究區進行網格化處理后,通過隨機選區的采樣策略完成初始數據集的創建。從研究區劃分網格單元中隨機選取若干單元作為目標采樣區。通過野外地質調查在選區內開展地質填圖,獲取區內巖性單元分布情況。模型試驗將通過從已有地質圖中直接提取選區內的巖性單元分布來代替野外實際填圖工作。
通過距離反比權重法(IDW)對試驗區3200個地球化學數據點進行空間插值,得到的柵格數據作為模型試驗的預測數據。將初始選區的巖性填圖結果與對應的地球化學數據進行標簽化整合,完成初始數據集的創建。最后,通過模型訓練建立巖性算法分類模型,根據模型評價標準實施迭代填圖,預測全區巖性分布結果,進而探索基于GBDT算法的巖性填圖方法。
梯度提升決策樹 (Gradient boosting decision tree, GBDT)算法(Friedman, 2001)是一種采用集成學習思想的迭代決策樹算法。所謂集成學習,即通過對多個學習器(如決策樹)的組合得到比單一學習器性能更好的算法模型訓練策略。一般情況下,GBDT以決策樹(Quinlan, 1986)為基礎分類器,并利用損失函數的負梯度作為提升樹殘差的近似值進行算法實現。其中,提升樹fM(x)可表示為:

其中,Tm(x)為弱學習器,即決策樹;γm為每個弱學習器最優擬合的權重;M為樹的個數,即迭代次數。
模型的訓練過程是損失函數L的最小化過程。假設訓練樣本數據量為N,第i條數據的變量與真值分別為xi和yi,則參數調優的目標函數為:

其中,表示訓練完成的預測模型;L為訓練過程中的損失函數;argmin則表示最小化損失函數L時f的取值;其他變量同公式(1)。
歸一化指數函數(Softmax)是邏輯函數在多分類任務上的一種推廣,其目的是將多分類結果以概率的形式展現出來。若以DT表示樣本訓練集,則DT={(xi,yi),i=1,…,nT}。 其中,xi是模型輸入的數據,如用來預測巖性單元的遙感、地球化學等數據;yi是對應地質目標名稱,如巖性單元標簽。假設訓練集巖性單元種類數為K,則一般情況下nT>K。在分類問題上,GBDT的作用是計算xi與yi之間的映射函數f:R15→RK。 對于輸入的x,輸出P維特征向量ν,并代入Softmax函數計算分類概率值:

其中,pk表示屬于第k類巖性的預測概率值。根據Softmax計算公式可知,對于任一數據x,各巖性預測概率之和必為1。
文中采用GBDT作為核心算法對區內地球化學數據與巖性單元的對應關系開展信息挖掘與擬合工作。針對小樣本數據集,特別是當前基礎預測數據小于104數量級的情況下,GBDT算法在訓練的過程中可能會出現過擬合問題。目標函數在機器學習過程中將會過度依賴訓練樣本集,將所有樣本(包括噪聲)都擬合到函數當中,從而只在訓練集中表現優異,對于未知樣本則無法正確預測。因此,為客觀判斷訓練參數對訓練集以外數據的符合程度,論文采用交叉驗證的思想對模型整體分類能力進行評估。將樣本數據集隨機分為F個不相交子集,從F個子集中逐次選取一個子集定義為測試集,其余F-1個子集定義為訓練集,基于訓練集進行訓練得到GBDT模型。利用測試集對模型進行分類器性能評價,將F次測試結果的均值定義為F折交叉驗證下模型性能指標,并以此來評估模型精度。此外,需要在交叉驗證基礎上進行多次參數調優,得到更為合理的模型參數,以保證訓練得到的GBDT模型具備較強的分類能力。
根據每次迭代過程中對模型進行多次訓練的結果(圖3)可知,經過300次訓練后模型表現趨于平穩,損失值基本穩定在0.2。這說明即使對于較為復雜的多分類問題,該模型仍具有較強的有效性和穩定性。

圖3 模型損失函數統計圖Fig.3 Statistical diagram of the loss function
從概率角度選定網格單元,將其作為目標填圖區進行針對性的迭代填圖,并逐步更新預測分類數據集是此次研究思路核心之一。迭代填圖這一思路作為整套方法流程中最主要的數據與知識補充過程,其準確性高低將對最終出圖結果造成直接影響。與傳統巖性填圖相結合,通過專家野外填圖的方式完成概率選區范圍內的信息采集工作,在保證結果精度的前提下減少傳統巖性填圖的野外實際工作量。在具體實施過程中,根據研究區預測概率分布結果(圖4),以網格為基本單元進行概率均值計算。按概率高低對全部單元進行排序,選取其中概率最低的若干網格單元(圖4中黑框位置)作為目標區域,開展野外局部實地填圖。將填圖區巖性分類結果與對應的地球化學數據進行整合,并更新至樣本數據庫。

圖4 概率分布選區示意圖Fig.4 Schematic diagram of probability distribution-based area selection
模型評價主要包括適用性和實用性評價兩個方面。模型適用性評價主要是從算法角度評價GBDT模型對地質問題的適用性。針對從區內網格中選取的野外填圖區,根據野外填圖獲得巖性分布,按比例劃分出預測評價區。訓練模型應用于預測區獲得相應的巖性分類結果。以地質圖為真值統計分類結果的準確率、宏平均F1分數等模型評價指標,并根據各類指標情況進行模型修正。
模型實用性評價主要是從預測概率角度評價預測結果對預期分類結果的滿意程度。預測概率值是將模型輸出值與各類巖性單元特征向量之間的殘差通過Softmax函數進行歸一化計算獲得。概率值高低代表當前地球化學數據分類結果與各巖性單元類型的相近程度。假設已知專家填圖區巖性單元種類集合為S,則概率分布高值區通常代表當前第i區域分類ki∈S,低值區表示當前分類范圍較大可能存在實際巖性單元種類ki?S的情況。基于以上原則,將模型預測概率與預期結果進行對比。若滿足,則將模型應用于全區地球化學數據并預測全區巖性單元分類,否則,進行迭代填圖,直至滿足預期。
文中采用準確率(Accuracy,簡記Ac)、宏平均精確率(Macro Average Precision,簡記Pr)、宏平均召回率(Macro Average Recall,簡記Re)以及宏平均F1分數 (Macro AverageF1,簡記F1)等指標對基于機器學習方法的巖性單元分類任務進行性能評估。其中,準確率表示正確預測的樣本比例,宏平均精確率表示預測為正樣本中正確的比例,宏平均召回率表示正樣本中預測正確的比例;宏平均F1分數是兼顧宏平均精確率與宏平均召回率的調和平均數。
假設混淆矩陣G:

其中,K表示巖性種類數。
在混淆矩陣G中準確率、宏平均精確率、宏平均召回率以及宏平均F1分數的計算公式:

其中,gaa表示a類巖性預測正確的數量;gab表示a類巖性預測為b類的數量。
利用上文所述方法在多龍礦集區開展巖性單元預測分類模型試驗,獲得了迭代過程各階段的模型評價指標。結果顯示,采用概率選區原則進行數據樣本逐步更新的思路具有良好表現,各指標隨迭代均保持遞增(表1)。以準確率為例,該指標表示當前分類結果與該區實際填圖獲得的巖性單元的匹配程度。模型經過7次迭代更新,準確率從初始47.3%增加至87%,性能提升近一倍。

表1 模型迭代性能統計表Table 1 Performance of model iteration
同時,結果顯示7次迭代后野外實際填圖的累計范圍占研究區面積的62.2%(表2),即,在全區約2/3范圍內開展野外填圖的情況下,獲得了與傳統填圖方法相近的巖性分類結果,說明文中提出的預測填圖方法在巖性填圖工作中的效率。根據已有地質圖可知,區內巖性單元種類數為20。巖性單元預測種類數在7次迭代過程中由13類增加至19類,覆蓋率達到95%。經統計發現,由于石英脈在研究區面積占比極少,僅為0.007%,缺少足夠的數據樣本,未能在研究中成功預測分類。由此可見,該研究方法從概率的角度定義迭代填圖范圍具有較高可行性。

表2 迭代分析結果信息統計表Table 2 Statistics table of iteration results
從7次迭代后的預測分類結果來看(圖5),在巖性單元分布較為復雜且多類型交替出現的場景下,相應的巖性單元邊界仍能被有效地劃分。該方法通過機器學習算法進行分類,提高了巖性單元填圖的工作效率。同時,與野外填圖結果對比發現具有較高的吻合度,體現了對巖性單元預測分類的準確性。

圖5 多龍礦集區巖石單元預測結果Fig.5 Prediction results of lithologic units in the Duolong mineral district
模型試驗經過7次迭代后,預測概率達到預期要求,分類結果涉及19類不同巖性單元。采用宏平均F1分數對各類單元進行精度評價(表3),模型分類精度整體表現優秀,各類預測精度均值達到0.845,其中5類超過0.9,僅有1類不足0.7,占比5%。最高為, 達到0.935,且該巖性單元僅占全區面積的0.47%,這說明該方法對于研究區面積占比較低的巖性單元仍具備較高的識別能力。最低為βμ,宏平均F1分數僅有0.683,但具備同等地球化學元素組成的,其宏平均F1分數達到了0.8,反映了該方法雖然對以巖石結構特征命名的地質單元無法有效區分,但對具備相同地球化學元素特征的巖性大類仍具備較高準確度。此外,通過預測結果與已知結果對比發現,第四紀區域預測與原地質圖有一定差別。經遙感查證,在排除第四系沖積扇區域之后,原1∶5萬圖幅的第四系分布范圍內局部顯示了露頭出露,表明該方法對已有填圖工作有部分修正作用。由于地球化學元素反演礦化蝕變的天然優勢,該方法對蝕變區域的有效識別,可產生重要的經濟價值。

表3 模型分類精度表Table 3 Table of classification accuracy of the current model
文中提出了一種基于GBDT算法的巖性單元預測分類方法,將西藏多龍礦集區作為試驗區,以1∶5萬勘查地球化學數據為例,對巖性填圖方法進行了有益的探索。研究強調了野外地質填圖與基于機器學習預測分類方法的深度融合,以及野地質調查工作在巖性預測填圖工作中的重要性和不可或缺性。在強調野外地質調查重要性的基礎上,將巖性填圖工作融入了機器學習方法。通過小范圍野外人工填圖迭代更新數據與知識庫,從而對全區進行巖性單元預測分類工作。該方法是對巖性單元填圖工作思路和流程的探索,是對現有工作模式的一種有益補充與輔助優化;體現了“基于大數據理論方法來促進地質問題的解決,并不意味著取代或摒棄地學傳統方法,而在于激活、提升和創新發展傳統方法”這一大數據科學范式在地質科學研究中的特點和優勢。
傳統巖性填圖方法通常要求對穿越地質體最多、地質構造復雜的路線進行復雜詳盡的野外調查工作,文中采用概率選區的方式來確定迭代填圖過程中的目標填圖區,使整個巖性填圖過程更具有針對性與高效性。根據試驗結果對比研究區
地質圖,該方法基于62.2%的已知研究區信息,有效實現了87%研究區范圍內的巖性單元分類。這一結果證明該方法不僅具有良好的填圖效果,而且能夠有效減輕野外填圖工作量。對在新疆、青海、西藏等野外環境條件艱苦地區的巖性填圖工作具有積極的參考作用。此外,為驗證該方法的通用性,未來可開展除化探數據以外其他數據資料,如遙感、航磁、航放、鉆井等數據資料的適用性研究,從而共同為地質資料相對匱乏或單一的研究區開展巖性填圖工作提供有效支撐。