胡　飛, 陳建國, 曾小華, 黎　蓉, 李　靖, 李　鵬

(1.中國地質(zhì)大學(xué) 資源學(xué)院,湖北 武漢　430074; 2.湖北省地質(zhì)調(diào)查院,湖北 武漢　430034; 3.湖北省地質(zhì)災(zāi)害防治中心,湖北 武漢　430034;4.湖北省地質(zhì)環(huán)境總站,湖北 武漢　430034)

0　引言

分形已被廣泛應(yīng)用到自然科學(xué)的眾多領(lǐng)域,在地球化學(xué)領(lǐng)域中,主要應(yīng)用分形技術(shù)進(jìn)行化探異常下限的確定,成秋明[1-3]指出化探數(shù)據(jù)具有多重分形結(jié)構(gòu)特征,其背景值往往服從正態(tài)分布或?qū)?shù)正態(tài)分布,而高低異常值服從多重分形分布。與統(tǒng)計方法相比,分形方法具有將樣品空間分布和統(tǒng)計特征隨空間尺度變化性納入化探數(shù)據(jù)的分形過程中,本質(zhì)是采用含量與頻數(shù)的關(guān)系來表達(dá),即統(tǒng)計大于某一含量的數(shù)據(jù)個數(shù),主要是含量-面積法(C-A模型)和含量-頻率法(C-S模型)。MAPGIS 6.7系統(tǒng)是中國地質(zhì)調(diào)查局指定使用的地質(zhì)工作軟件平臺,廣泛應(yīng)用于各地勘單位,利用MAPGIS平臺實現(xiàn)化探數(shù)據(jù)分形處理意義重大。

1　分形C-A模型

分形C-A模型是基于分形理論,研究含量與曲面面積之間關(guān)系的方法,并將此方法用于化探異常的研究。分形是指局部與總體具有某種相似性的形狀,或者是在不同尺度上看起來基本相似的形狀,即具有自相似性(self-similarity),而自相似性就是局部與整體在形態(tài)、功能和信息等方面具有統(tǒng)計意義上的相似[4]。定量描述這種自相似性的參數(shù)稱為“分維數(shù)”,用D表示。 分形分布的特點要求大于等于某一尺度的數(shù)目與物體大小之間存在冪指數(shù)關(guān)系,由此可得到含量-面積分形統(tǒng)計模型,即：

A(r)=Cr-D,r>0

①

式①中：r表示元素含量值；A(r)表示元素含量≥r的面積值;C>0表示比例常數(shù);D>0表示分維數(shù)。

分形C-A模型的關(guān)鍵在于確定分維數(shù),具體方法是將公式①兩端同時取對數(shù),得出線性回歸方程②。

lgA(r)=-Dlgr+lgC

②

把數(shù)據(jù)[A(r1),A(r2),…,A(rn)]和(r1,r2,…,rn)分別代入公式②,采用最小二乘法進(jìn)行直線擬合,計算斜率D,一般情況散點會分別在兩端或多段直線上,采用分段擬合方式,為提高精度避免人為確定分段拐點所造成的誤差較大,可采用最優(yōu)法確定分界點,使2個或多個擬合區(qū)間的直線與原始數(shù)據(jù)之間的殘差平方和Ei(i=1,2)在各區(qū)間的總和以公式③表示為最小。

③

式中：rj為分界點；D1和D2分別為相應(yīng)區(qū)間斜率,即分維數(shù),多段直線對應(yīng)D1,D2,D3…即多重分維,用維數(shù)確定異常下限后,把大于異常下限的數(shù)據(jù)刪除,利用所有異常點數(shù)據(jù)繪制異常圖。

2　MAPGIS 6.7平臺實現(xiàn)

MAPGIS地理信息系統(tǒng)是武漢中地信息工程有限公司研制的,具有自主版權(quán)的大型基礎(chǔ)地理信息系統(tǒng)軟件平臺,是集數(shù)字制圖、數(shù)據(jù)庫管理以及空間分析為一體的空間地理信息系統(tǒng)[5-6]。MAPGIS是目前國內(nèi)地勘單位廣泛使用的GIS平臺,在此平臺上實現(xiàn)化探數(shù)據(jù)多維分析處理顯得尤為重要。

在MAPGIS平臺上,用分析C-A模型確定化探異常下限的基本方法步驟,先將野外采集的化探數(shù)據(jù)進(jìn)行網(wǎng)格化,并繪制成平面等值線圖;然后利用MAPGIS屬性分析功能,統(tǒng)計出不同尺度的含量r值所對應(yīng)的面積A(r)值;最后利用GRAPHER軟件,將數(shù)據(jù)[A(r1),A(r2),…,A(rn)]和(r1,r2,…,rn)繪制在雙對數(shù)坐標(biāo)圖上,采用最小二乘法擬合成兩條或多條直線,從而求解出分維數(shù)D的估計值和對應(yīng)的拐點,拐點處所對應(yīng)的元素含量值即為該元素的異常下限[7]。

具體實現(xiàn)步驟如下:

(1) 化探數(shù)據(jù)按采樣坐標(biāo)X,Y,元素含量值保存為MAPGIS能導(dǎo)入的文本,如W.txt 。

(2) 利用MAPGIS空間分析功能DTM子模塊,對W.txt文件進(jìn)行網(wǎng)格化,選擇合適的網(wǎng)格化方法(克里格法、最小曲率、鄰近點法等),輸出W.grd網(wǎng)格文件;根據(jù)采樣數(shù)量設(shè)定合適的(r1,r2,…,rn)值,打開“平面等值線圖繪制”菜單,選擇W.grd文件,按r1,r2,…,rn設(shè)置參數(shù),生成區(qū)文件W.wp 。

(3) 利用MAPGIS空間分析子模塊,打開W.wp區(qū)文件,使用“屬性分析”的“雙屬性分類統(tǒng)計”功能將(r1,r2,…,rn)各個值段對應(yīng)的面積統(tǒng)計并導(dǎo)出為W.wb文件。

(4) 利用MAPGIS庫管理模塊轉(zhuǎn)換W.wb文件為W.dbf,在Excel軟件中導(dǎo)入的數(shù)據(jù)文件W.dbf,累加(r1,r2,…,rn)各個值段對應(yīng)的面積,求出A(r1),A(r2),…,A(rn),并利用Excel直接計算出lgri和lgA(ri)。

(5) 利用Excel或Grapher繪制lgri、lgA(ri)散點圖,利用最小二乘法公式③,求出分維數(shù)D1和D2及分界點含量值j,j即為元素異常下限值。

3　應(yīng)用實例

3.1　研究區(qū)地質(zhì)概況

3.1.1區(qū)域地質(zhì)特征

礦區(qū)位于桐柏—大別變質(zhì)核雜巖隆起亞帶中部。區(qū)域出露主要為TTG系列的英云閃長質(zhì)片麻巖、奧長花崗質(zhì)片麻巖、花崗閃長質(zhì)片麻巖組合;表殼巖系為沉積其上的含鐵碎屑建造、孔茲巖系及碎屑—碳酸鹽巖建造夾雙峰式火山建造(稱大別山(巖)群);而侵入巖則主要表現(xiàn)為晉寧期和燕山期的花崗巖類。

造山帶內(nèi)巖漿活動十分頻繁,發(fā)育多期巖漿巖,尤以燕山期后碰撞鈣堿性花崗巖最為強烈,工區(qū)北部有靈山、新縣和商城等巖體分布,規(guī)模較小的巖體主要形成于早白堊世,其中中酸性小巖體密布,巖性多為酸性富堿的花崗斑巖、似斑狀花崗巖、石英斑巖及花崗閃長斑巖等。

3.1.2礦床地質(zhì)特征

礦區(qū)出露地層屬秦嶺—大別地層區(qū)卡房—龜峰山小區(qū),出露主要巖性有黑云二長花崗質(zhì)片麻巖、斜長角閃巖、斜長角閃片巖等。廣泛出露的二長花崗質(zhì)片麻巖以及基性火山巖組,巖石化學(xué)性質(zhì)十分穩(wěn)定,巖石中節(jié)理、裂隙發(fā)育,為良好的成礦圍巖,也提供了儲礦空間。區(qū)內(nèi)巖漿巖活動頻繁,燕山期巖漿侵入尤為強烈,提供了豐富的巖漿熱液和成礦物質(zhì)來源。區(qū)內(nèi)構(gòu)造復(fù)雜,斷裂構(gòu)造極其發(fā)育,它們控制熱液活動、地表水和地下水運移,控制成礦流體的運移,提供有利的導(dǎo)礦通道及成礦空間,進(jìn)而控制礦體的空間分布。區(qū)內(nèi)構(gòu)造以斷裂活動為主。按照空間展布特征可劃分為北西西向、近東西向、北北東向三組,以北西西向較為發(fā)育,它們組成了區(qū)內(nèi)網(wǎng)格狀的斷裂構(gòu)造格架,基本與區(qū)域構(gòu)造配套。其中兩路口斷層、獅子垴斷層為工區(qū)主要斷層。兩路口斷層區(qū)內(nèi)出露長度約3 km,走向10°～30°,斷層產(chǎn)狀310°～340°∠45°～78°,形成的破碎帶一般寬5～15 m,帶內(nèi)發(fā)育強硅化巖、構(gòu)造角礫巖、碎裂巖及碎粒巖,明顯與圍巖面理斜交,常見斷層引起的拖曳褶曲,構(gòu)造巖具明顯分帶現(xiàn)象,從內(nèi)到外依次為碎粉巖、角礫巖、碎裂巖,部分有早期糜棱巖殘留,斷裂帶及其附近節(jié)理發(fā)育,沿斷層常常見石英細(xì)脈沿裂隙面充填及黃鐵礦化,該斷層是區(qū)內(nèi)規(guī)模最大的斷裂構(gòu)造,是商城—麻城斷裂的次級斷裂,且兩路口鎢礦床主礦體分布于其構(gòu)造破碎帶內(nèi)。

3.2　確定異常下限

以礦區(qū)土壤化探數(shù)據(jù)為例,有樣品2 600個,分析Ag、As、Bi、Mo、Pb、Sn、W、Zn元素,按C-A模型利用MAPGIS軟件進(jìn)行計算求出r和A(r),取對數(shù)部分統(tǒng)計結(jié)果如表1。

根據(jù)統(tǒng)計數(shù)據(jù)制作散點圖,如圖1。

以Ag元素為例,利用公式②③求擬合兩段直線方程為:

lgA(r)= 0.272lgr+6.859 2

(-3.91

lgA(r)=-2.627lgr-0.090 2

(-2.95

因此分維數(shù)為D1=0.272,D2=2.627,異常下限r(nóng)0=0.052 3,分界點r0的地質(zhì)意思是含量<0.052 3為背景數(shù)據(jù),>0.052 3為異常數(shù)據(jù)。經(jīng)分析表明:

(1) 在雙對數(shù)(lgr-lgA)圖中,元素含量與面積的投影點呈現(xiàn)曲線分布,可用兩條或多條直線去擬合,顯示出一種連續(xù)多重分型特征。其中第一段直線非常平坦,D1數(shù)值接近零,大致反應(yīng)了元素含量的低值背景,第二段直線斜率較陡,反應(yīng)元素含量的高值區(qū)及異常區(qū),兩條直線的拐點對應(yīng)元素含量可視為異常下限。

(2) 圖中分維值D定量反應(yīng)了元素含量在空間上的叢集集度和不均勻程度,D值越大,則元素空間分布均勻,且主要集中在低含量區(qū),在地質(zhì)上成礦性能極小,如As元素D2=8.48;D值越小,反應(yīng)元素低值含量點到高值含量點的頻率下降得越慢,元素在空間上

分布越不均勻,存在著較多高含量點,在地質(zhì)上有富集成礦的可能,如W元素D2=1.04、Mo元素D2=1.43,在兩路口地區(qū)主攻礦種就是鎢鉬礦,分形分維值真實反映了地質(zhì)事實,經(jīng)長期分形數(shù)據(jù)處理筆者認(rèn)為D2值<2,尤其是<1.8有極大可能性富集成礦[8]。

(3) 圖中有些元素多重分維,分維數(shù)D3明顯小于D2,在圖像上曲線明顯偏離第二條直線段,如Mo、Pb、Sn、Zn元素,說明這些元素在高含量點基礎(chǔ)上有二次局部富集[9],在地質(zhì)上極有可能是成礦元素的伴生元素,且本區(qū)Ag與Pb、Zn富集沒有相關(guān)性,Mo與Pb、Zn富集有一定相關(guān)性,W、Bi富集具有較強相關(guān)性,Mo與W、Bi富集也存在相關(guān)性。

表1　lgr(含量)與lgA(面積)統(tǒng)計結(jié)果Table 1　Statistical results of lgr (content) and lgA (area)

圖1　兩路口礦區(qū)土壤化探元素C-A分形雙對數(shù)圖Fig.1　Chart of double logarithm of C-A fractal dimension of soil elements from geochemical survey in Liangluokou Mining Area

3.3　繪制異常圖

依據(jù)分維確定的拐點處對應(yīng)的元素值,可直接繪制元素異常圖,以W元素為例,依據(jù)多重分維D1、D2、D3,確定區(qū)域異常下限為7.39 μg/g,局部異常(礦致異常)下限為40 μg/g繪制異常圖,如圖2。

圖2　W元素分形異常圖Fig.2　Fractal anomaly of W element

圖3　W元素異常圖Fig.3　Anomaly map of W element

對比傳統(tǒng)方法繪制的異常圖,如圖3。分形C-A模型處理方法保留了所有異常數(shù)據(jù),其區(qū)域異常能較好反應(yīng)礦區(qū)構(gòu)造特征,體現(xiàn)W元素富集與東西向構(gòu)造的密切關(guān)系;其局部異常面積小于傳統(tǒng)方法,更為有效地體現(xiàn)了最小面積、最大含礦率原則。

4　結(jié)論

基于分形C-A模型,確定化探異常下限主要有以下特點：

(1) 反映化探數(shù)據(jù)元素特征。傳統(tǒng)方法確定地球化學(xué)元素異常下限值的前提,是元素含量呈正態(tài)分布或者對數(shù)正態(tài)分布,所以傳統(tǒng)方法通常需要多次剔除特高值,不斷檢驗是否符合正態(tài)分布,即“檢驗—計算—再檢驗—再計算”的過程,分形C-A模型法確定異常下限值時,無需數(shù)據(jù)服從正態(tài)分布,也不用對數(shù)據(jù)進(jìn)行特高值剔除處理,大大降低了人為因素的干擾,能較客觀地反映研究區(qū)成礦元素的分布形態(tài)。

(2) 異常下限確定有理論依據(jù)。對服從多維分形的化探數(shù)據(jù)運用最小二乘法擬合曲線,斜率不同的直線段反映的是不同期次的地球化學(xué)場,所對應(yīng)的分界值往往可作為區(qū)分背景和異常的臨界值[10-11],這為研究地球化學(xué)元素的分布規(guī)律及其異常下限確定提供了理論支持,且元素分維值D2越小表示該元素越有富集成礦的可能,可利用分維值快速確定主攻礦種。

(3)C-A模型分析實現(xiàn)容易。僅利用MAPGIS和Excel就能實現(xiàn)化探數(shù)據(jù)分形C-A模型分析,操作簡單易行。結(jié)合傳統(tǒng)方法確定異常下限進(jìn)行對比分析,能取得較好的效果,提高礦產(chǎn)預(yù)測的可靠性和客觀性,對指導(dǎo)礦區(qū)下一步工作具有參考意義。

參考文獻(xiàn)：

[1]Cheng Q,Agterberg F P,Bonham-Carter G F.A spatial analysis method for geochemical anomaly separation[J].J Geochem Explor,1996，56:185-195.

[2]Cheng Q.Multifractality and spatial statistics[J].Computer & Geosciences,1999,5(9):949-962.

[3]Cheng Q,Agterberg F P,Bonhaml,Ballantyne S B.The separation of geochemical anomalies from background by fractal methods[J].J Geochem Explor,1994，51:109-130.

[4]成秋明.空間模式的廣義自相似分析與礦產(chǎn)資源評價[J].地球科學(xué),2004,29(6):733-743.

[5]張釗,韋龍明,陳三明.MAPGIS在地質(zhì)填圖及化探數(shù)據(jù)處理中的應(yīng)用[J].物探化探計算技術(shù),2010,32(2):119-122.

[6]吳信才.MAPGIS地理信息系統(tǒng)[M].北京:電子工業(yè)出版社,2007.

[7]楊震,梅紅波.基于MAPGIS的含量—面積法確定異常下限：以個舊西區(qū)為例[J].物探化探計算技術(shù),2013,35(4):477-480.

[8]李堃,胡光道,劉才澤.個舊地區(qū)地球化學(xué)元素的多重分型特征與找礦預(yù)測[J].礦產(chǎn)地質(zhì),2006，20(4/5):498-502.

[9]張慶敏,岳云娟.基于小波模極大值求取多重分形譜多重分形克里格算法探究[J].魯東大學(xué)學(xué)報：自然科學(xué)版,2013(1):22-25.

[10]陳建國,王仁鐸,陳永清.利用分形統(tǒng)計學(xué)提取化探數(shù)據(jù)中的隱蔽信息并圈定地球化學(xué)異常[J].地球科學(xué),1998,25(3):311-318.

[11]成秋明,趙鵬大,陳建國,等.奇異性理論在個舊錫銅礦產(chǎn)資源預(yù)測中的應(yīng)用:成礦弱信息提取和復(fù)合信息分解[J].地球科學(xué),2009,34(2):232-242.