基于隨機介質模型的GPR無單元法正演模擬

戴前偉，王洪華

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

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

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

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

1 隨機介質模型

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

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

假設空間隨機擾動φ具有零均值、一定方差及自相關函數的空間平穩隨機過程，則

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

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

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

3)用隨機數發生器生成0, 2π[]φ上服從均勻分布的獨立的二維隨機場

4)應用下式計算隨機功率譜：

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

8)將式(10)代入式(1)中，即可得到隨機介質模型的相對介電常數。隨機介質模型的電導率也可按上述方法構建。

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

2 GPR無單元法正演模擬

MAXWELL方程組描述了電磁場的運動學和動力學規律。由電磁波理論，高頻電磁波(GPR)在介質中的傳播規律也應服從MAXWELL方程組。以電場為例，根據文獻[21]的推導，GPR波滿足的波動方程為

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

首先定義如下形式的近似解函數：

圖1 采用混合型橢圓自相關函數產生的不同相關長度的隨機介質相對介電常數模型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)加權求和得：

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

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

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

式中：

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

式中：Ni(x)為節點i的形函數，它是坐標的函數。

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

式中：S為等效的電場源向量，為電場對時間的二次導數項，為電場對時間的一次導數項，M為質量矩陣，K’為阻尼矩陣，K為剛度矩陣。

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

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

3 數值模擬實例

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

3.1 模型1

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

圖2 隨機介質模型示意圖Fig.2 Sketch map of random medium model

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

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

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

4 結論

1)介紹了基于隨機過程的譜分解理論和混合型自相關函數的隨機介質模型構建方法，并詳細闡述了隨機介質模型的構造步驟，詳細推導了無單元法求解GPR波動方程的具體解法。

2)兩個隨機介質的GPR模型算例結果表明：與均勻介質中GPR波的傳播相比，隨機介質中GPR波的傳播由于受到小尺度異常的影響散射現象非常強烈，并且波形發生扭曲，隨機介質中異常體產生的反射波形非常不光滑，甚至發生同相軸不連續的現象，并且反射波能量較弱。認識GPR波在隨機介質中的傳播規律，更有利于指導雷達實測剖面的數據解譯。

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

[2]戴前偉, 王洪華, 馮德山, 陳德鵬.基于三角形剖分的復雜GPR模型有限元法正演模擬[J].中南大學學報: 自然科學版,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]黃真萍, 曹杰梅, 周成峰, 馮麗珍.探地雷達資料的高分辨去噪處理及應用[J].福州大學學報: 自然科學版, 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].石油地球物理勘探, 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].地球物理學進展, 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].中國地質大學學報: 地球科學, 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].石油物探, 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].石油天然氣學報, 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].地質學報, 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].勘探地球物理進展, 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].物探與化探, 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法探地雷達正演模擬[J].地球物理學報, 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]殷學鑫, 劉 洋.二維隨機介質模型正演模擬及其波場分析[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]底青云, 王妙月.雷達波有限元仿真模擬[J].地球物理學報,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].地球物理學進展, 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分解預處理的BICGSTAB算法在大地電磁二維正演模擬中的應用[J].中南大學學報: 自然科學版, 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.