基于三維構(gòu)造張量的位場邊界識(shí)別濾波器

周帥， 黃大年， 焦健

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

?

基于三維構(gòu)造張量的位場邊界識(shí)別濾波器

周帥， 黃大年*， 焦健

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

地質(zhì)體構(gòu)造邊界位置的確定是位場數(shù)據(jù)解釋中的一項(xiàng)重要工作，現(xiàn)有很多基于位場梯度張量數(shù)據(jù)的邊界檢測濾波器，但存在識(shí)別邊界位置模糊且無法均衡深淺地質(zhì)體異常的缺點(diǎn).本文定義了位場數(shù)據(jù)的三維構(gòu)造張量，并提出基于位場構(gòu)造張量的邊界濾波器.為了同時(shí)顯示不同振幅異常的邊界位置，對(duì)新定義的濾波器進(jìn)行歸一化處理.在高階均衡濾波器的計(jì)算中需要計(jì)算位場的垂向高階導(dǎo)數(shù)，本文引入一種計(jì)算的穩(wěn)定算法，基于拉普拉斯方程利用位場水平導(dǎo)數(shù)求解垂向?qū)?shù)，可減小垂向?qū)?shù)計(jì)算中產(chǎn)生的誤差.將定義的濾波器應(yīng)用到合成的重磁數(shù)據(jù)中證明了新方法相比傳統(tǒng)的濾波器能更加清晰、準(zhǔn)確地圈定邊界位置，而且針對(duì)同時(shí)含有正負(fù)異常的地質(zhì)情況，可避免產(chǎn)生額外的錯(cuò)誤邊界.最后將新的濾波器應(yīng)用到實(shí)測的重磁數(shù)據(jù)的解釋中，結(jié)果顯示基于構(gòu)造張量的濾波器可更準(zhǔn)確清晰地劃分出斷裂的邊界位置，發(fā)現(xiàn)更多的構(gòu)造細(xì)節(jié).

構(gòu)造張量； 邊界識(shí)別； 位場

1　引言

位場數(shù)據(jù)的邊界檢測是地球物理勘探中廣泛應(yīng)用的方法，它可以圈定地質(zhì)體的水平位置，指導(dǎo)下一步的解釋工作.目前，有很多方法用于地質(zhì)構(gòu)造邊界的檢測和拾取中.Cordell(1979)利用總水平導(dǎo)數(shù)的最大值圈定地質(zhì)體的邊界.Hood 和Teskey(1989)發(fā)現(xiàn)位場垂向?qū)?shù)的零等值線與地質(zhì)體的邊界相對(duì)應(yīng).Nabighian(1972)提出構(gòu)造邊界的位置可用解析信號(hào)振幅的最大值描述，后來，Hsu等(1996)利用高階解析信號(hào)圈定的磁性體的邊界位置更加清晰，提高了分辨率.但是上述方法均無法同時(shí)清晰地識(shí)別深部和淺部的重磁異常體的邊界.為了解決這一問題，人們開始致力于均衡邊界濾波器的研究，主要通過導(dǎo)數(shù)間的比值實(shí)現(xiàn).Wang等(2010)分析了簡單規(guī)則形體重力異常垂向?qū)?shù)零值位置的空間變化規(guī)律，垂向?qū)?shù)可圈定地質(zhì)體的邊界位置，且可用于均衡強(qiáng)弱振幅地質(zhì)體異常.Cooper和Cowan(2006)、馬國慶等(2012)利用垂向?qū)?shù)定義均衡濾波器，極大值位置對(duì)應(yīng)地質(zhì)體邊界位置.Miller和Singh(1994)利用位場總水平導(dǎo)數(shù)與垂向?qū)?shù)的比值定義了均衡邊界識(shí)別濾波器傾斜角(tilt angle).Verduzco和Fairhead(2004)計(jì)算傾斜角濾波器的總水平導(dǎo)數(shù)，并定義了新的高階邊界識(shí)別方法Thdr.Wijns等(2005)利用解析信號(hào)振幅均衡總水平導(dǎo)數(shù)，定義了新的濾波器Theta圖，它的最大值位置可有效的同時(shí)圈定深淺地質(zhì)體的邊界位置.王萬銀等(2010)對(duì)位場邊界識(shí)別的方法進(jìn)行綜述，指出位場傾斜角法和Theta圖法邊界檢測結(jié)果與垂向一階導(dǎo)數(shù)結(jié)果一致，且存在“解析奇點(diǎn)”的缺點(diǎn).Ma(2013)利用不同階導(dǎo)數(shù)的比值定義了新的邊界識(shí)別器，使得圈定的邊界位置更加清晰.Zhou等(2013)基于重力梯度張量矩陣的特征值定義了新的邊界識(shí)別器，解決了含正負(fù)重力異常的復(fù)雜地質(zhì)條件的邊界識(shí)別問題，避免產(chǎn)生錯(cuò)誤的邊界.

構(gòu)造張量被廣泛應(yīng)用于圖像處理與特征拾取中.Weickert(1999)提出構(gòu)造張量可用于確定圖像的邊界和角點(diǎn)位置.三維構(gòu)造張量可有效地描述視頻中的局部運(yùn)動(dòng)細(xì)節(jié)，并被廣泛應(yīng)用于視頻的剪切中(Wang and Ma, 2003).Jeong 等(2006)將構(gòu)造張量應(yīng)用到地球物理解釋中，利用構(gòu)造張量實(shí)現(xiàn)地震圖像斷層的識(shí)別.Sertcelik和Kafadar(2012)將高斯函數(shù)卷積二維構(gòu)造張量的特征值應(yīng)用于重力數(shù)據(jù)的邊界識(shí)別中，但此方法無法均衡深淺地質(zhì)體異常.Ma等(2015)利用構(gòu)造張量的水平和垂直方向?qū)?shù)定義新的均衡濾波器，解決了無法識(shí)別深部異常的問題.Yuan等(2014)指出與高斯函數(shù)進(jìn)行卷積而定義的基于構(gòu)造張量的邊界識(shí)別器會(huì)模糊圈定的邊界位置，他們定義新的濾波器，可更清晰地顯示地質(zhì)體位置.

本文定義了基于三維構(gòu)造張量的邊界識(shí)別濾波器，可清晰準(zhǔn)確地顯示不同振幅大小的邊界.而且針對(duì)同時(shí)含正負(fù)異常的復(fù)雜地質(zhì)情況，本文方法可避免產(chǎn)生多余的干擾邊界.

2　基于三維構(gòu)造張量的濾波器

2.1三維構(gòu)造張量(ST)

Sertcelik和Kafadar(2012)定義了含高斯函數(shù)的二維構(gòu)造張量，其表達(dá)式為：

(1)

(2)

(3)

位場三維構(gòu)造張量矩陣與位場梯度張量矩陣有相似的表達(dá)形式，現(xiàn)在有很多基于位場梯度張量的邊界識(shí)別方法(Pedersen and Rasmussen，1990；Hansen and deRidder，2006；Beiki and Pedersen，2010；Beiki et al.，2012；Pilkington and Beiki，2013).本文對(duì)位場的構(gòu)造張量進(jìn)行分析，并定義了基于位場構(gòu)造張量的邊界識(shí)別濾波器.

2.2新提出的邊界識(shí)別濾波器

我們利用位場構(gòu)造張量(公式(3))的每一列定義表達(dá)式：

(4)

(5)

(6)

其中，STx和STy最大值位置分別描述地質(zhì)體的N—S和E—W方向邊界，STz極小值位置與地質(zhì)體的邊界位置對(duì)應(yīng)，我們建立一棱柱體模型說明其效果.圖1i所示為棱柱體產(chǎn)生的重力異常，虛線為棱柱體的水平位置，棱柱體頂部埋深1 km,厚度0.5 km，剩余密度100 kg·m-3.圖1a—1b顯示STx和STy可圈定出N—S和E—W方向的邊界位置，STz的局部極小值位置描述出地質(zhì)體的邊界位置.利用STx和STy定義一種基于構(gòu)造張量的邊界濾波器：

(7)

SED的最大值對(duì)應(yīng)地質(zhì)體的邊界，如圖1g所示.對(duì)位場數(shù)據(jù)求垂向?qū)?shù)可提高分辨率(Hsu et al.，1996)，我們對(duì)公式(3)中的原始位場數(shù)據(jù)f求垂向一階導(dǎo)數(shù)，并定義對(duì)應(yīng)的邊界濾波器：

(8)

(9)

(10)

圖1　基于構(gòu)造張量的邊界濾波器效果(a) 異常的STx結(jié)果; (b) 異常的STy結(jié)果; (c) 異常的STz結(jié)果; (d) 異常的FSTx結(jié)果; (e) 異常的FSTy結(jié)果; (f) 異常的FSTz結(jié)果; (g) 異常的SED結(jié)果; (h) 異常的FSED結(jié)果; (i) 合成的重力異常.Fig.1　Results of edge detectors based on structural tensors(a) STx of anomaly; (b) STy of anomaly; (c) STz of anomaly; (d) FSTx of anomaly; (e) FSTy of anomaly; (f) FSTz of anomaly; (g) SED of anomaly; (h) FSED of anomaly; (i) Synthetic gravity anomaly.

圖2　不同方法邊界識(shí)別結(jié)果的剖面圖(a) 合成的重力異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果.Fig.2　Profiles showing edge detecting results using different methods(a) Synthetic gravity anomaly; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly.

圖3　不同方法邊界識(shí)別結(jié)果(a) 模型1合成的重力異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果.Fig.3　Edge detecting results using different methods(a) Synthetic gravity anomaly of model 1; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly.

(11)

圖1d—1f分別為FSTx、FSTy、FSTz邊界識(shí)別效果，F(xiàn)STx和FSTy相比低階的STx和STy可更清楚地圈定出N—S和E—W方向的邊界位置，對(duì)應(yīng)的邊界識(shí)別濾波器FSED(圖1h)相比SED具有更高的邊界識(shí)別分辨率.

2.3均衡的邊界識(shí)別濾波器

建立如圖2a所示的模型，埋深分別為1 km和2 km.圖2b和2c分別為新提出的邊界識(shí)別濾波器SED和FSED結(jié)果，利用垂向?qū)?shù)定義的濾波器FSED相比SED具有更高的分辨率，可更清晰地顯示邊界位置，但兩種方法均存在無法均衡深淺模型異常的缺點(diǎn).為了同時(shí)清晰地識(shí)別不同振幅的地質(zhì)體異常，我們提出改進(jìn)的基于位場構(gòu)造張量的邊界識(shí)別濾波器.STz的局部極小值與邊界位置對(duì)應(yīng)，我們采用STz對(duì)新定義的濾波器SED進(jìn)行歸一化處理，表達(dá)式為：

(12)

(13)

其中，R表示取實(shí)部運(yùn)算；常數(shù)k用于避免在正負(fù)異常同時(shí)存在情況下產(chǎn)生額外錯(cuò)誤的邊界，取值一般為0到0.1，當(dāng)只存在正異常或者負(fù)異常時(shí)，k取0，正負(fù)異常同時(shí)存在時(shí)，選擇一定數(shù)值的k可消除假異常邊界.k取值較大時(shí)會(huì)降低強(qiáng)弱振幅異常的均衡效果，對(duì)弱異常具有一定的壓制作用，造成識(shí)別的弱異常邊界位置不清晰；k取值較小時(shí)無法有效地消除在正負(fù)異常同時(shí)存在時(shí)帶來的錯(cuò)誤邊界信息.實(shí)際計(jì)算中需要解釋人員嘗試不同的k值完成邊界識(shí)別濾波器的計(jì)算，進(jìn)而確定最優(yōu)的邊界識(shí)別結(jié)果.圖2d和2e分別為NSED1和NSED2的邊界識(shí)別結(jié)果，SED和FSED濾波器均可平衡深淺地質(zhì)體異常.而且NSED2相比NSED1具有更高的分辨率，更清晰地顯示邊界位置.

邊界識(shí)別濾波器FSED的均衡的邊界濾波器表達(dá)式為：

(14)

NFSED2=

(15)

其中,

(16)

p用于平衡不同階導(dǎo)數(shù)間的強(qiáng)度，使均衡的邊界濾波器滿足數(shù)學(xué)運(yùn)算法則.為了說明本文方法的有效性，選擇常用的傾斜角總水平導(dǎo)數(shù)Thdr和Theta圖方法進(jìn)行效果對(duì)比.Verduzco和Fairhead(2004)定義了傾斜角的總水平導(dǎo)數(shù)，表達(dá)式為：

(17)

其中，

(18)

Wijns等(2005)利用解析信號(hào)的振幅均衡總水平導(dǎo)數(shù)，定義了均衡的邊界濾波器Theta圖：

(19)

其中，f表示原始位場異常.圖2h和2i分別為Thdr和Theta圖濾波器邊界識(shí)別效果，圖2f—2g表明NFSED1和NFSED2最大值顯示模型體的邊界.相比傳統(tǒng)的邊界識(shí)別器，新定義的方法具有更高的水平分辨率.

在均衡邊界識(shí)別器的計(jì)算中，需要計(jì)算位場f的高階導(dǎo)數(shù)，頻率域內(nèi)的導(dǎo)數(shù)計(jì)算會(huì)放大噪聲的干擾.為了減小噪聲干擾，f的各階水平導(dǎo)數(shù)在空間域采用有限差分的方式計(jì)算，而各階垂向?qū)?shù)采用下列變換進(jìn)行計(jì)算：

(20)

(21)

(22)

(23)

其中，fz的計(jì)算分兩步進(jìn)行，首先通過頻率域內(nèi)積分計(jì)算得到原始重磁異常f 的標(biāo)量位U，進(jìn)而利用空間域方法得到位U的二階水平導(dǎo)數(shù)Uxx和Uyy，然后利用拉普拉斯方程求得fz(FediandFlorio,2002).頻率域內(nèi)積分和空間域水平導(dǎo)數(shù)計(jì)算不會(huì)放大噪聲干擾，fz計(jì)算的表達(dá)式為：

(24)

3　理論模型試驗(yàn)

為了驗(yàn)證方法的有效性，建立如圖3a所示的模型1，虛線標(biāo)識(shí)地質(zhì)體的真實(shí)水平位置，棱柱體1和3的頂面埋深為1km，棱柱體2和4的埋深為0.5km，剩余密度均為100kg·m-3.采用上述的方法對(duì)該模型重力異常進(jìn)行處理，得到如圖3b—3i的邊界識(shí)別結(jié)果.SED和FSED濾波器可清晰地顯示淺部棱柱體的邊界位置，圈定的深部地質(zhì)體邊界位置模糊.傳統(tǒng)的濾波器Theta圖可均衡深淺異常，但是相比本文提出的均衡濾波器NSED1和NSED2，識(shí)別的邊界位置比較發(fā)散.對(duì)比同階導(dǎo)數(shù)的傳統(tǒng)濾波器Thdr和NFSED1、NFSED2，傳統(tǒng)的高階導(dǎo)數(shù)濾波器Thdr識(shí)別的邊界位置模糊，而本文新提出的方法可清晰準(zhǔn)確地顯示深淺異常，且NFSED2比NFSED1具有更高的邊界識(shí)別分辨率.

位場實(shí)測數(shù)據(jù)中不可避免的含有噪聲干擾，我們對(duì)模型2對(duì)應(yīng)的重力異常中加入3%的高斯噪聲.圖5a為加入噪聲后的重力異常，虛線為真實(shí)的地質(zhì)體位置.圖5h和5i所示傳統(tǒng)濾波器Thdr和Theta圖受噪聲影響較大，無法顯示棱柱體的邊界位置.SED和FSED受噪聲干擾影響較小，但無法均衡深淺地質(zhì)體異常.利用高階導(dǎo)數(shù)定義的濾波器NFSED1、NFSED2相比NSED1、NSED2放大了噪聲干擾，但可更清晰地顯示邊界位置，具有一定的抗噪性.并且NFSED2相比同階的NFSED1抗噪能力強(qiáng)，更加準(zhǔn)確地圈定出地下不同埋深異常體的水平邊界位置.在均衡濾波器計(jì)算中，因子k選擇為0.01.

為了檢驗(yàn)本文提出的均衡濾波器NSED1、NSED2、NFSED1、NFSED2在場源埋深差異較大情況下的邊界識(shí)別效果，我們將圖3a中的1和3號(hào)棱柱體的頂部埋深增大為2km，2和4號(hào)棱柱體頂部埋深仍為0.5km，剩余密度均為100kg·m-3，得到圖6a—6d所示四個(gè)均衡濾波器的邊界識(shí)別結(jié)果.改變1和3號(hào)棱柱體的剩余密度為-100kg·m-3，棱柱體其他參數(shù)不變，圖6e—6h為對(duì)應(yīng)的NSED1、NSED2、NFSED1、NFSED2邊界識(shí)別效果.結(jié)果顯示當(dāng)深淺場源埋深四倍差別時(shí)，正負(fù)異常同時(shí)存在時(shí)可獲得可靠的邊界識(shí)別結(jié)果.繼續(xù)增大1和3號(hào)棱柱體的頂部埋深為3km，剩余密度均為100kg·m-3，四個(gè)邊界識(shí)別器在只存在正密度異常時(shí)可圈定出清晰的邊界位置，可見圖6i—6l.但是由于1和3號(hào)棱柱體埋藏較深，識(shí)別的邊界位置與實(shí)際位置有一定的偏移.保持1和3號(hào)棱柱體的頂部埋深為3km，剩余密度改變?yōu)?100kg·m-3，圖6m—6n顯示NSED1和NSED2濾波器無法清晰地均衡深淺地質(zhì)體異常.NFSED1和NFSED2識(shí)別的1和3號(hào)棱柱體的邊界位置間斷，如圖6o和6p所示.當(dāng)場源埋深差異較大時(shí)，針對(duì)只存在正異常或負(fù)異常的地質(zhì)情況，均衡濾波器可獲得清晰的場源邊界位置，但圈定的邊界位置與真實(shí)位置有一定偏差.對(duì)于正負(fù)異常同時(shí)存在的地質(zhì)情況，新提出的邊界識(shí)別器在六倍場源埋深差異時(shí)，圈定的邊界位置出現(xiàn)間斷，場源之間的空間位置臨近關(guān)系也影響濾波器的邊界識(shí)別結(jié)果.

圖4　不同方法邊界識(shí)別結(jié)果(a) 模型2合成的重力異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果.Fig.4　Edge detecting results using different methods(a) Synthetic gravity anomaly of model 2; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly.

圖5　不同方法邊界識(shí)別結(jié)果(a) 模型2加入3%高斯噪聲合成的重力異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果.Fig.5　Edge detecting results using different methods(a) Synthetic gravity anomaly of model 2 with added 3% Gaussian noise; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly.

圖6　不同場源埋深均衡濾波器邊界識(shí)別結(jié)果埋深為2 km正異常時(shí)： (a)NSED1; (b) NSED2; (c) NFSED1; (d) NFSED2. 埋深為2 km負(fù)異常時(shí)：(e)NSED1; (f) NSED2; (g) NFSED1; (h) NFSED2. 埋深為3 km正異常時(shí)：(i)NSED1; (j) NSED2; (k) NFSED1; (l) NFSED2. 埋深為3 km負(fù)異常時(shí)：(m)NSED1; (n) NSED2; (o) NFSED1; (p) NFSED2.Fig.6　Edge detecting results using balanced filters with different source buried depthsPositive gravity anomaly with buried depth 2 km: (a) NSED1; (b) NSED2; (c) NFSED1; (d) NFSED2. Negative gravity anomaly with buried depth 2 km: (e) NSED1; (f) NSED2; (g) NFSED1; (h) NFSED2. Positive gravity anomaly with buried depth 3 km: (i) NSED1; (j) NSED2; (k) NFSED1; (l) NFSED2. Negative gravity anomaly with buried depth 3 km: (m) NSED1; (n) NSED2; (o) NFSED1; (p) NFSED2.

圖7　吉林某地區(qū)重力異常邊界識(shí)別結(jié)果(a) 實(shí)測重力異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果; (j) 異常的ETHD結(jié)果; (k) 異常的N_E結(jié)果; (l) 異常的N_HE圖結(jié)果.Fig.7　Edge detecting results of gravity anomalies in an area of Jilin using different methods(a) Measured gravity anomaly; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly; (j) ETHD of anomaly; (k) N_E of anomaly; (l) N_HE of anomaly.

圖8　朱日和地區(qū)化極磁異常邊界識(shí)別結(jié)果(a) 化極后磁異常; (b) 異常的SED結(jié)果; (c) 異常的FSED結(jié)果; (d) 異常的NSED1結(jié)果; (e) 異常的NSED2結(jié)果; (f) 異常的NFSED1結(jié)果; (g) 異常的NFSED2結(jié)果; (h) 異常的Thdr結(jié)果; (i) 異常的Theta圖結(jié)果; (j) 異常的ETHD結(jié)果; (k) 異常的N_E結(jié)果; (l) 異常的N_HE圖結(jié)果.Fig.8　Edge detecting results of magnetic anomalies reduced to the pole in the Zhurihe area(a) Magnetic anomaly reduced to the pole; (b) SED of anomaly; (c) FSED of anomaly; (d) NSED1 of anomaly; (e) NSED2 of anomaly; (f) NFSED1 of anomaly; (g) NFSED2 of anomaly; (h) Thdr of anomaly; (i) Theta of anomaly; (j) ETHD of anomaly; (k) N_E of anomaly; (l) N_HE of anomaly.

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

我們將本文所述方法應(yīng)用到實(shí)測的重力和磁異常數(shù)據(jù)中驗(yàn)證方法的實(shí)際應(yīng)用效果，對(duì)吉林某地區(qū)鐵礦的剩余重力異常和朱日和地區(qū)的磁異常數(shù)據(jù)進(jìn)行邊界識(shí)別.圖7a為吉林某地區(qū)的重力異常，采樣間距0.1km，同時(shí)包含正負(fù)密度異常，該區(qū)主要以中心位置磁鐵礦及黃鐵礦化的矽卡巖引起的正異常為主，圖中白線為根據(jù)前期工作劃定的成礦帶的大致范圍.Zhou和Huang(2015)基于張量矩陣特征值定義了均衡的邊界識(shí)別濾波器N_E和N_HE.王賽昕和劉林靜(2011)提出了均值歸一總水平導(dǎo)數(shù)邊界識(shí)別方法ETHD，可同時(shí)顯示深淺地質(zhì)體異常的邊界位置.我們引入N_E、N_HE和ETHD三種濾波器進(jìn)一步對(duì)比分析本文方法的邊界識(shí)別效果.從圖7b—7l不同方法的邊界識(shí)別結(jié)果可以看出，SED和FSED非均衡濾波器無法反映深部異常，識(shí)別的邊界位置模糊.在計(jì)算均衡的濾波器NSED1和NSED2時(shí)，選擇因子k等于0.0005，它們可同時(shí)顯示深淺異常的邊界位置，但分辨率差.高階導(dǎo)數(shù)定義的NFSED1和NFSED2可清晰準(zhǔn)確地圈定不同剩余密度的地質(zhì)體構(gòu)造邊界，在計(jì)算中因子k選擇為0.0005.傳統(tǒng)濾波器Thdr識(shí)別效果差，無法清晰顯示邊界位置，Theta圖可均衡深淺異常，但由于引入了額外的錯(cuò)誤邊界，造成邊界相互連接，給解釋帶來了很大的誤導(dǎo).N_E和ETHD可穩(wěn)定地圈定出異常體的邊界位置，但相比本文提出的均衡濾波器，它們識(shí)別的邊界位置模糊.當(dāng)數(shù)據(jù)中含有一定噪聲干擾時(shí)，高階濾波器N_HE由于需要計(jì)算位場的垂向高階導(dǎo)數(shù)，造成計(jì)算結(jié)果的不穩(wěn)定.

由于磁異常受磁化方向的干擾，為了獲得可靠的解釋結(jié)果，在應(yīng)用本文方法進(jìn)行磁異常處理時(shí)，需要對(duì)原始磁異常進(jìn)行化極處理.圖8a為朱日和地區(qū)化極后的磁異常圖，采樣間隔為20m，測量區(qū)域?yàn)樾略糯鷮有?，主要由大陸沉積物構(gòu)成，除一些含鐵的沙堰堤壩外基本不含鐵磁性物質(zhì).磁異常數(shù)據(jù)主要由堤壩產(chǎn)生的近線性異常和鐵異常構(gòu)成.本文提出的均衡的邊界濾波器相比傳統(tǒng)的濾波器Thdr和Theta圖具有更高的分辨率，且可避免產(chǎn)生由正負(fù)異常存在時(shí)的錯(cuò)誤邊界.N_E、N_HE和ETHD可大致顯示出構(gòu)造的邊界位置，但識(shí)別的邊界位置不清晰.其中，利用公式(12)—(15)計(jì)算時(shí)，k選擇為0.0001.

5　結(jié)論

本文定義了位場的構(gòu)造張量矩陣，對(duì)矩陣進(jìn)行分析并建立了新的邊界識(shí)別濾波器.為了克服濾波器無法均衡深淺異常的缺點(diǎn)，對(duì)其進(jìn)行歸一化處理，定義了均衡的邊界識(shí)別濾波器.在同時(shí)含有正負(fù)異常的復(fù)雜情況下，傳統(tǒng)的濾波器會(huì)產(chǎn)生錯(cuò)誤的干擾邊界信息，通過引入因子k消除了這一干擾，但當(dāng)場源埋深差別較大時(shí)會(huì)造成識(shí)別的深部場源邊界間斷.將基于位場構(gòu)造張量的邊界識(shí)別濾波器應(yīng)用到合成的和實(shí)測的重磁異常中，新定義的均衡濾波器相比傳統(tǒng)濾波器Thdr和Theta圖在均衡深淺異常的同時(shí)，更清晰準(zhǔn)確地圈定構(gòu)造邊界的水平位置，發(fā)現(xiàn)更多的小構(gòu)造信息，更有利于指導(dǎo)下一步的勘探工作.

致謝感謝吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院馬國慶副教授提供的實(shí)測數(shù)據(jù)及在論文寫作中給予的幫助.感謝審稿專家及編輯對(duì)本文提出的寶貴的修改意見.

Beiki M, Pedersen L B. 2010. Eigenvector analysis of gravity gradient tensor to locate geologic bodies.Geophysics, 75(6): 137-149.

Beiki M, Clark D A, Austin J R, et al. 2012. Estimating source location using normalized magnetic source strength calculated from magnetic gradient tensor data.Geophysics, 77(6): J23-J37. Cooper G R J, Cowan D R. 2006. Enhancing potential field data using filters based on the local phase.Computers&Geosciences, 32(10): 1585-1591. Cordell L. 1979. Gravimetric expression of graben faulting in Santa Fe Country and the Espanola Basin, New Mexico∥Ingersoll R V ed. Guidebook to Santa Fe Country, 30th Field Conference. New Mexico Geological Society Guidebook. Virginia, USA: 59-64.Fedi M, Florio G. 2002. A stable downward continuation by using the ISVD method.GeophysicalJournalInternational, 151: 146-156. Hansen R O, deRidder E. 2006. Linear feature analysis for aeromagnetic data.Geophysics, 71(6): L61-L67. Hood P J, Teskey D J. 1989. Aeromagnetic gradiometer program of the geological survey of Canada.Geophysics, 54(8): 1012-1022. Hsu S H, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential-field anomalies: An enhanced analytic signal technique.Geophysics, 61(2): 373-386.

Jeong W K, Whitaker R, Dobin M. 2006. Interactive 3D seismic fault detection on the Graphics Hardware.∥Proceedings of the 2006 International Workshop on Volume Graphics. Boston, USA: 111-118.

TCP采用網(wǎng)絡(luò)服務(wù)器和客戶端的方式進(jìn)行通信。服務(wù)器主要負(fù)責(zé)發(fā)送數(shù)據(jù)至多個(gè)客戶端或者從多個(gè)客戶端讀取數(shù)據(jù)，同樣，客戶端可以向所對(duì)應(yīng)的服務(wù)器發(fā)送請(qǐng)求并獲取所需數(shù)據(jù)。服務(wù)器與客戶端傳輸是雙向的，區(qū)別在于服務(wù)器可同時(shí)連接多個(gè)客戶端。

Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data.ChineseJournalofGeophysics(in Chinese), 55(12): 4288-4295, doi: 10.6038/j.issn.0001-5733.2012.12.040.

Ma G Q. 2013. Edge detection of potential field data using improved local phase filter.ExplorationGeophysics, 44(1): 36-41.

Ma G Q, Liu C, Huang D N. 2015. The removal of additional edges in the edge detection of potential field data.JournalofAppliedGeophysics, 114: 168-173.

Miller H G, Singh V. 1994. Potential field tilt—a new concept for location of potential field sources.JournalofAppliedGeophysics, 32(2-3): 213-217. Nabighian M N. 1972. The analytic signal of two-dimensional magnetic bodies with polygonal cross-section: Its properties and use for automated anomaly interpretation.Geophysics, 37(3): 507-517.

Pedersen L B, Rasmussen T M. 1990. The gradient tensor of potential field anomalies: Some implications on data collection and data processing of maps.Geophysics, 55(12): 1558-1566.

Pilkington M, Beiki M. 2013. Mitigating remanent magnetization effects in magnetic data using the normalized source strength.Geophysics, 78(3): J25-J32.

Sertcelik I, Kafadar O. 2012. Application of edge detection to potential field data using eigenvalue analysis of structure tensor.JournalofAppliedGeophysics, 84: 86-94.

Verduzco B, Fairhead J D. 2004. New insights into magnetic interpretation using the 3-D analytic signal.TheLeadingEdge, 23(2): 116-119.

Wang H Y, Ma K K. 2003. Automatic video object segmentation via 3D structure tensor.∥Proceedings of 2003 International Conference on Image Processing. Barcelona, Spain: IEEE.

Wang S X, Liu L J. 2011. Edge detection by equalized and normalized amplitude of total horizontal derivatives.ChineseJournalofEngineeringGeophysics(in Chinese), 8(6): 699-704. Wang W Y, Zhang G C, Liang J S. 2010. Spatial variation law of vertical derivative zero points for potential field data.AppliedGeophysics, 7(3): 197-209.

Wang W Y, Qiu Z Y, Yang Y, et al. 2010. Some advances in the edge recognition of the potential field.ProgressinGeophysics(in Chinese), 25(1): 196-210, doi: 10.3969/j.issn.1004-2903.2010.01.027.

Weickert J. 1999. Coherence-enhancing diffusion of color images.ImageandVisionComputing, 17(3-4): 201-212.

Wijns C, Perez C, Kowalczyk P. 2005. Theta map: Edge detection in magnetic data.Geophysics, 70(4): L39-L43.

Yuan Y, Huang D N, Yu Q L, et al. 2014. Edge detection of potential field data with improved structure tensor methods.JournalofAppliedGeophysics, 108: 35-42. Zhou S, Huang D N. 2015. Edge detection using directional eigenvalues of potential field gradient tensor data.GlobalGeology, 18(3): 188-195. Zhou W N, Du X J, Li J Y. 2013. The limitation of curvature gravity gradient tensor for edge detection and a method for overcoming it.JournalofAppliedGeophysics, 98: 237-242.

附中文參考文獻(xiàn)

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

王賽昕, 劉林靜. 2011. 均值歸一總水平導(dǎo)數(shù)邊界識(shí)別方法. 工程地球物理學(xué)報(bào), 8(6): 699-704.

王萬銀, 邱之云, 楊永等. 2010. 位場邊緣識(shí)別方法研究進(jìn)展. 地球物理學(xué)進(jìn)展, 25(1): 196-210, doi: 10.3969/j.issn.1004-2903.2010.01.027.

(本文編輯何燕)

A filter to detect edge of potential field data based on three-dimensional structural tensors

ZHOU Shuai, HUANG Da-Nian*， JIAO Jian

CollegeofGeoexplorationScienceandTechnology,JilinUniversity,Changchun130026,China

To detect the edge of potential field is an important task in identifying the horizontal position of geological structure. Although there are many methods based on potential field gradient tensors used to edge detection, the resolved edges remain blurred and inaccurate. This paper defines the structural tensors of a potential field and develops the related edge detectors to outline sources edges. In order to display different-amplitude anomalies simultaneously, we purpose two normalization strategies. When we calculate the new edge detectors, high-order vertical derivatives of potential field computation can increase the noise effect. The new calculation algorithm using horizontal derivatives to calculate vertical derivatives based on Laplace equation is applied to solve this problem. The presented filters tested on a synthetic potential field model show that they can outline the source edges clearly and accurately compared with traditional detectors. For complex situations with positive and negative anomalies, the new methods will not cause extra false edges. At last, we apply these new detectors to the real measured gravity and magnetic data. The results indicate that the new methods have ability to display edges of geological structures clearly and accurately and to discover more subtle information.KeywordsStructural tensor; Edge detection; Potential field

10.6038/cjg20161028.Zhou S, Huang D N， Jiao J. 2016. A filter to detect edge of potential field data based on three-dimensional structural tensors.ChineseJ.Geophys. (in Chinese),59(10):3847-3858,doi:10.6038/cjg20161028.

國家高技術(shù)研究發(fā)展計(jì)劃(863計(jì)劃)課題(2014AA06A613)資助.

周帥，男，1988年生，博士研究生，主要從事位場數(shù)據(jù)處理及解釋方面的研究.E-mail:zhoushuai14@mails.jlu.edu.cn

黃大年，男，1958年生，教授，博士生導(dǎo)師，主要從事移動(dòng)平臺(tái)探測數(shù)據(jù)處理與解釋及一體化軟件平臺(tái)開發(fā).E-mail:dnhuang@jlu.edu.cn

10.6038/cjg20161028

P631

2015-12-13，2016-06-27收修定稿

周帥， 黃大年， 焦健. 2016. 基于三維構(gòu)造張量的位場邊界識(shí)別濾波器. 地球物理學(xué)報(bào)，59(10):3847-3858,