重磁(梯度)張量數據邊界識別方法研究

趙建宇， 明彥伯， 孫成城， 張志東， 劉曉甲

(吉林大學 地球探測科學與技術學院，長春 130021)

重磁(梯度)張量數據邊界識別方法研究

趙建宇， 明彥伯， 孫成城， 張志東， 劉曉甲

(吉林大學 地球探測科學與技術學院，長春 130021)

邊界識別是重力數據解釋的常規任務之一，但地質體的邊界不能很好地與重力原始異常數據對應，而與重磁異常水平導數極大值、垂直導數零值相對應，因此大多利用該性質完成重磁邊界識別任務。研究發現，現有邊界識別方法會出現虛假邊界，針對現有邊界識別方法進行總結，并針對多余邊界問題進行改進，提出新的邊界識別，通過模型試驗和實際數據證明，該方法可有效地去除多余邊界，且降低了噪聲的干擾，為重磁數據邊界識別方法提供新的思路。

邊界識別； 重磁(梯度)張量數據； 重磁勘探

0 引言

重磁勘探是一種可以有效地圈定異常，劃分構造的物理勘探方法，由于它的經濟、快速、范圍廣等特點而得到廣泛應用。但是在實際生產工作中能夠發現，重磁原始異常數據和地質體的邊界不能很好地對應[1]，在異常邊界識別過程中可能會出現虛假的異常，給研究和生產造成不必要的麻煩。于是學者們開始對重磁異常數據的邊界識別進行研究，垂直導數[2-3]是最早被用來進行邊界識別的方法，后來又發現水平導數的極值位于密度和磁化率變化的位置[4-5]，從這點出發，同樣證明了水平導數對于邊界識別的作用。但是由于導數隨地質體埋深增加，其數值會迅速衰減[6]，因此，垂直導數和水平導數不能很好地識別深部異常。為了解決這一問題，學者將目光轉向對均衡濾波器的研究與開發，Miller等[7]提出了第一個均衡邊界識別濾波器—tilt angle；Verduzco等[8]提出利用tilt angle的總水平導數(THDR)進行邊界的識別；Wijns等[9]利用總水平導數與解析信號的比值(theta map)來進行邊界識別操作。這些方法有的可以均衡強弱異常之間的幅值，有的可以進行邊界識別。

重磁張量的研究是近年來興起的一種高精度勘探手段，目前國內、外只有針對原始異常邊界識別方法的研究，并未開展張量數據的研究，以及張量數據加權方式對邊界識別結果精度影響的討論，而利用重磁張量邊界識別技術可準確地獲得地質水平位置信息；利用延拓或者相關方法來進行張量數據的邊界識別，獲得地質體的空間位置信息。通過研究重磁張量數據的處理方法，進一步開發適用于全張量重磁數據的處理手段，并解決現有邊界識別方法所獲得的地質體邊界范圍較大，以及不能獲得深度信息的問題，提出新思路。在這一研究背景下，筆者分析了各種邊界識別方法，比較優缺點，并加以改進，整理出優化方案，最后得到新方法的成果圖像。

1 重磁張量數據的特征

張量表示為G的二階導數，可采用矩陣的形式表示為：

(1)

根據拉普拉斯方程可知gxy=gyx、gxz=gzx、gyz=gzy，張量測量僅有五個是獨立的。即：

(2)

根據棱柱體重力張量異常(圖1)分析張量屬性：gxx能有效地增強地質體的南北向邊界；gyy能增強異常的東西向邊界；gxy能描述地質體的角點信息；gzx的最大值能識別出地質體的南北邊界；gzy能識別地質體的東西邊界，不同的張量異常能反映出地質體的不同邊界信息。

2 邊界識別方法

常用重磁數據邊界識別方法如下：

1)Cordell[10]提出總水平導數(THD)。

圖1 棱柱體重力張量異常Fig.1 Gravity tensor abnormity of prismoid

(3)

2)Miller and Singh[7]傾斜角法。

(4)

3)Verduzco等[8]提出利用tilt angle的總水平導數(THDR)進行邊界的識別。

(5)

4)Wijns等[9]利用總水平導數與解析信號的比值(theta map)來進行邊界識別。

(6)

5)Cooper等[11]提出正則傾斜角法來進行異常體邊界的識別。

(7)

采用這五種邊界識別方法分別進行一正一負密度模型正演后再進行邊界識別處理(圖2)，通過比較可以看出這五種方法各有利弊。

圖2 不同邊界識別方法比較Fig.2 Compare with different boundary identification methods(a)原始異常數據;(b)總水平導數;(c)傾斜角法;(d)THDR;(e)theta map;(f)TDX

由圖2可以發現：圖2(a)原始異常數據對邊界識別效果不好，邊界非常模糊；圖2(b)總水平導數法較原始數據有所加強，但是邊界識別效果仍然有些模糊；圖2(c)傾斜角法、圖2(d)THDR 、圖2(e)theta map、圖2(f)正則化傾斜角法這四種方法會產生多余邊界。

3 改進的邊界識別方法

由上述結論發現：Cooper[11]提出正則傾斜角法來進行異常體邊界的識別(TDX)對于異常體的邊界呈現效果更好，但是卻在一正一負密度體模型成像時出現虛假異常邊界。所以，可以針對TXD方法進行改進，消除虛假異常邊界。

目前，對邊界識別方法的改進方面，馬國慶等[3]利用垂直導數的水平導數構成均衡邊界識別濾波器進行地質體的邊界識別；王萬銀等[2]提出歸一化總水平導數垂向導數進行地質體的邊界識別，并分析了位場總水平導數極值和垂直導數零值的空間位置變化，總結得到總水平導數極大值隨著地質體埋深的加大與地質體邊界偏差越大，利用垂直導數的總水平導數進行地質體的邊界識別可減小偏差。在這些成果的啟發下，需將TDX的公式進行適當的調整 (由于TDX只在一正一負密度模型中的邊界識別結果有虛假異常，而在同為正或負的密度模型中沒有，因此只討論前一種情況。)。

首先，分析TXD方法產生多余邊界的原因：

(7)

其中f為重磁異常值。如式(7)所示，該方法采用極大值法來進行地質體邊界的識別，即當垂直導數為零時取得極大值對應地質體邊界。而在一正一負密度模型中的垂直導數Vzz為“0”的點有5處(圖3)，其中四處在地質體邊界處，另外一處因為正密度地質體Vzz的負值和負密度地質體的Vzz正值相疊加產生的“0”，即為多余邊界產生的原因(其中Vzz為?f/?z，即垂直導數)。

圖3 重力垂直導數剖面圖Fig.3 Vertical derivative of gravity data profile

對TDX進行改進，提出一種基于TDX改進的新方法CTDX(Change TDX)：

(8)

由式(8)可以看出：CTDX是在TDX的基礎上在其分母加一個Δ值。使原有方法產生多余邊界處的總水平導數與式(8)分母的比值趨近于零值；而在真實的地質體邊界處總水平導數與式(8)分母的比值不為零值。這樣，可以使TDX中產生的多余邊界和地質體邊界分離開。

對TDX與CTDX方法進行比較(圖4)。

圖4 TDX和CTDX比較Fig.4 Comparison between TDX and CTDX(a)TDX ；(b)CTDX

從圖4可以看出，CTDX沒有TDX里的虛假異常邊界，并且其邊界識別能力也沒有得到影響。說明CTDX的方法可行(其中小數Δ取值0.01)，但是對于小數Δ的取值還有待商榷。

對于小變量或定量Δ，不同的取值會對邊界的識別產生不同的效果。因此對相同模型不同Δ取值進行了實驗：

Δ分別取0、0.000 1、0.001、0.01、0.1、1，以及總水平導數(圖5)。

從圖5可以得出結論,Δ太小依然會產生多余邊界，Δ增大會使結果逐漸趨近與總水平導數，而總水平導數對較深地質體效果不好。Δ的取值將受不同地質條件產生變化，所以Δ值的確定是目前研究的難點和重點。

以上Δ取值為具體數值，具有一定的經驗性，為了將Δ的取值問題應用到實際的數據處理過程中，進行了不同Δ計算方式的模型試算：

圖5 不同Δ取值結果比較圖Fig.5 Comparison with different Δ(a)Δ=0；(b)Δ=0.0001；(c)Δ=0.001；(d)Δ=0.01；(e)Δ=0.1；(f)Δ=1；(g)總水平導數

圖6 不同Δ計算方式的模型試算圖Fig.6 Comparison with different models of compute Δ(a)采用TDX對模型進行計算;(b)～(h)使用取垂直導數Vzz的最大值和最小值之差的1/2、1/5、1/10、1/20、1/30、1/40、1/50;(i)模型的垂直導數

將圖6(b)、圖6(c)、圖6(d)、圖6(e)和圖6(i)進行比較發現，若Δ取值太大，對于深處的地質體的邊界刻畫能力較弱，但是比較圖6(f)、圖6(g)、圖6(h)可以發現，Δ取值太小，對多余邊界的削弱程度不夠，仍然會有多余邊界出現。綜合上述對試算模型的結果分析，可以將Δ計算方式定義為式(9)。

Δ=|(Vzzmax-Vzzmin)|/30

(9)

此時的CTDX既可以對較深地質體邊界進行較好地識別，也沒有產生多余邊界，而對兩個正密度體模型也有著較好的效果(圖7)。

圖7 CTDX對不同模型的效果Fig.7 CTDX for different models(a)一正一負密度體模型； (b)兩個正密度體模型

4 噪聲壓制后的效果

通過正演理論模型試驗本文方法的抗噪性，在圖5正演模型中加入了隨機噪聲，得到了兩個正密度體加噪聲模型的處理結果(圖8)，一正一負密度體加噪聲模型的處理結果(圖9)。

圖8 兩個正密度體加噪聲Fig.8 Two positive density body models with noise(a)TDX ；(b)CTDX

圖9 一正一負密度體加噪聲Fig.9 One positive density body models and one negative density body models with noise(a)TDX ；(b)CTDX

通過對比處理結果可以得出結論：新方法與傳統方法相比，能較好地壓制噪聲和提高信噪比的效果。在隨機噪聲的干擾下，仍能較好地反映出異常體邊界。

5 實際數據的處理

根據研究區四川地區的重磁異常值(圖10、圖11)，可以發現，雖然該方法可以識別大體上的一些構造單元，但是對于有一些構造單元會產生多余的處理，加大實際解釋時的難度。將TDX方法與改進后的CTDX方法(圖12)行對比不難發現，TDX方法在處理磁異常數據時效果不是很好，對于異常的識別顯得非常雜亂，如果想要詳細劃分解釋必然要經過更多地處理。

分別經過TDX、CTDX邊界識別方法處理后，可以得到如圖12～圖15所示的結果。

圖10 研究區原始重力數據圖Fig.10 The original gravity data of the area

圖11 研究區磁數據圖Fig.11 The original magnetic data of the area

圖12 TDX方法處理研究區重力數據Fig.12 TDX of the gravity data in research area

圖13 TDX方法處理研究區重力數據Fig.13 TDX of the gravity data in research area

圖14 TDX方法處理研究區磁數據Fig.14 TDX of the magnetic data in research area

圖15 CTDX方法處理研究區磁數據Fig.15 CTDX of the magnetic data in research area

圖13相較圖12少了一些多余的邊界(紅色框中的陰影)，但對于已經可以識別的構造單元能更好地識別其邊界，并且有效地增強了深度場源體的邊界，使構造劃分工作變得簡單。

比較圖15與圖14，兩者的效果非常明顯，可以看出CTDX相較于TDX的處理結果可以更清楚地劃分出構造單元。從圖15與圖13結果結合來看，研究區內應存在一定數量的北東走向地質構造體，依據研究區的地質資料，推測是一系列的斷裂。且CTDX方法對于重力數據的改進效果有限，對磁數據的改進效果較為可觀。

以上從重、磁兩各方面印證了改進的正則化傾斜角法，對于邊界識別問題的效果明顯，說明改進的正則化傾斜角法可以在一定程度上優化邊界識別問題，去除多余邊界。

6 結論

總結了前人對異常體邊界識別的方法，分析了各個方法的優缺點，并在正則傾斜角法基礎上進行改進，從虛假異常邊界產生的原因入手，引入Δ=|(Vzzmax-Vzzmin)|/30，經過數據處理后，能夠很好地去除掉虛假異常邊界，并針對噪聲干擾進行模型實驗，發現該方法(CTDX)對噪聲的壓制效果也很可觀。通過對實際數據進行處理，可以發現改進后的方法可以去除多余邊界，增強深度場源體的邊界效果。

[1] 馬國慶, 黃大年, 于平, 等.改進的均衡濾波器在位場數據邊界識別中的應用[J].地球物理學報, 2012,55(12)：4288-4295.

MA G Q, HUANG D N, YU P, et al. Application of improved balancing filter to edge identification of potential field data[J].Chinese J Geophys, 2012,55 (12): 4288-4295.(In Chinese)

[2] 王萬銀，潘玉，邱之云.位場數據歸一化總水平導數垂向導數邊緣識別方法(英文)[J];Applied Geophysics,2009(03)：226 - 233.

WANG WY, PAN Y , QIU ZY . A new edge recognition technology based on the normalized vertical derivative of the total horizontal derivative for potential field data[J]. Applied Geophysics, 2009(3):226 - 233. (In Chinese)

[3] 馬國慶, 杜曉娟, 李麗麗. 利用水平與垂直導數的相關系數進行位場數據的邊界識別[J]. 吉林大學學報(地球科學版), 2011, 41(S1) 345-348.

MA G Q，DU X J，LI L L. Edge detection of potential filed data using correlation coefficients of horizontal and vertical derivatives[J].Journal of Jilin University( Earth Science Edition), 2011, 41(S1): 345-348.(In Chinese)

[4] 劉艷,李桐林,黃旭釗,等.航磁異常軸向特征提取算法研究及其應用[J].地球物理學進展,2002(1)： 141-144.

LIU Y， LI TL， HUANG X Z, et al. An automatic method of detection of aeromagnetic axes and its application[J].Progress in Geophysics, 2002 (01) :141-144.(In Chinese)

[5] 柏冠軍,吳漢寧,趙希剛,等.重力資料識別鄂爾多斯盆地線性構造方法研究[J].地球物理學進展,2007(05):1386-1392.

BAI G J,WU H N,ZHAO X G,et al.Research on recognition of linear structures using gravity data in ordos basin[J].Progress in Geophysics,2007(05): 1386-1392. (In Chinese)

[6] 郭華,吳燕岡,高鐵.重力斜導數方法在時間域中的理論模型與研究[J].吉林大學學報(地球科學版),2006(S2): 9-14.

GUO H, WU Y G, GAO T.The research of theories model in time area with the method of title derivative in gravity[J].Journal of Jilin University( Earth Science Edition),2006 (S2): 9-14.(In Chinese)

[7] MILLER H G, SINGH V. Potential field tilt—a new concept for location of potential field sources[J].Journal of Applied Geophysics, 1994, 32: 213-217.

[8] VERDUZCO B., FAIRHEAD, J. D. ,Green, C. M. New insights into magnetic derivatives for structural mapping[J].The Leading Edge, 2004, 23(2):116-119.

[9] EVJEN H M. The place of the vertical gradient in gravitational interpretations[J].Geophysics,1936(1):127-136.

[10] CORDELL L. Gravimetric expression of graben faulting in Santa Fe Country and the Espanola Basin, New Mexico[J]. New Mexico Geol. Soc. Guidebook, 30th Field Conf,1979:59 - 64.

[11] COOPER G R J, Cowan D R. Enhancing potential field data using filters based on the local phase[J]. Computers &Geosciences, 2006, 32: 1585-1591.

Investigationaboutedgerecognitionusing(gradient)tensordataofmagneticandgravity

ZHAO Jianyu , MING Yanbo , SUN Chengcheng, ZHANG Zhidong , LIU Xiaojia

(College of Geoexploration Science And Technology Jilin University , Changchun 130021, China)

The boundary identification is one of the regular tasks of gravity and magnetic data interpretation, but gravity and magnetic raw anomaly data does not well correspond to the boundary of geological body. The maximum value of the horizontal derivative of gravity and magnetic anomalies are corresponding to the zero value of vertical derivative, which is used mostly to complete the assignment of gravity and magnetic boundary identification. We find that false boundary is caused by using current boundary identification method. This article summarizes the current boundary identification method and corrects the problem of excessive boundaries to present new boundary identification method. This new method can effectively remove excessive boundaries and reduces the disturbance of noise, which provides a new thinking way for gravity and magnetic data boundary identification method.

boundary identification; (gradient) tensor data of magnetic and gravity; gravity and magnetic exploration

2016-10-21 改回日期： 2017-03-03

國家重點研發項目(2016YFC0600505)

趙建宇(1993-)，男，碩士，主要研究方向為固體地球物理學，E-mai：lzjy849891627@163.com。

明彥伯(1994-)，男，碩士，主要研究方向為固體地球物理學，E-mai：ybming1994@qq.com。

1001-1749(2017)06-0748-07

P 631.2

A

10.3969/j.issn.1001-1749.2017.06.06