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

(1. 中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083；2. 中南大學 地球科學與信息物理學院，長沙 410083)

0 引言

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

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

1 點源場的變分問題

點源二維穩(wěn)定電流場的邊值問題

(1)

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

(2)

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

2 網(wǎng)格剖分

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

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

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

3 有限單元法

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

K0(λr)dΓ

(3)

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

(4)

式(4)中，下標1、2和3是三角單元的三角節(jié)點號；σ=(σ1,σ2,σ3)T為節(jié)點上的電導率；V=(V1,V2,V3)T為三角單元節(jié)點上的波數(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ù)參考文獻[2]和[13]，得到變分問題式(2)的泛函F(V)的數(shù)值表達式：

(5)

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

KV=S

(6)

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

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

(7)

4 模型計算

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

參考文獻：

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

[2] 徐世浙. 地球物理中的有限單元法[M]. 北京：科學出版社, 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].桂林工學院學報, 2002,22: 11- 14.

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

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

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

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

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

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

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

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

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

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

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