基于方向tilt-Euler的三維磁數據快速反演

羅瀟， 王彥國， 葛坤朋， 鄧居智， 楊亞新

1 東華理工大學地球物理與測控技術學院， 南昌 330013 2 東華理工大學核資源與環(huán)境國家重點實驗室， 南昌 330013 3 核工業(yè)230研究所， 長沙 410011

0 引言

磁法數據處理與解釋是磁法勘探工作的核心環(huán)節(jié)，選取適當的方法技術是獲得良好地質解釋的關鍵.受磁化方向影響，磁異常往往與地質體分布之間存在一定偏差，直接進行異常解釋難度較大.化極處理可以將感應磁化磁異常轉換為垂直磁化異常，在一定程度上有利于異常的解釋和推斷.另外，許多數據處理方法建立在化極磁異常基礎上，如歸一化標準差(Cooper and Cowan, 2008)和 Tilt-depth法(Salem et al., 2007).然而，當磁化角度變化較大，研究區(qū)位于低緯度，或工區(qū)存在強剩磁時，化極很難獲得滿意的效果.解析信號是磁異常解釋的一類有效方法(Bastani and Pedersen, 2001; Srivastava and Agarwal, 2010; Cooper, 2015)，二維解析信號不受磁化方向影響(Atchuta et al., 1981)，三維磁異常解析信號受磁化方向影響較小，但仍不能無視斜磁化的影響(Huang and Guan, 1998; Li, 2006).Stavrev和Gerovsk(2000)在磁三分量和梯度張量基礎上給出了R、L、E、Q四個模量轉換量，并證實了這些轉換量在三維中對磁化方向敏感度較小，二維中不受磁化方向影響.Verduzco等(2004)在tilt梯度的基礎上提出了tilt梯度水平導數模，同樣在二維磁異常處理中不受斜磁化影響，在三維中受磁化角度影響較小，也被廣泛運用于磁異常識別中，不過該方法使用了磁異常二階導數，易受噪聲干擾影響，同時還會產生假異常.Wijns等(2005)提出了Theta map，這種方法可以有效識別低緯度磁異常的邊界信息，但在中高緯度地區(qū)的應用效果較差.梯度全張量能夠更好地描述場源的異常特征，近年來已被廣泛應用于磁法數據處理中(Oru?, 2010; Beiki et al., 2011; 馬國慶等，2012a)，其中張量不變量(Schmidt and Clark,2006)受磁化方向影響較小，但也依舊無法消除斜磁化的影響；歸一化磁源強度(Beiki et al., 2012)同樣是基于梯度張量的，該方法受磁化方向影響的程度更小，已被廣泛重視與使用(Pilkington and Beiki, 2013; 饒椿鋒等，2017).

磁異常的快速反演可以自動、快捷地獲取場源位置信息，在磁異常解釋中具有重要意義.歐拉反褶積(Thompson, 1982)是一種常用的方法，因其不需要先驗信息、計算效率高而得到廣泛的研究與應用(張季生等, 2011; Guo et al., 2014; 王林飛, 2016).Salem 和 Ravat(2003)基于歐拉反褶積和解析信號，提出了AN-EUL法，用于估計場源深度及構造指數，但該方法使用了三階導數，更易受干擾影響.Salem等(2008)將tilt梯度和歐拉反褶積相結合提出了tilt-Euler法，該方法不僅削弱了斜磁化對反演結果的影響，還無需事先給定構造指數，推動了歐拉反褶積的實用化發(fā)展.Zhang等(2000)提出的張量歐拉反褶積在不提高導數階次情況下，擴充了常規(guī)歐拉反褶積方程組個數，使得反演效果更佳，但該方法需要事先已知場源構造指數.馬國慶等(2012b)提出了張量局部波數法的歐拉反褶積，實驗顯示該方法明顯優(yōu)于常規(guī)局部波數法.陳國強和馬國慶(2016)在Theta圖和歐拉反褶積方法基礎上提出了Theta-Depth法，可有效反演場源的位置和深度，但這種方法需要事先對原始磁異常進行化極處理.

本文在方向tilt梯度(王彥國等, 2019)基礎上，推導出了兩個方向tilt梯度的不同方向導數，并同時將方向tilt梯度的導數與三個方向導數形式下的歐拉反褶積方程組相結合，推導出了方向tilt-Euler方程組，用于反演磁源參數.最后，通過模型實驗和實例應用證實了該方法在斜磁化磁異常解釋中的有效性及優(yōu)越性.

1 基本原理

1.1 方向tilt梯度及其導數

Tilt梯度(Miller and Singh, 1994)是一種可以均衡不同異常強度的位場數據處理方法，其理論公式為

(1)

方向tilt梯度(王彥國等，2019)的表達式為：

(2)

(3)

θx、θy分別為x，y方向上的方向tilt梯度.

對公式(2)、(3)分別求x，y，z三個方向的導數，并對分母進行均方根處理，則得：

(4)

公式(4)中的每一個表達式單位為nT·m-2或nT·km-2，與該式中使用的導數階次相一致.王彥國等(2019)指出，利用方向tilt梯度水平導數模可以進行磁異常識別，其表達式為：

(5)

1.2 方向tilt-Euler法

三維歐拉反褶積(Reid et al., 1990)的公式為：

(6)

其中，x，y，z是觀測點的坐標，x0，y0，z0是場源的位置坐標，N是構造指數.公式(6)對x，y，z三個方向求導，得：

(7)

(8)

(9)

(10)

(11)

將公式(4)中的6個表達式代入公式(10)和(11)，得：

(12)

公式(12)即為方向tilt-Euler反演方程組.給定一個滑動計算窗口D(文中均選為5)，通過求解超定方程組獲得窗口下的場源位置反演解(x0,y0,z0).獲得場源位置(x0,y0,z0)后，將公式(7)、(8)和(9)兩邊取平方，相加后取平方根，則得構造指數計算公式為：

(13)

通過逐步滑動窗口直至完成全區(qū)的數據覆蓋，則可以獲得大量的反演解.然而，遠離場源位置的反演解是不可靠的，需要進行有效的剔除處理.由于方向tilt梯度水平導數模THSθ極大值可作為磁性體的有效水平位置，因此可以剔除離THSθ極大值較遠的反演解(文中刪除了偏離THSθ極大值大于計算窗口寬度D的解)，另外，還需要剔除反演深度小于0及構造指數大于4的反演解.

2 模型試驗

2.1 單一模型

圖1 單一長方體模型磁異常Fig.1 Magnetic anomaly of a single cuboid model

2.2 組合模型

為了驗證本文方法對復雜磁異常的處理能力，設計了一個由球體、巖脈體、正方棱柱體、薄板及長方體5個不同類型模型構成的組合模型，各模型體具有不同的磁化強度、磁化方向、埋深及水平尺度等參數，具體參數見表1.圖4a為組合模型的原始磁異常，其中網格間距為0.1 km×0.1 km，可以看出，受磁化方向影響，磁異常與地質體并無明顯對應關系.圖4b是tilt梯度的水平導數模，可以看出，該方法能夠較好地識別出巖脈②的位置、長方體⑤的邊界位置及正方棱柱體③部分邊界，但是在球體①、巖脈②及薄板④附近存在著假異常.另外，極大值出現在薄板④邊界的兩側.圖4c是方向tilt梯度水平導數模，可以看出，該方法可以很好地反映出所有模型體的位置，且不存在多余異常.

表1 疊加模型參數Table 1 The parameters of the multisource model

圖2 單一長方體模型方向tilt梯度的改進導數Fig.2 The improved derivatives of the directional magnetic tilt-gradient of the single cuboid model

圖3 單一長方體模型磁異常的方向tilt-Euler反演結果(a) 方向tilt梯度水平導數模及反演解水平位置； (b) 深度反演解直方統(tǒng)計圖； (c) 構造指數反演解直方統(tǒng)計圖.Fig.3 The inversion results of magnetic anomaly generated by the single cuboid model using the directional tilt-Euler method(a) The total horizontal derivative of the directional tilt-gradient and the horizontal locations of inversion solutions; (b) Histogram of the estimated depth solutions; (c) Histogram of the estimated structural index solutions.

圖4 組合模型磁異常及識別結果(a) 理論磁異常； (b) Tilt梯度的水平導數模； (c) 方向tilt梯度的水平導數模.Fig.4 Magnetic anomaly generated by multisource model and the recognized results(a) Forward magnetic anomaly; (b) The total horizontal derivative of the tilt gradient; (c) The total horizontal derivative of the directional tilt-gradient.

圖5 Tilt-Euler及方向tilt-Euler反演解平面位置(a) Tilt-Euler; (b) 方向tilt-Euler.Fig.5 The horizontal locations of inversion solutions by using the tilt-Euler and the directional tilt-Euler methods(a) The tilt-Euler; (b) The directional tilt-Euler.

圖6 Tilt-Euler法反演結果(a1), (b1), (c1), (d1), (e1) 分別是模型體1—5的深度反演解直方圖；(a2), (b2), (c2), (d2), (e2) 分別是模型體1—5的構造指數反演解直方圖；(f1), (g1), (h1) 是三個虛假源S1—S3的深度反演解直方圖；(f2), (g2), (h2) 是三個虛假源S1—S3的構造指數反演解直方圖.Fig.6 Inversion results of the tilt-Euler method(a1), (b1), (c1), (d1), (e1) Histograms of the estimated depths of the models 1—5; (a2), (b2), (c2), (d2), (e2) Histograms of the estimated structural indices of the models 1—5; (f1), (g1), (h1) Histograms of the estimated depths of the three spurious sources; (f2), (g2), (h2) Histograms of the estimated structural indices of the three spurious sources.

圖7 方向tilt-Euler法反演結果(a1), (b1), (c1), (d1), (e1) 分別是模型體1—5的深度反演解直方圖；(a2), (b2), (c2), (d2), (e2) 分別是模型體1-5的構造指數反演解直方圖.Fig.7 Inversion results of the directional tilt-Euler method(a1), (b1), (c1), (d1), (e1) Histograms of the estimated depths of the models 1—5; (a2), (b2), (c2), (d2), (e2) Histograms of the estimated structural indices of the models 1—5.

圖8 含噪磁異常及識別結果(a) 含0.5%噪聲的磁異常； (b) 向上延拓2倍點距的tilt angle水平導數模； (c) 向上延拓2倍點距的方向tilt梯度水平導數模.Fig.8 The noisy magnetic anomaly and its recognized results(a) Magnetic anomaly added by 0.5% random noise; (b) and (c) The total horizontal derivatives of the tilt gradient and the directional tilt gradient, but the data smoothed by upward continuation of 2 intervals.

圖9 含噪磁異常向上延拓2倍點距后的tilt-Euler及方向tilt-Euler反演解平面位置(a) Tilt-Euler; (b) 方向tilt-Euler.Fig.9 The horizontal locations of inversion solutions by using the tilt-Euler and the directional tilt-Euler methods for the noisy magnetic anomaly, but the data smooth by upward continuation of 2 intervals(a) The tilt-Euler; (b) The directional tilt-Euler.

圖10 含噪磁異常的tilt-Euler法反演結果(a1), (b1), (c1), (d1) ,(e1) 分別是模型體1—5的深度反演解直方圖；(a2), (b2), (c2), (d2), (e2) 分別是模型體1—5的構造指數反演解直方圖；(f1), (g1)是兩個虛假源S1—S2的深度反演解直方圖；(f2), (g2)是兩個虛假源S1—S2的構造指數反演解直方圖.Fig.10 Inversion results of the tilt-Euler method for the noisy magnetic anomaly(a1), (b1), (c1), (d1), (e1) Histograms of the estimated depths of the models 1—5; (a2), (b2), (c2), (d2), (e2) Histograms of the estimated structural indices of the models 1—5; (f1), (g1) Histograms of the estimated depths of the two spurious sources; (f2), (g2) Histograms of the estimated structural indices of the two spurious sources.

圖11 含噪磁異常的方向tilt-Euler法反演結果(a1), (b1), (c1), (d1), (e1) 分別是模型體1-5的深度反演解直方圖；(a2), (b2), (c2), (d2), (e2) 分別是模型體1-5的構造指數反演解直方圖.Fig.11 Inversion results of the directional tilt-Euler method for the noisy magnetic anomaly (a1), (b1), (c1), (d1), (e1) Histograms of the estimated depths of the models 1-5; (a2), (b2), (c2), (d2), (e2) Histograms of the estimated structural indices of the models 1-5.

圖12 內蒙古塔木素地區(qū)地質圖Fig.12 Geological map of Tamusu area in Inner Mongolia

圖13 塔木素航磁異常及不同方法的異常識別結果(a) 航磁異常； (b) 解析信號振幅； (c) Tilt梯度水平導數模； (d) 方向tilt梯度水平導數模.Fig.13 Aeromagnetic anomaly of Tamusu area and the recognized results using different methods(a) Aeromagnetic anomaly; (b) Analytic signal amplitude; (c) The total horizontal derivative of tilt gradient; (d) The total horizontal derivative of the directional tilt-gradient.

圖14 塔木素地區(qū)航磁數據tilt-Euler及方向tilt-Euler法反演結果(a) Tilt-Euler法反演位置解； (b) 方向tilt-Euler法反演位置解； (c) Tilt-Euler反演構造指數解； (d) 方向tilt-Euler反演構造指數解.Fig.14 The inversion results of the tilt-Euler method and the directional tilt-Euler method for aeromagnetic data of Tamusu area(a) The inversion locations of the tilt-Euler method; (b) The inversion locations of the directional tilt-Euler method; (c) The inversion structural indices of the tilt-Euler method; (d) The inversion structural indices of the directional tilt-Euler method.

圖5是tilt-Euler及方向tilt-Euler法反演解的平面分布圖，從圖中可以看出，tilt-Euler及方向tilt-Euler法都能夠很好地反演出所有模型體的位置，但是tilt-Euler的反演解連續(xù)性及聚集度都沒有方向tilt-Euler的強，且在球體下方(S1)、巖脈左上方(S2)及薄板右方(S3)存在三組虛假反演解.圖6和圖7分別是場源深度及構造指數的tilt-Euler法和方向tilt-Euler反演解直方分布圖，可以看出，兩者都很好地反演出了所有模型的深度及構造指數，不過方向tilt-Euler法精度更高，其深度及構造指數的均方誤差分別為0.05 km和0.20，小于tilt-Euler法的0.09 km和0.32.且方向tilt-Euler法的反演解更接近正態(tài)分布，解的偏差也更小.另外，tilt-Euler法還獲得了埋深分別為1.21 km，0.19 km和0.72 km的三個假源，這顯然會給解釋工作帶來一定的干擾.

為了了解噪聲干擾對反演結果的影響，對疊加磁異常(圖4a)添加了0.5%的隨機噪聲(圖8a)，圖8b和8c分別是tilt梯度及方向tilt梯度的水平導數模，不過對含噪磁數據進行了2倍點距的向上延拓處理.從圖8b和8c可以看出，tilt梯度的水平導數模更易受噪聲影響，存在明顯的數據波動，另外巖脈模型②左上側的多余異常與巖脈異?；具B在一起，易被看作是有效異常；正方棱柱體③與長方體⑤之間存在一條明顯的假異常.方向tilt梯度的水平導數模THSθ仍能夠很好地識別出所有模型體的位置，不過巖脈和薄板的異常較為模糊.

圖9a和9b分別是tilt-Euler及方向tilt-Euler反演解的平面位置展布圖，對比可以看出，tilt-Euler反演解分布較為發(fā)散，尤其巖脈和正方棱柱體邊界位置；在薄板體左邊界，反演解呈現為兩條帶狀分布，分別位于薄板邊界的兩側；在巖脈左上側、正方棱柱體與長方體之間均存在一個條帶狀分布反演解，這也易給解釋工作帶來麻煩.而方向tilt-Euler反演解具有更強的聚集度，反演解也更接近真實位置，且不存在虛假解.

圖10與圖11分別是含噪聲磁異常的tilt-Euler及方向tilt-Euler反演磁源深度及構造指數解的直方分布圖，可以看出，方向tilt-Euler反演解的分布區(qū)間較小，解的均方差更小，也更接近于正態(tài)分布.方向tilt-Euler反演深度及構造指數的均方誤差分別為0.11 km 和0.43，也小于tilt-Euler的0.24 km 和0.47，尤其深度反演精度是tilt-Euler的2倍.

Tilt-Euler法還反演出了深度分別為0.05±0.06 km和0.35±0.31 km的兩個假源.

該試驗結果表明了方向tilt-Euler法相對于tilt-Euler法，具有反演解聚集度高、精度高、準確度高，及受噪聲干擾影響小和無虛假解等優(yōu)點.

3 實際資料應用

為了檢驗實際資料應用效果，選取了內蒙古塔木素地區(qū)的航磁數據進行測試.圖12為研究區(qū)的地質簡圖，工區(qū)西北側以火山巖為主，其他地區(qū)則主要是沉積巖分布，局部有火山巖出露，研究區(qū)內主構造為近NE向.圖中紅色實線為已知斷裂構造、紅色虛線為前人利用地震資料推斷的構造，紅色箭頭方向為地震資料推斷的構造傾向.

圖13a為塔木素地區(qū)航磁異常(網格間距1 km×1 km)，研究區(qū)北部及東南地區(qū)以高磁異常為主，磁異常高低與磁性高低的地質體分布并無明顯對應關系.圖13b、13c和13d分別是磁異常的解析信號、tilt梯度水平導數模和方向tilt梯度水平導數模.解析信號(圖13b)整體表現為西北高和中部及東南部低的異常特征，與地質圖中的火山巖和沉積巖分布均有良好的對應性，中部及東南部還存在著局部高值區(qū)，可能是含有火山巖成分的沉積巖或隱伏火山巖引起的，但解析信號對弱信號的識別能力較差，異常細節(jié)特征無法清晰識別.Tilt梯度水平導數模(圖13c)獲得了豐富的異常細節(jié)信息，異常圈閉走向與地質單元分布及斷裂構造走向具有一定的相關性，但該方法的計算結果中存在虛假信息(王彥國等，2019).方向tilt梯度水平導數模同樣獲得了較為豐富的異常細節(jié)信息，看起來像是對解析信號(圖13b)細節(jié)特征的突出.另外，在沉積巖分布區(qū)內，斷裂構造附近基本上都有方向tilt梯度水平總導數模的極大值響應，且這些極值異常恰位于斷裂構造傾向方向.

圖14是tilt-Euler法和方向tilt-Euler法的反演結果.從整體上來看，兩種方法計算得到的場源位置及構造指數均較為接近，反演解的水平位置基本上為條帶狀，反演深度值大多在0～2000 m之間，構造指數則主要分布在-0.6～1.4之間，尤其沉積巖分布區(qū)內的構造指數大部分在0值附近.這些信息表明了沉積巖分布區(qū)的磁異常主要是由斷裂構造產生的.從反演效果來看，方向tilt-Euler法的反演解個數更多、聚集度更強、連續(xù)性更好，更能有效地追蹤斷裂構造的延伸.如圖14中箭頭所指位置，方向tilt-Euler法反演解連續(xù)，條帶狀明顯，而tilt-Euler法反演解較為發(fā)散、連續(xù)性較差.另外，可根據方向tilt梯度水平導數模和方向tilt-Euler法反演結果，推斷隱伏斷裂構造及地下磁性體分布，為研究該地區(qū)地下地質結構特征提供了一種依據.如圖14b、d中勾勒出了部分(隱伏)斷裂帶和一個隱伏巖體，其中斷裂帶①基本上是巖漿巖與沉積巖的界限斷裂，應是東北向地表出露及前人推斷斷裂帶在西南側的延伸；斷裂帶②則將西南與東北側斷裂連接在一起；斷裂帶③推測是東側斷裂帶的向西延伸；而斷裂帶④、⑤、⑥推測為隱伏構造帶，這些斷裂帶在方向tilt梯度水平導數模(圖13d)及方向tilt-Euler(圖14b)都有著清晰的展示，而常規(guī)的tilt-Euler反演解(圖14a)偏少，連續(xù)性較差，不利于斷裂帶的精細劃分.在圈定的隱伏巖體Ⅰ范圍內，地表有花崗巖、閃長巖的零星分布，而解析信號(圖13b)和方向tilt梯度水平總梯度模(圖13d)在該區(qū)域存在大面積的強信號，方向tilt-Euler反演的深度較深及構造指數較大，故可推測地下存在大規(guī)模的隱伏巖體.

4 結論

在方向tilt梯度基礎上，推導出了方向tilt梯度的6個導數，構建了基于方向tilt梯度的歐拉反褶積(方向tilt-Euler法)，用于三維磁異常解釋工作.模型試驗表明，相對于tilt梯度水平導數模來說，方向tilt梯度水平導數模可以更好地反映不同磁化方向的磁源位置，且不會產生多余異常；而方向tilt-Euler法反演解的連續(xù)性、匯聚度及精度也均高于tilt-Euler法的.在內蒙古塔木素地區(qū)航磁資料應用中，方向tilt梯度水平總梯度模及方向tilt-Euler法獲得了豐富的異常信息和地下磁性體參數，反演信息揭示了沉積巖分布區(qū)內磁異常主要源自斷裂構造，且解釋結果與已知地質資料和地震解釋結果吻合較好，同時反演結果也為地下磁性體分布特征研究提供了依據.

致謝衷心感謝兩位匿名評審專家為本文提供的寶貴建議.