基于線性回歸的重力數據三維RMS成像

吉日嘎拉圖，魯光銀

(中南大學 地球科學與信息物理學院, 長沙 410083)

0 引言

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

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

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

1 方法原理

地下異常源在，水平測區的所產生的重力異常用公式(1)表示：

(1)

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

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

(2)

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

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

(3)

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

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

(4)

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

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

(5)

所以地上所測重力異常與地下第i個點所對應的異常源幾何構造函數存在著線性的對應關系，在進行線性回歸分析之前，必須進行相關性檢驗，驗證是否是存在線性相關性。

對任意一組觀測值與幾何構造函數，不管是否相關都可以建立一個線性回歸方程，這就產生了一個問題，這樣得到的線性回歸方程是否真的具有線性關系，而且所求的線性回歸方程是不是能夠反映它們之間的實際關系，因此進行線性之前，必須進行線性相關關系檢驗[6]。

根據統計學原理[21]，相關系數被定義為式(6)：

(6)

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

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

(7)

根據均方根殘差法(RMS)原理，對其進行修正，用相關系數R當作ΔGcalc的權，相關程度越高，擬合程度越接近，才能越接近真實觀測值，因此得到了等式(8)加入相關系數R修正后的RMSm，因為相關系數-1≤R≤1，故對于式(8)，只需要尋找局部極值(包括極大值和極小值)就可以準確找到異常源的所處位置。

RMSmodified=

(8)

圖1 利用線性回歸進行場源分析的流程圖Fig．1 The flow chart of processing 3-D RMSm imaging

2 基于線性回歸的重力數據三維RMS成像

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

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

3 算例分析

為了研究該方法的具體表現，我們利用兩個密度為0.5 g/cm3的異常體，設置在地面以下100 m處的位置如 圖 2所示。圖3(a)為測網100×100，縱橫網格間距各為30 m的平面網格上的合成重力異常在異常源正上方，即y=1500 m處的剖面數據，而圖3(c)則是在原始數據上添加了20%高斯隨機噪聲。利用式(6)、式(7)、式(8)對整個數據空間進行計算得到圖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)重力數據在y=1500 m處的剖面; (b)在y=1500 m處RMSm 深度切片 ;(c) 添加了20%高斯隨機噪聲的重力數據在y=1500 m處的剖面;(d) 添加噪聲之后的重力數據在y=1500 m處RMSm 深度切片，白色虛線框是實際異常源的真實輪廓

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

所以從以上算例分析中發現，本文所提出的成像方法對于異常源真實位置的反映是有效的且準確的，而且具有良好的抗噪性能。

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

4 實測資料計算

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

圖5 該實測異常計算后得到的相關系數的三維切片Fig．5 Three-dimensional slice maps of the intercept related to the measured bouguer gravity anomaly

5 結論

在概率成像和重磁異常對應分析的基礎上，提出了一種全新的基于線性回歸的重力數據三維RMS成像算法。該成像算法具有方法簡單，分辨率高，抗噪性能強等優點，而且在沒有任何先驗信息約束的情況下可以得到穩定的計算結果。模型算例和實測數據驗證分析表明，本文所提出的算法可以有效的定位地下重力異常源的位置和賦存范圍

參考文獻：

[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]. 地球物理學進展,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] 郭良輝，孟小紅，石磊，等. 重力和重力梯度數據三維相關成像[J]. 地球物理學報,2009, 52(4): 1098-1106.

[5] 曾華霖. 重力場與重力勘探[M]. 北京: 地質出版社, 2005.

[6] 姚長利，郝天珧，管志寧，等. 重磁遺傳算法三維反演中高速計算及有效存儲方法技術[J]. 地球物理學報,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] 管志寧. 地磁場與磁力勘探[M]. 北京: 地質出版社, 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三維相關成像[J]. 地球物理學報,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] 韓旭里，王家寶，陳亞力，等. 大學數學教程-概率論與數理統計[M]. 北京：科學出版社, 2006.