一種改進(jìn)的重力梯度全張量數(shù)據(jù)的邊界識(shí)別方法

杜威， 吳賀宇， 陳祥忠， 劉俊成， 張陽(yáng)陽(yáng)， 惠夢(mèng)琳

1 中國(guó)科學(xué)院地球化學(xué)研究所礦床地球化學(xué)重點(diǎn)研究實(shí)驗(yàn)室， 貴陽(yáng) 550081 2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 130026 3 北京桔燈地球物理勘探股份有限公司， 北京 102200 4 安徽省勘查技術(shù)院， 合肥 230031

0 引言

隨著重力梯度測(cè)量技術(shù)的發(fā)展，越來(lái)越多的重力梯度全張量數(shù)據(jù)被應(yīng)用到數(shù)據(jù)處理中.重力梯度全張量數(shù)據(jù)是重力位的二階偏導(dǎo)數(shù)，由9個(gè)張量數(shù)據(jù)組成，較重力數(shù)據(jù)包含有更高頻率的信息.對(duì)重力梯度全張量數(shù)據(jù)進(jìn)行處理，可以得到更全面，更準(zhǔn)確的地質(zhì)信息.

邊界識(shí)別在位場(chǎng)數(shù)據(jù)處理中有著重要的作用(Cooper and Cowan，2006；馬國(guó)慶等，2011；Archibald et al.，1999；魯寶亮等，2020)，其方法種類有很多，常用的有：垂向?qū)?shù)(Evjen，1936)、總水平導(dǎo)數(shù)法(THDR)(Verduzco et al.，2004)、解析信號(hào)(ASM)(Hsu et al.，1996)、傾斜角法(Miller and Singh，1994)、傾斜角水平導(dǎo)數(shù)法(王明等，2013)、Theta圖法(Wijns et al.，2005)以及其他方法(段杰翔和徐亞，2020；王丁丁等，2021)，這些方法都是在傳統(tǒng)的重力數(shù)據(jù)基礎(chǔ)上提出的.隨著重磁勘探方法的多樣化和高精度勘探儀器的發(fā)展，梯度儀可以測(cè)量所有的重磁張量分量.許多經(jīng)典的邊緣檢測(cè)方法被應(yīng)用于梯度張量數(shù)據(jù)(鄭強(qiáng)等，2019).Zhang等(2000)提出張量-歐拉反褶積方法，提高了地質(zhì)體水平位置的分辨率(Mikhailov et al.，2007).Beiki(2010a，b)利用重力張量數(shù)據(jù)的方向解析信號(hào)振幅來(lái)確定地質(zhì)體位置.朱自強(qiáng)等(2015)通過(guò)對(duì)重力梯度張量解析信號(hào)應(yīng)用歐拉反褶積，提高了反演結(jié)果的可靠性.Yuan和Yu(2015)定義了二階方向解析信號(hào)，提出了幾種基于水平方向解析信號(hào)和二階水平方向解析信號(hào)的新邊界識(shí)別方法.Oru?和Keskinsezer(2008)使用傾斜角法，利用重力梯度張量數(shù)據(jù)對(duì)簡(jiǎn)單模型的邊緣進(jìn)行檢測(cè).另外，重力梯度張量的特征值也被用來(lái)描繪地質(zhì)體的邊緣.Beiki和Pedersen(2010)利用重力梯度張量數(shù)據(jù)的特征向量來(lái)估計(jì)場(chǎng)源的位置和走向.Sertcelik和Kafadar(2012)提出了一種利用結(jié)構(gòu)張量特征值來(lái)識(shí)別地質(zhì)體邊緣的方法.Yuan等(2014a)提出了三種平衡結(jié)構(gòu)張量最大特征值的方法，并重新定義了結(jié)構(gòu)張量.Zhou等(2017)提出了基于三維結(jié)構(gòu)張量方向總水平導(dǎo)數(shù)的邊緣檢測(cè)新方法.Oru?等(2013)利用重力梯度張量數(shù)據(jù)曲率的最大特征值圈定了地質(zhì)構(gòu)造的邊緣.一些學(xué)者基于Oru? 等提出的方法進(jìn)行了相關(guān)改進(jìn)(Zhou et al.，2013；Wang et al.，2015；Zuo and Hu，2015).Yuan等(2016)提出了方向Theta方法，并利用水平方向Theta組合(ThetaX與ThetaY相加)的最大值定義了一種新的邊緣檢測(cè)方法(ED)來(lái)描繪重力梯度張量數(shù)據(jù)的邊界.當(dāng)?shù)刭|(zhì)體同時(shí)存在正異常和負(fù)異常時(shí)，該方法可以避免產(chǎn)生虛假邊緣，但該方法受地質(zhì)體走向的影響.本文針對(duì)其受地質(zhì)體走向影響的問(wèn)題，提出了一種改進(jìn)的邊界識(shí)別方法(IED)，替換原方法中的加法運(yùn)算，通過(guò)選擇合適的閾值來(lái)削弱虛假邊界的識(shí)別.模型試驗(yàn)表明，該方法不再受地質(zhì)體走向的影響，在實(shí)際數(shù)據(jù)處理中取得了良好的效果.

1 ED方法及存在問(wèn)題

重力梯度全張量(GGT)數(shù)據(jù)是重力位U的二階偏導(dǎo)數(shù)，其矩陣形式為

(1)

由位場(chǎng)理論可知，地質(zhì)體外的引力位滿足拉普拉斯方程(Beiki et al., 2010)：

gxx+gyy+gzz=0，

(2)

T為對(duì)稱矩陣，它只包括5個(gè)獨(dú)立張量，gxy、gxz、gyz、gxx和gzz.方向解析信號(hào)的振幅為

(3)

Yuan和Geng(2014b)、Yuan和Yu(2015)定義了方向總水平導(dǎo)數(shù).基于Theta方法的定義，提出了GGT數(shù)據(jù)的方向Theta方法(Yuan et al.，2016)，可以表示為

(4)

并在此基礎(chǔ)上，將ThetaX與ThetaY得到的結(jié)果進(jìn)行疊加，提出了ED方法，該方法用極大值識(shí)別邊界：

ED=ThetaX+ThetaY.

(5)

本文設(shè)計(jì)了模型試驗(yàn)(圖1).圖1a、b分別是模型的平面圖和三維視圖.測(cè)區(qū)范圍為100 km×100 km，點(diǎn)距為0.5 km，模型參數(shù)見表1.

圖1 模型位置圖 (a) 平面圖； (b) 立體圖.Fig.1 Positions of models (a) Planar view; (b) 3D view.

表1 圖1模型的參數(shù)Table 1 Parameters of the models in Fig.1

在模型試驗(yàn)中，當(dāng)重力梯度張量用于邊界識(shí)別時(shí)，可能會(huì)受到地質(zhì)體走向的影響，因此模型體設(shè)計(jì)成與坐標(biāo)軸成一定角度.此外，模型的設(shè)計(jì)考慮了不同深度地質(zhì)體及正負(fù)異常相間對(duì)邊界識(shí)別的影響.圖2和圖3分別為模型的重力異常圖和理論重力張量圖.

圖2 模型重力異常圖Fig.2 Gravity anomalies of models

圖3 模型理論張量 (a)gxx; (b) gyy; (c) gzz; (d) gxz; (e) gyz; (f) gxy.Fig.3 Synthetic gravity gradient tensors of models

本文將總水平導(dǎo)數(shù)(THDR)方法的邊界識(shí)別效果(圖4)與ED方法的邊界識(shí)別效果(圖5)進(jìn)行對(duì)比，其中THDR方法為重磁數(shù)據(jù)處理中常用的邊界識(shí)別方法，該方法受深度影響.圖4中，模型2較模型1、3深度大，所以模型2的邊界識(shí)別效果較模型1、3差，邊界模糊.

圖4 THDR方法邊界識(shí)別結(jié)果Fig.4 Edge detection resultby THDR

圖5是ED方法邊界識(shí)別效果，Yuan等(2016)等已經(jīng)驗(yàn)證該方法不受正負(fù)異常梯度帶的影響.與THDR方法識(shí)別結(jié)果進(jìn)行比較，其中模型2的邊界識(shí)別效果與模型1的相同，并且比THDR方法識(shí)別出的邊界清晰，說(shuō)明該方法不受深度的影響.圖5中模型3與模型1的深度相同，但是識(shí)別出的邊界清晰程度不同，考慮到只有模型3邊界垂直于x軸或y軸，所以推測(cè)該方法受地質(zhì)體走向的影響.

圖5 ED方法邊界識(shí)別結(jié)果Fig.5 Edge detection result by ED

從ThetaX和ThetaY的等值線圖(圖6)中可以看出，模型體1，2的邊界在ThetaX和ThetaY中分別對(duì)應(yīng)極大值，且大小相近.而模型體3，在ThetaX的圖中只能識(shí)別出垂直于x軸的邊界，在ThetaY的圖中只能識(shí)別出垂直于y軸的邊界，對(duì)ThetaX和ThetaY進(jìn)行疊加后，模型體1,2邊界的數(shù)值絕對(duì)值會(huì)達(dá)到模型體3邊界的2倍左右，因此模型體3的邊界較為模糊.因此，ED方法識(shí)別邊界的結(jié)果中，當(dāng)待識(shí)別區(qū)域內(nèi)存在垂直或平行于坐標(biāo)軸的地質(zhì)體和與坐標(biāo)軸成一定角度的地質(zhì)體時(shí)，與坐標(biāo)軸垂直或平行的地質(zhì)體的邊界較模糊，不利于邊界信息的提取.

圖6 邊界識(shí)別結(jié)果 (a) ThetaX； (b) ThetaY.Fig.6 Edge detection results

2 方法改進(jìn)

為了減小地質(zhì)體走向?qū)D方法的影響，本文提出IED方法，表達(dá)式為

IED(i,j)=max(ThetaX(i,j),ThetaY(i,j)),(6)

式中，i、j分別為點(diǎn)線號(hào)，將ThetaX和ThetaY逐點(diǎn)對(duì)比，取較大的值賦給IED.這樣可以避免與坐標(biāo)軸成一定角度的邊界疊加兩次.但這樣運(yùn)算會(huì)使ThetaX和ThetaY中的較弱的虛假異常凸顯出來(lái).為了減弱這種效應(yīng)，本文根據(jù)ThetaX和ThetaY各自的最小值TXmin和TYmin，定義了閾值T：

T=α·max(TXmin,TYmin).

(7)

由于ThetaX和ThetaY都為負(fù)值，T也為負(fù)值，α取值的范圍為0～1.設(shè)定判定條件：當(dāng)測(cè)點(diǎn)的ThetaX和ThetaY的值都小于T，則IED取二者較小的值；其他情況，IED取二者較大的值，即：

IED(i,j)=

(8)

將其應(yīng)用到模型試驗(yàn)中，令α分別為0.9、0.6、0.3，并與EI方法及常用的斜導(dǎo)數(shù)法、Theta法進(jìn)行對(duì)比(圖7).

由圖7，對(duì)比幾種常用方法及α分別為0.9、0.6、0.3的IED方法邊界識(shí)別結(jié)果，可以看出，斜導(dǎo)數(shù)法和Theta方法都會(huì)在正負(fù)異常之間產(chǎn)生虛假邊界，從而影響邊界識(shí)別效果.IED方法對(duì)三個(gè)模型體的邊界識(shí)別效果都較明顯，一定程度上改善了ED方法中模型體3相比于模型體1、2邊界模糊的現(xiàn)象，而且隨著α值的減小，虛假異常也逐漸減弱，但是當(dāng)α值減小到一定程度，邊界的有用信息也開始有一定的損失，通過(guò)實(shí)驗(yàn)對(duì)比分析得到當(dāng)α值取0.6～0.8時(shí)效果較好.根據(jù)圖1中虛線位置提取斜導(dǎo)數(shù)方法，Theta方法，ED方法及α值為0.6的IED方法邊界識(shí)別結(jié)果剖面圖如圖8所示.

圖7 幾種常用方法及不同α值的IED方法邊界識(shí)別結(jié)果 (a) 斜導(dǎo)數(shù)法； (b) Theta角法； (c) ED法； (d) IED of α=0.9； (e) IED of α=0.6； (f) IED of α=0.3.Fig.7 Edge detection results of several commonly used methods and IED with different values of α (a) Tilt derivative method； (b) Theta angle method； (c) ED method； (d) IED with α=0.9； (e) IED with α=0.6； (f) IED with α=0.3.

圖8a是斜導(dǎo)數(shù)方法邊界識(shí)別結(jié)果的剖面圖，該方法用零值線(圖中虛線)識(shí)別邊界，圖中零值點(diǎn)與地質(zhì)體邊界對(duì)應(yīng)較好，但是在正負(fù)異常之間有多余的零值點(diǎn)d，該點(diǎn)顯示為虛假邊界.圖8b是Theta角法邊界識(shí)別結(jié)果的剖面圖，與圖8a類似，在正負(fù)異常之間也有虛假邊界存在(f點(diǎn)).圖8c是ED方法邊界識(shí)別結(jié)果的剖面圖，從圖中可以看出該方法能夠較好的識(shí)別出地質(zhì)體邊界，且沒有虛假邊界產(chǎn)生，但是左側(cè)異常體邊界處的值是右側(cè)異常體邊界處值的2倍以上，如果地下地質(zhì)構(gòu)造情況較復(fù)雜，有些地質(zhì)體的邊界可能無(wú)法識(shí)別出來(lái).圖8d是IED方法邊界識(shí)別結(jié)果的剖面圖，可以看出該方法能夠較好的識(shí)別出地質(zhì)體邊界，沒有虛假異常產(chǎn)生，且異常體邊界處的值幾乎相等.為進(jìn)一步討論該方法的實(shí)用性，設(shè)計(jì)了加干擾的復(fù)雜模型，模型參數(shù)如表中所示.

圖8 幾種方法邊界識(shí)別結(jié)果剖面圖 (a) 斜導(dǎo)數(shù)法； (b) Theta角法； (c) ED法； (d) IED α=0.6.Fig.8 Edge detection resultsby methods shown by profile (a) Tilt derivative method； (b) Theta angle method； (c) ED method； (d) IED with α=0.6.

圖9為復(fù)雜模型的示意圖，表2為復(fù)雜模型參數(shù).在模型試驗(yàn)中加入異常最大幅值1%的隨機(jī)噪聲.用IED方法進(jìn)行邊界識(shí)別，結(jié)果如圖10所示.

圖9 復(fù)雜模型示意圖Fig.9 Schematic diagram of complex model

從圖10的處理結(jié)果可以看出，對(duì)于加入噪聲的復(fù)雜模型，IED方法也能較好地識(shí)別出邊界.邊界識(shí)別結(jié)果在模型體邊界位置受噪聲影響最小，則該方法在地質(zhì)體邊界處對(duì)干擾的壓制作用最強(qiáng)，這一特點(diǎn)在邊界識(shí)別中是有利的.

圖10 IED方法邊界識(shí)別結(jié)果(α=0.6)Fig.10 Edge detection result of IED(α=0.6)

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

本文實(shí)際數(shù)據(jù)是加拿大自然資源部(Natural Resources Canada)提供的圣喬治灣的(St.George′s Bay)全張量重力梯度測(cè)量數(shù)據(jù)(Miller and Singh，1994)，包含88條測(cè)線，線距為400 m，總飛行長(zhǎng)度為6602 km，測(cè)區(qū)覆蓋面積為3056 km2，飛行高度為80 m，聯(lián)絡(luò)線距為5000 m.該地區(qū)測(cè)得的重力全張量數(shù)據(jù)，如圖11所示圖，其中黑色實(shí)線是已知的地質(zhì)邊界.

圖11 圣喬治灣地區(qū)重力全張量 (a) gxx; (b) gyy; (c) gzz; (d) gxz; (e) gyz; (f) gxy.Fig.11 Fullgravity tensors in St. George′s Bay, Canada

測(cè)區(qū)位于圣勞倫斯在紐芬蘭島西南海岸的海灣，測(cè)區(qū)包含了圣喬治灣的大面積石炭系凹陷，石炭系和上泥盆統(tǒng)的沉積巖在北、東、南三個(gè)方向出露，覆蓋在前寒武紀(jì)基底上.測(cè)區(qū)中分布著不同時(shí)代的地層，可以通過(guò)重力梯度全張量數(shù)據(jù)對(duì)不同密度的地質(zhì)邊界進(jìn)行邊界識(shí)別.圖12為利用ED方法(圖12a、b)和α值為0.7時(shí)的IED方法(圖12c、d)進(jìn)行邊界識(shí)別的圖像.

對(duì)比圖12a、c，可以看出兩種方法都能識(shí)別出該地區(qū)大的構(gòu)造邊界.ED法和IED法在對(duì)斷裂F1、F2、F4、F6、F7的識(shí)別效果是相近的，因?yàn)檫@5條斷裂與水平方向存在一定的夾角，而對(duì)于斷裂F3與F5，對(duì)比圖中紅色虛線圈閉框(B和C)標(biāo)示出來(lái)的區(qū)域(即F5斷裂)，會(huì)發(fā)現(xiàn)ED方法一些邊界由于受走向的影響在圖中只有較弱的顯示，邊界不連續(xù)，點(diǎn)狀分布，而IED方法能清楚地分辨出邊界的位置；對(duì)比紅色虛線框(A)標(biāo)示出來(lái)的區(qū)域(即斷裂F3位置)，可以看出一些邊界通過(guò)ED方法甚至無(wú)法識(shí)別出來(lái)，但通過(guò)IED方法卻能得到清楚連續(xù)的邊界，可以對(duì)F3斷裂進(jìn)行更好的解釋.通過(guò)對(duì)實(shí)際數(shù)據(jù)的處理發(fā)現(xiàn)，與ED方法對(duì)比，IED方法在實(shí)際數(shù)據(jù)處理中邊界識(shí)別更加清晰，細(xì)節(jié)處更加明顯，處理結(jié)果更加可靠.

圖12 圣喬治灣重力梯度全張量數(shù)據(jù)處邊界識(shí)別效果 (a)和(b) ED方法； (c)和(d) IED方法.Fig.12 Edge detection results of full tensor gravity gradient data in the St. George′s Bay (a)and(b) ED method； (c)and(d) IED method.

4 結(jié)論

為了解決地質(zhì)構(gòu)造走向?qū)吔缱R(shí)別的影響，本文提出了一種改進(jìn)的ED方法(IED)，通過(guò)理論研究、模型試驗(yàn)及實(shí)際數(shù)據(jù)的處理驗(yàn)證了該方法的可靠性.得出以下結(jié)論：

(1)ED方法在邊界識(shí)別中，會(huì)受到地質(zhì)體走向的影響，對(duì)于平行或者垂直于坐標(biāo)軸方向的邊界識(shí)別效果較差，導(dǎo)致這種現(xiàn)象的原因是原方法中疊加項(xiàng)對(duì)于識(shí)別這種走向的地質(zhì)體較該坐標(biāo)軸成一定角度的地質(zhì)體的能力弱.

(2)改進(jìn)后IED方法，摒棄了原方法中的疊加項(xiàng)，選擇合適的閾值來(lái)減少虛假異常，從而使該方法不受地質(zhì)體走向的影響，進(jìn)一步提高邊界識(shí)別效果.

(3)在實(shí)際數(shù)據(jù)處理中IED方法比ED方法識(shí)別出的邊界細(xì)節(jié)更多，同時(shí)邊界也更連續(xù)，為地質(zhì)解釋提供了更詳細(xì)的信息.

致謝感謝加拿大自然資源部(Natural Resources Canada)提供圣喬治灣(St.George′s Bay)全張量重力梯度測(cè)量數(shù)據(jù)，并允許發(fā)表數(shù)據(jù)處理結(jié)果.