基于Delaunay三角形的非结构化有限元GPR正演

杜华坤，冯德山，汤井田



基于Delaunay三角形的非结构化有限元GPR正演

杜华坤，冯德山，汤井田

(中南大学地球科学与信息物理学院，有色资源与地质灾害探查湖南省重点实验室，湖南长沙，410083)

介绍传统Bowyer-Watson三角网逐点插入法的原理与实现步骤，并将固定边界限制、Laplacian光顺、边压缩、边分裂、点插入等拓扑变换技术应用于网格剖分的优化；为了使数值解的误差在全域内接近于均匀分布，通过间隔函数法实现点源、线源等网格渐变控制，结合局部粗化或细化技术，建立高质量Delaunay 三角形网格，实现自适应网格剖分。通过1个起伏地表与断层模型网格剖分实例验证非结构化网格对于物性参数分布复杂或几何特征不规则的地电模型的适应性。根据GPR有限元波动方程，应用三角形剖分、线性插值的Galerkin有限单元法进行求解。建立1个复杂GPR地电模型，利用Delaunay 三角形对该GPR地电模型进行自适应网格剖分。研究结果表明：非结构化网格对于物性参数分布复杂或几何特征不规则的地电模型都具有良好的适应性；非结构化三角形网格剖分质量好，单元密度易控制，易于实现自适应有限元，能提高复杂模型正演精度。

Delaunay三角形剖分；非结构化网格；有限单元法；正演模拟；探地雷达

对于复杂地球物理正、反演问题的求解，网格生成的重要性不容低估，它直接影响有限元的求解效率与计算精度。根据拓扑结构关系将网格生成技术分为结构化网格和非结构化网格两大类[1]。结构化网格具有统一的拓扑结构，节点关系可以通过网格编号规律推算出来，网格生成的速度快，数据结构简单[2]，但结构网格不再满足地球物理勘探日趋复杂化的要求。目前，非结构化网格取得了飞速发展，它突破了网格节点的结构性限制，节点和单元的分布是任意的，易于控制网格单元的大小、形状和网格节点的位置，删除或加入节点方便，具有更大的灵活性，可方便、合理分布网格的疏密并实现自适应算法[3−4]。非结构化网格的自动剖分方法有很多，目前流行的方法主要有Delaunay三角化方法、前沿生成法(advancing front technique，AFT)、有限四叉树(finite quad tree method)算法。AFT算法[5]采用递归方法对区域进行分割、搜寻和判断，在模拟区域的边界处能产生非常好的网格，对生成单元的控制能力强，但在网格结尾处时，大、小网格的结合是推进波前法的难点。四叉树[6]是一种递归的数据结构，可有效地控制平面中多等级的几何体，易于实现网格密度控制与自适应分析，但其缺点是生成网格与所选择的初始栅格及其取向有关，边界处单元质量较差，程序实现复杂，不利于并行处理。Delaunay三角剖分[7]是连接所有相邻的Voronoi 多边形的生长中心所形成的三角剖分，是目前最流行的全自动非结构化有限元网格生成方法，它的实现策略主要有逐点插入法、分冶法、三角网生长法共3种。Delaunay三角剖分实施过程中常遵循外接圆准则、最大最小三角形内角准则与Thiessen区域准则。目前，在地球物理领域，Delaunay非结构化剖分也得到了广泛应用，如：Key等[8]开展了基于非结构化网格的2D大地电磁自适应有限元模拟；Li等[9]将该自适应非结构化有限元应用到海洋可控源模拟中；汤井田等[10]采用非结构化三角形网格，开展了2.5D直流电阻率法自适应有限元数值模拟，提高了模拟复杂模型的灵活性；任政勇等[11]提出了一种新的非结构化局部节点加密策略，结合二次有限单元法开展了三维直流电阻率模拟，它具有计算量小、剖分精度高的特点；柳建新 等[12]采用自适应非结构网格剖分实现了起伏地形下的任意复杂地电模型MT二维正演模拟。在此，本文作者结合GPR模型的特点对传统Delaunay三角剖分算法进行优化，并将其用于线性插值FEM的探地雷达二维正演，有效地提高了复杂模型模拟的灵活性与 精度。

1 改进与细化的Bowyer-Watson逐点插入三角网生成算法

Bowyer-Watson逐点插入法是Lawson[13]提出的建立 Delaunay 三角网的算法。该算法在增加新点时，只对新插入的点周围三角形局部产生影响，无需对整个网格进行重构，思路简单，易于编程实现。其基本步骤如下。

步骤1：对所有的点进行遍历，求出点集的边界面网格，定义1个包含所有数据点的初始凸多边形，在实际计算中常选取最小外接矩形。得到辅助矩形框后，连接该矩形的任一条对角线。矩形一分为二，生成2个大三角形，这便是初始网格，这样就建立1个覆盖整个计算区域的原始网格与初始三角网。

步骤2：在已经建立的初始网格的基础上插入点集中1个新点。图1(a)所示是插入1个新的顶点1，然后形成4个新的三角形，再依次插入其他点。图1(b)所示是根据所有外接圆包含新点形成新三角形，删除这些三角形所形成的1个多边形Delaunay 空洞。该插入的新点对 Delaunay 空洞的各顶点是可见的空洞。图1(c)所示是连接插入的新点与Delaunay 空洞的各顶点所形成的新网格。

(a) 插入新项点形成4个三角形；(b) 外接圆包含新点形成Delaunay空洞；(c) 连接插入新点形成新网格；(d) 新点插入完毕；(e) 删除辅助凸多边形；(f) 最终剖分结果

步骤3：重复步骤2，直到所有的散点集合都插入完毕为止。图1(d)所示为所有新点插入完成后的网格。

步骤4：删除辅助凸多边形。图1(e)所示为当点集中全部的点都插入到三角网格后，删除辅助矩状凸多边形所得的网格；同时删除与辅助凸多边形有关的所有三角形，得到的最后三角剖分结果如图1(f)所示。

步骤5：删除所有包含初始多边形顶点(非原始数据点)的三角形。

在Bowyer-Watson逐点插入法实际应用中，当点集较大时，算法速度较慢，若顶点集范围是凹区域，则产生的网格会超出网格边界，可能得到的是非法三角网。因此，仍需要在网格单元的形状质量、单元的数目与稀疏分布、计算效率等方面进行局部优化，本文的优化主要包括以下几方面。

1) 以图2(a)为例，在经典的 Delaunay 三角形网格剖分后，由于固定边界，，，和的限制，导致局部优化过程中最小角最大准则不能得到满足，三角形存在过小的锐角，甚至存在生成的三角形与固定边界交叉的情况，如，和与固定边界交叉，为此采用图2 (b)所示的方法，将边界的中点或三角形的外接圆圆心添加到离散点集中，然后形成Delaunay三角网。

(a) 不考虑固定边界的Delaunay三角形部分；(b) 将边界中点添加进点集

2) 为了进一步提高网格质量，使网格变得更光滑、单元形状更接近等边三角形，需要开展细致的优化工作，实际优化处理主要有2类：第1类为不改变网格的拓扑结构，只对节点进行移动，提高网格质量。图3所示的Laplacian 光顺能使网格之间的变化更光滑。第2类称为拓扑变换，如图4所示，它包括边翻转、边分裂、边压缩、点插入、点删除，这些局部操作能在大部分网格不动的前提下，进行局部改变与更新，有效地节约了计算资源。网格优化操作中常依据如下次序：点的删除→局部调用 Laplacian 光顺→边的压缩→边翻转→边分裂与点插入。

(a) Laplacian 光顺前；(b) Laplacian光顺后

(a) 边翻转；(b) 边分裂；(c) 边压缩；(d) 点的删除

3) 网格渐变控制：开展复杂地球物理模型正演时，模型中常存在点源、线源、井孔、套管等点状或线状的源。为了兼顾计算效率与计算精度，需要在物性参数变化剧烈区域采用细网格，在缓变区采用粗网格，而在点插入的Delaunay三角网剖分生成中，考虑用间隔函数控制法来实现渐变网格的生成[14]。第1种方式为：图5(a)所示的点源疏密过渡控制方法，在点源附近的网格很密，由该点向四周网格逐渐变疏，离该点越近，网格越密；离点最远的点处，网格较粗，网格尺度呈线性变化。第2种方式为线源的疏密过渡控制方法。线源是指某一条直线或曲线附近的网格很密，由该直线或曲线向两边网格逐渐变疏。如图5(b)所示，线段附近网格较密，距较远处的点处网格尺寸最大。

(a) 点源A；(b) 线源AB

4) 网格细化与粗化方法。地球物理正演所用的网格都是在数值求解之前生成，因此，很难在首次计算时较好地解决局部地球物理特征，需要根据后验误差估计算法重新修改网格[15]以满足解的要求：若当前网格尺寸比预测的网格尺寸小，则应加大网格尺寸；若比预测的网格尺寸大，则需对网格进一步加密。因此，在网格剖分过程中，需要对网格进行细化和粗化。细化是当需要增加网格密度时，在现有网格的某个部位插入新点，提高数值解的品质。而粗化则是当要减少网格点的密度时，通过删除点来降低过高精度的解。1:4细化，找出单元格三边的中点处各插入1个新点，并将各个新点相连。粗化时可采用相反的操作来进行。1:4网格细化方法如图6所示。

图6 1:4网格细化方法

利用前面介绍的Delaunay 三角形非结构化网格剖分算法，对图7所示的起伏地表与断层模型进行非结构化Delaunay三角网剖分。图7(a)所示为逐点插入法采用494个节点、901个三角单元初次剖分所生成的Delaunay三角网。从图7可见：地表地形的起伏得到了很好拟合，由于网格剖分过程中实施了渐变控制，在地形变化剧烈区域采用了较小的网格、较多的单元，而在其他区域网格较稀疏。图7(b)所示为采用1:4插值加密后的剖分结果，加密后共计为1 888个节点， 3 604个三角单元，该网格可明显提高数值解的品质，但会占用更多的内存与计算时间。

(a) 初次剖分结果；(b) 1:4加密后的网格剖分

2 GPR波动方程有限元求解

根据电磁波理论，雷达波在媒质中的传播应遵循Maxwell方程组[16]，由Maxwell方程组结合媒质的本构关系，导出雷达波波动方程可以表示为[17]：

式中：为电场强度(V/m)；为磁场强度(A/m)；为磁感应强度(T)；为电位移(C/m2)；为电流密度(A/m2)；为电荷密度(C/m3)。

式(1)~(2)表明磁场和电场满足相同的微分方程。因此，仅以电场波动方程式(1)为例，结合初边值条件，应用三角单元剖分的FEM求解GPR满足的偏微分方程。而Galerkin有限元推导时设定未知函数的导数在边界上有已知值，其结果满足Neumann自然边界条件[18]。三角形单元剖分示意图如图8所示，假定二维求解区域已被剖分成许多个三角形单元，在三角单元内，3个顶点分别为P，P和P，节点按逆时针方向编号依次为，和，其坐标依次为(x，z)，(x，z)和(x，z)，节点函数值依次为E，E和E，这3点组成三角形单元为Δ，其面积用表示。三角形中的点(，)与，和这3点的连线，将Δ分割成3个小三角形Δ，Δ和Δ，其面积分别为Δ，Δ和Δ，点在单元中的位置，也可用下式表示：

式中：L，L和L为面积的比值，是量纲1的数，称为二维面积(自然)坐标，它们是单元上的局部坐标。

图8 三角形单元剖分示意图

其中：

；

根据Galerkin法实现方式，在求解区取任意单元，将式(1)两边同时乘以，并对单元求积分[19]，有

其中：为单元质量矩阵，

(10)

式(8)左边第2项展开为

式(8)左边第3项展开为

式(8)右边项展开为

根据式(9)，(11)，(13)和(15)，可将单元积分式 (8)改写为

式(19)便是时空域GPR的二维电场有限元波动方程。式中：为质量矩阵；为阻尼矩阵；为刚度矩阵；和为时间的一次导数；和为时间的二次导数[20]。对式(19)加载具体的边界条件，结合大型线性方程组高精度、快速求解算法即可求 解[21−22]。

3 GPR正演模拟实例

图9所示为复杂模型示意图。该模型为1个长×宽为10 m×6 m的矩形区域，中间是起伏分界面，上层介质的相对介电常数1为5.0，电导率1设置为0.000 1 S/m；下层介质的相对介电常数2为8.0，电导率2为0.001 0 S/m。在背景介质中存在3个圆状异常体，最上面的圆状异常体为电缆管线的金属外壳，其介电常数为1.0，电导率为106S/m，圆心位置为(6.5 m，2.0 m)，半径为0.05 m；左边圆为空洞，其介电常数为1.0，电导率为0 S/m，圆心位置为(5.5 m，2.8 m)，半径为0.10 m；右边为排水管线其介电常数为80.0，电导率为0.002 0 S/m，圆心位置为(7.5 m，2.8 m)，半径为0.15 m。在3圆的左边还有1个滑移断层，其形状、大小与空间位置如图9所示，其介电常数为20.0，电导率为0.010 0 S/m。采用图9所示的非结构化网格三角形自适应剖分方式，共计365个节点，644个单元。GPR波脉冲激励源的中心频率为100 MHz，模拟中波源为100 MHz的Ricker子波，置于界面下0.1 m，时间步长为0.2 ns，时窗长度为120 ns。

图9 复杂模型的非结构化网格剖分示意图

应用FEM对该模型进行正演，模拟所得的GPR剖面图如图10所示。从图10可见：上、下2层不平分界面处产生了1个能量较强反射界面，界面上起伏波动大、变化剧烈的角点、欠光滑的断棱处，产生了角点绕射波形，通过时深转换计算可知它与图9所示模型中的界面位置相吻合。

图10 复杂模型的正演剖面图

3个圆状异常体所在位置出现了双曲线绕射波，双曲线弧形两翼的宽度都比圆的直径大很多，但双曲线的弧顶能准确对应3个圆状异常体的顶点。对比3个圆产生的反射及绕射波，由于物性参数不一样，产生的反射波的能量强弱也是有区别的，最上边的电缆管线的金属外壳异常体由于产生全反射，故反射波能量最强；而左边的圆形空洞体能量最弱，排水管线的上、下界面的绕射波最清晰。而左边的滑移断层由于形状不规则，产生了较多的绕射波与反射波，很难了解其是哪条边所引起的。通过分析该复杂模型可知，基于非结构化三角剖分的FEM方法，能较好地拟合复杂GPR模型的物性参数分界面，其数值模拟结果能真实、客观地反映雷达波的实际探测波形特征。

4 结论

1) 在传统的Bowyer−Watson逐点插入法的基础上，将Laplacian光顺、固定边界限制、拓扑变换优化技术应用于网格剖分中，通过间隔函数网格渐变控制法和局部粗化或细化技术，建立了高质量Delaunay 三角形网格，实现自适应FEM网格剖分。

2) 非结构化网格对于物性参数分布复杂或几何特征不规则的模型都具有良好的适应性，网格剖分质量好，计算效率高，单元密度易控制。

3) 三角形非结构化剖分能很好地拟合物性参数分界面及任意复杂形体模型，数据光滑程度好，数值模拟结果能真实、客观地反映雷达波的实际探测波形特征。

[1] 戴阳豪. 基于非结构网格的二维水流数值模拟研究[D]. 威海: 哈尔滨工业大学船舶与海洋工程学院, 2012: 9−10. DAI Yanghao. Two-dimensional numerical simulation of flows based on unstructured mesh[D]. Weihai: Harbin Engineering University. School of Naval Architecture and Ocean Engineering, 2012: 9−10.

[2] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994: 135−137. XU Shize. Finite element method for geophysics[M]. Beijing: Science Press, 1994: 135−137.

[3] 张晓蒙, 张鉴, 陆忠华. 基于OpenMP的二维非结构网格生成算法Delaunay并行实现[J]. 科研信息化技术与应用, 2012, 3(5): 20−28.Zhang Xiaomeng, Zhang Jian, LU Zhonghua. Parallelization of Delaunay 2D unstructured mesh generation algorithm based on OpenMP[J]. E-Science Technology & Application, 2012, 3(5): 20−28.

[4] Jacobsen D W, Gunzburger M, Ringler T J, et al. Parallel algorithms for planar and spherical Delaunay construction with an application to centroidal Voronoi tessellations[J]. Geosci Model Dev, 2013, 6(4): 1353−1365.

[5] 宋超, 关振群, 顾元宪. 二维自适应网格生成的改进AFT 与背景网格法[J]. 计算力学学报, 2005, 22(6): 694−699. SONG Chao, GUAN Zhenqun, GU Yuanxian. Improved AFT and background-mesh methods for two-dimensional adaptive mesh generation[J]. Chinese Journal of Computational Mechanics, 2005, 22(6): 694−699.

[6] Mcmotris H, Kallinderis Y. Octree-advancing front method for generation of unstructured surface and volume meshes[J]. AIAA Journal, 1997, 35(6): 976−954.

[7] Du Q, Ju L, Gunzburger M. Adaptive finite element methods for elliptic PDE’s based on conforming centroidal Voronoi Delaunay triangulations[J]. SIAM J Sci Comput, 2006, 28(6): 2023−2053.

[8] Key K, Weiss C. Adaptive finite-element modeling using unstructured grids: The 2D magnetotelluric example[J]. Geophysics, 2006, 71(6): G291−G299.

[9] LI Yuguo, Kerry K. 2D marine controlled-source electromagnetic modeling. Part 1: An adaptive finite-element algorithm[J]. Geophysics, 2007, 72(2): WA51−WA62.

[10] 汤井田, 王飞燕, 任政勇. 基于非机构化网格的2.5-D直流电阻率自适应有限元数值模拟[J]. 地球物理学报, 2010, 53(3): 708−716. TANG Jingtian, WANG Feiyan, REN Zhongyong. 2.5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese Journal of Geophysics, 2010, 53(3): 708−716.

[11] 任政勇, 汤井田. 基于局部加密非结构化网格的三维电阻率法有限元数值模拟[J]. 地球物理学报, 2009, 52(10): 2627−2634. REN Zhongyong, TANG Jingtian. Finite element modeling of 3D DC resistivity using locally refined unstructured meshes[J]. Chinese Journal of Geophysics, 2009, 52(10): 2627−2634.

[12] 柳建新, 孙丽影, 童孝忠, 等. 基于自适应有限元的起伏地形MT二维正演模拟[J]. 地球物理学进展, 2012, 27(5): 2016−2023. LIU Jianxin, SUN Liying, TONG Xiaozhong, et al. Undulating terrain 2D MT forward modeling with adaptive finite element[J]. Progress in Geophys, 2012, 27(5): 2016−2023.

[13] Lawson C L. Generation of a triangular grid with applications to contour plotting[C]//Technical Memorandum. California: Institute of Technology. Jet Pollution Laboratory, 1972: 299.

[14] Remacle J F, Henrotte F, Carrier-Baudouin T, et al. A frontal delaunay quad mesh generator using the∞norm[J]. International Journal for Numerical Methods in Engineering, 2013, 94(5): 494−512.

[15] ShewchukJ R. Delaunay refinement algorithms for triangular mesh generation[J]. Computational Geometry: Theory and Applications, 2002, 22(1/2/3): 21−74.

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

[17] 冯德山, 陈承申, 王洪华. 基于混合边界条件的有限元法GPR正演模拟[J]. 地球物理学报, 2012, 55(11): 3774−3785. FENG Deshan, CHEN Chengshen, WANG Honghua. Finite element method GPR forward simulation based on mixed boundary condition[J]. Chinese Journal of Geophysics, 2012, 55(11): 3774−3785.

[18] 陈承申. 探地雷达二维有限元正演模拟[D]. 长沙: 中南大学地球科学与信息物理学院, 2011: 13−15. CHEN Chengshen. Two-dimensional finite element forward simulation for ground penetrating radar model[D]. Changsha: Central South University. School of Geosciences and Info-Physics, 2011: 13−15.

[19] 戴前伟, 王洪华, 冯德山, 等. 基于三角形剖分的复杂GPR模型有限元法正演模拟[J]. 中南大学学报(自然科学版), 2012, 43(7): 2668−2673. DAI Qianwei, WANG Honghua, FENG Deshan, et al. 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.

[20] 戴前伟, 王洪华, 冯德山, 等. 基于三角形剖分的复杂GPR模型有限元正演模拟[J]. 中南大学学报(自然科学版), 2012, 43(7): 2668−2673. DAI Qianwei, WANG Honghua, FENG Deshan, et al. 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.

[21] 柳建新, 蒋鹏飞, 童孝忠, 等. 不完全LU分解预处理的BICGSTAB算法在大地电磁二维正演模拟中的应用[J]. 中南大学学报(自然科学版), 2009, 40(2): 484−491. LIU Jianxin, JIANG Pengfei, TONG Xiaozhong, et al. 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.

[22] 薛东川, 王尚旭, 焦淑静. 起伏地表复杂介质波动方程有限元数值模拟方法[J]. 地球物理学进展, 2007, 22(2): 522−529. XUE Dongchuan, WANG shangxu, JIAO Shujing. Wave equation finite-element modeling including rugged topography and complicated medium[J]. Progress in Geophysics, 2007, 22(2): 522−529.

(编辑 陈灿华)

GPR simulation by finite element method of unstructured grid based on Delaunay triangulation

DU Huakun, FENG Deshan, TANG Jingtian

(Key Laboratory of Hunan, Nonferrous Resources and Geologic Disasters Prospecting, School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)

The theories and steps of traditional Bowyer-Watson triangulation incremental insertionmethod were introduced, and the fixed boundary limitation and topological transformation technology such as Laplacian smoothing, edge contraction, edge division and point insertion were applied to mesh generation optimization. In order to make the numerical solution error approach to the uniform distribution in the whole domain, mesh gradation controls such as point source and line source were realized by interval function method combined with local coarsening or local refinement technology, and high quality Delaunay triangle mesh was built to realize self-adaption mesh generation. A complex GPR geoelectric model was built, and the self-adaption mesh generation was made about this model by Delaunay triangle. The results show that the unstructured mesh has good adaption to the geoelectric model whose physical parameter distribution is complex or whose geometric features are irregular. Unstructured triangle mesh generation has high quality and easy controlling mesh density to realize self-adaption finite element algorithm, which can improve the complex model forward precision.

Delaunay triangulation; unstructured grid; finite element method; forward simulation; ground penetrating radar

10.11817/j.issn.1672-7207.2015.04.022

P631

A

1672−7207(2015)04−1326−09

2014−04−15；

2014−06−24

国家自然科学基金资助项目(41074085)；教育部新世纪优秀人才支持计划项目(NCET-12-0551)；湖湘青年创新创业平台培养对象基金资助项目(2013)；中南大学升华育英计划项目(2012年)(Project (41074085) supported by the National Natural Science Foundation of China; Project (NCET-12-0551) supported by the Program for New Century Excellent Talents in University; Project (2013) supported by the Youth Innovation Platform of Hunan Raise Object Fund; Project (2012) supported by the Yuying Plan Project of Central South University)

冯德山，博士，教授，从事地球物理数据处理与正反演研究；E-mail：fengdeshan@126.com