基于自然單元法的2.5維直流電阻率數值模擬

王佳新, 崔益安, 謝 靜, 羅議建, 柳建新

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.中南大學 有色資源與地質災害探查湖南省重點實驗室，長沙 410083；3.中南大學 教育部有色金屬成礦預測重點實驗室，長沙 410083)

0 引言

直流電阻率法數值模擬研究中，常用的數值算法包括有限單元法[1-3]及有限差分法[4]等網格法，但其在處理復雜、不規則地電模型過程中嚴格受網格約束，需要開展大量繁瑣的前期工作，影響計算效率、效果[5]。針對復雜模型適應性問題，麻昌英等[6]采用全局弱式無單元法對復雜地電條件下的直流電阻率模型開展數值模擬。該方法雖擺脫了網格的限制，但形函數構造復雜、計算量大、邊界條件施加困難。

自然單元法是一種新型的無網格數值方法，最早由Braun等提出[7]。該方法基于Voronoi圖[8]和Delaunay三角剖分[9]，采用Sibson[10]或non-Sibson[11]自然鄰點插值方法構造試函數，求解偏微分方程。自然單元法不僅具有無網格法和網格法的優點，還克服了它們的缺陷。與無網格法相比，自然單元法的計算量小，計算效率高；與網格法相比，自然單元法不受網格限制，只需要自然節點信息，這對復雜模型的剖分具有重要意義[5]。目前自然單元法主要應用于材料力學領域[12-14]，而很少見其在勘探地球物理領域中的應用報道[5, 15-16]。

這里首次成功將自然單元法應用于直流電阻率法數值模擬研究中，探究該算法在該領域的可行性及模擬效果。首先從2.5維點源場邊值問題出發，推導其滿足的變分及積分形式基本方程。隨后推導二維情況下的自然單元法形函數及其導數的基本方程。最后通過一系列數值算例驗證該算法的有效性，證明其在直流電阻率模型中的適用性。

1 2.5維點源場基本方程

對于點源位于地表的二維地電模型，其滿足的2.5維穩定直流場邊值問題可表示為[17]：

▽·(σ▽U)-k2σU=-I0δ(A)

(1)

(2)

(3)

筆者欲采用自然單元法求解該邊值問題，需先將其變換為變分問題。可表達為式(4)[17]。

(4)

在任意積分網格中，式(4)中的各項積分可表示為[17]：

(5)

(6)

源積分項及邊界積分項詳見文獻[17]。最后將各積分按節點編號累加至總剛度矩陣中，由變分為0可得到線性方程組：

KU=P

(7)

其中:K為總剛度矩陣;U為待求波數域電位向量；P為源向量。求解公式(7)即可得到各節點波數域電位值，而后再通過傅里葉反變換求解得到二維空間中任意一點的電位值：

(8)

其中采用的波數及其權值系數見表1[17]。

表1 波數k及其權值系數g

2 自然單元法

數值模擬算法的本質區別是構造形函數的方式不同，其影響數值模擬過程的計算效率、精度、復雜度等。自然單元法作為一種無網格法，僅需節點信息即可構造插值形函數，且整個過程可利用復合函數鏈式求導法則簡化為系列加減乘除計算，能有效保證計算效率。

Voronoi結構和Delaunay三角剖分是由一組離散點定義的最基本的幾何結構。如圖1所示[5]，假設在平面內給定7個離散節點，作任意兩節點的垂直平分線(圖1(a))，各半平面的交集構成多個凸多邊形區域，即為一階Voronoi單胞結構(圖1(b))；進而將所有具有共同邊界的Voronoi單胞對應的節點連接成三角形網絡以構成Delaunay三角剖分(圖1(c))。

圖1 Voronoi圖和Delaunay三角化Fig.1 Voronoi,Delaunay triangularization(a)7節點系統的Voronoi剖分；(b)7節點系統的Voronoi單元；(c)7節點系統的Delaunay三角剖分

基于Voronoi結構和Delaunay三角剖分，可以搜索各積分點的自然鄰點并構造插值函數。筆者選用簡便易行的non-Sibson (也稱為Laplace)自然鄰點插值格式[5, 18]。自然單元法形函數構造過程如圖2所示。插值形函數可寫作式(9)。

(9)

以圖2(b)中的積分點X為例，其自然鄰點為節點1、2、3、4。節點12、23、34、41分別為三角形1-2-X、2-3-X、3-4-X、4-1-X的外接圓圓心。對于任意三角形ABC，令其頂點坐標分別為(x1,y1)、(x2,y2)、(x3,y3)，其中C(x3,y3)與積分點X重合，則該三角形的外接圓圓心坐標可表示為：

圖2 Non-sibson插值法Fig.2 Non-sibson interpolation(a)假設點X為7節點系統的某積分點；(b)插值點X的自然鄰點及non-Sibson插值參數

(10)

(11)

(12)

(13)

(14)

(15)

其中：

D=2[(x1-x3)(y2-y3)-

(x2-x3)(y1-y3)]

(16)

α=(x2+x1)(x2-x1)+

(y2+y1)(y2-y1)

(17)

D,x=2(y1-y2)

(18)

D,y=2(x2-x1)

(19)

其中:vx和vy為外接圓圓心坐標；vx，x、vx，y、vy，x、vy，y為圓心坐標對坐標軸的偏導；α和D為中間變量；Dx和Dy表示D對坐標軸的偏導。

設任意Voronoi邊的兩個端點坐標為V1(vx1,vy1)和V2(vx2,vy2)，則si(x,y)可表達為：

(20)

進一步地，插值形函數的偏導可表示為：

(21)

(22)

其中：

(23)

(24)

si,x(x,y)=[(vx1-vx2)(vx1,x-vx2,x)+

(vy1-vy2)(vy1,x-vy2,x)]/

(25)

si,y(x,y)=[(vx1-vx2)(vx1,y-vx2,y)+

(vy1-vy2)(vy1,y-vy2,y)]/

(26)

(27)

(28)

(29)

將公式(9)～公式(29)帶入公式(5)、公式(6)，即可求解方程組(7)。其中，邊界條件的加載可類比點源二維電場的有限單元法處理過程[17]，即將邊界處相鄰的兩個自然節點看作網格的無窮遠邊界。最后由公式(8)計算得到二維模型空間的電位分布，進而計算地下視電阻率分布情況。

3 數值算例

為驗證自然單元法在2.5維穩定直流場數值模擬中的有效性，筆者對一系列直流電阻率地電模型進行算例分析。

3.1 水平層狀模型

采用水平層狀模型進行解析解驗證。其中上層電阻率為100 Ω·m，厚度為10 m，下層電阻率為10 Ω·m。采用固定點源的二極裝置進行測量，供電電極A位于(0, 0)處，測量電極M從(0, 1)逐漸移動至(0, 80)。161×161個自然節點均勻布設在二維地電模型中。視電阻率曲線如圖3所示。

圖3 水平兩層介質視電阻率曲線Fig.3 Apparent resistivity curves of a two-layered model

從數值結果可以看出：對于測區[0, 20]，數值計算結果與解析解吻合很好，而隨著測點M逐漸靠近邊界，自然單元法的計算精度呈現一定的誤差 (相對誤差在1%～10%之間)，但與常規有限單元法的計算結果[1]相比滿足精度需求。

3.2 均勻半空間中賦存塊狀異常體

均勻半空間中賦存3個塊狀電阻率異常體的驗證模型(圖4)。Xu等[19]采用該模型驗證了邊界元法，肖曉等[1]采用該模型驗證了有限元-無限元耦合法。筆者采用自然單元法計算該模型，并與參考文獻[1]及參考文獻[19]的計算結果相比較，以驗證自然單元法的有效性。

圖4 均勻半空間中賦存3個塊狀異常體Fig.4 Three blocks buried in homogeneous half-space

如圖4所示，3個異常體的電阻率分別為ρ1=0.1 Ω·m、ρ2=10 Ω·m、ρ3=50 Ω·m，均勻半空間電阻率為ρ4=1 Ω·m。整個異常體的邊長為1 m，頂部埋深為1 m。模型空間由321×321個自然節點均勻填充。采用二極AM裝置進行測量，點源A位于(0, 0)處且固定不動。測定M從-10 m移動到10 m處，電極間距為2 m。

圖5為邊界元法[19]、有限元法[1,19]以及自然單元法的計算結果。可見自然單元法的計算精度與邊界元法、有限單元法等相當，進一步驗證了算法的有效性和計算精度。

圖5 不同方法計算的電位對比Fig.5 A comparison of potential calculated with different methods

3.3 均勻半空間中賦存三棱柱體

采用均勻半空間中賦存三棱柱的地電模型，驗證自然單元法對復雜幾何體的有效剖分。模型如圖6所示，模型空間共布設61×31個自然節點，其中不規則目標體由自然節點靈活填充。三棱柱體電阻率為10 Ω·m，背景電阻率為100 Ω·m。三棱柱體的三角形面為高度等于10 m的等腰三角形，頂部埋深為10 m。采用對稱四極裝置進行測量。視電阻率斷面如圖7所示。數值結果表明，自然單元法能有效剖分復雜地電模型，且能夠較好地反映地下異常體的分布特征。

圖6 均勻半空間賦存三棱柱模型示意圖Fig.6 Infinite horizontal tri-prism buried in homogeneous half-space

圖7 均勻半空間中賦存三棱柱體視電阻率斷面圖Fig.7 Section of apparent resistivity for a triangular prism buried in homogeneous half-space

3.4 地質體模型

圖8所示地電模型中賦存兩個電阻率異常地質體。模型尺度為80 m×40 m，由81×41個自然節點均勻填充。背景電阻率為100 Ω·m；高阻體尺度為5 m×5 m，電阻率為1 000 Ω·m；低阻體尺度為5 m×5 m，電阻率為10 Ω·m。兩異常體相距10 m，且頂部埋深均為1 m。采用對稱四極裝置進行測量。圖9為視電阻率斷面，可見自然單元法能較好反映電阻率異常體的分布特征。

圖8 均勻半空間中賦存2個塊狀異常體Fig.8 Two blocks buried in homogeneous half-space

圖9 均勻半空間中賦存2個塊狀異常體視電阻率斷面圖Fig.9 Section of apparent resistivity for two blocks buried in homogeneous half-space

4 結論

首次將自然單元法引入直流電阻率模型數值模擬研究中，給出了詳細的公式推導，并通過4個地電模型驗證該算法的可行性，得到以下結論：

1) 自然單元法在直流電阻率法數值模擬領域是有效的，且計算精度滿足需求。

2) 自然單元法在構造形函數時只需要節點信息，且基于復合函數的鏈式求導法則使得計算過程簡單明了易編程實現。

3) 自然節點布設靈活，能較好地填充復雜異常體，適用于復雜模型數值模擬研究。