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

(吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130021)

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

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

(吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130021)

邊界識(shí)別是重力數(shù)據(jù)解釋的常規(guī)任務(wù)之一，但地質(zhì)體的邊界不能很好地與重力原始異常數(shù)據(jù)對(duì)應(yīng)，而與重磁異常水平導(dǎo)數(shù)極大值、垂直導(dǎo)數(shù)零值相對(duì)應(yīng)，因此大多利用該性質(zhì)完成重磁邊界識(shí)別任務(wù)。研究發(fā)現(xiàn)，現(xiàn)有邊界識(shí)別方法會(huì)出現(xiàn)虛假邊界，針對(duì)現(xiàn)有邊界識(shí)別方法進(jìn)行總結(jié)，并針對(duì)多余邊界問(wèn)題進(jìn)行改進(jìn)，提出新的邊界識(shí)別，通過(guò)模型試驗(yàn)和實(shí)際數(shù)據(jù)證明，該方法可有效地去除多余邊界，且降低了噪聲的干擾，為重磁數(shù)據(jù)邊界識(shí)別方法提供新的思路。

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

0 引言

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

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

1 重磁張量數(shù)據(jù)的特征

張量表示為G的二階導(dǎo)數(shù)，可采用矩陣的形式表示為：

(1)

根據(jù)拉普拉斯方程可知gxy=gyx、gxz=gzx、gyz=gzy，張量測(cè)量?jī)H有五個(gè)是獨(dú)立的。即：

(2)

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

2 邊界識(shí)別方法

常用重磁數(shù)據(jù)邊界識(shí)別方法如下：

1)Cordell[10]提出總水平導(dǎo)數(shù)(THD)。

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

(3)

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

(4)

3)Verduzco等[8]提出利用tilt angle的總水平導(dǎo)數(shù)(THDR)進(jìn)行邊界的識(shí)別。

(5)

4)Wijns等[9]利用總水平導(dǎo)數(shù)與解析信號(hào)的比值(theta map)來(lái)進(jìn)行邊界識(shí)別。

(6)

5)Cooper等[11]提出正則傾斜角法來(lái)進(jìn)行異常體邊界的識(shí)別。

(7)

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

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

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

3 改進(jìn)的邊界識(shí)別方法

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

目前，對(duì)邊界識(shí)別方法的改進(jìn)方面，馬國(guó)慶等[3]利用垂直導(dǎo)數(shù)的水平導(dǎo)數(shù)構(gòu)成均衡邊界識(shí)別濾波器進(jìn)行地質(zhì)體的邊界識(shí)別；王萬(wàn)銀等[2]提出歸一化總水平導(dǎo)數(shù)垂向?qū)?shù)進(jìn)行地質(zhì)體的邊界識(shí)別，并分析了位場(chǎng)總水平導(dǎo)數(shù)極值和垂直導(dǎo)數(shù)零值的空間位置變化，總結(jié)得到總水平導(dǎo)數(shù)極大值隨著地質(zhì)體埋深的加大與地質(zhì)體邊界偏差越大，利用垂直導(dǎo)數(shù)的總水平導(dǎo)數(shù)進(jìn)行地質(zhì)體的邊界識(shí)別可減小偏差。在這些成果的啟發(fā)下，需將TDX的公式進(jìn)行適當(dāng)?shù)恼{(diào)整 (由于TDX只在一正一負(fù)密度模型中的邊界識(shí)別結(jié)果有虛假異常，而在同為正或負(fù)的密度模型中沒(méi)有，因此只討論前一種情況。)。

首先，分析TXD方法產(chǎn)生多余邊界的原因：

(7)

其中f為重磁異常值。如式(7)所示，該方法采用極大值法來(lái)進(jìn)行地質(zhì)體邊界的識(shí)別，即當(dāng)垂直導(dǎo)數(shù)為零時(shí)取得極大值對(duì)應(yīng)地質(zhì)體邊界。而在一正一負(fù)密度模型中的垂直導(dǎo)數(shù)Vzz為“0”的點(diǎn)有5處(圖3)，其中四處在地質(zhì)體邊界處，另外一處因?yàn)檎芏鹊刭|(zhì)體Vzz的負(fù)值和負(fù)密度地質(zhì)體的Vzz正值相疊加產(chǎn)生的“0”，即為多余邊界產(chǎn)生的原因(其中Vzz為?f/?z，即垂直導(dǎo)數(shù))。

圖3 重力垂直導(dǎo)數(shù)剖面圖Fig.3 Vertical derivative of gravity data profile

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

(8)

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

對(duì)TDX與CTDX方法進(jìn)行比較(圖4)。

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

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

對(duì)于小變量或定量Δ，不同的取值會(huì)對(duì)邊界的識(shí)別產(chǎn)生不同的效果。因此對(duì)相同模型不同Δ取值進(jìn)行了實(shí)驗(yàn)：

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

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

以上Δ取值為具體數(shù)值，具有一定的經(jīng)驗(yàn)性，為了將Δ的取值問(wèn)題應(yīng)用到實(shí)際的數(shù)據(jù)處理過(guò)程中，進(jìn)行了不同Δ計(jì)算方式的模型試算：

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

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

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

Δ=|(Vzzmax-Vzzmin)|/30

(9)

此時(shí)的CTDX既可以對(duì)較深地質(zhì)體邊界進(jìn)行較好地識(shí)別，也沒(méi)有產(chǎn)生多余邊界，而對(duì)兩個(gè)正密度體模型也有著較好的效果(圖7)。

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

4 噪聲壓制后的效果

通過(guò)正演理論模型試驗(yàn)本文方法的抗噪性，在圖5正演模型中加入了隨機(jī)噪聲，得到了兩個(gè)正密度體加噪聲模型的處理結(jié)果(圖8)，一正一負(fù)密度體加噪聲模型的處理結(jié)果(圖9)。

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

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

通過(guò)對(duì)比處理結(jié)果可以得出結(jié)論：新方法與傳統(tǒng)方法相比，能較好地壓制噪聲和提高信噪比的效果。在隨機(jī)噪聲的干擾下，仍能較好地反映出異常體邊界。

5 實(shí)際數(shù)據(jù)的處理

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

分別經(jīng)過(guò)TDX、CTDX邊界識(shí)別方法處理后，可以得到如圖12～圖15所示的結(jié)果。

圖10 研究區(qū)原始重力數(shù)據(jù)圖Fig.10 The original gravity data of the area

圖11 研究區(qū)磁數(shù)據(jù)圖Fig.11 The original magnetic data of the area

圖12 TDX方法處理研究區(qū)重力數(shù)據(jù)Fig.12 TDX of the gravity data in research area

圖13 TDX方法處理研究區(qū)重力數(shù)據(jù)Fig.13 TDX of the gravity data in research area

圖14 TDX方法處理研究區(qū)磁數(shù)據(jù)Fig.14 TDX of the magnetic data in research area

圖15 CTDX方法處理研究區(qū)磁數(shù)據(jù)Fig.15 CTDX of the magnetic data in research area

圖13相較圖12少了一些多余的邊界(紅色框中的陰影)，但對(duì)于已經(jīng)可以識(shí)別的構(gòu)造單元能更好地識(shí)別其邊界，并且有效地增強(qiáng)了深度場(chǎng)源體的邊界，使構(gòu)造劃分工作變得簡(jiǎn)單。

比較圖15與圖14，兩者的效果非常明顯，可以看出CTDX相較于TDX的處理結(jié)果可以更清楚地劃分出構(gòu)造單元。從圖15與圖13結(jié)果結(jié)合來(lái)看，研究區(qū)內(nèi)應(yīng)存在一定數(shù)量的北東走向地質(zhì)構(gòu)造體，依據(jù)研究區(qū)的地質(zhì)資料，推測(cè)是一系列的斷裂。且CTDX方法對(duì)于重力數(shù)據(jù)的改進(jìn)效果有限，對(duì)磁數(shù)據(jù)的改進(jìn)效果較為可觀。

以上從重、磁兩各方面印證了改進(jìn)的正則化傾斜角法，對(duì)于邊界識(shí)別問(wèn)題的效果明顯，說(shuō)明改進(jìn)的正則化傾斜角法可以在一定程度上優(yōu)化邊界識(shí)別問(wèn)題，去除多余邊界。

6 結(jié)論

總結(jié)了前人對(duì)異常體邊界識(shí)別的方法，分析了各個(gè)方法的優(yōu)缺點(diǎn)，并在正則傾斜角法基礎(chǔ)上進(jìn)行改進(jìn)，從虛假異常邊界產(chǎn)生的原因入手，引入Δ=|(Vzzmax-Vzzmin)|/30，經(jīng)過(guò)數(shù)據(jù)處理后，能夠很好地去除掉虛假異常邊界，并針對(duì)噪聲干擾進(jìn)行模型實(shí)驗(yàn)，發(fā)現(xiàn)該方法(CTDX)對(duì)噪聲的壓制效果也很可觀。通過(guò)對(duì)實(shí)際數(shù)據(jù)進(jìn)行處理，可以發(fā)現(xiàn)改進(jìn)后的方法可以去除多余邊界，增強(qiáng)深度場(chǎng)源體的邊界效果。

[1] 馬國(guó)慶, 黃大年, 于平, 等.改進(jìn)的均衡濾波器在位場(chǎng)數(shù)據(jù)邊界識(shí)別中的應(yīng)用[J].地球物理學(xué)報(bào), 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] 王萬(wàn)銀，潘玉，邱之云.位場(chǎng)數(shù)據(jù)歸一化總水平導(dǎo)數(shù)垂向?qū)?shù)邊緣識(shí)別方法(英文)[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] 馬國(guó)慶, 杜曉娟, 李麗麗. 利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行位場(chǎng)數(shù)據(jù)的邊界識(shí)別[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 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] 劉艷,李桐林,黃旭釗,等.航磁異常軸向特征提取算法研究及其應(yīng)用[J].地球物理學(xué)進(jìn)展,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] 柏冠軍,吳漢寧,趙希剛,等.重力資料識(shí)別鄂爾多斯盆地線(xiàn)性構(gòu)造方法研究[J].地球物理學(xué)進(jìn)展,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] 郭華,吳燕岡,高鐵.重力斜導(dǎo)數(shù)方法在時(shí)間域中的理論模型與研究[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),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

國(guó)家重點(diǎn)研發(fā)項(xiàng)目(2016YFC0600505)

趙建宇(1993-)，男，碩士，主要研究方向?yàn)楣腆w地球物理學(xué)，E-mai：lzjy849891627@163.com。

明彥伯(1994-)，男，碩士，主要研究方向?yàn)楣腆w地球物理學(xué)，E-mai：ybming1994@qq.com。

1001-1749(2017)06-0748-07

P 631.2

A

10.3969/j.issn.1001-1749.2017.06.06