球坐標系下磁化強度矢量反演方法及其應用

楊 力， 周可法， 王金林

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

0 引言

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

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

1 磁化強度矢量反演方法

1.1 直角坐標系下的反演方法和原理

磁感應強度大小與觀測距離之間的關系定義為式(1)。

(1)

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

M=κH+MNRM

(2)

其中:κ是磁化率；H是磁場強度；MNRM是天然剩磁。磁感應強度觀測值與m的關系可通過算子F定義為：

b=F(m)

(3)

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

m=F(m)-1dobs

(4)

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

φ(m)=φd+βφm

(5)

其中:Wd是標準差的倒數矩陣；W是空間距離加權矩陣，β是權重系數。對φ(m)求關于m的偏導數為“0”，可得：

βWTW(m-mref)=0

(6)

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

(7)

其中:Gx、Gy、Gz是三個方向上梯度的有限差分算子，[Ws,Wi,i=x,y,z]和[Rs,Ri,i=x,y,z]是稀疏化對角矩陣。敏感度矩陣Wr的選取能夠對反演結果產生重要的影響。建立迭代更新的敏感度權重方程[11]，可以寫為：

(8)

目前，航空磁測數據仍然以總磁場強度數據(TMI)為主，為了能夠利用大量磁測數據，本文選取總場強度數據進行計算?？倛鰪姸瓤梢远x為：

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

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

(9)

其中:H0是地磁場強度；P是將磁場轉換為總場強度的投影矩陣。在迪卡爾直角坐標系中，磁感應強度大小可以寫為：

b=T[(HpHsHt]κ

(10)

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

κp=MP/H0

κs=MP/H0

κt=MP/H0

(11)

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

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

(12)

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

圖1 直角坐標系下的模型水平切片和模型垂直切片Fig.1 Plan section through the recovered magnetization vector model using the cartesian formulation(a)模型水平切片；(b)模型垂直切片

圖2 球坐標系下的模型水平切片和模型垂直切片Fig.2 Plan section through the recovered magnetization vector model using the spherical formulation(a)模型水平切片；(b)模型垂直切片

1.2 球坐標系下的反演方法和原理

在球坐標系下，磁場強度大小{α}和磁化方向{θ,φ}能夠較好地解耦。與直角坐標系下的反演相比具有兩大優(yōu)勢，①野外磁測數據和實驗室磁測數據都會給出球坐標系下的參數值，將其轉換到直角坐標系會增加運算的繁瑣程度；②球坐標系下能分別對磁化強度大小和磁化方向實施稀疏化約束。

對模型范數項求解關于{α,θ,φ}的偏微分，其形式可以寫成：

(13)

J=FS

(14)

其中:矩陣S可以寫成：

(15)

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

(16)

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

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

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

2.1 礦區(qū)地質特征

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

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

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

2.2 航磁ΔT異常特征

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

航磁ΔT磁異常值變化范圍大，磁異常主要沿北西西-近東西向和北西向呈帶狀分布，與區(qū)域構造方向一致。M1、M2、M4磁異常帶比M3高，對應于研究區(qū)中的中性-基性巖石。M3低磁異常可能是被破碎帶、碎屑巖覆蓋的區(qū)域，由于ΔT值不能反映磁化率的三維分布結構，因此限制了對磁異常體的進一步分析。

表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)直角坐標系下反演結果;(b)球坐標系下反演結果

3 反演結果

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

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

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

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

4 結論

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

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