區(qū)間B樣條小波有限元GPR模擬雙相隨機(jī)混凝土介質(zhì)

馮德山,王珣

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院,長沙　410083 2 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長沙　410083

?

區(qū)間B樣條小波有限元GPR模擬雙相隨機(jī)混凝土介質(zhì)

馮德山1,2,王珣1,2

1 中南大學(xué)地球科學(xué)與信息物理學(xué)院,長沙410083 2 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長沙410083

基于可分離小波理論，由一維區(qū)間B樣條小波尺度函數(shù)的張量積構(gòu)造二維B樣條小波基，并將它作為GPR波動(dòng)方程求解的插值函數(shù)，通過引入轉(zhuǎn)換矩陣，實(shí)現(xiàn)小波系數(shù)空間與雷達(dá)電磁場之間的轉(zhuǎn)換.應(yīng)用Galerkin算法，推導(dǎo)了二維區(qū)間B樣條小波有限元GPR波動(dòng)方程離散格式，求出了2階1尺度與2階2尺度BSWI尺度函數(shù)的積分值及聯(lián)系系數(shù)，給出了該算法的詳細(xì)求解過程.編制了BSWI的Matlab模擬程序，應(yīng)用該程序?qū)蓚€(gè)典型實(shí)例進(jìn)行了正演，結(jié)果表明：BSWI能采用較少的單元達(dá)到與FEM相似的精度，而BSWI算法尺度提升能提高解的精度，但耗時(shí)會(huì)急劇增加.最后，將BSWI算法應(yīng)用于雙相隨機(jī)混凝土模型，說明隨機(jī)介質(zhì)模型理論能靈活、有效地描述實(shí)際混凝土介質(zhì)的分布，正演剖面與實(shí)測剖面特征更相符，能更真實(shí)地模擬雷達(dá)波的傳播過程，可為提高GPR的探測效果和解釋準(zhǔn)確性提供理論基礎(chǔ).

探地雷達(dá)；區(qū)間B樣條；小波有限元；正演模擬；雙相隨機(jī)介質(zhì)

1　引言

混凝土是公路、橋梁、樁基、大壩、房屋等結(jié)構(gòu)工程的主要材料，由硬化的水泥漿、砂、碎石、礫石、卵石等骨料，再加上充水或空氣的孔隙等組成的非均勻多相復(fù)合凝聚體.簡單來看，它是水泥砂裝基體和骨料填充物組成的雙相復(fù)合材料(張劍等，2004).混凝土的質(zhì)量直接關(guān)系到生命及財(cái)產(chǎn)安全，始終是人們關(guān)注的焦點(diǎn).混凝土質(zhì)量無損檢測方法有許多，探地雷達(dá)(GPR)是重要的檢測方法之一.為了提高混凝土構(gòu)筑物的GPR檢測精度，需要開展GPR正演.傳統(tǒng)的GPR正演常將混凝土簡化為均勻介質(zhì)或各向同性介質(zhì)，這顯然不符合混凝土介質(zhì)的特點(diǎn).考慮到混凝土介質(zhì)中存在大量非均質(zhì)的卵石、碎石或砂漿微小異常體，加上混凝土中骨料顆粒的幾何尺寸不規(guī)則及位置隨機(jī)分布，雷達(dá)波在其中傳播時(shí)，必然會(huì)受到不均勻骨料的影響，造成大量小的不相干的擾動(dòng)，易產(chǎn)生波形畸變、能量衰減、傳播方向性差等現(xiàn)象，導(dǎo)致正演結(jié)果與實(shí)際探測結(jié)果在波場特征、幅值強(qiáng)弱、異常體可見性方面存在較大的偏差.需要開展以統(tǒng)計(jì)學(xué)理論為基礎(chǔ)的雙相隨機(jī)介質(zhì)模型來描述，使研究的對(duì)象更符合混凝土的真實(shí)情況(郭士禮，2013).

目前GPR正演算法主要有有限差分法(FDM，劉四新和曾昭發(fā)，2007；李靜等，2010；馮德山等，2010)、有限單元法(FEM，底青云和王妙月，1999；Di and Wang，2004；馮德山等，2012)，這些傳統(tǒng)算法在求解特定奇異性問題時(shí)，難以完全消除局部大梯度所引起的振蕩，影響求解精度.而小波有限元法利用小波函數(shù)或尺度函數(shù)作為插值函數(shù)構(gòu)造單元，使問題的求解在一個(gè)嵌套序列中進(jìn)行，根據(jù)實(shí)際需要任意改變分析尺度，在不改變網(wǎng)格剖分的前提下提高分辨率，使它可以在大梯度處采用小的分析尺度、高階單元以提高分析精度，而在小梯度處采用大的分析尺度、較低階單元以提高分析效率(何正嘉等，2006).在眾多的小波當(dāng)中，Daubechies小波(Daubechies，1998)因?yàn)榫哂芯o支撐性、正交性等諸多優(yōu)點(diǎn)，被廣泛應(yīng)用(Patton and Marks，1996；Chen et al.，2004，2006).Amaratunga等(1994)采用小波Galerkin法結(jié)合Dirichlet邊界條件求解一維Helmholtz方程及二維Green方程(Amaratunga and Williams，1993)，指出了小波嵌套空間能在不同尺度下求解的優(yōu)勢；Sarkar等(1994)將小波函數(shù)引入到傳統(tǒng)有限元插值函數(shù)中求解Maxwell方程，所得的系數(shù)矩陣呈對(duì)角線的稀疏分布，具有條件數(shù)不隨維數(shù)增加的優(yōu)點(diǎn)；Mishra和Sabina(2011)應(yīng)用小波Galerkin法求解一維諧波常微分方程及二維偏微分方程(Sabina and Mishra，2012).然而Daubechies小波沒有顯式表達(dá)，聯(lián)系系數(shù)求解困難，影響了它的進(jìn)一步推廣.而區(qū)間B樣條小波(B-spline wavelet on the interval，BSWI)是美國學(xué)者Chui &Quak(Chui and Quak，1992；Quak and Weyrich，1994)在求解邊值問題時(shí)，為避免定義在整個(gè)實(shí)數(shù)軸上經(jīng)典小波函數(shù)產(chǎn)生數(shù)值振蕩現(xiàn)象而提出的.與其他小波函數(shù)相比，BSWI具有良好的數(shù)值逼近性和明顯數(shù)學(xué)表達(dá)式，易于求導(dǎo)計(jì)算，能根據(jù)問題的求解精度要求靈活地選用小波基的尺度與階次.而且，由于BSWI在空間域具有良好的局部化性質(zhì)，能克服求解邊值問題在邊界上數(shù)值振蕩這一缺陷，可以采用較少單元獲得較高的精度，已引起國內(nèi)外學(xué)者的廣泛關(guān)注(Goswami et al.，1995；Xiang et al.，2007).

石鐘慈(1979)采用三次B樣條有限元解規(guī)則區(qū)域上板梁組合彈性結(jié)構(gòu)的平衡問題，具有比常規(guī)有限元計(jì)算量少、精確度高的特點(diǎn)；梁旭彪等(1988)以矩形單元的B樣條函數(shù)作為形函數(shù)，引入B樣條重節(jié)點(diǎn)理論，成功地將B 樣條有限元法應(yīng)用于多種媒質(zhì)場域內(nèi)電磁場求解問題；金堅(jiān)明等(2006)用張量積形式構(gòu)造二維最小支集樣條小波有限元，并將它用于求解彈性薄板小撓度問題；向家偉等(Xiang et al.，2008a；Xiang et al.，2008b)應(yīng)用構(gòu)造的一維、二維BSWI小波單元對(duì)桿、梁和板等結(jié)構(gòu)件及軸承轉(zhuǎn)子系統(tǒng)進(jìn)行變形和振動(dòng)分析；陳雪峰等(Chen et al.，2010)開展了二維BSWI多尺度小波自適應(yīng)有限元分析；孫惠香等(2014)構(gòu)造了三維BSWI小波單元，并將其應(yīng)用于爆炸荷載作用地下箱形結(jié)構(gòu)的動(dòng)力響應(yīng)的數(shù)值模擬.綜上所述，盡管對(duì)BSWI的研究取得了一些成果，但在地球物理領(lǐng)域仍處于起步階段.

目前開展隨機(jī)介質(zhì)GPR數(shù)值模擬的文獻(xiàn)有：戴前偉和王洪華(2013)基于隨機(jī)過程的譜分解和混合型自相關(guān)函數(shù)理論構(gòu)造隨機(jī)介質(zhì)模型，采用無單元法對(duì)其進(jìn)行了GPR波場正演；李靜(2014)通過在局部隨機(jī)位置引入高斯橢圓方程隨機(jī)干擾因子，結(jié)合混合型自相關(guān)函數(shù)，建立了多參數(shù)耦合隨機(jī)介質(zhì)模型，并應(yīng)用有限差分對(duì)該隨機(jī)介質(zhì)進(jìn)行了正演.本文借助工程材料學(xué)中較為成熟的“數(shù)值混凝土概念”技術(shù)，將瀝青混凝土視作隨機(jī)介質(zhì)，并將介電常數(shù)視為在空間上的隨機(jī)分布特征，使其在物質(zhì)組成，材料參數(shù)，骨料的形態(tài)及分布上與真實(shí)的混凝土在統(tǒng)計(jì)學(xué)意義上具有一致性(高政國和劉光廷，2003；余濤，2013)，建立更符合實(shí)際的雙相離散隨機(jī)混凝土模型，更加真實(shí)地模擬雷達(dá)波的傳播過程.

本文以雙相離散隨機(jī)介質(zhì)來建立“數(shù)值混凝土”模型(余濤，2013；郭士禮，2013)，采用BSWI算法對(duì)其進(jìn)行GPR數(shù)值模擬，能有效地指導(dǎo)混凝土介質(zhì)的GPR探測資料解釋，同時(shí)，為目前GPR正演算法提供一種新的算法和思路.

2　區(qū)間B樣條小波特性及二維張量積小波構(gòu)造

2.1區(qū)間B樣條小波性質(zhì)

通過連接節(jié)點(diǎn)間的分段多項(xiàng)式，構(gòu)成具有一定光滑性質(zhì)的函數(shù)成為B樣條函數(shù)，構(gòu)造[0,1]區(qū)間上m階j尺度B樣條函數(shù)空間，首先將任意區(qū)間[a,b]上函數(shù)f(x)通過坐標(biāo)變換ξ=(x-a)/(b-a)映射到標(biāo)準(zhǔn)[0,1]區(qū)間，再將[0,1]區(qū)間劃分為相等2j份(j∈z，表示尺度)，并在兩端點(diǎn)外側(cè)各添加m-1個(gè)重節(jié)點(diǎn)，節(jié)點(diǎn)總數(shù)為2j+m-1，節(jié)點(diǎn)序列xj用公式表示為

(1)

(2)

m階B樣條尺度函數(shù)表示為

(3)

其支撐區(qū)間為

(4)

(5)

(6)

(7)

2.2二維張量積小波構(gòu)造

(8)

式中，Φ1和Φ2分別為m階j尺度下的一維BSWI尺度函數(shù)，表示為

(9)

(10)

圖1為根據(jù)張量積構(gòu)造的二維張量積BSWI2尺度函數(shù).

圖1　二維張量積BSWI2尺度函數(shù)圖Fig.1　The map of two-dimensional tensor product scaling function of BSWI2

3　GPR波動(dòng)方程的區(qū)間B樣條小波有限元求解

由電磁波傳播理論可知(Yee，1966)，含衰減項(xiàng)的GPR波動(dòng)方程為

(11)

式(11)中，當(dāng)U表示電場值E時(shí)，S為激勵(lì)源Se，當(dāng)U表示磁場值H時(shí)，S為激勵(lì)源Sh，ε為介電常數(shù)，μ為磁導(dǎo)率，σ為電導(dǎo)率，t為時(shí)間，ΓU表示模擬邊界.本文采用滿足Neumann邊界條件的Galerkin原理(徐世浙，1994)來推導(dǎo)GPR小波有限元方程.首先，將待求解區(qū)域離散為若干個(gè)單元，對(duì)任一子單元，可以映射到標(biāo)準(zhǔn)BSWI求解單元.采用二維張量積BSWI尺度函數(shù)為插值函數(shù)構(gòu)造小波單元，求解域被分成n×n網(wǎng)格，其中n=2j+m-2,j>j0，j0為保證至少具有一個(gè)內(nèi)部小波的最小尺度，節(jié)點(diǎn)數(shù)為(n+1)×(n+1).本文采用的V1和V2尺度空間上的二階區(qū)間B樣條(BSWI2)小波單元的節(jié)點(diǎn)排列如圖2所示.

圖2　二維BSWI2小波單元Fig.2　Two-dimensional BSWI2 wavelet element

在小波單元內(nèi)未知場函數(shù)表示為

(12)

Ue為單元內(nèi)的電場值或磁場值，x,y為總體坐標(biāo)系，ζ,η為局部坐標(biāo)系.局部坐標(biāo)與總體坐標(biāo)的轉(zhuǎn)換關(guān)系為：x=a·ζ，y=b·η，a與b分別為1個(gè)單元的長與寬.其中Ae表示單元上待求的小波插值系數(shù)列向量，表示為

(13)

式中，

i=1,2,…,n+1.

定義物理自由度列向量:

(14)

Ue=ReAe,

(15)

式中，矩陣Re為

(16)

式中

(17)

聯(lián)立式(12)和式(15)，有

(18)

令矩陣Re的逆矩陣為轉(zhuǎn)換矩陣Te，即

(19)

則式(18)中形函數(shù)為

(20)

式(18)可表示為

(21)

利用Galerkin法，用形函數(shù)N做權(quán)函數(shù)，將(21)代入(11)式，兩邊同時(shí)乘以δU，在單元內(nèi)積分得

(22)

對(duì)式(22)中左邊第2項(xiàng)采用Green公式變換，得到

(23)

左邊第一項(xiàng)

(24)

得到

Me=∫eNTNdxdy=ab∫eNTNdξdη,

(25)

(26)

左邊第二項(xiàng)

(27)

其中

(28)

左邊第三項(xiàng)

(29)

(30)

右邊項(xiàng)

(31)

(32)

令

(33)

求出相關(guān)的尺度函數(shù)積分值和聯(lián)系系數(shù)值即可求得單元矩陣值，現(xiàn)分別列出1尺度及2尺度相關(guān)值.1尺度：

2尺度：

則

(34)

其中2階1尺度BSWI21的轉(zhuǎn)換矩陣T為9×9單位矩陣，2階2尺度BSWI22的轉(zhuǎn)換矩陣T為25×25單位矩陣.根據(jù)式(24)、(27)、(29)、(31)得到BSWI單元積分

δUTeMeU··e+δUTeKeUe+δUTeK'eU··e-δUTeSE=0.(35) 將各單元的積分相加.相加前,將單元的場值列向量U··e、U·e和Ue擴(kuò)展成全體節(jié)點(diǎn)的列向量U··、U·和U,U··=(U··1,U··2,U··3,…,U··ND)T,U·=(U·1,U·2,U·3,…,U·ND)T,U=(U1,U2,U3,…,UND)T,ND為節(jié)點(diǎn)總數(shù).將單元系數(shù)矩陣Me、Ke和K'e擴(kuò)展成ND×ND的矩陣M、K和K',將列向量SE擴(kuò)展成ND維列向量S. 0=δUTeMeU··e+δUTeKeUe+δUTeK'eU·e-δUTeSE=δUTe(MeU··e+KeUe+K'eU·e-SE). 由于δUT≠0,故有MeU··e+KeUe+K'eU·e=SE. 按照節(jié)點(diǎn)的總體序號(hào),將單元系數(shù)矩陣中的各元位置放置在總體的相應(yīng)行與列的交叉位置上,其余位置的元素為零,因?yàn)閱卧獢U(kuò)展后的U··、U·和U列向量是相同的,所以各單元積分相加時(shí),只要將對(duì)應(yīng)各單元系數(shù)相加即可,即∑NE1MeU··e+∑NE1KeUe+∑NE1K'eU·e=∑NE1SE,故有MU··+KU+K'U·=S.(36) 在求解方程組(36)時(shí),式中的一階及二階導(dǎo)數(shù)可采用中心差分來近似逼近(張新明,2005): U·t=12ΔtUt+Δt-Ut-Δt[],U··t=1(Δt)2Ut+Δt-2Ut+Ut-Δt[],(37)(36)式可化為(M/(Δt)2+K'/(2Δt))Ut+Δt=St+(2M/(Δt)2-K)Ut+(K'/(2Δt)-M/(Δt)2)Ut-Δt,(38)當(dāng)零時(shí)刻或-Δt時(shí)刻,場值為零,且St激勵(lì)源為已知值.因此,可以通過解上面方程逐步求得不同時(shí)間層位上的場值,故式(38)可化簡為Ax=b,(39)其中A=M/(Δt)2+K'/(2Δt),b=St+(2M/(Δt)2-K)Ut+(K'/(2Δt)-M/(Δt)2)Ut-Δt.至此,可以得到含衰減項(xiàng)的GPR波動(dòng)方程的數(shù)值解.由以上求解步驟可知,方程左端矩陣A的病態(tài)程度由M與Δt決定.因此,該GPR波動(dòng)方程的算法精度由M確定,數(shù)值穩(wěn)定性由Δt確定,若不考慮Δt而追求提高BSWI求解精度,可能會(huì)導(dǎo)致求解結(jié)果的不穩(wěn)定.圖3 雷達(dá)模型一示意圖Fig.3 ThesketchmapofGPRmodel14 BSWI算法實(shí)例分析4.1 模型一BSWI與FEM對(duì)比試驗(yàn)為了驗(yàn)證區(qū)間B樣條小波有限元算法的正確性,選取圖3所示的10.0m×10.0m兩層地電模型.上層介質(zhì)厚度為6m,上、下兩層相對(duì)介電常數(shù)分別為3.5與16.0,電導(dǎo)率分別為0.001S·m-1與0.1S·m-1.模型中部有一矩狀空洞異常體,長為2m,高為1m.波源為100MHz脈沖零相位Ricker子波,時(shí)窗長度為100ns,采樣間隔為0.1ns.采用BSWI及FEM兩種算法對(duì)該模型進(jìn)行正演.FEM采用100×100個(gè)邊長0.1m的正方形網(wǎng)格剖分,吸收邊界為20層.小波有限元選用圖1a所示的BSWI21以及圖1b所示的BSWI22尺度函數(shù)為基函數(shù),每個(gè)0.2m×0.2m小網(wǎng)格為1個(gè)小波單

(41)

式中，D0為篩孔直徑，Dmax為最大骨料粒徑，Pk為骨料體積占總體積的百分比.由(41)式可求得模擬區(qū)域內(nèi)截面上相應(yīng)粒徑骨料的顆粒數(shù)，即

(42)

由本文采用橢圓為基本骨料形狀，確定橢圓骨料形狀及其位置需要5個(gè)隨機(jī)參數(shù)：橢圓的長軸，短軸，傾角，形心坐標(biāo)(x0,y0).骨料生成與投放的算法如下所示(余濤，2013)：

(1)根據(jù)級(jí)配曲線計(jì)算各粒徑級(jí)顆粒數(shù)；

(2)依據(jù)粒級(jí)參數(shù)，依次隨機(jī)生成橢圓長軸，短軸，形心坐標(biāo)以及傾角參數(shù)；

(3)判斷是否與已有橢圓骨料干涉，若發(fā)生干涉返回步驟(2)；

(4)輸出該骨料參數(shù)；

(5)判斷該粒級(jí)是否完成，若未完成返回步驟(2)；

(6)判斷是否存在下一粒級(jí)，若存在，改變粒級(jí)參數(shù)并返回步驟(2)；

(7)橢圓骨料生成完畢.

判斷兩個(gè)橢圓骨料是否干涉，先考慮橢圓的定義：橢圓上的點(diǎn)到其兩個(gè)焦點(diǎn)的距離之和等于定長(長軸的長度).兩個(gè)平面上任意兩個(gè)橢圓，如果長軸較短的橢圓上的任意一點(diǎn)到另一橢圓兩個(gè)焦點(diǎn)的距離之和大于長軸較長的橢圓的長軸長度，則這兩個(gè)橢圓不會(huì)發(fā)生干涉.在編程實(shí)踐中，任意點(diǎn)是取橢圓上有限個(gè)點(diǎn)代替，即是用橢圓的內(nèi)接多邊形逼近橢圓.

5.2混凝土模型BSWI模擬算例

圖7所示為含有3個(gè)空洞的模型圖，大小空洞直徑分別為0.2 m和0.08 m，模擬區(qū)域?yàn)?.0 m×1.0 m，背景介質(zhì)都為混凝土，圖7a與圖7b不同之處為：前者將混凝土介質(zhì)當(dāng)成相對(duì)介電常數(shù)為5.45的均質(zhì)介質(zhì)，后者將混凝土看成雙相隨機(jī)介質(zhì)，其中水泥砂漿相對(duì)介電常數(shù)為5.45，骨料相對(duì)介電常數(shù)為8.2，骨料所占模型百分比為0.45，最大和最小骨料粒徑分別為0.05 m和0.008 m，兩者電導(dǎo)率均為0.001 S·m-1.波源為1000 MHz脈沖零相位Ricker子波，采樣時(shí)窗長度為12 ns.采用BSWI21模擬時(shí)，整個(gè)區(qū)域剖分為200×200個(gè)矩形網(wǎng)格，小波單元邊長為0.005 m，吸收邊界為20層，模擬過程中發(fā)射點(diǎn)與接收點(diǎn)同步移動(dòng)，采樣間隔為0.015 ns，共采樣201道雷達(dá)數(shù)據(jù).

圖8a與圖8b分別對(duì)應(yīng)背景為均質(zhì)介質(zhì)與雙相隨機(jī)介質(zhì)的GPR正演剖面圖，其中均質(zhì)介質(zhì)模擬耗時(shí)為8908.47 s.雙相隨機(jī)介質(zhì)正演耗時(shí)為9137.56 s.圖8a可見上部大圓與兩側(cè)小圓引起的雙曲繞射波非常圓滑、清晰，波形振幅強(qiáng)，雙曲線弧頂位置能準(zhǔn)確地指示空洞異常體的上界面.圖8b的雙相隨機(jī)介質(zhì)中，雷達(dá)波受礫石、卵石等骨料的影響散射非常嚴(yán)重，波形發(fā)生扭曲，上部大圓的上界面產(chǎn)生的繞射波弧形欠光滑、波形振幅較弱、畸變嚴(yán)重、同相軸連續(xù)性較差，剖面中出現(xiàn)了較多的干擾雜波.而下界面與兩側(cè)小圓所引起的雙曲線反射弧十分微弱，基本無法識(shí)別.避免了基于均勻介質(zhì)理論可探測與識(shí)別的GPR異常體，但在實(shí)際工作中卻無法探到的偏差，較好地解答了混凝土工程實(shí)測數(shù)據(jù)解釋中遇到的困惑，給實(shí)測資料的合理解譯提供了理論依據(jù).實(shí)例證明，基于雙相隨機(jī)介質(zhì)的混凝土模型所得的正演剖面與實(shí)測剖面更相符，更有利于指導(dǎo)雷達(dá)剖面的數(shù)據(jù)解譯.

圖8　模型三雷達(dá)正演剖面圖Fig.8　The section of model 3 GPR simulation

6　結(jié)論

(1)編制了二維區(qū)間B樣條小波有限元GPR正演程序，通過對(duì)比相同的剖分方式及節(jié)點(diǎn)數(shù)目條件下GPR單道波形圖，BSWI算法能以較少的計(jì)算時(shí)間達(dá)到FEM相似的精度，證明BSWI算法的正確性.而將BSWI的尺度提升后，模擬得到的雷達(dá)波

形更簡潔、能更清晰地體現(xiàn)異常體信息，對(duì)解的精度有較為明顯的提高，但計(jì)算時(shí)間大大增加.說明BSWI是區(qū)別于傳統(tǒng)網(wǎng)格法和多項(xiàng)式法的另一種提高精度的細(xì)化方法，能在不改變網(wǎng)格剖分的前提下提高分辨率，為GPR波動(dòng)方程求解提供新的思路.

(2)將混凝土視為隨機(jī)介質(zhì)，構(gòu)建了雙相離散隨機(jī)混凝土模型，使其在物質(zhì)組成材料參數(shù)、骨料形態(tài)及分布上與真實(shí)的混凝土更具一致性.采用BSWI方法模擬了雙相離散隨機(jī)模型，與均質(zhì)介質(zhì)GPR正演剖面圖相比，隨機(jī)介質(zhì)中GPR波形散射嚴(yán)重，波形扭曲，異常體繞射波弧形欠光滑、波形振幅較弱、剖面中出現(xiàn)了較多的干擾雜波.但隨機(jī)介質(zhì)模型的正演剖面能更真實(shí)地模擬雷達(dá)波的傳播過程，與實(shí)測剖面特征更相符，更有利于指導(dǎo)雷達(dá)剖面的數(shù)據(jù)解譯.

References

Amaratunga K,Williams J R.1993.Wavelet based Green's function approach to 2D PDES.Engineering Computations,10(4):349-367.

Amaratunga K,Williams J R,Qian S,et al.1994.Wavelet-Galerkin solutions for one-dimensional partial differential equations.International Journal for Numerical Methods in Engineering,37(16):2703-2716.

Chen X F,Yang S J,Ma J X,et al.2004.The construction of wavelet finite element and its application.Finite Elements in Analysis and Design,40(5-6):541-554.

Chen X F,He Z J,Xiang J W,et al.2006.A dynamic multiscale lifting computation method using Daubechies wavelet.Journal of Computational and Applied Mathematics,188(2):228-245.

Chen X F,Xiang J W,Li B,et al.2010.A study of multiscale wavelet-based elements for adaptive finite element analysis.Advances in Engineering Software,41(2):196-205.

Chui C K,Quak E.1992.Wavelets on a bounded interval.//Numerical Methods in Approximation Theory,Vol.9.Birkh?user Basel,53-75.Dai Q W,Wang H H.2013.Element free method forward modeling of GPR based on random medium model.The Chinese Journal of Nonferrous Metals (in Chinese),23(9):2436-2443.

Daubechies I.1998.Orthonormal bases of compactly supported wavelets.Communications on Pure and Applied Mathematics,41(7):909-996.

Di Q Y,Wang M Y.1999.2D finite element modeling for radar wave.Chinese J.Geophys.(in Chinese),42(6):818-825.

Di Q Y,Wang M Y.2004.Migration of ground-penetrating radar data with a finite-element method that considers attenuation and dispersion.Geophysics,69(2):472-477.

Feng D S,Chen C S,Dai Q W.2010.GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.Chinese J.Geophys.(in Chinese),53(10):2484-2496,doi:10.3969/j.issn.0001-5733.2010.10.022.

Feng D S,Chen C S,Wang H H.2012.Finite element method GPR forward simulation based on mixed boundary condition.Chinese J Geophys.(in Chinese),55(11):3774-3785,doi:10.6038/j.issn.0001-5733.2012.11.024.

Gao Z G,Liu G T.2003.Two-dimensional random aggregate structure for concrete.Journal of Tsinghua University(Science and Technology)(in Chinese),43(5):710-714.

Goswami J C,Chan A K,Chui C K.1995.On solving first-kind integral equations using wavelets on a bounded interval.IEEE Transactions on Antennas and Propagation,43(6):614-622.

Guo S L.2013.Research on GPR wave theory in random medium for detecting highway [Ph.D.thesis] (in Chinese).Wuhan:China University of Geosciences.

He Z J,Chen X F,Li B,et al.2006.Theory of the Wavelet Based Finite Element Methods and the Application in Engineering (in Chinese).Beijing:Science Press.

Jin J M,Xu Y X,Xue P X.2006.Minimum support spline wavelet finite element.Mathematica Numerica Sinica (in Chinese),28(1):89-112.

Li J,Zeng Z F,Wu F S,et al.2010.Study of three dimension high-order FDTD simulation for GPR.Chinese J.Geophys.( in Chinese),53(4):974-981,doi:10.3969/j.issn.0001-5733.2010.04.022.

Li J.2014.Ground penetrating radar detection and parameter inversion in stochastic effective medium [Ph.D.thesis] (in Chinese).Changchun:Jilin University.

Liang X B,Jian B D,Ni G Z.1988.The B-spline finite element method for electromagnetic field analysis.Journal of Zhejiang University (in Chinese),22(2):100-110.

Liu S X,Zeng Z F.2007.Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.Chinese J.Geophys.(in Chinese),50(1):320-326.

Mishra V,Sabina.2011.Wavelet Galerkin solutions of ordinary differential equations.Int.J.Math.Anal.,5(9-12):407-424.Patton R D,Marks P C.1996.One-dimensional finite elements based on the Daubechies family of wavelets.AIAA Journal,34(8):1696-1698.

Quak E,Weyrich N.1994.Decomposition and reconstruction algorithms for spline wavelets on a bounded interval.Applied and Computational Harmonic Analysis,1(3):217-231.

Sabina,Mishra V.2012.Wavelet-Galerkin solutions of one and two dimensional partial differential equations.Journal of Emerging Trends in Computing and Information Sciences,3(10):1373-1378.

Sarkar T K,Adve R S,García-Castillo L E,et al.1994.Utilization of wavelet concepts in finite elements for an efficient solution of Maxwell's equations.Radio Science,29(4):965-977.

Shih C T.1979.On spline finite element method.Mathematica Numerica Sinica (in Chinese),1(1):50-72.

Sun H X,Xu J Y,Kang T,et al.2014.The three-dimension spline wavelet finite element construction and application on antiknock of underground box structure.Journal of Xi′an University of Architecture &Technology (Natural Science Edition),46(6):816-821,832.

Walraven J C,Reinhardt H W.1981.Theory and experiments on the mechanical behaviour of cracks in plain and reinforced concrete subjected to shear loading.Heron,26(1A):26-35.

Xiang J W,Chen X F,He Z J,et al.2007.The construction of 1D wavelet finite elements for structural analysis.Computational Mechanics,40(2):325-339.

Xiang J W,Chen X F,He Z J,et al.2008a.A new wavelet-based thin plate element using B-spline wavelet on the interval.Computational Mechanics,41(2):243-255.

Xiang J W,Chen X F,Yang L F,et al.2008b.A class of wavelet-based flat shell elements using B-spline wavelet on the interval and its applications.CMES-Computer Modeling in Engineering and Sciences,23(1):1-12.Xu S Z.1994.The Finite Element Method in Geophysics (in Chinese).Beijing:Science Press.

Yee K S.1966.Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEE Transactions on Antennas and Propagation,14(3):302-307.

Yu T.2013.Numerical simulation of ultrasonic wave propagation in concrete [Master thesis] (in Chinese).Changsha:Central South University.

Zhang J,Jin N G,Jin X Y,et al.2004.Numerical simulation method for polygonal aggregate distribution in concrete.Journal of Zhejiang University (Engineering Science)(in Chinese),38(5):581-585.

Zhang X M,Liu K A,Liu J Q.2005.A wavelet finite element method for the 2-D wave equation in fluid-saturated porous media.Chinese J Geophys.(in Chinese),48(5):1156-1166.

Zhong G Q,Li L J,Liu F,et al.2008.Numerical simulation of two-dimensional random aggregate in concrete.Concrete (in Chinese),(9):70-73.

附中文參考文獻(xiàn)

戴前偉,王洪華.2013.基于隨機(jī)介質(zhì)模型的GPR無單元法正演模擬.中國有色金屬學(xué)報(bào),23(9):2436-2443.

底青云,王妙月.1999.雷達(dá)波有限元仿真模擬.地球物理學(xué)報(bào),42(6):818-825.

馮德山,陳承申,戴前偉.2010.基于UPML邊界條件的交替方向隱式有限差分法GPR全波場數(shù)值模擬.地球物理學(xué)報(bào),53(10):2484-2496,doi:10.3969/j.issn.0001-5733.2010.10.022.

馮德山,陳承申,王洪華.2012.基于混合邊界條件的有限單元法GPR正演模擬.地球物理學(xué)報(bào),55(11):3774-3785,doi:10.6038/j.issn.0001-5733.2012.11.024.

高政國,劉光廷.2003.二維混凝土隨機(jī)骨料模型研究.清華大學(xué)學(xué)報(bào)(自然科學(xué)版),43(5):710-714.

郭士禮.2013.基于隨機(jī)介質(zhì)的高速公路探地雷達(dá)檢測理論研究[博士論文].武漢:中國地質(zhì)大學(xué).

何正嘉,陳雪峰,李兵等.2006.小波有限元理論及其工程應(yīng)用.北京:科學(xué)出版社.

金堅(jiān)明,徐應(yīng)祥,薛鵬翔.2006.最小支集樣條小波有限元.計(jì)算數(shù)學(xué),28(1):89-112.

李靜,曾昭發(fā),吳豐收等.2010.探地雷達(dá)三維高階時(shí)域有限差分法模擬研究.地球物理學(xué)報(bào),53(4):974-981,doi:10.3969/j.issn.0001-5733.2010.04.022.

李靜.2014.隨機(jī)等效介質(zhì)探地雷達(dá)探測技術(shù)和參數(shù)反演[博士論文].長春:吉林大學(xué).

梁旭彪,簡柏敦,倪光正.1988.電磁場數(shù)值分析中的B樣條有限元法.浙江大學(xué)學(xué)報(bào),22(2):100-110.

劉四新,曾昭發(fā).2007.頻散介質(zhì)中地質(zhì)雷達(dá)波傳播的數(shù)值模擬.地球物理學(xué)報(bào),50(1):320-326.

石鐘慈.1979.樣條有限元.計(jì)算數(shù)學(xué),1(1):50-72.

孫惠香,許金余,康婷等.2014.三維樣條小波單元構(gòu)造及其在地下箱型結(jié)構(gòu)抗爆數(shù)值模擬中的應(yīng)用.西安建筑科技大學(xué)學(xué)報(bào)(自然科學(xué)版),46(6):816-821,832.

徐世浙.1994.地球物理中的有限單元法.北京:科學(xué)出版社.

余濤.2013.超聲波在混凝土中傳播的數(shù)值模擬[碩士論文].長沙:中南大學(xué).

張劍,金南國,金賢玉等.2004.混凝土多邊形骨料分布的數(shù)值模擬方法.浙江大學(xué)學(xué)報(bào)(工學(xué)版),38(5):581-585.

張新明,劉克安,劉家琦.2005.流體飽和多孔隙介質(zhì)二維彈性波方程正演模擬的小波有限元法.地球物理學(xué)報(bào),48(5):1156-1166.鐘根全,李麗娟,劉鋒等.2008.混凝土二維隨機(jī)骨料的數(shù)值模擬.混凝土,(9):70-73.

(本文編輯胡素芳)

The GPR simulation of bi-phase random concrete medium using finite element of B-spline wavelet on the interval

FENG De-Shan1,2，WANG Xun1,2

1 School of Geosciences and Info-Physics，Central South University，Changsha 410083，China 2 Key Laboratory of Non-ferrous Resources and Geological Detection,Ministry of Hunan Province,Changsha 410083,China

Based on the separable wavelet theory,the scaling functions of one-dimensional B-spline wavelet on the interval (BSWI)is employed to construct the two-dimensional B-spline wavelet bases on the interval.In order to solve the GPR wave equation,these constructed wavelet bases are used as the interpolation functions,and a transformation matrix is introduced to convert between the wavelet coefficient space and the physical space (radar electromagnetic field).In this study,the discrete format of GPR wave equation for two-dimensional interval B-spline wavelet finite element is derived using Galerkin algorithm.And the integral values and the connection coefficients of the second order one-scale and second order two-scale BSWI functions are calculated.The detailed process of the algorithm is given.Then two typical examples are forward modeled using the BSWI method by Matlab program.The result shows that the BSWI with fewer units has the same precision as FEM.Increasing the scale of BSWI algorithm,the precision of the results is also improved,but it is time-consuming.Finally,the BSWI algorithm is applied to model bi-phase random concrete medium,and it proves that the random medium model theory can describe the practical distribution of concrete medium flexibly and effectively.Furthermore,this study found that the forward modeling profile is more accordant with practical profile,which demonstrates that medium model theory can simulate the transmission process of the radar wave more accurately.It also provides a theoretical basis for improving the detection results and interpretation accuracy of GPR.

Ground Penetrating Radar;B-spline wavelet on the interval;Wavelet finite element method;Forward modeling;Bi-phase random medium

馮德山,王珣.2016.區(qū)間B樣條小波有限元GPR模擬雙相隨機(jī)混凝土介質(zhì).地球物理學(xué)報(bào)，59(8):3098-3109,

10.6038/cjg20160832.

Feng D S，Wang X.2016.The GPR simulation of bi-phase random concrete medium using finite element of B-spline wavelet on the interval.Chinese J.Geophys.(in Chinese),59(8):3098-3109,doi:10.6038/cjg20160832.

國家自然科學(xué)基金資助項(xiàng)目(41574116)，中南大學(xué)創(chuàng)新驅(qū)動(dòng)項(xiàng)目(2015CX008)，中南大學(xué)升華育英人才計(jì)劃，2014年度中南大學(xué)教師研究基金(2014JSJJ001)，湖湘青年創(chuàng)新創(chuàng)業(yè)平臺(tái)培養(yǎng)對(duì)象項(xiàng)目，中南大學(xué)研究生自主探索創(chuàng)新項(xiàng)目(2015ZZTS249)共同資助.

馮德山，男，博士，教授，從事地球物理數(shù)據(jù)處理與正反演研究.E-mail:fengdeshan@126.com

10.6038/cjg20160832

P631

2015-04-14，2016-06-14收修定稿