基于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