楊 力， 周可法， 王金林

(1.中國科學(xué)院 新疆生態(tài)與地理研究所，新疆礦產(chǎn)資源研究中心，烏魯木齊 830011；2.中國科學(xué)院大學(xué)，北京 100049)

0 引言

剩磁改變了總磁化強(qiáng)度的大小和方向，導(dǎo)致磁異常的幅值和形態(tài)產(chǎn)生畸變，從而影響磁測資料的反演解釋。磁化強(qiáng)度矢量反演方法因其能同時(shí)反演得到巖石的磁化強(qiáng)度的大小和方向，消除巖石剩磁的影響，成為研究的熱點(diǎn)。目前，研究磁異常強(qiáng)度矢量反演的方法主要有兩種：①直接反演磁化強(qiáng)度矢量(如王妙月等[1-2]通過反演磁化強(qiáng)度矢量的三個(gè)分量最終獲得磁化強(qiáng)度大小和方向，李澤林等[3-4]通過磁異常模量約束反演磁化強(qiáng)度矢量，Lelilèvre等[5]利用有限差分法、劉雙[6]利用有限元方法重建地下高磁化率磁性體的三維磁化率分布并取得較好效果)；②利用地面化極異常垂直梯度和總梯度的互相關(guān)關(guān)系來估計(jì)磁化強(qiáng)度的傾角和偏角[7]。磁場強(qiáng)度矢量反演方法適用于笛卡爾直角坐標(biāo)系和球坐標(biāo)系[5]，大多數(shù)的反演方法都是基于笛卡爾直角坐標(biāo)系展開的，但是磁化強(qiáng)度矢量反演將會(huì)產(chǎn)生三倍于常規(guī)磁化率反演方法的變量。另外，磁化強(qiáng)度矢量反演方法基于正則化方程計(jì)算磁化強(qiáng)度和磁化方向的三維分布，反演結(jié)果較為平滑，地質(zhì)體界限模糊，將反演結(jié)果直接用于地質(zhì)解釋會(huì)導(dǎo)致很多信息丟失。在球坐標(biāo)系中，地質(zhì)信息的稀疏化約束和正則化方程具有較好的耦合關(guān)系，反演結(jié)果能夠更加真實(shí)的反映地質(zhì)體分布信息。

近年來，我國加緊對(duì)老礦山深部及外圍的礦產(chǎn)資源進(jìn)行勘查，促進(jìn)了高精度航磁測量儀器和觀測技術(shù)的發(fā)展。航空高精度磁測具有探測精度高、受地形影響小、快速高效的特點(diǎn)[8]，因此在球坐標(biāo)系下對(duì)高精度磁測數(shù)據(jù)的反演解釋是必要的。筆者首先推導(dǎo)由笛卡爾直角坐標(biāo)系下磁化強(qiáng)度矢量反演的公式。然后將其轉(zhuǎn)換到球坐標(biāo)系下進(jìn)行計(jì)算并建立試驗(yàn)?zāi)Ｐ?，最后將該方法?yīng)用于對(duì)東天山卡拉塔格礦區(qū)，分析處理其1∶50 000高精度航磁測量數(shù)據(jù)[9]，對(duì)反演結(jié)果進(jìn)行實(shí)測數(shù)據(jù)驗(yàn)證，且將反演結(jié)果與笛卡爾直角系下的反演結(jié)果進(jìn)行對(duì)比，并進(jìn)一步分析該礦床磁化強(qiáng)度異常的展布特征，綜合地質(zhì)和巖石磁性資料，探討了紅海礦床的找礦潛力。

1 磁化強(qiáng)度矢量反演方法

1.1 直角坐標(biāo)系下的反演方法和原理

磁感應(yīng)強(qiáng)度大小與觀測距離之間的關(guān)系定義為式(1)。

(1)

將積分方程離散化后得到，b=TM，其中：T是二階張量，M=[(MxMyMz)]T，b=[(bxbybz)]T。每一個(gè)離散單元的磁化強(qiáng)度為：

M=κH+MNRM

(2)

其中:κ是磁化率；H是磁場強(qiáng)度；MNRM是天然剩磁。磁感應(yīng)強(qiáng)度觀測值與m的關(guān)系可通過算子F定義為：

b=F(m)

(3)

其中:m是三維空間的磁化強(qiáng)度矢量模型。磁化強(qiáng)度矢量反演問題轉(zhuǎn)化為求解m的問題，反問題可以寫為：

m=F(m)-1dobs

(4)

其中:dobs是觀測到的磁感應(yīng)場強(qiáng)度大小，F(xiàn)(m)-1是算子F(m)的逆。根據(jù)Tikhonov[10]正則化方程，磁化強(qiáng)度矢量反演的反問題可以轉(zhuǎn)化為求解正則化方程的形式：

φ(m)=φd+βφm

(5)

其中:Wd是標(biāo)準(zhǔn)差的倒數(shù)矩陣；W是空間距離加權(quán)矩陣，β是權(quán)重系數(shù)。對(duì)φ(m)求關(guān)于m的偏導(dǎo)數(shù)為“0”，可得：

βWTW(m-mref)=0

(6)

其中:雅克比矩陣J是預(yù)測觀測數(shù)據(jù)對(duì)m的偏導(dǎo)數(shù)。易得式中數(shù)據(jù)擬合差項(xiàng)依賴于雅克比矩陣，通過空間距離矩陣W加入先驗(yàn)信息對(duì)模型范數(shù)項(xiàng)約束可以得到更加聚焦的解，模型范數(shù)項(xiàng)可定義為：

(7)

其中:Gx、Gy、Gz是三個(gè)方向上梯度的有限差分算子，[Ws,Wi,i=x,y,z]和[Rs,Ri,i=x,y,z]是稀疏化對(duì)角矩陣。敏感度矩陣Wr的選取能夠?qū)Ψ囱萁Y(jié)果產(chǎn)生重要的影響。建立迭代更新的敏感度權(quán)重方程[11]，可以寫為：

(8)

目前，航空磁測數(shù)據(jù)仍然以總磁場強(qiáng)度數(shù)據(jù)(TMI)為主，為了能夠利用大量磁測數(shù)據(jù)，本文選取總場強(qiáng)度數(shù)據(jù)進(jìn)行計(jì)算。總場強(qiáng)度可以定義為：

dTMI=F(κ)=PTH0κ=Fκ

dTMI∈Rn,F∈Rn×n,κ∈Rn

(9)

其中:H0是地磁場強(qiáng)度；P是將磁場轉(zhuǎn)換為總場強(qiáng)度的投影矩陣。在迪卡爾直角坐標(biāo)系中，磁感應(yīng)強(qiáng)度大小可以寫為：

b=T[(HpHsHt]κ

(10)

其中:Hp、Hs、Ht是笛卡爾右手直角坐標(biāo)系中的磁場強(qiáng)度分塊矩陣；p軸指向磁感應(yīng)場方向；κ是等效磁化率；等效磁化率κe={κpκsκt} 可以定義為：

κp=MP/H0

κs=MP/H0

κt=MP/H0

(11)

式(9)中的總場強(qiáng)度可以寫成：

dTMI=Fpκp+Fsκs+Ftκt

(12)

通過一個(gè)簡單的立方體模型來說明直角坐標(biāo)系下磁化強(qiáng)度矢量反演方法對(duì)結(jié)果的影響。立方體的體積為15 m×15 m×15 m，磁性體模型的的磁化率為0.05，磁傾角為-45°，磁偏角為215°。背景磁場的磁場強(qiáng)度設(shè)定為5 000 nT，并且磁傾角為90°，磁偏角為0°。構(gòu)造的觀測數(shù)據(jù)為總磁場強(qiáng)度數(shù)據(jù)，并且數(shù)據(jù)均帶有1 nT的高斯白噪聲用于模擬實(shí)際情況。具體模型的反演結(jié)果見圖1。

圖1 直角坐標(biāo)系下的模型水平切片和模型垂直切片F(xiàn)ig.1 Plan section through the recovered magnetization vector model using the cartesian formulation(a)模型水平切片；(b)模型垂直切片

圖2 球坐標(biāo)系下的模型水平切片和模型垂直切片F(xiàn)ig.2 Plan section through the recovered magnetization vector model using the spherical formulation(a)模型水平切片；(b)模型垂直切片

1.2 球坐標(biāo)系下的反演方法和原理

在球坐標(biāo)系下，磁場強(qiáng)度大小{α}和磁化方向{θ,φ}能夠較好地解耦。與直角坐標(biāo)系下的反演相比具有兩大優(yōu)勢，①野外磁測數(shù)據(jù)和實(shí)驗(yàn)室磁測數(shù)據(jù)都會(huì)給出球坐標(biāo)系下的參數(shù)值，將其轉(zhuǎn)換到直角坐標(biāo)系會(huì)增加運(yùn)算的繁瑣程度；②球坐標(biāo)系下能分別對(duì)磁化強(qiáng)度大小和磁化方向?qū)嵤┫∈杌s束。

對(duì)模型范數(shù)項(xiàng)求解關(guān)于{α,θ,φ}的偏微分，其形式可以寫成：

(13)

J=FS

(14)

其中:矩陣S可以寫成：

(15)

在球坐標(biāo)系下，敏感度矩陣是非線性的[5]。在使用高斯-牛頓迭代方法進(jìn)行計(jì)算時(shí)，模型范數(shù)項(xiàng)、磁化強(qiáng)度大小和磁化方向收斂速度很快。所以，敏感度矩陣也需要適應(yīng)于這種變化。式(8)中的空間加權(quán)矩陣可以定義為：

(16)

由于磁傾角的取值范圍為[-π/2,π/2]，磁偏角的取值范圍為[-π,π]。磁化強(qiáng)度α的變化范圍很大，通過選取最大值約束進(jìn)行迭代更新。權(quán)重ωθ，ωφ可以定義為：

采用與直角坐標(biāo)系下反演時(shí)同樣的模型進(jìn)行反演，圖2(b)是球坐標(biāo)系下磁場強(qiáng)度矢量反演結(jié)果。對(duì)比直角坐標(biāo)系下磁化強(qiáng)度矢量反演得到模型的水平切片和垂直切片(圖1)和球坐標(biāo)系下的反演結(jié)果(圖2)，可以看出與直角坐標(biāo)系下的反演結(jié)果相比，在水平位置以及深度數(shù)值方面，球坐標(biāo)系下的反演結(jié)果都能得到更加聚焦的解。

2 航磁異常在勘查區(qū)的應(yīng)用

2.1 礦區(qū)地質(zhì)特征

紅海-黃土坡銅鋅多金屬礦床分布于卡拉塔格隆起區(qū)西南邊緣的斜坡地帶(圖3)。紅海銅鋅多金屬礦與黃土坡多金屬礦分屬不同公司的勘探區(qū)，但在地質(zhì)上屬于同一礦床。礦區(qū)發(fā)育的地層為晚奧陶-早志留統(tǒng)荒草坡群大柳溝組，為一套活動(dòng)陸源性質(zhì)的海相火山巖-火山碎屑巖建造，自下自上可分為四個(gè)巖性段。礦區(qū)主要出露第二及第四巖性段，第二巖性段巖性主要為英安巖、英安質(zhì)凝灰?guī)r、沉凝灰?guī)r、砂礫巖，第四巖性段巖性主要為流紋巖、流紋斑巖、含礫流紋質(zhì)沉凝灰?guī)r、角礫凝灰熔巖。礦區(qū)中發(fā)育一套形成于早古生代活動(dòng)大陸邊緣弧構(gòu)造環(huán)境的英云閃長巖-花崗閃長巖侵入組合，主要出露于礦區(qū)中南部，另有少量的北西-北北西向的石英閃長巖脈、二長花崗巖脈和流紋斑巖脈產(chǎn)出。地表發(fā)育北西向-北西西向-近東西向的斷裂系統(tǒng)，它們都是繼承早期基底構(gòu)造發(fā)展而成的長期多次火山活動(dòng)斷裂，銅鋅多金屬礦就位于早期構(gòu)造控制的火山構(gòu)造盆地中[14]。

圖3 卡拉塔格地區(qū)地質(zhì)簡圖 [12] Fig.3 Geological map of the Kalatage region

圖4 磁測數(shù)據(jù)ΔT等值線圖Fig.4 ΔT anomaly contour map of aeromagnetic in Kalatage district

2.2 航磁ΔT異常特征

通過磁異常ΔT等值線圖(圖4)和地質(zhì)圖(圖3)對(duì)比分析。其中，主要的航磁異常有M1、M2、M3、M4，其分布見圖4。區(qū)內(nèi)酸性巖(如流紋巖、花崗巖等)多具中低磁性，可引起低磁異常；中性巖(閃長巖)具有較高磁性，可引起較高的磁異常特征。研究區(qū)內(nèi)沉積碎屑巖磁性較低，火山巖類隨巖石基性程度有增強(qiáng)趨勢，且磁性變化范圍大(表1)。

航磁ΔT磁異常值變化范圍大，磁異常主要沿北西西-近東西向和北西向呈帶狀分布，與區(qū)域構(gòu)造方向一致。M1、M2、M4磁異常帶比M3高，對(duì)應(yīng)于研究區(qū)中的中性-基性巖石。M3低磁異?？赡苁潜黄扑閹?、碎屑巖覆蓋的區(qū)域，由于ΔT值不能反映磁化率的三維分布結(jié)構(gòu)，因此限制了對(duì)磁異常體的進(jìn)一步分析。

表1 紅海礦區(qū)巖礦石磁化率[13] Tab.1 Susceptibility parameter of rocks and minerals in Honghai deposit

圖5 磁化率模型Fig.5 Magnetic susceptibility model through the recovered MVI model(a)直角坐標(biāo)系下反演結(jié)果;(b)球坐標(biāo)系下反演結(jié)果

3 反演結(jié)果

磁化強(qiáng)度矢量反演方法能夠得到準(zhǔn)確的磁化率三維分布。在反演中，將地下半空間剖分為30 m×30 m×10 m的單元，為了減少邊界效應(yīng)，模型的邊界范圍比研究區(qū)大，模型核部包含183×178×80個(gè)長方體單元。針對(duì)該區(qū)域?qū)?22 m地下深部進(jìn)行切片分析。在第一節(jié)中通過模型試驗(yàn)驗(yàn)證了利用球坐標(biāo)系下磁化強(qiáng)度矢量反演方法的優(yōu)勢。對(duì)紅海-黃土坡礦區(qū)的實(shí)測數(shù)據(jù)進(jìn)行反演，顯示模型核部的反演結(jié)果。

直角坐標(biāo)系下，反演得到的結(jié)果總是發(fā)散的，不能真實(shí)地反映地質(zhì)體的展布特征(圖5(a))。而在球坐標(biāo)系下，能夠得到聚焦的解(圖5(b))。通過球坐標(biāo)系下三維磁化率模型的切片圖(圖5(b))，可以劃分出F1，F(xiàn)2，F(xiàn)3，F(xiàn)4四條斷層，與圖5a相比具有更清晰的邊界特征。對(duì)比ΔT等值線圖可以發(fā)現(xiàn)F2與F3之間的高磁異常區(qū)可能是磁性礦物的聚集區(qū)。紅海-黃土坡銅多金屬礦床屬于隱伏礦床，通過磁化率反演得出深部F2、F3斷裂控制了火山裂陷沉積中心的發(fā)育，并成為紅海-黃土坡銅鋅多金屬礦礦質(zhì)的熱液通道和主要富集沉淀地。

圖6 磁化方向的計(jì)算結(jié)果Fig.6 Magnetic inclination and declination through the recovered MVI model(a)磁偏角分布圖; (b)磁傾角分布圖

另外，通過磁化強(qiáng)度矢量反演得到的磁化方向分布，能夠指示巖體和沉積物的形成時(shí)代，為礦體的圈定提供依據(jù)。從磁化方向計(jì)算結(jié)果中可以看到(圖6(a),圖6(b))，M3磁異常帶具有強(qiáng)的偶極磁化異常，表明在F2與F3斷裂中間沉積了大量的碎屑物質(zhì)，而這些碎屑物質(zhì)為VMS礦床的形成提供了良好的就位空間。一般火山碎屑沉積越厚，銅鋅礦化越發(fā)育[15]。因此研究區(qū)深部有較大的找礦-成礦潛力，具有大型礦床規(guī)模遠(yuǎn)景。對(duì)比磁化方向和磁異常ΔT等值線圖發(fā)現(xiàn)，M1和M2雖然在地表表現(xiàn)出較強(qiáng)的磁異常，但在深部表現(xiàn)出低的磁化率，表明M1和M2不具有成礦潛力。M1、M2、M3磁異常帶，具有相同的磁化方向所以很有可能都是同一時(shí)期火山活動(dòng)的產(chǎn)物，而M4磁異常區(qū)，也具有較強(qiáng)的偶極異常，同時(shí)在觀測值中呈現(xiàn)出高的磁異常，也具有一定的找礦潛力。

4 結(jié)論

1)球坐標(biāo)系下磁化強(qiáng)度矢量反演不僅比笛卡爾直角坐標(biāo)系下反演的速度更快，而且得到的三維磁化率模型比笛卡爾直角坐標(biāo)系下反演得到的結(jié)果更加聚焦。

2)礦區(qū)發(fā)育的M1、M2、M3磁異常帶與區(qū)內(nèi)斷層發(fā)育方向一致，表明斷裂控制了該區(qū)域磁異常分布。F2和F3之間的磁異常M3是可能的成礦遠(yuǎn)景區(qū)。