基于隨機(jī)介質(zhì)模型的GPR無單元法正演模擬

戴前偉，王洪華

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

探地雷達(dá)(Ground penetrating radar, GPR)具有高效、快速、無損、抗干擾能力強(qiáng)的優(yōu)點(diǎn)，廣泛應(yīng)用于工程與環(huán)境地球物理探測的各個領(lǐng)域，成為淺部勘探的重要技術(shù)[1]。GPR正演模擬對數(shù)據(jù)解釋具有重要的指導(dǎo)作用。通過對GPR模型的正演模擬，可以加深對GPR探測剖面的認(rèn)識，提高解釋精度[2]。目前，GPR正演模擬大都以傳統(tǒng)的均勻介質(zhì)為基礎(chǔ)[3]，而在實際的探地雷達(dá)檢測中，地下介質(zhì)中存在大量分布不規(guī)則的微小異常，將會造成大量小的不相干的擾動，其實質(zhì)是源于介質(zhì)在小尺度上的非均勻性，在實際高分辨GPR探測中常常作為“噪聲”進(jìn)行處理[4]。為了精細(xì)地研究地下介質(zhì)的特征，必須對這種小尺度上的非均勻性產(chǎn)生的GPR波場有所了解[5]。顯然，沿用傳統(tǒng)的模型理論很難準(zhǔn)確描述這些不可忽略的微小異常，而以統(tǒng)計學(xué)理論為基礎(chǔ)的隨機(jī)介質(zhì)模型能較好地描述這種小尺度上的非均勻性特征[6]。

國內(nèi)外很多學(xué)者對此進(jìn)行了系統(tǒng)而深入的研究，IKELLE等[7]闡述了基于指數(shù)型橢圓自相關(guān)函數(shù)的二維隨機(jī)介質(zhì)模型的建立方法，并討論了隨機(jī)介質(zhì)對彈性波場的影響特征。VARADAN等[8]研究了隨機(jī)介質(zhì)中彈性波的散射和衰減特征；KORN[9]采用時域有限差分法對二維隨機(jī)介質(zhì)模型進(jìn)行了正演計算，并分析了隨機(jī)介質(zhì)對彈性波傳播的影響特征；KNEIB和KERNER[10]進(jìn)行了高精度、快速的隨機(jī)介質(zhì)模型的彈性波正演計算，并給出了隨機(jī)介質(zhì)模型的統(tǒng)計學(xué)描述方法；奚先和姚姚[11-13]系統(tǒng)研究了隨機(jī)介質(zhì)模型的建立方法及其特點(diǎn)，并對各種隨機(jī)介質(zhì)模型進(jìn)行了正演計算；郭乃川等[14]提出了一種拓展了小尺度非均勻性擇優(yōu)取向的隨機(jī)介質(zhì)建模新方法，并利用三維錐形函數(shù)表達(dá)式壓制計算誤差，使得建立的隨機(jī)介質(zhì)模型更具有可信度。朱生旺等[15]采用隨機(jī)介質(zhì)模型方法構(gòu)建了孔洞性油氣儲層模型；陳可洋[16]提出了三維隨機(jī)介質(zhì)建模方法，并對隨機(jī)介質(zhì)中的彈性波波場特征進(jìn)行了分析；王金山等[17]提出了一種局部隨機(jī)位置或固定位置任意形狀截取法并結(jié)合多尺度建模技術(shù)來共同構(gòu)造隨機(jī)介質(zhì)模型的新方法。目前，地球物理勘探中有關(guān)隨機(jī)介質(zhì)模型構(gòu)建及其正演模擬的研究主要集中在地震勘探中。

無單元法(Element free method, EFM)是GPR正演模擬的有效手段。無單元法的基本思想[18]是將計算區(qū)域離散成若干個點(diǎn)，由滑動最小二乘法(Moving least squares，MLS)來擬合場函數(shù)，從而擺脫了單元的束縛，具有更大的靈活性。無單元法由于拋棄了單元的概念，只需節(jié)點(diǎn)信息，及采用滑動最小二乘法構(gòu)造形函數(shù)，使得 EFM 具有前處理簡單、精度高、獨(dú)立變量解高次連續(xù)等優(yōu)點(diǎn)。無單元法已應(yīng)用于GPR正演模擬中，并取得了良好的效果。本文作者在上述理論的基礎(chǔ)上，根據(jù)隨機(jī)介質(zhì)模型構(gòu)建理論，構(gòu)建多尺度的二維GPR隨機(jī)介質(zhì)模型，然后采用無單元法進(jìn)行正演計算，研究了隨機(jī)介質(zhì)中GPR波場特征，并與均勻介質(zhì)的計算結(jié)果進(jìn)行了對比。

1 隨機(jī)介質(zhì)模型

隨機(jī)介質(zhì)模型主要由非均勻性大、小兩種尺度組成。大尺度描述介質(zhì)的背景特性，而小尺度則描述加在背景模型上的隨機(jī)擾動。隨機(jī)介質(zhì)模型通常用一個均值為零的二階平穩(wěn)隨機(jī)過程來描述[19]。在探地雷達(dá)(GPR)探測過程中，GPR波在地下介質(zhì)中傳播時主要受介電常數(shù)ε和介質(zhì)的電阻率σ的影響。以二維隨機(jī)介質(zhì)為例，根據(jù)隨機(jī)介質(zhì)模型的構(gòu)建理論，地下介質(zhì)的參數(shù)可以表示為

式中：ε0、σ0代表背景介質(zhì)參數(shù)，δε、δσ代表上述背景介質(zhì)上的非均勻擾動量，用來描述隨機(jī)介質(zhì)在小尺度上的非均勻性，

假設(shè)空間隨機(jī)擾動φ具有零均值、一定方差及自相關(guān)函數(shù)的空間平穩(wěn)隨機(jī)過程，則

構(gòu)造二階平穩(wěn)過程ε(x,z)，σ(x,z)的步驟如下[19-20]：

1)選擇自相關(guān)函數(shù)。目前，高斯型、指數(shù)型及Von KARMAN型自相關(guān)函數(shù)被廣泛地應(yīng)用于描述隨機(jī)介質(zhì)。高斯型相關(guān)函數(shù)能描述單尺度平滑的隨機(jī)介質(zhì)，指數(shù)型和Von KARMAN型相關(guān)函數(shù)能很好地描述具有多尺度平滑的隨機(jī)介質(zhì)。本文作者使用如下的混合型自相關(guān)函數(shù)：

式中：a、b分別是介質(zhì)在x、z方向上的自相關(guān)長度，p為粗糙度因子。當(dāng)p=0, 1時，φ(x,z)分別對應(yīng)高斯型自相關(guān)函數(shù)和指數(shù)型自相關(guān)函數(shù)。

3)用隨機(jī)數(shù)發(fā)生器生成0, 2π[]φ上服從均勻分布的獨(dú)立的二維隨機(jī)場

4)應(yīng)用下式計算隨機(jī)功率譜：

7)通過規(guī)范化產(chǎn)生均值為零、方差為2α，得到介電常數(shù)的小尺度相對擾動：

8)將式(10)代入式(1)中，即可得到隨機(jī)介質(zhì)模型的相對介電常數(shù)。隨機(jī)介質(zhì)模型的電導(dǎo)率也可按上述方法構(gòu)建。

圖1所示為按照上述方法，分別選擇不同的相關(guān)長度a、b所產(chǎn)生的4個不同特征的隨機(jī)介質(zhì)(相對介電常數(shù))，其中背景相對介電常數(shù)為3，隨機(jī)擾動量的標(biāo)準(zhǔn)差為 10%。從圖1可以看到：自相關(guān)長度a、b可以描述隨機(jī)介質(zhì)擾動的平均尺度，隨機(jī)介質(zhì)模型能靈活、有效地描述地下實際介質(zhì)分布。

2 GPR無單元法正演模擬

MAXWELL方程組描述了電磁場的運(yùn)動學(xué)和動力學(xué)規(guī)律。由電磁波理論，高頻電磁波(GPR)在介質(zhì)中的傳播規(guī)律也應(yīng)服從MAXWELL方程組。以電場為例，根據(jù)文獻(xiàn)[21]的推導(dǎo)，GPR波滿足的波動方程為

式中：ε為介電常數(shù)(F/m)，μ為磁導(dǎo)率(H/m)，σ為電導(dǎo)率(S/m)，E為電場強(qiáng)度(V/m)，SE為電場源，t為時間。

首先定義如下形式的近似解函數(shù)：

圖1 采用混合型橢圓自相關(guān)函數(shù)產(chǎn)生的不同相關(guān)長度的隨機(jī)介質(zhì)相對介電常數(shù)模型Fig.1 Random medium models of relative dielectric permittivity in different auto correlation length media produced with intermixed ellipsoidal autocorrelation function: (a)a=1, b=1; (b)a=1, b=5; (c)a=5, b=5; (d)a=1, b=20

對式(12)加權(quán)求和得：

式中：n為權(quán)函數(shù)w(x-xi)非零域內(nèi)的節(jié)點(diǎn)個數(shù)，在點(diǎn)x周圍一個有限的鄰域內(nèi)，權(quán)函數(shù)w(x-xi)＞0；在這個鄰域之外的權(quán)函數(shù)定義為0，該鄰域叫做點(diǎn)x的影響域[18]，w(x-xi)的大小隨著鄰域點(diǎn)xi遠(yuǎn)離中心點(diǎn)x而逐漸減小。本文中采用指數(shù)權(quán)函數(shù)：

式中：rinf決定影響域的大小，在二維情況下，影響半徑ri是點(diǎn)x與點(diǎn)xi之間的距離，c是一個控制相對權(quán)重的常數(shù)。

對式(13)求最小值可得：

式中：

將式(15)代入式(12)中可得：

式中：Ni(x)為節(jié)點(diǎn)i的形函數(shù)，它是坐標(biāo)的函數(shù)。

與有限元法類似，將式(19)通過變分原理用于式(11), 即得到空間上的離散方程。假設(shè)給定了基本邊界條件, 使用帶罰因子的 GALERKIN 方法[22]可得到如下離散方程：

式中：S為等效的電場源向量，為電場對時間的二次導(dǎo)數(shù)項，為電場對時間的一次導(dǎo)數(shù)項，M為質(zhì)量矩陣，K’為阻尼矩陣，K為剛度矩陣。

式中：?和Γ分別為所討論電場的域及其邊界，a為罰因子。以上各矩陣的求解都需要計算高斯數(shù)值積分。將式(22)在時間域中對加速度項E和速度項E采用中心差分法加以展開，

采用不完全LU分解預(yù)處理的BICGSTAB算法[23]進(jìn)行迭代求解式(26)。為了提高計算效率，采用集中質(zhì)量矩陣和集中阻尼矩陣使得方程組的求解無需對矩陣求逆。

3 數(shù)值模擬實例

設(shè)計了不含異常體的隨機(jī)介質(zhì)、包含矩形異常體隨機(jī)介質(zhì)2種GPR地電模型，應(yīng)用基于隨機(jī)介質(zhì)模型的無單元算法對設(shè)定的模型進(jìn)行了正演計算，分析了隨機(jī)介質(zhì)模型的GPR波散射特征，并與均勻介質(zhì)模型的計算結(jié)果進(jìn)行了對比分析。

3.1 模型1

圖2所示為一個不含異常體的隨機(jī)介質(zhì)模型示意圖，模型為一個10 m×10 m的矩形區(qū)域，其背景相對介電常數(shù)為ε為3.5，背景電導(dǎo)率σ為0.001 S/m，空間網(wǎng)格步長為0.1 m，網(wǎng)格總數(shù)為100×100。GPR波脈沖激勵源的中心頻率為100 MHz，時間步長為0.1 ns，時窗長度為50 ns。從圖2中可以看出，在隨機(jī)介質(zhì)模型中含有幾個小尺度的異常。應(yīng)用無單元法對該模型進(jìn)行正演計算，其模擬所得的GPR正演模擬剖面如圖3所示。由圖3可見，隨機(jī)介質(zhì)中的GPR波由于受到隨機(jī)介質(zhì)的影響散射現(xiàn)象非常強(qiáng)烈。

圖2 隨機(jī)介質(zhì)模型示意圖Fig.2 Sketch map of random medium model

為了更詳細(xì)地說明GPR波在隨機(jī)介質(zhì)中的傳播特征，將GPR波脈沖激勵源放置在模型的中間位置(5.0,5.0)，通過截取不同時刻GPR波場快照圖，觀測隨機(jī)介質(zhì)對GPR波傳播的影響。應(yīng)用無單元法對該隨機(jī)介質(zhì)模型進(jìn)行正演模擬，得到如圖4所示的波場快照，圖4(a)、(b)、(c)所示分別為隨機(jī)介質(zhì)模型中 GPR波10、20、30 ns時刻的波場快照，圖4(d)、(e)、(f)所示為相應(yīng)時刻均勻介質(zhì)(相對介電常數(shù)為 3.5)的波場快照。圖4(a)、(b)、(c)與(d)、(e)、(f)相比，GPR 波波前由于受到小尺度異常散射的影響，波形發(fā)生扭曲。此外，隨機(jī)介質(zhì)中 GPR波在傳播過程中由于受到小尺度異常的影響散射現(xiàn)象非常強(qiáng)烈，并呈無序狀。而均勻介質(zhì)中GPR波波前都是標(biāo)準(zhǔn)的圓形，并且GPR波傳播無散射現(xiàn)象出現(xiàn)。因此，在GPR數(shù)據(jù)處理過程中，準(zhǔn)確認(rèn)識地下介質(zhì)的小尺度異常產(chǎn)生的GPR波散射特征是非常有必要的。

圖3 模型1 GPR正演剖面圖Fig.3 GPR forward compose section of Model 1

圖4 隨機(jī)介質(zhì)模型與均勻介質(zhì)模型的波場快照圖Fig.4 Wave field snapshots of random medium model ((a), (b), (c))and homogeneous medium model ((d), (e), (f)): (a), (d)10 ns;(b), (e)20 ns; (c), (f)30 ns

3.2 模型 2

圖5 模型2示意圖Fig.5 Sketch map of random medium Model 2

圖6 均勻模型2 GPR正演剖面圖Fig.6 GPR forward compose section of Model 2

圖7 矩形異常體中相對介電常數(shù)不同時的GPR正演模擬剖面圖Fig.7 GPR forward compose section of rectangular anomalies with different relative dielectric permittivities: (a)6.0; (b)8.0;(c)10.0; (d)15.0

圖5所示為含有矩形異常體的隨機(jī)介質(zhì)模型示意圖，模型為一個10 m×6 m的矩形區(qū)域，其背景相對介電常數(shù)ε為3.5，背景電導(dǎo)率σ為0.001 S/m，隨機(jī)介質(zhì)模型在背景相對介電常數(shù)和背景電導(dǎo)率上進(jìn)行擾動。坐標(biāo)(5.0, 2.0)的位置有一個大小為0.6 m×0.4 m的矩形異常體，其相對介電常數(shù)ε1分別為6.0、8.0、10.0和15.0；電導(dǎo)率σ1為0.01 S/m?？臻g網(wǎng)格步長為0.1 m，網(wǎng)格總數(shù)為100×60。GPR波脈沖激勵源的中心頻率為100 MHz，時間步長為0.01 ns，時窗長度為50 ns, 具體參數(shù)如圖5所示。應(yīng)用無單元法分別對均勻介質(zhì)模型(ε1為 8.0)和隨機(jī)介質(zhì)模型(ε1為 6.0、8.0、10.0和15.0)進(jìn)行正演模擬，其模擬所得的GPR正演模擬剖面如圖6和圖7所示。由圖6可見：矩形異常體的上界面為一水平反射界面，矩狀異常體的兩個棱角位置出現(xiàn)繞射波，由于矩形的頂邊很短，導(dǎo)致矩形異常體的上邊在雷達(dá)剖面圖中近似為雙曲線形弧形，非常清晰，矩形異常體的下界面與上界面類似。圖7(a)、(b)、(c)和(d)所示分別為隨機(jī)介質(zhì)模型中矩形異常體的相對介電常數(shù)為6.0、8.0、10.0和15.0時得到的剖面圖。由圖7可知：雷達(dá)剖面圖GPR波散射非常嚴(yán)重，但是，矩形異常體產(chǎn)生的雙曲線反射弧形還是可見；由于GPR波散射現(xiàn)象的出現(xiàn)，使得雙曲線弧形非常不光滑，甚至發(fā)生同相軸不連續(xù)的現(xiàn)象。對比圖7(a)、(b)、(c)和(d)可知，隨著矩形異常體與背景介電常數(shù)差異的增大，雙曲線波形也越來越明顯。但是，對比圖6與圖7(b)可以看到，當(dāng)矩形異常體的相對介電常數(shù)保持不變時，矩形異常體在隨機(jī)介質(zhì)中產(chǎn)生的雙曲線波形比在均勻介質(zhì)中產(chǎn)生的雙曲線波形更微弱。隨機(jī)介質(zhì)模型模擬所得的正演剖面與實測剖面更相符，具有更高的模擬精度，更有利于指導(dǎo)雷達(dá)剖面的數(shù)據(jù)解譯。

4 結(jié)論

1)介紹了基于隨機(jī)過程的譜分解理論和混合型自相關(guān)函數(shù)的隨機(jī)介質(zhì)模型構(gòu)建方法，并詳細(xì)闡述了隨機(jī)介質(zhì)模型的構(gòu)造步驟，詳細(xì)推導(dǎo)了無單元法求解GPR波動方程的具體解法。

2)兩個隨機(jī)介質(zhì)的GPR模型算例結(jié)果表明：與均勻介質(zhì)中GPR波的傳播相比，隨機(jī)介質(zhì)中GPR波的傳播由于受到小尺度異常的影響散射現(xiàn)象非常強(qiáng)烈，并且波形發(fā)生扭曲，隨機(jī)介質(zhì)中異常體產(chǎn)生的反射波形非常不光滑，甚至發(fā)生同相軸不連續(xù)的現(xiàn)象，并且反射波能量較弱。認(rèn)識GPR波在隨機(jī)介質(zhì)中的傳播規(guī)律，更有利于指導(dǎo)雷達(dá)實測剖面的數(shù)據(jù)解譯。

[1]曾昭發(fā), 劉四新.探地雷達(dá)原理與應(yīng)用[M].北京: 電子工業(yè)出版社, 2010: 168-169.ZHENG Zhao-fa, LIU Si-xin.Ground penetrating radar theory and applications [M].Beijing: Electronic Industry Press, 2010:168-169.

[2]戴前偉, 王洪華, 馮德山, 陳德鵬.基于三角形剖分的復(fù)雜GPR模型有限元法正演模擬[J].中南大學(xué)學(xué)報: 自然科學(xué)版,2012, 43(7): 2668-2673.DAI Qian-wei, WANG Hong-hua, FENG De-shan, CHEN De-peng.Finite element method forward simulation for complex geoelectricity GPR model based on triangle mesh dissection [J].Journal of Central South University: Science and Technology,2012, 43(7): 2668-2673.

[3]JAMES I, ROSEMARY K.Numerical modeling of groundpenetrating radar in 2-D using MATLAB [J].Computers &Geosciences, 2006, 32(9): 1247-1258.

[4]黃真萍, 曹杰梅, 周成峰, 馮麗珍.探地雷達(dá)資料的高分辨去噪處理及應(yīng)用[J].福州大學(xué)學(xué)報: 自然科學(xué)版, 2012, 40(4):521-526.HUANG Zhen-ping, CAO Jie-mei, ZHOU Cheng-feng, FENG Li-zhen.The high resolution de-noise processing of groundpenetrating radar data and its engineering application [J].Journal of Fuzhou University: Natural Science Edition, 2012, 40(4):521-526.

[5]奚 先, 姚 姚.二維隨機(jī)介質(zhì)及波動方程正演模擬[J].石油地球物理勘探, 2001, 36(5): 546-552.XI Xian, YAO Yao.2-D random media and wave equation forward modeling [J].Oil Geophysical Prospecting, 2001, 36(5):546-552.

[6]MICHAEL K, HARUO S.Synthesis of plane vector wave envelopes in two dimensional random elastic media based on the Markov approximation and comparison with finite-difference simulations [J].Geophysical Journal International, 2005, 161(3):839-848.

[7]IKELLE L, YUNG S, DAUBE F.2-D random media with ellipsoidal autocorrelation function [J].Geophysics, 1993, 58(4):576-588.

[8]VARADAN V K, MA Y, VARADAN V V.Scattering and attenuation of elastic waves in random media [J].Pageoph,131(4): 577-603.

[9]KORN M.Seismic wave in random media [J].Journal of Applied Geophysics, 1993, 29(3, 4): 247-269.

[10]KNEIB G, KERNER C.Accurate and efficient seismic modeling in random media [J].Geophysics, 1993, 58(4): 576-588.

[11]奚 先, 姚 姚.二維彈性隨機(jī)介質(zhì)中的波場特征[J].地球物理學(xué)進(jìn)展, 2005, 20(1): 147-154.XI Xian, YAO Yao.The wave field characteristics in 2-D elastic random media [J].Progress in Geophysics, 2005, 20(1):147-154.

[12]奚 先, 姚 姚.隨機(jī)介質(zhì)模型的模擬與混合型隨機(jī)介質(zhì)[J].中國地質(zhì)大學(xué)學(xué)報: 地球科學(xué), 2002, 27(1): 67-71.XI Xian, YAO Yao.Simulations of random medium model and intermixed random medium [J].Journal of China University of Geosciences: Earth Science, 2002, 27(1): 67-71.

[13]奚 先, 姚 姚.隨機(jī)介質(zhì)模型正演模擬及其地震波場分析[J].石油物探, 2002, 41(1): 31-36.XI Xian, YAO Yao.Modeling in random medium and its seismic wave field analysis [J].Geophysical Prospecting for Petroleum, 2002, 41(1): 31-36.

[14]郭乃川, 王尚旭, 郭 銳, 啜曉宇.地震勘探中三維小尺度非均勻性隨機(jī)介質(zhì)模型的建立及其特點(diǎn)分析[J].石油天然氣學(xué)報, 2012 34(7): 62-67.GUO Nai-chuan, WANG Shang-xu, GUO Rui, CHUAI Xiao-yu.Construction and feature analysis of three dimensional small scale inhomogeneities in seismic prospecting [J].Journal of Oil and Gas Technology, 2012, 34(7): 62-67.

[15]朱生旺, 魏修成, 曲壽利, 劉春園, 吳開龍.用隨機(jī)介質(zhì)模型方法描述孔洞型油氣儲層[J].地質(zhì)學(xué)報, 2008, 82(3): 420-427.ZHU Sheng-wang, WEI Xiu-cheng, QU Shou-li, LIU Chun-yuan,WU Kai-long.Description of the carbonate karst reservoir with random media model [J].Acta Geologica Sinica, 2008, 82(3):420-427.

[16]陳可洋.三維隨機(jī)建模方法及其波場分析[J].勘探地球物理進(jìn)展, 2009, 32(5): 315-320.CHEN Ke-yang.3-D random modeling scheme and wave field simulation analysis [J].Progress in Exploration Geophysics,2009, 32(5): 315-320.

[17]王金山, 陳可洋, 吳清嶺, 楊 微.隨機(jī)介質(zhì)模型的一種構(gòu)造方法[J].物探與化探, 2010, 34(2): 191-194.WANG Jin-shan, CHEN Ke-yang, WU Qing-ling, YANG Wei.A construction method for random medium model [J].Geophysical & Geochemical Exploration, 2010, 34(2): 191-194.

[18]馮德山, 王洪華, 戴前偉.基于無單元法Galerkin法探地雷達(dá)正演模擬[J].地球物理學(xué)報, 2013, 56(1): 298-308.FENG De-shan, WANG Hong-hua, DAI Qian-wei.Forward simulation of ground penetrating radar based on the element free Galerkin method [J].Chinese Journal of Geophysics, 2013, 56(1):298-308.

[19]殷學(xué)鑫, 劉 洋.二維隨機(jī)介質(zhì)模型正演模擬及其波場分析[J].石油地球物理勘探, 2011, 46(6): 862-872.YIN Xue-xin, LIU Yang.Modeling in 2-D random medium and its wave field analysis [J].Oil Geophysical Prospecting, 2011,46(6): 862-872.

[20]JIANG Zi-meng, ZENG Zhao-fa, LI Jing, LIU Feng-shan, WU Feng-shou.Simulation and analysis of GPR signal based on stochastic media model [C]// Proceedings of the 14th International Conference on Ground Penetrating Radar.Shanghai: IEEE, 2012: 214-219.

[21]底青云, 王妙月.雷達(dá)波有限元仿真模擬[J].地球物理學(xué)報,1999, 42(6): 818-825.DI Qing-yun, WANG Miao-yue.2D finite element modeling for radar wave [J].Chinese Journal of Geophysics, 1999, 42(6):818-825.

[22]賈曉峰, 胡天躍, 王潤秋.無單元法用于地震波波動方程模擬與成像[J].地球物理學(xué)進(jìn)展, 2006, 21(1): 11-17.JIA Xiao-feng, HU Tian-yue, WANG Run-qiu.Wave equation modeling and imaging by using element free method [J].Progress in Geophysics, 2006, 21(1): 11-17.

[23]柳建新, 蔣鵬飛, 童曉忠, 徐凌華, 謝 維, 王 浩.不完全LU分解預(yù)處理的BICGSTAB算法在大地電磁二維正演模擬中的應(yīng)用[J].中南大學(xué)學(xué)報: 自然科學(xué)版, 2009, 40(2):484-491.LIU Jian-xin, JIANG Peng-fei, TONG Xiao-zhong, XU Liang-hua, XIE Wei, WANG Hao.Application of BICGSTAB algorithm with incomplete LU decomposition preconditioning to two-dimensional magnetotelluric forward modeling [J].Journal of Central South University: Science and Technology, 2009,40(2): 484-491.