, ,,
(1.陜西省煤田地質局一八五隊,陜西 榆林 719000; 2.中煤科工集團西安研究院,陜西 西安 710077; 3.中國地質大學(武漢)資源學院,湖北 武漢 430074; 4.中國地質大學地質過程與礦產資源國家重點實驗室,湖北 武漢 430074)
現代礦床學研究表明,眾多礦床的形成具有多來源、多階段和多成因的特點 (翟裕生等,1996 )。故而可將形成一個礦床的過程看作是多個總體的混合或疊加,在數學上這種混合或疊加可以用混合總體篩分的數學模型來描述。研究勘查數據的混合分布特征可以幫助揭示多期次、多重疊加的地質作用與地質過程,是勘查數據統計分布特征研究的重要內容之一。為了深入研究地質體(如巖體、礦體)的成因或形成過程,需要從混合分布的樣本數據中將各成分總體分離出來(或估計出來),對成分總體的有關參數進行估計,這個過程稱為混合分布的篩分?;旌峡傮w的篩分方法可分為3類:解析法(McLachlan et al,2000)、圖解法(辛克萊,1976; 趙鵬大等,1983; 王琳,2005)、數值法。數值法是近年來的主流方法,它運用各種數學最優化理論方法(如最小二乘法、最大似然法等)進行混合總體的篩分。EM算法目前是數值法中發展最快的一種(劉向沖等,2011; 張鋒,2012)。EM算法及其擴展算法如MML-EM算法等為勘查地球化學數據、金屬礦床品位分析數據的混合分布篩分以及解釋地質成因提供了一種快捷優良的定量化工具,然而在煤田地質中的應用還比較少見。
掌握煤層厚度變化規律,對于礦區煤炭資源估算及生產建設具有重要意義。前人對煤層厚度(簡稱“煤厚”)空間變化性總體上認為沉積環境控制煤厚區域變化,構造主要引起煤層厚度的局部變化(琚宜文等, 2002; 高榮斌等2011)。實際上煤厚往往呈現更加復雜多樣的變化,它往往不是單一因素造成的,而常常是幾種因素綜合作用的結果,基底的沉降幅度和速度、沉積物的補給及植物遺體的堆積速度、成煤期后的構造活動對煤層厚度及其變化都具有一定影響(李恒堂等,1995)。在同一礦井中的不同煤層或同一煤層的不同部位,引起煤厚變化的原因也可能不同,有些地方以原生變化為主,而另一些地方則可能因后發生沖蝕或構造擠壓引起。因此,煤厚變化通常是多種因素復合作用的結果,并且不同因素所起的作用還有著強弱上的差別,這些因素共同作用的信息可以體現在煤厚勘查數據中。因此,應用混合分布理論與方法來研究煤厚主要控制因素是可行的。
本研究以陜北侏羅紀煤田榆神礦區中雞勘查區83個鉆孔中測量所得的3個主要煤層(2-2(2-2上)煤、3-1煤、4-3煤)的煤厚數據為例,研究煤厚的混合分布特征,定量描述煤層的變化規律并進一步分析析混合分布中各子總體的地質意義。
地質勘查數據的多總體性特點具有普遍性,在開展統計分析中不容忽視。有限混合模型提供了一種用簡單分布擬合復雜分布的強勁靈活的概率分析數學工具。不失一般性,假定樣本的概率密度模型為(McLachlan et al,2000):
(1)

有限混合分布中,當模型分支數為1時,混合分布總體是單一分布,反映數據是同質的;當分支數多于1時,混合分布總體是多成分的,反映數據是異質的?;旌戏植伎傮w通??捎啥鄠€單一分布以某種權重疊加來模擬,這些單一分布與某一類型的端元組分相聯系,可能指示著某種內在的控制因素。因此,如何有效獲取混合分布數據中的分布類型及其參數、分支個數、混合權重,既是有限混合分布模型研究的核心問題,也是地質勘查數據統計分布研究的重要內容。
高斯混合模型以其形式簡單、計算方便等特點,受到普遍采用。盡管地質勘查數據統計大多具有非線性、非高斯特性,然而,相對于勘查地球化學微量元素數據而言,煤厚數據的非高斯特性通常表現不是很強烈。另外,混合高斯概率密度函數模型具有廣泛適用性,原理在于:當樣本足夠大時,混合分布中各子分布漸進分布為正態分布,選擇合適的模型參數,混合高斯概率密度可對任意形狀概率密度函數進行精確擬合(柳貴東等,2011)。因此,可以簡化研究的復雜度,采用高斯混合模型來開展本研究,即利用一定數量的正態分布概率密度函數,通過線性疊加混合來逼近煤厚變量的概率密度函數。在有限混合分布模型現有的參數估計方法中,基于極大似然估計原理的EM算法近年來最受關注。
EM算法(Expectation-Maximization Algorithm,期望最大化算法)是一種從“不完全數據”中求解模型參數的極大似然估計的方法。它提供了一個高效的迭代算法用來計算這些數據的最大似然估計(張士峰,2004)。EM總體算法流程為:初始化分布參數,然后反復迭代直到滿足收斂條件,終止迭代。EM算法的每一次迭代過程都分為2個步驟:(1) E步:估計未知參數的期望值,給出當前的參數估計;(2) M步:重新估計分布參數,以使得數據的似然性最大,給出未知變量的期望估計。當迭代結束,最終可求出使極大似然值最大的1組參數,此組參數即作為混合分布概率密度函數未知參數的極大似然估計解。
EM算法由Dempster等(1977)提出,現已成為混合模型擬合的主流方法,并出現了基于EM算法的多種變體方法。經典的EM算法解決了如何估計混合模型參數的難題,然而混合模型分支數量確定問題尚未得到徹底解決。目前最典型的方法是Figueiredo等(2002)提出的基于MML(Minimum Message Length Criterion,最小信息編碼)準則和EM框架的算法,通常稱為MML-EM算法。該算法將參數估計和模型選擇緊密結合到同一處理流程中,能同時處理分支數和模型參數這2個問題,其算法思想為:在EM算法迭代過程的E步驟和M步驟之后增加L步驟,即利用前2步計算參數計算相應的信息長度,重復上述第3步,直到求出在信息長度最小條件下的最優分支數和分布參數。詳細公式理論推導和算法實現見文獻(Figueiredo et al, 2002)。此外,還有其他角度對EM算法進行改進,尋求最優分支數,如基于遺傳算法的GA-EM算法等(連軍艷, 2006)。本研究將采用MML-EM算法進行陜北中雞勘查區煤厚數據混合分布特征研究。
中雞勘查區位于陜北侏羅紀煤田的北部,與神府礦區相鄰,是鄂爾多斯盆地中生代中侏羅世含煤沉積——陜北侏羅紀煤田的一部分。其行政區劃隸屬于榆林市神木縣中雞鎮所管轄,北部靠近內蒙古自治區。該地區煤質優良,煤層多且厚,煤炭資源儲量豐富,開采地質條件簡單,有廣闊的發展前景和巨大的開發潛力,是陜北能源化工基地的重要組成部分(陜西省煤田地質局一八五隊,2012)。
中雞勘查區所處的鄂爾多斯盆地是以中生代陸相沉積為主體的大型內陸沉積盆地,構造單元處于鄂爾多斯寬緩的東翼——陜北斜坡上。盆地基底是堅固的前震旦系結晶巖系,故成煤前后的整個地質發展過程繼承了深部基底的穩定性。中生代以來,地質史上歷次構造運動對本區影響甚微,以垂向運動為主,形成了一系列假整合面,沒有發現火成巖,發現少量斷層。
中雞勘查區內地表大部分被現代風積沙、第四系薩拉烏蘇組所覆蓋,基巖中北部出露。地層總體為走向北東、傾向北西、傾角<1°的單斜構造,未發現明顯的褶皺構造,未發現落差>50 m的斷層,也無巖漿活動。地層由老到新依次有:三疊系上統永坪組(T3y),侏羅系中統延安組(J2y)、直羅組(J2z)、安定組(J2a),白堊系下統洛河組(K1l),新近系上新統保德組(N2b),第四系中更新統離石組(Qpl)、第四系上更新統薩拉烏蘇組(Qps)、第四系全新統風積層(Qheol)和沖積層(Qhal)。其中含煤地層延安組(J2y)巖性以灰白色-淺灰白色粗、中、細粒長石石英砂巖、巖屑長石砂巖及鈣質砂巖為主,次為灰-灰黑色粉砂巖、砂質泥巖、泥巖及煤層,少量炭質泥巖,局部地段夾有透鏡狀泥灰巖及黃鐵礦結核。與下伏永坪組為假整合接觸,與上覆直羅組呈假整合接觸,厚度147.02~235.89 m,平均厚度201.20 m,總體變化趨勢由中部向四周逐漸增厚。
延安組可劃分為5個煤組,煤層編號見表1,其中2-2上及2-2下是2-2煤層分岔后上下分層的編號。本區可采煤層有9層,主要可采煤層為1-2下、2-2(2-2上)、3-1、4-3、5-2上,次要可采煤層為1-2下、2-2下、4-2、5-2下。以勘查區內83個鉆孔的2-2煤(2-2上)、3-1煤、4-3煤的煤厚數據為例,運用MML-EM算法,對煤厚的混合分布特征進行研究,其中2-2煤分岔后煤厚數據取自2-2上煤,為方便敘述,代號統一簡記為2-2煤。

表1 中雞勘查區延安組分段及煤層編號
傳統統計分析方法是以數據總體滿足正態假設為依據,并在此基礎上建立模型和統計推斷。為“讓數據說話”,采用探索性數據分析(EDA)方法來展示煤厚數據的統計特征及空間分布特點,在試探出數據的統計特點后再對數據進行進一步分析。注意到本區83個鉆孔中有若干鉆孔煤層不出現,2-2煤、3-1煤、4-3煤不出現煤層的鉆孔數分別為:2,3,9個,于是,首先剔除2-2煤、3-1煤、4-3煤的煤厚為0的數據。圖1是該3層煤的煤厚數據的矩陣散點圖-直方圖-箱線圖的疊加綜合統計圖。由3個煤厚變量的直方圖可見:2-2煤、3-1煤、4-3煤均呈現多峰型混合分布特征,從箱線圖可見3個煤厚變量的箱線圖均反映出煤厚分布是非正態的,并且中位數偏向高值端;3-1煤的煤厚具有若干異常值,特別是特低值成分更為明顯。圖1中的散點圖反映了任意2個煤厚變量之間的相關關系,從中可見,2-2煤與4-3煤之間存在一定的相關性,并且按高值和低值可以分成2組點群;另外2組變量之間的相關性則不明顯。

圖1 2-2煤、3-1煤、4-3煤煤厚數據的探索性數據分析
將探索性數據分析技術數據用于空間統計中,產生了探索性空間數據分析(ESDA)技術。將2-2煤、3-1煤、4-3煤煤厚數據所在的鉆孔空間位置以及3個變量之間的煤厚柱狀圖繪制投影到勘查區地形圖上。柱狀圖中煤厚為一水平線者表示對應煤層未出現,在圖中對應的鉆孔符號顏色設置為紅色(圖2)。由西北向東南,2-2下煤的出現意味著2-2煤開始出現分岔,該分岔界線呈近北東向延伸(見圖2中黃綠色虛線AB)。由鉆孔資料分析可知,分岔煤層間距約 0.93~32.75 m,平均16.57 m,總的趨勢是由北向南間距逐漸增大,到本區東南角又開始有復合趨勢。由圖2還可見:在該分岔界線AB兩側,2-2煤、3-1煤、4-3煤煤厚具有良好的空間變化規律。2-2煤在本區西北部與中部(AB線左上側)煤厚普遍較大,且明顯厚于3-1煤和4-3煤;而在東南部(AB線右下側),2-2煤厚度明顯降低,該區以3-1煤煤厚相對占優。4-3煤主要出現在本區西北部和中部;東南部則厚度較薄,甚至未見,未見4-3煤的9個鉆孔均出現在該區。在西北部與中部,4-3煤厚與3-1煤厚的比值以大于1為主,而在東南部則小于1;2-2煤厚與3-1煤厚的比值也有類似的結果。
由探索性數據分析結果可知,中雞勘查區2-2煤、3-1煤、4-3煤煤厚數據具有混合總體性,其空間分布特征比較復雜但又有一定規律可循。應用混合分布理論來開展煤厚控制因素研究,為對混合分布形態及分支數有個經驗性的認識,繪制了正態概率圖來展示2-2煤、3-1煤、4-3煤煤厚數據的混合總體累積概率特征(圖3)。按煤礦生產對煤厚進行分類結果表明:2-2煤以厚煤層為主,薄煤層為次,中厚煤層相對較少;3-1煤與4-3煤以中厚煤層為主,其次為薄煤層,部分位置還出現極薄煤層;4-3煤的中厚煤層的層厚以2 m以內占主導,而3-1煤的中厚煤層的層厚在2~3 m區間上有個較明顯的集中區。由圖3展示的正態概率圖可見,2-2煤、3-1煤、4-3煤煤厚數據具有混合分布特征,不能以單一分布如正態分布來描述,混合分支數量可能為2~3個,總體而言分支數量不會太多,各分支的概率密度函數可以用正態分布概率密度函數來刻畫,采用高斯有限混合模型來開展混合分布篩分。
MML-EM算法是EM算法的改良版,采用它來估計子總體的最優個數及正態概率密度函數參數。設迭代收斂時的精度為10-8,分支數的搜索范圍為[1,3],經過迭代計算得到最終參數估值。2-2煤、3-1煤厚數據均篩分出2個分支分布,4-3煤篩分出3個分支分布(表2、圖4)。為便于討論,對2個分支分布情形,按其子總體均值大小分別稱為高值子總體和低值子總體,對3個分支情形,則均值位于中間者所對應的第2個分支分布稱為中值子總體。

圖2 中雞勘查區探煤鉆孔位置及煤厚柱狀圖

圖3 2-2煤、3-1煤、4-3煤的煤厚數據正態概率圖

表2 層厚數據高斯混合總體篩分結果

圖4 層厚數據高斯混合總體篩分結果
應用GeoDA軟件(Anselin,2005)對煤層厚度與其底板高程進行了探索性空間數據分析,以期發現煤厚混合分布子總體樣本與底板高程之間的統計關系及空間展布特點。圖5展示了由GeoDA軟件繪制的粒狀示意地圖及兩者的散點圖,
圖5a示意地圖上,各點符號位置近似于鉆孔空間位置,點符號大小表示煤厚大小,填充顏色表示煤層底板高程,通過點符號輪廓的差異區分不同的子分布樣品。對比2-2煤、3-1煤、4-3煤3層煤底板高程的變化趨勢,可見它們有相同的變化趨勢,即底板高程由北西向東南總體上逐漸增高,這與“陜北斜坡”的構造形態基本吻合,值得注意的是,3層煤底板高程在不同位置有不同程度的凹陷。2-2煤以分岔界線AB(圖5a中的AB線,與圖2中的AB線相同)為界正好對應篩分出的2個子總體,西北及中部鉆孔所見2-2煤厚屬于高值正態總體,分岔后在東南部的鉆孔所見2-2煤厚屬于低值正態總體,這些低值總體樣品總體上分布于勘查區東南部的凹陷區,呈條帶狀、近北東向展布;在AB線附近鉆孔所見的3-1煤厚明顯偏低,對應該煤厚混合分布中的低值正態總體,這些低值總體樣品分布于勘查區中部3-1煤的凹陷區,呈串珠狀、近北東向展布;4-3煤的凹陷區分布在東南角上,呈條帶狀、近東西向展布,這些位置的煤厚數據服從低值正態總體分布,AB線穿越的區域對應篩分所得子總體為中值子總體。
圖5b散點圖所展示的煤厚與底板高程顯示出一定的規律性。在PQ線左上側范圍內,2-2煤、4-3煤的煤厚與底板高程呈負相關關系,即底板高程越高的地方,煤層越?。?-1煤煤厚則受底板高程影響較小。在PQ線右下范圍內,3-1煤煤厚與底板高程呈正相關關系,2-2煤煤厚也具有類似關系,2-2煤在該區域由于樣品數較少,與底板高程的關系不易識別??傮w而言,這3個煤層的煤厚混合數據在PQ線右側與高程的正相關關系還是相對較清晰的,2-2煤以極薄煤層形式出現。
綜上所述,底板高程與煤厚有密切關系,高程的變化可能是影響煤厚變化的主要因素。值得指出的是,圖5b所示PQ線對于劃分2-2煤、4-3煤煤厚樣品是有效的,而3-1煤煤厚與其底板高程的關系相對其他2層煤而言表現形式要復雜一些,PQ線所劃分的樣品與本次研究篩分所得的3個子總體并不對應。若按PQ線將3-1煤樣品劃分成2個子總體,則煤厚數據近似服從“相交雙峰型”混合分布,即其中一個子總體完全重疊在變化范圍更寬的另一個子總體之上,所采用的一維EM算法會失效。采用煤厚單變量來研究混合分布,根據圖4b中3-1煤煤厚直方圖特征,對該煤厚數據理解為近似服從“非相交雙峰型”混合分布。在進一步研究中可應用二維EM算法對煤厚及其底板高程2個變量開展混合總體篩法,并探討其地質意義。

圖5 煤層底板高程與煤層厚度變化
根據以往資料,整個陜北侏羅紀煤田構造簡單、穩定,沒有大的褶皺和斷裂,在此大背景下,可從沉積與剝蝕的角度來分析煤層厚度與其底板高程間關系的形成原因。由上述分析可知,在凹陷區及附近煤厚突然變薄是2-2煤、3-1煤、4-3煤共同的變化趨勢特點之一。不同的煤層,凹陷區分布位置有所不同。3-1煤凹陷區靠近勘查區中部AB線附近,2-2煤為東南部,4-3煤在東南角上,因此導致煤厚趨勢變化隨空間位置而不同,煤厚與煤層底板總體上表現為正相關關系。在沉積基底低凹的地區,由于地勢突然降低,沉積物的補給量和沉積速度很快(黃克興等,1991),因此在沉積過程中,泥沙及礫石等沉積物的沉積速度和補給量總體上大于植物遺體,搶占了植物遺體的補給空間,于是出現現存煤層底板低的區域煤層反而薄的情形,2-2煤的凹陷區出現2-2下煤可能就是這種原因造成的。隨著底板地勢的逐漸增高且沉積物的補給量和速度有所減緩,植物遺體獲得了最佳的沉積機會,于是煤層隨著其底板的增高不斷變厚,因此可形成本區中的中厚煤層。
在本區更廣泛出現的是煤厚與煤層底板呈負相關關系的情形,所涉及煤層含煤性更好。負相關關系的原因可能是由于在區域上,隨著底板高程的持續增高,植物遺體來源不斷減少而沉積速率逐漸降低,于是導致隨著煤層底板增高煤層逐漸減薄的狀況。圖5b中PQ線左側,2-2煤、4-3煤底板高程變化幅度可達150 m,煤厚變化幅度也可達3 m左右。在植物遺體來源充足的情況下,結合本區良好的沉積-構造條件,易于形成中厚煤層和厚煤層。
通過圖5a所展示的混合分布正態子總體所屬樣品的空間位置,發現每個子總體在空間上都具有良好的聚集性。2-2煤厚以AB線為界(左上側為高值子總體樣品,右下側為低值子總體樣品);4-3煤厚以AB線附近(低值子總體樣品)與遠處(高值子總體樣品)來劃分,4-3煤在本研究中篩分出3個子總體,在空間上由北西向東南依次為高值子總體—中值子總體—低值子總體,其中中值子總體在空間上處于中間過渡部位,AB線出現在區域內部??傮w而言,2-2煤、3-1煤、4-3煤煤厚數據經混合分布篩分都獲得了一個極薄—薄層煤厚的低值子總體以及中厚—厚層煤厚的高值子總體,兩者可能代表了2種不同的沉積環境,即低值總體指示了沉積物的沉積速度和補給量大于植物遺體的沉積速度和補給量;高值總體正好相反,代表著低值總體代表沉積物的沉積速度和補給量小于植物遺體的沉積速度和補給量。
中雞勘查區2-2煤以厚煤層為主;3-1煤以中厚煤層為主,局部不可采;4-3煤以中厚煤層為主,局部不可采。2-2煤、3-1煤、4-3煤的煤厚數據具有混合分布特征,采用一維MML-EM算法獲得了該3個煤層的煤厚數據的子總體個數及其參數,其中,2-2煤、3-1煤的煤厚數據近似服從由2個子分布組成的混合正態分布,4-3煤厚數據近似服從由3個子分布組成的混合正態分布。
EDA/ESDA技術為發現與理解煤厚數據統計特征、空間分布特征提供了有效的分析手段,特別是良好的可視化手段。初步認為煤厚變化的主要控制因素是底板高程的變化,其中既有正相關關系,也有負相關關系,以前者在本區占主導。正相關關系所指示的煤層厚度達到中厚—厚層級別,煤層主要分布在本區的西北和中部區域;負相關關系所指示的煤層厚度通常較薄,在地理位置上受控于局部凹陷。煤厚數據混合總體篩分獲得了低值與高值2個子正態子總體,它們可能分別對應于2種不同的沉積環境,即低值總體指示泥沙及礫石等沉積物的沉積速度大于植物遺體的沉積環境,高值總體則反之,植物遺體的沉積速度和補給量占主導。
中雞勘查區煤厚數據具有多峰型混合總體分布特征,2-2煤、4-3煤表現為典型的非相交雙峰型混合分布;3-1煤可能為非相交雙峰型混合分布,也可能為相交雙峰型混合分布??煽紤]煤厚及其底板高程的二維混合總體分布的篩分,進一步開展關于混合總體分布的研究。
杜文鳳,彭蘇萍.2010.利用地質統計學預測煤層厚度[J].巖石力學與工程學報,29(增刊1):2762-2767.
高榮斌,賀志強,來爭武,等.2011.豫西新安煤田煤層厚度變化規律及其控制因素[J].煤田地質與勘探,39(4):13-15,19.
黃克興,夏玉成.1991.構造控煤概論[M].北京:煤炭工業出版社.
琚宜文,王桂梁,胡超.2002.海孜煤礦構造變形及其對煤厚變化的控制作用[J].中國礦業大學學報,31(4):374-379.
李恒堂,呂志發.1995.鄂爾多斯盆地延安組控煤古構造趨勢分析[J].煤田地質與勘探,24(5):5-8.
連軍艷.2006.EM算法及其改進在混合模型參數估計中的應用研究[D].西安:長安大學.
劉向沖,侯翠霞,申維,等.2011.MML-EM方法及其在化探數據混合分布中的應用[J].地球科學:中國地質大學學報, 36(2):355-359.
柳貴東,山拜·達拉拜.2011.基于EM算法的非高斯噪聲參數估計[J].通信技術,44(1):151-153.
苗霖田等.2008.神木北部礦區5—2煤層厚度及其底板高程趨勢分析[J].煤田地質與勘探,36(3):12-15.
陜西省煤田地質局一八五隊.2012.陜西省陜北侏羅紀煤田榆神礦區中雞勘查區詳查報告[R].榆林:陜西省煤田地質局一八五隊.
王琳.2005.可視化概率圖解法軟件的研制與應用[D].北京:中國地質大學(北京).
辛克萊.1976.概率圖在礦床勘探中的應用[M].北京:地質出版社.
翟裕生,姚書振,崔彬.1996.成礦系列研究[M].北京:中國地質大學出版社.
張鋒.2012.閃鋅礦礦石標本便攜式XRF測量數據混合總體篩分及其地質意義[D].武漢:中國地質大學(武漢).
張士峰.2004.混合正態分布參數極大似然估計的EM算法[J].飛行器測控學報,23(4):47-52.
張展適,吳信民.1998.概率圖法在茅排金礦的應用[J].華東地質學院學報,21(3):254-256.
趙鵬大,胡旺亮,李紫金.1983.礦床統計預測[M].北京:地質出版社.
ANSELIN L.2005-03-06.Exploring Spatial Data with GeoDaTM:A Workbook[M/OL]. http://geodacenter.asu.edu/system/files/geodaworkbook.pdf.
DEMPSTER A P, LAIRD N M, RUBIN D B.1977.Maximum likelihood from imcomplete data via the EM alogorithm[J].Journal of the Royal Statistical Society:Series B,39:1-38.
FIGUEIREDO M A T,JAIN A K.2002.Unsupervised learning of finite mixture models[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,24(3):381-396.
MCLACHLAN G,PEELD.2000.Finite Mixture Models[M].New York,USA:John Wiley & Sons Inc.