任意密度分布復雜地質體重力異常快速三維正演方法

陳輕蕊， 陳龍偉， 戴世坤， 張錢江， 歐陽芳

(1.中南大學 地球科學與信息物理學院，長沙 410083; 2.桂林理工大學 地球科學學院，桂林 541006)

任意密度分布復雜地質體重力異常快速三維正演方法

陳輕蕊1， 陳龍偉2， 戴世坤1， 張錢江1， 歐陽芳1

(1.中南大學 地球科學與信息物理學院，長沙 410083; 2.桂林理工大學 地球科學學院，桂林 541006)

解決任意密度分布復雜地質體重力異常三維正演快速、高精度計算問題，是實現重力三維反演、人機交互解釋建模的關鍵。針對該問題，從積分方程出發，提出一種波數域重力異常三維正演方法，其關鍵環節包括三個方面：①將研究區域剖分成許多規則小棱柱體，每個小棱柱體密度值可以任意給定，以此刻畫任意密度分布和起伏地形條件下的復雜地質體；②給出一種新的高精度均勻棱柱體重力異常二維波數域的計算公式，用于計算組合棱柱體模型的重力異常；③采用Gauss-FFT法將重力異常從波數域轉換到空間域，保證計算效率的同時，有效克服了傳統FFT法引起的邊界效應問題。模型算例檢驗結果表明，該算法計算速度快、精度高，對于剖分為百萬個棱柱體的模型，耗時只需幾秒。

重力異常； 三維正演； Gauss-FFT； 積分方程

0 引言

隨著觀測數據的急劇增加，重力反演已由近似的二維反演發展到符合實際情況的三維反演。正演是反演的基礎，任意密度分布復雜地質體重力異常快速、高精度三維正演是實現重力三維反演、人機交互解釋建模的關鍵。解決重力異常三維正演問題，既可以從偏微分方程邊值問題入手，也可以從三維積分入手，兩者理論上是等價的，但由兩種思路得到的數值方法在計算效率和計算精度兩方面的表現是不同的。

從偏微分方程的邊值問題角度來求解三維正演問題，主要采用的數值方法包括有:①限單元法;②有限差分法;③有限體元法。微分方法的關鍵環節主要包括:①邊值條件給定;②稀疏矩陣方程求解。CAI Y[1]提出了一種新的邊界條件，提高了數值結果精度，并通過采用一種加權函數策略來直接求解節點上的重力值，避免了傳統有限單元法中先計算高斯點重力值再通過插值計算節點重力值，計算效率得到提高；C G Farquharson[2]從重力場高斯定理出發，采用有限差分法求解泊松方程，由于滿足其次邊界條件，使用時計算區域邊界需要遠于實際異常區域邊界，提高計算精度的同時降低了計算效率；May D A[3]對比研究了直接累加法、有限元法和快速多偶極子法；Jahandari H[4]采用了非結構化網格有限元法求解重力偏微分方程；Guzman S[5]采用有限體元法直接對重力場進行求解。偏微分方法的一個優勢是能夠一次計算所有節點上的場值，生成的系數矩陣為稀疏矩陣，存儲要求低，但目前這類方法的計算效率還是較低。

用積分方程求解重力三維正演問題的研究相對更多一些，大致可分為：①空間域方法；②波數域方法。

空間域法主要是推導規則幾何形體的解析解，對復雜結構形態和非均勻密度分布的地質體采用已有解析解的幾何體的模型剖分或近似，通過組合求和得到整個復雜形體的場值[6]。利用積分求解析解的優點是可精確地求出空間任意點的場值，缺點在于不同的模型的解析式往往較為復雜，推導過程比較繁瑣，對于密度連續變化的模型需要剖分很多段才能算準，當形體復雜、需要計算大量異常數據時，計算量大，速度較慢。

波數域方法將重力積分進行二維或三維傅里葉變換[7]，借助FFT算法實現場的高效計算。由于離散采樣等原因，傳統FFT存在邊界效應問題，導致計算精度低，通常需要進行擴邊來減小邊界效應，提高計算精度，但這使得計算量增大。Wuly[8]在偏移抽樣理論[9]的基礎上，提出一種實現傅里葉變換高精度數值計算的Gauss-FFT方法。該方法能有效抑制邊界效應，且能避免零波數計算時的奇異問題，十分適合于波數域方法。

筆者從積分方程出發，提出一種波數域三維正演方法。該方法給出了一種新的組合棱柱體模型波數域重力異常表達式，借助Gauss-FFT方法，很好解決了任意密度分布復雜地質體快速、高精度重力異常三維正演問題，為解決大規模重力異常三維反演和人機交互解釋提供了有力的方法支持。

1 新方法原理

為了模擬任意密度分布三維復雜地質體，我們采用的策略是將研究區域剖分成很多個規則小棱柱體，每個棱柱體的密度為常值，剖分棱柱體個數越多，模型刻畫越精細。這種策略十分適合于重力異常三維反演[10-12]和人機交互解釋建模[13-14]。依據重力異常的疊加原理，復雜地質體的重力異常計算問題，可以轉化為組合棱柱體的重力異常的計算問題(圖1)。

圖1 組合棱柱體示意圖Fig.1 Illustration of combined prism model

快速、準確計算組合棱柱體產生的重力異常，是我們研究的重點。組合棱柱體的重力異常gz可以表示為:

(1)

筆者提出的波數域重力異常快速、高精度三維正演方法，關鍵環節之一是對式(1)兩邊施加二維傅里葉變換，推導得到一種新的重力異常波數域表達式，將式(1)左右兩邊進行二維傅里葉變換，可得式(2)。

(2)

(3)

式(3)在形式上類似于二維離散傅里葉變換，在數值計算實現時，可利用快速傅里葉算法(FFT)實現。

在式(2)的推導過程中，利用了式(4)與式(5)兩個中間結論。

(4)

式中：sign(·)表示符號函數。

(5)

從理論上來說，空間域重力異常可通過二維傅里葉反變換得到。

(6)

2 數值算法

這里給出基于Gauss-FFT法的重力異常數值計算公式。對式(6)離散過程中，需要對空間域坐標和波數域坐標分別進行離散化，我們采用文獻[15]中的方法，給出任意采樣點情況下離散空間域坐標和離散波數值(以采樣點個數為偶數為例，奇數情況離散波數)如下:

xm=x1+(m-1)Δxm=1,2,…,Nx

(7)

yn=y1+(n-1)Δyn=1,2,…,Ny

(8)

kp=p·Δkx

(9)

kq=q·Δky

(10)

(11)

(12)

(13)

(14)

式中：ta、tb是區間[-1,1]上高斯點[16]的坐標；Aa、Ab為相應的高斯點加權系數。

根據Gauss-FFT法[8]離散傅里葉變換計算公式，可得波數域重力異常計算式為式(15)。

(15)

式中:

(16)

式中：大括號的部分利用FFT算法實現。

空間域重力異常數值計算表達式為式(17)。

(17)

分析式(17)可以發現，公式中大括號的部分可用FFT算法實現，這樣保證了計算效率，也提高了計算精度。由于不用計算零波數，避免了零波數造成的奇異問題。但是Gauss-FFT需要用Mx·My次FFT，時間代價增大，可用通過采用并行計算技術解決。

在實際解決具體的復雜地質體的重力異常三維正演時，實現本方法的具體流程如下：

1)給定剖分區域，確定剖分棱柱體個數Nx、Ny、Nz，水平方向均勻剖分，垂直方向可可非均勻剖分；在研究含有起伏地形的情況時，起伏地形要完全包含在剖分區域內；給定Gauss-FFT反變換時高斯點個數Mx、My，高斯坐標ta、tb和高斯系數Aa、Ab。

2)根據地質體密度分布給每個剖分的小棱柱體的剩余密度進行賦值。

4)根據式(16)，利用FFT算法計算密度波數域值。

6)利用式(17)和FFT算法，計算空間域重力異常。

7)計算起伏觀測面(圖2)或者起伏地形(圖3)上的重力異常時，采用如下策略。

采用步驟1)～步驟7)的方法，計算起伏觀測面最高點和最低點之間的多個平面的異常值，然后采用三次樣條差值方法，插值得到起伏觀測面上的場值。

圖2 起伏觀測面和插值平面示意圖Fig.2 Rugged observation surface and interpolation planes

圖3 起伏地形和插值平面示意圖Fig.3 Rugged topography and interpolation planes

3 模型算例

為了檢驗筆者提出的算法在解決任意密度分布情況下復雜地質體重力異常三維正演問題的計算速度和計算精度，設計了棱柱體異常模型，分別計算平面和起伏面上的重力異常值。兩種觀測面上的重力異常都可通過解析公式計算得到，用筆者提出的算法計算得到的重力異常數值解與解析解進行對比，以此驗證算法的計算精度,并用式(18)給出的均方誤差進行衡量。

(18)

模型的剖分區域：x方向-1 000 m～1 000 m,y方向-1 000 m～1 000 m,z方向0 m～500 m；棱柱體異常的位置為：x方向-300 m～300 m,y方向為-300 m～300 m,z方向100 m～300 m；圍巖密度均為0 kg/m3，棱柱異常體密度為1 000 kg/m3；高斯FFT算法中x方向和y方向選用高斯點均為4個。

筆者提出的數值算法均用Fortran語言實現，并在CPU主頻2.80GHz、內存為32GB的個人筆記本電腦運行。

3.1 重力異常

圖4顯示了觀測面z=0 m 的平面重力異常均方誤差和計算時間隨棱柱體剖分個數的變化，從圖中可以看出，計算精度隨著剖分個數的增加越來越高，達到一定的剖分個數時，計算精度提高緩慢；計算時間隨剖分個數的增加而增多，在剖分為100*100*100時，均方根誤差為0.000 13 mGal，計算時間約為1.3 s,計算效率高。兩條曲線的交點可認為是計算精度與計算時間之間一個很好的折中點，所對應的剖分個數可認為是最佳剖分個數。本算例中最優的剖分個數約為100*100*50。

圖4 均方誤差和計算時間隨棱柱體剖分個數變化Fig.4 Curves of RMSE and calculation time with the varied number of subdivided prisms

圖5和圖6分別給出的是剖分個數為100*100*50時，z=0 m平面重力異常的解析解和數值解。對比圖5與圖6，兩者形態幾乎一致，解析解與數值解作差，誤差如圖7所示。從圖7中可看出，總體誤差很小(誤差與解析解相差3個～4個)，相比采用傳統FFT的波數域方法，筆者基于Gauss-FFT法的波數域方法邊界誤差很小，具有很高的計算精度。

圖5 平面重力異常解析解Fig.5 Analytical solution of gravity anomalies on a plane

圖6 平面重力異常數值解Fig.6 Numerical solution of gravity anomalies on a plane

圖7 平面重力異常誤差圖Fig.7 Computation errors of gravity anomalies on a plane

圖8 起伏觀測面Fig.8 Rugged observation surface

圖9 均方誤差和計算時間隨插值平面個數變化圖Fig.9 Curves of RMSE and calculation time with the varied number of interpolation planes

3.2 觀測重力異常

如圖8所示，筆者采用簡化的起伏觀測面，該觀測面由z=-100 m、z=0 m 、 z=100 m三個平面組合而成，插值平面在z=-110 m～110 m之間。

圖9顯示了模型剖分個數為100*100*50時，起伏觀測面的重力異常均方誤差和計算時間隨差值平面個數的變化，從圖9中可以看出，計算精度隨著剖分個數的增加而提高，但當插值平面個數達到一定時，計算精度幾乎不變；同時，計算時間隨插值平面個數的增加逐漸增多，在插值平面為15個時，均方根誤差為0.000 14 mGal，計算時間約為4.7 s,計算效率高。兩條曲線的交點可認為是計算精度與計算時間之間一個很好的折中點，所對應的插值平面個數可認為是最佳個數。本算例中最優的插值平面個數約為15。

當插值平面個數為10時，圖10和圖11分別給出起伏觀測面上重力異常的解析解和數值解。對比圖10與圖11，兩者形態幾乎一致，解析解與數值解作差，誤差如圖12所示。從圖12中可看出，總體誤差很小(誤差與解析解相差3個～4個數量級)，大的誤差主要集中在異常體區域以及起伏觀測凹凸面的邊界，說明了本文算法在計算起伏觀測面的重力異常時也具有很高的精度。對比圖11和圖6可以看出，起伏地形上的重力異常比起平面上的場存在變化，在起伏凹凸面上變化形態明顯，凹面的異常值變大，凸面的異常場變小。

圖10 起伏觀測面重力異常解析解Fig.10 Analytical solution of gravity anomalies on a rugged surface

圖11 起伏觀測面重力異常數值解Fig.11 Numerical solution of gravity anomalies on a rugged surface

圖12 起伏觀測面重力異常誤差圖Fig.12 Computation errors of gravity anomalies on a rugged surface

4 結論

針對任意密度分布復雜地質體重力異常正演問題，筆者提出一種基于高斯FFT法的波數域算法，在數值實現過程中，采用了Gauss-FFT法，保證計算效率的同時，有效克服了傳統FFT法引起的邊界效應和零波數點的奇異問題。

筆者提出的算法對于剖分規模為128*128*64的算例，觀測點為128*128個時，在CPU主頻為1.90 GHz、內存為4 GB的電腦上運行的計算時間約為38 s。筆者提出的算法的計算速度在沒有并行的前提下，在主頻更低的計算機上的計算速度，比前人提出的方法[16]快了約2.6倍。由此可以說明我們提出的算法速度快，精度高。此外，本算法具有很好的并行性，結合并行計算技術，可以進一步提高本方法的計算效率，特別適合應用于三維重力反演和人機交互建模。

[1] CAI Y,WANG C.Fast finite-element calculation of gravity anomaly in complex geological regions [J].Geophysical Journal of International,2005,162:696-708.

[2] C G FARQUHARSON, C R W MOSHER. Three-dimensional modelling of gravity data using finite differences[J]. Journal of Applied Geophysics, 2009, 68(3): 417-422.

[3] MAY D A, KNEPLEY M G. Optimal, scalable forward models for computing gravity anomalies[J]. Geophysical Journal International, 2011, 187(1): 161-177.

[4] JAHANDARI H, FARQUHARSON C G. Forward modeling of gravity data using finite-volume and finite-element methods on unstructured grids[J]. Geophysics, 2013, 78(3): 69-G80.

[5] GUZMAN S. Forward modeling and inversion of potential field data using partial differential equations [D]. MS thesis of Colorado School of Mines, 2015.

[6] 何昌禮, 鐘本善. 復雜形體的高精度重力異常正演方法[J]. 物探化探計算技術, 1988,10(2): 121-128. HE C L, ZHONG B S, A high accuracy forward method for gravity anomaly of complex body [J]. Computing Techniques for Geophysical and Geochemical Exploration, 1988，10(2): 121-128.(In Chinese)

[7] TONTINI F C,COCCHI L,CARMISCIANO C.Rapid 3-D forward model of potential fields with application to the Palinuro Seamount magnetic anomaly (southern Tyrrhenian Sea, Italy)[J].Journal of Geophysical Research Solid Earth,2009,114(B2):1205-1222.

[8] WU L Y , TIAN G. High precision Fourier forward modeling of potential fields[J]. Geophysics, 2014, 79(5): G59-G68.

[9] 柴玉璞.偏移抽樣理論及其應用[M]. 北京: 石油工業出版社，1997. CHAI Y P.Shift Sampling Theory and its Application[M].Beijing: China Petroleum Industry Press,1997. (In Chinese)

[10]LI Y G, OLDENBURG DW. 3-D inversion of gravity data[J]. Geophysics, 1998, 63(1): 109-119.

[11]ZHDANOV M S, R ELLIS, S MUKHERJEE. Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 2004. 69(4): 925-937.

[12]周宇軒, 雷宛. 重力異常物性正則化反演研究[J]. 物探化探計算技術, 2016,38(5): 579-583. ZHOU Y X, LEI W, Study on physical properties regularization inversion anomaly[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2016,38(5): 579-583. (In Chinese)

[13]XIA H. Interactive modeling of potential fields in three dimensions[J]. SEG Expanded Abstracts of the 63th Annual Meeting,1993:403-404.

[14]吳文鸝, 管志寧, 高艷芳, 等. 重磁異常數據三維人機聯作模擬[J]. 物探化探計算技術, 2005, 27(3): 227-232. WU W L, GUAN Z N, GAO Y F, et al, Interactiv man/computer moding 3D body of gravity and magnetic anomaly data[J]. Computing Techniques for Geophysical and Geochemical Exploration, 2005, 27(3): 227-232. (In Chinese)

[15]陳龍偉, 戴世坤, 吳美平. 應用任意采樣點數FFT算法時離散頻率計算[J]. 地球物理學進展, 2016(1):164-169. CHEN L W, DAI S K，WU M P. Computation of discrete frequency when using FFT algorithm with random sampling points[J]. Progress in Geophysics, 2016(1):164-169. (In Chinese)

[16]徐世浙.地球物理中的有限單元法[M].北京：科學出版社,1994. XU S Z.Finite element method in Geophysics[M].Beijing:Science Press,1994. (In Chinese)

[17]陳召曦, 孟小紅, 郭良輝,等. 基于GPU并行的重力、重力梯度三維正演快速計算及反演策略[J]. 地球物理學報, 2012, 55(12):4069-4077. CHEN Z X, MENG X H, GUO L H, et al. Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU[J]. Chinese Journal of Geophysics- Chinese Edition, 2012, 55(12):4069-4077.(In Chinese)

第39卷 第2期2017年3月物探化探計算技術COMPUTING TECHNIQUES FOR GEOPHYSICAL AND GEOCHEMICAL EXPLORATIONVol.39 No.2Mar. 2017

Three-dimensional forward modeling method for gravity anomalies of complex geological bodies with arbitrary density distribution

CHEN Qingrui1, CHEN Longwei2, DAI Shikun1, ZHANG Qianjiang1, OUYANG fang1

(1.School of Geosciences and Info-Physics, Central South University, Changsha 410083, China; 2.College of Earth Sciences, Guilin University of Technology, Guilin 541006,China )

Rapid three-dimensional (3D) forward modeling of gravity anomalies of complex geologic bodies with arbitrary density distribution plays an important role on 3D inversion of gravity and human-computer interaction modeling and interpretation. This paper proposes a wavenumber domain 3D forward modeling method. This new method has three key aspects: (1) prism based subdivision strategy. The research area is divided into many small regular prisms, density for each small prism can be given arbitrarily, in order to simulate complex geological body with arbitrary density distribution and undulating topography; (2) new wavenumber domain formula for computing gravity anomaly of combined prism model; (3) Gauss-FFT method. This paper uses Gauss-FFT method to transform the gravity anomaly from the wave number domain to the spatial domain. The Gauss-FFT method not only has high computational efficiency, but also suppresses the boundary effect and avoids the singular problem at zero wavenumber caused by the traditional FFT method. Numerical tests show that the proposed method is fast and accurate. It takes only a few seconds to calculate gravity anomaly for the combined 3D model with millions of prisms.

gravity anomalies; three-dimensional forward modeling; Gauss-FFT; integral equation

2016-11-29 改回日期：2017-01-09

國家自然科學基金項目(41404106)

陳輕蕊(1992-)，女，碩士，從事重磁、電磁數值模擬研究，E-mail：1532550316@qq.com。

陳龍偉(1982-)，男，博士，從事重磁勘探、電磁勘探研究，E-mail：longweichen_glut@glut.edu.cn。

1001-1749(2017)02-0176-07

P 631.1

A

10.3969/j.issn.1001-1749.2017.02.04