膠西北金礦集區地球化學多重分形模型與異常識別

毛先成， 張顯洋， 劉占坤， 陳 進， 鄧 浩

(中南大學 a.地球科學與信息物理學院；b.有色金屬成礦預測與地質環境監測教育部重點實驗室，長沙 410083)

0 引言

膠西北礦集區是我國重要的金礦產出地，累計探明金儲量超過4 500 t[1-2]。地球化學勘查是尋找熱液金礦床的有效手段，但受埋深深度和地表覆蓋層的影響，能夠反映深部成礦信息的弱小地球化學異常，在膠西北金礦集區常被掩蓋，地球化學異常模式分解難度大；同時膠西北金礦床的形成，是伸展背景下拆離斷層活動、核雜巖隆升、熱液流體滲流-交代復雜耦合作用的產物[3-4]，金和其他元素在空間中的分布規律具有復雜的非線性特征。此外，膠西北金礦床礦化作用還呈現多階段性[5]，致使出現地球化學元素在空間中的疊加套合，影響后續地球化學異常的推斷。因此，如何從復雜地球化學背景中有效地識別出元素異常，是膠西北地區勘查地球化學應用中面臨的主要問題。

傳統地球化學數據處理方法需要數據服從正態分布或者近似正態分布[6]，但膠西北地區地球化學數據自身復雜程度高，難以滿足該條件。由于復雜的地質環境、成礦機理往往會導致復雜的空間以及頻率特性，而能譜密度-面積多重分形模型[7-9]，是在綜合考慮地球化學數據的空間屬性、形態特征以及尺度特征的基礎上，通過傅里葉變換將成礦信息從空間域轉化到頻率域，探索能譜密度與面積大小之間的冪律分布關系，采用最小二乘方法擬合的方式來構建不同類型的分形濾波器，通過濾波器將地球化學元素組合圖分解為地球化學噪聲、背景、異常圖[10-12]。因此，能譜密度-面積模型能夠有效地降低變化背景的影響，識別出弱小的地球化學異常。此外，該模型采用冪律分布的策略，使其能夠全面地刻畫成礦過程的非線性[13]，從而實現分解復雜地球化學模式。因此，能譜密度-面積多重分形模型常應用于地球化學勘查中[14-15]，但目前在膠西北地區尚未開展。

筆者采用基于對數比變換的因子分析方法，對水系沉積物地球化學數據進行預處理，提取其中能夠反映金礦化富集帶的元素組合[16]，運用能譜密度-面積多重分形模型進行地球化學模式的分解與異常識別工作，結合對區內構造與礦床分布的認識，對膠西北地區的成礦遠景區進行評價。

1 研究區地質特征

膠西北礦集區位于膠東半島的西北部，大地構造位置屬于華北板塊東南緣的膠北斷隆。區內發現的金礦形成大體上受前寒武紀變質基底巖系、中生代燕山期構造—巖漿活動、NE-NNE向韌—脆性構造三大成礦因素聯合控制(圖1)[17]。

圖1 膠西北地質特征圖[22-23]Fig.1 Geological map of the northwestern Jiaodong peninsula

1.1 地層

研究區自太古宙到新生代地層均有分布，主要的地層單元為新太古代膠東巖群和下元古界荊山群。前者分布較為廣泛，主要巖性為角閃變粒巖夾磁鐵石英、斜長角閃巖、黑云變粒巖和英云閃長巖、奧長花崗巖等。后者主要分布在膠西北地區的東南部，主要巖性為大理巖、透輝巖、黑云變粒巖、麻粒巖、石墨片麻巖、石榴矽線黑云片巖等。

1.2 控礦構造

研究區構造活動復雜，構造變形相互疊加，北北東向和北東向斷裂構造是研究區的重要控礦構造，自西向東為：三山島斷裂、焦家斷裂、招平斷裂(圖1)。三山島斷裂分布于膠西北金礦集區的西段，總體走向為35°～45°，局部為70°～85°，傾向南東，呈舒緩波狀延展，控制了三山島、北部海域、倉上、新立等大型-超大型金礦；焦家斷裂分布于朱橋、新城、辛莊一帶，平面上呈“S”形，走向35°～40°間，傾角為30°～50°，控制新城、焦家、寺莊等大型-特大型金礦；招遠-平度斷裂南起上莊，經招遠轉為北東東向，延伸至龍口市顏家溝一帶，控制了玲瓏、大尹格莊、夏甸等大型-超大型金礦的產出。

1.3 巖漿巖

區內侵入巖主要為中生代花崗巖和鎂鐵質脈巖，中生代花崗巖依據侵位時間劃分為晚侏羅世玲瓏花崗巖和欒家河花崗巖、早白堊世郭家嶺二長花崗巖和艾山花崗巖。玲瓏和欒家河單元花崗巖是分布最廣泛的巖漿單元，常見前寒武變質巖捕虜體，為地殼重融型花崗巖。玲瓏花崗巖常具中細粒結構，片麻狀構造，主要的礦物為石英、黑云母、斜長石、鉀長石，副礦物有鋯石、榍石、磷灰石、褐簾石、磁鐵礦等。區域斷裂主要沿玲瓏花崗巖與前寒武變質巖地層發育，多數金礦體發育在玲瓏花崗巖內。

郭家嶺花崗巖主要分布在焦家礦田北東部，發育巨大的鉀長石斑晶，主要成分為鉀長石、石英、黑云母和角閃石，副礦物為鋯石、磁鐵礦、榍石、磷灰石、綠簾石等。被認為是下地殼重融和幔源巖漿混合作用的產物。區內艾山期巖體發育較少。在巖漿侵位后，膠西北地區經歷了持續的整體抬升過程，發生了持續剝蝕作用，未經歷明顯沉積，目前區內整體的剝蝕深度在 5 km左右[18-20]。

2 數據與方法

2.1 實驗數據

本研究使用的數據集為 1：200 000比例采集的水系沉積物地球化學數據[21]、膠西北礦床礦點分布數據、膠西北斷裂數據。

水系沉積物采樣密度為1 km2采集一件樣品，通常采集的重量大于120 g，然后將每4 km2組合成一件樣品，進行測試分析。主量元素和微量元素以X射線熒光光譜(XRF)和等離子體質譜(ICP-MS)分析為主結合其他方法得到，重復的樣品數據分析顯示誤差不超過5%，詳細的元素分析方法和檢出見表1。本次研究獲取的地球化學數據包含Au、Ag、Cu、Sb、Co、Ni、As等在內的39種元素(表2)。

表1 元素測試分析方法[21]Tab.1 Element test and analysis method

表2 膠西北地區地球化學元素數據特征統計Tab.2 Statistics of geochemical data from the northwestern Jiaodong peninsula

膠西北礦床(礦點)分布數據包含了礦床礦點的名稱、大小、經緯度坐標、金品位平均值(Au)等相關信息，為了便于后續的結果分析，將已知礦床礦點數據中的經緯度坐標，通過逆地理編碼轉化為具體的地理位置坐標，礦床礦點在地質圖上的分布是結果解釋的事實依據。礦床礦點分布數據與地球化學數據處理結果的耦合關系圖，為進一步識別地球化學異常提供了重要的數據支撐。

膠西北斷裂數據是將山東省萊州-招遠地區區域地質圖進行配準、矢量化得到的，并且將斷裂數據的空間坐標系進行轉化，加入最終的多重分形結果圖中，斷裂帶數據為結果的解釋提供了地質依據。

2.2 能譜密度-面積模型(S-A)

復雜的地球化學模式可視為是由多種不同頻率的信號共同作用的，基于此，成秋明等[8-9]基于頻率域建立了“能譜密度-面積”(S-A)多重分形模型：

A(S>e)∝e-β

(1)

式中：A(S>e)為能譜密度大于閾值e的面積大?。弧貫檎?；β為在研究區域A(S>e)內的能譜密度的變化程度。在A(S>e)和e的雙對數結果中，A(S>e)和e呈現分段式的線性關系并且可通過最小二乘法進行擬合，β可通過不同的線段的擬合斜率值進行估算，不同線段的交點可作為分解地球化學背景和地球化學異常的閾值。

2.3 對數比變換算法原理

地球化學數據是典型的成分數據[24-25]，具有“閉合效應”。因此，在數據處理之前需要對原始地球化學數據進行變換，否則可能會導致錯誤的分析結果或者得到元素之間的偽相關關系，這會給后續研究結果的準確性帶來很大的干擾，目前，對于原始地球化學數據進行預處理變換，有效的方法是對數比變換(ILR)[26-28]，其算法實質是根據成分分量比值的對數值通常符合正態分布的特點，利用標準正交基對成分數據進行投影變換，實現原理如下。

為了便于對數比變換算法原理的介紹，將原始地球化學數據空間抽象成一個二維數組(m為總的地球化學樣品個數，n為總的地球化學元素個數，i=1,…,n-1)。

(2)

3 結果與討論

3.1 基于ILR的因子分析處理

考慮到地球化學數據是成分數據，在進一步探討礦區的地球化學異常之前，對膠西北礦集區地球化學數據進行基于ILR變換的因子分析，通過研究區因子分析特征值和累計方差貢獻率表(表3)，以解釋總方差70%以及特征值大于“1”為基準選取八個主因子。

表3 地球化學數據因子分析結果Tab.3 Factor analysis results of geochemical data

為了深入探討ILR在解決地球化學閉合效應上的有效性，以第二主因子為例，第二主因子是研究區主要研究元素Au的主要荷載主因子，對于研究區的研究意義重大，通過對比采用ILR變換進行數據預處理后因子分析的結果(表4)與采用對數變換(Ln)預處理后因子分析的結果(表5)可知，使用兩種不同的方法,對膠西北礦集區地球化學數據進行預處理后得到的第二主因子元素組合存在較明顯的差異，表4分析得出的第二主因子為：Au、Ag、Cu、Pb、Zn、Hg、B、Ba、Mo等，而表5分析得出的第二主因子為：Au、Ag、Cu、Pb、Zn、P、Hg、Bi、Mo。在近礦特征元素：Au、Ag、Cu、Pb、Zn等元素上得到的元素組合基本上是相同的，但是通過ILR方法處理后Au元素的因子荷載系數是第二主因子中的最大值，更加突出了Au元素在第二主因子中主導性，更加符合膠西北實際地質情況，而與之相反通過對數變換(Ln)得到的則是Ag元素；相較于采用對數變換(Ln)，ILR能有效提升其他近礦特征元素在第二主因子中的因子荷載系數，進一步提升預測結果的準確性。同時，在前緣暈特征指示元素組合： As、Sb、Hg、B上也存在明顯的差別，通過ILR數據處理的As、Hg、B元素異常明顯，而通過對數變換處理(Ln)的結果中As、B的元素異常被掩蓋，同時ILR數據處理后對Hg、Sb元素的荷載系數值有所提升。結合膠西北礦集區的地球化學特征分析，在成礦帶具有Au異常的基礎上，如果還存在著As、Sb、Hg、B等前緣暈元素異常，則指示該區域深部可能有盲礦存在，所以元素組合的選擇對下一步的靶區圈定，以及靶區的成礦概率的預測上面都起著重要的作用，因此，在對地球化學數據進行因子分析之前進行ILR對數處理是必要的，它可能間接地決定了因子分析結果的準確性。

表4 地球化學數據旋轉成分矩陣(ILR變換)Tab.4 Rotated component matrix of geochemical data (ILR transform)

表5 地球化學數據旋轉成分矩陣(Ln變換)Tab.5 Rotated component matrix of geochemical data (Ln transform)

為了突出研究的重點，筆者在對研究區旋轉成分矩陣(表4)進行分析的基礎上，以前兩個主因子的元素組合為例進行討論，第一主因子：Cr、Cu、Ti、V、Co、Fe2O3，占總方差的25.408%；第二主因子：Au、Ag、Cu、Pb、Zn、Sb、Hg、B、Mo，占總方差的18.834%。由膠西北地球化學特征可知，Cr、Cu、Co、V等元素的富集與基性巖、超基性巖有關，第二主因子中Pb、Zn、W可能指示了多金屬硫化物階段的疊加，第二主因子作為研究區重要研究元素Au的荷載主因子，是膠西北金礦集區礦床預測分析的重要要素，這一主因子可以看作是膠西北金礦集區最佳指示元素組合。

3.2 多重分形模型異常識別

3.2.1 S-A多重分形建模

本研究將因子分析得到的第二主因子得分，采用距離平方反比插值得到地球化學元素組合插值圖，為了進一步探討膠西北礦集區地球化學異常，對得到的元素組合插值圖進行了S-A多重分形模型分析，將其映射到頻率域，繪制出能譜密度與面積的散點圖。在此基礎上，通過最小二乘法進行擬合，得到能譜密度和面積的雙對數擬合圖(圖2)。圖2中使用三條直線來進行擬合，代表了不同分形維數，得到第二主因子的分形方程為：

圖2 S-A模型分析最小二乘擬合圖Fig.2 Least square fitting of the S-A model analysis

S(e)=4.867e-0.3510

(3)

S(e)=5.787e-0.6403.451

(4)

S(e)=8.275e-1.0536.714

(5)

3.2.2 S-A多重分形綜合異常識別

通過圖2構建的能譜密度-面積多重分形方程以能譜密度對數值：3.451、6.714、7.907為分割閾值，構建了三個能譜密度-面積模型濾波器:噪聲、異常、背景，基于背景、異常濾波器識別出研究區地球化學背景圖、地球化學異常圖(圖3、圖4)。

圖3 S-A模型提取的地球化學背景分布圖Fig.3 Geochemical background distribution map extracted by S-A model

圖4 S-A模型提取的地球化學異常分布圖Fig.4 Geochemical anomaly distribution map extracted by S-A model

由圖2構建的分形方程可知，在不同的濾波器對應的能譜密度區間內，能譜密度與面積都成冪律分布，說明這三個區間的能譜密度分布都具有自相似性，并且各個能譜密度區間的擬合優度也有所差別，如地球化學噪聲區(圖2中噪聲區域)擬合優度R2為0.971，說明S-A模型對噪聲的擬合效果明顯，為了防止對地球化學背景與地球化學異常的分解產生干擾，需要將擬合的噪聲剔除；地球化學異常區(圖2中異常區域)，最小二乘法擬合優度R2為0.999，體現這一區間的自相似性強，說明采用S-A模型在膠西北金礦集區地球化學異常的分解與識別效果顯著；地球化學背景區(圖2中背景區域)擬合優度R2為0.955，是三個區間中的擬合最低值，可能受膠西北地區地球化學背景值跳越變換的影響。同時，由分形方程可知，冪指數大小在三個能譜密度區間存在顯著的差異，分別為-0.351、-0.640、-1.053，可以通過這種差異對地球化學場進行分解。

由圖3可知，膠西北金礦集區的主要的背景值區，從左下角到右上角大致分為東宋鎮-平度、萊州-云山鎮-萊西、龍口-焦家-招遠，三個背景值區的背景值有跳躍變化的現象，由圖1可知，這三個區域的巖體類型存在變換，斷裂帶分布差異明顯，背景值的變換可能受此影響。其中又以龍口-焦家-招遠背景值最高，萊州-云山鎮-萊西背景值最低，如果在整個膠西北金礦集區進行地球化學異常的識別的時候，采用龍口-焦家-招遠背景值區作為整個區域的背景值，會造成東宋鎮-平度、萊州-云山鎮-萊西這兩部分區域的異常，被高背景掩蓋而造成地球化學異常識別不完全，影響下一步的結果分析。由圖4可知，在龍口-焦家-招遠區中進一步縮小了地球化學異常的范圍，而且礦床礦點基本都分布在該區的中心部分，在萊州-云山鎮-萊西區消除了背景差的影響，并且增強了這一區域的地球化學異常強度，讓原本被掩蓋的異常變得清晰，和該區礦床礦點的分布耦合度高，符合膠西北金礦集區的地質事實，從而驗證了圖2建立的S-A多重分形模型的準確性。

從圖4可知，在消除背景差異影響的基礎上，地球化學異常高值區基本上分布在焦家主斷裂以及一系列次級斷裂周圍，且與斷裂的走向延展趨勢一致。此外，膠西北地區金礦床主要發育在區域斷裂面下盤的黃鐵絹英巖化蝕變帶內[29-30]。

筆者通過S-A模型識別出的地球化學異常高值區域整體呈現出以斷裂面為界，向斷層傾向方向延伸的趨勢(招平斷裂東傾、焦家斷裂西傾)，這種分布規律與已知的礦化富集帶的套合性較好，暗示了區域斷裂或構造蝕變帶是膠西北金礦床的主要控礦地質體。

進一步的統計顯示，90%的膠西北金礦集區的礦床礦點落于本文提取的地球化學異常高值區(圖4)，說明研究區S-A模型識別的地球化學異常與礦體賦存位置對應性較高，可通過地球化學異常進行成礦潛力區預測。

此外，紗嶺礦區深部標高范圍在-904 m～-2 030 m內取得重大找礦突破，探獲金資源量389 t，紗嶺礦區是焦家金礦田的重要組成部分，是目前國內平均深度最大的金礦床[31-32]。由圖4中可知，此次紗嶺礦區的深部找礦突破位置正好位于焦家斷裂帶周圍的地球化學異常高值區內。

3.3 成礦潛力區評價

筆者通過觀察S-A多重分形模型識別的地球化學異常，對比已知礦化區異常分布特征，發現膠西北地區仍存在多處成礦潛力區(圖4)。

I號靶區位于玲瓏礦田NW方向，在異常圖中具有最高的異常值。重要的是，該區位于招平斷裂帶NE方向的“拐彎”處，這種部位在膠西北地區普遍發育較大規模的蝕變和礦化，如焦家成礦帶的上莊礦區、新城礦區。

II號靶區位于焦家成礦帶NE延伸方向，該區同樣呈現出高的異常值，異常中心位于焦家主裂面的NW方向，這種規律與現有焦家成礦帶的礦體分布規律相一致，暗示了該區存在隱伏礦體。

III號靶區位于招平斷裂帶中段，本次研究結果中，招平斷裂帶中段和南部異常與焦家礦田、玲瓏礦田相比，相對較弱，但整體呈現出沿斷裂帶分布的趨勢，且與已知的礦床對應較好。靶區位置異常與南部礦床出露位置異常相似，且具東傾特征，因此可能具有找礦潛力。

IV號靶區位于招平斷裂帶北部，具有較好的異常連續性和較高的異常值。兩個異常中心分別位于兩條NE方向的次級斷裂中，可能代表了次級斷裂控制的礦化部位。

4 結論

筆者以膠西北金礦集區為研究對象，利用因子分析確定成礦作用密切相關的元素組合，運用多重分形模型進行地球化學異常的分解與提取，圈定了膠西北地區成礦遠景區。主要取得以下結論：

1)相較于對數變換(Ln)數據預處理方法，采用對數比變換(ILR)處理效果更好，能有效提高因子分析結果的準確性。

2)S-A模型能有效識別膠西北金礦集區不均一的地球化學背景，在提取弱小地球化學異常上效果顯著，本文超過92%已知礦床(礦點)落在提取的地球化學異常高值區，說明能譜密度-面積多重分形方法在膠西北金礦集區應用的可靠性，能夠為膠西北金礦集區下一步的找礦工作提供重要依據。

3)膠西北金礦集區地球化學異常分布，受到研究區礦床(礦點)和控礦構造的共同影響。

4)招平斷裂段中段、北段以及焦家斷裂帶的北東延展方向顯示出較高的地球化學異常值，具有較好的找礦前景。