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

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

任意密度分布復(fù)雜地質(zhì)體重力異?？焖偃S正演方法

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

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

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

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

0 引言

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

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

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

空間域法主要是推導(dǎo)規(guī)則幾何形體的解析解，對(duì)復(fù)雜結(jié)構(gòu)形態(tài)和非均勻密度分布的地質(zhì)體采用已有解析解的幾何體的模型剖分或近似，通過組合求和得到整個(gè)復(fù)雜形體的場(chǎng)值[6]。利用積分求解析解的優(yōu)點(diǎn)是可精確地求出空間任意點(diǎn)的場(chǎng)值，缺點(diǎn)在于不同的模型的解析式往往較為復(fù)雜，推導(dǎo)過程比較繁瑣，對(duì)于密度連續(xù)變化的模型需要剖分很多段才能算準(zhǔn)，當(dāng)形體復(fù)雜、需要計(jì)算大量異常數(shù)據(jù)時(shí)，計(jì)算量大，速度較慢。

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

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

1 新方法原理

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

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

快速、準(zhǔn)確計(jì)算組合棱柱體產(chǎn)生的重力異常，是我們研究的重點(diǎn)。組合棱柱體的重力異常gz可以表示為:

(1)

筆者提出的波數(shù)域重力異?？焖佟⒏呔热S正演方法，關(guān)鍵環(huán)節(jié)之一是對(duì)式(1)兩邊施加二維傅里葉變換，推導(dǎo)得到一種新的重力異常波數(shù)域表達(dá)式，將式(1)左右兩邊進(jìn)行二維傅里葉變換，可得式(2)。

(2)

(3)

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

在式(2)的推導(dǎo)過程中，利用了式(4)與式(5)兩個(gè)中間結(jié)論。

(4)

式中：sign(·)表示符號(hào)函數(shù)。

(5)

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

(6)

2 數(shù)值算法

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

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是區(qū)間[-1,1]上高斯點(diǎn)[16]的坐標(biāo)；Aa、Ab為相應(yīng)的高斯點(diǎn)加權(quán)系數(shù)。

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

(15)

式中:

(16)

式中：大括號(hào)的部分利用FFT算法實(shí)現(xiàn)。

空間域重力異常數(shù)值計(jì)算表達(dá)式為式(17)。

(17)

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

在實(shí)際解決具體的復(fù)雜地質(zhì)體的重力異常三維正演時(shí)，實(shí)現(xiàn)本方法的具體流程如下：

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

2)根據(jù)地質(zhì)體密度分布給每個(gè)剖分的小棱柱體的剩余密度進(jìn)行賦值。

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

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

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

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

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

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

3 模型算例

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

(18)

模型的剖分區(qū)域：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方向?yàn)?300 m～300 m,z方向100 m～300 m；圍巖密度均為0 kg/m3，棱柱異常體密度為1 000 kg/m3；高斯FFT算法中x方向和y方向選用高斯點(diǎn)均為4個(gè)。

筆者提出的數(shù)值算法均用Fortran語言實(shí)現(xiàn)，并在CPU主頻2.80GHz、內(nèi)存為32GB的個(gè)人筆記本電腦運(yùn)行。

3.1 重力異常

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

圖4 均方誤差和計(jì)算時(shí)間隨棱柱體剖分個(gè)數(shù)變化Fig.4 Curves of RMSE and calculation time with the varied number of subdivided prisms

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

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

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

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

圖8 起伏觀測(cè)面Fig.8 Rugged observation surface

圖9 均方誤差和計(jì)算時(shí)間隨插值平面?zhèn)€數(shù)變化圖Fig.9 Curves of RMSE and calculation time with the varied number of interpolation planes

3.2 觀測(cè)重力異常

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

圖9顯示了模型剖分個(gè)數(shù)為100*100*50時(shí)，起伏觀測(cè)面的重力異常均方誤差和計(jì)算時(shí)間隨差值平面?zhèn)€數(shù)的變化，從圖9中可以看出，計(jì)算精度隨著剖分個(gè)數(shù)的增加而提高，但當(dāng)插值平面?zhèn)€數(shù)達(dá)到一定時(shí)，計(jì)算精度幾乎不變；同時(shí)，計(jì)算時(shí)間隨插值平面?zhèn)€數(shù)的增加逐漸增多，在插值平面為15個(gè)時(shí)，均方根誤差為0.000 14 mGal，計(jì)算時(shí)間約為4.7 s,計(jì)算效率高。兩條曲線的交點(diǎn)可認(rèn)為是計(jì)算精度與計(jì)算時(shí)間之間一個(gè)很好的折中點(diǎn)，所對(duì)應(yīng)的插值平面?zhèn)€數(shù)可認(rèn)為是最佳個(gè)數(shù)。本算例中最優(yōu)的插值平面?zhèn)€數(shù)約為15。

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

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

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

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

4 結(jié)論

針對(duì)任意密度分布復(fù)雜地質(zhì)體重力異常正演問題，筆者提出一種基于高斯FFT法的波數(shù)域算法，在數(shù)值實(shí)現(xiàn)過程中，采用了Gauss-FFT法，保證計(jì)算效率的同時(shí)，有效克服了傳統(tǒng)FFT法引起的邊界效應(yīng)和零波數(shù)點(diǎn)的奇異問題。

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

[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] 何昌禮, 鐘本善. 復(fù)雜形體的高精度重力異常正演方法[J]. 物探化探計(jì)算技術(shù), 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] 柴玉璞.偏移抽樣理論及其應(yīng)用[M]. 北京: 石油工業(yè)出版社，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]. 物探化探計(jì)算技術(shù), 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]吳文鸝, 管志寧, 高艷芳, 等. 重磁異常數(shù)據(jù)三維人機(jī)聯(lián)作模擬[J]. 物探化探計(jì)算技術(shù), 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]陳龍偉, 戴世坤, 吳美平. 應(yīng)用任意采樣點(diǎn)數(shù)FFT算法時(shí)離散頻率計(jì)算[J]. 地球物理學(xué)進(jìn)展, 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].北京：科學(xué)出版社,1994. XU S Z.Finite element method in Geophysics[M].Beijing:Science Press,1994. (In Chinese)

[17]陳召曦, 孟小紅, 郭良輝,等. 基于GPU并行的重力、重力梯度三維正演快速計(jì)算及反演策略[J]. 地球物理學(xué)報(bào), 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月物探化探計(jì)算技術(shù)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

國(guó)家自然科學(xué)基金項(xiàng)目(41404106)

陳輕蕊(1992-)，女，碩士，從事重磁、電磁數(shù)值模擬研究，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