Daubechies小波有限元求解GPR波動方程

馮德山, 楊炳坤, 王珣, 杜華坤

1 中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 長沙　410083 2 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室, 長沙　410083 3 福建省建筑科學(xué)研究院， 福州　350025

?

Daubechies小波有限元求解GPR波動方程

馮德山1,2, 楊炳坤1,3, 王珣1,2, 杜華坤1,2

1 中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 長沙410083 2 中南大學(xué)有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點(diǎn)實(shí)驗(yàn)室, 長沙410083 3 福建省建筑科學(xué)研究院， 福州350025

摘要基于可分離小波理論，由一維Daubechies尺度函數(shù)的張量積構(gòu)造二維Daubechies小波基，并將它作為GPR波動方程求解的插值函數(shù)，導(dǎo)出了二維Daubechies小波有限元GPR方程離散格式；通過引入轉(zhuǎn)換矩陣，實(shí)現(xiàn)小波系數(shù)空間與雷達(dá)場值之間轉(zhuǎn)換.引入自由度凝聚技術(shù)，有效解決了小波有限元求解中小波單元內(nèi)部自由度過多的問題，節(jié)約了計(jì)算量并方便與傳統(tǒng)有限元法耦合.然后，詳細(xì)闡述了Daubechies小波有限元聯(lián)系系數(shù)計(jì)算方法，有效解決了小波有限元求解偏微分方程的難點(diǎn)與核心問題.最后，以兩個典型GPR模型為例，對比了Daubechies小波有限元與傳統(tǒng)有限元的雷達(dá)正演剖面圖與單道波形圖，結(jié)果表明：在相同的剖分方式及節(jié)點(diǎn)數(shù)目條件下，Daubechies小波有限元的緊支性與正交性一定程度上提高了求解效率，它與有限元法求解結(jié)果能較好地吻合，驗(yàn)證了Daubechies小波有限元算法的正確性.

關(guān)鍵詞探地雷達(dá)； Daubechies小波有限元； 自由度凝聚技術(shù)； 聯(lián)系系數(shù)； 波動方程； 正演模擬

1引言

GPR正演傳統(tǒng)算法主要有：有限差分法(劉四新和曾昭發(fā)，2007；李靜等，2010；馮德山等，2010)、有限元法(底青云和王妙月，2010；馮德山等，2012)，它們理論體系日趨成熟完善，但都是選用多項(xiàng)式作為基函數(shù)，是全局化的函數(shù)，當(dāng)開展簡單的模型計(jì)算時具有優(yōu)勢.然而，隨著GPR工程勘探日益復(fù)雜化、精細(xì)化，如仍以低階多項(xiàng)式的基函數(shù)去逼近場函數(shù)，難以完全消除局部大梯度問題所引起的振蕩，影響求解精度.以有限單元法(FEM)為例，它在求解斷層斷點(diǎn)、地質(zhì)裂隙等局部大梯度、奇異性問題時，計(jì)算誤差取決于對模型的一次性剖分，當(dāng)計(jì)算結(jié)果存在較大誤差時，需要加密網(wǎng)格或提高插值多項(xiàng)式的階次，原來形成的剛度矩陣不能夠在網(wǎng)格細(xì)化后的計(jì)算中被繼承，勢必增加大量計(jì)算成本，浪費(fèi)計(jì)算資源，所以尋找更優(yōu)的插值函數(shù)也成為當(dāng)前有限元發(fā)展的新研究方向之一.

小波分析是調(diào)和分析發(fā)展史上里程碑式的進(jìn)展(程正興，2006)，它比傳統(tǒng)Fourier分析能更好地處理局部奇異性問題，具有多分辨特性(Qian and Weiss，1993；Ko et al，1995).近年來，小波分析與有限元方法相結(jié)合，產(chǎn)生了小波有限元，它作為一種偏微分方程(PDE)的潛在高效求解方法被提出，能將分析對象依次放入一個逐級擴(kuò)大、互相嵌套的函數(shù)空間序列…，V-1,V0,V1…中進(jìn)行分析，而且作為空間Vj+1的子空間Vj，存在相應(yīng)的補(bǔ)空間Wj，當(dāng)需要提高求解精度時，通過增加Wj空間以擴(kuò)大單元容許空間至Vj+1.因此，小波有限元能根據(jù)實(shí)際需要任意改變分析尺度，在不改變網(wǎng)格剖分的前提下提高分辨率，使其可以在大梯度處采用小的分析尺度、高階單元以提高分析精度，而在小梯度處采用大的分析尺度、較低階單元以提高分析效率.在眾多的小波當(dāng)中，Daubechies小波(Daubechies，1988)因?yàn)榫哂芯o支撐性、正交性等諸多優(yōu)點(diǎn)，已引起國內(nèi)外學(xué)者的廣泛關(guān)注.Daubechies小波的緊支性保證了在沒有截?cái)嗾`差的條件下，能以最小的單元自由度最大限度地逼近待求函數(shù)，通過閾值運(yùn)算，達(dá)到待求系數(shù)最少；而正交性使得小波有限元的剛度矩陣是帶狀稀疏矩陣，從而可以減少代數(shù)方程組的求解運(yùn)算量；小波函數(shù)的消失矩特性指明Daubechies小波函數(shù)φN(x)可以精確地表征出不大于N-1階的冪級數(shù)(Daubechies，1992).

Amaratunga等(1994)采用小波Galerkin法結(jié)合Dirichlet邊界條件求解一維Helmholtz方程及二維Green方程(Amaratunga and Williams，1993)，指出了小波嵌套空間能在不同尺度下求解的優(yōu)勢；Sarkar等(1994)將小波函數(shù)引入到傳統(tǒng)有限元插值函數(shù)中求解Maxwell方程，所得的系數(shù)矩陣呈對角線的稀疏分布，具有條件數(shù)不隨維數(shù)增加的優(yōu)點(diǎn)；Patton和Marks(1996)構(gòu)造了一維Daubechies小波單元；西安交通大學(xué)機(jī)械自動化研究所在何正嘉等(2006)的帶領(lǐng)下，涌現(xiàn)出一大批小波有限元數(shù)值模擬(Chen et al.，2004，2006)、裂紋織別(Dong et al.，2009)等應(yīng)用研究優(yōu)秀成果；楊仕友與倪光正(1999)將無窮區(qū)間Daubechies小波有限元法應(yīng)用到了電磁場數(shù)值計(jì)算中；石陸魁等(2001)采用小波-Galerkin法求解了二維多介質(zhì)靜態(tài)電磁場邊值問題；張新明等(2005)把小波有限元法引入到流體飽和多孔隙介質(zhì)二維波動方程的正演模擬中；Chen等(2006)應(yīng)用Daubechies小波開展了動態(tài)多尺度提升計(jì)算，并對小波有限元聯(lián)系系數(shù)、剛度矩陣與載荷列陣的計(jì)算進(jìn)行了詳細(xì)的探討；陳雅琴(2008)在其博士論文中提出一種基于廣義變分原理的Daubechies條件小波法，使小波Galerkin法和小波Ritz法的求解精度得到了一定的提高；Mishra和Sabina(2011)應(yīng)用小波Galerkin法求解一維諧波常微分方程及二維偏微分方程(Sabina and Mishra，2012)；Suk-In和Schulz (2013)應(yīng)用小波Galerkin法求解非線性黏稠Burgers偏微分方程.綜上所述，Daubechies小波有限元的研究取得了一些成果，但在地球物理領(lǐng)域仍處于起步階段，有待進(jìn)一步探索與發(fā)展.本文應(yīng)用Daubechies小波有限元開展GPR正演，能對目前GPR算法形成有益補(bǔ)充.

2Daubechies小波特性及二維張量積小波構(gòu)造

2.1Daubechies小波性質(zhì)

Daubechies小波自問世以來，引起了眾多學(xué)者的廣泛關(guān)注，其理論和應(yīng)用研究異?；钴S，主要性質(zhì)有(Daubechies，1992)：

(1) 緊支性：緊支性反映尺度函數(shù)與小波函數(shù)在時域的局部化能力，支撐區(qū)間越小，局部化能力越強(qiáng).對給定的正整數(shù)N≥2，尺度函數(shù)φN(x)和小波函數(shù)ψN(x)的支撐域分別為：

(1)

(2) 消失矩特性：Daubechies小波的光滑性可以用消失矩來刻畫，消失矩階數(shù)越大則小波函數(shù)越光滑，追求好的光滑性，必然以擴(kuò)大支撐區(qū)間，降低時域局部化能力為代價的.Daubechies小波的消失矩為N，即

(2)

(3) 正交性：Daubechies小波是正交小波，其尺度函數(shù)φN(x)與小波函數(shù)ψN(x)滿足如下正交條件：

(3)

(4)

式(3)中

(5)

(4) 歸一性：在Daubechies小波的構(gòu)造中，常要求尺度函數(shù)φN(x)滿足歸一性條件

(6)

(5) 兩尺度方程：由于Daubechies小波沒有明確的數(shù)學(xué)表達(dá)式，Daubechies小波的尺度函數(shù)φN(x)和小波函數(shù)ψN(x)由兩尺度方程構(gòu)造而來：

(7)

其中pN(k)稱為小波逼近空間低通濾波系數(shù)，qN(k)稱為小波空間帶通濾波系數(shù).

(8)

對于給定的正整數(shù)N，Daubechies小波僅有2N個pN(k)不等于零.為了便于敘述，將具有2N個非零濾波系數(shù)的Daubechies小波簡稱為DBN小波.

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

(9)

φ(x,y)=φ1(x)φ2(y)=φj,k,l(x,y)

=2jφjk(2jx-k)φjl(2jy-l),

(10)

式(10)中，-(2jL+2N-1)≤k,l≤0.L為x和y的定義域[0,L].則{φj,k,l(x,y)}是Vj的基底，所以{Vj}形成L2(R2)中的一個多分辨分析，而φ(x,y)是相應(yīng)的尺度函數(shù).圖1為根據(jù)張量積構(gòu)造的二維DB3小波尺度函數(shù)圖.

圖1　DB3小波二維張量積尺度函數(shù)圖

3GPR波動方程的Daubechies小波有限元求解

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

(11)

×φ(η-l),

(12)

式(12)中ak,l(t)為待求的小波系數(shù)，它僅與時間有關(guān)；φ為Daubechies小波尺度函數(shù)，僅與空間有關(guān)，Ue為單元內(nèi)的電場值或磁場值；x,y為總體坐標(biāo)系，ξ,η為局部坐標(biāo)系；以矩陣的形式表示為：

(13)

其中，Φ=[Φ(ξ)?Φ(η)]，Ae為待求小波系數(shù)組成的列向量，其SuppφN(x)=[0,2N-1].以V0空間上的DB3小波為例，1個二維小波單元自由度總數(shù)為：(2N-1)×(2N-1)，即25個單元自由度，則構(gòu)造的具有25個節(jié)點(diǎn)的小波單元如圖2所示.

圖2構(gòu)造的小波單元的局部坐標(biāo)取值范圍是[-1,1]，結(jié)合Daubechies小波尺度函數(shù)的支撐區(qū)間SuppφN(x)=[0,2N-1]，可以求得k,l的取值范圍為：

圖2　25個自由度的零尺度二維DB3小波單元

2-2N≤k,l≤0,

(14)

則局部坐標(biāo)與總體坐標(biāo)的轉(zhuǎn)換關(guān)系如下：

(15)

其中，x0、y0為總體坐標(biāo)下子單元的中心坐標(biāo)，a、b是總體坐標(biāo)下子單元的邊長.-1≤ξ,η≤1，則有如下微分關(guān)系：

(16)

則由總體坐標(biāo)與局部坐標(biāo)關(guān)系可知，若設(shè)在總體坐標(biāo)中滿足：

(17)

利用Galerkin法，將(17)式代入(11)式，兩邊同時乘以二維Daubechies小波基函數(shù)：

(18)

在單元內(nèi)積分得：

(19)

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

=?ΩSφk2,l2(ξ,η)dΩ,

(20)

再將(20)式中的U用(17)式代入，可得到：

= ?ΩSφ(ξ-k2)φ(η-l2)dΩ,

(21)

進(jìn)一步展開得：

= ?ΩSφ(ξ-k2)φ(η-l2)dΩ,

(22)

?ΩS·φk2,l2(ξ,η)dΩ=F1.

(23)

則(22)式可化為：

(24)

將系數(shù)ak,l(t)按照一定的順序重新組合形成一維向量A，即

A=(a2-2N,2-2N,a2-2N,3-2N,…,a2-2N,0,a3-2N,2-2N,…,a3-2N,0,…,a0,1-2N,a0,2-2N,…,a0,0)T.

(25)

再作如下假定：

(26)

(27)

則(24)式可改寫為：

(28)

(29)

在求解方程組(29)時，式中的一階及二階導(dǎo)數(shù)可采用中心差分來近似逼近：

(30)

(29)式可化為：

(31)

當(dāng)零時刻或-Δt時刻，電場值為零，故ak,l(0)=0，ak,l(-Δt)=0，且激勵源F為已知值.因此，可以通過解上面方程逐步求得不同時間層位上的小波系數(shù)，故式(31)可化簡為

(32)

(33)

由于小波有限元單元平衡方程建立于小波空間，此時的單元自由度是小波系數(shù)而非電場值，為了實(shí)現(xiàn)小波系數(shù)空間與GPR波動方程電場(或磁場)空間之間的變換，可以構(gòu)造快速小波變換(張新明等，2005；何正嘉等，2006)：

(34)

E為單元電場向量，Φ為小波系數(shù)變換矩陣，A為小波系數(shù)向量.仍以V0尺度的25節(jié)點(diǎn)DB3小波單元為例，以節(jié)點(diǎn)電場作為單元自由度，則

(35)

(36)

(37)

式(37)中的Φ元素為二維DB3尺度函數(shù)整數(shù)點(diǎn)及二分點(diǎn)上的函數(shù)值.

4自由度凝聚技術(shù)

采用DBN小波尺度函數(shù)構(gòu)造的二維小波單元時，每個單元有(2N-1)×(2N-1)個自由度，相對傳統(tǒng)FEM法增加了很多單元內(nèi)部自由度，單元矩陣將變得比較龐大，給單元離散和單元組合帶來一定困難.因此，有必要引入自由度凝聚技術(shù)將多余的單元自由度消去，剩下需要的自由度(MehraeenandChen，2006).假如記矩陣Φ的逆矩陣為P，則由(34)式可得：

(38)

將(38)式代入(32)式，并在等式兩邊同乘以PT，可得

(39)

(40)

將(40)式的矩陣重新排列并進(jìn)行分塊，得

(41)

由式(41)第2個等式可知

(42)

將(42)式再代回式(41)第1個等式，就消去了除角點(diǎn)之外的自由度，得

(43)

化簡(43)式可得

(44)

5Daubechies小波聯(lián)系系數(shù)計(jì)算

圖3　DB3小波尺度函數(shù)一階導(dǎo)數(shù)圖

圖4　DB3小波尺度函數(shù)一階導(dǎo)數(shù)內(nèi)積圖

下面著重介紹式(23)中的Daubechies小波聯(lián)系系數(shù)的計(jì)算.由于式(23)中的積分區(qū)間為[-1,1]，而大部分聯(lián)系系數(shù)求解方法采用的區(qū)間為[0,1]，為此可將任意長度為L的區(qū)域投射到[0,1]的區(qū)間，譬如：

(45)

本文中僅介紹x的積分區(qū)間為[0,1]，結(jié)合Daubechies小波尺度函數(shù)的支撐區(qū)間，可求得k1、k2的取值范圍為：

(46)

(47)

令ξ=2x，則

(48)

根據(jù)式(7)的尺度函數(shù)表達(dá)式，可得

(49)

式(49)中s取值范圍與k1,k2一致.將式(49)兩邊同求d次導(dǎo)，則

(50)

本文采用Vj尺度下聯(lián)系系數(shù)通用式求解過程作為范例，V0尺度聯(lián)系系數(shù)可類似得到.

根據(jù)式(50)可得

(51)

式(51)中，s,t的取值范圍與k1,k2一致，且必須保證2-2N≤s-2k1，t-2k2，s-2k1+2j，t-2k2+2j≤2j-1.將式(51)寫成矩陣形式，可得如下方程組：

(52)

式(52)中：

(53)

(54)

等式(54)左邊矩陣為奇異陣，不能唯一確定聯(lián)系系數(shù)的值.目前常用的求解方法是：首先由求解矩陣的性質(zhì)得出聯(lián)系系數(shù)的解空間，引入與解空間自由度個數(shù)相等的附加方程，從而組成恰定方程組，求解后可唯一確定準(zhǔn)確的聯(lián)系系數(shù)值.對于添加的附加方程，Latto等(1999)、Yang等(2000)、陳雪峰等(2006)給出如下的計(jì)算式：

(55)

(56a)

(56b)

求解公式(56a)秩為24，待求聯(lián)系系數(shù)總數(shù)為25，僅需添加式(55)中1個方程式.求解公式(56b)秩為22，待求聯(lián)系系數(shù)總數(shù)為25，需添加式(55)中3個方程式.2個聯(lián)系系數(shù)計(jì)算結(jié)果分別參見表1與表2，利用該數(shù)據(jù)繪出的聯(lián)系系數(shù)如圖5a與圖5b所示.其中表1中的計(jì)算結(jié)果若保留8位有效位數(shù)字，與Landragin-Frassati(2009)文章中計(jì)算結(jié)果完全一致，證明程序的正確性.

表1　DB3小波聯(lián)系系數(shù)k1,k2值

表2　DB3小波聯(lián)系系數(shù)k1,k2值

圖5　DB3小波尺度函數(shù)聯(lián)系系數(shù)計(jì)算結(jié)果圖(a) 聯(lián)系系數(shù)的計(jì)算結(jié)果； (b) 聯(lián)系系數(shù)的計(jì)算結(jié)果.

6數(shù)值算例

6.1矩狀模型

選取圖6所示矩狀地電模型，模擬區(qū)域?yàn)?.0 m×2.0 m，上層混凝土介質(zhì)的相對介電常數(shù)為8.0，電導(dǎo)率為0.0001 S·m-1，厚度為0.8 m，下層介質(zhì)介電常數(shù)為15.0，電導(dǎo)率為0.01 S·m-1.上層介質(zhì)中埋有矩狀異常體，其介電常數(shù)為3.0，電導(dǎo)率為0.002 S·m-1，異常體長0.4 m，寬為0.2 m，異常體頂邊距地面0.4 m，波源為脈沖Ricker子波，中心頻率取900 MHz，采樣時窗長度為30 ns，采樣時間間隔為0.02 ns.

基于FEM算法對該模型進(jìn)行正演，整個區(qū)域由單元邊長0.025 m的正方形剖分為160×80的網(wǎng)格空間.采用DB3小波有限元進(jìn)行正演，共計(jì)40×20個小波單元，每個小波單元網(wǎng)格被DB3小波插入3×3個內(nèi)部節(jié)點(diǎn)，每個小波單元細(xì)分為4×4區(qū)間，則最終的小波區(qū)間由邊長0.025 m的正方形網(wǎng)格組成，剖分總網(wǎng)格與正方形FEM算法相當(dāng)為160×80個網(wǎng)格.

圖6　矩狀雷達(dá)地電模型示意圖

圖7a為160道數(shù)據(jù)的FEM正演剖面圖，圖中可見，8 ns位置為矩狀異常體的上界面，由于繞射波的影響，異常體的兩個棱角處出現(xiàn)了繞射波，10 ns左右可見異常體下界面反射.圖7b為DB3小波有限元法正演剖面圖，觀察圖7a與圖7b，兩圖非常類似，都能反映異常體的形態(tài)，但是小波有限元多次波波形更為離散一些，推斷與小波有限元聯(lián)系系數(shù)的計(jì)算精度有關(guān)聯(lián).

圖8為兩種方法第80道雷達(dá)波形對比圖，圖8a為0～24 ns時的顯示，說明兩種算法能很好地?cái)M合.圖8b為了突出波場細(xì)微的差別，僅顯示12～24 ns時的單道雷達(dá)波形，圖中可見，小波有限元與FEM計(jì)算結(jié)果之間仍存在細(xì)小的差別，但總體擬合趨勢很好，說明了小波有限元計(jì)算結(jié)果的可靠性.

6.2組合模型

為了說明小波有限元對復(fù)雜模型模擬的有效性，選取圖9所示的組合雷達(dá)地電模型驗(yàn)證.模擬區(qū)域、背景介質(zhì)、網(wǎng)格剖分方式及雷達(dá)頻率、采樣參數(shù)與上例一致.在上層介質(zhì)中左邊為矩狀空洞異常體，其長與寬都為0.2 m；中間為菱形素填土介質(zhì)，其介電常數(shù)為3.0，電導(dǎo)率為0.002 S·m-1，菱形長0.6 m，高為0.3 m，其上頂點(diǎn)距模擬區(qū)域上邊界0.3 m；右邊為半徑為0.05的金屬球體.在Intel(R) CoreTMi7-4500U CPU@1.80 GHz，8.00 GB的內(nèi)存物理地址擴(kuò)展，Window 8操作系統(tǒng)IBM X240s筆記本上計(jì)算該模型雷達(dá)剖面160道數(shù)據(jù)，采用DB3小波有限元法計(jì)算時間為6903.25 s，而采用傳統(tǒng)有限元計(jì)算時間為7845.84 s.可見，盡管小波有限元增加了小波有限元聯(lián)系系數(shù)與凝聚技術(shù)的轉(zhuǎn)換運(yùn)算，但是由于這些計(jì)算能事先計(jì)算好再調(diào)用，并不增加太多計(jì)算工作量，相反小波有限元具有的正交性，使得求解的剛度矩陣為極度稀疏矩陣，可提高計(jì)算效率.

圖7　矩狀模型雷達(dá)正演剖面圖(a) FEM方法； (b) 小波有限元.

圖8　矩狀模型第80道波形數(shù)據(jù)對比(a) 0～24 ns；(b) 12～24 ns.

圖9　組合雷達(dá)模型示意圖

圖10a與圖10b分別為應(yīng)用FEM及小波有限元正演所得的160道波形的正演剖面圖，從圖中可見，橫坐標(biāo)2.0 m，縱坐標(biāo)6 ns處為菱形的上界面，菱形的上面兩條邊也較好地對應(yīng)了圖中的反射波；左邊橫坐標(biāo)1.0 m，縱坐標(biāo)9.5 ns及11 ns處兩處繞射波較好地對應(yīng)了空洞異常體上下界面；而右邊橫坐標(biāo)3.0 m縱坐標(biāo)10 ns處為金屬球體上界面繞射波，由于金屬球體的全反射，導(dǎo)致其下界面基本見不到繞射波.對比圖10a與圖10b，兩圖波形大致類似，僅多次波的體現(xiàn)上存在細(xì)微差別，小波有限元波形更為豐富，在增益倍數(shù)均相同條件下，小波有限元幅值稍強(qiáng)于有限元法.

圖11a為0～24 ns兩種方法第40道波形對比， 黑線表示的FEM計(jì)算結(jié)果與紅色點(diǎn)表示的小波有限元計(jì)算結(jié)果能較好地吻合.圖11b為僅顯示4～24 ns時的單道雷達(dá)波形，由于去掉了1.5 ns左右直達(dá)波的大振幅值，波場細(xì)節(jié)得已體現(xiàn)，圖中可見，在10.0 ns處的異常體波形幅值稍大.

圖10　組合模型雷達(dá)正演剖面圖(a) FEM方法；(b) 小波有限元.

圖11　組合模型的第40道雷達(dá)波形對比(a) 0～24 ns；(b) 4～24 ns.

圖12a為第120道0～24 ns的雷達(dá)波形數(shù)據(jù)，兩種方法計(jì)算結(jié)果同樣能較好地?cái)M合，10.0 ns處為圓形球體的上界面反射.圖12b為僅顯示4～24 ns時的單道雷達(dá)波形，圖中可見，小波有限元與FEM計(jì)算結(jié)果僅在16.0～20.0 ns之間的多次反射波及界面反射之間存在細(xì)微的差別，但兩者總體擬合很好，說明了小波有限元計(jì)算結(jié)果的可靠性.

圖12　組合模型的第120道雷達(dá)波形數(shù)據(jù)對比(a) 0～24 ns；(b) 4～24 ns.

7結(jié)論

(1) 詳細(xì)闡述了Daubechies小波有限元聯(lián)系系數(shù)計(jì)算方法，有效解決了小波有限元求解偏微分方程的難點(diǎn)與核心問題.通過引入自由度凝聚技術(shù)消除單元內(nèi)部的自由度，剩下四個角點(diǎn)的自由度，解決了小波單元自由度過多的問題，該技術(shù)能在小波有限元與有限元耦合計(jì)算中得到廣泛應(yīng)用.

(2) 編制了二維張量積Daubechies小波有限元GPR正演程序，通過對比相同的剖分方式及節(jié)點(diǎn)數(shù)目條件下GPR正演剖面圖與單道波形圖，表明小波有限元與傳統(tǒng)有限元模擬結(jié)果能較好地吻合，證明小波有限元算法的正確性.由于Daubechies小波的正交性、緊支性使得Daubechies小波有限元剛度矩陣為條帶狀稀疏矩陣，提高了求解效率，為GPR波動方程求解提供了一種新的思路.

References

Amaratunga K, Williams J R. 1993. Wavelet based Green′s function approach to 2D PDES.EngineeringComputations, 10(4): 349-367.Amaratunga K, Williams J R, Qian S, et al. 1994. Wavelet-Galerkin solutions for one-dimensional partial differential equations.InternationalJournalforNumericalMethodsinEngineering, 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.FiniteElementsinAnalysisandDesign, 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.JournalofComputationalandAppliedMathematics, 188(2): 228-245. Chen Y Q. 2008. Study on Generalized Variational Principle in bridge structural analysis-Daubechies conditional wavelet [Ph. D. thesis] (in Chinese). Xi′an: Chang′an University.

Cheng Z X. 2006. Wavelet Analysis and its Application Examples (in Chinese). Xi′an: Xi′an Jiaotong University Press. Daubechies I. 1992. Ten Lectures on Wavelets (CBMS-NSF Regional Conference Series in Applied Mathematics, 61. Philadelphia: Society for Industrial and Applied Mathematics.Daubechies I. 1988. Orthonormal bases of compactly supported wavelets.CommunicationsonPureandAppliedMathematics, 41(7): 909-996.

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

Dong H B, Chen X F, Li B, et al. 2009. Rotor crack detection based on high-precision modal parameter identification method and wavelet finite element model.MechanicalSystemsandSignalProcessing, 23(3): 869-883.

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.ChineseJ.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.ChineseJ.Geophys. (in Chinese), 55(11): 3774-3785, doi: 10.6038/j.issn.0001-5733.2012.11.024.

Gao Y L. 2009. Method of numerical integration with Daubechies wavelet.JournalofJiangnanUniversity(NaturalScienceEdition) (in Chinese), 8(1): 122-125.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.

Ko J, Kurdila A J, Pilant M S. 1995. A class of finite element methods based on orthonormal, compactly supported wavelets.ComputationalMechanics, 16(4): 235-244. Landragin-Frassati A, Bonnet S, Silva A D, et al. 2009. Application of a wavelet-Galerkin method to the forward problem resolution in fluorescence Diffuse Optical Tomography.OpticsExpress, 17(21): 18433-18448.Latto A, Resnikoff L, Tenenbaum E. 1999. The evaluation of connection coefficients of compactly supported wavelets. Aware Inc., Technical Report AD910708.Li J, Zeng Z F, Wu F S, et al. 2010. Study of three dimension high-order FDTD simulation for GPR.ChineseJ.Geophys. (in Chinese), 53(4): 974-981, doi: 10.3969/j.issn.0001-5733.2010.04.022.Liu S X, Zeng Z F. 2007. Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.ChineseJ.Geophys. (in Chinese), 50(1): 320-326.Mehraeen S, Chen J S. 2006. Wavelet Galerkin method in multi-scale homogenization of heterogeneous media.InternationalJournalforNumericalMethodsinEngineering, 66(3): 381-403.

Mishra V, Sabina. 2011. Wavelet Galerkin solutions of ordinary differential equations.Int.JournalofMath.Analysis, 5(9): 407-424.

Patton R D, Marks P C. 1996. One-dimensional finite elements based on the Daubechies family of wavelets.AIAAJournal, 34(8): 1696-1698.

Qian S, Weiss J. 1993. Wavelets and the numerical solution of partial differential equations.JournalofComputationalPhysics, 106(1): 155-175.Restrepo J M, Leaf G K. 1997. Inner product computations using periodized Daubechies wavelets.InternationalJournalforNumericalMethodsinEngineering, 40(19): 3557-3578.Sabina, Mishra V. 2012. Wavelet-Galerkin solutions of one and two dimensional partial differential equations.JournalofEmergingTrendsinComputingandInformationSciences, 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.RadioScience, 29(4): 965-977.

Shi L K, Shen X Q, Yan W L, et al. 2001. A wavelet interpolation Galerkin method for the solution of boundary value problems in 2D electrostatic field.JournalofHebeiUniversityofTechnology(in Chinese), 30(1): 62-66.Suk-In S, Schulz E. 2013. Wavelet-Galerkin solution of a partial differential equation with nonlinear viscosity.AppliedMathematicalSciences, 7(38): 1849-1880.Xu S Z. 1994. The Finite Element Method in Geophysics (in Chinese). Beijing: Science Press.

Yang S Y, Ni G Z. 1999. Wavelet-Galerkin method for the numerical calculation of electromagnetic fields.ProceedingsoftheCSEE(in Chinese), 19(1): 56-61.

Yang S Y, Ni G Z, Ho S L, et al. 2000. Wavelet-Galerkin method for computations of electromagnetic fields-computation of connection coefficients.IEEETransactionsonMagnetics, 36(4): 644-648.Yee K S. 1966. Numerical solution of initial boundary value problems involving Maxwell′s equations in isotropic media.IEEETransactionsonAntennasandPropagation, 14(3): 302-307.

Zhang P W, Liu F Q, Zhang Y. 1995. The numerical computation of wavelets.ComputationalMathematics(in Chinese), (2): 173-185.

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.ChineseJ.Geophys. (in Chinese), 48(5): 1156-1166.

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

陳雅琴. 2008. 橋梁結(jié)構(gòu)分析的廣義變分原理-Daubechies條件小波法研究[博士論文]. 西安: 長安大學(xué).

程正興. 2006. 小波分析與應(yīng)用實(shí)例. 西安: 西安交通大學(xué)出版社.

底青云, 王妙月. 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.

高友蘭. 2009. 數(shù)值積分的Daubechies小波方法. 江南大學(xué)學(xué)報(bào)(自然科學(xué)版), 8(1): 122-125.

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

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

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

石陸魁, 沈雪勤, 顏威利等. 2001. 小波插值Galerkin法解二維靜電場中的邊值問題. 河北工業(yè)大學(xué)學(xué)報(bào), 30(1): 62-66.

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

楊仕友, 倪光正. 1999. 小波-伽遼金有限元法及其在電磁場數(shù)值計(jì)算中的應(yīng)用. 中國電機(jī)工程學(xué)報(bào), 19(1): 56-61.

張平文, 劉法啟, 張宇. 1995. 小波函數(shù)值的計(jì)算. 計(jì)算數(shù)學(xué), (2): 173-185.

張新明, 劉克安, 劉家琦. 2005. 流體飽和多孔隙介質(zhì)二維彈性波方程正演模擬的小波有限元法. 地球物理學(xué)報(bào), 48(5): 1156-1166.

(本文編輯何燕)

基金項(xiàng)目國家自然科學(xué)基金項(xiàng)目(41574116,41074085)，中南大學(xué)創(chuàng)新驅(qū)動項(xiàng)目(2015CX008)，中南大學(xué)升華育英人才計(jì)劃，教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃(NCET-12-0551)，中南大學(xué)教師研究基金(2014JSJJ001)，湖湘青年創(chuàng)新創(chuàng)業(yè)平臺培養(yǎng)對象項(xiàng)目共同資助.

作者簡介馮德山，男，1978年生，漢族，湖南祁陽人，博士，教授，從事地球物理數(shù)據(jù)處理與正反演研究. E-mail:fengdeshan@126.com

doi:10.6038/cjg20160129 中圖分類號P631

收稿日期2014-10-31，2015-06-30收修定稿

Daubechies wavelet finite element method for solving the GPR wave equations

FENG De-Shan1,2， YANG Bing-Kun1,3， WANG Xun1,2， DU Hua-Kun1,2

1SchoolofGeosciencesandInfo-Physics，CentralSouthUniversity，Changsha410083，China2Keylaboratoryofmetallogenicpredictionofnon-ferrousmetalsandGeologicalEnvironmentMonitor(CentralSouthUniversity)，MinistryofEducation,Changsha410083，China3FujianAcademyofBuildingResearch,Fuzhou350025，China

AbstractBased on the separable wavelet theory, we construct the two-dimensional Daubechies wavelet bases by means of one-dimensional Daubechies scaling functions, which is used for interpolation functions of solving the GPR wave equation, thus present the discrete format of two-dimensional Daubechies wavelet finite element GPR equation. By introducing a transformation matrix, the transformation between the wavelet coefficient space and the GPR electromagnetic field is implemented. By introducing the degree of freedom condensation technique, it effectively solves the problem of too much freedom in internal wavelet unit during the solution process of the wavelet finite element, reducing the amount of calculation and can be coupled easily with traditional finite element method. Then the calculation formulas of connection coefficient used in Daubechies wavelet finite element are elaborated, which effectively resolve the difficulty and core problem in solving partial differential equations by wavelet finite element. Finally, with two typical GPR models as example, comparing the radar forward sections and the single waveforms between Daubechies wavelet finite element method and the traditional finite element method, and the result shows that under the conditions of the same dividing method and the number of nodes, the compact support and orthogonality of Daubechies wavelet finite element improves the solving efficiency to some extent, and it can be fitted well with the solving result of finite element method, validating the correctness of the Daubechies wavelet finite element method, which provides a new idea for solving the GPR wave equation.

KeywordsGround Penetrating Radar； Daubechies wavelet finite element method; Degree of freedom condensation technique; Connection coefficient; Wave equation; Forward modeling

馮德山, 楊炳坤, 王珣等. 2016. Daubechies小波有限元求解GPR波動方程.地球物理學(xué)報(bào)，59(1):342-354,doi:10.6038/cjg20160129.

Feng D S, Yang B K, Wang X, et al. 2016. Daubechies wavelet finite element method for solving the GPR wave equations.ChineseJ.Geophys. (in Chinese),59(1):342-354,doi:10.6038/cjg20160129.