一種改進的重力梯度全張量數據的邊界識別方法

杜威， 吳賀宇， 陳祥忠， 劉俊成， 張陽陽， 惠夢琳

1 中國科學院地球化學研究所礦床地球化學重點研究實驗室， 貴陽 550081 2 吉林大學地球探測科學與技術學院， 長春 130026 3 北京桔燈地球物理勘探股份有限公司， 北京 102200 4 安徽省勘查技術院， 合肥 230031

0 引言

隨著重力梯度測量技術的發(fā)展，越來越多的重力梯度全張量數據被應用到數據處理中.重力梯度全張量數據是重力位的二階偏導數，由9個張量數據組成，較重力數據包含有更高頻率的信息.對重力梯度全張量數據進行處理，可以得到更全面，更準確的地質信息.

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

1 ED方法及存在問題

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

(1)

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

gxx+gyy+gzz=0，

(2)

T為對稱矩陣，它只包括5個獨立張量，gxy、gxz、gyz、gxx和gzz.方向解析信號的振幅為

(3)

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

(4)

并在此基礎上，將ThetaX與ThetaY得到的結果進行疊加，提出了ED方法，該方法用極大值識別邊界：

ED=ThetaX+ThetaY.

(5)

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

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

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

在模型試驗中，當重力梯度張量用于邊界識別時，可能會受到地質體走向的影響，因此模型體設計成與坐標軸成一定角度.此外，模型的設計考慮了不同深度地質體及正負異常相間對邊界識別的影響.圖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

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

圖4 THDR方法邊界識別結果Fig.4 Edge detection resultby THDR

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

圖5 ED方法邊界識別結果Fig.5 Edge detection result by ED

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

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

2 方法改進

為了減小地質體走向對ED方法的影響，本文提出IED方法，表達式為

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

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

T=α·max(TXmin,TYmin).

(7)

由于ThetaX和ThetaY都為負值，T也為負值，α取值的范圍為0～1.設定判定條件：當測點的ThetaX和ThetaY的值都小于T，則IED取二者較小的值；其他情況，IED取二者較大的值，即：

IED(i,j)=

(8)

將其應用到模型試驗中，令α分別為0.9、0.6、0.3，并與EI方法及常用的斜導數法、Theta法進行對比(圖7).

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

圖7 幾種常用方法及不同α值的IED方法邊界識別結果 (a) 斜導數法； (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，該點顯示為虛假邊界.圖8b是Theta角法邊界識別結果的剖面圖，與圖8a類似，在正負異常之間也有虛假邊界存在(f點).圖8c是ED方法邊界識別結果的剖面圖，從圖中可以看出該方法能夠較好的識別出地質體邊界，且沒有虛假邊界產生，但是左側異常體邊界處的值是右側異常體邊界處值的2倍以上，如果地下地質構造情況較復雜，有些地質體的邊界可能無法識別出來.圖8d是IED方法邊界識別結果的剖面圖，可以看出該方法能夠較好的識別出地質體邊界，沒有虛假異常產生，且異常體邊界處的值幾乎相等.為進一步討論該方法的實用性，設計了加干擾的復雜模型，模型參數如表中所示.

圖8 幾種方法邊界識別結果剖面圖 (a) 斜導數法； (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為復雜模型的示意圖，表2為復雜模型參數.在模型試驗中加入異常最大幅值1%的隨機噪聲.用IED方法進行邊界識別，結果如圖10所示.

圖9 復雜模型示意圖Fig.9 Schematic diagram of complex model

從圖10的處理結果可以看出，對于加入噪聲的復雜模型，IED方法也能較好地識別出邊界.邊界識別結果在模型體邊界位置受噪聲影響最小，則該方法在地質體邊界處對干擾的壓制作用最強，這一特點在邊界識別中是有利的.

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

3 實際數據處理

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

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

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

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

圖12 圣喬治灣重力梯度全張量數據處邊界識別效果 (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 結論

為了解決地質構造走向對邊界識別的影響，本文提出了一種改進的ED方法(IED)，通過理論研究、模型試驗及實際數據的處理驗證了該方法的可靠性.得出以下結論：

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

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

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

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