譚曉迪，李麗麗，馬國慶，張代磊

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

0 前言

在位場(chǎng)數(shù)據(jù)處理過程中，對(duì)場(chǎng)源進(jìn)行邊界識(shí)別是不可或缺的，其實(shí)質(zhì)是根據(jù)目標(biāo)體邊界兩側(cè)存在明顯的密度及磁性差異來劃分邊界，從而確定目標(biāo)體的水平位置。場(chǎng)源邊界識(shí)別方法分為數(shù)理統(tǒng)計(jì)方法和數(shù)值計(jì)算方法。

數(shù)理統(tǒng)計(jì)方法主要有小子域?yàn)V波[1]和歸一化標(biāo)準(zhǔn)差法[2]，這類方法都是采用滑動(dòng)窗口的方式對(duì)異常數(shù)據(jù)進(jìn)行統(tǒng)計(jì)運(yùn)算，優(yōu)點(diǎn)在于能夠利用窗口來控制噪聲，但窗口的尺寸通常憑經(jīng)驗(yàn)給定，因此存在主觀因素，會(huì)影響計(jì)算結(jié)果。

數(shù)值計(jì)算方法主要包括垂向?qū)?shù)、總水平導(dǎo)數(shù)(THDR)、解析信號(hào)振幅(ASM)、傾斜角(TILT)和Theta圖法。這5種方法是最基本的數(shù)值計(jì)算方法，在此之后所提出的數(shù)值計(jì)算方法都是以這5種方法進(jìn)行改進(jìn)得到。Hood等[3]提出利用垂向一階導(dǎo)數(shù)確定鉛錘臺(tái)階邊緣位置，垂向?qū)?shù)是通過零值點(diǎn)的位置來識(shí)別邊界，該方法分辨能力強(qiáng)，但易受高頻干擾。Cordell[4]提出總水平導(dǎo)數(shù)法，通過對(duì)重力異常求總水平導(dǎo)數(shù)來劃分目標(biāo)體邊界，但該方法極易受場(chǎng)源埋深影響，不能有效識(shí)別出深部的、小的地質(zhì)體的邊緣位置。Nabighian[5]提出解析信號(hào)振幅法，該方法受磁異常分量和磁化方向影響較小，不足之處是分辨率低。上述3種方法均不能很好地均衡不同振幅異常，尤其是振幅值存在較大差異時(shí)，邊界識(shí)別效果不理想。Miller等[6]提出第一個(gè)邊界識(shí)別均衡濾波器——傾斜角法，但用該方法進(jìn)行邊界識(shí)別不能突出邊界的異常值，劃分邊界不夠準(zhǔn)確。Theta圖法由Wijins等[7]提出，該方法能夠均衡高低幅值的重磁異常數(shù)據(jù)，但存在解析奇點(diǎn)。馬國慶等[8]提出利用總水平導(dǎo)數(shù)與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行地質(zhì)體邊界的識(shí)別，但該方法對(duì)小范圍異常體的識(shí)別效果不是很好。馬國慶等[9]進(jìn)而提出增強(qiáng)型濾波器圈定邊界，得到了較好的識(shí)別效果。綜上，數(shù)值計(jì)算方法便于實(shí)現(xiàn)且計(jì)算速度快，但在處理實(shí)際數(shù)據(jù)時(shí)，求導(dǎo)會(huì)放大噪聲，而且識(shí)別出的地質(zhì)體邊界過于寬泛，不能對(duì)邊界進(jìn)行精確的劃分；因此，如何解決噪聲干擾和邊界收斂問題需要進(jìn)一步研究。

小波變換由于具有多尺度分析的特點(diǎn)，能夠有效地提取目標(biāo)體特征，同時(shí)具有較強(qiáng)的抗噪聲能力，在數(shù)據(jù)處理中得到了廣泛的應(yīng)用。徐亞等[10]總結(jié)了小波變換在位場(chǎng)數(shù)據(jù)應(yīng)用中的理論方法及發(fā)展現(xiàn)狀，認(rèn)為小波變換可應(yīng)用于位場(chǎng)的分離、去噪、地質(zhì)體邊緣提取等。

本文采用小波結(jié)合傳統(tǒng)邊界識(shí)別方法對(duì)數(shù)據(jù)進(jìn)行處理，驗(yàn)證小波對(duì)噪聲的壓制能力。針對(duì)傳統(tǒng)方法在進(jìn)行邊界識(shí)別時(shí)存在的邊界收斂問題，本文提出冪次變換的方法對(duì)邊界進(jìn)行收斂，并將這兩種方法相結(jié)合應(yīng)用于模型和實(shí)際位場(chǎng)數(shù)據(jù)的邊界識(shí)別，同時(shí)與傳統(tǒng)方法識(shí)別效果進(jìn)行對(duì)比分析，證實(shí)改進(jìn)的方法獲得了更好的邊界識(shí)別效果。

1 小波變換

1.1 離散小波變換

設(shè)g(t)為母小波函數(shù)，且g(t)∈L2(R)，L(R)為平方可積實(shí)數(shù)空間，對(duì)g(t)做伸縮和平移后得到[11]：

(1)

式中：t∈R；ga,b(t)為小波函數(shù)，簡稱小波；a是尺度因子,表征函數(shù)的伸縮程度；b是位置因子，反映函數(shù)的平移位置。

對(duì)a,b進(jìn)行離散化，即：

b=kb0,b0∈Z,k∈Z。

則離散小波變?yōu)?/p>

(2)

信號(hào)s(t)∈L2(R)的離散小波變換定義為

式中：g*(t)表示g(t)的復(fù)共軛。

其反變換[11]定義為

(4)

1.2 多分辨率分析理論和Mallat算法

多分辨率分析是將L2(R)空間分解為不同分辨率的子空間序列，其子空間序列的極限就是L2(R)。其中任意函數(shù)s(t)∈L2(R)，s(t)投影到每一個(gè)子空間得到一系列近似函數(shù)，這一系列近似函數(shù)逼近極限就是s(t)。通過研究各個(gè)子空間的投影，從而可分析s(t)在不同分辨率子空間上的性質(zhì)和特征[12]。

在多分辨率分析中，Vj稱為尺度空間。通常，不同的尺度空間對(duì)應(yīng)尺度不同函數(shù)，從而可用來研究總空間L2(R)的不同分辨率特征。設(shè)φ是尺度函數(shù)，則存在系數(shù)序列h(k),k∈Z[12]，得到

(5)

對(duì)任意的s(t)∈Vj，首先將其分解為細(xì)節(jié)部分和近似部分，然后再將近似部分進(jìn)一步分解；再重復(fù)以上兩步，直到得到任意尺度上的細(xì)節(jié)部分和近似部分，這就是Mallat算法的原理。s(t)在Vj空間的展開式是

(6)

式中：k∈Z;j∈Z；cj,k為j尺度上的近似系數(shù)。

將s(t)分解一次得到

(7)

式中：dj-1,k是j-1尺度上的細(xì)節(jié)系數(shù)。將cj-1,k進(jìn)一步分解下去，可以分別得到系數(shù)cj-2,k和dj-2,k；再繼續(xù)分解，從而達(dá)到任意空間尺度的分解系數(shù)。其重構(gòu)算法表達(dá)式如下：

(8)

其中：n∈Z；h(k)(k∈Z)相當(dāng)于低通濾波器；g(k)(k∈Z)相當(dāng)于高通濾波器。由式(8)可以看出，對(duì)信號(hào)進(jìn)行分解得到的近似系數(shù)，是信號(hào)與低通濾波器作用的結(jié)果，而細(xì)節(jié)系數(shù)是信號(hào)和高通濾波器作用得到的；因此近似系數(shù)相對(duì)細(xì)節(jié)系數(shù)包含更多的低頻信息和更少的高頻信息，從而達(dá)到壓制噪聲的目的。

圖1給出了第一層的小波分解和重構(gòu)計(jì)算過程，首先將信號(hào)分解為低頻信息L與高頻信息H；然后再將低頻信息L分解為低頻信息LL1與其高頻信息LH1，同時(shí)將高頻信息H分解為HL1和HH1[13]。其中LL1為第一層近似系數(shù)，其余成分為細(xì)節(jié)系數(shù)，這就是第一層的分解過程；第二層的分解是以LL1為原始信號(hào)重復(fù)相同的操作[14]，以此類推，如圖2所示。

圖1 Mallat算法第一層分解重構(gòu)圖示Fig.1 First layer of Mallat algorithm decomposition and reconstruction

圖2 Mallat算法塔式分解示意圖Fig.2 Sketch map of Mallat algorithm pyramid decomposition

當(dāng)對(duì)所分解的小尺度上的系數(shù)進(jìn)行重構(gòu)時(shí)，我們只選取近似系數(shù)，這樣可以有效壓制噪聲。層數(shù)越高，對(duì)噪聲的壓制作用越明顯，但會(huì)丟失細(xì)節(jié)信息，導(dǎo)致邊界識(shí)別結(jié)果發(fā)散；因此，本文選取小波變換的第一層和第二層近似系數(shù)進(jìn)行重構(gòu)來驗(yàn)證試驗(yàn)效果。

2 冪次變換原理

針對(duì)傳統(tǒng)邊界識(shí)別方法在處理數(shù)據(jù)時(shí)邊界不夠收斂、識(shí)別效果不理想的問題，本文提出冪次變換的方法來提高收斂。該方法實(shí)質(zhì)上是在傳統(tǒng)數(shù)值計(jì)算方法上加以改進(jìn)，提高整體計(jì)算的冪級(jí)數(shù)。本文選取了總水平導(dǎo)數(shù)、傾斜角和Theta圖這3種邊界識(shí)別中應(yīng)用廣泛且效果較好的數(shù)值計(jì)算方法進(jìn)行冪次變換，并與原方法的識(shí)別效果進(jìn)行對(duì)比，驗(yàn)證收斂效果。

2.1 數(shù)值計(jì)算方法

總水平導(dǎo)數(shù)法是通過極大值來判定地質(zhì)構(gòu)造的邊界。該方法對(duì)噪聲不敏感，但當(dāng)場(chǎng)源埋深增大時(shí)，識(shí)別出的邊界會(huì)與實(shí)際邊界位置存在較大偏差[15]。假設(shè)V表示位場(chǎng)異常數(shù)據(jù)，總水平導(dǎo)數(shù)計(jì)算公式為

(9)

傾斜角法將垂向?qū)?shù)和水平導(dǎo)數(shù)做比值計(jì)算，可有效克服導(dǎo)數(shù)計(jì)算隨深度增大迅速衰減的問題。傾斜角計(jì)算公式為

(10)

該方法受場(chǎng)源埋深影響較小，可實(shí)現(xiàn)在場(chǎng)源內(nèi)部為正值，邊界附近為零值，場(chǎng)源外部為負(fù)值，從而識(shí)別深部和淺部場(chǎng)源體的邊界[16]；但這種識(shí)別方法效果不理想，并不能直觀地顯示出目標(biāo)體的邊緣位置。為了突出目標(biāo)體的邊界，本文對(duì)傾斜角法做進(jìn)一步的改進(jìn)，稱其為可視化傾斜角法，可視化傾斜角ETILT的公式為

ETILT=|TILT2-C|。

(11)

式中：C為常數(shù)，通常情況下C≥2。這種方法能夠提高邊界的異常幅值，降低周圍的異常值，能夠更好地突出邊界的位置信息。

Theta圖法是求取總水平導(dǎo)數(shù)比上總水平導(dǎo)數(shù)與垂直導(dǎo)數(shù)之和的夾角θ的余弦，公式為

(12)

Theta圖的數(shù)值在[0, 1]范圍內(nèi)，利用極大值進(jìn)行場(chǎng)源體的邊界識(shí)別，并可增強(qiáng)弱異常、降低場(chǎng)源埋深的影響。

這幾種常用的數(shù)值計(jì)算方法在邊界識(shí)別中都得到了較好的結(jié)果，同時(shí)也都存在明顯的缺點(diǎn)，因此本文提出冪次變換的方法加以改進(jìn)。

2.2 冪次變換

對(duì)以上幾種數(shù)值計(jì)算方法進(jìn)行冪次變換，給出表達(dá)式，總水平導(dǎo)數(shù)法為

(13)

可視化傾斜角法為

ETILTn=|TILT2-C|n；

(14)

Theta圖法為

(15)

式中：n表示級(jí)數(shù)(n≥1)。當(dāng)n=1時(shí)，與原數(shù)值計(jì)算方法相同。n的取值決定于原方法對(duì)地質(zhì)體邊界的識(shí)別能力，根據(jù)邊界識(shí)別效果擇最優(yōu)來選?。籲越大，對(duì)邊界的收斂程度越高，但級(jí)數(shù)過高會(huì)丟失細(xì)節(jié)信息，通常n≤8。

3 模型試驗(yàn)

重力異常模型包含3個(gè)剩余密度均勻的棱柱體，且在水平面上的投影位置沒有疊加，將棱柱體的起點(diǎn)(距離空間直角坐標(biāo)系0點(diǎn)最近)坐標(biāo)和尺寸(長寬高)分為2個(gè)數(shù)組，與剩余密度一起在表1中給出，它們的平面分布及空間分布分別如圖3、圖4所示。計(jì)算這個(gè)模型的重力異常響應(yīng)，并在正演數(shù)據(jù)中加入隨機(jī)噪聲，信噪比為35。正演數(shù)據(jù)和加噪聲數(shù)據(jù)如圖5所示。

表1 棱柱體位置、尺寸和剩余密度

圖3 重力異常模型棱柱體的平面分布Fig.3 Planar distribution of the prisms

圖4 重力異常模型棱柱體的空間分布Fig.4 Spatial distribution of the prisms

此模型正演數(shù)據(jù)均為正值，在棱柱體的邊界位置存在接觸區(qū)域，因此在進(jìn)行邊界識(shí)別時(shí)存在一定

難度。添加隨機(jī)噪聲可以有效檢驗(yàn)小波壓制噪聲的能力。將小波與冪次變換結(jié)合，與3種傳統(tǒng)的位場(chǎng)邊界識(shí)別數(shù)值計(jì)算方法總水平導(dǎo)數(shù)、傾斜角、Theta圖共同應(yīng)用于模型數(shù)據(jù)。模型數(shù)據(jù)及加噪聲數(shù)據(jù)小波處理后邊界識(shí)別結(jié)果分別如圖6、圖7所示，再進(jìn)行冪次變換，邊界識(shí)別效果如圖8、圖9所示，可以看出單獨(dú)小波處理以及與冪次變換結(jié)合2種方法在模型數(shù)據(jù)的應(yīng)用中都取得了較好的邊界識(shí)別效果。

圖6a、圖7a為可視化傾斜角法識(shí)別結(jié)果，圖6d、圖7d為總水平導(dǎo)數(shù)法識(shí)別結(jié)果，圖6g、圖7g為Theta圖法識(shí)別結(jié)果?？梢钥闯?，3種方法都基本給出了3個(gè)棱柱體真實(shí)的邊界信息，并且可視化傾斜角法的邊界識(shí)別效果與Theta圖法十分相似，均為高值表示目標(biāo)體邊界，能夠更好地突出邊緣的位置信息。圖6b和圖7b、圖6e和圖7e、圖6h和圖7h是分別在相應(yīng)的數(shù)值計(jì)算方法的基礎(chǔ)上進(jìn)行小波變換，選取第一層的近似系數(shù)進(jìn)行重構(gòu)得到，圖6c和圖7c、圖6f和圖7f、圖6i和圖7i是應(yīng)用小波變換第二層近似系數(shù)重構(gòu)得到，但因?yàn)閳D6是基于無噪聲的模型數(shù)據(jù)進(jìn)行處理，因此是否采用小波變換沒有明顯差別。由圖7可以看出，通過對(duì)分解得到的近似系數(shù)進(jìn)行重構(gòu)可以有效地起到壓制噪聲的作用，且層數(shù)越高，壓制作用越明顯，同時(shí)也會(huì)導(dǎo)致更多細(xì)節(jié)信息的損失。經(jīng)過試驗(yàn)，本文選擇兩層近似系數(shù)進(jìn)行重構(gòu)。圖8、圖9分別在圖6、圖7的基礎(chǔ)上，基于小波處理結(jié)果再進(jìn)行冪次變換得到，其中，可視化傾斜角和Theta圖進(jìn)行了8次冪變換，總水平導(dǎo)數(shù)進(jìn)行了4次冪變換。將圖6與圖8進(jìn)行對(duì)比，可以看出在無噪聲數(shù)據(jù)中，冪次變換能夠?qū)δ繕?biāo)體的邊界進(jìn)行有效的收斂，使得邊界劃分更加清晰；由圖7與圖9的對(duì)比可以看出，對(duì)于含噪聲的模型數(shù)據(jù)，小波變換能夠有效壓制噪聲，再結(jié)合冪次變換可以達(dá)到更好的邊界收斂效果，從而進(jìn)一步提高邊界識(shí)別精度。

圖5 模型重力異常(a)和含有噪聲的模型重力異常(b)Fig.5 Gravity anomaly of the model(a) and gravity anomaly contaminated with random noise(b)

a. ETILT邊界識(shí)別結(jié)果；b. ETILT第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；c. ETILT第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；d.THDR邊界識(shí)別結(jié)果；e. THDR第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；f. THDR第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；g. Theta圖邊界識(shí)別結(jié)果；h. Theta圖第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；i；Theta圖第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果。其中黑色邊框表示棱柱體在測(cè)量平面的水平位置投影。圖6 模型數(shù)據(jù)小波變換邊界識(shí)別結(jié)果Fig.6 Edge detection results of wavelet transfromed model data

a. ETILT邊界識(shí)別結(jié)果；b. ETILT第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；c. ETILT第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；d. THDR邊界識(shí)別結(jié)果；e. THDR第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；f. THDR第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；g. Theta圖邊界識(shí)別結(jié)果；h. Theta圖第一層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果；i. Theta圖第二層近似系數(shù)重構(gòu)邊界識(shí)別結(jié)果。其中黑色邊框表示棱柱體在測(cè)量平面的水平位置投影。圖7 含噪聲模型數(shù)據(jù)小波變換邊界識(shí)別結(jié)果Fig.7 Edge detection results of wavelet transfromed model data with noise

a. ETILT邊界識(shí)別結(jié)果；b. ETILT第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；c. ETILT第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；d. THDR邊界識(shí)別結(jié)果；e. THDR第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；f. THDR第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；g. Theta圖邊界識(shí)別結(jié)果；h. Theta圖第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；i. Theta圖第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果。其中黑色邊框表示棱柱體在測(cè)量平面的水平位置投影。圖8 模型數(shù)據(jù)小波結(jié)合冪次變換邊界識(shí)別結(jié)果Fig.8 Edge detection results of wavelet and power transfromed model data

a. ETILT邊界識(shí)別結(jié)果；b. ETILT第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；c. ETILT第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；d. THDR邊界識(shí)別結(jié)果；e. THDR第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；f. THDR第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；g. Theta圖邊界識(shí)別結(jié)果；h. Theta圖第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；i. Theta圖第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果。其中黑色邊框表示棱柱體在測(cè)量平面的水平位置投影。圖9 含有噪聲的模型數(shù)據(jù)小波結(jié)合冪次變換邊界識(shí)別結(jié)果Fig.9 Edge detection results of wavelet and power transfromed model data with noise

a. 四川盆地重力異常；b. ETILT邊界識(shí)別結(jié)果；c. ETILT第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；d. ETILT第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；e. THDR邊界識(shí)別結(jié)果；f. THDR第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；g. THDR第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；h. Theta圖邊界識(shí)別結(jié)果；i. Theta圖第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；j. Theta圖第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果。。圖10 四川盆地重力異常數(shù)據(jù)邊界識(shí)別結(jié)果Fig.10 Edge detection results of gravity anomaly of Sichuan basin

a. 經(jīng)過化磁極的磁異常數(shù)據(jù)；b. ETILT邊界識(shí)別結(jié)果；c. ETILT第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；d. ETILT第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；e. THDR邊界識(shí)別結(jié)果；f. THDR第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；g. THDR第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；h. Theta圖邊界識(shí)別結(jié)果；i. Theta圖第一層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果；j. Theta圖第二層近似系數(shù)重構(gòu)加冪次邊界識(shí)別結(jié)果。圖11 朱日和地區(qū)磁異常數(shù)據(jù)邊界識(shí)別結(jié)果Fig.11 Edge detection results of magnetic anomaly of Zhurihe region

4 實(shí)際數(shù)據(jù)應(yīng)用

為驗(yàn)證本文方法在實(shí)際位場(chǎng)數(shù)據(jù)中的應(yīng)用效果，我們分別處理了四川盆地重力異常數(shù)據(jù)和朱日和地區(qū)的磁測(cè)數(shù)據(jù)。

首先對(duì)四川盆地布格重力異常數(shù)據(jù)進(jìn)行邊界識(shí)別處理，結(jié)果見圖10，其中圖10a為布格重力異常圖。從圖10可以看出：總水平導(dǎo)數(shù)法對(duì)實(shí)測(cè)數(shù)據(jù)的識(shí)別效果不太理想，僅能識(shí)別出大的斷裂，對(duì)于一些細(xì)小斷裂不能清晰地顯示出來(圖10e)；可視化傾斜角和Theta圖的邊界識(shí)別結(jié)果有一定相似性，均可大體識(shí)別出研究區(qū)域內(nèi)存在的斷裂構(gòu)造，但構(gòu)造分布范圍較廣，未實(shí)現(xiàn)精確劃分(圖10b、h)；在總水平導(dǎo)數(shù)的基礎(chǔ)上對(duì)其進(jìn)行小波變換并分別選取第一層和第二層近似系數(shù)進(jìn)行重構(gòu)、再進(jìn)行二次冪變換后，噪聲均得到了壓制，同時(shí)識(shí)別出的斷裂位置更收斂(圖10f、g)；對(duì)可視化傾斜角和Theta圖做小波變換第一層近似系數(shù)重構(gòu)和冪次變換的處理，冪級(jí)數(shù)分別為3、4，則噪聲影響明顯降低，且斷裂水平位置更清晰，提高了識(shí)別精度(圖10c、i)；對(duì)可視化傾斜角和Theta圖做小波變換第二層系數(shù)重構(gòu)及冪次變換的處理，且冪級(jí)數(shù)不變，則噪聲被進(jìn)一步壓制，收斂程度提高，整體效果與第一層系數(shù)重構(gòu)處理結(jié)果大致相同(圖10d、j)。

將本文方法應(yīng)用于磁異常數(shù)據(jù)的邊界識(shí)別，對(duì)朱日和地區(qū)的磁數(shù)據(jù)進(jìn)行化極處理得到異常如圖11所示。從圖11a可明顯看出成礦帶的大致范圍；從總水平導(dǎo)數(shù)法的識(shí)別結(jié)果(圖11e)可以看出該方法能夠大致識(shí)別出斷裂的水平位置，但對(duì)斷裂的邊界識(shí)別效果不理想，而且對(duì)于細(xì)小的地質(zhì)構(gòu)造不能清晰劃分；圖11f是在總水平導(dǎo)數(shù)的基礎(chǔ)上做小波第一層近似系數(shù)重構(gòu)及冪次變換的結(jié)果，冪級(jí)數(shù)為2，該方法能夠清晰地識(shí)別出大斷裂的邊界位置，識(shí)別效果較為理想；圖11g是在總水平導(dǎo)數(shù)的基礎(chǔ)上做第二層系數(shù)重構(gòu)及冪次變換的結(jié)果，冪級(jí)數(shù)也為2，可以看出與第一層相比較，第二層系數(shù)重構(gòu)能夠進(jìn)一步壓制噪聲，但會(huì)損失細(xì)節(jié)信息，導(dǎo)致邊界發(fā)散；圖11b、h分別為可視化傾斜角法、Theta圖法的識(shí)別結(jié)果，2種方法對(duì)磁異常數(shù)據(jù)進(jìn)行處理的邊界識(shí)別效果相似，都能夠清晰識(shí)別出大斷裂和一些細(xì)小斷裂的邊界，但是邊界范圍不夠收斂；圖11c、i是分別對(duì)其做第一層近似系數(shù)重構(gòu)結(jié)合冪次變換的處理結(jié)果，冪級(jí)數(shù)分別為3和4，由圖看出該方法能夠有效對(duì)邊界進(jìn)行收斂，更加精確地劃分出斷裂邊界的位置；圖11d、j是采用第二層系數(shù)重構(gòu)結(jié)合冪次的結(jié)果，冪級(jí)數(shù)不變，可見，經(jīng)第二層重構(gòu)之后，會(huì)丟失細(xì)節(jié)信息，導(dǎo)致識(shí)別效果發(fā)散，因此在處理朱日和地區(qū)磁異常數(shù)據(jù)時(shí)，第一層系數(shù)重構(gòu)結(jié)合冪次變換能夠更好地劃分測(cè)區(qū)內(nèi)斷裂構(gòu)造的邊界。

綜上所述，在實(shí)際數(shù)據(jù)應(yīng)用中，可視化傾斜角法和Theta圖法相較于總水平導(dǎo)數(shù)法能夠更好地進(jìn)行地質(zhì)體的邊界識(shí)別，小波變換結(jié)合冪次變換能夠很好地起到壓制噪聲以及收斂邊界的作用。

5 結(jié)論

1)本文針對(duì)數(shù)值計(jì)算方法在位場(chǎng)數(shù)據(jù)邊界識(shí)別中受噪聲影響大、識(shí)別邊界不夠收斂的問題，提出應(yīng)用小波變換結(jié)合冪次變換的方法在原有方法上進(jìn)行改進(jìn)，并選取了3種傳統(tǒng)數(shù)值計(jì)算方法，將改進(jìn)方法應(yīng)用于模型數(shù)據(jù)和實(shí)際重磁數(shù)據(jù)，得到了較好的結(jié)果。

2)傳統(tǒng)傾斜角法不能直觀地顯示出目標(biāo)體的邊緣位置，為了突出目標(biāo)體的邊界，本文對(duì)傾斜角法做進(jìn)一步的改進(jìn)，即可視化傾斜角法，該方法能夠提高邊界異常的幅值，從而提高邊界識(shí)別的精度。

3)小波變換方法是將數(shù)據(jù)分解成近似系數(shù)和細(xì)節(jié)系數(shù)，再對(duì)近似系數(shù)進(jìn)行重構(gòu)，由于近似系數(shù)包含較少的噪聲信息，從而減少噪聲對(duì)邊界識(shí)別的影響。

4)本文提出冪次變換方法對(duì)識(shí)別出的邊界進(jìn)行有效收斂，通常冪級(jí)數(shù)越大，收斂效果越好，但由于冪級(jí)數(shù)過高會(huì)丟失細(xì)節(jié)信息，因此冪級(jí)數(shù)選取通常不大于8。

[1] 楊高印. 位場(chǎng)數(shù)據(jù)處理的一項(xiàng)新技術(shù)：小子域?yàn)V波法[J]. 石油地球物理勘探，1995, 30(2): 240-244.

Yang Gaoyin. A New Technique for Potential-Field Data Processing: Small Sub-Domain Filtering[J]. Oil Geophysical Prospecting, 1995, 30(2): 240-244.

[2]Cooper G R J, Cowan D R. Edge Enhancement of Potential-Field Data Using Normalized Statistics[J]. Geophysics, 2008, 73(3): 1-4.

[3]Hood P, McClure D J. Gradient Measurements in Ground Magnetic Prospecting[J]. Geophysics, 1965, 30(3): 403-410.

[4] Cordell L. Gravimetric Expression of Graben Faulting in Santa Fe Country and the Espanola Basin, New Mexico[C]//New Mexico Geological Society Guidebook, 30th Annual Fall Field Conference. New Mexico: New Mexico Geological, Society, 1979: 59-64.

[5] Nabighian M N. The Analytic Signal of Two-Dimen-sional Magnetic Bodies With Polygonal Cross-Section: Its Properties and Use for Automated Anomaly Interpretation[J]. Geophysics, 1972, 37(3): 507-517.

[6] Miller H G, Singh V. Potential Field Tilt: A New Concept for Location of Potential Field Sources[J]. Journal of Applied Geophysics, 1994, 32(2/3): 213-217.

[7] Wijins C, Perez C, Kowalczyk P. Theta Map: Edge Detection in Magnetic Data[J]. Geophysics, 2005, 70(4): 39-43.

[8] 馬國慶，杜曉娟，李麗麗. 利用水平與垂直導(dǎo)數(shù)的相關(guān)系數(shù)進(jìn)行位場(chǎng)數(shù)據(jù)的邊界識(shí)別[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2011, 41(增刊1): 345-348.

Ma Guoqing, Du Xiaojuan, Li Lili. Edge Detection of Potential Filed Data Using Correlation Coefficients of Horizontal and Vertical Derivatives[J]. Journal of Jilin University (Earth Science Edition), 2011, 41(Sup.1): 345-348.

[9] 馬國慶，黃大年，于平，等. 改進(jìn)的均衡濾波器在位場(chǎng)數(shù)據(jù)邊界識(shí)別中的應(yīng)用[J]. 地球物理學(xué)報(bào)，2012, 55(12): 4288-4295.

Ma Guoqing, Huang Danian, Yu Ping, et al. Application of Improved Balancing Filters to Edge Identification of Potential Field Data[J]. Chinese Journal of Geophysics, 2012, 55(12): 4288-4295.

[10] 徐亞，郝天珧，周立宏，等. 位場(chǎng)小波變換研究進(jìn)展[J]. 地球物理學(xué)進(jìn)展，2006, 21(4): 1132-1138.

Xu Ya, Hao Tianyao, Zhou Lihong, et al. Review of Wavelet Transform in Potential Field[J]. Progress in Geophysics, 2006, 21(4): 1132-1138.

[11] 劉彩云. 基于小波變換的位場(chǎng)場(chǎng)源識(shí)別與異常分離方法研究[D]. 北京：中國地質(zhì)大學(xué)，2014.

Liu Caiyun. The Research on Methods of Edge Detection and Anomaly Separation of Potential Field Source Based on Wavelet Transform[D]. Beijing: China University of Geosciences, 2014.

[12] 張郝. 基于小波變換的圖像去噪方法研究[D].北京： 北京交通大學(xué)，2008.

Zhang Hao. The Research on Image Denoising Based on Wavelet Transform[D].Beijing: Beijing Jiaotong University, 2008.

[13] 嚴(yán)奉霞. 小波理論及其在圖像邊緣檢測(cè)中的應(yīng)用[D]. 長沙：國防科學(xué)技術(shù)大學(xué)，2003.

Yan Fengxia. Wavelet Theory and Its Application to The Edge Detection of Images[D].Changsha: National University of Defense Technology, 2003.

[14] 張代磊，黃大年，張沖. 基于遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡(luò)在密度界面反演中的應(yīng)用[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2017, 47(2): 580-588.

Zhang Dailei, Huang Danian, Zhang Chong. Application of BP Neural Network Based on Genetic Algorithm in the Inversion of Density Interface[J]. Journal of Jilin University (Earth Science Edition), 2017, 47(2): 580-588.

[15] 王明，郭志宏，何輝，等. 基于反雙曲正切的位場(chǎng)邊界識(shí)別技術(shù)[J]. 物探與化探，2013, 37(4): 655-663.

Wang Ming, Guo Zhihong, He Hui, et al. Edge Detection of Potential Field Data Using Inverse Hyperbolic Tangent[J]. Geophysical and Geochemical Exploration, 2013, 37(4): 655-663.

[16] 冷眉. 重磁資料異常分離與構(gòu)造邊界識(shí)別方法的研究及應(yīng)用[D]. 成都：成都理工大學(xué)，2014.

Leng Mei. Research and Application of the Method of Separating Anomalous and the Structure Boundary Identification with the Data of Gravity and Magnetic[D].Chengdu: Chengdu University of Technology, 2014.