磁總場(chǎng)梯度數(shù)據(jù)三維反演方法及應(yīng)用

孫石達(dá),張青杉,陳 超,王浩然,陳海弟,戴繼舒

(1. 中國(guó)冶金地質(zhì)總局地球物理勘查院，河北保定 071000；2. 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074)

?

磁總場(chǎng)梯度數(shù)據(jù)三維反演方法及應(yīng)用

孫石達(dá)1,2,張青杉1,2,陳 超2,王浩然2,陳海弟1,戴繼舒1

(1. 中國(guó)冶金地質(zhì)總局地球物理勘查院，河北保定 071000；2. 中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430074)

磁總場(chǎng)三方位梯度數(shù)據(jù)相對(duì)于總場(chǎng)數(shù)據(jù)包含更豐富的異常信息，將梯度數(shù)據(jù)應(yīng)用于三維磁化率反演中，可以更準(zhǔn)確地描述異常體。本文采用最小模型結(jié)構(gòu)反演方法進(jìn)行三維磁化率反演成像，并采用對(duì)數(shù)障礙法對(duì)磁化率的反演取值范圍進(jìn)行約束。通過(guò)模型試驗(yàn)，對(duì)磁總場(chǎng)異常數(shù)據(jù)及其三軸梯度數(shù)據(jù)進(jìn)行單獨(dú)反演、聯(lián)合反演，結(jié)果表明三軸梯度聯(lián)合反演結(jié)果可以更好地刻畫(huà)異常體形態(tài)，更有效地分辨鄰近異常體，反演物性更合乎實(shí)際。將此反演技術(shù)應(yīng)用于大冶鐵礦高精度航磁數(shù)據(jù)反演解釋?zhuān)〉昧溯^好的應(yīng)用效果。

磁總場(chǎng)梯度 磁化率 最小模型結(jié)構(gòu)反演 對(duì)數(shù)障礙法

Sun Shi-da, Zhang Qing-shan, Chen Chao, Wang Hao-ran, Chen Hai-di, Dai Ji-shu. 3D inversion of gradients of total field magnetic anomalies and its application[J]. Geology and Exploration, 2015, 51(6): 1016-1024.

1 引言

磁法測(cè)量的對(duì)象除地磁場(chǎng)的總場(chǎng)強(qiáng)度外，還包括總場(chǎng)空間變化率(梯度)、總場(chǎng)三分量及其空間變化率(梯度張量)。其中，總場(chǎng)強(qiáng)度的梯度測(cè)量始于20世紀(jì)80年代，并已取得了較好的應(yīng)用效果(Hood,1981；Hoodetal.,1989)。由于地磁場(chǎng)三分量測(cè)量受測(cè)量設(shè)備精度及姿態(tài)影響嚴(yán)重，目前尚未得到較為理想的測(cè)量結(jié)果，Schmidtetal.(1997,1998)證實(shí)由觀測(cè)磁總場(chǎng)強(qiáng)度進(jìn)而推算其三軸分量更加切實(shí)可行。磁梯度張量數(shù)據(jù)由于其包含了更為豐富的信息而被認(rèn)為是磁法勘探的下一個(gè)突破點(diǎn)(Christensenetal.,2000)，Schmidtetal.(2006)介紹了梯度張量測(cè)量相比分量測(cè)量和總場(chǎng)梯度測(cè)量的優(yōu)勢(shì)，并討論了梯度張量數(shù)據(jù)的應(yīng)用前景。目前，國(guó)內(nèi)航空磁測(cè)仍以觀測(cè)磁總場(chǎng)強(qiáng)度為主，梯度測(cè)量方面垂直梯度測(cè)量開(kāi)展相對(duì)較早(姚長(zhǎng)利等，1996)，并且近幾年也實(shí)現(xiàn)了航磁三方向梯度測(cè)量工作(李曉祿等，2009；郭華等，2013)。中國(guó)冶金地質(zhì)總局物勘院于20世紀(jì)80年代開(kāi)始使用核旋磁力儀開(kāi)展直升機(jī)吊艙式垂向磁總場(chǎng)梯度測(cè)量工作，效果較為理想；2007年后，該院使用光泵磁力儀開(kāi)展三方位磁總場(chǎng)梯度測(cè)量工作，工作模式為直升機(jī)吊艙式及固定翼軟補(bǔ)償式，取得了較好成果。

針對(duì)國(guó)內(nèi)磁法勘探現(xiàn)狀，本文將磁總場(chǎng)梯度數(shù)據(jù)作為主要研究對(duì)象。磁總場(chǎng)梯度測(cè)量較總場(chǎng)測(cè)量而言，其原始資料具備如下主要優(yōu)點(diǎn)(Hoodetal.,1989；秦葆瑚，1985；張昌達(dá)，2006；Schmidtetal.,2006；張青杉等，2010)：

(1) 數(shù)據(jù)質(zhì)量較高，各探頭受外界磁干擾比較一致，且梯度觀測(cè)儀器的精度要求高于總場(chǎng)觀測(cè)，一定程度上提高了磁異常圖的質(zhì)量；

(2) 數(shù)據(jù)類(lèi)別及數(shù)量有較大幅度提升，而且梯度數(shù)據(jù)包含垂直測(cè)線方向上一定距離的信息，更利于欠采樣觀測(cè)數(shù)據(jù)的解釋?zhuān)岣吡舜艤y(cè)效能，同時(shí)不同梯度數(shù)據(jù)間可相互印證，更多的數(shù)據(jù)處理方法更利于推斷解釋工作；

(3) 梯度測(cè)量可以消除區(qū)域磁場(chǎng)和隨時(shí)間變化磁場(chǎng)(地磁日變及磁暴)的影響，因此可以不設(shè)置日變站及進(jìn)行日變校正。

就數(shù)據(jù)解譯而言，總場(chǎng)梯度數(shù)據(jù)尤其水平梯度對(duì)于場(chǎng)源體的邊界位置具有更高的分辨能力(張青杉，2009)，許多處理方法如歐拉反褶積法(Thompson,1982；王明等，2012)、總梯度模法(Roestetal.,1992；黃臨平等，1998；張季生，2000)均需使用梯度信息推算磁性體邊界位置或估算頂面邊緣埋深等，而且總梯度模法還能夠在一定程度上消除剩磁影響；石磊等(2012)對(duì)磁總場(chǎng)及其垂向梯度數(shù)據(jù)采用相關(guān)成像技術(shù)進(jìn)行處理，結(jié)果表明梯度數(shù)據(jù)效果更好。

本文針對(duì)磁化率參數(shù)研究適用最小模型結(jié)構(gòu)反演方法，并通過(guò)對(duì)數(shù)障礙法對(duì)磁化率范圍施加約束。通過(guò)多組模型試驗(yàn)，對(duì)(單或多)梯度數(shù)據(jù)及總場(chǎng)數(shù)據(jù)的反演結(jié)果進(jìn)行對(duì)比，可以看出兩者均有準(zhǔn)確可靠的表現(xiàn)。在磁性體形態(tài)描述的精細(xì)性、鄰近磁性體的分離能力以及反演磁化率與其真值的偏離度等方面，梯度(尤其是多梯度)數(shù)據(jù)反演結(jié)果更為優(yōu)越。將上述反演技術(shù)應(yīng)用于大冶地區(qū)高精度航磁數(shù)據(jù)的反演解釋?zhuān)Ч^為理想。

2 方法原理

2.1 磁總場(chǎng)及其梯度正演

磁總場(chǎng)及其梯度正演公式采用駱遙等(2007)提出的長(zhǎng)方體磁場(chǎng)及其梯度無(wú)解析奇點(diǎn)表達(dá)式。假設(shè)地下空間剖分為M個(gè)規(guī)則網(wǎng)格排列的長(zhǎng)方體，觀測(cè)面上存在N個(gè)數(shù)據(jù)觀測(cè)點(diǎn)，則第j個(gè)長(zhǎng)方體單元在第i個(gè)觀測(cè)點(diǎn)處的磁總場(chǎng)或其梯度異常均可簡(jiǎn)化為：

(1)

其中，mj為第j個(gè)單元的磁化率，Gii為根據(jù)正演公式計(jì)算得到的核系數(shù)，其只與正常場(chǎng)的磁化強(qiáng)度及方向、長(zhǎng)方體單元的剩磁以及與觀測(cè)點(diǎn)的相對(duì)位置有關(guān)。

由于位場(chǎng)的疊加性，整個(gè)模型空間對(duì)第i個(gè)觀測(cè)點(diǎn)的異?？赏ㄟ^(guò)所有單元在該點(diǎn)產(chǎn)生異常值的線性疊加，因此可將所有數(shù)據(jù)點(diǎn)的正演計(jì)算寫(xiě)為矩陣形式：

T=Gm

(2)

其中G為N×M維矩陣，可稱(chēng)為核矩陣或敏感度矩陣，其中第i行第j列元素即為式(1)中核系數(shù)，T和m分別為正演結(jié)果和模型磁化率構(gòu)成的N維向量和M維向量。

2.2 最小模型結(jié)構(gòu)反演

磁化率反演的目的是根據(jù)磁測(cè)數(shù)據(jù)求取比較符合實(shí)際地質(zhì)情況的地下空間磁化率分布。為使反演結(jié)果盡量精細(xì)，一般情況下地下空間剖分單元個(gè)數(shù)遠(yuǎn)多于觀測(cè)數(shù)據(jù)點(diǎn)數(shù)，反演過(guò)程需要求解欠定方程組，同時(shí)由于觀測(cè)數(shù)據(jù)中存在噪聲和誤差，僅僅通過(guò)精確擬合觀測(cè)數(shù)據(jù)所獲得的反演結(jié)果未必理想，因此需要加入某些確定的模型約束來(lái)控制解模型空間的結(jié)構(gòu)，改善反演的病態(tài)程度，降低多解性，這種反演方法統(tǒng)稱(chēng)為最小模型結(jié)構(gòu)反演，目前已被廣泛應(yīng)用于地球物理數(shù)據(jù)反演中(Lietal.,1998；Farquharson,2008；Sunetal.,2014)，而加入模型結(jié)構(gòu)約束的過(guò)程也被稱(chēng)為正則化(Tikhonovaetal.,1977)。

最終的反演目標(biāo)函數(shù)通常包括兩部分：數(shù)據(jù)擬合差函數(shù)Φd和模型目標(biāo)函數(shù)Φm，前者體現(xiàn)的是反演結(jié)果正演數(shù)據(jù)與觀測(cè)數(shù)據(jù)的擬合程度，后者則表示反演結(jié)果的模型復(fù)雜程度和光滑程度，兩者相對(duì)重要性由正則化參數(shù)β來(lái)確定。因此目標(biāo)函數(shù)形式可寫(xiě)為：

(3)

目前常用的確定最優(yōu)正則化參數(shù)的方法主要有L曲線法和廣義交叉驗(yàn)證法(Asteretal.,2005)，并且均被較為廣泛的應(yīng)用到了重磁單數(shù)據(jù)反演中(Li and Oldenburg,2003；Farquharson and Oldenburg,2004)，但是由于后者目前較難應(yīng)用于多數(shù)據(jù)聯(lián)合反演的情況，因此本文中選用L曲線法來(lái)確定最優(yōu)正則化參數(shù)，具體方法可參見(jiàn)Asteretal.(2005)，這里不再贅述。

本文以磁化率為反演參數(shù)，可將Φd和Φm表述為常規(guī)的計(jì)算矢量長(zhǎng)度形式，式(3)演化為：

(4)

其中dobs為觀測(cè)數(shù)據(jù)，G為依據(jù)所剖分地下空間及數(shù)據(jù)觀測(cè)點(diǎn)位正演計(jì)算得到的核矩陣，m為磁化率模型空間，Gm即為所求取數(shù)據(jù)的正演表達(dá)，Wd為依據(jù)觀測(cè)數(shù)據(jù)標(biāo)準(zhǔn)差確定的對(duì)角矩陣；Ws和Wi(i=x,y,z)分別是測(cè)算模型空間m模型復(fù)雜度的單位矩陣和三方向粗糙度的差分矩陣(Lietal.,1996)；mref為參考模型。Φg表示測(cè)度矢量長(zhǎng)度的函數(shù)，不同的測(cè)度函數(shù)具有不同的性質(zhì)(Farquharson,2008；杜勁松，2014；Sunetal.,2014)，可導(dǎo)致反演結(jié)果呈現(xiàn)不同的分布特征，本文采用常規(guī)的L2范數(shù)測(cè)度形式，式(4)推演如下：

(5)

這是應(yīng)用較多的一種形式，得到的物性結(jié)果較為平滑，并且由于磁正演是線性問(wèn)題，因此采用二階范數(shù)求解時(shí)反演過(guò)程也是線性的，最終求解公式為：

(6)

上式適用于針對(duì)磁總場(chǎng)或磁梯度單一觀測(cè)參數(shù)的三維反演，對(duì)于三方位磁總場(chǎng)梯度數(shù)據(jù)的聯(lián)合反演問(wèn)題，反演目標(biāo)函數(shù)可表述為將針對(duì)三個(gè)參量的反演目標(biāo)函數(shù)加權(quán)疊加(由于三個(gè)方位磁梯度正演核函數(shù)隨深度衰減特性較為一致，本文采用相同的深度加權(quán)函數(shù)，以便于目標(biāo)函數(shù)的實(shí)現(xiàn))，即將式(6)中數(shù)據(jù)擬合差部分?jǐn)U展為三梯度數(shù)據(jù)擬合差即可：

式中μl為各梯度數(shù)據(jù)反演權(quán)重，其余參數(shù)同式(6)。

(7)

2.3 物性約束

為避免出現(xiàn)有悖于地質(zhì)理論的物性反演結(jié)果，提高反演結(jié)果的可靠性，合理描述地質(zhì)體形態(tài)，對(duì)反演模型的物性范圍進(jìn)行合理約束是必要的。本文采用對(duì)數(shù)障礙法(Logarithmic barrier method)來(lái)實(shí)現(xiàn)物性約束。該方法最早用于解決不等值約束條件下的線性和二次規(guī)劃問(wèn)題(e.g.Gill et al.,1991；Rooset al.,1997)，是一種近似全局精確的范圍約束方法(白富生等，2000)，也被廣泛應(yīng)用在多種重磁位場(chǎng)數(shù)據(jù)反演中并取得了較好的物性范圍約束效果(e.g.LiandOldenburg,2003；Farquharson,2008；Martinezet al.,2013)。假設(shè)反演磁化率值約束在(m-,m+)的范圍內(nèi)，按照對(duì)數(shù)障礙法把該不等值約束轉(zhuǎn)換為式(9)所示的障礙函數(shù)并將其作為反演目標(biāo)函數(shù)的一部分，因此式(3)擴(kuò)展為：

(8)

其中λ為障礙參數(shù)，Xbf為障礙函數(shù)：

(9)

(10)

使用對(duì)數(shù)障礙法約束反演物性時(shí)，給定的初始模型值m(0)需全部處于物性約束范圍(m-,m+)內(nèi)，而初始障礙參數(shù)λ(0)可通過(guò)式(10)求得。由于加入障礙函數(shù)使得式(5)的線性反演問(wèn)題非線性化，因此本文采用加權(quán)迭代最小二乘方法(Iteratively Reweighted Least Squares，IRLS)求解此非線性問(wèn)題，求解思路為用第(n-1)次迭代的觀測(cè)數(shù)據(jù)擬合差計(jì)算第n次迭代的物性模型修改量。求解物性模型修改量Δm的公式推導(dǎo)如下(以單參量反演為例，對(duì)于多參量聯(lián)合反演，只需比照式(7)將數(shù)據(jù)擬合差部分?jǐn)U展為多參量即可)。

設(shè)初始物性模型為m(0)；設(shè)第(n-1)次迭代的物性模型為m(n-1)，其正演結(jié)果為d(n-1)(等同于Gm(n-1))，則第n次迭代可視為用物性模型變化量Δm擬合dobs與d(n-1)的差值；第n次迭代得到的物性模型即為m(n)=m(n-1)+Δm；因此第n次迭代反演目標(biāo)函數(shù)由(8)式演化為：

(11)

在僅考慮單一參量反演的情況下，求取第n次Δm的最終表達(dá)式由式(6)變?yōu)椋?/p>

(12)

式中，X1和X2均為對(duì)角矩陣，分別由障礙函數(shù)對(duì)物性模型向量的一階和二階求導(dǎo)所得，λ(n)為第n次障礙參數(shù)，通過(guò)上式即可求得第n次Δm。

為使物性模型保持在約束范圍內(nèi)，第n次物性結(jié)果可表示為：

(13)

其中γ為處于(0,1)區(qū)間內(nèi)近于1的常數(shù)，η為最大允許步長(zhǎng)，用于將m(n)約束在(m-,m+)內(nèi)，其值可由(14)式求得。

(14)

為使反演過(guò)程穩(wěn)定收斂，應(yīng)使障礙函數(shù)對(duì)反演目標(biāo)函數(shù)的影響逐步降低，故將下一次迭代的障礙參數(shù)設(shè)定為λ(n+1)=[1-min(η,γ)]λ(n)，可見(jiàn)障礙參數(shù)隨著迭代次數(shù)越來(lái)越小，而障礙函數(shù)對(duì)整個(gè)目標(biāo)函數(shù)的作用也相應(yīng)越來(lái)越弱，當(dāng)達(dá)到規(guī)定的終止條件時(shí)(如反演目標(biāo)函數(shù)變化量小于整體1%、障礙函數(shù)小于反演目標(biāo)函數(shù)1%)即可終止迭代，求得最終物性模型結(jié)果。

3 模型試驗(yàn)

為驗(yàn)證上述反演方法的可靠性，測(cè)試磁總場(chǎng)及其梯度數(shù)據(jù)的反演效果，以下采用傾斜體模型及兩組組合模型進(jìn)行試驗(yàn)。通過(guò)對(duì)磁總場(chǎng)及其梯度數(shù)據(jù)進(jìn)行單參量反演、多參量聯(lián)合反演，對(duì)反演結(jié)果在異常體形態(tài)描述準(zhǔn)確性、鄰近異常體區(qū)分能力以及與真實(shí)物性值接近程度等方面的表現(xiàn)特征進(jìn)行分析。

圖1 傾斜體模型三維顯示圖及其正演磁異常圖Fig.1 Perspective view of tilting model and its magnetic anomalies of forward modelinga-三維模型；b-磁總場(chǎng)異常ΔT；c-磁總場(chǎng)北向梯度Tx；d-磁總場(chǎng) 東向梯度Ty；e-磁總場(chǎng)垂向梯度Tza-3D model;b-total-field magnetic anomalies ΔT;c-Tx,d-Ty;e-Tzare the total magnetic gradient in x,y and z directions, respectively

3.1 傾斜體模型

該模型因其傾斜特性多被應(yīng)用于磁物性反演試驗(yàn)中(Lietal.,1996；Boulangeretal.,2001)。將地下空間剖分為21×21×10=4410個(gè)網(wǎng)格單元，網(wǎng)格單元邊長(zhǎng)均為50 m，傾斜體模型三維立體顯示如圖1a，其傾斜角度為45°，傾向東，傾斜體磁化率為0.1 SI，背景空間磁化率為0 SI(本文中所有模型試驗(yàn)均假定磁性體無(wú)剩磁，自退磁效應(yīng)忽略不計(jì))；地磁場(chǎng)傾角45°，偏角10°，強(qiáng)度為57000 nT；正演計(jì)算得到磁總場(chǎng)異常ΔT及其梯度異常如圖1所示。為計(jì)算反演結(jié)果正演后與原始模型正演數(shù)據(jù)的偏差，在此未對(duì)異常數(shù)據(jù)加入任何噪聲。

對(duì)上述ΔT及其梯度異常數(shù)據(jù)進(jìn)行磁化率反演，反演過(guò)程中僅加入物性范圍約束，約束范圍為0～0.1 SI，正則化參數(shù)通過(guò)L曲線方法確定，得到最終磁化率反演結(jié)果(圖2)?？梢钥闯?，ΔT及其三軸梯度單獨(dú)反演結(jié)果大體相當(dāng)，Tx、Ty反演結(jié)果分別在x及y方向表現(xiàn)略好；而三軸梯度聯(lián)合反演結(jié)果則表現(xiàn)較好，既繼承了單梯度數(shù)據(jù)反演結(jié)果的優(yōu)點(diǎn)，能夠?qū)Ξ惓ｓw邊界進(jìn)行較準(zhǔn)確的刻畫(huà)，同時(shí)對(duì)深部描述也更好。從反演模型及原始模型分別正演后標(biāo)準(zhǔn)偏差統(tǒng)計(jì)結(jié)果(表1)看，ΔT反演結(jié)果不能很好地?cái)M合梯度數(shù)據(jù)，單梯度反演結(jié)果也不能很好地?cái)M合ΔT及其它梯度數(shù)據(jù)，而三軸梯度聯(lián)合反演結(jié)果對(duì)所有數(shù)據(jù)均擬合較好，具備真實(shí)刻畫(huà)實(shí)際異常體的能力。

圖2 傾斜體模型正演數(shù)據(jù)反演結(jié)果切片對(duì)比圖Fig.2 Comparison of slices for forward and inversion results recovered from the anomalies in Fig.1a-理論模型；b-ΔT反演結(jié)果；c-Tx反演結(jié)果；d-Ty反演結(jié)果；e-Tz 反演結(jié)果；f-三方向梯度聯(lián)合反演結(jié)果Slices through (a)the model and results recovered from (b)ΔT,(c) Tx,(d) Ty,(e) Tz and (f) joint three gradients

3.2 組合模型1

為進(jìn)一步測(cè)試上述反演方法的效果，采用如圖3a的組合模型，使水平x、y方向均存在高-高磁與高-低磁組合異常體，以便考察三軸梯度反演結(jié)果對(duì)鄰近異常體的區(qū)分分辨能力。將地下空間剖分為31×31×10= 9610個(gè)網(wǎng)格單元，網(wǎng)格單元邊長(zhǎng)均為50 m；1、2、3號(hào)異常體磁化率均為0.1 SI，4號(hào)異常體磁化率為0.05 SI，背景空間磁化率0 SI；假定地磁場(chǎng)傾角45°，偏角10°，強(qiáng)度為57000 nT；正演計(jì)算得到磁總場(chǎng)異常ΔT及其梯度異常(均加入異常幅值2%的高斯噪聲)如圖3b～e所示。

表1 傾斜體下反演模型與原始模型正演結(jié)果標(biāo)準(zhǔn)差Table 1 Standard deviations between inversion and forward modeling for a tilting body

圖3 組合模型1三維顯示圖及其正演磁異常圖Fig. 3 Same as Fig. 1 but on combined model 1a-三維模型；b-磁總場(chǎng)異常ΔT；c-磁總場(chǎng)北向梯度Tx；d- 磁總場(chǎng)東向梯度Ty；e-磁總場(chǎng)垂向梯度Tza-3D model;b-total-field magnetic anomaly ΔT;c-Tx;d-Ty;e-Tz are the total magnetic gradient in x,y and z directions,re- spectively

反演過(guò)程同3.1節(jié)。為便于查看橫向分辨力，本文給出深度200 m處切片圖(圖4)?？梢钥闯?，ΔT反演結(jié)果(圖4b)的異常體邊緣刻畫(huà)能力弱于梯度數(shù)據(jù)反演結(jié)果，當(dāng)異常體靠近時(shí)，ΔT反演結(jié)果中兩者之間的邊界或空隙漸漸模糊，容易形成兩者連為一體的假象；Tx、Ty單獨(dú)反演結(jié)果(圖4c～d)在各自對(duì)應(yīng)方向上能夠較為清晰地區(qū)分異常體；Tz單獨(dú)反演(圖4e)與三梯度聯(lián)合反演(圖4f)結(jié)果在切片深度具有較為一致的異常體分辨能力，能夠清晰地分辨所有異常體。

圖4 組合模型1正演數(shù)據(jù)反演結(jié)果切片對(duì)比圖Fig. 4 Same as Fig. 2 but on combined model 1a-理論模型；b-ΔT反演結(jié)果；c-Tx反演結(jié)果；d-Ty反演結(jié)果；e-Tz 反演結(jié)果；f-三方向梯度聯(lián)合反演結(jié)果Slices through (a)the model and results recovered from (b)ΔT,(c) Tx,(d) Ty,(e) Tz and (f) joint three gradients

3.3 組合模型2

為檢測(cè)反演技術(shù)的綜合能力，本文模擬某地實(shí)際地質(zhì)情況建立模型(圖5a)，地下空間剖分為31×31×10 = 9610個(gè)網(wǎng)格單元，網(wǎng)格單元邊長(zhǎng)均為50 m；1號(hào)異常體磁化率0.05 SI， 2號(hào)異常體磁化率0.1 SI， 3號(hào)異常體磁化率0.15 SI；地磁場(chǎng)傾角45°，偏角10°，強(qiáng)度為57000 nT；令異常觀測(cè)面為0.1 m，正演計(jì)算得到ΔT及其梯度異常(均加入異常幅值2%的高斯噪聲)如圖5b～e所示。

圖5 組合模型2三維顯示圖及其正演磁異常圖Fig. 5 Same as Fig.3 but on combined model 2a-三維模型；b-磁總場(chǎng)異常ΔT；c-磁總場(chǎng)北向梯度Tx；d-磁 總場(chǎng)東向梯度Ty；e-磁總場(chǎng)垂向梯度Tza-3D model;b-total-field magnetic anomaly ΔT;c-Tx;d-Ty;e-Tz are the total magnetic gradient in x,y and z direction respec tively

在此僅分析比對(duì)ΔT和三梯度聯(lián)合反演結(jié)果(圖6b～c)，可以看出，后者能量更為集中，也對(duì)異常體邊界描述更準(zhǔn)確，反演磁化率更接近實(shí)際。

圖6 組合模型2正演數(shù)據(jù)反演結(jié)果切片對(duì)比圖Fig. 6 Same as Fig.4 but on combined model 2a-理論模型；b-ΔT反演結(jié)果；c-三方向梯度聯(lián)合反演結(jié)果Slices through (a)the model and results recovered from (b)ΔT and (c) joint three gradients

通過(guò)上述理論模型試驗(yàn)，三梯度聯(lián)合反演結(jié)果與ΔT反演結(jié)果相比較，能夠更為準(zhǔn)確地刻畫(huà)異常體形態(tài)，橫向上可更好地區(qū)分鄰近異常體，反演磁化率更接近實(shí)際。

4 實(shí)例應(yīng)用

大冶鐵(銅)礦床位處鐵山巖體南緣中段與三疊系下統(tǒng)大冶組大理巖接觸帶，頂?shù)装鍑鷰r分別為中細(xì)粒含石英閃長(zhǎng)斑巖及黑云母透輝石閃長(zhǎng)巖。據(jù)巖礦石物性資料，礦區(qū)鐵礦石實(shí)測(cè)磁化率平均值約0.8 SI，黑云母透輝石閃長(zhǎng)巖和閃長(zhǎng)斑巖磁化率稍低，大理巖磁化率最低(朱永剛等，2006；高寶龍，2010；陶德益等，2011)。礦區(qū)1∶1萬(wàn)航磁ΔT數(shù)據(jù)經(jīng)匹配濾波剔除區(qū)域場(chǎng)后用于反演(圖7a)，圖7b～d為推算的三方位梯度異常圖。

圖7 大冶礦區(qū)航磁總場(chǎng)異常及換算三方位梯度異常圖Fig.7 Total-field magnetic anomalies in the Daye Mine and its derived gradients a-磁總場(chǎng)ΔT局部異常圖；b-磁總場(chǎng)北向梯度Tx；c-磁總場(chǎng)東向梯度Ty；d-磁總場(chǎng)垂向梯度Tza-total-field magnetic anomaly ΔT；b-Tx;c-Ty;d-Tz are the total mag-netic gradient in x,y and z direction respectively derived from ΔT

對(duì)上述數(shù)據(jù)進(jìn)行ΔT反演及三梯度聯(lián)合反演(不考慮剩磁影響)，將地下模型空間剖分為56×36×20個(gè)單元，每個(gè)單元邊長(zhǎng)均為100 m；礦區(qū)地磁場(chǎng)傾角45.6，偏角-3；磁化率約束范圍0～0.8 SI；反演結(jié)果如圖8。ΔT反演及三梯度聯(lián)合反演結(jié)果中的高值異常自西北向東南分別與鐵門(mén)坎、龍洞、象鼻山、獅子山四個(gè)礦區(qū)對(duì)應(yīng)良好，并且對(duì)四個(gè)異常體都能較為清晰的刻畫(huà)。相比較而言，后者結(jié)果聚焦更好，且反演磁化率值更接近實(shí)際，可以更好的描述異常體淺部以及深部的形態(tài)，能夠更為清晰的標(biāo)示出可能存在礦體的高磁性異常體的空間分布范圍。結(jié)果中(圖8b)磁化率高值區(qū)南西側(cè)存在北西西-南東東走向的連續(xù)低值帶，由于該區(qū)礦產(chǎn)主要賦存于燕山期含石英閃長(zhǎng)巖與三疊系下統(tǒng)大冶群灰?guī)r的接觸帶上，灰?guī)r在接觸帶附近已變質(zhì)為大理巖及白云質(zhì)大理巖(高寶龍等，2010)，因此推測(cè)該低值帶為無(wú)磁或低磁性的大理巖或灰?guī)r，同時(shí)三維立體圖(圖8f)中顯示高磁性體的走向與當(dāng)?shù)刂饕獦?gòu)造方向即北西西向相同，其傾向以及與接觸帶的接觸面整體向南西方向傾斜，與薛清潑等(2006)結(jié)論相同。

圖8 ΔT反演結(jié)果及三梯度聯(lián)合反演結(jié)果水平切片圖、斷面圖及三維顯示圖Fig. 8 Horizontal and vertical sections and perspective view of the inversion results of ΔT and joint inversion of three gradientsΔT反演結(jié)果：a-550 m深度切片圖；c-AA’斷面；e-三維圖；三梯度聯(lián)合反演結(jié)果：b-550 m深度切片圖；d-AA'斷面；f- 三維圖ΔT inversion result: a-horizontal sections at 550m depth; c-vertical section AA’; e-3D result; three-dimensional joint inversion results: b-horizontal sections at 550m depth; d-ertical section AA'; e-vertical section AA'; f-3D result

5 結(jié)論

本文將基于對(duì)數(shù)障礙法進(jìn)行物性范圍約束的最小模型結(jié)構(gòu)反演方法應(yīng)用于磁總場(chǎng)及其梯度數(shù)據(jù)反演中，通過(guò)模型數(shù)據(jù)驗(yàn)證了方法的可靠性，并將該反演方法應(yīng)用于大冶礦區(qū)航磁數(shù)據(jù)反演，反演結(jié)果與實(shí)際地質(zhì)情況較為一致?？傮w而言，ΔT數(shù)據(jù)反演結(jié)果對(duì)異常體邊界刻畫(huà)稍顯模糊，異常體橫向區(qū)分能力偏弱，反演磁化率與實(shí)際偏差稍大；水平單梯度數(shù)據(jù)反演結(jié)果偏重于刻畫(huà)異常體相應(yīng)邊界，具備較好的橫向分辨能力；垂向梯度數(shù)據(jù)反演結(jié)果在深部表現(xiàn)較好；三軸梯度數(shù)據(jù)聯(lián)合反演結(jié)果保持了單梯度反演結(jié)果的優(yōu)點(diǎn)，能夠更準(zhǔn)確地描述異常體的整體形態(tài)。值得指出的是，梯度數(shù)據(jù)反演結(jié)果會(huì)出現(xiàn)反演物性值遠(yuǎn)大于實(shí)際值情況，需要更準(zhǔn)確的物性約束以確保反演結(jié)果真實(shí)可靠。

本文中模型試驗(yàn)及實(shí)際資料處理均未考慮剩磁，由于總梯度模對(duì)磁性體總磁化方向不敏感，后續(xù)研究中將補(bǔ)充梯度數(shù)據(jù)與總梯度模數(shù)據(jù)的聯(lián)合反演，使本文反演技術(shù)能夠適用于強(qiáng)剩磁情況。

Aster R C,Borchers B,Thurber C H.2011.Parameter Estimation and Inverse Problems [M].The second edition.Elsevier: 89-109

Bai Fu-sheng,Zhang Liang-sheng.2000.Approximate global exact bnarrier method for logarithmic barrier function [J].OR Transactions,4(3): 13-18

Boulanger O,Chouteau M.2001.Constraints in 3D gravity inversion [J].Geophysical Prospecting,49: 265-280

Christensen A,Rajagopalan S.2000.The magnetic vector and gradient tensor in mineral and oil exploration [J].Preview,84: 77

Du Jin-song.2014.Study on processing,forward modeling and inversion algorithms of satellite magnetic anomaly data in spherical coordinate system [D].Wuhan: China University of Geosciences(Wuhan): 84-85(in Chinese)

Farquharson C G.2008.Constructing piecewise-constant models in multidimensional minimum-structure inversions [J].Geophysics,73(1): K1-K9

Farquharson C G,Oldenburg D W.2004.A comparison of automatic techniques for estimating the regularization parameter in non-Linear inverse problems [J].Geophysical Journal International,156(3):411-425

Gao Bao-long,Tao De-yi,Zhan Ying-lin,Fan Zhi-xiong,Zhou Kui.2010.Application of aero-surface and borehole magnetic exploration to the prospecting of exhausted mines in the Daye iron mine [J].Geology and Exploration,46(3): 0483-0490(in Chinese with English abstract)

Gill P E,Murray W,Poncele D B.1991.Solving reduced KKT systems in barrier methods for linear and quadratic programming[Z].Technical Report SOL 91-7,Stanford University,Stanford,CA: 1-23

Guan Zhi-ning,Hou Jun-sheng,Yao Chang-li.1996.Application of aeromagnetic gradient data in geological mapping and metallogenetic prognosis of gold deposits [J].Geoscience Journal of Graduate School,China University of Geosciences,10(02): 239-249(in Chinese with English abstract)

Guo Hua,Wang Ping,Xie Ru-kuan.2013.A study of geological interpretation with the tri-axial aeromagnetic gradients [J].Progress in Geophys,28(05): 2688-2692(in Chinese with English abstract)

Huang Lin-ping,Guan Zhi-ning.1998.The determination of magnetic causative boundaries using total gradient modules of magnetic anomalies [J].Journal of East China Geological Institute,21(2): 44-51(in Chinese with English abstract)

Hood PJ. 1980. Aeromagnetic gradiometry: A superior geological mapping tool for mineral exploration programs. In: Weinstock H, and Overton WC. (eds.), SQUID applications to geophysics. Society of Exploration Geophysicists: 72-77

Hood P J,Teskey D J.1989.Aeromagnetic gradiometer program of the Geological Survey of Canada [J].Geophysics,54(8): 1012-1022

Li Xiao-lu,Chang Shu-shuai.2009.Aeromagnetic gradient survey and elementary application in sandstone type uranium deposits prospecting [J].Uranium Geology,25(06): 355-360(in Chinese with English abstract)

Li Y G,Oldenburg D W.1996.3-D inversion of magnetic data[J].Geophysics,61(2): 394-408

Li Y G,Oldenburg D W.1998.3-D inversion of gravity data[J].Geophysics,63(1): 109-119

Li Y G,Oldenburg D W.2003.Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method [J].Geophysical Journal International,152(2): 251-265

Luo Yao,Yao Chang-li.2007.Theoretical study on cuboid magnetic field and its gradient expression without analytic singular point [J].OGP,42(6): 714-719(in Chinese with English abstract)

Martinez C,Li Y,Krahenbuhl R.2013.3D inversion of airborne gravity gradiometry data in mineral exploration: A case study in the Quadrilátero Ferrífero,Brazil [J].Geophysics,78(1): B1-B11

Roest W R,Verhoef J,Pilkington M.1992.Magnetic interpretation using the 3-D analytic signal [J].Geophysics,57(1):116-125

Roos C,Terlaky T,Vial J.1997.Theory and algorithms for linear optimization: An interior point approach [M].Wiley Chichester:109-144

Schmidt P W,Clark D A.1997.Directions of magnetization and vector anomalies derived from total field surveys [J].Preview,70: 30-32

Schmidt P W,Clark D A.1998.The calculation of magnetic components and moments from TMI: A case study from the Tuckers igneous complex,Queensland [J].Exploration Geophysics,29: 609-614

Schmidt P W,Clark D A.2006.The magnetic gradient tensor: Its properties and uses in source characterization [J].The Leading Edge,25(1): 75

Shi Lei,Guo Liang-hui,Meng Xiao-hong.2012.3D correlation imaging of the vertical gradient of magnetic total field anomaly [J].Progress in Geophys,27(4): 1609-1614(in Chinese with English abstract)

Sun J J,Li Y G.2014.Adaptive Lp inversion for simultaneous recovery of both blocky and smooth features in a geophysical model [J].Geophysical Journal International,197: 882-899

Tao De-yi,Yang Hai-yan,Xiao Ming-yao,Gao Bao-long,Shu Xiu-feng.2011.The application effect of borehole magnetic exploration in the prospecting of the Daye iron mine [J].Resources Environment & Engineering,25(4): 358-363(in Chinese with English abstract)

Thompson D T.1982.EULDPH: A new technique for making computer-assisted depth estimated from magnetic data [J].Geophysics,47(1): 31-37

Tikhonov A V,Arsenin V Y.1977.Solution of Ill-Posed Problems[M].John Wiley and Sons:93-127

Wang Ming,Luo Yao,Luo Feng,Tian Song.2012.The application and development of Euler deconvolution in gravity and magnetic field [J].Geophysical & Geochemical Exploration,36(5): 834-841(in Chinese with English abstract)

Zhang Chang-da.2006.Airborne tensor magnetic gradiometry—The latest progress of airborne magnetometric technology [J].Chinese Journal of Engineering Geophysics,3(5): 354-361(in Chinese with English abstract)

Zhang Ji-sheng.2000.The application of 3D analytical signal technique to the treatment of aeromagnetic anomalies in south China [J].Geophysical & Geochemical Exploration,24(03): 190-196(in Chinese with English abstract)

Zhang Qing-shan,Ma Feng-lin,Xu Li-yun.2010.Study on the 3D gradient aeromagnetic survey [J].Geology and Exploration,46(6): 1087-1091(in Chinese with English abstract)

Zhu Yong-gang,Yu Chang-chun.2007.Analysis of magnetic properties of Daye iron deposit in Hubei province [J].Contributions to Geology and Mineral Resources Research,21(sup): 155-159(in Chinese with English abstract)

[附中文參考文獻(xiàn)]

白富生,張連生.2000.近似全局精確障礙方法：對(duì)數(shù)障礙函數(shù)法[J].運(yùn)籌學(xué)學(xué)報(bào),4(3): 13-18

杜勁松.2014.基于球坐標(biāo)系的衛(wèi)星磁異常數(shù)據(jù)處理與正反演方法研究[D].武漢：中國(guó)地質(zhì)大學(xué)(武漢)：84-85

高寶龍,陶德益,詹應(yīng)林,范志雄,周 逵,2010.大冶鐵礦接替資源勘查項(xiàng)目中“空、地、井”磁法測(cè)量的應(yīng)用[J].地質(zhì)與勘探,46(03):483-490

管志寧,侯俊勝,姚長(zhǎng)利.1996.航磁梯度資料在金礦地質(zhì)填圖和成礦預(yù)測(cè)中的應(yīng)用[J].現(xiàn)代地質(zhì),10(2): 239-249

郭 華,王 平,謝汝寬.2013.航磁全軸梯度數(shù)據(jù)地質(zhì)解釋優(yōu)勢(shì)研究[J].地球物理學(xué)進(jìn)展,28(05): 2688-2692

黃臨平,管志寧.1998.利用磁異?？偺荻饶４_定磁源邊界位置[J].華東地質(zhì)學(xué)院學(xué)報(bào),21(2): 44-51

李曉祿,常樹(shù)帥.2009.航磁梯度測(cè)量及其在砂巖型鈾礦勘查中的應(yīng)用初探[J].鈾礦地質(zhì),25(06): 355-360

駱 遙,姚長(zhǎng)利.2007.長(zhǎng)方體磁場(chǎng)及其梯度無(wú)解析奇點(diǎn)表達(dá)式理論研究[J].石油地球物理勘探,42(6): 714-719

秦葆瑚.1985.國(guó)外磁場(chǎng)梯度測(cè)量的新進(jìn)展.國(guó)外地質(zhì)勘探技術(shù)[J].(5): 18-24

石 磊,郭良輝,孟小紅.2012.磁總場(chǎng)異常垂直梯度三維相關(guān)成像[J].地球物理學(xué)進(jìn)展,27(4): 1609-1614

陶德益,楊海燕,肖明堯,高寶龍,舒秀峰.2011.井中磁測(cè)在大冶鐵礦深部勘查中的應(yīng)用效果[J].資源環(huán)境與工程,25(4): 358-363

王 明,駱 遙,羅 鋒,田 嵩.2012.歐拉反褶積在重磁位場(chǎng)中應(yīng)用與發(fā)展[J].物探與化探,36(5): 834-841

薛清潑,徐九華,陳偉,石教波,謝玉玲,2006.大冶鐵礦尖林山—獅子山礦段接觸帶形態(tài)及控礦規(guī)律分析[J].有色金屬(礦山部分),58(05): 14-16

張昌達(dá).2006.航空磁力梯度張量測(cè)量——航空磁測(cè)技術(shù)的最新進(jìn)展[J].工程地球物理學(xué)報(bào),3(5): 354-361

張季生.2000.用三維解析信號(hào)技術(shù)處理華南航磁異常[J].物探與化探,24(3): 190-196

張青杉.2009.發(fā)展三維磁梯度勘查技術(shù)，提升深部鐵礦勘查能力[C].中國(guó)地質(zhì)學(xué)會(huì)2009年學(xué)術(shù)年會(huì)論文摘要匯編: 594-596

張青杉.2010.三維磁梯度測(cè)量的高精度方案及其應(yīng)用前景[C].中國(guó)冶金地質(zhì)總局科技大會(huì)論文集：412-423

張青杉,麻豐林,許麗云.2010.航空三維磁梯度測(cè)量方案研究[J].地質(zhì)與勘探,46(5): 1087-1091

朱永剛,于長(zhǎng)春.2007.湖北省大冶鐵礦區(qū)內(nèi)礦石磁性特征分析[J].資源環(huán)境與工程,21(sup): 155-159

3D Inversion of Gradients of Total Field Magnetic Anomalies and its Application

SUN Shi-da2,ZHANG Qing-shan1,2,CHEN Chao2,WANG Hao-ran2,CHEN Hai-di1,DAI Ji-shu1

(1.GeophysicalExplorationAcademyofChinaMetallurgicalGeologyBureau,Baoding,Hebei071000;2.HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan,Hubei430074)

The 3-direction gradients of total-field magnetic anomalies (ΔT) contain more information compared with ΔT.Applying such data to 3D inversion of susceptibility allows us to describe anomaly bodies more accurately.This work employs the minimum-structure inversion method to conduct such inversion and uses the logarithmic barrier method to constrain the susceptibility limits.We compare the recovered models from inverting data ΔT,their gradients independently and jointly,separately.The results indicate that the inversion gradients,especially by the joint inversion,are capable of distinguishing nearby sources and can describe the sources more completely and exactly,and the inverted physical properties are closer to the true ones.And we apply this method to the aeromagnetic data acquired in Daye for inversion and interpretation,and obtain good results.

magnetic total-field gradients,susceptibility,minimum structure inversion,logarithmic barrier method

2015-03-28；

2015-07-09；[責(zé)任編輯]郝情情。

國(guó)家“973”計(jì)劃課題(編號(hào)：2012CB416805)和科研院所基本科研資金資助項(xiàng)目(編號(hào)： WHS201211)聯(lián)合資助。

孫石達(dá)(1990年-)，男，在讀博士生，主要從事重磁數(shù)據(jù)處理、反演與解釋方面研究。E-mail: shidasun.cug@gmail.com。

張青杉(1968年-)，男，中國(guó)冶金總局物勘院副總工程師，主要從事地球物理勘探與研究工作。E-mail: qingshan-zhang@163.com。

P631

A

0495-5331(2015)06-1016-9