地球物理学报  2016, Vol. 59 Issue (1): 342-354   PDF    
Daubechies小波有限元求解GPR波动方程
冯德山1,2, 杨炳坤1,3, 王珣1,2, 杜华坤1,2    
1. 中南大学地球科学与信息物理学院, 长沙 410083;
2. 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙 410083;
3. 福建省建筑科学研究院, 福州 350025
摘要: 基于可分离小波理论,由一维Daubechies尺度函数的张量积构造二维Daubechies小波基,并将它作为GPR波动方程求解的插值函数,导出了二维Daubechies小波有限元GPR方程离散格式;通过引入转换矩阵,实现小波系数空间与雷达场值之间转换.引入自由度凝聚技术,有效解决了小波有限元求解中小波单元内部自由度过多的问题,节约了计算量并方便与传统有限元法耦合.然后,详细阐述了Daubechies小波有限元联系系数计算方法,有效解决了小波有限元求解偏微分方程的难点与核心问题.最后,以两个典型GPR模型为例,对比了Daubechies小波有限元与传统有限元的雷达正演剖面图与单道波形图,结果表明:在相同的剖分方式及节点数目条件下,Daubechies小波有限元的紧支性与正交性一定程度上提高了求解效率,它与有限元法求解结果能较好地吻合,验证了Daubechies小波有限元算法的正确性.
关键词: 探地雷达     Daubechies小波有限元     自由度凝聚技术     联系系数     波动方程     正演模拟    
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    
1. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
2. Key Laboratory of Metallogenic Prediction of Non-Ferrous Metals and Geological Environment Monitor(Central South University), Ministry of Education, Changsha 410083, China;
3. Fujian Academy of Building Research, Fuzhou 350025, China
Abstract: Based 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.
Key words: Ground Penetrating Radar     Daubechies wavelet finite element method     Degree of freedom condensation technique     Connection coefficient     Wave equation     Forward modeling    
1 引言

GPR正演传统算法主要有:有限差分法(刘四新和曾昭发,2007李静等,2010冯德山等,2010)、有限元法(底青云和王妙月,2010;冯德山等,2012),它们理论体系日趋成熟完善,但都是选用多项式作为基函数,是全局化的函数,当开展简单的模型计算时具有优势.然而,随着GPR工程勘探日益复杂化、精细化,如仍以低阶多项式的基函数去逼近场函数,难以完全消除局部大梯度问题所引起的振荡,影响求解精度.以有限单元法(FEM)为例,它在求解断层断点、地质裂隙等局部大梯度、奇异性问题时,计算误差取决于对模型的一次性剖分,当计算结果存在较大误差时,需要加密网格或提高插值多项式的阶次,原来形成的刚度矩阵不能够在网格细化后的计算中被继承,势必增加大量计算成本,浪费计算资源,所以寻找更优的插值函数也成为当前有限元发展的新研究方向之一.

小波分析是调和分析发展史上里程碑式的进展(程正兴,2006),它比传统Fourier分析能更好地处理局部奇异性问题,具有多分辨特性(Qian and Weiss,1993Ko et al,1995).近年来,小波分析与有限元方法相结合,产生了小波有限元,它作为一种偏微分方程(PDE)的潜在高效求解方法被提出,能将分析对象依次放入一个逐级扩大、互相嵌套的函数空间序列…,V-1V0V1…中进行分析,而且作为空间Vj+1的子空间Vj,存在相应的补空间Wj,当需要提高求解精度时,通过增加Wj空间以扩大单元容许空间至Vj+1.因此,小波有限元能根据实际需要任意改变分析尺度,在不改变网格剖分的前提下提高分辨率,使其可以在大梯度处采用小的分析尺度、高阶单元以提高分析精度,而在小梯度处采用大的分析尺度、较低阶单元以提高分析效率.在众多的小波当中,Daubechies小波(Daubechies,1988)因为具有紧支撑性、正交性等诸多优点,已引起国内外学者的广泛关注.Daubechies小波的紧支性保证了在没有截断误差的条件下,能以最小的单元自由度最大限度地逼近待求函数,通过阈值运算,达到待求系数最少;而正交性使得小波有限元的刚度矩阵是带状稀疏矩阵,从而可以减少代数方程组的求解运算量;小波函数的消失矩特性指明Daubechies小波函数ΦN(x)可以精确地表征出不大于N-1阶的幂级数(Daubechies,1992).

Amaratunga等(1994)采用小波Galerkin法结合Dirichlet边界条件求解一维Helmholtz方程及 二维Green方程(Amaratunga and Williams,1993),指出了小波嵌套空间能在不同尺度下求解的优势;Sarkar等(1994)将小波函数引入到传统有限元插值函数中求解Maxwell方程,所得的系数矩阵呈对角线的稀疏分布,具有条件数不随维数增加的优点;Patton和Marks(1996)构造了一维Daubechies小波单元;西安交通大学机械自动化研究所在何正嘉等(2006)的带领下,涌现出一大批小波有限元数值模拟(Chen et al.,20042006)、裂纹织别(Dong et al.,2009)等应用研究优秀成果;杨仕友与倪光正(1999)将无穷区间Daubechies小波有限元法应用到了电磁场数值计算中;石陆魁等(2001)采用小波-Galerkin法求解了二维多介质静态电磁场边值问题;张新明等(2005)把小波有限元法引入到流体饱和多孔隙介质二维波动方程的正演模拟中;Chen等(2006)应用Daubechies小波开展了动态多尺度提升计算,并对小波有限元联系系数、刚度矩阵与载荷列阵的计算进行了详细的探讨;陈雅琴(2008)在其 博士论文中提出一种基于广义变分原理的Daubechies 条件小波法,使小波Galerkin法和小波Ritz法的求解精度得到了一定的提高;Mishra和Sabina(2011)应用小波Galerkin法求解一维谐波常微分方程及二维偏微分方程(Sabina and Mishra,2012);Suk-In和Schulz(2013)应用小波Galerkin法求解非线性黏稠Burgers偏微分方程.综上所述,Daubechies小波有限元的研究取得了一些成果,但在地球物理领域仍处于起步阶段,有待进一步探索与发展.本文应用Daubechies小波有限元开展GPR正演,能对目前GPR算法形成有益补充.

2 Daubechies小波特性及二维张量积小波构造 2.1 Daubechies小波性质

Daubechies小波自问世以来,引起了众多学者的广泛关注,其理论和应用研究异常活跃,主要性质有(Daubechies,1992):

(1)紧支性:紧支性反映尺度函数与小波函数在时域的局部化能力,支撑区间越小,局部化能力越强.对给定的正整数N≥2,尺度函数N(x)和小波函数ψN(x)的支撑域分别为:

(2)消失矩特性:Daubechies小波的光滑性可以用消失矩来刻画,消失矩阶数越大则小波函数越光滑,追求好的光滑性,必然以扩大支撑区间,降低时域局部化能力为代价的.Daubechies小波的消失矩为N,即

(3)正交性:Daubechies小波是正交小波,其尺度函数 ΦN(x)与小波函数ψN(x)满足如下正交条件:

式(3)中

(4)归一性:在Daubechies小波的构造中,常要求尺度函数ΦN(x)满足归一性条件

(5)两尺度方程:由于Daubechies小波没有明确的数学表达式,Daubechies小波的尺度函数ΦN(x)和小波函数ψN(x)由两尺度方程构造而来:

其中pN(k)称为小波逼近空间低通滤波系数,qN(k)称为小波空间带通滤波系数.

对于给定的正整数N,Daubechies小波仅有2N个pN(k)不等于零.为了便于叙述,将具有2N个非零滤波系数的Daubechies小波简称为DBN小波.

2.2 二维张量积小波构造

求解二维GPR偏微分方程时,电场 E 及磁场 H 是关于时间变量t和空间变量x,y的函数,需要采用二维多分辨分析,根据可分离小波理论,二维小波基可以用一维小波基的张量积来表示(Restrepo and Leaf,1997).现设一元尺度函数Φ1(x)与Φ2(y)分别生成一个多分辨分析{Vj1}与{Vj2},则空间{Vj1}与{Vj2}的张量积空间为:

式(9)中,ⓧ是Kronecker符号.Vj1的基底是{2j/2Φ1(2jx-k)},而Vj2的基底是{2j/2Φ2(2jy-k)}.因此,对于二元函数Φ(x,y)引入记号

式(10)中,-(2jL+2N-1)≤kl≤0.Lxy的定义域[0,L].则{Φj,k,l(x,y)}是Vj的基底,所以{Vj}形成L2(R2)中的一个多分辨分析,而Φ(x,y)是相应的尺度函数.图 1为根据张量积构造的二维DB3小波尺度函数图.

图 1 DB3小波二维张量积尺度函数图Fig. 1 The map of two-dimensional tensor product scaling function of DB3 wavelet
3 GPR波动方程的Daubechies小波有限元求解

由电磁波传播理论(Yee,1996)可知,含衰减项的GPR波动方程为:

式(11)中,当 U 表示电场值 E 时,S为激励源Se,当 U 表示磁场值 H 时,S为激励源Shε为介电常数,μ为磁导率,σ为电导率,t为时间,Γ U 表示模拟边界.本文采用满足Neumann边界条件的Galerkin原理来推导GPR小波有限元方程(徐世浙,1994). 首先,将待求解区域离散为若干个单元,在单元内有:

式(12)中ak,l(t)为待求的小波系数,它仅与时间有关;Φ为Daubechies小波尺度函数,仅与空间有关,U e为单元内的电场值或磁场值;x,y为总体坐标系,ξ,η为局部坐标系;以矩阵的形式表示为:

其中,Φ =[Φ(ξ)Φ(η)],A e为待求小波系数组成的列向量,其SuppφN(x)=[0,2N-1].以V0空间上的DB3小波为例,1个二维小波单元自由度总数为:(2N-1)×(2N-1),即25个单元自由度,则构造的具有25个节点的小波单元如图 2所示.

图 2 25个自由度的零尺度二维DB3小波单元Fig. 2 Two-dimensional DB3 wavelet element with 25 degrees of freedom under zero scale

图 2构造的小波单元的局部坐标取值范围是[-1,1],结合Daubechies小波尺度函数的支撑区间SuppφN(x)=[0,2N-1],可以求得k,l的取值范围为:

则局部坐标与总体坐标的转换关系如下:

其中,x0y0为总体坐标下子单元的中心坐标,a、b是总体坐标下子单元的边长.-1≤ξ,η≤1,则有如下微分关系:

则由总体坐标与局部坐标关系可知,若设在总体坐标中满足:

利用Galerkin法,将(17)式代入(11)式,两边同时乘以二维Daubechies小波基函数:

在单元内积分得:

对式(19)中左边第2项采用Green公式变换,可得到

再将(20)式中的 U 用(17)式代入,可得到:

进一步展开得:

则(22)式可化为:

将系数ak,l(t)按照一定的顺序重新组合形成一维向量 A,即

再作如下假定:

则(24)式可改写为:

M k1k2l1l2T k1k2l1l2 K k1k2l1l2是(2N-1)×(2N-1)阶矩阵,À为时间的一阶导数,Ä 为时间的二阶导数.当将模拟区域划分为n个小波单元时,将相应的系数矩阵 M k1k2l1l2T k1k2l1l2 K k1k2l1l2 F 1拓展成n×n单元的矩阵 M,T,K,F,而每个单元的矩阵为(2N-1)×(2N-1)阶矩阵.

在求解方程组(29)时,式中的一阶及二阶导数可采用中心差分来近似逼近:

(29)式可化为:

当零时刻或-Δt时刻,电场值为零,故ak,l(0)=0,ak,l(-Δt)=0,且激励源 F 为已知值.因此,可以通过解上面方程逐步求得不同时间层位上的小波系数,故式(31)可化简为

求得小波系数值后,再代入公式(33)就可以求出任意点上的函数值,即雷达波场值.

由于小波有限元单元平衡方程建立于小波空间,此时的单元自由度是小波系数而非电场值,为了实现小波系数空间与GPR波动方程电场(或磁场)空间之间的变换,可以构造快速小波变换(张新明等,2005何正嘉等,2006):

E 为单元电场向量,Φ 为小波系数变换矩阵,A 为小波系数向量.仍以V0尺度的25节点DB3小波单元为例,以节点电场作为单元自由度,则

式(37)中的Φ元素为二维DB3尺度函数整数点及二分点上的函数值.

4 自由度凝聚技术

采用DBN小波尺度函数构造的二维小波单元时,每个单元有(2N-1)×(2N-1)个自由度,相对传统FEM法增加了很多单元内部自由度,单元矩阵将变得比较庞大,给单元离散和单元组合带来一定困难.因此,有必要引入自由度凝聚技术将多余的单元自由度消去,剩下需要的自由度(Mehraeen and Chen,2006).假如记矩阵 Φ 的逆矩阵为 P,则由(34)式可得:

将(38)式代入(32)式,并在等式两边同乘以 P T,可得

若记: K = P TKPB = P T B,则(39)式可简记为:

将(40)式的矩阵重新排列并进行分块,得

式(41)中 E i=[E1 E5 E21 E25]为单元内4个角点的场值; E j=[E2 E3 E4 E6E20 E22 E23 E24] 为单元内部节点的场值; B iBj分别为4个角点及内部节点对应小波系数经矩阵转换后的值;

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

将(42)式再代回式(41)第1个等式,就消去了除角点之外的自由度,得

化简(43)式可得

式(44)中 K= Kii-KijKjj-1K ji B = Bi-KijKjj-1Bj.得到转换后的单元小波有限元方程后,即可和传统有限元类似,建立整体平衡方程,加载边界条件.

5 Daubechies小波联系系数计算

本文计算过程中需要用到的Daubechies小波的尺度函数φN(k)、尺度函数的导数值φ(d)N(k)的求解,可参照文章(陈雅琴,2008高友兰,2009张平文等,1995)的算法,在此不再赘述.本文以DB3小波为例,根据该算法求得的二分点{k=i/2n,i=0,1,2,…,2N-1,nZ +}处尺度函数的导数值如图 3所示.根据二维小波可分离理论及图 3中求得的尺度函数的导数值φ(d1)N(k)φ(d2)N(k),很容易求得尺度函数的内积值.图 4为DB3小波尺度函数一阶导数内积值.

图 3 DB3小波尺度函数一阶导数图Fig. 3 Scale function one order derivative of DB3 wavelet

图 4 DB3小波尺度函数一阶导数内积图Fig. 4 Scale function one derivative inner product of DB3 wavelet

下面着重介绍式(23)中的Daubechies小波联 系系数的计算.由于式(23)中的积分区间为[-1,1],而大部分联系系数求解方法采用的区间为[0,1],为此可将任意长度为L的区域投射到[0,1]的区间,譬如:

本文中仅介绍x的积分区间为[0,1],结合Daubechies小波尺度函数的支撑区间,可求得k1k2的取值范围为:

为了方便公式推导,引入特征函数 χ[0,1](x)

ξ=2x,则

根据式(7)的尺度函数表达式,可得

式(49)中s取值范围与k1k2一致.将式(49)两边同求d次导,则

本文采用Vj尺度下联系系数通用式求解过程作为范例,V0尺度联系系数可类似得到.

根据式(50)可得

式(51)中,s,t的取值范围与k1k2一致,且必须保证2-2Ns-2k1t-2k2s-2k1+2jt-2k2+2j≤2j-1. 将式(51)写成矩阵形式,可得如下方程组:

式(52)中:

由于 Λk1k2j,d1d2≠ 0,根据齐次方程组有非零解的条件,可知:

等式(54)左边矩阵为奇异阵,不能唯一确定联系系数的值.目前常用的求解方法是:首先由求解矩阵的性质得出联系系数的解空间,引入与解空间自由度个数相等的附加方程,从而组成恰定方程组,求解后可唯一确定准确的联系系数值.对于添加的附加方程,Latto等(1999)Yang等(2000)陈雪峰等(2006)给出如下的计算式:

用此方程组作为附加方程组时,需要求解系数pj,kqpj,kw的值,但是因为q、w取值范围比较广,计算结果会随附加方程个数的不同而改变,难以确定哪组解是最优的联系系数值.实际计算中可通过选取更严格的限定条件:d1+d2<q=wN-1,通常可以保证添加的补充方程对于求解联系系数严格有效,从而可以获得准确的计算值.根据上述算法,利用Matlab编写计算程序,V0尺度下的DB3小波有限元联系系数如下:

求解公式(56a)秩为24,待求联系系数总数为25,仅需添加式(55)中1个方程式.求解公式(56b)秩为22,待求联系系数总数为25,需添加式(55)中3个方程式.2个联系系数计算结果分别参见表 1表 2,利用该数据绘出的联系系数如图 5a图 5b所示.其中 表 1中的计算结果若保留8位有效位数字,与Frassati(2009)L and ragin- 文章中计算结果完全一致,证明程序的正确性.

表 1 DB3小波联系系数Λk1,k20,0,0 Table 1 Value of DB3 wavelet coefficient Λk1,k20,0,0

表 2 DB3小波联系系数Λk1,k20,1,1 Table 2 Value of DB3 wavelet coefficient Λk1,k20,1,1

图 5 DB3小波尺度函数联系系数计算结果图 (a) 联系系数Λk1,k20,0,0的计算结果; (b) 联系系数Λk1,k20,1,1的计算结果. Fig. 5 Maps of calculating values of connection coefficient with DB3 scaling function (a) Results of connection coefficients Λk1,k20,0,0;(b) Results of connection coefficients Λk1,k20,1,1.
6 数值算例 6.1 矩状模型

选取图 6所示矩状地电模型,模拟区域为4.0 m×2.0 m,上层混凝土介质的相对介电常数为8.0,电导 率为0.0001 S·m-1,厚度为0.8 m,下层介质介电常数为15.0,电导率为0.01 S·m-1. 上层介质中埋有矩状异常体,其介电常数为3.0,电导率为0.002 S·m-1,异常体长0.4 m,宽为0.2 m,异常体顶边距地面 0.4 m,波源为脉冲Ricker子波,中心频率取900 MHz,采样时窗长度为30 ns,采样时间间隔为0.02 ns.

图 6 矩状雷达地电模型示意图Fig. 6 The sketch map of rectangle GPR model

基于FEM算法对该模型进行正演,整个区域由单元边长0.025 m的正方形剖分为160×80的网格空间.采用DB3小波有限元进行正演,共计40× 20个小波单元,每个小波单元网格被DB3小波插入 3×3个内部节点,每个小波单元细分为4×4区间,则最终的小波区间由边长0.025 m的正方形网格组成,剖分总网格与正方形FEM算法相当为160×80 个网格.

图 7a为160道数据的FEM正演剖面图,图中可见,8 ns位置为矩状异常体的上界面,由于绕射波的影响,异常体的两个棱角处出现了绕射波,10 ns 左右可见异常体下界面反射.图 7b为DB3小波有 限元法正演剖面图,观察图 7a图 7b,两图非常类似,都能反映异常体的形态,但是小波有限元多次波波形更为离散一些,推断与小波有限元联系系数的计算精度有关联.

图 7 矩状模型雷达正演剖面图 (a) FEM方法; (b) 小波有限元. Fig. 7 The section of GPR simulation (a) FEM method; (b) Wavelet finite element method.

图 8为两种方法第80道雷达波形对比图,图 8a为0~24 ns时的显示,说明两种算法能很好地拟合.图 8b为了突出波场细微的差别,仅显示12~24 ns时的单道雷达波形,图中可见,小波有限元与FEM计算结果之间仍存在细小的差别,但总体拟合趋势很好,说明了小波有限元计算结果的可靠性.

图 8 矩状模型第80道波形数据对比 (a) 0~24 ns;(b) 12~24 ns. Fig. 8 Contrasting between FEM and Daubechies wavelet FEM of the 80th waveform of rectangle GPR model
6.2 组合模型

为了说明小波有限元对复杂模型模拟的有效性,选取图 9所示的组合雷达地电模型验证.模拟区域、背景介质、网格剖分方式及雷达频率、采样参数与上例一致.在上层介质中左边为矩状空洞异常体,其长与宽都为0.2 m;中间为菱形素填土介质,其介 电常数为3.0,电导率为0.002 S·m-1,菱形长0.6 m,高为0.3 m,其上顶点距模拟区域上边界0.3 m;右 边为半径为0.05的金属球体.在Intel(R)CoreTMi7-4500U CPU@1.80 GHz,8.00 GB的内存物理地址扩展,Window 8操作系统IBM X240s笔记本上计算该模型雷达剖面160道数据,采用DB3小波有限元法计算时间为6903.25 s,而采用传统有限元计算时间为7845.84 s.可见,尽管小波有限元增加了 小波有限元联系系数与凝聚技术的转换运算,但是由于这些计算能事先计算好再调用,并不增加 太多计算工作量,相反小波有限元具有的正交性,使得求解的刚度矩阵为极度稀疏矩阵,可提高计算效率.

图 9 组合雷达模型示意图Fig. 9 The sketch map of GPR combination model

图 10a图 10b分别为应用FEM及小波有限元正演所得的160道波形的正演剖面图,从图中可见,横坐标2.0 m,纵坐标6 ns处为菱形的上界面,菱形的上面两条边也较好地对应了图中的反射波;左边横坐标1.0 m,纵坐标9.5 ns及11 ns处两处绕射波较好地对应了空洞异常体上下界面;而右边横坐标3.0 m纵坐标10 ns处为金属球体上界面绕射波,由于金属球体的全反射,导致其下界面基本见不到绕射波.对比图 10a图 10b,两图波形大致类似,仅多次波的体现上存在细微差别,小波有限元波形更为丰富,在增益倍数均相同条件下,小波有限元幅值稍强于有限元法.

图 10 组合模型雷达正演剖面图 (a) FEM方法;(b) 小波有限元. Fig. 10 The GPR simulation section of combination model (a) FEM method; (b) Wavelet finite element method.

图 11a为0~24 ns两种方法第40道波形对比,黑线表示的FEM计算结果与红色点表示的小波有限元计算结果能较好地吻合.图 11b为仅显示4~24 ns时的单道雷达波形,由于去掉了1.5 ns左 右直达波的大振幅值,波场细节得已体现,图中可见,在10.0 ns处的异常体波形幅值稍大.

图 11 组合模型的第40道雷达波形对比 (a) 0~24 ns;(b) 4~24 ns. Fig. 11 Contrasting between FEM and wavelet FEM of the 40th GPR waveform of combination model

图 12a为第120道0~24 ns的雷达波形数据,两种方法计算结果同样能较好地拟合,10.0 ns处为圆形球体的上界面反射.图 12b为仅显示4~24 ns时的单道雷达波形,图中可见,小波有限元与FEM计算结果仅在16.0~20.0 ns之间的多次反射波及界面反射之间存在细微的差别,但两者总体拟合很好,说明了小波有限元计算结果的可靠性.

图 12 组合模型的第120道雷达波形数据对比 (a) 0~24 ns;(b) 4~24 ns. Fig. 12 Contrasting between FEM and wavelet FEM of the 120th GPR waveform of combination model
7 结论

(1)详细阐述了Daubechies小波有限元联系系数计算方法,有效解决了小波有限元求解偏微分方程的难点与核心问题.通过引入自由度凝聚技术消除单元内部的自由度,剩下四个角点的自由度,解决了小波单元自由度过多的问题,该技术能在小波有限元与有限元耦合计算中得到广泛应用.

(2)编制了二维张量积Daubechies小波有限元GPR正演程序,通过对比相同的剖分方式及节点数目条件下GPR正演剖面图与单道波形图,表明小波有限元与传统有限元模拟结果能较好地吻合,证明小波有限元算法的正确性.由于Daubechies小波的正交性、紧支性使得Daubechies小波有限元刚度矩阵为条带状稀疏矩阵,提高了求解效率,为GPR波动方程求解提供了一种新的思路.

参考文献
[1] Amaratunga K, Williams J R. 1993. Wavelet based Green's function approach to 2D PDES. Engineering Computations, 10(4):349-367.
[2] 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.
[3] 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.
[4] 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.
[5] 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.
[6] Cheng Z X. 2006. Wavelet Analysis and its Application Examples(in Chinese). Xi'an:Xi'an Jiaotong University Press.
[7] Daubechies I. 1992. Ten Lectures on Wavelets(CBMS-NSF Regional Conference Series in Applied Mathematics, 61. Philadelphia:Society for Industrial and Applied Mathematics.
[8] Daubechies I. 1988. Orthonormal bases of compactly supported wavelets. Communications on Pure and Applied Mathematics, 41(7):909-996.
[9] Di Q Y, Wang M Y. 1999. 2D finite element modeling for radar wave. Chinese J. Geophys.(in Chinese), 42(6):818-825.
[10] 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. Mechanical Systems and Signal Processing, 23(3):869-883.
[11] 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.
[12] 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.
[13] Gao Y L. 2009. Method of numerical integration with Daubechies wavelet. Journal of Jiangnan University(Natural Science Edition)(in Chinese), 8(1):122-125.
[14] 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.
[15] Ko J, Kurdila A J, Pilant M S. 1995. A class of finite element methods based on orthonormal, compactly supported wavelets. Computational Mechanics, 16(4):235-244.
[16] 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. Optics Express, 17(21):18433-18448.
[17] Latto A, Resnikoff L, Tenenbaum E. 1999. The evaluation of connection coefficients of compactly supported wavelets. Aware Inc., Technical Report AD910708.
[18] 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.
[19] 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.
[20] Mehraeen S, Chen J S. 2006. Wavelet Galerkin method in multi-scale homogenization of heterogeneous media. International Journal for Numerical Methods in Engineering, 66(3):381-403.
[21] Mishra V, Sabina. 2011. Wavelet Galerkin solutions of ordinary differential equations. Int. Journal of Math. Analysis, 5(9):407-424.
[22] Patton R D, Marks P C. 1996. One-dimensional finite elements based on the Daubechies family of wavelets. AIAA Journal, 34(8):1696-1698.
[23] Qian S, Weiss J. 1993. Wavelets and the numerical solution of partial differential equations. Journal of Computational Physics, 106(1):155-175.
[24] Restrepo J M, Leaf G K. 1997. Inner product computations using periodized Daubechies wavelets. International Journal for Numerical Methods in Engineering, 40(19):3557-3578.
[25] 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.
[26] 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.
[27] 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. Journal of Hebei University of Technology(in Chinese), 30(1):62-66.
[28] Suk-In S, Schulz E. 2013. Wavelet-Galerkin solution of a partial differential equation with nonlinear viscosity. Applied Mathematical Sciences, 7(38):1849-1880.
[29] Xu S Z. 1994. The Finite Element Method in Geophysics(in Chinese). Beijing:Science Press.
[30] Yang S Y, Ni G Z. 1999. Wavelet-Galerkin method for the numerical calculation of electromagnetic fields. Proceedings of the CSEE(in Chinese), 19(1):56-61.
[31] Yang S Y, Ni G Z, Ho S L, et al. 2000. Wavelet-Galerkin method for computations of electromagnetic fields-computation of connection coefficients. IEEE Transactions on Magnetics, 36(4):644-648.
[32] 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.
[33] Zhang P W, Liu F Q, Zhang Y. 1995. The numerical computation of wavelets. Computational Mathematics(in Chinese),(2):173-185.
[34] 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.
[35] 陈雅琴. 2008. 桥梁结构分析的广义变分原理-Daubechies条件小波法研究[博士论文]. 西安:长安大学.
[36] 程正兴. 2006. 小波分析与应用实例. 西安:西安交通大学出版社.
[37] 底青云, 王妙月. 1999. 雷达波有限元仿真模拟. 地球物理学报, 42(6):818-825.
[38] 冯德山, 陈承申, 戴前伟. 2010. 基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟. 地球物理学报, 53(10):2484-2496, doi:10.3969/j.issn.0001-5733.2010.10.022.
[39] 冯德山, 陈承申, 王洪华. 2012. 基于混合边界条件的有限单元法GPR正演模拟. 地球物理学报, 55(11):3774-3785, doi:10.6038/j.issn.0001-5733.2012.11.024.
[40] 高友兰. 2009. 数值积分的Daubechies小波方法. 江南大学学报(自然科学版), 8(1):122-125.
[41] 何正嘉, 陈雪峰, 李兵等. 2006. 小波有限元理论及其工程应用. 北京:科学出版社.
[42] 李静, 曾昭发, 吴丰收等. 2010. 探地雷达三维高阶时域有限差分法模拟研究. 地球物理学报, 53(4):974-981, doi:10.3969/j.issn.0001-5733.2010.04.022.
[43] 刘四新, 曾昭发. 2007. 频散介质中地质雷达波传播的数值模拟. 地球物理学报, 50(1):320-326.
[44] 石陆魁, 沈雪勤, 颜威利等. 2001. 小波插值Galerkin法解二维静电场中的边值问题. 河北工业大学学报, 30(1):62-66.
[45] 徐世浙. 1994. 地球物理中的有限单元法. 北京:科学出版社.
[46] 杨仕友, 倪光正. 1999. 小波-伽辽金有限元法及其在电磁场数值计算中的应用. 中国电机工程学报, 19(1):56-61.
[47] 张平文, 刘法启, 张宇. 1995. 小波函数值的计算. 计算数学,(2):173-185.
[48] 张新明, 刘克安, 刘家琦. 2005. 流体饱和多孔隙介质二维弹性波方程正演模拟的小波有限元法. 地球物理学报, 48(5):1156-1166.