位場數(shù)據(jù)解釋的Theta-Depth法

陳國強， 馬國慶

1 吉林大學地球科學學院， 長春　130061　2 吉林大學地球探測科學與技術學院， 長春　130021

?

位場數(shù)據(jù)解釋的Theta-Depth法

陳國強1， 馬國慶2*

1 吉林大學地球科學學院， 長春1300612 吉林大學地球探測科學與技術學院， 長春130021

摘要Theta圖是利用位場(重磁)數(shù)據(jù)識別邊界的常用方法，其表達式為重磁異常水平變化與垂直變化的比值函數(shù).該方法計算淺源地質體邊界的效果較好，而由于深源位場數(shù)據(jù)在換算過程中會產(chǎn)生趨同效應，在深源地質體識別應用中計算結果不準確，為此，本文提出Theta-Depth法并進行地質體埋深的計算.首先給出直接利用Theta圖像進行場源體深度估算的方法，然后推導出基于Theta導數(shù)的線性方程來自動估算場源位置參數(shù)，本文方法可有效地利用Theta圖像的特征為約束條件來提高反演結果的精度.理論模型試驗證明本文提出的Theta-Depth法能有效地計算出場源體位置和深度.將該方法應用于滿都拉地區(qū)實測磁數(shù)據(jù)的解釋，幫助圈定了礦脈的分布.

關鍵詞Theta； 位場； 深度； 線性方程

1引言

邊界識別是位場數(shù)據(jù)解釋的常規(guī)任務之一，現(xiàn)今有很多方法來完成這一任務(Evjen,1936; Cordell and Grauch, 1982; Blakely, 1995).重磁異常水平導數(shù)的最大值、垂直導數(shù)的零值與地質體的邊界相對應(Roest et al., 1992; Miller and Singh, 1994; Blakely, 1995)，因而邊界識別方法大多是水平與垂直導數(shù)的函數(shù).前期邊界識別方法僅能給出較淺地質體的邊界，而不能給出深部地質體的邊界.為了改善這一缺點，人們開始致力于均衡邊界識別濾波器的研究(Miller and Singh,1994; Verduzco et al., 2004; Wijns et al.,2005; Cooper and Cowan, 2006, 2008; Cooper, 2009; 馬國慶等， 2012； Ma, 2013)，該類方法能同時識別出不同深度地質體的邊界.傾斜角法(Miller and Singh,1994)是第一個提出的均衡邊界識別方法，但該方法所得到的邊界并不十分清晰.Theta圖(Wijns et al.,2005)是一種常用的均衡邊界識別濾波器，其能清晰地給出場源體的水平位置.Cooper和Cowan (2006) 將該方法與傾斜角法進行比較，并將其應用于實際數(shù)據(jù)的邊界識別任務.馬國慶等(2012)將該方法應用于四川盆地線性構造的劃分，斷裂結果與實際地質資料吻合較好.Ma和Li (2012) 總結Theta圖法的應用特征，并用于中國內蒙古航磁數(shù)據(jù)的構造劃分.

地質體位置反演是重磁數(shù)據(jù)解釋的主要任務，隨著勘探數(shù)據(jù)量的增加，人們更加傾向于選擇自動解釋方法(Debeglia and Corpel, 1997; Thurston and Smith, 1997; Hsu et al., 1998; Smith et al., 1998; Kirkham, 2001; Salem et al., 2004; Salem and Smith, 2005)來進行計算，該類方法在給定異常后可直接計算出地質體的位置參數(shù)，主要包括主要包括歐拉反褶積法(Thompson, 1982; Reid et al., 1990)、局部波數(shù)法(Salem et al., 2008)和解析信號法(Kirkham, 2001)，通過原始異常構建相應函數(shù)于位置參數(shù)的對應方程，進而完成計算，在實際數(shù)據(jù)解釋中應用較為廣泛.歐拉反褶積法是由Thompson在1982年提出應用于磁法數(shù)據(jù)的解釋，隨后被應用于重力數(shù)據(jù)的解釋，解析信號和局部波數(shù)法在進行磁異常解釋時具有不受傾斜磁化干擾的優(yōu)點.眾所周之，方程求解會造成結果的發(fā)散效應，因而以上方法在求解后需設定結果篩選準則來獲得更加準確的結果，由于以上方法的函數(shù)本身與地質體水平位置之間不具有直接的對應關系，往往需要借助其他函數(shù)來進行，從而造成結果之間缺乏對應性.

本文提出Theta-Depth法計算異常體的深度信息，首先推導出一種直接利用Theta進行深度計算的方法，通過理論推導證明Theta反正弦函數(shù)兩條45°等值線之間距離的一半等于地質體的埋深，利用此關系進行深度的估算.此外，本文還推導出一種基于Theta的自動解釋方法，其利用Theta導數(shù)組成的線性方程來完成場源體位置和深度的計算，并根據(jù)Theta函數(shù)與地質體位置之間的關系設定相應的篩選準則.

2Theta-Depth法

Theta是由Wijns于2005年提出的一種進行邊界識別的方法，其表達式為

(1)

Theta圖的最大值對應著地質體的邊界，該方法具有能夠同時識別出不同深度地質體邊界的優(yōu)勢，在實際數(shù)據(jù)解釋中應用較為廣泛.

(2)

(3)其中，K為磁化率差，F(xiàn)為當?shù)氐牡卮艌鰪姸?，c=1-cos2isin2A，A為h軸與磁北方向的夾角，i為磁傾角，tanI=tani/cosA，d為地質體的傾斜角度.

在垂直磁化情況下式(2)和式(3)可以改寫為

(4)

(5)

將式(4)和式(5)帶入到式(1)，整理后可以得到

(6)

通過式(6)可以看出，當h=0時式(6)的反正弦函數(shù)取得最大值90°，因此可通過式(6)反正弦的最大值判斷出地質體的水平位置；當h=z時式(6)反正弦函數(shù)為45°，隨著h的增大式(6)的反正弦函數(shù)繼續(xù)減小直至接近最小值0.通過以上的分析可知，利用Theta反正弦函數(shù)兩條45°等值線之間距離的1/2可估算出地質體的埋藏深度.需要指出的是，利用45°等值線之間距離計算深度只是其中一種計算深度的方法，任何的角度都可以建立起水平距離與埋深的關系.

本文還推導出一種利用Theta導數(shù)組成的線性方程進行位場異常反演的方法.常規(guī)歐拉反褶積法的表達式為

-N(f-B)，

(7)

其中，x、y和z為觀測點坐標，x0、y0和z0為待求的場源體的空間位置坐標，B為背景場值，N為表征場源體類型的構造指數(shù).計算式(7)在x和y方向上的導數(shù)，由于背景異常變化較為平緩，其導數(shù)較小可忽略，因此式(7)導數(shù)的表達式為

(8)

(9)

(10)

Huang等(1995) 證明解析信號的歐拉反褶積方程的表達式為

(11)

(12)

計算式(1)在x,y和z方向上的導數(shù)為

(13)

(14)

(15)

將式(13)、(14)和(15)帶入到式(12)可以得到

(16)

根據(jù)式(16)可知，利用Theta的導數(shù)所組成的線性方程可估算出場源體的位置坐標x0、y0和z0.方程求解往往會造成結果的發(fā)散，因此利用Theta函數(shù)極大值與場源體之間的對應關系設定如下的結果篩選準則.

(1) 距離Theta極大值超過窗口長度的解去掉.

(2) 當求解結果與窗口中心點之間的距離大于窗口半徑時該解無效，因為當窗口在場源附近時解的準確性較高.

(3) 經(jīng)過以上的篩選后，解應該已經(jīng)比較集中，對于一定范圍內解的總數(shù)小于給定值的異常進行濾除，因為當解的個數(shù)小于一定數(shù)量時不認為是有效異常.

3模型試驗

為驗證上述算法的精確性與實用性，本文進行了模型試驗,主要分為原始異常和含噪異常來進行試驗，并分別針對重磁異常進行了試驗來驗證方法的穩(wěn)定性和適應性.

由于Theta法的結果會受到傾斜磁化的干擾，因此利用理論模型產(chǎn)生的磁異常和重力異常來試驗Theta-Depth方法的應用效果.圖1a垂直磁化條件下埋深分別為5 m和7 m正方體所產(chǎn)生的磁異常.圖1b為磁異常Theta計算結果的反正弦，其中45°等值線采用虛線表示，并用灰色進行填充之間的區(qū)域.通過公式(6)的推導證明地質體的埋深等于Theta反正弦結果的兩條45°等值線間距離的一般，其區(qū)域極大值標識地質體的水平范圍.根據(jù)圖1b利用兩條45°等值線(虛線)的垂直距離可計算出地質體的埋深分別為4.86 m和6.72 m，與地質體理論埋深差距較小，說明該方法具有良好的應用效果.

利用公式(16)計算地質體的埋深，其窗口大小為7×7，并采用本文給出的結果篩選準則來獲得更加準確的結果.圖1c為利用Theta導數(shù)組成線性方程所計算得到的異常體的水平位置分布，其結果與第一種方法所得到的異常體的水平位置相一致，且與地質體的理論位置吻合較好.圖1d為Theta導數(shù)線性方程反演結果的三維顯示，地質體的埋深分別為4.91 m和6.87 m，與理論埋深相接近.通過該試驗證明Theta-Depth法能準確地計算出地下異常體的深度.

圖1　(a) 垂直磁化情況下正方體所產(chǎn)生的磁力異常， (b) 異常Theta的反正弦計算結果，(c) 利用式(16)反演得到的場源體的水平位置， (d) 反演結果的三維顯示Fig.1　(a) Magnetic anomalies of a cub in the vertical magnetization， (b) Arcsine calculation of anomalies using Theta method，(c) Horizontal position of source by inversion， (d) Three-dimensional display of inversion results

下面試驗一下Theta-Depth方法在含噪重力異常解釋中的應用效果.圖2a為埋深分別為7 m和10 m正方體所產(chǎn)生的重力異常，并加入了均值為0，方差為1 mGal的隨機噪聲.圖2b為重力異常的Theta結果，其最大值能準確地標識出地質體的水平位置.圖2c為Theta的反正弦計算結果，并利用綠色填充45°等值線間的區(qū)域，根據(jù)45°等值線之間的距離可以得到異常體的埋深分別為6.81 m和9.66 m.

圖2d為利用公式(16)計算得到的結果，其窗口大小為7×7.從結果中可以看出由于噪聲的加入產(chǎn)生了較多的干擾異常.采用上述篩選準則來獲得更加準確的結果，圖2e為方程反演結果經(jīng)過篩選后得到的結果篩，可以看出該篩選準則很好地去除了隨機噪聲的干擾.圖2f為篩選后反演結果的三維顯示，異常體的深度為6.87 m和9.73 m，與地質體的真實埋深相接近.

圖2　(a) 正方體產(chǎn)生的重力異常并加入隨機噪聲后異常; (b) 圖2a所示異常的Theta結果； (c)Theta的反正弦結果; (d) Theta-Depth法反演得到的異常體的水平位置; (e) 篩選后的反演結果; (f) 反演結果的三維顯示Fig.2　(a) Gravity anomalies of the cube with added random noise; (b) Theta results of anomalies in (a);(c)Arcsine results by Theta method; (d) Horizontal position of abnormal body obtained by inversion of the Theta-Depth method; (e) Inversion results after screening; (f) Three-dimensional display of inversion results

4實際數(shù)據(jù)應用效果

將該方法應用于滿都拉地區(qū)實測磁異常數(shù)據(jù)的解釋中，磁數(shù)據(jù)的采集點距為40 m，該地區(qū)主要是為了尋找地下鐵礦場的位置，周圍礦產(chǎn)的開采深度在100 m左右.圖3a為化極后磁異常，可以看出該地區(qū)構造呈東南向展布.圖3b為磁異常的Theta計算結果，其最大值給出了了地下礦體及地層之間的界線.圖3c為Theta的反正弦結果，并采用綠色對45°等值線間進行填充，根據(jù)距離可估算礦脈的深度約為140 m.圖3d為利用式(16)反演得到的異常體的分布，并對結果進行篩選.從結果中可以看出異常體的深度大多集中在80～200 m之間，與另一種Theta-Depth方法計算得到的結果相一致，說明本文提出的Theta-Depth法具有較好的實際應用效果，且反演深度與以往礦產(chǎn)深度相接近.

圖3　(a)滿都拉地區(qū)化極后磁異常; (b) 磁異常的Theta結果; (c) Theta反正弦結果; (d) 式(16)反演得到的結果Fig.3　(a) Magnetic anomalies after reduction to the pole in Mandula area; (b) Theta results of magnetic anomalies;(c) Arcsine results by Theta method; (d) Inversion results by Eq.16

5結論

地質體的埋深確定是重磁數(shù)據(jù)處理與解釋的主要任務，在混成性地質空間分析中，確定地質體的水平邊界與垂向邊界對于地質體數(shù)字特征研究十分重要，尤其是有關含礦地質體埋深推斷事關資源預測精度水平，本文提出的Theta-Depth法較好的完成了上述任務.Theta-Depth法首次將Theta函數(shù)進行擴展來完成異常體深度的計算，通過推導證明根據(jù)Theta反正弦函數(shù)可直接給出地質體的埋深，此外還推導出基于Theta導數(shù)的線性方法來估算地質體的水平位置和深度，且該方法在進行反演前不需要給定任何的先驗信息.通過理論重力和磁數(shù)據(jù)證明本文提出的Theta-Depth方法能有效地完成重磁異常的解釋工作，且反演結果與理論結果相一致.利用噪聲重力異常試驗可以看出本文方法依舊可以獲得較為準確的結果，且對結果進行篩選后所得到的結果與理論值相一致.將本文方法應用于滿都拉地區(qū)實際磁異常的解釋，獲得了巖脈的大致深度.

ReferencesBlakely R J. 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press. Cooper G R J. 2009. Balancing images of potential-field data.Geophysics, 74(3): L17-L20.

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.

Cooper G R J, Cowan D R. 2008. Edge enhancement of potential-field data using normalized statistics.Geophysics, 73(3): H1-H4.

Cordell L, Grauch V J S. 1982. Mapping basement magnetization zones from aeromagnetic data in the San Juan basin, New Mexico.∥SEG Technical Program Expanded Abstracts. SEG, 246-247.

Debeglia N, Corpel J. 1997. Automatic 3-D interpretation of potential field data using analytic signal derivatives.Geophysics, 62(1): 87-96. Evjen H M. 1936. The place of the vertical gradient in gravitational interpretations.Geophysics, 1(1): 127-136.

Hsu S K, Coppens D, Shyu C T. 1998. Depth to magnetic source using the generalized analytic signal.Geophysics, 63(6): 1947-1957.

Huang D, Gubbins D, Clark R A, et al. 1995, Combined Study of Euler′s Homogeneity Equation for Gravity and Magnetic field . 57th EAGE conference. Glasgow UK Extended Abstracts, 144.

Kirkham K. 2001. Investigations of a high-resolution aeromagnetic survey over the southeastern portion of the Illinois Basin [Master′s thesis]. Carbondale: Southern Illinois University of Carbondale. Ma G Q. 2013. Edge detection of potential field data using improved local phase filter.ExplorationGeophysics, 44(1): 36-41.

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, Li L L. 2012. Edge detection in potential fields with the normalized total horizontal derivative.Computers&Geosciences, 41: 83-87. 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:507-517.Reid A B, Allsop J M, Granser H, et al. 1990. Magnetic interpretation in three dimensions using Euler deconvolution.Geophysics, 55(1): 80-91.

Roest W R, Verhoef J, Pilkington M. 1992. Magnetic interpretation using the 3-D analytic signal.Geophysics, 57(1): 116-125. Salem A, Ravat D, Mushayandebvu M F, et al. 2004. Linearized least-squares method for interpretation of potential-field data from sources of simple geometry.Geophysics, 69(3): 783-788.

Salem A, Smith R S. 2005. Depth and structural index from normalized local wavenumber of 2D magnetic anomalies.GeophysicalProspecting, 53(1): 83-89.

Salem A, Williams S, Fairhead D, et al. 2008. Interpretation of magnetic data using tilt-angle derivatives.Geophysics, 73(2): L1-L10.

Smith R S, Thurston J B, Dai T F, et al. 1998. iSPITM—The improved source parameter imaging method.GeophysicalProspecting, 46(2): 141-151. Thompson D T. 1982. EULDPH: A new technique for making computer-assisted depth estimates from magnetic data.Geophysics, 47(1): 31-37. Thurston J B, Smith R S. 1997. Automatic conversion of magnetic data to depth, dip, and susceptibility contrast using the SPI (TM) method.Geophysics, 62(3): 807-813.

Verduzco B, Fairhead J D, Green C M, et al. 2004. New insights into magnetic derivatives for structural mapping.TheLeadingEdge, 23(2): 116-119.

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

附中文參考文獻

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

(本文編輯張正峰)

基金項目國家重點基礎研究發(fā)展計劃(973計劃)(2015CB453000)、國家自然科學基金(41404089)和中國地質調查局地質礦產(chǎn)調查評價專項項目(GZH003-07-03)聯(lián)合資助.

作者簡介陳國強，男，1985年生，博士研究生，主要從事礦產(chǎn)資源定量與定位預測研究. E-mail：chenguoqiang0506@163.com *通訊作者馬國慶，男，1984年生，講師，主要從事重磁數(shù)據(jù)解釋方法研究.E-mail：maguoqing@jlu.edu.cn

doi:10.6038/cjg20160625 中圖分類號P631

收稿日期2015-06-02，2016-03-16收修定稿

Theta-Depth method for the interpretation of potential field data

CHEN Guo-Qiang1, MA Guo-Qing2*

1CollegeofEarthScience,JilinUniversity,Changchun130061,China2CollegeofGeoexplorationScienceandTechnology,JilinUniversity,Changchun130021,China

AbstractThe Theta map is a commonly used tool to detect edges of subsurface geological bodies using potential field (gravity and magnetic) data, which is the ratio of function of horizontal derivative to vertical derivative of gravity or magnetic anomalies. Although the effect of calculation of shallow geological body edges by this method is good, it is not successful in probing the deep geological body edges process due to the convergence of the conversion process that leads inaccurate calculation. To solve this problem, we proposed the Theta-Depth method to compute the depth of the buried geological body. Frist, we directly used the Theta map to estimate the depth of the source. Then we derived a linear equation based on the derivatives of Theta to compute the location and depth of the source. The presented method can use the Theta map as a constraint to obtain more accurate results effectively. The tests show that this Theta-Depth method can effectively calculate the location and depth of the source. We applied the Theta-Depth method to interpretation of measured magnetic data in the Mandula area, which helped delineate the distribution of ore veins.

KeywordsTheta; Potential field data; Depth; Linear equation

陳國強， 馬國慶. 2016. 位場數(shù)據(jù)解釋的Theta-Depth法. 地球物理學報，59(6):2225-2231,doi:10.6038/cjg20160625.

Chen G Q, Ma G Q. 2016. Theta-Depth method for the interpretation of potential field data.ChineseJ.Geophys. (in Chinese),59(6):2225-2231,doi:10.6038/cjg20160625.