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

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

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

0 引言

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

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

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

1 2.5維點源場基本方程

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

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

(1)

(2)

(3)

筆者欲采用自然單元法求解該邊值問題，需先將其變換為變分問題?？杀磉_(dá)為式(4)[17]。

(4)

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

(5)

(6)

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

KU=P

(7)

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

(8)

其中采用的波數(shù)及其權(quán)值系數(shù)見表1[17]。

表1 波數(shù)k及其權(quán)值系數(shù)g

2 自然單元法

數(shù)值模擬算法的本質(zhì)區(qū)別是構(gòu)造形函數(shù)的方式不同，其影響數(shù)值模擬過程的計算效率、精度、復(fù)雜度等。自然單元法作為一種無網(wǎng)格法，僅需節(jié)點信息即可構(gòu)造插值形函數(shù)，且整個過程可利用復(fù)合函數(shù)鏈?zhǔn)角髮?dǎo)法則簡化為系列加減乘除計算，能有效保證計算效率。

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

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

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

(9)

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

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

(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為外接圓圓心坐標(biāo)；vx，x、vx，y、vy，x、vy，y為圓心坐標(biāo)對坐標(biāo)軸的偏導(dǎo)；α和D為中間變量；Dx和Dy表示D對坐標(biāo)軸的偏導(dǎo)。

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

(20)

進一步地，插值形函數(shù)的偏導(dǎo)可表示為：

(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]，即將邊界處相鄰的兩個自然節(jié)點看作網(wǎng)格的無窮遠(yuǎn)邊界。最后由公式(8)計算得到二維模型空間的電位分布，進而計算地下視電阻率分布情況。

3 數(shù)值算例

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

3.1 水平層狀模型

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

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

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

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

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

圖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個自然節(jié)點均勻填充。采用二極AM裝置進行測量，點源A位于(0, 0)處且固定不動。測定M從-10 m移動到10 m處，電極間距為2 m。

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

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

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

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

圖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 地質(zhì)體模型

圖8所示地電模型中賦存兩個電阻率異常地質(zhì)體。模型尺度為80 m×40 m，由81×41個自然節(jié)點均勻填充。背景電阻率為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 結(jié)論

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

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

2) 自然單元法在構(gòu)造形函數(shù)時只需要節(jié)點信息，且基于復(fù)合函數(shù)的鏈?zhǔn)角髮?dǎo)法則使得計算過程簡單明了易編程實現(xiàn)。

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