朱姣， 殷長春， 任秀艷， 劉云鶴， 惠哲劍， 谷宇

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

0 引言

直流電阻率法是電法勘探領(lǐng)域一個重要的分支已在環(huán)境和工程領(lǐng)域獲得廣泛應(yīng)用.該方法利用接地電極向地下供電、建立地下穩(wěn)定電流場，通過觀測測量電極間的電位差以達到推斷和解釋地下目標(biāo)體的目的.在野外勘探中我們常常采用剖面或者測深裝置進行高密度電阻率測量.為此，需要在地面布設(shè)數(shù)十至數(shù)百個測量電極，通過多芯電纜連接到觀測設(shè)備，從而實現(xiàn)多種電極裝置的分布式測量(Loke and Barker, 1996).由于環(huán)境工程領(lǐng)域勘探目標(biāo)大多為三維地質(zhì)體，采用這種陣列測量方法可以準(zhǔn)確地重構(gòu)出三維地質(zhì)體的空間分布，因此高密度電阻率法可獲得豐富的地質(zhì)信息(Loke et al., 2013).

隨著高密度電阻率法在環(huán)境、工程、地下水資源等勘探領(lǐng)域的廣泛應(yīng)用，三維反演和解釋技術(shù)研究受到廣泛關(guān)注(Udphuay et al., 2011; Kneisel et al., 2014; Neyamadpour, 2019; Nthaba et al., 2020).同時，隨著電阻率法向三維精細化探測方向發(fā)展，能夠提取更多的細節(jié)信息，為數(shù)據(jù)解釋打下良好的基礎(chǔ).然而，對精細結(jié)構(gòu)的刻畫將導(dǎo)致龐大的三維數(shù)據(jù)體，對這些三維數(shù)據(jù)體進行反演解釋需要高效和高精度的三維正演模擬技術(shù)(Günther et al., 2006; Rücker et al., 2006).常用的電阻率三維模擬算法包括積分方程法(IE)、有限差分法(FD)和有限元法(FE)(Wang and Signorelli, 2004; Yoon et al., 2016; Ren et al., 2018a).積分方程法僅對異常體進行剖分，計算效率高，但不適合復(fù)雜異常體的數(shù)值計算(Sheng et al., 2006; Ueda and Zhdanov, 2006).有限差分算法簡單、求解速度快，但由于要求計算域被剖分成規(guī)則單元，對復(fù)雜模型的適用性較差(Jaysaval et al., 2014; Davydycheva et al., 2003).有限元法由于網(wǎng)格剖分靈活、計算精度高，是目前最主流的方法(Song et al., 2017；谷宇等，2020).傳統(tǒng)有限元法僅在單元端點或棱邊上產(chǎn)生未知數(shù)，計算精度對網(wǎng)格的依賴性較為嚴重，為達到較高的精度通常需要增加單元內(nèi)未知數(shù)的個數(shù)或者網(wǎng)格數(shù)量(Ren et al., 2013).對于規(guī)則的八叉樹網(wǎng)格，通過引入懸掛點對局部區(qū)域加密，在控制網(wǎng)格數(shù)量的同時提高精度(Grayver and Bürg, 2014).非結(jié)構(gòu)網(wǎng)格及自適應(yīng)網(wǎng)格加密技術(shù)目前已發(fā)展得較為成熟，它能夠根據(jù)后驗誤差不斷迭代得到最優(yōu)的計算網(wǎng)格(Key and Ovall, 2011).然而，其結(jié)果是未知數(shù)大幅增加，并且網(wǎng)格的迭代過程也會帶來時間損耗.為此，人們嘗試通過增加插值基函數(shù)階數(shù)來提高數(shù)值精度，目前在電法勘探領(lǐng)域大多只用到二階插值基函數(shù)(Chen et al., 2018).

本文研究非結(jié)構(gòu)譜元法(Spectral Element Method, SEM)進行電阻率三維各向異性正演模擬問題.譜元法屬于高階有限元算法的一種.傳統(tǒng)高階有限元法采用等間距插值，由于龍格現(xiàn)象的存在，數(shù)值解在靠近端點位置會產(chǎn)生震蕩，因此隨著基函數(shù)階數(shù)的升高，精度不一定會得到提升(Runge,1901).譜元法采用非線性正交基函數(shù)描述場的擴散，在端點附近插值點更加密集，因此能夠有效提高數(shù)值精度(Maday et al., 1993).譜元法結(jié)合了有限元法與譜方法的雙重優(yōu)點.它不僅具有有限元法的網(wǎng)格靈活性，還具有譜方法良好的數(shù)值精度和指數(shù)收斂性，可以通過細化網(wǎng)格或提高插值基函數(shù)階數(shù)達到提升計算精度的目的.Komatitsch 等(2002)將譜元法用于全球地震波場的精確模擬，Huang 等(2019)驗證了譜元法在航空電磁正演模擬中的有效性.然而，所有這些基于譜元法的三維正演均采用六面體網(wǎng)格，這種幾何結(jié)構(gòu)清晰的網(wǎng)格限制了譜元法對于地形及復(fù)雜異常體的刻畫能力，給地球物理探測用于解決復(fù)雜工程問題帶來困難(Yin et al., 2016b; Wang et al., 2013).相比之下，基于非結(jié)構(gòu)四面體網(wǎng)格的譜元法可以很好地解決這些問題.它除了繼承傳統(tǒng)譜元法通過細化網(wǎng)格或提高插值函數(shù)階數(shù)達到很高的計算精度外，還可利用非結(jié)構(gòu)網(wǎng)格模擬起伏地形和地下復(fù)雜構(gòu)造.Liu等(2014)分別利用三角形網(wǎng)格有限元法和譜元法進行二維時域彈性波場模擬，對比結(jié)果驗證了譜元法能夠達到更高的計算精度.Zhu等(2020)將四面體網(wǎng)格譜元法引入電法勘探領(lǐng)域，在與自適應(yīng)加密的有限元算法比較中發(fā)現(xiàn)譜元法只需求解較少的未知數(shù)，就能夠達到更高的精度要求.由于四面體網(wǎng)格坐標(biāo)之間存在的復(fù)雜耦合關(guān)系，無法使用六面體網(wǎng)格對應(yīng)的基函數(shù)和插值節(jié)點(Cohen, 2002)，因此，本文我們將利用Proriol-Koornwinder-Dubiner(PKD)正交多項式(Proriol, 1957; Dubiner, 1991; Koornwinder, 1975)構(gòu)建基函數(shù).另外，我們還使用Warp & Blend插值節(jié)點集，保證良好的插值屬性，改善算法的穩(wěn)定性(Warburton, 2006).

各向異性是環(huán)境與工程問題中經(jīng)常遇到的地質(zhì)現(xiàn)象.從微觀上，各向異性是由于溫度壓力等因素造成礦物結(jié)構(gòu)發(fā)生變化，形成了頁巖和板巖等層理發(fā)育的巖石；從宏觀上，由于構(gòu)造運動造成地下介質(zhì)變質(zhì)分層或者產(chǎn)生裂隙，這些各向同性的薄層相互疊加或者定向排列的構(gòu)造裂隙充水后，導(dǎo)致電阻率呈現(xiàn)各向異性特征(Herwanger et al., 2004; Aissaoui et al., 2019；蔡義宇等，2020).在這些層理或裂隙發(fā)育地區(qū)利用傳統(tǒng)的各向同性模型進行數(shù)據(jù)解釋會產(chǎn)生錯誤結(jié)果(Yin and Weidelt, 1999).因此，利用各向異性模型進行數(shù)據(jù)反演解釋非常必要.雖然針對由介質(zhì)不均勻性造成的各向異性本質(zhì)上可以利用三維精細模型進行模擬，然而受計算條件限制，利用不均勻介質(zhì)模擬這些精細結(jié)構(gòu)難度很大.因此，當(dāng)這些介質(zhì)的層理或裂隙尺寸比采用的電法勘探技術(shù)的有效尺度小很多時，我們可以利用各向異性描述這種不均勻性，從而大大減少計算量.需要注意的是，各向異性與電性不均勻性在數(shù)學(xué)描述上存在本質(zhì)差異.不均勻性指示了巖礦石電阻率隨位置發(fā)生變化，與方向無關(guān)，因此描述各向同性不均勻電導(dǎo)率模型僅需要一個隨位置變化的標(biāo)量函數(shù).然而，各向異性是同一位置上電導(dǎo)率隨電流流動方向發(fā)生變化，此時電流密度與不同方向的電場存在耦合關(guān)系，因此描述各向異性電導(dǎo)率是一個3×3的張量(Yin, 2000).在直流電阻率法各向異性數(shù)值模擬方面，Zhou等(2009)使用高斯正交格網(wǎng)(GQG)方法處理復(fù)雜各向異性介質(zhì)模型，能夠達到與譜方法近似的收斂速度.Li和Spitzer(2005)結(jié)合有限元正演模擬和Cholesky預(yù)處理共軛梯度法，研究了各向異性模型的電阻率分布特征.Ren等(2018b)針對各向異性電阻率模型研究了不同的非結(jié)構(gòu)網(wǎng)格自適應(yīng)加密策略，取得了很好的應(yīng)用效果.本文針對直流電阻率法驗證高階譜元法模擬各向異性模型的有效性.

在本文后續(xù)討論中，我們首先對各向異性模型進行數(shù)學(xué)描述，并以此為基礎(chǔ)建立任意各向異性介質(zhì)電阻率法三維正演理論，進而我們利用正交多項式構(gòu)建譜元法基函數(shù)，并結(jié)合插值節(jié)點與數(shù)值積分節(jié)點完成有限元線性方程中系數(shù)矩陣的建立，最后利用直接求解器實現(xiàn)方程的快速求解.在對本文譜元法進行精度驗證后，我們系統(tǒng)分析各向異性對電阻率響應(yīng)的影響特征.最后，我們重點探討各項異性特征識別方法以及地形效應(yīng)和地形校正技術(shù).

1 各向異性模型

針對任意各向異性模型，我們用一個對稱正定的電導(dǎo)率張量對其進行描述(Yin, 2000)，即：

(1)

其中，σxy=σyx,σxz=σzx,σyz=σzx.電導(dǎo)率張量和電阻率張量之間滿足σ=ρ-1.另外，電導(dǎo)率張量σ可由主軸電導(dǎo)率張量經(jīng)三次歐拉旋轉(zhuǎn)獲得，即：

σ=DσpDT，

(2)

其中，主軸電導(dǎo)率張量σp為

(3)

而σx、σy、σz為主軸電導(dǎo)率.

式(2)中的總旋轉(zhuǎn)矩陣D可表示為

D=DxDyDz，

(4)

其中：

(5)

分別為繞x、y、z軸的旋轉(zhuǎn)矩陣，1、2和3為對應(yīng)的歐拉旋轉(zhuǎn)角.圖1給出任意各向異性介質(zhì)的主軸電導(dǎo)率張量旋轉(zhuǎn)示意圖.其中，實線和虛線分別表示沿不同坐標(biāo)軸歐拉旋轉(zhuǎn)前后的介質(zhì)模型，而陰影面表示同一個面旋轉(zhuǎn)前后的位置.圖1a中(x,y,z)為測量坐標(biāo)系, 主軸電導(dǎo)率繞x軸旋轉(zhuǎn)1之后，坐標(biāo)系變?yōu)?(x1,y1,z1)；參見圖1b，(x1,y1,z1)再繞y軸旋轉(zhuǎn)2，主軸坐標(biāo)系變?yōu)?x2,y2,z2)；同理，(x2,y2,z2)繞z軸旋轉(zhuǎn)3之后，可得最終的主軸坐標(biāo)系(x3,y3,z3)，如圖1c所示.由此完成了沿三個方向的歐拉旋轉(zhuǎn)，此時三個主軸電導(dǎo)率不再沿測量坐標(biāo)系方向，電導(dǎo)率張量變?yōu)?×3的對稱正定矩陣.

圖1 任意各向異性旋轉(zhuǎn)示意圖(參考Yin et al., 2016a)(a) 繞x軸旋轉(zhuǎn)； (b) 繞y軸旋轉(zhuǎn)； (c) 繞z軸旋轉(zhuǎn).Fig.1 Sketch of arbitrary anisotropic rotation(refer to Yin et al., 2016a)(a) Rotation around x axis; (b) Rotation around y axis; (c) Rotation around z axis.

由上述討論可知，我們只需要三個主軸電導(dǎo)率σx、σy、σz和三個歐拉角1、2、3即可對地下各向異性介質(zhì)進行定量描述.參見圖2，按照主軸電導(dǎo)率及定向不同可以將地下介質(zhì)劃分為垂直方向的橫向各向同性(Vertically Transverse Isotropy，VTI)、水平方向的橫向各向同性(Horizontally Transverse Isotropy，HTI)、傾斜各向異性 (Dipping anisotropy)、三軸各向異性(Triple-axis anisotropy)及任意各向異性(Arbitrary anisotropy)(劉云鶴等，2018).VTI介質(zhì)(圖2a)層理面水平，介質(zhì)沿層理面的電導(dǎo)率不隨方向變化，其電導(dǎo)率張量為一個對角陣，且σx=σy≠σz；而HTI介質(zhì)(圖2b)的層理面直立，可由VTI介質(zhì)繞一個水平軸旋轉(zhuǎn)90°得到，其電導(dǎo)率沿直立的層理面同樣不隨方向變化，電導(dǎo)率張量也為一個對角陣，且有σy=σz≠σx或者σx=σz≠σy；三軸各向異性(圖2c)對應(yīng)的電導(dǎo)率張量也是對角陣，但是其三個元素不相同σx≠σz≠σy.傾斜各向異性介質(zhì)(圖2d)是由VTI介質(zhì)繞水平軸旋轉(zhuǎn)一個任意角度，其電導(dǎo)率張量包含5個分量，其中旋轉(zhuǎn)軸對應(yīng)的電導(dǎo)率不隨旋轉(zhuǎn)發(fā)生變化；任意各向異性(圖2e)由VTI介質(zhì)或三軸各向異性繞多個坐標(biāo)軸旋轉(zhuǎn)得到，其電導(dǎo)率張量為3×3的滿秩矩陣.考慮到地球物理中各向異性大多是由層理或定向排列的裂隙引起的，人們通常引入沿層理方向上電導(dǎo)率為σL及垂直層理方向上為σT(滿足σT<σL)，并定義平均電導(dǎo)率和各向異性參數(shù)來描述介質(zhì)的各項異性特征.當(dāng)λ越大時，各項異性差異越明顯.

圖2 各向異性模型及電導(dǎo)率張量(a) VTI介質(zhì)； (b) HTI介質(zhì)； (c) 三軸各向異性介質(zhì)； (d) 傾斜各向異性介質(zhì)； (e) 任意各向異性介質(zhì).Fig.2 Anisotropic model and conductivity tensor(a) VTI medium; (b) HTI medium; (c) Triple-axis anisotropic medium; (d) Dipping anisotropic medium; (e) Arbitrarily anisotropic medium.

2 高階譜元計算方法

為求解三維直流電阻率正演問題，我們首先建立如下坐標(biāo)系：x、y軸位于水平地表，z軸垂直向下.假設(shè)電流強度為I的電流源位于r0=(x0,y0,z0)，則各向異性介質(zhì)中任意位置r=(x,y,z)處的電位U滿足如下定解問題 (Dey and Marrison, 1979)：

(6)

其中，σ為大地電導(dǎo)率張量，δ為狄拉克函數(shù)，n是與邊界相關(guān)的法向量，而：

B=σxx(x-x0)2+2σxy(x-x0)(y-y0)

+2σxz(x-x0)(z-z0)+σyy(y-y0)2

+2σyz(y-y0)(z-z0)+σzz(z-z0)2.

(7)

對于電勢U，我們引入n階正交基函數(shù)進行譜展開，即：

(8)

則根據(jù)伽遼金有限元理論，并選取相同的n階譜元基函數(shù)與權(quán)函數(shù)φi(i=1,…,Nt)，Nt=(n+1)(n+2)(n+3)/6為每個四面體單元的節(jié)點數(shù)，我們可將公式(6)轉(zhuǎn)化為積分弱形式后得到如下線性方程組：

Κx=b,

(9)

其中，左端項包括全局電勢x及系數(shù)矩陣K，右端項b為源項.在每個單元上我們有：

(10)

(11)

其中，Ωe表示四面體單元，Γe表示單元邊界.公式(10)中的第二項僅在包含外邊界的四面體單元中出現(xiàn).

在每個單元內(nèi)，φ=[φ1,…,φNt]是一組正交完備的多項式基，可表示為

(12)

其中，r1,r2,…rNt表示插值節(jié)點，V表示對稱正定的范德蒙德矩陣，而內(nèi)部元素ψj(ri)為[-1,1]正交區(qū)間內(nèi)的Proriol-Koornwinder-Dubiner(PKD)多項式，即：

(13)

在正演模擬中，合理的插值節(jié)點分布是避免數(shù)值解震蕩的關(guān)鍵問題所在.前人研究發(fā)現(xiàn)Gauss-Lobatto-Legendre(GLL)點集是最優(yōu)的六面體網(wǎng)格插值節(jié)點(Pasquetti and Rapetti, 2010, 圖3a)，然而對于四面體網(wǎng)格，簡單地將GLL節(jié)點映射到四面體單元中會在頂點位置產(chǎn)生大量重合的節(jié)點(圖3b)，因此我們需要在四面體中構(gòu)建合理的插值點集.本文選用具有良好插值屬性的Warp & Blend點集，其主要思想是將等距節(jié)點拉伸到最優(yōu)位置，使得節(jié)點分布具有向端點聚攏的特性(Warburton, 2006).圖3c展示了4階譜插值節(jié)點分布.

圖3 不同網(wǎng)格四階插值節(jié)點分布(a) 傳統(tǒng)譜元法六面體網(wǎng)格節(jié)點(每個單元有125個GLL節(jié)點)； (b) 將六面體GLL節(jié)點直接映射到四面體網(wǎng)格上(頂點位置上產(chǎn)生大量重合的節(jié)點)； (c) 四面體網(wǎng)格W&B節(jié)點(每個單元僅有25個節(jié)點，數(shù)量遠小于六面體).Fig.3 Distribution of fourth-order interpolation nodes on different grids(a) Traditional spectral element method with hexahedral mesh (each unit has 125 GLL nodes); (b) Hexahedral GLL nodes are directly mapped to the tetrahedral mesh (a large number of overlapping nodes are generated at vertex positions); (c) W&B nodes in tetrahedral mesh (each element has only 25 nodes, much less than that of hexahedron).

在得到插值節(jié)點和譜插值基函數(shù)后，我們可將公式(10)中每個單元的系數(shù)矩陣計算出來.需要強調(diào)的是，由于基函數(shù)梯度的三重積分項沒有解析形式，故采用高斯數(shù)值積分計算，即：

(14)

式中，我們已將任意電導(dǎo)率張量展開.同時，為了對基函數(shù)進行積分，我們將物理域Ωe轉(zhuǎn)換到參考坐標(biāo)域Ωs.Nf為參考域內(nèi)高斯數(shù)值積分節(jié)點數(shù)，wi為權(quán)重，J3D為雅克比矩陣，表示了物理坐標(biāo)系與參考坐標(biāo)系之間的轉(zhuǎn)換關(guān)系，可表示為

(15)

綜上所述，本文針對任意各向異性介質(zhì)的電導(dǎo)率張量模型，從直流電阻率法滿足的泊松方程出發(fā)，采用靈活的四面體網(wǎng)格對復(fù)雜模型進行精細剖分，結(jié)合高精度的譜插值基函數(shù)與插值節(jié)點建立了伽遼金法有限元控制方程，使用多波前直接求解器MUMPS進行線性方程求解(Amestoy et al.，2000)，并利用公式(8)中的插值方法計算任意接收點處的電勢.

3 算例分析

3.1 精度驗證

圖4 測線沿y方向的視電阻率及與解析解的相對誤差(圖中的n代表譜插值基函數(shù)階數(shù))Fig.4 Apparent resistivities and relative errors of analytical solution along survey line in y direction (n represents the order of the spectral interpolation basis function)

3.2 各向異性識別與效率分析

圖5 立方體異常各向異性模型及測量裝置示意圖Fig.5 Sketch of anisotropic cube model and survey configuration

我們采用圖5中的二級裝置研究兩組不同收發(fā)距 (AM1：r=40 m和AM2：r=150 m)條件下的視電阻變化特征.為獲得視電阻率極性圖，我們進行36組測量：發(fā)射源點位置A固定，接收點繞A點依次順時針旋轉(zhuǎn)10°，完成一個閉合環(huán)路觀測.結(jié)合二極裝置極距與探測深度的關(guān)系，我們知道當(dāng)極距r=40 m時，能夠探測到異常體的電阻率分布，而當(dāng)r=150 m時，能夠探測到圍巖的電阻率.圖6a給出當(dāng)異常體電阻率繞z軸旋轉(zhuǎn)0°、30°、45°、60°、90°、135°時的視電阻率極性圖(圖中從原點到極性圖端點的長度表征了視電阻率大小，而其方向表征了二極裝置的定向)，圖6b給出圍巖電阻率繞x軸旋轉(zhuǎn)0°、30°、45°、60°、90°、135°的視電阻率極性圖.由圖可以看出，當(dāng)r=40m時，隨著沿z軸方向旋轉(zhuǎn)角度的變化，視電阻率曲線形態(tài)基本不變，均為軸對稱圖形，并且有兩個正交的對稱軸.然而，曲線繞中心發(fā)生了旋轉(zhuǎn)，其中長軸方向角與異常體的旋轉(zhuǎn)角吻合，即視電阻率極性圖中長軸方向與走向方向一致.當(dāng)接收半徑r=150 m時，隨著沿x軸方向旋轉(zhuǎn)角度的變化視電阻率曲線呈現(xiàn)明顯的不對稱.長軸的方向基本不變，但是短軸的一側(cè)沿著y軸正方向逐漸被拉伸.當(dāng)旋轉(zhuǎn)角達到90°時曲線為正圓形.當(dāng)旋轉(zhuǎn)角為45°和135°時視電阻率極性圖關(guān)于x軸對稱，但兩者短軸的拉伸方向發(fā)生了翻轉(zhuǎn).另外，從圖可以看出短軸的大頭指示了地層傾向，這與圍巖電阻率的變化規(guī)律相反，呈現(xiàn)電各向異性反?，F(xiàn)象.由此可以得出結(jié)論，利用視電阻率極性圖我們可以確定各向異性主軸方向與傾向，為未來各向異性電阻率數(shù)據(jù)三維反演提供重要的旋轉(zhuǎn)角度參數(shù)信息，減少多解性.

圖6 各向異性圍巖和異常體主軸電阻率發(fā)生旋轉(zhuǎn)時視電阻率極性圖(a) 圍巖旋轉(zhuǎn)角1=0°,0°,0°，異常體旋轉(zhuǎn)角2=0°,0°,0°/30°/45°/60°/90°/135°，接收半徑r=40 m; (b) 圍巖旋轉(zhuǎn)角1=0°/30°/45°/60°/90°/135°,0°,0°,異常體旋轉(zhuǎn)角2=0°,0°,0°，接收半徑r=150 m.Fig.6 Polar plots of apparent resistivity when principal resistivity axes of anomaly body and anisotropic host rock rotate(a) Rotation angle of host rock 1=0°,0°,0°, rotation angle of abnormal body 2=0°,0°,0°/30°/45°/60°/90°/135°, receiving radius is r=40 m; (b) Rotation angle of host rock 1=0°/30°/45°/60°/90°/135°,0°,0°, rotation angle of abnormal body 2=0°,0°,0°, receiving radius is r=150 m.

下面進一步研究多個各向異性參數(shù)變化時視電阻率極性圖分布特征.圖7中我們給出兩組極距的模擬結(jié)果.對比圖7a、d中小極距極性圖發(fā)現(xiàn)，異常體繞x軸旋轉(zhuǎn)45°和135°時，兩者曲線關(guān)于x軸對稱，大頭分別指向y軸負方向和正方向.大極距視電阻率極性圖長軸對應(yīng)的方位角與背景電阻率繞z軸旋轉(zhuǎn)角基本相同(分別為135°和45°)，兩支曲線關(guān)于y軸對稱.由于受異常體各向異性的影響，大極距曲線在立方體傾向方向也發(fā)生了一定的拉伸.圖7b中圍巖與異常體分別繞x軸旋轉(zhuǎn)30°和60°，大極距與小極距對應(yīng)視電阻率極性圖的大頭方向相反，分別指向y軸正、負方向，而圖7e中除繞x方向旋轉(zhuǎn)以外，兩者均繞z軸再旋轉(zhuǎn)45°，因此兩支曲線長軸也發(fā)生了旋轉(zhuǎn)，旋轉(zhuǎn)角接近45°，并且曲線大頭的方向也發(fā)生了改變，大致指示走向與傾向的變化.同樣地，我們將模型主軸電導(dǎo)率先繞z軸旋轉(zhuǎn)30°和60°(圖7c)，再繞x軸旋轉(zhuǎn)45°(圖7f)后計算視電阻率極性圖.從圖7c可以看出兩支曲線分別關(guān)于30°和60°呈現(xiàn)軸對稱，而從圖7f可以看出兩只曲線均發(fā)生了拉伸，但拉伸方向相反.由此，我們可以得出結(jié)論：當(dāng)?shù)叵陆橘|(zhì)各向異性屬性較為復(fù)雜時，通過分析視電阻率極性圖的分布形態(tài)，我們?nèi)钥纱_定各向異性主軸的走向和傾向.這些參數(shù)的確定對未來電阻率數(shù)據(jù)三維反演中設(shè)定合理反演參數(shù)、減少各向異性數(shù)據(jù)反演多解性至關(guān)重要.

圖7 各向異性圍巖和異常體主軸電阻率發(fā)生旋轉(zhuǎn)時視電阻率極性圖(a) 圍巖旋轉(zhuǎn)角1=0°,0°,135°, 異常體旋轉(zhuǎn)角2=45°,0°,0°； (b) 圍巖旋轉(zhuǎn)角1=30°,0°,0°, 異常體旋轉(zhuǎn)角2=60°,0°,0°； (c)圍巖旋轉(zhuǎn)角1=0°,0°,30°, 異常體旋轉(zhuǎn)角2=0°,0°,60°； (d) 圍巖旋轉(zhuǎn)角1=0°,0°,45°,異常體旋轉(zhuǎn)角2=135°,0°,0°； (e)圍巖旋轉(zhuǎn)角1=30°,0°,45°, 異常體旋轉(zhuǎn)角2=60°,0°,45°；(f)圍巖旋轉(zhuǎn)角1=45°,0°,30°, 異常體旋轉(zhuǎn)角2=45°,0°,60°.Fig.7 Polar plots of apparent resistivities when the principal resistivity axes of anomaly body and anisotropic host rock rotate(a) Rotation angle of host rock 1=0°,0°,135°, rotation angle of abnormal body 2=45°,0°,0°; (b) Rotation angle of host rock 1=30°,0°,0°, rotation angle of abnormal body 2=60°,0°,0°; (c) Rotation angle of host rock 1=0°,0°,30°, rotation angle of abnormal body 2=0°,0°,60°; (d) Rotation angle of host rock 1=0°,0°,45°, rotation angle of abnormal body 2=135°,0°,0°; (e) Rotation angle of host rock 1=30°,0°,45°, rotation angle of abnormal body 2=60°, 0°, 45°; (f) Rotation angle of host rock 1=45°,0°,30°, rotation angle of abnormal body 2=45°,0°,60°.

為檢驗本文算法的計算效率，我們還對比了針對不同網(wǎng)格的譜元和有限元算法的計算結(jié)果.為此，我們將圖5中的立方體主軸電阻率繞x軸旋轉(zhuǎn)90°，并在圖8中給出了三種不同情況的視電阻率極性圖.圖8a為多次自適應(yīng)網(wǎng)格加密的有限元計算結(jié)果 (Ren and Tang, 2010，開源代碼https:∥software.seg.org/2010/0002/index.html).其中，Ap_0為初始網(wǎng)格極性圖，Ap_1、Ap_2、Ap_3為網(wǎng)格自適應(yīng)加密1、2及3次的極性圖.可以看出隨著網(wǎng)格加密極性圖越來越光滑，視電阻率計算精度越來越高.圖8b、c分別給出粗、細網(wǎng)格上不同階數(shù)譜元法的極性圖.有圖可以看出，當(dāng)譜元基函數(shù)階數(shù)為n=2時，粗、細網(wǎng)格上的視電阻率極性圖均很光滑.進一步對比圖8b、c可以發(fā)現(xiàn)，在細網(wǎng)格上n=3與n=4的視電阻率計算結(jié)果已很好地收斂.

圖8 不同網(wǎng)格與不同算法視電阻率極性圖(a) 自適應(yīng)加密網(wǎng)格的有限元算法加密不同次數(shù)； (b) 粗網(wǎng)格譜元法不同階數(shù)； (c) 細網(wǎng)格譜元法不同階數(shù).Fig.8 Polar plots of apparent resistivities with different grids and different algorithms(a) Adaptive finite element algorithm with different grid refinement； (b) Coarse grid spectral element method with different orders； (c) Fine grid spectral element method with different orders.

如前所述，加密網(wǎng)格與提高階數(shù)均能提升譜元法數(shù)值模擬精度.本節(jié)我們對不同計算方法之間的計算效率進行對比.以細網(wǎng)格4階譜元模擬結(jié)果作為參考, 我們計算針對不同方法及不同網(wǎng)格模擬結(jié)果的相對誤差.表1給出了對比結(jié)果.圖8a、b、c分別對應(yīng)表中網(wǎng)格加密的有限元法、粗網(wǎng)格譜元法和細網(wǎng)格譜元法.由表可以看出，通過簡單加密網(wǎng)格有限元法相對誤差逐漸減小，但多次網(wǎng)格加密之后的相對誤差仍然較大.相比之下，譜元法能夠在保持初始網(wǎng)格不變的條件下，通過提高階數(shù)使得模擬精度得到明顯的提升.自適應(yīng)加密一次網(wǎng)格(AP_1)與2階插值(n=2)未知數(shù)數(shù)量相近，但譜元法的相對誤差遠小于有限元法，在計算時間上也略少于后者.對比網(wǎng)格自適應(yīng)加密三次(AP_3)有限元與4階譜元(n=4)的計算結(jié)果可知，譜元法只需求解較少的未知數(shù)就能達到更高的精度.進一步對比粗、細網(wǎng)格下不同階次譜元法的網(wǎng)格數(shù)、自由度、求解時間、計算誤差等參數(shù)發(fā)現(xiàn)，采用更加精細的網(wǎng)格比粗網(wǎng)格能夠獲得更高的計算精度.由此可以得出結(jié)論：基于譜元法我們可以通過采用雙重策略獲得高精度的數(shù)值計算結(jié)果.

表1 自適應(yīng)加密有限元法與譜元法參數(shù)對比Table 1 Forward modelling data comparison between adaptive finite element method and spectral element method

3.3 復(fù)雜地形各向異性響應(yīng)

為檢驗四面體網(wǎng)格譜元法對于復(fù)雜各向異性模型的適應(yīng)性，我們設(shè)計了一個如圖9所示的高度差為20 m，直徑為100 m的三維山脊地形模型.山脊頂部中心點的坐標(biāo)為(0,0,-20).在地形下方10 m處埋藏一個半徑為15 m的球體，模型區(qū)域為3000 m×3000 m×3000 m，共有165793個四面體網(wǎng)格.測線沿x軸布設(shè)，間距為2m.由于在源點附近、電阻率分界面及地表面等區(qū)域會產(chǎn)生急劇變化的場，為保證算法穩(wěn)定，我們對這些區(qū)域進行了局部網(wǎng)格加密(圖9b).

圖9 三維山脊地形模型及其剖分網(wǎng)格Fig.9 3D ridge terrain model with a sphere embedded and grid

(16)

(17)

即可計算出任意點的電流密度矢量J.由此我們繪制了如圖10所示的xz面上電流密度等值線切片圖.其中，圖10a—d分別對應(yīng)于圍巖主軸電導(dǎo)率繞y軸旋轉(zhuǎn)0°、45°、90°、135°的結(jié)果.圖中球形異常體的輪廓用白線圈出.觀察不同旋轉(zhuǎn)軸對應(yīng)的電流場分布，我們可以看出電流向?qū)永砻鎯A斜的方向集聚，這是由電流向低阻方向聚焦產(chǎn)生的通道效應(yīng).隨著旋轉(zhuǎn)角的增大，電流等值線由垂向逐漸發(fā)生傾斜，直至旋轉(zhuǎn)角為90°時電流等值線變?yōu)榻跛?旋轉(zhuǎn)角為135°時電流傾斜方向與45°正好相反.由于圍巖各向異性的影響，低阻異常體內(nèi)部的電流密度也發(fā)生了畸變.

圖10 圍巖為傾斜各向異性介質(zhì)時xz平面上電流密度分布主軸電阻率繞x軸的旋轉(zhuǎn)角(a)1=0°,0°,0°； (b) 1=0°,45°,0°； (c) 1=0°,90°,0°； (d) 1=0°,135°,0°.Fig.10 Current density distribution on the xz- plane when host rock is dipping anisotropicThe rotation angle of theprincipal resistivity axis around the x axis are(a)1=0°,0°,0°； (b) 1=0°,45°,0°； (c) 1=0°,90°,0°； (d) 1=0°,135°,0°.

針對聯(lián)合剖面裝置形式，我們計算了上述四個模型的視電阻率剖面曲線，并與純地形(不存在異常體)的視電阻率曲線進行對比，如圖11a—d所示.每組聯(lián)合剖面測量產(chǎn)生A、B兩支曲線，每對曲線均有三個公共交點，其中山脊頂部對應(yīng)高阻交點，在兩側(cè)山腳處出現(xiàn)低阻交點.由于地形效應(yīng)遠大于異常體響應(yīng)，為消除地形效應(yīng)我們利用如下地形校正公式(Fox et al., 1980)：

(18)

利用式(18)，我們對上述四個模型的視電阻率進行地形校正，結(jié)果如圖11e—h所示.從圖可以看出，四組校正后的數(shù)據(jù)均很好地展示球形良導(dǎo)體上的聯(lián)合剖面曲線.曲線左右支交點均位于球體中心點位置(x=0).各向異性圍巖在地形校正后的視電阻率值遠大于各向異性異常體的響應(yīng).對比圖11e、f和圖11g、h，我們還發(fā)現(xiàn)各向異性的球體使得視電阻率的極小值更小.必須指出，由于地形和各向異性之間的復(fù)雜耦合關(guān)系，導(dǎo)致地下電流場分布異常復(fù)雜(參見圖10)，因此簡單應(yīng)用地形校正無法完全消除地形效應(yīng)的影響.然而，從圖11 給出的算例可以看出，即使大地存在起伏地形和復(fù)雜各向異性，視電阻率響應(yīng)經(jīng)地形校正后仍能很好地確定異常體位置與導(dǎo)電性特征.

圖11 各向異性山脊模型聯(lián)合剖面裝置視電阻率曲線(模型參見圖8)(a) 圍巖為HTI介質(zhì)； (b) 圍巖與球體異常均為HTI介質(zhì)； (c) 球體異常體為HTI介質(zhì)； (d) 圍巖與球體異常均為各向同性介質(zhì)； (e)—(h) 地形校正后視電阻率曲線.Fig.11 Apparent resistivity of the composite profile device for the anisotropic ridge model shown in Fig.8(a) Host rock is HTI； (b) Host rock and sphere are both HTI； (c) Sphere is HTI； (d) Host rock and sphere are isotropic； (e)—(h) Apparent resistivity curves after terrain correction.

4 結(jié)論

本文將四面體網(wǎng)格的譜元法成功應(yīng)用于各向異性介質(zhì)直流電阻率法正演模擬之中，并通過層狀模型的解析解驗證了算法的正確性.相比于傳統(tǒng)有限元算法，譜元法在提高數(shù)值模擬精度上具有獨特優(yōu)勢.它可以靈活利用加密網(wǎng)格或提高譜插值基函數(shù)階數(shù)，改善計算精度.通過對各向異性模型模擬及響應(yīng)特征分析發(fā)現(xiàn)，利用視電阻率極性圖的長短軸走向以及極性圖不對稱性，可以很好地確定地下介質(zhì)各向異性主軸電阻率及其走向和傾向.針對復(fù)雜地形條件下各向異性大地模型，采用傳統(tǒng)的地形校正方法仍然可以有效地消除地形效應(yīng)，獲得地下導(dǎo)電體的可靠信息.這將有助于對各向異性環(huán)境下地下三維電性結(jié)構(gòu)進行精確反演和解釋.

致謝感謝中南大學(xué)邱樂穩(wěn)博士、陳煌博士以及吉林大學(xué)殷長春教授研究團隊其他成員對本文研究提供的幫助，感謝兩位審稿人提出的寶貴意見.