基于線性回歸的重力數(shù)據(jù)三維RMS成像

吉日嘎拉圖，魯光銀

(中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 長(zhǎng)沙 410083)

0 引言

重力勘探是最早應(yīng)用在油氣勘探領(lǐng)域的地球物理方法，雖然現(xiàn)在逐漸被地震勘探所取代[1]，但是在深部勘探等方面的優(yōu)勢(shì)，是其他地球物理勘探方法所無(wú)法取代的。而且重力勘探方法能在各種平臺(tái)上獲得大量觀測(cè)數(shù)據(jù)，從而在大地構(gòu)造分區(qū)、礦產(chǎn)資源勘查以及基礎(chǔ)地質(zhì)研究提供了重要的地球物理依據(jù)[2]，同時(shí)也作為其他地球物理勘探方法反演算法中的約束發(fā)揮著重要作用[1]。在任何地球物理勘探方法的應(yīng)用過(guò)程中，必須對(duì)實(shí)際測(cè)量異常進(jìn)行反演解釋，重力勘探也不例外，重力數(shù)據(jù)解釋的基本目標(biāo)就是要得到一個(gè)清晰的異常源[3]，而且這也是應(yīng)用重力勘探的重要環(huán)節(jié)，常規(guī)重力異常反演方法是基于反演理論的、在最小二乘意義下，使目標(biāo)函數(shù)達(dá)到最小的線性或者非線性反演[4]。根據(jù)重力異常，利用一些統(tǒng)計(jì)分析的方法(如地層視密度的解釋方法)，來(lái)求取異常源的位置和物性參數(shù)，大體上也是屬于反演問(wèn)題[5]。但前一種反演方法必須計(jì)算非常龐大的數(shù)據(jù)量，盡管姚常利[6]等學(xué)者都提出了相關(guān)的解決方法，但是仍然要耗費(fèi)大量的計(jì)算時(shí)間。

基于概率統(tǒng)計(jì)的數(shù)據(jù)解釋及處理方法，最早是由Spector和Grant提出的[7]，經(jīng)過(guò)不斷地發(fā)展，在1979年利用統(tǒng)計(jì)學(xué)的滑窗法處理磁場(chǎng)剖面數(shù)據(jù)的程序用來(lái)進(jìn)行場(chǎng)源深度估計(jì)被發(fā)表[8-9]；1981年Chandler[10]利用滑動(dòng)窗口的方式進(jìn)行重磁數(shù)據(jù)的線性回歸計(jì)算，這也是日后被大家所熟知的重磁異常對(duì)應(yīng)分析[11]；概率成像，最初是由自然電位法的解釋方法發(fā)展而來(lái)，經(jīng)Patella和Mauriello[12-14]兩位學(xué)者不斷發(fā)展，已將該方法應(yīng)用位場(chǎng)勘探等領(lǐng)域，而且取得了令人滿意的結(jié)果；Domenico Patella[15]利用概率密度成像方法對(duì)維蘇威火山地區(qū)進(jìn)行概率密度成像研究，而且提出了聯(lián)合概率成像；郭良輝，孟曉紅[4, 16-19]提出基于相同匹配濾波技術(shù)的三維相關(guān)成像[20]，并且已將該方法推廣到重力梯度數(shù)據(jù)應(yīng)用領(lǐng)域。由于梯度數(shù)據(jù)本身所具有高分辨率的性質(zhì)，使得重力梯度數(shù)據(jù)的成像結(jié)果要明顯好于直接利用重力數(shù)據(jù)進(jìn)行成像的結(jié)果。對(duì)于以上這些基于統(tǒng)計(jì)學(xué)的反演方法，普遍都具有的一個(gè)優(yōu)勢(shì)就是即使沒(méi)有任何先驗(yàn)信息，都可以得到穩(wěn)定的解釋結(jié)果[20]。

作者正是受到重磁異常對(duì)應(yīng)分析和重磁概率成像(相關(guān)成像)算法的啟發(fā)，提出了一種基于線性回歸的重力數(shù)據(jù)均方根殘差(RMS)成像方法，對(duì)地下重力異常源的分布情況進(jìn)行研究和探索方法。該方法對(duì)于重力異常源的位置較為敏感，所以可以用來(lái)進(jìn)行成像分析。

1 方法原理

地下異常源在，水平測(cè)區(qū)的所產(chǎn)生的重力異常用公式(1)表示：

(1)

其中 G為萬(wàn)有引力常量 G=6.67×10-11m3/(kg·s2)。

借助相關(guān)成像[4]的原理，將式(1)進(jìn)行推廣，設(shè)置笛卡爾坐標(biāo)的(x,y)平面在基準(zhǔn)面上，z取垂直向下為正，假設(shè)測(cè)區(qū)地下任意點(diǎn)源坐標(biāo)為(xi,yi,zi)、體積為vi的第i個(gè)點(diǎn)質(zhì)量的剩余密度為σi，則它在測(cè)區(qū)上任意位置(x,y,z)處的重力異常為：

(2)

根據(jù)式(2)，我們將重力異常抽象概括為地下異常體物性參數(shù)與幾何構(gòu)造函數(shù)[6]的乘積：

Δgi=GσvQ(xi,yi,zi)

(3)

Q(xi,yi,zi)表示地下異常體的幾何構(gòu)造函數(shù)：

yi)2+(z-zi)2]3/2

(4)

利用一元線性回歸的原理，將重力異常表示為：

Δg(x,y,z)=A·Q(xi,yi,zi)+b

(5)

所以地上所測(cè)重力異常與地下第i個(gè)點(diǎn)所對(duì)應(yīng)的異常源幾何構(gòu)造函數(shù)存在著線性的對(duì)應(yīng)關(guān)系，在進(jìn)行線性回歸分析之前，必須進(jìn)行相關(guān)性檢驗(yàn)，驗(yàn)證是否是存在線性相關(guān)性。

對(duì)任意一組觀測(cè)值與幾何構(gòu)造函數(shù)，不管是否相關(guān)都可以建立一個(gè)線性回歸方程，這就產(chǎn)生了一個(gè)問(wèn)題，這樣得到的線性回歸方程是否真的具有線性關(guān)系，而且所求的線性回歸方程是不是能夠反映它們之間的實(shí)際關(guān)系，因此進(jìn)行線性之前，必須進(jìn)行線性相關(guān)關(guān)系檢驗(yàn)[6]。

根據(jù)統(tǒng)計(jì)學(xué)原理[21]，相關(guān)系數(shù)被定義為式(6)：

(6)

其中 cov(Δg(x,y,z),Q(xi,yi,zi))為二者的協(xié)方差；D(Δg(x,y,z))和D(Q(xi,yi,zi))表示其相應(yīng)的方差。

當(dāng)滿足線性相關(guān)性之后，利用最小二乘原理，求解線性回歸方程，得到等式(7)所表示的線性回歸方程斜率和截距公式：

(7)

根據(jù)均方根殘差法(RMS)原理，對(duì)其進(jìn)行修正，用相關(guān)系數(shù)R當(dāng)作ΔGcalc的權(quán)，相關(guān)程度越高，擬合程度越接近，才能越接近真實(shí)觀測(cè)值，因此得到了等式(8)加入相關(guān)系數(shù)R修正后的RMSm，因?yàn)橄嚓P(guān)系數(shù)-1≤R≤1，故對(duì)于式(8)，只需要尋找局部極值(包括極大值和極小值)就可以準(zhǔn)確找到異常源的所處位置。

RMSmodified=

(8)

圖1 利用線性回歸進(jìn)行場(chǎng)源分析的流程圖Fig．1 The flow chart of processing 3-D RMSm imaging

2 基于線性回歸的重力數(shù)據(jù)三維RMS成像

如圖1所示的流程已經(jīng)完整的闡述了利用線性回歸方法進(jìn)行重力異常場(chǎng)源分析的工作流程。需要注意的是，此處的重力異常數(shù)據(jù)是經(jīng)過(guò)布格校正和網(wǎng)格化的。在經(jīng)過(guò)位場(chǎng)延拓和相關(guān)的濾波算法，可以在一定程度上提高沿深度方向上的分辨率[4]。

在計(jì)算時(shí)，首先要對(duì)地下空間進(jìn)行剖分計(jì)算，為了方便計(jì)算，一般選擇利用等間距的網(wǎng)格剖分方法，將地下空間剖分成均勻網(wǎng)格。此時(shí)假設(shè)地下空間中必定存在某一些點(diǎn)是滿足等式(5)的，所以利用公式(6)、公式(7)、公式(8)由淺到深的逐層進(jìn)行計(jì)算，最后利用RMSm的數(shù)據(jù)進(jìn)行成像。

3 算例分析

為了研究該方法的具體表現(xiàn)，我們利用兩個(gè)密度為0.5 g/cm3的異常體，設(shè)置在地面以下100 m處的位置如 圖 2所示。圖3(a)為測(cè)網(wǎng)100×100，縱橫網(wǎng)格間距各為30 m的平面網(wǎng)格上的合成重力異常在異常源正上方，即y=1500 m處的剖面數(shù)據(jù)，而圖3(c)則是在原始數(shù)據(jù)上添加了20%高斯隨機(jī)噪聲。利用式(6)、式(7)、式(8)對(duì)整個(gè)數(shù)據(jù)空間進(jìn)行計(jì)算得到圖3(b)、(d)。

圖2 組合模型立體圖Fig．2 Perspective of 3D model consisting of two blocks

圖3 模型RMSm切片圖及模型的重力異常Fig．3 The slice maps of the synthetic model and the profiles of its gravity anomalies(a)重力數(shù)據(jù)在y=1500 m處的剖面; (b)在y=1500 m處RMSm 深度切片 ;(c) 添加了20%高斯隨機(jī)噪聲的重力數(shù)據(jù)在y=1500 m處的剖面;(d) 添加噪聲之后的重力數(shù)據(jù)在y=1500 m處RMSm 深度切片，白色虛線框是實(shí)際異常源的真實(shí)輪廓

圖3(d)是y=1500 m處RMSm成像結(jié)果。從圖3可以看出，本文方法所估算出的異常源位置與白色虛線框所標(biāo)示出的異常源的實(shí)際位置非常吻合，所以本文所提出的方法在理論上是有效的。為了說(shuō)明本算法對(duì)于噪聲的穩(wěn)健性，利用添加20%隨機(jī)高斯噪聲的合成重力數(shù)據(jù)進(jìn)行驗(yàn)證，如圖3(c)所示。圖3(d)即為添加噪聲之后所得的y=1500 m處RMSm成像結(jié)果，從圖中可知，該算法具有良好的抗噪性能。

所以從以上算例分析中發(fā)現(xiàn)，本文所提出的成像方法對(duì)于異常源真實(shí)位置的反映是有效的且準(zhǔn)確的，而且具有良好的抗噪性能。

圖4 某地實(shí)測(cè)的布格重力異常的等值線填充圖Fig．4 The contour map of the measured bouguer gravity anomaly in a certain region

4 實(shí)測(cè)資料計(jì)算

圖4是某地地面重力測(cè)量所得的布格重力異常的等值線圖，測(cè)區(qū)范圍為73 km×73 km，平面數(shù)據(jù)網(wǎng)格為100×100，橫縱網(wǎng)格間距為730 m，對(duì)其進(jìn)行線性回歸的分析并計(jì)算RMSm，成像深度大概為0 km～18 km，深度步長(zhǎng)為720 m，對(duì)其RMSm結(jié)果進(jìn)行三維切片成圖，如圖5所示，圖中白色標(biāo)記為異常位置，而且極大值正好對(duì)應(yīng)了負(fù)重力異常，而圖中兩處極小值則對(duì)應(yīng)了兩處正重力異常，圖5的計(jì)算結(jié)果不僅定性地標(biāo)記異常源的賦存范圍，而且也確定了異常源的深度，因此基于線性回歸的重力場(chǎng)源分析的成像結(jié)果，將為測(cè)區(qū)的地質(zhì)解釋提供參考。

圖5 該實(shí)測(cè)異常計(jì)算后得到的相關(guān)系數(shù)的三維切片F(xiàn)ig．5 Three-dimensional slice maps of the intercept related to the measured bouguer gravity anomaly

5 結(jié)論

在概率成像和重磁異常對(duì)應(yīng)分析的基礎(chǔ)上，提出了一種全新的基于線性回歸的重力數(shù)據(jù)三維RMS成像算法。該成像算法具有方法簡(jiǎn)單，分辨率高，抗噪性能強(qiáng)等優(yōu)點(diǎn)，而且在沒(méi)有任何先驗(yàn)信息約束的情況下可以得到穩(wěn)定的計(jì)算結(jié)果。模型算例和實(shí)測(cè)數(shù)據(jù)驗(yàn)證分析表明，本文所提出的算法可以有效的定位地下重力異常源的位置和賦存范圍

參考文獻(xiàn)：

[1] NABIGHIAN M N, ANDER M E, GRAUCH V J S, et al. 75th Anniversary: Historical development of the gravity method in exploration[J]. Geophysics,2005, 70(6): 63N-89N.

[2] 管志寧，郝天珧. 21 世紀(jì)重力與磁法勘探的展望[J]. 地球物理學(xué)進(jìn)展,2002, 17(2): 237-244.

[3] SILVA DIAS F J, BARBOSA V C, SILVA J B. 3D gravity inversion through an adaptive-learning procedure[J]. Geophysics,2009, 74(3): 19-121.

[4] 郭良輝，孟小紅，石磊，等. 重力和重力梯度數(shù)據(jù)三維相關(guān)成像[J]. 地球物理學(xué)報(bào),2009, 52(4): 1098-1106.

[5] 曾華霖. 重力場(chǎng)與重力勘探[M]. 北京: 地質(zhì)出版社, 2005.

[6] 姚長(zhǎng)利，郝天珧，管志寧，等. 重磁遺傳算法三維反演中高速計(jì)算及有效存儲(chǔ)方法技術(shù)[J]. 地球物理學(xué)報(bào),2003, 46(2):252-258.

[7] NABIGHIAN M N, GRAUCH V, HANSEN R O, et al. The historical development of the magnetic method in exploration[J]. Geophysics,2005, 70(6): 33N-61N.

[8] BLAKELY R J, HASSANZADEH S. Estimation of depth to magnetic source using maximum entropy power spectra, with application to the Peru-Chile Trench[J]. Nazca Plate. 1981: 667-682.

[9] PHILLIPS J D. ADEPT: a program to estimate depth to magnetic basement from sampled magnetic profiles[M]. US Geological Survey, 1978.

[10] CHANDLER V W, KOSKI J S, HINZE W J, et al. Analysis of multisource gravity and magnetic anomaly data sets by moving-window application of Poisson’s theorem[J]. Geophysics,1981, 46(1): 30-39.

[11] 管志寧. 地磁場(chǎng)與磁力勘探[M]. 北京: 地質(zhì)出版社, 2005.

[12] MAURIELLO P, PATELLA D. Localization of magnetic sources underground by a data adaptive tomographic scanner[J]. arXiv preprint physics/0511192,2005.

[13] MAURIELLO P, PATELLA D. Imaging polar and dipolar sources of geophysical anomalies by probability tomography. Part II: Application to the Vesuvius volcanic area[J]. arXiv preprint physics/0602057,2006.

[14] ALAIA R, PATELLA D, MAURIELLO P. Imaging multipole gravity anomaly sources by 3D probability tomography[J]. Journal of Geophysics and Engineering,2009, 6(3): 298-310.

[15] IULIANO T, MAURIELLO P, PATELLA D. Looking inside Mount Vesuvius by potential fields integrated probability tomographies[J]. Journal of Volcanology and Geothermal Research,2002, 113(3): 363-378.

[16] CHEN Z, MENG X, GUO L, et al. GICUDA: A parallel program for 3D correlation imaging of large scale gravity and gravity gradiometry data on graphics processing units with CUDA[J]. Computers & Geosciences,2012，46:119-128.

[17] 郭良輝，孟小紅，石磊. 磁異常 ΔT三維相關(guān)成像[J]. 地球物理學(xué)報(bào),2010, 53(2):435-441.

[18] GUO L, MENG X, SHI L. 3D correlation imaging of the vertical gradient of gravity data[J]. Journal of Geophysics and Engineering, 2011, 8(1): 6-12.

[19] GUO L, SHI L, MENG X. 3D correlation imaging of magnetic total field anomaly and its vertical gradient[J]. Journal of Geophysics and Engineering, 2011, 8(2): 287-293.

[20] FEDI M, PILKINGTON M. Understanding imaging methods for potential field data[J]. Geophysics, 2012, 77(1): G13-G24.

[21] 韓旭里，王家寶，陳亞力，等. 大學(xué)數(shù)學(xué)教程-概率論與數(shù)理統(tǒng)計(jì)[M]. 北京：科學(xué)出版社, 2006.