圣安陳，楊悅，李亞彬，鄒宗霖，翁愛華

吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026

0 引言

花牛山鉛鋅礦床位于甘、新交界部位的北山地區(qū)，是該區(qū)一處中型的多金屬礦床。該區(qū)礦產(chǎn)資源豐富，已發(fā)現(xiàn)有鉛、鋅、金、銀和鎢等多處金屬礦床。礦田大地構(gòu)造位置位于北山塔里木板塊敦煌地塊北緣活動(dòng)帶北緣，花牛山--黑山--雙鷹山早古生代裂谷裂陷帶中[1]。

礦區(qū)北鄰方山口--營(yíng)毛沱--牛圈子早古生代裂陷帶，南接紅十井--柳園--大奇山晚古生代陸內(nèi)裂谷帶[1]。礦區(qū)內(nèi)褶皺和斷裂發(fā)育。褶皺常伴隨近EW向斷裂發(fā)育，并控制著礦區(qū)內(nèi)主要礦體的分布。礦區(qū)內(nèi)巖漿活動(dòng)頻繁且強(qiáng)烈，分布較為廣泛，從侵入巖到噴發(fā)巖，從超基性巖到酸性巖均有出露，巖漿活動(dòng)明顯受EW向與NW向構(gòu)造帶控制[2]。

甘肅省地質(zhì)局物探隊(duì)、西安地質(zhì)礦產(chǎn)研究院及廊坊物化探所先后在花牛山鉛鋅礦區(qū)開展了勘查研究工作，發(fā)現(xiàn)礦體賦存有利部位是碳酸鹽巖與碎屑巖接觸帶。相關(guān)文獻(xiàn)表明，花牛山鉛鋅礦的有利成礦地段為深大斷裂帶邊緣，成礦主要圍巖為碳酸鹽巖[3]。楊建國(guó)認(rèn)為花牛山鉛鋅礦無論是以碎屑巖和碳酸鹽巖為主的容礦區(qū)，還是以中基性火山夾碎屑巖和碳酸鹽巖為主的容礦區(qū)，均賦存于震旦系上統(tǒng)三巖組第二巖性段千枚巖、千枚狀板巖和第三巖性段大理巖夾千枚巖過渡部位，礦體主要產(chǎn)于大理巖與千枚巖接觸帶附近及其附近層間破碎帶[4]。

可控源音頻大地電磁法能克服天然場(chǎng)源信號(hào)的隨機(jī)性和信號(hào)弱的缺點(diǎn)，已成為解決地球深部資源勘探問題的一種有效手段[5]?；ㄅＩ姐U鋅礦找礦模型過去一直從地面地質(zhì)和地球化學(xué)方向出發(fā)，較少結(jié)合物探方法。為進(jìn)一步研究鉛鋅礦區(qū)的三維電性結(jié)構(gòu)，劃分成礦有利區(qū)域，廊坊物化探所在花牛山地區(qū)開展了以可控源音頻大地電磁法(CSAMT)為主的地球物理勘探工作。近年來，可控源方法的一維、二維反演技術(shù)已日趨成熟。但考慮到地下電性結(jié)構(gòu)的復(fù)雜性，一維、二維反演可靠性較差，采用三維反演將是一個(gè)更有效的手段。目前三維反演方法主要有Occam反演法[6]、擬牛頓法[7]、非線性共軛梯度法[8--9]和快速松弛反演法[10]等。筆者采用的有限內(nèi)存擬牛頓法[11]，是由Broyden[12--14]等在1970年提出并發(fā)展起來的，是目前三維電磁法反演中最有效的擬牛頓方法之一，由于該方法不需要存儲(chǔ)大型海森矩陣，因此克服了牛頓法計(jì)算量大、產(chǎn)生非下降方向的缺點(diǎn)。2006年Avdeeva[15--16]首次將該方法應(yīng)用到大地電磁(MT)一維反演中，并進(jìn)一步推廣到三維反演，非常適合求解大規(guī)模的最優(yōu)化問題。筆者利用該反演方法獲得花牛山鉛鋅礦測(cè)區(qū)的三維電性模型和界面形態(tài)，并結(jié)合前人的激電法測(cè)量結(jié)果，推測(cè)出可能的鉛鋅成礦區(qū)。

1 礦區(qū)地質(zhì)及地球物理特征

花牛山鉛鋅礦田由4個(gè)礦區(qū)組成，出露地層主要為震旦系洗腸井群二、三、四巖組，礦體主要賦存于震旦系洗腸井群碎屑巖--碳酸鹽巖巖系中，屬火山沉積碎屑巖--碳酸鹽巖系有關(guān)的噴流沉積型鉛鋅成礦系列[4]。成礦有利地層位于千枚巖與大理巖的接觸帶。礦床的形成受震旦系地層、巖漿熱液和構(gòu)造作用影響。二巖組(Zxb)為一套淺變質(zhì)的淺海相泥質(zhì)巖建造，在其內(nèi)分布有大量酸性斑巖和石英脈，整體巖性較為單一；三巖組(Zxc)主要為一套淺海相碳酸鹽巖與泥質(zhì)巖建造，巖性主要為含粒狀石英的大理巖夾絹云千枚巖及粉砂質(zhì)板巖，鉛鋅礦體產(chǎn)于該巖組中；四巖組(Zxd)主要為板巖和角巖等[17]。礦田南部有EW向的花西灘--花牛山壓(扭)性斷裂(F1)，北部有EW向的五井河隱伏斷裂(F3)，F(xiàn)1斷裂規(guī)模較大，控制著巖體的侵位，是主要的控礦構(gòu)造[18]。

測(cè)區(qū)位于一礦區(qū)，北部為震旦系洗腸井群三巖組大理巖夾千枚巖，南部為震旦系洗腸井群四巖組千枚巖和角巖。兩者的接觸帶為二礦帶賦存位置，位于三巖組第三巖性段大理巖夾千枚與三巖組第二巖性段千枚巖過度部位。東部、西部均有印支期花崗巖侵入。測(cè)區(qū)南端為三礦帶，測(cè)點(diǎn)未覆蓋該礦帶。測(cè)區(qū)內(nèi)有NW向的花黑灘--雙峰山斷裂(F2)與二礦帶構(gòu)造走向交叉，是后期次級(jí)斷裂帶，與礦體產(chǎn)狀和走向有一定關(guān)系(圖1)。

物性差異是地球物理工作的前提。測(cè)區(qū)內(nèi)主要的巖石為千枚巖、含炭千枚巖、花崗閃長(zhǎng)巖、大理巖和鉛鋅磁黃鐵礦。5種主要巖石電阻率的差異較大，其中花崗閃長(zhǎng)巖的電阻率最高，千枚巖、大理巖的電阻率相對(duì)較高，含炭千枚巖、鉛鋅磁黃鐵礦的電阻率最低(表1)，這為可控源音頻大地電磁反演異常的解釋提供了基礎(chǔ)。此外，測(cè)區(qū)內(nèi)鉛鋅磁黃鐵礦礦石、含炭千枚巖具有低阻高極化的特點(diǎn)。因此，含炭千枚巖是電法找礦的重要地質(zhì)干擾因素。

表1 測(cè)區(qū)主要巖石電性特征[19]Table 1 Electrical characteristics of main rocks in survey area

2 數(shù)據(jù)及反演

可控源音頻大地電磁法(CSAMT)是在大地電磁(MT)的基礎(chǔ)上發(fā)展起來的一種人工源電磁測(cè)深方法，主要適用于油氣、地?zé)岬荣Y源勘探。目前，由于其勘探深度較大，逐漸成為深部找礦的重要地球物理手段。

1.中奧陶統(tǒng)花牛山群；2.震旦系洗腸井群四巖組；3.震旦系洗腸井群三巖組；4.震旦系洗腸井群二巖組；5.印支期正長(zhǎng)花崗巖；6.印支期二長(zhǎng)花崗巖；7.華力西期花崗閃長(zhǎng)巖；8.橄欖巖；9.斷裂；10.測(cè)區(qū)；11.CSAMT測(cè)線；F1.花西灘—花牛山斷裂；F2.花黑灘—雙峰山斷裂；F3.五井河斷裂。圖1 花牛山礦區(qū)地質(zhì)概況圖Fig.1 Geological map of Huaniushan ore deposit

2.1 數(shù)據(jù)采集

本次研究采用的可控源音頻大地電磁數(shù)據(jù)是廊坊物化探所獨(dú)立研制的多功能電磁工作站采集的野外實(shí)測(cè)數(shù)據(jù)。測(cè)區(qū)面積約3 km2，觀測(cè)方式采用赤道偶極標(biāo)量方式，觀測(cè)垂直測(cè)線方向的磁場(chǎng)分量Hy與沿測(cè)線方向的電場(chǎng)分量Ex，發(fā)射極距為1.5 km，接收極距為50 m，最大供電電流為 20 A，最小收發(fā)距為 8 km，共布設(shè)13條測(cè)線，分別命名為L(zhǎng)1、L2到L13，測(cè)線L1距離發(fā)射源最遠(yuǎn)，測(cè)線L13距離發(fā)射源最近，線距為 200 m，點(diǎn)距為50 m，每條測(cè)線24個(gè)測(cè)點(diǎn)，一共312個(gè)測(cè)點(diǎn)。將收集到的數(shù)據(jù)進(jìn)行處理，主要掃描時(shí)間序列，根據(jù)電磁干擾程度選擇合適的頻率段參與電場(chǎng)分量Ex估計(jì)。

數(shù)據(jù)質(zhì)量對(duì)反演結(jié)果至關(guān)重要。圖2為76號(hào)點(diǎn)、124號(hào)點(diǎn)、172號(hào)點(diǎn)和220號(hào)點(diǎn)本次反演的數(shù)據(jù)電場(chǎng)分量Ex。原始數(shù)據(jù)質(zhì)量均較好，高頻部分出現(xiàn)的跳點(diǎn)、飛點(diǎn)不明顯，說明數(shù)據(jù)受到的噪聲干擾較小，反演可以得到地下可靠的信息。

2.2 有限內(nèi)存擬牛頓反演方法

根據(jù)Tikhonov正則化理論，地球物理反演問題的目標(biāo)函數(shù)為：

Φ=Φd+λΦm

(1)

式中：Φd為觀測(cè)數(shù)據(jù)和理論模型響應(yīng)的擬合差；Φm為模型范數(shù)；λ為正則化子。

本文的正則化目標(biāo)函數(shù)定義為：

(2)

式中：m是模型參數(shù)向量；d是數(shù)據(jù)向量；Cd是數(shù)據(jù)方差矩陣；Cm是模型方差函數(shù)；F是正演算子。

L--BFGS方法可以直接求取目標(biāo)函數(shù)的最優(yōu)化問題[20]。將模型的搜索方向定義為：

(3)

利用充分下降條件和曲率條件來約束搜索步長(zhǎng)αk[21]。由于求解式(3)需要很長(zhǎng)的耗時(shí)和很大的存儲(chǔ)空間，L--BFGS在牛頓法的基礎(chǔ)上，將Hessian矩陣的逆近似改進(jìn)為：

(4)

第一步：令k=1，給定初始模型mk，初始正定矩陣Hk，允許誤差范圍ε>0；

第二步：計(jì)算梯度rk=-▽U(λ,mk)，當(dāng)‖rk‖≤ε，輸出最終解mk；

第三步：否則,利用(3)式進(jìn)行步長(zhǎng)αk的線性搜索。當(dāng)步長(zhǎng)αk小于闕值，減小正則化因子λ；如果λ小于闕值，輸出最終解mk，反演終止；

第四步：更新模型mk+1=mk+αkPk，若‖rk+1‖≤ε，得到最終反演結(jié)果mk+1；

第五步：否則，利用(4)式更新正定矩陣Hk+1；

第六步：令k=k+1，返回第二步。

2.3 理論模型模擬

為研究有限內(nèi)存擬牛頓法三維反演的可靠性，設(shè)計(jì)了理論模型(圖3)。模型為均勻半空間，背景電阻率為100 Ω·m，交錯(cuò)放置8個(gè)300 m×300 m×300 m的立方異常體，紅色代表高阻體為500 Ω·m，藍(lán)色代表低阻體為10 Ω·m。發(fā)射源AB使用1 000 m的長(zhǎng)導(dǎo)線源，最小發(fā)射極距為7 500 m，發(fā)射頻率為2-2、2-1、…、213Hz，共計(jì)16個(gè)頻點(diǎn)。測(cè)區(qū)(虛線框)布設(shè)11條測(cè)線，線距設(shè)置為100 m，點(diǎn)距設(shè)置為50 m，每條測(cè)線21個(gè)測(cè)點(diǎn)，共計(jì)231個(gè)測(cè)點(diǎn)，觀測(cè)參數(shù)為電場(chǎng)Ex復(fù)分量。

圖3 理論模型及工作裝置示意圖Fig.3 Theoretica model and working device diagram

理論模型數(shù)值模擬時(shí)，測(cè)區(qū)范圍內(nèi)將網(wǎng)格剖分為20×20×25，網(wǎng)格長(zhǎng)度均為50 m，由于模型的范圍應(yīng)略大于測(cè)區(qū)范圍，因此模型以2.0倍的系數(shù)各擴(kuò)邊5個(gè)網(wǎng)格，最終生成的網(wǎng)格為30×30×30，共計(jì)27 000個(gè)。反演部分參數(shù)的設(shè)置，背景電阻率為100 Ω·m，正則化因子初始值為1，疊加5%的噪聲到理論數(shù)據(jù)。

圖4三維反演結(jié)果深度為200 m和900 m的水平切片。圖中可以看出淺部異常體，高、低阻均能清晰的刻畫出異常體的區(qū)域，反演的電阻率大小與理論模型異常體電阻率比較接近。深部異常體的低阻的輪廓有所展現(xiàn)，高阻不明顯，與理論模型電阻率異常體電阻率相差很大。說明有限內(nèi)存擬牛頓法淺部反演結(jié)果是比較可靠的。

a.深度200 m切片；b.深度900 m切片。圖4 Ex分量復(fù)數(shù)據(jù)三維反演結(jié)果Fig.4 3D model inversed results from Ex component

2.4 實(shí)測(cè)數(shù)據(jù)反演

雖然L--BFGS反演方法相比于非線性共軛梯度法(NLCG)在反演速度和反演結(jié)果準(zhǔn)確性方面有一定的優(yōu)勢(shì)[11]，但在實(shí)際應(yīng)用中，三維模型模擬的計(jì)算量很大，計(jì)算時(shí)間很長(zhǎng)。為減少模型網(wǎng)格數(shù)量，提高反演速度，對(duì)測(cè)線進(jìn)行了16°的逆時(shí)針旋轉(zhuǎn)，建立x軸指向正北，y軸指向正東，z軸指向地心的空間坐標(biāo)系，其中坐標(biāo)原點(diǎn)為發(fā)射源的中點(diǎn)。實(shí)測(cè)數(shù)據(jù)的頻率范圍為0.28～8 000 Hz,共計(jì)45個(gè)頻點(diǎn)，本次反演使用的頻率為0.2 Hz、55 Hz、22 Hz、640 Hz、1 280 Hz 5個(gè)頻點(diǎn)。反演模型沿x方向網(wǎng)格30個(gè)，網(wǎng)距40 m，兩側(cè)以2.0倍系數(shù)各擴(kuò)邊5個(gè)網(wǎng)格。y方向網(wǎng)格30個(gè)，網(wǎng)距100 m，兩側(cè)以2.0倍系數(shù)各擴(kuò)邊5個(gè)網(wǎng)格。z方向首層厚度20 m，其下各層厚度按系數(shù)1.10遞增，邊界以1.5倍系數(shù)擴(kuò)邊5個(gè)網(wǎng)格，最終生成的反演網(wǎng)格為40(東西)×40(南北)×30(垂向，不包含空氣層)，反演初始模型采用2 000 Ω·m的均勻半空間模型，擴(kuò)邊網(wǎng)格的電阻率同樣設(shè)定為2 000 Ω·m，誤差門限設(shè)置為5%的噪聲值。

圖5給出了迭代反演過程中參數(shù)λ和參數(shù)rms的變化情況，其中λ和rms在收斂時(shí)表征觀察數(shù)據(jù)的擬合程度。從圖中可以看出經(jīng)過651次的反演迭代，λ在收斂時(shí)<10-5；擬合差降到了6.86。雖然該擬合差>1，但初始模型擬合差為700，通過反演迭代擬合差已經(jīng)下降了近2個(gè)數(shù)量級(jí)，說明了反演結(jié)果的可靠性。

圖5 反演參數(shù)λ(a)和rms(b)曲線圖Fig.5 Curves of inverse parameter λ(a)and rms(b)

3 結(jié)果與討論

3.1 反演結(jié)果

為了反演結(jié)果成圖效果更好，將電阻率取對(duì)數(shù)后進(jìn)行三維成圖，圖6給出了地下1 000 m±的三維立體反演結(jié)果顯示。從圖中可以看出，地表3個(gè)條帶狀低阻異常帶J-D1、J-D2和J-D3，其中J-D1范圍最大，位于測(cè)區(qū)西北角，下部可能賦有礦體；J-D2、J-D3圍巖的電阻率較高，根據(jù)測(cè)區(qū)成礦規(guī)律，礦體產(chǎn)于大理巖與千枚巖接觸帶部位，與J-D2、J-D3邊緣呈現(xiàn)的中、低阻過渡帶較吻合，下部也可能賦有礦體。同時(shí)，李建華等[2]在該測(cè)區(qū)開展了TDIP法測(cè)量，其視極化率剖面平面圖和視電阻率剖面平面圖的結(jié)果圈定的異常區(qū)域與筆者的三維反演結(jié)果地表圈定的異常帶吻合，證實(shí)了采用有限內(nèi)存擬牛頓法三維反演的淺部信息是可靠的。

圖6 可控源音頻大地電磁測(cè)深反演結(jié)果三維圖Fig.6 3D inversion results of CSAMT

由于大理巖和千枚巖兩種巖石巖性的極化率相差不大，而電阻率相差較大，因此僅靠TDIP測(cè)量達(dá)不到找礦的可靠效果，仍需結(jié)合可控源音頻大地電磁測(cè)深結(jié)果來綜合探討礦體賦存的有利部位。

3.2 控礦斷裂分析

斷裂是花牛山鉛鋅礦區(qū)主要控礦因素,斷裂帶的發(fā)育部位，常伴隨礦體的產(chǎn)出。為了分析地下斷裂F2(花黑灘--雙峰山斷裂)走向，將電阻率模型沿測(cè)線方向切片。根據(jù)測(cè)區(qū)地質(zhì)資料，斷裂F2位于測(cè)區(qū)西部，只選取了測(cè)線L2--L8的切片圖。圖7中L2(y=10 200 m)、L3(y=10 000 m)、L4(y=9 800 m)、L6(y=9 400 m)、L7(y=9 200 m)、L8(y=9 000 m)測(cè)線切片平面圖可以明顯看出地下地層為二層電性結(jié)構(gòu)，淺部地層電阻率較低，深部地層電阻率較高，說明測(cè)區(qū)的鉛鋅礦帶埋深較淺。圖中斷裂F2均有所呈現(xiàn)，L2斷面圖中斷裂F2地表露出位置為x=200 m處，L4斷面圖中斷裂F2地表露出位置為x=0 m處，L6斷面圖中斷裂F2地表露出位置為x=-300 m處，L8斷面圖中斷裂F2地表露出位置為x=-400 m處，綜合判斷分析，斷裂F2走向?yàn)楸蔽鞣较?，傾角45°，向南西傾斜，延伸較淺(<500 m)，為次級(jí)斷裂帶。

圖7 三維反演電阻率模型沿測(cè)線的垂直切片圖Fig.7 Vertical slice of 3D inversion resistivity model along survey line

3.3 成礦有利區(qū)圈定

由于斷裂F2位于測(cè)區(qū)西部，礦體的產(chǎn)出常在斷裂帶邊緣，因此取出L1、L8、L13測(cè)線切片三維圖來劃分成礦區(qū)域。圖8中L1測(cè)線的北端有明顯的低阻異常，異常區(qū)延伸至地下400 m±，根據(jù)測(cè)區(qū)地質(zhì)資料，該區(qū)域地層為低電阻率的含炭千枚巖地層，是主要的找礦干擾地層。L8測(cè)線中部的地下低阻異常區(qū)明顯，異常帶南北方向長(zhǎng)約500 m，東西方向長(zhǎng)約800 m，延伸至地下300 m±，與圖6地表的低阻異常帶J-D2，J-D3吻合。測(cè)區(qū)地質(zhì)資料顯示該區(qū)域地層為中電阻率大理巖夾千枚巖地層，說明淺部異?？赡苁堑妥桡U鋅礦引起的，為找礦有利區(qū)域。L13測(cè)線下方地層均為高阻地層，非找礦有利區(qū)域。

圖8 三維反演模型電阻率切片圖Fig.8 Resistivity slice of 3D inversion model

為了更直觀地展示鉛鋅礦帶的成礦區(qū)域和成礦地層接觸帶，將三維反演結(jié)果進(jìn)行電阻率劃分，圖9顯示了測(cè)區(qū)不同巖性的三維立體圖像。結(jié)合圖8找礦的預(yù)測(cè)結(jié)果和測(cè)區(qū)地面地質(zhì)資料，圖中顯示的黃色區(qū)域?yàn)橹须娮杪蕩r層，可能表征千枚巖，綠色區(qū)域?yàn)橹须娮杪蕩r層，可能表征大理巖，測(cè)區(qū)左上方有明顯的低阻層覆蓋，推測(cè)為含炭千枚巖地層，與圖6的三維反演結(jié)果吻合。同時(shí)，可以推測(cè)出測(cè)區(qū)的西部地區(qū)地層淺部表征低阻含炭千枚巖，深部表征中阻震旦系洗腸井群大理巖；中部地區(qū)地層表征中阻震旦系洗腸井群大理巖；東部地區(qū)地層表征中阻震旦系洗腸井群千枚巖，淺部可能有高阻印支期花崗巖侵入。圖中白色區(qū)域代表兩種巖石的巖性接觸帶，紫色的低阻異常位于大理巖與千枚巖接觸帶附近，L8測(cè)線之下，地下約200 m±。同時(shí)李健華等[2]的結(jié)果地表圈定該部位的異常也表現(xiàn)為高極化率，與圖6中J-D2條帶狀低阻異常區(qū)對(duì)應(yīng)，推測(cè)該低阻異常為鉛鋅礦帶，走向?yàn)榻鼥|西方向，與已知二礦帶對(duì)應(yīng)。斷裂F2與鉛鋅礦帶走向相近，位于礦體邊緣，證實(shí)了該斷裂是測(cè)區(qū)鉛鋅礦帶的控礦構(gòu)造。

圖9 三維地層分布和成礦有利區(qū)域電性結(jié)構(gòu)圖Fig.9 3D stratigraphic distribution and electrical structure map of favorable regional mineralization

4 結(jié)論

(1)使用L--BFGS法進(jìn)行反演得到的三維電性模型，模型清晰地展示了測(cè)區(qū)的電性特征，顯示測(cè)區(qū)西北角大范圍的低阻異常主要是由含炭千枚巖地層引起的，鉛鋅礦找礦地層電性特征標(biāo)志為低、中阻過渡帶部位。

(2)根據(jù)三維反演結(jié)果，探明了斷裂F2的情況，圈定出測(cè)區(qū)內(nèi)鉛鋅礦帶的成礦有利區(qū)域，認(rèn)為該礦體產(chǎn)于大理巖與千枚巖接觸帶附近，斷裂F2為該礦帶的控礦構(gòu)造。

(3)采用可控源方法探測(cè)鉛鋅礦是可行的，可以將其應(yīng)用到其他類似的鉛鋅礦勘探工作中。

致謝本文采用的可控源數(shù)據(jù)由科研任務(wù)(2011YQ05006010)所屬國(guó)家重大科學(xué)儀器研發(fā)計(jì)劃開發(fā)的儀器采集。資料處理和建?；谌S電磁數(shù)據(jù)處理平臺(tái)EMDesk。感謝吉林金太地球探測(cè)技術(shù)有限公司提供的三維反演計(jì)算資源。