李澤林, 姚長(zhǎng)利, 鄭元滿, 孟小紅, 張聿文

地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室, 地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，中國(guó)地質(zhì)大學(xué)(北京)， 北京 100083

?

數(shù)據(jù)空間磁異常模量三維反演

李澤林, 姚長(zhǎng)利*, 鄭元滿, 孟小紅, 張聿文

地下信息探測(cè)技術(shù)與儀器教育部重點(diǎn)實(shí)驗(yàn)室, 地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，中國(guó)地質(zhì)大學(xué)(北京)， 北京 100083

強(qiáng)剩磁的存在通常導(dǎo)致了總磁化強(qiáng)度方向未知，進(jìn)而影響了磁異常的反演和解釋.磁異常模量是一種受磁化方向影響小的轉(zhuǎn)換量，可以在強(qiáng)剩磁條件下通過(guò)反演三維磁化強(qiáng)度大小分布來(lái)推測(cè)場(chǎng)源分布狀態(tài).我們提出了一種數(shù)據(jù)空間磁異常模量反演算法來(lái)減少剩磁的影響.與標(biāo)準(zhǔn)的模型空間L2范數(shù)正則化反演方法相比，我們的方法有兩個(gè)優(yōu)點(diǎn)：一是無(wú)需搜索正則化參數(shù)(需要反復(fù)求解非線性反演問(wèn)題)，因而可以減少計(jì)算時(shí)間；二是反演結(jié)果更加聚焦，深度分辨率更高，我們對(duì)此進(jìn)行了原因分析.通過(guò)模型和實(shí)測(cè)數(shù)據(jù)測(cè)試證明了該算法的有效性和更好的反演效果.

數(shù)據(jù)空間； 磁異常模量； 剩磁； 三維反演

1 引言

三維磁化率反演(Li and Oldenburg，1996，2003；Pilkington，1997，2009；Portniaguine and Zhdanov，2002；姚長(zhǎng)利等，2003，2007)逐步成為磁異常定量解釋中的一種重要方法，已經(jīng)成功地應(yīng)用到了礦產(chǎn)和油氣資源勘探以及區(qū)域地質(zhì)構(gòu)造解釋中.傳統(tǒng)的三維磁化率反演方法通常假設(shè)場(chǎng)源的磁化方向已知且與地磁場(chǎng)方向一致.然而，由于實(shí)際地質(zhì)情況的復(fù)雜性，當(dāng)強(qiáng)剩磁存在時(shí)，場(chǎng)源的磁化方向往往與地磁場(chǎng)方向大不相同，此時(shí)采用傳統(tǒng)的反演方法會(huì)得到錯(cuò)誤的結(jié)果.

近年來(lái)，針對(duì)上述問(wèn)題，國(guó)內(nèi)外學(xué)者已進(jìn)行了相當(dāng)深入的研究.在已有的研究成果中，減少反演中剩磁影響的方法主要有三種：第一種方法是首先估計(jì)磁異常的磁化方向(Lourenco and Morrison，1973；Phillips，2005；Dannemiller and Li，2006；Gerovska et al.，2009)，然后將估計(jì)的磁化方向用于三維反演，條件理想時(shí)，這種方法的反演精度較高，但不足之處是其僅適用于孤立場(chǎng)源情況.第二種方法是直接反演磁化強(qiáng)度矢量.王妙月等(2004)提出了二維層狀模型的磁化強(qiáng)度矢量反演方法；Lelièvre和Oldenburg(2009)以及Ellis等(2012)提出了直接反演三維磁化強(qiáng)度矢量的方法；劉雙等(2013)提出了基于磁異常模量的井中二維磁化強(qiáng)度矢量反演方法.這些方法可以同時(shí)反演磁化強(qiáng)度的大小和方向，但通常需要添加特定的約束來(lái)減少反演參數(shù)增多帶來(lái)的更加明顯的非唯一性(Lelièvre and Oldenburg，2009；Li S L and Li Y G，2014)，目前的應(yīng)用還很少，研究還需進(jìn)一步深化.第三種方法是首先將磁異常數(shù)據(jù)轉(zhuǎn)化為受磁化方向影響小的轉(zhuǎn)換量，然后再對(duì)這個(gè)轉(zhuǎn)換量進(jìn)行反演.Shearer(2005)和Li等(2010)通過(guò)反演磁異常模量來(lái)得到磁化強(qiáng)度大小的三維分布；Pilkington和Beiki(2013)通過(guò)反演歸一化磁源強(qiáng)度來(lái)減少剩磁的影響.這種基于轉(zhuǎn)換量的反演方法可以用于多場(chǎng)源復(fù)雜磁異常的反演.但是，進(jìn)一步的分析也發(fā)現(xiàn)，由于轉(zhuǎn)換量消除了原有磁異常中包含的相位信息，因此有時(shí)難以恢復(fù)地質(zhì)體的準(zhǔn)確產(chǎn)狀，只能反演出大概位置(Li et al.，2010)，這方面的研究工作還需要不斷深化.

本文主要的工作屬于第三種方法中磁異常模量反演的范疇.因?yàn)楸容^而言，這種方法需要最少的先驗(yàn)信息，且可以用于多場(chǎng)源復(fù)雜磁異常的反演(Li S L and Li Y G，2014).磁異常模量反演已經(jīng)成功應(yīng)用到盆地火山巖勘探(Li et al.，2012b)、金屬礦勘探(Santos et al.，2012；Ribeiro et al.，2013)以及海洋基底構(gòu)造研究(Li et al.，2012a)，具有廣泛的應(yīng)用.由于磁異常模量與地下場(chǎng)源磁化率之間是非線性關(guān)系，Shearer(2005)和Li等(2010)提出的反演方法需要通過(guò)反復(fù)求解非線性反演問(wèn)題來(lái)搜索正則化因子，這個(gè)過(guò)程非常耗時(shí)；此外，由于該算法屬于L2范數(shù)正則化反演方法，因此反演結(jié)果較為模糊.針對(duì)上述兩個(gè)問(wèn)題，我們將數(shù)據(jù)空間反演算法(Pilkington，2009)擴(kuò)展到了磁異常模量反演中，在數(shù)據(jù)空間中，無(wú)需搜索正則化因子，可大幅提高計(jì)算效率；此外由于模型范數(shù)定義為加權(quán)磁化率平方根向量的平方和，反演結(jié)果中的高磁化率部分和相鄰單元較大的磁化率變化將會(huì)被保留(Lelièvre and Oldenburg，2006)，反演結(jié)果更加聚焦，可在一定程度上改善深部分辨率.

本文首先回顧磁異常模量和數(shù)據(jù)空間反演方法的基本原理，然后再推導(dǎo)基于數(shù)據(jù)空間的磁異常模量反演方法的迭代公式，最后通過(guò)模型試驗(yàn)和實(shí)測(cè)數(shù)據(jù)驗(yàn)證本方法的有效性和計(jì)算效率.

2 數(shù)據(jù)空間磁異常模量反演方法

2.1 磁異常模量

由于實(shí)際勘探所用的磁力儀測(cè)量的是磁性體產(chǎn)生的磁異常的某個(gè)投影分量，其結(jié)果受磁化方向影響很大.前人的研究發(fā)現(xiàn)磁異常模量具有弱敏感于磁化方向和與場(chǎng)源平面位置對(duì)應(yīng)關(guān)系更好的特點(diǎn)，可以減少磁測(cè)數(shù)據(jù)解釋中的剩磁影響(Stavrev and Gerovska，2000；Shearer，2005；Li et al.，2010；Liu et al.，2013；Li S L and Li Y G，2014).在三維情況下，磁異常模量的定義如下：

(1)

其中Ta是磁性體產(chǎn)生的磁異常矢量Ta的模量，Hax、Hay、Za分別為磁異常矢量Ta在x、y、z三個(gè)方向的投影分量.

由于現(xiàn)在磁力儀(如質(zhì)子磁力儀)實(shí)測(cè)數(shù)據(jù)通常為總場(chǎng)異常ΔT(ΔT相當(dāng)于是磁異常矢量Ta在地磁場(chǎng)方向上的投影分量)，因此需要利用ΔT換算磁異常模量.對(duì)于平面網(wǎng)格數(shù)據(jù)，可通過(guò)頻率域轉(zhuǎn)換的方法(Pedersen，1978)將總場(chǎng)異常轉(zhuǎn)化為異常三分量，然后投影合成得到模量數(shù)據(jù).對(duì)于起伏地形觀測(cè)數(shù)據(jù)，常規(guī)分量轉(zhuǎn)換不易實(shí)現(xiàn)，可以采用空間域的等效源方法(Dampney，1969)，進(jìn)一步得到模量數(shù)據(jù)，然后再進(jìn)行基于模量數(shù)據(jù)的反演.由于磁異常模量與磁異常的量級(jí)一致，因此在轉(zhuǎn)換的過(guò)程中不會(huì)放大噪聲，保持了長(zhǎng)波長(zhǎng)信息.

2.2 數(shù)據(jù)空間的優(yōu)勢(shì)

在磁異常模量反演中，假設(shè)模型和數(shù)據(jù)個(gè)數(shù)分別為M和N，由于M通常是遠(yuǎn)大于N的，特別是三維反演情況.所以，該反演屬于純欠定問(wèn)題.在模型空間需要求解一個(gè)M×M的方程組，而在數(shù)據(jù)空間中只需要求解一個(gè)N×N的方程組(Tarantola，1987)，因此，數(shù)據(jù)空間反演算法將在一定程度上提高計(jì)算效率，減少計(jì)算時(shí)間.此外，由于模型空間中包含龐大的零空間(至少M(fèi)-N維)，在模型空間中求解需要施加正則化來(lái)改善系數(shù)矩陣的病態(tài)程度，而正則化因子的搜索是一個(gè)非常耗時(shí)的過(guò)程(通常是需要反復(fù)求解非線性反演的問(wèn)題)；在數(shù)據(jù)空間中則不存在零空間，因此大大提高了系數(shù)矩陣的穩(wěn)定性，同時(shí)因?yàn)闊o(wú)需搜索正則化參數(shù)，所以，計(jì)算量得到進(jìn)一步減少.為此，這里我們立足于數(shù)據(jù)空間進(jìn)行磁異常模量反演(關(guān)于模型空間和數(shù)據(jù)空間的差異，在附錄A-1中有介紹).

另外，我們經(jīng)過(guò)深入分析認(rèn)為，模型空間中改善系數(shù)矩陣的病態(tài)程度的正則化過(guò)程，也有其副作用，就是增加了反演結(jié)果圖像的模糊性.而這恰恰是基于數(shù)據(jù)空間不需要的，因而，基于數(shù)據(jù)空間的反演具有更好的模型聚焦度，也即反演結(jié)果圖像會(huì)更加清晰.

2.3 模型定義與磁化率正約束

在實(shí)際礦產(chǎn)資源勘探中，絕大多數(shù)巖石的磁化率均為正值，在反演中加入磁化率正約束可以有效地恢復(fù)地質(zhì)體的產(chǎn)狀信息，得到更有地質(zhì)意義的反演結(jié)果(Li and Oldenburg，1996).首先定義有效磁化率κ=μ0M/T0，其中，μ0=4π×10-7H·m-1為真空中的磁導(dǎo)率，M為有效磁化強(qiáng)度，單位為A·m-1，T0為地磁場(chǎng)強(qiáng)度，單位為T.在數(shù)據(jù)空間反演算法中，利用平方根算子m=κ1/2對(duì)模型進(jìn)行投影變換實(shí)現(xiàn)正約束(Lelièvre and Oldenburg，2006)是一種簡(jiǎn)單的技巧，其中m為模型向量，κ為有效磁化率向量.即模型向量可以在(-∞,∞)變化，但有效磁化率向量只能在(0,∞)之間變化.這樣做的優(yōu)勢(shì)有兩個(gè)：一是反演所需的雅可比矩陣與常規(guī)方式比較，具有一種更簡(jiǎn)單的形式；二是反演結(jié)果中的高磁化率部分和相鄰單元較大的磁化率變化將會(huì)被保留，反演結(jié)果更加聚焦(Lelièvre and Oldenburg，2006)，這個(gè)優(yōu)點(diǎn)也是我們需要重視的.

2.4 總目標(biāo)函數(shù)

磁異常模量反演和磁異常分量反演的目標(biāo)函數(shù)形式上是一樣的，帶參考模型約束的總目標(biāo)函數(shù)有如下形式：

S(m)= (d-F(m))TD-1(d-F(m))

+(m-m0)TW-1(m-m0),

(2)

2.5 目標(biāo)函數(shù)的優(yōu)化

對(duì)(2)式的求解通常是在模型空間進(jìn)行的.這里我們參考Tarantola(1987)和Pilkington(2009)的處理方式，建立起磁異常模量基于數(shù)據(jù)空間反演的目標(biāo)函數(shù)求解公式.經(jīng)過(guò)推導(dǎo)在數(shù)據(jù)空間中極小化(2)式可以得到模型m的不動(dòng)點(diǎn)迭代方程(推導(dǎo)過(guò)程見(jiàn)附錄A-1)：

m=m0+WJT(D+JWJT)-1(d-F(m)

+J(m-m0)),

(3)

其中J代表F(m)的雅可比矩陣(維數(shù)為N×M，具體形式見(jiàn)附錄A-2).由于(3)式直接對(duì)模型進(jìn)行迭代更新，其求解的穩(wěn)定性和收斂性會(huì)很差，所以，我們采取將每次迭代得到的模型m作為下次迭代的參考模型，得到第k次迭代時(shí)數(shù)據(jù)空間的近似迭代公式：

Δmk=mk+1-mk

(4)

(4)式對(duì)模型的修正量進(jìn)行迭代更新，穩(wěn)定性和收斂性較好.其中α為迭代步長(zhǎng)，實(shí)際計(jì)算中，初始選擇1，計(jì)算擬合差，若不減小，則逐次減少為原來(lái)的1/3，直至擬合差減小為止.

將(4)式改寫(xiě)為緊湊形式：

(5)

fk=Akbk,

(6)

2.6 實(shí)際應(yīng)用中需要注意的問(wèn)題

在實(shí)際應(yīng)用中，還需要注意以下幾個(gè)方面的問(wèn)題：

首先，在實(shí)測(cè)數(shù)據(jù)中，可能存在異常不完整或分布于數(shù)據(jù)體邊緣的情況，需要在模型剖分時(shí)對(duì)水平方向進(jìn)行適當(dāng)?shù)臄U(kuò)邊，數(shù)值模擬表明，這樣做不但可以減少邊界效應(yīng)，還可以加速算法的收斂.

其次，在反演迭代的過(guò)程中，建議將共軛梯度迭代的殘差設(shè)置為期望的噪聲水平(Pilkington，2009)，這種不完全共軛梯度迭代有兩個(gè)優(yōu)點(diǎn)：一是通過(guò)降低迭代的次數(shù)來(lái)減少反演時(shí)間，二是可以起到類似正則化的作用，提高反演算法的收斂性和穩(wěn)定性.

另外，因?yàn)檠趴杀染仃囉?jì)算的要求，初始模型m0不能給0，通?？梢越o一個(gè)較小的值，例如0.001.

3 模型實(shí)驗(yàn)

我們采用兩個(gè)簡(jiǎn)單模型和一個(gè)復(fù)雜模型來(lái)測(cè)試反演算法的有效性.簡(jiǎn)單模型包括立方體模型和傾斜板狀體模型，復(fù)雜模型為組合傾斜板狀體模型.這三種模型已被多位學(xué)者用于重磁異常物性反演的測(cè)試(Li and Oldenburg，1996，1998；Portniaguine and Zhdanov，2002；姚長(zhǎng)利等，2007).

觀測(cè)數(shù)據(jù)均為21×21的網(wǎng)格數(shù)據(jù)，網(wǎng)格間距均為50 m.我們首先利用最小曲率法對(duì)磁異常進(jìn)行擴(kuò)邊，然后利用頻率域方法將其轉(zhuǎn)化為模量，最后分別采用模型空間反演方法(Li et al., 2010)和數(shù)據(jù)空間反演方法對(duì)模量進(jìn)行反演.在反演中，模型剖分為立方體單元，邊長(zhǎng)與觀測(cè)數(shù)據(jù)的網(wǎng)格間距一致，為了減少邊界效應(yīng)，模型剖分為31×31×10的立方體，同時(shí)為了便于對(duì)比，只顯示中間的21×21×10的數(shù)據(jù).在模型空間反演方法中，為了減少模型的光滑效應(yīng)，我們將αx、αy、αz均設(shè)置為0，使模型保持最高的分辨率.在模型空間和數(shù)據(jù)空間反演方法中，初始模型均設(shè)置為0.001.在主頻為2.90 GHz的計(jì)算機(jī)上，利用兩種反演算法對(duì)上述三個(gè)模型進(jìn)行了反演測(cè)試.

3.1 簡(jiǎn)單模型

模型一是一個(gè)邊長(zhǎng)為200 m的立方體，立方體的中心位于地下250 m，其有效磁化率為0.05SI，地磁場(chǎng)強(qiáng)度為50000 nT，地磁場(chǎng)傾角和偏角分別為90°和0°.磁化傾角和偏角分別為60°和30°.立方體模型產(chǎn)生的磁異常如圖1a所示，其中包含了標(biāo)準(zhǔn)差為2%數(shù)據(jù)絕對(duì)值加上0.5 nT的高斯噪聲.由磁異常換算得到的模量如圖1b所示，可以看出，模量與立方體在水平面上的投影位置對(duì)應(yīng)關(guān)系良好.

圖2a、2b和2c、2d分別為采用模型空間反演算法和數(shù)據(jù)空間反演算法的反演結(jié)果切片圖.可以看出，模型空間的反演結(jié)果較為光滑，場(chǎng)源邊界模糊一些，且磁化率值整體偏低.而數(shù)據(jù)空間的反演結(jié)果聚焦程度高，場(chǎng)源邊界清晰且磁化率值更接近模型的真實(shí)磁化率.數(shù)據(jù)空間反演共需185次共軛梯度(以下簡(jiǎn)稱CG)迭代，耗時(shí)2.2 s.模型空間反演共需802次CG迭代，耗時(shí)14.3 s.

需要指出的是，由于三維磁性體磁異常模量依然受到磁化方向的影響，所以兩種方法的反演結(jié)果都略微偏離了真實(shí)模型的位置.但是，這是模量反演中在磁化方向未知的情況下得到的結(jié)果.同樣條件下，如果是磁異常的反演，則磁化方向的選擇誤差會(huì)造成很大的結(jié)果誤差，甚至是完全錯(cuò)誤的結(jié)果，這樣的情況，我們?cè)谙旅娴慕M合模型例子中得到了更明顯的體現(xiàn).

圖1 (a)立方體模型產(chǎn)生的磁異常；(b)由磁異常轉(zhuǎn)換得到的模量Fig.1 (a) The total field anomaly produced by the cube model； (b) The amplitude data computed from the total field anomaly in Fig.1a

圖2 圖1b所示的磁異常模量的反演結(jié)果，黑框代表真實(shí)模型的位置(僅顯示有效磁化率為0.0015SI以上的值) (a)和(b)利用模型空間反演算法的Z=-225 m和N=500 m的反演結(jié)果切片圖； (c)和(d)利用數(shù)據(jù)空間反演算法的Z=-225 m和N=500 m的反演結(jié)果切片圖.Fig.2 Inversion results of the amplitude data in Fig.1b. The black lines indicate the position of the true model. Only values above 0.0015SI are displayed(a，b) Model-space inversion results, slice at 225 m depth and 500 m north； (c，d) Data-space inversion results, slice at 225 m depth and 500 m north.

模型二是一個(gè)45°傾角的傾斜板狀體，其有效磁化率為0.05SI，地磁場(chǎng)傾角和偏角分別為65°和-25°，磁化傾角和偏角分別為45°和75°.傾斜板狀體模型產(chǎn)生的磁異常如圖3a所示，其中包含了標(biāo)準(zhǔn)差為2%數(shù)據(jù)絕對(duì)值加上2 nT的高斯噪聲.由磁異常換算得到的模量如圖3b所示，可以看出，模量與板狀體在水平面上的投影位置對(duì)應(yīng)關(guān)系良好.

圖4a、4b和4c、4d分別為采用模型空間反演算法和數(shù)據(jù)空間反演算法的反演結(jié)果切片圖.可以看出，模型空間的反演結(jié)果場(chǎng)源邊界模糊，產(chǎn)狀不明顯，但磁化率極大值接近真實(shí)場(chǎng)源磁化率.而數(shù)據(jù)空間的反演結(jié)果場(chǎng)源邊界清晰，產(chǎn)狀明顯，不足之處是磁化率極大值大于模型的真實(shí)磁化率，這也是聚焦性體現(xiàn)的一個(gè)附帶現(xiàn)象.

本例的數(shù)據(jù)空間反演共需89次CG迭代，耗時(shí)1.6 s.對(duì)應(yīng)的模型空間反演共需832次CG迭代，耗時(shí)15.2 s.

3.2 復(fù)雜模型

模型三為組合傾斜板狀體模型，由兩個(gè)長(zhǎng)寬、磁化率、延伸長(zhǎng)度、磁化方向不同，傾向相反但走向相同的傾斜板狀體組成，主要用于測(cè)試干擾場(chǎng)源的影響.地磁場(chǎng)傾角和偏角分別為65°和-25°.西側(cè)板狀體的有效磁化率為0.04SI，磁化傾角和偏角分別為-25°和-25°.東側(cè)板狀體的有效磁化率為0.05SI，磁化傾角和偏角分別為45°和75°.該組合模型產(chǎn)生的磁異常如圖5a所示，其中包含了標(biāo)準(zhǔn)差為2%數(shù)據(jù)絕對(duì)值加上2 nT的高斯噪聲.由磁異常換算得到的模量如圖5b所示.

圖3 (a) 傾斜板狀體模型產(chǎn)生的磁異常；(b) 由磁異常轉(zhuǎn)換得到的模量Fig.3 (a) The total field anomaly produced by the dipping slap model； (b) The amplitude data computed from the total field anomaly in Fig.3a

圖4 圖3b所示的磁異常模量的反演結(jié)果，黑框代表真實(shí)模型的位置(僅顯示有效磁化率為0.0066SI以上的值)(a)和(b)利用模型空間反演算法的Z=-225 m和N=500 m的反演結(jié)果切片圖；(c)和(d)利用數(shù)據(jù)空間反演算法的Z=-225 m和N=500 m的反演結(jié)果切片圖.Fig.4 Inversion results of the amplitude data in Fig.3b. The black lines indicate the position of the true model. Only values above 0.0066SI are displayed(a，b) Model-space inversion results, slice at 225 m depth and 500 m north； (c，d) Data-space inversion results, slice at 225 m depth and 500 m north.

對(duì)于疊加異常，為了體現(xiàn)模量反演和常規(guī)磁異常反演的差異，我們首先假設(shè)磁化方向與地磁場(chǎng)方向一致，利用標(biāo)準(zhǔn)磁異常三維反演算法對(duì)圖5a所示磁總場(chǎng)異常進(jìn)行了反演.圖6a、6b和6c、6d分別為采用模型空間反演算法(Li and Oldenburg，1996，2003)和數(shù)據(jù)空間反演算法(Pilkington，2009)的反演結(jié)果切片圖.可以看出，標(biāo)準(zhǔn)磁異常三維反演算法的反演結(jié)果幾乎是完全錯(cuò)誤的.

然后我們利用磁異常模量反演算法對(duì)圖5b所示磁異常模量進(jìn)行了反演，圖7a、7b和7c、7d分別為采用模型空間反演算法和數(shù)據(jù)空間反演算法的反演結(jié)果切片圖.可以看出，磁異常模量反演算法的反演結(jié)果較好.比較而言，模型空間的反演結(jié)果難以區(qū)分兩個(gè)板狀體，而數(shù)據(jù)空間的反演結(jié)果場(chǎng)源邊界則更加清晰.當(dāng)然，兩種反演方法均不能反演出此模型的深部構(gòu)造特征，原因主要是深部構(gòu)造產(chǎn)生的磁異常較小，且被西側(cè)板狀體產(chǎn)生的異常所掩蓋.

圖5 (a)組合傾斜板狀體模型產(chǎn)生的磁異常；(b) 由磁異常轉(zhuǎn)換得到的模量Fig.5 (a) The total field anomaly produced by the model composed of two dipping slaps；(b) The amplitude data computed from the total field anomaly in Fig.5a

圖6 圖5a所示的磁總場(chǎng)異常的反演結(jié)果，黑框代表真實(shí)模型的位置(僅顯示磁化率為0.0046SI以上的值)(a)和(b)利用模型空間反演算法的Z=-125 m和N=500 m的反演結(jié)果切片圖；(c)和(d)利用數(shù)據(jù)空間反演算法的Z=-125 m和N=500 m的反演結(jié)果切片圖.Fig.6 Inversion results of the total field anomaly data in Fig.5a. The black lines indicate the position of the true model. Only values above 0.0046SI are displayed(a，b) Model-space inversion results, slice at 125 m depth and 500 m north；(c，d) Data-space inversion results, slice at 125 m depth and 500 m north.

反演計(jì)算對(duì)比中，數(shù)據(jù)空間反演共需189次CG迭代，耗時(shí)2.5 s.模型空間反演共需1435次CG迭代，耗時(shí)20.1 s.

從上述模型的反演結(jié)果來(lái)看，數(shù)據(jù)空間反演算法會(huì)得到一個(gè)更加聚焦的反演結(jié)果，同時(shí)計(jì)算量和計(jì)算時(shí)間更少.

4 實(shí)際資料試驗(yàn)

圖8a為澳大利亞某地區(qū)的實(shí)測(cè)航磁數(shù)據(jù)，數(shù)據(jù)大小為53×71，網(wǎng)格距為50 m.該地區(qū)地磁場(chǎng)強(qiáng)度為58240 nT，地磁傾角和偏角分別為-64.7°和2.1°.

圖7 圖5b所示的磁異常模量的反演結(jié)果，黑框代表真實(shí)模型的位置(僅顯示磁化率為0.0042SI以上的值)(a)和(b)利用模型空間反演算法的Z=-125 m和N=500 m的反演結(jié)果切片圖；(c)和(d)利用數(shù)據(jù)空間反演算法的Z=-125 m和N=500 m的反演結(jié)果切片圖.Fig.7 Inversion results of the amplitude data in Fig.5b. The black lines indicate the position of the true model. Only values above 0.0042SI are displayed(a，b) Model-space inversion results, slice at 125 m depth and 500 m north;(c，d) Data-space inversion results, slice at 125 m depth and 500 m north.

圖8 (a) 實(shí)測(cè)航磁數(shù)據(jù)；(b) 由磁異常轉(zhuǎn)換得到的模量Fig.8 (a) The measured aeromagnetic data；(b) The amplitude data computed from the total field anomaly in Fig.8a

可以看出，圖中存在三個(gè)磁化方向明顯不同的疊加異常，采用單一磁化方向?qū)υ摂?shù)據(jù)進(jìn)行反演顯然是不合適的.由磁異常換算得到的模量如圖8b所示.

在模量反演中，為了減少邊界效應(yīng)，地下剖分為63×81×27的立方體，即模型的邊界范圍比采集的數(shù)據(jù)邊界范圍稍大.同時(shí)為了便于對(duì)比，僅顯示中間53×71×27的反演結(jié)果.圖9a、9b和9c、9d分別為采用模型空間反演算法和數(shù)據(jù)空間反演算法的反演結(jié)果切片圖.從圖9a和9c可以看出，兩種方法均顯示了三個(gè)場(chǎng)源的平面位置，數(shù)據(jù)空間反演的結(jié)果要更加聚焦一些.對(duì)于北側(cè)異常而言，模型空間顯示為一個(gè)場(chǎng)源，而數(shù)據(jù)空間則顯示為兩個(gè).由于沒(méi)有鉆孔數(shù)據(jù)，還無(wú)法進(jìn)行直接驗(yàn)證，但該反演結(jié)果更細(xì)致地展示了場(chǎng)源分布的細(xì)節(jié).圖9b和9d均顯示了北側(cè)場(chǎng)源的產(chǎn)狀，模型空間反演的結(jié)果模糊一些，而數(shù)據(jù)空間反演的結(jié)果則更加聚焦，產(chǎn)狀也更加明顯.

該實(shí)測(cè)數(shù)據(jù)的數(shù)據(jù)空間反演共需187次CG迭代，耗時(shí)600.7 s.對(duì)應(yīng)的模型空間反演共需2402次CG迭代，耗時(shí)4509.3 s.

圖9 實(shí)測(cè)數(shù)據(jù)的反演結(jié)果，虛線代表切片的位置(僅顯示磁化率為0.005SI以上的值)(a)和(b)利用模型空間反演算法的Z=-425 m和N=6493825 m的反演結(jié)果切片圖； (c)和(d)利用數(shù)據(jù)空間反演算法的Z=-425 m和N=6493825 m的反演結(jié)果切片圖.Fig.9 Inversion results of the amplitude data in Fig.8b. The dash lines indicate the position of the slices. Only values above 0.005SI are displayed(a，b) Model-space inversion results, slice at 425 m depth and 6493825 m north； (c，d) Data-space inversion results, slice at 425 m depth and 6493825 m north.

5 結(jié)論

強(qiáng)剩磁的存在改變了磁化強(qiáng)度的方向，進(jìn)而影響了磁異常的反演和解釋.磁異常模量是一種受磁化方向影響小的轉(zhuǎn)換量，且由磁異常轉(zhuǎn)換模量的過(guò)程中能達(dá)到不放大噪聲，保持了長(zhǎng)波長(zhǎng)信息.我們將數(shù)據(jù)空間反演算法擴(kuò)展到了磁異常模量反演中，模型試驗(yàn)表明，與模型空間反演算法相比，其具有計(jì)算量小，成像結(jié)果聚焦的優(yōu)勢(shì).但該算法同時(shí)也存在反演結(jié)果的磁化率極大值可能偏大的缺點(diǎn)，這是在實(shí)際應(yīng)用中需要注意的地方.

致謝 本文研究工作得到了加拿大地質(zhì)調(diào)查局Mark Pilkington的幫助，使用的實(shí)測(cè)數(shù)據(jù)來(lái)源于澳大利亞地質(zhì)調(diào)查局網(wǎng)站公布的資料，兩位審稿專家提出了寶貴的修改建議，在此一并感謝.

附錄A 數(shù)據(jù)空間磁異常模量反演算法

1 數(shù)據(jù)空間反演算法迭代公式的推導(dǎo)

磁性反演計(jì)算，就是對(duì)目標(biāo)函數(shù)(2)式求解極小值的過(guò)程.在目標(biāo)函數(shù)的極小值處，(2)式的梯度一定為0，故可得到：

=0,

(A1)

其中，J為F(m)的雅可比矩陣.我們對(duì)(A1)式進(jìn)一步轉(zhuǎn)化，可以得到適合計(jì)算的具體公式.由(A1)式整理得W-1(m-m0)=JTD-1(d-F(m))，兩邊加上JTD-1J(m-m0)則得到：

(JTD-1J+W-1)(m-m0)=

JTD-1(d-F(m)+J(m-m0)),

(A2)

再由(A2)式兩邊乘上(JTD-1J+W-1)-1，得到：

m=m0+(JTD-1J+W-1)-1JTD-1(d-F(m)

+J(m-m0)).

(A3)

在(A3)式中，JTD-1J+W-1為M×M階矩陣，即存在于模型空間，該矩陣的計(jì)算如求逆等，計(jì)算量通常會(huì)是非常大的.為此，設(shè)法對(duì)其進(jìn)一步轉(zhuǎn)換.

由于(JTD-1J+W-1)和JTD-1滿足下列恒等式：

JTD-1(D+JWJT) =JT+JTD-1JWJT

=(JTD-1J+W-1)WJT,

(A4)

且(D+JWJT)和(JTD-1J+W-1)均可逆，由(A4)式則可以得出：

(JTD-1J+W-1)-1JTD-1=WJT(D+JWJT)-1,

(A5)

進(jìn)一步將(A5)式代入(A3)式中進(jìn)行等式替換即得到了方程在數(shù)據(jù)空間中的不動(dòng)點(diǎn)形式：m=m0+WJT(D+JWJT)-1(d-F(m)+J(m-m0)),

(A6)

其中，D+JWJT為N×N階矩陣，即存在于數(shù)據(jù)空間.與(A3)式比較，該矩陣求逆等計(jì)算的計(jì)算量將大大減少.

2 雅可比矩陣J的推導(dǎo)

由于磁異常模量與地下模型的磁化率(或磁化強(qiáng)度)之間是非線性關(guān)系，因此需要在迭代的過(guò)程中計(jì)算其雅可比矩陣.首先將磁異常三分量的正演過(guò)程表達(dá)為：

Hax=Gxκ,

Hay=Gyκ,

Za=Gzκ,

(A7)

其中，Hax、Hay、Za分別為磁異常向量的x、y、z分量.Gx、Gy、Gz是磁異常三分量正演的靈敏度矩陣.κ=m2為有效磁化率向量.

在迭代的過(guò)程中，Ta的第i個(gè)分量Tai對(duì)第j個(gè)模型mj的偏導(dǎo)數(shù)為：

×2mj,

(A8)

由(A8)式可以得到雅可比矩陣J：

×diag(2m),

(A9)

其中，diag(Hax/Ta)表示主對(duì)角由Hax/Ta向量構(gòu)成的對(duì)角矩陣.diag(Hay/Ta)，diag(Za/Ta)以及diag(2m)的定義類似.

這樣在反演迭代的過(guò)程中，僅需存儲(chǔ)Gx，Gy，Gz三個(gè)矩陣即可快速完成相關(guān)的矩陣向量運(yùn)算.計(jì)算這三個(gè)矩陣需要正確的磁化方向，但真實(shí)的磁化方向是未知的，由于磁異常模量受磁化方向的影響小，因此在計(jì)算這三個(gè)矩陣時(shí)，可以采用地磁場(chǎng)方向作為假設(shè)的磁化方向，這是模量反演中需要的(Li et al.，2010).顯然，這個(gè)磁化方向的選擇已經(jīng)不再像常規(guī)磁異常分量反演那么重要了.

3 數(shù)據(jù)空間反演步驟

為便于了解在數(shù)據(jù)空間進(jìn)行反演的計(jì)算細(xì)節(jié)，這里將具體計(jì)算步驟列出如下：

(6)mi+1=mi+Δmi；

Dampney C N G. 1969. The equivalent source technique.Geophysics, 34(1): 39-53.Dannemiller N, Li Y G. 2006. A new method for determination of magnetization direction.Geophysics, 71(6): L69-L73.

Ellis R G, de Wet B, Macleod I N. 2012. Inversion of magnetic data from remanent and induced sources. ∥ 22nd International Geophysical Conference and Exhibition. Brisbane, Australia.

Gerovska D, Aráuzo-Bravo M J, Stavrev P. 2009. Estimating the magnetization direction of sources from southeast Bulgaria through correlation between reduced-to-the-pole and total magnitude anomalies.GeophysicalProspecting, 57(4): 491-505.

Lelièvre P G, Oldenburg D W. 2006. Magnetic forward modelling and inversion for high susceptibility.GeophysicalJournalInternational, 166(1): 76-90.

Lelièvre P G, Oldenburg D W. 2009. A 3D total magnetization inversion applicable when significant, complicated remanence is present.Geophysics, 74(3): L21-L30.

Li S L, Li Y G, Meng X H. 2012a. The 3D magnetic structure beneath the continental margin of the northeastern South China Sea.AppliedGeophysics, 9(3): 237-246.

Li S L, Li Y G. 2014. Inversion of magnetic anomaly on rugged observation surface in the presence of strong remanent magnetization.Geophysics, 79(2): J11-J19.

Li Y G, Oldenburg D W. 1996. 3-D inversion of magnetic data.Geophysics, 61(2): 394-408.

Li Y G, Oldenburg D W. 1998. 3-D inversion of gravity data.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.GeophysicalJournalInternational, 152(2): 251-265.

Li Y G, Shearer S E, Haney M M, et al. 2010. Comprehensive approaches to 3D inversion of magnetic data affected by remanent magnetization.Geophysics, 75(1): L1-L11.

Li Y G, He Z X, Liu Y X. 2012b. Application of magnetic amplitude inversion in exploration for volcanic units in a basin environment.Geophysics, 77(5): B219-B225.

Liu S, Feng J, Gao W L, et al. 2013. 2D inversion for borehole magnetic data in the presence of significant remanence and demagnetization.ChineseJ.Geophys. (in Chinese), 56(12): 4297-4309, doi: 10.6038/cjg20131232.

Liu S, Hu X Y, Liu T Y, et al. 2013. Magnetization vector imaging for borehole magnetic data based on magnitude magnetic anomaly.Geophysics, 78(6): D429-D444.

Lourenco J S, Morrison H F. 1973. Vector magnetic anomalies derived from measurements of a single component of the field.Geophysics, 38(2): 359-368.

Nocedal J, Wright S. 2006. Numerical Optimization. New York: Springer.

Pedersen L B. 1978. Wavenumber domain expressions for potential fields from arbitrary 2-, 21/2-, and 3-dimensional bodies.

Geophysics, 43(3): 626-630.

Phillips J D. 2005. Can we estimate total magnetization directions from aeromagnetic data using Helbig′s integrals?Earth,PlanetsandSpace, 57(8): 681-689. Pilkington M. 1997. 3-D magnetic imaging using conjugate gradients.Geophysics, 62(4): 1132-1142. Pilkington M. 2009. 3D magnetic data-space inversion with sparseness constraints.Geophysics, 74(1): L7-L15. Pilkington M, Beiki M. 2013. Mitigating remanent magnetization effects in magnetic data using the normalized source strength.Geophysics, 78(3): J25-J32.

Portniaguine O, Zhdanov M S. 2002. 3D magnetic inversion with data compression and image focusing.Geophysics, 67(5): 1532-1541.

Santos M L, Li Y G, Moraes R. 2012. Application of 3D magnetic amplitude inversion to Fe oxide-Cu-Au deposits at low magnetic latitudes: A case study from Carajás Mineral Province, Brazil.∥ 2012 SEG Annual Meeting. Society of Exploration Geophysicists. 1-5.

Shearer S E. 2005. Three-dimensional inversion of magnetic data in the presence of remanent magnetization [Master′s thesis]. Colorado: Colorado School of Mines.

Stavrev P, Gerovska D. 2000. Magnetic field transforms with low sensitivity to the direction of source magnetization and high centricity.GeophysicalProspecting, 48(2): 317-340.

Tarantola A. 1987. Inverse Problem Theory. Methods for Data Fitting and Model Parameter Estimation. New York: Elsevier.

Wang M Y, Di Q Y, Xu K, et al. 2004. Magnetization vector inversion equations and 2D forward and inversed model study.ChineseJ.Geophys. (in Chinese), 47(3): 528-534, doi: 10.3321/j.issn:0001-5733.2004.03.025.

Yao C L, Hao T Y, Guan Z N, et al. 2003. High-speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms.ChineseJ.Geophys. (in Chinese), 46(2): 252-258.

Yao C L, Zheng Y M, Zhang Y W. 2007. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces.ChineseJ.Geophys. (in Chinese), 50(5): 1576-1583.

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

劉雙, 馮杰, 高文利等. 2013. 強(qiáng)剩磁強(qiáng)退磁條件下的二維井中磁測(cè)反演. 地球物理學(xué)報(bào), 56(12): 4297-4309, doi: 10.6038/cjg20131232.

王妙月, 底青云, 許琨等. 2004. 磁化強(qiáng)度矢量反演方程及二維模型正反演研究. 地球物理學(xué)報(bào), 47(3): 528-534, doi: 10.3321/j.issn:0001-5733.2004.03.025.

姚長(zhǎng)利, 郝天珧, 管志寧等. 2003. 重磁遺傳算法三維反演中高速計(jì)算及有效存儲(chǔ)方法技術(shù). 地球物理學(xué)報(bào), 46(2): 252-258.

姚長(zhǎng)利, 鄭元滿, 張聿文. 2007. 重磁異常三維物性反演隨機(jī)子域法方法技術(shù). 地球物理學(xué)報(bào), 50(5): 1576-1583.

(本文編輯 何燕)

3D data-space inversion of magnetic amplitude data

LI Ze-Lin, YAO Chang-Li*, ZHENG Yuan-Man, MENG Xiao-Hong, ZHANG Yu-Wen

KeyLaboratoryofGeo-detection(ChinaUniversityofGeosciences,Beijing),MinistryofEducation,StateKeyLaboratoryofGeologicalProcessesandMineralResources,Beijing100083,China

3D magnetic inversion for susceptibility distribution is a powerful tool in quantitative interpretation of magnetic data and has been successfully applied to exploration of mineral and oil gas resources and to interpretation of regional geologic structure. The traditional inversion algorithms require knowledge of magnetization direction, which means that one should assume there is no remanent magnetization and self-demagnetization. Consequently, the direction of magnetization is assumed to be the same as the current geomagnetic field direction. However, strong remanent magnetization always exists and the total magnetization direction can be significantly different from that of the geomagnetic field direction due to complicated geological conditions, and in this case the traditional inversion algorithms become ineffective and the inversion result may be wrong.Magnetic amplitude data are less sensitive to the total magnetization direction and can be used to invert for 3D magnetization strength distribution to delineate the positions of causative bodies in the presence of strong remanent magnetization. We present a data-space magnetic amplitude inversion algorithm to reduce the effects of remanent magnetization. We also give a detail formula derivation of the proposed algorithm. In the data space, the matrix dimensions are equal toN×N(the number of data) rather thanM×M(the number of model parameters), whereNis far less thanM. So the computational efficiency is improved. The computational time is further reduced because this method does not need to search for a regularization parameter by using an incomplete conjugate gradient method. Moreover, a square root operator is used to impose a positivity constraint on the effective susceptibility and likewise to focus the inversion results.Three synthetic data examples and a real data example are used to verify the proposed data-space algorithm. Tests on these examples show that the proposed method can focus the inversion results and likewise increase solution speed by up to an order of magnitude compared with the standard model-space, L2-norm regularized inversion.

Data-space; Magnetic amplitude data; Remanence; 3D inversion

10.6038/cjg20151030.

Li Z L, Yao C L, Zheng Y M, et al. 2015. 3D data-space inversion of magnetic amplitude data.ChineseJ.Geophys. (in Chinese),58(10):3804-3814,doi:10.6038/cjg20151030.

國(guó)家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)項(xiàng)目(2014AA06A613)和國(guó)家自然科學(xué)基金項(xiàng)目(41304104)資助．

李澤林，主要從事重磁處理反演研究. E-mail：zelin.lee@foxmail.com

*通訊作者 姚長(zhǎng)利，從事重磁勘探理論與方法技術(shù)研究.E-mail：clyao@cugb.edu.cn

10.6038/cjg20151030

P631

2014-10-24，2015-03-05收修定稿

≤≥? ?李澤林, 姚長(zhǎng)利, 鄭元滿等. 2015. 數(shù)據(jù)空間磁異常模量三維反演.地球物理學(xué)報(bào)，58(10):3804-3814,