任意各向異性介質三維非結構譜元法直流電阻率正演模擬研究

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

吉林大學地球探測科學與技術學院， 長春 130026

0 引言

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

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

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

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

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

1 各向異性模型

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

(1)

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

σ=DσpDT，

(2)

其中，主軸電導率張量σp為

(3)

而σx、σy、σz為主軸電導率.

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

D=DxDyDz，

(4)

其中：

(5)

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

圖1 任意各向異性旋轉示意圖(參考Yin et al., 2016a)(a) 繞x軸旋轉； (b) 繞y軸旋轉； (c) 繞z軸旋轉.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.

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

圖2 各向異性模型及電導率張量(a) VTI介質； (b) HTI介質； (c) 三軸各向異性介質； (d) 傾斜各向異性介質； (e) 任意各向異性介質.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 高階譜元計算方法

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

(6)

其中，σ為大地電導率張量，δ為狄拉克函數，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階正交基函數進行譜展開，即：

(8)

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

Κx=b,

(9)

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

(10)

(11)

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

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

(12)

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

(13)

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

圖3 不同網格四階插值節點分布(a) 傳統譜元法六面體網格節點(每個單元有125個GLL節點)； (b) 將六面體GLL節點直接映射到四面體網格上(頂點位置上產生大量重合的節點)； (c) 四面體網格W&B節點(每個單元僅有25個節點，數量遠小于六面體).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).

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

(14)

式中，我們已將任意電導率張量展開.同時，為了對基函數進行積分，我們將物理域Ωe轉換到參考坐標域Ωs.Nf為參考域內高斯數值積分節點數，wi為權重，J3D為雅克比矩陣，表示了物理坐標系與參考坐標系之間的轉換關系，可表示為

(15)

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

3 算例分析

3.1 精度驗證

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

圖6 各向異性圍巖和異常體主軸電阻率發生旋轉時視電阻率極性圖(a) 圍巖旋轉角1=0°,0°,0°，異常體旋轉角2=0°,0°,0°/30°/45°/60°/90°/135°，接收半徑r=40 m; (b) 圍巖旋轉角1=0°/30°/45°/60°/90°/135°,0°,0°,異常體旋轉角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.

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

圖7 各向異性圍巖和異常體主軸電阻率發生旋轉時視電阻率極性圖(a) 圍巖旋轉角1=0°,0°,135°, 異常體旋轉角2=45°,0°,0°； (b) 圍巖旋轉角1=30°,0°,0°, 異常體旋轉角2=60°,0°,0°； (c)圍巖旋轉角1=0°,0°,30°, 異常體旋轉角2=0°,0°,60°； (d) 圍巖旋轉角1=0°,0°,45°,異常體旋轉角2=135°,0°,0°； (e)圍巖旋轉角1=30°,0°,45°, 異常體旋轉角2=60°,0°,45°；(f)圍巖旋轉角1=45°,0°,30°, 異常體旋轉角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°.

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

圖8 不同網格與不同算法視電阻率極性圖(a) 自適應加密網格的有限元算法加密不同次數； (b) 粗網格譜元法不同階數； (c) 細網格譜元法不同階數.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.

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

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

3.3 復雜地形各向異性響應

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

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

(16)

(17)

即可計算出任意點的電流密度矢量J.由此我們繪制了如圖10所示的xz面上電流密度等值線切片圖.其中，圖10a—d分別對應于圍巖主軸電導率繞y軸旋轉0°、45°、90°、135°的結果.圖中球形異常體的輪廓用白線圈出.觀察不同旋轉軸對應的電流場分布，我們可以看出電流向層理面傾斜的方向集聚，這是由電流向低阻方向聚焦產生的通道效應.隨著旋轉角的增大，電流等值線由垂向逐漸發生傾斜，直至旋轉角為90°時電流等值線變為近乎水平.旋轉角為135°時電流傾斜方向與45°正好相反.由于圍巖各向異性的影響，低阻異常體內部的電流密度也發生了畸變.

圖10 圍巖為傾斜各向異性介質時xz平面上電流密度分布主軸電阻率繞x軸的旋轉角(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°.

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

(18)

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

圖11 各向異性山脊模型聯合剖面裝置視電阻率曲線(模型參見圖8)(a) 圍巖為HTI介質； (b) 圍巖與球體異常均為HTI介質； (c) 球體異常體為HTI介質； (d) 圍巖與球體異常均為各向同性介質； (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 結論

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

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