李澤林, 姚長利, 鄭元滿, 孟小紅, 張聿文

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

?

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

李澤林, 姚長利*, 鄭元滿, 孟小紅, 張聿文

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

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

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

1 引言

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

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

本文主要的工作屬于第三種方法中磁異常模量反演的范疇.因為比較而言，這種方法需要最少的先驗信息，且可以用于多場源復(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)用.由于磁異常模量與地下場源磁化率之間是非線性關(guān)系，Shearer(2005)和Li等(2010)提出的反演方法需要通過反復(fù)求解非線性反演問題來搜索正則化因子，這個過程非常耗時；此外，由于該算法屬于L2范數(shù)正則化反演方法，因此反演結(jié)果較為模糊.針對上述兩個問題，我們將數(shù)據(jù)空間反演算法(Pilkington，2009)擴展到了磁異常模量反演中，在數(shù)據(jù)空間中，無需搜索正則化因子，可大幅提高計算效率；此外由于模型范數(shù)定義為加權(quán)磁化率平方根向量的平方和，反演結(jié)果中的高磁化率部分和相鄰單元較大的磁化率變化將會被保留(Lelièvre and Oldenburg，2006)，反演結(jié)果更加聚焦，可在一定程度上改善深部分辨率.

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

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

2.1 磁異常模量

由于實際勘探所用的磁力儀測量的是磁性體產(chǎn)生的磁異常的某個投影分量，其結(jié)果受磁化方向影響很大.前人的研究發(fā)現(xiàn)磁異常模量具有弱敏感于磁化方向和與場源平面位置對應(yīng)關(guān)系更好的特點，可以減少磁測數(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三個方向的投影分量.

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

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

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

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

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

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

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)化

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

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

+J(m-m0)),

(3)

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

Δmk=mk+1-mk

(4)

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

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

(5)

fk=Akbk,

(6)

2.6 實際應(yīng)用中需要注意的問題

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

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

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

另外，因為雅可比矩陣計算的要求，初始模型m0不能給0，通常可以給一個較小的值，例如0.001.

3 模型實驗

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

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

3.1 簡單模型

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

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

需要指出的是，由于三維磁性體磁異常模量依然受到磁化方向的影響，所以兩種方法的反演結(jié)果都略微偏離了真實模型的位置.但是，這是模量反演中在磁化方向未知的情況下得到的結(jié)果.同樣條件下，如果是磁異常的反演，則磁化方向的選擇誤差會造成很大的結(jié)果誤差，甚至是完全錯誤的結(jié)果，這樣的情況，我們在下面的組合模型例子中得到了更明顯的體現(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é)果，黑框代表真實模型的位置(僅顯示有效磁化率為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.

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

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

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

3.2 復(fù)雜模型

模型三為組合傾斜板狀體模型，由兩個長寬、磁化率、延伸長度、磁化方向不同，傾向相反但走向相同的傾斜板狀體組成，主要用于測試干擾場源的影響.地磁場傾角和偏角分別為65°和-25°.西側(cè)板狀體的有效磁化率為0.04SI，磁化傾角和偏角分別為-25°和-25°.東側(cè)板狀體的有效磁化率為0.05SI，磁化傾角和偏角分別為45°和75°.該組合模型產(chǎn)生的磁異常如圖5a所示，其中包含了標(biāo)準(zhǔn)差為2%數(shù)據(jù)絕對值加上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é)果，黑框代表真實模型的位置(僅顯示有效磁化率為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.

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

然后我們利用磁異常模量反演算法對圖5b所示磁異常模量進行了反演，圖7a、7b和7c、7d分別為采用模型空間反演算法和數(shù)據(jù)空間反演算法的反演結(jié)果切片圖.可以看出，磁異常模量反演算法的反演結(jié)果較好.比較而言，模型空間的反演結(jié)果難以區(qū)分兩個板狀體，而數(shù)據(jù)空間的反演結(jié)果場源邊界則更加清晰.當(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所示的磁總場異常的反演結(jié)果，黑框代表真實模型的位置(僅顯示磁化率為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.

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

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

4 實際資料試驗

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

圖7 圖5b所示的磁異常模量的反演結(jié)果，黑框代表真實模型的位置(僅顯示磁化率為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ù)據(jù)；(b) 由磁異常轉(zhuǎn)換得到的模量Fig.8 (a) The measured aeromagnetic data；(b) The amplitude data computed from the total field anomaly in Fig.8a

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

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

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

圖9 實測數(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é)論

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

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

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

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

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

=0,

(A1)

其中，J為F(m)的雅可比矩陣.我們對(A1)式進一步轉(zhuǎn)化，可以得到適合計算的具體公式.由(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階矩陣，即存在于模型空間，該矩陣的計算如求逆等，計算量通常會是非常大的.為此，設(shè)法對其進一步轉(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)

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

(A6)

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

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

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

Hax=Gxκ,

Hay=Gyκ,

Za=Gzκ,

(A7)

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

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

×2mj,

(A8)

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

×diag(2m),

(A9)

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

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

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

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

(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.

附中文參考文獻

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

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

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

姚長利, 鄭元滿, 張聿文. 2007. 重磁異常三維物性反演隨機子域法方法技術(shù). 地球物理學(xué)報, 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.

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

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

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

10.6038/cjg20151030

P631

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

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