柯頌頌，宋 滔，劉 云

(1.中國科學(xué)院 地球化學(xué)研究所礦床地球化學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 貴陽 550002；2.中國科學(xué)院大學(xué)，北京 100049)

0 引言

直流電阻率法各向同性的正反演技術(shù)已經(jīng)比較成熟，并且在工程、找礦等領(lǐng)域有了廣泛應(yīng)用[1-2]。隨著數(shù)值模擬技術(shù)的發(fā)展，研究的熱點(diǎn)聚集到了更符合實(shí)際情況的連續(xù)介質(zhì)和各向異性介質(zhì)。

對連續(xù)介質(zhì)的研究，徐世浙[3]使用有限單元法矩形單元剖分;阮百堯[4]使用三角單元剖分，實(shí)現(xiàn)了對連續(xù)介質(zhì)的數(shù)值模擬，取得較高精度;劉云[5]在阮百堯的基礎(chǔ)上，使用矩形內(nèi)剖分四個三角形的剖分方式實(shí)現(xiàn)了對連續(xù)介質(zhì)、復(fù)雜地形以及復(fù)雜模型的數(shù)值模擬。

近年來，對各向異性直流電法的研究也越來越多[6-9]。但關(guān)于直流電阻率法2.5維各向異性正演模擬的研究相對較少。徐世浙等[3]使用有限單元法對二維各向異性的直流電阻率法進(jìn)行了模擬;Zhou等[9]使用高斯正交網(wǎng)格(Gaussian quadrature grids)實(shí)現(xiàn)了對2.5維復(fù)雜各向異性介質(zhì)的數(shù)值模擬，取得較好的精度。由于三維各向異性參數(shù)太多，導(dǎo)致基于三維各向異性正演的反演研究工作進(jìn)展緩慢，所以研究直流電阻率法2.5維的正反演方法成為探索各向異性反演工作的橋梁。在前人的研究中2.5維正演均是基于總場法的，因?yàn)楫惓龇ㄐ枰蠼庖淮螆?，而電性各向異性介質(zhì)點(diǎn)源傅氏空間中的解析解較難求得，所以關(guān)于2.5維各向異性介質(zhì)異常場法的研究較少。

這里給出各向異性介質(zhì)2.5維總場和異常場的變分問題。在實(shí)現(xiàn)異常場法時，因?yàn)閷Ⅻc(diǎn)源各向異性介質(zhì)空間電位轉(zhuǎn)換到傅里葉空間中具有一定困難，所以筆者進(jìn)行了簡化處理，假設(shè)點(diǎn)源附近的介質(zhì)為主軸各向異性，從而實(shí)現(xiàn)2.5維各向異性介質(zhì)異常場法的模擬，通過算例證明該方法的正確性。

1 變分問題

一般定義各向異性電阻率為張量形式[11]，如式(1)所示。

(1)

選取如圖1所示的觀測坐標(biāo)系，z方向?yàn)榇怪狈较?，x、y方向?yàn)樗椒较?，假設(shè)介質(zhì)構(gòu)造為x方向，即沿x方向介質(zhì)沒有變化，設(shè)介質(zhì)電性主軸的平面x′y′與坐標(biāo)軸xy平面的夾角為α，此時電阻率張量表達(dá)式(1)中旋轉(zhuǎn)矩陣D為：

(2)

圖1 二維各向異性

將式(2)代入式(1)中得到介質(zhì)的電阻率張量為：

(3)

相應(yīng)的電導(dǎo)率為：

(4)

(5)

根據(jù)Zhou[9]、嚴(yán)波[11]等的推導(dǎo)，傅里葉空間中的電位U滿足的邊值問題如式(6)所示。

(6)

對應(yīng)的變分問題如式(7)所示。

(7)

當(dāng)采用異常場法模擬時，將總場分為背景場(一次場)和異常場(二次場)，以消除源的影響，提高模擬精度。

與各向同性的邊值問題和變分問題類似，我們給出異常場滿足的變分問題如式(8)所示。

(8)

2 有限單元法

2.1 區(qū)域剖分

首先將整個區(qū)域剖分成矩形單元，然后再將每個矩形剖分成兩個三角形，如圖 2所示。

圖2 區(qū)域剖分

在三角單元內(nèi)假設(shè)電位是線性變化的，在單元內(nèi)任意位置的電位u可以通過形函數(shù)和三角形三個節(jié)點(diǎn)的電位表示，如式(9)所示。

u=Niui+Njuj+Nmum=NTu=uTN

(9)

其中:NT=(Ni,Nj,Nm)為形函數(shù)；uT=(ui,uj,um)為三角形節(jié)點(diǎn)的電位。

其中：Δ為三角形的面積，且有：

2.2 單元分析

將式(8)中的積分在區(qū)域離散化，表示成所有單元的線性組合如式(10)所示。

(10)

式(10)中的積分項(xiàng)依次記為積分1、2、3、4、5和6。

根據(jù)嚴(yán)波[11]等的推導(dǎo)，以及積分1和積分4的相似性，得到：

(11)

(12)

單元積分2和積分5為：

(13)

(14)

單元積分3和積分6為：

(15)

(16)

將單元矩陣添加到整體系數(shù)矩陣中的相應(yīng)位置，得到式(17)。

(17)

令式(17)的變分為0，得到線性方程組(18)[3]。

(18)

解線性方程組得到波數(shù)域中異常場的電位，進(jìn)行傅里葉反變換得到空間場中的異常場電位，最后加上一次場電位得到總場電位。

2.3 歐拉角的簡化

觀察式(17)和式(18)，要得到方程組還需計(jì)算波數(shù)域中的一次場電位U0，點(diǎn)源均勻半空間各向異性介質(zhì)電位表達(dá)式為式(19)[12]。

(19)

其中B=(r-r0)T·ρ·(r-r0)，將B展開后，直接使用該表達(dá)式進(jìn)行傅里葉變換較為困難，因此筆者采用簡化歐拉角的方法進(jìn)行處理。假設(shè)二維構(gòu)造下點(diǎn)源附近的介質(zhì)(背景介質(zhì))電性主軸與觀測坐標(biāo)系的夾角為零得到：

ρ=diag(ρx,ρy,ρz)

(20)

對式(20)進(jìn)行傅里葉變換，得到傅氏空間中電位表達(dá)式為式(21)。

(21)

得到波數(shù)域中的一次場后，代入有限元公式進(jìn)行計(jì)算，將角度信息也作為異常來處理，由有限元完成計(jì)算異常場的工作。

此處我們假設(shè)了背景電阻率的電性主軸與觀測坐標(biāo)系的夾角α為零，所以該方法對α≠0的模型并不適合。

圖3 模型1：兩層含VTI介質(zhì)模型

3 算例分析

3.1 兩層含垂直各向異性(VTI)介質(zhì)的地電模型

設(shè)計(jì)一個兩層單界面模型，其中第一層為VTI介質(zhì)，第二層為各向同性介質(zhì)，模型和電極布置如圖3所示。發(fā)射電極為A點(diǎn)，模擬供入1A的電流，1～10為接收電極且電極之間距離為1 m。

3.1.1 算法驗(yàn)證

分別使用總場法和異常場法進(jìn)行模擬，選用的波數(shù)為徐世浙[3]計(jì)算的最優(yōu)波數(shù)。異常場法中背景場設(shè)為ρx=ρy=0.5 Ω·m，ρz=2.0 Ω·m產(chǎn)生的電場，也即點(diǎn)源處的電阻率作為背景電阻率。將總場法與異常場法的數(shù)值模擬結(jié)果與解析解分別進(jìn)行對比，如表1所示。

從模擬結(jié)果可以看出，使用異常場法精度較高，誤差均在1%以內(nèi)，而總場法在點(diǎn)源附近誤差較大，并且整體的誤差也較異常場法較高,也證明了本文算法的正確性和準(zhǔn)確性。

3.1.2 波數(shù)的討論

在以上計(jì)算中，如果采用宋滔等[14]計(jì)算的7點(diǎn)波數(shù)，總場法和異常場法的計(jì)算結(jié)果分別與解析解對比，相對誤差如表2所示。

從結(jié)果中可以看出，異常場法模擬的結(jié)果誤差非常小，但是總場法采用7點(diǎn)波數(shù)的模擬結(jié)果與解析解的相對誤差較大。這是因?yàn)樵诓〝?shù)域中點(diǎn)源附近有限元解誤差較大，所以變換到空間域時誤差也較大；7點(diǎn)波數(shù)本來的精度是較高的，即如果有限元解越精確，得到空間域的電位也越準(zhǔn)確，相反如果解的誤差較大，使用七點(diǎn)波數(shù)反而會使空間域的電位誤差變大。

因?yàn)樵诘匦螚l件下無法使用異常場法，所以在各向異性的正演模擬中用總場法時采用5點(diǎn)波數(shù)；采用異常場法時選用7點(diǎn)波數(shù)，可以取得相對更高的精度。

表1 數(shù)值解與解析解對比

表2 采用7點(diǎn)波數(shù)的模擬結(jié)果對比

表3 TTI介質(zhì)模擬結(jié)果對比

圖4 解析解與異常場法數(shù)值解曲線

圖5 解析解與異常場法數(shù)值解相對誤差

圖6 含異常體模型

3.1.3 背景電阻率的取值

采用異常場法簡化歐拉角時，假設(shè)點(diǎn)源附近的介質(zhì)與選定坐標(biāo)系的夾角為零，通過模型來計(jì)算當(dāng)點(diǎn)源處介質(zhì)為TTI時異常場的結(jié)果。

假設(shè)均勻半空間，電阻率為ρx=ρy=0.5 Ω·m，ρz=2.0 Ω·m，采用異常場法模擬α=0°、30°、60°、90°的結(jié)果，與解析解進(jìn)行對比驗(yàn)證文中的假設(shè)；設(shè)置電極為距離點(diǎn)源1 m～10 m且電極距為1 m，10個測量點(diǎn)。

表3給出了α=30°時，總場法與異常場法的解，以及它們與解析解的相對誤差，從表3中可以看出，總場法在臨近點(diǎn)源的第二個點(diǎn)的誤差較大，達(dá)到了2.571 1 %，在其他測點(diǎn)處的誤差均在1%以內(nèi)；但異常場法的解誤差均較大，在源附近已經(jīng)達(dá)到了63.012 1 %。因?yàn)榇藭r背景電阻率的電性主軸與坐標(biāo)系有α=30°的夾角，不滿足文中關(guān)于歐拉角的假設(shè)，所以此處采用異常場法處理的結(jié)果不準(zhǔn)確。

設(shè)α=0°、30°、60°、90°，解析解與異常場法求解結(jié)果的曲線見圖4，以及相對誤差的曲線見圖5。

從以上模擬可以發(fā)現(xiàn)點(diǎn)源處的介質(zhì)如果為TTI介質(zhì)(背景介質(zhì))，只能采用總場法進(jìn)行正演，文中給出的異常場法不再適用。

3.2 含異常體模型

設(shè)計(jì)如圖6所示的含異常體模型，異常體距離地面3 m，大小為3 m×3 m，背景介質(zhì)為各向同性介質(zhì)電阻率為ρo=100 Ω·m，在地表進(jìn)行測量設(shè)置40個電極，異常體位于測量區(qū)域的中心部分。

設(shè)置三個模型，mod1異常體為各向同性電阻率為ρ=10 Ω·m，mod1異常體為各向異性電阻率為ρx=ρy=10 Ω·m，ρz=100 Ω·m，mod3異常體為各向異性，電阻率為ρx=ρy=100 Ω·m，ρz=10 Ω·m。

模擬對稱四極裝置的響應(yīng)，取極距AB為3 m、7 m、11 m和21 m的視電阻率曲線進(jìn)行對比，結(jié)果如圖7所示。

圖7中MD-10、MD-10-100、MD-100-10分別代表mod1、mod2和mod3。

圖7 含低阻異常體模型極距為3 m、7 m、11 m和21 m模擬結(jié)果曲線圖

從圖7可以清晰得看出，當(dāng)極距較小時，對稱四極法反映的是較淺介質(zhì)的電性，所以三個模型的結(jié)果相近，均接近100 Ω·m，當(dāng)極距變大時，受到不同異常體的影響，三個模型的視電阻率曲線出現(xiàn)分離。圖7中mod3在四個不同的極距下，視電阻率的值均與背景電阻率較為接近，為100 Ω·m，其異常體的橫向電阻率為100 Ω·m，縱向電阻率為10 Ω·m，而mod2模型的模擬結(jié)果與異常體為各向同性的 的模擬結(jié)果較為接近。對稱四極裝置、偶極偶極裝置以及二極裝置的響應(yīng)，如圖8、圖9和圖10所示。

圖8 對稱四極裝置視電阻率剖面

圖9 偶極偶極裝置視電阻率剖面

圖10 二極裝置視電阻率剖面

通過以上模擬，發(fā)現(xiàn)介質(zhì)橫向電阻率的變化對測量的結(jié)果影響較大，而縱向電阻率對結(jié)果影響相對較小。二極和偶極偶極裝置相較于三極和對稱四極裝置，對縱向電阻率的反映更加靈敏，并且對異常的位置和形態(tài)反映也更加準(zhǔn)確，其中偶極偶極裝置對縱向電阻率的變化最為靈敏。

4 結(jié)論

我們給出了點(diǎn)源2.5維各向異性異常場法的邊值問題和變分問題，并用有限元實(shí)現(xiàn)正演模擬。使用異常場法求解時，假設(shè)點(diǎn)源附近的介質(zhì)為VTI介質(zhì)，簡化空間電位解析解的歐拉角，使得異常場法可以進(jìn)行。

通過算例分析，表明異常場法的精度更高。在模擬計(jì)算中，建議地形模型采用5點(diǎn)波數(shù)使用總場法，平地形模型使用異常場法7點(diǎn)波數(shù)，可以獲得更高的精度。

因?yàn)槭褂卯惓龇〞r，假設(shè)點(diǎn)源附近的介質(zhì)為VTI介質(zhì)，所以文中所提出的簡化方法并不適合點(diǎn)源附近的介質(zhì)為TTI的情況，即介質(zhì)的電性主軸與觀測坐標(biāo)夾角不為0時，該方法并不適合。

通過模擬，發(fā)現(xiàn)對于二維介質(zhì)，直流電阻率法對橫向電阻率的變化更為靈敏，而縱向電阻率對結(jié)果影響相對較小。