陳國(guó)強(qiáng)， 馬國(guó)慶

1 吉林大學(xué)地球科學(xué)學(xué)院， 長(zhǎng)春　130061　2 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春　130021

?

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

陳國(guó)強(qiáng)1， 馬國(guó)慶2*

1 吉林大學(xué)地球科學(xué)學(xué)院， 長(zhǎng)春1300612 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春130021

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

關(guān)鍵詞Theta； 位場(chǎng)； 深度； 線性方程

1引言

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

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

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

2Theta-Depth法

Theta是由Wijns于2005年提出的一種進(jìn)行邊界識(shí)別的方法，其表達(dá)式為

(1)

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

(2)

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

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

(4)

(5)

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

(6)

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

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

-N(f-B)，

(7)

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

(8)

(9)

(10)

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

(11)

(12)

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

(13)

(14)

(15)

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

(16)

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

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

(2) 當(dāng)求解結(jié)果與窗口中心點(diǎn)之間的距離大于窗口半徑時(shí)該解無效，因?yàn)楫?dāng)窗口在場(chǎng)源附近時(shí)解的準(zhǔn)確性較高.

(3) 經(jīng)過以上的篩選后，解應(yīng)該已經(jīng)比較集中，對(duì)于一定范圍內(nèi)解的總數(shù)小于給定值的異常進(jìn)行濾除，因?yàn)楫?dāng)解的個(gè)數(shù)小于一定數(shù)量時(shí)不認(rèn)為是有效異常.

3模型試驗(yàn)

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

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

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

圖1　(a) 垂直磁化情況下正方體所產(chǎn)生的磁力異常， (b) 異常Theta的反正弦計(jì)算結(jié)果，(c) 利用式(16)反演得到的場(chǎng)源體的水平位置， (d) 反演結(jié)果的三維顯示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

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

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

圖2　(a) 正方體產(chǎn)生的重力異常并加入隨機(jī)噪聲后異常; (b) 圖2a所示異常的Theta結(jié)果； (c)Theta的反正弦結(jié)果; (d) Theta-Depth法反演得到的異常體的水平位置; (e) 篩選后的反演結(jié)果; (f) 反演結(jié)果的三維顯示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í)際數(shù)據(jù)應(yīng)用效果

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

圖3　(a)滿都拉地區(qū)化極后磁異常; (b) 磁異常的Theta結(jié)果; (c) Theta反正弦結(jié)果; (d) 式(16)反演得到的結(jié)果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結(jié)論

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

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.

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

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

(本文編輯張正峰)

基金項(xiàng)目國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(973計(jì)劃)(2015CB453000)、國(guó)家自然科學(xué)基金(41404089)和中國(guó)地質(zhì)調(diào)查局地質(zhì)礦產(chǎn)調(diào)查評(píng)價(jià)專項(xiàng)項(xiàng)目(GZH003-07-03)聯(lián)合資助.

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

doi:10.6038/cjg20160625 中圖分類號(hào)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

陳國(guó)強(qiáng)， 馬國(guó)慶. 2016. 位場(chǎng)數(shù)據(jù)解釋的Theta-Depth法. 地球物理學(xué)報(bào)，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.