復(fù)雜地形下高密度激電法2.5維有限單元法數(shù)值模擬

麻昌英，柳建新*1,，劉海飛1,，郭榮文1,

(1. 中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083；2. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083)

0 引言

激發(fā)極化法是勘探金屬、非金屬礦產(chǎn)等具有激發(fā)極化效應(yīng)礦產(chǎn)的主要物探方法之一，在工程物探中也有著重要的應(yīng)用。1971年Coggon[1]首先將有限元法應(yīng)用到電磁場(chǎng)模擬中，自此有限單元法在地球物理模擬計(jì)算中得到廣泛的應(yīng)用。在我國(guó)，徐世哲[2]對(duì)有限單元法進(jìn)行了深入研究。近年來(lái)，對(duì)于激發(fā)極化法正演模擬研究取得了較大的成果，主要是對(duì)激發(fā)極化法應(yīng)有優(yōu)先發(fā)進(jìn)行三維模擬。A.Dey等[3]對(duì)三維電阻率任意形狀模型展開(kāi)了研究；黃俊革等[4-8]對(duì)有限元法三維地電斷面電阻率正演模擬展開(kāi)了深入研究，實(shí)現(xiàn)了齊次條件下及介質(zhì)參數(shù)連續(xù)條件下三維電阻率有限元模擬；熊彬等[9-11]對(duì)復(fù)雜地形下電阻率/極化率有限單元法三維數(shù)值模擬進(jìn)行了研究，實(shí)現(xiàn)了復(fù)雜地形三維電阻率/極化率有限元模擬。

雖然電阻率/極化率三維有限元模擬取得了很好的成果，但用有限元法解決三維問(wèn)題是比較麻煩的，網(wǎng)格剖分量大，導(dǎo)致計(jì)算時(shí)間較慢，對(duì)計(jì)算機(jī)內(nèi)存要求較高。在2.5維問(wèn)題中，介質(zhì)為二維構(gòu)造電性特征，場(chǎng)源為三維場(chǎng)，通過(guò)傅里葉變換，將三維電場(chǎng)問(wèn)題轉(zhuǎn)化成帶參數(shù)的二維問(wèn)題求解，從而降低網(wǎng)格剖分量，節(jié)省計(jì)算時(shí)間。目前，2.5維有限元模擬在直流電阻率模擬中應(yīng)用已較成熟，徐世哲[2,12]對(duì)2.5維穩(wěn)定電流場(chǎng)問(wèn)題進(jìn)行了深入研究；阮百堯等[13]實(shí)現(xiàn)了三角元部分電導(dǎo)率分塊連續(xù)變化2.5維電場(chǎng)有限元模擬；湯井田等[14-15]實(shí)現(xiàn)了基于非結(jié)構(gòu)化網(wǎng)格的2.5維直流電阻率法自適應(yīng)有限元模擬。對(duì)于高密度激電法模擬，一條剖面上一般模擬點(diǎn)數(shù)較大，因此對(duì)模擬計(jì)算效率要求較高，2.5維模擬相對(duì)于三維模擬計(jì)算量要小很多，因而計(jì)算效率相對(duì)于三維模擬較明顯。本次研究在參考上述的基礎(chǔ)上，利用等效電阻率法實(shí)現(xiàn)高密度激電法2.5維有限元數(shù)值模擬。

1 點(diǎn)源場(chǎng)的變分問(wèn)題

點(diǎn)源二維穩(wěn)定電流場(chǎng)的邊值問(wèn)題

(1)

采用有限元法求解點(diǎn)源2.5D穩(wěn)定電流場(chǎng)的邊值問(wèn)題，首先需將方程(1)轉(zhuǎn)化為等價(jià)的變分問(wèn)題

(2)

然后再利用有限元法解變分問(wèn)題式(2)。

2 網(wǎng)格剖分

對(duì)傳統(tǒng)的二維矩形剖分，其不適合于對(duì)復(fù)雜地形的模擬，對(duì)于非結(jié)構(gòu)化網(wǎng)格，其雖然可以很好地模擬復(fù)雜地形，正演計(jì)算精度也較高，但其網(wǎng)格不規(guī)則，剖分難度較大，計(jì)算量較大，不易應(yīng)用于反演計(jì)算。采用結(jié)構(gòu)化三角形網(wǎng)格，可以較好地解決網(wǎng)格剖分復(fù)雜度及復(fù)雜地形模擬問(wèn)題，在高密度激電法中，采用規(guī)則的網(wǎng)格剖分有利于模擬高密度激電法的滾動(dòng)測(cè)量。在縱向方向上，從上到下網(wǎng)格剖分間距由小到大，在保證計(jì)算精度的情況下，減少網(wǎng)格剖分?jǐn)?shù)，可節(jié)省計(jì)算時(shí)間。

在具體實(shí)現(xiàn)過(guò)程中，根據(jù)同一條測(cè)線上所有測(cè)點(diǎn)的相對(duì)位置構(gòu)成橫向上的網(wǎng)格剖分，并根據(jù)測(cè)深點(diǎn)所能達(dá)到的勘探深度，確定研究區(qū)域的縱向網(wǎng)格剖分，再根據(jù)每個(gè)測(cè)深點(diǎn)的供電和測(cè)量極距關(guān)系進(jìn)行網(wǎng)格剖分，構(gòu)成正演模擬網(wǎng)格(如圖1中的粗實(shí)線網(wǎng)格)。為保證正演模擬的精度，在網(wǎng)格的左、右和下邊緣做一定程度的網(wǎng)格外延，如圖1所示。

圖1 網(wǎng)格剖分示意圖Fig．1 Mesh sketch map

3 有限單元法

將積分區(qū)域剖分成許多三角單元，總節(jié)點(diǎn)數(shù)為N(圖1)，則式(2)中對(duì)區(qū)域Ω和邊界?！薜姆e分可分界為對(duì)各三角單元e和Γe的積分之和則有：

K0(λr)dΓ

(3)

假設(shè)三角單元內(nèi)轉(zhuǎn)換電位V和電導(dǎo)率σ線性變化，即在每個(gè)三角單元有：

(4)

式(4)中，下標(biāo)1、2和3是三角單元的三角節(jié)點(diǎn)號(hào)；σ=(σ1,σ2,σ3)T為節(jié)點(diǎn)上的電導(dǎo)率；V=(V1,V2,V3)T為三角單元節(jié)點(diǎn)上的波數(shù)域電位；N=(N1,N2,N3)T，Nl=(alx+blz+cl)/2Δ為x和z的線性函數(shù)(l=1,2,3)，其中Δ=(a1b2-a2b1)為三角單元的面積，其中a1=z2-z3、a2=z3-z1、a3=z1-z2；b1=x3-x2、b2=x1-x3、b3=x2-x1；c1=x2z3-x3z2、c2=x3z1-x1z3、c3=x1z2-x2z1。

根據(jù)上面所述的網(wǎng)格剖分，根據(jù)參考文獻(xiàn)[2]和[13]，得到變分問(wèn)題式(2)的泛函F(V)的數(shù)值表達(dá)式：

(5)

式中K是由全部三角單元和邊界單元的(Ke1+Ke2+Ke3)相加組成的N×N階對(duì)稱帶寬矩陣，V是由所有N個(gè)三角網(wǎng)格節(jié)點(diǎn)上的轉(zhuǎn)換電位組成的列矢量，S是與源有關(guān)的列矢量。式(5)泛函F(V)對(duì)V求變分，并令其為零，就得到求波數(shù)域中各節(jié)點(diǎn)上電位的線性方程組:

KV=S

(6)

通過(guò)變帶寬存儲(chǔ)喬里斯基法解線性方程組(6)，便得各節(jié)點(diǎn)的波數(shù)域電位V。

對(duì)波數(shù)域電位V進(jìn)行付氏逆變換，付氏逆變換公式為

(7)

4 模型計(jì)算

為了驗(yàn)證程序的正確性，計(jì)算均勻半空間ρ0=100 Ω·m的電阻率測(cè)深曲線。如圖2所示，實(shí)點(diǎn)為本文程序模擬結(jié)果，實(shí)線為解析解，從圖2中可看出，除了前三個(gè)點(diǎn)由于電源點(diǎn)影響外，其余各點(diǎn)擬合很好。

圖2 均勻半空間測(cè)深曲線Fig．2 Sounding curve of homogeneous half space

1)模型一。水平均勻半空間下，圍巖電阻率ρ0=100 Ω·m，極化率η0=1%，存在一個(gè)電阻率ρ1=100 Ω·m，極化率η2=10%，沿Y方向無(wú)限延伸的4 m×4 m的正方形極化異常體，其與圍巖極化率在接觸帶處逐漸變化(線性)。異常體頂部埋深6 m，用對(duì)稱四極裝置進(jìn)行斷面測(cè)量，點(diǎn)距2 m，極化率正演結(jié)果如圖3所示。

圖3 模型一：正演視極化率擬斷面圖Fig．3 Model one: polarization forward profile

2)模型二。山脊地形起伏高度4 m，左右各延伸6 m，圍巖電阻率ρ0=100 Ω·m，極化率η0=1%，在山脊正下方存在一個(gè)電阻率ρ1=100 Ω·m，極化率η2=10%，沿Y方向無(wú)限延伸的4 m×4 m的正方形極化異常體，其與圍巖極化率在接觸帶處逐漸變化(線性)。異常體頂部距離水平地面6 m，用對(duì)稱四極裝置進(jìn)行斷面測(cè)量，點(diǎn)距2 m，極化率正演結(jié)果如圖4所示。

圖4 模型二正演視極化率擬斷面圖Fig．4 Model two: polarization forward profile

3)模型三。山谷地形起伏高度4 m，左右各延伸6 m，圍巖電阻率ρ0=100 Ω·m，極化率η0=1%，在山谷正下方存在一個(gè)電阻率ρ1=100 Ω·m，極化率η2=10%，沿Y方向無(wú)限延伸的4 m×4 m的正方形極化異常體，其與圍巖極化率在接觸帶處逐漸變化(線性)。異常體頂部距離水平地面6 m，用對(duì)稱四極裝置進(jìn)行斷面測(cè)量，點(diǎn)距2 m，極化率正演結(jié)果如圖5所示。

圖5 模型三正演視極化率擬斷面圖Fig．5 Model three: polarization forward profile

4)模型四。山脊-山谷地形，起伏高度均為4 m，左右各延伸均為6 m，山脊與山谷之間為4 m水平地形連接，在此段水平地形正下方存在一個(gè)電阻率ρ1=100 Ω·m，極化率η2=10%，沿Y方向無(wú)限延伸的4 m×4 m的正方形極化異常體，其與圍巖極化率在接觸帶處逐漸變化(線性)，圍巖電阻率ρ0=100 Ω·m，極化率η0=1%。異常體頂部埋深6 m，用對(duì)稱四極裝置進(jìn)行斷面測(cè)量，點(diǎn)距2 m，極化率正演結(jié)果如圖6所示。

圖6 模型四正演視極化率擬斷面圖Fig．6 Model four: polarization forward profile

5 結(jié)論

在實(shí)際地質(zhì)環(huán)境中，巖石、礦物的物性參數(shù)多為變化的，本研究基于矩形-三角形網(wǎng)格，設(shè)計(jì)單元內(nèi)極化率及電阻率為雙線性變化，實(shí)現(xiàn)了連續(xù)介質(zhì)的激發(fā)極化法2.5維有限元數(shù)值模擬。計(jì)算了激發(fā)極化模型，通過(guò)三角元模擬了地形的起伏對(duì)激發(fā)極化法的影響，從上述模型計(jì)算結(jié)果可以分析，當(dāng)?shù)匦纹鸱捌鸱匦蜗麓嬖跇O化異常體時(shí)，山脊地形與山谷山谷都影響視極化率的觀測(cè)，從模型二與模型三模擬結(jié)果分析比較可知，山谷地形影響較大，山脊地形影響較小。在非結(jié)構(gòu)網(wǎng)格系統(tǒng)中，網(wǎng)格的的不規(guī)則性質(zhì)導(dǎo)致其較難應(yīng)用于反演計(jì)算，而在矩形-三角形網(wǎng)格系統(tǒng)中，利用三角元進(jìn)行正演，可以模擬復(fù)雜的地形起伏，當(dāng)應(yīng)用于反演計(jì)算時(shí)，可以將其矩形網(wǎng)格用于反演，使得反演問(wèn)題相對(duì)簡(jiǎn)化。

參考文獻(xiàn)：

[1] COGGON J H. Electromagnetic and eletrical modeling by the finite element method[J].Geophysics, 1971, 36(1): 132-155.

[2] 徐世浙. 地球物理中的有限單元法[M]. 北京：科學(xué)出版社, 1994.

[3] DEY A, MORRISION H F. Resistivity modeling for arbitrarily shaped 3-D structures[J].Geophysics, 1979, 44( 4): 753-780.

[4] 阮百堯, 黃俊革, 鮑光淑. 齊次邊界條件下三維地電斷面電阻率有限元數(shù)值模擬法[J].桂林工學(xué)院學(xué)報(bào), 2002,22: 11- 14.

[5] 黃俊革, 阮百堯,鮑光淑. 基于有限單元法的三維地電斷面電阻率反演[J]. 中南大學(xué)學(xué)報(bào)：自然科學(xué)版, 2004, 35(2): 295-299.

[6] 阮百堯, 熊彬, 徐世浙. 三維地電斷面電阻率測(cè)深有限元數(shù)值模擬[J]. 地球科學(xué)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào), 2001, 26(1): 73-77.

[7] 阮百堯, 熊彬. 電導(dǎo)率連續(xù)變化的三維電阻率測(cè)深有限元模擬[J]. 地球物理學(xué)報(bào),2002, 45(1): 131-138.

[8] 黃俊革. 三維電阻率/極化率有限元正演模擬與反演成像[D]. 長(zhǎng)沙：中南大學(xué),2003.

[9] 熊彬, 阮百堯, 羅延鐘. 復(fù)雜地形條件下直流電阻率異常三維數(shù)值模擬研究[J]. 地質(zhì)與勘探, 2003, 30(4): 60-64.

[10] 林家勇, 湯井田, 丁茂斌. 復(fù)雜地形條件下激發(fā)極化有限單元法三維數(shù)值模擬[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版, 2010, 40(5): 1183-1187.

[11] 林家勇, 湯井田, 李開(kāi)壁，等. 雙頻激電法有限元數(shù)值模擬研究[J]. 物探化探計(jì)算技術(shù), 2009, 31(3): 189-192.

[12] 徐世浙. 點(diǎn)電源二維電場(chǎng)問(wèn)題中付氏變換的波數(shù)k的選擇[J]. 物探化探計(jì)算技術(shù), 1988, 10(3): 235-239.

[13] 阮百堯. 三角單元部分電導(dǎo)率分塊連續(xù)變化點(diǎn)源二維電場(chǎng)有限元數(shù)值模擬[J]. 廣西科學(xué), 2001, 30(5): 1-3.

[14] 湯井田, 王飛燕. 基于非結(jié)構(gòu)化網(wǎng)格的2.5D直流電阻率模擬[J]. 物探化探計(jì)算技術(shù), 2008, 8(1): 413-418.

[15] 王飛燕. 基于非結(jié)構(gòu)化網(wǎng)格的2.5-D直流電阻率法自適應(yīng)有限元數(shù)值模擬[D]. 長(zhǎng)沙：中南大學(xué), 2009.