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

冯德山, 杨炳坤, 王珣, 杜华坤

1 中南大学 地球科学与信息物理学院, 长沙　410083 2 中南大学有色金属成矿预测与地质环境监测教育部重点实验室, 长沙　410083 3 福建省建筑科学研究院， 福州　350025



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小波有限元； 自由度凝聚技术； 联系系数； 波动方程； 正演模拟

1引言

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

小波分析是调和分析发展史上里程碑式的进展(程正兴，2006)，它比传统Fourier分析能更好地处理局部奇异性问题，具有多分辨特性(Qian and Weiss，1993；Ko et al，1995).近年来，小波分析与有限元方法相结合，产生了小波有限元，它作为一种偏微分方程(PDE)的潜在高效求解方法被提出，能将分析对象依次放入一个逐级扩大、互相嵌套的函数空间序列…，V-1,V0,V1…中进行分析，而且作为空间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.，2004，2006)、裂纹织别(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算法形成有益补充.

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

2.1Daubechies小波性质

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

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

(1)

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

(2)

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

(3)

(4)

式(3)中

(5)

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

(6)

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

(7)

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

(8)

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

2.2二维张量积小波构造

(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)是相应的尺度函数.图1为根据张量积构造的二维DB3小波尺度函数图.

图1　DB3小波二维张量积尺度函数图

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

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

(11)

×φ(η-l),

(12)

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

(13)

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

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

图2　25个自由度的零尺度二维DB3小波单元

2-2N≤k,l≤0,

(14)

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

(15)

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

(16)

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

(17)

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

(18)

在单元内积分得：

(19)

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

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

(20)

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

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

(21)

进一步展开得：

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

(22)

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

(23)

则(22)式可化为：

(24)

将系数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)时，式中的一阶及二阶导数可采用中心差分来近似逼近：

(30)

(29)式可化为：

(31)

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

(32)

(33)

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

(34)

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

(35)

(36)

(37)

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

4自由度凝聚技术

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

(38)

将(38)式代入(32)式，并在等式两边同乘以PT，可得

(39)

(40)

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

(41)

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

(42)

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

(43)

化简(43)式可得

(44)

5Daubechies小波联系系数计算

图3　DB3小波尺度函数一阶导数图

图4　DB3小波尺度函数一阶导数内积图

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

(45)

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

(46)

(47)

令ξ=2x，则

(48)

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

(49)

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

(50)

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

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

(55)

(56a)

(56b)

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

表1　DB3小波联系系数k1,k2值

表2　DB3小波联系系数k1,k2值

图5　DB3小波尺度函数联系系数计算结果图(a) 联系系数的计算结果； (b) 联系系数的计算结果.

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.

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

图6　矩状雷达地电模型示意图

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

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

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.可见，尽管小波有限元增加了小波有限元联系系数与凝聚技术的转换运算，但是由于这些计算能事先计算好再调用，并不增加太多计算工作量，相反小波有限元具有的正交性，使得求解的刚度矩阵为极度稀疏矩阵，可提高计算效率.

图7　矩状模型雷达正演剖面图(a) FEM方法； (b) 小波有限元.

图8　矩状模型第80道波形数据对比(a) 0～24 ns；(b) 12～24 ns.

图9　组合雷达模型示意图

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

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

图10　组合模型雷达正演剖面图(a) FEM方法；(b) 小波有限元.

图11　组合模型的第40道雷达波形对比(a) 0～24 ns；(b) 4～24 ns.

图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.

7结论

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

(2) 编制了二维张量积Daubechies小波有限元GPR正演程序，通过对比相同的剖分方式及节点数目条件下GPR正演剖面图与单道波形图，表明小波有限元与传统有限元模拟结果能较好地吻合，证明小波有限元算法的正确性.由于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.

附中文参考文献

陈雅琴. 2008. 桥梁结构分析的广义变分原理-Daubechies条件小波法研究[博士论文]. 西安: 长安大学.

程正兴. 2006. 小波分析与应用实例. 西安: 西安交通大学出版社.

底青云, 王妙月. 1999. 雷达波有限元仿真模拟. 地球物理学报, 42(6): 818-825.

冯德山, 陈承申, 戴前伟. 2010. 基于UPML边界条件的交替方向隐式有限差分法GPR全波场数值模拟. 地球物理学报, 53(10): 2484-2496, doi: 10.3969/j.issn.0001-5733.2010.10.022.冯德山, 陈承申, 王洪华. 2012. 基于混合边界条件的有限单元法GPR正演模拟. 地球物理学报, 55(11): 3774-3785, doi: 10.6038/j.issn.0001-5733.2012.11.024.

高友兰. 2009. 数值积分的Daubechies小波方法. 江南大学学报(自然科学版), 8(1): 122-125.

何正嘉, 陈雪峰, 李兵等. 2006. 小波有限元理论及其工程应用. 北京: 科学出版社.

李静, 曾昭发, 吴丰收等. 2010. 探地雷达三维高阶时域有限差分法模拟研究. 地球物理学报, 53(4): 974-981, doi: 10.3969/j.issn.0001-5733.2010.04.022.

刘四新, 曾昭发. 2007. 频散介质中地质雷达波传播的数值模拟. 地球物理学报, 50(1): 320-326.

石陆魁, 沈雪勤, 颜威利等. 2001. 小波插值Galerkin法解二维静电场中的边值问题. 河北工业大学学报, 30(1): 62-66.

徐世浙. 1994. 地球物理中的有限单元法. 北京: 科学出版社.

杨仕友, 倪光正. 1999. 小波-伽辽金有限元法及其在电磁场数值计算中的应用. 中国电机工程学报, 19(1): 56-61.

张平文, 刘法启, 张宇. 1995. 小波函数值的计算. 计算数学, (2): 173-185.

张新明, 刘克安, 刘家琦. 2005. 流体饱和多孔隙介质二维弹性波方程正演模拟的小波有限元法. 地球物理学报, 48(5): 1156-1166.

(本文编辑何燕)

基金项目国家自然科学基金项目(41574116,41074085)，中南大学创新驱动项目(2015CX008)，中南大学升华育英人才计划，教育部新世纪优秀人才支持计划(NCET-12-0551)，中南大学教师研究基金(2014JSJJ001)，湖湘青年创新创业平台培养对象项目共同资助.

作者简介冯德山，男，1978年生，汉族，湖南祁阳人，博士，教授，从事地球物理数据处理与正反演研究. 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波动方程.地球物理学报，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.